LOJ138 类欧几里得算法

类欧几里得算法

给出 (T) 组询问，每组用 (n, a, b, c, k_1, k_2) 来描述。对于每组询问，请你求出

[sum_{x = 0} ^ {n} x ^ {k_1} {left lfloor frac{ax + b}{c} ight floor} ^ {k_2} ]

对 (1000000007) 取模。

对于 (100 \%) 的数据，(T = 1000, 1 le n, a, b, c le {10} ^ 9, 0 le k_1 + k_2 le 10) 。

数形结合

这里只计算

[sum_{i=0}^{N}lfloorfrac{ai+b}{c} floor ]

首先通过简单转化把 (a,b) 都放到 (mod c) 的意义下。

考虑计算式 (sum_{i=0}^{N}lfloorfrac{ai+b}{c} floor) 的含义。它要求的其实是 (y=frac{a}{c}x+frac{b}{c}~(0 leq x leq N)) 这条线段下方（包括线段上）的整点个数。要求整点 ((xgeq 0,y > 0))。

使用差分，转化为用 ((0,0) ightarrow (N,frac{aN+b}{c})) 这个矩形里的整点个数减去 (y=frac{a}{c}x+frac{b}{c}~(0 leq x leq N)) 这条线段上方（不包括线段上）的整点个数。要求整点 ((xgeq 0,y > 0))。

然后翻转坐标轴， 转化成求 (y=frac{c}{a}x-frac{b}{a}~(0 < x leq frac{aN+b}{c})) 这条线段下方（不包括线段上）的整点个数。要求整点 ((x>0,ygeq 0))。

化归原问题，求 (y=frac{c}{a}(x+1)-frac{b}{a}-frac{1}{a}~(0 leq x < frac{aN+b}{c})) 这条线段下方（包括线段上）的整点个数。要求整点 ((xgeq 0,y geq 0))。(y=0) 的那些点要特殊处理掉，由于 (c-b-1geq 0)，所以每个可能的横坐标都对应了一个这样的点。

此时相当于求

[(N+1)lfloorfrac{aN+b}{c} floor - sum_{i=0}^{lfloorfrac{aN+b}{c} floor-1} (lfloorfrac{ci+c-b-1}{a} floor+1) ]

此时 (c>a)，所以可以再次取模，进行递归。注意到在递归中的参数变化：((a,c) ightarrow (amod c,c) ightarrow (c,amod c))，所以递归层数为 (O(log V))。

推公式

https://blog.csdn.net/qq_39972971/article/details/94618394

是哪个天才发明的把 (x^k) 变成 (sum_{i=0}^{x-1}(i+1)^k-i^k)？

以下考虑实现函数 (func(N,a,b,c)) ，计算 (0leq k_1+k_2leq10) 的情况下所求式子的值，即

[sum_{i=0}^{N}i^{k_1}lfloorfrac{ai+b}{c} floor^{k_2} ]

递归基

若 (a=0) 或 (lfloorfrac{aN+b}{c} floor=0) ，那么， (lfloorfrac{ai+b}{c} floor) 的取值始终相同，答案即为

[lfloorfrac{aN+b}{c} floorsum_{i=0}^{N}i^{k_1} ]

可以用插值解决。

(a geq c)

若 (ageq c) ，则令 (q=lfloorfrac{a}{c} floor,r=amod c) ，所求式子即为

[sum_{i=0}^{N}i^{k_1}(qi+lfloorfrac{ri+b}{c} floor)^{k_2}\ =sum_{i=0}^{N}i^{k_1}sum_{j=0}^{k_2}inom{k_2}{j}(qi)^jlfloorfrac{ri+b}{c} floor^{k_2-j}\ =sum_{j=0}^{k_2}inom{k_2}{j}q^jsum_{i=0}^{N}i^{k_1+j}lfloorfrac{ri+b}{c} floor^{k_2-j} ]

递归计算 (func(N,r,b,c)) 即可。

(bgeq c)

若 (bgeq c) ，则令 (q=lfloorfrac{b}{c} floor,r=bmod q)，所求式子即为

