快速傅里叶变换#

插值多项式#

通过点值表达式还原多项式的操作就是插值，这等价于求解一个线性方程组：

表示成矩阵形式就是

这里系数矩阵 $\begin{bmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^{n-1}\\ 1 & x_1 & x_1^2 & \cdots & x_1^{n-1}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & x_{n-1} & x_{n-1}^2 & \cdots & x_{n-1}^{n-1}\\ \end{bmatrix}$ 是范德蒙德矩阵，该矩阵的行列式值是 $\prod_{0\le i\lt j\lt n}(x_j-x_i)$，因为 $x_i(i=0,1,\cdots,n-1)$ 各不相同，所以该矩阵的行列式值 $\neq0$，即系数矩阵满秩，该线性方程组只能有唯一解。并且可以计算出多项式系数向量：

点值表示法下的多项式乘法#

假设多项式 $C(x)=A(x)B(x)$，其中 $A(x)$ 的度数（最高次幂）为 $n$，$B(x)$ 的度数为 $m$，此时 $C(x)$ 的度数 $\deg C = \deg A + \deg B = n+m$。如果要通过点值表示法插值出原多项式 $C(x)$，显然至少需要 $n+m+2$ 个点。因此，多项式 $A(x),B(x)$ 的点值表示都需要 $n+m+2$ 个点，设点集分别为 $S_A,S_B$：

对于任意点 $x_k$，$C(x_k)=A(x_k)B(x_k)$。因此，多项式 $C(x)$ 的点值表示为

利用点值表示法求多项式乘法的流程实际上就是：

1. 将给定的多项式 $A(x),B(x)$（设度数都是 $O(N)$ 量级）分别转变为点值表示
2. 通过 $A(x),B(x)$ 的点值表示计算出 $C(x)$ 的点值表示
3. 通过 $C(x)$ 的点值表示插值还原出系数多项式

其中第2步的时间复杂度显然为 $O(N)$，如果我们可以通过合理的点集选取，使得第一步和第三步都能够在 $O(N\log N)$ 的复杂度下实现，那么我们就得到了一个总复杂度为 $O(N\log N)$ 的多项式乘法算法。

单位复根#

$n$ 次单位复根 是满足 $\omega^n=1$ 的复数 $\omega$，这样的复根恰好有 $n$ 个：$\omega = e^{\frac{2\pi k i}{n}}(k=0,1,\cdots,n-1)$。根据欧拉公式可知，这些复根均匀地分布在复平面的单位圆上。

一些引理#

离散傅里叶变换(DFT)#

回想一下我们的问题：计算多项式乘法 $C(x)=A(x)B(x)$。为了 $O(N\log N)$ 解决这个问题，我们需要将系数表示转换为点值表示，这里我们对于点值表示法的选取就是 $n$ 次单位复根，即我们希望计算多项式

在 $\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}$ 处的值（即 $n$ 个 $n$ 次单位复根处）。假设 $A(x)$ 的系数表示用向量 $\vec a = (a_0,a_1,\cdots,a_{n-1})$ 表示，定义：

则向量 $\vec y = (y_0,y_1,\cdots,y_{n-1})$ 就是系数向量 $\vec a = (a_0,a_1,\cdots,a_{n-1})$ 的离散傅里叶变换（DFT），记作 $y=\text{DFT}_n(a)$。

这里的公式可以与傅里叶变换文中的离散傅里叶变换公式作比较，不难发现本质上是相同的。

快速傅里叶变换(FFT)#

本章节中，假设 $n$ 是 $2$ 的整数次幂。

**快速傅里叶变换（FFT）**实际上就是利用 $n$ 次单位复根的性质，在 $O(n\log n)$ 的时间复杂度下计算出 $n-1$ 阶多项式 $A(x) = \sum_{j=0}^{n-1} a_jx^j$ 的点值表示。按照下标的奇偶性，分别定义两个多项式 $A^{[0]}(x),A^{[1]}(x)$：

原多项式 $A(x) = A^{[0]}(x^2) + xA^{[1]}(x^2)$。

比如 $7$ 阶多项式 $f(x) = a_0+a_1x+a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6+a_7x^7$ 就处理为：

$$f(x) = (a_0+a_2x^2+a_4x^4+a_6x^6)+(a_1x+a_3x^3+a_5x^5+a_7x^7)$$

其中 $f^{[0]}(x) = a_0 + a_2x+a_4x^2+a_6x^3,f^{[1]}(x) = a_1+a_3x+a_5x^2+a_7x^3$。于是 $f(x) = f^{[0]}(x^2)+xf^{[1]}(x^2)$。

