多项式学习笔记（1）

Phartial · 2026-05-07 09:59:16 · 算法·理论

前情提要。

本文又名《李白天信息学竞赛中的生成函数计算理论框架》[^5]学习笔记（但行文结构可能不太一样，而且会加一些里面没讲的东西）。理解可能会较浅显，欢迎交流批评。

因为一次写完不可能（其实是因为我是懒狗），所以我把文章拆成了多篇，这是第一篇，主要讲解一些基础理论。（第二篇绝赞咕咕咕中！！！）

默认在 \mathbb{F}_{\mathsf{NTT}} 域中做运算，也即模数为形如 m\times 2^k+1 的质数时的模域。

感谢 Gemini 老师解答了我的很多疑惑。

概述

使用生成函数解决问题时，我们关心的无非两点：

怎么用形式幂级数描述问题；
怎么从形式幂级数中提取我们需要的答案。

描述问题的方式被称为符号化方法，本文作为第一篇并不会介绍这一部分的理论。

而提取时，难点会出现在以下对象：

复合：当描述的问题出现子结构时，复合便自然出现。得益于 noshi91 与 EI 在 2024 年的研究[^3]，多项式复合在如今几乎不会成为计算时的瓶颈；
方程：当描述的问题出现递归/自指时，方程便自然出现。由于方程形式的灵活性，其求解有时会相当困难。但我们还是有三种强大的工具能处理幂级数方程：在线算法、牛顿迭代与拉格朗日反演；
远处求值：当数据范围过大时，远处求值便自然出现。这类问题常常要求我们在提取 n 次项系数时在 o(n) 时间内完成计算。目前仅有两类远处求值的理论是较为成熟的：线性递推与整式递推。

本文期望能对上述问题与工具给出一个梳理性的介绍。

这篇文章会讲什么？

前置知识：大概地过一下多项式的基础定义与算法；
在线算法：讲解在线问题、（半）在线卷积与一般的在线算法设计思路；
拉格朗日反演；
牛顿迭代；
远处求值：讲解线性递推与整式递推的远处求值；
再谈远处求值：
1. 线性递推：讲解有理函数重建与 Half-GCD 算法；
2. 整式递推：讲解 D-Finite 相关理论。

1 前置知识

设 \mathbb{A} 为一个特征为 0 的域。

定义 1.1（多项式） 定义 \mathbb{A} 上的一个以 z 为形式变元的多项式 F(z)（无歧义时可简记为 F，下同）为：
F(z)=\sum_{i=0}^n f_iz^i
满足 f_n\ne 0。记 [z^k]F=f_k 为 F 的 k 次项系数，记 \deg F=n 为 F 的度数。

特别的，我们单独定义零多项式 F(z)=0，其度数被看作是 -\infty。

容易定义其上的加法与乘法，并验证由 \mathbb{A} 上的所有多项式所形成的代数结构是交换环。我们记其为多项式环 \mathbb{A}[z]。

定理 1.1（多项式乘法） 计算两个度数不超过 n 的多项式的乘积可以在 \mathsf{M}(n)=\Theta(n\log n) 的时间内完成。

值得一提的是本文介绍的大部分工具都不会强依赖于 FFT，也就是说就算你只会暴力卷积，你也可以简单地把下文中的 \mathsf{M}(n) 替换为 \Theta(n^2) 来得到一系列基于暴力卷积的多项式算法（不过直接这么做可能会多 \log，建议还是手推一下）。

定理 1.2（多项式多点求值） 已知多项式 F 与 n 个数 x_1,\dots,x_n，则我们可以在 \Theta(\mathsf{M}(n)\log n) 的时间内计算出 F(x_1),\dots,F(x_n)。
定理 1.3（多项式多点插值） 已知 n 对点值 (x_1,y_1),\dots,(x_n,y_n) 满足 x_i 互不相同，则唯一存在 \deg F<n 的多项式 F 使得 \forall i,F(x_i)=y_i，且我们可以在 \Theta(\mathsf{M}(n)\log n) 的时间内计算出 F 的系数。

通过令多项式的度数取到无穷大，我们得到形式幂级数：

定义 1.2（形式幂级数） 定义 \mathbb{A} 上的一个以 z 为形式变元的形式幂级数（无歧义时可简称幂级数）F(z) 为：
F(z)=\sum_{i\ge 0} f_iz^i
记 [z^k]F=f_k。

容易验证形式幂级数仍形成一个环（因为只有 [z^0]F 可逆的幂级数可以求逆，所以还不是域），我们将其记作形式幂级数环 \mathbb{A}[[z]]。

这里讨论一下稍微复杂一点的复合与复合逆。

定义 1.3（幂级数复合） 两个形式幂级数 F,G 的复合 F\circ G（或 F(G)）被定义为 \displaystyle\sum_{i=0}^\infty ([z^i]F)G^i。
定义 1.4（幂级数复合逆） 我们称形式幂级数 G 是形式幂级数 F 的复合逆，当且仅当 F\circ G=z，记作 G=F^{\lang -1\rang}。

首先可以发现的是复合与复合逆的合法性并不那么显然了，为了避免讨论 \mathbb{A} 上的极限，我们在这里只规定两种复合是合法的：F 只有有限项非零，或是 [z^0]G=0。

并且我们只讨论 [z^0]F=0\land [z^1]F\ne 0 的 F 的复合逆，此时可以证明复合逆存在且唯一，证明思路是我们可以直接递推出来。

接下来有两种方法能将多项式扩张成一个域，一种方法是直接将对象记作两多项式之比，这给出有理函数：

定义 1.5（有理函数） 定义 \mathbb{A} 上的一个以 z 为形式变元的有理函数 F(z) 为：
F(z)=\frac{P(z)}{Q(z)}
满足 P(z) 与 Q(z) 均为 \mathbb{A} 上的多项式，且 Q(z)\ne 0。

容易验证有理函数形成一个域，我们记其为分式域 \mathbb{A}(z)。

另一种方法是添加负次项，这给出形式 Laurent 级数：

定义 1.6（形式 Laurent 级数） 定义 \mathbb{A} 上的一个以 z 为形式变元的形式 Laurent 级数 F(z) 为：
F(z)=\sum_{i\ge n_0}f_iz^i\\
或等价表达：
F(z)=z^{n_0}\sum_{i\ge 0}f_{i+n_0}z^i
满足 f_{n_0}\ne 0，n_0 可以是负数。记 [z^k]F=f_k，\operatorname{ord} F=n_0 为 F 的阶数，\operatorname{res} F=[z^{-1}]F 为 F 的形式留数。

我们同样特殊定义 \operatorname{ord}0=+\infty。

容易验证这给出了形式 Laurent 级数域 \mathbb{A}((z))。并且形式 Laurent 级数上复合的收敛性与形式幂级数一致。

假定读者知道一些基础的生成函数知识。

2 什么是在线算法？

在离线算法中，假设我们想计算 G(F(z))\bmod z^{n}，则必须一次性给出 F(z)\bmod z^n。然而在幂级数方程中输入可能由输出决定，这给出在线算法。

定义 2.1（后效变换） 一个变换 \Psi:\mathbb{A}[[ z]]\times\cdots\times\mathbb{A}[[ z]]\to \mathbb{A}[[ z]] 被称为是后效的，当且仅当对任意 i 与一组参数 F_1(z),\dots,F_n(z)，记 \psi(z)=\Psi(F_1(z),\dots,F_n(z))，则 \psi(z)\bmod z^i 仅依赖于 F_1(z)\bmod z^i,\dots,F_n(z)\bmod z^i。
引理 2.1 卷积变换 \Psi(F,G)=F\times G 是后效的。
引理 2.2 左复合一个多项式 F 的变换 \Psi(G)=F(G) 是后效的。
定义 2.2（在线算法） 对于一个有 n 个输入的后效变换 \Psi 与一个黑盒，一个算法被称为是在线的，当且仅当其能实现这一过程：
1. 令输入 F_1(z),\dots,F_n(z) 均为 0；

这就给了我们一种解决问题的通用方法：

引理 2.3 设问题规模为 n，且有一系列其在线算法用时分别为 T_1(n),\dots,T_k(n) 的后效变换 \Psi_1,\dots,\Psi_k 所构成的表达式树，则其整体所形成的变换 \Phi 是后效的，且存在用时为 T_1(n)+\dots+T_k(n) 的在线算法。

正如多项式算法的基础是卷积，在线算法的基础自然是在线卷积。我们先考虑单输入的半在线卷积：

定理 2.1（半在线卷积）[^1] 给定多项式 G，定义后效变换 \Psi(F)=F\times G，则计算这一变换的在线算法可以做到 \mathsf{S}(n)=\Theta(n\log ne^{\sqrt{2\log 2\log\log n}}\sqrt{\log\log n}) 的时间复杂度。

但事实上 OI 里较为 pratical 的半在线卷积是 \Theta(\frac{n\log^2 n}{\log\log n}) 的。它基于一个 \Theta(n\log^2 n) 的算法。其流程为：

