很多人小时候应该都玩过这样一个游戏给你一串数字让你猜下一个是什么。比如 1, 1, 2, 3, 5, 8, 13……稍微有点经验的人都能看出这是斐波那契数列下一项是 21。但如果把这串数字拉长到几千项或者背后隐藏的规律根本不是“后一项等于前两项之和”这么简单单纯靠肉眼猜就行不通了。这时候就需要一套系统的数学工具最短线性递推式求解以及跟它紧密相关的有理函数重建。这两个问题看起来是纯代数领域的理论活儿但实际应用非常广。信号处理里要从观测序列恢复系统模型密码学里要分析线性反馈移位寄存器LFSR的结构编码理论里 BCH/Reed-Solomon 码的译码要解关键方程数值计算里的 Padé 逼近本质上也在做类似的事。我当年是在算法竞赛里第一次被这个东西震撼到一道题要求根据数列前若干项求通项公式标准做法就是先跑一遍 Berlekamp-Massey 算法把最短递推式求出来再由递推式还原生成函数。整个过程优雅得不像话。这篇文章我就想把这条线完整梳理一遍最短线性递推式怎么求Berlekamp-Massey 算法的原理和代码怎么落地以及递推式如何进一步还原成有理函数。同时我会把我实际踩过的一些坑、验证结果的小技巧也一起写出来尽量让读者看完之后能直接照着写一套能跑的代码。1. 先把两个问题的数学模型一起摊开1.1 线性递推式到底在描述什么所谓线性递推式就是序列里每一项都能由前面固定数量的若干项线性组合出来。形式化写是这样的s[n] c1·s[n-1] c2·s[n-2] ... ck·s[n-k]这个 k 就叫递推阶数整数系数 c1 到 ck 是递推系数。斐波那契数列就是 k2c11c21 的特例。当然系数不一定是整数可以是实数、复数、有理数甚至模某个素数 p 意义下的整数。满足某个 k 阶线性递推的序列全体构成了一个 k 维线性空间。这个观点特别重要因为“最短”这个概念其实就是在问给定序列的前 N 项能不能找到一个最小的 k以及对应的系数使得从第 k1 项开始每一项都严格符合这个递推关系。注意这里说的是前 N 项内部的约束不是要求递推对未来所有项都成立——虽然在实际场景中一旦某个递推对整个序列成立它就对未来所有项成立。在工程背景里这个模型最常见的化身就是线性反馈移位寄存器。一个 k 级的 LFSR 输出的序列本质上就是一个满足 k 阶线性递推的序列。流密码分析里有个经典攻击手段就是拿到足够长的密钥流之后用 Berlekamp-Massey 算法反推出 LFSR 的结构参数。这就是最短线性递推式最直接的一个应用场景。1.2 有理函数重建是什么有理函数重建是反过来的视角。给定一个形式幂级数的前若干项我们希望找到一个有理函数F(x) P(x) / Q(x)其中 P(x) 和 Q(x) 都是多项式并且 Q(0) ≠ 0一般还约定 Q(0)1使得 F(x) 的泰勒展开前若干项恰好和给定序列吻合。也就是说如果我们把序列看成某个生成函数的系数那么这个生成函数很可能就是一个有理函数我们现在要根据截断的系数把它找回来。这里有个非常漂亮的结论一个形式幂级数是有理函数的充分必要条件就是它的系数序列满足某个有限阶线性递推。必要性很好理解如果 F P/Q那么把等号两边同时乘以 Q比较系数就能得到一个线性递推关系递推阶数不超过 Q 的次数。反过来如果系数序列满足 k 阶递推那就能构造一个分母次数为 k 级的有理函数。所以“最短线性递推式求解”和“有理函数重建”根本就是一枚硬币的两面。前者求的是序列背后的最小递推结构后者求的是生成函数的最小分母表示。某个意义下Padé 逼近就是在干这件事给定截断幂级数找一对次数尽量低的多项式 (P, Q) 来逼近它。当我们要求误差是 O(x^N) 级别时Padé 逼近给出的正是有理函数重建的一个合法解。1.3 为什么要把两者打通把这两个问题放在一起讨论不是因为它们长得像而是因为实际操作中可以串成一条完整流水线。第一步拿到序列的前 N 项用 Berlekamp-Massey 算法求出最短线性递推得到递推阶数 k 和系数 c1..ck第二步把递推系数直接翻译成分母多项式 Q(x)1-c1·x-c2·x^2-...-ck·x^k第三步利用序列前 k 项和 Q(x) 的卷积反推出分子多项式 P(x)。三步走完一个完整的有理函数就出来了。可能有人会问为什么不直接用 Padé 逼近一步到位因为 Berlekamp-Massey 算法的数值稳定性通常更好而且它对“最短”的把握非常精确。Padé 逼近虽然也能做但涉及构造 Hankel 矩阵、解线性方程组繁琐不说在模素数环境下还没有 BM 算法那么干净利落。所以很多高质量代码库内部都是先用 BM 求出递推再用递推恢复有理函数而不是直接去解线性方程组。2. Berlekamp-Massey 算法从“猜规律”到严格求解2.1 算法直觉像修水管一样增量修正我第一次接触 BM 算法的时候第一反应是这玩意儿怎么能叫“算法”它看起来更像是“猜”出来的。后来才明白它其实是把“猜”的过程系统化了每一步都在用当前最优的递推去预测下一位预测错了就修正修正完之后再继续往下走。整个过程有点像修水管哪漏水就补哪里补完继续加压测试再漏再补。用一个具体的场景来说。假设现在已经处理完序列的前 i-1 项手里有一个长度为 L 的递推多项式 C(x)它能完美预测前 i-1 项。现在轮到第 i 项我们用 C(x) 去预测得到一个预测值。如果预测值和真实值一致那就继续往下走如果不一致就说明当前递推在局部“漏了”需要更新 C(x)。关键的技巧在于更新 C(x) 的时候不是随便修修补补而是利用上一次失配时保存的修正信息构造一个新的修正项加到旧多项式上。这个修正项要保证两个性质第一修完之后前 i-1 项依然满足第二第 i 项也能被修正到正确值。只要算法能始终维护这两条那么递推长度 L 就会一直保持最小。为什么不是暴力枚举所有可能的递推因为暴力枚举是 O(N^3) 甚至更高而 BM 算法巧妙地利用了修正项的结构每一步只做一个多项式加法总复杂度只有 O(N^2)。N 是序列长度。这在竞赛和工程中都是可以接受的。2.2 关键公式失配时怎么更新形式化描述 BM 算法需要引入连接多项式。记当前递推多项式为C(x) 1 c1·x c2·x^2 ... cL·x^L注意这里常数项是 1所以“递推关系”其实是 C(x) 与序列卷积的前若干项为 0。也就是对任意 n ≥ L有s[n] c1·s[n-1] ... cL·s[n-L] 0在算法第 i 步我们用当前的 C(x) 计算一个误差项delta s[i] c1·s[i-1] ... cL·s[i-L]如果 delta 0万事大吉如果 delta ≠ 0就需要更新。更新公式长这样C_new(x) C(x) - (delta / delta_last) · x^(i - i_last) · C_last(x)其中 C_last(x) 是上一次失配时保存的旧连接多项式delta_last 是那次失配时的误差值i_last 是那次失配的索引。这个公式的直观含义是把旧修正多项式搬到当前位置再乘以一个比例系数用来抵消当前的 delta。为什么能保证修正后前面依然正确因为旧修正多项式在 i_last 之前的卷积结果全是 0乘以 x 的幂次之后它对更早位置没有任何影响而在当前位置上它正好贡献出 delta 的相反数把误差归零。这套逻辑几乎是几何级的巧妙。长度 L 的更新也有规则只有在新长度超过当前长度的时候才更新 L并且更新 L 的瞬间要保存 C_last 和 delta_last。具体来说是当 L i 1 - L 时L 变成 i 1 - L。这个条件保证算法求出的长度在所有可能递推中是最小的标准的证明用的是反证法这里不展开但结论可以直接信任BM 算法给出的递推长度一定是最短的。2.3 可直接运行的 Python 实现说了这么多理论直接上代码。我用 Python 写一个最经典的版本支持模素数 p 运算这样比特权重的浮点实现更稳也更容易嵌进各种算法题或密码学代码里。def berlekamp_massey(s, mod): # s: 序列前 N 项list 或数组 # mod: 素数模数所有运算在这个模下进行 # 返回 (L, C)L 是递推阶数C 是连接多项式系数C[0]1 C [1] # 当前连接多项式 B [1] # 上次失配时的连接多项式 L 0 # 当前递推长度 m 1 # 距上次失配的步数 b 1 # 上次失配时的误差 delta_last for i in range(len(s)): # 计算当前误差 delta d s[i] % mod for j in range(1, L 1): d (d C[j] * s[i - j]) % mod if d 0: m 1 continue # 需要修正 if 2 * L i: # 保存当前 C 作为新的 B T C[:] # 计算修正比例 coef d * pow(b, mod - 2, mod) % mod # 扩展 C 到足够长度 if len(C) len(B) m: C [0] * (len(B) m - len(C)) # C_new C - coef * x^m * B for j in range(len(B)): C[j m] (C[j m] - coef * B[j]) % mod L i 1 - L B T b d m 1 else: # 长度不变只修正系数 coef d * pow(b, mod - 2, mod) % mod if len(C) len(B) m: C [0] * (len(B) m - len(C)) for j in range(len(B)): C[j m] (C[j m] - coef * B[j]) % mod m 1 # 裁掉可能多余的高位系数 C C[:L 1] return L, C这段代码有几个细节值得注意。一是模逆元用了费马小定理 pow(b, mod-2, mod)要求 mod 是素数二是两个分支里都做了 C[jm] 的更新区别只在于是否更新 L、B、b、m三是我最后把 C 裁到 L1 长度保证输出干净。如果序列本身是整数序列且不需要模运算可以把所有运算放在有理数域里做用 Fraction 类型代替整数。但一般不建议这么做因为分数运算会让时间复杂度膨胀得很厉害而且容易溢出大整数的范围。实际工程里要么用模大素数比如 998244353 这种 NTT 素数兜底要么用浮点 BM 算法后面我在避坑部分会细说。2.4 手推一个例子斐波那契数列光看代码可能还是有点抽象我们手动推一遍最简单的斐波那契序列。设 s[1, 1, 2, 3, 5, 8]取模数 1000000007。初始C[1], B[1], L0, m1, b1。i0s[0]1。计算 d 1。2L0 i0所以进入长度更新分支。coef1C 变成 [1,-1]模意义下是 [1, mod-1]L 01-01B[1]b1m1。现在 C 表示的递推是 s[n]s[n-1] 吗其实是常数数列的递推。i1s[1]1。d s[1] C[1]*s[0] 1 (-1)*1 0。没有失配m2。i2s[2]2。d s[2] C[1]s[1] 2 - 1 1非零。2L2 i2所以又进入长度更新分支。coef d * b^{-1} 1。当前 C[1, -1]B[1]m2。更新后 C [1, -1] - x^2[1] [1, -1, -1]也就是 1 - x - x^2。新长度 L 21-2 1等等这里按公式 L i 1 - L 3 - 1 2。对是 L2。B[1,-1]b1m1。i3s[3]3。用 C[1,-1,-1] 预测d s[3] (-1)*s[2] (-1)*s[1] 3 - 2 - 1 0正确。i4s[4]5。d 5 - 3 - 2 0正确。i5s[5]8。d 8 - 5 - 3 0正确。最终得到 L2C[1, mod-1, mod-1]对应递推 s[n] s[n-1] s[n-2]完美还原斐波那契。整个过程和预期完全一致。从这个例子能看出BM 算法一开始可能会给出一个“错的”常数递推但后面会逐步修正成正确结果而且一旦修正到位后续就一路零误差说明递推找对了。3. 从递推系数到有理函数分子分母怎么拼出来3.1 分母直接抄分子靠卷积反推拿到最短递推多项式 C(x) 1 c1·x c2·x^2 ... ck·x^k 之后有理函数的分母就直接对应Q(x) C(x) 1 c1·x c2·x^2 ... ck·x^k注意这个符号选择和前面递推方向一致递推关系写成 s[n] a1·s[n-1] ... ak·s[n-k] 的话那么 C(x) 1 - a1·x - ... - ak·x^k。所以 BM 输出的连接多项式和生成函数的分母是同一个东西只是中间那些系数的符号要按上面的规则对应好。分子怎么求设生成函数 F(x) P(x)/Q(x)那么 P(x) F(x)·Q(x)。如果只看截断幂级数记 A(x) s[0] s[1]·x ... s[N-1]·x^(N-1)那么P(x) A(x) · Q(x) mod x^{k1}这里只要取到 x^k 这一项就够了因为分子多项式最多 k 次。原因很简单P 的次数 Q 的次数 k而 Q 的常数项是 1所以 F 的幂级数展开是唯一的。我们拿 A(x) 和 Q 卷积然后把高于 k 次的项全部丢掉剩下的就是 P。实际操作里有个更简单的等价做法用序列前 k 项和 Q 的系数逐个卷积。设 P(x) p0 p1·x ... pk·x^k那么它有递推关系p[n] sum_{j0}^{n} s[j] · C[n-j], for n 0..k直接算就行复杂度 O(k^2)。3.2 最少需要多少项2k 原则这里必须强调一个边界条件要唯一确定一个 k 阶递推序列至少需要给多少项答案是需要大约 2k 项。为什么因为 k 阶递推有 k 个未知系数而递推关系对序列的第 k1 到第 N 项都要成立这给出了 N-k 个约束方程。只有当 N-k ≥ k也就是 N ≥ 2k 时方程个数才不少于未知数个数递推系数才可能被唯一确定。BM 算法的输出在样本不足的情况下可能会不稳定多个不同的低阶递推都能适配当前样本算法返回哪一个取决于内部细节。这个道理也可以从生成函数角度理解。一个分母 k 次的有理函数分子也是 k-1 次多项式所以一共有 2k 个自由度。只看前 N 项只有在 N ≥ 2k 的时候才能把这 2k 个自由度的信息全部包含进来。所以做实验的时候我会先跑一遍 BM看看递推阶数 L 大概是多少。如果 L 已经接近 N/2那说明给的数据量有点危险还得再多采集一些样本才能下结论。这个检查在时序分析和系统辨识里非常重要。3.3 端到端小实验从有理函数出发再绕回来为了验证整套流程我建议做一个闭环实验先随便构造一个有理函数展开成序列再用 BM 重建。这里构造一个稍复杂的例子F(x) (x 1) / (1 - 2x - 3x^2)展开前 8 项试试。按长除法F(x) 1 3x 9x^2 27x^3 81x^4 243x^5 729x^6 2187x^7 ...这个序列看似简单前 8 项就是 3 的幂但其实它被有理函数约束着递推阶数是 2。跑一遍 BM应该得到 C(x) 1 - 2x - 3x^2也就是 L2。然后再利用前 2 项算分子p0 s0 1p1 s1 C1·s0 3 (-2)·1 1。所以 P(x) 1 x和初始构造完全一致。我写了一段简单的 Python 验证代码思路是这样的def seq_from_rational(P, Q, n): # 通过递推关系生成前 n 项 k max(len(P), len(Q)) - 1 P P [0] * (k 1 - len(P)) Q Q [0] * (k 1 - len(Q)) s [] for i in range(n): val P[i] if i len(P) else 0 for j in range(1, i 1): val - Q[j] * s[i - j] s.append(val) return s def rational_from_seq(s, mod): L, C berlekamp_massey(s, mod) k L # 分子 P A * C mod x^{k1} P [0] * (k 1) for n in range(k 1): acc 0 for j in range(n 1): acc (acc s[j] * C[n - j]) % mod P[n] acc return P, C实测下来从 8 项序列能稳定恢复到原来的 (1x)/(1-2x-3x^2)。这套闭环比单纯跑递推更有说服力能直接检验分子分母整体是否恢复正确。3.4 一个容易忽略的细节Q(0) 归一化严格说Padé 逼近或者有理函数重建的通常约定是 Q(0)1。BM 算法输出的连接多项式常数项就是 1天然满足这个要求。但如果你的输入数据来自别的地方比如线性方程组解出来的分母常数项不是 1那一定要先做归一化分子分母同时除以 Q(0)。这个看起来很蠢的坑在实际代码里非常常见。尤其是当你用浮点数运算时归一化能显著提升数值稳定性而在模运算下归一化就相当于用 Q(0) 的模逆乘一遍分子分母。要是忘了做后面所有比较、代入验证、画图都会错得非常诡异。4. 实战中的坑与调优清单4.1 边界情况零序列、长度不足、分母退化第一个边界序列全零。这时候 BM 算法会得到 L0C[1]对应的有理函数是什么0/1也就是恒零函数。这个结果是对的。但代码实现里要特别小心C 的长度只有 1后续所有依赖 C[j] 的地方都要防止越界。第二个边界序列长度 N 小于 2k。刚刚说了这种数据量不足的情况下结果不唯一。BM 算法本身会返回一个递推但它可能只是“过拟合”当前样本不代表序列的真实结构。比如说我随手取一个周期为 4 的序列只给前 3 项BM 可能会返回一个 1 阶递推但这个递推显然不能外推。处理办法是拿到 BM 结果后一定做外推验证用递推生成接下来的若干项跟真实新增的序列对比看看是否持续一致。第三个边界分母有重根或者 Q(x) 和 P(x) 有公因子。从算法角度来说BM 不关心这些它只负责找最短短递推。但如果在重建有理函数之后做部分分式分解重根会导致分解式出现 x - r 的幂次项这属于后续处理的范畴。需要提醒的是如果 Q 和 P 共享公因子那么这个有理函数实际上可以被约分你可以看到此时 BM 给出的递推阶数其实等于约分后的分母次数而不是约分前的所以输出始终是最简形式。4.2 浮点运算 vs 模素数运算怎么选这是一个非常现实的问题。我是两者都用过最后形成了一个经验判断除非应用场景明确要求实数系数比如系统辨识、控制论、连续信号否则优先用模素数。我整理了一张对比表维度浮点 BM模素数 BM有理数 BM数值稳定性低误差会随长度累积高完全无舍入误差高但分数膨胀严重速度快快慢实现难度中低中适用场景实数序列、数值分析竞赛、密码学、精确代数小规模精确计算浮点 BM 最容易踩的坑就是“毁灭性的精度灾难”。特别是递推阶数超过 20 之后误差会指数级放大经常出现递推系数算出小数点后好几位却依然无法正确外推的情况。这时候用模素数就非常舒服所有运算都是整数只要不溢出结果精确无比。选模数的时候尽量选大素数比如 998244353 这种它同时是 NTT 素数后面如果要优化卷积还能直接衔接。如果用浮点 BM额外建议每步对 delta 做一个阈值判断当 abs(delta) 小于某个 epsilon比如 1e-10就强制当作 0 处理。否则微小的舍入误差会被当成失配导致算法疯狂修正、递推长度快速膨胀最后输出一个高得离谱的“最短递推”。4.3 验证正确性的三个小实操做完一次求解怎么知道自己没写错我总结了三招从易到难。第一招外推验证。用求出的递推生成未来 10 项跟真实数据逐项比较。如果递推阶数是 k外推 10 项足够发现大部分问题。这个方法适合任何场景5 秒钟就能跑完。第二招随机自检。构造一个已知递推比如 s[n] 3·s[n-1] - 2·s[n-2] s[n-3]随机生成初始项然后生成 200 项喂给 BM看输出是否还原出原递推。这个自检能暴露绝大多数实现 bug尤其是初始化顺序和更新分支写反的问题。我当年第一次写 BM就是靠这个自检发现 C 和 B 更新先后顺序搞反了。第三招与有理函数闭环校验。用 rational_from_seq 恢复出 P 和 Q 之后直接计算 P/Q 的展开序列和原始序列对比。如果完全一致说明整条流水线没问题。这一步适合在正式处理数据前跑一次“冒烟测试”。4.4 性能边界与优化空间BM 算法的时间复杂度是 O(N^2)其中 N 是序列长度。如果递推阶数 k 远小于 N实际计算量大概在 O(N·k) 左右还是可以接受的。但 N 到十万、百万级别的时候O(N^2) 就比较吃力了。有几个可行的优化思路第一如果数据是在模素数环境中且递推阶数 k 比较大可以考虑把 BM 中的卷积部分用 NTT 加速。不过这会显著增加代码复杂度市面上的实现也少通常只有在 N 特别大、k 也特别大时才值得。第二如果只关系递推阶数而不关心系数可以用随机投影技巧随便取几个随机向量跟序列做内积然后在这些低维序列上跑 BM再用概率论保证结果的正确性。这种方法适合快速估计阶数但不适合精确恢复系数。第三在不行的情况下可以改用分块 BM 算法。这个我其实没有在实际项目中用过只是调研时见过相关论文。一般应用场景下老老实实写个 O(N^2) 的 BM配合合理的模数和数据预处理已经完全够用了。我个人在实际操作中的一点体会是最短线性递推式和有理函数重建这两个问题越是深入用越觉得它们不只是“算法题”而是一套理解序列和生成函数关系的语言。以前我遇到一串数字第一反应是去拟合多项式或者做最小二乘后来学会了先跑一遍 BM看看数据底层是不是真的存在一个低阶线性结构。这个习惯帮我避免了好多次无谓的过拟合。另外写 BM 的时候千万不要贪快一步一步把状态变量理清楚尤其是 B、b、m 这些保存的历史信息它们在更新分支里的使用顺序很容易出错。建议任何生产环境里的代码都要先挂一个随机自检再放到真实数据上跑。这套流程我重复了几十次几乎成了信条拿到序列先 BM再重建有理函数最后外推验证。下次你手里碰上一串来路不明的数字不妨也先试试这条路。