【数学】【多项式】快速傅里叶变换（FFT）
引入
求两个多项式的卷积
Description
给定两个多项式 \(F\left(x\right), G\left(x\right)\) 的系数表示法，求两个多项式的卷积。
如：
得到的结果即为：
一种计算方法是将两个多项式中的每一项两两相乘，最后再将结果相加。
可惜，这样实现的复杂度不够优秀，为 \(\mathcal O\left(n^2\right)\) 级别。
因此我们需要寻找更快的计算方法。
多项式的点值表示和系数表示
我们知道，一个 \(n\) 次多项式，可以被 \(n + 1\) 个横坐标互不相同的点唯一确定。这样用 \(n + 1\) 个点的集合来表示一个多项式的方法被称为点值表示。对应的，用每一项的系数构成的集合表示一个多项式的方法就被称为系数表示。
例如方才的多项式 \(G\left(x\right)\)，系数表示： \(\{1, 2, 1\}\)，点值表示（不唯一）：\(\{\left(-1, 0\right), \left(0, 1\right), \left(1, 4\right)\}\)
这样，设 \(n\) 为 \(F\left(x\right)G\left(x\right)\) 的次数，下同，也就意味着我们可以通过找 \(n + 1\) 个 \(x\) 坐标，求出 \(F\left(x\right)\) 和 \(G\left(x\right)\) 在这，\(n + 1\) 个位置上的函数值，再将这 \(n - 1\) 个位置上的函数值分别相乘，就得到了 \(F\left(x\right)G\left(x\right)\) 的一个点值表示。
计算的复杂度是 \(\mathcal O\left(n\right)\) 级别，非常优秀。
求值和插值
实际应用上，我们很少通过点值表示来表示一个多项式。
因此，若想快速计算两个多项式的卷积，需要将两个多项式由系数表示转换为点值表示，然后在点值表示下计算两个多项式的卷积，再将点值表示转换为系数表示。
计算卷积的过程已经解决，现在需要的是实现将系数表示转换为点值表示，和将点值表示转换为系数表示的方法。这两个过程分别被称为求值和插值。
考虑暴力计算。找到 \(n + 1\) 个两两不同的 \(x\) 值，然后分别代入两个多项式里去，求出对应的点值，进行完乘法后，再用待定系数法，根据得到的点列出 \(n + 1\) 个方程，然后消元，或者使用插值法。复杂度？
消元复杂度：\(\mathcal O\left(n^3\right)\)；插值复杂度：单个 \(\mathcal O\left(n^2\right)\)，在 \(x\) 取值连续时可以优化到 \(\mathcal O\left(n\right)\)，由于 \(x\) 值可以任取，可以认为其复杂度为 \(\mathcal O\left(n^2\right)\)。没有起到优化作用。
快速傅里叶变换（FFT）
一种思路
考虑一个多项式 \(F\left(x\right) = x^2\)。发现因为它的函数图像具有对称性，因此只需要取 \(\dfrac{n}{2}\) 个 \(x\) 值，就可以得到 \(n\) 个点值。
那对于 \(F\left(x\right) = x^3\) 呢？显然还是有的。
那对于一个任意多项式 \(P\left(x\right)\) 呢？
（接下来我们默认多项式的项数为 \(2^n\) 形式，即便不是，也可以简单地将不足的项的系数视为 \(0\)。）
考虑将多项式 \(P\left(x\right)\)：
将其按奇偶次项分成两个多项式 \(P_e\left(x\right), P_o\left(x\right)\)：
则显然有：
看，发现了什么？我们只需要求出 \(P_e\left(x\right)\) 和 \(P_o\left(x\right)\)，就可以得到 \(P\left(x\right)\) 和 \(P\left(-x\right)\) 两处的多项式的点值！
而 \(P_e\left(x\right), P_o\left(x\right)\) 又是两个关于 \(x^2\) 的 \(\frac{n}{2}\) 次多项式，这似乎意味着我们可以在 \(\mathcal O\left(n \log n\right)\) 的时间复杂度级别下，递归地来实现这个算法！
但是，注意到这里面有什么问题吗？因为这个算法能够实现，是基于求 \(P\left(x\right)\) 的点值时 \(x\) 的取值能够正负配对的假设。而 \(x^2, x^4\) 等偶次项，不是正负配对的。这样递归就无法进行。
解决方案
都进行到这里了，有没有解决方案呢？答案是肯定的。
在实数域内，显然是不能使得 \(x\) 的取值总是能够两两配对的了。因此考虑把 \(x\) 的取值扩展到复数域。
我们知道， \(\mathbb i^2 = -1\)。因此若 \(x\) 的取值集合是 \(\{1， -1， \mathbb i, -\mathbb i\}\) 的话，\(x^2\) 就也满足两两配对了。
而这四个数又分别是什么呢？没错，就是 \(4\) 次单位根。
这便是 FFT 的基本思想。
快速傅里叶变换（DFT）
通过系数表示求点值表示的方式称为 DFT。
\(n\) 次单位根（\(\omega_n\)）可以通过如下的公式求得，其中 \(\alpha\) 是幅角，采用弧度制，证明自行搜索复数相关内容：
或者如下的公式求得，其中 \(e\) 是自然对数的底数：
根据单位根的性质有：
这样，想要求 \(P\left(x\right)\) 在 \(\{1, \omega_n, \omega_n^2, \cdots \omega_n^{\left(n - 1\right)}\}\) 时的点值，转而求 \(P_e\left(x\right)\) 和 \(P_o\left(x\right)\) 分别在 \(\{1, \omega_{\frac{n}{2}}, \omega_{\frac{n}{2}} ^ 2, \cdots, \omega_{\frac{n}{2}} ^ {\frac{n}{2} - 1}\}\) 时的点值，这样递归下去，在要求的多项式只有一项的时候返回系数，合并的时候把得到的 \(P_o\left(x\right)\) 的结果乘上一个当前的单位根，然后再分别一加一减计算出 \(P\left(x\right)\) 和 \(P\left(-x\right)\)。
下面会给出递归版的伪代码，不需要考虑怎么实现，只是用于理解 FFT 的基本思想即可。因为接下来还会介绍常数更小且更好写的迭代写法。
伪代码：
DFT(P{p_0, p_1, p_2, ..., p_n-1}):
if n = 1
返回 P
定义两个求值点子序列 Pe{p_0, p_2,...} Po{p_1, p_3,...}
DFT(Pe) DFT(Po)
计算 n 次单位根 wn
枚举 i
将 Pe_i + wn ^ i * Po_i 和 Pe_i - wn ^ i * Po_i 填入结果序列对应的位置
快速傅里叶逆变换（IDFT）
求值的问题解决了，我们现在来解决插值。
不难想到，插值应该和求值是基本相同的问题。
考虑 DFT 前的系数矩阵，让它乘上 DFT 矩阵后，变为点值矩阵，则 DFT 矩阵即运算为：
奇妙的，若找到一个矩阵，使得它乘上点值矩阵能够得到系数矩阵，则这个矩阵就是 IDFT 矩阵。
不难想到 IDFT 矩阵是 DFT 矩阵的逆矩阵。
那么 IDFT 矩阵是什么样子的呢？非常奇妙：
我自己认为的一个简单而且直观的证明：
证明过程可根据复数运算规则、单位根的性质、矩阵乘法得出。
也就是说，要实现 IDFT，只需要将原来 DFT 的代码中，改变函数名称，然后将单位根改成单位根的 \(-1\) 次方，最后把得出的结果除以 \(n\) 即可！
更快的实现
如果我们能够找到一种规律，使得每次计算的求值点对象有迹可循，我们就能够通过迭代的方式实现 FFT，从而省去递归调用的时间并减小常数，提高运行效率。
这种规律就在接下来要介绍的 FFT 的迭代实现中。
蝴蝶变换
考虑一个包含 \(8\) 个求值点的序列，手模一下不断奇偶分组的过程，可以得到这样的最终序列：
考虑下标的变化关系，没错，变换后序列的下标，就是原序列中对应数的二进制翻转后的下标，这种变换方式被称之为“蝴蝶变换”。
位逆序置换
考虑归纳法递推。
当要求 \(i\) 的位逆序置换时，假定 \(0\) 到 \(i - 1\) 的位逆序置换都已求出。
那么 \(\lfloor\dfrac{i}{2}\rfloor\) 的位逆序置换也必定求出。因此可以利用 \(\lfloor \dfrac{i}{2} \rfloor\) 的位逆序置换结果，得到 \(i\) 除了最低位上的位逆序置换。然后再判断最后一位上是否是 \(1\)，如果是，就在当前计算结果的最高位上补 \(1\) 即可。
边界条件：\(0\) 的位逆序置换恒为 \(0\)。
这样就可以在 \(\mathcal O\left(n\right)\) 的时间复杂度下求出所有位置的位逆序置换。
代码
int rev[Maxn];
inline void Init(int len) {
for(int i = 0; i < len; ++i) {
rev[i] = rev[i >> 1] >> 1;
if(i & 1) rev[i] |= len >> 1;
}
}
FFT 的迭代实现
至此，FFT 的迭代实现步骤已经很清楚了。
首先利用开始求出的位逆序置换结果，对整个序列进行蝴蝶变换，再通过循环的结构通过倍增模拟递归实现的步骤即可。
代码：
void FFT(CP a[], int len, int type) {//type = 1: DFT; type = -1: IDFT.
for(int i = 0; i < len; ++i) if(i < rev[i]) swap(a[i], a[rev[i]]);
for(int h = 2; h <= len; h <<= 1) {
CP wn(cos(2.0 * Pi / h), sin(2.0 * Pi / h * type));
for(int j = 0; j < len; j += h) {
CP wk(1.0, 0.0);
for(int k = j; k < j + h / 2; ++k) {
CP e = a[k], o = wk * a[k + h / 2];
a[k] = e + o, a[k + h / 2] = e - o;
wk = wk * wn;
}
}
}
if(type == -1) for(int i = 0; i < len; ++i) a[i].R /= len;
}
其它
完整代码
完整代码：（洛谷 P3803 多项式乘法（FFT））
#include <bits/stdc++.h>
#define LF long double
template <typename Temp> inline void read(Temp & res) {
Temp fh = 1; res = 0; char ch = getchar();
for(; !isdigit(ch); ch = getchar()) if(ch == '-') fh = -1;
for(; isdigit(ch); ch = getchar()) res = (res << 3) + (res << 1) + (ch ^ '0');
res = res * fh;
}
using namespace std;
const int Maxn = 2097200;
const LF Pi = 3.1415926535897932384626;
struct CP {
LF R, I;
CP() {R = I = 0;}
CP(LF A, LF B) {R = A; I = B;}
CP operator + (CP B) {return CP(R + B.R, I + B.I);}
CP operator - (CP B) {return CP(R - B.R, I - B.I);}
CP operator * (CP B) {return CP(R * B.R - I * B.I, R * B.I + I * B.R);}
};
int rev[Maxn];
inline void Init(int len) {
for(int i = 0; i < len; ++i) {
rev[i] = rev[i >> 1] >> 1;
if(i & 1) rev[i] |= len >> 1;
}
}
void FFT(CP a[], int len, int type) {
for(int i = 0; i < len; ++i) if(i < rev[i]) swap(a[i], a[rev[i]]);
for(int h = 2; h <= len; h <<= 1) {
CP wn(cos(2.0 * Pi / h), sin(2.0 * Pi / h * type));
for(int j = 0; j < len; j += h) {
CP wk(1.0, 0.0);
for(int k = j; k < j + h / 2; ++k) {
CP e = a[k], o = wk * a[k + h / 2];
a[k] = e + o, a[k + h / 2] = e - o;
wk = wk * wn;
}
}
}
if(type == -1) for(int i = 0; i < len; ++i) a[i].R /= len;
}
void polymul(CP a[], CP b[], int lenA, int lenB) {
int L = lenA + lenB, len = 1; while(L) {len <<= 1; L >>= 1;}
Init(len); FFT(a, len, 1); FFT(b, len, 1);
for(int i = 0; i < len; ++i) a[i] = a[i] * b[i];
FFT(a, len, -1);
}
int n, m, x;
CP A[Maxn], B[Maxn];
int main() {
read(n); read(m);
for(int i = 0; i <= n; ++i) {read(x); A[i].R = x;}
for(int i = 0; i <= m; ++i) {read(x); B[i].R = x;}
polymul(A, B, n, m);
for(int i = 0; i <= n + m; ++i) {printf("%d ", (int)round(A[i].R));}
return 0;
}
鸣谢
排名不分先后
远航之曲（图片）
B 站上某讲的特别清楚的视频