对 f_{0\dots\frac{n}{2}-1} 与 g_{0\dots\frac{n}{2}-1} 做半在线卷积；
对 f_{0\dots\frac{n}{2}-1} 与 g_{\frac{n}{2}\dots n-1} 做离线卷积；
对 f_{\frac{n}{2}\dots n-1} 与 g_{0\dots\frac{n}{2}-1} 做半在线卷积。
对 f_{\frac{n}{2}\dots n-1} 与 g_{\frac{n}{2}\dots n-1} 做离线卷积。

如图所示[^2]：

该算法用时为 \mathsf{S}(n)=2\mathsf{S}(\frac{n}{2})+2\mathsf{M}(\frac{n}{2})，根据主定理有 \mathsf{S}(n)=\Theta(n\log^2 n)。

而我们可以将其改为多叉的，将 f 和 g 以 b 为块长分块，则有这样的算法流程：

对 i 从 0 到 \frac{n}{b}-1，执行：
1. 对 f_{ib\dots (i+1)b-1} 与 g_{0\dots b-1} 做半在线卷积；
2. 对 f_{ib\dots (i+1)b-1} 与 g_{b\dots 2b-1}、g_{2b\dots 3b-1}、……各做一次离线卷积。

那么我们就只需要做 2\cdot\frac{n}{b}-1 次规模为 2b 的 DFT，随后是 (\frac{n}{b})^2 次规模为 b 的点乘，最后我们不需要做 (\frac{n}{b})^2 次规模为 b 的 IDFT，因为我们可以把贡献到的下标范围相同的点值加起来再 IDFT 回去，这样只要做 2\cdot\frac{n}{b}-2 次规模为 2b 的 IDFT。

于是该算法用时为 \mathsf{S}(n)=\frac{n}{b}\mathsf{S}(b)+\Theta(n\log n+\frac{n^2}{b})，取 b=\frac{n}{\log n} 可得 \mathsf{S}(n)=\Theta(\frac{n\log^2 n}{\log\log n})。

而可以证明的是，半在线卷积和在线卷积是等难的。

定理 2.2（在线卷积） 定义后效变换 \Psi(F,G)=F\times G，则计算这一变换的在线算法可以做到 \mathsf{R}(n)=\Theta(\mathsf{S}(n)) 的时间复杂度。

证明：考虑如下算法流程：
1. 对 f_{0\dots\frac{n}{2}-1} 与 g_{0\dots\frac{n}{2}-1} 做在线卷积；
2. 并行地对 f_{0\dots\frac{n}{2}-1}\times g_{\frac{n}{2}\dots,n-1} 与 f_{\frac{n}{2}\dots n-1}\times g_{0\dots\frac{n}{2}-1} 做半在线卷积；
3. 对 f_{\frac{n}{2}\dots n-1} 与 g_{\frac{n}{2}\dots n-1} 做离线卷积。
不难发现并行的半在线卷积与一般的半在线卷积是等难的。所以该算法用时为 \mathsf{R}(n)=\mathsf{R}(\frac{n}{2})+\mathsf{S}(\frac{n}{2})+\mathsf{M}(\frac{n}{2})，根据主定理，有 \mathsf{R}(n)=\Theta(\mathsf{S}(n))。

以在线卷积为基础，根据引理 2.3，我们就可以得到一系列多项式初等函数的在线算法了。

逆元：在线地给出 F，求 G=\frac{1}{F}。

考虑 \displaystyle g_n=-\frac{1}{f_0}\sum_{i=0}^{n-1}g_if_{n-i}，记 h_i=f_{i+1}，则 g_n=-\frac{[z^{n-1}]H(z)G(z)}{f_0}，这就是一个在线卷积的形式了。
直接求导，有 $G'=GF'$，同样是一个在线卷积的形式。
直接求导，有 $G'=\frac{F'}{F}$，把在线逆元和在线卷积拼起来即可。

同时在线算法还能用来解幂级数方程，我们以一道例题来讲解其过程。

例 2.1（无标号有根树计数） 解方程 \displaystyle F(z)=A(z)\exp\left(\sum_{i\ge 1}\frac{F(z^i)}{i}\right)，其中 [z^0]A(z)=0。

分析常数项可以发现 [z^0]F(z)=0，而在知道 [z^n]F(z) 后，我们可以更新 \displaystyle [z^n]G(z)=[z^n]\sum_{i\ge 1}\frac{F(z^i)}{i}，进而更新 [z^n]\exp G(z)，从而在乘上 A(z) 后更新 [z^{n+1}]F(z)。这是一个在线的过程，将刚才提到的算法都替换为对应的在线算法即可。

3 什么是拉格朗日反演？

借助拉格朗日反演，我们得以建立原函数与其复合逆系数的关系。而这基于形式 Laurent 级数上特殊的形式留数，它的特殊之处在于求导，因为 z^0 求导后系数不会传递给 z^{-1}。

引理 3.1 对任意形式 Laurent 级数 F(z)，\operatorname{res}F'(z)=0。
引理 3.2 对任意阶数为 1 的形式 Laurent 级数 F(z)，\operatorname{res}F'(z)F^k(z)=[k=-1]。

证明：当 k\ne -1 时，F'(z)F^k(z)=(\frac{1}{k+1}F^{k+1}(z))'，根据引理 3.1，其 z^{-1} 项系数为 0。

