LOJ138 类欧几里得算法
类欧几里得算法
给出 \(T\) 组询问，每组用 \(n, a, b, c, k_1, k_2\) 来描述。对于每组询问，请你求出
对 \(1000000007\) 取模。
对于 \(100 \%\) 的数据，\(T = 1000, 1 \le n, a, b, c \le {10} ^ 9, 0 \le k_1 + k_2 \le 10\) 。
数形结合
这里只计算
首先通过简单转化把 \(a,b\) 都放到 \(\bmod c\) 的意义下。
考虑计算式 \(\sum_{i=0}^{N}\lfloor\frac{ai+b}{c}\rfloor\) 的含义。它要求的其实是 \(y=\frac{a}{c}x+\frac{b}{c}~(0 \leq x \leq N)\) 这条线段下方（包括线段上）的整点个数。要求整点 \((x\geq 0,y > 0)\)。
使用差分，转化为用 \((0,0) \rightarrow (N,\frac{aN+b}{c})\) 这个矩形里的整点个数减去 \(y=\frac{a}{c}x+\frac{b}{c}~(0 \leq x \leq N)\) 这条线段上方（不包括线段上）的整点个数。要求整点 \((x\geq 0,y > 0)\)。
然后翻转坐标轴， 转化成求 \(y=\frac{c}{a}x-\frac{b}{a}~(0 < x \leq \frac{aN+b}{c})\) 这条线段下方（不包括线段上）的整点个数。要求整点 \((x>0,y\geq 0)\)。
化归原问题，求 \(y=\frac{c}{a}(x+1)-\frac{b}{a}-\frac{1}{a}~(0 \leq x < \frac{aN+b}{c})\) 这条线段下方（包括线段上）的整点个数。要求整点 \((x\geq 0,y \geq 0)\)。\(y=0\) 的那些点要特殊处理掉，由于 \(c-b-1\geq 0\)，所以每个可能的横坐标都对应了一个这样的点。
此时相当于求
此时 \(c>a\)，所以可以再次取模，进行递归。注意到在递归中的参数变化：\((a,c)\rightarrow (a\bmod c,c)\rightarrow (c,a\bmod 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)\) ，计算 \(0\leq k_1+k_2\leq10\) 的情况下所求式子的值，即
递归基
若 \(a=0\) 或 \(\lfloor\frac{aN+b}{c}\rfloor=0\) ，那么， \(\lfloor\frac{ai+b}{c}\rfloor\) 的取值始终相同，答案即为
可以用插值解决。
\(a \geq c\)
若 \(a\geq c\) ，则令 \(q=\lfloor\frac{a}{c}\rfloor,r=a\bmod c\) ，所求式子即为
递归计算 \(func(N,r,b,c)\) 即可。
\(b\geq c\)
若 \(b\geq c\) ，则令 \(q=\lfloor\frac{b}{c}\rfloor,r=b\bmod q\)，所求式子即为
递归计算 \(func(N,a,r,c)\) 即可。
一般情况
否则，即 \(a < c,b < c,\lfloor\frac{aN+b}{c}\rfloor > 0\)，令 \(M=\lfloor\frac{aN+b}{c}\rfloor\) 。
注意到 \(\lfloor\frac{ai+b}{c}\rfloor^{k_2}\) 较难处理，需要将其变形，有
注意这里规定 \(0^0=0\)。
上述变换的好处是提供了交换求和顺序的可能，由此，原式可变为
靠前的求和符号可以直接计算，考虑靠后的求和符号。
\((j+1)^{k_2}-j^{k_2}\) 显然是关于 \(j\) 的 \(k_2-1\) 次多项式，而 \(\sum_{i=0}^{\lfloor\frac{cj+c-b-1}{a}\rfloor}i^{k_1}\) 也是关于 \(\lfloor\frac{cj+c-b-1}{a}\rfloor\) 的 \(k_1+1\) 次多项式，不妨令为 \(A(x),B(x)\)。
显然 \(A\) 是组合数，\(B\) 可以在一开始就用拉格朗日插值法预处理出来。
那么
递归计算 \(func(M-1,c,c-b-1,a)\) 即可。
注意到在递归中的参数变化：\((a,c)\rightarrow (a\bmod c,c)\rightarrow (c,a\bmod c)\)，所以递归层数为 \(O(\log V)\)。
时间复杂度 \(O(TK^4\log 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\n",solve(n,a,b,c).v[k1][k2]);
}
return 0;
}
