多项式全家桶

基础工具

快速傅立叶变换

引入

多项式：\(A(x) = \sum_{i = 0}^{n}\limits a_ix^i\)。

考虑多项式运算（设 \(A(x), B(x)\) 长度均为 \(n\)，因为长度不同时可以在后面补充若干系数为 \(0\) 的项）。

多项式加法： 求 \(C(x) = A(x) + B(x)\)。

容易发现 \(\forall i \in[0, n], c_i = a_i + b_i\)。

多项式减法：求 \(C(x) = A(x) - B(x)\)。

容易发现 \(\forall i \in[0, n], c_i = a_i - b_i\)。

多项式乘法：求 \(C(x) = A(x)⋅B(x)\)。

此时并不满足 \(\forall i \in [0, n], c_i = a_i\times b_i\)。

容易发现 \(\forall i \in [0, 2n], c_i = \sum_{j=0}^{min(i, n)}\limits a_jb_{i-j}\)。\((1)\)

这是因为所有满足 \(i + j = k\) 的 \(a_ix^i, b_jx^j\) 之积的 \(x\) 的次数都为 \(k\)，都对 \(c_kx^k\) 作出了贡献。

此时我们通过枚举 \((1)\) 中的 \(i\) 和 \(j\) 来计算 \(C(x)\) 的时间复杂度为 \(O(n^2)\)。

快速傅立叶变换就是一种可以在 \(O(n\log n)\) 复杂度内计算出 \(C(x)\) 的优异算法。

快速傅立叶变换 FFT

多项式的点值表示法：通过 \(n + 1\) 组形如 \((x_i, y_i)\) 的不同的点可以唯一确定一个 \(n\) 次多项式 \(A(x)\)，其中满足 \(\forall i\in[0, n], y_i = A(x_i)\)。

可以唯一确定的原因是我们可以将这个 \(n\) 次多项式的每个系数看作一个 \(n + 1\) 元方程的未知数，如此我们的 \(n + 1\) 个点就可以看成 \(n + 1\) 个方程式。此时无解的情况为：其中存在未知数系数线性相关的若干组方程式。但因为多项式的常数项没有乘 \(x\)，因此任意两个方程式线性相关时的比必须为 \(1\)，也就意味着 \(x_i\) 相等，与题设不符。

容易想到一个利用点值表示法的做法：任意取 \(2n + 1\) 个互不相同的数 \(x_i\)，分别代入两个多项式 \(O(n)\) 求出对应的 \(y_{i, A}, y_{i, B}\)，然后将所有的 \(y_{i, A}\times y_{i, B}\) 变为新的 \(y_i\)，再根据这些 \((x_i, y_i)\) 推导出 \(C(x)\)。

我们先不谈如何通过 \((x_i, y_i)\) 推导出 \(C(x)\)。这里前面分别代入的复杂度高达 \(O(n^2)\)，与暴力没有区别。

但这种方法给了我们很大的自由：代入的互不相同的 \(x_i\) 由我们选择。如果我们选择恰到好处的 \(x_i\)，那么通过特定的优化手法，是否可以将复杂度降低？答案是肯定的。

我们称 \(n\) 次单位根为满足 \(x^n = 1\) 的 \(x\)。

因为 \(n\) 次方程一定有 \(n\) 个解，所以在复数域上一定存在 \(n\) 个单位根。

我们设幅角为正且最小的向量对应的复数为 \(\omega_n\)。

那么这 \(n\) 个单位根分别是 \(\omega_n^0, \omega_n^1...\omega_n^{n-1}\)。

根据欧拉公式，我们知道：

\(\omega_n^k = \cos (k\times \dfrac{2\pi}{n}) +i\sin (k\times \dfrac{2\pi}{n})\)

接下来我们研究单位根的性质：

\(1.\) \(\omega_{2n}^{2k} = \omega_n^k\)

