多项式乘法 快速傅里叶变换（FFT）

前置知识1

1.多项式：一个以\(x\)为变量的多项式定义在一个代数域\(F\)上，将函数\(A(x)\)表示为形式和:

\[ A(x) \ =\ \sum _{i = 0} ^ {n - 1} a_i x^i \]

2.多项式的系数表示法；即由多项式的系数组成的向量 \(a\) $ = (a_0, \ a_1, \ ... \ a_{n - 1})$

2.多项式的点值表示法：即一个由\(n\)的点组成的集合

\[\{ \ (x_0, y_0), \ (x_1, y_1), \ ..., \ (x_{n - 1}, y_{n - 1}) \ \} \]

其中满足对于\(\forall \ i \ \in [0, \ n - 1]\)

\[ y_i = A(x_i) \]

3.点值表示法的多项式乘法：
已知两个多项式的点值表示：

\[\{ \ (x_0, y_0), \ (x_1, y_1), \ ..., \ (x_{n - 1}, y_{n - 1}) \ \} \\ \{ \ (x_0, y{'}_0), \ (x_1, y{'}_1), \ ..., \ (x_{n - 1}, y{'}_{n - 1}) \ \} \]

则两式相乘后的点值表示为：

\[ \{ \ (x_0, y_0 \cdot y{'}_0), \ (x_1, y_1 \cdot y{'}_1), \ ..., \ (x_{n - 1}, y_{n - 1} \cdot y{'}_{n - 1}) \ \} \]

4.；离散傅里叶变换（DFT）：我们称向量 \(y\) \(\ = \ (y_0, \ y_1, \ ... \ y_{n - 1})\) 为向量 \(a\) \(= (a_0, \ a_1, \ ... \ a_{n - 1})\)的离散傅里叶变换，记为 \(y\) \(= DFT(\) \(a\) \()\), 同时，\(a\) \(= DFT^{-1}(\) \(y\) \()\)

快速傅里叶变换（FFT）：

可以在\(O(n \ log \ n)\)的时间复杂度内求解两个次数界为\(n\)的多项式的乘法。算法主要由三步构成：

1.求值\(O(n \ log \ n)\)：将多项式由系数表示法转化为点值表示法，即\(DFT\)。

2.点值乘法\(O(n)\)：将两个多项式的点值相乘，得到所求多项式的点值表达式。

3.插值\(O(n \ log \ n)\)：将所求多项式由点值表示法转变为系数表示法，即\(DFT^{-1}\)。

一个次数界为\(n\)和\(m\)的多项式相乘，需要\(n + m\)个点的点值表达式，朴素的求法为\(O(n^2)\)，显然\(x\)的取值是任意的，我们可以通过选取特殊的\(x\)来降低求值的复杂度。

n次单位复数根：

是指满足\(w^n = 1\)的复数\(w\)，共有n个，分别为

\[w_n^{k} = e^{2 \pi ik / n} \qquad (k \in [0, n - 1]) \]

根据欧拉公式可得：

\[ w_n^k = cos\frac{2k \pi}{n} + i \cdot sin\frac{2k \pi}{n} \]

n次单位复数根有几个性质：

性质1： \(w_{pk}^{pn} = w_{k}^{n}\)

性质2： \(w_n^{k + \frac{n}{2}} = -w_n^k\)

性质3： 对任意\(n \geq 1\)和不能被\(n\)整除的非负整数\(k\)，有：

\[ \sum_{i = 0} ^ {n - 1} \ (w_n^k)^i = 0 \]

FFT的过程：

1.求值

(以下假定n为2的整数次幂)

对于多项式

\[A(x) = a_0 + a_1x + a_2x^2 + ... + a_{n - 1}x^{n - 1} \]

将其按照奇偶项拆分：

\[A_1(x) = a_0 + a_2x + a_4x^2 + ... + a_{n - 2}x^{\frac{n}{2} - 1} \\ A_2(x) = a_1 + a_3x + a_5x^2 + ... + a_{n - 1}x^{\frac{n}{2} - 1} \]

则：

\[A(x) = A_1(x^2) + xA_2(x^2) \]

求\(A(x)\)在\(w_n^0, w_n^1, ..., w_n^{n - 1}\)处的值，即为求\(A_1(x)\)和\(A_2(x)\)在\((w_n^0)^2, (w_n^1)^2, ..., (w_n^{n - 1})^2\)处的值。

根据性质1可得：\((w_n^k)^2 = w_{n/2}^k\)，因此\((w_n^0)^2, (w_n^1)^2, ..., (w_n^{n - 1})^2\)这n个单位复数根仅是由\(n / 2\)个不同的值组成的，我们可以先求出\(w_{n / 2}^0, ..., w_{n / 2}^{n / 2 - 1}\)这\(n / 2\)个值，再根据性质2，求出后一半的值，于是问题的规模缩小了一半，使用分治的方法即可在\(O(n \ log \ n)\)的时间内求值。

2.插值

将\(DFT\)的过程写成矩阵乘法的形式，则可得

\[y = V_n a \]

其中\(V_n\)为\(x = w_n\)时的范德蒙德矩阵。则：

\[a = V_n^{-1}y \]

根据范德蒙德矩阵的性质，易得\(V_n^{-1}\)的\((i, j)\)处的元素为\(w_n^{-ij} \ / \ n\) 则可知：

\[a_j = \frac{1}{n} \cdot \sum_{k = 0}^{n - 1} \ y_kw_n^{-kj} \]

式子与上面求值的\(DFT\)相似，通过类似方法也可以在\(O(n \ log \ n)\)的时间内插值。

详细证明

模板

代码

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
const int N = 1000005;
const double PI = acos(-1);
typedef complex <double> comp;
void FFT(int M, vector <comp> &a, int typ) {
    if(M == 1) return ;
    vector <comp> a1(M / 2 + 2), a2(M / 2 + 2);
    int cnt = 0;
    for(int i = 0; i <= M; i += 2) {
        a1[cnt] = a[i]; a2[cnt++] = a[i + 1];
    }
    FFT(M >> 1, a1, typ); FFT(M >> 1, a2, typ);
    comp wn(cos(PI * 2 / M), typ * sin(PI * 2 / M)), w(1, 0);
    for(int i = 0; i < (M >> 1); i++) {
        a[i] = a1[i] + w * a2[i];
        a[i + (M >> 1)] = a1[i] - w * a2[i];
        w = w * wn;
    }
}
int n, m, M;
int main() {
    // freopen("data.in", "r", stdin);
    cin >> n >> m;
    M = 1; while(M <= n + m) M <<= 1;
    vector <comp> a(M + 2), b(M + 2);
    for(int i = 0; i <= n; i++) {
        double x; cin >> x; 
        a[i] = complex <double> (x, 0);
    }
    for(int i = 0; i <= m; i++) {
        double x; cin >> x;
        b[i] = complex <double> (x, 0);
    }
    int M = 1; while(M <= n + m) M <<= 1;
    
    FFT(M, a, 1); FFT(M, b, 1);
    vector <comp> c(M + 2);
    for(int i = 0; i <= M; i++) {
        c[i] = a[i] * b[i];
    }
    FFT(M, c, -1);
    for(int i = 0; i <= n + m; i++) {
        cout << (int)(c[i].real() / M + 0.5) << " ";
    }
    return 0;
}