于是，求 $A(x)$ 在 $\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}$ 处的值就转化为：求两个新的多项式 $A^{[0]}(x),A^{[1]}(x)$ 在 $(\omega_n^0)^2,(\omega_n^1)^2,\cdots,(\omega_n^{n-1})^2$ 处的值，根据折半引理可知这 $n$ 个点的点值实际上只有 $\frac n 2$ 个不同的取值（因为 $(\omega_{n}^{k+n/2})^2 = (\omega_n^k)^2$），这就使得子问题的量级减半了，即求 $\text{DFT}_n$ 转化为了两个 $\text{DFT}{\frac n 2}$ 的子问题，时间复杂度如下：

快速傅里叶逆变换(Inverse FFT)#

求出多项式的点值表示后，我们还需要将点值表示逆运算为常见的系数表示，这一步就是离散傅里叶逆变换(IDFT)。前文中，我们已经提到了IDFT就是插值，即解方程

在FFT中，点值表示选用了 $n$ 次单位复根，上方的方程可以进一步写成

这个单位复根构成的逆矩阵可以表示为以下形式：

我们记 $W=\begin{bmatrix} \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0\\ \omega_n^0 & \omega_n^1 & \omega_n^2 & \cdots & \omega_n^{n-1}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ \omega_n^0 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)(n-1)}\\ \end{bmatrix},W^{-1}=\frac 1 n\begin{bmatrix} \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0\\ \omega_n^0 & \omega_n^{-1} & \omega_n^{-2} & \cdots & \omega_n^{-(n-1)}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ \omega_n^0 & \omega_n^{-(n-1)} & \omega_n^{-2(n-1)} & \cdots & \omega_n^{-(n-1)(n-1)}\\ \end{bmatrix}$，只需要验证 $WW^{-1}=I_n$，这里略去计算过程。

因此有

与公式 $y_k = \sum_{j=0}^{n-1}a_j\omega_n^{kj}$ 相比较，不难发现这两个问题几乎是一样的（计算 $\text{IFFT}_n(x)$ 只需要将向量 $y,a$ 互换，单位复根 $\omega_n$ 取逆）。

分治#

假设我们已知 $y^{[0]}=\text{DFT}(A^{[0]}),y^{[1]}=\text{DFT}(A^{[1]})$（这里 $y^{[0]},y^{[1]}$ 分别是两个长度为 $\frac n 2$ 的向量），求解 $y=\text{DFT}(A)$。也就是我们已知 $y^{[0]} = (A^{[0]}(\omega_{\frac n 2}^{0}), A^{[0]}(\omega_{\frac n 2}^{1}),\cdots, A^{[0]}(\omega_{\frac n 2}^{\frac n 2 - 1}))$ 和 $y^{[1]} = (A^{[1]}(\omega_{\frac n 2}^{0}), A^{[1]}(\omega_{\frac n 2}^{1}),\cdots, A^{[1]}(\omega_{\frac n 2}^{\frac n 2 - 1}))$，求 $y=(A(\omega_n^0),A(\omega_n^1),\cdots,A(\omega_n^{n-1}))$。

根据前文推导的公式 $A(x) = A^{[0]}(x^2) + xA^{[1]}(x^2)$ 可知 $0\le k\lt \frac n 2$ 时：

而后一半的计算则略有不同（$0\le k\lt \frac n 2$）：

然后就可以写出FFT的代码了：

1
using cp = complex<double>;
2
void fft(vector<cp>& a, int inv) {
3
    int n = a.size();
4
    if (n == 1) return;
5
    vector<cp> a0(n / 2), a1(n / 2);
6
    for (int i = 0; i * 2 < n; i++) {
7
        a0[i] = a[i * 2];
8
        a1[i] = a[i * 2 + 1];
9
    }
10
    fft(a0, inv);
11
    fft(a1, inv);
12

13
    double angle = 2 * pi / n * (inv ? -1 : 1);
14
    cp w(1, 0), wn(cos(angle), sin(angle));
15
    for (int i = 0; i * 2 < n; i++) {
16
        a[i] = a0[i] + w * a1[i];
17
        a[i + n / 2] = a0[i] - w * a1[i];
18
        if (inv) {
19
            a[i] /= 2;
20
            a[i + n / 2] /= 2;
21
        }
22
        w *= wn;
23
    }
24
}

上方函数中 inv 取 $0$ 时就是一次FFT的正变换；否则是FFT的逆变换。

下面是模板题SPOJ - POLYMUL的一个AC代码：

1
#include <bits/stdc++.h>
2
using namespace std;
3
using db = double;
4
const db pi = acos(-1);
5
using cp = complex<db>;
6

