【學習筆記】快速離散傅里葉變換（FFT）（遞歸版）

本文講述的是快速離散傅里葉變換的遞歸版，並非倍增版。

---

零、前言

參考：

具體學習並實現快速傅里葉變換 - 鶴翔萬里

洛谷日報 71：傅里葉變換 (FFT) 學習筆記 - command_block

在這里特別感謝。

代碼中的 ll 是 long long，有在代碼之前加上 typedef long long ll;。

1. 概念

快速傅里葉變換（Fast Fourier Transform)，全稱快速離散傅里葉變換（Fast Discrete Fourier Transform），即利用計算機計算離散傅里葉變換（DFT）的高效、快速計算方法的統稱，簡稱 FFT。快速傅里葉變換是 1965 年由 J.W.庫利 和 T.W.圖基 提出的。采用這種算法能使計算機計算離散傅里葉變換所需要的乘法次數大為減少，特別是被變換的抽樣點數 \(N\) 越多，FFT 算法計算量的節省就越顯著。

以上改自 百度百科 - 快速傅里葉變換。

快速傅里葉變換不是傅里葉變換的優化，而是離散傅里葉變換的。

2. 前置知識

為了學習快速傅里葉變換，我們先來學習以下前置知識：

1. 復數
2. 單位根
3. 一元多項式

現在開始講一下前置知識。

一、復數

1. 概念

復數（復雜的數，complex number），是數的概念擴展。

我們把形如 \(z=a+bi\)（\(a\)、\(b\)均為實數）的數稱為復數。其中，\(a\) 稱為實部，\(b\) 稱為虛部，\(i\) 稱為虛數單位。當 \(z\) 的虛部 \(b=0\) 時，則 \(z\) 為實數；當 \(z\) 的虛部 \(b\neq0\) 時，實部 \(a=0\) 時，常稱 \(z\) 為純虛數。復數域是實數域的代數閉包，即任何復系數多項式在復數域中總有根。

以上摘自 百度百科 - 復數。

虛數單位 \(i=\sqrt{-1}\)。

2. 四則運算

復數加法

\[(a+bi)+(c+di)=a+bi+c+di=(a+c)+(bi+di)=(a+c)+(b+d)i \]

復數減法

\[(a+bi)-(c+di)=a+bi-c-di=(a-c)+(bi-di)=(a-c)+(b-d)i \]

復數乘法

\[(a+bi)(c+di)=ac+adi+bci+bi\cdot di=ac+adi+bci-bd=(ac-bd)+(ad+bc)i \]

復數除法

\[\dfrac{a+bi}{c+di}=\dfrac{(a+bi)(c-di)}{(c+di)(c-di)}=\dfrac{(ac+bd)+(-ad+bc)i}{(c^2+d^2)+(-cd+cd)i}=\dfrac{ac+bd}{c^2+d^2}+\dfrac{bc-ad}{c^2+d^2}i \]

3. 歐拉定理

復變函數中，\(e^{ix}=\cos{x}+i\sin{x}\) 稱為歐拉公式，\(e\) 是自然對數的底，\(i\) 是虛數單位。

當 \(x=\pi\) 時

\[e^{i\pi}=\cos{\pi}+i\sin{\pi}=-1+0\\e^{i\pi}+1=0 \]

即歐拉恆等式。

4. 三角形式

把 \(x\) 換一下，變成 \(\theta\)，方便一點。

明顯地，\(e^{i\theta}=\cos{\theta}+i\sin{\theta}\)，把一個復數寫成了 \(e\) 的指數的形式，於是一個復數也可以寫成 \(z=re^{i\theta}\)。

像這樣的表示形式，被稱為三角形式。

其中 \(r\) 是 \(z\) 的模，即 \(r = |z|\)；\(θ\) 是 \(z\) 的輻角，在區間 \([-\pi,\pi]\) 內的輻角稱為輻角主值，記作 \(\arg(z)\)。

三角形式下的復數乘法

根據乘方的運算法則，我們可以得到：

\[ae^{i\alpha}\cdot be^{i\beta}=abe^{i(\alpha+\beta)} \]