[sum_{i=0}^{N}i^{k_1}(q+lfloorfrac{ai+r}{c} floor)^{k_2}\ =sum_{i=0}^{N}i^{k_1}sum_{j=0}^{k_2}inom{k_2}{j}q^jlfloorfrac{ai+r}{c} floor^{k_2-j}\ =sum_{j=0}^{k_2}inom{k_2}{j}q^jsum_{i=0}^{N}i^{k_1}lfloorfrac{ai+r}{c} floor^{k_2-j} ]

递归计算 (func(N,a,r,c)) 即可。

一般情况

否则，即 (a < c,b < c,lfloorfrac{aN+b}{c} floor > 0)，令 (M=lfloorfrac{aN+b}{c} floor) 。

注意到 (lfloorfrac{ai+b}{c} floor^{k_2})​ 较难处理，需要将其变形，有

[lfloorfrac{ai+b}{c} floor^{k_2}=sum_{j=0}^{lfloorfrac{ai+b}{c} floor-1}((j+1)^{k_2}-j^{k_2}) ]

注意这里规定 (0^0=0)。

上述变换的好处是提供了交换求和顺序的可能，由此，原式可变为

[sum_{i=0}^{N}i^{k_1}sum_{j=0}^{lfloorfrac{ai+b}{c} floor-1}((j+1)^{k_2}-j^{k_2})\ =sum_{j=0}^{M-1}((j+1)^{k_2}-j^{k_2})sum_{i=0}^{N}i^{k_1}[lfloorfrac{ai+b}{c} floorgeq j+1]\ =sum_{j=0}^{M-1}((j+1)^{k_2}-j^{k_2})sum_{i=0}^{N}i^{k_1}[i>lfloorfrac{cj+c-b-1}{a} floor]\ =sum_{j=0}^{M-1}((j+1)^{k_2}-j^{k_2})sum_{i=0}^{N}i^{k_1}-sum_{j=0}^{M-1}((j+1)^{k_2}-j^{k_2})sum_{i=0}^{lfloorfrac{cj+c-b-1}{a} floor}i^{k_1} ]

靠前的求和符号可以直接计算，考虑靠后的求和符号。

((j+1)^{k_2}-j^{k_2})​ 显然是关于 (j) 的 (k_2-1) 次多项式，而 (sum_{i=0}^{lfloorfrac{cj+c-b-1}{a} floor}i^{k_1})​ 也是关于 (lfloorfrac{cj+c-b-1}{a} floor) 的 (k_1+1) 次多项式，不妨令为 (A(x),B(x))。

显然 (A) 是组合数，(B) 可以在一开始就用拉格朗日插值法预处理出来。

那么

[sum_{j=0}^{M-1}((j+1)^{k_2}-j^{k_2})sum_{i=0}^{lfloorfrac{cj+c-b-1}{a} floor}i^{k_1}\ =sum_{i=0}^{M-1}sum_{j=0}^{k_2-1}A_ji^jsum_{k=0}^{k_1+1}B_klfloorfrac{ci+c-b-1}{a} floor^{k}\ =sum_{j=0}^{k_2-1}A_jsum_{k=0}^{k_1+1}B_ksum_{i=0}^{M-1}i^jlfloorfrac{ci+c-b-1}{a} floor^{k} ]

递归计算 (func(M-1,c,c-b-1,a)) 即可。

注意到在递归中的参数变化：((a,c) ightarrow (amod c,c) ightarrow (c,amod c))，所以递归层数为 (O(log V))。

时间复杂度 (O(TK^4log V))。

