线性递推求值问题

常系数线性递推关系的定义

对于满足
$$a_n=\sum _{i=1}^k{c}_i{a}_{n-i}+f(n)$$
一类的递推关系（其中 $c,k$ 为常数），我们称之为常系数线性递推关系。

只要给定 $a_0,a_1,\dots ,a_{k-1}$ 以及 $c_1,c_2,\dots ,c_k$ 以及 $f(n)$ ，我们就能求出 $a$ 的任意一项。

这个问题的一般形式是不能快速求解的，比如令 $f(n)=\ln n$ 。但有一些特殊形式我们可以快速求解。我们先从齐次常系数线性递推关系来研究。

齐次常系数线性递推关系的求值

齐次常系数线性递推关系即 $f(n)=0$ 的情况，此时 $a_n={\sum }_{i=1}^k{c}_i{a}_{n-i}$ 。

这看起来像线性变换，因此考虑从线性变换的角度解决这个问题。

矩阵优化

我们构造初始向量

$$\vec{x}=\left[\begin{array}{c}{a}_{k-1}\\ a_{k-2}\\ ⋮\\ a_1\\ a_0\end{array}\right]$$

以及转移矩阵

$$A=\left[\begin{array}{cccccc}{c}_1& c_2& \cdots & c_{k-2}& c_{k-1}& c_k\\ 1& 0& \cdots & 0& 0& 0\\ 0& 1& \cdots & 0& 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& \cdots & 1& 0& 0\\ 0& 0& \cdots & 0& 1& 0\end{array}\right]$$

于是

$$A^n\cdot \vec{x}=\left[\begin{array}{c}{a}_{n+k-1}\\ a_{n+k-2}\\ ⋮\\ a_{n+1}\\ a_n\end{array}\right]$$

这样我们就可以用矩阵快速幂来求 $a_n$ 了，复杂度 $O(k^3\log n)$ 。

但这个复杂度依然不够优秀，考虑从线性代数的角度切入。

特征多项式优化

事实上，我们可以把 $A^n$ 表示为 $A^0,A^1,\cdots ,A^{k-1}$ 的线性组合，并且我们有：

$$A^n=\sum _{i=1}^k{c}_i{A}^{n-i}$$

为什么呢？
我们可以从 $Caylay-Camilton$ 定理¹来考虑：

对于矩阵 $A$ 的特征多项式 $f(\lambda )=\det \left|\lambda I-A\right|$ ，我们有 $f(A)=0$ 。

而对于上述转移矩阵，我们有：

$$f(\lambda )={\lambda }^k-\sum _{i=1}^k{c}_i{\lambda }^{n-i}$$

将 $A$ 代入上式即可。

当然我们也可以从其他角度考虑。

$$A^n\cdot \vec{x}=\left[\begin{array}{c}{a}_{n+k-1}\\ a_{n+k-2}\\ ⋮\\ a_{n+1}\\ a_n\end{array}\right]$$

$$\sum _{i=1}^k{c}_i{A}^{n-i}\cdot \vec{x}=\left[\begin{array}{c}{\sum }_{i=1}^k{c}_i{a}_{n-i+k-1}\\ \\ {\sum }_{i=1}^k{c}_i{a}_{n-i+k-2}\\ \\ ⋮\\ \\ {\sum }_{i=1}^k{c}_i{a}_{n-i+1}\\ \\ {\sum }_{i=1}^k{c}_i{a}_{n-i}\\ \end{array}\right]=\left[\begin{array}{c}{a}_{n+k-1}\\ \\ a_{n+k-2}\\ \\ ⋮\\ \\ a_{n+1}\\ \\ a_n\\ \end{array}\right]=A^n\cdot \vec{x}$$

也即：

$$(A^n-\sum _{i=1}^k{c}_i{A}^{n-i})\cdot \vec{x}=0$$

由于任意向量 $\vec{x}$ 都满足这个关系，因此我们就得到了 $A^n-{\sum }_{i=1}^k{c}_i{A}^{n-i}=0$ 。

取 $n=k$ 就有 $f(A)=0$ 。

无论如何，有了上述关系，我们只要求出 ${\lambda }^n\mathrm{m}\mathrm{o}\mathrm{d}f(\lambda )$ 即可。

设 $g(\lambda )={\lambda }^n\mathrm{m}\mathrm{o}\mathrm{d}f(\lambda )$ ，则 $A^n=f(A)P(A)+g(A)=0\cdot P(A)+g(A)=g(A)$ 。

再设

$$g(\lambda )=\sum _{i=0}^{k-1}{h}_i{\lambda }^i$$