7
void fft(vector<cp>& a, int inv) {
8
    int n = a.size();
9
    if (n == 1) return;
10
    vector<cp> a0(n / 2), a1(n / 2);
11
    for (int i = 0; i * 2 < n; i++) {
12
        a0[i] = a[i * 2];
13
        a1[i] = a[i * 2 + 1];
14
    }
15
    fft(a0, inv);
16
    fft(a1, inv);
17

18
    db angle = 2 * pi / n * (inv ? -1 : 1);
19
    cp w(1, 0), wn(cos(angle), sin(angle));
20
    for (int i = 0; i * 2 < n; i++) {
21
        a[i] = a0[i] + w * a1[i];
22
        a[i + n / 2] = a0[i] - w * a1[i];
23
        if (inv) {
24
            a[i] /= 2;
25
            a[i + n / 2] /= 2;
26
        }
27
        w *= wn;
28
    }
29
}
30

31
vector<long long> multiply(vector<int>& a, vector<int>& b) {
32
    int n = 1;
33
    vector<cp> fa(a.begin(), a.end()), fb(b.begin(), b.end());
34
    while (n < a.size() + b.size()) n <<= 1;
35
    fa.resize(n);
36
    fb.resize(n);
37

38
    fft(fa, 0), fft(fb, 0);
39

40
    for (int i = 0; i < n; i++)
41
        fa[i] *= fb[i];
42

43
    fft(fa, 1);
44
    vector<long long> res(n);
45
    for (int i = 0; i < n; i++)
46
        res[i] = round(fa[i].real());
47
    return res;
48
}
49

50
int main() {
51
    ios::sync_with_stdio(false);
52
    cin.tie(nullptr);
53
    cout.tie(nullptr);
54
    int t;
55
    cin >> t;
56
    while (t--) {
57
        int n;
58
        cin >> n;
59
        vector<int> a(n + 1), b(n + 1);
60
        for (auto& i : a)
61
            cin >> i;
62
        for (auto& i : b)
63
            cin >> i;
64
        auto c = multiply(a, b);
65
        for (int i = 0; i < n * 2 + 1; i++)
66
            cout << c[i] << " \n"[i == n * 2];
67
    }
68
    return 0;
69
}

优化（倍增FFT）#

下图展示了一个 $7$ 阶多项式的系数在每一轮递归后的位置：

注意每个系数的下标在最终状态（树的叶子结点）时的位置，$(0,1,2,3,4,5,6,7)\rightarrow (0,4,2,6,1,5,3,7)$。我们从二进制的角度找找规律：$(000,001,010,011,100,101,110,111)\rightarrow(000,100,010,110,001,101,011,111)$。观察到：把原始序列下标的二进制翻转对称一下，就是最终那个位置的下标。这个性质被称为是位逆序置换（bit-reversal permutation）。

位逆序置换显然可以在 $O(n\log n)$ 的复杂度下实现，但也有更优的 $O(n)$ 做法，设 $R(x)$ 表示二进制数 $x$ 翻转后的数，从小到大递推 $R(x)$。

• 首先 $R(0)=0$。
• 因为我们从小到大递推 $R(x)$ 的值，因此在求 $R(x)$ 时，$R(\lfloor \frac x2\rfloor)$ 的值已知。只要将 $x$ 右移一位，然后翻转，再右移一位，就得到了 $x$ 除了二进制最高位之外所有位的翻转结果，而 $R(x)$ 的最高位可以通过 $x\bmod 2$ 求出。于是 $$R(x) = \bigg\lfloor\frac{R(\lfloor \frac x2\rfloor)}{2}\bigg\rfloor + (x\bmod 2)\times 2^k$$ 这里 $2^k$ 表示将 $x\bmod 2$ 移动至最高位的偏移量。

代码如下：

1
int bit_reorder(int n, vector<cp>& a) {
2
    int len = __builtin_ctz(n);
3
    if ((int)rev.size() != n) {
4
        rev.assign(n, 0);
5
        for (int i = 0; i < n; i++) {
6
            rev[i] = (rev[i >> 1] >> 1) + ((i & 1) << (len - 1));
7
        }
8
    }
9
    for (int i = 0; i < n; i++) {
10
        if (i < rev[i]) {
11
            swap(a[i], a[rev[i]]);
12
        }
13
    }
14
    return len;
15
}

也就是说我们能够直接找到FFT递归树最后一层的状态，然后自下而上地合并两个叶子，每一次合并实际上都是在做“已知 $y^{[0]} = (A^{[0]}(\omega_{\frac n 2}^{0}), A^{[0]}(\omega_{\frac n 2}^{1}),\cdots, A^{[0]}(\omega_{\frac n 2}^{\frac n 2 - 1}))$ 和 $y^{[1]} = (A^{[1]}(\omega_{\frac n 2}^{0}), A^{[1]}(\omega_{\frac n 2}^{1}),\cdots, A^{[1]}(\omega_{\frac n 2}^{\frac n 2 - 1}))$，求 $y=(A(\omega_n^0),A(\omega_n^1),\cdots,A(\omega_n^{n-1}))$”这个子问题，从前文已经推出了这个问题的解法：

