快速傅里叶变换（FFT）学习笔记

#前置知识

#多项式

#系数表示与点值表示

多项式的系数表示，设 $A(x)$ 表示一个 $n - 1$ 次多项式，所有项的系数组成的 $n$ 维向量 $(a_0, a_1, a_2, \dots, a_{n - 1})$ 唯一确定了这个多项式。

$$\begin{align*} A(x) & = \sum_{i = 0}^{n - 1} a_i x^i \\ & = a_0 + a_1 x + a_2 x^2 + \dots + a_{n - 1} x^{n - 1} \\ \end{align*}$$

多项式的点值表示，将一组 互不相同 的 $n$ 个 $x$ 带入多项式，得到的 $n$ 个值。设它们组成的 $n$ 维向量分别为 $(x_0, x_1, x_2, \dots, x_{n - 1})$、$(y_0, y_1, y_2, \dots, y_{n - 1})$。

$$\begin{align*} y_i &= A(x_i) \\ &= \sum_{j = 0}^{n - 1} a_j x_i^j \\ \end{align*}$$

#求值与插值

定理：一个 $n - 1$ 次多项式在 $n$ 个不同点的取值唯一确定了该多项式。

证明：假设命题不成立，存在两个不同的 $n - 1$ 次多项式 $A(x)$、$B(x)$，满足对于任何 $i \in [0, n - 1]$，有 $A(x_i) = B(x_i)$。

令 $C(x) = A(x) - B(x)$，则 $C(x)$ 也是一个 $n - 1$ 次多项式。对于任何 $i \in [0, n - 1]$，有 $C(x_i) = 0$。

即 $C(x)$ 有 $n$ 个根，这与代数基本定理（一个 $n - 1$ 次多项式在复数域上有且仅有 $n - 1$ 个根）相矛盾，故 $C(x)$ 并不是一个 $n - 1$ 次多项式，原命题成立，证毕。

如果我们按照定义求一个多项式的点值表示，时间复杂度为 $O(n^2)$。

已知多项式的点值表示，求其系数表示，可以使用 插值。使用待定系数法的插值算法时间复杂度为 $O(n^3)$。

#多项式乘法

我们定义两个多项式 $A(x) = \sum_{i = 0}^{n - 1} a_i x^i$ 与 $B(x) = \sum_{i = 0}^{n - 1} b_i x^i$ 相乘的结果为 $C(x)$（假设两个多项式次数相同，若不同可在后面补零）。

$$C(x) = A(x) \times B(x) = \sum_{k = 0}^{2n - 2} (\sum_{k = i + j} a_i b_j) x^k$$

两个 $n - 1$ 次多项式相乘，得到的是一个 $2n - 2$ 次多项式，时间复杂度为 $O(n^2)$。

如果使用两个多项式在 $2n - 1$ 个点处取得的点值表示，那么

$${y_3}_i = (\sum_{j = 0}^{2n - 1} a_j x_i^j) \times (\sum_{j = 0}^{2n - 1} b_j x_i^j) = {y_1}_i \times {y_2}_i$$

时间复杂度为 $O(n)$。

#复数

设 $a$、$b$ 为实数，$i^2 = -1$，形如 $a + bi$ 的数叫做 复数，其中 $i$ 被称为 虚数单位。复数域是已知最大的域。

关于复数的更多内容详见人民教育出版社出版的《普通高中教科书 数学（必修 第二册）》。

#复平面

在复平面中，$x$ 轴代表实数、$y$ 轴（除原点外的所有点）代表虚数。每一个复数 $a + bi$ 对应复平面上一个从 $(0, 0)$ 指向 $(a, b)$ 的向量。

该向量的长度 $\sqrt{a^2 + b^2}$ 叫做模长。表示从 $x$ 轴正半轴到该向量的转角的有向（以逆时针为正方向）角叫做幅角。

复数相加遵循平行四边形定则。

复数相乘时，模长相乘，幅角相加。

#单位根

下文中，如不特殊指明，均设 $n$ 为 $2$ 的正整数次幂。

在复平面上，以原点为圆心，$1$ 为半径作圆，所得的圆叫做单位圆。以原点为起点，单位圆的 $n$ 等分点为终点，作 $n$ 个向量。设所得的 幅角为正且最小 的向量对应的复数为 $\omega_n$，称为 $n$ 次单位根。

由复数乘法的定义（模长相乘，幅角相加）可知，其与的 $n - 1$ 个向量对应的复数分别为 $\omega_n^2, \omega_n^3, \dots, \omega_n^n$，其中 $\omega_n^n = \omega_n^0 = 1$。

单位根的幅角为周角的 $1 \over n$，这为我们提供了一个计算单位根及其幂的公式：

$$\omega_n^k = \cos k \frac{2 \pi}{n} + i\sin k \frac{2 \pi}{n}$$