则

$$A^n\cdot \vec{x}=g(A)\cdot \vec{x}=\sum _{i=0}^{k-1}{h}_i{A}^i\cdot \vec{x}=\left[\begin{array}{c}{\sum }_{i=0}^{k-1}{h}_i{a}_{k-1+i}\\ \\ {\sum }_{i=0}^{k-1}{h}_i{a}_{k-2+i}\\ \\ ⋮\\ \\ {\sum }_{i=0}^{k-1}{h}_i{a}_{1+i}\\ \\ {\sum }_{i=0}^{k-1}{h}_i{a}_i\\ \end{array}\right]$$

这实际上是一个卷积的形式。我们只要把 $h$ 和 $a$ 的前 $l+k-1$ 项做多项式乘法就能求出连续的 $l$ 项，如果只要求出单项就只需计算 ${\sum }_{i=0}^{k-1}{h}_i{a}_i$。

求 $g(\lambda )$ 我们可以用倍增+多项式取模的方式。

如果使用 $FFT$ 那么总复杂度就是 $O(k\log k\log n)$ ，暴力取模就是 $O(k^2\log n)$ 。

以及还有一个问题：怎么求出 $a$ 的前 $l$ 项？暴力求是 $O(kl)$ 的。

由于这是一个序列问题，要求的项数不多，我们考虑用生成函数法解决。
设

$$f(x)=\sum _{i=0}^{\infty }{a}_i{x}^i,C(x)=\sum _{i=1}^k{c}_i{x}^i$$

由于将 $f,C$ 卷积在前 $k$ 项是错误的，因此我们构造

$$g(x)=f(x)-f(x)C(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^k)$$

于是我们有 $f(x)=f(x)C(x)+g(x)$ 。

因此 $f(x)=\frac{g(x)}{1-C(x)}$ 。

值得一提的是 $1-C(x)$ 和转移矩阵特征多项式的关系。以及逆推这个过程我们会发现：若序列 $a$ 的生成函数满足 $f(x)=\frac{g(x)}{1-C(x)}$ 且 $\mathrm{deg}g<\mathrm{deg}C=k$ 且 $C$ 的常数项为 $0$ （实际上不为 $1$ 即可），我们就可以把该序列与一个 $k$ 阶线性递推方程相对应，这一点对非齐次线性递推也同样适用。

于是我们就能在 $O(k\log k\log n)$ 的时间内解决齐次常系数线性递推的求值问题了。

一类特殊非齐次常系数线性递推关系的求值

注：以下部分的一些多项式算法，在我的这一篇博文中均有给出解法。

我们现在已经掌握齐次线性递推关系的求值算法了。

现在来讨论 $f(n)$ 为多项式函数的情况。

我们尝试用化归法转化为齐次线性线性递推关系求解。

我们知道：
$$a_n=\sum _{i=1}^k{c}_i{a}_{n-i}+f(n)$$
$$a_{n-1}=\sum _{i=1}^k{c}_i{a}_{n-1-i}+f(n-1)$$
$$a_n-a_{n-1}=\sum _{i=1}^k{c}_i{a}_{n-i}-\sum _{i=1}^k{c}_i{a}_{n-1-i}+f(n)-f(n-1)=\sum _{i=1}^k{c}_i{a}_{n-i}-\sum _{i=1}^k{c}_i{a}_{n-1-i}+f_1(n)$$

显然 $\mathrm{deg}{f}_1=\mathrm{deg}f-1$ ，设 $f$ 的次数为 $m$ ，我们只需要对其进行 $m+1$ 次对减就能转化为齐次关系。
现在问题在于如何求出新的递推关系以及 $a_k,\cdots ,a_{m+k}$ 。

求 $a_k,\cdots ,a_{m+k}$ 的值我们可以仿照上面的生成函数法：
设

$$F(x)=\sum _{i=0}^{\infty }{a}_i{x}^i,C(x)=\sum _{i=1}^k{c}_i{x}^i,G(x)=F(x)-F(x)C(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^k),H(x)=\sum _{i=k}^{\infty }f(i)x^i$$

则 $F(x)=F(x)C(x)+G(x)+H(x),F(x)=\frac{G(x)+H(x)}{1-C(x)}$

不过这样我们需要使用一次多点求值算法，复杂度 $O(m{\log }^{2}m)$ 。

至于新的递推关系，由于对减的特点，设原递推关系的转移矩阵特征多项式为 $f(\lambda )$ ，新的递推关系的转移矩阵的特征多项式就是 $f(\lambda )(x-1{)}^{m+1}$ ，这样就得到了新的递推关系式。

