上个月我在做流密码序列的线性复杂度分析时同时撞上了两个看起来八竿子打不着的需求一边要用 Berlekamp-Massey 算法求序列的最短线性递推式一边要从一个幂级数的前若干项恢复它的分式表达。做着做着我才反应过来这两个问题在代数结构上根本是同一件事——递推式本质上就是一个有理函数的分母而有理函数重建本质上就是在求一个最短的、能解释已知观测的递推模型。这篇文章就是想把这条线彻底捋清楚并给出可以直接抄走的实现代码和踩坑记录。这篇内容适合以下几类读者做密码学/编码理论的经常跟 LFSR、BCH 码、伴随式打交道的人做信号处理、语音建模想理解线性预测背后代数原理的人还有写序列预测工具、想从少量数据里恢复生成函数的人。我会从数学原理讲到代码实现中间穿插完整的手动演算和验证方法保证看完能直接上手。1. 两个名字一个本质递推式与有理函数互为表里1.1 线性递推到底在描述什么先明确一个基本概念。所谓线性递推式指的是序列中的每一项都能由前面若干项的线性组合表示。比如 Fibonacci 序列0, 1, 1, 2, 3, 5, 8, 13, ...满足s[n] s[n-1] s[n-2]这就是一个二阶线性递推式。如果写成更一般的形式s[n] c1 * s[n-1] c2 * s[n-2] ... cL * s[n-L] 0那这个L就叫递推阶数也叫线性复杂度。最短线性递推式就是让这个L尽可能小同时又能在给定前缀上完全成立的递推关系。这个最短在工程里不是洁癖。一个序列如果能被一个短递推式完整描述意味着它的信息量其实很小——你可以用极少的参数存储它、预测它、甚至判断它是否安全。流密码里的 LFSR 序列就是典型的线性递推序列如果它的最短递推阶数太低攻击者只要收集 2L 个连续比特就能完全重建整个序列这在密码学里是致命的弱点。1.2 有穷递推的背后是一个有理函数现在把序列看成一个幂级数也就是生成函数S(x) s0 s1*x s2*x^2 s3*x^3 ...如果序列满足上面的 L 阶递推式我定义一个连接多项式C(x) 1 c1*x c2*x^2 ... cL*x^L然后把C(x)和S(x)相乘。稍微展开一下就能发现当次数大于等于 L 时每一项的系数恰好就是递推式的左端也就是 0。换句话说C(x) * S(x) P(x)其中P(x)是一个次数严格小于 L 的多项式。于是S(x) P(x) / C(x)这是一个实实在在的有理函数——分子是次数小于 L 的多项式分母就是递推式的连接多项式。拿 Fibonacci 序列举例0, 1, 1, 2, 3, 5, 8, 13的生成函数是S(x) x / (1 - x - x^2)它对应的递推式s[n] s[n-1] s[n-2]连接多项式C(x) 1 - x - x^2分子P(x) x一切严丝合缝。1.3 为什么最短和重建在工程里都重要一旦你接受了递推式 有理函数分母这个视角最短线性递推式求解和有理函数重建就成了同一枚硬币的两面给定序列的前 N 项求最短递推式本质上是找一个分母次数最小的有理函数去解释已知数据给定幂级数的前 N 项恢复分子分母都有次数约束的有理函数本质上也是在找一个足够短的模型。这种用最少参数解释观测数据的诉求在编码理论里是纠错在密码学里是安全度量在信号处理里就是线性预测。理解了这层关系后面所有算法都是在同一张地图上找路径。2. Berlekamp-Massey 求解最短线性递推式从手动演算到可运行实现2.1 算法维护的到底是什么Berlekamp-Massey后面统一叫 BM算法的核心是维护一个当前最优的连接多项式C(x)并随着新数据逐项进入不断修正它。算法涉及的变量有五个我按直觉解释一下C当前的连接多项式C[0]恒为 1递推关系是sum(C[i] * s[n-i]) 0L当前递推阶数也就是C的实际长度减 1B上一次发生过长度更新时的连接多项式可以理解为旧模型的备份b上次长度更新时对应的误差值m从上次长度更新到现在走过的步数用来记录旧模型落后了多远。每次读入一个新的s[n]算法先算当前模型在这一项上的预测误差d C[0]*s[n] C[1]*s[n-1] ... C[L]*s[n-L]如果d为 0说明当前递推式靠得住什么都不改只要m 1。如果d不为 0说明模型失配了要修正C修正的方式是拿旧模型做一个补偿C_new C - (d / b) * x^m * B这里的x^m * B表示把旧模型整体向高次方向移动 m 位。至于是否需要增大长度看一个关键条件2L n——如果当前长度已经不够覆盖观察到的模式算法就把长度更新为n 1 - L同时把旧C存到B、当前误差存到b、m重置为 1否则只修正系数长度不变m 1。2.2 手动演算用 Fibonacci 序列跑一遍 BM纸上谈兵没有用我拿 Fibonacci 序列0, 1, 1, 2, 3, 5, 8完整走一遍保证你看完就懂每一步在干嘛。初始状态C[1]B[1]L0m1b1。ns[n]误差 d是否更新更新后 CL00C[0]*s[0] 0否[1]011C[0]*s[1] 1是2Ln[1, 0, -1]221C[0]*s[2]C[1]*s[1]C[2]*s[0] 1是2Ln[1, -1, -1]2322 - 1 - 1 0否[1, -1, -1]2433 - 2 - 1 0否[1, -1, -1]2555 - 3 - 2 0否[1, -1, -1]2688 - 5 - 3 0否[1, -1, -1]2注意看n1到n2的变化。n1时因为当前长度 L0 根本解释不了非零序列触发长度翻倍的更新C变成[1, 0, -1]这相当于在猜下一项等于上上项。n2时这个猜测错了误差 d1但因为2L n不需要扩大长度只是用修正项把C调整成[1, -1, -1]。之后连续 4 项误差都是 0说明这个递推式已经稳定生效。最终结果C [1, -1, -1]即递推式s[n] s[n-1] s[n-2]完全正确。2.3 可直接复现的 Python 实现下面给两版实现。第一版用fractions.Fraction做精确有理数运算不涉及求逆适合学习和验证第二版是模素数域版本竞赛和密码学场景更常用。有理数域版本from fractions import Fraction def berlekamp_massey(seq): C [Fraction(1)] B [Fraction(1)] L 0 m 1 b Fraction(1) for n in range(len(seq)): d Fraction(0) for i in range(L 1): d C[i] * seq[n - i] if d 0: m 1 continue T C[:] coef d / b if len(C) len(B) m: C [Fraction(0)] * (len(B) m - len(C)) for j in range(len(B)): C[j m] - coef * B[j] if 2 * L n: L n 1 - L B T b d m 1 else: m 1 return C[:L 1], L模素数域版本def berlekamp_massey_mod(seq, p): C [1] B [1] L 0 m 1 b 1 for n in range(len(seq)): d 0 for i in range(L 1): d (d C[i] * seq[n - i]) % p if d 0: m 1 continue T C[:] coef d * pow(b, -1, p) % p 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]) % p if 2 * L n: L n 1 - L B T b d m 1 else: m 1 return C[:L 1], L验证方法很简单。拿前面的 Fibonacci 序列跑一下seq [0, 1, 1, 2, 3, 5, 8, 13] C, L berlekamp_massey(seq) print(C, L) # [Fraction(1, 1), Fraction(-1, 1), Fraction(-1, 1)] 2然后把C还原成递推公式算下一个值def predict_next(seq, C, L): return -sum(C[i] * seq[len(seq) - i] for i in range(1, L 1)) print(predict_next(seq, C, L)) # 21提示pow(b, -1, p)需要 Python 3.8 以上版本。如果用的是老版本自己写一个扩展欧几里得求逆函数就行。3. 有理函数重建用扩展欧几里得实现 Padé 逼近3.1 模 x^N 视角与 Padé 逼近有理函数重建的问题可以这么描述已知幂级数f(x) a0 a1*x a2*x^2 ...的前 N 项系数想找一个有理函数p(x) / q(x)使得它和f(x)在前 N 项完全一致同时分子、分母的次数都受到约束。这个问题的经典名字叫 Padé 逼近和整数上的有理重建结构完全一致整数版本是已知a mod m找p/q且|p|, |q|有界这里只是把整数换成了多项式环把模m换成了模x^N。为什么能用扩展欧几里得做因为扩展欧几里得在每一步都维护一个组合等式r u * x^N v * f两边同时模掉x^N就得到r ≡ v * f (mod x^N)也就是说每一步的余式r和组合系数v天然构成一个合法的重建候选令q v、p r就有q * f ≡ p (mod x^N)。随着欧几里得算法不断进行r的次数越来越低v的次数越来越高我们只需要找到一对同时满足次数约束的(v, r)就行。3.2 扩展欧几里得重建法的实现和自检我直接给出一个可以运行的 Python 实现。这里用低次到高次的列表表示多项式所有系数统一用Fraction保证精确运算。from fractions import Fraction as F def poly_trim(p): while p and p[-1] 0: p.pop() return p def poly_deg(p): p poly_trim(p[:]) return len(p) - 1 if p else -1 def poly_sub(a, b): n max(len(a), len(b)) res [F(0)] * n for i in range(len(a)): res[i] a[i] for i in range(len(b)): res[i] - b[i] return poly_trim(res) def poly_mul(a, b): if not a or not b: return [] res [F(0)] * (len(a) len(b) - 1) for i, ai in enumerate(a): for j, bj in enumerate(b): res[i j] ai * bj return poly_trim(res) def poly_divmod(a, b): a [F(v) for v in a] b poly_trim([F(v) for v in b]) if not b: raise ZeroDivisionError(zero divisor) q [F(0)] * max(0, len(a) - len(b) 1) while len(a) len(b): coef a[-1] / b[-1] deg_diff len(a) - len(b) q[deg_diff] coef for i in range(len(b)): a[deg_diff i] - coef * b[i] a poly_trim(a) return poly_trim(q), a def rational_reconstruct(f, N, deg_q_max, deg_p_max): f [F(v) for v in f] if len(f) N: f [F(0)] * (N - len(f)) else: f f[:N] xN [F(0)] * (N 1) xN[N] F(1) r0, r1 xN, f v0, v1 [], [F(1)] while poly_deg(r1) 0: q, r poly_divmod(r0, r1) v_next poly_sub(v0, poly_mul(q, v1)) r0, r1 r1, r v0, v1 v1, v_next if poly_deg(r0) deg_p_max and poly_deg(v0) deg_q_max: prod poly_mul(v0, f) diff poly_sub(prod, r0) if poly_deg(diff) N: q_out, p_out v0, r0 if q_out and q_out[0] ! 0: inv F(1) / q_out[0] q_out poly_trim([c * inv for c in q_out]) p_out poly_trim([c * inv for c in p_out]) return q_out, p_out return None用 Fibonacci 序列验证一下。前面提到它的生成函数是x / (1 - x - x^2)于是设分子次数界为 1、分母次数界为 2seq [0, 1, 1, 2, 3, 5, 8, 13] q, p rational_reconstruct(seq, N8, deg_q_max2, deg_p_max1) print(q) # [Fraction(1, 1), Fraction(-1, 1), Fraction(-1, 1)] print(p) # [Fraction(0, 1), Fraction(1, 1)]输出q [1, -1, -1]即1 - x - x^2p [0, 1]即x与理论完全一致。代码里已经内置了自检也就是验证q * f - p的前 N 项是否全部为零这一步很关键能避免把某个不满足条件的中间余式当成结果输出。提示如果q[0]为 0说明重建出来的有理函数在x0处有极点和幂级数本身有定义矛盾。正常条件下不会出现真遇到就要回头检查输入数据和次数约束是不是给错了。3.3 两条路线殊途同归BM 也能直接做有理重建细心的读者会发现既然 BM 求出的连接多项式C直接就是生成函数的分母那我根本不需要跑扩展欧几里得也能完成有理函数重建。具体做法是用 BM 从序列前 N 项求出C(x)和阶数L取序列的前 L 项作为「初始值」把C(x)和S(x)的截断多项式相乘只保留次数小于 L 的部分得到分子P(x)。拿前面的结果手动验证C [1, -1, -1]序列前 2 项是[0, 1]乘积(1 - x - x^2) * (0 x ...)的低次部分正好是x所以P(x) x。BM 和 Padé 重建的关系可以这样概括BM 只约束分母次数即递推阶数不显式约束分子次数而扩展欧几里得版本的 Padé 重建可以同时控制分子、分母的次数界。当分子次数确实小于分母次数时两者给出相同结果但应用场景各有侧重。4. 实战中怎么选序列预测、纠错码与密码学里的真实用途4.1 序列预测一个完整的实例假设我拿到了一个自称随机的序列前几项是2, 4, 6, 8, 10, 12。直觉上这就是等差数列但我不想用肉眼判断直接用 BMseq [2, 4, 6, 8, 10, 12] C, L berlekamp_massey(seq) print(C, L) # 输出 [Fraction(1, 1), Fraction(-2, 1), Fraction(1, 1)] 2C [1, -2, 1]意味着s[n] - 2*s[n-1] s[n-2] 0也就是二阶差分全零。用这个递推预测下一个值完全没问题。如果某个序列混入了一个错误项比如把12改成13BM 给的递推阶数会立刻变大因为一个非均匀的错误破坏了低阶结构。这一点在异常检测里非常实用短递推阶数的突然上升往往意味着数据里混入了异常。4.2 纠错码与密码学中的 BMBM 最经典的应用场景是 BCH 码和 Reed-Solomon 码的译码。这类码的译码过程会计算伴随式序列而这些伴随式恰好满足一个由错误位置多项式决定的线性递推关系。用 BM 求出最短递推式本质上就是求出了错误位置多项式后续只需找根就能定位错误、再做纠错。这也是为什么 BM 在编码理论里地位这么高——它把找错误位置这个看似困难的问题转化成了求一个序列的最短递推这个精确可解的代数问题。在密码学里BM 是衡量序列线性复杂度的标准工具。给定一个二元序列线性复杂度就是能够生成它的最短 LFSR 的级数。如果这个数字太小说明序列可以用一个很短的线性模型完全复现那它在流密码里就不能用。反过来设计者也需要通过 BM 验证自己设计的序列线性复杂度是否足够高这是流密码安全性分析的基本功。4.3 BM 还是 Padé 重建选型对比需求推荐方案原因已知序列求最短递推阶数并预测后续项BM只关心分母次数实现简单复杂度是 O(N^2)已知幂级数恢复分子分母都有次数约束的有理函数扩展欧几里得 / Padé 逼近可以同时控制分子和分母的次数界两者似乎都行看下游需求需要分子信息用 Padé 更直接只做递推预测用 BM 更快我个人的习惯是如果是纯序列预测或者递推阶数分析无脑 BM如果是要还原一个有理函数本身比如从牛顿插值、有理插值或者某些密码协议中间量中恢复分式那就用扩展欧几里得版本。它们的数学底层是相通的但在工程接口上各有方便之处。5. 必须避开的坑从数值精度到索引错位的调试记录5.1 浮点数是最大的陷阱BM 算法对数值误差极度敏感。我在早期版本里图省事用float算Fibonacci 短序列前几步看起来没问题但误差会随着递推逐步放大到了第 20 项左右d不再是精确的 0而是一个1e-10级别的残差算法会认为递推式失配强行拉长递推阶数最后输出一个完全错误的连接多项式。更隐蔽的是如果你的序列本身来自浮点计算哪怕只有1e-12的噪声BM 都会把这个噪声解释成新的结构。所以我的原则是能用整数就绝不用浮点能用 Fraction 就绝不用小数。如果输入数据确实来自浮点计算要么先做有理化处理要么换用针对近似序列设计的稳定算法而不是硬上精确 BM。5.2 假递推与过拟合必须留出验证集BM 有个很迷惑人的性质给定任意一段长度有限的序列它总能给出一个递推式把这 N 项全部解释清楚。换句话说如果你只给 6 项算法一定能找出一个长度不超过 3 的递推式把 6 项完美拟合——但这不代表第 7 项能预测对。我踩过的坑是这样的给了一段平方数的前四项1, 4, 9, 16BM 很快就返回了一个长度 2 的递推式看着很美预测第五项却完全错了。原因很简单平方数序列真正满足的是三阶递推它的生成函数分母含(1-x)^3四项数据不足以识别出三阶结构算法只能拿一个更高阶的短模型硬拟合。从那以后我每次做序列预测都会留出最后若干项做交叉验证用前 N 项求递推预测后面几项如果预测值和真实值不一致这个递推式就是不可信的。提示验证时至少要留出 L 个连续项因为短递推式可能在很短窗口内碰巧正确连续多步验证才能暴露问题。5.3 递推式不唯一与 Padé 表中的方块现象在一般的整环上BM 给出的最短递推式不一定唯一。特别是模合数域或者特征 2 的有限域上可能同时存在多个长度相同的递推式它们都能解释已知数据。BM 返回的是算法内部选择的那一个不一定是你期望的那一个。这种情况下最稳妥的做法是增加观测项数数据越多不唯一性越容易被消除。有理函数重建也有类似问题在 Padé 逼近里叫方块现象。当分子分母次数约束周围存在多个候选时扩展欧几里得可能会停在不同的余式上。我会用自检条件 次数约束双重过滤但真正干净的办法还是加长观测窗口让有理函数的结构充分暴露出来。5.4 代码实现里的高频 Bug索引错位和边界判断写过 BM 的人大概率都在索引上翻过车。最典型的是修正公式里的x^m * B的位移方向B的第j个系数对应x^j乘上x^m之后应该存到新数组的jm位置。我第一次写的时候顺手写成了C[j-m]结果递推式完全错乱而且前几步因为下标越界还不会立刻报错只在某个边界上静默爆炸。另一个容易错的是多项式除法里的停止条件。如果poly_deg函数没有正确去掉尾部零系数除法循环会一直跑不停或者余式次数判断错位。我的建议是给每个 poly 工具函数写独立的单元测试用简单的小多项式先验证一遍再进主算法。排查 BM 问题时把每一步的C、L、m、d全部打印出来和手动演算表逐行对一遍比盯代码盯半小时有效得多。最后再分享一个我自己的习惯写完 BM 或者 Padé 重建第一件事不是看返回的C和p/q长得好不好看而是先往后多预测几项或者做一次q*f-p的零校验。只有自检通过我才会把这个递推式交给下游的代码。这种验证习惯救过我很多次因为算法在数学上正确和实现正确是两回事而数据本身有没有问题唯有验证能告诉你。