当 k=-1 时，\displaystyle\frac{F'(z)}{F(z)}=\frac{\sum_{i\ge 0} (i+1)f_{i+1}z^{i}}{z\sum_{i\ge 0}f_{i+1}z^i}，容易验证该式的 z^{-1} 项系数是 1。

事实上观察上述证明过程可以发现这一引理完全可以扩展到 \operatorname{ord}F\ne 1 的情况，此时有 \operatorname{res}F'(z)F^k(z)=[k=-1]\operatorname{ord}F，不过和接下来要干的事情没啥关系，所以不管了。

到这里你大概能猜出来我们接下来要做什么了：多项式复合能同时对不同的 k 提供 F^k(z)，那么借助引理 3.2，我们只需要乘上 F'(z) 再取 [z^{-1}] 就能筛选出满足 [k=-1] 的那一项，从而做到定向提取系数。

定理 3.1（拉格朗日反演） 对任意阶数为 1 的形式 Laurent 级数 F(z)，对任意整数 n,k，有：
n[z^n]F^k(z)=k[z^{-k}]\left(F^{\lang -1\rang}(z)\right)^{-n}
证明：记 G=F^{\lang -1\rang}。我们从 F(G(z))=z 开始，先将两侧求 k 次方得到：
F^k(G(z))=z^k
记 F_0=F^k，再求导，有：
G'(z)F_0'(G(z))=kz^{k-1}
展开得到：
G'(z)\sum_i i([z^i]F_0) G^{i-1}=kz^{k-1}
此时已经有了引理 3.2 的形式，为了锁定 i=n 这一项，我们希望得到 [i-1-n=-1] 的艾弗森括号，于是将等式两边同时乘上 G^{-n} 再取形式留数：
\begin{aligned} \sum_i i([z^i]F_0)\operatorname{res}G'G^{i-1-n}&=\operatorname{res}kz^{k-1}G^{-n}\\ \sum_i i([z^i]F_0)[i-1-n=-1]&=k[z^{-k}]G^{-n}\\ n[z^n]F^k&=k[z^{-k}]G^{-n}\\ \end{aligned}
即证。

根据线性性，我们还可以得到拉格朗日反演的扩展版本：

推论 3.1（扩展拉格朗日反演） 对任意形式 Laurent 级数 H(z) 与阶数为 1 的形式 Laurent 级数 F(z)，对任意整数 n，有：
n[z^n]H(F(z))=[z^{-1}]H'(z)\left(F^{\lang -1\rang}(z)\right)^{-n}
或更常用的形式：
n[z^n]H(F(z))=[z^{n-1}]H'(z)\left(\frac{z}{F^{\lang -1\rang}(z)}\right)^{n}

通过修改证明过程，我们还能得到拉格朗日反演的变体：

定理 3.2（另类拉格朗日反演） 对任意阶数为 1 的形式 Laurent 级数 F(z) 与任意整数 n,k，设 G=F^{\lang -1\rang}，则有：
[z^n]F^k=[z^{-k-1}]G'G^{-n-1}
证明：设 F_0=F^k，我们仍然从 F_0(G(z))=z^k 开始，但此时不再求导，而是直接乘上 G' 得到：
G'\sum_i ([z^i]F_0)G^i=z^kG'
构造 [i-n-1=-1]，在等式两侧乘上 G^{-n-1} 并提取形式留数得到：
\begin{aligned} \sum_i ([z^i]F_0)[i-n-1=-1]&=\operatorname{res}z^kG'G^{-n-1}\\ [z^n]F^k&=[z^{-1}]z^kG'G^{-n-1}\\ \end{aligned}

它同样也有其扩展版本：

推论 3.2 对任意形式 Laurent 级数 H(z) 与阶数为 1 的形式 Laurent 级数 F(z)，对任意整数 n，设 G=F^{\lang -1\rang}，有：
[z^n]H(F(z))=[z^{n}]HG'\left(\frac{z}{G}\right)^{n+1}

这里应该有一些例题，但是我咕了。

4 什么是牛顿迭代？

我们熟知实数域上的泰勒展开：

f(z)=\sum_{i\ge 0}\frac{f^{(i)}(z_0)}{i!}(z-z_0)^i

借助泰勒展开，我们得到了一种在方程求解领域相当强大的工具：牛顿迭代。我们自然会希望得到它在幂级数环上的自然扩展。

一般的幂级数方程描述起来相当麻烦，因此我们这里只考虑一类叫做代数方程的幂级数方程。

定义 4.1（代数方程） 一个代数方程基于一个系数 \in\mathbb{A}[[z]] 的幂级数 G(t)\in\mathbb{A}[[z]][[t]]，其解为满足 [z^0]y=0\land G(y)=0 的幂级数 y\in\mathbb{A}[[z]]。

（这里的 G(t) 也能看成一个二元幂级数 G(z,t)。于是我们可能会混用一些 G(z,t) 与 G(t)，请读者自行鉴别。）

而我们可以证明，实数上的泰勒展开在 \mathbb{A}[[z]][[t]] 上仍然成立。

定理 4.1（泰勒展开） 已知幂级数 G(t)\in\mathbb{A}[[z]][[t]] 与 t_0\in\mathbb{A}[[z]] 满足 [z^0]t_0=0，则：
G(t)=\sum_{i\ge 0}\frac{G^{(i)}(t_0)}{i!}(t-t_0)^i
证明：展开 G(t+t_0) 可得：
\begin{aligned} G(t+t_0)&=\sum_{i\ge 0}([t^i]G)(t+t_0)^i\\ &=\sum_{i\ge 0}([t^i]G)\sum_{j=0}^i {i\choose j}t^jt_0^{i-j}\\ &=\sum_{j\ge 0}\frac{t^j}{j!}\sum_{i\ge 0}(i+j)^{\underline{j}}([t^{i+j}]G)t_0^i\\ &=\sum_{j\ge 0}\frac{t^j}{j!}G^{(j)}(t_0)\\ \end{aligned}

猜测牛顿迭代也是成立的。

定理 4.2（牛顿迭代） 已知幂级数 G(t)\in\mathbb{A}[[z]][[t]] 与幂级数 y_0\in\mathbb{A}[[z]] 满足 [z^0]G(y_0)=0（若 G(t) 有无穷项非零，则为了使 G(y_0) 收敛，必须有 [z^0]y_0=0）且 G'(y_0) 可逆（也即 [z^0]G'(y_0)\ne 0）。则存在唯一的 y 使得 G(y)=0 且 [z^0]y=[z^0]y_0。

证明：先证明 y 存在。考虑构造序列 \{y_i\}，构造方法正是我们熟知的牛顿迭代 y_{n+1}=y_n-\dfrac{G(y_n)}{G'(y_n)}。设 \delta_n=-\dfrac{G(y_n)}{G'(y_n)}，我们需要证明的只有如下几点：
3. 序列 \{y_i\} 收敛，我们实际上可以证明 \delta_n\equiv 0\pmod{z^{2^n}}，所以这一条是显然的。
考虑归纳，对 G(y_{n+1}) 在 y_n 处执行泰勒展开，那么有：
\begin{aligned} G(y_{n+1}) &\equiv\sum_{i\ge 0}\frac{G^{(i)}(y_n)}{i!}\delta_n^i\pmod{z^{2^{n+1}}}\\ &\equiv G(y_n)+G'(y_n)\delta_n\pmod{z^{2^{n+1}}}\\ &\equiv 0\pmod{z^{2^{n+1}}}\\ \end{aligned}
第二个等号是根据 \delta_n\equiv 0\pmod{z^{2^n}} 得到的 \delta_n^2\equiv 0\pmod{z^{2^{n+1}}}，因此 i\ge 2 的项都被消去。同时也因为 [z^0]\delta_n=0，所以 [z^0]y_{n+1}=[z^0]y_n=[z^0]y_0，于是 [z^0]G'(y_n)=[z^0]G'(y_0)\ne 0。

那么 \delta_{n+1} 是良定义的，又因为 G(y_{n+1})\equiv 0\pmod{z^{2^{n+1}}}，所以 \delta_{n+1} 也 \equiv 0\pmod{z^{2^{n+1}}}，这就证明了 y 的存在性。
对 $G(y^*)$ 在 $y$ 处泰勒展开得到 $0=G(y^*)=G(y)+G'(y)\Delta+O(\Delta^2)$，因为 $G(y)=0$，所以我们有 $\Delta(G'(y)+O(\Delta))=0$。因为 $[z^0]G'(y)\ne 0$ 但 $[z^0]O(\Delta)=0$，所以一定有 $\Delta=0$，于是 $y=y^*$，这就证明了 $y$ 的唯一性。

和在线卷积相比，牛顿迭代在方程求解的时间复杂度上更胜一筹。这是因为从 y_n 推到 y_{n+1} 只需要在 {\!}\bmod z^{2^{n+1}} 意义下运算，因此计算 y\bmod z^{N} 的用时会满足 T(N)=T(\frac{N}{2})+R(N)，而绝大多数情况下 R(N)=\Omega(N)，所以有 T(N)=\Theta(R(N))。

但是牛顿迭代的常数可能会更差，这是因为计算 \delta_n 时我们可能会调用其他过程，比如另一次牛顿迭代，而常数便会随之翻倍。一些精细实现可以大幅优化其常数，本文不提。

我们简单地提一下常见函数怎么使用牛顿迭代计算：

逆元：给出 F，求 G=\frac{1}{F}。

