快速傅里葉變換簡介

F a s t F o u r i e r T r a n s f o r m Fast \ Fourier \ Transform Fast Fourier Transform ，簡稱 F F T FFT FFT，這東西似乎在數字通信領域有大用處：（EDA表示很感動）

計算量小的顯著的優點，使得FFT在信號處理技術領域獲得了廣泛應用，結合高速硬體就能實作對信號的實時處理，例如，對語音信號的分析和合成，對通信系統中實作全數字化的時分制與頻分制(TDM/FDM)的復用轉換，在頻域對信號濾波以及相關分析，通過對雷達、聲納、振動信號的頻譜分析以提高對目標的搜索和跟蹤的解析度等等，都要用到FFT，可以說FFT的出現，對數字信號處理學科的發展起了重要的作用，——百度百科（快速傅里葉變換詞條）

但是本博客只介紹其在演算法競賽中的應用——快速計算兩個多項式的乘積，

假設現在有兩個多項式： A ( x ) = x 2 + 3 x + 2 A(x)=x^2+3x+2 A(x)=x2+3x+2， B ( x ) = 2 x 2 + 1 B(x)=2x^2+1 B(x)=2x2+1，要求 C ( x ) = A ( x ) × B ( x ) C(x)=A(x) \times B(x) C(x)=A(x)×B(x)，樸素做法是按照乘法分配律，將其中一個多項式的每一項分別與另一個多項式的每一項相乘，如果這兩個多項式最高次數為 N N N，這個演算法的時間復雜度將會是 O ( n 2 ) O(n^2) O(n2)，然而我們可以用 F F T FFT FFT在 O ( n l o g n ) O(nlogn) O(nlogn)的時間內計算出 C ( x ) C(x) C(x)， F F T FFT FFT演算法的核心，是四個十分精妙的想法，發明這個演算法的人應該是天選者，

下面將介紹這四個巧妙的想法，以及 F F T FFT FFT的代碼和實作細節，

1. 多項式的另一種表示

一般怎么用計算機存盤多項式？

開一個陣列 A [ i ] A[i] A[i]，表示多項式 x i x^i xi 項的系數，例如 F ( x ) = 4 x 4 + 6 x 3 + 5 x 2 + 3 x + 2 F(x)=4x^4+6x^3+5x^2+3x+2 F(x)=4x4+6x3+5x2+3x+2，那么 A [ ] = { 4 ， 6 ， 5 ， 3 ， 2 } A[]=\{4，6，5，3，2\} A[]={4，6，5，3，2}，這完全沒有問題，形如 F ( x ) = ∑ a i x i F(x)=\sum a_ix^i F(x)=∑ai?xi 的多項式表示方法，叫做“系數表示法”，

然而現在為了更快計算多項式乘積，我們使用多項式的另一種表示法，要理解這種表示方法，可以從最簡單的多項式開始，設 A ( x ) = a 0 + a 1 x A(x)=a_0+a_1x A(x)=a0?+a1?x，一次函式就是一種簡單的多項式，它的最高冪次為1，

我們知道 A ( x ) = a 0 + a 1 x A(x)=a_0+a_1x A(x)=a0?+a1?x 不但是一個多項式，還是一條直線，由于兩點能確定一條直線，我們在平面中隨便取兩點 ( x 1 ， y 1 ) (x_1，y_1) (x1?，y1?)和 ( x 2 ， y 2 ) (x_2，y_2) (x2?，y2?)，那么系數 a 0 a_0 a0?， a 1 a_1 a1? 也就隨之確定了，而我們現在在講的“多項式的另一種表示方法”，就是用平面上的點來確定多項式，

高次多項式的影像是一條曲線，顯然，無論你在平面上取多少點，也無法唯一確定一條曲線，但卻可以確定一個多項式，因為多項式“并不是普通的（或者說任意的）曲線”，這里有一個結論：平面上的 N + 1 N+1 N+1 個點，能唯一確定一個 N N N 次多項式，例如：你在平面上隨機取四個點，那么有且只有一個三次函式同時穿過這四個點，這個結論的正確性需要用線性代數的知識來證明，詳細寫可能需要一篇新的博客，所以我只能簡單說明一下正確性，其實是不會