于是求解该问题的总复杂度为 $O(m{\log }^{2}m+(m+k)\log (m+k)\log n)$ 。

对该问题的进一步优化

实际上除了求 $a_0,a_1,\cdots ,a_{m+k}$ 的值以外，我们仅需要求得 $x^n\mathrm{m}\mathrm{o}\mathrm{d}f(x)(x-1{)}^{m+1}$ 。

若 $f(x)$ 与 $(x-1{)}^{m+1}$ 互质，我们可以求出 $x^n\mathrm{m}\mathrm{o}\mathrm{d}f(x)$ ，再求出 $x^n\mathrm{m}\mathrm{o}\mathrm{d}(x-1{)}^{m+1}$ ，最后使用中国剩余定理合并。不互质的话我们只需提取出 $f(x)$ 中的 $(x-1)$ 因式即可。

对于 $x^n\mathrm{m}\mathrm{o}\mathrm{d}(x-1{)}^{m+1}$ ，如果我们使用倍增取模，那么复杂度就是 $O(m\log m\log n+k\log k\log n)$ ，并没有真正的优化复杂度。

但显然由于所模函数的特殊性，我们可以快速取模。我们只要求出 $f(x+1)\mathrm{m}\mathrm{o}\mathrm{d}{x}^{m+1}$ ，再代入 $(x-1)$ 就能求出 $f(x)\mathrm{m}\mathrm{o}\mathrm{d}(x-1{)}^{m+1}$ 。

由于被取模函数的特殊性甚至可以不需要用到多项式算法，而直接使用组合数学的方法 $O(m)$ 算出。当然这不是重点。

根据中国剩余定理我们需要求解这样的方程：
$$f(x)P(x)+(x-1{)}^{m+1}Q(x)=1$$

我们可以直接套用 $exgcd$ 算法，复杂度是 $O((m+k)\log (m+k))$ ，但常数较大。

一个较好的做法是转化为 $f(x+1)P(x+1)+x^{m+1}Q(x+1)=1$ 。这样就是多项式求逆了。

求 $F(x+q)$ 是一个经典的问题，可见于上文提到的另一篇博文。

求出 $P(x)$ 之后，设 $A(x)=x^n\mathrm{m}\mathrm{o}\mathrm{d}f(x),B(x)=x^n\mathrm{m}\mathrm{o}\mathrm{d}(x-1{)}^{m+1}$ ，那么就有：

$$x^n\equiv A(x)(x-1{)}^{m+1}Q(x)+B(x)f(x)P(x)\equiv A(x)(1-f(x)P(x))+B(x)f(x)P(x)(\mathrm{m}\mathrm{o}\mathrm{d}f(x)(x-1{)}^{m+1})$$

这样求解的总复杂度是 $O(m{\log }^{2}m+k\log k\log n)$ 。

更一般的形式

对于

$$f(n)=\sum _{i=1}^t{F}_i(n)q_i^n$$

其中 $F_i(n)$ 为多项式函数， $q_i$ 两两不相等。

我们依然可以仿照上面的方式进行求值。

每次设新的数列 $b_n=\frac{a_n}{q_i^n}$ 并逐次差分即可。

设 $F_i$ 的次数为 $m_i$ ，总差分次数为 $M=\sum _{i=1}^t(m_i+1)$ ，由于每次消去一个多项式时其余多项式的次数并不会发生变化，因此我们最终会将其化为齐次的 $k+M$ 阶方程。

考虑如何求解新递推关系的特征方程。

设我们上一次消去的函数为第 $l$ 项即 $F_l(n)q_l^n$ ，则新的序列 $b_n$ 满足 $b_n=\frac{a_n}{q_l^n}$ 。

则 $b_n$ 的递推方程的特征多项式：
$$(x^k-\sum _{i=1}^k\frac{c_i}{q_l^i}{x}^{n-i})\prod _{i=1}^l(x-\frac{q_i}{q_l}{)}^{m_i+1}$$

而此时的函数形式为

$$f^{\prime }(n)=\sum _{i=l+1}^t{F}_i^{\prime }(n)(\frac{q_i}{q_l}{)}^n$$

且各多项式函数的次数不变。

