多项式算法

Published: July 09, 2019 by Daltao

• Categories:
• algorithm 101

• Tags:

• 快速数论变换
• 多项式求逆
• 多项式除法取模
• 多项式多点求值
• 等差卷积
• 求序列等幂和
• 多项式差与多项式前缀和
• Chirp-Z Transform
• 多元多项式卷积
• 参考资料

快速数论变换

了解快速傅里叶变换的人应该知道，快速傅里叶变换之所以能成立源于对任意正整数n，存在这样的一个单位根$w(n)$，满足四个条件：

1. $w(n)^0,w(n)^1, \ldots, w(n)^{n-1}$各不相同
2. $w(n)^n = 1$
3. $w(dn)^{dk}=w(n)^k$。
4. $w(2n)^n=-1$

在快速傅里叶变换中，我们选取$w(n)=e^{(2\pi /n)i}$作为单位根。

快速傅里叶变换会存在一定的精度损失问题，我们可以使用快速数论变换来替代。但是快速数论变换实际上还是沿用了快速傅里叶变换的整体框架，只是替换了单位根。我们需要找一个素数p，并找到它的一个原根g。可以发现当我们取$w(n)=g^{\frac{p-1}{n}}$时，可以满足上面所有条件。不是所有p都适合用于快速数论变换，我们需要找一个比较特殊的p，使得$n|p-1$。由于n一般是2的幂次，因此我们希望$p-1$包含尽可能多的因子2，下面介绍两个常用的选择：

$$998244353=2^{23}\cdot17 \cdot 7 +1\\ 1004535809=2^{21}\cdot 479+1$$

剩下的和快速傅里叶变换没有什么区别，直接套用就好了。

多项式求逆

给出多项式A，找到A的逆多项式B，满足$AB=1\pmod{x^n}$。

当n为1时，问题非常简单，此时A的逆多项式一定是$1/A[0]$。假设我们已经找到了某个多项式C，使得$AC=1\pmod{x^{n/2}}$。那么，可以推导出：

$$AC=AB=1\pmod{x^{n/2}}\\ \Rightarrow A(B-C)=0\pmod{x^{n/2}}\\ \Rightarrow (B-C)=0\pmod{x^{n/2}}\\ \Rightarrow (B-C)^2=0 \pmod{x^n}$$

最后一步需要说明下，对于$B-C$的后n项都是0，因此对于$B-C$的$x^i$的系数，若$i<n/2$，那么系数一定是$B-C$的后$n/2$项中两项的乘积的累加，一定是0，若$i\geq n/2$，那么一定是后$n/2$项和前$n/2$项的系数乘积的累和，也是0。

将二次方程展开：

$$(B-C)^2=B^2+C^2-2BC=0\pmod{x^n}\\ \Rightarrow B^2=2BC-C^2\pmod{x^n}$$

之后两端都乘上A，借助$AB=1$可知：

$$AB^2=2ABC-AC^2\pmod{x^n}\\ \Rightarrow B=2C-AC^2\pmod{x^n}$$

这样我们就得到了多项式B，并且可以发现，整个流程我们仅假设$A[0]$可逆，即一个多项式可逆当且仅当其$x^0$项的系数可逆。

多项式求逆的时间复杂度为$T(n)=T(\frac{n}{2})+\alpha n\log_2n=O(2\alpha n\log_2n)=O(n\log_2n)$，其中$\alpha$是固定常数。

多项式除法取模

对于阶数为n的多项式$f(x)$，记$f^R(x)=x^nf(\frac{1}{x})$，换句话说$f^R(x)$是$f(x)$的系数前后交换后得到的。

接下来考虑阶数为$n$的多项式$A$和阶数为$m$的多项式$B$，以及阶数小于$m$的多项式$R$，我们希望得到如下形式等式$A(x)=B(x)C(x)+R(x)$。推一下公式：

