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

---

零、前言

參考：

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

洛穀日報 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}{ll} \\ 1 & \operatorname{FFT}(\boldsymbol{a}, limit) \\ 2 & \qquad \textbf{if } limit = 1 \\ 3 & \qquad \qquad \textbf{return } \\ 4 & \qquad \boldsymbol{a}^{[0]} = (a_0,a_2,a_4,\cdots,a_{n-2}) \\ 5 & \qquad \boldsymbol{a}^{[1]} = (a_1,a_3,a_5,\cdots,a_{n-1}) \\ 6 & \qquad \operatorname{FFT}(\boldsymbol{a}^{[0]}, \frac{limit}{2}) \\ 7 & \qquad \operatorname{FFT}(\boldsymbol{a}^{[1]}, \frac{limit}{2}) \\ 8 & \qquad \omega_n = e^{\frac{2\pi}{n}i} = \cos(\frac{2\pi}{n})+i\sin(\frac{2\pi}{n}) \\ 9 & \qquad \omega = 1 \\ 10 & \qquad \textbf{for } k = 0 \ldots \frac{n}{2} - 1 \\ 11 & \qquad \qquad \boldsymbol{a}_k = \boldsymbol{a}_k^{[0]} + \omega \boldsymbol{a}_k^{[1]} \\ 12 & \qquad \qquad \boldsymbol{a}_{k+\frac{n}{2}} = \boldsymbol{a}_k^{[0]} - \omega \boldsymbol{a}_k^{[1]} \\ 13 & \qquad \qquad \omega = \omega \omega_n \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。