設 F ( x ) = ∑ i = 0 i = n a i x i F(x)=\sum_{i=0}^{i=n}a_ix^i F(x)=∑i=0i=n?ai?xi，我們選取的點為 { x 0 ， x 1 ， x 2 ， . . . ， x n } \{x_0，x_1，x_2，...，x_n\} {x0?，x1?，x2?，...，xn?}，那么：
[ F ( x 0 ) F ( x 1 ) ? F ( x n ) ] = [ 1 x 0 x 0 2 ? x 0 n 1 x 1 x 1 2 ? x 1 n ? ? ? ? ? 1 x n x n 2 ? x n n ] [ a 0 a 1 ? a n ] \left[ \begin{matrix} F(x_0) \\ F(x_1) \\ \vdots \\F(x_n) \end{matrix} \right]= \left[ \begin{matrix} 1 & x_0 &x^2_0 & \cdots & x_0^n \\ 1 & x_1 &x^2_1 & \cdots & x_1^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n &x^2_n & \cdots & x_n^n \\ \end{matrix} \right] \left[ \begin{matrix} a_0 \\ a_1 \\ \vdots \\ a_n \end{matrix} \right] ??????F(x0?)F(x1?)?F(xn?)???????=??????11?1?x0?x1??xn??x02?x12??xn2???????x0n?x1n??xnn??????????????a0?a1??an????????
易證中間矩陣的行列式不為0，所以方程組只有一組解，好了證畢，

這種表示多項式的方法叫做“點值表示法”，假設對于兩個多項式 f ( x ) f(x) f(x)和 g ( x ) g(x) g(x)，我們選取相同的 x x x 序列 { x 0 , x 1 . . . . , x n } \{ x_0, x_1....,x_n\} {x0?,x1?....,xn?}，那么兩個多項式可以表示為：
f ( x ) = ( ( x 0 , f ( x 0 ) ) , ( x 1 , f ( x 1 ) ) , . . . , ( x n , f ( x n ) ) ) g ( x ) = ( ( x 0 , g ( x 0 ) ) , ( x 1 , g ( x 1 ) ) , . . . , ( x n , g ( x n ) ) ) f(x)=(\ (x_0,f(x_0)),(x_1,f(x_1)),...,(x_n,f(x_n))\ ) \\ g(x)=(\ (x_0,g(x_0)),(x_1,g(x_1)),...,(x_n,g(x_n))\ ) f(x)=( (x0?,f(x0?)),(x1?,f(x1?)),...,(xn?,f(xn?)) )g(x)=( (x0?,g(x0?)),(x1?,g(x1?)),...,(xn?,g(xn?)) )
設 F ( x ) = f ( x ) × g ( x ) F(x)=f(x) \times g(x) F(x)=f(x)×g(x)，則 F ( x ) = { ( x 0 , f ( x 0 ) g ( x 0 ) ) , ( x 1 , f ( x 1 ) g ( x 1 ) ) , . . . , ( x n , f ( x n ) g ( x n ) ) } F(x)=\{(x_0,f(x_0)g(x_0)),(x_1,f(x_1)g(x_1)),...,(x_n,f(x_n)g(x_n))\} F(x)={(x0?,f(x0?)g(x0?)),(x1?,f(x1?)g(x1?)),...,(xn?,f(xn?)g(xn?))}

容易發現，如果已知兩個多項式的點值表示，求兩個多項式乘積多項式的點值表示，可以在線性時間內處理完，如果能找到一種方法，能快速將一個多項式在系數表示和點值表示之間轉化，就能實作快速多項式乘法，

2.系數到點值的轉換（DFT）

既然任意 N + 1 N+1 N+1 個點就可以確定一個 N N N 階多項式，那么如果我直接暴力隨機選 N + 1 N+1 N+1 個橫坐標，帶入多項式求值會怎樣呢？容易想到這樣做的時間復雜度是 O ( n 2 ) O(n^2) O(n2)（找 N + 1 N+1 N+1個點 O ( N ) O(N) O(N)，每個點求值 O ( N ) O(N) O(N)），這樣就又回到了原點，沒有達到加速的目的，所以需要找另外的方法，

依舊考慮最簡單的情況：現在，我想在二次函式 f ( x ) = x 2 f(x)=x^2 f(x)=x2上，任意取 8 個點，并求出這8個點的函式值，有沒有什么比較快的方法？容易想到，其實我只需要任意取四個 x x x 就行了，因為如果我知道了 x = x 0 x=x_0 x=x0? 的函式值，我立刻就知道了 x = ? x 0 x=-x_0 x=?x0? 的函式值，因為 f ( x ) = x 2 f(x)=x^2 f(x)=x2 是一個偶函式，奇函式同理，只不過我需要在函式值上再加一個負號，總之：我們可以利用奇偶性來減少選取點的數目，

