浅谈算法——多项式乘法相关

从多项式乘法到FFT

这一段大部分是复制以前我写的这篇博客：

https://blog.csdn.net/wang3312362136/article/details/79510933

这篇博客有详细的代码，但是阅读体验布星，因此我把它的大部分证明过程复制下来。

如果会fft，请跳过这一段。

多项式乘法

我们知道，多项式可以表示成：

A=∑i=0naixi$$A=\sum _{i=0}^n{a}_i{x}^i$$

的形式。

对于两个多项式A(x)$A(x)$和B(x)$B(x)$，我们可以计算乘积A⋅B$A\cdot B$：

A⋅B=∑i=0|A|∑j=0|B|aibjxi+j$$A\cdot B=\sum _{i=0}^{|A|}\sum _{j=0}^{|B|}{a}_i{b}_j{x}^{i+j}$$

但是，这样算是O(|A|×|B|)$O(|A|\times |B|)$的，太慢了，怎么办？

我们需要换一条思路。

首先，我们得知道一个东西：多项式的点值表示法。

我们把上面的称为多项式的系数表示法，而点值表示法就是：

若A$A$多项式的次数为n$n$，则任取n$n$个不相同的x0,x1,⋯,xn$x_0,x_1,\cdots ,x_n$，求出A$A$多项式的A(x0),A(x1),⋯,A(xn)$A(x_0),A(x_1),\cdots ,A(x_n)$。记为：

<(x0,A(x0)),(x1,A(x1)),⋯,(xn,A(xn))>$$<(x_0,A(x_0)),(x_1,A(x_1)),\cdots ,(x_n,A(x_n))>$$

显然，一个有n+1$n+1$个点的点对唯一表示一个n$n$次多项式。

对于一个点值表示法下多项式

<(x0,A(x0)),(x1,A(x1)),⋯,(xn,A(xn))>$$<(x_0,A(x_0)),(x_1,A(x_1)),\cdots ,(x_n,A(x_n))>$$

和

<(x0,B(x0)),(x1,B(x1)),⋯,(xn,B(xn))>$$<(x_0,B(x_0)),(x_1,B(x_1)),\cdots ,(x_n,B(x_n))>$$

它们的乘积是

<(x0,A(x0)⋅B(x0)),(x1,A(x1)⋅B(x1)),⋯,(xn,A(xn)⋅B(xn))>$$<(x_0,A(x_0)\cdot B(x_0)),(x_1,A(x_1)\cdot B(x_1)),\cdots ,(x_n,A(x_n)\cdot B(x_n))>$$

可以看出点值表示法的多项式相乘是O(n)$O(n)$的。

等等，我们好像找到了一个突破口！

为啥不把原来的系数表示法下的多项式转化成点值表示法呢？

仔细想一想：系数表示法与点值表示法互相转换，这个步骤好像是O(n2)$O(n^2)$的。

而FFT（快速傅里叶变换）就是为了优化这个O(n2)$O(n^2)$。

PS：对于O(n2)$O(n^2)$的点值表示法转化成系数表示法可以看百度百科中关于插值法的介绍。

FFT（快速傅里叶变换）

如果未特别说明，那么下面的多项式次数将是2k−1$2^k-1$的形式。

如果不是关键部分的公式或定理，不提供证明，自己出门右转百度。

首先介绍两个概念：

DFT（离散傅里叶变换）是将多项式由系数表示法转化成点值表示法；

IDFT（离散傅里叶逆变换）是将多项式由点值表示法转化成系数表示法；

而FFT就是上述两种变换的优化。

DFT部分

前置技能：

下面的内容将会提到复数，不会的可以参考百度百科中关于复数的介绍；

为了介绍FFT中的DFT部分，首先要介绍的是一个概念：单位根。

单位根：若有

zn=1$$z^n=1$$

此时将z$z$称为n$n$次单位根。

若有z∈ℝ$z\in \mathbb{R}$，显然，z$z$可以等于1$1$，如果n$n$是偶数，则z$z$还可以等于−1$-1$。

我们把范围扩大到z∈ℂ$z\in \mathbb{C}$，那么，我们可以得到n$n$个复数，它们将复平面上的单位圆等分成n$n$份。

为了表示n$n$次单位根，我们引入一个公式。

欧拉公式：

ein=cosn+isinn$$e^{in}=\cos n+i\sin n$$

