多项式相关算法及模板(未完待续)

  • FFT模板
const int N=3e5+10;
const double PI=acos(-1.0);
struct cd{
    double x,y;
    cd(double a=0,double b=0):x(a),y(b){}
};
inline cd operator +(cd a,cd b){return cd(a.x+b.x,a.y+b.y);}
inline cd operator -(cd a,cd b){return cd(a.x-b.x,a.y-b.y);}
inline cd operator *(cd a,cd b){return cd(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
inline cd conj(cd a){return cd(a.x,-a.y);}
struct FastFourierTransform{
    int n,rev[N];
    cd w[2][N];
    void init(int lim){
        n=1;int k=0;
        while (n<lim) n<<=1,k++;
        for (int k=0;k<n;k++){
            w[0][k]=cd(cos(2*PI/n*k),sin(2*PI/n*k));
            w[1][k]=conj(w[0][k]);
        }
    }
    void fft(cd *a,int v){
        for (int i=0,j=0;i<n;i++){
            if(i<j) swap(a[i],a[j]);
            for (int l=n>>1;(j^=l)<l;l>>=1);
        }
        for (int l=2;l<=n;l<<=1){
            int m=l>>1;
            for (int i=0;i<n;i+=l){ 
                for (int k=0;k<m;k++){
                    cd t=w[v][n/l*k]*a[i+k+m];
                    a[i+k+m]=a[i+k]-t;
                    a[i+k]=a[i+k]+t;
                }
            }
        }
        if (!v) return;
        for (int i=0;i<n;i++) a[i].x/=n;
    }
    void mul(cd *a,cd *b,int m){
        init(m);
        fft(a,0),fft(b,0);
        for(int i=0;i<n;i++) a[i]=a[i]*b[i];
        fft(a,1);
    }
}f;
  • NTT模板
const int P=998244353;//3是998244353的原根
const int N=800000+10;
int fexp(int a,int n){
    int res=1;
    for (;n;n>>=1,a=(ll)a*a%P) if(n&1)res=(ll)res*a%P;
    return res;
}
struct NumberTheoreticTransform{
    int n,w[2][N];
    void init(int lim){
        for (n=1;n<lim;n<<=1);
        w[0][0]=w[0][n]=1;
        int g=fexp(3,(P-1)/n);
        for (int i=1;i<n;i++) w[0][i]=(ll)w[0][i-1]*g%P;
        for (int i=0;i<=n;i++) w[1][i]=w[0][n-i];
    }
    void ntt(int* a,int v){
        for (int i=0,j=0;i<n;i++){
            if(i>j) swap(a[i],a[j]);
            for (int l=n>>1;(j^=l)<l;l>>=1);
        }
        for (int l=2;l<=n;l<<=1){
            int m=l>>1;
            for (int i=0;i<n;i+=l){
                for (int k=0;k<m;k++){
                    int t=(ll)w[v][n/l*k]*a[i+k+m]%P;
                    a[i+k+m]=(a[i+k]-t>=0)?(a[i+k]-t):(a[i+k]-t+P);
                    a[i+k]=(a[i+k]+t<P)?(a[i+k]+t):(a[i+k]+t-P);
                }
            }
        }
        if (!v) return;
        int rv=fexp(n,P-2);
        for (int i=0;i<n;i++) a[i]=(ll)a[i]*rv%P;
    }
    void work(int* a,int* b,int m){
        init(m);
        ntt(a,0),ntt(b,0);
        for (int i=0;i<n;i++) a[i]=(ll)a[i]*b[i]%P;
        ntt(a,1);
    }
}NTT;
  • 任意模数FFT
#define double long double
const int N=5e5+10;
const double PI=acos(-1.0);
struct cd{
    double x,y;
    cd(double a=0,double b=0):x(a),y(b){}
};
inline cd operator +(cd a,cd b){return cd(a.x+b.x,a.y+b.y);}
inline cd operator -(cd a,cd b){return cd(a.x-b.x,a.y-b.y);}
inline cd operator *(cd a,cd b){return cd(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
inline cd conj(cd a){return cd(a.x,-a.y);}
int rev[N];
cd w[2][N];
void init(int n){
    for (int k=0;k<n;k++){
        w[0][k]=cd(cos(2*PI/n*k),sin(2*PI/n*k));
        w[1][k]=conj(w[0][k]);
    } 
}
void fft(cd *a,int n,int v){
    for (int i=0,j=0;i<n;i++){
        if(i<j) swap(a[i],a[j]);
        for (int l=n>>1;(j^=l)<l;l>>=1);
    }
    for (int l=2;l<=n;l<<=1){
        int m=l>>1;
        for (int i=0;i<n;i+=l){ 
            for (int k=0;k<m;k++){
                cd t=w[v][n/l*k]*a[i+k+m];
                a[i+k+m]=a[i+k]-t;
                a[i+k]=a[i+k]+t;
            }
        }
    }
    if (!v) return;
    for (int i=0;i<n;i++) a[i].x/=n;
}
int n,m,mod,i,f[N],g[N],h[N];
void mtt(int deg,int*lhs,int*rhs,int*ret){
    static cd A[N],B[N],C[N],D[N],
              E[N],F[N],G[N],H[N];
    int len,i;
    for (len=1;len<=deg;len<<=1);
    init(len);
    for (i=0;i<len;i++){
        lhs[i]%=mod,rhs[i]%=mod;
        A[i].x=lhs[i]>>15;B[i].x=lhs[i]&0x7fff;
        C[i].x=rhs[i]>>15;D[i].x=rhs[i]&0x7fff;
    }
    fft(A,len,0),fft(B,len,0),fft(C,len,0),fft(D,len,0);
    for (i=0;i<len;i++){
        E[i]=A[i]*C[i];F[i]=B[i]*C[i];
        G[i]=A[i]*D[i];H[i]=B[i]*D[i];
    }
    fft(E,len,1),fft(F,len,1),fft(G,len,1),fft(H,len,1);
    for (i=0;i<len;i++){
        ret[i]=((((ll)(E[i].x+0.5)%mod)<<30)%mod+(((ll)(F[i].x+0.5)%mod)<<15)%mod
            +(((ll)(G[i].x+0.5)%mod)<<15)%mod+(ll)(H[i].x+0.5))%mod;
    }
}
int main(){
    read(n),read(m),read(mod);
    for (i=0;i<=n;i++) read(f[i]);
    for (i=0;i<=m;i++) read(g[i]);
    mtt(n+m,f,g,h);
    for (i=0;i<=n+m;i++) printf("%d%c",h[i],i==n+m?'\n':' ');
    return 0;
}

发表评论

电子邮件地址不会被公开。 必填项已用*标注