形式幂级数复合|复合逆

数学 / 多项式

本地源文件：docs/math__poly__comp-rev.md

形式幂级数复合|复合逆

形式幂级数的复合和复合逆也是常见的形式幂级数操作，对于没有特殊性质的 𝑓f 之前我们一直使用的多是 𝑂(𝑛2)O(n2) 的算法来计算 𝑓(𝑔)mod𝑥𝑛f(g)modxn 其中 𝑓 ∈ℂ[[𝑥]],𝑔 ∈𝑥ℂ[[𝑥]]f∈C[[x]],g∈xC[[x]]，但是因为效率较低应用较少．我们介绍 Kinoshita–Li 的 𝑂(𝖬(𝑛)log⁡𝑛)O(M(n)log⁡n) 的算法，其中 𝑂(𝖬(𝑛))O(M(n)) 为两个次数为 𝑂(𝑛)O(n) 的多项式相乘的时间．

形式幂级数/多项式的复合

若要计算 𝑓(𝑔(𝑥))mod𝑥𝑛f(g(x))modxn 那么需要 𝑓(𝑔(𝑥))f(g(x)) 的每一项系数都是有限项之和，所以之前要求 𝑓(𝑥) ∈ℂ[[𝑥]],𝑔(𝑥) ∈𝑥ℂ[[𝑥]]f(x)∈C[[x]],g(x)∈xC[[x]]，而如果 𝑓(𝑥),𝑔(𝑥) ∈ℂ[𝑥]f(x),g(x)∈C[x] 也可以满足这个条件．因为我们需要将 𝑓(𝑔(𝑥))f(g(x)) 的系数截断，不妨直接考虑 𝑓(𝑥),𝑔(𝑥)f(x),g(x) 都是多项式的情况．对于 𝑓(𝑥) =∑𝑛−1𝑗=0𝑓𝑗𝑥𝑗f(x)=∑j=0n−1fjxj，有

𝑓(𝑔(𝑥))=𝑛−1∑𝑗=0𝑓𝑗𝑔(𝑥)𝑗f(g(x))=∑j=0n−1fjg(x)j

我们考虑环 ℂ𝑥)Cx) 上的有理函数

𝑓(𝑦−1)1−𝑦⋅𝑔(𝑥)=∑𝑗≥0(⋯+𝑓𝑗𝑦−𝑗+⋯)𝑔(𝑥)𝑗𝑦𝑗𝑓(𝑔(𝑥))=[𝑦0]𝑓(𝑦−1)1−𝑦⋅𝑔(𝑥)f(y−1)1−y⋅g(x)=∑j≥0(⋯+fjy−j+⋯)g(x)jyjf(g(x))=[y0]f(y−1)1−y⋅g(x)

根据 常系数齐次线性递推 中提到的 Bostan–Mori 算法，Kinoshita 和 Li 指出可以将其修改为二元形式：

𝑃(𝑦)𝑄(𝑥,𝑦)mod𝑥𝑛=(𝑃(𝑦)𝑄(𝑥,𝑦)𝑄(−𝑥,𝑦)mod𝑥𝑛)𝑄(−𝑥,𝑦)mod𝑥𝑛=(𝑃(𝑦)𝑉(𝑥2,𝑦)mod𝑥𝑛)𝑄(−𝑥,𝑦)mod𝑥𝑛=(𝑃(𝑦)𝑉(𝑧,𝑦)mod𝑧⌈𝑛/2⌉)∣𝑧=𝑥2𝑄(−𝑥,𝑦)mod𝑥𝑛P(y)Q(x,y)modxn=(P(y)Q(x,y)Q(−x,y)modxn)Q(−x,y)modxn=(P(y)V(x2,y)modxn)Q(−x,y)modxn=(P(y)V(z,y)modz⌈n/2⌉)|z=x2Q(−x,y)modxn

这样递归的计算在 𝑛 =1n=1 时我们只需计算

𝑃(𝑦)𝑄(𝑥,𝑦)mod𝑥=𝑃(𝑦)𝑄(0,𝑦)∈ℂ((𝑦))P(y)Q(x,y)modx=P(y)Q(0,y)∈C((y))

