专栏文章

[笔记] 多项式

jjr_zch
算法·理论参与者 2已保存评论 1

文章操作

快速查看文章及其快照的属性，并进行相关操作。

当前评论
1 条
当前快照
1 份
快照标识符
@mip1ll1n
此快照首次捕获于
2025/12/03 04:40
3 个月前
此快照最后确认于
2025/12/03 04:40
3 个月前
查看原文

多项式

前言

多项式全家桶就是分治倍增分治倍增分治倍增倍增分治分治分治倍增倍增倍增分治分治。
下文中所有系数默认在  mod  p\bmod~p 意义下存在，pp 默认为 998244353998244353。

多项式乘法

  • 加法卷积
  • 位运算卷积

多项式求逆

一些数学细节

给定 n−1n-1 次多项式 A(x)=a0+a1x+a2x2+a3x3+⋯+an−1xn−1A(x)=a_0+a_1x+a_2x^2+a_3x^3+\cdots+a_{n-1}x^{n-1}，求多项式 B(x)B(x) 使得 A(x)⋅B(x)=1A(x) \cdot B(x)=1。
令 A(x)⋅B(x)=C(x)A(x) \cdot B(x)=C(x)。
显然有 a0b0≡c0≡1(modp)a_0b_0 \equiv c_0 \equiv 1 \pmod p，所以当且仅当 a0a_0 在  mod  p\bmod~p 意义下不存在逆元时，A(x)A(x) 不存在逆元。
当 a0≠0a_0 \not= 0 时，易得 b0=1a0b_0=\frac 1 {a_0}。
此时将 b0b_0 乘 A(x)A(x) 的每一项，得到目前 C(x)C(x) 的各项系数，根据系数可以推出使得 c1=0c_1=0 的 b1b_1 的值，然后重复上述步骤。
用这样递推的方法可以得到 B(x)B(x)，但是 B(x)B(x) 的次数显然有可能是无穷大的，并且这样的 B(x)B(x) 在 xx 取大多数值的情况下都是发散的，讨论这个无穷次的多项式本身在 OI\text{OI} 中没有意义。
所以一般我们只需要求出 B(x)B(x) 的前 nn 项，即 A(x)⋅B(x)=1(modxn)A(x) \cdot B(x)=1 \pmod {x^n}。
直接按上述流程暴力模拟的时间复杂度为 O(n2)O(n^2)。

倍增法

考虑对于 22 的幂 mm，m−1m-1 次多项式 Bm(x)B_m(x) 满足 A(x)⋅Bm(x)=1(modxm)A(x) \cdot B_m(x)=1 \pmod{x^m}。
显然，Bm(x)B_m(x) 就是上述无穷次多项式的前 mm 项。
考虑已知 Bm(x)B_m(x)，如何求出 B2m(x)B_{2m}(x)，列出已知条件：
B2m(x)⋅A(x)≡1(modx2m)B2m−Bm≡0(modxm)\begin{aligned} &B_{2m}(x) \cdot A(x) \equiv 1 \pmod {x^{2m}}\\ &B_{2m}-B_m \equiv 0 \pmod{x^m} \end{aligned}
将二式左右两边同时平方并同时乘 A(x)A(x) 后将一式代入：
B2m2(x)+Bm2(x)−2B2m(x)Bm(x)≡0(modx2m)B2m(x)+Bm2(x)A(x)−2Bm(x)≡0(modx2m)B2m(x)≡Bm(x)(2−Bm(x)A(x))(modx2m)\begin{aligned} &B^2_{2m}(x)+B^2_m(x)-2B_{2m}(x)B_m(x) \equiv 0 \pmod {x^{2m}}\\ &B_{2m}(x)+B_m^2(x)A(x)-2B_m(x) \equiv 0 \pmod {x^{2m}}\\ &B_{2m}(x) \equiv B_m(x)(2-B_m(x)A(x)) \pmod {x^{2m}} \end{aligned}
两遍多项式乘法即可完成一次迭代，不难分析出总复杂度为 O(nlog⁡2n)O(n \log_2 n)。
于是我们得到了多项式全家桶的一部分！

多项式带余除法

