[整理]FFT（快速傅立葉變換）隨記
0.前言
本文對資訊學競賽中的 FFT 演算法及其推導作出一點粗淺的探討，希望能夠幫到大家一點點。
另：作者預設大家已經有了一些高中數學基礎，對於課本內基礎知識部分可能會涉及較少。建議初中生先學完有關三角函式、複數、向量的部分之後再來學習 FFT 。
1.前置知識
多項式
形如\(A(x)=\sum_{i=0}^n a_ix^i\)的式子稱作多項式（Polynomial），多項式可以進行加減乘（除法在此不涉及）等運算。
這種表示方法稱作多項式的係數表示（Coefficient representation），多項式還有一種點值表示（Point-value representation）法：
眾所周知，平面上的\(n\)
點值表示的優點：將多項式乘法的複雜度由係數表示的\(O(n^2)\)減少到了\(O(n)\)（只需要將點值相乘）。
缺點：係數表示與點值表示的互化還是需要\(O(n^2)\)，不過它是可以優化的（這也正是 FFT 做的事情）。
複數及其運算
首先我們有\(\sqrt{-1}=i\)，並將形如\(z=a+bi\ (a,b\in\mathbb{R})\)的數稱為複數，並定義其共軛\(\bar{z}=a-bi\)。
其中，\(\Re z=a\)稱為實部
以上是複數的代數表示，其運算同多項式運算，運算時要注意\(i^2=-1\)。（例如\((a+bi)(c+di)=(ac-bd)+(ad+bc)i\)）
複數還可以與平面上的向量一一對應，我們把這種表示叫做複數的幾何表示，該平面叫做複平面。
對於一個複數，它的幾何表示由模長（向量的長度）和幅角（從\(x\)軸正半軸逆時針轉到該向量的角度）組成。
複數幾何表示的加減法運算同向量運算，而乘法的規則是模長相乘、幅角相加。
尤拉公式：\(e^{i\theta}=\cos\theta+i\sin\theta\)。
單位根及其性質
接下來車速加快，請仔細閱讀！
定義方程\(z^n=1\)
如圖，這是方程\(z^3=1\)的三個複數根\(1,-\dfrac 12+\dfrac{\sqrt{3}}2i,-\dfrac 12-\dfrac{\sqrt{3}}2i\)：
如果一個單位根的\(0,\cdots,n-1\)次方恰好構成互不相同的所有單位根，我們就將其稱為本原單位根。
顯然\(\omega_n=e^{\frac{2\pi i}{n}}\)是一個\(n\)次本原單位根（為方便，下文中“本原單位根”即特指\(\omega_n\)）。
單位根有一些美妙的性質，我們來看一下。（假設以下的\(n\)都是正偶數）
\(\text{Theorem 1: }\omega_n^{2k}=\omega_{\frac{n}2}^k\)
\(\text{Proof: }\)利用複數幾何表示的相乘法則證明。
\(\text{Theorem 2: }\omega_n^{\frac{n}2+k}=-\omega_n^k\)
\(\text{Proof: }\)可以理解為多轉了半圈之後恰與原來相反。
習題一：寫出方程\(z^6=1\)的所有複數根。
2.Cooley-Tukey演算法
該演算法分為兩個過程： DFT（Discrete Fourier Transform） 和 IDFT（Inverse Discrete Fourier Transform） ，主要思想是分治。
鋪墊了這麼多，我們要想，如何利用單位根的一些性質來減少運算量呢？
兩位巨佬 Cooley 和 Tukey 想出了一種方法：將每一項按照指數奇偶分類。
我們要計算一個多項式的 DFT ，也就是計算其點值表示，需要將每一個單位根代入計算。
即\(A(\omega_n^k)=\sum_{j=0}^{n-1}a_i\omega_n^{kj}=\sum_{j=0}^{\frac{n}2-1}a_{2j}\omega_n^{2kj}+\omega_n^k\sum_{j=0}^{\frac{n}2-1}a_{2j+1}\omega_n^{2kj}\)。
根據剛剛提到的\(\text{Theorem 1}\)，我們可以將其改寫成\(\sum_{j=0}^{\frac{n}2-1}a_{2j}\omega_{\frac{n}2}^{kj}+\omega_n^k\sum_{j=0}^{\frac{n}2-1}a_{2j+1}\omega_{\frac{n}2}^{kj}\)。
此時再根據\(\text{Theorem 2}\)，我們可以寫出\(A(\omega_n^k)\)和\(A(\omega_n^{k+\frac{n}2})\)的表示：
我們成功了！現在只需要分別求奇數項和偶數項的 DFT 就能求出來兩個值了！運算量減少了一半！（可以證明這個過程的確是\(O(n\log n)\)的）
但是先別高興得太早，我們還有一步：將運算完的點值表示轉化回係數表示（畢竟大多數時候我們要的是乘積多項式的係數）。這時候應該怎麼辦呢？
Notice:由於證明較複雜且超出了高中的知識範圍，建議從此處開始跳過直接看結論。
注意到我們求值的過程實際上就是計算了一個矩陣乘法：
而我們的任務就是求係數矩陣（假設叫做\(V\)）的逆矩陣。
注意到這個矩陣是一個 Vandermonde 矩陣，對此，我們可以很容易地求出它的逆矩陣。
兩個比較繁複的證明：
利用行列式及伴隨矩陣
利用拉格朗日插值
然而最終的結果非常優美：設\(V=\{v_{ij}\}\)，則\(V^{-1}=\frac 1n\{v_{ij}^{-1}\}\)。
Notice:現在你可以在此停止了，結論就在下方。
所以，我們只需要將 DFT 中\(\omega_n^k\)換成\(\omega_n^{-k}\)（體現為虛部取相反數），再將答案乘以\(\frac 1n\)就行了！
由此，我們就可以用一個函式來實現 DFT 和 IDFT 的過程了。
3.遞迴實現 FFT
一層層向下遞迴，求出 FFT ：
void FFT(C *f,int len,int opt){//f是要進行DFT的陣列，len是長度，opt是要進行的操作DFT/IDFT
if(len==1)return;
C *f0=f,*f1=f+(len>>1);
for(rg int i=0;i<len;i++)tmp[i]=f[i];
for(rg int i=0;i<(len>>1);i++){//按奇偶性分成兩個多項式
f0[i]=tmp[i<<1],f1[i]=tmp[i<<1|1];
}
FFT(f0,len>>1,opt),FFT(f1,len>>1,opt);//分別進行FFT
C wn=C(cos(2*Pi/len),opt*sin(2*Pi/len)),w=C(1,0);//wn是本原單位根，w是wn的冪
for(rg int i=0;i<(len>>1);i++){
tmp[i]=f0[i]+w*f1[i];//按照上面推出來的式子進行計算
tmp[i+(len>>1)]=f0[i]-w*f1[i];
w=w*wn;
}
for(rg int i=0;i<len;i++)f[i]=tmp[i];
}
C++ STL 為我們提供了複數類，但是據說常數巨大，而手寫複數類也不難寫，所以還是建議大家手寫：
struct C {
double re,im;
C(double _re=0.0,double _im=0.0){
re=_re,im=_im;
}
}f[N],g[N],tmp[N];
il C operator + (C a,C b){
return C(a.re+b.re,a.im+b.im);
}
il C operator - (C a,C b){
return C(a.re-b.re,a.im-b.im);
}
il C operator * (C a,C b){
return C(a.re*b.re-a.im*b.im,a.re*b.im+a.im*b.re);
}
而最終實現可以這麼寫：
FFT(f,len,1),FFT(g,len,1);
...//進行計算
FFT(f,len,-1);
4.迭代實現 FFT
注意到上面的實現方式需要進行大量遞迴，常數巨大（據說在洛谷會只有77分，但我實測能很寬鬆地通過），我們考慮如何優化常數。
這時候就需要拿出這張著名的圖片了（注意下面最後兩個框順序反了）：
這是我們進行遞迴的順序，大家觀察一下二進位制有什麼發現？
可以發現，我們其實就是按照反轉二進位制排序之後依次分組的！如\(0000,1000,0100,1100\)倒過來就是\(0000,0001,0010,0011\)。
所以，我們可以將要進行 FFT 的陣列按照反轉二進位制排序，然後直接迭代即可！
for(rg int i=1;i<len;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));//求出i的反轉二進位制（其中l是二進位制位數）
}
而我們不需要遞迴的 FFT 函式變成了這樣：
void FFT(C *f,int len,int opt){
for(rg int i=0;i<len;i++){
if(i<r[i])swap(f[i],f[r[i]]);//排序
}
for(rg int i=2;i<=len;i<<=1){
int mid=i>>1;
C wn=C(cos(2*Pi/i),opt*sin(2*Pi/i));
for(rg int j=0;j<len;j+=i){
C w=C(1,0);
for(rg int k=0;k<mid;k++){
C z=f[j+k+mid]*w;
f[j+k+mid]=f[j+k]-z,f[j+k]=f[j+k]+z;
w=w*wn;
}
}
}
}
用時和記憶體都有了優化。
另外，關於 FFT 還有一個優化：三次變兩次。實現方法是將多項式\(g\)放到多項式\(f\)的虛部，這時\(f^2\)的虛部除以二就是我們要的答案。這個方法也可以省去不少常數。
5.應用
FFT 的最大應用是優化多項式卷積的計算。所謂卷積其實就是多項式乘法，它的定義是：對於兩個函式\(f(n),g(n)\)，它們的卷積\((f*g)(n)=\sum_{i+j=n}f(i)g(j)\)。（就是個聽起來挺nb的名字）
所以，當我們碰到一個噁心的式子時，可以嘗試將其化為卷積的形式。
例題：洛谷P3338 [ZJOI2014]力
題目讓我們求\(E_j=\sum_{i=1}^{j-1}\dfrac{q_i}{(i-j)^2}-\sum_{i=j+1}^{n}\dfrac{q_i}{(i-j)^2}\)的值，我們想辦法將它化一化。
為了好看，我們設\(f_i=q_i,g_i=\dfrac{1}{i^2}\)。這時候原式變成了\(\sum_{i=1}^{j-1}f_ig_{j-i}-\sum_{i=j+1}^nf_ig_{i-j}\)。
我們優化一下上下界，令\(f_0=g_0=0\)，有\(E_j=\sum_{i=0}^jf_ig_{j-i}-\sum_{i=j}^nf_ig_{i-j}\)。
發現第一項已經是卷積形式，我們來著重推一推第二項\(\sum_{i=j}^nf_ig_{i-j}\)。
運用和式中變數代換的技巧，這個式子就等於\(\sum_{i=0}^{n-j}f_{i+j}g_i\)。
這裡還有一個常用的技巧，就是引入反轉函式\(f'_i=f_{n-i}\)。此時原式變成了\(\sum_{i=0}^{n-j}f'_{n-j-i}g_i\)，再代換一次就是\(\sum_{i=0}^mf'_{m-i}g_i\)。
我們把它化成了卷積形式，現在可以用 FFT 做了！
核心程式碼（程式碼裡h[i]就是我們說的\(f'_i\)）：
int main(){
Read(n);
for(rg int i=1;i<=n;i++){
cin>>f[i].re;h[n-i].re=f[i].re;
g[i].re=1.0/(1.0*i*i);
}
while(len<n+n+1)len<<=1,l++;
for(rg int i=1;i<len;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
}
FFT(f,len,1),FFT(g,len,1),FFT(h,len,1);
for(rg int i=0;i<=len;i++)f[i]=f[i]*g[i],h[i]=h[i]*g[i];
FFT(f,len,-1),FFT(h,len,-1);
for(rg int i=1;i<=n;i++){
printf("%.3lf\n",(f[i].re-h[n-i].re)/len);
}
return 0;
}
6.總結
FFT 是多項式的基礎演算法，也是許多新手學多項式要過的第一關。理解了此處的知識，才能為之後的學習打下良好基礎。
習題二：洛谷P3803 【模板】多項式乘法（FFT）
習題三：洛谷P1919 【模板】A*B Problem升級版（FFT快速傅立葉）