在计算 𝑃(𝑦)𝑉(𝑧,𝑦)mod𝑧⌈𝑛/2⌉ ∈ℂ𝑧)P(y)V(z,y)modz⌈n/2⌉∈Cz) 时我们不需要保留所有 𝑦y 的系数，因为最后我们只需要提取 𝑦0y0 的系数，所以 𝑦>0y>0 的系数是不需要的，而因为求出前者之后需要将其乘以若干个形如 𝑄( −𝑥,𝑦) ∈ℂ[𝑥,𝑦]Q(−x,y)∈C[x,y] 的「多项式 」，所以只需要保留对于 𝑦0y0 有贡献的系数即可．我们准备好给出伪代码：

𝐀𝐥𝐠𝐨𝐫𝐢𝐭𝐡𝐦 𝖢𝗈𝗆𝗉⁡(𝑃(𝑦),𝑄(𝑥,𝑦),𝑛,𝑚):𝐈𝐧𝐩𝐮𝐭: 𝑃=∑0≤𝑗<𝑛𝑝𝑗𝑦−𝑗∈ℂ((𝑦)),𝑄∈ℂ[𝑥,𝑦],𝑛,𝑚∈ℕ>0.𝐎𝐮𝐭𝐩𝐮𝐭: [𝑦(−𝑚,0]]𝑃(𝑦)𝑄(𝑥,𝑦)mod𝑥𝑛.𝐑𝐞𝐪𝐮𝐢𝐫𝐞: [𝑥0𝑦0]𝑄=1.1𝐢𝐟 𝑛=1 𝐭𝐡𝐞𝐧 𝐫𝐞𝐭𝐮𝐫𝐧 ([𝑦−𝑚+1]𝑃(𝑦)𝑄(0,𝑦),…,[𝑦0]𝑃(𝑦)𝑄(0,𝑦))2𝑉(𝑥2,𝑦)←𝑄(𝑥,𝑦)𝑄(−𝑥,𝑦)mod𝑥𝑛3𝑑←deg𝑦⁡𝑄(−𝑥,𝑦)4(𝑡−(𝑚+𝑑)+1,…,𝑡0)←𝖢𝗈𝗆𝗉⁡(𝑃(𝑦),𝑉(𝑥,𝑦),⌈𝑛/2⌉,𝑚+𝑑)5𝑇(𝑥,𝑦)←∑0𝑗=−(𝑚+𝑑)+1𝑡𝑗𝑦𝑗6𝑈(𝑥,𝑦)=∑𝑑𝑗=−(𝑚+𝑑)+1𝑢𝑗𝑦𝑗←𝑇(𝑥2,𝑦)𝑄(−𝑥,𝑦)mod𝑥𝑛7𝐫𝐞𝐭𝐮𝐫𝐧 (𝑢−𝑚+1,…,𝑢0)Algorithm Comp⁡(P(y),Q(x,y),n,m):Input: P=∑0≤j<npjy−j∈C((y)),Q∈C[x,y],n,m∈N>0.Output: [y(−m,0]]P(y)Q(x,y)modxn.Require: [x0y0]Q=1.1if n=1 then return ([y−m+1]P(y)Q(0,y),…,[y0]P(y)Q(0,y))2V(x2,y)←Q(x,y)Q(−x,y)modxn3d←degy⁡Q(−x,y)4(t−(m+d)+1,…,t0)←Comp⁡(P(y),V(x,y),⌈n/2⌉,m+d)5T(x,y)←∑j=−(m+d)+10tjyj6U(x,y)=∑j=−(m+d)+1dujyj←T(x2,y)Q(−x,y)modxn7return (u−m+1,…,u0)

那么我们有

𝑓(𝑔(𝑥))mod𝑥𝑛=𝖢𝗈𝗆𝗉⁡(𝑓(𝑦−1),1−𝑦⋅𝑔(𝑥),max{1+deg⁡𝑓,𝑛},1)mod𝑥𝑛f(g(x))modxn=Comp⁡(f(y−1),1−y⋅g(x),max{1+deg⁡f,n},1)modxn

注意第三个参数是因为 𝑔(0)g(0) 可能不为零，如果 deg⁡𝑓 ≥𝑛deg⁡f≥n 此时不能截断 𝑓(𝑥)f(x) 来计算 𝑓(𝑔(𝑥))f(g(x))，我们也可以选择计算 𝑓(𝑔) =𝑓 ∘(𝑥+𝑔(0)) ∘(𝑔−𝑔(0))f(g)=f∘(x+g(0))∘(g−g(0))，此时可以取 𝐹 :=𝑓(𝑥+𝑔(0))mod𝑥𝑛F:=f(x+g(0))modxn 和 𝐺 :=𝑔 −𝑔(0)G:=g−g(0) 转而计算 𝖢𝗈𝗆𝗉⁡(𝐹(𝑦−1),1−𝑦⋅𝐺(𝑥),𝑛,1)Comp⁡(F(y−1),1−y⋅G(x),n,1)．