5. 形式互換

\(a+bi=r(\cos{\theta}+i\sin{\theta})=re^{i\theta}\)，其中 \(r=\sqrt{a^2+b^2}\)，\(\cos{\theta}=\dfrac{a}{r}\)，\(\sin{\theta}=\dfrac{b}{r}\)

6. 幾何意義

復數的幾何意義

由於有實部和虛部，我們可以在一個平面直角坐標系中表示一個復數，\(y\) 軸作為虛軸，\(x\) 作為實軸。而這個平面，叫做復數平面，簡稱復平面。

若用三角形式表示，復數所對應的向量長度稱為復數的模，該向量與實軸正方向的夾角為復數的輻角。

於是我們可以更好地理解形式互換公式了。

復數加法的幾何意義

設 \(A=1+i,B=2+i\)，

\(C=A+B=3+2i,D=A+C=4+3i,E=B+D=6+4i\)。

我們把 \(A,B,C,D,E\) 在復平面上表示出來。

發現規律了嗎？

同向量加法，我們把每一個復數都看成向量，然后做向量加法。

對於每一個復數，加上另一個復數，相當於以第一個加數為原點重新畫復平面，然后在這個新的復平面上畫第二個加數，然后回來，即為和。

上述例子是虛部實部都為正數的情況，對於其他的情況也是一樣的，大家可以自行探究。

三角形式下的復數乘法的幾何意義

設 \(A=1+i,B=2+i\)，

\(C=AB=1+3i,D=AC=-2+4i\)。

我們把 \(A,B,C,D\) 在復平面上表示出來。

我們把復數換成三角形式，然后就會發現，復數乘積的模，即為乘數的模之積，復數乘積的輻角，即為乘數的輻角之和。

誒，你問我為什么不給出上述復數的三角形式？

你可以自己算，根據形式互換公式（其實是因為我懶，這里空間太小了，我寫不下）。

7. 代碼實現

C++ 有自帶的復數頭文件 complex，屬於 STL，可以使用但是慢，所以還是建議手寫復數。

const double Pi = acos(-1.0);
struct Complex_Number {
    double real, imag;
    Complex_Number() {}
    Complex_Number(double num1, double num2) : real(num1) , imag(num2)  {}
};
Complex_Number operator + (const Complex_Number &num1, const Complex_Number &num2) {
    return Complex_Number(num1.real + num2.real, num1.imag + num2.imag);
}
Complex_Number operator - (const Complex_Number &num1, const Complex_Number &num2) {
    return Complex_Number(num1.real - num2.real, num1.imag - num2.imag);
}
Complex_Number operator * (const Complex_Number &num1, const Complex_Number &num2) {
    return Complex_Number(num1.real * num2.real - num1.imag * num2.imag, num1.real * num2.imag + num1.imag * num2.real);
}

二、單位根

1. 概念

數學上，\(n\) 次單位根（unit root）是 \(n\) 次冪為 \(1\) 的復數。它們位於復平面的單位圓上，構成正 \(n\) 邊形的頂點，其中一個頂點是 \(1\)。

以上內容改編自 百度百科 - 單位根。

簡單來說單位根就是方程 \(z^n=1 \ (n=1,2,3,\cdots)\) 在復數范圍內的 \(n\) 個根。

這方程的復數根 \(z\) 為 \(n\) 次單位根。

\[z^n=r^ne^{n\theta i}=1 \]

因為單位根的 \(n\) 次方為 \(1\)，所以單位根的模一定為 \(1\)。

因為單位根的 \(n\) 次方為 \(1\)，所以單位根的輻角的 \(n\) 倍都是 \(2\pi\) 的倍數。

