类欧几里得算法

标签：步骤   upd   tchar   cpp   log   tps   display   作用   define   

类欧几里得算法

作用

比较快速地算出下面的式子
\[ 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;
}

类欧几里得算法

标签：步骤   upd   tchar   cpp   log   tps   display   作用   define   

原文地址：https://www.cnblogs.com/OYJason/p/11806416.html