另外因为调用的限制最后递归终止时的 𝑄(0,𝑦)−1Q(0,y)−1 是可以直接导出的，不需要使用形式幂级数的乘法逆元算法来计算，我们只需计算一次乘法然后提取需要的系数．

常见的特殊形式复合

我们常用的 多项式初等函数 都可以通过复合计算：

𝑔(0)=1, 𝑔−1=1+(1−𝑔)+(1−𝑔)2+⋯𝑔(0)=1, log⁡𝑔=−1−𝑔1−(1−𝑔)22−(1−𝑔)33−⋯𝑔(0)=0, exp⁡𝑔=1+𝑔1!+𝑔22!+𝑔33!+⋯𝑔(0)=1, 𝑔𝑒=1+𝑒1(𝑔−1)+𝑒(𝑒−1)2(𝑔−1)2+⋯g(0)=1, g−1=1+(1−g)+(1−g)2+⋯g(0)=1, log⁡g=−1−g1−(1−g)22−(1−g)33−⋯g(0)=0, exp⁡g=1+g1!+g22!+g33!+⋯g(0)=1, ge=1+e1(g−1)+e(e−1)2(g−1)2+⋯

在复合逆的计算中我们也会用到求幂函数．

Kronecker 代换

在分析时间复杂度之前我们先考虑如何作二元多项式乘法，一种想法是将系数「打包」，这一方法由 Kronecker 在 1882 年通过 𝑦 ↦𝑥𝑁y↦xN 将 𝑅[𝑥,𝑦]R[x,y] 上的乘法缩减为 𝑅[𝑥]R[x] 上的乘法，但是要求 𝑁N 足够大．

不妨设 deg𝑥⁡(𝐴𝐵) <𝑁degx⁡(AB)<N，那么我们计算 𝐴(𝑥,𝑥𝑁)𝐵(𝑥,𝑥𝑁)A(x,xN)B(x,xN) 之后仍然可以还原出 𝐴(𝑥,𝑦)𝐵(𝑥,𝑦)A(x,y)B(x,y) 且「打包」和「拆包」的时间为线性．

我们使用 Kronecker 代换再计算一元多项式乘法即可，不难发现在 𝑛n 为二的幂时上述算法可以在 𝑂(𝖬(𝑛)log⁡𝑛)O(M(n)log⁡n) 时间完成，因为每一次递归中 𝑦y 的次数翻倍，但是 𝑥x 的次数减半．

模板（P5373【模板】多项式复合函数）

代码相对于原算法作了一些简化及修改，使得代码更短．

---|---

## 形式幂级数的复合逆

现给出 𝑓 ∈𝑥ℂ[[𝑥]]f∈xC[[x]]![](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7) 且 𝑓′(0) ≠0f′(0)≠0![](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)，求出 𝑔(𝑥)mod𝑥𝑛g(x)modxn![](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7) 满足 𝑓(𝑔) ≡𝑔(𝑓) ≡𝑥(mod𝑥𝑛)f(g)≡g(f)≡x(modxn)![](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)．

根据 [Lagrange 反演](../lagrange-inversion/)，对于 𝑛 >1,𝑘 ≥0n>1,k≥0![](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7) 我们有

[𝑥𝑛−1]𝑓(𝑥)𝑘=𝑘𝑛−1[𝑥𝑛−1−𝑘](𝑔(𝑥)𝑥)−(𝑛−1)[xn−1]f(x)k=kn−1[xn−1−k](g(x)x)−(n−1)![](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)

也就是我们如果能对 𝑘 =0,1,…,𝑛 −1k=0,1,…,n−1![](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7) 求出 [𝑥𝑛−1]𝑓(𝑥)𝑘[xn−1]f(x)k![](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)，那么就可以求出其复合逆．

Kinoshita 和 Li 指出我们可以考虑二元有理函数

11−𝑦⋅𝑓(𝑥)=∑𝑗≥0𝑓(𝑥)𝑗𝑦𝑗11−y⋅f(x)=∑j≥0f(x)jyj![](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)