这并不难推导出来。当我们需要消去 $F_{l+1}^{\prime }$ 时，新的序列为 $b^{\prime }(n)=\frac{b_n}{\frac{q_{l+1}^n}{q_l^n}}=b_n\frac{q_l^n}{q^{l+1}}=\frac{a_n}{q_{l+1}^n}$ ，则特征方程中每一个因式从高到低各项分别乘上 $(\frac{q_l}{q_{l+1}}{)}^j$ ，且所有 $F_i^{\prime }(n)(\frac{q_i}{q_l}{)}^n$ 乘上 $\frac{q_l^n}{q^{l+1}}$ ，变为 $F_i^{\prime }(n)(\frac{q_i}{q_{l+1}}{)}^n$ 。

此时 $F_{l+1}^{\prime }$ 的指数为 $1$ ，其余函数指数不为 $1$ ，差分 $m_{l+1}+1$ 次会导致特征多项式乘上 $(x-1{)}^{m_{l+1}+1}$ 并消去 $F_{l+1}^{\prime }$，且其余函数的多项式次数不变。

归纳证明即可。

通过这个结论，我们可以推导出 $a_n$ 最终化为齐次递推方程后方程的特征多项式为

$$(x^k-\sum _{i=1}^k{c}_i{x}^{n-i})\prod _{i=1}^t(x-q_i{)}^{m_i+1}$$

我们依然可以用原来的方法计算 $x^n\mathrm{m}\mathrm{o}\mathrm{d}(x-q_i{)}^{m_i+1}$ ，并用类似解同余方程的方法以及 $exgcd$ 算法合并，注意要使用分治合并的方法，最坏情况下复杂度为 $O(M\log M\log t)$ 。

同时我们要求出 $a$ 的 $O(k+M)$ 项。使用生成函数法我们还是要求出每个函数的点值，单次求点值的复杂度为 O ( min ⁡ { m i log ⁡ 2 m i + M log ⁡ m i , m i M } ) O( \min\{m_i\log^2 m_i+M\log m_i,m_iM\}) O(min{mi​log2mi​+Mlogmi​,mi​M}) ，（注意多项式求值也可以当做线性递推问题求连续点的点值），并且我们需要求 $t$ 次。当 $t$ 较小时复杂度是 $O(tM\log M+M{\log }^{2}M)$ 的， $t$ 较大时使用暴力求值法更优秀，复杂度为 $O(M^2)$ 。

总复杂度 O ( min ⁡ { t M log ⁡ M + M log ⁡ 2 M , M 2 } + k log ⁡ k log ⁡ n ) O(\min\{tM\log M+M\log^2 M,M^2\}+k\log k \log n) O(min{tMlogM+Mlog2M,M2}+klogklogn) 。

对一般形式更进一步的优化

在上述算法中，最劣情况下复杂度含有一个 $M^2$ ，显然是非常不优秀的。

考虑什么地方出现了问题。

显然瓶颈在于如何求 $a_k,a_{k+1},\cdots ,a_{k+M}$ ，而当 $t$ 很大时也就是底数不同的指数函数过多时暴力求出点值用以表示序列的生成函数是不优秀的。

类似于序列 $1^n+2^n+\cdots +m^n$ 的生成函数为 $\sum _{i=1}^m\frac{1}{1-ix}$ ，显然可以使用分治 $FFT$ 在 $O(m{\log }^{2}m)$ 时间内求出，但是若是使用上面的点值法表示生成函数复杂度就会达到 $m^2$ 。

现在我们考虑如何快速求出序列的通项形如 $a_n=F(n)q^n$ 的生成函数。

首先我们依然可以设 $b_n=\frac{a_n}{q^n}$ 将该问题化归为求解 $F(n)$ 的生成函数。

多项式生成函数问题在上文提到的博文中有给出解法，这里不再赘述。

对于线性递推式：

$$a_n=\sum _{i=1}^k{c}_i{a}_{n-i}+f(n)$$

$$f(n)=\sum _{i=1}^t{F}_i(n)q_i^n$$

我们只需设 $$H(x)=\sum _{i=k}^{\infty }f(i)x^i=x^k\sum _{i=1}^t\frac{G_i(x)}{(1-q_{i}x{)}^{m_i+1}}$$
即可。

只需 $O(\sum _{i=1}^t{m}_i{\log }^2{m}_i)$ 求出每一个 $G_i$ 。

然后使用与上面使用生成函数法时一致的做法，即：

$$F(x)=\frac{G(x)+H(x)}{1-C(x)}$$

就可以求出 $a_k,a_{k+1},\cdots ,a_{M+k}$ 的值了。

类似的，依然有总复杂度为 $O(M{\log }^{2}M+k\log k\log n)$ 。

---

1. $Caylay-Camilton$ 定理可以从线性代数的相关资料以及网络上找到证明。 ↩︎