不妨假设 [z^0]F=1，考虑设 H(t)=\frac{1}{1-t}-F，则方程 H(t)=0 的解就是 1-G。这个古怪的形式是为了方便我们将其表示成幂级数的形式 H(t)=(1-F)+t+t^2+\cdots。
$$ \begin{aligned} y_{n+1} &=y_n-\frac{H(y_n)}{H'(y_n)}\\ &=y_n-\frac{(1-y_n)^{-1}-F}{(1-y_n)^{-2}}\\ &=y_n-((1-y_n)-F(1-y_n)^2)\\ &=1-(1-y_n)(2-F(1-y_n))\\ \end{aligned} $$ 实现时直接对 $1-y_n$ 做递推即可，最后得到的就是 $G$ 了。用时满足 $T(N)=R(N)=\Theta(\mathsf{M}(N))$。
开方：给出 F，求 G=\sqrt{F}。

设 H(t)=t^2-F，则方程 H(t)=0 的解就是 G。
$$ \begin{aligned} y_{n+1} &=y_n-\frac{H(y_n)}{H'(y_n)}\\ &=y_n-\frac{y_n^2-F}{2y_n}\\ &=\frac{y_n^2+F}{2y_n}\\ \end{aligned} $$ 递推时需要求一次逆，不过因为求逆的用时也是 $\Theta(\mathsf{M}(N))$ 的，故总用时仍满足 $T(N)=R(N)=\Theta(\mathsf{M}(N))$。
对数/指数：只要能计算其中一者就能计算另一者。因为 \ln' F=\dfrac{F'}{F} 相对容易计算，我们只看指数函数怎么算。也就是给出 F 满足 [z^0]F=0，求 G=\exp F。

设 H(t)=\ln (1+t)-F，则方程 H(t)=0 的解就是 G-1。
$$ \begin{aligned} y_{n+1} &=y_n-\frac{H(y_n)}{H'(y_n)}\\ &=y_n-(1+y_n)(\ln(1+y_n)-F)\\ &=-1+(1+y_n)(1+F-\ln(1+y_n))\\ \end{aligned} $$ 实现时直接对 $1+y_n$ 做递推即可，最后得到的就是 $G$。用时仍满足 $T(N)=R(N)=\Theta(\mathsf{M}(N))$。

而只要对牛顿迭代做一些扩展，我们还能处理更复杂的一类方程，比如含有下标偏移的方程：

例 2.1' 解方程 \displaystyle F(z)=A(z)\exp\left(\sum_{i\ge 1}\frac{F(z^i)}{i}\right)，其中 [z^0]A(z)=0。

我们在第二章里给出了一个基于在线算法的做法，现在我们希望给出一个基于牛顿迭代的做法。

这个方程不是代数方程，因为其中出现了下标偏移 F(z^i)。但是重新考察序列 \{y_i\} 的含义，实质上是逐步求出 y_n\equiv F\pmod{z^{2^n}}。而我们在计算 y_{n+1} 时可以在 {\!}\bmod z^{2^{n+1}} 意义下计算原方程，并且对任意 k\ge 2，F(z^k)\equiv y_n(z^k)\pmod {z^{2^{n+1}}}，所以我们在递推时总是可以把 F(z^i) 当作常数！

这样一来方程就被重新书写为 F\equiv A\exp\left(F+B\right)\pmod{z^{2^{n+1}}}，这就是一个普通的代数方程了，后续步骤略。

5 什么是远处求值？

有时，当我们提取某一幂级数 F(z) 或数列 \lang f_i\rang 的 n 次项系数时，希望得到一个 o(n) 的做法，这被称为远处求值问题，目前较成熟的两类远处求值题型分别为线性递推与整式递推。

线性递推

定义 5.1（常系数齐次线性递推） 我们称数列 a_n 是一个 k 阶常系数齐次线性递推数列（或简称为 k 阶线性递推数列），当且仅当存在 b_1,\dots,b_k（被称为 a 上的线性递推式），满足 \displaystyle \forall n\ge k,a_n=\sum_{i=1}^k a_{n-i}b_i。
定理 5.1 已知 k 阶线性递推数列 a_n 的第 0\sim k-1 项与其上的线性递推式 b_1,\dots,b_k，则计算 a_n 可以做到 \Theta(\mathsf{M}(k)\log n) 的时间复杂度。

过去较为普及的线性递推算法由 Fiduccia 于 1985 年提出，这一算法常数较差，本文将介绍常数更优秀且更易实现的算法，其由 Bostan 和 Mori 于 2020 年提出[^4]，被称为 Bostan-Mori 算法。

该算法分为两部分。首先，我们会希望把线性递推转化为有理函数。

引理 5.1 对上述线性递推数列 a，设 \displaystyle Q(z)=1-\sum_{i=1}^{k}b_iz^i，则存在多项式 P(z) 满足 \deg P<k，使得 \displaystyle\frac{P(z)}{Q(z)}=\sum_{i\ge 0}a_iz^i。

证明：将 Q(z) 移到等式右边，整理得：
Q(z)\sum_{i\ge 0}a_iz^i=\sum_{i\ge 0}z^i\left(a_i-\sum_{j=1}^kb_ja_{i-j}\right)
根据线性递推数列的定义，我们知道 \displaystyle a_i-\sum_{j=1}^kb_ja_{i-j} 在 i\ge k 时恒为 0，因此剩余部分为一个次数小于 k 的多项式，其正是 P(z)。

这样一来，我们就把线性递推数列的远处求值问题转化成了有理函数 \dfrac{P(z)}{Q(z)} 的远处求值问题。这一引理同时还指出，这一转化可以做到 \Theta(\mathsf{M}(k))，不会成为瓶颈。

定理 5.2（Bostan-Mori 算法） 已知次数不超过 k 的多项式 P(z),Q(z)，计算 [z^n]\frac{P(z)}{Q(z)} 可以做到 \Theta(\mathsf{M}(k)\log n) 的时间复杂度。

证明：该算法基于这样的观察：
\frac{P(z)}{Q(z)}=\frac{P(z)Q(-z)}{Q(z)Q(-z)}=\frac{U_0(z^2)+zU_1(z^2)}{V(z^2)}
可以证明 Q(z)Q(-z) 只有偶数次项，所以我们可以将其写成 V(z^2) 的形式。U_0 与 U_1 就是把 P(z)Q(-z) 奇偶分组所得到的两个多项式。

那么我们就可以根据 n 的奇偶性确定 z^n 一项在 \frac{U_0(z^2)}{V(z^2)} 还是 \frac{zU_1(z^2)}{V(z^2)} 中，也就是：
[z^n]\frac{P(z)}{Q(z)}=[z^{\lfloor n/2\rfloor}]\frac{U_{n\bmod 2}(z)}{V(z)}
因为 Q(z)Q(-z) 与 P(z)Q(-z) 次数均不超过 2k，所以 U_0,U_1,V 次数均不超过 k。于是该算法用时满足递归式 T(n,k)=T(\frac{n}{2},k)+\Theta(\mathsf{M}(k))，有 T(n,k)=\Theta(\mathsf{M}(k)\log n)。

虽说 Bostan-Mori 算法的常数优于 Fiduccia 算法，但 Fiduccia 算法也是有其优势的：其算法核心是先计算一个形如 z^n\bmod \Gamma(z) 的式子，再计算一个形如 \lang z^n\bmod \Gamma(z),v_0\rang 的内积。而这里的 \Gamma(z) 只和 a 的线性递推式 b_1,\dots,b_k 有关，v_0 只和 a 的初值 a_0,\dots,a_{k-1} 有关。这允许我们对同一线性递推式预处理出 z^n\bmod\Gamma(z)，以此 \Theta(k) 地计算单组初值所给出的 a_n。

我们希望对 Bostan-Mori 做类似的改造。定义 \displaystyle\mathcal{F}_{n,k}(A)=\sum_{i=0}^{k-1}z^i[z^{n-k+1+i}]A，则对 \deg P<k 的多项式 P，就有 [z^n]\frac{P(z)}{Q(z)}=[z^{k-1}]P(z)\mathcal{F}_{n,k}(\frac{1}{Q(z)})。

定理 5.3（Bostan-Mori 算法变体） 已知次数不超过 k 的多项式 Q(z)，计算 \mathcal{F}_{n,k}(\frac{1}{Q(z)}) 可以做到 \Theta(\mathsf{M}(k)\log n) 的时间复杂度。

证明：仍然考虑上下同乘 Q(-z)，设 V(z^2)=Q(z)Q(-z)。
\begin{aligned} \mathcal{F}_{n,k}\left(\frac{1}{Q(z)}\right) &=\mathcal{F}_{n,k}\left(\frac{Q(-z)}{V(z^2)}\right)\\ &=\mathcal{F}_{n,k}\left(Q(-z)z^{n-2k+1}\mathcal{F}_{n,2k}\left(\frac{1}{V(z^2)}\right)\right)\\ &=\mathcal{F}_{2k-1,k}\left(Q(-z)\mathcal{F}_{n,2k}\left(\frac{1}{V(z^2)}\right)\right)\\ \end{aligned}
可以发现 A(z)=\mathcal{F}_{n,2k}\left(\frac{1}{V(z^2)}\right) 与 B(z)=\mathcal{F}_{\lfloor \frac{n}{2}\rfloor,k}\left(\frac{1}{V(z)}\right) 有相当大的关联，具体的，有 A(z)=z^{(n+1)\bmod 2}B(z^2)。故该算法用时满足递归式 T(n,k)=T(\frac{n}{2},k)+\Theta(\mathsf{M}(k))，同样有 T(n,k)=\Theta(\mathsf{M}(k)\log n)。

对多组不同初值提取系数的一个特例是单组初值提取一段区间内的系数。

推论 5.1 已知次数 <k 的多项式 P(z) 与次数 \le k 的多项式 Q(z)，计算 [z^l]\frac{P(z)}{Q(z)},[z^{l+1}]\frac{P(z)}{Q(z)},\dots,[z^{r}]\frac{P(z)}{Q(z)} 可以做到 \Theta(\mathsf{M}(r-l+k)\log r)。

证明：只需要求出 \mathcal{F}_{r,r-l+k}(\frac{1}{Q(z)})。

接下来我们稍微加强一下问题。

定理 5.4（常系数非齐次线性递推） 已知数列 a_n 的第 0\sim k-1 项、其上的线性递推式 b_1,\dots,b_k 与一个 \deg C=m 的多项式 C(z)，满足 \displaystyle\forall n\ge k,a_n=C(n)+\sum_{i=1}^k a_{n-i}b_i。则计算 a_n 可以做到 \Theta(\mathsf{M}(m+k)(\log n+\log(m+k))) 的时间复杂度。

证明：考虑将其转化成齐次的线性递推。

去掉 C(n) 高次项的一个好方法是差分，此时有：
a_{n+1}=\Delta C(n)+\sum_{i=1}^{k+1} (b_i-b_{i-1})a_{n+1-i}
（假设 b_0=-1，b_{k+1}=0）

那么做 m+1 次差分后有 \Delta^{m+1} C(n)=0，原式变成了 k+m+1 阶的齐次线性递推。

新的问题是我们此时要求出 a_k,\dots,a_{k+m} 作为初值。以生成函数的视角看待线性递推，设 \displaystyle\mathcal{A}(z)=\sum_{i=0}^\infty a_iz^i、\displaystyle\mathcal{C}(z)=\sum_{i=0}^{k-1} a_iz^i+\sum_{i=k}^\infty C(i)z^i、\displaystyle\mathcal{B}(z)=\sum_{i=1}^k b_iz^i，则有：
\begin{aligned} \mathcal{A}(z)&=\mathcal{C}(z)+\mathcal{A}(z)\mathcal{B}(z)\\ &=\frac{\mathcal{C}(z)}{1-\mathcal{B}(z)}\\ \end{aligned}
因为我们只要求 \mathcal{A}(z)\bmod{z^{k+m+1}}，所以只需要知道 \mathcal{C}(z)\bmod{z^{k+m+1}}。用多项式多点求值求出 C(k),\dots,C(k+m) 即可。

整式递推

定义 5.2（整式递推） 我们称数列 a_n 是一个 k 阶整式递推数列，当且仅当存在整式递推式 P_0,\dots,P_k，满足 \forall n\ge k,P_0(n)\ne 0，且 \displaystyle \forall n\ge k,a_n=-\frac{1}{P_0(n)}\sum_{i=1}^k a_{n-i}P_i(n)。
定理 5.5 已知 k 阶整式递推数列 a_n 的第 0,\dots,k-1 项与其上的整式递推式 P_0,\dots,P_k，设 d=\max\limits_{i=0}^k \deg P_i，则计算 a_n 可以做到 \Theta(k^{2.5}\sqrt{nd}\log nd) 或 \Theta(\sqrt{nd}(k^3+k^2\log nd)) 的时间复杂度。

为了证明这一定理，我们需要一些前置知识。
引理 5.2（拉格朗日插值） 已知 n 次多项式 f 的 n+1 个点值 (x_0,y_0),\dots,(x_n,y_n)，则有：
f(x)=\sum_{i=0}^n y_i\prod_{j=0\land j\ne i}^n \frac{x-x_j}{x_i-x_j}
证明略。
引理 5.3 已知 n 次多项式 f 的 n+1 个点值 y_0=f(0),\dots,y_n=f(n)，已知正整数 m>n，则计算 f(m),\dots,f(m+n) 可以做到 \Theta(\log p+\mathsf{M}(n)) 的时间复杂度，此处 \log p 的来源是 \mathbb{F}_p 下的逆元。

证明：考虑直接应用引理 5.2，则有：
\begin{aligned} f(m+k) &=\sum_{i=0}^n y_i\prod_{j=0\land j\ne i}^n\frac{m+k-j}{i-j}\\ &=\sum_{i=0}^n y_i\frac{(m+k)^{\underline{n+1}}}{(m+k-i)i!(n-i)!(-1)^{n-i}}\\ &=\sum_{i=0}^n \frac{(m+k)^{\underline{n+1}}}{(m+k-i)}\frac{y_i}{i!(n-i)!(-1)^{n-i}}\\ \end{aligned}
设 \displaystyle F(z)=\sum_{i=0}^n \frac{y_i}{i!(n-i)!(-1)^{n-i}}z^i、\displaystyle G(z)=\sum_{i=m-n}^{m+n} \frac{z^i}{i}，则 f(m+k)=(m+k)^{\underline{n+1}}[z^{m+k}]F(z)\times G(z)。m^{\underline{n+1}},\dots,(m+n)^{\underline{n+1}} 不难 \Theta(n+\log p) 计算。
定理 5.5 的证明（第一部分） 设 \displaystyle b_{n,m}=a_n\prod_{i=k}^{m} P_0(i)，则 b 满足：
- 对任意 m>n\ge k，b_{n,m}=b_{n,m-1}P_0(m)
- 对任意 n\ge k，\displaystyle b_{n,n}=-\sum_{i=1}^k b_{n-i,n-1}P_i(n)。
这可以使用矩阵的形式刻画：
\begin{bmatrix} -P_1(n)&-P_2(n)&-P_3(n)&\cdots&-P_{k-1}(n)&-P_k(n)\\ P_0(n)&0&0&\cdots&0&0\\ 0&P_0(n)&0&\cdots&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&\cdots&P_0(n)&0 \end{bmatrix} \begin{bmatrix} b_{n-1,n-1}\\ b_{n-2,n-1}\\ b_{n-3,n-1}\\ \vdots\\ b_{n-k,n-1} \end{bmatrix} = \begin{bmatrix} b_{n,n}\\ b_{n-1,n}\\ b_{n-2,n}\\ \vdots\\ b_{n-k+1,n} \end{bmatrix}
于是设 M_n 为上面给出的转移矩阵，\bm{v_0}=\begin{bmatrix}a_{k-1}\\\vdots\\ a_0\end{bmatrix}，则 M_n\cdots M_{k}\bm{v_0} 的第一项即为 b_{n,n}，随后我们只需要再求出 \displaystyle\prod_{i=k}^n P_0(i) 即可从 b_{n,n} 还原 a_n。

既然如此，我们不妨先考察一下如何计算 \displaystyle\prod_{i=k}^n P_0(i)，并期望能从做法中得到原问题的启示。
引理 5.4 已知正整数 n 与 \deg P=d 的多项式 P，则计算 \displaystyle\prod_{i=0}^{n-1} P(i) 可以做到 \Theta(\sqrt{nd}\log nd) 的时间复杂度。

证明：考虑分块，设 \displaystyle F_t(z)=\prod_{i=0}^{t-1} P(z+i)，则我们只需要设置一个参数 B，并求出 h_{B,0},\dots,h_{B,\lfloor\frac{n}{B}\rfloor-1} 满足 h_{i,j}=F_i(Bj)，则答案即为 \displaystyle\prod_{i=0}^{\lfloor\frac{n}{B}\rfloor-1}h_{B,i}\prod_{i=\lfloor\frac{n}{B}\rfloor B}^{n-1}P(i)，后一乘积可以做到 \Theta(Bd)，我们重点关注前一乘积。

考虑倍增地计算 h，这一想法的动机是注意到 F_{2i}(j)=F_i(j)F_i(i+j)，因此 h_{2i,j}=h_{i,j}h_{i,j+\Delta}，其中 \Delta=\frac{i}{B}（下标可以不是整数）。这给了我们应用引理 5.3 的可能。

另一方面，h_{i,j} 可以看成是多项式 H_i(z)=F_i(Bz)（该多项式的次数不超过 id）的点值，所以这确实满足引理 5.3 的应用条件。我们便知道：
2. 如果 h_{i,0},\dots,h_{i,id} 已知，则根据引理 5.3，我们可以在 \Theta(id\log id)（忽略 \log p）的时间内计算出 h_{i,id+1},\dots,h_{i,2id}；
3. 如果 h_{i,0},\dots,h_{i,2id} 已知，根据引理 5.3，我们可以在 \Theta(id\log id) 的时间内计算出 h_{i,\Delta},\dots,h_{i,2id+\Delta}；
4. 如果上述均已知，则我们可以在 \Theta(id) 的时间内计算出 h_{2i,0},\dots,h_{2i,2id}。
为了精确地得到 h_{B,j}，我们还需要支持从 h_{i,j} 得到 h_{i+1,j}。而因为 h_{i+1,j}=h_{i,j}h_{1,j+\Delta}，所以我们同样可以先 \Theta(id\log id) 地求出 h_{i,id+1},\dots,h_{i,(i+1)d} 与 h_{1,\Delta},\dots,h_{1,(i+1)d+\Delta}，再由它们推出 h_{i+1,0},\dots,h_{i+1,(i+1)d}。

求出 h_{m,0},\dots,h_{m,md} 的用时 T(m) 满足 T(m)=T(\frac{m}{2})+\Theta(md\log md)，解得 T(m)=\Theta(md\log md)。同时最后我们还要把 h_{B,0},\dots,h_{B,Bd} 扩张到 \lfloor\frac{n}{B}\rfloor-1，于是总用时为 \Theta((\frac{n}{B}+Bd)\log Bd+\frac{n}{B})，取 B=\sqrt{\frac{n}{d}} 平衡得 \Theta(\sqrt{nd}\log nd)。
定理 5.5 的证明（第二部分） 我们发现引理 5.4 的证明过程几乎可以原封不动地用来计算 M_n\cdots M_k。具体过程高度相似，此处留做习题。

用时 T(m) 满足 T(m)=T(\frac{m}{2})+\Theta(md(k^3+k^2\log md))，解得 T(m)=\Theta(md(k^3+k^3\log md))，总用时为 \Theta(Bdk^3+(Bd+\frac{n}{B})k^2\log Bd)，取 B=\sqrt{\frac{n}{dk}} 平衡得 \Theta(k^{2.5}\sqrt{nd}\log nd)，也可以选择取 B=\sqrt{\frac{n}{d}} 平衡得 \Theta(\sqrt{nd}(k^3+k^2\log nd))。

（分析的有点混乱，如果有问题可以联系我。）

6 你说得对，但是什么是远处求值？

上一章中，我们讨论了线性递推与整式递推提取远处系数的方法，但求出递推式本身也是一个难题，本章将着重探讨这方面的内容，并揭示两种递推背后更深层的结构。

线性递推与有理函数重建

碎碎念：在撰写本节时，作者不明不白地跑去学了一下多项式上的逼近理论，然后发现其实并不需要，，，你可以在这里看到被我删掉的这部分内容。

根据引理 5.1，任意线性递推的结构都可以由有理函数 \frac{P(z)}{Q(z)} 刻画。我们自然会好奇一个问题（这被称为有理函数重建问题）：如果给定数列的前若干项 \displaystyle A(z)=\sum_{i=0}^{n-1} a_iz^i，能不能求出一个 \deg P(z)<\deg Q(z) 且 \deg Q(z) 最小的有理函数 \frac{P(z)}{Q(z)}，使得 A(z)\equiv\frac{P(z)}{Q(z)}\pmod{z^n}？

我们先将问题适当泛化：

定义 6.1（有理函数重建、有理逼近） 已知多项式 F(z),M(z) 满足 \deg F<\deg M，已知度数限制 d_P,d_Q。若存在多项式 P(z),Q(z) 满足 P\equiv QF\pmod M 且 \deg P\le d_P、\deg Q\le d_Q、Q\ne 0，则我们称 (P,Q) 是 F 的一组 (d_P,d_Q) 有理逼近。

那么还原线性递推式相当于在找到 M=z^n 时 A 的一组有理逼近（d_P 与 d_Q 的取值我们先按下不表），要求 \deg P<\deg Q 且 \deg Q 尽量小。

我们可以把这一同余方程改写成 P=QF+RM，这启示我们将其刻画为 Euclid 过程：

令 R_{-1}=M、R_{0}=F；
对所有 i\ge 0，令：
重复上述过程，直到首次出现 R_i=0。

容易在这一过程中同时维护 S_i 与 T_i，使得对任意 i 都有 R_i=S_iF+T_iM，它们满足与渐近分数 \frac{P_i}{Q_i} 类似的递推关系：

对所有 i\ge 0，S_{i+1}=S_{i-1}-A_iS_{i}、T_{i+1}=T_{i-1}-A_iT_{i}。

经过一些观察，可以发现：

引理 6.1 对任意 i\ge 0，\deg R_{i-1}+\deg S_i=\deg M。

证明：因为 A_i 是 R_{i-1} 除以 R_{i} 的商，且其不等于 0，所以 \deg R_{i-1}=\deg R_{i}+\deg A_i。而另一侧，我们不难根据 \deg A_i>0 归纳证明 \deg S_i<\deg S_{i+1}，所以也有 \deg S_{i+1}=\deg S_i+\deg A_i。

因此 \deg R_{i}+\deg S_{i+1}=(\deg R_{i-1}-\deg A_i)+(\deg S_i+\deg A_i)=\deg R_{i-1}+\deg S_i，而初始时 \deg R_{-1}+\deg S_0=\deg M，故引理得证。

而 R_i=S_iF+T_iM 与 P=QF+RM 中，左式的 R_i 与 S_i 扮演的角色与右式的 P 与 Q 一样，所以实际上我们可以直接根据这一点来找有理逼近。

定理 6.1（有理逼近的存在性） 如果 d_P+d_Q\ge \deg M-1，则 (P,Q) 一定存在。

证明：我们找到首个满足 \deg R_i\le d_P 的 R_i，则有 d_P<\deg R_{i-1}，根据引理 6.1，\deg S_i=\deg M-\deg R_{i-1}\le\deg M-1-d_P\le d_Q，故 (R_i,S_i) 为一组满足条件的有理逼近。
定理 6.2（有理逼近的唯一性） 如果 d_P+d_Q<\deg M，则 (P,Q) 一定唯一。

证明：若存在两组有理逼近 (P_1,Q_1) 与 (P_2,Q_2)，说明 \frac{P_1}{Q_1}\equiv\frac{P_2}{Q_2}\pmod M，交叉相乘后因为 \deg P_1+\deg Q_2\le d_P+d_Q<\deg M 且右侧也是如此，故 P_1Q_2=P_2Q_1，得证。

定理 6.1 与定理 6.2 实际上保证了 d_P+d_Q=\deg M-1 的有理逼近是唯一存在的，所以我们可以只关心这些解。而回到线性递推，我们还要求 \deg P<\deg Q 且 \deg Q 尽量小，而定理 6.1 实际上保证了每一组 (R_i,S_i) 都是 d_P+d_Q=\deg M-1 时的一组有理逼近（且是唯一一组），找到第一组满足 \deg R_i<\deg S_i 的即可。

上半部分中我们将有理函数重建问题完全解决，而实现上我们会需要支持快速对一组 (F,M) 运行 Euclid 过程（下面设 n=\deg M）。因为计算 \lfloor\frac{F}{G}\rfloor 与 F\bmod G 可以在 \Theta((\deg F-\deg G)\deg G) 的时间内完成，因此朴素地进行模拟就能做到 \Theta(\sum (\deg R_{i-1}-\deg R_i)\deg R_i)=\Theta(n^2) 的时间。部分场景已足够使用，不过我们确实能在给力一点，因为有 Half-GCD 算法：

定理 6.3（Half-GCD 算法） 存在算法能在 \Theta(n\log^2 n) 的时间内计算两个多项式 A,B 的 gcd。

下面将介绍这一算法，并给出其正确性证明。

gcd 就是 R 的最后一项非零项，于是我们希望找到一种方法快速求出 R 的值。可以发现 R 上的递推关系可以使用矩阵刻画为 \begin{bmatrix}R_i\\ R_{i+1}\end{bmatrix}=\begin{bmatrix}0&1\\ 1&-Q_i\end{bmatrix}\begin{bmatrix}R_{i-1}\\ R_i\end{bmatrix}，此处 Q_i=\lfloor\frac{R_{i-1}}{R_i}\rfloor。我们记 \lang Q\rang=\begin{bmatrix}0&1\\ 1&-Q\end{bmatrix}，则 \begin{bmatrix}R_k\\ R_{k+1}\end{bmatrix}=\lang Q_k\rang\cdots\lang Q_0\rang\begin{bmatrix}R_{-1}\\ R_0\end{bmatrix}。

而我们知道 \begin{bmatrix}R_{-1}\\ R_0\end{bmatrix}=\begin{bmatrix}A\\ B\end{bmatrix}，现在问题转化为求出 \lang Q_k\rang\cdots\lang Q_0\rang，这里的 Q_i 依赖于 \lang Q_{i-1}\rang\cdots\lang Q_0\rang，所以我们的算法应该类似于一个在线算法。

在线算法的流程通常都和分治有关，但这里朴素算法的用时瓶颈在于矩阵系数的度数可能过大，并且我们还不知道 k 的具体值，所以基于下标的分治方式并不可行，我们需要从度数切入去进行分治。

假设我们已经实现了这样一个叫做 \operatorname{hgcd}(A,B) 的算法，它输入两个多项式 A,B 满足 \deg A>\deg B\ge 0，并输出 A 与 B 之间的 Euclid 过程中的一步 \mathbf{M}=\lang Q_k\rang\cdots\lang Q_0\rang，使得令 \begin{bmatrix}C\\ D\end{bmatrix}=\mathbf{M}\begin{bmatrix}A\\ B\end{bmatrix} 后有 \deg C\ge\lceil\frac{\deg A}{2}\rceil>\deg D。显然 \operatorname{hgcd} 的输出是唯一的。

那么在计算 \gcd(A,B) 时，我们只要再求出 Q_{k+1}=\lfloor\frac{C}{D}\rfloor 并将 \begin{bmatrix}C\\ D\end{bmatrix} 左乘上 \lang Q_{k+1}\rang，则新的向量就满足两维的度数均 <\lceil\frac{\deg A}{2}\rceil，也就可以递归下去了。于是接下来我们只需要考虑 \operatorname{hgcd} 的设计。

（不直接使得 \lceil\frac{\deg A}{2}\rceil>\deg C>\deg D 是为了方便 \operatorname{hgcd} 算法自身的设计。）

剖析一下带余除法的结构，设 Q=\lfloor\frac{A}{B}\rfloor，我们发现当 \deg A\ge \deg B 时总有 \deg Q=\deg A-\deg B，而计算 Q 的 \deg Q-i 次项时我们只要知道 A 的 \deg A-i 次项与 B 的最高次项，而为了得到正确的 A 的 [\deg A-\deg Q,\deg A] 次项，我们还要知道 B 的 [\deg B-\deg Q,\deg B] 次项。这就是我们为了计算 Q 需要知道的所有信息了！

这启示我们在进行 Euclid 过程时，前半部分可能不会涉及到 A 与 B 的更低项，于是我们可以分治。基于这一启示，我们就可以设计出 \operatorname{hgcd} 算法的一个大致框架了：

设 m=\lceil\frac{\deg A}{2}\rceil；
若 \deg B<m，则单位矩阵就满足 \operatorname{hgcd} 的输出要求，直接返回即可；
否则，我们将 A 和 B 分别按 m 次项为界划分为两个多项式 A_0z^m+A_1 与 B_0z^m+B_1，满足 \deg A_1,\deg B_1<m；
只考虑 A_0 与 B_0，在其上运行 \operatorname{hgcd} 算法，我们猜测这一步对低次项的影响是极其有限的，也就是仍有 \deg A=\deg A_0+m、\deg B=\deg B_0+m；
此时 \deg B_0<\frac{\deg A_0}{2}\le\frac{\deg A}{4}，于是辗转相除一次后 A' 和 B 的次数都小于 \frac{3}{4}\deg A（忽略一些取整）；
再取 m 为 \frac{\deg A}{4} 重新划分 A'_0,A'_1,B_0,B_1，并重新对 A'_0,B_0 运行 \operatorname{hgcd} 算法；
此时 \deg B_0<\frac{\deg A'_0}{2}\le\frac{\deg A}{4}，于是 \deg B<\frac{\deg A}{2}，成功满足了输出要求。

这一算法用时 T(n) 满足递归式 T(n)=2T(\frac{n}{2})+\Theta(\mathsf{M}(n))，也就有 T(n)=\Theta(\mathsf{M}(n)\log n)。接下来我们补全一些算法和正确性证明的细节。

定义 6.2（正则矩阵） 若一个矩阵 \mathbf{M} 可以表示成 \lang Q_1\rang\lang Q_2\rang\cdots 的形式，我们就称其为正则（regular）的。

回忆一下，\lang Q\rang=\begin{bmatrix}0&1\\ 1&-Q\end{bmatrix}。

观察 6.1 若矩阵 \mathbf{M} 是正则的，则 \det\mathbf{M}=\pm 1；
观察 6.2 在 Euclid 过程中用到的 \lang Q\rang 总是满足 \deg Q\ge 1；
观察 6.3 若矩阵 \mathbf{M}=\begin{bmatrix}P&Q\\ R&S\end{bmatrix} 是正则的且不是单位矩阵，则 \deg P\le \min\{\deg Q,\deg R\}\le\max\{\deg Q,\deg R\}\le\deg S，且 \deg P<\deg S；
观察 6.4 若矩阵 \mathbf{M}=\begin{bmatrix}P&Q\\ R&S\end{bmatrix} 是正则的，则 \mathbf{M}^{-1}=\pm\begin{bmatrix}S&-Q\\ -R&P\end{bmatrix}，此处的正负号由 \deg\mathbf{M} 唯一确定。

确认一下我们要验证的事情：只需要验证算法的第四步中的猜测是对的，那么 A_0 和 B_0 上的 \operatorname{hgcd} 算法所计算的商的次数就始终不会超过 \frac{\deg A}{4}，于是也就不会受到低次项的影响，可以直接作用到 A 和 B 上。

引理 6.2 已知多项式 A,B、阈值 m 与正则矩阵 \mathbf{M}，设 A=A_0z^m+A_1、B=B_0z^m+B_1 满足 \deg A_1,\deg B_1<m。

我们定义 \mathbf{M}\begin{bmatrix}A_0\\ B_0\end{bmatrix}=\begin{bmatrix}A'_0\\ B'_0\end{bmatrix}、\mathbf{M}\begin{bmatrix}A_1\\ B_1\end{bmatrix}=\begin{bmatrix}A'_1\\ B'_1\end{bmatrix}，A'=A'_0z^m+A'_1、B'=B'_0z^m+B'_1。若 \deg A_0'>\deg B_0' 且 \deg A_0'\ge\frac{\deg A}{2}（这是我们在执行 \operatorname{hgcd}(A_0,B_0) 时能保证的条件），则 \deg A'=\deg A'_0+m 且 \deg B'\le \max\{\deg B'_0,\deg A_0-\deg A'_0-1\}+m.

证明：设 \mathbf{M}=\begin{bmatrix}P&Q\\ R&S\end{bmatrix}，则根据观察 6.4，\mathbf{M}^{-1}=\pm\begin{bmatrix}S&-Q\\ -R&P\end{bmatrix}，于是 A_0=\pm(A'_0S-B'_0Q)，同时 \deg A'_0>\deg B'_0 且 \deg S\ge \deg Q（根据观察 3.3），故 \deg A_0=\deg A'_0+\deg S。再结合 \deg A'_0\ge\frac{\deg A}{2}，有 \deg S\le\frac{\deg A}{2}\le\deg A'_0。

再看 \deg A'_1，有 A'_1=A_1P+B_1Q，于是有：
\begin{aligned} \deg A'_1&\le \max\{\deg A_1+\deg P,\deg B_1+\deg Q\}\\ &\le\max\{\deg A_1,\deg B_1\}+\max\{\deg P,\deg Q\}\\ &<m+\deg S\\ &\le m+\deg A'_0 \end{aligned}
因此 \deg A'_0z^m>\deg A'_1，前半句话得证。

对于后半句话，我们也类似地分析一下，但是就不需要这么精细了。有 B'_1=A_1R+B_1S，因此：
\begin{aligned} \deg B'_1&\le \max\{\deg A_1+\deg R,\deg B_1+\deg S\}\\ &\le m-1+\deg S\\ &=m-1+\deg A_0-\deg A'_0 \end{aligned}
于是有 \deg B'\le\max\{\deg B'_0,\deg A_0-\deg A'_0-1\}+m。

补全一下算法细节（框架中唯一一处需要修正的是第八行）：

设 m=\lceil\frac{\deg A}{2}\rceil；
若 \deg B<m，直接返回单位矩阵；
否则按 m 次项划分 A 与 B。设 \mathbf{R}=\operatorname{hgcd}(A_0,B_0)；
设 \begin{bmatrix}C\\ D\end{bmatrix}=\mathbf{R}\begin{bmatrix}A\\ B\end{bmatrix}；
若 \deg D<m，直接返回 \mathbf{R}；
否则设 Q=\lfloor\frac{C}{D}\rfloor、E=C\bmod D；
若 \deg E<m，直接返回 \lang Q\rang\mathbf{R}；
否则设 \bm {k=2m-\deg D}，并依此划分 \bm D 与 \bm E；
设 \mathbf{S}=\operatorname{hgcd}(D_0,E_0)，返回 \mathbf{S}\lang Q\rang\mathbf{R}。

定理 6.3' 上述算法是正确的。

证明：我们只需证明两部分：返回的矩阵确实是 A 和 B 之间的 Euclid 过程中的某一步，且作用到 \begin{bmatrix}A\\ B\end{bmatrix} 上后满足度数限制。

第二行的返回显然正确。

对于第五行的返回，我们先证明 \mathbf{R} 里的辗转相除作用到 A,B 上仍然正确。显然 \mathbf{R} 里每一项的度数都不会超过 \lfloor\frac{\deg A_0}{2}\rfloor，同时根据引理 6.2，A 与 B 之间的最后一次辗转相除前的 \deg A>\deg B=\deg C=\deg C_0+m\ge\lceil\frac{\deg A_0}{2}\rceil+m，于是 \mathbf{R} 在作用到其上后不会涉及到 A_1 与 B_1，我们便验证了这一部分的正确性。度数限制不难验证，此处略。

第七行的返回显然正确。

对于第九行的返回，首先我们有 \deg D\le \max\{\deg D_0,\deg A_0-\deg C_0-1\}+m\le\lceil\frac{\deg A_0}{2}\rceil-1+m，因此 k=2m-\deg D\ge m-\lceil\frac{\deg A_0}{2}\rceil-1\ge 0。接下来与第五行同理，我们发现 \mathbf{S} 作用到 D 和 E 上后最后一次辗转相除前有 \deg D>\deg E=\deg D'=\deg D'_0+k\ge\lceil\frac{\deg D_0}{2}\rceil+k，而 \mathbf{S} 里系数的度数不超过 \lfloor\frac{\deg D_0}{2}\rfloor，于是 \mathbf{S} 是正确的作用。

对于度数限制，我们知道 \deg D_0=\deg D-k=2\deg D-2m，因此 \deg D'\ge \deg D-m+k=m；而 \deg E'\le \max\{\deg E'_0,\deg D_0-\deg D'_0-1\}+k\le\deg D-1-m+k=m-1，于是 (D',E') 确实是满足度数限制的返回值。故算法正确性得证。

有了强大的 Half-GCD 算法，我们就不难计算下列值了，细节留作习题：

对多项式 A,B 计算 \gcd(A,B)；
对多项式 A,B 计算多项式 P,Q 满足方程 AP+BQ=\gcd(A,B)；
有理函数重建。

整式递推与 D-Finite

与线性递推不同，整式递推强大之处在于其表达力，以至于相当多的数列都能用整式递推刻画，但另一方面，我们又难以得到整式递推的一个简短的封闭形式。于是我们并不会那么关心怎么用有限项数列“还原”其整式递推式，而是关心怎么从一个完整的表达式来“找到”其整式递推式。

我们针对数列的工具并不多，一个明智的想法是搬到形式幂级数环上工作。

定义 6.3（微分有限） 已知幂级数 F(z)，若存在一系列多项式 C_0,\dots,C_k 使得 F 满足：
\sum_{i=0}^k C_iF^{(i)}=0
（此处 F^{(i)} 意为 i 阶导数），则我们称 F 是微分有限（D-Finite）的。

这实际上在说 F 满足某个线性常微分方程。可以发现整式递推与微分有限是一体两面：

定理 6.4 数列 \lang f_i\rang 是整式递推数列，当且仅当其生成函数 F 是微分有限的。

证明：从 F 微分有限推 f 整式递推是显然的，提取系数即可。而如果 f 整式递推，我们就可以用线性算子构造 F 的一个微分方程。

记算子 \vartheta=\frac{z\mathrm{d}}{\mathrm{d}z}，其作用到 F 上后系数会发生 f_n\gets nf_n 的变化。那么对于 f 的整式递推的一项系数 \displaystyle\sum_n f_{n-i}(n-i)^jz^n（我们先把所有 P_i(n) 改写为 P_i(n-i) 的形式），其就可以表示成 z^i\vartheta^j F=z^{i+j}(\frac{\mathrm{d}}{\mathrm{d}z})^j F。因此我们可以将完整的递推式改写成一个关于 F 的线性常微分方程，故 F 微分有限。

值得注意的是 \lang f_i\rang 的初值信息并没有蕴含在 C_i 中，也就是说一组 C_i 可以对应多个 F，需要人为限定初值，这些初值也就是 [z^0]F,[z^0]F',\dots,[z^0]F^{(k-1)}。

接下来我们就可以讨论更易操作的微分有限函数了。

为了判定 F 是不是微分有限的，我们将工作分为两部分：一些常见的微分有限函数与其上可以出现的构造方法。

很神奇的是常见的微分有限函数其实只有两种：一种是幂函数 F=z^k，其满足微分方程 kF-zF'=0。

而另一种就是超几何函数 \displaystyle F=F\left(\begin{matrix}a_1,\dots,a_m\\ b_1,\dots,b_n\end{matrix}\middle| z\right)=\sum_{i\ge 0}\frac{a_1^{\overline{i}}\dots a_m^{\overline{i}}}{b_1^{\overline{i}}\dots b_n^{\overline{i}}}\frac{z^i}{i!}。不难验证它满足微分方程：

\left(\frac{\mathrm{d}}{\mathrm{d}z}\prod_{i=1}^n (\vartheta+b_i-1)\right)F=\left(\prod_{i=1}^m (\vartheta+a_i)\right)F

证明：左式在求导前等于 \displaystyle\sum_{i\ge 0}\frac{a_1^{\overline{i}}\dots a_m^{\overline{i}}}{b_1^{\overline{i-1}}\dots b_n^{\overline{i-1}}}\frac{z^i}{i!}，右式等于 \displaystyle\sum_{i\ge 0}\frac{a_1^{\overline{i+1}}\dots a_m^{\overline{i+1}}}{b_1^{\overline{i}}\dots b_n^{\overline{i}}}\frac{z^i}{i!}。

而很多常见幂级数实际上都只是超几何函数的一个特例：

\displaystyle F\left(\begin{matrix}1,1\\2\end{matrix}\middle|{-z}\right)=z^{-1}\ln(1+z)

更一般地，所有常数项为 1 且相邻项系数比值为“关于次数的有理函数”的幂级数均为超几何函数，证明略。

接下来我们再看构造方法，为方便验证一种构造方法是否真的能生成微分有限的幂级数，我们有引理：

引理 6.3 形式幂级数 F 微分有限当且仅当 F,F',F^{(2)},\dots 在域 \mathbb{A}(z) 上所张成的空间维数有限。

证明：若前者成立。则我们不难根据定义式将任意 F^{(n)} 表示为 F,F',\dots,F^{(k)} 的线性组合，因此 F,F',\dots 所张成的空间维数不超过 k+1。

若后者成立，则说明存在一个 k（且这个 k 可以取到维数）使得 F,F',\dots,F^{(k)} 线性相关，这可以直接导出 F 微分有限。

为什么是 \mathbb{A}(z) 而不是等价的 \mathbb{A}[z]（只需要乘上/除以所有分母的 \operatorname{lcm}）？这是因为前者是域而后者只是环，我们自然会希望在性质更好的结构上工作。

也可以用数列的语言写出与其等价的命题。

引理 6.3' 定义作用在数列上的位移算子 S 满足 (Sf)_n=f_{n+1}。则数列 \lang f_i\rang 可整式递推当且仅当由 f,Sf,S^2f,\dots 在域 \mathbb{A}(n) 上所张成的空间维数有限。

值得注意的是，你不应该将数列看成一串数 f_0,f_1,\dots，而应该看成一个以下标 n 为自变量的函数 f(n)，这样一来引理 6.3' 的表述应该就很自然了。

定理 6.5（加法） 若幂级数 F,G 微分有限，则 H=F+G 微分有限。

证明：设 F,\dots,F^{(n-1)} 线性相关，G,\dots,G^{(m-1)} 线性相关（以下将沿用这一假设），因此任意的 H^{(k)} 都可以表示成上述 n+m 项的线性组合，于是 H,\dots,H^{(n+m-1)} 必然线性相关。
定理 6.6（乘法） 若幂级数 F,G 微分有限，则 H=FG 微分有限。

证明：不难发现 H^{(k)} 总是可以表示成 F^{(i)}G^{(j)} 的线性组合，于是 H,\dots,H^{(nm-1)} 必然线性相关。
定理 6.7（点积） 若数列 \lang f_i\rang,\lang g_i\rang 可整式递推，则数列 \lang h_i\rang 满足 h_i=f_ig_i 也可整式递推。

证明：不难发现 S^kh(n) 总是可以表示成 \{S^if(n)\cdot S^jg(n)\mid i<n,j<m\} 内元素的线性组合，这和定理 6.6 是一样的。
定理 6.8（复合） 若幂级数 F 微分有限，G 代数（指的是，存在多项式 P(z,t)\ne 0 使得 P(z,G)=0），则 H=F\circ G 微分有限。

我们将证明分为两部分：第一部分是证明将域 \mathbb{A}(z) 扩张到 \mathbb{A}(z,G) 后 H,H',\dots 张成的空间维数有限，第二部分是证明域 \mathbb{A}(z,G) 上的有限维空间在 \mathbb{A}(z) 中仍然是有限维的。
引理 6.4 G'\in\mathbb{A}(z,G)。

证明：不妨设 P(z,t) 是 G 的最小多项式。

我们已知 P(z,G)=0，对 z 一元求导可得：
\frac{\mathrm{d}P(z,G(z))}{\mathrm{d}{z}}=\frac{\partial P}{\partial z}(z,G)+\frac{\partial P}{\partial t}(z,G)G'=0
于是有 G'=-\dfrac{\frac{\partial P}{\partial z}(z,G)}{\frac{\partial P}{\partial t}(z,G)}，分子与分母均为 \mathbb{A}[z,G] 中的多项式，且因为 P 最小，所以 \frac{\partial P}{\partial t} 不可能以 G 为根，因此分母不为零，故 G'\in\mathbb{A}(z,G)。
引理 6.5 对任意幂级数 H，若 H\in\mathbb{A}(z,G)，则 H'\in\mathbb{A}(z,G)。

证明：设 H=\dfrac{A(z,G)}{B(z,G)}，其中 A,B\in\mathbb{A}[z,t]。则我们对 H 求导可得：
\begin{aligned} H'&=\left(\frac{A(z,G)}{B(z,G)}\right)'\\ &=\frac{\frac{\mathrm{d}A(z,G)}{\mathrm{d}z}B(z,G)-\frac{\mathrm{d}B(z,G)}{\mathrm{d}z}A(z,G)}{B^2(z,G)}\\ \end{aligned}
而 \dfrac{\mathrm{d}{A(z,G)}}{\mathrm{d}z} 展开后和引理 6.4 中的形式一样，因为我们已知 G'\in\mathbb{A}(z,G)，所以 \dfrac{\mathrm{d}{A(z,G)}}{\mathrm{d}z}\in\mathbb{A}(z,G)，故 H'\in\mathbb{A}(z,G)。
定理 6.8（第一部分） H,H',\dots 在 \mathbb{A}(z,G) 上张成的空间维数有限。

证明：根据链式求导法则，我们可以把 H^{(k)} 展开成如下形式：
H^{(k)}=\sum_i (F^{(i)}\circ G)B_{k,i}(G,G',\dots)
其中 B_{k,i} 是多项式。根据引理 6.5，我们知道 B_{k,i}(G,G',\dots) 一定 \in\mathbb{A}(z,G)，因此接下来我们只需要考察 F^{(i)}\circ G。

而根据 F 的定义，我们知道：
\sum_{i=0}^{n-1}(C_i\circ G)\cdot(F^{(i)}\circ G)=0
显然 C_i\circ G\in\mathbb{A}(z,G)，所以 F\circ G,\dots,F^{(n-1)}\circ G 在 \mathbb{A}(z,G) 上线性相关，第一部分得证。
引理 6.6 [\mathbb{A}(z,G):\mathbb{A}(z)] 有限，也就是说 \mathbb{A}(z,G) 在 \mathbb{A}(z) 上的维数是有限的。

证明：设 d=\deg_t P，则直接有 1,G,\dots,G^{d} 在 \mathbb{A}(z) 上线性相关，故引理得证。
定理 6.8（第二部分） H 微分有限。

证明：结合定理 6.8.1 与引理 6.6 即证，这是因为域的扩张次数满足乘法原理。

这就给了我们构造性地得到某个幂级数的 ODE 的能力。当然手算时还有一种可行的方法是直接大力求导并对比系数，可参见我的一篇魔怔文。

不过一个可能更 pratical 的做法是在分析出 F 能整式递推后直接用其他方式求出 F 的前若干项，这样我们就能得到关于 C_i 系数的若干条方程，直接求解即可。因为整式递推式通常都不会很复杂，所以这个做法在用时上是可接受的。

[^1]: Hoeven, J. Faster relaxed multiplication. (HAL,2012), http://hal.archives-ouvertes.fr/hal-00687479

[^2]: Hoeven, J. New algorithms for relaxed multiplication. https://www.texmacs.org/joris/newrelax/newrelax.html#dichrelax-fig

[^3]: Kinoshita, Y. \& Li, B. Power Series Composition in Near-Linear Time. (2024), https://arxiv.org/abs/2404.05177

[^4]: Bostan, A. \& Mori, R. A Simple and Fast Algorithm for Computing the N-th Term of a Linearly Recurrent Sequence. CoRR. abs/2008.08822 (2020), https://arxiv.org/abs/2008.08822

[^5]: 李白天, 信息学竞赛中的生成函数计算理论框架. IOI2021 国家集训队论文集 (2021)