很多人小时候应该都玩过这样一个游戏：给你一串数字，让你猜下一个是什么。比如 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 是递推系数。斐波那契数列就是 k=2，c1=1，c2=1 的特例。当然，系数不一定是整数，可以是实数、复数、有理数，甚至模某个素数 p 意义下的整数。

满足某个 k 阶线性递推的序列，全体构成了一个 k 维线性空间。这个观点特别重要，因为“最短”这个概念其实就是在问：给定序列的前 N 项，能不能找到一个最小的 k，以及对应的系数，使得从第 k+1 项开始，每一项都严格符合这个递推关系。注意这里说的是前 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[j+m] 的更新，区别只在于是否更新 L、B、b、m；三是我最后把 C 裁到 L+1 长度，保证输出干净。

如果序列本身是整数序列，且不需要模运算，可以把所有运算放在有理数域里做，用 Fraction 类型代替整数。但一般不建议这么做，因为分数运算会让时间复杂度膨胀得很厉害，而且容易溢出大整数的范围。实际工程里，要么用模大素数（比如 998244353 这种 NTT 素数）兜底，要么用浮点 BM 算法，后面我在避坑部分会细说。

2.4 手推一个例子：斐波那契数列

光看代码可能还是有点抽象，我们手动推一遍最简单的斐波那契序列。设 s=[1, 1, 2, 3, 5, 8]，取模数 1000000007。

初始：C=[1], B=[1], L=0, m=1, b=1。

i=0，s[0]=1。计算 d = 1。2L=0 <= i=0，所以进入长度更新分支。coef=1，C 变成 [1,-1]（模意义下是 [1, mod-1]），L = 0+1-0=1，B=[1]，b=1，m=1。现在 C 表示的递推是 s[n]=s[n-1] 吗？其实是常数数列的递推。

i=1，s[1]=1。d = s[1] + C[1]*s[0] = 1 + (-1)*1 = 0。没有失配，m=2。

i=2，s[2]=2。d = s[2] + C[1]s[1] = 2 - 1 = 1，非零。2L=2 <= i=2，所以又进入长度更新分支。coef = d * b^{-1} = 1。当前 C=[1, -1]，B=[1]，m=2。更新后 C = [1, -1] - x^2[1] = [1, -1, -1]，也就是 1 - x - x^2。新长度 L = 2+1-2 = 1？等等，这里按公式 L = i + 1 - L = 3 - 1 = 2。对，是 L=2。B=[1,-1]，b=1，m=1。

i=3，s[3]=3。用 C=[1,-1,-1] 预测：d = s[3] + (-1)*s[2] + (-1)*s[1] = 3 - 2 - 1 = 0，正确。

i=4，s[4]=5。d = 5 - 3 - 2 = 0，正确。

i=5，s[5]=8。d = 8 - 5 - 3 = 0，正确。

最终得到 L=2，C=[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^{k+1}

这里只要取到 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_{j=0}^{n} s[j] · C[n-j], for n = 0..k

直接算就行，复杂度 O(k^2)。

3.2 最少需要多少项：2k 原则

这里必须强调一个边界条件：要唯一确定一个 k 阶递推，序列至少需要给多少项？

答案是需要大约 2k 项。为什么？因为 k 阶递推有 k 个未知系数，而递推关系对序列的第 k+1 到第 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，也就是 L=2。然后再利用前 2 项算分子：p0 = s0 = 1；p1 = 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^{k+1} 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 项序列能稳定恢复到原来的 (1+x)/(1-2x-3x^2)。这套闭环比单纯跑递推更有说服力，能直接检验分子分母整体是否恢复正确。

3.4 一个容易忽略的细节：Q(0) 归一化

严格说，Padé 逼近或者有理函数重建的通常约定是 Q(0)=1。BM 算法输出的连接多项式常数项就是 1，天然满足这个要求。但如果你的输入数据来自别的地方，比如线性方程组解出来的分母常数项不是 1，那一定要先做归一化：分子分母同时除以 Q(0)。

这个看起来很蠢的坑，在实际代码里非常常见。尤其是当你用浮点数运算时，归一化能显著提升数值稳定性；而在模运算下，归一化就相当于用 Q(0) 的模逆乘一遍分子分母。要是忘了做，后面所有比较、代入验证、画图都会错得非常诡异。

4. 实战中的坑与调优清单

4.1 边界情况：零序列、长度不足、分母退化

第一个边界：序列全零。这时候 BM 算法会得到 L=0，C=[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 最容易踩的坑就是“毁灭性的精度灾难”。特别是递推阶数超过 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，再重建有理函数，最后外推验证。下次你手里碰上一串来路不明的数字，不妨也先试试这条路。