一种线性求和黑魔法

Credit: Entropy Increaser.

接下来我们介绍一种由 Entropy Increaser (Baitian Li) 等人发明的一种针对多项式和形式幂级数的线性求和算法。这个算法推导简单，并且它统一了很多关于多项式的求和的问题中 ad hoc 的推导，是一个值得一学的新技术。

算法. 设 \(F(x)\) 是一个微分有限 (D-Finite) 的函数。令 \(G(x)\) 为一个生成函数，\((a_i)\) 为一个未知数列。设对每个 \(0\le k\le n\) 给定了

\[s_k = \sum_{j=0}^n a_j[x^j]G(x)^k, \]

那么存在一个线性(\(O(n)\))时间的算法计算

\[\sum_{j=0}^{n} a_j [x^j]F(G(x)). \]

注. 所谓微分有限，是指 \(F\) 满足一个微分方程：

\[\sum_{i=0}^{t} p_i(z) F^{(i)}(z) = 0, \]

其中 \(p_i(z)\) 为一系列多项式，\(t\) 是一个有限的常数。

在很多的应用中，我们会取 \(G(x) = \exp(x)\)，此时 \([x^j]\exp(kx) =\frac{k^j}{j!}\). 在此情形下，上述算法自然对应了一类多项式求和问题。

推论. 假设 \(F(z) = \sum_{i\ge 0} f_i z^i\) 微分有限。令 \(P(x)=\sum_{i=0}^D c_i\cdot x^i\) 为一个 \(D\) 次多项式。给定 \(P(0),P(1),\dots,P(D)\)，存在一个算法在 \(O(D)\) 时间内计算

\[\sum_{j=0}^D c_j \left( \sum_{i=0}^{\infty} f_i \cdot i^j\right) =\sum_{i=0}^{\infty} f_i\sum_{j=0}^D c_j \cdot i^j = \sum_{i=0}^{\infty} f_i \cdot P(i). \]

(如果觉得这个推论不好理解，请参考下面 [例1. 插值] 中的内容。)

例子

在叙述算法之前，我们先代入几个例子看看这个算法的通用性。

例1. 插值

问题. 假设 \(P(x)\) 是一个度数为 \(D\) 的多项式，已经给出 \(P(0),P(1),\dots, P(D)\)，我们想在 \(O(D)\) 的时间内计算出 \(P(n)\)。

设 $P(x)=\sum_{j=0}^D \frac{a_j}{j!} x^j $ 。我们取 \(G(x)\) 为 \(\exp(x)\)，那么对于任意 \(k\in [0,D]\)，我们有

\[\sum_{j=0}^D a_j[x^j]G(x)^k = \sum_{j=0}^{D} a_j [x^j]\exp(kx) = \sum_{j=0}^D \frac{a_j}{j!} k^j = P(k). \]

