cdcq的密码学教程四——BM算法

发表于2022-08-28更新于2023-04-30阅读次数

一、前言

BM 算法是一个神奇的算法，它不仅是序列密码的基础算法，还能在 ACM 比赛中起到出人意料的效果。它的做法不复杂，正确性也容易理解，但你真的仔细考虑过它的最优性吗？这篇文章就将带你理解 BM 算法的原理，并将证明它的正确性和最优性。

文章中大多数内容都来自 Massey （也就是 BM 中的 M）的原文《Shift-Register Synthesis and BCH Decoding》，原文中对于算法的描述和证明已经很详细了，我这篇文章也仅仅是在不好懂的地方做了一些补充。原文除介绍 BM 算法本身，还给出了有关 LFSR 的一些更强性质，如解的唯一性等，以及 BM 算法在解码 BCH 码上的应用，还是推荐有能力的读者去看原文。虽然原文是英文，但是看清楚每一个细节的话，还是很好懂的。

另外是我在每篇数学内容比较多的文章里都会强调的建议：文章里会出现较多的数学推导，它们看上去很难，但往往只是形状复杂，而读者对符号不甚熟悉。只要能跟着文章自己动手写一遍，即使是普通高中生也能理解推导的过程。因此推荐大家一定要动手写，实在懒得动手或不方便的话也要认真看清每个符号，否则是看不懂的。

二、常系数齐次线性递推式

形如

sj=−∑i=1Lcisj−i,j=L,L+1,L+2,⋯s_j=-\sum\limits_{i=1}^Lc_is_{j-i}, \quad j=L, L+1, L+2, \cdots

的递推式被称为常系数齐次线性递推式，其中 L 为此递推式的长度。这个公式生成的序列中从第 L 项开始，每一项都由之前的 L 项加权求和得到。只要定义了前 L 项 s0,s1,⋯ ,sL−1s_0, s_1, \cdots, s_{L-1} 和系数 c1,c2,⋯ ,cLc_1, c_2, \cdots, c_L ，就可以根据此式不断生成其余的项。这个公式还有另一个常用的等式形式

sj+∑i=1Lcisj−i=0,j=L,L+1,L+2,⋯ .s_j+\sum\limits_{i=1}^Lc_is_{j-i}=0, \quad j=L, L+1, L+2, \cdots.

这个形式看上去与递推形式没有大的差别，但在后文中可以看到，等式形式在 BM 算法的分析中是很方便的。

如果将这个公式应用在 GF(2)GF(2)（即模 2 域）下，就被称为是线性反馈移位寄存器（LFSR）。为了行文方便，下文都以 LFSR 代替常系数齐次线性递推式进行分析，得到的结果可以很轻易地应用到实数域、模质数域等其他域上。

BM 算法的主要作用就是根据一个比较长的序列，推测出产生此序列的 LFSR 。

三、需要用到的定理

要想理解 BM 算法，需要先推导出两个小小的结论。

定理 1

如果一个长度为 L 的 LFSR 可以生成序列 s0,s1,⋯ ,sN−1s_0, s_1, \cdots, s_{N-1} 但不能生成 s0,s1,⋯ ,sNs_0, s_1, \cdots, s_N ，那么任意可以生成第二个序列的 LFSR 的长度 L′L' 都满足

L′≥N+1−L.L'\geq N+1-L.

证明：对于 L≥NL\geq N 的情况，答案是显然的，所以只考虑 L<NL<N 的情况。令 c1,c2,⋯ ,cLc_1, c_2, \cdots, c_L 和 c1′,c2′,⋯ ,cL′′c'_1, c'_2, \cdots, c'_{L'} 分别表示只能生成定理中第一个序列和可以生成全部两个序列的 LFSR 的反馈系数。即

−∑i=1Lcisj−i=sj,j=L,L+1,⋯ ,N−1≠sj,j=N,-\sum\limits_{i=1}^Lc_is_{j-i} \quad \begin{aligned} =s_j, \quad j&=L, L+1, \cdots, N-1 \\ \neq s_j, \quad j&=N, \end{aligned}

以及

−∑i=1L′c′isj−i=sj,j=L′,L′+1,⋯ ,N.-\sum\limits_{i=1}^{L'}{c'}_is_{j-i}=s_j, \quad j=L', L'+1, \cdots, N.

用反证法，假设 L′≤N−LL'\leq N-L ，那么当 j=Nj=N 时 −∑i=1LcisN−i-\sum\limits_{i=1}^Lc_is_{N-i} 中的 sN−is_{N-i} 均可以用第二个 LFSR 展开，得到