#单位根的性质

性质一：$\omega_{2n}^{2k} = \omega_n^k$

从几何意义上看，在复平面上，二者表示的向量终点相同。

更具体的，有

$$\cos 2k \frac{2 \pi}{2n} + i\sin 2k \frac{2 \pi}{2n} = \cos k \frac{2 \pi}{n} + i \sin k \frac{2 \pi}{n}$$

性质二：$\omega_n^{k + \frac{n}{2}} = -\omega_n^k$

等式左边相当于 $\omega_n^k$ 乘上 $\omega_n^{\frac{n}{2}}$，考虑其值

$$\begin{align*} \omega_n^{ \frac{n}{2} } &= \cos \frac{n}{2} \cdot \frac{2 \pi}{n} + i\sin \frac{n}{2} \cdot \frac{2 \pi}{n} \\ &= \cos \pi + i\sin \pi \\ &= -1 \end{align*}$$

#离散傅里叶变换（DFT）

离散傅里叶变换（Discrete Fourier Transform，DFT）是一种将系数表达转换为点值表达的变换。这一变换可以求出一个 $n$ 次多项式在每个 $n$ 次单位根下的点值。

将 $n$ 次单位根的 $0$ 到 $n - 1$ 次幂带入多项式的系数表示，所得点值向量 $(A(\omega_n^0), A(\omega_n^1), \dots, A(\omega_n^{n - 1}))$ 称为其系数向量 $(a_0,\ a_1,\ \dots,\ a_{n - 1})$ 的离散傅里叶变换。

这个过程是 $O(n^2)$ 的。

#快速傅里叶变换（FFT）

FFT 是一种高效实现 DFT 的算法，其时间复杂度为 $O(n \log n)$。

#递归分治

考虑将多项式按照系数下标的奇偶分为两部分：

$$A(x) = (a_0 + a_2 x^2 + a_4 x^4 + \dots + a_{n - 2} x^{n - 2}) + (a_1 x + a_3 x^3 + a_5 x^5 + \dots + a_{n - 1} x^{n - 1})$$

令

$$\begin{align*} A_1(x) &= a_0 + a_2 x + a_4 x^2 + \dots + a_{n - 2} x^{\frac{n}{2} - 1} \\ A_2(x) &= a_1 + a_3 x + a_5 x^2 + \dots + a_{n - 1} x^{\frac{n}{2} - 1} \\ \end{align*}$$

则有

$$A(x) = A_1(x^2) + x A_2(x^2)$$

假设 $k < \frac{n}{2}$，现在要求 $A(\omega_n^k)$

$$\begin{align*} A(\omega_n^k) &= A_1(\omega_n^{2k}) + \omega_n^{k} A_2(\omega_n^{2k}) \\ &= A_1(\omega_{\frac{n}{2}}^{k}) + \omega_n^{k} A_2(\omega_{\frac{n}{2}}^{k}) \\ \end{align*}$$

这一步转化利用了单位根的性质一。

对于 $A(\omega_n^{k + \frac{n}{2}})$

$$\begin{align*} A(\omega_n^{k + \frac{n}{2}}) &= A_1(\omega_n^{2k + n}) + \omega_n^{k + \frac{n}{2}} A_2(\omega_n^{2k + n}) \\ &= A_1(\omega_n^{2k} \times \omega_n^n) - \omega_n^{k} A_2(\omega_n^{2k} \times \omega_n^n) \\ &= A_1(\omega_n^{2k}) - \omega_n^{k} A_2(\omega_n^{2k}) \\ &= A_1(\omega_{\frac{n}{2}}^{k}) - \omega_n^{k} A_2(\omega_{\frac{n}{2}}^{k}) \\ \end{align*}$$

这一步转化除性质一外，还利用到了性质二和 $\omega_n^n = 1$ 这一显然的结论。

注意到，当 $k$ 取遍 $[0, \frac{n}{2} - 1]$ 时，$k$ 和 $k + \frac{n}{2}$ 取遍了 $[0, n - 1]$。

也就是说，如果已知 $A_1(x)$ 和 $A_2(x)$ 在 $\omega_{\frac{n}{2}}^0, \omega_{\frac{n}{2}}^1, \dots, \omega_{\frac{n}{2}}^{\frac{n}{2} - 1}$ 处的点值，就可以在 $O(n)$ 的时间内求得 $A(x)$ 在 $\omega_n^0, \omega_n^1, \dots, \omega_n^{n - 1}$ 处的取值。而关于 $A_1(x)$ 和 $A_2(x)$ 的问题都是相对于原问题规模缩小了一半的子问题，分治的边界为一个常数项 $a_0$。

根据主定理，该分治算法的时间复杂度为

$$T(n) = 2T(\frac{n}{2}) + O(n) = O(n \log n)$$

