组合数学学习笔记（三）：生成函数 1（OGF）

orange_new · 2024-11-20 14:35:49 · 算法·理论

\text{A generating function is a device somewhat similar to a bag.}

\text{Instead of carrying many little objects detachedly, which could be embarrassing,}

\text{we put them all in a bag.}

\text{And then, we have only one object to carry, the bag.}

\text{——George Polya}

其实，一个数列可以表示一个函数，特别是在离散意义下。但有些时候，用一个函数来计数要方便一些，这时我们就需要用到生成函数。生成函数又叫做形式幂级数。

定义一个无穷数列 h_0, h_1, h_2, \dots, h_n, \dots 的生成函数 H(x) 为 h_0 + h_1x + h_2x^2 + \dots + h_nx^n + \dots。

而一个有限的数列，可以通过在后面补 0 的方式转化为无穷数列，这样就都可以用生成函数了。

我觉得，生成函数就像一个压缩包，我们用离散微积分或者特殊的数将一个无穷多项式压缩。而我们之前提到过的泰勒展开，它可以把一个函数重新拆成一个无穷多项式，相当于解压操作，此时我们就具备了操作生成函数的能力，也就拥有了操作原数列的能力。

生成函数的计算

其实生成函数的某些计算和暴力多项式计算一模一样，只是有些时候，我们不想暴力用泰勒展开，因此我们需要了解生成函数的操作对原数列的影响。

加法

生成函数的加法和减法同普通多项式一样，也就是 a\displaystyle\sum_{i \geq 0}f_i x^i \pm b\sum_{i \geq 0}g_i x^i = \sum_{i \geq 0}(f_i \pm g_i) x^i，此时我们得到了数列 \langle af_n + bg_n \rangle 的生成函数。

平移

我们现在要将数列 \langle f_0, f_1, f_2, \dots \rangle 向右平移 m 位，变成 \langle 0, \dots, 0, f_0, f_1, f_2, \dots \rangle，此时原来第 n 项的系数就变成了第 n + m 项的系数，于是这个数列的生成函数就变成了 \displaystyle\sum_{i \geq 0} f_i x^{i + m} = x^m \sum_{i \geq 0} f_i x^i，于是想让一个生成函数向右平移 m 位，只用乘上一个 x^m。对应到递推式中，如果出现 f_i = g_{i - k} 的形式，那么对应到生成函数上，相当于将数列 g 右平移了 k 次，再加到了 f 上，因此 F(x) = x^k G(x)

向左平移也是一样的，这里不再赘述。

乘法

这里先来介绍一下卷积。严谨定义的卷积我没有看懂，只能说一下我感性理解的卷积，这是一种特殊的乘法，就是如果两个多项式函数两个下标满足通过某种运算后变成一个定值，那么我们把它们乘起来，比如之前见过的范德蒙德卷积 \displaystyle\sum_{i + j = n} \binom ni \binom nj（此时 i + j 是定值），以及狄利克雷卷积 \displaystyle\sum_{i \times j = n} f(i) g(j)（此时 i \times j 是定值）。

生成函数的乘法可以通过乘法分配律来推导，可以发现生成函数的乘法就是这两个数列的卷积 \displaystyle\sum_{i \geq 0} \left(\sum_{j = 0}^i a_j b_{i - j} \right)x^i。因此，也可以发现生成函数的单位元就是 \displaystyle\sum_{i \geq 0}[i = 0]x^i，也就是 1。卷积对于原数列的影响比较复杂，我们放在后面再说。

求逆

考虑生成函数求逆，也就是我们在已知生成函数 F(x) 的情况下要求出另一个生成函数 G(x) 满足 F(x)G(x) = 1，于是我们要求 f_0g_0 = 1, \displaystyle\sum_{j = 0}^i f_j g_{i - j}= 0(i > 0)，将 g_i 这一项单独提出来可以得到 g_i f_0 = -\displaystyle\sum_{j = 1}^i f_j g_{i - j}，最终可以化简成 g_i = -g_0 \displaystyle\sum_{j = 1}^i f_j g_{i - j}，这样就可以通过 O(n^2) 递推求出 g_n，如果套用牛顿迭代，则可以在 O(n \log_2 n) 的时间复杂度内求出 g_n。此时，我们也知道了生成函数 F(x) 存在逆的充要条件是 g_0 \neq 0。不过目前似乎没有用过求逆去操作数列。

除法

现在我们就知道如何对生成函数做除法了，考虑 \displaystyle\frac{F(x)}{G(x)}，也就是 F(x)G(x)^{-1}，假设 G(x)^{-1} = H(x)，那么 H(x) = \displaystyle\sum_{i \geq 0} \left(- h_0\sum_{j = 1}^i g_j h_{i - j} \right) x^i，于是答案 R(x) 就是 \displaystyle\sum_{i \geq 0} \left(r_i = \begin{cases} f_0g_0^{-1} & i = 0 \\ \left(f_i - \displaystyle\sum_{j = 1}^{i}g_j r_{i - j} \right) \times g_0^{-1} & i < d \\ -g_0^{-1} \displaystyle\sum_{j = 1}^{d}g_j r_{i - j} & i \geq d\end{cases} \right) x^i。不过目前似乎也没有用过除法去操作数列。

求导与积分

生成函数的求导与普通多项式一样，就是 F'(x) = \displaystyle\sum_{i \geq 1} i f_i x^{i - 1}，如果我们再将其向右平移一位，可以得到 x F'(x) = \displaystyle\sum_{i \geq 1} i f_i x^i，此时我们给原数列的每一项都乘以了 i，如果求导与平移合理，我们可以用任何 i 的多项式乘以 x^i 项。

积分和求导互为逆运算，因此 \displaystyle\int F(x) \mathrm dx = \sum_{i \geq 1} \frac 1i f_{i - 1} x^i，如果我们再向左平移一位，就可以得到 \displaystyle\int \frac{F(x)}{x} \mathrm dx = \sum_{i \geq 1} \frac 1i f_{i} x^i，那么求导、积分、平移结合起来，我们可以用任何 i 的分式乘以 x^i 项。

求前缀和

这是生成函数一个超级妙用，建议先看一下下面 \langle 1, 1, 1, \dots \rangle 数列的生成函数的封闭形式，我们用 \displaystyle\sum_{i \geq 0} x^i 乘以某个多项式函数 F(x)，可以得到 \displaystyle\frac{1}{1 - x} F(x) = \sum_{i \geq 0} \left(\sum_{j \geq 0} f_{i - j} \right) x^i = \sum_{i \geq 0} \left(\sum_{j \leq i} f_i \right) x^n，可以发现这对原数列求了一次前缀和。推广一下，当 F(x) 乘以 \displaystyle\frac{1}{(1 - x)^k}，就是对原数列做 k 次前缀和，此时数列中第 i 个数就变成了 \displaystyle\binom{i + k - 1}{k - 1} f_i。

我们再将其推广，还是建议去看一下下面 \langle 1, \underbrace{0, \dots, 0}_{p - 1 个 0}, 1, \underbrace{0, \dots, 0}_{p - 1 个 0}, 1, \dots\rangle 数列的生成函数的封闭形式，我们用 \displaystyle\sum_{i \geq 0} x^{ip} 乘以某个多项式函数 F(x)，可以得到 \displaystyle\frac{1}{1 - x^p} F(x) = \sum_{i \geq 0} \left(\sum_{j \geq 0} f_{i - jp} \right) x^i，它对于原数列的影响就是，每个位置往前递减，每减 p 次做一次前缀和。此时序列中第 i 个数就变成了 \displaystyle\sum_{j = 0}^{\lfloor \frac ij \rfloor} f_{i - jp}。此做法依然可以拓展，当 F(x) 乘以 \displaystyle\frac{1}{(1 - x^p)^k}，就是对原数列做了 k 次这样的前缀和。

普通数列的生成函数

这些数列一般只用掌握一些式子的变换或离散微积分就可以求出该数列的生成函数的封闭形式，比如：

我们发现这个式子不太好求和，考虑到 1 + x + x^2 + \dots = \displaystyle\frac{1}{1 - x}。

对两边求导，可以得到 \displaystyle\frac{d}{dx} (1 + x + x^2 + \dots) = \frac{d}{dx}(\frac{1}{1 - x})。

求出两边的式子，根据多项式求导公式（见多项式学习笔记（一）：多项式基础知识、快速傅里叶变换）f'_i = (i + 1)f_{i + 1}，可以得到 1 + 2x + 3x^2 + \dots = \displaystyle\frac{1}{(1 - x)^2}。

那么原数列的生成函数就是 \displaystyle\frac{1}{(1 - x)^2}。

此时我们发现，我们用一个较为简单的函数表示了一个复杂的无穷数列，但是，就第 2 个数列的生成函数而言，当我们取 x \geq 1 时，发现两边的式子根本不相等，只是通过一系列运算法则和公式推导出来的，这时，这个生成函数就被叫做形式幂级数，化简后的结果就叫做这个生成函数的封闭形式。

练习一

求 a_n = f(n) 的生成函数的封闭形式，其中 f(x) = \displaystyle\sum_{i = 0}^{k - 1} f_i x^i。

解：这就是一个典型的利用求导与积分算生成函数的例子。

原式的生成函数 = \displaystyle\sum_{i \geq 0} x^i \sum_{j = 0}^{k - 1} f_j i^j

根据斯特林数的展开式，可以得到：\displaystyle\sum_{i \geq 0} x^i \sum_{j = 0}^{k - 1} f_j \sum_{t = 0}^j {j \brace t} i^{\underline t}