证：\(\omega_{2n}^{2k} = \cos (2k\times \dfrac{2\pi}{2n}) + i\sin (2k\times \dfrac{2\pi}{2n}) = \cos (k\times \dfrac{2\pi}{n}) + i\sin (k\times \dfrac{2\pi}{n})=\omega_n^k\)

\(2.\) \(\omega_{n}^{k+\frac{n}{2}} = -\omega_n^k\)

证：\(\omega_n^{k+\frac{n}{2}} = \omega_{n}^k \times \omega_n^{\frac{n}{2}}=\omega_n^k \times(\cos (\dfrac{n}{2}\times \dfrac{2\pi}{n})+i\sin (\dfrac{n}{2}\times \dfrac{2\pi}{n}))=\omega_n^k \times (-1)=-\omega_n^k\)

在这里，“恰到好处”的 \(x_i\) 的取值恰好是这 \(n\) 个单位根。

这里说明一下：为了方便计算，我们将 \(n\) 变成严格大于 \(2n + 1\) 的 \(2\) 的整次幂，因此我们只用将它看作 \(n - 1\) 次多项式处理即可。

考虑多项式：\(A(x) = \sum_{i = 0}^{n - 1}\limits a_ix^i\)。

将其系数按奇偶分组：\(A(x) = \sum_{i = 0}^{\frac{n}{2} - 1}\limits a_{2i}x^{2i} + \sum_{i = 0}^{\frac{n}{2} - 1}\limits a_{2i+1}x^{2i+1}\)。

将两组看成两个新的多项式：\(A_1(x) = \sum_{i = 0}^{\frac{n}{2} - 1}\limits a_{2i}x^{2i}, A_2(x) = \sum_{i = 0}^{\frac{n}{2} - 1}\limits a_{2i+1}x^{2i+1}\)。

由于两个多项式次数不同，将第二个多项式除以 \(x\)，此时 \(A_1(x) = \sum_{i = 0}^{\frac{n}{2} - 1}\limits a_{2i}x^{2i}, A_2(x) = \sum_{i = 0}^{\frac{n}{2} - 1}\limits a_{2i+1}x^{2i}\)。

此时 \(A(x) = A_1(x) + xA_2(x)\)。

注意到 \(A_1\) 和 \(A_2\) 都不是正常多项式的形式，于是将每一项的 \(x\) 都开根。

此时 \(A_1(x) = \sum_{i = 0}^{\frac{n}{2} - 1}\limits a_{2i}x^i, A_2(x) = \sum_{i = 0}^{\frac{n}{2} - 1}\limits a_{2i+1}x^i, A(x) = A_1(x^2) + xA_2(x^2)\)。

将 \(n\) 个单位根分两类代入：\(\omega_n^k(0 \le k \le \dfrac{n}{2} - 1), \omega_n^k (\dfrac{n}{2} \le k \le n-1)\)。

显然，第二类等价于 \(\omega_n^{k + \frac{n}{2}}(0\le k\le\dfrac{n}{2} - 1)\)。

此时，对于第一类：

\(A(\omega_n^k) = A_1(\omega_n^{2k}) + \omega_n^k A_2(\omega_n^{2k})\)

\(A(\omega_n^k) = A_1(\omega_{\frac{n}{2}}^k) + \omega_n^kA_2(\omega_{\frac{n}{2}}^k)\)

对于第二类：

\(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(\omega_n^{k+\frac{n}{2}}) = A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k})\)

\(A(\omega_n^{k + \frac{n}{2}}) = A_1(\omega_{\frac{n}{2}}^k) - \omega_n^kA_2(\omega_{\frac{n}{2}}^k)\)

至此，我们就通过利用单位根的奇妙性质将所有 \(A(\omega_n^k)\) 的值转化为 \(A_1(\omega_{\frac{n}{2}}^k)\) 和 \(A_2(\omega_{\frac{n}{2}}^k)\) 运算的结果。