−∑i=1LcisN−i=∑i=1Lci∑k=1L′c′ksN−i−k=−∑k=1L′c′k∑i=1LcisN−k−i=−∑k=1L′c′ksN−k=sN.\begin{aligned} -\sum\limits_{i=1}^Lc_is_{N-i}&=\sum\limits_{i=1}^{L}{c}_i\sum\limits_{k=1}^{L'}{c'}_ks_{N-i-k} \\ &=-\sum\limits_{k=1}^{L'}{c'}_k\sum\limits_{i=1}^Lc_is_{N-k-i} \\ &=-\sum\limits_{k=1}^{L'}{c'}_ks_{N-k} \\ &=s_N. \end{aligned}

这显然与第一个 LFSR 的定义不符，所以假设失败，必有 L′≥N+1−LL'\geq N+1-L 。

现在定义 s\mathbf{s} 为无限序列 s0,s1,s2,⋯s_0, s_1, s_2, \cdots ，LN(s)L_N(\mathbf{s}) 为生成 s\mathbf{s} 的前 N 项 s0,s1,⋯ ,sN−1s_0, s_1, \cdots, s_{N-1} 所需的最小 LFSR 长度。特别地，如果前 N 项均为 0 ，那么定义 LN(s)=0L_N(\mathbf{s})=0 。关于 LN(s)L_N(\mathbf{s}) ，有如下结论：

引理 1

如果一个长度为 LN(s)L_N(\mathbf{s}) 的 LFSR 可以生成 s\mathbf{s} 的前 N 项，而不能生成前 N+1N+1 项，那么有

LN+1(s)≥max[LN(s),N+1−LN(s)].L_{N+1}(\mathbf{s})\geq max[L_N(\mathbf{s}), N+1-L_N(\mathbf{s})].

证明：首先 LN(s)L_N(\mathbf{s}) 一定是递增的，因为对于任意能生成前 N+1N+1 项的 LFSR ，显然也可以生成前 N 项，只需要将第 N+1N+1 项忽略就可以了。而可能存在能生成前 N 项却不能生成前 N+1N+1 项的 LFSR ，因此有 LN+1(s)≥LN(s)L_{N+1}(\mathbf{s})\geq L_N(\mathbf{s}) 。其次根据定理 1 可以得到 LN+1(s)≥N+1−LN(s)L_{N+1}(\mathbf{s})\geq N+1-L_N(\mathbf{s}) ，因此结论成立。

引理 1 看上去很简单，证明也很显然，但对接下来 BM 算法的推导有很大作用。后文也将证明，这里的大于等于号其实可以替换为等号。

四、BM 算法的推导

由于序列的前 N 项已经给定了，所以还原 LFSR 只需找出反馈系数。为了方便表示，需要定义一个多项式来表示 LFSR 的反馈系数：

C(D)=1+c1D+c2D2+⋯+cLDL.C(D)=1+c_1D+c_2D^2+\cdots+c_LD^L.

特别地，当 LFSR 的长度为 0 时，C(D)=1C(D)=1 。

需要指出的是，这个多项式和 LFSR 以及生成的序列都没有关系，它只是用来方便地表示 LFSR 的长度以及系数的位置而已。

对于给定的 s\mathbf{s} ，定义

C(N)(D)=1+c1(N)D+⋯+cLN(s)(N)DLN(s)C^{(N)}(D)=1+c_1^{(N)}D+\cdots+c_{L_N(\mathbf{s})}^{(N)}D^{L_N(\mathbf{s})}

表示一个能生成此序列前 N 项的 LFSR 。其中上标 (N)(N) 用括号表明这是一个标号，而不是指数。

BM 的求解是一个迭代的过程，现在假设已经得到了对于 N=1,2,⋯ ,nN=1, 2, \cdots, n 的所有 C(N)(D)C^{(N)}(D) ，下一步的目标是找出 Ln(s)L_n(\mathbf{s}) 以及一个任意的 C(n)(D)C^{(n)}(D) 。

令 dnd_n 为 C(n)(D)C^{(n)}(D) 在第 n+1n+1 项的误差，即

