OI-wiki 的 Berlekamp–Massey 算法完全指南从最短递推式到稀疏矩阵全套黑科技【免费下载链接】OI-wiki:star2: Wiki of OI / ICPC for everyone. 某大型游戏线上攻略内含炫酷算术魔法项目地址: https://gitcode.com/GitHub_Trending/oi/OI-wikiBerlekamp–MasseyBM算法是 OI / ICPC 中求解数列最短线性递推式的经典算法给定一个长为 $n$ 的数列当最短递推式阶数为 $m$ 时它能在 $O(nm)$ 时间内求出每个前缀的最短递推式。本指南以 OI-wiki 仓库的 docs/math/berlekamp-massey.md 为主体结合仓库中 常系数齐次线性递推、随机化技巧、特征多项式 等关联文档完整覆盖 BM 算法的定义、增量构造流程、参考实现以及其在矩阵快速幂加速、最小多项式、稀疏矩阵行列式/秩、稀疏方程组求解中的全套实战用法并给出每类问题的复杂度结论。读完本文你将掌握如何仅凭一个数列的前缀即可自动恢复递推式并把这一能力推广到矩阵与线性方程组问题。算法概述与适用场景BM 算法解决的核心问题是给定一个数列求出它的最短线性递推式。这里的最短指递推阶数最小。从仓库的 docs/math/poly/linear-recurrence.md 可知常系数齐次线性递推数列C-finite / C-recursive 数列广泛出现在递推题、矩阵快速幂题与生成函数题中BM 算法的价值在于当递推系数未知时可以只根据数列的前若干项自动学习出递推式再交由 线性递推求第 k 项 的快速算法求出远处项。适用范围与前提由于 BM 算法的数值稳定性较差处理实数问题时一般很少使用。以下讨论均假定在某个质数 $p$ 的剩余系 $\mathbb{F}_p$ 下进行运算代码中的取模运算也默认以模 $p$ 为前提。复杂度若最短递推式阶数为 $m$算法复杂度为 $O(nm)$最坏情况下 $m O(n)$因此最坏复杂度为 $O(n^2)$由于每次调整递推系数时只用到上次调整时的递推系数若只需整个数列的最短递推式空间复杂度为 $O(n)$仅需保存当前递推系数与上次调整时的递推系数。递推式的定义定义一个数列 ${a_0, \dots, a_{n-1}}$ 的递推式为满足下式的序列 ${r_0, \dots, r_m}$$$ \sum_{j0}^{m} r_j a_{i-j} 0, \quad \forall i \ge m $$其中 $r_0 1$$m$ 称为该递推式的阶数。数列 ${a_i}$ 的最短递推式就是阶数最小的递推式。在实现中通常采用与上述定义稍有差别的形式定义新的递推系数 ${f_0, \dots, f_{m-1}}$满足$$ a_i \sum_{j0}^{m-1} f_j a_{i-j-1}, \quad \forall i \ge m $$容易看出 $f_i -r_{i1}$且两种表述中的阶数 $m$ 相同。代码中最终返回的递推式即第一种定义以 $1$ 开头、整体满足 $\sum_j a_{i-j} r_j 0$的形式。增量式构造算法核心流程BM 算法按顺序逐个考虑 ${a_i}$ 的每一位在递推结果出错时对递推系数进行调整。记前 $i$ 位的最短递推式为 $F_i {f_{i,j}}$。初始时显然有 $F_0 {}$空递推式。假设 $F_{i-1}$ 对数列前 $i-1$ 项均成立处理第 $i$ 项时有两种情况递推系数对 $a_i$ 也成立无需调整直接令 $F_i F_{i-1}$递推系数对 $a_i$ 不成立需要调整得到新的 $F_i$。定义差值即预测误差$$ \Delta_i a_i - \sum_{j0}^{m} f_{i-1,j}, a_{i-j-1} $$首次修正如果这是第一次修改递推系数说明 $a_i$ 是序列中第一个非零项此时直接令 $F_i$ 为 $i$ 个 $0$ 即可——这显然是一个合法的最短递推式。一般修正否则设上一次修改时已考虑的项数为 $k$。若存在序列 $G {g_0, \dots, g_{m-1}}$ 满足$$ \sum_{j0}^{m-1} g_j a_{i-j-1} 0, \quad \forall i \in [m, i) $$且 $\sum_{j0}^{m-1} g_j a_{i-j-1} \Delta_i$那么将 $F_k$ 与 $G$ 按位相加即得到一个合法的递推系数 $F_i$。构造 G 的一种可行方案$$ G {, \underbrace{0, 0, \dots, 0}{i-k-1\text{ 个}}, \ \frac{\Delta_i}{\Delta_k},\ -\frac{\Delta_i}{\Delta_k} F{k-1} ,} $$其中最后的 $-\frac{\Delta_i}{\Delta_k}F_{k-1}$ 表示把 $F_{k-1}$ 的每一项都乘以 $-\frac{\Delta_i}{\Delta_k}$ 后接在序列末尾。验证$$ \sum_{j0}^{m-1} g_j a_{i-j-1} \Delta_k \cdot \frac{\Delta_i}{\Delta_k} \Delta_i $$因此这样构造的 $G$ 是合法的将 $F_i$ 赋值为 $F_k$ 与 $G$ 逐项相加的结果即可。符号约定如果要求的是最开头定义$r_0 1$ 形式的递推式 ${r_i}$则把 ${f_j}$ 全部取相反数后在最前面插入 $r_0 1$。参考实现解析OI-wiki 在 docs/math/berlekamp-massey.md 给出了完整可运行的 C 参考实现核心逻辑如下以下p为模数power为快速幂求逆元vectorint berlekamp_massey(const vectorint a) { vectorint v, last; // v 是答案0-basedp 是模数 int k -1, delta 0; for (int i 0; i (int)a.size(); i) { int tmp 0; for (int j 0; j (int)v.size(); j) tmp (tmp (long long)a[i - j - 1] * v[j]) % p; if (a[i] tmp) continue; // 情况 1预测正确无需调整 if (k 0) { // 第一次修正第一个非零项 k i; delta (a[i] - tmp p) % p; v vectorint(i 1); continue; } vectorint u v; int val (long long)(a[i] - tmp p) * power(delta, p - 2) % p; // 即 Δ_i / Δ_k if (v.size() last.size() i - k) v.resize(last.size() i - k); (v[i - k - 1] val) % p; // 放置 Δ_i / Δ_k for (int j 0; j (int)last.size(); j) { // 拼接 -Δ_i/Δ_k * F_{k-1} v[i - k j] (v[i - k j] - (long long)val * last[j]) % p; if (v[i - k j] 0) v[i - k j] p; } // 只有当前递推式“变差”时才更新 last 与 k if ((int)u.size() - i (int)last.size() - k) { last u; k i; delta a[i] - tmp; if (delta 0) delta p; } } // 转换为 r_0 1 的表述取相反数并在最前面插入 1 for (auto x : v) x (p - x) % p; v.insert(v.begin(), 1); return v; // 满足 ∀i, Σ_{j0..m} a_{i-j} v_j 0 }实现要点解读tmp是当前递推式 $F_{i-1}$ 对 $a_i$ 的预测值a[i] tmp时直接跳过。更新last与k时使用条件u.size() - i last.size() - k这保证last始终记录上一次修正时最优的递推系数对应项数 $k$是算法能保持 $O(nm)$ 复杂度且求出最短递推式的关键。由于每次修正都是把 $F_k$ 与构造的 $G$ 逐项相加且 $G$ 的前 $i-k-1$ 项为 $0$代码中先 resize 到last.size() i - k再在i-k-1处累加val、随后逐位减去val * last[j]正是对上述数学构造的直接翻译。最后统一取模并转换为 $r_0 1$ 的标准形式返回。关于取前 $2m$ 项的实用结论朴素的 BM 算法求解的是有限项数列的最短递推式。如果待求递推式的序列有无限项但已知最短递推式的阶数上界为 $m$那么只需取出序列的前 $2m$ 项即可求出整个序列的最短递推式证明从略。这个性质是所有打表 → 跑 BM → 线性递推加速流派的根基。应用一求向量列或矩阵列的最短递推式BM 算法本身只能处理标量数列但通过随机线性投影可以推广到向量列与矩阵列向量列要求 $n$ 维向量列 $\boldsymbol{v}_i$ 的最短递推式可随机一个 $n$ 维行向量 $\boldsymbol{u}^T$计算标量序列 ${\boldsymbol{u}^T \boldsymbol{v}_i}$ 的最短递推式。由 Schwartz–Zippel 引理二者最短递推式至少有 $1 - \frac{n}{p}$ 的概率相同。矩阵列设矩阵大小为 $n \times m$随机一个 $1 \times n$ 行向量 $\boldsymbol{u}^T$ 和一个 $m \times 1$ 列向量 $\boldsymbol{v}$计算标量序列 ${\boldsymbol{u}^T A_i \boldsymbol{v}}$ 的最短递推式相同概率至少为 $1 - \frac{nm}{p}$。理论依据仓库 docs/misc/rand-technique.md 中给出了 Schwartz–Zippel 引理的完整陈述——设 $f \in F[z_1,\dots,z_k]$ 是域 $F$ 上的 $k$ 元 $d$ 次非零多项式$S$ 是 $F$ 的有限子集则至多有 $d \cdot |S|^{k-1}$ 组取值使 $f 0$。推论是若 $z_1,\dots,z_k$ 在 $S$ 中等概率独立随机选取则 $\Pr[f 0] \le d / |S|$。这里取 $|S| p$$d$ 对应向量维数 $n$或 $nm$即得上式。这解释了为何 BM 与随机化结合的错误率可以做到可接受的 $\frac{n}{p}$ 级别通常取 $p$ 为约 $10^9$ 的大质数。应用二优化矩阵快速幂设 $\boldsymbol{f}_i$ 是 $n$ 维列向量且转移满足 $\boldsymbol{f}i A \boldsymbol{f}{i-1}$则 ${\boldsymbol{f}_i}$ 是一个阶数不超过 $n$ 的线性递推向量列证明从略。做法暴力求出 $\boldsymbol{f}0, \dots, \boldsymbol{f}{2n-1}$共 $2n$ 项随机投影后调用 BM 求出最短递推式再调用 常系数齐次线性递推 的快速算法求出远处的项。复杂度对比若要求 $\boldsymbol{f}_m$$O(n^3 n\log n \log m)$第一部分为暴力递推 $2n$ 项做 $O(n)$ 次稀疏/稠密矩阵乘向量第二部分为线性递推求第 $m$ 项若 $A$ 是只有 $k$ 个非零项的稀疏矩阵复杂度降至 $O(nk n\log n \log m)$由于算法至少需要 $O(nk)$ 时间预处理压力不大的场景下也可直接用 $O(n^2 \log m)$ 的线性递推算法复杂度同样可以接受。与仓库线性递推文档的衔接仓库 docs/math/poly/linear-recurrence.md 详细介绍了两种在已知递推系数后求第 $k$ 项的算法——Fiduccia 算法构造特征多项式 $\Gamma(x) : x^d - \sum_{j0}^{d-1} c_{d-j} x^j$计算 $a_k \langle x^k \bmod \Gamma(x), A(x) \rangle$复杂度 $O(\mathsf{M}(d)\log k)$与 Bostan–Mori 算法基于 Graeffe 迭代与偶函数化分治同样 $O(\mathsf{M}(d)\log k)$。BM 算法负责猜递推这两者负责快速求值三者合起来构成完整的解题流水线。应用三求矩阵的最小多项式方阵 $A$ 的最小多项式是次数最小的满足 $f(A) 0$ 的多项式 $f$。关键观察最小多项式就是 ${A^i}$ 的最小递推式因此直接调用 BM 即可。若 $A$ 是 $n$ 阶方阵最小多项式次数显然不超过 $n$。避免 $O(n^4)$ 的优化直接对 ${A^i}$ 做矩阵乘法逐项计算会达到 $O(n^4)$。但求矩阵列最短递推式时实际求的是 ${\boldsymbol{u}^T A^i \boldsymbol{v}}$ 的最短递推式因此只需计算 $A^i \boldsymbol{v}$向量迭代而非矩阵乘法。若 $A$ 有 $k$ 个非零项复杂度为 $O(kn n^2)$。关联佐证仓库 docs/math/linear-algebra/char-poly.md 指出特征多项式可与常系数齐次线性递推联系并可与 Cayley–Hamilton 定理、多项式取模结合来加速域上求矩阵幂次的算法这与本节最小多项式 ⇒ 递推式 ⇒ BM的思路互为印证。应用四求稀疏矩阵的行列式如果能求出方阵 $A$ 的特征多项式则常数项乘上 $(-1)^n$ 就是行列式。但最小多项式不一定等于特征多项式。随机化修正把 $A$ 乘上一个随机对角阵 $B$则 $AB$ 的最小多项式至少有 $1 - \frac{2n^2 - n}{p}$ 的概率就是特征多项式。求出后再除以 $\det B$ 即可还原。设 $A$ 为 $n$ 阶方阵且有 $k$ 个非零项复杂度为 $O(kn n^2)$。这里同样用到了随机线性投影 Schwartz–Zippel 的错误率分析框架随机对角阵 $B$ 的引入把最小多项式恰好是特征多项式变成一个概率事件而 $A$ 的稀疏性保证了每次矩阵乘向量迭代 $A^i \boldsymbol{v}$只需 $O(k)$ 时间。应用五求稀疏矩阵的秩设 $A$ 是 $n \times m$ 矩阵先随机一个 $n \times n$ 对角阵 $P$ 和一个 $m \times m$ 对角阵 $Q$然后计算 $QAP A^T Q$ 的最小多项式即可。技巧不必真正调用矩阵乘法——求最小多项式时需要 $QAP A^T Q$ 乘一个向量所以依次把这几个矩阵乘到向量上即可每次对角阵乘法是 $O(n)$中间两个矩阵乘向量各为 $O(k)$。答案就是最小多项式除掉所有 $x$ 因子后剩下的次数。设 $A$ 有 $k$ 个非零项且 $n \le m$复杂度为 $O(kn n^2)$。应用六解稀疏方程组问题已知 $A\mathbf{x} \mathbf{b}$其中 $A$ 是 $n \times n$ 的满秩稀疏矩阵$\mathbf{b}$、$\mathbf{x}$ 是 $n$ 维列向量$A, \mathbf{b}$ 已知需要在低于 $n^\omega$$\omega$ 为矩阵乘法指数的复杂度内解出 $\mathbf{x}$。做法显然 $\mathbf{x} A^{-1}\mathbf{b}$。若能求出 ${A^i \mathbf{b}}$$i \ge 0$的最小递推式 ${r_0, \dots, r_{m-1}}$$m \le n$则有结论证明从略$$ A^{-1}\mathbf{b} -\frac{1}{r_{m-1}} \sum_{i0}^{m-2} A^i \mathbf{b} \cdot r_{m-2-i} $$因为 $A$ 稀疏直接按定义递推出 $\mathbf{b}, A\mathbf{b}, \dots, A^{2n-1}\mathbf{b}$ 即可共 $2n$ 个向量每个向量一次稀疏矩阵乘向量。设 $A$ 有 $k$ 个非零项复杂度为 $O(kn n^2)$。OI-wiki 在 docs/math/berlekamp-massey.md 给出了配套参考实现solve_sparse_equationsvectorint solve_sparse_equations(const vectortupleint, int, int A, const vectorint b) { int n (int)b.size(); // 0-based vectorvectorint f({b}); for (int i 1; i 2 * n; i) { // 迭代计算 b, A b, ..., A^{2n-1} b vectorint v(n); auto u f.back(); for (auto [x, y, z] : A) // 稀疏矩阵乘向量[行, 列, 值] v[x] (v[x] (long long)u[y] * z) % p; f.push_back(v); } vectorint w(n); mt19937 gen; for (auto x : w) x uniform_int_distributionint(1, p - 1)(gen); // 随机投影成标量序列 { w · A^i b } vectorint a(2 * n); for (int i 0; i 2 * n; i) for (int j 0; j n; j) a[i] (a[i] (long long)f[i][j] * w[j]) % p; auto c berlekamp_massey(a); // 求出最小递推式r_0 1 形式 int m (int)c.size(); // 按公式 A^{-1} b -1/r_{m-1} * Σ A^i b * r_{m-2-i} 累加 vectorint ans(n); for (int i 0; i m - 1; i) for (int j 0; j n; j) ans[j] (ans[j] (long long)c[m - 2 - i] * f[i][j]) % p; int inv power(p - c[m - 1], p - 2); // 除以 -r_{m-1} for (int i 0; i n; i) ans[i] (long long)ans[i] * inv % p; return ans; }实现细节说明A以稀疏三元组列表(x, y, z)存储表示 $A[x][y] z$递推 $A^i \mathbf{b}$ 时对每个非零元做一次乘加单次迭代 $O(k)$随机向量 $\boldsymbol{w}$ 的分量取自 $[1, p-1]$避开 0配合 docs/misc/rand-technique.md 中 Schwartz–Zippel 引理的错误率分析保证投影不丢失递推信息最后一步利用 BM 返回系数$r_0 1$ 形式直接组装答案避免了求逆矩阵这是整个算法低于 $n^\omega$ 的关键。例题练习OI-wiki 为该主题收录了两道配套练习原文档列于 docs/math/berlekamp-massey.md 文末LibreOJ #163. 高斯消元 2——用于训练 BM 与稀疏方程组/高斯消元结合的应用ICPC2021 台北 Gym103443E. Composition with Large Red Plane, Yellow, Black, Gray, and Blue——综合性更强的竞赛原题。建议按暴力递推前 $2m$ 项 → 随机投影 → BM 求递推式 → 线性递推求值的流水线进行演练并在实现时注意模质数 $p$ 下逆元与负数取模的正确性。总结与复杂度速查应用场景核心手法复杂度$A$ 有 $k$ 个非零项$n \le m$标量数列最短递推式增量修正 上次修正记录$O(nm)$最坏 $O(n^2)$空间 $O(n)$向量列 / 矩阵列随机投影降为标量序列正确率 $\ge 1 - \frac{nm}{p}$优化矩阵快速幂前 $2n$ 项 线性递推求值$O(nk n\log n \log m)$矩阵最小多项式最小多项式 ${A^i}$ 最小递推式$O(kn n^2)$稀疏矩阵行列式随机对角阵扰动后求最小多项式$O(kn n^2)$稀疏矩阵的秩随机对角阵 最小多项式去 $x$ 因子$O(kn n^2)$稀疏方程组求解${A^i\mathbf{b}}$ 递推式反解 $\mathbf{x}$$O(kn n^2)$整体上BM 算法的思想可以概括为三步1暴力生成足够长的前缀通常 $2m$ 项2随机投影消除向量/矩阵维度3用 BM 求出标量最短递推式。其正确性由 Schwartz–Zippel 引理保证效率由稀疏矩阵乘向量与线性递推求值算法共同支撑。若要进一步研究已知递推式后的快速求值可继续阅读仓库的 常系数齐次线性递推含 Fiduccia 与 Bostan–Mori 两种算法与 随机化技巧含 Schwartz–Zippel 引理完整证明与错误率分析。【免费下载链接】OI-wiki:star2: Wiki of OI / ICPC for everyone. 某大型游戏线上攻略内含炫酷算术魔法项目地址: https://gitcode.com/GitHub_Trending/oi/OI-wiki创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考
OI-wiki 的 Berlekamp–Massey 算法完全指南：从最短递推式到稀疏矩阵全套黑科技
OI-wiki 的 Berlekamp–Massey 算法完全指南从最短递推式到稀疏矩阵全套黑科技【免费下载链接】OI-wiki:star2: Wiki of OI / ICPC for everyone. 某大型游戏线上攻略内含炫酷算术魔法项目地址: https://gitcode.com/GitHub_Trending/oi/OI-wikiBerlekamp–MasseyBM算法是 OI / ICPC 中求解数列最短线性递推式的经典算法给定一个长为 $n$ 的数列当最短递推式阶数为 $m$ 时它能在 $O(nm)$ 时间内求出每个前缀的最短递推式。本指南以 OI-wiki 仓库的 docs/math/berlekamp-massey.md 为主体结合仓库中 常系数齐次线性递推、随机化技巧、特征多项式 等关联文档完整覆盖 BM 算法的定义、增量构造流程、参考实现以及其在矩阵快速幂加速、最小多项式、稀疏矩阵行列式/秩、稀疏方程组求解中的全套实战用法并给出每类问题的复杂度结论。读完本文你将掌握如何仅凭一个数列的前缀即可自动恢复递推式并把这一能力推广到矩阵与线性方程组问题。算法概述与适用场景BM 算法解决的核心问题是给定一个数列求出它的最短线性递推式。这里的最短指递推阶数最小。从仓库的 docs/math/poly/linear-recurrence.md 可知常系数齐次线性递推数列C-finite / C-recursive 数列广泛出现在递推题、矩阵快速幂题与生成函数题中BM 算法的价值在于当递推系数未知时可以只根据数列的前若干项自动学习出递推式再交由 线性递推求第 k 项 的快速算法求出远处项。适用范围与前提由于 BM 算法的数值稳定性较差处理实数问题时一般很少使用。以下讨论均假定在某个质数 $p$ 的剩余系 $\mathbb{F}_p$ 下进行运算代码中的取模运算也默认以模 $p$ 为前提。复杂度若最短递推式阶数为 $m$算法复杂度为 $O(nm)$最坏情况下 $m O(n)$因此最坏复杂度为 $O(n^2)$由于每次调整递推系数时只用到上次调整时的递推系数若只需整个数列的最短递推式空间复杂度为 $O(n)$仅需保存当前递推系数与上次调整时的递推系数。递推式的定义定义一个数列 ${a_0, \dots, a_{n-1}}$ 的递推式为满足下式的序列 ${r_0, \dots, r_m}$$$ \sum_{j0}^{m} r_j a_{i-j} 0, \quad \forall i \ge m $$其中 $r_0 1$$m$ 称为该递推式的阶数。数列 ${a_i}$ 的最短递推式就是阶数最小的递推式。在实现中通常采用与上述定义稍有差别的形式定义新的递推系数 ${f_0, \dots, f_{m-1}}$满足$$ a_i \sum_{j0}^{m-1} f_j a_{i-j-1}, \quad \forall i \ge m $$容易看出 $f_i -r_{i1}$且两种表述中的阶数 $m$ 相同。代码中最终返回的递推式即第一种定义以 $1$ 开头、整体满足 $\sum_j a_{i-j} r_j 0$的形式。增量式构造算法核心流程BM 算法按顺序逐个考虑 ${a_i}$ 的每一位在递推结果出错时对递推系数进行调整。记前 $i$ 位的最短递推式为 $F_i {f_{i,j}}$。初始时显然有 $F_0 {}$空递推式。假设 $F_{i-1}$ 对数列前 $i-1$ 项均成立处理第 $i$ 项时有两种情况递推系数对 $a_i$ 也成立无需调整直接令 $F_i F_{i-1}$递推系数对 $a_i$ 不成立需要调整得到新的 $F_i$。定义差值即预测误差$$ \Delta_i a_i - \sum_{j0}^{m} f_{i-1,j}, a_{i-j-1} $$首次修正如果这是第一次修改递推系数说明 $a_i$ 是序列中第一个非零项此时直接令 $F_i$ 为 $i$ 个 $0$ 即可——这显然是一个合法的最短递推式。一般修正否则设上一次修改时已考虑的项数为 $k$。若存在序列 $G {g_0, \dots, g_{m-1}}$ 满足$$ \sum_{j0}^{m-1} g_j a_{i-j-1} 0, \quad \forall i \in [m, i) $$且 $\sum_{j0}^{m-1} g_j a_{i-j-1} \Delta_i$那么将 $F_k$ 与 $G$ 按位相加即得到一个合法的递推系数 $F_i$。构造 G 的一种可行方案$$ G {, \underbrace{0, 0, \dots, 0}{i-k-1\text{ 个}}, \ \frac{\Delta_i}{\Delta_k},\ -\frac{\Delta_i}{\Delta_k} F{k-1} ,} $$其中最后的 $-\frac{\Delta_i}{\Delta_k}F_{k-1}$ 表示把 $F_{k-1}$ 的每一项都乘以 $-\frac{\Delta_i}{\Delta_k}$ 后接在序列末尾。验证$$ \sum_{j0}^{m-1} g_j a_{i-j-1} \Delta_k \cdot \frac{\Delta_i}{\Delta_k} \Delta_i $$因此这样构造的 $G$ 是合法的将 $F_i$ 赋值为 $F_k$ 与 $G$ 逐项相加的结果即可。符号约定如果要求的是最开头定义$r_0 1$ 形式的递推式 ${r_i}$则把 ${f_j}$ 全部取相反数后在最前面插入 $r_0 1$。参考实现解析OI-wiki 在 docs/math/berlekamp-massey.md 给出了完整可运行的 C 参考实现核心逻辑如下以下p为模数power为快速幂求逆元vectorint berlekamp_massey(const vectorint a) { vectorint v, last; // v 是答案0-basedp 是模数 int k -1, delta 0; for (int i 0; i (int)a.size(); i) { int tmp 0; for (int j 0; j (int)v.size(); j) tmp (tmp (long long)a[i - j - 1] * v[j]) % p; if (a[i] tmp) continue; // 情况 1预测正确无需调整 if (k 0) { // 第一次修正第一个非零项 k i; delta (a[i] - tmp p) % p; v vectorint(i 1); continue; } vectorint u v; int val (long long)(a[i] - tmp p) * power(delta, p - 2) % p; // 即 Δ_i / Δ_k if (v.size() last.size() i - k) v.resize(last.size() i - k); (v[i - k - 1] val) % p; // 放置 Δ_i / Δ_k for (int j 0; j (int)last.size(); j) { // 拼接 -Δ_i/Δ_k * F_{k-1} v[i - k j] (v[i - k j] - (long long)val * last[j]) % p; if (v[i - k j] 0) v[i - k j] p; } // 只有当前递推式“变差”时才更新 last 与 k if ((int)u.size() - i (int)last.size() - k) { last u; k i; delta a[i] - tmp; if (delta 0) delta p; } } // 转换为 r_0 1 的表述取相反数并在最前面插入 1 for (auto x : v) x (p - x) % p; v.insert(v.begin(), 1); return v; // 满足 ∀i, Σ_{j0..m} a_{i-j} v_j 0 }实现要点解读tmp是当前递推式 $F_{i-1}$ 对 $a_i$ 的预测值a[i] tmp时直接跳过。更新last与k时使用条件u.size() - i last.size() - k这保证last始终记录上一次修正时最优的递推系数对应项数 $k$是算法能保持 $O(nm)$ 复杂度且求出最短递推式的关键。由于每次修正都是把 $F_k$ 与构造的 $G$ 逐项相加且 $G$ 的前 $i-k-1$ 项为 $0$代码中先 resize 到last.size() i - k再在i-k-1处累加val、随后逐位减去val * last[j]正是对上述数学构造的直接翻译。最后统一取模并转换为 $r_0 1$ 的标准形式返回。关于取前 $2m$ 项的实用结论朴素的 BM 算法求解的是有限项数列的最短递推式。如果待求递推式的序列有无限项但已知最短递推式的阶数上界为 $m$那么只需取出序列的前 $2m$ 项即可求出整个序列的最短递推式证明从略。这个性质是所有打表 → 跑 BM → 线性递推加速流派的根基。应用一求向量列或矩阵列的最短递推式BM 算法本身只能处理标量数列但通过随机线性投影可以推广到向量列与矩阵列向量列要求 $n$ 维向量列 $\boldsymbol{v}_i$ 的最短递推式可随机一个 $n$ 维行向量 $\boldsymbol{u}^T$计算标量序列 ${\boldsymbol{u}^T \boldsymbol{v}_i}$ 的最短递推式。由 Schwartz–Zippel 引理二者最短递推式至少有 $1 - \frac{n}{p}$ 的概率相同。矩阵列设矩阵大小为 $n \times m$随机一个 $1 \times n$ 行向量 $\boldsymbol{u}^T$ 和一个 $m \times 1$ 列向量 $\boldsymbol{v}$计算标量序列 ${\boldsymbol{u}^T A_i \boldsymbol{v}}$ 的最短递推式相同概率至少为 $1 - \frac{nm}{p}$。理论依据仓库 docs/misc/rand-technique.md 中给出了 Schwartz–Zippel 引理的完整陈述——设 $f \in F[z_1,\dots,z_k]$ 是域 $F$ 上的 $k$ 元 $d$ 次非零多项式$S$ 是 $F$ 的有限子集则至多有 $d \cdot |S|^{k-1}$ 组取值使 $f 0$。推论是若 $z_1,\dots,z_k$ 在 $S$ 中等概率独立随机选取则 $\Pr[f 0] \le d / |S|$。这里取 $|S| p$$d$ 对应向量维数 $n$或 $nm$即得上式。这解释了为何 BM 与随机化结合的错误率可以做到可接受的 $\frac{n}{p}$ 级别通常取 $p$ 为约 $10^9$ 的大质数。应用二优化矩阵快速幂设 $\boldsymbol{f}_i$ 是 $n$ 维列向量且转移满足 $\boldsymbol{f}i A \boldsymbol{f}{i-1}$则 ${\boldsymbol{f}_i}$ 是一个阶数不超过 $n$ 的线性递推向量列证明从略。做法暴力求出 $\boldsymbol{f}0, \dots, \boldsymbol{f}{2n-1}$共 $2n$ 项随机投影后调用 BM 求出最短递推式再调用 常系数齐次线性递推 的快速算法求出远处的项。复杂度对比若要求 $\boldsymbol{f}_m$$O(n^3 n\log n \log m)$第一部分为暴力递推 $2n$ 项做 $O(n)$ 次稀疏/稠密矩阵乘向量第二部分为线性递推求第 $m$ 项若 $A$ 是只有 $k$ 个非零项的稀疏矩阵复杂度降至 $O(nk n\log n \log m)$由于算法至少需要 $O(nk)$ 时间预处理压力不大的场景下也可直接用 $O(n^2 \log m)$ 的线性递推算法复杂度同样可以接受。与仓库线性递推文档的衔接仓库 docs/math/poly/linear-recurrence.md 详细介绍了两种在已知递推系数后求第 $k$ 项的算法——Fiduccia 算法构造特征多项式 $\Gamma(x) : x^d - \sum_{j0}^{d-1} c_{d-j} x^j$计算 $a_k \langle x^k \bmod \Gamma(x), A(x) \rangle$复杂度 $O(\mathsf{M}(d)\log k)$与 Bostan–Mori 算法基于 Graeffe 迭代与偶函数化分治同样 $O(\mathsf{M}(d)\log k)$。BM 算法负责猜递推这两者负责快速求值三者合起来构成完整的解题流水线。应用三求矩阵的最小多项式方阵 $A$ 的最小多项式是次数最小的满足 $f(A) 0$ 的多项式 $f$。关键观察最小多项式就是 ${A^i}$ 的最小递推式因此直接调用 BM 即可。若 $A$ 是 $n$ 阶方阵最小多项式次数显然不超过 $n$。避免 $O(n^4)$ 的优化直接对 ${A^i}$ 做矩阵乘法逐项计算会达到 $O(n^4)$。但求矩阵列最短递推式时实际求的是 ${\boldsymbol{u}^T A^i \boldsymbol{v}}$ 的最短递推式因此只需计算 $A^i \boldsymbol{v}$向量迭代而非矩阵乘法。若 $A$ 有 $k$ 个非零项复杂度为 $O(kn n^2)$。关联佐证仓库 docs/math/linear-algebra/char-poly.md 指出特征多项式可与常系数齐次线性递推联系并可与 Cayley–Hamilton 定理、多项式取模结合来加速域上求矩阵幂次的算法这与本节最小多项式 ⇒ 递推式 ⇒ BM的思路互为印证。应用四求稀疏矩阵的行列式如果能求出方阵 $A$ 的特征多项式则常数项乘上 $(-1)^n$ 就是行列式。但最小多项式不一定等于特征多项式。随机化修正把 $A$ 乘上一个随机对角阵 $B$则 $AB$ 的最小多项式至少有 $1 - \frac{2n^2 - n}{p}$ 的概率就是特征多项式。求出后再除以 $\det B$ 即可还原。设 $A$ 为 $n$ 阶方阵且有 $k$ 个非零项复杂度为 $O(kn n^2)$。这里同样用到了随机线性投影 Schwartz–Zippel 的错误率分析框架随机对角阵 $B$ 的引入把最小多项式恰好是特征多项式变成一个概率事件而 $A$ 的稀疏性保证了每次矩阵乘向量迭代 $A^i \boldsymbol{v}$只需 $O(k)$ 时间。应用五求稀疏矩阵的秩设 $A$ 是 $n \times m$ 矩阵先随机一个 $n \times n$ 对角阵 $P$ 和一个 $m \times m$ 对角阵 $Q$然后计算 $QAP A^T Q$ 的最小多项式即可。技巧不必真正调用矩阵乘法——求最小多项式时需要 $QAP A^T Q$ 乘一个向量所以依次把这几个矩阵乘到向量上即可每次对角阵乘法是 $O(n)$中间两个矩阵乘向量各为 $O(k)$。答案就是最小多项式除掉所有 $x$ 因子后剩下的次数。设 $A$ 有 $k$ 个非零项且 $n \le m$复杂度为 $O(kn n^2)$。应用六解稀疏方程组问题已知 $A\mathbf{x} \mathbf{b}$其中 $A$ 是 $n \times n$ 的满秩稀疏矩阵$\mathbf{b}$、$\mathbf{x}$ 是 $n$ 维列向量$A, \mathbf{b}$ 已知需要在低于 $n^\omega$$\omega$ 为矩阵乘法指数的复杂度内解出 $\mathbf{x}$。做法显然 $\mathbf{x} A^{-1}\mathbf{b}$。若能求出 ${A^i \mathbf{b}}$$i \ge 0$的最小递推式 ${r_0, \dots, r_{m-1}}$$m \le n$则有结论证明从略$$ A^{-1}\mathbf{b} -\frac{1}{r_{m-1}} \sum_{i0}^{m-2} A^i \mathbf{b} \cdot r_{m-2-i} $$因为 $A$ 稀疏直接按定义递推出 $\mathbf{b}, A\mathbf{b}, \dots, A^{2n-1}\mathbf{b}$ 即可共 $2n$ 个向量每个向量一次稀疏矩阵乘向量。设 $A$ 有 $k$ 个非零项复杂度为 $O(kn n^2)$。OI-wiki 在 docs/math/berlekamp-massey.md 给出了配套参考实现solve_sparse_equationsvectorint solve_sparse_equations(const vectortupleint, int, int A, const vectorint b) { int n (int)b.size(); // 0-based vectorvectorint f({b}); for (int i 1; i 2 * n; i) { // 迭代计算 b, A b, ..., A^{2n-1} b vectorint v(n); auto u f.back(); for (auto [x, y, z] : A) // 稀疏矩阵乘向量[行, 列, 值] v[x] (v[x] (long long)u[y] * z) % p; f.push_back(v); } vectorint w(n); mt19937 gen; for (auto x : w) x uniform_int_distributionint(1, p - 1)(gen); // 随机投影成标量序列 { w · A^i b } vectorint a(2 * n); for (int i 0; i 2 * n; i) for (int j 0; j n; j) a[i] (a[i] (long long)f[i][j] * w[j]) % p; auto c berlekamp_massey(a); // 求出最小递推式r_0 1 形式 int m (int)c.size(); // 按公式 A^{-1} b -1/r_{m-1} * Σ A^i b * r_{m-2-i} 累加 vectorint ans(n); for (int i 0; i m - 1; i) for (int j 0; j n; j) ans[j] (ans[j] (long long)c[m - 2 - i] * f[i][j]) % p; int inv power(p - c[m - 1], p - 2); // 除以 -r_{m-1} for (int i 0; i n; i) ans[i] (long long)ans[i] * inv % p; return ans; }实现细节说明A以稀疏三元组列表(x, y, z)存储表示 $A[x][y] z$递推 $A^i \mathbf{b}$ 时对每个非零元做一次乘加单次迭代 $O(k)$随机向量 $\boldsymbol{w}$ 的分量取自 $[1, p-1]$避开 0配合 docs/misc/rand-technique.md 中 Schwartz–Zippel 引理的错误率分析保证投影不丢失递推信息最后一步利用 BM 返回系数$r_0 1$ 形式直接组装答案避免了求逆矩阵这是整个算法低于 $n^\omega$ 的关键。例题练习OI-wiki 为该主题收录了两道配套练习原文档列于 docs/math/berlekamp-massey.md 文末LibreOJ #163. 高斯消元 2——用于训练 BM 与稀疏方程组/高斯消元结合的应用ICPC2021 台北 Gym103443E. Composition with Large Red Plane, Yellow, Black, Gray, and Blue——综合性更强的竞赛原题。建议按暴力递推前 $2m$ 项 → 随机投影 → BM 求递推式 → 线性递推求值的流水线进行演练并在实现时注意模质数 $p$ 下逆元与负数取模的正确性。总结与复杂度速查应用场景核心手法复杂度$A$ 有 $k$ 个非零项$n \le m$标量数列最短递推式增量修正 上次修正记录$O(nm)$最坏 $O(n^2)$空间 $O(n)$向量列 / 矩阵列随机投影降为标量序列正确率 $\ge 1 - \frac{nm}{p}$优化矩阵快速幂前 $2n$ 项 线性递推求值$O(nk n\log n \log m)$矩阵最小多项式最小多项式 ${A^i}$ 最小递推式$O(kn n^2)$稀疏矩阵行列式随机对角阵扰动后求最小多项式$O(kn n^2)$稀疏矩阵的秩随机对角阵 最小多项式去 $x$ 因子$O(kn n^2)$稀疏方程组求解${A^i\mathbf{b}}$ 递推式反解 $\mathbf{x}$$O(kn n^2)$整体上BM 算法的思想可以概括为三步1暴力生成足够长的前缀通常 $2m$ 项2随机投影消除向量/矩阵维度3用 BM 求出标量最短递推式。其正确性由 Schwartz–Zippel 引理保证效率由稀疏矩阵乘向量与线性递推求值算法共同支撑。若要进一步研究已知递推式后的快速求值可继续阅读仓库的 常系数齐次线性递推含 Fiduccia 与 Bostan–Mori 两种算法与 随机化技巧含 Schwartz–Zippel 引理完整证明与错误率分析。【免费下载链接】OI-wiki:star2: Wiki of OI / ICPC for everyone. 某大型游戏线上攻略内含炫酷算术魔法项目地址: https://gitcode.com/GitHub_Trending/oi/OI-wiki创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考