将所有求和外移，可以得到：\displaystyle\sum_{i \geq 0} \sum_{j = 0}^{k - 1} \sum_{t = 0}^j x^i f_j {j \brace t} i^{\underline t}

将求和顺序调整，可以得到：\displaystyle\sum_{t = 0}^{k - 1 } \sum_{j = t}^{k - 1} f_j {j \brace t} \sum_{i \geq 0} x^i i^{\underline t}

那么现在就需要对 \displaystyle\sum_{i \geq 0} x^i i^{\underline t} 求生成函数。

由于当 i < t 时，下降幂会乘到 0，对答案没有贡献，那么可以将原式改写成 \displaystyle\sum_{i \geq t} x^i i^{\underline t}。

我们已知 f'_i = (i + 1)f_{i + 1}，那么对 1 + x + x^2 + x^3 + \dots + x^t + \dots 求一次导就等于 1 + 2x + 3x^2 + \dots + tx^t + \dots，求两次导就等于 2 + 6x + 12x^2 + \dots + t(t + 1)x^t，可以发现求 t 次导就变成 1^{\overline t} + 2^{\overline t}x + 3^{\overline t}x^2 + \dots + t^{\overline t}x^t，将其乘以 x^t，并将上升幂转为下降幂，可以得到 t^{\underline t}x^t + (t + 1)^{\underline t}x^{t + 1} + (t + 2)^{\underline t}x^{t + 2} + \dots + (2t)^{\underline t}x^{2t} + \dots = \displaystyle\sum_{i \geq t} x^i i^{\underline t}。有点凑的感觉，但是如果掌握了多项式求导，那么还是可以算出来的。根据之前求出的 1 + x + x^2 + x^3 + \dots + x^t + \dots 的生成函数，\displaystyle\sum_{i \geq t} x^i i^{\underline t} 的生成函数就是 x^t(\displaystyle\frac{1}{1 - x})^{(t)}。

现在再切换回复合函数的求导，那么该生成函数就为 \displaystyle\frac{t!x^t}{(1 - x)^{t + 1}} = \frac{(-1)^{\underline t}x^t}{(1 - x)^{t + 1}}。

那么原式就可以成 \displaystyle\sum_{t = 0}^{k - 1 } \frac{x^t(-1)^{\underline t}}{(1 - x)^{t + 1}} \sum_{j = t}^{k - 1} f_j {j \brace t}，此时式子已经与 i 无关，是原函数的封闭形式。

递推数列的生成函数（解递归式）

一个经典的递推函数就是斐波那契数列，它满足 f_i = \begin{cases} 1 & i = 1 或 2 \\ f_{i - 1} + f_{i - 2} & i > 2\end{cases}，它的通项公式是可以求出来的，主要有生成函数与线性代数两种解法，这里只讨论生成函数的做法。

首先可以写出斐波那契数列的前几项：\{1, 1, 2, 3, 5, 8, \dots\}，然后写出它的生成函数：1 + x + 2x^2 + 3x^3 + 5x^4 + 8x^5 + \dots。

设这个生成函数为 A，那么 xA = x + x^2 + 2x^3 + 3x^4 + 5x^5 + 8x^6 + \dots，x^2A = x^2 + x^3 + 2x^4 + 3x^5 + 5x^6 + 8x^7 + \dots，将两个式子加在一起可以得到 xA + x^2A = x + 2x^2 + 3x^3 + 5x^4 + 8x^5 + \dots，再用 A 减掉这个式子就可以得到 A - xA - x^2A = 1，因此 A = \displaystyle\frac{1}{1 - x - x^2}。

这个式子特别复杂，考虑能否把这个生成函数转为之前求出的几个简单的生成函数的和，因为这些生成函数所对应数列的通项公式很容易知道。

将分母转化为 (1 - \phi_1 x)(1 - \phi_2 x)，其中 \phi_1 = \displaystyle\frac{1 + \sqrt 5}{2}, \phi_2 = \displaystyle\frac{1 - \sqrt 5}{2}，将原式裂项，变成 \displaystyle\frac{\displaystyle\frac{\phi_1}{\sqrt 5}}{(1 - \phi_1 x)} - \frac{\displaystyle\frac{\phi_2}{\sqrt 5}}{(1 - \phi_2 x)}，这时，我们发现 \displaystyle\frac{1}{(1 - \phi_1 x)} 所对应的数列是 \{1, \phi_1, \phi_1^2, \phi_1^3, \dots\}，这个数列的第 n 项为 \phi_1^{n - 1}；\displaystyle\frac{1}{(1 - \phi_2 x)} 所对应的数列是 \{1, \phi_2, \phi_2^2, \phi_2^3, \dots\}，这个数列的第 n 项为 \phi_2^{n - 1}，那么原数列，也就是斐波那契数列的第 n 项就为 \displaystyle\frac{\phi_1}{\sqrt 5} \phi_1^{n - 1} - \displaystyle\frac{\phi_2}{\sqrt 5} \phi_2^{n - 1} = \frac{1}{\sqrt 5}(\phi_1^n - \phi_2^n) = \frac{1}{\sqrt 5}[(\frac{1 + \sqrt 5}{2})^n - (\frac{1 - \sqrt 5}{2})^n]。

我们将其进行推广，考虑常系数齐次线性递推的式子：f_n = \displaystyle\sum_{i = 1}^k c_i f_{n - i}，我们将 f 数列的生成函数 F(x) 表示出来

\begin{aligned} F(x) &= \displaystyle\sum_{i \geq 0} f_i x^i \\&= \sum_{i \geq 0} \sum_{j = 1}^k c_j f_{i - j} x^i \\&= \sum_{j = 1}^k c_j \sum_{i \geq 0} f_{i - j} x^i \\&= \sum_{j = 1}^k c_j \sum_{i \geq 0} f_i x^{i + j} \\&= \sum_{j = 1}^k c_j x^j F(x) \end{aligned}

这是一个方程，我们把它解出来，可以得到 F(x) = \displaystyle\frac{1}{1 - \sum\limits_{i = 1}^k c_ix^i}，这个是 F(x) 的封闭形式，它是一个分式的形式，因此我们再次把问题推广，如果给定我们一个有理函数 R(x) = \displaystyle\frac{P(x)}{Q(x)}，此时我们怎么确定 x^n 项的系数呢？

和求斐波那契数列通项公式一样，我们希望能将 R(x) 转换成已知的数列的生成函数。这时，需要我们再次发扬凑的方法，其实数学也不老是严谨的证明，很多时候需要凭一点感觉。

我们考虑到生成函数一般是拿来做排列组合用的，因此我们把 \displaystyle\frac{1}{(1 - x)^k} = \sum_{i \geq 0} \binom{i + k - 1}{k - 1} x^i 拿来做基本框架，我们在两边同时乘以 a，可以得到：

\displaystyle\frac{a}{(1 - x)^k} = \sum_{i \geq 0} \binom{i + k - 1}{k - 1} a x^i

等式两边再同时乘以 \rho^n，根据第 5 个数列的生成函数的封闭形式，等式又可以变成:

\displaystyle\frac{a}{(1 - \rho x)^k} = \sum_{i \geq 0} \binom{i + k - 1}{k - 1} a \rho^i x^i

最后，再将组合数简化一下，最终可以得到：

\displaystyle\frac{a}{(1 - \rho x)^{m + 1}} = \sum_{i \geq 0} \binom{m + i}{i} a \rho^i x^i

这么一大坨式子。像是所有我们已知的简单数列生成函数的融合，我们再设 S(x) = \displaystyle\sum_{i = 0}^l \sum_{j \geq 0} \binom{m_i + j}{j} a_i \rho_i^j x^j = \sum_{i = 1}^n \frac{a_i}{(1 - \rho_i x)^{m_i + 1}}，其中 [x^n] = \displaystyle\sum_{i = 1}^n a_i \binom{m_i + n}{m_i} \rho_i^n。

为了能将我们求出的生成函数 R(x) 第 n 项的系数，我们希望能将 R(x) 表示成 S(x) + T(x) 的形式，其中 T(x) 是关于 x 的多项式函数。由于 T(x) 和 S(x) 的 x^n 项的系数都是好求的，那么我们可以快速得到 R(x) 的 x^n 项的系数。

部分分式展开法

我们这里先来回忆一下我们在平时数学考试中是如何裂项的，以此来证明 T(x) 并不存在。

我们考虑一个分式 F(x) = \displaystyle\frac{G(x)}{\prod\limits_{i = 1}^n (x - a_i)}（G(x) 的次数小于 n），其中分母已经被因式分解了，我们假设它裂项之后变成了 \displaystyle\sum_{i = 1}^n \frac{b_i}{x - a_i}。

如果我们把 b_i 看成未知数，这其实就是一个 n 元方程，我们只需要带入 n 个不同的 x，就可以解一个 n 元 n 次方程了。不过我们肯定希望这些方程越简单越好，套用拉格朗日插值和中国剩余定理的推导经验，我们希望第 i 个方程只和 b_i 有关系，于是我们同时将等式两边乘以 (x - a_1)，得到：\displaystyle\frac{G(x)}{\prod\limits_{i = 2}^n (x - a_i)} = b_1 + \sum_{i = 2}^n \frac{b_i}{x - a_i} (x - a_1)。

我们再取 x \rightarrow a_1，那么可以得到 b_1 = \displaystyle\frac{G(a_1)}{\prod\limits_{i = 2}^n (a_1 - a_i)}。此时对于所有 b_i 都进行如下操作，那么可以得到 b_i = G(a_i) \displaystyle\prod_{j = 1 \wedge j \neq i}^n \frac{1}{a_i - a_j}，此时我们发现 b_i 一定存在，且是一个常数，那么 R(x) 一定可以表示成 \displaystyle\sum_{i = 1}^n \frac{b_i}{x - a_i}。我们将分式上下同时除以 -a_i，可以得到 R(x) = \displaystyle\sum_{i = 1}^n \frac{-\frac{b_i}{a_i}}{1 - \frac{1}{a_i} x}，我们发现这就是 S(x) 的形式，因此一个有理函数 P(x) 一定可以表示成 S(x) 的形式，从而 T(x) 并不存在。