给定 nn 次多项式 A(x)A(x) 和 mm 次多项式 B(x)B(x)，若 A(x)=B(x)⋅D(x)+R(x)A(x)=B(x) \cdot D(x) + R(x)，求 n−mn-m 次多项式 D(x)D(x) 和 m−1m-1 次多项式 R(x)R(x)。
定义下面四个新的多项式：
A′(x)=A(1x)xnB′(x)=B(1x)xmD′(x)=D(1x)xn−mR′(x)=R(1x)xm−1\begin{aligned} &A'(x)=A(\frac 1 x) x^n\\ &B'(x)=B(\frac 1 x) x^m\\ &D'(x)=D(\frac 1 x) x^{n-m}\\ &R'(x)=R(\frac 1 x) x^{m-1}\\ \end{aligned}
即将系数翻转后得到的新多项式。
由带余除法的定义式容易得到：
A(x)xn=B(x)D(x)xn+R(x)xnA(x)xn=B(x)xm⋅D(x)xn−m+R(x)xm−1⋅xn−m+1A′(x)=B′(x)⋅D′(x)+R′(x)xn−m+1\begin{aligned} &A(x) x^n=B(x)D(x)x^n+R(x)x^n\\ &A(x)x^n=B(x)x^m \cdot D(x)x^{n-m}+R(x)x^{m-1}\cdot x^{n-m+1}\\ &A'(x)=B'(x) \cdot D'(x)+R'(x)x^{n-m+1} \end{aligned}
由于 D′(x)D'(x) 的次数为 n−mn-m，所以 D′(x)=D′(x) mod xn−m+1D'(x)=D'(x) \bmod {x^{n-m+1}}。
我们将上述式子左右两边同时模 xn−m+1x^{n-m+1}。
A′(x)≡B′(x)⋅D′(x)(modxn−m+1)D′(x)=A′(x)B′(x) mod xn−m+1R(x)=A(x)−B(x)⋅D(x)\begin{aligned} &A'(x) \equiv B'(x) \cdot D'(x) \pmod {x^{n-m+1}}\\ &D'(x) = \frac {A'(x)}{B'(x)} \bmod {x^{n-m+1}}\\ &R(x)=A(x)-B(x) \cdot D(x) \end{aligned}
一次求逆和两次乘法即可解决，时间复杂度 O(nlog⁡2n)O(n \log_2 n)。
成功通关四则运算！

多项式开方

这道模板题还有一道加强版。
给定 n−1n-1 次多项式 A(x)A(x)，求 n−1n-1 次多项式 B(x)B(x)，使得 B(x)2≡A(x)(modxn)B(x)^2 \equiv A(x) \pmod {x^n}。
下面 aia_i 表示 A(x)A(x) 的 ii 次项系数，bib_i 同理。

朴素做法

显然有 b02≡a0(modp)b_0^2 \equiv a_0 \pmod p，若 a0a_0 在  mod  p\bmod~p 意义下没有二次剩余，则 B(x)B(x) 无解。
用 Cipolla\text{Cipolla} 或 BSGS\text{BSGS} 等求解二次剩余的算法得出 b0b_0，然后考虑递推。
ak=∑i=0kbibk−ibk=ak−∑i=0k−1bibk−ib0\begin{aligned} &a_k=\sum_{i=0}^k b_i b_{k-i}\\ &b_k=\frac {a_k-\sum_{i=0}^{k-1}b_ib_{k-i}} {b_0} \end{aligned}
根据递推式得出有解的另一个条件是 b0b_0 在  mod  p\bmod~p 意义下要有逆元。若有解，直接暴力即可做到 O(n2)O(n^2) 的复杂度。

倍增法

类比多项式乘法逆的倍增法，设 m−1m-1 次多项式 Bm(x)B_m(x) 满足 Bm(x)2≡A(x)(modxm)B_m(x)^2 \equiv A(x) \pmod {x^{m}}，则有：
B2m(x)2≡A(x)(modx2m)B2m(x)−Bm(x)≡0(modxm)\begin{aligned} &B_{2m}(x)^2 \equiv A(x) \pmod {x^{2m}}\\ &B_{2m}(x)-B_m(x) \equiv 0 \pmod {x^m} \end{aligned}
那么将二式平方就易得：
B2m(x)2+Bm(x)2−2B2m(x)⋅Bm(x)≡0(modx2m)A(x)+Bm(x)2−2B2m(x)⋅Bm(x)≡0(modx2m)B2m(x)≡A(x)2Bm(x)+Bm(x)2(modx2m)\begin{aligned} &B_{2m}(x)^2+B_m(x)^2-2B_{2m}(x) \cdot B_m(x) \equiv 0 \pmod {x^{2m}}\\ &A(x)+B_m(x)^2-2B_{2m}(x) \cdot B_m(x) \equiv 0 \pmod {x^{2m}}\\ &B_{2m}(x) \equiv \frac {A(x)} {2B_m(x)}+ \frac {B_m(x)} 2 \pmod {x^{2m}} \end{aligned}
一次求逆和一次多项式乘法就可以迭代一轮，不难分析总复杂度为 O(nlog⁡2n)O(n \log_2 n)。

多项式多点求值