注意到，如果我们把向量 $y^{[0]}$ 和 $y^{[1]}$ 依次排列，即 $(y^{[0]},y^{[1]})=(y_0^{[0]},y_1^{[0]},\cdots,y_{\frac n2-1}^{[0]},y_0^{[1]},y_1^{[1]},\cdots,y_{\frac n2-1}^{[1]})$，对比所求向量 $y=(y_0,y_1,\cdots,y_\frac n2,y_{\frac n2+1},\cdots,y_{n-1})$，不难发现 $y_{k}^{[0]},y_{k}^{[1]}$ 和 $y_k,y_{k+\frac n2}$ 的位置两两对应，因此我们可以直接在原数组上计算傅里叶变换，无需递归。这个优化被称为是蝴蝶变换（butterfly transform）。

1
void fft(vector<cp>& a, int inv) {
2
    int n = a.size();
3
    int len = bit_reorder(n, a);
4

5
    for (int k = 1; k < n; k <<= 1) {
6
        db angle = 2 * pi / (k * 2) * (inv ? -1 : 1);
7
        cp wn(cos(angle), sin(angle));
8
        for (int i = 0; i < n; i += k * 2) {
9
            cp w(1, 0);
10
            for (int j = i; j < i + k; j++) {
11
                cp yk0 = a[j], yk1 = a[j + k] * w;
12
                a[j] = yk0 + yk1;
13
                a[j + k] = yk0 - yk1;
14
                w *= wn;
15
            }
16
        }
17
    }
18
    if (inv) {
19
        for (auto& i : a)
20
            i /= n;
21
    }
22
}

下面是模板题SPOJ - POLYMUL的倍增法AC代码：

1
#include <bits/stdc++.h>
2
using namespace std;
3
using db = double;
4
const db pi = acos(-1);
5
using cp = complex<db>;
6
vector<int> rev;
7
int bit_reorder(int n, vector<cp>& a) {
8
    int len = __builtin_ctz(n);
9
    if ((int)rev.size() != n) {
10
        rev.assign(n, 0);
11
        for (int i = 0; i < n; i++) {
12
            rev[i] = (rev[i >> 1] >> 1) + ((i & 1) << (len - 1));
13
        }
14
    }
15
    for (int i = 0; i < n; i++) {
16
        if (i < rev[i]) {
17
            swap(a[i], a[rev[i]]);
18
        }
19
    }
20
    return len;
21
}
22

23
void fft(vector<cp>& a, int inv) {
24
    int n = a.size();
25
    int len = bit_reorder(n, a);
26

27
    for (int k = 1; k < n; k <<= 1) {
28
        db angle = 2 * pi / (k * 2) * (inv ? -1 : 1);
29
        cp wn(cos(angle), sin(angle));
30
        for (int i = 0; i < n; i += k * 2) {
31
            cp w(1, 0);
32
            for (int j = i; j < i + k; j++) {
33
                cp yk0 = a[j], yk1 = a[j + k] * w;
34
                a[j] = yk0 + yk1;
35
                a[j + k] = yk0 - yk1;
36
                w *= wn;
37
            }
38
        }
39
    }
40
    if (inv) {
41
        for (auto& i : a)
42
            i /= n;
43
    }
44
}
45

46
vector<long long> multiply(vector<int>& a, vector<int>& b) {
47
    int n = 1;
48
    vector<cp> fa(a.begin(), a.end()), fb(b.begin(), b.end());
49
    while (n < a.size() + b.size()) n <<= 1;
50
    fa.resize(n);
51
    fb.resize(n);
52

53
    fft(fa, 0), fft(fb, 0);
54

55
    for (int i = 0; i < n; i++)
56
        fa[i] *= fb[i];
57

58
    fft(fa, 1);
59
    vector<long long> res(n);
60
    for (int i = 0; i < n; i++)
61
        res[i] = round(fa[i].real());
62
    return res;
63
}
64

65
int main() {
66
    ios::sync_with_stdio(false);
67
    cin.tie(nullptr);
68
    cout.tie(nullptr);
69
    int t;
70
    cin >> t;
71
    while (t--) {
72
        int n;
73
        cin >> n;
74
        vector<int> a(n + 1), b(n + 1);
75
        for (auto& i : a)
76
            cin >> i;
77
        for (auto& i : b)
78
            cin >> i;
79
        auto c = multiply(a, b);
80
        for (int i = 0; i < n * 2 + 1; i++)
81
            cout << c[i] << " \n"[i == n * 2];
82
    }
83
    return 0;
84
}