感謝 路人黑的紙巾， 理論部分來源於地址
FFT原理：將多項式的系數表示轉換為點值表示，從而進行卷積運算，理論上從\(O(n^2)\)降低到\(O(nlogn)\)。
\(f(x)\)為\(n\)次多項式，\(g(x)\)為\(m\)次多項式，\(f(x)\cdot g(x)\)為\(m+n\)次多項式。
將多項式從系數表示法轉為點值表示法更加方便兩個多項式相乘（只要將他們的對應項的做乘積即可），同時為了方便對遞歸算法進行改寫，將點值擴展到它們的和向上最近的那個2的冪次方個：
這里\(r\)大於\(m+n\)最小的2的冪次方。然后選擇\(r\)個插值點，使得每個點的函數計算的時間花費最小。
如上圖當\(r=8\)時可以在一維復空間上這樣等分的選擇這8個插值點（每個點是單位圓上的等分點），是因為這些點有以下的幾個性質可以用來減少計算量
性質
-
每個插值點：\(w_r^k=cos\frac{k}{r}2\pi + isin\frac{k}{r}2\pi,(k=0,...,r-1)\)
-
\(w_{2r}^{2k}=w_r^k\)
證：\(w_{2r}^{2k}=cos\frac{2k}{2r}2\pi+isin\frac{2k}{2r}2\pi=cos\frac{k}{r}2\pi+isin\frac{k}{r}2\pi=w_r^k\)
-
\(w_r^{k+\frac{r}{2}}=-w_r^{k}\)
證：\(w_r^{k+\frac{r}{2}}=cos\frac{k+\frac{r}{2}}{r}2\pi+isin(\frac{k+\frac{r}{2}}{r}2\pi)=cos(\frac{k}{r}2\pi+\pi)+isin(\frac{k}{r}2\pi+\pi)=-w_r^k\)
-
\(w_n^0=w_n^n\)
-
\(w_n^iw_n^j=w_n^{i+j}\)（復數乘法運算規律）
過程
首先對多項式進行整理有：
這里有\(t_1(x) = \sum_{i=0}^{i=\frac{r}{2}-1}a_{2i}x^i\)， \(t_2(x) = \sum_{i=0}^{i=\frac{r}{2}-1}a_{2i+1}x^i\)。
設\((k<\frac{r}{2})\)，\(w^k_r\)帶入有：
將\(w_r^{k+\frac{n}{2}}\)帶入有：
可以發現兩個問題：
- 計算\(w_r^k\)時只需要在計算前一半的時候將符號進行改變，就能得到關於原點對稱的另外一個值，就不需要單獨的在計算一次計算。
- 計算\(w_r^k\)的過程實際上是將該問題分解成為兩個不相交的子問題，每個子問題可以繼續向下遞歸求解。
以\(w_8^0\)的計算過程為例：
- 每個公示的上面的數字表示的是當前多項式的系數
- 如果當求解的數值變成了\(w_8^k\)，只需要對每個黑框中的數值進行相應的替換
如果考慮到不同的\(w_8^k\)的求解，那么上面的求解過程實際上是在下面的表中遍歷了一棵二叉樹
- 第一個\(+0\)表示的是左邊的分支（\(t_1(w_4^0)\)） ，加上（也就是”+“的意思）\(w_8^0\)（“0”代表的是右上標）乘上右邊的分支（\(t_2(w_4^0)\)）。
- 由於兩個關於原點對稱的單位復數的符號相反，所以后一部分第二項前面乘的系數可以\(+w_8^{4+k}\)寫成\(-w_8^k\)（簡寫成為“\(-k\)"）
- 中間的黑色豎杠將區間在不同深度上分成了不相交的子區間
下面是幾種情況的求解過程，可以理解下
可以看出所有的多項式的求解都是在這個\(4\times8\)的表格上進行的，從下到上維護整個表格便可以得到傅里葉變換，每一層維護的時間復雜度是\(O(n)\)，從\(w_2\)到\(w_8\)一共三層（每層都是類似開方的降系數），所以層數是\(O(logn)\)，因此總的時間復雜度是\(O(nlogn)\)
傅里葉變化的優化實際上利用\(w_{2r}^{2k}=w_r^k\) 這個性質，巧妙的將計算的復雜度降了下來。如果遞歸求解，實際上也是要遞歸到最后一層，將該層的信息整體維護好供上層使用。
逆傅里葉變換（IFFT）（最近有點忙，有時間再看）未完待續。。。。
模板：
代碼的一些細節，比如rev[i]的作用和如何實現的，可以參考別人的文章。
typedef complex<double> cpx;
const double pi = acos(-1.0);
void fft_init(cpx *coe1, cpx *coe2, int n, int &bit, int &bit_len) {
bit = 0;
while ((1<<bit) < n) bit++;
bit_len = 1<<bit; // 找到最小的2的冪次方數
for (int i = 0; i < bit_len; i++) coe1[i] = coe2[i] = 0;
for (int i = 0; i < bit_len; i++) //初始化每個元素對應的最后的位置
rev[i] = (rev[i>>1]>>1) | ((i&1)<<(bit-1));
}
void fft(cpx *coe, int bitn, int bit, int inv) {
for (int i = 0; i < bitn; i++)
if (i < rev[i]) swap(coe[i], coe[rev[i]]);
// mid枚舉隔板的長度
for (int mid = 1; mid < bitn; mid <<= 1) {
cpx tmp(cos(pi/double(mid)), inv*sin(pi/mid));
// 枚舉當前在哪一個隔板那里
for (int i = 0; i < bitn; i += (mid<<1)) {
cpx omega(1.0, 0.0);
for (int j = 0; j < mid; j++, omega *= tmp) {
// 計算兩邊函數的值，隔板的左邊是x+y，隔板右邊是x-y
cpx x = coe[i+j], y = omega*coe[i+j+mid];
coe[i+j] = x + y, coe[i+j+mid] = x - y;
}
}
}
// 如果是逆映射的話要除掉整體的個數
if (inv < 0) for (int i = 0; i < bitn; i++) coe[i] /= double(bitn);
}