sj+∑i=1Ln(s)ci(n)sj−i={0,j=Ln(s),⋯ ,n−1dn,j=n,s_j+\sum\limits_{i=1}^{L_n(\mathbf{s})}c^{(n)}_is_{j-i}= \left\{ \begin{aligned} 0, \quad j&=L_n(\mathbf{s}), \cdots, n-1 \\ d_n, \quad j&=n, \end{aligned} \right.

如果 dn=0d_n=0 ，那么 C(n)(D)C^{(n)}(D) 也可以生成第 n+1n+1 项，因此不必做出变化。接下来考虑 dn≠0d_n\neq 0 的情况，尝试消除误差。

假设 m 是最近一次可以引起最短 LFSR 长度变短的位置，即

Lm(s)<Ln(s)Lm+1(s)=Ln(s).\begin{aligned} L_m(\mathbf{s})<L_n(\mathbf{s}) \\ L_{m+1}(\mathbf{s})=L_n(\mathbf{s}). \end{aligned}

因为从第 m 项到第 m+1m+1 项时 LFSR 的长度增加了，因此一定有

sj+∑i=1Lm(s)ci(m)sj−i={0,j=Lm(s),⋯ ,m−1dm≠0,j=m,s_j+\sum\limits_{i=1}^{L_m(\mathbf{s})}c^{(m)}_is_{j-i}= \left\{ \begin{aligned} 0, \quad j&=L_m(\mathbf{s}), \cdots, m-1 \\ d_m\neq 0, \quad j&=m, \end{aligned} \right.

可以利用此处的 dmd_m 来消除 dnd_n 的误差。

定义

C′(D)=C(n)(D)−dndm−1Dn−mC(m)(D).C'(D)=C^{(n)}(D)-d_nd_m^{-1}D^{n-m}C^{(m)}(D).

式中 Dn−mD^{n-m} 一项的含义是将 C(m)(D)C^{(m)}(D) 偏移到指数为 n−m,n−m+1,⋯ ,n−m+Lm(s)n-m, n-m+1, \cdots, n-m+L_m(\mathbf{s}) 的位置再减去，C′(D)C'(D) 的长度实际上就变成 L′=max[Ln(s),n−m+Lm(s)]L'=max[L_n(\mathbf{s}), n-m+L_m(\mathbf{s})] 。对于 L′L' 到 n−1n-1 的位置，由于同时满足

sj+∑i=1Lci(n)sj−i=0,j=L′,L′+1,⋯ ,n−1s_j+\sum\limits_{i=1}^Lc^{(n)}_is_{j-i}=0, \quad j=L', L'+1, \cdots, n-1

和

sj+∑i=1Lci(m)sj−i=0,j=L′,L′+1,⋯ ,n−1,s_j+\sum\limits_{i=1}^Lc^{(m)}_is_{j-i}=0, \quad j=L', L'+1, \cdots, n-1,

所以 C′(D)C'(D) 可以生成 L′L' 到 n−1n-1 的项，即前 n 项。而对于 j=nj=n 的情况，有

sn+∑i=1Ln(s)ci(n)sn−i−dndm−1[sn−(n−m)+∑i=1Lm(s)ci(m)sn−(n−m)−i]=sn+∑i=1Ln(s)ci(n)sn−i−dndm−1[sm+∑i=1Lm(s)ci(m)sm−i]=dn−dndm−1dm=0.\begin{aligned} &s_n+\sum\limits_{i=1}^{L_n(\mathbf{s})}c^{(n)}_is_{n-i}-d_nd^{-1}_m\left[s_{n-(n-m)}+\sum\limits_{i=1}^{L_m(\mathbf{s})}c^{(m)}_is_{n-(n-m)-i}\right] \\ &\quad =s_n+\sum\limits_{i=1}^{L_n(\mathbf{s})}c^{(n)}_is_{n-i}-d_nd^{-1}_m\left[s_m+\sum\limits_{i=1}^{L_m(\mathbf{s})}c^{(m)}_is_{m-i}\right] \\ &\quad =d_n-d_nd^{-1}_md_m \\ &\quad =0. \end{aligned}

因此可以说 C′(D)C'(D) 也可以生成第 n+1n+1 项。这样，我们就找到了一个可以生成前 n+1n+1 项的 LFSR 。

这里有一个容易迷惑的点，有的读者可能会认为从 n−mn-m 到 Lm(s)L_m(\mathbf{s}) 的项无法生成，因为这一段位置只对 C(m)(D)C^{(m)}(D) 的部分系数进行求和。产生这种误解的原因是混淆了 n 和 Ln(s)L_n(\mathbf{s}) 。事实上，n−m+Lm(s)n-m+L_m(\mathbf{s}) 不一定小于 Ln(s)L_n(\mathbf{s}) ，而 C′(D)C'(D) 的长度 L′L' 已经将 C(m)(D)C^{(m)}(D) 的所有包含在内了。小于 L′L' 的项都是作为初始项而不必生成的。

找到可行解之后，接下来就要证明 C′(D)C'(D) 是生成前 n 项的最短 LFSR 。根据引理 1 ，知道

L′=Ln+1(s)≥max[Ln(s),n+1−Ln(s)],L'=L_{n+1}(\mathbf{s})\geq max[L_n(\mathbf{s}), n+1-L_n(\mathbf{s})],

因此只要能证明 LN+1(s)=max[LN(s),N+1−LN(s)]L_{N+1}(\mathbf{s})=max[L_N(\mathbf{s}), N+1-L_N(\mathbf{s})] ，就可以说明 C′(D)C'(D) 是最短 LFSR 。

证明方法为归纳证明，先假设此条件对于 N=mN=m 成立，即 Lm+1(s)=max[Lm(s),m+1−Lm(s)]L_{m+1}(\mathbf{s})=max[L_m(\mathbf{s}), m+1-L_m(\mathbf{s})] 。根据 m 的定义，有

Ln(s)=Lm+1(s)=max[Lm(s),m+1−Lm(s)]Ln(s)>Lm(s).\begin{aligned} L_n(\mathbf{s})&=L_{m+1}(\mathbf{s})=max[L_m(\mathbf{s}), m+1-L_m(\mathbf{s})] \\ L_n(\mathbf{s})&>L_m(\mathbf{s}). \end{aligned}

所以

Ln(s)=m+1−Lm(s).L_n(\mathbf{s})=m+1-L_m(\mathbf{s}).

即

m=Ln(s)+Lm(s)−1.m=L_n(\mathbf{s})+L_m(\mathbf{s})-1.

所以

Ln+1(s)=max[Ln(s),n−m+Lm(s)]=max[Ln(s),n+1−Ln(s)].\begin{aligned} L_{n+1}(\mathbf{s})&=max[L_n(\mathbf{s}), n-m+L_m(\mathbf{s})] \\ &=max[L_n(\mathbf{s}), n+1-L_n(\mathbf{s})]. \end{aligned}

而对于初始情况，即 n=0n=0 时，显然有

L0(s)=0L1(s)=1,\begin{aligned} L_0(\mathbf{s})&=0 \\ L_1(\mathbf{s})&=1, \\ \end{aligned}

即

L1(s)=max[L0(s),0+1−L0(s)].L_1(\mathbf{s})=max[L_0(\mathbf{s}), 0+1-L_0(\mathbf{s})].

因此归纳成立，LN+1(s)=[LN(s),N+1−LN(s)]L_{N+1}(\mathbf{s})=[L_N(\mathbf{s}), N+1-L_N(\mathbf{s})] 。这样就证明了 C′(D)C'(D) 是最短的可以生成前 n+1n+1 项的 LFSR 。

五、BM 算法的应用

根据公式可以看到，BM 算法的全部计算都可以作用在任意域上，因此文章开头就说，BM 算法不止可以用来求解 LFSR ，还可以用来求解在模质数域、实数域等其他域上的常系数齐次线性递推式。

因为我自己的了解也有限，BM 算法在密码学和编码学的应用中就不多介绍了。这里只给出一个 BM 算法在 ACM 比赛中的神奇应用：

有一些题目其实是结论题，它需要选手用较高的数学水平或思维能力推出一个公式，然后利用此公式快速计算。虽然推导公式的过程是困难的，但是最后的结论一般是简单的。如果某道题的结论是一个常系数齐次线性递推式，那么就可以先用暴力算法跑出答案前几项，然后用 BM 算法得到递推式。我和队友在 CCPC2020 绵阳站就通过这种方式做出了一道题，最终拿了银牌。

最后给出 BM 算法的 c++ 实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
using ll=long long int;
ll mo=2;  // 模数可以换任意质数
ll qcp(ll x,ll y){  // 快速幂
	ll z=x;
	for(;y;y>>=1){
		if(y&1)  z=z*x%mo;
		x=x*x%mo;
	}
	return z;
}
namespace BM{
	vector<ll> BM(vector<ll> s){
		vector<ll> C(1,1),B(1,1);
		ll L=0,m=1,b=1;
		for(int n=0;n<(int)s.size();++n){
			ll d=0;
			for(int i=0;i<L+1;++i)
				d=(d+(ll)C[i]*s[n-i])%mo;
			if(d==0)  ++m;
			else if(2*L<=n){
				vector<ll> T=C;
				ll c=mo-d*qcp(b,mo-2)%mo;
				while(C.size()<B.size()+m)
					C.push_back(0);
				for(int i=0;i<(int)B.size();++i)
					C[i+m]=(C[i+m]+c*B[i])%mo;
				L=n+1-L;
				B=T;
				b=d;
				m=1;
			}
			else{
				ll c=mo-d*qcp(b,mo-2)%mo;
				while(C.size()<B.size()+m)
					C.push_back(0);
				for(int i=0;i<(int)B.size();++i)
					C[i+m]=(C[i+m]+c*B[i])%mo;
				++m;
			}
		}
		return C;
	}
	bool check(vector<ll> u,vector<ll> v){
        // 注意这里长度小于递推式长度的序列直接判否了
        // 有些场景可能不适用
		if(u.size()>v.size())  return false;
		for(int i=u.size()-1;i<(int)v.size();++i){
			ll d=0;
			for(int j=0;j<(int)u.size();++j)
				d=(d+u[j]*v[i-j])%mo;
			if(d!=0)  return false;
		}
		return true;
	}
}