快速傅里叶变换（FFT）
liserver
·
·
算法·理论
快速傅里叶变换
前言
关于此算法，应用广泛。但是在 OI 算法竞赛中，我们基本只关注它“加速多项式乘法” 这一用途。
本算法在《NOI 大纲（2023 年修订版）》中难度评级为 10，请酌情观看。
本文适用于未接触过此算法的初学者。
对于本文用词不当、概念错误等问题，请发布在讨论区，我看到会及时修改，力争全文的每一句话都可以被引用而无误。
规范
算法的名称为“快速傅里叶变换（Fast Fourier Transform）”。在一些文章中，它被翻译为“快速傅立叶变换”。其原因为作者 Fourier 的音译问题。本文章按照“傅里叶” 翻译。所有此名称的衍生名词的翻译也遵守此规定。
用途
快速傅里叶变换可以用来加速多项式乘法，将 O(n^2) 的暴力算法优化为 O(n \log n)。
并且，我们可以将形如 1234 这样的十进制数字看作 1x^3+2x^2+3x^1+4x^0 这个多项式在 x=10 情况下的值，这也是高精度的原理。两个数字相乘即是像这样的两个多项式相乘在 x=10 情况下的值（进位另算），所以快速傅里叶变换还可以加速高精度乘法。
FFT 还有更多的用途，如信号处理、图像处理、量子计算等领域均有广泛应用，但是本文并不深入研究。
基础知识
本算法大量依赖数学，大部分为高中数学内容，掺杂初中知识，本文一并讲解。本文未进行解释的，即默认你阅读本文章所要掌握的知识，你可以在其它途径查询。
多项式
既然快速傅里叶变换是用来加速多项式乘法，我们首先要理解明白什么是多项式。
下文所讨论的多项式默认指只包含一个字母 x 的多项式，符号可能与一些数学教材有些许出入。
多项式的定义
多项式的定义即是几个单项式的和。因为次数相同的单项式可以合并同类项，所以多项式的最简形式应是一些次数各不相同的单项式加法形式。
\begin{aligned}
F(x)&=\sum_{i=0}^na_ix^i\\
&=a_nx^n+a_{n-1}x^{n-1}+...+a_0x^0
\end{aligned}
此多项式 F 是多项式的一般形式。若单项式中最高次数为 n，那么多项式的次数定义为 n。该多项式的次数界定义为所有 m 满足 m\gt n。
多项式的表示
既然我们是 OIer ，那我们必然会联想到一个问题：如何存储多项式？
观察多项式，发现对于 n 次多项式，每个单项式所不同的只有每个单项式的系数，所以我们可以开辟一个数组，存储多项式。数组的容量，就是可以存储的多项式的最高次数。
这种方法叫做系数表示法。
如果我们用这种方法，那么多项式 A 与 B 相乘就可以写作下列形式（伪代码）：
for i from 0 to n
for j from 0 to n
plus A[i] and B[j] to C[i + j]
不难发现时间复杂度 O(n^2)，在 n 较大时，就会非常慢。
那么能否有其它的方法能存储多项式，并且乘法乘起来还很快呢？
当然会有，但是在介绍这种方法之前，我们需要一条引理：
对于一个 n 次多项式，可以用 n+1 个点来确定。
思考此引理的含义，我们先从特殊情况出发，如果多项式是一次的，那么多项式的图像当然是一条直线。
而两点确定一条直线，我们可以用两个点确定这一条直线，此引理成立。
可以证明，对于普遍情况，此引理也成立。
有了这个非常棒的性质，我们就可以通过两个数组 X 和 Y 存储一个多项式 F。
如果点的 $x$ 默认确定，那么可以只用一个数组存储。
这种方法叫做点值表示法。
在此情况下，两个选点相同的多项式相乘可以写作：
```
for i from 0 to n
plus A[i] and B[i] to C[i]
```
不难发现 $O(n)$ 的时间复杂度。
那么这样是否就解决了多项式乘法速度问题？
然而，我们忽略了一个问题：在转换和计算点值时，选点是 $O(n)$，计算是 $O(n)$，最后累计 $O(n^2)$。
是的，我们的思路是对的，它确实可以将 $O(n^2)$ 的计算转化为 $O(n)$ 的计算。但是题目一般都不会是直接提供点值表示法，我们需要转换，而**转换才是多项式乘法计算中的大头**。
转换，正是快速傅里叶变换的用武之地，它可以将多项式从系数表示法转化为点值表示法，并且时间复杂度只有 $O(n \log n)$。
并且，快速傅里叶逆变换还可以将点值表示法转化为系数表示法，时间复杂度同样只有 $O(n \log n)$。
实现这样的转换，我们则需要借助复数。
### 复数
复数是高中数学的内容，此文简略介绍，只保证能够看懂本算法，更多内容可以参见 OI Wiki（链接参见文末）。
#### 复数的引入
在实数域内，根号下的数字只能是正数，不能是负数。
原因很简单，在实数域内，平方满足如下性质：
$$
a^2\ge0
$$
当且仅当 $a = 0$ 时，$a^2=0$。
而在我们高中学习时，会知道有一个数字 $i$，有性质 $i^2=-1$，称为**虚数单位**。
虚数没有现实意义，它也不能被任何实数以及它们之间的运算所表示，但是它却有一个不错的性质。
形如 $bi$ 的形式的数字叫做**纯虚数**。
一切数字都可以写作 $a+bi$ 的形式，其中 $a$ 是实数，$bi$ 是纯虚数，这种形式的数字称作**复数**。
其中，$a$ 称为**实部**，$b$ 称为**虚部**。
如此，我们就可以得到目前数学上最大的数域：**复数域**。
复数也可以作为方程的解，并且延伸出一个重要的定理 —— **代数基本定理**：简单概括，$n$ 次方程有且仅有 $n$ 个根。
#### 复数的运算
复数运算的推导很简单，利用分配律把原式化为 $a+bi$ 的形式即可。
有复数 $z_1=a_1+b_1i$ 以及 $z_2=a_2+b_2i$，下文讨论两者的四则运算。
##### 复数加法
$$
\begin{aligned}
z_1+z_2&=(a_1+b_1i)+(a_2+b_2i)\\
&=a_1+b_1i+a_2+b_2i\\
&=(a_1+a_2)+(b_1+b_2)i
\end{aligned}
$$
##### 复数减法
$$
\begin{aligned}
z_1-z_2&=(a_1+b_1i)-(a_2+b_2i)\\
&=a_1+b_1i-a_2-b_2i\\
&=(a_1-a_2)+(b_1-b_2)i
\end{aligned}
$$
##### 复数乘法
$$
\begin{aligned}
z_1\times z_2&=(a_1+b_1i)\times(a_2+b_2i)\\
&=a_1(a_2+b_2i)+b_1i(a_2+b_2i)\\
&=a_1a_2+a_1b_2i+a_2b_1i-b_1b_2\\
&=(a_1a_2-b_1b_2)+(a_1b_2+a_2b_1)i
\end{aligned}
$$
##### 复数除法
$$
\begin{aligned}
z_1 \div z_2&=\frac{a_1+b_1i}{a_2+b_2i}\\
&=\frac{(a_1+b_1i)(a_2-b_2i)}{(a_2+b_2i)(a_2-b_2i)}\\
&=\frac{(a_1a_2+b_1b_2)+(a_2b_1-a_1b_2)i}{a_2^2+b_2^2}\\
&=\frac{a_1a_2+b_1b_2}{a_2^2+b_2^2}+\frac{a_2b_1-a_1b_2}{a_2^2+b_2^2}i
\end{aligned}
$$
#### 复数的几何意义
同实数一样，复数也有几何意义。实数与数轴上的点一一对应，复数则与平面上的一点一一对应。
但是，坐标轴需要改造一下：
复数的基本形式是 $a+bi$，我们按照此规定了复数的坐标轴。
横轴叫做实轴，代表复数的实部 $a$。纵轴是虚轴，代表复数的虚部 $b$。
而此坐标轴形成的平面，叫做**复平面**。
复数在复平面下的加减乘除也有特殊的意义，但是涉及到向量，超出了本文章的讨论范围。
#### 共轭复数
复数 $z=a+bi$ 的共轭复数定义为 $a-bi$ 记作 $\bar z$。
#### 欧拉公式
数学家欧拉发现了一个等式，就是著名的欧拉公式：
$$
e^{i\pi}+1=0
$$
这个公式，有自然常数 $e$ 和圆周率 $\pi$，有乘法的单位元 $1$，有加法的单位元 $0$，堪称数学最美公式。
但其实此式子只是欧拉公式的一个特殊情况，真正的欧拉公式长这样：
$$
e^{i\theta}=\cos(\theta)+i\sin(\theta)
$$
此公式将指数函数延伸到复数域，并与三角函数建立起了桥梁。
#### 单位根
##### 单位根的定义
考虑一个这样的方程：
$$
\omega^n=1
$$
根据当文提到的代数基本定理，此方程有 $n$ 个根。
此方程的 $n$ 个根，叫做 $n$ 次单位根，记作 $w_n$。
而这个方程还有特殊的性质，我们可以通过复平面来窥其究竟。
下图是当 $n=8$ 时的方程 8 个根在复平面上的分布。
是否发现，这 $n$ 个解在以原点为圆心，以 1 为半径的一个圆上均匀分布？
这个圆叫做单位圆。
可以证明，$\omega_n$ 是单位圆的 $n$ 等分点。
根据欧拉公式，我们还可以直接求出单位根：
$$
\begin{aligned}
w_n&=e^{\frac{2\pi i}{n}}\\
&=\cos(\frac{2\pi}{n})+i\sin(\frac{2\pi}{n})
\end{aligned}
$$
##### 单位根的性质
$n$ 次单位根有 $n$ 个，都是单位圆的 $n$ 等分点，此性质已经在上文提到。
$n$ 个单位根分别是欧拉公式求出的 $\omega_n$ 的 $m$ 次方，$m$ 满足 $0 \le m \le n$。
例如， 4 次单位根分别是 $\omega_4^1, \omega_4^2, \omega_4^3, \omega_4^4$。
同时，$\omega_n^0=\omega_n^n=1$。
现在提出一条引理（消去引理）：
$$
w_{dn}^{dk}=w_n^k
$$
关于此引理的证明，可以直接带入公式：
$$
w_{dn}^{dk}=e^{\frac{2\pi dki}{dn}}=e^{\frac{2\pi ki}{n}}=\omega_n^k
$$
根据消去引理，我们可以得到（$n$ 为偶数）：
$$
w_n^{\frac{n}{2}}=w_{2}=-1
$$
此推论的证明，也很简单：
$$
w_n^{\frac{n}{2}}=\omega_{2n}^n=\omega_2=e^{\frac{2\pi i}{2}}=e^{i\pi}=-1
$$
根据此结论，还可以进一步得到：
$$
w_n^{m+\frac{n}{2}}=\omega_n^m\times w_n^\frac{n}{2}=-\omega_n^m
$$
第二条引理（折半定理）：
如果 $n\gt0$ 为偶数，那么 $n$ 个 $n$ 次单位根的平方集合就是 $\frac{n}{2}$ 个 $\frac{n}{2}$ 次单位复数根集合。
证明：
根据消去引理：
$$
(w_n^m)^2=w_n^{2m}=w_{\frac{n}{2}}^m
$$
此时集合的数量仍然是 $n$，而引理说明应该是 $\frac{n}{2}$。
$$
(\omega_n^{m+\frac{n}{2}})^2=(-\omega_n^m)^2=(\omega_n^m)^2
$$
上述式子证明了 $(w_n^{m+\frac{n}{2}})^2$ 与 $(w_n^m)^2$ 相等，也就是集合少了一半。
## 多项式转换的优化
### 特殊情况下的优化
如前文所说，将系数表示法转化为点值表示法正是我们迫切需要解决的。
我们只考虑一个简单的多项式：
$$
\begin{aligned}
F(x)&=4x^2+0x^1+0x^0\\
&=4x^2
\end{aligned}
$$
根据前文的讨论，我们需要带入至少三个值。
现在思考一个问题：
我们能否少计算几个点？
换句话说，我们能否在取到一些点后，通过 $O(1)$ 算出另一些点的值？
答案是存在的，我们只需要把取的点分组，每一组的两个点互为相反数，那么一组中一个点得到了值，另一个点的值也就知道了。
原因就是，这个多项式是一个偶函数，即
$$
F(-x)=(-4x)^2=4x=F(x)
$$
在这个例子中，我们只是简单的优化了选点的方式，就使得计算量减半。
那么，再看一个例子：
$$
F(x)=7x^3
$$
这个例子是否适用？
当然适用，选点同样是相反数，但是一组中得到的两个值不再相等，而是互为相反数。
原因也很简单，$F$ 是一个奇函数：
$$
F(-x)=(-7x)^3=-(7x)^3=-F(x)
$$
### 由特殊到一般
现在，在两个例子中，我们都利用选点的艺术减少了时间复杂度，现在我们将它从特殊推向一般。
即，对于一个 $n$ 次多项式 $F$，我们考虑其如何优化。
$$
F(x) = \sum_{i=0}^{n}a_ix^i
$$
那么在这个情况下我们不能以普通的奇函数、偶函数来优化，因为这个多项式可能非奇非偶。
那我们只好凑出奇函数和偶函数。
将原多项式拆解为：
$$
\begin{aligned}
F_1(x)&=a_1x^0+a_3x^1+a_5x^2+...+a_{n-1}x^{\frac{n}{2}-1}\\
F_2(x)&=a_0x^0+a_2x^1+a_4x^2+...+a_nx^{\frac{n}{2}}
\end{aligned}
$$
我们按照原多项式的系数奇偶拆除了两个新的多项式。
注意，我们默认 $F$ 的次数 $n$ 是偶数，如果不是偶数，可以为 $F$ 添加一个次数为 $n+1$、系数为 0 的单项式，这样 $F$ 就变成了次数为偶数的多项式。
由于新多项式与原多项式的次数不同，所以不能直接暴力相加。
但是很容易得到：
$$
F(x)=F_1(x^2)+xF_2(x^2)
$$
这样，我们就将求 $F$ 的点值变成了求 $F_1$ 与 $F_2$ 这两个函数的点值。
而这又变成了两个新的求值问题，因为多项式是偶函数，按照第一个例子，我们要求的值变成了一半。
对于剩下需要求的值，我们仍然可以通过上述策略拆分函数。
如果一切顺利，那么时间复杂度就应该是 $O(n\log n)$。
但注意到 $x^2 \ge 0$。 所以，等到 $F_1$ 和 $F_2$ 求值时，所有的 $x$ 都变成了正数，与互为相反数的要求不符，递归不成立。
## 快速傅里叶变换
### 思想
如果你读到现在对于前文复数的知识有些遗忘，请立即返回温习，下文将会高频使用前文的引理、推理。
梳理一下，我们的点子是对的。但是选点上仍然有问题。
我们所需要的，是平方后的数字仍然以相反数成对出现。
既然平方后还能出现负数，则一定需要复数的参与才能完成。
恰好，**单位根满足这样的性质**。
所以，我们求的 $n$ 个点，就恰好可以选择 $n$ 次单位根，正恰好 $n$ 次单位根有 $n$ 个。
这正是快速傅里叶变换的精髓所在。
那么我们就将上文中的式子替换为 $n$ 次单位根。
$$
\begin{aligned}
F(\omega_n^m)&=F_1((\omega_n^m)^2)+\omega_n^mF_2((\omega_n^m)^2)\\
&=F_1(\omega_n^{2m})+\omega_n^mF_2(\omega_n^{2m})\\
&=F_1(\omega_{\frac{n}{2}}^m)+\omega_n^mF_2(\omega_{\frac{n}{2}}^m)
\end{aligned}
$$
更加惊奇的是，我们将 $\omega_n^{m+\frac{n}{2}}$ 带入：
$$
\begin{aligned}
F(\omega_n^{m+\frac{n}{2}})&=F_1((\omega_n^{m+\frac{n}{2}})^2)+w_n^{m+\frac{n}{2}}F_2((\omega_n^{m+\frac{n}{2}})^2)\\
&=F_1((-w_n^m)^2)-\omega_n^mF_2((-\omega_n^m)^2)\\
&=F_1(\omega_n^{2m})-\omega_n^mF_2(\omega_n^{2m})\\
&=F_1(\omega_{\frac{n}{2}}^m)-\omega_n^mF_2(\omega_\frac{n}{2}^m)
\end{aligned}
$$
发现带入 $\omega_n^m$ 与 $\omega_n^{m+\frac{n}{2}}$ 后两者只差了一个单项式系数的符号，所以我们在求出左半部分时，可以 $O(1)$ 求出右半部分。
这再次缩减了求值的次数。
并且，根据折半定理，求 $(\omega_n^m)^2$ 的集合实际上就是求所有 $\omega_{\frac{n}{2}}^m$ 的集合，所以 $F_1$ 与 $F_2$ 这两个子问题就变成了求 $\frac{n}{2}$ 个点值的问题，与父问题大抵相同，只不过是求值的个数少了一半，满足递归的条件。
这样，我们就可以递归地求解 $F_1$ 与 $F_2$。
不过这次我们没有用上偶函数的思想，而是在拆解的时候做了文章，同样缩小了一半的范围。
不难发现这是分治的思想，时间复杂度 $O(n\log n)$。
### 实现
关于这个算法的实现，一共有两种方法，一种是按照上面所讲述的递归实现，另一种是继续研究，将“拆分” 改为“合并” 的非递归实现。
#### 递归实现
递归实现没有什么新的理论知识，完全按照上文的实现，请直接阅读伪代码：
```
func FFT(a)
n = a.length // a 必须是 2 的 n 次幂
if n == 1
return a // 只有一个常量 代入值必定是本身
wn = complex(cos(2 * pi / n), sin(2 * pi / n))
w = complex(1, 0)
a1 = {a[1], a[3], a[5], ..., a[n - 1]} // 根据奇偶性对系数分组
a2 = {a[0], a[2], a[4], ..., a[n]}
y1 = FFT(a1) // 递归地调用、分解
y2 = FFT(y2)
y = {} // 返回值
for m from 0 to n / 2 - 1
y[m] = y1[m] + w * y2[m]
y[m + n / 2] = y1[m + n / 2] - w * y2[m + n / 2]
w = w * wn // 得到下一个单位根
return y
```
#### 非递归实现
观察图片中系数在每轮递归后的位置所形成的树，注意叶节点的位置以及它的原位置。
可以证明，对于每一个系数，它最后位置位于自身二进制的“倒数”（从 0 开始计数）。
所谓“倒数” 即将原二进制反转，得到新的二进制数字。
例如，处于第 $6=(110)_2$ 位置的系数最终位于 $3=(011)_2$ 位置。
那么，我们可以先将所有系数调整到最终的位置，随后合并。
我们按照这棵树来进行遍历，第一层循环遍历树的每一层，第二层循环遍历目前的合并位置，第三层循环则开始进行合并。
需要注意，我们表面上是在遍历一棵树，实际上我们是在遍历一个数组。
伪代码：
```
func FFT(a)
change_pos(a) // 将每个系数移动到最终位置
n = a.length
for layer from 1 to log2(n) // 遍历树的层数
len = power(2, layer) // 每一组的长度（一组即两个节点）
wn = complex(cos(2 * pi / len), sin(2 * pi / len)) // 计算单位根
for pos from 0 to n - 1 by len // 合并到的位置
w = complex(1, 0)
for i from 0 to len / 2 - 1 // 开始合并
t = w * a[pos + i + len / 2]
u = a[pos + i]
a[pos + i] = u + t
a[pos + i + len / 2] = u - t
w = w * wn
return a
```
## 快速傅里叶逆变换
我们将多项式转化为点值表示后，需要再将它转换回来，此时就需要快速傅里叶逆变换。
我们在快速傅里叶变换时，只需要将代入的单位根取倒数，并将得到的结果都除以 $n$，即是快速傅里叶逆变换。
而单位根的倒数，即是单位根的共轭复数。
相应的代码，由于思维难度不大，留作课后练习。
## 一些例题
[P1919 【模板】高精度乘法 | A\*B Problem 升级版 - 洛谷](https://www.luogu.com.cn/problem/P1919)
[P3803 【模板】多项式乘法（FFT） - 洛谷](https://www.luogu.com.cn/problem/P3803)
## 参考与致谢
本文章参考了 OI Wiki 、《算法导论》的相关部分，以及洛谷的部分文章。感谢前人所作的贡献，让我学会了这个算法并得以写出此文章。
文中提到的资料的链接：
- [OI Wiki](https://oi.wiki/)
- [《NOI 大纲（2023 年修订版）》](https://www.noi.cn/upload/resources/file/2023/03/15/1fa58eac9c412e01ce3c89c761058a43.pdf)
## 更新日志
2024/12/5：
- 使用 manim 更新了单位根的图片。
- 把引用链接移动到了文章最后。
- 修改了一些用词、修正了一些公式。
2024/12/6：
- 修改了文章格式。