且这个问题有一个更一般的形式即 Power Projection 问题：我们考虑计算

𝑢:=[𝑥𝑛−1]𝑃(𝑥,𝑦)𝑄(𝑥,𝑦)mod𝑦𝑚u:=[xn−1]P(x,y)Q(x,y)modym![](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)

当 𝑛 −1 =0n−1=0![](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7) 时显然有 𝑢 =𝑃(0,𝑦)𝑄(0,𝑦)mod𝑦𝑚u=P(0,y)Q(0,y)modym![](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)，否则我们有

𝑃(𝑥,𝑦)𝑄(𝑥,𝑦)=𝑃(𝑥,𝑦)𝑄(−𝑥,𝑦)𝑄(𝑥,𝑦)𝑄(−𝑥,𝑦)=𝑈𝑒(𝑥2,𝑦)+𝑥𝑈𝑜(𝑥2,𝑦)𝑉(𝑥2,𝑦)P(x,y)Q(x,y)=P(x,y)Q(−x,y)Q(x,y)Q(−x,y)=Ue(x2,y)+xUo(x2,y)V(x2,y)![](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)

那么

𝑢=⎧{ { {⎨{ { {⎩[𝑥𝑛−1]𝑈𝑒(𝑥2,𝑦)𝑉(𝑥2,𝑦) if 𝑛−1 is even,[𝑥𝑛−1]𝑥𝑈𝑜(𝑥2,𝑦)𝑉(𝑥2,𝑦) if 𝑛−1 is odd.=⎧{ { {⎨{ { {⎩[𝑥⌈𝑛/2⌉−1]𝑈𝑒(𝑥,𝑦)𝑉(𝑥,𝑦) if 𝑛−1 is even,[𝑥⌈𝑛/2⌉−1]𝑈𝑜(𝑥,𝑦)𝑉(𝑥,𝑦) if 𝑛−1 is odd.u={[xn−1]Ue(x2,y)V(x2,y) if n−1 is even,[xn−1]xUo(x2,y)V(x2,y) if n−1 is odd.={[x⌈n/2⌉−1]Ue(x,y)V(x,y) if n−1 is even,[x⌈n/2⌉−1]Uo(x,y)V(x,y) if n−1 is odd.![](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)

我们给出其伪代码：

𝐀𝐥𝐠𝐨𝐫𝐢𝐭𝐡𝐦 𝖯𝗈𝗐𝖯𝗋𝗈𝗃⁡(𝑃(𝑥,𝑦),𝑄(𝑥,𝑦),𝑛,𝑚):𝐈𝐧𝐩𝐮𝐭: 𝑃,𝑄∈ℂ[𝑥,𝑦],𝑛,𝑚∈ℕ>0.𝐎𝐮𝐭𝐩𝐮𝐭: [𝑥𝑛−1]𝑃(𝑥,𝑦)𝑄(𝑥,𝑦)mod𝑦𝑚.𝐑𝐞𝐪𝐮𝐢𝐫𝐞: [𝑥0𝑦0]𝑄=1.1𝐰𝐡𝐢𝐥𝐞 𝑛>1 𝐝𝐨2𝑈(𝑥,𝑦)←𝑃(𝑥,𝑦)𝑄(−𝑥,𝑦)mod𝑥𝑛mod𝑦𝑚3𝐢𝐟 𝑛−1 is even 𝐭𝐡𝐞𝐧4𝑃(𝑥,𝑦)←∑⌈𝑛/2⌉+1𝑗=0([𝑥2𝑗]𝑈(𝑥,𝑦))𝑥𝑗5𝐞𝐥𝐬𝐞6𝑃(𝑥,𝑦)←∑⌈𝑛/2⌉+1𝑗=0([𝑥2𝑗+1]𝑈(𝑥,𝑦))𝑥𝑗7𝐞𝐧𝐝 𝐢𝐟8𝑉(𝑥,𝑦)←𝑄(𝑥,𝑦)𝑄(−𝑥,𝑦)mod𝑥𝑛mod𝑦𝑚9𝑄(𝑥,𝑦)←∑⌈𝑛/2⌉+1𝑗=0([𝑥2𝑗]𝑉(𝑥,𝑦))𝑥𝑗10𝑛←⌈𝑛/2⌉11𝐞𝐧𝐝 𝐰𝐡𝐢𝐥𝐞12𝐫𝐞𝐭𝐮𝐫𝐧 (𝑃(0,𝑦)𝑄(0,𝑦)mod𝑦𝑚)Algorithm PowProj⁡(P(x,y),Q(x,y),n,m):Input: P,Q∈C[x,y],n,m∈N>0.Output: [xn−1]P(x,y)Q(x,y)modym.Require: [x0y0]Q=1.1while n>1 do2U(x,y)←P(x,y)Q(−x,y)modxnmodym3if n−1 is even then4P(x,y)←∑j=0⌈n/2⌉+1([x2j]U(x,y))xj5else6P(x,y)←∑j=0⌈n/2⌉+1([x2j+1]U(x,y))xj7end if8V(x,y)←Q(x,y)Q(−x,y)modxnmodym9Q(x,y)←∑j=0⌈n/2⌉+1([x2j]V(x,y))xj10n←⌈n/2⌉11end while12return (P(0,y)Q(0,y)modym)![](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)

同样的我们也可以直接导出 𝑄(0,𝑦)−1Q(0,y)−1![](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7) 而不需要计算形式幂级数的乘法逆元，那么复合逆的算法就是

𝐀𝐥𝐠𝐨𝐫𝐢𝐭𝐡𝐦 𝖱𝖾𝗏⁡(𝑓(𝑥),𝑛):𝐈𝐧𝐩𝐮𝐭: 𝑓∈𝑥ℂ[[𝑥]],𝑓′(0)≠0,𝑛∈ℕ≥2.𝐎𝐮𝐭𝐩𝐮𝐭: 𝑔(𝑥)mod𝑥𝑛 such that 𝑓(𝑔)≡𝑔(𝑓)≡𝑥(mod𝑥𝑛).1𝑡←𝑓′(0)2𝐹(𝑥)←𝑓(𝑡−1𝑥)3∑𝑛−1𝑘=0𝑎𝑘𝑦𝑘←𝖯𝗈𝗐𝖯𝗋𝗈𝗃⁡(1,1−𝑦⋅𝐹(𝑥),𝑛,𝑛)4𝐺(𝑥)←∑𝑛−1𝑘=1𝑛−1𝑘𝑎𝑘𝑥𝑛−1−𝑘5𝐻(𝑥)←(𝐺(𝑥)1/(𝑛−1))−1mod𝑥𝑛−16𝐫𝐞𝐭𝐮𝐫𝐧 ((𝑡−1𝑥)∘(𝑥⋅𝐻))Algorithm Rev⁡(f(x),n):Input: f∈xC[[x]],f′(0)≠0,n∈N≥2.Output: g(x)modxn such that f(g)≡g(f)≡x(modxn).1t←f′(0)2F(x)←f(t−1x)3∑k=0n−1akyk←PowProj⁡(1,1−y⋅F(x),n,n)4G(x)←∑k=1n−1n−1kakxn−1−k5H(x)←(G(x)1/(n−1))−1modxn−16return ((t−1x)∘(x⋅H))![](data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7)模板（[P5809【模板】多项式复合逆](https://www.luogu.com.cn/problem/P5809)）

代码相对于原算法作了一些简化及修改，使得代码更短．

---|---

由转置原理导出

Power Projection 问题是 Modular Composition 的转置，Kinoshita 和 Li 指出我们前文的复合算法可以由 Power Projection 算法直接转置得到．同样的，如果优化可以应用于 Power Projection 算法，其也可以应用于 Modular Composition 算法．我们省略细节．

参考文献

  1. Yasunori Kinoshita, Baitian Li.Power Series Composition in Near-Linear Time. FOCS 2024.
  2. Alin Bostan, Ryuhei Mori.A Simple and Fast Algorithm for Computing the N-th Term of a Linearly Recurrent Sequence. SOSA 2021: 118-132
  3. R. P. Brent and H. T. Kung. 1978.Fast Algorithms for Manipulating Formal Power Series. J. ACM 25, 4 (Oct. 1978), 581–595.
  4. Daniel J. Bernstein. "Fast multiplication and its applications." Pages 325–384 in Algorithmic number theory: lattices, number fields, curves and cryptography, edited by Joe Buhler, Peter Stevenhagen, Cambridge University Press, 2008, ISBN 978-0521808545.
本页面最近更新： 2026/1/7 08:56:54，更新历史 发现错误？想一起完善？在 GitHub 上编辑此页！ 本页面贡献者：Tiphereth-A, c-forrest, hly1204 本页面的全部内容在CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用