洛谷 P4245 [模板]任意模数NTT —— 三模数NTT / 拆系数FFT(MTT)

题目：https://www.luogu.org/problemnew/show/P4245

用三模数NTT做，需要注意时间和细节；

注意各种地方要取模！传入 upt() 里面的数一定要不超过2倍 mod！

乘法会爆 long long 时用快速乘！

两次合并的模数，第一次是 (ll) p1*p2，第二次直接对题目的模数取模即可！

注意局部开 (ll)！

合并时用到的逆元每次都一样，所以要先处理好而不是现场快速幂算！！

然而为什么时间还是 Narh 的两倍！

一晚上的心血...

代码如下：

#include<iostream>

#include<cstdio>

#include<cstring>

#include<algorithm>

using namespace std;

typedef long long ll;

int const xn=(<<);

int n,m,lim,rev[xn],a[][xn],b[][xn],p[]={,,,};

ll mod;

int rd()

{

  int ret=,f=; char ch=getchar();

  while(ch<''||ch>''){if(ch=='-')f=; ch=getchar();}

  while(ch>=''&&ch<='')ret=ret*+ch-'',ch=getchar();

  return f?ret:-ret;

}

ll upt(ll x,int md){while(x>=md)x-=md; while(x<)x+=md; return x;}

ll mul(ll a,ll b,int md)

{

  ll ret=; a=a%md; b=b%md;

  if(a<)a+=md; if(b<)b+=md;//

  for(;b;b>>=1ll,a=(a+a)%md)if(b&)ret=(ret+a)%md;

  return ret;

}

ll pw(ll a,int b,int md)

{

  ll ret=; a=a%md;

  for(;b;b>>=,a=mul(a,a,md))if(b&)ret=mul(ret,a,md);//mul!!

  return ret;

}

void ntt(int *a,int tp,int md)

{

  for(int i=;i<lim;i++)

    if(i<rev[i])swap(a[i],a[rev[i]]);

  for(int mid=;mid<lim;mid<<=)

    {

      int wn=pw(,(md-)/(mid<<),md);

      if(tp==-)wn=pw(wn,md-,md);

      for(int j=,len=(mid<<);j<lim;j+=len)

    {

      int w=;

      for(int k=;k<mid;k++,w=(ll)w*wn%md)

        {

          int x=a[j+k],y=(ll)w*a[j+mid+k]%md;

          a[j+k]=upt(x+y,md); a[j+mid+k]=upt(x-y,md);

        }

    }

    }

  if(tp==)return; int inv=pw(lim,md-,md);

  for(int i=;i<lim;i++)a[i]=(ll)a[i]*inv%md;

}

ll uni(ll r1,ll r2,ll m1,int m2,int tp,int inv)

{

  ll k=mul(r2-r1,inv,m2);

  if(!tp)return (r1+k*m1)%(m1*m2);

  return upt((r1+mul(k,m1,mod))%mod,mod);//%mod!!

}

int main()

{

  n=rd(); m=rd(); mod=rd();

  for(int i=;i<=n;i++)a[][i]=a[][i]=a[][i]=rd();

  for(int i=;i<=m;i++)b[][i]=b[][i]=b[][i]=rd();

  lim=; int l=;

  while(lim<=n+m+)lim<<=,l++;//+2!

  for(int i=;i<lim;i++)

    rev[i]=((rev[i>>]>>)|((i&)<<(l-)));

  for(int i=;i<=;i++)

    {

      ntt(a[i],,p[i]); ntt(b[i],,p[i]);

      for(int j=;j<lim;j++)a[i][j]=(ll)a[i][j]*b[i][j]%p[i];

      ntt(a[i],-,p[i]);

    }

  int inv1=pw(p[],p[]-,p[]);

  int inv2=pw((ll)p[]*p[],p[]-,p[]);//inv!!

  for(int i=;i<=n+m;i++)

    {

      ll ans=uni(a[][i],a[][i],p[],p[],,inv1);

      ans=uni(ans,a[][i],(ll)p[]*p[],p[],,inv2);//(ll)!!!

      printf("%lld ",ans);

    }

  puts("");

  return ;

}

关于拆系数FFT，这篇博客说得十分清晰：https://blog.csdn.net/lvzelong2014/article/details/80156989

而且代码也十分简洁优美，所以就模仿着写了；

注意：

1. 读入的初始数组要先取模；