但是，如果存在 a_i = a_j，那么此时 b_i 就求不出来，那么我们分解的方式就有所变化。这在后面讲解一般的有理展开定理会讲到。我们在这里只用知道 R(x) 依然可以裂项成 S(x) 的形式就可以了。

此时我们得到了一个重要的结论，那就是 R(x) = S(x)，但是每次都将 R(x) 裂项太麻烦了，能不能有更加简便的方法呢？这就是下面我们要来介绍的解递归式中最重要，也是最难的有理展开定理。

不同根的有理展开定理

这里我们先讲解不同根的有理展开定理，我们依然给出有理函数 R(x) = \displaystyle\frac{P(x)}{Q(x)}，我们将 Q(x) 因式分解成 q_0 \displaystyle\prod_{i = 0}^l (1 - \rho_i x)（其中 \rho_i 互不相等，也就是说，S(x) 中所有 m_i 都 = 0），若 P(x) 是次数小于 l 的多项式，那么有结论 [x^n]R(x) = \displaystyle\sum_{i = 1}^l a_i \rho_i^n，其中 a_i = \displaystyle\frac{- \rho_i P \left( \frac{1}{\rho_i}\right)}{Q'\left( \frac{1}{\rho_i}\right)}。

证明：

由于此时 S(x) 中所有 m_i 都为 0，那么 S(x) 可以简化成 \displaystyle\sum_{i = 1}^l \frac{a_i}{1 - \rho_i x}。

设 \alpha_i = \displaystyle\frac{1}{\rho_i}。我们可以利用部分分式展开的思想，将 R(x) = \displaystyle\frac{P(x)}{Q(x)} 上下同乘 1 - \rho_i x = \alpha_i - x，可以得到 (\alpha_i - x) R(x) = (\alpha_i - x) S(x)，那么 (x - \alpha_i) R(x) = (x - \alpha_i) S(x)。

我们在两边同时取极限 x \rightarrow \alpha_i，得到 \displaystyle\lim_{x \rightarrow \alpha_i} (x - \alpha_i) R(x) = \lim_{x \rightarrow \alpha_i} (x - \alpha_i) S(x)。

我们先来计算 \displaystyle\lim_{x \rightarrow \alpha_i} (x - \alpha_i) S(x)，对于 j \neq i，\displaystyle\lim_{x \rightarrow \alpha_i} (x - \alpha_i) \frac{a_j}{1 - \rho_i x} 一定为 0，而由于 (1 - \rho_i x) = -\rho_i(x - \alpha_i)，那么对于 j = i 这一项，取极限以后这一项为 \displaystyle\frac{-a_i}{\rho_i}。

我们再来计算 \displaystyle\lim_{x \rightarrow \alpha_i} (x - \alpha_i) R(x)。我们将 R(x) 写成 \displaystyle\frac{P(x)}{Q(x)}，于是我们要求 \displaystyle\lim_{x \rightarrow \alpha_i} (x - \alpha_i) \frac{P(x)}{Q(x)}，我们将这个式子改写成 P(\alpha_k) \displaystyle\lim_{x \rightarrow \alpha_i} \frac{x - \alpha_i}{Q(x)}，根据洛必达法则，\displaystyle\lim_{x \rightarrow \alpha_i} \frac{x - \alpha_i}{Q(x)} = \frac{1}{Q'(\alpha_i)}，那么右边的极限就是 \displaystyle\frac{P(\alpha_i)}{Q'(\alpha_i)}。

于是我们现在知道了 \displaystyle\frac{-a_i}{\rho_i} = \frac{P(\alpha_i)}{Q'(\alpha_i)}，于是 a_i = \displaystyle\frac{- \rho_i P \left( \frac{1}{\rho_i}\right)}{Q'\left( \frac{1}{\rho_i}\right)}。而由于所有 m_i = 0，因此 \displaystyle\binom{m_i + n}{m_i} = 1，那么最终 [x^n] R(x) = \displaystyle\sum_{i = 1}^l a_i \rho_i^n，证明成立。

一般的有理展开定理

现在开始讲一般化的有理展开定理。对于有理函数 R(x) = \displaystyle\frac{P(x)}{Q(x)}，其中 Q(x) = q_0 \displaystyle\prod_{i = 1}^l (1 - \rho_i x)^{d_i}，\rho_1, \rho_2, \dots, \rho_l 互不相等，而 P(x) 是一个次数小于 d_1 + d_2 + \dots + d_l 的多项式，那么 [x^n]R(x) = \displaystyle\sum_{i = 1}^l f_i(n) \rho_i^n，其中 f_i(n) 是一个次数为 d_i - 1 且首项系数为 a_i = \displaystyle\frac{(-\rho_i)^{d_i} P \left(\frac{1}{\rho_i} \right) d_i}{Q^{(d_i)} \left(\frac{1}{\rho_i} \right)} = \frac{P \left(\frac{1}{\rho_i} \right)}{(d_i - 1)! q_0 \prod\limits_{j = 1 \wedge j \neq i}^l \left(1 - \frac{\rho_j}{\rho_i} \right)^{d_j}}。

证明：

考虑数学归纳法。

我们现在已经知道 \max(d_1, d_2, \dots, d_l) = 1 时该定理成立。假设对于 \max(d_1, d_2, \dots, d_l) \leq m，该定理都成立，我们要证明对于 \max(d_1, d_2, \dots, d_l) = m + 1 时，该定理依然成立。

首先，如果有某个 d_i \leq m，那么根据归纳假设，f_i(n) 一定满足以上条件，因此我们只考虑 d_i = m + 1 的 i。不失一般性，假设 d_1 = m + 1，那我们将 (1 - \rho_1 x)^{d_1} 提出来，那么 Q(x) = (1 - \rho_1 x)^{d_1} \widetilde Q(x)，其中 \widetilde Q(x) = q_0 \displaystyle\prod_{i = 2}^l (1 - \rho_i x)^{d_i}。

我们将 R(x) 从新写成 \displaystyle\frac{P(x)}{(1 - \rho_1 x)^{d_1} \widetilde Q(x)}，再仿照部分分式展开的求解方法，假设我们将分母中 (1 - \rho_1 x)^{d_1} 提出来后，原式变成了 \displaystyle\sum_{i = 1}^{d_1} \frac{c_i}{(1 - \rho_1 x)^k} + \frac{\widetilde P(x)}{\widetilde Q(x)}（\widetilde P(x) 的次数小于 \widetilde Q(x)），现在我们要证明 c_i 是存在的。

我们在等式两边同乘 (1 - \rho_1 x)^{d_1}，我们可以得到 \displaystyle\frac{P(x)}{\widetilde Q(x)} = \sum_{i = 1}^{d_1} c_i (1 - \rho_1 x)^{d_1 - i} + (1 - \rho_i x)^{d_i} \frac{\widetilde P(x)}{\widetilde Q(x)}。此时，我们取 x \rightarrow \displaystyle\frac{1}{\rho_1}，那么等式右边除了 c_{d_1} 项以外，所有项都趋近于 0，因此 c_{d_1} = \displaystyle\frac{P \left(\frac{1}{\rho_1} \right)}{\widetilde Q \left(\frac{1}{\rho_1} \right)}。此时，我们再对于原式两边同时求导，将最高此项系数不断降低，就可以求出所有 c_i。

我们现在希望求出 \displaystyle\sum_{i = 1}^{d_1} \frac{c_i}{(1 - \rho_1 x)^k} + \frac{\widetilde P(x)}{\widetilde Q(x)} 的 x^n 项的系数。根据 S(x) 的 x^n 项的系数，此时 \displaystyle\sum_{i = 1}^{d_1} \frac{c_i}{(1 - \rho_1 x)^k} 的 x^n 项的系数为 \displaystyle\sum_{i = 1}^{d_1} c_i \binom{n + i - 1}{i - 1} \rho_i^n，此时我们发现 c_i \displaystyle\binom{n + i - 1}{i - 1} 是一个 i - 1 次多项式，那么 \displaystyle\sum_{i = 1}^{d_1} c_i \binom{n + i - 1}{i - 1} 就是一个 d_1 - 1 次多项式，这就是 f_1(n)。而根据归纳假设，\displaystyle\frac{\widetilde P(x)}{\widetilde Q(x)} 的 x^n 项的系数为 \displaystyle\sum_{i = 2}^l f_i(n) \rho_i^n。那么 [x^n] R(x) = \displaystyle\sum_{i = 1}^l f_i(n) \rho_i^n。

此时还差最后一个任务，那就是求出 f_1(n) 的首项系数。因此我们将 \displaystyle\frac{c_{d_1}}{(1 - \rho_1)^{d_1}} 写成 S(x) 的形式 c_{d_1} \displaystyle\binom{n + d_1 - 1}{d_1 - 1} \rho_1^{n}，它的首项为 c_{d_1} \displaystyle\frac{n^{d_1 - 1}}{(d_1 - 1)!} \rho_1^n，而 c_{d_1} = \displaystyle\frac{P \left(\frac{1}{\rho_1} \right)}{\widetilde Q \left(\frac{1}{\rho_1} \right)} = \frac{P \left(\frac{1}{\rho_1} \right)}{q_0 \prod\limits_{j = 2}^l \left(1 - \frac{\rho_j}{\rho_1}\right)^{d_j}}。由于我们已经将 \rho_i^n 提出去了，于是 f_1(n) 的首项系数为 \displaystyle\frac{c_{d_1}}{(d_1 - 1)!} = \frac{P \left(\frac{1}{\rho_1} \right)}{(d_1 - 1)! q_0 \prod\limits_{j = 2}^l \left(1 - \frac{\rho_j}{\rho_1} \right)^{d_j}}。

所有证明均符合数学归纳法，定理成立。

P4451 [国家集训队] 整数的lqp拆分

这道题目相对简单，因为我们已经求出了斐波那契数列的通项公式，这可以大大简化此题。

我们设 dp_{i, j} 表示 m = i 且 n = j 时的答案，于是我们可以很轻松得到递推式 dp_{i, j} = \displaystyle\sum_{k = 0}^j dp_{i - 1, k} \times F(j - k)，那么最终的答案可以写成 \displaystyle\sum_{i = 1}^n dp_{i, n}。

我们发现 k + (j - k) 是一个定值，因此这是一个卷积的形式，于是 dp_i 可以看做 dp_{i - 1} 与 F 数列的卷积，那么从 dp_0 开始，一共卷了 i 次，于是答案就是 \displaystyle\sum_{i = 1}^n [x^n] F(x)^i = [x^n] \sum_{i = 1}^n F(x)^i。

继续发扬人类智慧，我们考虑当 i > n 时，此时除了常数项，F 其他项至少都是从 n + 1 开始，因此 x^n 项的系数一定为 0，于是我们直接更改求和的上界，变成 [x^n] \displaystyle\sum_{i = 1}^{\infty} F(x)^i，由于 F(x)^0 = 1，不会影响 x^n 项的系数，因此我们再次改写成 [x^n] \displaystyle\sum_{i = 0}^{\infty} F(x)^i，我们发现这又是一个生成函数的形式，可以写成 [x^n] \displaystyle\frac{1}{1 - F(x)}，而 F(x) 的封闭形式我们已经求出了，于是式子再次化简成 \displaystyle\frac{1 - x - x^2}{1 - 2x - x^2}。

直接通过有理展开定理，将 Q(x) 拆成 (1 - (1 + \sqrt 2)x)(1 - (1 - \sqrt 2)x)，那么此时 x^n 项的系数就是 a_1 (1 + \sqrt 2)^n + a_2 (1 - \sqrt 2)^n，其中 a_1 = \displaystyle\frac{-(1 + \sqrt 2) P \left(\frac{1}{1 + \sqrt 2} \right)}{Q' \left(\frac{1}{1 + \sqrt 2} \right)} = \frac{-(1 + \sqrt 2) \times \frac{3 + 2 \sqrt 2}{7 + 5 \sqrt 2}}{\frac{-4 - 2 \sqrt 2}{1 + \sqrt 2}} = \frac{\sqrt 2 + 1}{2 \sqrt 2}，a_2 = \displaystyle\frac{-(1 - \sqrt 2) P \left(\frac{1}{1 - \sqrt 2} \right)}{Q' \left(\frac{1}{1 - \sqrt 2} \right)} = \frac{-(1 - \sqrt 2) \times \frac{3 + 8 \sqrt 2}{7 - 5 \sqrt 2}}{\frac{-4 + 4 \sqrt 2}{1 - \sqrt 2}} = \frac{\sqrt 2 - 1}{2 \sqrt 2}，于是答案就是 \displaystyle\frac{\sqrt 2 - 1}{2 \sqrt 2}(1 + \sqrt 2)^n + \frac{\sqrt 2 + 1}{2 \sqrt 2}(1 - \sqrt 2)^n。

现在，我们用求出 2 在模 10^9 + 7 的二次剩余 ~~（其实我不会，直接贺题解）~~，是 59713600，那么原式直接化简成 485071604 \times 940286408^n + 514928404 \times 59713601^n。根据费马小定理，只要边读边对 10^9 + 6 取模即可，然后这到题目就做完了。

其实这到题目的式子比较简单，有更简单的方法求出该式子，这个放到后面再讲。

:::info[点击查看代码]

#include <bits/stdc++.h>
using namespace std;
#define int long long
const int MOD = 1e9 + 7;
int read(){
    int x = 0, f = 1;
    char ch = getchar();
    while(!isdigit(ch)){
        if(ch == '-') 
            f = -1;
        ch = getchar();
    }
    while(isdigit(ch)){
        x = ((x << 1) + (x << 3) + (ch ^ 48)) % (MOD - 1);
        ch = getchar();
    }
    return x * f;
}
void write(int x){
    if(x < 0){
        putchar('-');
        x = -x;
    }
    if(x > 9)
        write(x / 10);
    putchar(x % 10 + '0');
}
int qpow(int a, int b){
    int res = 1;
    while(b > 0){
        if(b & 1)
            res = res * a % MOD;
        a = a * a % MOD;
        b >>= 1;
    }
    return res;
}
int n;
signed main(){
    n = read();
    write((485071604 * qpow(940286408, n - 1) % MOD + 514928404 * qpow(59713601, n - 1) % MOD) % MOD);
    return 0;
}

:::

还原递推式

这里介绍一下如何将生成函数转化回递推式。首先先来讲解一下不是封闭形式的情况。

我们先考虑一个最简单的形式，那就是告诉了你数列 g 的生成函数 G(x) = \displaystyle\sum_{i \geq 0} g_i x^i，让你求 g_n 的递推式。这太显然了，递推式就是 g_n = \displaystyle\sum_{i \geq 0} [n = i] g_i。

现在，我们将 G(x) 再乘以一个多项式 F(x) = \displaystyle\sum_{i = 0}^m f_i x^i，那么此时 g_n 的递推式又是什么呢？我们将 G(x) 也看成一个多项式，根据多项式的乘法，我们知道了 x^n 项的系数等于将 n 拆成两部分 n_1 和 n_2，将 F 中 x^{n_1} 项的系数乘以 G 中 x^{n_2} 项的系数，于是我们最终得到 f_n = \displaystyle\sum_{i = 0}^m f_i g_{n - i}。

我们再将其扩展，如果我们已知 G'(x) 与 F(x) 的乘积，那么此时 g_n 的递推式又是什么呢？根据多项式求导法则，G'(x) = \displaystyle\sum_{i \geq 0} (i + 1) g_{i + 1} x^i，再根据我们刚才推导的式子，最终我们会得到 g_n = \displaystyle\sum_{i = 0}^m f_i (n - i + 1) g_{n - i + 1}。

接着讲一下封闭形式的情况，设有理函数 R(x) = \displaystyle\frac{P(x)}{Q(x)}，我们先将其看成 P(x) 和 \displaystyle\frac{1}{Q(x)} 的卷积，而我们之前又求出一个常系数其次线性递推的生成函数为 \displaystyle\frac{1}{1 - \displaystyle\sum_{i = 1}^k c_ix^i}，那么我们可以很轻松求出分子对应的递推式，再根据上面讲解的整式的情况求解即可。

生成函数的应用（卷积）

考虑数列 h_0, h_1, h_2, \dots, h_n, \dots，其中 h_n 表示 e_1 + e_2 + \dots + e_k = n 的非负整数解的数目，那么这就变成了将 n 个不可区分小球放入 k 个可区分的盒子，每个盒子可以为空的方案数。

考虑到先给每个盒子里装 1 个小球，那么每个盒子就不能为空了，此时有 n + k 个物品了，考虑插板法，那么就是在 n + k - 1 个空中选 k - 1 个，方案数为 \displaystyle\binom{n + k - 1}{k - 1}，此时这个数列的生成函数为 \displaystyle\sum_{i \geq 0} \binom{i + k - 1}{k - 1}x^i = \frac{1}{(1 - x)^k}。

此时你会发现，\displaystyle\frac{1}{(1 - x)^k} = \underbrace{\frac{1}{1 - x} \times \frac{1}{1 - x} \times \dots \times \frac{1}{1 - x}}_{k - 1 个} = (1 + x + x^2 + \dots)(1 + x + x^2 + \dots)(1 + x + x^2 + \dots) = (\displaystyle\sum_{e_1 \geq 0} x^{e_1})(\sum_{e_2 \geq 0} x^{e_2}) \dots (\sum_{e_k \geq 0} x^{e_k})，于是，单个物品能取到的数所形成数列的生成函数之积等于所有物品能取到的数所形成数列的生成函数。

这个性质再求某些带限制的背包问题中有非常大的用处，下面举两个例子。

P10780 BZOJ3028 食物

各种食物的生成函数如下：

承德汉堡：\{1, 0, 1, 0, 1, 0, \dots \} = 1 + x^2 + x^4 + x^6 + \dots = \displaystyle\frac{1}{1 - x^2}；
可乐：\{1, 1, 0, 0, 0, 0, \dots \} = 1 + x；
鸡腿：\{1, 1, 1, 0, 0, 0, \dots \} = 1 + x + x^2；
蜜桃多：\{0, 1, 0, 1, 0, 1, \dots \} = x + x^3 + x^5 + x^7 + \dots = \displaystyle\frac{x}{1 - x^2}；
鸡块：\{1, 0, 0, 0, 1, 0, \dots \} = 1 + x^4 + x^8 + x^{12} + \dots = \displaystyle\frac{1}{1 - x^4}；
包子：\{1, 1, 1, 1, 0, 0, \dots \} = 1 + x + x^2 + x^3；
土豆片炒肉：\{1, 1, 0, 0, 0, 0, \dots \} = 1 + x；
面包：\{1, 0, 0, 1, 0, 0, \dots \} = 1 + x^3 + x^6 + x^9 + \dots = \displaystyle\frac{1}{1 - x^3}。

全部乘起来为 \displaystyle\frac{x}{1 - x^4}，它是 \displaystyle\sum_{i \geq 0} x^{i + 1} \binom{i + 3}{3} = \sum_{i \geq 1} x^i \binom{i + 2}{3} 的封闭形式，于是 \displaystyle\binom{n + 2}{3} 就是这个问题的答案。

:::info[点击查看代码]

#include <bits/stdc++.h>
using namespace std;
#define int long long
const int MOD = 10007;
int extend_gcd(int a, int b, int &x, int &y){
    if(b == 0){
        x = 1;
        y = 0;
        return 0;
    }
    int d = extend_gcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}
int mod_inverse(int a, int m){
    int x, y;
    extend_gcd(a, m, x, y);
    return (x % m + m) % m;
}
int n;
char c;
signed main(){
    while(isdigit(c = getchar()))
        n = (n * 10 + c - '0') % MOD;
    int inv = mod_inverse(6, MOD);
    printf("%lld", n * (n + 1) * (n + 2) * inv % MOD);
    return 0;
}

:::

P4389 付公主的背包

假设每个物品体积为 1，于是每个物品只能取 v_i 倍数个，于是每个物品的生成函数就是 \displaystyle\frac{1}{1 - x^{v_i}}，总的方案数序列的生成函数就是 \displaystyle\prod_{i = 1}^n \frac{1}{1 - x^{v_i}}。

下面就考虑各位的推式子能力了，由于求积不太能计算，于是考虑先多项式 \ln 再多项式 \exp，那么求积就变成了求和：\exp \left[\displaystyle\sum_{i = 1}^n \ln \left(\frac{1}{1 - x^{v_i}} \right) \right]。

将 \ln \left(\displaystyle\frac{1}{1 - x^{v_i}} \right) 先求导再做积分，根据复合函数求导，原式变成了 \exp \left[\displaystyle\sum_{i = 1}^n \int \left(\frac{1}{1 - x^{v_i}} \right)'(1 - x^{v_i}) \mathrm d x \right]。

将式子化简，可以得到 \exp \left[\displaystyle\sum_{i = 1}^n \int \frac{v_i x^{v_i - 1}}{1 - x^{v_i}} \mathrm d x \right]。

发现 \displaystyle\frac{v_i x^{v_i - 1}}{1 - x^{v_i}} 和等比数列的和很像，展开后变成 \exp \left(\displaystyle\sum_{i = 1}^n \int \sum_{j \geq 1}v_i x^{jv_i - 1} \mathrm d x \right)。

根据幂函数的积分，原式最终等于 \exp \left(\displaystyle\sum_{i = 1}^n \sum_{j \geq 0} \frac{x^{jv_i}}{j} \right)。

里面的多项式是好求的，再套个 \exp 就可以了。

:::info[点击查看代码]

#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N = 1e6 + 9, MOD = 998244353;
int rev[N], in[N];
void init(){
    in[1] = 1;
    for(int i = 2; i <= N; i++)
        in[i] = (MOD - MOD / i) * in[MOD % i] % MOD;
}
int qpow(int x, int y){
    int res = 1;
    while(y){
        if(y & 1)
            res = res * x % MOD;
        x = x * x % MOD;
        y >>= 1;
    }
    return res;
}
void NTT(int *f, int n, int opt){
    for(int i = 0; i < n; i++){
        rev[i] = rev[i >> 1] >> 1;
        if(i % 2 == 1)
            rev[i] |= n >> 1;
    }
    for(int i = 0; i < n; i++)
        if(i < rev[i])
            swap(f[i], f[rev[i]]);
    for(int h = 2; h <= n; h <<= 1){
        int mi = h >> 1;
        int gn = qpow(3, (MOD - 1) / h);
        for(int j = 0; j < n; j += h){
            int g = 1;
            for(int k = 0; k < mi; k++){
                int tmp = f[j + k + mi] * g % MOD;
                f[j + k + mi] = (f[j + k] - tmp + MOD) % MOD;
                f[j + k] = (f[j + k] + tmp) % MOD;
                g = g * gn % MOD;
            }
        }
    }
    if(opt == -1){
        reverse(f + 1, f + n);
        int inv = qpow(n, MOD - 2);
        for(int i = 0; i < n; i++)
            f[i] = f[i] * inv % MOD;
    }
}
int tmp[N << 1];
void inv(int *f, int *g, int n){
    if(n == 1){
        g[0] = qpow(f[0], MOD - 2);
        return;
    }
    inv(f, g, (n + 1) >> 1);
    int lim = 1;
    while(lim < (n << 1))
        lim <<= 1;
    for(int i = 0; i < n; i++)
        tmp[i] = f[i];
    for(int i = n; i < lim; i++)
        tmp[i] = 0;
    NTT(tmp, lim, 1);
    NTT(g, lim, 1);
    for(int i = 0; i < lim; i++)
        g[i] = (2 - g[i] * tmp[i] % MOD + MOD) % MOD * g[i] % MOD;
    NTT(g, lim, -1);
    for(int i = n; i < lim; i++)
        g[i] = 0;
}
int f2[N << 1];
void ln(int *f, int *g, int n){
    for(int i = 0; i < (n << 2); i++)
        g[i] = 0;
    inv(f, g, n);
    int lim = 1;
    while(lim < (n << 1))
        lim <<= 1;
    for(int i = 0; i < n - 1; i++)
        f2[i] = f[i + 1] * (i + 1) % MOD;
    for(int i = n - 1; i < lim; i++)
        f2[i] = 0;
    NTT(f2, lim, 1);
    NTT(g, lim, 1);
    for(int i = 0; i < lim; i++)
        g[i] = g[i] * f2[i] % MOD;
    NTT(g, lim, -1);
    for(int i = n - 1; i > 0; i--)
        g[i] = g[i - 1] * in[i] % MOD;
    for(int i = n; i < lim; i++)
        g[i] = 0;
    g[0] = 0;
}
int lng[N << 1];
void exp(int *f, int *g, int n){
    if(n == 1){
        g[0] = 1;
        return;
    }
    exp(f, g, (n + 1) >> 1);
    ln(g, lng, n);
    int lim = 1;
    while(lim < (n << 1))
        lim <<= 1;
    for(int i = 0; i < n; i++){
        if(f[i] >= lng[i])
            lng[i] = f[i] - lng[i];
        else
            lng[i] = f[i] - lng[i] + MOD;
    }
    for(int i = n; i < lim; i++)
        lng[i] = g[i] = 0;
    lng[0]++;
    NTT(lng, lim, 1);
    NTT(g, lim, 1);
    for(int i = 0; i < lim; i++)
        g[i] = g[i] * lng[i] % MOD;
    NTT(g, lim, -1);
    for(int i = n; i < lim; i++)
        g[i] = 0;
}
int v[N], cnt[N], f[N], g[N], n, m;
signed main(){
    init();
    scanf("%lld%lld", &n, &m);
    for(int i = 1; i <= n; i++){
        scanf("%lld", &v[i]);
        cnt[v[i]]++;
    }
    for(int i = 1; i <= m; i++)
        if(cnt[i])
            for(int j = i; j <= m; j += i)
                f[j] = (f[j] + cnt[i] * in[j / i]) % MOD;
    exp(f, g, m + 1);
    for(int i = 1; i <= m; i++)
        printf("%lld\n", g[i]);
    return 0;
}

:::

生成函数技巧

有理展开定理的式子真是一大坨，因此对于某些特别的生成函数，我们有简单的方法解递归式，我们来看一看。

特征多项式

对象：分母容易因式分解的有理函数

P4451 [国家集训队] 整数的lqp拆分

这种方法也就是我们之前求斐波那契数列时用的方法，也算是有理展开定理的一种简单应用。

我们目前求出了答案序列的生成函数为 \displaystyle\frac{1 - x - x^2}{1 - 2x - x^2}。根据还原递推式的方法，我们知道分母对应的递推式为 a_n = 2 a_{n - 1} + a_{n - 2}，再将其和分子乘在一起。得到最终的递推式为 ans = a_n - a_{n - 1} - a_{n - 2} = a_{n - 1}。

接下来我们将转移写成矩阵形式

\begin{bmatrix} 2 & 1\\ 1 & 0 \end{bmatrix} \begin{bmatrix} a_{n - 1}\\ a_{n - 2} \end{bmatrix} = \begin{bmatrix} a_n\\ a_{n - 1} \end{bmatrix}

我们求解出该矩阵的两个特征值 v_1 = \sqrt 2 + 1，v_2 = -\sqrt 2 + 1，根据我们在线性代数学习笔记（一）：线性代数基础知识中的推导，a_n 的表达式一定为 c_1 (\sqrt 2 + 1)^n + c_2 (-\sqrt 2 + 1)^n，随后带入 n = 0 以及 n = 1 可以求出 c_1 = \displaystyle\frac{\sqrt 2 - 1}{2 \sqrt 2}, c_2 = \frac{\sqrt 2 + 1}{2 \sqrt 2}，此时的通项公式和我们用有理展开定理推出来的一毛一样。

拆分法

对象：分子次数大于分母次数的有理函数

AT_abc241_h Card Deck Score

我们写出每一种卡牌 i 的生成函数 \displaystyle\sum_{j = 0}^{b_i} a_i^j x^j = \frac{1 - a_i^{b_i + 1} x^{b_i + 1}}{1 - a_i x}，那么我们将每种卡牌的生成函数乘在一起，根据之前对于卷积的讲解，此时该函数的 x^m 项的系数就是答案。

但是我们发现，对于整体的生成函数 \displaystyle\prod_{i = 1}^n \frac{1 - a_i^{b_i + 1} x^{b_i + 1}}{1 - a_i x}，它分母的次数竟然大于分母的次数，无法使用有理展开定理。此时我们有了一个想法，那就是将 \displaystyle\prod_{i = 1}^n 1 - a_i^{b_i + 1} x^{b_i + 1} 和 \displaystyle\prod_{i = 1}^n \frac{1}{1 - a_i x} 分别看成两个生成函数 F(x) 和 G(x)，那我们要求的就是 \displaystyle\sum_{j = 0}^m [x^j] F(x) [x^{m - j}] G(x)。

此时我们发现由于 n 非常的小，那么 F(x) 中只有很少的项，仿照二项式定理的推导，我们可以花费 \mathcal O(2^n) 的时间复杂度，枚举每一个因式是选 - a_i^{b_i + 1} x^{b_i + 1} 这一项还是 1 这一项，然后就可以算出贡献了。而 G(x) 就可以非常容易地用有理展开定理求出第 m 项的系数，然后就做完了。

:::info[点击查看代码]

#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N = 19, MOD = 998244353;
int qpow(int a, int b){
    int res = 1;
    while(b > 0){
        if(b & 1)
            res = res * a % MOD;
        a = a * a % MOD;
        b >>= 1;
    }
    return res;
}
int a[N], b[N], f[N], n, m, ans;
map <int, int> mp;
map <int, bool> flag;
signed main(){
    scanf("%lld%lld", &n, &m);
    for(int i = 1; i <= n; i++)
        scanf("%lld%lld", &a[i], &b[i]);
    for(int s = 0; s < (1 << n); s++){
        int c = 1, d = 0;
        for(int i = 0; i < n; i++){
            if(s & (1 << i)){
                c = c * (MOD - qpow(a[i + 1], b[i + 1] + 1)) % MOD;
                d += b[i + 1] + 1;
            }
        }
        mp[d] += c;
    }
    for(int i = 1; i <= n; i++){
        int p = 1;
        for(int j = 1; j <= n; j++){
            if(j == i)
                continue;
            p = p * (1 - a[j] * qpow(a[i], MOD - 2) % MOD + MOD) % MOD;
        }
        f[i] = qpow(p, MOD - 2);
    }
    for(int s = 0; s < (1 << n); s++){
        int d = 0;
        for(int i = 0; i < n; i++)
            if(s & (1 << i))
                d += b[i + 1] + 1;
        if(!flag[d] && d <= m){
            flag[d] = true;
            int sum = 0;
            for(int i = 1; i <= n; i++)
                sum = (sum + f[i] * qpow(a[i], m - d) % MOD) % MOD;
            ans = (ans + mp[d] * sum % MOD) % MOD;
        }
    }
    printf("%lld", ans);
    return 0;
}

:::

拉格朗日插值法

对象：分母次数较高的有理函数

AT_aising2020_f Two Snuke

我们先设

\begin{cases} a_1 + a_2 = n_1 & a_2 - a_1 = m_1\\ b_1 + b_2 = n_2 & b_2 - b_1 = m_2\\ c_1 + c_2 = n_3 & c_2 - c_1 = m_3\\ d_1 + d_2 = n_4 & d_2 - d_1 = m_4\\ e_1 + e_2 = n_5 & e_2 - e_1 = m_5 \end{cases}

（其中的 a, b, c, d, e 分别对应原问题中的 s, n, u, k, e）

我们再设 f_{i, j} 表示前 i 组数 \displaystyle\sum_{k = 1}^i n_k = j 时，所有可能的 \displaystyle\prod_{k = 1}^i m_k 的和，那么最终答案就是 \displaystyle\sum_{i = 1}^n f_{5, i}。

现在我们考虑第 i 组数，如果它们是 0 和 n_i，那么 m_i = n_i；如果它们是 1 和 n_i - 1，那么 m_i = n_i - 2……由此看来，m_i 构成了一个等差数列。根据乘法分配律，我们现在需要求出所有可能的 m_i 的和，因此直接等差数列求和，可以得到当 n_i 为奇数时，所有可能 m_i 的和为 \displaystyle\frac{(i + 1)^2}{4}；当 n_i 为偶数时，所有可能 m_i 的和为 \displaystyle\frac{i(i + 2)}{4}，我们设其为 g_{n_i}。

此时我们就可以得出 f_{i, j} 的递推式了：f_{i, j} = \displaystyle\sum_{k = 0}^j f_{i - 1, j - k} \times g_k，而且 f_{1, j} = g_{j}，那么我们仿照P4451 [国家集训队] 整数的lqp拆分的思路，将其看成两个生成函数的卷积，那么 f_{5, i} = [x^i] g^5。

我们现在写出 g 的生成函数。分开奇偶考虑，分别求出它们的生成函数，再乘在一起，可以得到最终的生成函数为 \displaystyle\frac{x}{(1 - x)^2 (1 - x^2)}，因此 g^5 = \displaystyle\frac{x^5}{(1 - x)^{10} (1 - x^2)^5}。而我们最终要求的是一段前缀和，因此还需要乘上一个 \displaystyle\frac{1}{(1 - x)}，变成 \displaystyle\frac{x^5}{(1 - x)^{11} (1 - x^2)^5}。

这里肯定可以使用有理展开定理，但是这可是相同根，式子特别复杂，因此我们仿照上一道题目的思路，将分子分母拆开，分子相当于做了 5 次前缀和，不用管它，只用求分母就可以了，此时再套用有理展开定理的式子，就会简化很多。于是最终 x^n 项的系数为 \displaystyle\sum_{k = 0}^{\lfloor \frac n2 \rfloor} \binom{k + 4}{k} \times \binom{n - 2k + 10}{10}。

此时我们发现 \displaystyle\binom{k + 4}{k} \times \binom{n - 2k + 10}{10} 是一个关于 k 的 14 次多项式，套上一个 \displaystyle\sum 就变成了 15 次多项式，因此直接拉格朗日差值即可。

:::info[点击查看代码]

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N = 19, MOD = 1e9 + 7;
int fac[N], sum[N], T, n;
int qpow(int a, int b){
    int res = 1;
    while(b > 0){
        if(b & 1)
            res = res * a % MOD;
        a = a * a % MOD;
        b >>= 1;
    }
    return res;
}
int binom(int n, int m){
    int res = 1;
    for(int i = 1; i <= m; i++)
        res = res * (n - i + 1) % MOD;
    return res * qpow(fac[m], MOD - 2) % MOD;
}
signed main(){
    fac[0] = 1;
    for(int i = 1; i <= 16; i++)
        fac[i] = fac[i - 1] * i % MOD;
    scanf("%lld", &T);
    while(T--){
        scanf("%lld", &n);
        if(n < 5){
            printf("0\n");
            continue;
        }
        n -= 5;
        int tmp = min(16ll, n / 2), res = 0;
        for(int j = 0; j <= tmp; j++)
            sum[j] = binom(j + 4, 4) * binom(n - 2 * j + 10, 10) % MOD;
        for(int j = 1; j <= tmp; j++)
            sum[j] = (sum[j] + sum[j - 1]) % MOD;
        n /= 2;
        if(n <= tmp){
            printf("%lld\n", sum[n]);
            continue;
        }   
        for(int i = 0; i <= tmp; i++){
            int s1 = sum[i], s2 = 1;
            for(int j = 0; j <= tmp; j++)
                if(j != i){
                    s1 = s1 * (n - j) % MOD;
                    s2 = s2 * (i - j) % MOD;
                }
            res = (res + s1 * qpow(s2, MOD - 2)) % MOD;
        }
        printf("%lld\n", res);
    }
    return 0;
}

:::

暴力方法

对象：相互递推数列

P6112 直接自然溢出啥事没有加强版

真的是特别暴力的方法，我从么有见过这么长的式子（但是总比有理展开定理简单）。

我们这里先来推一个式子。现在我们有方程 G^2 + F_1 G + F_2 = 0（G, F_1, F_2 均为连续可导函数），我们现在希望得到 G 函数的表达式（其实 5 次以下方程皆可，但那个太难写了，就不讲解了）。

我们希望能将 G^2 降次，来求出 G 的值，因此我们将原等式两边求导，根据求导法则 [f(x) g(x)]' = f'(x) g(x) + f(x) g'(x)，2GG' + F_1 G' + F_1' G + F_2'，于是 G' = -\displaystyle\frac{F_1 G' + F_2'}{2G + F_1}，因此 G' = -\left(\displaystyle\frac{F_1'}{2} + \frac{F_2' - \frac{F_1 F_1'}{2}}{2G + F_1} \right)。

