类欧几里得算法

Posted 2020-11-30 oyjason

类欧几里得算法

作用

比较快速地算出下面的式子
[ 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 floorq=lfloorfrac bc floorA=amod cB=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\Updownarrowax+b geq cw+c\Updownarrowaxgeq cw+c-b\Updownarrowax> cw+c-b-1\Updownarrowx>lfloorfrac{ cw+c-b-1}{a} floorF(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;
}