任意模数 NTT

2020-09-12 13:00:21 阅读：234 来源： 互联网

标签：MOD2 MOD3 int ll NTT 模数任意 inv mod

#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

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

ICode9

任意模数 NTT