类欧几里得算法&洛谷P5170题解

类欧几里得算法&洛谷P5170题解

好久没更博客了说，不能这样颓废了(flag

Problem

给定\(n,a,b,c\),求\(f(a,b,c,n)=\sum\limits_{i=0}^n\lfloor\frac{ai+b}c\rfloor,g(a,b,c,n)=\sum\limits_{i=0}^n\lfloor\frac{ai+b}c\rfloor^2,h(a,b,c,n)=\sum\limits_{i=0}^ni\lfloor\frac{ai+b}c\rfloor\)

Algorithm

先考虑看上去最可做的第一个柿子。

首先可以通过一些trivial的讨论把\(a,b\)对\(c\)取模。

\[\begin{align} f(a,b,c,n)&=\sum\limits_{i=0}^n\lfloor\frac{ai+b}c\rfloor\\ &=\sum\limits_{i=0}^n\lfloor\frac{(a\%c)i+b\%c}c\rfloor+\frac{(n+1)*n}2\lfloor\frac ac\rfloor+(n+1)\lfloor\frac bc\rfloor\\ &=f(a\%c,b\%c,c,n)+\frac{(n+1)*n}2\lfloor\frac ac\rfloor+(n+1)\lfloor\frac bc\rfloor \end{align} \]

我们接下来考虑\(a,b<c\)的情况,考虑此时的\(\lfloor\frac{ai+b}c\rfloor\)取值不会超过\(n\),可以枚举取值。

\[\begin{align} f(a,b,c,n) &=\sum\limits_{i=0}^n\lfloor\frac{ai+b}c\rfloor\\ &=\sum\limits_{i=0}^n\sum\limits_{d=1}^{\lfloor\frac{an+b}{c}\rfloor}[\lfloor\frac{ai+b}c\rfloor\ge d]\\ 两边大力乘上a^{-1}c &=\sum\limits_{i=0}^n\sum\limits_{d=1}^{\lfloor\frac{an+b}{c}\rfloor}[a^{-1}c\lfloor\frac{ai+b}c\rfloor\ge a^{-1}cd\gt a^{-1}(cd-1)]\\ &=\sum\limits_{i=0}^n\sum\limits_{d=1}^{\lfloor\frac{an+b}{c}\rfloor}[i\gt \frac{(cd-b-1)}a]\\ &=\sum\limits_{i=0}^n\sum\limits_{d=0}^{\lfloor\frac{an+b}{c}\rfloor-1}[i\gt \frac{(cd+c-b-1)}a]\\ &=\sum\limits_{d=0}^{\lfloor\frac{an+b}{c}\rfloor-1}(n-\lfloor\frac{(cd+c-b-1)}a\rfloor)\\ &=n\lfloor\frac{an+b}{c}\rfloor-\sum\limits_{d=0}^{\lfloor\frac{an+b}{c}\rfloor-1}\lfloor\frac{(cd+c-b-1)}a\rfloor\\ &=n\lfloor\frac{an+b}{c}\rfloor-f(c,c-b-1,a,\lfloor\frac{an+b}{c}\rfloor-1)\\ \end{align} \]

于是我们就可以在\(log(n)\)的渐进时间复杂度下求出\(f\)了

好的，相信\(g,h\)不是问题。

于是我们看看\(g\)

\[g(a,b,c,n)=\sum\limits_{i=0}^n\lfloor\frac{ai+b}c\rfloor^2 \]

套路一波。

\[\begin{aligned}g(a,b,c,n)&=\sum\limits_{i=0}^n(\lfloor\frac{(a\%c)i+b\%c}c\rfloor+i\lfloor\frac ac\rfloor+\lfloor\frac bc\rfloor)^2\\大力展开,上式&=\sum\limits_{i=0}^n(\lfloor\frac{(a\%c)i+b\%c}c\rfloor^2+i^2\lfloor\frac ac\rfloor^2+\lfloor\frac bc\rfloor^2+2i\lfloor\frac ac\rfloor\lfloor\frac{(a\%c)i+b\%c}c\rfloor+2i\lfloor\frac bc\rfloor\lfloor\frac{(a\%c)i+b\%c}c\rfloor+2i\lfloor\frac ac\rfloor\lfloor\frac bc\rfloor)^2\\设S1(n)&=\sum\limits_{i=0}^ni=\frac{i*(i+1)}{2},S2(n)=\sum\limits_{i=0}^ni^2=\frac {i*(i+1)*(2i+1)}6\\整理,原式&=g(a\%c,b\%c,c,n)+S2(n)\lfloor\frac ac\rfloor^2+(n+1)\lfloor\frac bc\rfloor^2+2\lfloor\frac ac\rfloor h(a\%c,b\%c,c,n)+2\lfloor\frac bc\rfloor f(a\%c,b\%c,c,n)+2S1(n)\lfloor\frac ac\rfloor\lfloor\frac bc\rfloor\end{aligned} \]

你发现事情变得辣手了起来，我们似乎并不知道这个\(h\)怎么搞。

但是你先不管它,假设你知道了\(h\)怎么算。

\[\begin{aligned}g(a,b,c,n)&=\sum\limits_{i=0}^n\lfloor\frac{ai+b}c\rfloor^2\\&=\sum\limits_{i=0}^n(2\sum\limits_{d=1}^{\lfloor\frac{ai+b}c\rfloor}d-\lfloor\frac{ai+b}c\rfloor)\\&=\sum\limits_{i=0}^n(2\sum\limits_{d=0}^{\lfloor\frac{ai+b}c\rfloor-1}(d+1)-\lfloor\frac{ai+b}c\rfloor)\\&=2\sum_{d=0}^{\lfloor\frac {an+b}c\rfloor-1}(d+1)\sum\limits_{i=0}^n[\lfloor\frac{ai+b}c\rfloor\ge d+1]-f(a,b,c,n)\\&=2\sum_{d=0}^{\lfloor\frac {an+b}c\rfloor-1}(d+1)(n-\lfloor\frac{cd+c-b-1}a\rfloor)-f(a,b,c,n)\\&=n*\lfloor\frac {an+b}c\rfloor*(\lfloor\frac {an+b}c\rfloor+1)-2\sum_{d=0}^{\lfloor\frac {an+b}c\rfloor-1}(d+1)\lfloor\frac{cd+c-b-1}a\rfloor-f(a,b,c,n)\\&=n*\lfloor\frac {an+b}c\rfloor*(\lfloor\frac {an+b}c\rfloor+1)-2h(c,c-b-1,a,\lfloor\frac{an+b}c\rfloor-1)-2f(c,c-b-1,a,\lfloor\frac{an+b}c\rfloor-1)-f(a,b,c,n)\\\end{aligned} \]

我们最后考虑一下\(h\),首先大力展开

\[\begin{aligned}h(a,b,c,n)=S2(n)\lfloor\frac ac\rfloor+S1(n)\lfloor\frac bc\rfloor+h(a\%c,b\%c,c,n)\end{aligned} \]

好的,继续

\[\begin{aligned} h(a,b,c,n)&=\sum\limits_{i=0}^ni\lfloor\frac {ai+b}c\rfloor\\ &=\sum\limits_{i=0}^ni\sum\limits_{d=1}^{\lfloor\frac{an+b}{c}\rfloor}[a^{-1}c\lfloor\frac{ai+b}c\rfloor\ge a^{-1}cd\gt a^{-1}(cd-1)]\\ &=\sum\limits_{i=0}^ni\sum\limits_{d=0}^{\lfloor\frac{an+b}{c}\rfloor-1}[i\gt\frac{cd+c-b-1}{a}]\\ &=S1(n)\lfloor\frac {an+b}c\rfloor-\sum\limits_{d=0}^{\lfloor\frac{an+b}{c}\rfloor-1}\sum\limits_{i=0}^ni[\frac{cd+c-b-1}a\ge i]\\ &=S1(n)\lfloor\frac {an+b}c\rfloor-\sum\limits_{d=0}^{\lfloor\frac{an+b}{c}\rfloor-1}\frac{(\frac{cd+c-b-1}a)*(\frac{cd+c-b-1}a+1)}2\\ &=S1(n)\lfloor\frac {an+b}c\rfloor-\frac12f(c,c-b-1,a,\lfloor\frac{an+b}{c}\rfloor-1)-\frac12\\g(c,c-b-1,a,\lfloor\frac{an+b}{c}\rfloor-1) \end{aligned} \]

如此递归，我们就可以在\(O(logn)\)的渐进时间复杂度内把\(f,g,h\)都求出来了。

代码

/*
@Date    : 2020-02-13 08:16:26
@Author  : Adscn (adscn@qq.com)
@Link    : https://www.cnblogs.com/LLCSBlog
*/
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define int ll
#define IL inline
#define RG register
#define gi geti<int>()
#define gl geti<ll>()
#define gc getchar()
#define File(a) freopen(a".in","r",stdin);freopen(a".out","w",stdout)
template<typename T>IL bool chkmax(T &x,const T &y){return x<y?x=y,1:0;}
template<typename T>IL bool chkmin(T &x,const T &y){return x>y?x=y,1:0;}
template<typename T>
IL T geti()
{
	RG T xi=0;
	RG char ch=gc;
	bool f=0;
	while(!isdigit(ch))ch=='-'?f=1:f,ch=gc;
	while(isdigit(ch))xi=xi*10+ch-48,ch=gc;
	return f?-xi:xi;
}
template<typename T>
IL void pi(T k,char ch=0)
{
	if(k<0)k=-k,putchar('-');
	if(k>=10)pi(k/10);
	putchar(k%10+'0');
	if(ch)putchar(ch);
}
const int P=998244353,i2=(P+1)>>1;
inline void exgcd(int a,int b,int &x,int &y){
	if(!b)return x=1,y=0,void(); 
	exgcd(b,a%b,y,x);
	y=(y-a/b*x)%P;
}
inline int inv(int a){
	int x,y;
	exgcd(a,P,x,y);
	return (x+P)%P;
}
const int i6=inv(6);
struct data{ll a,b,c;};
inline data solve(int a,int b,int c,int n)
{
	ll n1=(n+1)%P,S1=n1*n%P*i2%P,S2=n*n1%P*(n+n1)%P*i6%P;
	ll m=(1ll*a*n+b)/c,ac=a/c,bc=b/c,acq=ac*ac%P,bcq=bc*bc%P,nm=n*m%P;
	if(!a)return (data){n1*bc%P,n1*bcq%P,bc*S1%P};
	if(a>=c||b>=c){
		data ret=solve(a%c,b%c,c,n);
		return (data){
			(ret.a+S1*ac%P +n1*bc%P)%P,
			(ret.b
				+S2*acq%P
				+n1*bcq%P
				+2*ac*ret.c%P
				+2*bc*ret.a%P
				+2*S1*ac%P*bc%P)%P,
			(ret.c+S2*ac%P +S1*bc%P)%P};
	}
	auto lst=solve(c,c-b-1,a,m-1);
	data ret;
	ret.a=(nm-lst.a)%P;
	ret.b=(nm*(m+1)%P-2ll*(lst.a+lst.c)%P-ret.a)%P;
	ret.c=(S1*m%P-1ll*i2*(lst.a+lst.b)%P)%P;
	return ret;
}
signed main(void)
{
	#ifndef ONLINE_JUDGE
//	File("");
	#endif
	int T=gi;
	while(T--)
	{
		int n=gi,a=gi,b=gi,c=gi;
		auto ret=solve(a,b,c,n);
		printf("%lld %lld %lld\n",(ret.a+P)%P,(ret.b+P)%P,(ret.c+P)%P);
	}
	return 0;
}