这就是最常用的 FFT 算法 —— Cooley-Tukey 算法。

#迭代优化

递归实现的 FFT 效率不高，实际中一般采用迭代实现。

#二进制位翻转

考虑递归 FFT 分治到边界时，每个数的顺序，及其二进制位。

000 001 010 011 100 101 110 111
 0   1   2   3   4   5   6   7
 0   2   4   6 - 1   3   5   7
 0   4 - 2   6 - 1   5 - 3   7
 0 - 4 - 2 - 6 - 1 - 5 - 3 - 7
000 100 010 110 001 101 011 111

发现规律，分治到边界后的下标等于原下标的二进制位翻转。

代码实现，枚举每个二进制位即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14

// assert(a.size() == 1 << std::__lg(a.size()));
int k = std::__lg(a.size());

for (int i = 0; i < a.size(); i++) {
    int t = 0;

    for (int j = 0; j < k; j++) {
        if (i & (1 << j)) {
            t |= 1 << (k - j - 1);
        }
    }

    if (i < t) std::swap(a[i], a[t]);
}

#蝴蝶操作

考虑合并两个子问题的过程，假设 $A_1(\omega_{ \frac{n}{2} }^k)$ 和 $A_2(\omega_{ \frac{n}{2} }^k)$ 分别存在 $a_k$ 和 $a_{\frac{n}{2} + k}$ 中，$A(\omega_n^{k})$ 和 $A(\omega_n^{k + \frac{n}{2} })$ 将要被存放在 $b_k$ 和 $b_{\frac{n}{2} + k}$ 中，合并的单位操作可表示为

$$\begin{align*} b_k &= a_k + \omega_n^k \times a_{\frac{n}{2} + k} \\ b_{\frac{n}{2} + k} &= a_k - \omega_n^k \times a_{\frac{n}{2} + k} \\ \end{align*}$$

考虑加入一个临时变量 $t$，使得这个过程可以在原地完成，不需要另一个数组 $b$，也就是说，将 $A(\omega_n^{k})$ 和 $A(\omega_n^{k + \frac{n}{2} })$ 存放在 $a_k$ 和 $a_{\frac{n}{2} + k}$ 中，覆盖原来的值

$$\begin{align*} t &= \omega_n^k \times a_{\frac{n}{2} + k} \\ a_{\frac{n}{2} + k} &= a_k - t \\ a_k &= a_k + t \\ \end{align*}$$

这一过程被称为 蝴蝶操作。

#离散傅里叶逆变换（IDFT）

将点值表示的多项式转化为系数表示，同样可以使用快速傅里叶变换，这个过程叫做 傅里叶逆变换（Inverse Discrete Fourier Transform，IDFT）。

设 $(y_0, y_1, y_2, \dots, y_{n - 1})$ 为 $(a_0, a_1, a_2, \dots, a_{n - 1})$ 的傅里叶变换。考虑另一个向量 $(c_0, c_1, c_2, \dots, c_{n - 1})$ 满足

$$c_k = \sum_{i=0}^{n-1} y_i (\omega_n^{-k})^i$$

即多项式 $B(x) = y_0 + y_1 x + y_2 x^2 + \dots + y_{n - 1} x^{n - 1}$ 在 $\omega_n^0, \omega_n^{-1}, \omega_n^{-2}, \dots, \omega_n^{-(n - 1)}$ 处的点值表示。

将上式展开，得

$$\begin{align*} c_k &= \sum_{i = 0}^{n - 1} y_i (\omega_n^{-k})^i \\ &= \sum_{i = 0}^{n - 1} (\sum_{j = 0}^{n - 1} a_j (\omega_n^i)^j) (\omega_n^{-k})^i \\ &= \sum_{i = 0}^{n - 1} (\sum_{j = 0}^{n - 1} a_j (\omega_n^j)^i) (\omega_n^{-k})^i \\ &= \sum_{i = 0}^{n - 1} (\sum_{j = 0}^{n - 1} a_j (\omega_n^j)^i (\omega_n^{-k})^i) \\ &= \sum_{i = 0}^{n - 1} \sum_{j = 0}^{n - 1} a_j (\omega_n^{j - k})^i \\ &= \sum_{j = 0}^{n - 1} a_j (\sum_{i = 0}^{n - 1} (\omega_n^{j - k})^i) \\ \end{align*}$$

考虑一个式子

$$S(\omega_n^k) = 1 + \omega_n^k + (\omega_n^k)^2 + \dots + (\omega_n^k)^{n - 1}$$

当 $k \neq 0$ 时，两边同时乘上 $\omega_n^k$ 得

$$\omega_n^k S(\omega_n^k) = \omega_n^k + (\omega_n^k)^2 + (\omega_n^k)^3 + \dots + (\omega_n^k)^n$$

两式相减，整理后得

