【学习笔记】BM算法和线性递推

$\text{BM}$ 算法

一、概述

全称是 $\text{Berlekamp-Massey}$ 算法，用于求出一个序列 $⟨a_0,a_1,a_2,\dots ,a_n⟩$ 的最短 线性递推式。

线性递归式：一个有限非空数列 $⟨r_0,r_1,\dots ,r_m⟩$ 使得任意 $p⩾m$ 有 ${\sum }_{i=0}^m{a}_{p-i}{r}_i=0$，则 $⟨r⟩$ 就是 $⟨a⟩$ 的线性递归式。从生成函数的角度来看，令 $⟨a⟩$ 对应 $f(x)$，而 $⟨r⟩$ 对应 $A(x)$，则存在次数不超过 $(m-1)$ 的多项式 $R(x)$ 使得 $A(x)f(x)=R(x)$ 。

线性递推式：满足 $r_0=1$ 的线性递归式。或者说 $[x^0]A(x)=1$ 。然并无甚区别。

记号约定：记 $|f(x)|$ 为形式幂级数 $f(x)$ 的次数，即最大的 $k$ 使得 $[x^k]f(x)\ne 0$ 。

二、过程

看上去可怕，过程好像也并不复杂。假设我们求出了前 $n$ 项的线性递推式
$$A_n(x)f_n(x)\equiv R_n(x)(\mathrm{mod}{x}^n)$$

现在试试能否据此找到前 $(n+1)$ 项的线性递推式。如果 $A_n(x)f_n(x)\equiv R_n(x)(\mathrm{mod}{x}^{n+1})$，那么 $A_n(x)$ 就是。所以我们只讨论
$$A_n(x)f_n(x)\equiv R_n(x)+d{x}^n(\mathrm{mod}{x}^{n+1})\text{\hspace{1em}(1)}$$

假如上次在 $p(p<n)$ 处，也有一次 “失配” 的情况，即
$$A_p(x)f(x)\equiv R_p(x)+c{x}^p(\mathrm{mod}{x}^{p+1})\text{\hspace{1em}(2)}$$

并且 满足 $|A_{p+1}(x)|>|A_p(x)|$ 。那么 $(2)$ 式两边同时乘 $x^{n-p}d{c}^{-1}$ 得
$$x^{n-p}d{c}^{-1}{A}_p(x)f(x)\equiv x^{n-p}d{c}^{-1}{R}_p(x)+d{x}^n(\mathrm{mod}{x}^{n+1})$$

然后与 $(1)$ 式相减得
$$\left[A_n(x)-x^{n-p}d{c}^{-1}{A}_p(x)\right]f(x)\equiv R_n(x)-x^{n-p}d{c}^{-1}{R}_p(x)(\mathrm{mod}{x}^{n+1})$$