给定一个 nn 次多项式 F(x)F(x) 和一个长度为 mm 的数组 aa，要求对于 ∀i∈[1,n]\forall i \in [1,n]，求出 f(ai)f(a_i)。

缩减多项式次数

暴力算法显然是 O(nm)O(nm) 的，先考虑一种 mm 较小且 O(m2)O(m^2) 的算法能够满足需求的情形。
尝试将 F(x)F(x) 进行降幂，构造如下的 mm 次多项式 G(x)G(x)。
G(x)=∏i=1m(x−ai)\begin{aligned} G(x)=\prod_{i=1}^m (x-a_i) \end{aligned}
然后令 H(x)=F(x) mod G(x)H(x)=F(x) \bmod G(x)，将 a1…ama_1 \dots a_m 暴力代入 H(x)H(x) 中求答案。
由于对于 ∀i∈[1,n]\forall i \in [1,n]，当 x=aix=a_i 时，G(x)G(x) 的值均为 00，所以 F(x) mod G(x)F(x) \bmod G(x) 本质上是减去了若干个 00，不影响最终的答案，但是用于暴力计算的 H(x)H(x) 的次数就缩减到了 m−1m-1。
G(x)G(x) 可以 O(m2)O(m^2) 暴力计算，也可以分治 NTT\text{NTT} O(mlog⁡22m)O(m \log_2^2 m) 计算，总体能轻松做到 O(m2)O(m^2) 的复杂度。

分治优化

抓住上述优化中模 00 这个有趣的点，如果将 F(x)F(x) 变为 F(x) mod (x−a1)F(x) \bmod (x-a_1)，忽略取模的复杂度，就可以 O(1)O(1) 求出 F(a1)F(a_1) 的值。
考虑结合分治，求解 i∈[l,r]i \in [l,r] 时 F(ai)F(a_i) 的值。
构造一系列与上文中的 G(x)G(x) 和 F′(x)F'(x) 类似的多项式：
Gl,r(x)=∏i=lr(x−ai)Fl,r′(x)=F(x) mod Gl,r(x)\begin{aligned} &G_{l,r}(x)=\prod_{i=l}^r (x-a_i)\\ &F'_{l,r}(x)=F(x) \bmod G_{l,r}(x) \end{aligned}
利用原先分治 NTT\text{NTT} 时求出的所有 Gl,r(x)G_{l,r}(x)，在第二次分治向下递归时，由 Fl,r′(x)F'_{l,r}(x) 去推出 Fl,mid′(x)F'_{l,mid}(x) 和 Fmid+1,r′(x)F'_{mid+1,r}(x) 即可。
这里注意当区间长度较小时，采用暴力而不是大量的 NTT\text{NTT} 会更快。
不难分析出总时间复杂度为 O(nlog⁡22n)O(n \log_2^2 n)。

转置原理

咕咕咕。

优化常系数齐次线性递推

概念

fi=∑j=1kajfi−jf_i=\sum_{j=1}^k a_j f_{i-j}
形如上式这样子的递推式，学名叫做常系数齐次线性递推。
现在我们需要求 fnf_n 的值，其中 n≤1018,k≤32000n \leq 10^{18},k \leq 32000。

朴素算法

在我们的幼年时期，面对需要递推的题目，我们只会 O(n)O(n) 的暴力，更准确地说是 O(nk)O(nk)，不过那时脑子里并没有 kk 这个概念，因为面对的递推式是这样的：
fi=fi−1+fi−2f_i=f_{i-1}+f_{i-2}
被玩烂的斐波那契数列。
后来了解了线性代数初步，其中的矩阵被用来优化这类递推，简单的快速幂可以做到 O(k3log⁡2n)O(k^3 \log_2 n) 的复杂度，具备了解决 kk 较小，但 nn 可以很大的问题的能力。
从矩阵中可以看出后面更优算法的影子。
在用矩阵求解斐波那契数列第 nn 项时，有一种常数较小的做法，它不需要显式地体现矩阵，只需要实现一个“二元组快速幂”。
因为所有的 fif_i 都可以被表示成 pif0+qif1p_if_0+q_if_1，考虑下面的流程。
pp 和 qq 都是广义斐波那契数列，但这与接下来要介绍的无关。
已知 fi=af1+bf0f_i=af_1+bf_0，fi+1=cf1+df0f_{i+1}=cf_1+df_0，推算出用 f0f_0 和 f1f_1 表示 f2if_{2i} 和 f2i+1f_{2i+1} 的系数。
因为 f0f_0 到 fif_i 与 fif_i 到 f2if_{2i} 中间的项数相同，所以有以下推导：
f2i=afi+1+bfi=a(cf1+df0)+b(af1+bf0)=(ac+ab)f1+(ad+b2)f0\begin{aligned} f_{2i}&=af_{i+1}+bf_i\\ &=a(cf_1+df_0)+b(af_1+bf_0)\\ &=(ac+ab)f_1+(ad+b^2)f_0 \end{aligned}
同理易得：
f2i+1=cfi+1+dfi=(ad+c2)f1+(bd+cd)f0\begin{aligned} f_{2i+1}&=cf_{i+1}+df_i\\ &=(ad+c^2)f_1+(bd+cd)f_0 \end{aligned}
由于该做法没有很好的合并性，但是由 fi,fi+1f_i,f_{i+1} 可以推出 fi+2f_{i+2}，所以只能用二分快速幂实现。