$$\begin{align*} \omega_n^k S(\omega_n^k) - S(\omega_n^k) &= (\omega_n^k)^n - 1 \\ S(\omega_n^k) &= \frac{(\omega_n^k)^n - 1}{\omega_n^k - 1} \end{align*}$$

分子为零，分母不为零，所以

$$S(\omega_n^k) = 0$$

当 $k = 0$ 时，显然 $S(\omega_n^k) = n$。

继续考虑上式

$$\begin{align*} c_k &= \sum_{j = 0}^{n - 1} a_j (\sum_{i = 0}^{n - 1} (\omega_n^{j - k})^i) \\ &= \sum_{j = 0}^{n - 1} a_j S(\omega_n^{j - k}) \end{align*}$$

当 $j = k$ 时，$S(\omega_n^{j - k}) = n$，否则 $S(\omega_n^{j - k}) = 0$，即

$$\begin{align*} c_i &= n a_i \\ a_i &= \frac{1}{n} c_i \end{align*}$$

所以，使用单位根的倒数代替单位根，做一次类似快速傅里叶变换的过程，再将结果每个数除以 $n$，即为傅里叶逆变换的结果。

但是这样实现起来有点麻烦，而 $w_n^{-k}$ 可以看做 $n$ 次本原单位根每次逆时针旋转本原单位根幅角的弧度，因此 $\omega_n^{-k}$ 和 $\omega_n^k$ 是一一对应的。具体的，$w_n^{-k} = w_n^{k + n}$。因此我们只需要使用 FFT 的方法，求出 $B(x)$ 在 $\omega_n$ 各个幂次下的值，然后数组反过来，即令 $a_k=\frac{1}{n} \sum_{i = 0}^n B(w_n^{n - k})$ 即可。

这一步快速计算插值的过程叫做 快速傅里叶逆变换（Inverse Fast Fourier Transform，IFFT）。

#代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

#include <cmath>
#include <complex>
#include <valarray>

const double PI = std::acos(-1);

void fast_fourier_transform(std::valarray<std::complex<double>>& a) {
    if (a.size() == 1) return;

    int m = a.size() >> 1;

    std::valarray<std::complex<double>>
        even = a[std::slice(0, m, 2)],
        odd = a[std::slice(1, m, 2)];

    fast_fourier_transform(even);
    fast_fourier_transform(odd);

    for (int i = 0; i < m; i++) {
        auto w = std::polar(1.0, -2 * PI * i / a.size()) * odd[i];
        a[i] = even[i] + w;
        a[i + m] = even[i] - w;
    }
}

void dft(std::vector<std::complex<double>>& a) {
    fast_fourier_transform(a);
}

void idft(std::vector<std::complex<double>>& a) {
    fast_fourier_transform(a);
    std::reverse(a.begin() + 1, a.end());
    std::transform(a.begin(), a.end(), a.begin(), [&](std::complex<double> x) {
        return static_cast<int>(std::round(x.real() / a.size()));
    });
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

#include <cmath>
#include <complex>
#include <vector>

const double PI = std::acos(-1);

void fast_fourier_transform(std::vector<std::complex<double>>& a) {
    if (a.size() == 1) return;

    // assert(a.size() == 1 << std::__lg(a.size()));
    int k = std::__lg(a.size());

    for (int i = 0; i < a.size(); i++) {
        int t = 0;

        for (int j = 0; j < k; j++) {
            if (i & (1 << j)) {
                t |= 1 << (k - j - 1);
            }
        }

        if (i < t) std::swap(a[i], a[t]);
    }

    for (int len = 2; len <= a.size(); len <<= 1) {
        std::complex<double> wlen(std::cos(2 * PI / len), std::sin(2 * PI / len));

        for (int i = 0; i < a.size(); i += len) {
            std::complex<double> w(1);

            for (int j = 0; j < len / 2; j++) {
                std::complex<double> u = a[i + j],
                                     v = a[i + j + len / 2] * w;

                a[i + j] = u + v;
                a[i + j + len / 2] = u - v;
                w *= wlen;
            }
        }
    }
}

void dft(std::vector<std::complex<double>>& a) {
    fast_fourier_transform(a);
}

void idft(std::vector<std::complex<double>>& a) {
    fast_fourier_transform(a);
    std::reverse(a.begin() + 1, a.end());
    std::transform(a.begin(), a.end(), a.begin(), [&](std::complex<double> x) {
        return static_cast<int>(std::round(x.real() / a.size()));
    });
}

#参考资料

本文的大部分内容转载自 FFT 学习笔记 - Menci’s OI Blog，在此表示感谢。此外，本文还参考了以下文章：

1. P3803 【模板】多项式乘法（FFT） 题解，—扶苏—，2019 年 12 月 25 日。
2. 快速傅里叶变换(FFT)详解，自为风月马前卒，2018 年 2 月 11 日。