$$\begin{aligned} A(x)&=B(x)C(x)+R(x)\\ \Rightarrow A(\frac{1}{x})&=B(\frac{1}{x})C(\frac{1}{x})+R(\frac{1}{x})\\ \Rightarrow x^nA(\frac{1}{x})&=(x^mB(\frac{1}{x}))(x^{n-m}C(\frac{1}{x}))+x^nR(\frac{1}{x})\\ \Rightarrow A^R(x)&=B^R(x)C^R(x)+x^{n-m+1}R^R(x)\\ \Rightarrow A^R(x)&=B^R(x)C^R(x)\pmod{x^{n-m+1}}\\ \Rightarrow C^R(x)&=A^R(x)(B^R(x))^{-1}\pmod{x^{n-m+1}} \end{aligned}$$

最后的结果可以通过对多项式求逆得到。

而取模可以通过$R(x)=A(x)-B(x)C(x)$得到。

多项式除法和取模都仅涉及到一次的多项式求逆和乘法，因此时间复杂度为$O(n\log_2n)$。

多项式多点求值

给出$n$阶多项式$f$，以及$m$个不同的数$x_1,\ldots, x_m$，要求计算$f(x_1),\ldots,f(x_m)$。

很显然可以设计一套暴力算法，$O(nm)$时间复杂度解决。但是这个问题存在丧心病狂的$O(n\log_2n+m(\log_2m)^2)$的解法。

我们可以$m$个数值分成相同大小的两部分。第一部分是$X_1={x_1,\ldots,x_{\frac{m}{2}}}$，第二部分是$X_2={x_{\frac{m}{2}+1},\ldots, x_m}$。

