FFT——从入门到入土

FFT 是一种可以在 $O(n \log n)$ 的时间复杂度将多项式转为点值表达的算法。实际上， FFT 只是在求解方法上优化了 DFT（离散傅里叶变换）的过程，并没有提出新的理论。但是其高效的复杂度使得它被广泛使用。

阅读此文章前请先对复数有基本的了解。

定义#

系数表示法#

就是使用一个多项式的系数序列来表达这个多项式。

$$f(x) = \sum_{i=0}^{n-1} a_ix^i \Leftrightarrow f(x) = \{a_0, a_1, \dots , a_{n-1}\}$$

点值表示法#

一个 $n$ 次多项式可以由 $n+1$ 个点来确定。

设 $y_i = f(x_i)$ 则

$$f(x) = \sum_{i=0}^{n-1} a_ix^i \Leftrightarrow f(x) = \{(x_0,y_0),(x_1,y_1),\dots,(x_{n-1},y_{n-1})\}$$

多项式卷积#

定义多项式 $A(x) = \sum_{i=0} ^{n-1} a_ix^i$，$B(x)=\sum_{i=0}^{m}b_ix^i$，它们的乘积为

$$(A * B)(x) = \sum_{i=0}^{n+m-2} x^i \sum_{k=0}^ia_kb_{i-k}$$

实际上就是简单的直接相乘

单位复根#

定义#

$n$ 次单位为满足 $\omega^n = 1$ 的解 $\omega$，记作 $\omega_n$。

这 $n$ 个单位根均匀分布在复平面的单位圆上。

自然地，$w_n^k$ 的辐角为 $\dfrac{2k\pi}{n}$。因此，有

$$\omega_n^k = \cos\dfrac{2k\pi}{n} + i\cdot\sin \dfrac{2k\pi}{n},k\in\Z$$

单位根的性质#

1. $\omega_n^k = \omega_{pn}^{pk}$

   等分成 $pn$ 块相当于先分成 $n$ 块后再将每个块分成 $p$ 块，边界自然不变，

2. $\omega_n^{k+\frac{n}{2}} = -\omega_n^k$。

   指数加 $\frac{n}{2}$ 相当于绕原点旋转 $\pi$ 。

3. $\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)$，即为

$$F(x)=G(x^2)+x \cdot H(x^2)$$

这样拆分利用了单位根的性质，我们尝试将单位根代入求值：

证明：

$$\begin{aligned} F(\omega_n^k)&=G[(\omega_n^k)^2]+\omega_n^k \cdot H[(\omega_n^k)^2] \\ &= G(\omega_n^{2k}) + \omega_n^k \cdot H(\omega_n^{2k}) \\ &= G(\omega_{n/2}^k)+\omega_n^k \cdot H(\omega_{n/2}^k) \end{aligned}$$

由性质二可得

$$F(\omega_n^{k+n/2})=G(\omega_{n/2}^k)-\omega_n^k \cdot H(\omega_{n/2}^k)$$

这样，我们就利用单位根，将原式在 $\omega_n^{k}$ 的值转化成了两个规模更小的多项式的值，递归处理即可。

由递归式 $T(n)=2T(\frac{n}{2})+n$ 得到 $T(n)=O(n\log n)$。

注意到这里我们每次都要对每个系数进行点乘，所以必须将多项式长度拓展到二的整数幂上。

按频域抽取#

从另一个角度，我们直接把原多项式分为前半和后半：

$$\begin{aligned} F(\omega_n^k)&= \sum_{i=0}^{n/2 -1} a_i \cdot \omega_n^{ik} + a_{i+n/2} \cdot \omega_{n}^{k(i+n/2)} \\ &=\sum_{i=0}^{n/2-1} (a_i+a_{i+n/2}\cdot \omega_n^{kn/2})\cdot \omega_n^{ik} \end{aligned}$$

注意到 $\omega_n^{n/2}=-1$，即 $\omega_n^{kn/2}=(-1)^k$，我们根据 $k$ 的奇偶性分类，可以得到两个 $\frac{n}{2}$ 项的多项式，递归处理即可。

这样，我们将问题规模成功缩小了一半，同样可以达到 $O(n\log n)$ 的复杂度。

傅里叶逆变换#

将点值表达式转回系数表达式的过程被称作傅里叶逆变换。

考虑 $\text{DFT}$ 本质上是一个线性变换，我们可以给得到的新矩阵乘上一个逆矩阵，来得到初始矩阵。

$$\begin{bmatrix} F(\omega_n^0)\\ F(\omega_n^1) \\ \vdots\\ F(\omega_n^{n-1}) \end{bmatrix} = \begin{bmatrix} 1 &1&\cdots&1 \\ 1&\omega_n^1&\cdots&\omega_n^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n-1} &\cdots &\omega_{n}^{(n-1)^2} \end{bmatrix} \begin{bmatrix} a_0\\a_1\\\vdots\\a_{n-1} \end{bmatrix}$$

观察上面的矩阵，不难看出其逆矩阵 $(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 在开头同样要做一次蝴蝶变换，而做两次蝴蝶变换可以抵消，因此可以直接省略。

更详细的内容可以自行搜索转置原理。

三次变两次#

不知道有啥用，咕。

原根#

呃呃

快速数论变换#

呃呃