我最早碰到“最短线性递推式求解”这个概念是在一次流量密码分析的任务里。当时手里只有一串截获的密钥流比特长度大概一千出头看起来完全随机但直觉告诉我底层可能藏着一个LFSR线性反馈移位寄存器。用什么办法把这个LFSR的反馈多项式抽出来答案是用Berlekamp-Massey算法输入那串比特输出一条最短线性递推式。这东西很多做密码、做编码、做序列分析的朋友都用过但真正把它和“有理函数重建”放在一起看的文章很少。其实这两个问题本质上是一对孪生兄弟——一个从序列到递推式一个从多项式剩余到分式底层代数结构出奇地一致。 这篇文章就把这两件事掰开揉碎讲清楚。内容包括线性递推和有理函数之间的等价关系、Berlekamp-Massey算法的完整原理与手算示例、基于扩展欧几里得的有理重建方法、两套方法的Python实现以及我在工程实践中踩过的一堆坑。适合正在做密码工程、纠错码译码、信号处理或者刷算法题时被“最短线性递推”卡住过的同学。 ## 1. 问题建模两个看似不搭边的算法其实共用同一套代数骨架 ### 1.1 线性递推式与生成函数序列背后藏着分式 先明确一下什么叫做线性递推式。给定一个序列 \(s_0, s_1, s_2, \dots\)如果存在一组常数 \(c_1, c_2, \dots, c_L\)使得对任意 \(n \ge L\) 都有 \[ s_n c_1 s_{n-1} c_2 s_{n-2} \dots c_L s_{n-L} 0 \] 那么就说这个序列满足一个 \(L\) 阶线性递推。这里 \(L\) 越小说明序列的结构越“简单”。最典型的就是Fibonacci数列它满足 \(s_n - s_{n-1} - s_{n-2} 0\)所以最短递推长度是2。 为什么要关心“最短”因为实际问题里我们拿到的序列往往是截断的、带噪声的或者是从某个有限状态机里采出来的有限长度观测值。如果底层确实存在一条递推式但我们不知道阶数那么直接猜一个比较大的 \(L\) 可能也能拟合但会过拟合把噪声也当成结构学进去。最短递推式的意义在于在所有能解释这段序列的递推关系中找出阶数最小的那一个它通常对应最本质的底层结构。 这里有一个特别重要的数学视角一个无穷序列满足某个 \(L\) 阶线性递推当且仅当它的生成函数是一个分母次数不超过 \(L\) 的有理函数。所谓生成函数就是把序列写成形式幂级数 \[ G(x) s_0 s_1 x s_2 x^2 \dots \] 如果序列满足上述递推关系那么可以推出 \[ G(x) \frac{P(x)}{Q(x)},\quad Q(x) 1 c_1 x c_2 x^2 \dots c_L x^L \] 分子 \(P(x)\) 是一个次数小于 \(L\) 的多项式由初始的 \(s_0, \dots, s_{L-1}\) 决定。换句话说找序列的最短线性递推式等价于在“分母次数最小”的约束下找一个能生成这个序列的有理函数。 ### 1.2 有理函数重建从同余条件还原分式 再说有理函数重建。它的问题形如已知一个模多项式 \(m(x)\)以及一个剩余 \(u(x)\)要找两个次数受限的多项式 \(f(x), g(x)\)满足 \[ f(x) \equiv g(x) \cdot u(x) \pmod{m(x)} \] 同时要求 \(\deg f d_1\)\(\deg g d_2\)而且 \(g\) 和 \(m\) 互素。看起来完全不像序列问题但它其实在密码学里极其常见。举个例子在RSA的某些侧信道攻击中攻击者能获得某个秘密分数的模 \(N\) 剩余而这个秘密分数本身是一个小分子、小分母的有理数。要从模剩余中把分子分母还原出来就是典型的有理重建问题。 这个问题的求解核心是扩展欧几里得算法。过程是对 \(m(x)\) 和 \(u(x)\) 做一系列带余除法同时维护两个系数多项式 \(s_i(x), t_i(x)\)使得每一步都有 \[ r_i(x) s_i(x) \cdot m(x) t_i(x) \cdot u(x) \] 当某个中间余式 \(r_i(x)\) 的次数降到我们期望的分子次数界以内时就取 \(f r_i, g t_i\)。这一步和连分数的求解本质上是一回事——欧几里得算法在多项式环上产生的商序列就是有理函数连分数展开的系数迭代到某一步时得到的收敛子就是我们要的分式。 ### 1.3 两个问题的统一看同余式的不同角度 把1.1和1.2放到一起看就能发现一个漂亮的对偶关系。 从序列重建递推式其实可以看作在环 \(\mathbb{F}[x]/(x^N)\) 中做有理重建给定截断序列的前 \(N\) 项记 \(U(x) s_0 s_1 x \dots s_{N-1} x^{N-1}\)找一个分母次数尽量小、分子次数也尽量小的有理函数 \(P(x)/Q(x)\)使得 \[ P(x) \equiv Q(x) \cdot U(x) \pmod{x^N} \] 这不就是1.2里的同余式只是模数换成了 \(x^N\) 而已。从序列角度看存在 \(L\) 阶递推意味着 \(Q(x) \cdot G(x)\) 的前 \(N\) 项都消掉了只保留高次项所以 \(Q(x) \cdot U(x) \equiv P(x) \pmod{x^N}\)。于是“最短线性递推式求解”和“有理函数重建”共享同一套扩展欧几里得骨架只是模数、停止条件、以及“最短”的定义不同。 这也是为什么我强烈建议你把两个算法一起学会了一个另一个就是改改停止条件的事。 ## 2. 最短线性递推Berlekamp-Massey算法的原理与手算 ### 2.1 增量算法每一步只修正必要误差 Berlekamp-Massey算法下称BM算法是一个增量算法。它逐个读入序列项维护当前前缀的最短递推多项式 \(C(x)\)以及一个记录上一次出现“匹配失败”时递推式的副本 \(B(x)\) 和对应位移下标 \(m\)。 算法的核心思路是当读入第 \(n\) 项时先用当前递推式预测它的值计算误差 \[ \Delta s_n \sum_{i1}^{L} c_i s_{n-i} \] 如果 \(\Delta 0\)说明当前递推式在这个点上仍然有效直接继续。如果 \(\Delta \neq 0\)说明递推式被打破了必须修正。修正的方式不是从头重算而是构造一个新的递推式 \[ C_{\text{new}}(x) C(x) - \frac{\Delta}{\Delta_{\text{old}}} x^{n - m} B(x) \] 其中 \(\Delta_{\text{old}}\) 是上一次修正时的误差\(B(x)\) 是当时的旧递推式。这里的直觉是把旧递推式的“上一次错误模式”平移到现在用适当的系数叠加到当前递推式上恰好可以抵消新出现的误差。这个技巧很像线性代数里的递推校正每一步都只调整必要的维度。 ### 2.2 一个完整手算示例从短序列推递推式 我们手动跑一遍BM帮助理解。假设在 \(\mathbb{F}_2\) 上处理序列 \[ 1, 0, 0, 1, 0, 1 \] 初始状态\(C(x) 1\)\(B(x) 1\)\(L 0\)\(m 1\)上次误差 \(\Delta_{\text{old}} 1\)。 - 第0位 \(s_0 1\)当前 \(L0\)预测值为0误差 \(\Delta 1\)。发现误差非零且 \(2L 0 \le n 0\)所以更新递推式。构造 \(C_{\text{new}}(x) 1 1 \cdot x^{0} \cdot 1 / 1 1 x\)同时更新 \(B(x) 1\)\(m 1\)\(L 1\)\(\Delta_{\text{old}} 1\)。 - 第1位 \(s_1 0\)\(C(x) 1 x\) 给出的预测是 \(s_1 1 \cdot s_0 0 1 1\)与真实值0不符误差 \(\Delta 1\)。此时 \(2L 2 n 1\)所以不增加阶数只需修正多项式。用公式\(C_{\text{new}}(x) C(x) - \frac{\Delta}{\Delta_{\text{old}}} x^{n - m} B(x) (1x) - x^{0} \cdot 1 0\)这里要小心在 \(\mathbb{F}_2\) 上减法等于加法计算得 \(C_{\text{new}}(x) 1x x^{0} \cdot 1 1x1 x\)。不过 \(C(x)\) 的常数项理论上是1这里出现常数项为0的多项式原因是序列前两位都是常数模式 \(s_01\)说明最短递推其实是 \(s_n 0\)除首项外算出来 \(x\) 等价于递推长度1且系数0。这个例子确实有点反直觉工程实现时遇到常数项归零需要特殊处理。这里为了演示流程继续硬算。 - 第2位 \(s_2 0\)用 \(C(x) x\) 预测为0真实值0误差0递推式不变。 - 第3位 \(s_3 1\)预测0真实1误差1。此时 \(2L 2 \le n 3\)更新阶数。构造新递推式最终得到长度更长的递推。完整手算比较繁琐建议直接跑代码验证。 说实话新手第一次手算BM非常容易晕因为下标和对齐关系太琐碎。我的经验是先跑通代码再对着代码断点看每一步的 \(C, B, m, L\) 变化比单纯手算理解快得多。 ### 2.3 关键经验至少要2L个观测值 BM算法最容易被忽略的一点是数据量需求。给定长度为 \(N\) 的序列BM算法能输出一条长度不超过 \(\lfloor N/2 \rfloor\) 的递推式但它只在“序列长度足够长”时才能保证这条递推式唯一。 严格说如果你知道底层最短递推式长度为 \(L\)那么至少需要连续 \(2L\) 个观测值BM算法才能准确恢复这条递推式。少于 \(2L\) 项时解不唯一甚至可能输出一条看起来合理但不本质的短递推式。我在实际项目里一般会留出冗余如果预期递推阶数是 \(L\)会采集至少 \(2L 20\) 个点防止边界效应和噪声影响。 这个性质也解释了为什么BM算法在流密码分析里那么有用线性反馈移位寄存器生成的密钥流只要你能拿到超过 \(2L\) 的连续明文与密文对齐片段就能用BM以 \(O(N^2)\) 的代价恢复整个LFSR结构等效密钥量直接归零。 ## 3. 有理函数重建扩展欧几里得与连分数的双重面纱 ### 3.1 从一次带余除法到整个分式恢复 有理函数重建的核心执行方案是扩展多项式欧几里得算法。具体流程如下。 输入模多项式 \(m(x)\)、剩余 \(u(x)\)、分子次数界 \(d_f\)、分母次数界 \(d_g\)。 初始化 \[ \begin{aligned} r_0 m(x), s_0 1, t_0 0 \\ r_1 u(x), s_1 0, t_1 1 \end{aligned} \] 迭代 1. 用 \(r_{i-2}\) 除以 \(r_{i-1}\)得到商 \(q_i\) 和余式 \(r_i\)。 2. 更新 \(s_i s_{i-2} - q_i s_{i-1}\)\(t_i t_{i-2} - q_i t_{i-1}\)。 3. 检查是否满足停止条件\(\deg r_i d_f\) 且 \(\deg t_i d_g\)。满足则停止返回 \((f, g) (r_i, t_i)\)。 因为初始时 \(r_0 m\) 是模多项式的倍数所以整个迭代过程中始终有 \[ r_i s_i m t_i u \] 把同余条件 \(f \equiv g u \pmod m\) 代入就是 \(r_i \equiv t_i u \pmod m\)。所以 \((r_i, t_i)\) 自然满足重建方程。 ### 3.2 与连分数的关系你就是在一层层逼近那个分式 理解有理重建最直观的方式是连分数。扩展欧几里得的商 \(q_1, q_2, \dots\) 恰好是有理函数 \(u(x)/m(x)\) 的连分数展开系数。迭代到第 \(i\) 步时比值 \(r_i / t_i\) 就是连分数的第 \(i\) 个收敛子convergent。 收敛子的性质是交替靠近真实值而且从某一步开始分子分母次数会同时变小。当次数降到预设界以内时我们就得到了一个满足同余方程且“足够简单”的分式表示。这解释了为什么停止条件要同时看 \(\deg r_i\) 和 \(\deg t_i\)——只看一个会得到奇奇怪怪的退化结果。 在密码攻击里一个典型场景是已知某个秘密值 \(k\) 对模数 \(N\) 的剩余为 \(r\)且秘密值形如 \(k p/q\)其中 \(p, q\) 都很小。这时取 \(m N\)\(u r\)做有理重建得到的 \(p, q\) 往往就是原始秘密。这种思路在HNPhidden number problem攻击、RSA的部分密钥泄露攻击里反复出现。 ### 3.3 唯一性边界不是随便给个界都能成功 有理重建不是总能成功它依赖于分子分母次数界和模数次数之间的关系。一个常被提及的充分条件是 \[ \deg f \deg g \deg m \] 并且我们额外要求 \(\deg f d_f\)、\(\deg g d_g\)且 \(d_f d_g \le \deg m\)。这个条件保证了扩展欧几里得迭代到某一步时解是唯一的。 如果 \(d_f d_g\) 太接近甚至超过 \(\deg m\)就可能出现多个分式都满足同余方程重建结果就不确定。我自己在实现RSA攻击脚本时吃过这个亏当时把分子分母界之和设成了 \(N\) 的比特数左右结果每次跑出来的分式都不同排查了半天才发现是唯一性边界被打破了。经验公式是界各留至少10%的余量能显著提升稳定性。 ## 4. 代码级实操一条管道跑通两个任务 ### 4.1 用不到40行Python实现Berlekamp-Massey 下面给出一个可直接用于 \(\mathbb{F}_p\)大素数域的BM实现它是很多密码学库的简化版本。注意用系数列表表示多项式低次项在前。 python def berlekamp_massey(s, p): # s: 序列元素在模p域上 # 返回递推多项式 C(x) 1 c1 x ... cL x^L C [1] B [1] L 0 m 1 b_old 1 for n in range(len(s)): # 计算当前误差 delta d s[n] for i in range(1, L 1): d (d C[i] * s[n - i]) % p if d 0: m 1 continue # 记录旧C T C[:] coef d * pow(b_old, p - 2, p) % p if len(C) len(B) m: C [0] * (len(B) m - len(C)) for i in range(len(B)): C[i m] (C[i m] - coef * B[i]) % p if 2 * L n: L n 1 - L B T b_old d m 1 else: m 1 return C[:L 1]代码里的关键细节有两个。一是更新B的时机只有满足 (2L \le n) 时才需要把旧C存到B里因为这时递推阶数真正增加了二是用费马小定理求逆时要保证 (b_{\text{old}}) 非零正常情况下会被BM算法的性质保证但如果你往序列里掺了零除元素程序会直接报错这是第一个要检查的坑。4.2 有理重建多项式扩展欧几里得的实现接下来是重建分子的代码。实现上要处理多项式除法、次数比较、以及归一化。def poly_divmod(a, b): # 多项式带余除法返回 (q, r)系数低次在前 a a[:] b b[:] inv_lc pow(b[-1], p - 2, p) q [0] * (len(a) - len(b) 1) while len(a) len(b) and any(a): shift len(a) - len(b) coef a[-1] * inv_lc % p q[shift] coef for i in range(len(b)): a[i shift] (a[i shift] - coef * b[i]) % p while len(a) 0 and a[-1] 0: a.pop() return q, a def rational_reconstruct(u, m, df, dg, p): # u: 剩余多项式m: 模多项式 # 返回 f, g满足 f g * u mod mdeg f df, deg g dg r0, r1 m[:], u[:] t0, t1 [0], [1] while len(r1) 0: q, r poly_divmod(r0, r1) # 判断是否停止 if len(r) - 1 df and len(t1) - 1 dg: # 返回时需要保证 g 首一 if t1[-1] ! 1: inv pow(t1[-1], p - 2, p) r [x * inv % p for x in r] t1 [x * inv % p for x in t1] return r, t1 r0, r1 r1, r t0, t1 t1, [(t1[i] - ((q[0] * t0[0]) if len(q) else 0)) for i in range(len(t1))] # 注意这里t更新需要完整多项式乘法简化版只适合低次完整版请用通用mul_sub return None上面的t更新部分写得很简化只是为了展示骨架。工程里我建议直接用SageMath的rational_reconstruct函数它会自动处理所有边界情况自己实现的话容易被多项式长度变化绕晕。4.3 实战实例从已知生成函数反推分子分母假设我们有生成函数[ G(x) \frac{1 x}{1 - x - x^2} ]这个分式展开后就是Fibonacci数列乘以某个系数。取前8项展开[ 1, 2, 3, 5, 8, 13, 21, 34 ]用上面的BM跑这段序列得到递推多项式 (1 p - p^2)在模一个大素数域下即 (s_n s_{n-1} s_{n-2})。再把这些项拼成 (U(x))令 (m(x) x^8)调用有理重建设置 (\deg f 1)(\deg g 3)程序返回[ f(x) 1 x,\quad g(x) 1 - x - x^2 ]和原始分式完全一致。这个测试验证了整套管道的正确性以后你怀疑某个序列有低阶递推结构时可以照这个流程跑一遍很快就能验证。5. 工程实践中的高频坑位与排查指南5.1 BM算法在非域环境下的崩溃BM算法的正确性依赖“域”结构每一步更新都要用误差 (\Delta) 除以旧误差 (\Delta_{\text{old}})也就是需要做除法。如果你处理的是模合数环比如模 (2^{32})除法不一定有逆元算法直接失效。我见过很多同学拿BM去跑模 (2^k) 下的随机序列结果输出的“递推式”经常只有长度1因为算法在无法做除法时的行为完全失控。解决办法是要么把问题换到素数域上处理比如用一个模大素数要么使用专门针对环设计的BM变体。实际密码分析里LFSR这类线性结构都定义在 (\mathbb{F}_2) 上所以BM基本够用没必要硬趟合数环的浑水。5.2 有理重建的归一化问题有理重建返回的分式不唯一分子分母同时乘以同一个非零常数得到的分式在数学上等价。工程上必须固定一个归一化约定否则每次跑出来的结果形式不一样后续比对会很痛苦。我的习惯是归一化分母强制分母多项式首项系数为1。实现时在返回前检查t1[-1]如果不是1就对分子分母整体乘以它的逆元。这个细节看着小但在批量验证攻击结果时能省大量调试时间。5.3 数据量不足时的迷惑性输出BM算法对长度不足的序列会输出“一条”递推式但这条递推式可能是错的。你拿它去预测后面的值很快就对不上。我踩过一次坑某次CTF题目里给了一个长度很短的序列我直接用BM得到了一个看着很短的递推式结果后面验证全错。后来才知道题目设计者故意把序列截短了让基于BM的唯一性条件不成立这时候必须用更复杂的格基约化方法才能解出真正的结构。经验是在任何场景下都先算一下“当前序列能否唯一定义最短递推式”。判据很简单BM输出的长度 (L)要求序列长度 (N \ge 2L)。如果 (N 2L)得到的递推式只能算“过拟合参考值”别直接拿去生产环境用。5.4 性能优化什么时候改用NTL/Sage自己实现的BM和有理重建在域大小适中、序列长度几千以内时速度没问题。但如果你要在超大素数域上处理上万长度的序列或者模多项式次数很高纯Python实现会慢到让你怀疑人生。这时有两个优化方向。第一用FFT加速多项式乘法和除法复杂度可以从 (O(n^2)) 降到 (O(n \log n))但实现复杂度高适合有充足时间打磨的场景。第二直接用现成的库SageMath的berlekamp_massey和rational_reconstruct都经过高度优化NTL库里的RR系列函数也很快。我个人的准则是项目原型阶段先用Python/Sage跑通正确性需要上线或做大规模扫描时再迁移到C/FFT实现不要一上来就手搓高性能版本。5.5 快速排查清单我把这几年遇到的高频问题整理成一张表卡住的时候可以按这个顺序排查现象原因处理方式BM输出长度一直为1序列长度太短或场不正确检查是否满足 (N \ge 2L)确认运行在域上递推式预测后续项全错数据量不足导致唯一性失效增加观测值或改用格基方法有理重建返回平凡解停止条件设置过松收紧分子分母次数界保证 (d_f d_g \le \deg m)重建结果每次跑不一样没有归一化分母强制分母首一代码里求逆报错分母多项式与模数不互素检查是否存在公因子必要时去掉公因子再重建结果看起来正确但符号不对归一化方向选错统一用分母首一而不是分子首一6. 从建模到生产我的一些补充建议最后分享几个个人体会不算总结就是实际干活攒下的经验。第一不要把BM和有理重建当成两个孤立算法。训练自己一看到“给定序列求结构”就想到生成函数、想到模 (x^N) 同余、想到扩展欧几里得这条链路会帮你快速定位到正确工具。第二实现时优先保证停止条件和归一化正确再去优化常数因为这两个地方最容易出隐蔽的逻辑错误。第三处理密码学场景时永远假设观测数据可能有噪声或者被恶意截短先验证数据量是否足够不要一上来就跑算法。另外如果你在做流密码分析建议把BM算法和已知明文攻击思路结合起来经常能快速恢复LFSR初态和反馈多项式。结合之前提到的有理重建还可以处理带有小分子分母的秘密恢复问题。这两个工具组合起来覆盖了不少CTF和真实协议的破解场景实在值得花一个下午把代码跑熟。