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

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

矩阵快速幂

已知 \(x_n=a_1x_{n-1}+a_2x_{n-2}+\dots+a_kx_{n-k}\) 以及 \(x_0,x_1,\dots,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_1x_{n-1}+a_2x_{n-2}+\dots+a_kx_{n-k}\) 以及 \(x_0,x_1,\dots,x_{k-1}\),求 \(x_n\)

接下来将通过某种方法将 \(x_{n}\) 的值表示为 \(\alpha_1x_{k-1}+\alpha_2x_{k-2}+\dots+\alpha_kx_0\),记为 \(\vec x_n=(\alpha_1, \alpha_2,\dots,\alpha_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\) 个向量的线性组合 (\(\vec x_{n+1}=a_1\vec x_{n-1}+a_2\vec x_{n-2}+\dots+a_k\vec x_{n-k}\))。如果得到了 \(x_n\) 的这种表示,那么求出它的值也就轻而易举了。接下来要做的就是加速这个过程。

正文

\(x_n\) 的值为 \(\lambda^n\),问题转化为利用 \({x_n}\) 的特征方程 \(f(\lambda)=-\lambda^k+a_1\lambda^{k-1}+a_2\lambda^{k-2}+\dots+a_k\lambda^0=0\) 以及 \(\lambda^i(0\le i\lt k)\) 的值求 \(\lambda^n\)

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

如何求 \(r(\lambda)\) 呢?只需借助多项式快速幂求 \(\lambda^n \mod f(\lambda)\) 即可。

复杂度 \(O(k^2\log n)\),用 FFT 或 NTT+CRT 可以优化到 \(O(k\log k \log n)\)

补充

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

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

已知 \(x_0,x_1,\dots,x_{m}\),求递推式 \(x_n=a_1x_{n-1}+a_2x_{n-2}+\dots+a_kx_{n-k}\) 在满足数列的同时,\(k\) 尽可能小。

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

\(f(n)=-\lambda^n+a_1\lambda^{n-1}+a_2\lambda^{n-2}+\dots+a_k\lambda^{n-k}(n\ge k)\) (含义与上同),目标就是使 \(\forall n\in[k,m],f(n) =0\)

当前值与递推式 \(f(n)\) 矛盾意味着 \(f(u)=d \neq 0\),若能找到 \(g(n)\) 使 \(g(u)=-1\)\(\forall v\lt u,g(v)=0\),那么将 \(f(n)\) 修正为 \(f(n) + d\cdot g(n)\) 即可。为了使 \(k\) 尽可能小,\(g\) 的项数要么不超过 \(f\) 的项数,要么尽可能小。

\(g(n)\) 来自于之前的 \(f(n)\) 。对于一次矛盾,得到了 \(f(v)=d\neq 0\) 那么满足 \(g(u)=-1\)\(\forall v\lt 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\lambda^n\) 即可,\(C\) 可以是任意数字。

写在最后

代码见模板库。

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