快速傅里叶变换(FFT)模板
//有多少项N取多至少8~10倍 #include <algorithm> #include <cmath> #include <complex> #include <cstdlib> using namespace std; typedef complex<double> cn; const int N=1e7+1; const double Pi=3.1415926535;//Pi通常取10位即可(随意) cn a[N],b[N]; int r[N]; int mxb,bit; int n,m; //二进制翻转 void Get_Rev(int bit) { for (int i=0;i<(1<<bit);i++) r[i]=(r[i>>1]>>1)|((i&1)<<(bit-1)); } //傅里叶变换 void DFT(cn *a,int mx,int inot) {//inot为1就是DFT,为-1就是IDFT for (int i=0;i<mx;i++) if (i<r[i]) swap(a[i],a[r[i]]); for (int mlen=1;mlen<mx;mlen<<=1) { cn wn=exp(cn(0,inot*Pi/mlen)); for (int len=mlen<<1,l=0;l<mx;l+=len) { cn wnk=cn(1,0); for (int k=l;k<l+mlen;k++,wnk*=wn) { cn x=a[k],y=wnk*a[k+mlen]; a[k]=x+y;a[k+mlen]=x-y;//蝴蝶操作 } } } if (inot==-1)//IDFT的除n操作 for (int i=0;i<mx;i++) a[i]/=mx; } void FFT() { mxb=1; while (mxb<=n+m) mxb<<=1,bit++; Get_Rev(bit); DFT(a,mxb,1);DFT(b,mxb,1); for (int i=0;i<=mxb;i++) a[i]*=b[i]; DFT(a,mxb,-1); }
在日渐沉没的世界里,我发现了你。