容易发现，这时 \(A_1\) 与 \(A_2\) 的问题与 \(A\) 完全相同：我们需要的代入 \(A_1\) 和 \(A_2\) 的单位根恰好都是 \(\dfrac{n}{2}\) 次单位根，与 \(A_1\) 和 \(A_2\) 的长度相同。于是，我们分治 \(A_1\) 和 \(A_2\)，然后再根据它们的值算出 \(A\) 即可，时间复杂度为 \(O(n\log n)\)。

值得注意的是，我们在过程中利用 \(\omega_n^k = \omega_n^{k-1} \times \omega_n\) 计算 \(\omega_n^k\)。

快速傅立叶逆变换 IFFT

容易发现，我们刚才只是计算出了 \(A\) 在代入 \(n\) 个 \(n\) 次单位根时的点值表示。

我们用相同方法可以计算出 \(B\) 在代入 \(n\) 个 \(n\) 次单位根时的点值表示。

接下来，我们将这些对应的点值表示相乘得出 \(C(x)\) 在代入 \(n\) 个 \(n\) 次单位根时的点值表示。

我们的问题转化为：给出一个多项式在取 \(n\) 个单位根的值时的点值表示，求出这个多项式。

设 \(y_0, y_1, ..., y_{n-1}\) 为 \(C(x)\) 在代入 $n $ 单位根时的点值表示。

设 \(D(x) = \sum_{i = 0}^{n-1}\limits y_ix^i\)，\(z_0, z_1, ..., z_{n-1}\) 为 \(D(x)\) 在代入 \(\omega_n^0, \omega_n^{-1}, ..., \omega_n^{-n+1}\) 时的点值表示。

那么有：

\[z_k = \sum_{i=0}^{n-1}\limits y_i(\omega_n^{-k})^i\\ z_k = \sum_{i=0}^{n-1}\limits (\sum_{j=0}^{n-1} c_j(\omega_n^i)^j)(\omega_n^{-k})^i\\ z_k = \sum_{i=0}^{n-1}\limits (\sum_{j=0}^{n-1} c_j(\omega_n^j)^i)(\omega_n^{-k})^i\\ z_k = \sum_{i=0}^{n-1}\limits (\sum_{j=0}^{n-1} c_j(\omega_n^j)^i(\omega_n^{-k})^i)\\ z_k = \sum_{i=0}^{n-1}\limits (\sum_{j=0}^{n-1} c_j(\omega_n^{j-k})^i)\\ z_k = \sum_{i=0}^{n-1}\limits \sum_{j=0}^{n-1} c_j(\omega_n^{j-k})^i\\ z_k = \sum_{j=0}^{n-1}\limits \sum_{i=0}^{n-1} c_j(\omega_n^{j-k})^i\\ z_k = \sum_{j=0}^{n-1}\limits c_j (\sum_{i=0}^{n-1}(\omega_n^{j-k})^i) \]

设 \(E(x) = \sum_{i=0}^{n-1}\limits x^i\)。

显然 \(\omega_n^kE(\omega_n^k) - E(\omega_n^k) = (\omega_n^k)^n - 1\)。

那么 \(E(\omega_n^k)(\omega_n^k - 1) = (\omega_n^k)^n-1\)。\((2)\)

分类讨论：

当 \(k = 0\) 时，显然 \(\omega_n^0 = 1\)，所以 \(E(\omega_n^0) = n\)。

当 \(k \ne 0\) 时，因为 \(\omega_n^k \ne 1\)，所以将 \((2)\) 两侧同除 \((\omega_n^k-1)\)，得到：

\(E(\omega_n^k) = \dfrac{(\omega_{n}^k)^n-1}{\omega_n^k-1}\)

\(E(\omega_n^k) = \dfrac{(\omega_{n}^n)^k-1}{\omega_n^k-1}\)

\(E(\omega_n^k) = \dfrac{1-1}{\omega_n^k-1}=0\)

再回到 \(z_k = \sum_{j=0}^{n-1}\limits c_j (\sum_{i=0}^{n-1}\limits(\omega_n^{j-k})^i)\)。显然，此时 \(z_k = \sum_{j=0}^{n-1}\limits c_j E(\omega_n^{j-k})\)。

