FFT（快速傅里叶变换）摘要
怎么说呢。。。这次的代码码风有点。。。
从这篇博客开始，我终于会用LATEX了！撒花
注：以下涉及多项式的n均表示多项式的次数
FFT用途
很简单，多项式相乘。
FFT原理
如果暴力乘，复杂度是$O(n^{2})$的，显然不行
所以考虑点值法。点值表示下的多项式相乘是$O(n)$的。这就是DFT（离散傅里叶变换）。
但很显然，随机带入n个值的复杂度是$O(n^{2})$的。
因此我们考虑以下推导。
设$$F(n)=a_{0}+a_{1}x+a_{2}x^{2}+...+a_{n}x^{n}$$那么我们一定可以将其拆成两个函数，即$$G(n)=a_{0}+a_{2}x+a_{4}x^{2}+a_{6}x^{3}+...$$$$H(n)=a_{1}+a_{3}x+a_{5}x^{2}+a_{7}x^{3}+...$$所以$$F(n)=G(x^{2})+xH(x^{2})$$从这个式子我们可以看出，如果我们取这么一个点值集合，使得其中一半的值为所有值的平方，比如${-1,1}$，那么我们就可以使计算规模减半（对于$G(n)$和$H(n)$，我们只需带入其中一半的值即可）。看来这是个很完美的策略。
很不幸，如果我们只在实数中这么思考，那复杂度仍为$O(n^{2})$，可以说没什么用。但这给了我们一个启发，如果每次对函数进行拆分递归并使点值集合规模减半，我们就可以在$O(nlogn)$的复杂度内求出一个多项式的点值表示。毫无疑问，问题的关键是这个点值集合。那是否存在这样一个点值集合呢？实数域当然没有，但复数域有，也就是单位根！（用复数做自变量真的让我很长时间无法接受）至于单位根么，出门右转问百度
当然，既然用了FFT，自然要把多项式项数补成2的幂，具体做法为高次项补零。
递归式FFT的很好打，但常数大，且容易爆栈。因此我们考虑迭代。
易发现，每次拆分多项式的时候，我们总是按下标的奇偶性来拆分。这就导致了一个结果：原多项式经过递归拆分后的新多项式下标是原多项式下标的二进制反转。
1 void unit() { 2 lg[1]=0; 3 fui(i,2,mxle,1)lg[i]=lg[i/2]+1; 4 bz[0]=1; 5 fui(i,1,mxlg,1)bz[i]=bz[i-1]<<1; 6 n=max(n,m); 7 n=bz[lg[n]+2]-1; 8 root[0].x=1; 9 root[1].x=cos((2*PI)/(n+1)); 10 root[1].y=sin((2*PI)/(n+1)); 11 fui(i,2,n,1)root[i]=root[i-1]*root[1]; //求单位根 12 frot[0]=root[0]; 13 fui(i,1,n,1)frot[i]=root[n-i+1]; //求逆单位根（后面用 14 for(int i=0; i<=n; i++)fz[i]=(fz[i>>1]>>1)|((i&1)<<(lg[n+1]-1)) ; //二进制反转 15 }
然后就是FFT了
1 void pre_FFT() { 2 fui(i,0,n,1)if(i<fz[i]) { 3 swp=fft[i]; 4 fft[i]=fft[fz[i]]; 5 fft[fz[i]]=swp; 6 } 7 } 8 void FFT() { 9 pre_FFT(); //对原多项式进行二进制反转 10 int zl=(n+1)>>1; //单位根增量 11 for(int i=1; i<=n; i=(i<<1)|1,zl>>=1) { //每次处理的多项式长度 12 for(int j=0; j<n; j+=i+1) { //每次处理的多项式左端点 13 int mid=(j+j+i)>>1; 14 for(int k1=j,k2=mid+1,nw=0; k1<=mid; nw+=zl,k1++,k2++) { //对该多项式的点值进行迭代处理 15 x1=fft[k2]*root[nw]+fft[k1],x2=fft[k2]*root[nw+((n+1)>>1)]+fft[k1]; 16 fft[k1]=x1; 17 fft[k2]=x2; 18 } 19 } 20 } 21 }
记住，千万不要用stl里的复数，会被坑死的，强烈建议手打
IFFT（快速傅里叶逆变换）
FFT求出了点值，然后相乘，求出了结果多项式的点值表示。但这不是结果，我们还要将其“翻译”回系数表示。所以就有了IFFT，可以叫其插值。
插值么，具体做法就是对结果求一遍类似FFT的东西，只不过把单位根改成{${{\omega}^{0},{\omega}^{-1},{\omega}^{-2},{\omega}^{-3},...,{\omega}^{-n}}$}，即上面预处理的时候求过的那个所谓逆单位根带入即可。证明么，我不会，反正作为一个OIer记结论就好了。具体证明参见自为风月马前卒大佬的博客。
1 void IFFT() { 2 pre_FFT(); 3 int zl=(n+1)>>1; 4 for(int i=1; i<=n; i=(i<<1)|1,zl>>=1) { 5 for(int j=0; j<n; j+=i+1) { 6 int mid=(j+j+i)>>1; 7 for(int k1=j,k2=mid+1,nw=0; k1<=mid; nw+=zl,k1++,k2++) { 8 x1=fft[k2]*frot[nw]+fft[k1],x2=fft[k2]*frot[nw+((n+1)>>1)]+fft[k1]; //改成逆单位根带入 9 fft[k1]=x1; 10 fft[k2]=x2; 11 } 12 } 13 } 14 }
模板题
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define INF 0x7fffffff 4 #define ME 0x7f 5 #define FO(s) freopen(s".in","r",stdin);freopen(s".out","w",stdout) 6 #define fui(i,a,b,c) for(int i=(a);i<=(b);i+=(c)) 7 #define fdi(i,a,b,c) for(int i=(a);i>=(b);i-=(c)) 8 #define fel(i,a) for(register int i=h[a];i;i=b[i].ne) 9 #define ll long long 10 #define ld double 11 #define MEM(a,b) memset(a,b,sizeof(a)) 12 #define maxn 1000010 13 #define mxlg 21 14 #define mxle 4194310 15 const ld PI=3.1415926535897932384626; 16 int n,m,mn; 17 ld a0[mxle],b0[mxle]; 18 ld ans[mxle]; 19 int lg[mxle],bz[mxlg+3]; 20 int fz[mxle]; 21 struct Complex{ 22 ld x,y; 23 Complex(){} 24 Complex(ld __x,ld __y){x=__x;y=__y;} 25 Complex operator *(Complex xx)const{return (Complex){x*xx.x-y*xx.y,x*xx.y+y*xx.x};} 26 Complex operator *(ld xx)const{return (Complex){x*xx,y*xx};} 27 Complex operator +(Complex xx)const{return (Complex){x+xx.x,y+xx.y};} 28 }root[mxle],frot[mxle],a[mxle],b[mxle],fft[mxle],c[mxle],swp,x1,x2; 29 template<class T> 30 inline T read(T &n){ 31 n=0;int t=1;double x=10;char ch; 32 for(ch=getchar();!isdigit(ch)&&ch!='-';ch=getchar());(ch=='-')?t=-1:n=ch-'0'; 33 for(ch=getchar();isdigit(ch);ch=getchar()) n=n*10+ch-'0'; 34 if(ch=='.') for(ch=getchar();isdigit(ch);ch=getchar()) n+=(ch-'0')/x,x*=10; 35 return (n*=t); 36 } 37 template<class T> 38 T write(T n){ 39 if(n<0) putchar('-'),n=-n; 40 if(n>=10) write(n/10);putchar(n%10+'0');return n; 41 } 42 template<class T> 43 T writeln(T n){write(n);putchar('\n');return n;} 44 int inc(int &x){int in=(n+1)>>1;for(;x∈in>>=1)x^=in;x|=in;return x;} 45 void unit(){ 46 lg[1]=0;fui(i,2,mxle,1)lg[i]=lg[i/2]+1; 47 bz[0]=1;fui(i,1,mxlg,1)bz[i]=bz[i-1]<<1; 48 n=max(n,m);n=bz[lg[n]+2]-1; 49 root[0].x=1; 50 root[1].x=cos((2*PI)/(n+1));root[1].y=sin((2*PI)/(n+1)); 51 fui(i,2,n,1)root[i]=root[i-1]*root[1]; 52 frot[0]=root[0]; 53 fui(i,1,n,1)frot[i]=root[n-i+1]; 54 fui(i,0,n,1)a[i]=root[0]*a0[i]; 55 fui(i,0,n,1)b[i]=root[0]*b0[i]; 56 for(int i=0;i<=n;i++)fz[i]=(fz[i>>1]>>1)|((i&1)<<(lg[n+1]-1)) ; 57 } 58 void pre_FFT(){ 59 fui(i,0,n,1)if(i<fz[i]){ 60 swp=fft[i];fft[i]=fft[fz[i]];fft[fz[i]]=swp; 61 } 62 } 63 void FFT(){ 64 pre_FFT();int zl=(n+1)>>1; 65 for(int i=1;i<=n;i=(i<<1)|1,zl>>=1){ 66 for(int j=0;j<n;j+=i+1){ 67 int mid=(j+j+i)>>1; 68 for(int k1=j,k2=mid+1,nw=0;k1<=mid;nw+=zl,k1++,k2++){ 69 x1=fft[k2]*root[nw]+fft[k1],x2=fft[k2]*root[nw+((n+1)>>1)]+fft[k1]; 70 fft[k1]=x1;fft[k2]=x2; 71 } 72 } 73 } 74 // pre_FFT(); 75 } 76 void IFFT(){ 77 pre_FFT();int zl=(n+1)>>1; 78 for(int i=1;i<=n;i=(i<<1)|1,zl>>=1){ 79 for(int j=0;j<n;j+=i+1){ 80 int mid=(j+j+i)>>1; 81 for(int k1=j,k2=mid+1,nw=0;k1<=mid;nw+=zl,k1++,k2++){ 82 x1=fft[k2]*frot[nw]+fft[k1],x2=fft[k2]*frot[nw+((n+1)>>1)]+fft[k1]; 83 fft[k1]=x1;fft[k2]=x2; 84 } 85 } 86 } 87 } 88 int main(){ 89 read(n);read(m);mn=m+n; 90 fui(i,0,n,1)read(a0[i]);fui(i,0,m,1)read(b0[i]);unit(); 91 fui(i,0,n,1)fft[i]=a[i];FFT();fui(i,0,n,1)a[i]=fft[i]; 92 fui(i,0,n,1)fft[i]=b[i];FFT();fui(i,0,n,1)b[i]=fft[i]; 93 fui(i,0,n,1)c[i]=a[i]*b[i]; 94 fui(i,0,n,1)fft[i]=c[i];IFFT();fui(i,0,n,1)c[i]=fft[i]; 95 fui(i,0,mn,1)ans[i]=(c[i].x+0.5)/((n+1)); 96 fui(i,0,mn,1)printf("%.0f ",ans[i]); 97 return 0; 98 }
常数太大，导致最后两个点2s多。。。
