类欧几里得算法
类欧几里得算法
作用
比较快速地算出下面的式子
\[F(n,a,b,c,k_1,k_2)=\sum\limits_{x=0}^n x^{k_1} \lfloor\frac{ax+b}{c}\rfloor ^{k_2}
\]
步骤
不妨假设现在 \(a\geq c\) 或 \(b \geq c\) ，那么
\[\sum\limits_{x=0}^n x^{k_1} \lfloor\frac{ax+b}{c}\rfloor ^{k_2} \\
=\sum\limits_{x=0}^n x^{k_1} (\lfloor\frac ac\rfloor x+\lfloor\frac bc\rfloor+\lfloor\frac{(a \mod c)x+(b\mod c)}{c}\rfloor) ^{k_2}
\]
设
\[p=\lfloor\frac ac\rfloor\\
q=\lfloor\frac bc\rfloor\\
A=a\mod c\\
B=b\mod c\\
\]
那么
\[\sum\limits_{x=0}^n x^{k_1} (\lfloor\frac ac\rfloor x+\lfloor\frac bc\rfloor+\lfloor\frac{(a \mod c)x+(b\mod c)}{c}\rfloor) ^{k_2}\\
=\sum\limits_{x=0}^n x^{k_1} (p x+q+\lfloor\frac{Ax+B}{c}\rfloor) ^{k_2}\\
=\sum\limits_{i+j+k=k_2}\frac{k_2!}{i!j!k!} p^iq^j\sum\limits_{x=0}^nx^{k_1+i}\lfloor\frac{Ax+B}{c}\rfloor^{k}\\
=\sum\limits_{i+j+k=k_2}\frac{k_2!}{i!j!k!} p^iq^jF(n,A,B,c,k_1+i,k)
\]
注意到这一步我们令 \(a\) 变成了 \(a\mod c\)
于是，我们现在只需讨论 \(a,b<c\) 的情况了，注意到此时 \(\lfloor \frac{a(x+1)+b}{c}\rfloor-\lfloor \frac{ax+b}{c}\rfloor \leq 1\)
我们设 \(m=\lfloor \frac{an+b}{c}\rfloor\) ,
\[F(n,a,b,c,k_1,k_2)\\
=\sum\limits_{x=0}^n x^{k_1} \lfloor\frac{ax+b}{c}\rfloor ^{k_2}\\
=\sum\limits_{x=0}^n x^{k_1} (\sum\limits_{w=0}^{m-1} \left[ \lfloor\frac{ax+b}{c}\rfloor >w \right] )^{k_2}\\
=\sum\limits_{x=0}^n x^{k_1} \sum\limits_{w=0}^{m-1} \left[ \lfloor\frac{ax+b}{c}\rfloor >w \right] ((w+1)^{k_2}-w^{k_2})\\
\lfloor \frac{ax+b}{c}\rfloor > w \\
\Updownarrow\\
\lfloor \frac{ax+b}{c}\rfloor \geq w+1\\
\Updownarrow\\
ax+b \geq cw+c\\
\Updownarrow\\
ax\geq cw+c-b\\
\Updownarrow\\
ax> cw+c-b-1\\
\Updownarrow\\
x>\lfloor\frac{ cw+c-b-1}{a}\rfloor\\
F(n,a,b,c,k_1,k_2)\\
=\sum\limits_{w=0}^{m-1}((w+1)^{k_2}-w^{k_2})\sum\limits_{x=0}^n x^{k_1} \left[ x>\lfloor\frac{ cw+c-b-1}{a}\rfloor \right] \\
=(\sum\limits_{w=0}^{m-1}((w+1)^{k_2}-w^{k_2})\sum\limits_{x=0}^n x^{k_1})- (\sum\limits_{w=0}^{m-1}((w+1)^{k_2}-w^{k_2})\sum\limits_{x=0}^{\lfloor\frac{ cw+c-b-1}{a}\rfloor} x^{k_1} )\\
=m^{k_2}\sum\limits_{x=0}^n x^{k_1}- \sum\limits_{w=0}^{m-1}((w+1)^{k_2}-w^{k_2})\sum\limits_{x=0}^{\lfloor\frac{ cw+c-b-1}{a}\rfloor} x^{k_1}\\
\]
注意此时左半部分可以用拉格朗日差值快速求得。
\(\sum\limits_{w=0}^{m-1}((w+1)^{k_2}-w^{k_2})\) 是关于 \(w\) 的 \(k_2-1\) 次多项式，设为 \(A\) ， \(i\) 次项系数为 \(A_i\)
\(\sum\limits_{x=0}^{\lfloor\frac{ cw+c-b-1}{a}\rfloor} x^{k_1}\) 是关于 \(\lfloor\frac{ cw+c-b-1}{a}\rfloor\) 的 \(k_1+1\) 次多项式，设为 \(B\) ， \(i\) 次项系数为 \(B_i\)
\[\sum\limits_{w=0}^{m-1}((w+1)^{k_2}-w^{k_2})\sum\limits_{x=0}^{\lfloor\frac{ cw+c-b-1}{a}\rfloor} x^{k_1}\\
=\sum\limits_{w=0}^{m-1}\sum\limits_{i=0}^{k_2-1}A_iw^i\sum\limits_{j=0}^{k_1+1}B_j\lfloor\frac{ cw+c-b-1}{a} \rfloor^{j}\\
=\sum\limits_{i=0}^{k_2-1}\sum\limits_{j=0}^{k_1+1}A_iB_j\sum\limits_{w=0}^{m-1}w^i\lfloor\frac{ cw+c-b-1}{a}\rfloor^{j}\\
=\sum\limits_{i=0}^{k_2-1}\sum\limits_{j=0}^{k_1+1}A_iB_jF(m-1,c,c-b-1,a,i,j)
\]
此时我们本质上在接下来的运算中交换了 \(a,c\)
因此总的迭代次数是 \(O(log)\) 的，并且可以根据迭代的层数作为下标记忆化，并且 \(k_1,k_2\) 的和不会增大。
如果 \(a=0\) 或 \(m=0\) 或 \(k_2=0\) 之类的边界情况可以直接拉格朗日插值快速求得。
LOJ138 AC代码如下
#include<bits/stdc++.h>
using namespace std;
#define LL long long
#define debug(x) cerr<<#x<<" = "<<x
#define sp <<" "
#define el <<endl
#define fgx cerr<<" ---------------------------------------------- "<<endl
#define uint unsigned int
#define ULL unsigned long long
#define DB long double
#define pii pair<int,int>
#define mp make_pair
#define pb push_back
inline int read(){
int nm=0; bool fh=true; char cw=getchar();
for(;!isdigit(cw);cw=getchar()) fh^=(cw=='-');
for(;isdigit(cw);cw=getchar()) nm=nm*10+(cw-'0');
return fh?nm:-nm;
}
#define mod 1000000007
namespace CALC{
inline int add(int x,int y){return (x+y>=mod)?(x+y-mod):(x+y);}
inline int mns(int x,int y){return (x-y<0)?(x-y+mod):(x-y);}
inline int mul(int x,int y){return (LL)x*(LL)y%mod;}
inline void upd(int &x,int y){x=((x+y>=mod)?(x+y-mod):(x+y));}
inline void dec(int &x,int y){x=((x-y<0)?(x-y+mod):(x-y));}
inline int qpow(int x,int sq){int res=1;for(;sq;sq>>=1,x=mul(x,x))if(sq&1)res=mul(res,x);return res;}
}using namespace CALC;
int fac[50],ifac[50],C[50][50],G[80][11][11];
namespace Lagrange{
const int n=11;
#define M 14
int B[M][M],A[M][M];
inline int _C(int N,int K){if(N<K||K<0)return 0;return mul(fac[N],mul(ifac[N-K],ifac[K]));}
inline void solve(int *f,int m){
static int t[M][M];
for(int i=0;i<=n;i++){
C[i][0]=1;
for(int j=1;j<=i;j++) C[i][j]=add(C[i-1][j],C[i-1][j-1]);
}
for(int i=0;i<=m;i++){
t[i][0]=1,t[i][m+1]=f[i+1];
for(int j=1;j<=m;j++) t[i][j]=mul(t[i][j-1],i+1);
}
for(int i=0,k;i<=m;i++){
for(k=i;!t[k][i];++k);
if(k^i){for(int j=k;j<=m+1;j++)swap(t[k][j],t[i][j]);}
int bs=qpow(t[i][i],mod-2);
for(int j=i;j<=m+1;j++) t[i][j]=mul(t[i][j],bs);
for(int w=0;w<=m;w++) if(w!=i&&t[w][i]>0)
for(int tmp=t[w][i],j=i;j<=m+1;++j) dec(t[w][j],mul(tmp,t[i][j]));
}
for(int i=0;i<=m;i++) f[i]=t[i][m+1];
}
inline int calc_B(int X,int k){
if(!k) return X+1; int res=0,bs=1;
for(int i=0;i<=k+1;i++,bs=mul(bs,X)) upd(res,mul(B[k][i],bs));
return res;
}
void init(){
fac[0]=ifac[0]=1;
for(int i=1;i<=20;i++) fac[i]=mul(fac[i-1],i);
for(int i=1;i<=20;i++) ifac[i]=qpow(fac[i],mod-2);
for(int i=1;i<=n;i++) for(int j=1;j<=i+2;j++) B[i][j]=add(B[i][j-1],qpow(j,i));
for(int i=1;i<=n;i++) for(int j=1;j<=i;j++) A[i][j]=mns(qpow(j+1,i),qpow(j,i));
for(int i=1;i<=n;i++) solve(A[i],i-1),solve(B[i],i+1); B[0][0]=B[0][1]=1;
}
}
using Lagrange::A;
using Lagrange::B;
using Lagrange::calc_B;
namespace __Euclid{
inline int F(int k1,int k2,int a,int b,int c,int n,int lev=0){
if(G[lev][k1][k2]!=-1) return G[lev][k1][k2]; int &res=G[lev][k1][k2]; res=0;
if(!k2) return res=Lagrange::calc_B(n,k1);
if((LL)a*(LL)n+b<(LL)c) res=0;
if(!k2||!a){
res=mul(calc_B(n,k1),qpow(b/c,k2));
return res;
}
if(a>=c||b>=c){
int ac[12],bc[12]; ac[0]=bc[0]=1;
for(int i=1;i<=k2;i++) ac[i]=mul(ac[i-1],a/c),bc[i]=mul(bc[i-1],b/c);
for(int i=1;i<=k2;i++) ac[i]=mul(ac[i],ifac[i]),bc[i]=mul(bc[i],ifac[i]);
for(int i=0;i<=k2;++i) for(int j=0;i+j<=k2;++j){
int k=k2-i-j,tmp=mul(mul(ac[i],bc[j]),mul(fac[k2],ifac[k]));
upd(res,mul(F(k1+i,k,a%c,b%c,c,n,lev+1),tmp));
}
return res;
}
int m=((LL)a*(LL)n+(LL)b)/c; res=mul(qpow(m,k2),calc_B(n,k1));
for(int p=0;p<k2;++p) for(int q=0;q<=k1+1;++q)
dec(res,mul(mul(C[k2][p],B[k1][q]),F(p,q,c,c-b-1,a,m-1,lev+1)));
return res;
}
} using __Euclid::F;
int main(){
Lagrange::init();
for(int Cas=read();Cas;--Cas){
int n=read(),a=read(),b=read(),c=read(),k1=read(),k2=read();
memset(G,-1,sizeof(G)),printf("%d\n",F(k1,k2,a,b,c,n));
} return 0;
}