poly operator*(CO poly&a,CO poly&b){
	int n=a.size()-1,m=b.size()-1;
	poly ans(n+m+1);
	for(int i=0;i<=n;++i)for(int j=0;j<=m;++j)
		ans[i+j]=add(ans[i+j],mul(a[i],b[j]));
	return ans;
}
poly operator/(poly a,CO poly&b){
	int n=a.size()-1,m=b.size()-1;
	poly ans(n-m+1);
	int inv=fpow(b[m],mod-2);
	for(int i=n;i>=m;--i){
		ans[i-m]=mul(a[i],inv);
		for(int j=i;j>=i-m;--j) a[j]=add(a[j],mod-mul(ans[i-m],b[j-(i-m)]));
	}
	return ans;
}
poly lagrange(int n,CO poly&x,CO poly&y){
	poly ans(n+1),p={1};
	for(int i=0;i<=n;++i) p=p*(poly){mod-x[i],1};
	for(int i=0;i<=n;++i){
		poly q=p/(poly){mod-x[i],1};
		int c=1;
		for(int j=0;j<=n;++j)if(j!=i) c=mul(c,add(x[i],mod-x[j]));
		c=fpow(c,mod-2);
		for(int j=0;j<=n;++j) ans[j]=add(ans[j],mul(c,mul(y[i],q[j])));
	}
	return ans;
}

CO int N=11;
poly sum[N];
int binom[N][N];
struct info {int v[N][N];};

info solve(int64 n,int64 a,int64 b,int64 c){
	info ans={};
	if(a==0 or (a*n+b)/c==0){
		for(int k1=0;k1<=10;++k1){
			int s=0;
			for(int i=0,t=1;i<=k1+1;++i,t=mul(t,n%mod))
				s=add(s,mul(sum[k1][i],t));
			int q=(a*n+b)/c%mod;
			for(int k2=0,t=1;k1+k2<=10;++k2,t=mul(t,q))
				ans.v[k1][k2]=mul(s,t);
		}
		return ans;
	}
	if(a>=c){
		info tmp=solve(n,a%c,b,c);
		for(int k1=0;k1<=10;++k1)for(int k2=0;k1+k2<=10;++k2){
			int q=a/c%mod;
			for(int i=0,t=1;i<=k2;++i,t=mul(t,q))
				ans.v[k1][k2]=add(ans.v[k1][k2],
				mul(binom[k2][i],mul(t,tmp.v[k1+i][k2-i])));
		}
		return ans;
	}
	if(b>=c){
		info tmp=solve(n,a,b%c,c);
		for(int k1=0;k1<=10;++k1)for(int k2=0;k1+k2<=10;++k2){
			int q=b/c%mod;
			for(int i=0,t=1;i<=k2;++i,t=mul(t,q))
				ans.v[k1][k2]=add(ans.v[k1][k2],
				mul(binom[k2][i],mul(t,tmp.v[k1][k2-i])));
		}
		return ans;
	}
	int64 m=(a*n+b)/c;
	info tmp=solve(m-1,c,c-b-1,a);
	for(int k1=0;k1<=10;++k1){
		int s=0;
		for(int i=0,t=1;i<=k1+1;++i,t=mul(t,n%mod))
			s=add(s,mul(sum[k1][i],t));
		for(int k2=0;k1+k2<=10;++k2){
			ans.v[k1][k2]=mul(fpow(m%mod,k2),s);
			for(int i=0;i<=k2-1;++i)for(int j=0;j<=k1+1;++j)
				ans.v[k1][k2]=add(ans.v[k1][k2],mod-
				mul(binom[k2][i],mul(sum[k1][j],tmp.v[i][j])));
		}
	}
	return ans;
}

int main(){
	for(int i=0;i<=10;++i){
		poly x(i+2),y(i+2);
		x[0]=0,y[0]=fpow(0,i);
		for(int j=1;j<=i+1;++j) x[j]=j,y[j]=add(y[j-1],fpow(j,i));
		sum[i]=lagrange(i+1,x,y);
	}
	for(int i=0;i<=10;++i){
		binom[i][0]=binom[i][i]=1;
		for(int j=1;j<i;++j) binom[i][j]=add(binom[i-1][j-1],binom[i-1][j]);
	}
	for(int T=read<int>();T--;){
		int64 n=read<int64>(),a=read<int64>(),b=read<int64>(),c=read<int64>();
		int k1=read<int>(),k2=read<int>();
		printf("%d
",solve(n,a,b,c).v[k1][k2]);
	}
	return 0;
}