数学杂谈 #21

初探整式递推和 ODE

约定

记 \(D=\frac{\text d}{\text dx},\vartheta=xD\)​​​。前者是微分算子，后者可以看作对于形式幂级数，将每一项的指数乘到它的系数里去。

记 \(f^{(k)}=D^kf\)。

ODE

ODE 和 D-finite

ODE 的意思是“常微分方程（Ordinary Differential Equation）”。

Note.

查了一下百科，Ordinary 指的是“只有一个变元”，也就是未知函数是一元函数的微分方程。

它的一般形式是：

\[F(x,y,y',y'',\dots,y^{(n)})=0 \]

其中 \(n\) 是某个正整数。（据说）OI 中常用的形式为：

\[\sum_{k=0}^np_k(x)y^{(k)}(x)=0 \]

（据说）\(p_k(x)\)​​ 一般是常数次数或者常数项数的有理分式，\(n\)​ 一般是常数。特别地，我们要求 \(p_n(x)\neq 0\)，否则这个微分方程不是阶数最小的。

对于一个函数 \(f(x)\)​，若它满足一个上述常微分方程，那么就称它为“微分有限（D-finite）”的。

常见的 D-finite 的函数

有限多项式

如果 \(A(x)=\sum_{k=0}^ma_kx^k\)​，显然有 \(A^{(m+1)}=0\)​。

这个东西没有 \(\Delta^{m+1}A=0\) 好用 qwq。

当然，如果差分对于下降幂来说是平凡的位移效果，而在普通幂上面可以建立方程，我们也许可以考虑微分作用在下降幂上的结果。

幂函数

如果 \(f(x)=x^a,a\neq 0\)​​​，显然有 \(af-xf'=0\)​。

指数函数

如果 \(f(x)=e^x\)​​​，显然有​​ \(f-f'=0\)​。其它的 \(f=a^x\)​ 可以通过复合来表达。

对数函数

如果 \(f(x)=\ln x\)，显然有 \(xf'-1=0\)。其它的 \(f=\log_ax\) 可以通过复合来表达。

注：虽然这个东西不是一般的线性组合形式，但是它仍然满足一般的 ODE 的形式。

三角函数和反三角函数

三角的 ODE 很简单，但是没啥用，不写了。

反三角的 ODE 可以通过查导数得到，并且更加没用，也不写了。

超几何级数

注：以下内容均来自于《具体数学》和 slide.pdf，我都是抄的。

超几何级数的定义是：

\[F\left(\left.\begin{matrix}a_1,a_2,\dots,a_m\\b_1,b_2,\dots,b_n\end{matrix}\right|z\right)=\sum_{k\ge 0}\frac{z^k\prod_{j=1}^ma_j^{\overline{k}}}{k!\prod_{t=1}^nb_t^{\overline{k}}} \]

其中任意的 \(b\)​ 都不可以是非正整数。

根据 slide.pdf，我们可以得到超几何级数的一般的 ODE 为：

\[\left(D\prod_{t=1}^n(\vartheta+b_t-1)\right)F=\left(\prod_{t=1}^n(\vartheta +a_t)\right)F \]

证明不会。但是超几何级数相当之强悍，它涵盖了超级多的常见的幂级数。比如几何级数 \(\frac{1}{1-x}\)​、指数函数 \(e^x\)、（几乎是）对数函数 \(x^{-1}\ln(1+x)\)​​ 等等。所以记一下也无妨。

Note.

多说几句。《具体数学》上说，超几何级数涵盖了所有“相邻两项的系数之比是 \(k\)​​ 的有理分式”的级数。简单来说，如果令 \(c_k\) 为：

\[\begin{aligned} c_k=\frac{\prod_{j=1}^ma_{j}^{\overline k}}{k!\prod_{t=1}^nb_{t}^{\overline k}} \end{aligned} \]

比一下就是：

\[\frac{c_{k+1}}{c_k}=\frac{k!}{(k+1)!}\frac{\prod_{j=1}^ma_j^{\overline{k+1}}}{\prod_{j=1}^ma_j^{\overline{k}}}\frac{\prod_{t=1}^nb_t^{\overline{k}}}{\prod_{t=1}^nb_t^{\overline{k+1}}}=\frac{\prod_{j=1}^m(a_j+k)}{(k+1)\prod_{t=1}^n(b_t+k)} \]

ODE 与函数运算

我们已经得到了五类基本函数的 ODE，之后只需要理解在若干函数运算下 ODE 的相应变化，我们（理论上）就可以给出任意初等函数的 ODE。

在此之前，我们先来简单看一条比较重要的引理：

Lemma.

形式幂级数 \(u(x)\) 微分有限当且仅当 \(\dim_{\mathbb K(x)}\operatorname{span}_{\mathbb K(x)}\{u,u',u'',\dots\}\) 有限。

从右到左是极其容易的。从左到右感性理解就好，如果要深入理解证明思路，看一看下面“加法”部分就明白了。

加法

Theorem.

如果 \(f,g\) 是微分有限的，则 \(f+g\)​ 也是微分有限的。

---

Proof.

假设 \(f,g\)​ 各自满足微分方程：

\[\sum_{k=0}^np_k(x)f^{(k)}=0\\ \sum_{k=0}^mq_k(x)g^{(k)}=0 \]

则我们只需要构建一个关于 \((f+g)\)​ 的微分方程即可完成证明。

方程本身给出了关系：

\[f^{(n)}=-p_n(x)^{-1}\sum_{k=0}^{n-1}p_k(x)f^{(k)} \]

两边同时对 \(x\)​ 求导并用上式替换 \(f^{(n)}\)​，我们即可得到 \(f^{(n+1)}\)​ 关于 \(f,f',f'',\dots,f^{(n-1)}\)​ 的线性表示。

重复这个过程，我们可以得到 \(f^{(n+k)}\) 关于 \(f,f',f'',\dots,f^{(n-1)}\) 的线性表示，其中 \(k\in \mathbb N\)。

\(g\)​ 类似。于是任意的 \((f+g)^{(k)}\)​ 都可以被 \(f,f',f'',\dots,f^{(n-1)},g,g',g'',\dots,g^{(m-1)}\) 线性表示。

从而 \((f+g),(f+g)',(f+g)^{''},\dots,(f+g)^{(n+m)}\) 必然线性相关，所以存在一个阶数不超过 \(n+m\) 的 ODE。实现的时候只需要维护 \((f+g)^{(k)}\) 的线性基即可。

乘法

Theorem.

如果 \(f,g\) 是微分有限的，则 \(fg\) 也是微分有限的。

---

Proof.

使用公式：

\[(fg)^{(n)}=\sum_{k=0}^{n}\binom{n}{k}f^{(k)}g^{(n-k)} \]

所以任意的 \((fg)^{(n)}\)​​ 均可被 \(f^{(k)}g^{(j)}\) 线性表示，其中 \(0\le k<n,0\le j<m\)​。故必然存在一个阶数不超过 \(nm\) 的 \(fg\) 的 ODE。

复合

Definition.

一个形式幂级数 \(A(x)\)​​​ 是“代数形式幂级数”，当且仅当 \(A(x)\)​​ 可以作为某个系数为 \(\mathbb K(x)\)​​​​ 的多项式方程的根。其中 \(\mathbb K\) 为某个域，\(\mathbb K(x)\) 表示系数在 \(\mathbb K\) 内的有理分式域。

我们将要指出：

Theorem.

如果​ \(f\)​ 微分有限，而 \(g\)​​ 是一个代数形式幂级数，则 \(f\circ g\) 也微分有限。

---

Explanation.

考虑 \(g\)​​ 是有理分式的情况。反复使用链式法则我们可以得到：

\[(f\circ g)^{(k)}=\sum_{j=0}^{k}c_{k,j}(x)\left(f^{(j)}\circ g\right) \]

\(c_j(x)\)​​​ 是关于 \(g\)​​​ 及其多阶微分的线性组合，但是 \(g\)​ 是有理分式。

而 \(f^{(j)}\circ g\) 的 ODE 可以直接由 \(f\)​ 的 ODE 生成：

\[\sum_{k=0}^np_k(g)(f^{(k)}\circ g)=0 \]

因此 \(f^{(k)}\circ g\) 构成一组基，其中 \(0\le k<n\)​。进而必然存在一个阶数不超过 \(n\) 的 ODE。

---

假设 \(g\) 满足代数方程 \(P(x,g)=0\)​，其中 \(P(x,g)=\sum_{k=0}^{m-1}b_k{(x)}g^k\)。

仍然试图得到使用链式法则，但是这个时候我们就必须直面 \(g\) 的导数了。

先考虑 \(g'\)，我们来通过 \(P(x,g)\) 生成关系。在方程两侧同时对于 \(x\) 求导：

\[\frac{\text d}{\text dx}P(x,g)=\frac{\partial}{\partial x}P(x,g)+\frac{\partial}{\partial g}P(x,g)\cdot g'=0 \]

所以有 \(g'=-\frac{\frac{\partial}{\partial x}P(x,g)}{\frac{\partial}{\partial g}P(x,g)}\)。上下都是关于 \(g\) 的多项式，我们尝试将 \(g'\) 也写成关于 \(g\) 的多项式。

我们知道 \(P(x,g)=0\)，因此对于 \(P(x,g)\) 做截取。也即，我们只需要求：

\[\frac{\partial}{\partial x}P(x,g)+\frac{\partial}{\partial g}P(x,g)u(x,g)\equiv 0\pmod{P(x,g)} \]

“可以证明这个总有解，为了防止多解，我们再假设 \(P(x,g)\)​ 关于 \(g\)​ 在 \(\mathbb K(x)\)​ 上不可约”。解出来 \(u(x,g)\)​​​ 即可。据说这个过程需要多项式欧几里德，不过考虑到平时大家都用手算，无所谓了。

另一方面，我们回头看一下：

\[\sum_{k=0}^np_k(g)(f^{(k)}\circ g)=0 \]

怎么算 \(a_k(g)\)​​？“可以使用多项式欧几里德”。我觉得大家手算的时候还是直接复合好了。

如果 \(a_k(g)\)​ 的分子分母都比较短，我们把所有 \(a_k\) 通分，之后去掉分母即可将有理分式变成有限多项式。

最后补充说明一下（不是泼冷水），下面很多 ODE 都可以通过反复求导建立，所以不要忘了这个基本策略！

整式递推

整式递推是一类并不特殊的递推数列：

Definition.

对于无穷数列 \(\{a_k\}_{k\ge 0}\)​​​ 和有限非空多项式列 \(\{r_0,r_1,r_2,\dots,r_{m}\}\)​​​，如果对于任意的 \(n\ge m\)​，都有：

\[\sum_{k=0}^ma_{n-k}r_{k}(n)=0 \]

则我们称 \(r\)​ 为 \(a\)​ 的整式递推式。我们称存在整式递推式的无限数列为整式递推数列。

对于有穷数列 \(\{a_k\}_{k=0}^{p}\)，我们可以类似地定义。此处略去不谈。

整式递推涵盖的范围很广。比如，如果所有的 \(r\) 均取到常数，则这个递推式退化成了线性递推。

下面这个定理很好地指出了整式递推和数列的 GF \(A(x)\) 的 ODE 的关系（之一）：

Theorem.

一个数列 \(\{a_k\}_{k\ge 0}\) 是整式递推数列，当且仅当数列的 GF \(A(x)\) 是 D-finite 的。

显然，我们只需要给出一个构造性证明即可在证明定理的同时，给出从“整式递推”到“ODE”的转化途径。

从 ODE 到整式递推是极其容易的——我们只需要对于 ODE 取它的 \([x^n]\) 即可建立一个方程，怎么怎么移位一下就可以得到递推式了。

从整式递推到 ODE 反而不太容易。考虑到这一部分内容在 OI 中的应用基本都是前者，所以我并不准备展开讲。

例题

短分式幂

给定 \(F(x)=\frac{P(x)}{Q(x)}\)，满足 \(P(x),Q(x)\) 均为多项式且 \(\max\{\deg P,\deg Q\}=m\)。现在请求出 \(\ln F(x),F(x)^\alpha,\exp F(x)\) 的前 \(n\) 项。你可以认为 \(m\ll n\)，所以这里要求你给出一个 \(O(nm)\) 的算法。运算在某个域上进行。

---

经典问题了。

首先解决 \(\ln F(x)\)。显然有 \(\ln F(x)=\ln P(x)-\ln Q(x)\)，所以我们这里仅考察 \(U(x)=\ln P(x)\) 的情况。

记得上面我们提到的 \(\ln x\) 的微分方程吗？这里我们还暂时用不到它（笑。不过我们可以以此获得启发，在等式两侧同时微分：

\[P(x)U'(x)=P'(x) \]

这样我们实际上就得到了一个 \(U(x)\) 的微分方程。作为例子，我们再详细一点——两侧同时取 \([x^n]\)：

\[\sum_{k=0}^np_k(n+1-k)u_{n+1-k}=(n+1)p_{n+1} \]

整理一下即可得到：

\[u_{n+1}=\frac{1}{n+1}\left((n+1)p_{n+1}-\sum_{k=1}^nku_kp_{n+1-k}\right) \]

这相当于是建立了 \(\{u\}\) 的一个整式递推数列。由于 \(P(x)\) 比较短，因此 \(k\) 只需要枚举 \(O(m)\) 项。递推的复杂度就是 \(O(nm)\) 的。

其它两类照猫画虎。对于 \(V(x)=F(x)^\alpha\)，两侧同时微分：

\[\begin{aligned} F(x)V'(x)&=\alpha F'(x)V(x)\\ \alpha (P'(x)Q(x)-Q'(x)P(x))V(x)&=P(x)Q(x)V'(x) \end{aligned} \]

这里 \(\max\{\deg P'(x)Q(x),\deg Q'(x)P(x), \deg P(x)Q(x)\}\le 2m\)。因此我们可以先 \(O(m^2)\) 暴力卷积，之后再做递推。

对于 \(W(x)=\exp F(x)\)，两侧同时微分：

\[\begin{aligned} W'(x)&=F'(x)W(x)\\ (P'(x)Q(x)-Q'(x)P(x))W(x)&=Q(x)^2W'(x) \end{aligned} \]

类似处理即可。

---

我们在上面伏了一笔。上面的微分方程都是由普通的微分建立的，如何用上之前 ODE 与运算的结论呢？

举个例子，设 \(g(x)=x^\alpha\)，那么 \(V=g\circ F\)。我们知道 \(g(x)\) 满足 \(\alpha g-xg'=0\)，所以有 \(g'=\frac{\alpha g}{x}\)。

怎么变成 \(V\) 的微分方程？根据结论，\(V\) 及其任意阶微分必然都可以由 \(V=g\circ F\) 自身来线性表示。所以我们单独考察 \(V'\) 即可生成 ODE。

链式法则告诉我们 \(V'=(g\circ F)'=F'\cdot (g'\circ F)=F'\cdot \frac{\alpha V}{F}\)，因此有 \(F(x)V'(x)=\alpha F'(x)V(x)\)。这和我们直接微分的结果（当然）是一致的。

「UOJ514」通用测评号

开门见山，我们需要求：

\[n(n-1)\sum_{k\ge 0}k![x^k]\left(\frac{x^{b-1}n^{-b}}{(b-1)!}\left(e^{xn^{-1}}-A(x)\right)\left(e^{xn^{-1}}-B(x)\right)^{n-2}\right) \]

其中 \(A(x)=\sum_{k=0}^{a-1}\frac{x^kn^{-k}}{k!},B(x)=\sum_{k=0}^{b-1}\frac{x^kn^{-k}}{k!}\)。

数据满足 \(1\le n\le 250,1\le b<a\le 250\)，你只需要给出一个 \(O(n^3)\) 的算法（认为 \(n,a,b\) 同阶）。

---

首先发现把式子内层不完全展开之后可以得到 \(x^{\alpha}e^{\beta xn^{-1}}\) 的线性组合。而 \(\sum_{k\ge 0}k![x^k]x^{\alpha}e^{\beta xn^{-1}}\) 是容易计算的。因此，我们的目标就是算出来每一个单项之前的系数。

使用二项式定理展开 \((e^{xn^{-1}}-B(x))^{n-2}\)，我们相当于只需要快速算出 \(U_k(x)=A(x)B(x)^k,V_k(x)=B(x)^k\) 的所有系数。考虑到它们的项数之和就是 \(O(n^3)\) 级别的，我们能想到的一种比较好的方法，就是递推系数。

以较复杂的 \(U(x)\) 为例，我们在两侧同时微分：

\[U_k'(x)=A'(x)B(x)^k+kA(x)B'(x)B(x)^{k-1} \]

注意到 \(A'(x),B'(x)\) 都是对于 \(e^{xn^{-1}}\) 进行截断之后的结果，这直接导致它们的 ODE 极其简单：

\[\begin{aligned} A(x)-A'(x)=\frac{x^{a-1}n^{-a+1}}{(a-1)!}\\ B(x)-B'(x)=\frac{x^{b-1}n^{-b+1}}{(b-1)!} \end{aligned} \]

之后我们会详细讨论一下截断对于 ODE 的影响。

下面我们就记 \(\iota_A(x)=A(x)-A'(x),\iota_B(x)=B(x)-B'(x)\)。所以，我们可以将 ODE 改写成：

\[\begin{aligned} U_k'(x)&=(A(x)-\iota_A(x))B(x)^k+kA(x)(B(x)-\iota_B(x))B(x)^{k-1}\\ &=(k+1)U_k(x)-\iota_A(x)V_k(x)-k\iota_B(x)U_{k-1}(x) \end{aligned} \]

此外，我们显然还有边界值 \(U_0=A,V_0=B\)。因此，我们可以按照 \(k\) 从小到大递推，复杂度即为 \(O(n^3)\)。

「CTS2019」珍珠

有 \(n\) 个在范围 \([1,D]\) 内的整数均匀随机变量。

求至少能选出 \(m\) 个瓶子，使得存在一种方案，选择一些变量，并把选出来的每一个变量放到一个瓶子中，满足每个瓶子都恰好装两个值相同的变量的概率。

请输出概率乘上 \(D^n\) 后对 \(998244353\) 取模的值。

对于 \(100\%\) 的数据，满足 \(1\le D\le 10^5,0\le m\le 10^9,1\le n\le 10^9\)。

---

注意到，我们实际上只需要讨论一下有多少种值出现了奇数次即可。假如有 \(k\) 种值出现了奇数次，则我们要求 \(\lfloor\frac{n-k}{2}\rfloor\ge m\) 和 \(k\equiv n\pmod 2\)。

当我们写出 EGF 的时候，后面的同余限制实际上不用管，因为不满足的项根本不会贡献。假设第一条给出限制为 \(k\le K\)，答案为：

\[n!2^{-D}[x^n]\sum_{k=0}^K\binom{D}{k}(e^x-e^{-x})^k(e^x+e^{-x})^{D-k} \]

注意到和式内的其实是 \(e^x\) 的生成函数。因此我们可以换元得到 \(F(z)\)：

\[F(z)=z^{-D}\sum_{k=0}^K\binom{D}{k}(z^2-1)^k(z^2+1)^{D-k} \]

现在我们就可以直接计算 \(F(z)\) 了，据说可以 \(O(D\log^2 D)\) 算。如果将和式展开我们还可以得到 \(O(D\log D)\) 的做法。

知道了 \(F(z)=\sum_j f_jz^j\) 之后，我们相当于要计算 \(n!2^{-D}\sum_j f_j[x^n]e^{jx}\)，这个就可以直接计算了。

---

但是这样还不够好。我们将 \(F(z)\) 处理成：

\[F(z)=z^{-D}(z^2+1)^D\sum_{k=0}^K\binom{D}{k}\left(\frac{z^2-1}{z^2+1}\right)^k \]

再次换元。设 \(H(z)=\frac{z^2-1}{z^2+1},G(t)=\sum_{k=0}^K\binom{D}{k}t^k\)，现在就有 \(F(z)=z^{-D}(z^2+1)^DG(H(z))\)。

压力来到了 \((z^2+1)^DG(H(z))\) 身上，我们尝试给出这样一个 GF 的 ODE，这样就可以递推了。

先从最简单的 \(G(t)\) 开始。我们发现有 \(G(t)=(1+t)^D\bmod t^{K+1}\)。这能带来什么信息呢？

没有截断的 ODE 我们是知道的。那么，引入截断的 ODE 必然可以表示成：

\[DG(t)-(1+t)G'(t)=\mathscr D(t) \]

这个 \(D(t)\) 的性状取决于 LHS 的抵消情况。考察 LHS 的某一项系数：如果它的系数全由 \(G\) 中（在微分前的）原先指数 \(\le K\) 的项提供，那么这一项系数不会受到截断的影响；如果它的系数全由 \(G\) 中（在微分前的）原先指数 \(>K\) 的项提供，那么这一项系数必然全是 \(0\)。所以 \(D(t)\) 包含的就是 LHS 中系数同时由 \(\le K\) 和 \(>K\) 的项提供的那些项。

依据这样的分析，我们可以发现 \(\mathscr D(t)=(K+1)\binom{D}{K+1}t^K\)，因为 \(1\times G'(t)\) 恰好在 \(t^K\) 的位置缺了一项。因此：

\[DG(t)-(1+t)G'(t)=(K+1)\binom{D}{K+1}t^K \]

现在我们可以快速地递推 \(G\) 的系数了。

下一步，考虑 \(G(H(z))\) 的 ODE。从这里开始，我们设 \(P(z)=G(H(z))\)。

然而我不想动脑筋了，直接使用 ODE 机械化构造的方法。

首先，给出 \(G'(t)=\frac{DG(t)-\mathscr D(t)}{1+t}\)。根据结论，我们考察 \((G\circ H)'\) 即可得到一个 ODE。有 \((G\circ H)'=H'\cdot (G'\circ H)\)，所以有：

\[\begin{aligned} P'(z)&=H'(z)\cdot \frac{DP(z)-\mathscr D(H(z))}{1+H(z)}\\ (1+H(z))(H'(z))^{-1}P'(z)&=DP(z)-(K+1)\binom{D}{K+1}H(z)^K \end{aligned} \]

这个其实也可以线性推系数。

手算得到 \(H'(z)=\frac{4z}{(z^2+1)^2},1+H(z)=\frac{2z^2}{z^2+1}\)，二者恰好抵消得到 \(\frac{1}{2}z(z+1)\)，因此两个系数都很短。

\(H(z)^K\) 的系数虽然不方便直接算出来，但是我们可以使用短分式幂的方法和 \(P(z)\) 同步递推。

最后，我们只需要给出 \((z^2+1)^DP(z)\) 的 ODE 就好啦！仍然使用结论，\((z^2+1)^D\) 和 \(P(z)\) 各自都有一个一阶 ODE，因此表示 \(F(z)\) 及其任意次微分的结果也只需要用 \((z^2+1)^D\cdot P(z)\) 一项。因此直接考察 \(F'(z)\) 肯定可以得到一个 ODE：

\[F'(z)=D\frac{(z^2+1)^D}{z^2+1}P(z)+(z^2+1)^D\cdot \frac{2}{z(z+1)}(DP(z)-\mathscr D(H(z))) \]

整理一下即可得到最终的 ODE：

\[z(z+1)(z^2+1)F'(z)=(Dz(z+1)+2D(z^2+1))F(z)-2(K+1)\binom{D}{K+1}(z^2-1)^K(z^2+1)^{D+1-K} \]

递推系数都比较短。后面的 \(R(z)=(z^2-1)^K(z^2+1)^{D+1-K}\) 可以建立另外的微分方程：

\[(z^2-1)(z^2+1)R'(z)=(K(z^2+1)+(D+1-K)(z^2-1))R(z) \]

同步递推即可。这样我们就可以正儿八经地做到 \(O(D)\) 了。

初探载谭

先让我们从一道例题开始。

「CF932E」Binomial Sum

给定 \(n,m\)，求出：

\[\sum_{k=0}^n\binom{n}{k}k^m \bmod (10^9+7) \]

请给出一个 \(O(m+\log n)\) 的算法。

---

我们可以敏锐地觉察到，这个问题要求的其实就是：

\[m^n \]

这样通过基础多项式运算我们可以得到一个 \(O(m\log m)\) 的算法，已经相当不错了。

我们依然可以注意到，这个 GF 是一个 \(F(G(x))\) 的复合形式，其中 \(F(z)=(1+z)^n,G(x)=e^x\)。

回忆一下，复合有两种良好定义：

对于 \(F(G(x))\)，当 \(F\) 有限或者 \([x^0]G(x)=0\) 的时候，该复合是良定义的。

这里虽然 \((1+x)^n\) 已经满足了“\(F\) 有限”的条件，但是它太长长长长长了，不方便我们处理。我们尝试将内层函数变成 \([x^0]G(x)=0\) 的形式。

这很简单，由于 \(G(0)=1\)，我们只需要在 \(1\) 处展开 \(F(z)\) 即可：

\[F(z)=\sum_{k\ge 0}f_z(z-1)^k \]

关键的一步来了：由于 \([x^0](e^x-1)=0\)，所以我们的求和上界只需要保留到 \(m\)。也就是说，当我们求 \([x^m]F(G(x))\) 的时候，我们用 \(\mathscr F(z)=F(z)\bmod (z-1)^{m+1}\) 不会导致结果出错，还可以降低次数！

现在，我们可以先算出来 \(F(z+1)\bmod z^{m+1}\)，然后做一个多项式平移得到 \(\mathscr F(z)\)。这个算法仍然是 \(O(m\log m)\) 的。

但是 \(F(z),\mathscr F(z)\) 的形式都很简洁，我们尝试给出一个 ODE。先给出 \(\mathscr F(z+1)\) 的 ODE——我们之前已经研究过的那一个：

\[n\mathscr F(z+1)-(2+z)\mathscr F'(z+1)=2^{n-m}(m+1)\binom{n}{m+1}z^m \]

现在再平移就 OK 辣！

\[n\mathscr F(z)-(1+z)\mathscr F'(z)=2^{n-m}(m+1)\binom{n}{m+1}(z-1)^m \]

我们还需要解决 \([z^0]\mathscr F(z)\) 的边界情况。实际上只需要带进去算一算：

\[\mathscr F(0)=\sum_{k=0}^K\binom{n}{k}(-1)^k2^{n-k} \]

就好了。之后系数是可以 \(O(m)\) 递推的。而答案就是 \(m![x^m]\mathscr F(e^x)\)，线性筛整数幂就好了。

EI 的载谭

我们把上面的方法再扩展一下。

先考虑朴素的情况。假如我们有两个 GF \(F(z),G(x)\)，满足 \(F(z)\) 微分有限，它的 ODE 最短阶数为 \(n\)。我们现在就要求一个 \([x^{m-1}]F(G(x))\)。

首先，我们把复合构造成良复合形式。假设 \(c=[x^0]G(x)\)，我们可以在 \(c\) 处展开 \(F(z)=\sum_{k\ge 0}f_k(z-c)^k\)。这样的话，在复合 \(G(x)\) 的时候我们就可以做一个截断 \(\mathscr F(z+c)=F(z+c)\bmod z^m\)。显然，该截断满足 \([x^{m-1}]\mathscr F(G(x))=[x^{m-1}]F(G(x))\)。展开之后，我们有 \([x^{m-1}]=[x^{m-1}]\sum_{k=0}^{m-1}([z^k] \mathscr F(z))G(x)^k\)。问题变成了如何快速计算 \(G(x)^k\) 的单项系数和 \(\mathscr F(z)\) 的所有系数。

接下来，我们可以用多种方式给出 \(\mathscr F\) 的系数。首先，如果 \(F\) 自己有好算的封闭形式，那么我们可以先使用多项式工业算出来 \(F(z+c)\bmod z^{m}\)，然后做一个平移，这是朴素的用不着 ODE 的方法。

其次，我们用上之前对于截断的讨论。如果原先的 \(F\) 的 ODE 为：

\[\sum_{k=0}^np_k(z)F^{(k)}(z)=0 \]

我们可以通过某种方法找到一个修正因子 \(\mathscr D(z)\)：

\[\sum_{k=0}^np_k(z+c)\mathscr F^{(k)}(z+c)=\mathscr D(z) \]

再平移回来得到：

\[\sum_{k=0}^np_k(z)\mathscr F^{(k)}(z)=\mathscr D(z-c) \]

现在就可以递推 \(\mathscr F\) 的系数了。

快速计算 \(G(x)^k\) 的单项系数的话......如果 \(G(x)\) 本身性质比较好，比如 \(G(x)=e^x\) 或者 \(\frac{1}{1-x}\) 之类的，我们可以很快地算出来单项系数。退一步，如果存在复合逆的话，我们也可以使用拉格朗日反演算 \([x^{m-1}]\frac{1}{1-yG(x)}\bmod y^m\)。其它情况的方法暂时不明。

实际上，我们可以将它推广到 \(G(x)^k\) 的线性组合，这也是 EI 在原博客中给出的形式：

如果对于任意的 \(0\le k<m\) 的整数 \(k\)，下面的：

\[\sum_{j=0}^ma_j[x^j]G(x)^k \]

都已知，那么我们可以通过上述方法快速地算出：

\[\sum_{j=0}^ma_j[x^j]F(G(x)) \]

从我们的结果到这一步是容易的，只需要交换求和号即可快速计算。

---

所以这个做法有一定的局限性：

1. 由于我还不会什么 \(\mathscr D\) 的机械化算法，所以一旦 \(n\) 很大或者 \(p_k(z)\) 很长，这个方法就会变得尤其繁琐。

2. 这个方法只能快速算一项结果。当需要一次求多项结果的时候，它会显得有点鸡肋，甚至可能直接失效......