这似乎降了个寂寞，于是我们考虑能否将分母也乘上 2G + F_1，这样就可以将 G 和 G' 分解开了。

于是我们将原等式左右两边同时除以 2G + F_1，然后大力推式子：

\begin{aligned} \displaystyle\frac{G^2 + F_1 G + F_2}{2G + F_1} &= 0\\ \frac G2 + \frac{F_1}{4} + \frac{F_2 - \frac{F_1^2}{4}}{2G + F_1} &= 0\\ \frac{1}{2G + F_1} &= -\frac{\frac G2 + \frac{F_1}{4}}{F_2 - \frac{F_1^2}{4}} = - \frac{2G + F_1}{4F_2 - F_1^2} \end{aligned}

此时我们成功将 2G + F_1 翻到了分子上，那么就可以带入式子 G' = -\left(\displaystyle\frac{F_1'}{2} + \frac{F_2' - \frac{F_1 F_1'}{2}}{2G + F_1} \right)，得到：

\begin{aligned} G' &= -\displaystyle\frac{F_1'}{2} + \frac{(F_2' - \frac{F_1 F_1'}{2})(2G + F_1)}{4F_2 - F_1^2}\\ (4F_2 - F_1^2) G' &= -\frac{F_1' (4F_2 - F_1^2)}{2} + (F_2' - \frac{F_1 F_1'}{2})(2G + F_1)\\ (4F_2 - F_1^2) G' &= -2F_1' F_2 + \frac{F_1' F_1^2}{2} + 2GF_2' + F_1 F_2' - G F_1 F_1' - \frac{F_1' F_1^2}{2}\\ (4F_2 - F_1^2) G' + (F_1 F_1' - 2F_2')G + 2F_1' F_2 - F_1 F_2' &= 0 \end{aligned}

我们将其扩展一下，设 F_1 G^2 + F_2 G + F_3 = 0，于是 G^2 + \displaystyle\frac{F_2}{F_1} G + \frac{F_3}{F_1} = 0，运用刚才推出来的式子，可以得到：

\begin{aligned} (\displaystyle\frac{4 F_3}{F_1} - \frac{F_2^2}{F_1^2})G' + (\frac{F_2}{F_1} \times \frac{F_2' F_1 - F_2 F_1'}{F_1^2} - 2 \times \frac{F_3' F_1 - F_3 F_1'}{F_1^2}) G + 2 \times \frac{F_2' F_1 - F_2 F_1'}{F_1^2} \times \frac{F_3}{F_1} - \frac{F_2}{F_1} \times \frac{F_3' F_1 - F_3 F_1'}{F_1^2} &= 0\\ (4 F_3 F_1^2 - F_1 F_2^2)G' + (F_1 F_2 F_2' - F_1' F_2^2 - 2 F_3' F_1^2 + 2 F_1 F_1' F_3)G + 2 F_1 F_2' F_3 - F_1' F_2 F_3 - F_1 F_2 F_3' &= 0 \end{aligned}

似乎感觉还是转换了个寂寞，还是存在 G' 和 G 两个未知多项式，但是在生成函数中，求导有组合意义，这在这道题中很有用处。

我们设 a_i, b_i, c_i, d_i, e_i, f_i 表示长度为 i 的序列是「语句」、「程序片段」、「语句块」、「函数」、「值」、「空串」的数量，现在开始逐个分析每个 DP 函数的递推式：

最简单的是 f_i，只有当字符串长度为 0 的时候这个串才是「空串」，因此 f_i = [i = 0]；
然后是 a_i。题目定义 “单个分号 ; 是一个「语句」”，我们其实可以将其看成一个「空串」加上一个 ;，这就和“如果字符串 A 是「值」，那么 A; 是「语句」”的定义一致了，由于增加一个 ; 会使字符串长度增加 1，那么 a_i = e_{i - 1} + f_{i - 1}，而又由于「语句块」是语句，因此 a_i = c_i + e_{i - 1} + f_{i - 1}；
然后是 c_i。它只有一个定义“如果字符串 A 是「程序片段」，则 {A} 是「语句块」”，由于增加一个 {} 会使字符串长度增加 2，因此 c_i = b_{i - 2}；
然后是 d_i。题目定义“如果字符串 A 是「语句块」，则 []A 和 []()A 都是「函数」”，“如果字符串 A 是「函数」，则 (A) 是「函数」”，由于增加一个 []、()、[]() 会使字符串长度分别增加 2, 2, 4，因此 d_i = c_{i - 2} + c_{i - 4} + d_{i - 2}；
然后是 b_i。题目定义“如果字符串 A 是「程序片段」，字符串 B 是「语句」，则 AB 是「程序片段」”，那么我们可以将当前长度为 i 的字符串拆成 j 和 i - j 两部分，前一部分为「程序片段」，后一部分为「语句」，再将「空串」考虑进来，因此 b_i = \displaystyle\sum_{j = 0}^i b_j a_{i - j} + f_i；
最后是 e_i。题目定义“如果字符串 A 是「函数」，A 和 A() 都是「值」”，因此首先 e_i = d_{i - 2} + d_i。接着考虑“如果字符串 A 是「值」，则 (A) 是「值」”，此时我们不能直接加上 e_{i - 2}，因为假设该递推式成立，那么 e_{i} = d_{i} + 2 \times d_{i - 2} + d_{i - 4} + e_{i - 4}，此时我们发现计算了两边 d_{i - 2} 的贡献，但是事实上 d_{i - 2} 只应该造成一次贡献，因此我们在递推式中还要再减掉一个 d_{i - 2}，于是 e_i = d_i + e_{i - 2}。

我们运用之前讲解的生成函数平移和卷积的含义，将所有递推式写成生成函数的形式（f_i 的生成函数 F(x) 就是 1，因此不单独列出来）：

\begin{cases} A = C + x E + x\\ B = AB + 1\\ C = x^2 B\\ D = x^2 C + x^4 C + x^2 D\\ E = D + x^2 E \end{cases}

我们最终要求的是 [x^n] B(x)，因此我们试着将 B 表示出来。

首先，C 可以完全被 B 表示出来，我们将其带入 D 的表达式，得到 D = x^6 B + x^4 B + x^2 D，D = \displaystyle\frac{x^6 + x^4}{1 - x^2} B。我们再将 D 带入 E 的表达式，得到 E = \displaystyle\frac{x^6 + x^4}{1 - x^2} B + x^2 E，E = \displaystyle\frac{x^4 + x^6}{x^4 - 2x^2 + 1} B。

我们把 E 再代入 A 的表达式，于是 A = x^2 B + \displaystyle\frac{x^7 + x^5}{x^4 - 2x^2 + 1} B + x，最终将 A 带入 B 的表达式，得到 B = x^2 B^2 + \displaystyle\frac{x^7 + x^5}{x^4 - 2x^2 + 1} B^2 + xB + 1。

我们将等式两边同时乘以 x^4 - 2x^2 + 1，得到 (x^4 - 2x^2 + 1) B = (x^6 - 2x^4 + x^2) B^2 + (x^7 + x^5) B^2 + (x^5 - 2x^3 + x)B + x^4 - 2x^2 + 1，整理可得 (x^7 + x^6 + x^5 - 2x^4 + x^2) B^2 + (x^5 - x^4 - 2x^3 + 2x^2 + x - 1) B + x^4 - 2x^2 + 1。

此时我们就可以运用刚才推出来的式子了，此处 F_1 = x^7 + x^6 + x^5 - 2x^4 + x^2, F_2 = x^5 - x^4 - 2x^3 + 2x^2 + x - 1, F_3 = x^4 - 2x^2 + 1。那么我们就可以得到以下超长等式：

\begin{aligned} &(4x^{18} + 7x^{17} + 5x^{16} - 20x^{15} - 33x^{14} + 5x^{13} + 55x^{12} + 42x^{11} - 67x^{10} - 63x^9 + 59x^8 + 40x^7 - 31x^6 - 13x^5 + 9x^4 + 2x^3 - x^2)B'\\ &+ (6x^{17} + 8x^{16} - 4x^{15} - 30x^{14} - 63x^{13} + 29x^{12} + 111x^{11} + 27x^{10} - 94x^9 - 78x^8 + 70x^7 + 64x^6 - 39x^5 - 23x^4 + 15x^3 + 3x^2 - 2x)B\\ &+ (-x^{15} + 3x^{14} + 11x^{13} - 15x^{12} - 30x^{11} + 22x^{10} + 40x ^9 - 6x^8 - 37x^7 - 9x^6 + 27x^5 + 5x^4 - 12x^3 + 2x) = 0 \end{aligned}

通过惊人的注意力以及高超的因式分解技巧，我们发现等式两边可以同时除以 x(x - 1)(x + 1)，得到：

\begin{aligned} &(4x^{15} + 7x^{14} + 9x^{13} - 13x^{12} - 24x^{11} - 8x^{10} + 31x^9 + 34x^8 - 36x^7 - 29x^6 + 23x^5 + 11x^4 - 8x^3 - 2x^2 + x)B'\\ &+ (6x^{14} + 8x^{13} + 2x^{12} - 22x^{11} - 61x^{10} + 7x^9 + 50x^8 + 34x^7 - 44x^6 - 44x^5 + 26x^4 + 20x^3 - 13x^2 - 3x + 2)B\\ &+ (-x^{12} + 3x^{11} + 10x^{10} - 12x^9 - 20x^8 + 10x^7 + 20x^6 + 4x^5 - 17x^4 - 5x^3 + 10x^2 - 2) = 0 \end{aligned}

根据我们之前讲解的还原递推式的技巧，我们求出 b_i 第 n 项的递推式为：

\begin{aligned} &n b_n - 2 (n - 1) b_{n - 1} - 8 (n - 2) b_{n - 2} + 11 (n - 3) b_{n - 3} + 23 (n - 4) b_{n - 4} - 29 (n - 5) b_{n - 5} - 36 (n - 6) b_{n - 6}\\ &+ 34 (n - 7) b_{n - 7} + 31 (n - 8) b_{n - 8} - 8 (n - 9) b_{n - 9} - 24 (n - 10) b_{n - 10} - 13 (n - 11) b_{n - 11} + 9 (n - 12) b_{n - 12}\\ &+ 7 (n - 13) b_{n - 13} + 4 (n - 14) b_{n - 14} + 2 b_n - 3 b_{n - 1} - 13 b_{n - 2} + 20 b_{n - 3} + 26 b_{n - 4} - 44 b_{n - 5} - 44 b_{n - 6} + 34 b_{n - 7}\\ &+ 50 b_{n - 8} + 7 b_{n - 9} - 61 b_{n - 10} - 22 b_{n - 11} + 2 b_{n - 12} + 8 b_{n - 13} + 6 b_{n - 14} - [n = 12] + 3 [n = 11] + 10 [n = 10]\\ &- 12 [n = 9] - 20 [n = 8] + 10 [n = 7] + 20 [n = 6] + 4 [n = 5] - 17 [n = 4] - 5[n = 3] + 10 [n = 2] - 2[n = 0] = 0 \end{aligned}

现在就可以递推出 b_n 了。

:::info[点击查看代码]

#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N = 1e7 + 9, MOD = 998244353;
int inv[N];
void init(){
    inv[1] = 1;
    for(int i = 2; i < N; i++)
        inv[i] = (MOD - MOD / i) * inv[MOD % i] % MOD;
}
struct Variable{
    int x;
    inline Variable operator + (const Variable &b) const{
        Variable res = Variable{x + b.x > MOD ? x + b.x - MOD : x + b.x};
        return res;
    }
    inline Variable operator + (const int &b) const{
        Variable res = Variable{x + b > MOD ? x + b - MOD : x + b};
        return res;
    }
    inline Variable operator - (const Variable &b) const{
        Variable res = Variable{x - b.x < 0 ? x - b.x + MOD : x - b.x};
        return res;
    }
    inline Variable operator - (const int &b) const{
        Variable res = Variable{x - b < 0 ? x - b + MOD : x - b};
        return res;
    }
    inline Variable operator * (const Variable &b) const{
        Variable res = Variable{x * b.x % MOD};
        return res;
    }
    inline Variable operator * (const int &b) const{
        Variable res = Variable{x * b % MOD};
        return res;
    }
    inline Variable operator / (const int &b) const{
        Variable res = Variable{x * inv[b] % MOD};
        return res;
    }
} f[N];
int n;
inline Variable F(int x){
    return x >= 0 ? f[x] : Variable{0};
}
signed main(){
    init();
    scanf("%lld", &n);
    for(int i = 0; i <= n; i++){
        f[i] = F(i - 1) * 2 * (i - 1) + F(i - 2) * 8 * (i - 2) - F(i - 3) * 11 * (i - 3) - F(i - 4) * 23 * (i - 4) + F(i - 5) * 29 * (i - 5) + F(i - 6) * 36 * (i - 6) - F(i - 7) * 34 * (i - 7) - F(i - 8) * 31 * (i - 8) + F(i - 9) * 8 * (i - 9) + F(i - 10) * 24 * (i - 10) + F(i - 11) * 13 * (i - 11) - F(i - 12) * 9 * (i - 12) - F(i - 13) * 7 * (i - 13) - F(i - 14) * 4 * (i - 14) 
        + F(i - 1) * 3 + F(i - 2) * 13 - F(i - 3) * 20 - F(i - 4) * 26 + F(i - 5) * 44 + F(i - 6) * 44 - F(i - 7) * 34 - F(i - 8) * 50 - F(i - 9) * 7 + F(i - 10) * 61 + F(i - 11) * 22 - F(i - 12) * 2 - F(i - 13) * 8 - F(i - 14) * 6 
        + (i == 12) - 3 * (i == 11) - 10 * (i == 10) + 12 * (i == 9) + 20 * (i == 8) - 10 * (i == 7) - 20 * (i == 6) - 4 * (i == 5) + 17 * (i == 4) + 5 * (i == 3) - 10 * (i == 2) + 2 * (i == 0);
        f[i] = f[i] / (i + 2);
    }
    printf("%lld\n", f[n].x);
    return 0;
}

:::

伪多元生成函数

这种生成函数看起来是二元的，但是可以通过一元生成函数来做（真正的多元生成函数放在[]()最后），比如我们之前讲解的P4451 [国家集训队] 整数的lqp拆分和AT_abc241_h Card Deck Score。这里，我们再讲解一道这样的题目LOJ6077 「2017 山东一轮集训 Day7」逆序对。

首先，设 f_{i, j} 表示 1 到 i 的排列中逆序对数为 j 的方案数，那么我们可以很轻松写出递推式 f_{i, j} = \displaystyle\sum_{k = 0}^{i - 1} f_{i - 1, j - k}，朴素的前缀和优化 DP 只能做到 \mathcal O(n^2)，因此我们考虑生成函数。

由于递推式是二维的，感觉用普通的一元生成函数不可做，但是我们注意到 f_{i, j} 只由上一行转移而来，此时依然可以用一元生成函数解决。

我们设第 i 行所有 f_{i, j} 的生成函数为 F_i(x)，再将递推式转化成 f_{i, j} = \displaystyle\sum_{k \geq 0} f_{i - 1, j - k} [k \leq i - 1]，如果我们设 H_i(x) 表示 [k \leq i] 这个序列的生成函数，那么递推式又可以写成 F_i(x) = F_{i - 1}(x) \times H_{i - 1}(x)，因此 F_i 可以看成 \displaystyle\sum_{j = 0}^{i - 1} H_i(x)，而 H_i(x) = (1 + x + x^2 + \dots + x^i) = \displaystyle\frac{1 - x^{i + 1}}{1 - x}，因此 F_i(x) = \displaystyle\prod_{j = 1}^i \frac{1 - x^j}{1 - x} = \frac{\prod\limits_{j = 1}^i 1 - x^j}{(1 - x)^i}。

我们发现这又是一个分子次数大于分母次数的有理函数，因此我们将分子分母分开来看，分母相当于做了 i 次前缀和，不用管它。现在我们考虑如何算出分子的第 x^k 项的系数。此时式子已经推不下去了，于是考虑能否将其转化成 DP 或者图论问题，那么一个很显然的转化就是，现在有 n 个物品，每个物品的重量分别是 1, 2, \dots, n，求装满这个背包的方案数。

参考我们在DP 学习笔记（一）：DP 基础知识，基础 DP 类型中LOJ6089 小 Y 的背包计数问题讲解的方法，我们将加物品的操作拆分成以下两种操作：

让背包中所有物品的重量增加 1；
让背包中所有物品的重量增加 1，并将一个重量为 1 的物品加入背包。

我们设 dp_{i, j} 表示放入了 i 个物品，总重量为 j 的方案数，那么根据上面的两种操作，我们可以得到 dp_{i, j} = dp_{i, j - i} + dp_{i - 1, j - i}。但是考虑到 k 最大可以达到 \displaystyle\binom n2，但是每个物品的重量却都小于等于 n，因此操作的次数必须小于等于 n。不过考虑我们生成物品序列的方法，因此如果有物品的重量超过了 n，它一定会首先到达 n + 1，我们刨去 n + 1 这个数字，剩下的方案数为 dp_{i - 1, j - n - 1}，于是我们最终的 DP 转移方程为 dp_{i, j} = dp_{i, j - i} + dp_{i - 1, j - i} - [j > n] dp_{i - 1, j - n - 1}。由于每个数字互不相同，那么 i 依然是 \mathcal O(\sqrt n) 的，因此这部分的时间复杂度为 \mathcal O(n \sqrt n) 的。我们再对其做 k 次前缀和，可以发现 \displaystyle\sum_{i = 1}^n dp_{i, k} 会乘上一个组合数，于是再花费 \mathcal O(k) 的时间复杂度就做完了这道题目。

:::info[点击查看代码]

#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N = 2e5 + 9, MOD = 1e9 + 7;
int dp[2][N], ans[N], fac[N], inv[N], n, k, len, cur, res;
int qpow(int a, int b){
    int res = 1;
    while(b){
        if(b & 1)
            res = res * a % MOD;
        a = a * a % MOD;
        b >>= 1;
    }
    return res;
}
void init(){
    fac[0] = 1;
    for(int i = 1; i < N; i++)
        fac[i] = fac[i - 1] * i % MOD; 
    inv[N - 1] = qpow(fac[N - 1], MOD - 2);
    for(int i = N - 2; i >= 0; i--)
        inv[i] = inv[i + 1] * (i + 1) % MOD; 
}
int binom(int n, int m) {
    if(n < m)
        return 0;
    return fac[n] * inv[m] % MOD * inv[n - m] % MOD;
}
signed main(){
    init();
    scanf("%lld%lld", &n, &k);
    int len = sqrt(2 * k);
    dp[0][0] = ans[0] = 1;
    for(int i = 1; i <= len; i++){
        cur ^= 1;
        for(int j = 0; j <= (i + 1) * i / 2 - 1; j++)
            dp[cur][j] = 0;
        for(int j = (i + 1) * i / 2; j <= k; j++){
            dp[cur][j] = (dp[cur ^ 1][j - i] + dp[cur][j - i]) % MOD;
            if(j >= n + 1)
                dp[cur][j] -= dp[cur ^ 1][j - n - 1];
            ans[j] = (ans[j] + (cur ? MOD - 1 : 1) * dp[cur][j] % MOD) % MOD;
        }
    }
    for(int j = 0; j <= k; j++) 
        res = (res + binom(n + j - 1, j) * ans[k - j] % MOD) % MOD;
    printf("%lld", res);
    return 0;
}

:::

参考资料

付 wj、黄 jy、张 hr 的博客；
Introductory Combinatorics(Fifth Edition) Richard A. Brualdi
Concrete Mathmatics(Second Edition) Donald L.Graham, Donald K.Knuth and Oren Patashnik
生成函数封闭形式的正确理解 Sunlight-zero
利用生成函数求斐波那契数列通项公式 zwfymqz
生成函数 RandomShuffle
部分分式展开法 ljxtt
有理展开定理与递推数列通项公式 lcw
算法学习笔记：生成函数浅讲 Re_star
题解 P6112 【直接自然溢出啥事没有加强版】 01190220csl
多项式计数杂谈 command_block