特征方程和多项式

一个常系数齐次线性递推式有唯一的特征方程，下面是一个例子：
fi=afi−1+bfi−2+cfi−3+dfi−4f_i=af_{i-1}+bf_{i-2}+cf_{i-3}+df_{i-4}
其特征方程为：
x4=ax3+bx2+cx+dx^4=ax^3+bx^2+cx+d
对应的特征多项式为：
x4−ax3−bx2−cx−dx^4-ax^3-bx^2-cx-d
说人话其实就是把下标映射成 xx 的次数。
显然，特征多项式是等于 00 的。

O(k2log⁡2n)O(k^2 \log_2n) 算法

以斐波那契数列为例，其特征方程为 x2=x+1x^2=x+1，易得下列式子：
x3=x2⋅x=(x+1)⋅x=x2+x=2x+1x4=3x+2x5=5x+3x6=8x+5⋯\begin{aligned} x^3&=x^2 \cdot x\\ &=(x+1)\cdot x\\ &=x^2+x\\ &=2x+1\\ x^4&=3x+2\\ x^5&=5x+3\\ x^6&=8x+5\\ &\cdots \end{aligned}
注意到这等价于原本的递推，相当于用 xx 和 11 表示出了 xx 的若干次方，但多项式的形式给予了更大的操作空间。
沿用快速幂的思想，已知 xi=ax+bx^i=ax+b，试求 x2ix^{2i}。
x2i=(xi)2=(ax+b)2=a2x2+2abx+b2=a2(x+1)+2abx+b2=(a2+2ab)x+(a2+b2)\begin{aligned} x^{2i}&=(x^i)^2\\ &=(ax+b)^2\\ &=a^2x^2+2abx+b^2\\ &=a^2(x+1)+2abx+b^2\\ &=(a^2+2ab)x+(a^2+b^2)\\ \end{aligned}
这与上文中二元组快速幂的最大区别是仅用 xix^i 就能得出 x2ix^{2i}，且 xx 的指数相加可以多项式相乘，所以用倍增快速幂实现。
单轮中将 kk 次多项式暴力平方复杂度为 O(k2)O(k^2)，平方后出现的 k−1k-1 个次数 ≥k\geq k 的项的降幂过程复杂度也为 O(k2)O(k^2)，故总复杂度为 O(k2log⁡2n)O(k^2 \log_2 n)。

多项式技巧的进一步优化

上面算法中，将多项式平方可以使用 NTT\text{NTT} 简单地做到 O(klog⁡2k)O(k \log_2 k)，瓶颈在于降幂的过程。
由于特征多项式值为 00，且次数恰好是 kk，所以把平方后的多项式对特征多项式取模，既没有改变多项式的值，又实现了降幂的目的。
这样就做到了总时间复杂度 O(klog⁡2klog⁡2n)O(k \log_2 k \log_2 n)，但常数较大。

多项式牛顿迭代

