算法 - 齐次线性递推求解和优化及 BM 算法

• k 项线性递推求第 n 项的 O(k²log n) 的算法。
• 已知线性递推数列的前若干项，求递推式。
• 结合以上两者可以由线性递推数列的起始若干项得到任意一项。

矩阵快速幂

已知 x_n = a₁x_n − 1 + a₂x_n − 2 + … + a_kx_n − k 以及 x₀, x₁, …, x_k − 1，求 x_n。

$$ \begin{eqnarray} \left( \begin{array}{c} a_1 & a_2 & a_3 & \ldots & a_{k-1} & a_k \\ 1 & 0 & 0 & \ldots & 0 & 0\\ 0 & 1 & 0 & \ldots & 0 & 0\\ 0 & 0 & 1 & \ldots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \ldots & 1 & 0 \end{array} \right) \cdot \left( \begin{array}{c} x_{n-1} \\ x_{n-2} \\ x_{n-3} \\ x_{n-4} \\ \vdots \\ x_{n-k} \end{array} \right) = \left( \begin{array}{c} x_{n} \\ x_{n-1} \\ x_{n-2} \\ x_{n-3} \\ \vdots \\ x_{n-k+1} \end{array} \right) \end{eqnarray} $$

剩下的就是矩阵快速幂了。

特征多项式

这部分可略过

已知 x_n = a₁x_n − 1 + a₂x_n − 2 + … + a_kx_n − k 以及 x₀, x₁, …, x_k − 1，求 x_n。

接下来将通过某种方法将 x_n 的值表示为 α₁x_k − 1 + α₂x_k − 2 + … + α_kx₀，记为 x⃗_n = (α₁, α₂, …, α_k)。

首先显然有 $$ \begin{eqnarray} \vec x_{k}&=&(a_1,a_2,\dots,a_{k-1},a_k) \\ \vec x_{k-1}&=&(1,0,\dots,0,0) \\ \vec x_{k-2}&=&(0,1,\dots,0,0) \\ &\vdots& \\ \vec x_1&=&(0,0,\dots,1,0) \\ \vec x_0&=&(0,0,\dots,0,1) \\ \end{eqnarray} $$ 从 k + 2 项开始，利用递推式能够将其表示为前 k 个向量的线性组合 (x⃗_n + 1 = a₁x⃗_n − 1 + a₂x⃗_n − 2 + … + a_kx⃗_n − k)。如果得到了 x_n 的这种表示，那么求出它的值也就轻而易举了。接下来要做的就是加速这个过程。

正文

记 x_n 的值为 λⁿ，问题转化为利用 x_n 的特征方程 f(λ) = −λ^k + a₁λ^k − 1 + a₂λ^k − 2 + … + a_kλ⁰ = 0 以及 λⁱ(0 ≤ i < k) 的值求 λⁿ。

若将 λⁿ 表示为 g(λ) × f(λ) + r(λ)，其中 g(λ) 是关于 λ 的多项式，r(λ) 是关于 λ 的一个最高次数小于 k 的多项式。由于 f(λ) = 0，所以 λⁿ = r(λ)，而 r(λ) 的计算就很容易了。

如何求 r(λ) 呢？只需借助多项式快速幂求 λⁿ mod  f(λ) 即可。

复杂度 O(k²log n)，用 FFT 或 NTT+CRT 可以优化到 O(klog klog n)。

补充

代码量差不多，但效率更差的矩阵快速幂的方法就没有用了吗？如果是双递推（{a_n} 的递推式中有 {b_n}， {b_n} 的递推式中有 {a_n} ）之类的，可能还是要借助矩阵来做。

Berlekamp-Massey (BM) 算法求递推式

已知 x₀, x₁, …, x_m，求递推式 x_n = a₁x_n − 1 + a₂x_n − 2 + … + a_kx_n − k 在满足数列的同时，k 尽可能小。

总体思路就是依次考虑每一项，如果发现与当前递推式矛盾则进行修正，同时保证递推项数 k 尽可能小。

记 f(n) = −λⁿ + a₁λ^n − 1 + a₂λ^n − 2 + … + a_kλ^n − k(n ≥ k) （含义与上同），目标就是使 ∀n ∈ [k, m], f(n) = 0 。

当前值与递推式 f(n) 矛盾意味着 f(u) = d ≠ 0，若能找到 g(n) 使 g(u) = −1 且 ∀v < u, g(v) = 0，那么将 f(n) 修正为 f(n) + d ⋅ g(n) 即可。为了使 k 尽可能小，g 的项数要么不超过 f 的项数，要么尽可能小。

g(n) 来自于之前的 f(n) 。对于一次矛盾，得到了 f(v) = d ≠ 0 那么满足 g(u) = −1 且 ∀v < u, g(v) = 0 的一个函数为 $g(n)=-\frac{f(n)}{d \lambda^{u-v}}$。

g(n) 可以重复利用，只是随着 u 的增长，g(n) 的项数越来越多。当 f(n) 的项数小于 g(n) 时，利用 f(n) 将 g(n) 替换掉，这样可以保证 k 尽可能小。

一开始没有 g(n) 怎么办呢，将 g(n) 设成 g(n) = Cλⁿ 即可，C 可以是任意数字。

写在最后

代码见模板库。

看了很多网上的资料，几乎都是晦涩难懂。特征多项式这部分我没有涉及到矩阵，应该更好理解了。BM 算法这部分算是对 wiki 的一些补充，但记号可能不太一样。最后还是希望我写的这些对你有所帮助。