FFT，即快速傅里葉變換，是離散傅里葉變換的快速方法，可以在很低復雜度內解決多項式乘積的問題（兩個序列的卷積）
卷積
卷積通俗來說就一個公式（本人覺得卷積不重要）
$$C_k=\sum_{i+j=k}A_i*B_i$$
那么這個表達式是啥意思了：
有兩個序列A和B，其中$A=\left\{A_1,A_2,...\right\}，B=\left\{B_1,B_2,...\right\}$
A序列有a個元素，B序列有b個元素。於是，由這兩個序列可以推出另一個序列C，C序列如何確定了？確定方法就按照卷積的公式得來的，即$$C=\left\{\sum_{i+k=0}A_i*B_j\,,\,\sum_{i+k=1}A_i*B_j\,,...,\sum_{i+k=2^a+2^b}A_i*B_j\right\}$$
這里只介紹一下卷積，下面來探究卷積和多項式之間的關系。
多項式（預備知識）
多項式的定義
形如
$$f(x)=a_0+a_1x+a_2x^2+...+a_nx^n$$
的函數關系式叫做關於x的多項式，其中多項式系數可以構成一個序列，即
$$A=\left\{a_0,a_1,a_2,...,a_n\right\}$$
多項式乘法與序列的卷積
假如有兩個多項式，其中
$f(x)=a_0+a_1x+a_2x^2+...+a_nx^n$
$g(x)=b_0+b_1x+b_2x^2+...+b_mx^m$
現在要求f(x)*g(x)的序列，很明顯
$$f(x)*g(x)=\sum_{i+j=0}a_i*b_j+\sum_{i+j=1}a_i*b_jx+\sum_{i+j=2}a_i*b_jx^2+...+\sum_{i+j=m+n}a_i*b_jx^{m+n}$$
現在把$f(x),g(x)$以及$f(x)*g(x)$三個多項式的系數拿出來寫成三個序列，可得:
$A=\left\{a_0,a_1,a_2,...a_n\right\}$
$B=\left\{b_0,b_1,b_2,...b_m\right\}$
$C=\left\{\sum_{i+j=0}a_i*b_j,\sum_{i+j=1}a_i*b_j,\sum_{i+j=2}a_i*b_j,...\sum_{i+j=m+n}a_i*b_j\right\}$
於是驚訝的發現，兩個序列的卷積相當於兩個多項式乘積后得到的多項式的系數序列。而FFT算法的目的，就是求兩個多項式乘積后得到的多項式的系數序列。
多項式的表示方法
多項式有兩種表示方法，即系數表示法和點值表示法。
1.系數表示法是我們常用的表示方法，即
$$f(x)=a_0+a_1x+a_2x^2+...+a_nx^n$$
的形式
2.點值表示法，顧名思義，就是在這個多項式上任意選取不重復的n+1個點，即
$$f(x)=\left\{(x_0,y_0),(x_1,y_1),...,(x_2,y_2)\right\}$$
可以證明：任何n+1個點可以確定唯一一個最高次項為n次方的多項式（下面是證明，可以看看，也可以忽略）
將一個多項式的系數寫成一個系數矩陣（當然，這些系數我們是不知道的）
$$\begin{bmatrix}a_0\\a_1\\a_2\\\vdots\\a_n\end{bmatrix}$$
然后將剛剛選取的n+1個點寫成下面兩個矩陣
$$\begin{bmatrix}1&x_0&x_0^2&x_0^3&\cdots&x_0^n\\1&x_1&x_1^2&x_1^3&\cdots&x_1^n\\1&x_2&x_2^2&x_2^3&\cdots&x_2^n\\\vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\1&x_n&x_n^2&x_n^3&\cdots&x_n^n\end{bmatrix}\;\begin{bmatrix}y_0\\y_1\\y_2\\\vdots\\y_n\end{bmatrix}$$
可以得出以下三個矩陣的關系
$$\begin{bmatrix}1&x_0&x_0^2&x_0^3&\cdots&x_0^n\\1&x_1&x_1^2&x_1^3&\cdots&x_1^n\\1&x_2&x_2^2&x_2^3&\cdots&x_2^n\\\vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\1&x_n&x_n^2&x_n^3&\cdots&x_n^n\end{bmatrix}\;\begin{bmatrix}a_0\\a_1\\a_2\\\vdots\\a_n\end{bmatrix}=\begin{bmatrix}y_0\\y_1\\y_2\\\vdots\\y_n\end{bmatrix}$$
解出系數矩陣
$$\begin{bmatrix}a_0\\a_1\\a_2\\\vdots\\a_n\end{bmatrix}=\begin{bmatrix}1&x_0&x_0^2&x_0^3&\cdots&x_0^n\\1&x_1&x_1^2&x_1^3&\cdots&x_1^n\\1&x_2&x_2^2&x_2^3&\cdots&x_2^n\\\vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\1&x_n&x_n^2&x_n^3&\cdots&x_n^n\end{bmatrix}^{-1}\begin{bmatrix}y_0\\y_1\\y_2\\\vdots\\y_n\end{bmatrix}$$
只能唯一確定一個系數矩陣，即只能唯一確定一組系數，即只能唯一確定一個多項式，證明完畢
多項式乘法與FFT算法
為了實現多項式乘法並得到系數序列，我們可以考慮實現的方法，如果直接暴力（通過定義直接算系數）是O(n2)復雜度，肯定會超時，於是我們這樣考慮
首先可以選取n+1個點，把兩個多項式轉化為點值表示法。
然后將兩個點值表示法的多項式相乘（復雜度為O(n)）,然后將新得到的多項式的點值表示法轉化為系數表示法。
注意：假設f(x)最高有n次方，g(x)最高有m次方，所以f(x)*g(x)最高有m+n次方，但是m+n可能不是2的冪次方，如果不是，則需要選取點的數量應該是大於m+n的2的冪次方個，假設這個值為1<<k，所以從系數表示法到點值表示法的轉化過程中，我們要在f(x)和g(x)內選擇1<<k的點，才能保證點值表示法的f(x)和g(x)之后，至少仍有m+n個點，才能確定唯一一個f(x)*g(x)的多項式。所以在選取點的數量的時候，要保證點的個數能確定f(x)*g(x)后的多項式。
如果任意的選取n+m個點然后轉化，復雜度還是太高，所以我們需要巧妙的選取1<<k個點，從而使得系數表示法與點值表示法之間轉化的復雜度降為O(nlogn)，這就是FFT算法。
虛數的乘積及其表示（預備知識）
一個虛數，是由實部和虛部組成，表示為$(a+bi)$，虛數的幾何表示可以在坐標系中體現，如下圖所示
如圖所示，$r$為此虛數的模長，$\theta$為此虛數的幅角，於是虛數還能表示為$r(cos\theta,isin\theta)$的形式
若$z_1=r_1(cos\theta_1,isin\theta_1),z_2=r_2(cos\theta_2,isin\theta_2)$
於是$$z_1*z_2=r_1(cos\theta_1,isin\theta_1)*r_2(cos\theta_2,isin\theta_2)=r_1r_2(cos(\theta_1+\theta_2),isin(\theta_1+\theta_2))$$
可以得到虛數相乘的幾何意義：模長相乘，幅角相加
n次單位根（預備知識）
定義
n次單位根$w_n$，即$x^n=1$在虛數范圍內的解，即$w_n^n=1$，故n次單位根$w_n$也是一個復數，且模長為1，所以n次單位根在乘方的時候模長不變，幅角等差增大
如圖所示，將一個單位圓分成n塊，幅角為正且最小的那個虛數即為n次單位根（$w_n^1$，紅色點），下圖是n=8的情況
還要注意，$w_n^0=w_n^n=1$
由圖及n次單位根的性質（n次單位根模長為1，幅角為$\frac{2\pi}{n}$）可得
$$w_n^k=cos\frac{2k\pi}{n}+i\,sin\frac{2k\pi}{n}$$
性質
比較簡單的性質：
1.$w_n^{m+n}=w_n^m*w_n^n=w_n^m$
2.$(w_n^m)^n=(w_n^n)^m=1,m\in [0,n-1]$
比較復雜（重要）的性質：
1.$w_{2n}^{2k}=w_n^k$，如圖所示
2.$w_n^k=-w_n^{k+\frac n{2}}$
通過上面那個圖也可以看出來，注意一點：虛數乘-1，那么虛數的實部與虛部都要乘上-1，所以$-w_n^k$在單位圓上的點和$w_n^k$在單位圓上的點關於原點對稱
3.$(w_n^k)^2=w_{\frac n{2}}^k$
推導:$(w_n^k)^2=(-w_n^{k+\frac n{2}})^2=w_n^{2*k}=w_{\frac n{2}}^k$
更加復雜（重要）的性質
$$\sum_{k=0}^{n-1}w_n^{ik}=\begin{cases} 0 & i\bmod n \ne 0 \\ n & i\bmod n=0 \end{cases}$$
證明：
當$i\bmod n \ne 0$時，相當如等比數列前n項和求和，故$$\sum_{k=0}^{n-1}w_n^{ik}=\frac {w_n^0(1-(w_n^i)^n)}{1-w_n^i}=\frac {1-(w_n^n)^i}{1-w_n^i}=0$$
當$i\bmod n = 0$時，$w_n^{ik}=1$，故$$\sum_{k=0}^{n-1}w_n^{ik}=n*1=n$$
至此，你應該知道，多項式從系數表示法轉化為點值表示法時，選取的n個點的x值即為n次單位根序列$\left\{w_n^0,w_n^1,....,w_n^{n-1}\right\}$
注意，此處的n是f(x)*g(x)的多項式最高次項的次數，不是f(x)或者g(x)最高此項的系數，前面（多項式乘法與FFT）講過要選取1<<k（再次強調這個值比f(x)*g(x)的多項式最高次的次數m+n還要大，且是2的冪次方）個點
FFT
假設現在多項式為
$f(x)=a_0+a_1x+a_2x^2+...+a_{k-1}x^{k-1}$，注意，多項式最高次項為k-1次方，共k項，k如果不等於2的n次冪，就湊成2的n次冪（加系數為0的項）
$g(x)=b_0+b_1x+b_2x^2+...+b_{h-1}x^{h-1}$，注意，多項式最高次項為h-1次方，共h項，h如果不等於2的n次冪，就湊成2的n次冪（加系數為0的項）
那么現在，我們從多項式中選取了這樣一些點
${(w_n^0,f(w_n^0)),(w_n^1,f(w_n^1)),(w_n^2,f(w_n^2)),...,(w_n^{n-1},f(w_n^{n-1}))}$，共n個點
注意n是大於k+h的最小2次冪，而不是等於k+h！！！！！
此時選取得點數n大於k+h，故可以確定f(x)*g(x)這個多項式，現在，只需確定$f(w_n^0),f(w_n^1),...,f(w_n^{n-1})$的值即可。
如何算值？我們可以觀察一下$f(w_n^i)$的值$(0\le i<n)$
$f(w_n^i)=a_0+a_1w_n^i+a_2(w_n^i)^2+a_3(w_n^i)^3+...+a_{k-1}(w_n^i)^{k-1}$
$=\color{Red}{a_0+a_2(w_n^i)^2+a_4(w_n^i)^4+...+a_{k-2}(w_n^i)^{k-2}}+\color{LimeGreen}{w_n^i}\color{Blue}{(a_1+a_3(w_n^i)^2+...+a_{k-1}(w_n^i)^{k-2})}$
根據$(w_n^k)^2=w_{\frac n{2}}^k$可得
$=\color{Red}{a_0+a_2w_{\frac n{2}}^i+a_4(w_{\frac n{2}}^i)^2+...+a_{k-2}(w_{\frac n{2}}^i)^{\frac {k-2}2}}+\color{LimeGreen}{w_n^i}\color{Blue}{(a_1+a_3w_{\frac n{2}}^i+...+a_{k-1}(w_{\frac n{2}}^i)^{\frac {k-2}2})}$
$=\color{Red}{f_1(w_n^i)}+\color{LimeGreen}{w_n^i}\color{Blue}{f_2(w_n^i)}$
再算$f_1(w_n^i)$的值和$f_2(w_n^i)$的值的時候，我們又可以按照奇偶分開，像算$f(w_n^i)$一樣變成兩個問題
同時我們發現
$f(w_n^{i+\frac n{2}})=f(-w_n^i)$
$=\color{Red}{a_0+a_2(-w_n^i)^2+a_4(-w_n^i)^4+...+a_{k-2}(-w_n^i)^{k-2}}-\color{LimeGreen}{w_n^i}\color{Blue}{(a_1+a_3(-w_n^i)^2+...+a_{k-1}(-w_n^i)^{k-2})}$
$=\color{Red}{a_0+a_2w_{\frac n{2}}^i+a_4(w_{\frac n{2}}^i)^2+...+a_{k-2}(w_{\frac n{2}}^i)^{\frac {k-2}2}}-\color{LimeGreen}{w_n^i}\color{Blue}{(a_1+a_3w_{\frac n{2}}^i+...+a_{k-1}(w_{\frac n{2}}^i)^{\frac {k-2}2})}$
$=\color{Red}{f_1(w_n^i)}-\color{LimeGreen}{w_n^i}\color{Blue}{f_2(w_n^i)}$
又發現，如果我們要算$f(w_n^i)$的值，我們會先算$f_1(w_n^i)$和$f_2(w_n^i)$的值，但是算出$f_1(w_n^i)$和$f_2(w_n^i)$的值之后，不僅能算出$f(w_n^i)$的值，還能同時算出$f(w_n^{i+\frac n{2}})$的值
總的來說，如果n=8算出了$f(w_n^0)$的值，就算出了$f(w_n^4)$的值；算出了$f(w_n^1)$，就算出來$f(w_n^5)$的值；...；算出了$f(w_n^3)$的值，就算出了$f(w_n^7)$的1值
同理，算$f_1(w_n^i)$和$f_2(w_n^i)$的值的時候，也是算出一半，得到另一半的值------每次解決問題，都會變成解決一半問題然后直接得到另一半的答案，在解決這一半問題的時候，也變成了解決一半的一半的問題，另一半的一半的問題又被直接得到答案。
這樣，就可以把$f(w_n^0),f(w_n^1),f(w_n^2),...,f(w_n^{n-1})$的值全部算出來，得到了n個點的值；同理算出$g(w_n^0),g(w_n^1),f(w_n^2),...,g(w_n^{n-1})$的值。得到了f(x)和g(x)的點值表示法，點值表示法相乘，O(n)復雜度就得到了f(x)*g(x)的點值表示法。
代碼如下：
#include <complex>//c++自帶復數模板 typedef complex<double> Complex;//重定義數據類型，記得要用double類型 void FFT(Complex *a,int n,int inv){ //a是你要進行變換的數組（多項式系數數組），n是數組大小，inv暫時不管，你就當做它等於1 if(n==1) return;//如果分治后的序列大小為1，就直接返回，不用繼續分治了 int mid=n>>1;//如果n不等於1，那就繼續分治，mid是分治后變成兩個序列的長度 static Complex b[1000];//創建一個輔助用的b數組，后面體現用處 for(int i=0;i<mid;i++) b[i]=a[2*i],b[i+mid]=a[2*i+1];//將a數組奇數位置的值和偶數位置的值分開，b數組前mid個位置存a數組偶數位置值，后mid個位置存a數組奇數位置值 for(int i=0;i<n;i++) a[i]=b[i];//將b數組值重新賦給a數組 //對分治后的兩個序列（一個是a數組前mid個元素組成的序列，另一個數a數組后mid元素組成的序列） 進行FFT變換 FFT(a,mid,inv); FFT(a+mid,mid,inv); //算出一半，得到另一半的值，a數組前mid元素相當於之前說的f1，a數組后mid元素相當於之前說的f2 for(int i=0;i<mid;i++){ Complex com(cos(2*pi*i/n),inv*sin(2*pi*i/n)); b[i]=a[i]+com*a[i+mid];b[i+mid]=a[i]-com*a[i+mid]; } for(int i=0;i<n;i++) a[i]=b[i];//重新將b數組賦值給a數組 return; }
這里再來詳細解釋一下，上面代碼中的n值
加入開始兩個序列是n次多項式和m次多項式，那么開始n值就等於大於等於$(m+n)$的最小2次冪，如下代碼所示
cin>>n>>m; int k=n+m,lim=1; while(k>>=1) lim<<=1; if(lim<n+m) lim<<=1; FFT(a,lim,1);//對a序列進行變換 FFT(b,lim,1);//對b序列進行變換
此時lim值就是n的初值
強調一遍$$w_n^k=cos\frac{2k\pi}{n}+i\,sin\frac{2k\pi}{n}$$
FFT逆變換
FFT逆變換則是將點值表示法轉化為系數表示法的步驟。在逆變換之前，我們已經對兩個多項式系數序列進行了正變換，即
FFT(a,lim,1);//正變換 FFT(b,lim,1); for(int i=0;i<lim;i++) a[i]=a[i]*b[i];//將正變換的兩個序列乘在一起，得到新的多項式的點值表示法序列。
如上所示，只要將a序列進行逆變換，就可以得到新多項式的系數序列。在FFT正變換之中，我們一直在求$f(w_n^i)$的值。
假設所求多項式系數序列為$k_0,k_1,...k_{n-1}$
現在對於新的多項式a（新得到的序列，不是原來的a序列）來求卷積$$c_h=\sum_{i=0}^{n-1}a_i(w_n^{-h})^i$$
相當於對a序列又進行FFT變換，只不過，我們現在以新序列a為系數的多項式在$x=w_n^{-i}$的一系列值
相當於求$C(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}$在$x=w_n^0,w_n^-1,w_n^{-2},...w_n^{-(n-1)}$的點值表示法
而我們知道$a_x=k_0+k_1w_n^x+k_2(w_n^k)^2+...=\sum_{j=0}^{n-1}k_j(w_n^x)^j$
$$\sum_{i=0}^{n-1}a_i(w_n^{-h})^i$$
$$=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}k_j*(w_n^i)^j)(w_n^{-h})^i$$
$$=\sum_{j=0}^{n-1}k_j(\sum_{i=0}^{n-1}(w_n^{j-h})^i)$$
$$=\sum_{j=0}^{n-1}k_j(\sum_{i=0}^{n-1}w_n^{(j-h)*i})$$
由於當且只當$j=h$時，$(j-h)$才為i的倍數，其余時候為0，根據n次單位根更加復雜的性質可得（不知道的看看n次單位根預備性質）
$$\sum_{j=0}^{n-1}k_j(\sum_{i=0}^{n-1}w_n^{(j-h)*i})$$
$$=n*k_h$$
所以
$$k_h=\frac {\sum_{i=0}^{n-1}a_i(w_n^{-h})^i}{h}=\frac {c_h}{n}$$
說明：$w_n^{-i}$與$w_n^i$實部相同，虛部相反。現在知道代碼中inv的作用了吧
當inv=1，表示FFT變換；當inv=-1，表示FFT逆變換
FFT(a,lim,1); FFT(b,lim,1); for(int i=0;i<lim;i++) a[i]=a[i]*b[i]; FFT(a,lim,-1); //經過FFT逆變換之后，a序列虛部都已經為0，只有實部有值，且為整數 for(int i=0;i<n+m-1;i++) printf("%d ",(int)(a[i].real()/lim+0.5));//實部根據前面推導要除以n（這里除以lim）
為啥要加0.5捏：首先我們確定答案肯定是整數，但是在代碼實現過程中，由於有$\pi$介入計算，導致答案可能變小了（誤差很小），所以加一個0.5，然后強制類型轉換切掉小數部位
比如說本來答案是1，但是代碼實現的結果可能是0.999999999999，此時加上了0.5，然后轉換為int，答案就是1了。
強調一遍$$w_n^k=cos\frac{2k\pi}{n}+i\,sin\frac{2k\pi}{n}$$
FFT代碼
#include <complex>//c++自帶復數模板 typedef complex<double> Complex;//重定義數據類型，記得要用double類型 void FFT(Complex *a,int n,int inv){ //a是你要進行變換的數組（多項式系數數組），n是數組大小，inv=1表示正變換，inv=-1表示逆變換 if(n==1) return;//如果分治后的序列大小為1，就直接返回，不用繼續分治了 int mid=n>>1;//如果n不等於1，那就繼續分治，mid是分治后變成兩個序列的長度 static Complex b[1000];//創建一個輔助用的b數組，后面體現用處 for(int i=0;i<mid;i++) b[i]=a[2*i],b[i+mid]=a[2*i+1];//將a數組奇數位置的值和偶數位置的值分開，b數組前mid個位置存a數組偶數位置值，后mid個位置存a數組奇數位置值 for(int i=0;i<n;i++) a[i]=b[i];//將b數組值重新賦給a數組 //對分治后的兩個序列（一個是a數組前mid個元素組成的序列，另一個數a數組后mid元素組成的序列） 進行FFT變換 FFT(a,mid,inv); FFT(a+mid,mid,inv); //算出一半，得到另一半的值，a數組前mid元素相當於之前說的f1，a數組后mid元素相當於之前說的f2 for(int i=0;i<mid;i++){ Complex com(cos(2*pi*i/n),inv*sin(2*pi*i/n)); b[i]=a[i]+com*a[i+mid];b[i+mid]=a[i]-com*a[i+mid]; } for(int i=0;i<n;i++) a[i]=b[i];//重新將b數組賦值給a數組 return; }
FFT系數儲存注意
在進行FFT之前，我們有兩個多項式系數序列，用數組存系數的不要用普通的數組存，而要用復數數組，因為在計算的時候系數涉及復數計算，所以這樣
#include <complex>//c++自帶復數模板 Complex a[1000],b[1000]; for(int i=0;i<n;i++) scanf("%d",&t),a[i]=Complex(t,0);//構造函數 for(int i=0;i<m;i++) scanf("%d",&t),b[i]=Complex(t,0);
這樣，復數的實部存系數的值，虛部為0
一個算多項式相乘后系數序列的程序
//開頭說明一下，最下面有輸入輸出樣例 #include <iostream> #include <cmath> #define pi acos(-1) #include <complex>//c++自帶復數模板 using namespace std; typedef complex<double> Complex;//重定義數據類型，記得要用double類型 Complex a[1000],b[1000]; int n,m; void FFT(Complex *a,int n,int inv){ //a是你要進行變換的數組（多項式系數數組），n是數組大小，inv=1表示正變換，inv=-1表示逆變換 if(n==1) return;//如果分治后的序列大小為1，就直接返回，不用繼續分治了 int mid=n>>1;//如果n不等於1，那就繼續分治，mid是分治后變成兩個序列的長度 static Complex b[1000];//創建一個輔助用的b數組，后面體現用處 for(int i=0;i<mid;i++) b[i]=a[2*i],b[i+mid]=a[2*i+1];//將a數組奇數位置的值和偶數位置的值分開，b數組前mid個位置存a數組偶數位置值，后mid個位置存a數組奇數位置值 for(int i=0;i<n;i++) a[i]=b[i];//將b數組值重新賦給a數組 //對分治后的兩個序列（一個是a數組前mid個元素組成的序列，另一個數a數組后mid元素組成的序列） 進行FFT變換 FFT(a,mid,inv); FFT(a+mid,mid,inv); //算出一半，得到另一半的值，a數組前mid元素相當於之前說的f1，a數組后mid元素相當於之前說的f2 for(int i=0;i<mid;i++){ Complex com(cos(2*pi*i/n),inv*sin(2*pi*i/n)); b[i]=a[i]+com*a[i+mid];b[i+mid]=a[i]-com*a[i+mid]; } for(int i=0;i<n;i++) a[i]=b[i];//重新將b數組賦值給a數組 return; } int main(){ cin>>n>>m;//輸入兩個多項式的項數 int k=n+m,lim=1,t; for(int i=0;i<n;i++) scanf("%d",&t),a[i]=Complex(t,0);//第一個多項式系數 for(int i=0;i<m;i++) scanf("%d",&t),b[i]=Complex(t,0);//第二個多項式系數 while(k>>=1) lim<<=1; if(lim<n+m) lim<<=1; FFT(a,lim,1); FFT(b,lim,1); for(int i=0;i<lim;i++) a[i]=a[i]*b[i]; FFT(a,lim,-1); for(int i=0;i<n+m-1;i++) printf("%d ",(int)(a[i].real()/lim+0.5)); } /* //輸入 2 3//第一個多項式有兩項 ，第二個多項式有三項 1 2//第一個多項式為(1+2x) 2 1 2//第二個多項式為(2+x+2x^2) //輸出 2 5 4 4//(1+2x)*(2+x+2x^2)=2+5x+4x^2+4x^3,系數為2 5 4 4 */
測試用例
//輸入 2 2//兩個多項式都有兩項 1 1//第一個多項式為(1+x) 1 1//第二個多項式為(1+x) //輸出 1 2 1//(1+x)*(1+x)=1+2*x+x^2,系數為1 2 1
FFT變換的示例演示
在講FFT迭代版本之前，我打算先用一個例子來展示FFT變換，加深認識
假設
多項式$f(x)=3+2x+x^2$
多項式$g(x)=2+x+2x^2$
在草稿紙上算一下，這兩個多項式相乘的結果為$h(x)=6+7x+10x^2+5x^3+2x^4$
先得到$f(x)$與$g(x)$多項式的系數序列，注意應該寫成復數形式的序列（復數形式為(real,image)）
多項式$f(x)$的系數序列為$A=\left\{(3,0),(2,0),(1,0)\right\}$
多項式$g(x)$的系數序列為$B=\left\{(2,0),(1,0),(2,0)\right\}$
然后判斷兩個多項式相乘后有多少項，應該不超過3+3=6項。而大於6的最小2次冪應該為8，所以要把兩個多項式序列補成8項。
故（用(0,0)來補項）
多項式$f(x)$的系數序列為$A=\left\{(3,0),(2,0),(1,0),(0,0),(0,0),(0,0),(0,0),(0,0)\right\}$
多項式$g(x)$的系數序列為$B=\left\{(2,0),(1,0),(2,0),(0,0),(0,0),(0,0),(0,0),(0,0)\right\}$
對於每一個$f(w_8^i)$，可以尋得規律
$f(w_8^i)=A_0+A_1w_8^i+A_2(w_8^i)^2+A_3(w_8^i)^3+A_4(w_8^i)^4+A_5(w_8^i)^5+A_6(w_8^i)^6+A_7(w_8^i)^7$
$=A_0+A_2(w_8^i)^2+A_4(w_8^i)^4+A_6(w_8^i)^6+w_8^i(A_1+A_3(w_8^i)^2+A_5(w_8^i)^4+A_7(w_8^i)^6)$
$=A_0+A_2w_4^i+A_4(w_4^i)^2+A_6(w_4^i)^3+w_8^i(A_1+A_3w_4^i+A_5(w_4^i)^2+A_7(w_4^i)^3)$
$=A_0+A_4(w_4^i)^2+w_4^i(A_2+A_6(w_4^i)^2)+w_8^i(A_1+A_5(w_4^i)^2+w_4^i(A_3+A_7(w_4^i)^2))$
$=A_0+A_4w_2^i+w_4^i(A_2+A_6w_2^i)+w_8^i(A_1+A_5w_2^i+w_4^i(A_3+A_7w_2^i))$
然后
$f(w_8^{i+4})=f(-w_8^i)=A_0+A_1(-w_8^i)+A_2(-w_8^i)^2+A_3(-w_8^i)^3+A_4(-w_8^i)^4+A_5(-w_8^i)^5+A_6(-w_8^i)^6+A_7(-w_8^i)^7$
$=A_0+A_2(w_8^i)^2+A_4(w_8^i)^4+A_6(w_8^i)^6-w_8^i(A_1+A_3(w_8^i)^2+A_5(w_8^i)^4+A_7(w_8^i)^6)$
只要算出來$A_0+A_2w_4^i+A_4(w_4^i)^2+A_6(w_4^i)^3$和$A_1+A_3w_4^i+A_5(w_4^i)^2+A_7(w_4^i)^3$的值，就可以算出$f(w_8^i)$和$f(w_8^{i+4})$的值
同理，算出$A_0+A_4(w_4^i)^2$和$A_2+A_6w_2^i$的值，不僅可以得到$A_0+A_2w_4^i+A_4(w_4^i)^2+A_6(w_4^i)^3$的值，還可以得到$A_0+A_2w_4^{i+2}+A_4(w_4^{i+2})^2+A_6(w_4^{i+2})^3$
然后將每一個序列分成長度為2的小序列(下面是用1~8來編號的)
后面的是重點，對於每一個長度為2的小序列（設為(x,y)）:
第一個元素$x$重新賦值為$x+w_1^0*y$
第二個元素$y$重新賦值為$x-w_1^0*y$
(這個$w_1^1$這樣來的：$w_{當前序列一半長度l}^{一半長度l的第i個元素,i從0到l-1}$)
重新賦值完之后就可以合並了，每兩個長度為2的相鄰序列合並為一個長度為4的序列
對於每一個長度為4的序列（設為(a,b,c,d)）:
第一個元素$a$重新賦值為$a+w_2^0*c$
第三個元素$c$重新賦值為$a-w_2^0*c$
第二個元素$b$重新賦值為$b+w_2^1*d$
第四個元素$d$重新賦值為$b-w_2^1*d$
重新賦值完之后就可以合並了，每兩個長度為4的相鄰序列合並為一個長度為8的序列（此時全部小序列已經合並為一個序列了）
對於整個序列（長度為8），按照上面的賦值方法重新賦值（可以參考代碼），於是FFT變換就ok了
兩個序列就變成了
$A=\left\{A_0',A_1',A_2',A_3',A_4',A_5',A_6',A_7'\right\}$
$B=\left\{B_0',B_1',B_2',B_3',B_4',B_5',B_6',B_7'\right\}$
將兩個序列對應元素相乘，得到新序列
$C=\left\{A_0'*B_0',A_1'*B_1',A_2'*B_2',A_3'*B_3',A_4'*B_4',A_5'*B_5',A_6'*B_6',A_7'*B_7'\right\}$
然后將C序列進行FFT逆變換，就是多項式相乘后的系數序列了。
ps：FFT逆變換就是把當前序列再來一次FFT變換，只不過在重新賦值環節的時候是乘上$w_{mid}^{-1}$，而不是$w_{mid}^i$，最后得到的序列虛部為0，將實部全部除以8即可（這里不懂參考代碼或前面逆變換的講解）
FFT迭代版本
FFT的序列更新規律
綜上所述，FFT的序列更新是有規律的，先把序列分解為序列長度為2的小段，每一段更新完畢后就合並小段，繼續一組一組更新，再合並，最后合並成一個序列，然后再最后一次更新完畢。
如圖所示：
最后序列$\left\{A_0''',A_1''',...,A_{n-1}'''\right\}$為原序列進過FFT變換后的序列。
不像遞歸版本，迭代版本是直接從每一個小區間開始更新，然后合並，然后更新...，但是我們需要找到原系數序列經過奇偶分離的一系列操作后每個值的位置發生了什么變化。
假設原系數序列為$\left\{A_0,A_1,A_2,A_3,A_4,A_5,A_6,A_7\right\}$，下圖是奇偶分離的示意圖
沒錯，這不是巧合，序列經過奇偶分離之后，前后位置數的二進制是對稱的，根據這一點我們就不需要用遞歸來實現，可以直接來得到奇偶分離最終序列，然后再進行FFT
FFT代碼
FFT迭代版比遞歸版快的不是一點。。
關於序列奇偶分離可以用下面方法
int rev[(n+2)>>1]; for(i=0;i<n;i++){ rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1)); if (i<rev[i])swap(a[i],a[rev[i]]); }
於是就可以實現迭代版本FFT
void FFT(Complex *a,int n,int inv){ int bit=0,i,rev[(n+2)>>1]; while((1<<bit)<n) bit++; //交換序列 for(i=0;i<n;i++){ rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1)); if (i<rev[i])swap(a[i],a[rev[i]]);//if來保證只存在小的和大的交換，不會發生大的和小的交換 } for (int mid=1;mid<n;mid*=2){//mid枚舉區間長度的一半 Complex com(cos(pi/mid),inv*sin(pi/mid));//單位根 for (int i=0;i<n;i+=mid*2){//i來枚舉每個區間的開頭位置 Complex omega(1,0); for (int j=0;j<mid;j++,omega*=com){//j枚舉每個區間的前一半元素 Complex x=a[i+j],y=omega*a[i+j+mid]; a[i+j]=x+y,a[i+j+mid]=x-y;//蝴蝶效應 } } } }
總接代碼
FFT遞歸版本
void FFT(Complex *a,int n,int inv){ //a是你要進行變換的數組（多項式系數數組），n是數組大小，inv=1表示正變換，inv=-1表示逆變換 if(n==1) return;//如果分治后的序列大小為1，就直接返回，不用繼續分治了 int mid=n>>1;//如果n不等於1，那就繼續分治，mid是分治后變成兩個序列的長度 static Complex b[1000];//創建一個輔助用的b數組，后面體現用處 for(int i=0;i<mid;i++) b[i]=a[2*i],b[i+mid]=a[2*i+1];//將a數組奇數位置的值和偶數位置的值分開，b數組前mid個位置存a數組偶數位置值，后mid個位置存a數組奇數位置值 for(int i=0;i<n;i++) a[i]=b[i];//將b數組值重新賦給a數組 //對分治后的兩個序列（一個是a數組前mid個元素組成的序列，另一個數a數組后mid元素組成的序列） 進行FFT變換 FFT(a,mid,inv); FFT(a+mid,mid,inv); //算出一半，得到另一半的值，a數組前mid元素相當於之前說的f1，a數組后mid元素相當於之前說的f2 for(int i=0;i<mid;i++){ Complex com(cos(2*pi*i/n),inv*sin(2*pi*i/n));// b[i]=a[i]+com*a[i+mid];b[i+mid]=a[i]-com*a[i+mid]; // cout<<n<<" "<<i<<" "<<b[i].real()<<" "<<b[i].imag()<<endl; } for(int i=0;i<n;i++) a[i]=b[i];//重新將b數組賦值給a數組 return; }
FFT迭代版本
void FFT(Complex *a,int n,int inv){ //a是你要進行變換的數組（多項式系數數組），n是數組大小，inv=1表示正變換，inv=-1表示逆變換 int bit=0,i,rev[(n+2)>>1]; while((1<<bit)<n) bit++; //交換序列 for(i=0;i<n;i++){ rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1)); if (i<rev[i])swap(a[i],a[rev[i]]);//if來保證只存在小的和大的交換，不會發生大的和小的交換 } for (int mid=1;mid<n;mid*=2){//mid枚舉區間長度的一半 Complex com(cos(pi/mid),inv*sin(pi/mid));//單位根 for (int i=0;i<n;i+=mid*2){//i來枚舉每個區間的開頭位置 Complex omega(1,0); for (int j=0;j<mid;j++,omega*=com){//j枚舉每個區間的前一半元素 Complex x=a[i+j],y=omega*a[i+j+mid]; a[i+j]=x+y,a[i+j+mid]=x-y;//蝴蝶效應 } } } }