任意模数 NTT

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll MAXN=1<<21,MOD1=998244353,MOD2=1004535809,MOD3=469762049;
inline ll fpow(ll a,ll x,ll mod) { ll ans=1; for(;x;x>>=1,a=a*a%mod) if(x&1) ans=ans*a%mod; return ans; }
inline ll inv(ll a,ll mod) { return fpow(a,mod-2,mod); }
ll L,Rev[MAXN],invL1,invL2,invL3;
inline void preNTT(){
    Rev[1]=L>>1;
    for(int i=2;i<L;i++) Rev[i]=(Rev[i>>1]>>1)|Rev[i&1];
    invL1=inv(L,MOD1);
    invL2=inv(L,MOD2);
    invL3=inv(L,MOD3);
}
inline void NTT(int f[MAXN][4],ll op){
    ll g1=(op==1?3:inv(3,MOD1)),g2=(op==1?3:inv(3,MOD2)),g3=(op==1?3:inv(3,MOD3));
    for(int i=1;i<L;i++) if(i<Rev[i])
        for(int j=1;j<=3;j++) swap(f[i][j],f[Rev[i]][j]);
    for(int bas=2,mid=1,lb=1;bas<=L;bas<<=1,mid<<=1,lb++)
        for(int l=0,omega1=fpow(g1,MOD1>>lb,MOD1),omega2=fpow(g2,MOD2>>lb,MOD2),omega3=fpow(g3,MOD3>>lb,MOD3);l<L;l+=bas)
            for(int k=l,buf1=1,buf2=1,buf3=1,tmp1,tmp2,tmp3;
            tmp1=1ll*buf1*f[k+mid][1]%MOD1,tmp2=1ll*buf2*f[k+mid][2]%MOD2,tmp3=1ll*buf3*f[k+mid][3]%MOD3,k<mid+l;
            k++,buf1=1ll*buf1*omega1%MOD1,buf2=1ll*buf2*omega2%MOD2,buf3=1ll*buf3*omega3%MOD3)
                    f[k+mid][1]=(0ll+f[k][1]-tmp1)%MOD1,f[k][1]=(0ll+f[k][1]+tmp1)%MOD1,
                    f[k+mid][2]=(0ll+f[k][2]-tmp2)%MOD2,f[k][2]=(0ll+f[k][2]+tmp2)%MOD2,
                    f[k+mid][3]=(0ll+f[k][3]-tmp3)%MOD3,f[k][3]=(0ll+f[k][3]+tmp3)%MOD3;
    if(op==1){
        for(int i=0;i<L;i++)
            f[i][1]=(0ll+f[i][1]+MOD1)%MOD1,f[i][2]=(0ll+f[i][2]+MOD2)%MOD2,f[i][3]=(0ll+f[i][3]+MOD3)%MOD3;
    }
    else{
        for(int i=0;i<L;i++)
            f[i][1]=(0ll+f[i][1]+MOD1)*invL1%MOD1,f[i][2]=(0ll+f[i][2]+MOD2)*invL2%MOD2,f[i][3]=(0ll+f[i][3]+MOD3)*invL3%MOD3;
    }
}
inline void mul(int f[MAXN][4],int g[MAXN][4],int h[MAXN][4],int N){
    static int p[MAXN][4],q[MAXN][4];
    L=1;
    while(L<(N<<1)) L<<=1;
    preNTT();
    for(int i=0;i<L;i++)
        p[i][1]=f[i][1],q[i][1]=g[i][1],
        p[i][2]=f[i][2],q[i][2]=g[i][2],
        p[i][3]=f[i][3],q[i][3]=g[i][3];
    NTT(p,1);
    NTT(q,1);
    for(int i=0;i<L;i++)
        h[i][1]=1ll*p[i][1]*q[i][1]%MOD1,
        h[i][2]=1ll*p[i][2]*q[i][2]%MOD2,
        h[i][3]=1ll*p[i][3]*q[i][3]%MOD3;
    NTT(h,-1);
}
inline int crt(ll a1,ll m1,ll a2,ll m2,ll a3,ll m3,ll p){
    __int128 mod=1,sum=0; mod=mod*m1*m2*m3;
    sum+=mod/m1*inv(mod/m1%m1,m1)*a1;
    sum+=mod/m2*inv(mod/m2%m2,m2)*a2;
    sum+=mod/m3*inv(mod/m3%m3,m3)*a3;
    return sum%mod%p;
}

int F[MAXN][4],G[MAXN][4],H[MAXN][4],N,M,P;
inline void work(){
    cin>>N>>M>>P;
    for(int i=0;i<=N;i++) cin>>F[i][0],F[i][3]=F[i][2]=F[i][1]=F[i][0];
    for(int i=0;i<=M;i++) cin>>G[i][0],G[i][3]=G[i][2]=G[i][1]=G[i][0];
    mul(F,G,H,N+M+2);
    for(int i=0;i<L;i++) H[i][0]=crt(H[i][1],MOD1,H[i][2],MOD2,H[i][3],MOD3,P);
}
int main(){
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    work();
    for(int i=0;i<=N+M;i++) cout<<H[i][0]<<" ";
    cout<<endl;
    return 0;
}

标签：MOD2,MOD3,int,ll,NTT,模数,任意,inv,mod
来源： https://www.cnblogs.com/JustinRochester/p/13656635.html