英文名是 Newton’s method\text{Newton's method}，很直观了。
这个方法原本是用于求近似方程解，那么列出关于要求的多项式的方程，就可以用牛顿迭代来解出多项式了。

迭代原理

下面为了避免歧义，统一将自变量为 xx 的多项式省略字母后的 (x)(x)。
现在有一个已知多项式 F(P)F(P)，其自变量是一个多项式（注意不表示多项式复合），F(P)F(P) 就是我们所求多项式需要满足的条件。
具体而言，我们要求 n−1n-1 次多项式 BB，使得 F(B)≡0(modxn)F(B) \equiv 0 \pmod {x^n}，当多项式 F(A)F(A) 的次数不为 00 和 11 时，我们就可以采用牛顿迭代法求出 BB。
该方法的核心思想仍为倍增，所以和之前一样，令 m−1m-1 次多项式 BmB_m 满足 F(Bm)≡0(modxm)F(B_m) \equiv 0 \pmod {x^m}。
考虑 F(B)F(B) 在 BmB_m 处的无穷阶泰勒展开：
F(B)=F(Bm)+∑k=1∞F(k)(Bm)k!(B−Bm)kF(B)=F(B_m)+\sum_{k=1}^\infty \frac {F^{(k)}(B_m)} {k!} (B-B_m)^k
因为 B−BmB-B_m 的最低次项次数为 mm，所以 (B−Bm)k(B-B_m)^k 的最低次项次数为 kmkm，故当 k≥2k \geq 2 时有：
(B−Bm)k≡0(modx2m)(B-B_m)^k \equiv 0 \pmod {x^{2m}}
代入展开式可得：
F(B)≡F(Bm)+F′(Bm)(B−Bm)(modx2m)F(B) \equiv F(B_m)+F'(B_m) (B-B_m) \pmod {x^{2m}}
进一步地，根据 BB 的定义可得：
F(Bm)+F′(Bm)(B−Bm)≡0(modx2m)B≡Bm−F(Bm)F′(Bm)(modx2m)B2m=Bm−F(Bm)F′(Bm) mod x2m\begin{aligned} &F(B_m)+F'(B_m) (B-B_m) \equiv 0 \pmod {x^{2m}}\\ &B \equiv B_m-\frac {F(B_m)}{F'(B_m)} \pmod {x^{2m}}\\ &B_{2m}=B_m-\frac {F(B_m)}{F'(B_m)} \bmod {x^{2m}} \end{aligned}
求是求出来了，但感觉很怪。
关于多项式 F(P)F(P) 的求导等操作需要放在具体问题中具体分析，这样会更好理解一些。

实际应用

接下来尝试用牛顿迭代解决多项式求逆。
已知多项式 AA，求它的逆多项式 BB，考虑构造符合题意的多项式 F(P)F(P)。
F(P)=A⋅P−1F′(P)=A\begin{aligned} &F(P)=A \cdot P -1\\ &F'(P)=A \end{aligned}
注意这里的 AA 已知，相当于常量，当 F(B)=0F(B)=0 时，BB 显然是 AA 的逆元，尝试套公式：
B2m=Bm−A⋅Bm−1AB_{2m}=B_m - \frac {A \cdot B_m-1}{A}
发现迭代过程中需要求 AA 的逆元，陷入循环求解，这是因为 F(P)F(P) 的次数为 11，不符合使用牛顿迭代的前提。
F(P)=1P−AF′(P)=−1P2\begin{aligned} &F(P)= \frac 1 P -A\\ &F'(P)=-\frac 1 {P^2}\\ \end{aligned}
如上构造 F(P)F(P)，当 F(B)F(B) 为 00 时，同样能得到 A,BA,B 互为逆元，且次数为 −1-1，符合要求，直接代入牛顿迭代公式。
B2m=Bm−F(Bm)F′(Bm) mod x2mB2m=2Bm−A⋅Bm2\begin{aligned} &B_{2m}=B_m-\frac {F(B_m)}{F'(B_m)} \bmod {x^{2m}}\\ &B_{2m}=2B_m-A \cdot B_m^2 \end{aligned}
这和使用倍增法得到的结果是一样的！
当然，多项式开方也可以用牛顿迭代推式子。

多项式对数函数

给定 n−1n-1 次多项式 A(x)A(x)，其中 [x0]A(x)≠1[x^0]A(x) \not= 1，求 n−1n-1 次多项式 B(x)B(x)，满足 B(x)≡ln⁡A(x)(modxn)B(x) \equiv \ln A(x) \pmod {x^n}。

求导再积分

先对同余式左右两边进行求导：
B′(x)≡A′(x)A(x)(modxn)B'(x) \equiv \frac {A'(x)}{A(x)} \pmod {x^n}
先对 A(x)A(x) 求导得到 A′(x)A'(x)，然后通过多项式求逆算出 B′(x)B'(x)，再对 B′(x)B'(x) 进行积分就能得到 B(x)B(x)，如下式：
B(x)=∫B′(x)dx+CB(x) = \int B'(x)dx + C
注意到需要求 A(x)A(x) 的逆元，所以需要 [x0]A(x)≠0[x^0]A(x) \not= 0，事实上也可以证明当且仅当 [x0]A(x)=1[x^0]A(x)=1 时，存在多项式 B(x)=ln⁡A(x)B(x) = \ln A(x)，但这里不展开讲。
根据微积分基本定理，还能得出积分常数 C=ln⁡[x0]A(x)=ln⁡1=0C=\ln {[x^0]A(x)}=\ln 1=0。
综上，单次求多项式 ln⁡\ln 的时间复杂度为 O(nlog⁡2n)O(n \log_2 n)。

化标

称常数项为 11 的多项式为标准多项式形式，设 B(x)B(x) 为标准多项式，将任意多项式 A(x)A(x) 变为 vxk⋅B(x)vx^k \cdot B(x) 的过程称为多项式化标。
具体地，找到最低的系数不为 00 的项 vxkvx^k，然后把它提出来，就完成了化标。

多项式指数函数

给定 n−1n-1 次多项式 A(x)A(x)，保证 [x0]A(x)=0[x^0]A(x)=0，求 n−1n-1 次多项式 B(x)B(x)，满足 B(x)≡eA(x)(modxn)B(x) \equiv e^{A(x)} \pmod {x^n}。
这一次两边同时求导的方法不太行得通，因为 eA(x)e^{A(x)} 求导后是 eA(x)⋅A′(x)e^{A(x)} \cdot A'(x)，该消失的没有消失，反而多了些东西。
注意到 B(x)B(x) 的常数项 [x0]B(x)=e[x0]A(x)=1[x^0]B(x)=e^{[x^0]A(x)}=1。
考虑牛顿迭代，设 F(P)=ln⁡P−A(x)F(P)=\ln P -A(x)，则 F′(P)=1PF'(P)=\frac 1 P，直接套公式（为避免引起歧义，下面同样省略 (x)(x)）：
B2m≡Bm−F(Bm)F′(Bm)≡Bm−ln⁡Bm−A1Bm≡Bm+Bm(A−ln⁡Bm)(modx2m)\begin{aligned} B_{2m} &\equiv B_m - \frac {F(B_m)}{F'(B_m)}\\ &\equiv B_m-\frac{\ln B_m-A}{\frac 1 {B_m}}\\ &\equiv B_m+ B_m(A-\ln B_m) \pmod {x^{2m}} \end{aligned}
单次迭代复杂度为 O(mlog⁡2m)O(m \log_2 m)，总复杂度显然为 O(nlog⁡2n)O(n \log_2 n)。

多项式快速幂

这道题也有加强版。
有 n−1n-1 次多项式 A(x)A(x) 和整数 kk，求 n−1n-1 次多项式 B(x)B(x)，使得 B(x)≡A(x)k(modxn)B(x) \equiv A(x)^k \pmod {x^n}。
对于同余式两边求 ln⁡\ln 可得：
ln⁡B(x)≡kln⁡A(x)(modxn)\ln {B(x)}\equiv k\ln {A(x)} \pmod {x^n}
求出 ln⁡B(x)\ln B(x) 再 exp⁡\exp 回去即可，但是若 [x0]A(x)≠1[x^0]A(x) \not= 1，ln⁡A(x)\ln A(x) 就不存在，但是这并不意味着 A(x)kA(x)^k 不存在，所以考虑对 A(x)A(x) 进行化标，然后将系数和标准多项式分别 kk 次方即可求解，时间复杂度 O(nlog⁡2n)O(n \log_2 n)。

多项式快速插值

给定 nn 个点 (xi,yi)(x_i,y_i)，求 n−1n-1 次多项式 F(x)F(x)，满足对于 ∀i,F(xi)=yi\forall i,F(x_i)=y_i。
通过拉格朗日插值，可以轻松做到 O(n2)O(n^2) 的复杂度，考虑结合多项式技巧进行优化，先观察拉格朗日插值的式子：
F(x)=∑i=1nyi∏j≠ix−xjxi−xjF(x)=\sum_{i=1}^n y_i \prod_{j \not= i} \frac {x-x_j}{x_i-x_j}
注意到 ∏\prod 里面的分子似乎可以用分治 NTT\text{NTT} 计算，但分母部分与 i,ji,j 都有关，不好处理，将其提出来作为系数，构造下述多项式：
Wi(x)=∏j≠ix−xjW_i(x)=\prod_{j \not= i} x-x_j
注意到对于 ∀k≠i,Wi(xk)=0\forall k \not= i,W_i(x_k)=0，于是可以将所有的 WiW_i 加起来，利用新多项式计算系数，并不影响 x=xix=x_i 时算出的系数大小。
W(x)=∑i=1nWi(x)=∑i=1n∏j≠ix−xjW(x)=\sum_{i=1}^n W_i(x)=\sum_{i=1}^n \prod_{j \not= i} x-x_j
如果已知 W(x)W(x)，对其进行多点求值即可得到所有系数，直接暴力分治 NTT\text{NTT} 可以求得 W(x)W(x)，但每次合并需要两次多项式乘法，常数稍大。
考虑构造 G(x)G(x) 并利用导数的乘法原则对其求导：
G(x)=∏i=1nx−xiG′(x)=∑i=1n(x−xi)′∏j≠ix−xj=∑i=1n∏j≠ix−xj=W(x)\begin{aligned} G(x)&=\prod_{i=1}^n x-x_i\\ G'(x)&=\sum_{i=1}^n (x-x_i)' \prod_{j \not= i} x-x_j\\ &=\sum_{i=1}^n \prod_{j \not= i} x-x_j\\ &=W(x) \end{aligned}
直接对 G(x)G(x) 做分治 NTT\text{NTT} 并求导得到 W(x)W(x) 即可减小常数。
那么现在的式子就变成了：
F(x)=∑i=1nyiW(xi)∏j≠ix−xjF(x)=\sum_{i=1}^n \frac {y_i} {W(x_i)} \prod_{j \not= i} x-x_j
利用刚才的 G(x)G(x) 直接分治 NTT\text{NTT} 求解，具体地，设 Fl,r(x)F_{l,r}(x) 为区间 [l,r][l,r] 的答案，Gl,r(x)G_{l,r}(x) 同理，根据下式合并即可。
Fi,i(x)=yiW(xi)Fl,r(x)=Fl,mid(x)⋅Gmid+1,r(x)+Fmid+1,r(x)⋅Gl,mid(x)\begin{aligned} &F_{i,i}(x)=\frac {y_i}{W(x_i)}\\ &F_{l,r}(x)=F_{l,mid}(x) \cdot G_{mid+1,r}(x)+F_{mid+1,r}(x) \cdot G_{l,mid}(x) \end{aligned}
具体实现上，多点求值和求 F(x)F(x) 可以在同一次分治里分别在递归和回溯的位置实现，且每个区间的 W(x)W(x) 和 F(x)F(x) 都可以用同一个数组在递归和回溯时分别存储。
分治到区间长度较短时可以考虑暴力，速度会快很多。
总时间复杂度 O(nlog⁡22n)O(n \log_2^2 n)。

例题

咕咕咕。

附件

多项式全家桶

使用结构体封装四则运算，采用 vector<int> 存多项式，禁止用 #define int long long！
C
const int mod=998244353;
vector<int> _b,_c,_g,_fac,_inv;

inline int _mod(int x){
	return x>=mod?x-mod:x;
}

inline int qpow(int x,int y){
	int res=1;
	while(y){
		if(y&1) res=(ll)res*x%mod;
		x=(ll)x*x%mod,y>>=1;
	}
	return res;
}

void inv_table(int lim){
	_fac.resize(lim+1);
	_inv.resize(lim+1);
	_fac[0]=1;
	for(int i=1;i<=lim;i++) _fac[i]=(ll)_fac[i-1]*i%mod;
	int now=qpow(_fac[lim],mod-2);
	for(int i=lim;i;i--) _inv[i]=(ll)_fac[i-1]*now%mod,now=(ll)now*i%mod;
	return ;
}

void ntt(int len,vector<int> &f,int typ){
	for(int i=0,j=0;i<len;i++){
		if(i<j) swap(f[i],f[j]);
		for(int k=len>>1;(j^=k)<k;k>>=1);
	}
	_g.resize(len),_g[0]=1;
	for(int step=1;step<len;step<<=1){
		int now=qpow(3,(mod-1)/(step<<1));
		if(typ==-1) now=qpow(now,mod-2);
		for(int i=step-1<<1;i>=0;i-=2){
			_g[i]=_g[i>>1];
			_g[i^1]=(ll)_g[i]*now%mod;
		}
		for(int i=0;i<len;i+=step<<1){
			for(int j=i;j<i+step;j++){
				int x=f[j],y=(ll)f[j+step]*_g[j-i]%mod;
				f[j]=_mod(x+y),f[j+step]=_mod(x-y+mod);
			}
		}
	}
	if(typ==-1){
		int inv=qpow(len,mod-2);
		for(int i=0;i<len;i++) f[i]=(ll)f[i]*inv%mod;
	}
	return ;
}

struct poly{
	int len;
	vector<int> a;
	poly(int _len=0){
		len=_len;
		a.clear();
		a.resize(len);
	}
	poly operator *(const poly &x) const{
		int flen=len+x.len-1;
		flen=1<<(__lg(flen)+(__builtin_popcount(flen)>1));
		poly f=poly(flen);
		_b.resize(flen);
		_c.resize(flen);
		for(int i=0;i<flen;i++){
			_b[i]=i<len?a[i]:0;
			_c[i]=i<x.len?x.a[i]:0;
		}
		ntt(flen,_b,1);
		ntt(flen,_c,1);
		for(int i=0;i<flen;i++) f.a[i]=(ll)_b[i]*_c[i]%mod;
		ntt(flen,f.a,-1);
		f.len=len+x.len-1;
		f.a.resize(f.len);
		return f;
	}
	poly operator ~() const{
		int lim=1<<(__lg(len-1)+2);
		poly f=poly(lim);
		f.a[0]=qpow(a[0],mod-2);
		_b.resize(lim);
		_c.resize(lim);
		for(int step=1;step<len;step<<=1){
			for(int i=0;i<step<<2;i++){
				_b[i]=i<step?f.a[i]:0;
				_c[i]=i<min(step<<1,len)?a[i]:0;
			}
			ntt(step<<2,_b,1);
			ntt(step<<2,_c,1);
			for(int i=0;i<step<<2;i++) _b[i]=(ll)_b[i]*_b[i]%mod*(ll)_c[i]%mod;
			ntt(step<<2,_b,-1);
			for(int i=step;i<step<<1;i++) f.a[i]=(_b[i]>0?mod-_b[i]:0);
		}
		f.len=len;
		f.a.resize(len);
		return f;
	}
	poly operator /(const poly &b) const{
		int flen=len-b.len+1;
		poly x=poly(flen),y=poly(flen);
		for(int i=0;i<flen;i++) x.a[i]=a[len-i-1];
		for(int i=0;i<flen;i++) y.a[i]=i<b.len?b.a[b.len-i-1]:0;
		poly f=x*(~y);
		f.len=flen;
		f.a.resize(flen);
		reverse(f.a.begin(),f.a.end());
		return f;
	}
	poly operator %(const poly &b) const{
		poly x=poly(len);
		for(int i=0;i<len;i++) x.a[i]=a[i];
		if(len<b.len) return x;
		poly d=(x/b)*b;
		poly f=poly(b.len-1);
		for(int i=0;i<f.len;i++) f.a[i]=_mod(a[i]-d.a[i]+mod);
		return f;
	}
};

poly sqrt(poly x){
	int lim=1<<(__lg(x.len-1)+2);
	poly f=poly(lim),g,h;
	f.a[0]=1;
	for(int step=1;step<x.len;step<<=1){
		g=poly(step<<1),h=poly(step<<1);
		for(int i=0;i<step;i++) g.a[i]=_mod(f.a[i]+f.a[i]);
		for(int i=0;i<step<<1;i++) h.a[i]=i<x.len?x.a[i]:0;
		g=h*(~g);
		for(int i=step;i<step<<1;i++) f.a[i]=g.a[i];
	}
	f.len=x.len;
	f.a.resize(x.len);
	return f;
}

poly ln(poly x){
	poly f=poly(x.len-1);
	for(int i=0;i<f.len;i++) f.a[i]=(ll)x.a[i+1]*(i+1)%mod;
	f=f*(~x);
	f.len=x.len;
	f.a.resize(x.len);
	for(int i=f.len-1;i;i--) f.a[i]=(ll)f.a[i-1]*_inv[i]%mod;
	f.a[0]=0;
	return f;
}

poly exp(poly x){
	int lim=1<<(__lg(x.len-1)+2);
	poly f=poly(lim),g,h;
	f.a[0]=1;
	for(int step=1;step<x.len;step<<=1){
		g=poly(step),h=poly(step<<1);
		for(int i=0;i<step;i++) h.a[i]=g.a[i]=f.a[i];
		h=ln(h);
		for(int i=0;i<step<<1;i++) h.a[i]=_mod((i<x.len?x.a[i]:0)-h.a[i]+mod);
		g=g*h;
		for(int i=step;i<step<<1;i++) f.a[i]=g.a[i];
	}
	f.len=x.len;
	f.a.resize(x.len);
	return f;
}

poly qpow(poly x,ll y){
	int pos=0,val=0,ym=y%mod;
	poly b=poly(x.len);
	for(int i=0;i<x.len;i++){
		if(x.a[i]){
			pos=i,val=x.a[i];
			break;
		}
	}
	if(pos&&(y>=x.len||y*pos>=x.len)) return b;
	int inv=qpow(val,mod-2);
	for(int i=0;i<b.len-pos;i++) b.a[i]=(ll)x.a[i+pos]*inv%mod;
	b=ln(b);
	for(int i=0;i<b.len;i++) b.a[i]=(ll)b.a[i]*ym%mod;
	b=exp(b);
	val=qpow(val,(int)(y%(mod-1))),pos*=y;
	for(int i=b.len-1;i>=pos;i--) b.a[i]=(ll)b.a[i-pos]*val%mod;
	for(int i=0;i<pos;i++) b.a[i]=0;
	return b;
}

多项式初等变换依赖图

该图根据上面代码阐述了运算的依赖关系。
由于求逆元过程中的多项式乘法直接调用了 NTT\text{NTT}，所以不依赖乘法。

Finally\text{Finally}

对多项式的评价是又丑又长，实用性较低，建议别学。
Thanks for your reading.\text{Thanks for your reading.}

评论

共 1 条评论，欢迎与作者交流。

正在加载评论...