下面的問題就是，如何將一個一般的多項式，轉換成可以利用奇偶性進行優化的多項式，我們可以將原多項式，按照每一項 x x x 的冪次分組，偶數次冪一組，奇數次冪一組，然后再把奇數次冪組中的每一項提出一個 x x x 來，就像下面這樣：
F ( x ) = 3 x 5 + 2 x 4 + x 3 + 7 x 2 + 5 x + 1 ? F ( x ) = ( 2 x 4 + 7 x 2 + 1 ) + x ( 3 x 4 + x 2 + 5 ) F(x)=3x^5+2x^4+x^3+7x^2+5x+1 \Rightarrow \\ F(x)=(2x^4+7x^2+1)+x(3x^4+x^2+5) F(x)=3x5+2x4+x3+7x2+5x+1?F(x)=(2x4+7x2+1)+x(3x4+x2+5)
我們可以把左右兩部分視為新的，關于 x 2 x^2 x2 的多項式，
F 1 ( x 2 ) = 2 x 4 + 7 x 2 + 1 F 2 ( x 2 ) = 3 x 4 + x 2 + 5 F_1(x^2) = 2x^4+7x^2+1 \\ F_2(x^2) = 3x^4+x^2+5 F1?(x2)=2x4+7x2+1F2?(x2)=3x4+x2+5
我們只需要解決這兩個小多項式的求值，再利用 F ( x ) = F 1 ( x 2 ) + x F 2 ( x 2 ) F(x)=F_1(x^2)+xF_2(x^2) F(x)=F1?(x2)+xF2?(x2) 和 F ( ? x ) = F 1 ( x 2 ) ? x F 2 ( x 2 ) F(-x)=F_1(x^2)-xF_2(x^2) F(?x)=F1?(x2)?xF2?(x2) 這兩個公式代回即可，最終取的點的形式是： { F ( x 1 ) , F ( ? x 1 ) , F ( x 2 ) , . . . . , F ( x n / 2 ) , F ( ? x n / 2 ) } \{F(x_1), F(-x_1), F(x_2),....,F(x_{n/2}),F(-x_{n/2})\} {F(x1?),F(?x1?),F(x2?),....,F(xn/2?),F(?xn/2?)}，而兩個小多項式都只有原多項式一半的規模，所以這個方法的復雜度是 O ( n log ? n ) O(n\log{n}) O(nlogn)，但是這個方法還存在一個問題：我們想利用奇偶性和每次選取 x = x 0 x=x_0 x=x0? 和 x = ? x 0 x=-x_0 x=?x0? 兩個點，來減少選取點的次數，也就是說我們每次取的是成相反數的兩個點，然而當問題變為關于 F 1 ( x 2 ) F_1(x^2) F1?(x2) 和 f 2 ( x 2 ) f_2(x^2) f2?(x2) 時，可以發現： x 2 x^2 x2 永遠為正，我們不再能取一對相反數了，而解決這個問題的方法，被人們視為 F F T FFT FFT的精髓，

3.單位復根與快速傅里葉變換

可以看出，問題出在多項式的定義域上，一旦問題的規模縮小，原多項式的定義域會失效，那么如果我們把多項式的定義域擴大到復數域呢？這次不從最簡單的情況入手了，這次需要一個逆向思維，下面結合圖來說明，

假設遞回已經進行到了終點，也就是現在多項式只剩下一次項了，也就是說我們需要選取一個點來求值，那么這個點如何選取？顯然我可以選 x = 1 x=1 x=1，然后考慮上一層遞回，我要選兩個值，我希望他們倆是相反數，還希望他們倆的平方等于 1，顯然我可以取 x = 1 x=1 x=1 和 x = ? 1 x=-1 x=?1，再往上一層，左邊是 1，和第一層一樣，關鍵是右邊的 -1，我同樣希望上一層的兩個數，他們互為相反數，而他們的平方等于 -1，這兩個數顯然是 i i i 和 ? i -i ?i，