那么只有当 \(j\) 取到 \(k\) 时，\(E(\omega_n^{j-k})\ne 0\)。

此时 \(c_jE(\omega_n^{j-k}) = nc_k\)。

故 \(z_k = nc_k, c_k = \dfrac{z_k}{n}\)。

根据这个深刻的性质，我们可以想到一个做法：

先通过 \(FFT\) 将 \(D(x)\) 在代入 \(\omega_n^0, \omega_n^{-1}, ..., \omega_n^{-n+1}\) 时的点值表示得出，然后再将它们除以 \(n\) 得到 \(c_k\)。

这里有一个十分显然的性质：\(\omega_n^{-k} = \cos (k\times \dfrac{2\pi}{n}) - i\sin (k\times \dfrac{2\pi}{n})\)。

于是我们在 \(FFT\) 中额外传入一个参数 \(coe\)，表示虚部的系数。

此时我们利用 \(\omega_n^{-k} = \omega_n^{-k+1}\times \omega_n^{-1}\) 来计算 \(\omega_n^{-k}\)。

到这里，我们就用 \(O(n\log n)\) 的时间复杂度计算出了两个多项式的卷积。

代码

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <queue>
#include <map>
#include <set>
#include <cstdlib>
#include <cmath>
#include <deque>
using namespace std;
#define pii pair<int,int>
#define mp make_pair
const int N = 2.1e6 + 10;
const double PI = acos(-1.0);
int read()
{
    int x = 0, f = 1;
    char c = getchar();
    while (c < '0' || c > '9')
    {
        if (c == '-') f = -1;
        c = getchar();
    }
    while (c >= '0' && c <= '9')
    {
        x = (x << 1) + (x << 3) + (c ^ 48);
        c = getchar();
    }
    return x * f;
}
struct complex
{
    double x, y;
    complex () {x = 0; y = 0;}
    complex (double a, double b) {x = a; y = b;}
    complex operator + (const complex &a) {return complex(a.x + x, a.y + y);}
    complex operator - (const complex &a) {return complex(x - a.x, y - a.y);}
    complex operator * (const complex &a) {return complex(x * a.x - y * a.y, y * a.x + x * a.y);}
};
namespace tool
{
    void FFT(int len, complex * a, int coe)
    {
        if (len == 1) return ;
        complex a1[len >> 1], a2[len >> 1];
        for (int i = 0; i < len; i += 2)
        {
            a1[i >> 1] = a[i];
            a2[i >> 1] = a[i + 1];
        }
        FFT(len >> 1, a1, coe);
        FFT(len >> 1, a2, coe);
        complex Wn = complex(cos(2.0 * PI / len), coe * sin(2.0 * PI / len)), W = complex(1, 0);
        for (int i = 0; i < (len >> 1); i++, W = Wn * W)
        {
            a[i] = a1[i] + W * a2[i];
            a[i + (len >> 1)] = a1[i] - W * a2[i];
        }
    }
}
using namespace tool;
void mul(int len, complex * a, complex * b, complex * c)
{
    FFT(len, a, 1); FFT(len, b, 1);
    for (int i = 0; i <= len; i++)
        a[i] = a[i] * b[i];
    FFT(len, a, -1);
    for (int i = 0; i <= len; i++)
        c[i].x = int(a[i].x / len + 0.5);
}
complex a[N], b[N], res[N];
int main()
{
    int n, m, len;
    n = read(), m = read();
    for (int i = 0; i <= n; i++) a[i].x = read();
    for (int i = 0; i <= m; i++) b[i].x = read();
    if (n == m && m == 0)
    {
        printf("%d", (int)(a[0].x * b[0].x));
        return 0;
    }
    len = 1;
    while (len <= n + m) len <<= 1;
    mul(len, a, b, res);
    for (int i = 0; i <= n + m; i++) printf("%d ", (int)res[i].x);
    return 0;
}