之后我们定义$g(x)=\prod_{i=1}^{\frac{m}{2}} (x-x_i)$，很显然对于任意$x'\in X_1$，满足$g(x')=0$。我们利用多项式除法可以得到：$f(x)=g(x)A(x)+R(x)$，其中$R$的阶数小于$\frac{m}{2}$。利用$x'\in X_1 \Rightarrow f(x')=R(x')$，因此我们现在只需要计算$X_1$在某个不超过$m$阶的多项式下的取值。重复这个过程，直到点集不可再分。

分治的时间复杂度为：

$$T(m)=2T(\frac{m}{2})+\alpha m\log_2m = \alpha m(\log_2m)^2=O(m(\log_2m)^2)。$$

总的时间复杂度为$O(n\log_2n+m(\log_2m)^2)$。

提供一道题目。

等差卷积

考虑计算下面多项式：

$$H(x)=\sum_{i=0}^nx^i\sum_{j=i}^nf(j)g(j-i)$$

它非常类似于一般卷积，只是内部是等差，而一般卷积是等和。我们定义$f'(x)=f(n-x)$。那么代入得到

$$H'(x)=\sum_{i=0}^nx^i\sum_{j=0}^if'(j)g(i-j)$$

很显然$H'$可以利用一般的快速卷积在$O(n\log_2n)$时间复杂度内得到，下面我们来看$H'$和$H$的关系。

$H'_i=\sum_{j=0}^if'(j)g(i-j)=\sum_{j=0}^if(n-j)g(i-j)$ 这里包含了所有差为$n-i$的项的和。因此可以得到：

$$H_i'=H_{n-i}$$

即最终我们可以通过翻转$H'$得到$H$。

注意这里我们要在$f$前面填充0之前翻转，且$H'$也仅翻转后面的$n$项目。

求序列等幂和

假设提供一个序列${a_i}$，记$f(k)=\sum_{i=1}^na_i^k$，要求我们计算出$f(1),f(2),\ldots, f(m)$。

这个问题可以用生成函数$G(x)$来求解。

$$G(x)=\sum_{i=0}\sum_{j=1}^na_{j}^ix^i\\ =\sum_{j=1}^n\sum_{i=0}(a_{j}x)^i\\ =\sum_{j=1}^n\frac{1}{1-a_jx}\\ =\frac{\sum_{j=1}^n\prod_{j\neq i}(1-a_jx)}{\prod_{j}(1-a_jx)}$$

假如我们记$D(x)=\prod_{j}(1-a_jx)$，记$U(x)=\sum_{j=1}^n\prod_{j\neq i}(1-a_jx)$，那么我们可以发现$D_i$和$U_i$存在某种联系。事实上，考虑所有对$x^j$有贡献的系数一定是$j$个${-a_i}$序列中的数的乘积，而这个乘积会在$U_i$中出现贡献$n-j$次。因此我们可以直接断定$U_i=(n-i)D_i$。

于是乎我们得到了$U(x)$和$D(x)$后，就可以用多项式求逆的技术得出$G(x)=U(x)D^{-1}(x)$。

这里我们计算$D(x)$的时候，可以利用分治+快速卷积的技术将时间复杂度优化到$O(n(\log_2n)^2 )$。

多项式差与多项式前缀和

定义1：对于一个任意多项式$f(x)$，我们记其前缀和函数为$g(x)=\sum_{i=0}^xf(i)$。

命题1：对于任意$n$阶多项式$A(x)$，我们始终可以找到另外一个$n+1$阶多项式$B(x)$，满足$B(x)-B(x-1)=A(x)$

证明：

我们可以利用归纳法进行证明，当$A(x)$是0阶多项式的时候，即$A(x)=c$，那么我们可以得到$B(x)=cx$，很显然$B(x)-B(x-1)=cx-c(x-1)=c=A(x)$。

假设对于任意不大于$k-1$阶的多项式$A(x)$，命题都成立，那么当$A(x)$为$k$阶多项式的时候，记$A(x)=a_kx^k+\ldots+a_0$，注意到$x^k-(x-1)^k$是$k-1$阶多项式，因此我们可以通过为$B(x)$选择合适的最高项系数$b_{k+1}=\frac{a_k}{k}$，将问题化简为$(B(x)-b_{k+1}x^{k+1})-(B(x-1)-b_{k+1}(x-1)^{k+1})=A(x)-(b_{k+1}x^{k+1}-b_{k+1}(x-1)^{k+1})$，我们重新标记后得到公式：$B'(x)-B'(x-1)=A'(x)$，其中我们要找的$B'(x)$是$k$阶多项式，而$A'(x)$是$k-1$阶多项式，而之前的归纳法已经证明了$B'(x)$的存在性，因此$B(x)=B'(x)+b_{k+1}x^{k+1}$也是存在的。

命题2：如果多项式$f(x)$是$n$阶多项式，那么其前缀和函数$g(x)$一定是$n+1$阶多项式

证明：

首先注意到$g(x)-g(x-1)=f(x)$，按照命题1我们一定可以刻画出这样一个$n+1$阶多项式$g(x)$，现在只需要保证$g(0)=f(0)$即可。注意到修改$g(x)$的常数项不会影响$g(x)-g(x-1)=f(x)$等式的成立，因此我们可以通过修改$g(x)$的常数项为$f(0)$即可找到与前缀和函数完全相同的一个$n+1$阶多项式。

Chirp-Z Transform

考虑这样一个问题：给定一个$n$阶多项式$f(x)=\sum_{i=0}^na_ix^i$，以及某个数$c$和$m$，要求计算$f(c^0),f(c^1),\ldots,f(c^m)$。其中$1\leq n, m\leq 3\times 10^5$。

暴力做法$n^2$显然不可行。

可以发现

$$f(c^j)=\sum_{i=0}^na_ic^{ij}$$

由于$ij={i+j\choose 2}-{i\choose 2}-{j\choose 2}$，因此有

$$\begin{aligned} f(c^j)&=\sum_{i=0}^na_ic^{ij}\\ &=\sum_{i=0}^na_ic^{ {i+j\choose 2}-{i\choose 2}-{j\choose 2} }\\ \end{aligned}$$

稍加转换可以得到：

$$\begin{aligned} &f(c^j)c^{j\choose 2}=\sum_{i=0}^na_ic^{ {i+j\choose 2}-{i\choose 2} S}\\ \Rightarrow &f(c^j)c^{j\choose 2}=\sum_{i=0}^n[c^{i+j\choose 2}][a_ic^{-{i\choose 2}}] \end{aligned}$$

可以发现右边是两个多项式的等差卷积，因此可以直接用快速卷积算法$O((n+m)\log_2(n+m))$快速求解。

提供一道题目，卡常题，慎做。

多元多项式卷积

非常神奇的技术，从这个博客学到的。

考虑给定两个$k$元多项式$f(x_1,\ldots,x_k)$和$g(x_1,\ldots,x_k)$，求它们的卷积对$(x_1^{n_1},\ldots,x_k^{n_k})$取模的结果（注意这里取模并不是不进位的意思，而是舍弃发生进位的项的意思）。记$N=\prod_{i=1}^kn_i$

我们可以用一个$0$到$N-1$的下标来表示多项式中的一项。

考虑如何移除进位项带来的影响，我们可以引入一个新的哈希函数$\chi(i)=\sum_{j=1}^{k-1} \lfloor \frac{i}{\prod_{t\leq j} n_t} \rfloor$。考虑到$0\leq \lfloor \frac{a+b}{\prod_{t\leq j} n_t}\rfloor -\lfloor \frac{a}{\prod_{t\leq j} n_t} \rfloor - \lfloor \frac{b}{\prod_{t\leq j} n_t} \rfloor\leq 1$，因此有$0\leq \chi(a+b)-\chi(a)-\chi(b)<k$，其中如果$\chi(a+b)-\chi(a)-\chi(b)=0$，则表示$a+b$过程中不发生进位。

之后我们将$f$重新写作如下的形式$f(x,z)=\sum_{i=0}^{N-1}f_ix^iz^{\chi(i)}$，同理$g(x,z)=\sum_{i=0}^{N-1}g_ix^iz^{\chi(i)}$。我们可以从它们中还原出原先的$k$元多项式。

考虑$f$和$g$中两项的乘积：

$$f_ix^iz^{\chi(i)}\cdot g_jx^jz^{\chi(j)}=f_ig_jx^{i+j}z^{\chi(i)+\chi(j)}[\chi(i+j)=\chi(i)+\chi(j)]$$

其中它们的乘积会被统计到我们结果中，当且仅当$\chi(i)+\chi(j)=\chi(i+j)$。由于$0\leq \chi(a+b)-\chi(a)-\chi(b)<k$以及$\chi(i)$可能比较大，因此我们可以让结果对$z^{k}-1$取模，可以发现信息是不会丢失的。

接下来我们只需要计算$f*g$，并舍弃那些$x$项指数$i$与$t$项指数$j$不匹配的项，即$\chi(i)\neq j$，这些项是发生了进位的。

最后我们要考虑如何快速计算两个二元多项式的卷积。事实上我们可以按$t$的指数对两边进行展开，令$F_i=\sum_{i=0}^{N-1}f_ix^i[\chi(i)=t]$，$G_i=\sum_{i=0}^{N-1}g_ix^i[\chi(i)=t]$：

$$\begin{aligned} f*g=&(\sum_{i=0}^{N-1}f_ix^iz^{\chi(i)}) * (\sum_{i=0}^{N-1}g_ix^iz^{\chi(i)})\\ =&(\sum_{t=0}^{k-1}z^t\sum_{i=0}^{N-1}f_ix^i[\chi(i)=t]) * (\sum_{t=0}^{k-1}z^t\sum_{i=0}^{N-1}g_ix^i[\chi(i)=t])\\ =&(\sum_{t=0}^{k-1}z^tF_t)*(\sum_{t=0}^{k-1}z^tG_t)\\ =&\sum_{t=0}^{k-1}\sum_{p=0}^{k-1}z^{t+p}F_t*G_p \end{aligned}$$

现在问题就变成了$k^2$个多项式的乘积了，一种方案就是直接暴力FFT，时间复杂度为$O(k^2N\log_2N)$，当然我们可以先对每个$F_i$和$G_i$预处理它们的DFT形式，那么公式中卷积就可以简单的用点积来替代，之后把同类项合并，在一同取逆即可，时间复杂度可以减少到$O(kN\log_2N+k^2N)$，考虑到$k\leq \log_2 N$，因此时间复杂度可以简写为$O(kN\log_2N)$。

提供几道题目：

• Multivariate Convolution

话说子集卷积其实也是多元多项式的卷积，对$(x_1^2,\ldots,x_k^2)$进行取模而已。所以同样可以用这个算法解子集卷积，不过由于FFT操作会慢于大部分FWT，因此最好还是用FWT做子集卷积。

参考资料

• 多项式求逆litble博客
• https://www.cnblogs.com/cjyyb/p/10129563.html
• 快速计算数列a的多点等幂和
• 多项式多点求值快速插值
• 多项式牛顿迭代
• [Tutorial] Generalized Chirp Z-Transform
• 集训队互测 2021 Round #1 题解