所以我们找到了新的解是
$$\{\begin{array}{l}{A}_{n+1}(x)=A_n(x)-\frac{d}{c}{x}^{n-p}{A}_p(x)\\ R_{n+1}(x)=R_n(x)-\frac{d}{c}{x}^{n-p}{R}_p(x)\end{array}$$

由 $|R_n(x)|<|A_n(x)|\wedge |R_p(x)|<|A_p(x)|$ 可知 $|R_{n+1}(x)|<|A_{n+1}(x)|$，而且 $[x^0]A_{n+1}(x)=1$，符合定义。时间复杂度 $\mathcal{O}(n^2)$，其中 $n$ 是最短线性递推式长度；我们后面也会说到，$n$ 跟序列长度是同级的。

三、最短性

可行的线性递推式一抓一大把；唯独 $\text{BM}$ 抓住了最短的那个！下文记 $A_n(x)$ 为前 $n$ 项的最短线性递推式。

首先 $|A_n(x)|⩾|A_{n-1}(x)|$ 是显然的。再结合下面的引理，就可以说明最短性了。

Lemma：若 $|A_n(x)|>|A_{n-1}(x)|$，则最短线性递推式长度 $|A_n(x)|⩾n-|A_{n-1}(x)|$ 。

插嘴：这个 $\mathrm{Lemma}$ 使用 “线性递推式” 的定义其实很好证明；我已写了生成函数的证明，不愿换掉罢了。

证明：反证法。若 $|A_n(x)|⩽n-1-|A_{n-1}(x)|$，由 $[x^0]A_n(x)=1$ 可知 $[1-A_n(x)]$ 是 $x$ 的倍数，同乘于 $A_{n-1}(x)$ 的定义式方程
$$A_{n-1}(x)f(x)\equiv R_{n-1}(x)(\mathrm{mod}{x}^n)$$

就可以让模数乘 $x$ 了。得
$$[1-A_n(x)]A_{n-1}(x)f(x)\equiv [1-A_n(x)]R_{n-1}(x)(\mathrm{mod}{x}^n)$$

对等式左侧进行处理，把 $f(x)$ 乘进括号内，用 $A_n(x)f(x)\equiv R_n(x)$ 作代换得
$$[f(x)-R_n(x)]A_{n-1}(x)\equiv R_{n-1}(x)-A_n(x)R_{n-1}(x)$$

再把括号拆开，然后移项。
$$A_{n-1}(x)f(x)\equiv R_{n-1}(x)+A_{n-1}(x)R_n(x)-A_n(x)R_{n-1}(x)(\mathrm{mod}{x}^n)\text{\hspace{1em}(3)}$$

我们敏锐地注意到后两项是类似叉积的形式，联想到克莱默法则。或许它是某个方程组的解？正是如此。
$$\begin{gathered} \{\begin{array}{l}{A}_{n-1}(x)f(x)\equiv R_{n-1}(x)\\ A_n(x)f(x)\equiv R_n(x)\end{array}(\mathrm{mod}{x}^{n-1}) \\ ⟹A_{n-1}(x)R_n(x)-A_n(x)R_{n-1}(x)\equiv 0(\mathrm{mod}{x}^{n-1}) \end{gathered}$$

联系我们的前提条件
$$|A_{n-1}(x)R_n(x)-A_n(x)R_{n-1}(x)|<|A_n(x)|+|A_{n-1}(x)|⩽n-1$$

可知此处取模毫无意义，二者就是相等。代入 $(3)$ 式即知 $A_{n-1}(x),R_{n-1}(x)$ 是取模 $x^n$ 意义下的合法解，与前提条件 $|A_n(x)|>|A_{n-1}(x)|$ 矛盾。证毕。$■$

证明 $\mathrm{Lemma}$ 之后，我们就有 $|A_{n+1}(x)|⩾\max \{|A_n(x)|,n+1-|A_n(x)|\}$ 。用归纳法证明这个下界恰好被取到，则是比较简单的事情了：

“失配” 时，递推式告诉我们 $|A_{n+1}(x)|=\max \{|A_n(x)|,n-p+|A_p(x)|\}$ 。只需看后者。由于 $p$ 是上一次变长的位置，即 $i\in (p,n]$ 的 $|A_i(x)|$ 都相同，再由归纳假设，得 $n-p+|A_p(x)|=n-p+p+1-|A_{p+1}(x)|=n+1-|A_n(x)|$，证毕。

四、无穷项序列

无穷项序列也可能存在线性递推式。若我们已知其线性递推式满足 $|A(x)|⩽m$，那么只需要求出前 $2m$ 项的线性递推式。

证明蕴含在 $\mathrm{Lemma}$ 中：若以后还想更改线性递推式，就至少是 $n+1-|A(x)|⩾2m+1-m=m+1$ 长度的线性递推式了，不满足我们的已知信息。

注意：此方法不能逆用；仅在已知 $m$ 的范围时可行。

五、初值

$\text{BM}$ 算法的赋初值，真的不是那么平凡的。比如 $A_1(x)$ 是什么？或者从 $A_0(x)$ 开始么？

值得注意的是，$|A_n(x)|$ 不一定小于 $n$ 。虽然这对于线性递推式的成立是毫无帮助的，但是它可以让 $|R_n(x)|<|A_n(x)|$ 恒成立。说白了还是妥协之计。

那么，我们其实可以从 $A_0(x)$ 开始赋值。就设定为 $\{\begin{array}{l}{A}_0(x)=1\\ R_0(x)=0\end{array}$ 就好。其实 $|R_0(x)|=-\infty$，这是为了使得 $|f(x)g(x)|=|f(x)|+|g(x)|$ 恒成立；所以此时 $|A_0(x)|>|R_0(x)|$ 成立，当然也满足 $[x^0]A_0(x)=1$ 的定义。

这样赋值以后，可能会让 $A_n(x)$ 及以前都保留该设定：在 $a_0=a_1=\cdots =a_{n-1}=0$ 的时候就会发生。然后如果 $A_{n+1}(x)$ 失配了呢？即 $A_n(x)f(x)\equiv R_n(x)+a_n{x}^n(\mathrm{mod}{x}^{n+1})$，那么根据我们的调整方式，我们会得到
$$\{\begin{array}{l}{A}_{n+1}(x)=A_n(x)-x^{n-p}\frac{a_n}{c}{A}_p(x)\\ R_{n+1}(x)=R_n(x)-x^{n-p}\frac{a_n}{c}{R}_p(x)\end{array}$$

而事实上唯一可行的解是
$$\{\begin{array}{ll}{A}_{n+1}(x)=1+\lambda x^{n+1}& (\lambda \ne 0)\\ R_{n+1}(x)=a_n{x}^n& \end{array}$$

毕竟我们需要满足 $|A_{n+1}(x)|=n+1-|A_n(x)|=n+1$ 嘛。对比一下系数可知，
$$\{\begin{array}{ll}p=0& \\ R_p(x)=-c& (c\ne 0)\\ A_p(x)=x& \end{array}$$

就是我们需要的初始转移点。其中 $c$ 是任取的。由于在代码实现时，$A_p,R_p$ 与 $A_n,R_n$ 是分别存储的，此处 $A_0,R_0$ 矛盾的定义在代码实现中是可以被接受的。

六、示例代码

int a[MAXN], r[MAXN];
int BM(const int v[],int n){
	static int ap[MAXN], rp[MAXN], p, invc;
	static int an[MAXN], rn[MAXN]; // new
	int lenp = 2, lenn = 1; p = 0;
	rp[0] = ap[1] = 1, rp[1] = ap[0] = 0;
	a[0] = 1, r[0] = 0; invc = MOD-1;
	for(int i=0,d; i!=n; ++i){
		for(int j=d=0; j!=lenn; ++j)
			d = int((d+llong(a[j])*v[i-j])%MOD);
		if(d == 0) continue; // keep it
		const llong f = MOD-llong(d)*invc%MOD;
		for(int j=i-p; j!=lenp+i-p; ++j){
			an[j] = int((a[j]+f*ap[j-i+p])%MOD);
			rn[j] = int((r[j]+f*rp[j-i+p])%MOD);
		}
		if(lenp+i-p > lenn){ // update p
			memcpy(ap,a,lenn<<2), memcpy(rp,r,lenn<<2);
			memcpy(a+i-p,an,lenp<<2), memcpy(r+i-p,rn,lenp<<2);
			swap(lenp,lenn), lenn += i-p, p = i;
			invc = int(qkpow(d,MOD-2));
		}
		else memcpy(a+i-p,an,lenp<<2), memcpy(r+i-p,rn,lenp<<2);
	}
	return lenn;
}

比 $\text{BM}$ 更快

可见于 $\text{EntropyIncreaser}$ 的博客。因为我还没搞懂，因此讲不了。

常系数齐次线性递推

零、前言

在学废这个特征多项式的做法之后，我发现了 四、更好的做法 。推荐大家只阅读那个，而不需要在这个渣做法上浪费时间了。

Comment. 该方法在求解远处的连续值时倒是确实无可替代。

一、定义

形如
$$x_n=\sum _{i=1}^k{a}_i{x}_{n-i}$$

也就是上面的所谓线性递推式。

• 常系数：系数 $a_i$ 不随 $n$ 的改变而改变。
• 齐次线性：$x^q$ 等于前面几项的 $x^q$ 的线性和。

那么，定长（一定由前 $k$ 项确定下一项）的含义在哪里呢 😂

二、特征方程

想必大家都看过数竞的小蓝本。其中有一种求递推通项公式的方法：特征方程。在这里，我也不便引入数竞的那一套；我将不加证明地说，下面的 $f(x)$ 是特征多项式，$f(x)=0$ 就是特征方程。
$$f(x)=x^k-\sum _{i=1}^k{a}_i{x}^{k-i}$$

特征方程 $f(x)=0$ 的含义就是，以 $x$ 为公比的等比数列满足该递推式。

巧妙的是，如果我们用 矩阵 改写这个方程式，得到
$$f(A)=A^k-\sum _{i=1}^k{a}_i{A}^{k-i}=0$$

这里 $A$ 是一个矩阵，$0$ 代表全零矩阵。你会发现，这个方程的一个解就是转移矩阵。

简单证明一下。如果您对这一套有深入的了解，晓得 $f(x)$ 真的是特征多项式，又晓得特征多项式是零化多项式，就不用管下面的证明了。下面的证明比较初等，但是正确。

不妨设状态向量 $v_i=(x_i,x_{i-1},x_{i-2},\dots ,x_{i-k+1})$ 为列向量。转移矩阵的第一行是 $a_1,a_2,a_3,\dots ,a_k$，以后的每一行都只有 $A_{i(i-1)}=1$ 。

此时考虑矩阵的幂的意义，就是在有向图上走特定的步数，两个点之间有多少条路径。考虑从点 $1$ 出发，如果你第一步从 $1$ 走到 $r$，那么第二步只能是 $r$ 到 $(r-1)$，第三步到 $(r-2)$，直到第 $r$ 步回到 $1$ 。这一路径的方案数是 $a_r$，由第一步带来的系数。

那么，我可以大胆地说，这就像 势能。那么枚举第一步（如果起点不是 $1$，那么第一步是走回 $1$ 之后的第一步）获得的势能，一定满足
$$A^k=\sum _{i=1}^k{a}_i{A}^{k-i}$$

所以 $f(A)=0$ 得证。这个证明感觉不太严谨，仔细想想，又挺正确的……

接下来，根据多项式取模，可以得到 $g(x),r(x)$，满足
$$x^{n-k}=g(x)f(x)+r(x)$$

此处可以使用快速幂，就是一个带取模的快速幂，计算
$$r(x)\equiv x^{n-k}(\mathrm{mod}f(x))$$

复杂度是 $\mathcal{O}(k\log k\log n)$ 的，还不算太差。令 $x=A$ 则
$$A^{n-k}=g(A)f(A)+r(A)$$

直接改变了定义域，看上去是胆大妄为，实际上是很自然的。因为关于 $A$ 矩阵的整式也构成环。而我们已知 $f(A)=0$，所以我们得到重要等式
$$A^{n-k}=r(A)$$

不妨设 $r(x)$ 的系数为 { b i }    ( 0 ⩽ i < k ) \{b_i\}\;(0\leqslant i<k) {bi​}(0⩽i<k)，我们想要计算的是 $A^{n-k}\overrightarrow{{\mathbf{v}}_k}$ 即
$$r(A)\overrightarrow{{\mathbf{v}}_k}=\sum _{i=0}^{k-1}{b}_i{A}^i\overrightarrow{{\mathbf{v}}_k}$$

注意到 $A^i\overrightarrow{{\mathbf{v}}_k}=\overrightarrow{{\mathbf{v}}_{k+i}}=(x_{k+i},x_{k+i-1},x_{k+i-2},\dots ,x_{i+1})$，所以我们实际上要求的就是
$$r(A)\overrightarrow{{\mathbf{v}}_k}=\sum _{i=0}^{k-1}{b}_i\overrightarrow{{\mathbf{v}}_{k+i}}$$

而且我们只需要向量的首项，所以右侧实际上是 ${\sum }_{i=0}^{k-1}{b}_i{x}_{k+i}$ 。计算这些项，可以用卷积，所以复杂度 $\mathcal{O}(k\log k)$ 。

总时间复杂度是 $\mathcal{O}(k\log k\log n)$ 的，吊打 $\mathcal{O}(k^3\log n)$ 矩阵快速幂优化啊！

三、总结

我回看自己的博客，发现我不知道代码咋打 😢 赶紧总结一下。

对于递推 $x_n={\sum }_{i=1}^k{a}_i{x}_{n-i}$ 而言，

• 第一步：令 $f(x)=x^k-{\sum }_{i=0}^{k-1}{a}_{k-i}{x}^i$ 。
• 第二步：令 $g(x)=x^{n-k}\mathrm{mod}f(x)$ ，求出 $g(x)={\sum }_{i=0}^{k-1}{b}_i{x}^i$ 。
• 第三步：$x_n={\sum }_{i=0}^{k-1}{b}_i{x}_{k+i}$ 。

四、更好的方法

毒瘤读毒瘤论文：以毒攻毒。$\mathsf{EI}$ 搞了个新做法，据 $\text{Tiw-Air-OAO}$ 称，好写并且常数小（即便无下文中的常数优化），值得信赖！

在 $\text{BM}$ 算法里讲过，线性递推等价于 $A(x)f(x)=R(x)$，所以可以解出 $f(x)=\frac{R(x)}{A(x)}$ 。更一般地，我们就是求
$$[x^{\tau }]\frac{A(x)}{B(x)}([x^0]B(x)\ne 0)$$

上下同时乘以 $B(-x)$，然后按照指数的奇偶性分类。即，存在形式幂级数 $F(x),G(x),D(x)$ 使得
$$\{\begin{array}{l}A(x)B(-x)=xF(x^2)+G(x^2)\\ B(x)B(-x)=D(x^2)\end{array}$$

毕竟 $B(x)B(-x)$ 是偶函数，无奇数次项。于是有
$$\frac{A(x)}{B(x)}=x\cdot \frac{F(x^2)}{D(x^2)}+\frac{G(x^2)}{D(x^2)}$$

若 $2\nmid n$ 则递归到求 $[x^{\frac{n-1}{2}}]\frac{F(x)}{D(x)}$，否则递归到 $[x^{\frac{n}{2}}]\frac{G(x)}{D(x)}$ 。若 $A(x),B(x)$ 次数不超过 $n$，则 $F(x),D(x)$ 次数也不超过 $n$，复杂度 $\mathcal{O}(n\log n\log \tau )$ 没问题。

常数优化一

计算 $B(x)B(-x)$ 时，无需两次正变换，因为 $B(-x)$ 的对应位置的点值 $B(-{\omega }_n^k)=B({\omega }_n^{k+\frac{n}{2}})$，即 $B(x)$ 的点值的第 $(k+\frac{n}{2})$ 项。

常数优化二

求出 $A(x)B(-x)$ 的点值表达式之后，我们只需要提取 $F(x)$ 或 $G(x)$ 其中之一。是否存在更快的逆变换方法？

记 $\delta (x)=A(x)B(-x)$，那么只需求出 $F(x^2)=\frac{\delta (x)+\delta (-x)}{2}$ 或 $G(x^2)=\frac{\delta (x)-\delta (-x)}{2}$ 。看来得到 $F(x)$ 的点值 $F({\omega }_n^k)=\frac{\delta ({\omega }_{2n}^k)+\delta ({\omega }_{2n}^{k+n})}{2}$ 毫不费力，因为 $\delta (x)$ 的点值是已知的。对 $G(x^2)$ 同理。然后逆变换得到 $F(x)$ 或 $G(x)$ 就行了。逆变换的长度减小到了 $\frac{1}{2}$，可喜可贺！

当然，在 $D(x^2)=B(x)B(-x)$ 中也可以使用这个技巧。

常数优化三

在 常数优化二 中，我们成功算出了 $F(x)$ 或 $G(x)$ 的点值表达式。可否压根儿就不做那逆变换，直接递归？同样地，我们也不做 $Q(x)$ 的逆变换，只保留点值？

仍然是求解 $[x^n]\frac{A(x)}{B(x)}$，但是 $A(x),B(x)$ 被表示为 ${\omega }_n$ 的 $\text{DFT}$ 点值结果。为了求得 $A(x)B(-x)$，需要得到 ${\omega }_{2n}$ 的 $\text{DFT}$ 点值结果。似乎这并不需要从头再算，因为 $A({\omega }_{2n}^{2k})=A({\omega }_n^k)$ 可直接得到。

想办法在 $n$ 的长度内解决 ${\omega }_{2n}^{2k+1}$ 。则必然要递归到 ${\omega }_n^k$ 的点值。提出一个 ${\omega }_{2n}$ 就行了！所以我们得到
$$A({\omega }_{2n}^{2k+1})=\sum _i({\omega }_{2n}^i{a}_i)\cdot {\omega }_n^{ik}$$

其中 $a_i=[x^i]A(x)$ 。只需要先做长为 $n$ 的 $\text{IDFT}$，将系数 $a_i$ 乘上 ${\omega }_{2n}^i$，然后再做一次 $\text{DFT}$ 就结束了。

审视一下目前得到的方法：每次递归只需要做长度为 $n$ 的 $\text{IDFT}$ 和 $\text{DFT}$ 各两次（对于 $A(x)$ 和 $B(x)$ 都需要点值表达式的长度翻倍），然后是两次长度为 $2n$ 的点值乘法和加法，再没别的！多么清爽的做法！

注记

其实多项式取模和多项式求逆（多项式除法）的本质相同，因为 $\mathcal{O}(n\log n)$ 的多项式取模不就是先求逆么？