洛谷 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=(1<<19);
int n,m,lim,rev[xn],a[5][xn],b[5][xn],p[5]={0,469762049,998244353,1004535809};
ll mod;
int rd()
{
  int ret=0,f=1; char ch=getchar();
  while(ch<'0'||ch>'9'){if(ch=='-')f=0; ch=getchar();}
  while(ch>='0'&&ch<='9')ret=ret*10+ch-'0',ch=getchar();
  return f?ret:-ret;
}
ll upt(ll x,int md){while(x>=md)x-=md; while(x<0)x+=md; return x;}
ll mul(ll a,ll b,int md)
{
  ll ret=0; a=a%md; b=b%md;
  if(a<0)a+=md; if(b<0)b+=md;//
  for(;b;b>>=1ll,a=(a+a)%md)if(b&1)ret=(ret+a)%md;
  return ret;
}
ll pw(ll a,int b,int md)
{
  ll ret=1; a=a%md;
  for(;b;b>>=1,a=mul(a,a,md))if(b&1)ret=mul(ret,a,md);//mul!!
  return ret;
}
void ntt(int *a,int tp,int md)
{
  for(int i=0;i<lim;i++)
    if(i<rev[i])swap(a[i],a[rev[i]]);
  for(int mid=1;mid<lim;mid<<=1)
    {
      int wn=pw(3,(md-1)/(mid<<1),md);
      if(tp==-1)wn=pw(wn,md-2,md);
      for(int j=0,len=(mid<<1);j<lim;j+=len)
    {
      int w=1;
      for(int k=0;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==1)return; int inv=pw(lim,md-2,md);
  for(int i=0;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=0;i<=n;i++)a[1][i]=a[2][i]=a[3][i]=rd();
  for(int i=0;i<=m;i++)b[1][i]=b[2][i]=b[3][i]=rd();
  lim=1; int l=0;
  while(lim<=n+m+2)lim<<=1,l++;//+2!
  for(int i=0;i<lim;i++)
    rev[i]=((rev[i>>1]>>1)|((i&1)<<(l-1)));
  for(int i=1;i<=3;i++)
    {
      ntt(a[i],1,p[i]); ntt(b[i],1,p[i]);
      for(int j=0;j<lim;j++)a[i][j]=(ll)a[i][j]*b[i][j]%p[i];
      ntt(a[i],-1,p[i]);
    }
  int inv1=pw(p[1],p[2]-2,p[2]);
  int inv2=pw((ll)p[1]*p[2],p[3]-2,p[3]);//inv!!
  for(int i=0;i<=n+m;i++)
    {
      ll ans=uni(a[1][i],a[2][i],p[1],p[2],0,inv1);
      ans=uni(ans,a[3][i],(ll)p[1]*p[2],p[3],1,inv2);//(ll)!!!
      printf("%lld ",ans);
    }
  puts("");
  return 0;
}

关于拆系数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=(1<<18)+5;//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=0,f=1; char ch=getchar();
  while(ch<'0'||ch>'9'){if(ch=='-')f=0; ch=getchar();}
  while(ch>='0'&&ch<='9')ret=ret*10+ch-'0',ch=getchar();
  return f?ret:-ret;
}
void init()
{
  lim=1; int l=0;
  while(lim<=n+m)lim<<=1,l++;
  for(int i=0;i<lim;i++)rev[i]=((rev[i>>1]>>1)|((i&1)<<(l-1)));
}
int upt(int x){while(x>=P)x-=P; while(x<0)x+=P; return x;}
void fft(com *a,int tp)
{
  for(int i=0;i<lim;i++)
    if(i<rev[i])swap(a[i],a[rev[i]]);
  for(int mid=1;mid<lim;mid<<=1)
    {
      com wn=(com){cos(Pi/mid),tp*sin(Pi/mid)};
      for(int j=0,len=(mid<<1);j<lim;j+=len)
    {
      com w=(com){1,0};
      for(int k=0;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==1)return;
  for(int i=0;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=0;i<lim;i++)A[i]=upt(A[i]%P),B[i]=upt(B[i]%P);//
  int M=(1<<15)-1;
  for(int i=0;i<lim;i++)a[i]=(com){A[i]&M,A[i]>>15};
  for(int i=0;i<lim;i++)b[i]=(com){B[i]&M,B[i]>>15};
  fft(a,1); fft(b,1);
  a[lim]=a[0]; b[lim]=b[0];//
  for(int i=0,j=lim;i<lim;i++,j--)
    {
      com da,db,dc,dd;
      da=(a[i]+conj(a[j]))*(com){0.5,0};
      db=(a[i]-conj(a[j]))*(com){0,-0.5};
      dc=(b[i]+conj(b[j]))*(com){0.5,0};
      dd=(b[i]-conj(b[j]))*(com){0,-0.5};
      Da[i]=da*dc; Db[i]=da*dd; Dc[i]=db*dc; Dd[i]=db*dd;
    }
  a[lim]=b[lim]=(com){0,0};
  for(int i=0;i<lim;i++)a[i]=Da[i]+Db[i]*(com){0,1};
  for(int i=0;i<lim;i++)b[i]=Dc[i]+Dd[i]*(com){0,1};
  fft(a,-1); fft(b,-1);
  for(int i=0;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)<<15)+((ll)dd<<30))%P;//(ll)
    }
}
int main()
{
  n=rd(); m=rd(); P=rd(); init();
  for(int i=0;i<=n;i++)f[i]=rd();
  for(int i=0;i<=m;i++)g[i]=rd();
  mul(f,g,ans);
  for(int i=0;i<=n+m;i++)printf("%d ",upt(ans[i])); puts("");
  return 0;
}