\[\left\{\begin{matrix} r = 1 \\ n\theta = 2\pi k \end{matrix}\right. \]

單位的 \(n\) 次根有 \(n\) 個，我們可以找到這些點：

\[\omega_n^k=e^{i\frac{2\pi k}{n}}=\cos{\frac{2\pi k}{n}}+i\sin{\frac{2\pi k}{n}} \ \ \ (k = 0,1,2,\cdots,n-1) \]

每一個單位根都均勻地落在單位圓上，如圖（\(8\) 次單位根）：

同時每一個單位根都可以看作 \(\omega_n = e^{i\frac{2\pi}{n}}\) 的冪，所以這個單位根也被稱作主 \(n\) 次單位根，記作 \(\omega_n\)。

單位根有它的性質，這里有 \(3\) 個重要的性質

2. 消去引理

基本

\[\omega_{dn}^{dk}=\omega_n^k \]

證明：

\[\omega_{dn}^{dk}=(e^{i\frac{2\pi}{dn}})^{dk}=(e^{i\frac{2\pi}{n}})^{k}=\omega_n^k \]

幾何意義

\(z\) 是 \(4\) 次單位根，\(a\) 是 \(8\) 次單位根。

應該可以看出來吧...

3. 折半引理

基本

\[(\omega_n^{k+\frac{n}{2}})^2=(\omega_n^k)^2=\omega_\frac{n}{2}^k \]

證明:

\[\omega_n^{k+\frac{n}{2}}=\omega_n^k\omega_n^\frac{n}{2}=-\omega_n^k,\ (-\omega_n^k)^2=\ (\omega_n^k)^2=\omega_n^{2k}=\omega_{\frac{n}{2}}^k \]

幾何意義

還是拿這張圖：

\(z\) 是 \(4\) 次單位根，\(a\) 是 \(8\) 次單位根。

\((\omega_n^{k+\frac{n}{2}})^2\) 相當於，先繞半圈，再繞半圈，但是會多一位。

\((\omega_n^k)^2\) 相當於 \(\omega_n^k\) 的下一個單位根。

然后使用消去引理。

4. 求和引理

基本

\[\sum_{i = 0}^{n - 1}{(\omega_n^k)^i} = 0 \]

證明：

\[\sum_{i = 0}^{n - 1}{(\omega_n^k)^i} = \dfrac{(\omega_n^k)^n - 1}{\omega_n^k - 1} = \dfrac{(\omega_n^n)^k - 1}{\omega_n^k - 1} = \dfrac{1^k - 1}{\omega_n^k - 1} = 0 \]

三、一元多項式

1. 概念

由數和字母的積組成的代數式叫做單項式，單獨的一個數或一個字母也叫做單項式（例：\(0\) 可看做 \(0\) 乘 \(a\)，\(1\) 可以看做 \(1\) 乘指數為 \(0\) 的字母，\(b\) 可以看做 \(b\) 乘 \(1\)），分數和字母的積的形式也是單項式。

在數學中，由若干個單項式相加組成的代數式叫做多項式（若有減法：減一個數等於加上它的相反數）。多項式中的每個單項式叫做多項式的項，這些單項式中的最高項次數，就是這個多項式的次數。其中多項式中不含字母的項叫做常數項。

以上摘自 百度百科 - 單項式 和 百度百科 - 多項式。

• 單項式次數：所有變數字母的指數之和

• 多項式次數：多項式 \(F\) 的次數為其最高項的次數，記作 \(\operatorname{degree}(F)\)

• 多項式次數界：多項式 \(F\) 的次數界為任意一個嚴格大於\(F\) 的次數的整數

對於 FFT，我們只研究一元多項式，即只有一個未知數的多項式。

為了方便，對於多項式的項的系數數列編號從 \(0\) 開始。

2. 系數表示

如果用一個 \(n\) 項數列 \(a\) 表示多項式 \(A\) 的每一項數，\(x\) 表示字母部分，則：

\[A(x) = a_0+a_1x+a_2x^2+a_3x^3+\cdots+a_{n-1}x^{n-1} \]

或者：

\[A(x) = \sum_{i = 0}^{n-1}{a_ix^i} \]

這就是多項式的系數表示法。

明顯地，計算一個系數表示的多項式的值為 \(O(n)\)。

一元多項式在系數表示下的加法

整式加法，對應次數項加起來即可。

設多項式 \(A(x) = \sum_{i = 0}^{n-1}{a_ix^i}\) 和多項式 \(B(x) = \sum_{i = 0}^{n-1}{b_ix^i}\)。

則：

\[C(x)=A(x)+B(x)=\sum_{i = 0}^{n-1}{a_ix^i}+\sum_{i = 0}^{n-1}{b_ix^i}=\sum_{i = 0}^{n-1}{(a_i + b_i)x^i} \]

若我們運用計算機做整式加法，明顯地，時間復雜度為 \(O(n)\)。

一元多項式在系數表示下的乘法

整式乘法，用一個多項式的每一項去乘另一個多項式的每一項，將積相加。

設多項式 \(A(x) = \sum_{i = 0}^{n - 1}{a_ix^i}\) 和多項式 \(B(x) = \sum_{i = 0}^{n-1}{b_ix^i}\)。

則：

\[C(x)=A(x)\cdot B(x)=\sum_{i = 0}^{n - 1}{a_ix^i}\cdot \sum_{i = 0}^{n - 1}{b_ix^i}=\sum_{i = 0}^{\operatorname{degree}(C)}{c_ix^i} \]

其中 \(c_i = \sum_{j=0}^{i}{a_jb_{i-j}}\)，\(\operatorname{degree}(C)=\operatorname{degree}(A)+\operatorname{degree}(B)=n - 1+n- 1=2n - 2=2(n-1)\)。

\(c\) 為 \(a,b\) 的卷積，記作 \(c=a\otimes b\)。

若我們運用計算機做整式乘法，明顯地，時間復雜度為 \(O(n^2)\)，時間復雜度過高，不優秀。

3. 點值表示

平時我們的多項式表示方法一般都是系數表示，因為較為直觀，為多項式的概念。

即為單項式之和的形式。

多項式還有一種表示方法，叫做點值表示。

即用至少 \(n\) 個多項式上的點表示。

\[\{(x_0,A(x_0)),(x_1,A(x_1)),(x_2,A(x_2),\cdots,(x_n,A(x_n))\} \]

求這個多項式的值，被稱為插值，可以使用拉格朗日插值公式進行求解。

由於我自己也不會，所以直接給出公式。

\[A(x)=\sum_{i=0}^{n}{A(x_i)\dfrac{{\textstyle \prod_{j\neq i}{(x-x_j)}}}{{\textstyle \prod_{j\neq i}{(x_i-x_j)}}}} \]

時間復雜度為 \(O(n^2)\)，不夠優秀。

一元多項式在點值表示下的加法

對於每個對應的 \(x\) 的點，直接相加即可

\[C(x_i)=A(x_i)+B(x_i) \]

明顯地，時間復雜度為 \(O(n)\)。

一元多項式在點值表示下的乘法

對於每個對應的 \(x\) 的點，直接相加即可

\[C(x_i)=A(x_i)B(x_i) \]

明顯地，時間復雜度為 \(O(n)\)，效率遠遠高於在系數表示下的乘法。

四、快速傅里葉變換

我們看到，一元多項式在點值表示下的乘法是非常快的，所以我們想把系數表示的多項式轉化成點值表示的多項式，這樣就可以快速求出兩個多項式的卷積了。但是我們轉化的過程是很慢的，而快速傅里葉變換（FFT）就是使用單位根優化這個過程的算法。

1. 離散傅里葉變換

離散傅里葉變換（Discrete Fourier Transform，DFT），是這樣一個東西：

\[A(x) = \sum_{i = 0}^{n - 1}{a_ix^i} \ \text{在} \ \omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1} \ \text{處的值} \ y_0,y_1,\cdots,y_{n-1} \]

\[y_i = A(\omega_n^i) = \sum_{j = 0}^{n - 1}{\omega_n^{ij} a_j} \]

記作：\(\boldsymbol{y}=\text{DFT}_n(\boldsymbol{a})\) 或 \(\boldsymbol {y}=\mathcal{F} \boldsymbol{a}\)

看的可能有點迷惑。

簡單來說，就是傅里葉非常無聊，把單位根帶到一個多項式里面去了...

然后這個東西可以在其他的地方發揮作用。

2. 快速傅里葉變換

離散傅里葉變換並沒有降低時間復雜度，霍納法則（秦九韶算法）求值時間復雜度為 \(O(n^2)\)。

但是我們發現離散傅里葉變換帶入了單位根，我們可以使用單位根的性質進行優化。

首先我們有一個 \(n-1\) 次多項式，\(A(x)\)，有系數向量 \(\boldsymbol{a} = [a_0,a_1,a_2,\cdots,a_{n-1}]^T\)。

注意 FFT 的 \(n\) 為 \(2\) 的冪。

我們現在把這個向量拆開，分為偶數項和奇數項。

記為 \(\boldsymbol{a}^{[0]}\) 和 \(\boldsymbol{a}^{[1]}\)。這兩個向量所對應的兩個多項式我們分別記為 \(A^{[0]}(x)\) 和 \(A^{[1]}(x)\)。

\[\boldsymbol{a}^{[0]} = [a_0,a_2,a_4,\cdots,a_{n-2}]^T \rightarrow A^{[0]}(x) \\ \boldsymbol{a}^{[1]} = [a_1,a_3,a_5,\cdots,a_{n-1}]^T \rightarrow A^{[1]}(x) \]

現在把三個多項式都帶入 \(x\)。

\[A(x) = a_0+a_1x+a_2x^2+a_3x^3+\cdots+a_{n-1}x^{n-1} \\ A^{[0]}(x) = a_0+a_2x+a_4x^2+a_6x^3+\cdots+a_{n-2}x^{\frac{n}{2}-1} \\ A^{[1]}(x) = a_1+a_3x+a_5x^2+a_7x^3+\cdots+a_{n-1}x^{\frac{n}{2}-1} \]

把后兩個多項式自變量換成 \(x^2\)。

\[A(x) = a_0+a_1x+a_2x^2+a_3x^3+\cdots+a_{n-1}x^{n-1} \\ A^{[0]}(x^2) = a_0+a_2x^2+a_4x^4+a_6x^6+\cdots+a_{n-2}x^{n-2} \\ A^{[1]}(x^2) = a_1+a_3x^2+a_5x^4+a_7x^6+\cdots+a_{n-1}x^{n-2} \]

把 \(A^{[1]}\) 乘上 \(x\)。

\[A(x) = a_0+a_1x+a_2x^2+a_3x^3+\cdots+a_{n-1}x^{n-1} \\ A^{[0]}(x^2) = a_0+a_2x^2+a_4x^4+a_6x^6+\cdots+a_{n-2}x^{n-2} \\ xA^{[1]}(x^2) = a_1x+a_3x^3+a_5x^5+a_7x^7+\cdots+a_{n-1}x^{n-1} \]

然后就會發現：

\[A(x) = A^{[0]}(x^2) + xA^{[1]}(x^2) \]

我們現在就將原問題轉化成了求次數界為 \(\dfrac{n}{2}\) 的兩個多項式在每個單位根的平方上的值然后相加的值。

比較拗口但是還是很好理解的。

我們把兩個具體的單位根帶入看看。

\[A(\omega_n^k) = A^{[0]}((\omega_n^k)^2) + \omega_n^kA^{[1]}((\omega_n^k)^2) \\ A(\omega_n^{k+\frac{n}{2}}) = A^{[0]}((\omega_n^{k+\frac{n}{2}})^2) + \omega_n^{k+\frac{n}{2}}A^{[1]}((\omega_n^{k+\frac{n}{2}})^2) \]

再根據消去引理：\(\omega_{dn}^{dk}=\omega_n^k\)，

和折半引理：\((\omega_n^{k+\frac{n}{2}})^2=(\omega_n^k)^2=\omega_\frac{n}{2}^k\)，

以及折半引理證明中的 \(\omega_n^{k+\frac{n}{2}} = -\omega_n^k\)。

帶入化簡得到：

\[A(\omega_n^k) = A^{[0]}(\omega_\frac{n}{2}^k) + \omega_n^kA^{[1]}(\omega_\frac{n}{2}^k) \\ A(\omega_n^{k+\frac{n}{2}}) = A^{[0]}(\omega_\frac{n}{2}^k) - \omega_n^kA^{[1]}(\omega_\frac{n}{2}^k) \]

發現了嗎？除了符號不同之外都是一樣的。

所以我們只要知道了前一半，后一半的值我們也能知道了。

這里寫一下偽代碼。

\[\begin{array}{l} \operatorname{Fast-Fourier-Transform}(\boldsymbol{a}, \textit{limit}) \\ \begin{array}{ll} 1 & \textbf{if } \textit{limit} = 1 \\ 2 & \qquad \textbf{return } \\ 3 & \boldsymbol{a}^{0} \gets \{ a_0,a_2,a_4,\cdots,a_{n-2} \} \\ 4 & \boldsymbol{a}^{1} \gets \{ a_1,a_3,a_5,\cdots,a_{n-1} \} \\ 5 & \operatorname{Fast-Fourier-Transform}(\boldsymbol{a}^{0}, \frac{\textit{limit}}{2}) \\ 6 & \operatorname{Fast-Fourier-Transform}(\boldsymbol{a}^{1}, \frac{\textit{limit}}{2}) \\ 7 & \omega_n \gets e^{\frac{2\pi}{n}i} \\ 8 & \omega \gets 1 \\ 9 & \textbf{for } k \text{ from } 0 \text{ to } \frac{n}{2} - 1 \\ 10 & \qquad \boldsymbol{a}_k \gets \boldsymbol{a}_k^{0} + \omega \boldsymbol{a}_k^{1} \\ 11 & \qquad \boldsymbol{a}_{k+\frac{n}{2}} \gets \boldsymbol{a}_k^{0} - \omega \boldsymbol{a}_k^{1} \\ 12 & \qquad \omega \gets \omega \omega_n \end{array} \end{array} \]

3. 代碼實現

例：洛谷 - P3803 【模板】多項式乘法（FFT）。

void FFT(Complex_Number * p, ll len, bool IDFT) {
    if(len == 1) return;
    for(ll i = 0; i < len; i++) poly[i] = p[i];
    Complex_Number * lp = p, * rp = p + (len >> 1);
    for(ll i = 0; i < (len >> 1); i++) {
        lp[i] = poly[i << 1]; rp[i] = poly[i << 1 | 1];
    }
    FFT(lp, len >> 1, IDFT);
    FFT(rp, len >> 1, IDFT);
    Complex_Number unitroot(cos(2.0 * Pi / len), sin(2.0 * Pi / len));
    if(IDFT) unitroot.imag *= -1;
    Complex_Number w(1.0, 0.0), t;
    for(ll i = 0; i < (len >> 1); i++) {
        t = w * rp[i];
        poly[i] = lp[i] + t;
        poly[i + (len >> 1)] = lp[i] - t;
        w = w * unitroot;
    }
    for(ll i = 0; i < len; i++) p[i] = poly[i];
}
int main(){
    scanf("%lld%lld",&n,&m);
    for(ll i = 0; i <= n; i++) scanf("%lf", &poly1[i].real);
    for(ll i = 0; i <= m; i++) scanf("%lf", &poly2[i].real);
    for(m += n, n = 1; n <= m; n <<= 1);
    FFT(poly1, n, false); FFT(poly2, n, false);
    for(ll i = 0; i < n; i++) poly1[i] = poly1[i] * poly2[i];
    FFT(poly1, n, true);
    for(ll i = 0; i <= m; i++) printf("%lld ", ll(poly1[i].real / n + 0.45));
    return 0;
}

五、快速傅里葉變換的應用——大數乘法

1. 基本思路

把乘數分解為多項式即 \(a_1+a_2x+a_3x^2+a_4x^3+\cdots\)，其中 \(x = 10\) 的形式。

\[A = a_1+a_2x+a_3x^2+a_4x^3+\cdots \\ B = b_1+b_2x+b_3x^2+b_4x^3+\cdots \]

其中 \(x = 10\)。

然后使用 FFT 相乘，最后全加起來求值即可。

例：洛谷 - P1919 【模板】A*B Problem升級版（FFT快速傅里葉），洛谷 - P1303 A*B Problem。