因此套用算法所需要的输入已经齐全。我们构造 \(F(G)=G^n\)​，\(F\) 是一个微分有限的函数，它满足微分方程 \(z\cdot F'(z) - n\cdot F(z) = 0\). 最后注意到

\[\sum_{j=0}^D a_j[x^j] F(G(x)) = \sum_{j=0}^D a_j[x^j]\exp(nx) = \sum_{j=0}^D \frac{a^j}{j!} n^j = P(n), \]

这正是我们想要的答案。

例2. （清华集训）如何优雅地求和

问题. 给定一个 \(D\) 次多项式 \(f\) 在 \(f(0),\dots, f(D)\) 处的取值，再输入 \(n\le 10^9\) 和 \(\alpha\in (0, 1)\)，计算

\[Q(f,n,\alpha) = \sum_{k=0}^n f(k) \binom{n}{k} \alpha^k(1-\alpha)^{n-k}. \]

我们依然令 \(f\) 的系数为 \((a_i/i!)\)，取 \(G(x) = \exp(x)\)，那么对任意 \(k\in [0,D]\)，

\[\sum_{j=0}^D a_j[x^j] G(x)^k = \sum_{j=0}^D \frac{a_j}{j!} k^j = f(k). \]

我们构造 \(F\) 为 \(F(z) = \sum_{k=0}^n \binom{n}{k} (\alpha z)^k (1-\alpha)^{n-k} = (\alpha z + 1-\alpha)^n\)，那么

\[\begin{aligned} \sum_{j=0}^D a_j[x^j] F(G(x)) &= \sum_{j=0}^D a_j[x^j] \left( \sum_{k=0}^n (\alpha G(x))^k (1-\alpha)^{n-k}\binom{n}{k} \right) \\ &= \sum_{j=0}^D \frac{a_j}{j!} \sum_{k=0}^n \alpha^k k^j(1-\alpha)^{n-k}\binom{n}{k} \\ &= \sum_{k=0}^n\binom{n}{k} \alpha^k (1-\alpha)^{n-k}\sum_{j=0}^D \frac{a_j}{j!}\cdot k^j \\ &= \sum_{k=0}^n\binom{n}{k}\alpha^k(1-\alpha)^{n-k} f(k) \\ &= Q(f,n,\alpha). \end{aligned} \]

不难验证 \(F\) 是一个微分有限的幂级数，它满足微分方程 \((\alpha z + 1 - \alpha) F' - n F = 0\)。因此，套用上述算法可以在 \(O(D)\) 的时间内解决此问题，比传统解法的 \(O(D\log D)\) 更加优越。

之后我们会介绍更多的例子，下面我们来叙述算法。

算法

首先，我们把所求答案记为 \(ANS\)，也就是

\[ANS := \sum_{j=0}^{n} a_j [x^j]F(G(x)). \]

我们首先来研究 \(F(G(x))\)。设 \(F(z) =\sum_{i\ge 0}f_i z^i\)，那么我们可以自然地把 \(F(G(x))\) 写为

\[F(G(x)) = \sum_{i = 0}^{\infty} f_i \cdot G(x)^i. \]

这样的展开并不利于我们做计算，因为对于任意大的 \(i\)，都可能有 \([x^j]G(x)^i \ne 0\)。比较好的做法是，设 \(G(0) = c\)，我们考虑把 \(F\) 写为 \(F(z) = \sum_{i\ge 0} \tilde{f}_i \cdot (z-c)^i\)（也就是说，在常数 \(c\) 处展开 \(F\)），那么我们就得到了

\[F(G(x)) = \sum_{i=0}^{\infty} \tilde{f}_i \cdot (G(x)-c)^i. \]

这里我们注意到，当 \(i > j\) 时，一定有 \(x^{j+1}|(G(x)-c)^i\)，从而 \([x^j]G(x)^i = 0\)。由于我们只关心 \([x^{\le n}]F(G(x))\)，对于上式，我们只需要展开到第 \(n\) 层即可。于是，我们考虑定义 \(\hat{F}\) 为 \(F\) 的一个截断：

\[\hat{F}(z) = \sum_{i=0}^n \tilde{f}_i \cdot (z-c)^i = F(z) \bmod (z-c)^{n+1}. \]

那么我们有

\[\begin{aligned} ANS &= \left(\sum_{j=0}^n a_j[x^j] \right) F(G(x)) \\ &= \left(\sum_{j=0}^n a_j[x^j] \right) \left( \sum_{i=0}^\infty \tilde{f}_i \cdot (G(x) - c)^i \right) \\ &= \left(\sum_{j=0}^n a_j[x^j] \right) \left( \sum_{i=0}^{\color{red}{n}} \tilde{f}_i \cdot (G(x) - c)^i \right) \\ &= \left( \sum_{j=0}^n a_j[x^j] \right) \hat{F}(G(x)). \end{aligned} \]

现在注意到 \(\hat{F}(z)\) 是一个度数不超过 \(n\) 的多项式，我们知道它在 \(c\) 处的展开式为 \(\hat{F}(z) = \sum_{i=0}^n \tilde{f}_i\cdot (z-c)^i\)。假设它在 \(0\) 处的展开式为

\[\hat{F}(z) = \sum_{i=0}^n \hat{f}_i z^i. \]

假设我们能算出 \((\hat{f}_i)_{i=0}^n\)，我们可以依以下方法在 \(O(n)\) 时间内计算出答案：

\[\begin{aligned} ANS &= \left( \sum_{j=0}^n a_j[x^j]\right) \left( \sum_{i=0}^n \hat{f}_i \cdot G(x)^i \right) \\ &= \sum_{i=0}^n \hat{f}_i \left( \sum_{j=0}^n a_j [x^j] G(x)^i \right) \\ &= \sum_{i=0}^n \hat{f}_i \cdot s_i. \end{aligned} \]

于是，我们最后需要解决的问题是：如何在 \(O(n)\) 的时间内计算出 \((\hat{f}_i)_{i=0}^n\)？

我们需要利用微分方程。设 \(F(z)\) 满足的 \(t\) 阶微分方程为

\[\sum_{i=0}^t P_i(z)\cdot F^{(i)}(z) = 0. \]

注. 在大多数应用中，我们有 \(t=1\)，并且多项式 \(P_i(z)\) 的度数不超过 \(1\)。简单起见，在思考之后的推导时不妨先带入 \(t=1\)。

由于平移不变性，我们也有

\[\sum_{j=0}^t P_j(z + c) \cdot F^{(j)}(z+c) = 0. \]

注意到 \(F(z+c) = \sum_{i\ge 0} f_i\cdot (z+c)^i = \sum_{i\ge 0} \tilde{f}_i\cdot z^i\)，而上式给出了一个关于 \(\tilde{f}_i\) 的 \(O(t)\) 阶递推关系，不妨设为 \(\tilde{f}_i = \sum_{j=1}^{O(t)} \tilde{f}_{i-j} \cdot c_j(i)\).

回忆 \(\hat{F}(z+c)=\sum_{i=0}^{\color{red} n} \tilde{f}_i \cdot z^i\) 是 \(F(z+c)\) 的一个截断，对于 \(\hat{F}(z+c)\) 的系数 \((\tilde{f}_0,\tilde{f}_1,\dots, \tilde{f}_{n}, 0, 0,\dots)\) 来说，递推关系 \(\tilde{f}_i =\sum_{j=1}^{O(t)} \tilde{f}_{i-j}\cdot c_j(i)\) 在 \(i\le n\) 的时候依然成立，在 \(i-n\gg t\) 的时候由于递推的齐次性也成立，因此我们只需要在 \(|i-n|\lesssim t\) 的位置加一个扰动项即可。形式化地说，考虑如下幂级数

\[D(z):= \sum_{j=0}^t P_j(z+c)\cdot \left( F^{(j)}(z+c) -\hat{F}^{(j)}(z+c) \right) = \sum_{j=0}^t P_j(z+c) \left(\sum_{i=n+1}^\infty \tilde{f}_i z^i\right)^{(j)}. \]

当 \(i- n \gg t\) 时，\([z^i]D(z) = 0\)，当 \(n-i \gg t\) 时，同样也有 \([z^i]D(z) = 0\)。因此 \(D\) 是一个只包含 \(O(t)\) 个非零项的稀疏多项式。于是我们有

\[\sum_{j=0}^t P_j(z+c)\hat{F}^{(j)}(z+c) = -D(z)， \]

将其平移回去，得到

\[\sum_{j=0}^t P_j(z) \hat{F}^{(j)}(z)= -D(z-c). \]

这个式子蕴含了一个关于 \(\hat{f}_i\) 的递推关系，依此式进行递推即可。

总结

1. 要计算 \(ANS := \sum_{j=0}^{n} a_j [x^j]F(G(x))\)，如果依定义展开 \(F(G(x))=\sum_{i\ge 0} f_i \cdot G(x)^i\)，那么任意大的 \(G(x)^i\) 项对答案都会有贡献，无法快速计算。因此我们在 \(G(0)=c\) 处展开得到 \(\sum_{j=0}^n a_j[x^j] \sum_{i=0}^{\infty} \tilde{f}_i (G(x)-c)^i\).
2. 此时我们可以做截断 \(\sum_{j=0}^n a_j[x^j] \sum_{i=0}^{\color{red}n} \tilde{f}_i (G(x)-c)^i\)，令 \(\hat{F}=\sum_{i=0}^n \hat{f}_i z^i\) 为截断后的多项式。我们发现 \(ANS = \sum_{i=0}^n \hat{f}_i s_i\)，于是只需要计算 \((\hat{f}_i)_{i=0}^n\)。
3. 最后我们计算 \(\hat{f}_i\) . 考虑它是怎样得来的：将 \(F(z)\) 平移到 \(F(z+c)\)，截断取前 \(n\) 项得到 \(\hat{F}(z+c)\)，再平移得到 \(\hat{F}(z)\). 由 \(F\) 的微分方程得到 \(F(z+c)\) 的系数 \((\tilde{f}_i)\) 满足一个齐次递推，这个齐次递推蕴含了关于 \((\tilde{f}_0,\dots, \tilde{f}_n, 0, \dots)\) 的一个非齐次递推，把这个非齐次递推平移得到关于 \(\hat{f}_i\) 的递推即可。

更多例子

TJOI 2016 求和. 给定 \(n\le 2\times 10^5\)，计算 \(f(n) = \sum_{i=0}^n \sum_{j=0}^i S(i,j)\cdot 2^i\cdot j!\)，其中 \(S(i,j)\) 表示第二类斯特林数。

考虑此式的一个组合意义：准备 \(i\le n\) 个元素，划分成带标号的 \(j\le i\) 组，并对每一组黑白染色的方案数。

令 \(G(x)=e^x\), \(G'(x)=G(x)-1\), 那么划分成有顺序的 \(j\) 组，并对每一组染色的指数生成函数为 \((2G'(x))^j\)，对 \(j\) 求和得到划分成若干组的生成函数为

\[\frac{1}{1-2G'} = \frac{1}{3-2G}. \]

于是定义 \(F(z)=\frac{1}{3-2z}\)。我们要对 \(i\in \{0,1,\dots, n\}\) 求和，这也就是

\[\sum_{j=0}^{n}j!\cdot [x^j] F(G(x)). \]

要套用上述算法，我们需要能计算

\[s_k = \sum_{j=0}^n j![x^j] \exp(kx) = \sum_{j=0}^n k^j, \]

这只是一个简单的等比数列求和。

国王奇遇记加强版之再加强版. 给定 \(n\le 10^9, q\le 10^6\)，计算 \(\sum_{i=0}^n q^i i^q\)。

我们使用推论。构造多项式 \(P(x) = x^q\)。构造辅助函数 \(F(z) = \sum_{i=0}^n q^i z^i= \frac{1-(zq)^{n+1}}{1-zq}\)，可以验证 \(F\) 满足某个 \(1\) 阶微分方程。我们要计算的是 \(\sum_{i=0}^{n} f_i P(i)\).

可以在 \(O(q)\) 的时间内计算出 \(P(0),\dots, P(q)\)，套用算法即可。

一些常见生成函数的微分方程

来源：杜老师的博客。

\[\begin{aligned} g_1(x) = \sum_{i\ge 0} \frac{x^i}{i!} &\Rightarrow g_1 = g_1'\\ g_2(x) = \sum_{i=0}^k \frac{x^i}{i!} &\Rightarrow g_2 = g_2' + \frac{x^i}{i!} \\ g_3(x) = \sum_{i\ge 0} x^i i! &\Rightarrow g_3 = g_3'x^2 + g_3x + 1 \\ g_4(x)= \sum_{i\ge 0} \frac{x^i}{i!(i+k)!} &\Rightarrow g_4 = g_4''x + (k+1) g_4' \\ g_5(x) = \sum_{n\ge 0} \frac{1}{(k-1)n+1}\binom{kn}{n} x^{(k-1)n + 1} &\Rightarrow g_5 = \frac{kxg_5'}{1 + (k-1)g_5'}\\ g_6(x) = (1+x)^a(1-x)^b &\Rightarrow g_6'=\frac{ag_6}{1+x} +\frac{bg_6}{1-x} \\ g_7(x) = \sum_{i=0}^k \binom{n}{i}a^ib^{n-i} &\Rightarrow ng_7(a'+b')-g_7(a+b)=n\binom{n-1}{k}(a'a^kb^{n-k}-b'a^{k+1} b^{n-k-1}) \\ g_8(x) = g^n &\Rightarrow ng_8 f'=fg_8' \end{aligned} \]

注意这里我们讨论的是形式导数和形式幂级数，一般来说它们不能以微积分意义上的级数理解（举例来说，我们注意到 \(g_3(x)\) 在 \(x\ne 0\) 时甚至不收敛）。

习题：构造 \(F(z) = \frac{1}{3-2z}\) 所满足的微分方程。

拓展阅读

1. 关于整数拆分和五边形数定理的一个博文
2. 快速求和：Entropy Increaser 的讲稿
3. 微分法：杜老师的讲稿
4. 计算斯特林数的一些快速算法 (Entropy Increaser)