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 }

FFT

記住，千萬不要用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 }

IFFT

模板題

洛谷P3803 【模板】多項式乘法（FFT）

 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;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 }

AC代碼

常數太大，致使最後兩個點2s多。。。