类欧几里得算法
作用
比较快速地算出下面的式子
[F(n,a,b,c,k_1,k_2)=sumlimits_{x=0}^n x^{k_1} lfloorfrac{ax+b}{c}
floor ^{k_2}
]
步骤
不妨假设现在 (ageq c) 或 (b geq c) ，那么
[sumlimits_{x=0}^n x^{k_1} lfloorfrac{ax+b}{c}
floor ^{k_2} \
=sumlimits_{x=0}^n x^{k_1} (lfloorfrac ac
floor x+lfloorfrac bc
floor+lfloorfrac{(a mod c)x+(bmod c)}{c}
floor) ^{k_2}
]
设
[p=lfloorfrac ac
floor\
q=lfloorfrac bc
floor\
A=amod c\
B=bmod c\
]
那么
[sumlimits_{x=0}^n x^{k_1} (lfloorfrac ac
floor x+lfloorfrac bc
floor+lfloorfrac{(a mod c)x+(bmod c)}{c}
floor) ^{k_2}\
=sumlimits_{x=0}^n x^{k_1} (p x+q+lfloorfrac{Ax+B}{c}
floor) ^{k_2}\
=sumlimits_{i+j+k=k_2}frac{k_2!}{i!j!k!} p^iq^jsumlimits_{x=0}^nx^{k_1+i}lfloorfrac{Ax+B}{c}
floor^{k}\
=sumlimits_{i+j+k=k_2}frac{k_2!}{i!j!k!} p^iq^jF(n,A,B,c,k_1+i,k)
]
注意到这一步我们令 (a) 变成了 (amod c)
于是，我们现在只需讨论 (a,b<c) 的情况了，注意到此时 (lfloor frac{a(x+1)+b}{c} floor-lfloor frac{ax+b}{c} floor leq 1)
我们设 (m=lfloor frac{an+b}{c} floor) ,
[F(n,a,b,c,k_1,k_2)\
=sumlimits_{x=0}^n x^{k_1} lfloorfrac{ax+b}{c}
floor ^{k_2}\
=sumlimits_{x=0}^n x^{k_1} (sumlimits_{w=0}^{m-1} left[ lfloorfrac{ax+b}{c}
floor >w
ight] )^{k_2}\
=sumlimits_{x=0}^n x^{k_1} sumlimits_{w=0}^{m-1} left[ lfloorfrac{ax+b}{c}
floor >w
ight] ((w+1)^{k_2}-w^{k_2})\
lfloor frac{ax+b}{c}
floor > w \
Updownarrow\
lfloor frac{ax+b}{c}
floor geq w+1\
Updownarrow\
ax+b geq cw+c\
Updownarrow\
axgeq cw+c-b\
Updownarrow\
ax> cw+c-b-1\
Updownarrow\
x>lfloorfrac{ cw+c-b-1}{a}
floor\
F(n,a,b,c,k_1,k_2)\
=sumlimits_{w=0}^{m-1}((w+1)^{k_2}-w^{k_2})sumlimits_{x=0}^n x^{k_1} left[ x>lfloorfrac{ cw+c-b-1}{a}
floor
ight] \
=(sumlimits_{w=0}^{m-1}((w+1)^{k_2}-w^{k_2})sumlimits_{x=0}^n x^{k_1})- (sumlimits_{w=0}^{m-1}((w+1)^{k_2}-w^{k_2})sumlimits_{x=0}^{lfloorfrac{ cw+c-b-1}{a}
floor} x^{k_1} )\
=m^{k_2}sumlimits_{x=0}^n x^{k_1}- sumlimits_{w=0}^{m-1}((w+1)^{k_2}-w^{k_2})sumlimits_{x=0}^{lfloorfrac{ cw+c-b-1}{a}
floor} x^{k_1}\
]
注意此时左半部分可以用拉格朗日差值快速求得。
(sumlimits_{w=0}^{m-1}((w+1)^{k_2}-w^{k_2})) 是关于 (w) 的 (k_2-1) 次多项式，设为 (A) ， (i) 次项系数为 (A_i)
(sumlimits_{x=0}^{lfloorfrac{ cw+c-b-1}{a} floor} x^{k_1}) 是关于 (lfloorfrac{ cw+c-b-1}{a} floor) 的 (k_1+1) 次多项式，设为 (B) ， (i) 次项系数为 (B_i)
[sumlimits_{w=0}^{m-1}((w+1)^{k_2}-w^{k_2})sumlimits_{x=0}^{lfloorfrac{ cw+c-b-1}{a}
floor} x^{k_1}\
=sumlimits_{w=0}^{m-1}sumlimits_{i=0}^{k_2-1}A_iw^isumlimits_{j=0}^{k_1+1}B_jlfloorfrac{ cw+c-b-1}{a}
floor^{j}\
=sumlimits_{i=0}^{k_2-1}sumlimits_{j=0}^{k_1+1}A_iB_jsumlimits_{w=0}^{m-1}w^ilfloorfrac{ cw+c-b-1}{a}
floor^{j}\
=sumlimits_{i=0}^{k_2-1}sumlimits_{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
",F(k1,k2,a,b,c,n));
} return 0;
}