Berlekamp-Massey

感觉很有用。

算法

Berlekamp-Massey 可以求出一个数列的最短线性递推式。对于序列 \(a_1,a_2,\dots,a_n\)，若序列 \(r_0,r_1,\dots,r_k\) 满足 \(r_1\neq 0\)，且 \(\sum_{j=0}^k a_{i-j}r_j=0\) 对 \(i>k\) 都成立，则称 \(r\) 是 \(a\) 的一个线性递推式，且其长度为 \(k\)。注意最短递推式的长度可能会等于原序列的长度。

过程如下：

从左到右考虑每个前缀，维护当前的最短递推式 \(r\)。设加入 \(i\) 之后递推式求出的值为 \(\Delta_i\)。

将原序列和当前的递推式看作两个生成函数 \(A(x),R_i(x)\)，则余项 \(S_i(x)=A(x)R(x)\bmod x^{k+1}\) 应是次数不超过 \(k-1\) 的多项式。

• 若递推式仍满足，即 \(\Delta_i=0\)，则继续；

• 否则调整当前的递推式。设上一次调整是在 \(i=k\) 的时候，且误差为 \(\Delta_k\neq 0\)。（第一次调整特判掉）
  那么我们有：

  \[R_{i-1}(x)A(x)\equiv S_i(x)+\Delta_i x^i\pmod{x^{i+1}}\\ R_{k-1}(x)A(x)\equiv S_k(x)+\Delta_k x^k\pmod{x^{k+1}}\\ \]

  其中 \(S_i(x),S_k(x)\) 的次数分别不超过 \(i-1,k-1\)。
  把第二个式子乘上 \(\dfrac{\Delta_i}{\Delta_k}x^{i-k}\) 并作差，得到 \(A(x)(R_{i-1}(x)-x^{i-k}\dfrac{\Delta_i}{\Delta_k}R_{k-1}(x))\equiv S_i(x)-S_k(x)\pmod{x^{i+1}}\)，而右边的次数不超过 \(i-1\)，这样就得到了一个新的递推式。

为了让长度最短，调整的时候 \(k\) 应当取之前求出的递推式中 \(i-size\) 最大的一项。证明不会，就这吧。

设最短递推式长度为 \(r\)，复杂度 \(O(nr)\)，最坏为 \(O(n^2)\)。代码实现很简单（~700B），且常数不算大：

vector<int> BM(int *a, int n)
{
	vector<int> seq = {1}, lst;
	int k = -1, delta;
	for (int i = 1; i <= n; i++)
	{
		int sum = 0;
		for (int j = 0; j < seq.size(); j++) (sum += a[i - j] * seq[j]) %= MOD;
		if (!sum) continue;
		
		if (k == -1)
		{
			lst = seq;
			k = i, delta = sum;
			while (seq.size() < i + 1) seq.push_back(0);
			continue;
		}
		
		vector<int> tmp = seq;
		int c = sum * qpow(delta, MOD - 2) % MOD;
		seq.resize(max((int)(lst.size() + i - k), (int)seq.size()));
		for (int j = 0; j < lst.size(); j++)
			(seq[j + i - k] -= lst[j] * c) %= MOD;
		if ((int)tmp.size() - i < (int)lst.size() - k) lst = tmp, k = i, delta = sum;
	}
	return seq;
}

应用

求出线性递推式之后，如果要求第 \(k\) 项，通常可以直接用比较好写的 \(O(n^2\log k)\) 快速线性递推，毕竟 BM 自身就已经 \(O(n^2)\) 了。

找规律

有结论：如果一个无限长度的序列具有长度为 \(r\) 的最短递推式，则取前 \(2r\) 项一定可以求出来它。证明依然不会。

求向量序列的最短递推式

随机左乘/右乘一个向量后求递推式。可以证明求出的递推式有至少 \(1-n/p\) 的概率是对的。

求矩阵序列的最短递推式

随机左乘再右乘一个向量后求递推式。可以证明求出的递推式有至少 \(1-(n+m)/p\) 的概率是对的。

---

下面的算法通常用于稀疏矩阵。设矩阵中有 \(e\) 个非零位置。

优化稀疏矩阵快速幂

求 \(A^k v\)。可以证明 \(A^iv\) 是一个阶数不超过 \(n\) 的线性递推。

求出 \(A^0v,A^1v,A^2v,\dots,A^{2n}v\) 的最短递推式，然后跑一个线性递推即可。

如果有 \(q\) 次询问，复杂度是 \(O(ne+qn^2\log k)\) 的，用多项式做线性递推可以做到单次查询 \(O(n\log n\log k)\)。

有趣的是即使 \(e=O(n^2)\)，这个算法相较于朴素的矩阵快速幂（\(O(n^3\log k+qn^2\log k)\)）也少了一个 \(\log\)。

求稀疏矩阵的最小多项式

相当于求 \(A^0,A^1,A^2,\dots\) 的最短递推式。复杂度为 \(O(ne)\)。

求稀疏矩阵的行列式

特征多项式的常数项乘上 \((-1)^n\) 就是行列式。

随机乘上一个对角阵 \(B\) 之后，可以证明有至少 \(1-O(n^2/p)\) 的概率最小多项式等于特征多项式，所以做完了，记得最后除掉 \(\det B\)。

复杂度仍然是 \(O(ne)\)。

求稀疏矩阵的秩

按上述方法求 \(AB\) 的特征多项式。有多少个特征值非 \(0\)，\(rank\) 就是多少，因此考察最低次项即可。复杂度 \(O(ne)\)。

注意这两个 case 都是利用一些性质，只使用乘对角阵后的特征多项式获取答案；原矩阵的特征多项式可能不太好还原回去。

解稀疏方程组（矩阵求逆）

对于满秩矩阵 \(A\) 求解方程 \(Ax=b\)。

显然 \(x=A^{-1}b\)，而 \(A^{-1}\) 可以通过最短递推式得到：\(\sum_{i=0}^n r_iA^{n-i}=0\implies A^0=-\dfrac{1}{r_i}\sum_{i=0}^{n-1} r_iA^{n-i}\implies A^{-1}=-\dfrac{1}{r_i}\sum_{i=0}^{n-1} r_iA^{n-i-1}\)。

---

这样对于稀疏图的生成树计数（Laplacian）、最大匹配大小（Tutte 矩阵）、随机游走等等可以通过线性代数手法计算的问题，通常都能优化到 \(O(nm)=O(n^2)\)，相当于去掉了一个 \(n\)。