多项式算法1：多项式乘法

多项式简介

对于数域$\mathbb{F}$，若有∀i∈{1,2,3,⋯ ,n}\forall i\in\{1,2,3,\cdots,n \}∀i∈{1,2,3,⋯,n}，则
$$f(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_{n}x=\sum _{i=1}^n{a}_i{x}^i$$

为数域$\mathbb{F}$上的一个多项式.

多项式的度数

• 多项式最高项的次数称作多项式的次数

• $f(x)=\sum _{i=1}^n{a}_i{x}^i$，其中$a_n̸=0$

• 记作$\text{deg}f(x)=n$

多项式的系数表示

对于上述多项式$f(x)$，有系数向量$\vec{a}$与之映射
$$f(x)=\sum _{i=1}^n{a}_i{x}^i\iff \vec{a}=(a_1,a_2,a_3,\cdots ,a_n)$$

我们称$\vec{a}$为$f(x)$的系数表示

多项式卷积

首先定义
$$f(x)=\sum _{i=1}^n{a}_i{x}^i$$

$$g(x)=\sum _{i=1}^n{b}_i{x}^i$$

于是他们的卷积形式为
$$h(x)=(f\cdot g)(x)=\left(\sum _{i=1}^n{a}_i{x}^i\right)\left(\sum _{i=1}^n{b}_i{x}^i\right)$$

可得$\text{deg}h(x)=2n$，于是令
$$h(x)=\sum _{i=1}^n{c}_i{x}^i$$

$c_i$满足以下运算律
$$c_i=\sum _{j+k=i}{a}_j{b}_k=\sum _{j=i}^n{a}_j{b}_{i-j}$$

系数向量$c$为$a$与$b$的卷积，记作
$$c=a\ast b$$

多项式点值表示

设$\text{deg}f(x)=n$，$\forall x_i\in \mathbb{F}(i=(1,2,3,\cdots ,n))$
$$y_i=f(x_i)$$

则我们可以在数域$\mathbb{F}$的坐标系中绘制出$n$个点
$$(x_1,y_1),(x_2,y_2),(x_3,y_3),\cdots ,(x_n,y_n)$$
同样与$f(x)$函数有映射关系。

不妨将这$n$个点带入矩阵运算
$$\left[\begin{array}{c}1x_1{x}_1^2\cdots x_1^n\\ \\ 1x_2{x}_2^2\cdots x_2^n\\ \\ 1x_3{x}_3^2\cdots x_3^n\\ \\ ⋮⋮⋮\ddots ⋮\\ \\ 1x_n{x}_n^2\cdots x_n^n\end{array}\right]\left[\begin{array}{c}{a}_1\\ \\ a_2\\ \\ a_3\\ \\ ⋮\\ \\ a_n\end{array}\right]=\left[\begin{array}{c}{y}_1\\ \\ y_2\\ \\ y_3\\ \\ ⋮\\ \\ y_n\end{array}\right]$$

我们称最左边的矩阵为范德蒙矩阵.

多项式插值

对于多项式点值表示$(x_1,y_1),(x_2,y_2),(x_3,y_3),\cdots ,(x_n,y_n)$

求它的系数向量表达式$\vec{a}=(a_1,a_2,a_3,\cdots ,a_n)$

这个过程叫做插值.

朴素的插值法可以利用范德蒙矩阵进行高斯消元，但是时间复杂度是无法承受的.

我们可以用拉格朗日插值求解
$$f(x)=\sum _{i=1}^n{y}_i\prod _{i̸=j}\frac{x-x_j}{x_i-x_j}$$

复杂度优化到了$n^2$

点值求解多项式乘法

1. 选取$n$个$x_i$，带入计算点值，这一步是$O(n^2)$

2. 计算点值的卷积，这一步是$O(n)$

3. 插值计算系数向量，这一步是$O(n^2)$

所以最后的时间复杂度应该为$O(n^2)$，朴素的点值求解与暴力计算一致。

接下来我们会学习以下知识点优化多项式卷积

• 复数单位根

• 原根

复数单位根

考虑方程
$$z^n=1(z\in \mathbb{C})$$

我们按照极角排序，命第$k$个复数解为${\omega }_n^k$

由初等三角函数知识，我们解得
$${\omega }_n^k=\cos \frac{k}{n}2\pi +i\sin \frac{k}{n}2\pi (k\in \mathbb{Z})$$

由欧拉公式可得
$${\omega }_n^k=e^{\frac{2\pi ki}{n}}$$

可以$O(1)$求解。

由三角函数的周期性$(\tau =2\pi )$，不难考虑${\omega }_n^k$的周期性

• ${\omega }_n^k={\omega }_n^{k(\mathrm{m}\mathrm{o}\mathrm{d}n)}$

• ${\omega }_n^n={\omega }_n^0=1$

我们取k∈{0,2,3,⋯ ,n−1}k\in\{0,2,3,\cdots,n-1\}k∈{0,2,3,⋯,n−1}

有以下定理

相消定理

$${\omega }_{dn}^{dk}={\omega }_n^k$$

充分证明
$${\omega }_{dn}^{dk}=\cos \frac{dk}{dn}2\pi +i\sin \frac{dk}{dn}2\pi$$

$$=\cos \frac{k}{n}2\pi +i\sin \frac{k}{n}2\pi ={\omega }_n^k$$
证毕.

此外补充几种常见形式
$${\omega }_{2n}^n={\omega }_2^1=-1$$

$${\omega }_{2n}^{n+k}=-{\omega }_{2n}^k$$

折半定理

$n$次单位根的平方的集合等于$\frac{n}{2}$次单位根的集合
A={(ωnk)2∣k∈Z}A=\{\left(\omega_n^k\right)^2|k\in \mathbb Z\}A={(ωnk​)2∣k∈Z}

B={ωn2k∣k∈Z}B=\{\omega_{\frac{n}{2}}^k|k\in\mathbb Z\}B={ω2n​k​∣k∈Z}

$$A\subseteq B,B\subseteq A\iff A=B$$

$${\left({\omega }_n^{k+\frac{n}{2}}\right)}^2={\omega }_n^{n+2k}={\omega }_n^{2k}={\left({\omega }_n^k\right)}^2$$

求和引理

$$\sum _{i=0}^{n-1}{\left({\omega }_n^k\right)}^i=0(k̸=n)$$

由等比数列的求和公式可得
$$\sum _{i=0}^{n-1}{\left({\omega }_n^k\right)}^i=\frac{1-{\left({\omega }_n^k\right)}^n}{1-{\omega }_n^k}=\frac{1-1}{1-{\omega }_n^k}=0$$

对于$(k=n)$
$$\sum _{k=0}^{n-1}{\left({\omega }_n^k\right)}^i=n$$

单位根的带入可以优化多项式卷积.

注：以下出现的$n$均为$2$的整数幂.

离散傅里叶变换

对于多项式$f(x)=(a_0,a_2,a_3,\cdots ,a_{n-1})$

求解$\forall x_i={\omega }_n^k(i=(0,1,2,\cdots ,n-1))$处的点值.

可以令
$$y_k=f(x_k)=f({\omega }_n^k)$$

此过程读作长度为$n$的离散傅里叶变换，记作
$$y_k={\text{DFT}}_n(f)$$

快速傅里叶变换

考虑分治
$$f_0(x)=(a_0,a_2,a_4,\cdots ,a_{n-2})$$

$$f_1(x)=(a_1,a_3,a_5,\cdots ,a_{n-1})$$
于是有
$$f(x)=f_0(x^2)+x{f}_1(x^2)$$

对于$k<\frac{n}{2}$，我们有
$$f({\omega }_n^k)=f_0({\omega }_n^{2k})+{\omega }_n^k{f}_1({\omega }_n^{2k})=f_0({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k{f}_1({\omega }_{\frac{n}{2}}^k)$$

对于$k>\frac{n}{2}$，替换成$k+\frac{n}{2}$
$$f\left({\omega }_n^{k+\frac{n}{2}}\right)=f_0\left({\omega }_n^{2k+n}\right)+{\omega }_n^{k+\frac{n}{2}}{f}_1({\omega }_n^{2k+n})=f_0({\omega }_{\frac{n}{2}}^k)-{\omega }_n^k{f}_1({\omega }_{\frac{n}{2}}^k)$$

令
$$y_k=f({\omega }_n^k)y_k^{[0]}=f_0({\omega }_{\frac{n}{2}}^k)y_k^{[1]}=f_1({\omega }_{\frac{n}{2}}^k)$$

由上述结论，我们有
$$y_k=y_k^{[0]}+{\omega }_n^k{y}_k^{[1]}$$

$$y_{k+\frac{n}{2}}=y_k^{[0]}-{\omega }_n^k{y}_k^{[1]}$$
上述操作为蝴蝶操作，其中${\omega }_n^k$为旋转因子.

快速傅里叶变换时间复杂度分析

令$T(n)$为长度为$n$的快速傅里叶变换的时间复杂度

递归求得$y_k^{[0]},y_k^{[1]}$,规模均为$\frac{n}{2}$

即
$$T(n)=2T\left(\frac{n}{2}\right)+O(n)=O(n{\log }_{2}n)$$

快速傅里叶逆变换

设有向量$(c_0,c_1,c_2,\cdots ,c_{n-1})$，满足
$$c_k=\sum _{i=0}^{n-1}{y}_i({\omega }_n^{-k}{)}^i$$
展开$y_i$，我们有
$$c_k=\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j({\omega }_n^i{)}^j(w_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(w_n^{j-k}{)}^i$$

$$=\sum _{j=0}^{n-1}{a}_j\sum _{i=0}^{n-1}(w_n^{j-k}{)}^i$$

考虑到求和定理
$$\sum _{i=1}^{n-1}({\omega }_n^k{)}^i=0(k̸=n)$$

取$j=k$
$$c_k=n{y}_k$$

求解以下$n$的逆元就可以了
$$y_k=\frac{1}{n}{c}_k$$

于是快速傅里叶逆变换只需要取共轭复数就行了。

迭代求解

我们需要这样求解

对于$n=2^k$，考虑第$i$次分治

二进制第$i$位为$0$，放入左侧，最终位置的二进制的$k-i$位为$0$

二进制第$i$位为$1$，放入右侧，最终位置的二进制的$k-i$位为$1$

知识点讲完了，代码如下

#include <cstdio>
#include <complex>

using namespace std;
const int N = 1e+7 + 1;
const double Pi = acos(-1);
int n, m, rev[N];
complex<double> F[N], G[N], H[N];

int invert(int n) {
    int bit = 1;
    while((1 << bit) < n) bit++;
    return (1 << bit);
}
int getint() {
    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 - '0', c = getchar();
    return x * f;
}
void FFT(complex<double> *A, int n, int inv) {
    int bit = 1;
    while((1 << bit) < n) bit++;
    for(int i = 0; i < n; i++) {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
        if(i < rev[i]) swap(A[i], A[rev[i]]);
    }
    for(int mid = 1; mid < n; mid <<= 1) {
        complex<double> temp(cos(Pi / mid), inv * sin(Pi / mid));
        for(int i = 0; i < n; i += mid << 1) {
            complex<double> omega(1, 0);
            for(int j = 0; j < mid; j++, omega *= temp) {
                complex<double> x = A[i + j], y = omega * A[i + j + mid];
                A[i + j] = x + y;
                A[i + j + mid] = x - y;
            }
        }
    }
}
int main()
{
    scanf("%d %d", &n, &m);
    for(int i = 0; i <= n; i++) F[i].real(getint());
    for(int i = 0; i <= m; i++) G[i].real(getint());
    FFT(F, invert(n + m), 1);
    FFT(G, invert(n + m), 1);
    for(int i = 0; i <= invert(n + m); i++) H[i] = F[i] * G[i];
    FFT(H, invert(n + m), -1);
    for(int i = 0; i <= n + m; i++)
        printf("%d ", (int)(H[i].real() / invert(n + m) + 0.5));
    return 0;
}