[关闭]
@xzyxzy 2018-09-26T07:38:53.000000Z 字数 6063 阅读 982

FFT/NTT/MTT

数学

作业部落

评论地址


前言

这是网上的优秀博客
并不建议初学者看我的博客,因为我也不是很了解FFT的具体原理

一、概述

两个多项式相乘,不用,通过可以把复杂度优化到能够取模,可以对非模数取模,相对来说常数小些因为不要取模

二、我们来背板子(FFT)

先放一个板子(洛谷P3803 【模板】多项式乘法(FFT)

  1. #include<iostream>
  2. #include<cstdio>
  3. #include<cstdlib>
  4. #include<cmath>
  5. using namespace std;
  6. const int MAXN=3000005;
  7. const double pi=acos(-1);
  8. int N,M,r[MAXN],l;
  9. struct Complex
  10. {
  11. double rl,im;//real part / imaginary part
  12. Complex(){rl=im=0;}//以下是初始化的板子,虽然不懂为什么可以这样写
  13. Complex(double a,double b){rl=a,im=b;}
  14. Complex operator + (Complex B)
  15. {return Complex(rl+B.rl,im+B.im);}
  16. Complex operator - (Complex B)
  17. {return Complex(rl-B.rl,im-B.im);}
  18. Complex operator * (Complex B)
  19. {return Complex(rl*B.rl-im*B.im,rl*B.im+im*B.rl);}
  20. }A[MAXN],B[MAXN];//对A,B两个多项式进行乘法
  21. void FFT(Complex *P,int op)
  22. {
  23. for(int i=1;i<N;i++)//这个叫Rader排序
  24. /*
  25. 假设原来P[1...n].id=1..n
  26. 现在需要的序列是从1到n所对应的id分别为id[1..n],满足r[id[i]]是升序
  27. r[i]表示把i二进制上第1到l位的数反过来后的十进制数
  28. */
  29. if(i<r[i]) swap(P[i],P[r[i]]);
  30. //接下来的这个叫做蝴蝶操作,算法导论上有一张图较为清晰
  31. for(int i=1;i<N;i<<=1)//表示操作区间集的每个区间的长度
  32. {
  33. Complex W=(Complex){cos(pi/i),op*sin(pi/i)};
  34. for(int p=i<<1,j=0;j<N;j+=p)//表示每个区间集的最上端位置
  35. {
  36. Complex w=(Complex){1,0};//第0个单位复数根
  37. /*
  38. 转角公式:将一个点(x,y)绕原点逆时针旋转t后的点是(x*cost-y*sint,x*sint+y*cost)
  39. 用三角函数和差化积公式容易得证
  40. 单位复数根是把单位元分为若干等份,于是每次就要转一定角度
  41. 用w=w*W实现转角
  42. */
  43. for(int k=0;k<i;k++,w=w*W)//每个区间的最上端位置
  44. {
  45. Complex X=P[j+k],Y=w*P[j+k+i];//j+k+i便是每个区间下端位置
  46. P[j+k]=X+Y;P[j+k+i]=X-Y;//所谓蝴蝶操作
  47. }
  48. }
  49. }
  50. }
  51. int main()
  52. {
  53. cin>>N>>M;
  54. for(int i=0;i<=N;i++) cin>>A[i].rl;
  55. for(int i=0;i<=M;i++) cin>>B[i].rl;
  56. //读入实部,便是系数
  57. M+=N;//最终位数
  58. for(N=1;N<=M;N<<=1) l++;l--;//FFT必须是2^k项才能做,这里把他补全
  59. for(int i=0;i<N;i++) r[i]=(r[i>>1]>>1)|((i&1)<<l);//r是rader排序,将每个i的二进制位反过来
  60. FFT(A,1);FFT(B,1);//将AB化成点集形式
  61. //形如(w0,y0),(w1,y1)...(wn,yn)的这些点确定一条线
  62. for(int i=0;i<N;i++) A[i]=A[i]*B[i];//点集O(n)相乘
  63. FFT(A,-1);//再将点集转化为系数表示的形式
  64. for(int i=0;i<=M;i++) printf("%d ",(int)(A[i].rl/N+0.5));//这时虚部都是0了
  65. return 0;
  66. }

以下是预处理单位复数根的代码
代码长度会小些,精度也要高,建议使用这种写法
三角函数比乘法慢

  1. #include<iostream>
  2. #include<cstdio>
  3. #include<cstdlib>
  4. #include<cmath>
  5. #include<complex>
  6. using namespace std;
  7. const int MAXN=3e6+10;
  8. const double pi=acos(-1);
  9. int r[MAXN],N,M,l;
  10. complex<double>A[MAXN],B[MAXN],w[MAXN];
  11. void FFT(complex<double> *P,int op)
  12. {
  13. for(int i=1;i<N;i++) if(i>r[i]) swap(P[i],P[r[i]]);
  14. for(int i=1;i<N;i<<=1)
  15. for(int p=i<<1,j=0;j<N;j+=p)
  16. for(int k=0;k<i;k++)
  17. {
  18. complex<double> W=w[N/i*k];W.imag()*=op;//实际要得到的是cos(pi/i*k)
  19. complex<double> X=P[j+k],Y=W*P[j+k+i];//QAQ这里总是忘记乘W
  20. P[j+k]=X+Y;P[j+k+i]=X-Y;
  21. }
  22. }
  23. int main()
  24. {
  25. cin>>N>>M;
  26. for(int i=0;i<=N;i++) cin>>A[i].real();
  27. for(int i=0;i<=M;i++) cin>>B[i].real();
  28. M+=N;
  29. for(N=1;N<=M;N<<=1) l++;l--;
  30. for(int i=0;i<N;i++) r[i]=(r[i>>1]>>1)|((i&1)<<l);
  31. for(int i=0;i<N;i++) w[i].real()=cos(pi/N*i),w[i]=imag()=sin(pi/N*i);
  32. FFT(A,1);FFT(B,1);
  33. for(int i=0;i<N;i++) A[i]=A[i]*B[i];
  34. FFT(A,-1);
  35. for(int i=0;i<=M;i++) printf("%d ",(int)(A[i].real()/N+0.5));
  36. puts("");return 0;
  37. }

记忆方式:
循环的枚举当前处理的长度
枚举第几组(两组两组进行)
枚举位置
于是表示某组的第一小组的一个位置,是某组第二小组与第一小组对应的位置
然后先加再减,记得乘上

注意点:
1.最后要(int)(real()/N+0.5)
2.由于N要放大所以空间开两倍!!

三、我们再来背板子(NTT)

还是那道题

  1. #include<iostream>
  2. #include<cstdio>
  3. #include<cstdlib>
  4. #include<cmath>
  5. using namespace std;
  6. const int N=3000005;
  7. const int mod=998244353;
  8. int r[N],l,n,m,A[N],B[N],w[N];
  9. int ksm(int a,int k)
  10. {
  11. int s=1,b=a;
  12. for(;k;k>>=1,b=1ll*b*b%mod)
  13. if(k&1) s=1ll*s*b%mod;
  14. return s;
  15. }
  16. void NTT(int *P,int op)
  17. {
  18. for(int i=0;i<n;i++) if(i<r[i]) swap(P[i],P[r[i]]);
  19. for(int i=1;i<n;i<<=1)
  20. {
  21. int W=ksm(3,(mod-1)/(i<<1));//3是998244353的一个原根
  22. if(op<0) W=ksm(W,mod-2);w[0]=1;
  23. for(int j=1;j<i;j++) w[j]=1ll*w[j-1]*W%mod;
  24. for(int j=0,p=i<<1;j<n;j+=p)
  25. for(int k=0;k<i;k++)
  26. {
  27. int X=P[j+k],Y=1ll*w[k]*P[i+j+k]%mod;
  28. P[j+k]=(X+Y)%mod;P[i+j+k]=((X-Y)%mod+mod)%mod;
  29. }
  30. }
  31. }
  32. int main()
  33. {
  34. cin>>n>>m;
  35. for(int i=0;i<=n;i++) cin>>A[i];
  36. for(int i=0;i<=m;i++) cin>>B[i];
  37. m=n+m;for(n=1;n<=m;n<<=1) l++;l--;
  38. for(int i=0;i<n;i++) r[i]=(r[i>>1]>>1)|((i&1)<<l);
  39. NTT(A,1);NTT(B,1);
  40. for(int i=0;i<n;i++) A[i]=1ll*A[i]*B[i]%mod; NTT(A,-1);
  41. for(int i=0,inv=ksm(n,mod-2);i<n;i++) P[i]=1ll*P[i]*inv%mod;
  42. for(int i=0;i<=m;i++) printf("%d ",A[i]);
  43. return 0;
  44. }

一个数的原根满足各不相同且
对于且仅对于有原根存在
NTT的原根就代替了FFT中的单位复数根,要求形式是
常用的模数有

找质数的原根

最暴力的方法是枚举原根,然后判断是否相同
优化的话是检查的所有质因数中,是否存在一个质因子使得,若存在,则该数不是原根,否则是原根

证明(Thanks GXY)

首先可以明确的是,若对于属于,没有,则g是一个原根
因为,且,则一定有
利用反证法,假设存在一个使得

分两种情况讨论:
1.
为质数
,能够通过上述方法判定出来

2.

,由于,根据同余方程的EXGCD判断,可以在任意取值,都有符合条件的使得式子成立
根据欧拉定理/费马小定理得,使得所有的属于都模p余1,也会在之前的方法中判断出来

四、有个可以讲清的了(MTT)

处理任意模数问题
(这样子好像复杂度最优)
然后多项式的每一项拆成,于是都在之内就不会爆
所以两个数相乘就成为了

分别进行即可(一共8次,有些博客是7次,但是代码比我长)

Code

洛谷P4245 【模板】任意模数NTT

  1. #include<iostream>
  2. #include<cstdio>
  3. #include<cstring>
  4. #include<complex>
  5. #include<cmath>
  6. using namespace std;
  7. const double Pi=acos(-1);
  8. const int N=400100;
  9. const int M=30000;
  10. int n,m,p,F[N],G[N];
  11. int r[N],Ans[N],l,tt;
  12. complex<double> A1[N],B1[N],A2[N],B2[N],A[N],w[N];
  13. void FFT(complex<double> *P,int op)
  14. {
  15. for(int i=0;i<l;i++) if(r[i]<i) swap(P[i],P[r[i]]);
  16. for(int i=1;i<l;i<<=1)
  17. for(int p=i<<1,j=0;j<l;j+=p)
  18. for(int k=0;k<i;k++)
  19. {
  20. complex<double> W=w[l/i*k];W.imag()*=op;
  21. complex<double> X=P[j+k],Y=W*P[j+k+i];
  22. P[j+k]=X+Y;P[j+k+i]=X-Y;
  23. }
  24. }
  25. void Work(complex<double> *P1,complex<double> *P2,int base)
  26. {
  27. for(int i=0;i<l;i++) A[i]=P1[i]*P2[i];FFT(A,-1);
  28. for(int i=0;i<=m+n;i++) (Ans[i]+=(long long)(A[i].real()/l+0.5)%p*base%p)%=p;
  29. }
  30. int main()
  31. {
  32. scanf("%d%d%d",&n,&m,&p);
  33. for(int i=0,x;i<=n;i++) scanf("%d",&x),A1[i].real()=x/M,B1[i].real()=x%M;
  34. for(int i=0,x;i<=m;i++) scanf("%d",&x),A2[i].real()=x/M,B2[i].real()=x%M;
  35. for(l=1;l<=n+m;l<<=1) tt++;tt--;
  36. for(int i=0;i<l;i++) r[i]=(r[i>>1]>>1)|((i&1)<<tt);
  37. for(int i=0;i<l;i++) w[i].real()=cos(Pi/l*i),w[i].imag()=sin(Pi/l*i);
  38. FFT(A1,1);FFT(A2,1);FFT(B1,1);FFT(B2,1);
  39. Work(A1,A2,M*M%p); Work(A1,B2,M%p);
  40. Work(A2,B1,M%p); Work(B1,B2,1);
  41. for(int i=0;i<=m+n;i++) printf("%d ",Ans[i]);
  42. }

五、一些要点

这一部分还没玩成,待博主把这些算法完全弄懂后再来填坑~
NTT时(X-Y+mod)%mod时,Y为负数就可能爆int,可以不加mod然后最后输出的时候加
当乘起来不会超过mod(注意是乘后累加),那么NTT可以代替FFT,否则不行,例子见MTT
乘法通过原根变成加法再NTT
字符串匹配问题的两种做法
组合数公式给拆成可以NTT的形式

添加新批注
在作者公开此批注前,只有你和作者可见。
回复批注