2. 对 (x<<30) 开 (ll) 要写在括号里面而非外面；

3. IDFT中最后要 /lim，平常都之写 a[i].x/lim，但这里因为用到了 y，所以必须加上 a[i].y/lim！！

4. 卡精度，要开 long double

(注意数组范围，如果要用到 a[lim] 的话，再稍微把数组开大一点)

代码如下：

#include<iostream>

#include<cstdio>

#include<cstring>

#include<algorithm>

#include<cmath>

using namespace std;

typedef long long ll;

typedef long double db;//

int const xn=(<<)+;//a[lim]

db const Pi=acos(-1.0);

int n,m,lim,rev[xn],P,f[xn],g[xn],ans[xn];

struct com{db x,y;}a[xn],b[xn],Da[xn],Db[xn],Dc[xn],Dd[xn];

com operator + (com a,com b){return (com){a.x+b.x,a.y+b.y};}

com operator - (com a,com b){return (com){a.x-b.x,a.y-b.y};}

com operator * (com a,com b){return (com){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}

com conj(com a){return (com){a.x,-a.y};}

int rd()

{

  int ret=,f=; char ch=getchar();

  while(ch<''||ch>''){if(ch=='-')f=; ch=getchar();}

  while(ch>=''&&ch<='')ret=ret*+ch-'',ch=getchar();

  return f?ret:-ret;

}

void init()

{

  lim=; int l=;

  while(lim<=n+m)lim<<=,l++;

  for(int i=;i<lim;i++)rev[i]=((rev[i>>]>>)|((i&)<<(l-)));

}

int upt(int x){while(x>=P)x-=P; while(x<)x+=P; return x;}

void fft(com *a,int tp)

{

  for(int i=;i<lim;i++)

    if(i<rev[i])swap(a[i],a[rev[i]]);

  for(int mid=;mid<lim;mid<<=)

    {

      com wn=(com){cos(Pi/mid),tp*sin(Pi/mid)};

      for(int j=,len=(mid<<);j<lim;j+=len)

    {

      com w=(com){,};

      for(int k=;k<mid;k++,w=w*wn)

        {

          com x=a[j+k],y=w*a[j+mid+k];

          a[j+k]=x+y; a[j+mid+k]=x-y;

        }

    }

    }

  if(tp==)return;

  for(int i=;i<lim;i++)a[i].x/=lim,a[i].y/=lim;//y!! for use y

}

void mul(int *A,int *B,int *C)

{

  for(int i=;i<lim;i++)A[i]=upt(A[i]%P),B[i]=upt(B[i]%P);//

  int M=(<<)-;

  for(int i=;i<lim;i++)a[i]=(com){A[i]&M,A[i]>>};

  for(int i=;i<lim;i++)b[i]=(com){B[i]&M,B[i]>>};

  fft(a,); fft(b,);

  a[lim]=a[]; b[lim]=b[];//

  for(int i=,j=lim;i<lim;i++,j--)

    {

      com da,db,dc,dd;

      da=(a[i]+conj(a[j]))*(com){0.5,};

      db=(a[i]-conj(a[j]))*(com){,-0.5};

      dc=(b[i]+conj(b[j]))*(com){0.5,};

      dd=(b[i]-conj(b[j]))*(com){,-0.5};

      Da[i]=da*dc; Db[i]=da*dd; Dc[i]=db*dc; Dd[i]=db*dd;

    }

  a[lim]=b[lim]=(com){,};

  for(int i=;i<lim;i++)a[i]=Da[i]+Db[i]*(com){,};

  for(int i=;i<lim;i++)b[i]=Dc[i]+Dd[i]*(com){,};

  fft(a,-); fft(b,-);

  for(int i=;i<=n+m;i++)

    {

      int da=(ll)(a[i].x+0.5)%P;

      int db=(ll)(a[i].y+0.5)%P;

      int dc=(ll)(b[i].x+0.5)%P;

      int dd=(ll)(b[i].y+0.5)%P;

      C[i]=(da+((ll)(db+dc)<<)+((ll)dd<<))%P;//(ll)

    }

}

int main()

{

  n=rd(); m=rd(); P=rd(); init();

  for(int i=;i<=n;i++)f[i]=rd();

  for(int i=;i<=m;i++)g[i]=rd();

  mul(f,g,ans);

  for(int i=;i<=n+m;i++)printf("%d ",upt(ans[i])); puts("");

  return ;

}