如果我们令：

ωn=e2πi/n$${\omega }_n=e^{2\pi i/n}$$

那么，n$n$次单位根就可以表示成ω0n,ω1n,⋯,ωn−1n${\omega }_n^0,{\omega }_n^1,\cdots ,{\omega }_n^{n-1}$，它们的n$n$次方显然都是1$1$。

下面是关于ωn${\omega }_n$的两条性质：（都是在n$n$为偶数的情况下）

ωn/2n=e(2πi/n)⋅(n/2)=eπi=cosπ+isinπ=−1(1)$$\begin{array}{}\text{(1)}& {\omega }_n^{n/2}=e^{(2\pi i/n)\cdot (n/2)}=e^{\pi i}=cos\pi +i\sin \pi =-1\end{array}$$

ω2n=e2⋅2πi/n=ωn/2(2)$$\begin{array}{}\text{(2)}& {\omega }_n^2=e^{2\cdot 2\pi i/n}={\omega }_{n/2}\end{array}$$

下面，我们进入正题：DFT的求法。

在这里，我们令多项式次数为n−1$n-1$，那么我们可以用点值表示成

<(ω0n,A(ω0n),(ω1n,A(ω1n)),⋯,(ωn−1n,A(ωn−1n))>$$<({\omega }_n^0,A({\omega }_n^0),({\omega }_n^1,A({\omega }_n^1)),\cdots ,({\omega }_n^{n-1},A({\omega }_n^{n-1}))>$$

额……这时间复杂度好像并没有减少……
别急，我们来看A(ωkn)$A({\omega }_n^k)$能够表示成什么。

A(ωkn)====∑i=0n−1aiωkin∑i=0n/2−1a2iω2kin+∑i=0n/2−1a2i+1ω2ki+kn∑i=0n/2−1a2iωkin/2+∑i=0n/2−1a2i+1ωkin/2ωkn∑i=0n/2−1a2iωkin/2+ωn∑i=0n/2−1a2i+1ωkin/2(3)(4)(5)(6)$$\begin{array}{}\text{(3)}& A({\omega }_n^k)=& \sum _{i=0}^{n-1}{a}_i{\omega }_n^{ki}\text{(4)}& =& \sum _{i=0}^{n/2-1}{a}_{2i}{\omega }_n^{2ki}+\sum _{i=0}^{n/2-1}{a}_{2i+1}{\omega }_n^{2ki+k}\text{(5)}& =& \sum _{i=0}^{n/2-1}{a}_{2i}{\omega }_{n/2}^{ki}+\sum _{i=0}^{n/2-1}{a}_{2i+1}{\omega }_{n/2}^{ki}{\omega }_n^k\text{(6)}& =& \sum _{i=0}^{n/2-1}{a}_{2i}{\omega }_{n/2}^{ki}+{\omega }_n\sum _{i=0}^{n/2-1}{a}_{2i+1}{\omega }_{n/2}^{ki}\end{array}$$

我们来分别分析以下这几个神奇的步骤。

(3)$(3)$步骤就是将ωkn${\omega }_n^k$带入原来的A$A$多项式。

(4)$(4)$步骤就是将原多项式拆成两个部分，按奇偶分类。

(5)$(5)$步骤用到了上面提到的性质(2)$(2)$。

(6)$(6)$步骤就是上面式子的后半部分提出公因数。

有了这个等式，我们就可以分治+递归解决DFT了。

算法步骤：

• 对当前的多项式（一个数组）系数进行奇偶分类；
• 递归算出偶数部分的数组的anse$an{s}_e$和奇数部分的数组的anso$an{s}_o$；
• 这个多项式的ans=anse+ωnanso$ans=an{s}_e+{\omega }_{n}an{s}_o$

但是这个的常数好像很大啊？能不能减少一点呢？

上面的性质(1)$(1)$给了我们提示：

ωn/2+kn=ωn/2n⋅ωkn=−ωkn$${\omega }_n^{n/2+k}={\omega }_n^{n/2}\cdot {\omega }_n^k=-{\omega }_n^k$$

在算k<n2$k<\frac{n}{2}$时，可以顺便把k≥n2$k\ge \frac{n}{2}$的情况也算出来。

常数减小了一半！但是还是很大啊！

递归版的程序一般比非递归版慢，为啥不用非递归版呢？

算法核心就是奇偶分类，分来分去最后分到了哪里？我们来研究研究。

显然，一个序列原来是0,1,2,3,4,5,6,7$0,1,2,3,4,5,6,7$，最终变成0,4,2,6,1,5,3,7$0,4,2,6,1,5,3,7$。

把它们的二进制列出来：

000,001,010,011,100,101,110,111000,100,010,110,001,101,011,111$$\begin{gathered} 000,001,010,011,100,101,110,111 \\ 000,100,010,110,001,101,011,111 \end{gathered}$$

其中，上面是位置，下面是这个位置对应的数。

把上面的数翻转，好像就是下面的数！

没错，只需要计算一下每个数的二进制翻转后的结果，就能得到一个数最终对应的位置，也就能实现非递归版了。

IDFT部分

回顾上面的DFT部分，仔细思考一下，它本质就是在求：

⎧⎩⎨⎪⎪⎪⎪a0(ω0n)0+a1(ω0n)1+⋯+an−1(ω0n)n−1=A(ω0n)a0(ω1n)0+a1(ω1n)1+⋯+an−1(ω1n)n−1=A(ω1n)⋯a0(ωn−1n)0+a1(ωn−1n)1+⋯+an−1(ωn−1n)n−1=A(ωn−1n)$$\{\begin{array}{l}{a}_0({\omega }_n^0{)}^0+a_1({\omega }_n^0{)}^1+\cdots +a_{n-1}({\omega }_n^0{)}^{n-1}=A({\omega }_n^0)\\ a_0({\omega }_n^1{)}^0+a_1({\omega }_n^1{)}^1+\cdots +a_{n-1}({\omega }_n^1{)}^{n-1}=A({\omega }_n^1)\\ \cdots \\ a_0({\omega }_n^{n-1}{)}^0+a_1({\omega }_n^{n-1}{)}^1+\cdots +a_{n-1}({\omega }_n^{n-1}{)}^{n-1}=A({\omega }_n^{n-1})\end{array}$$

其中，给定了a0,a1,⋯,an−1$a_0,a_1,\cdots ,a_{n-1}$以及ω0n,ω1n,⋯,ωn−1n${\omega }_n^0,{\omega }_n^1,\cdots ,{\omega }_n^{n-1}$，求A(ω0n),A(ω1n),⋯,A(ωn−1n)$A({\omega }_n^0),A({\omega }_n^1),\cdots ,A({\omega }_n^{n-1})$的值。

用矩阵表示如下：

⎡⎣⎢⎢⎢⎢(ω0n)0(ω1n)0⋮(ωn−1n)0(ω0n)1(ω1n)1⋮(ωn−1n)1⋯⋯⋱⋯(ω0n)n−1(ω1n)n−1⋮(ωn−1n)n−1⎤⎦⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢a0a1⋮an−1⎤⎦⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢A(ω0n)A(ω1n)⋮A(ωn−1n)⎤⎦⎥⎥⎥⎥(7)$$\begin{array}{}\text{(7)}& \left[\begin{array}{cccc}({\omega }_n^0{)}^0& ({\omega }_n^0{)}^1& \cdots & ({\omega }_n^0{)}^{n-1}\\ ({\omega }_n^1{)}^0& ({\omega }_n^1{)}^1& \cdots & ({\omega }_n^1{)}^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ ({\omega }_n^{n-1}{)}^0& ({\omega }_n^{n-1}{)}^1& \cdots & ({\omega }_n^{n-1}{)}^{n-1}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{n-1}\end{array}\right]=\left[\begin{array}{c}A({\omega }_n^0)\\ A({\omega }_n^1)\\ ⋮\\ A({\omega }_n^{n-1})\end{array}\right]\end{array}$$

我们令：

V=⎡⎣⎢⎢⎢⎢(ω0n)0(ω1n)0⋮(ωn−1n)0(ω0n)1(ω1n)1⋮(ωn−1n)1⋯⋯⋱⋯(ω0n)n−1(ω1n)n−1⋮(ωn−1n)n−1⎤⎦⎥⎥⎥⎥$$V=\left[\begin{array}{cccc}({\omega }_n^0{)}^0& ({\omega }_n^0{)}^1& \cdots & ({\omega }_n^0{)}^{n-1}\\ ({\omega }_n^1{)}^0& ({\omega }_n^1{)}^1& \cdots & ({\omega }_n^1{)}^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ ({\omega }_n^{n-1}{)}^0& ({\omega }_n^{n-1}{)}^1& \cdots & ({\omega }_n^{n-1}{)}^{n-1}\end{array}\right]$$

那么IDFT的本质就是求V$V$矩阵的逆矩阵。

普通的矩阵求逆显然无法快速求出，这必须依靠人类智慧大力猜想逆矩阵。

人类智慧还真求出来了，考虑下面这个矩阵：

D=⎡⎣⎢⎢⎢⎢⎢(ω−0n)0(ω−1n)0⋮(ω−(n−1)n)0(ω−0n)1(ω−1n)1⋮(ω−(n−1)n)1⋯⋯⋱⋯(ω−0n)n−1(ω−1n)n−1⋮(ω−(n−1)n)n−1⎤⎦⎥⎥⎥⎥⎥$$D=\left[\begin{array}{cccc}({\omega }_n^{-0}{)}^0& ({\omega }_n^{-0}{)}^1& \cdots & ({\omega }_n^{-0}{)}^{n-1}\\ ({\omega }_n^{-1}{)}^0& ({\omega }_n^{-1}{)}^1& \cdots & ({\omega }_n^{-1}{)}^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ ({\omega }_n^{-(n-1)}{)}^0& ({\omega }_n^{-(n-1)}{)}^1& \cdots & ({\omega }_n^{-(n-1)}{)}^{n-1}\end{array}\right]$$

那么我们令E=D⋅V$E=D\cdot V$，则：

Ei,j===∑k=0n−1Di,kVk,j∑k=0n−1(ω−in)k(ωkn)j∑k=0n−1ωk(j−i)n$$\begin{array}{rl}{E}_{i,j}=& \sum _{k=0}^{n-1}{D}_{i,k}{V}_{k,j}\\ =& \sum _{k=0}^{n-1}({\omega }_n^{-i}{)}^k({\omega }_n^k{)}^j\\ =& \sum _{k=0}^{n-1}{\omega }_n^{k(j-i)}\end{array}$$

显然：

Ei,j={0(i≠j)n(i=j)$$E_{i,j}=\{\begin{array}{l}0(i\ne j)\\ n(i=j)\end{array}$$

因此

1nD⋅V=1nE=In$$\frac{1}{n}D\cdot V=\frac{1}{n}E=I_n$$

(7)$(7)$式两边同时左乘1nD$\frac{1}{n}D$，可得

⎡⎣⎢⎢⎢⎢a0a1⋮an−1⎤⎦⎥⎥⎥⎥=1n⎡⎣⎢⎢⎢⎢⎢(ω−0n)0(ω−1n)0⋮(ω−(n−1)n)0(ω−0n)1(ω−1n)1⋮(ω−(n−1)n)1⋯⋯⋱⋯(ω−0n)n−1(ω−1n)n−1⋮(ω−(n−1)n)n−1⎤⎦⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢A(ω0n)A(ω1n)⋮A(ωn−1n)⎤⎦⎥⎥⎥⎥$$\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{n-1}\end{array}\right]=\frac{1}{n}\left[\begin{array}{cccc}({\omega }_n^{-0}{)}^0& ({\omega }_n^{-0}{)}^1& \cdots & ({\omega }_n^{-0}{)}^{n-1}\\ ({\omega }_n^{-1}{)}^0& ({\omega }_n^{-1}{)}^1& \cdots & ({\omega }_n^{-1}{)}^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ ({\omega }_n^{-(n-1)}{)}^0& ({\omega }_n^{-(n-1)}{)}^1& \cdots & ({\omega }_n^{-(n-1)}{)}^{n-1}\end{array}\right]\left[\begin{array}{c}A({\omega }_n^0)\\ A({\omega }_n^1)\\ ⋮\\ A({\omega }_n^{n-1})\end{array}\right]$$

这就相当于把DFT中ωkn${\omega }_n^k$都换成ω−kn${\omega }_n^{-k}$。

IDFT过程中的小技巧

如果IDFT中逆矩阵使用的是V$V$（使用的是DFT的过程），要得到IDFT的正确结果只需要执行一遍std::reverse(a+1,a+n)，这是很容易证明的。

NTT（快速数论变换）

FFT要求采用复数和实数运算，如果数很大很容易产生精度误差。如果是模意义下呢？

前置技能

FFT，初级数论内容，原根。

一点微小的变化？

考虑原根的性质。

性质(1)$(1)$ 循环：gn+xn=gxn(modn)$g_n^{n+x}=g_n^x(\mathrm{mod}n)$

性质(2)$(2)$ 折半：gxn=gx/2n/2(modm)$g_n^x=g_{n/2}^{x/2}(\mathrm{mod}m)$

性质(3)$(3)$ ∑n−1i=0gin=0(modm)$\sum _{i=0}^{n-1}{g}_n^i=0(\mathrm{mod}m)$

可以发现原根的性质和单位根的性质差不多，因此可以用来做模意义下的FFT，就是NTT。

或者说，NTT就是把FFT的ω$\omega$替换成了g$g$。

由于使用的是gm−1i$g_i^{m-1}$，其中m$m$为模数，i$i$为2k$2^k$形式，因此模数是x⋅2k+1$x\cdot 2^k+1$的形式（NTT模数）。

常见的NTT模数：

998244353，原根是3

1004535809，原根也是3

FWT（快速沃尔什变换）

上面求的都是多项式卷积

F(n)=∑i=0na(i)b(n−i)$$F(n)=\sum _{i=0}^{n}a(i)b(n-i)$$

或者换种方式

F(n)=∑i,j,i+j=na(i)b(j)$$F(n)=\sum _{i,j,i+j=n}a(i)b(j)$$

但是如果是下面的方式

F(n)=∑i,j,i∘j=na(i)b(j)$$F(n)=\sum _{i,j,i\circ j=n}a(i)b(j)$$

其中∘$\circ$是种位运算呢？

依然考虑转化成点值表达式再相乘。类比FFT，考虑用一种变换tf$tf$和逆变换utf$utf$，只要tf(A×B)=tf(A)×tf(B)$tf(A\times B)=tf(A)\times tf(B)$，那么就可以得到A×B=utf(tf(A)×tf(B))$A\times B=utf(tf(A)\times tf(B))$

下面我们分类讨论∘$\circ$的位运算种类及tf,utf$tf,utf$的种类。

集合交卷积（and）

集合并卷积（or）

集合对称差卷积（xor）

子集卷积

其他的方式

多项式求逆

前置技能

多项式乘法（FFT/NTT）

多项式的度

一个多项式的度就是这个多项式最高此项的次数，例如4x2+3x+1$4x^2+3x+1$的度为2$2$。

多项式的逆元

现在有两个多项式F,G$F,G$，如果

F×G=1(modxn)(1)$$\begin{array}{}\text{(1)}& F\times G=1(\mathrm{mod}{x}^n)\end{array}$$

并且G$G$的度小于等于F$F$的度，那么G$G$称为F$F$在模xn$x^n$意义下的逆元。

怎么求？假设已经求出了F$F$在模x⌈n/2⌉$x^{\lceil n/2\rceil }$下的逆元，为G′$G^{\prime }$。

由定义得到

F×G′=1(modx⌈n/2⌉)(2)$$\begin{array}{}\text{(2)}& F\times G^{\prime }=1(\mathrm{mod}{x}^{\lceil n/2\rceil })\end{array}$$

将(1)$(1)$放到模x⌈n/2⌉$x^{\lceil n/2\rceil }$意义下

F×G=1(modx⌈n/2⌉)(3)$$\begin{array}{}\text{(3)}& F\times G=1(\mathrm{mod}{x}^{\lceil n/2\rceil })\end{array}$$

(2)−(3)$(2)-(3)$得到

F×G′−F×G=0(modx⌈n/2⌉)$$F\times G^{\prime }-F\times G=0(\mathrm{mod}{x}^{\lceil n/2\rceil })$$

两边平方，得到

F2×G2−2×F2×G′×G+F2×G′2=0(modxn)$$F^2\times G^2-2\times F^2\times G^{\prime }\times G+F^2\times G^{\prime 2}=0(\mathrm{mod}{x}^n)$$

同时除以F$F$并处理modxn$mod{x}^n$意义下为1$1$的项

G=2G′−FG′2(modxn)$$G=2G^{\prime }-F{G}^{\prime 2}(\mathrm{mod}{x}^n)$$

就是说，只要求出了模x⌈n/2⌉$x^{\lceil n/2\rceil }$意义下的逆元，就可以求出模xn$x^n$意义下的逆元。