【学习笔记】快速离散傅里叶变换（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。