FFT——从入门到入土
FFT 是一种可以在 \(O(n \log n)\) 的时间复杂度将多项式转为点值表达的算法。实际上， FFT 只是在求解方法上优化了 DFT（离散傅里叶变换）的过程，并没有提出新的理论。但是其高效的复杂度使得它被广泛使用。
阅读此文章前请先对复数有基本的了解。
定义
系数表示法
就是使用一个多项式的系数序列来表达这个多项式。
点值表示法
一个 \(n\) 次多项式可以由 \(n+1\) 个点来确定。
设 \(y_i = f(x_i)\) 则
多项式卷积
定义多项式 \(A(x) = \sum_{i=0} ^{n-1} a_ix^i\)，\(B(x)=\sum_{i=0}^{m}b_ix^i\)，它们的乘积为
实际上就是简单的直接相乘
单位复根
定义
\(n\) 次单位为满足 \(\omega^n = 1\) 的解 \(\omega\)，记作 \(\omega_n\)。
这 \(n\) 个单位根均匀分布在复平面的单位圆上。
自然地，\(w_n^k\) 的辐角为 \(\dfrac{2k\pi}{n}\)。因此，有
单位根的性质
-
\(\omega_n^k = \omega_{pn}^{pk}\)
等分成 \(pn\) 块相当于先分成 \(n\) 块后再将每个块分成 \(p\) 块，边界自然不变，
-
\(\omega_n^{k+\frac{n}{2}} = -\omega_n^k\)。
指数加 \(\frac{n}{2}\) 相当于绕原点旋转 \(\pi\) 。
-
\(\omega_n^p \cdot \omega_n^q = \omega_n^{p+q}\)
设 \(\alpha = \frac{2p\pi}{n}, \beta = \frac{2q\pi}{n}\)，则
\[\begin{aligned} \omega_{n}^{p} \times \omega_{n}^{q}&=(\cos \alpha \cos \beta-\sin \alpha \sin \beta)+i(\sin \alpha \cos \beta+\cos \alpha \sin \beta)\\ &=\cos (\alpha+\beta)+i \sin (\alpha+\beta)\\&=\omega_{n}^{p+q} \end{aligned} \]
离散傅里叶变换
我们发现，当两个多项式 \(A,B\) 同时取 \(x\) 时，得到的点分别为 \(y_A,y_B\)，\(A*B\) 取到的点即为 \((x,y_A \cdot y_B)\)。
这样做是 \(O(n)\) 的。
对于任意系数多项式转点值，我们可以随意取 \(x\) 的值。
但即使是这样，我们代入还是 \(O(n^2)\) 的复杂度。
快速傅里叶变换
傅里叶正变换
FFT 算法的基本思想是分治。
按时域抽取
具体的讲，我们将多项式分为奇数项和偶数项。
考虑对于多项式 \(F(x)\)，将其拆分成奇数项的多项式和偶数项的多项式 \(G(x)\) 和 \(H(x)\)，即为
这样拆分利用了单位根的性质，我们尝试将单位根代入求值：
证明：
由性质二可得
这样，我们就利用单位根，将原式在 \(\omega_n^{k}\) 的值转化成了两个规模更小的多项式的值，递归处理即可。
由递归式 \(T(n)=2T(\frac{n}{2})+n\) 得到 \(T(n)=O(n\log n)\)。
注意到这里我们每次都要对每个系数进行点乘，所以必须将多项式长度拓展到二的整数幂上。
按频域抽取
从另一个角度，我们直接把原多项式分为前半和后半：
注意到 \(\omega_n^{n/2}=-1\)，即 \(\omega_n^{kn/2}=(-1)^k\)，我们根据 \(k\) 的奇偶性分类，可以得到两个 \(\frac{n}{2}\) 项的多项式，递归处理即可。
这样，我们将问题规模成功缩小了一半，同样可以达到 \(O(n\log n)\) 的复杂度。
傅里叶逆变换
将点值表达式转回系数表达式的过程被称作傅里叶逆变换。
考虑 \(\text{DFT}\) 本质上是一个线性变换，我们可以给得到的新矩阵乘上一个逆矩阵，来得到初始矩阵。
观察上面的矩阵，不难看出其逆矩阵 \((b_{i,j})=\frac{(\overline{a_{i,j}})}{n}\)。
递归版代码实现
void fft(comp *x, int n, int type) { //时域
if (n == 1) return;
comp l[n >> 1], r[n >> 1];
for (int i = 0; i < n; i++) { //按奇偶分类
if (!(i & 1))
l[i >> 1] = x[i];
else
r[i >> 1] = x[i];
}
fft(l, n >> 1, type), fft(r, n >> 1, type);
comp wn1 = comp(cos(2 * pi / n), type * sin(2 * pi / n)), wnk = comp(1, 0);
for (int i = 0; i < (n >> 1); i++, wnk *= wn1;) {
x[i] = l[i] + wnk * r[i];
x[i + (n >> 1)] = l[i] - wnk * r[i]; // 计算左半部分
}
}
这份代码常数巨大，不推荐使用。我们需要一种常数更小的写法。
蝴蝶变换
我们每次递归的时候都要将系数分成两部分。为什么不可以在计算前就将其计算好呢？
以 \(n=8\) 时为例：
原本系数 \(F=[0,1,2,3,4,5,6,7]\)，变换后系数 \(\text{DFT}(F)=[0,4,2,6,1,5,3,7]\)。
写成二进制后发现，变换后的下标即为原数写成二进制再翻转后的数。
证明咕着。
for (int i = 0; i < len; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * (len >> 1));
使用时：
for (int i = 0; i < n; i++) if (i < rev[i]) swap(a[i], a[rev[i]]);
这样就免去了递归的大常数。
非递归代码实现
void fft(comp *x, int n, int tp) { // 时域
for (int i = 0; i <= n; i++) if (i < rev[i]) swap(x[i], x[rev[i]]);
for (int len = 1; len < n; len <<= 1) {
int sz = len * 2;
comp wn1 = comp(cos(PI / len), sin(PI / len) * tp);
for (int l = 0; l < n; l += sz) {
int r = l + len - 1;
comp wnk = 1;
for (int i = l; i <= r; i++) {
comp a = x[i], b = x[i + len];
x[i] = a + wnk * b, x[i + len] = a - wnk * b;
wnk *= wn1;
}
}
}
}
优化方法
省略蝴蝶变换
观察 按频域抽取 时进行的分类，手玩可以写出一个将偶数单位根放在左边，奇数单位根放在右边的算法（即按照 \(k\) 的奇偶性分类，偶数放左边奇数放右边）。
此时做一次蝴蝶变换即可得到按 \(\omega\) 升幂排序的点值表达。
但是考虑到按时域抽取的 IDFT 在开头同样要做一次蝴蝶变换，而做两次蝴蝶变换可以抵消，因此可以直接省略。
更详细的内容可以自行搜索转置原理。
三次变两次
不知道有啥用，咕。
原根
呃呃
快速数论变换
呃呃
浙公网安备 33010602011771号