到目前為止，我們得到了四個數：1，-1，i，-i ，如果我的原多項式是一個三階多項式，我只需要取 x x x 等于這四個數字就可以，但如果我的多項式階數更高，這四個數還不夠，我需要至少 N N N 個才行，這 N N N 個數需要滿足 x N = 1 x^N = 1 xN=1，而之前學過的代數知識告訴我們， x N = 1 x^N=1 xN=1 在復數域內恰好有 N N N 個解，也就是說我們前面得到的四個數字，其實是 x 4 = 1 x^4=1 x4=1 的四個解，

x n = 1 x^n=1 xn=1 在復數域上的 n n n 個解，均勻地分布在復平面的單位圓上，例如： x 8 = 1 x^8=1 x8=1 的八個解，就是單位圓的八個八等分點，讀者可以自行驗證，這樣的點叫做單位復根，將(1，0)點編號為0，其余點依次編號，我們定義 ω n k \omega_{n}^k ωnk? 是 x n = 1 x^n=1 xn=1 的第 k k k 號單位復根，
（本來想做個單位圓的圖，然而沒有比較好的幾何工具，就咕咕了）

簡單推導可得： ω n k = cos ? ( 2 π k n ) + i × sin ? ( 2 π k n ) \omega_{n}^k=\cos(\frac{2\pi k}{n})+i \times \sin(\frac{2\pi k}{n}) ωnk?=cos(n2πk?)+i×sin(n2πk?)

單位復根有三個重要性質，通過簡單推導就可以得到：

• ω n n = 1 \omega_{n}^n=1 ωnn?=1
• ω n k = ω 2 n 2 k \omega_{n}^k=\omega_{2n}^{2k} ωnk?=ω2n2k?
• ω 2 n k + n = ? ω 2 n k \omega_{2n}^{k+n}=-\omega_{2n}^k ω2nk+n?=?ω2nk?

那么原來的公式 F ( x ) = F 1 ( x 2 ) + x F 2 ( x 2 ) F(x)=F_1(x^2)+xF_2(x^2) F(x)=F1?(x2)+xF2?(x2)，利用單位復根的性質可得：

同理可得：

因此我們求出了 D F T ( G ( ω n / 2 k ) ) DFT(G(\omega^k_{n/2})) DFT(G(ωn/2k?)) 和 D F T ( H ( ω n / 2 k ) ) DFT(H(\omega^k_{n/2})) DFT(H(ωn/2k?)) 后，就可以同時求出 D F T ( f ( ω n k ) ) DFT(f(\omega^k_n)) DFT(f(ωnk?)) 和 D F T ( f ( ω n k + n / 2 ) ) DFT(f(\omega^{k+n/2}_n)) DFT(f(ωnk+n/2?)) ，于是對 G G G 和 H H H 分別遞回 DFT 即可，考慮到分治 DFT 能處理的多項式長度只能是二的整數次冪，否則在分治的時候左右不一樣長，右邊就取不到系數了，所以要在第一次 DFT 之前就把序列向上補成長度為 2 m 2^m 2m（高次系數補 0 0 0）、最高項次數為 2 m ? 1 2^m-1 2m?1 的多項式，

所以到現在我們可以發現，最初的原多項式代入的 n n n 個點，實際就是 { ω n 0 , ω n 1 , ω n 2 , . . . , ω n n ? 1 } \{\omega_n^0,\omega_n^1,\omega_n^2,...,\omega_n^{n-1}\} {ωn0?,ωn1?,ωn2?,...,ωnn?1?}，其中 n n n 是二的整數次冪，

得到兩個多項式乘積后，還需要把這個結果多項式從點值表示變回系數表示，

4.點值到系數的轉換（IDFT）

這里我真寫不下去了，因為 DFT 和 IDFT 怎么結合到一起的我真啥也沒看懂，就知道一個結論：做 DFT 和 IDFT 的區別就在于一個正負號，
1 ω k = ω k ? 1 = e ? 2 π i k = cos ? ( 2 π k ) + i × sin ? ( ? 2 π k ) \frac{1}{\omega_k}=\omega_k^{-1}=e^{-\frac{2\pi i}{k}}=\cos(\frac{2\pi}{k})+i \times \sin(-\frac{2\pi}{k}) ωk?1?=ωk?1?=e?k2πi?=cos(k2π?)+i×sin(?k2π?)
因此實際在做 F F T FFT FFT時，通常會傳入一個引數，來決定是做 D F T DFT DFT 還是 I D F T IDFT IDFT，

之后要是弄懂了可能會回來寫，

板子

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

P1919 【模板】A*B Problem升級版.