https://liam.page/https://liam.page/https://liam.page/

现代编程语言，大都在标准库中包含了随机库。例如，C++ 在 C++11 标准中添加了 random 头文件，提供了现代的随机库；Python 则有 random。C++11 的随机库将生成随机数的过程在逻辑上切分成了两个步骤：随机数生成引擎和分布。在学习 C++11 的 random 库时，std::mt19937 这一随机数生成引擎的名字看起来十分奇怪，成功吸引了我的注意力。

查询后得知，std::mt19937 中的 MT 是 Mersenne Twister 的缩写，这是伪随机数生成算法的名字(梅森旋转算法)；而 19937 则取自算法中用到的梅森素数 $2^{19937} - 1$。这里，梅森素数是算法生成伪随机数的循环长度(period)，而旋转则说的是算法内部对定长二进制串循环位移的过程。

此篇讲解梅森旋转算法的一些原理，并介绍对其的一个「爆破」方法。

伪随机数发生器质量的度量——$k$-维 $v$-比特准确度

梅森旋转算法(Mersenne Twister Algorithm，简称 MT)是为了解决过去伪随机数发生器(Pseudo-Random Number Generator，简称 PRNG)产生的伪随机数质量不高而提出的新算法。该算法由松本眞(Makoto Matsumoto)和西村拓士(Takuji Nishimura)在 1997 年提出，期间还得到了「算法之神」高德纳(Donald Ervin Knuth)的帮助。

既然 MT 是为了解决过去 PRNG 质量低下的问题，那么首先我们就必须要有一个能够度量 PRNG 质量的方法。否则，「公说公有理婆说婆有理」，我们就无法对 PRNG 作出统一的评价了。

这里介绍评价 PRNG 最严格指标[Tootill et al. 1973][Fushimi and Tezuka 1983][Couture et al. 1993][Tezuka 1995][Tezuka 1994][Tezuka and L'Ecuyer 1991][L'Ecuyer 1996]：$k$-维 $v$-比特准确度($k$-distributed to $v$-bit accuracy)。

假设某一 PRNG 能够产生周期为 $P$ 的 $w$-比特的随机数序列 ${vec x_i}$；同时，将 $w$-比特的随机数 $vec x$ 的最高 $v$ 位组成的数记作 $text{trunc}_v(vec x)$。现构造如下长度为 $kv$-比特的二进制数

$$text{PRNG}_{k, v}(i) overset{text{def}}{=}(text{trunc}_v(vec x_i),text{trunc}_v(vec x_{i + 1}),ldots,text{trunc}_v(vec x_{i + k - 1})).$$

由于 $text{PRNG}_{k, v}(i)$ 是长度为 $kv$-比特的二进制数，所以它可以有 $2^{kv}$ 中不同的取值。若当 $i$ 遍历 $[0, P)$，$text{PRNG}_{k, v}(i)$ 的取值在这 $2^{kv}$ 中均匀分布。具体来说，$text{PRNG}_{k, v}(i)$ 落在这 $2^{kv}$ 个取值上的数量完全相等($frac{P + 1}{2^{kv}}$)，除了全零的取值会少不多不少正好一次($frac{P + 1}{2^{kv}} - 1$)。

显而易见，对任意固定的 $v$，若 PRNG 是 $k$-维 $v$-比特准确的，那么必然也是 $(k - 1)$-维 $v$-比特准确的，但不一定是 $(k + 1)$-维 $v$-比特准确的。考虑 $k$ 是有上限的，因此，对于任意固定的 $v$，必然存在最大的 $k = k(v)$ 使得 PRNG 是 $k(v)$-维 $v$-比特准确的。那么，根据定义，显然有 $2^{k(v)v} - 1 leqslant P$。

$k$-维 $v$-比特准确度也可以有密码学角度的描述：若 PRNG 是 $k$-维 $v$-比特准确的，那么即使已知 PRNG 生成的 $l < k$ 个伪随机数的最高 $v$ 位，也无法推出 PRNG 生成的第 $l + 1$ 个伪随机数的最高 $v$ 位。

根据这样的定义，MT 算法具有非常优良的性能。首先 MT 算法(MT19937)的周期非常长。其周期 $P = 2^{19937} - 1$，要比预估的宇宙可观测的粒子总数($10^{87}$)还要高出数千个数量级。其次，作为一个 $32$ 比特随机数生成器，MT19937 是 $623$-维 $32$-比特准确的。考虑到 $bigllfloor frac{19937}{32}bigrrfloor = 623$，MT19937 在 $k$-维 $v$-比特准确度上的性能已达到理论上的最大值。因此，MT19937 具有非常优良的性能。

梅森旋转算法的描述

旋转

32 位的梅森旋转算法能够产生周期为 $P$ 的 $w$-比特的随机数序列 ${vec x_i}$；其中 $w = 32$。这也就是说，每一个 $vec x$ 是一个长度为 32 的行向量，并且其中的每一个元素都是二元数域 $mathbb{F}_2 overset{text{def}}{=} {0, 1}$ 中的元素。现在，我们定义如下一些记号，来描述梅森旋转算法是如何进行旋转(线性移位)的。

$n$：参与梅森旋转的随机数个数；

$r$：$[0, w)$ 之间的整数；

$m$：$(0, n]$ 之间的整数；

$mathbf{A}$：$w times w$ 的常矩阵；

$vec x^{(u)}$：$vec x$ 的最高 $w - r$ 比特组成的数(低位补零)；

$vec x^{(l)}$：$vec x$ 的最低 $r$ 比特组成的数(高位补零)。

梅森旋转算法，首先需要根据随机数种子初始化 $n$ 个行向量：

$$vec x_0, vec x_1, ldots, vec x_{n - 1}.$$

而后根据下式，从 $k = 0$ 开始依次计算 $vec x_{n}$：

begin{equation}vec x_{k + n} overset{text{def}}{=} vec x_{k + m}oplus bigl(vec x_{k}^{(u)}mid vec x_{k + 1}^{(l)}bigr)mathbf{A}.label{eq:twister}end{equation}

其中，$vec xmid vec x'$ 表示两个二进制数按位或；$vec xoplus vec x'$ 表示两个二进制数按位半加(不进位，也就是按位异或)；$vec xmathbf A$ 则表示按位半加的矩阵乘法。在 MT 中，$mathbf A$ 被定义为

$$begin{pmatrix}

& 1 \

& & 1 \

& & & ddots \

& & & & 1 \

a_{w - 1} & a_{w - 2} & a_{w - 3} & cdots & a_0

end{pmatrix}$$

因此

$$vec xmathbf A = begin{cases}vec x >> 1& text{if $x_0 = 0$} \ (vec x >> 1)oplusvec a& text{if $x_0 = 1$}end{cases}.$$

此处，若 $r = 0$，则 MT 退化为 TGFSR (Matsumoto and Kurita 1992, 1994)；若再加上 $mathbf A = mathbf I_{w}$，则又从 TGFSR 退化为 GFSR (Lewis and Payne 1973)。

因此，梅森旋转 ref{eq:twister} 式完全由位运算组成(移位、按位与、按位或、按位异或)。

线性反馈移位寄存器、旋转之名、周期

上一小节我们介绍了 MT 算法当中的「旋转」。但只凭抽象的数学公式(尤其是二进制的逻辑数学)，很难看出它为什么是「旋转」。这一节我们首先介绍线性反馈移位寄存器(Linear Feedback Shifting Register，简称 LFSR)，看到它是如何「旋转」的；最后再将 LFSR 和 MT 算法当中的旋转步骤统一起来。

反馈移位寄存器是对二进制序列进行等长变换的一种特殊函数。它包括两个部分：

级。等长変换的长度即是反馈移位寄存器的级。

反馈函数。若反馈函数是线性的，则称线性反馈移位寄存器；否则是非线性反馈移位寄存器。

一般来说，LFSR 是基于异或运算的。一个 LFSR 的工作步骤是这样的：

将原寄存器状态的最低位作为输出。

执行线性反馈函数；也就是选取其中若干位，从高位到低位迭代异或。

将元寄存器状态向低位移位 1 位，并以上述迭代异或的结果作为填充值填入最高位。

在不断的迭代异或、填充高位的过程中，二进制位在寄存器中循环旋转。这就是「旋转」的来由。

对于一个 4 级的 LFSR 来说，假设其反馈函数是 $f(x) = x^4 + x^2 + x + 1$。则 LFSR 每次从最低位取出结果，将最高位($x^4$)和倒数第二低位($x^2$)取异或后，再与最低位($x$)取异或后，填入移位后的最高位。

因此，若初始状态为 $(1000)_2$，则有

1

2

3

4

5

6

71000

1100

1110

0111

0011

0001

1000 # 复原

如此，我们就构建了一个循环长度为 6 的 LFSR。

考虑一个 $w$ 级的 LFSR，其最多共有 $2^w$ 种状态。又考虑对于异或来说，全零的状态是不可用的(因为不论如何运算都是全零)，因此一个 $w$ 级的 LFSR，其有效状态共有 $2^w - 1$ 个。因此，理论上，一个 LFSR 的循环长度最大为 $2^w - 1$。可以证明，当且仅当反馈函数是 $mathbb F_2$ 上的本原多项式(素多项式)时，LFSR 的循环长度达到最大值。

回过头来看 ref{eq:twister} 式，不难发现，这其实相当于一个 $nw - r$ 级的线性反馈移位寄存器(取 $vec x_k^{(u)}$ 的最高 $w - r$ 位与 $vec x_{k + 1}^{(l)}$ 的最低 $r$ 位进行迭代异或，再经过一个不影响周期的线性变换 $mathbf A$)。只不过，ref{eq:twister} 式每一次运算，相当于 LFSR 进行了 $w$ 轮计算。若 $w$ 与 $nw - r$ 互素，那么这一微小的改变是不会影响 LFSR 的周期的。考虑到 LFSR 的计算过程像是在「旋转」，这即是「梅森『旋转』」名字的来由。

对这个等效的 $nw - r$ 级 LFSR 来说，当且仅当反馈函数是 $mathbb F_2$ 上的本原多项式(素多项式)时，MT 的循环周期长度 $P$ 达到最大值($2^{nw - r} - 1$)。

提取(tempering)输出

MT19937 有两个主要特性。一是周期很长，达到 $2^{19937} - 1$，二是满足 $623$-维 $32$-比特准确性。上述「旋转」的过程，帮助我们达成了较长的周期。接下来，我们需要将每次旋转的结果提取(tempering)输出，以保证 MT 是 $623$-维 $32$-比特准确的。

提取的方法很简单，只需要将每次旋转得到的输出右乘一个可逆矩阵 $mathbf T$ 即可。将 $vec xmapsto vec xmathbf T$ 表述成位运算，则有

begin{align}

vec y &{}gets vec xoplus (vec x >> u) \

vec y &{}gets vec yoplus ((vec y << s) mathop{mathbf{AND}} vec b) \

vec y &{}gets vec yoplus ((vec y << t) mathop{mathbf{AND}} vec c) \

vec z &{}gets vec yoplus (vec y >> l)

end{align}

此处，$u$, $s$, $t$, $l$ 是整数参数；$vec b$ 和 $vec c$ 是 $w$-比特的整数，用作比特遮罩(bit mask)；最终能得到的 $vec z$ 即是当前轮次的随机数输出。

算法描述

这样一来，MT 算法的主要有两个部分：旋转、提取。

旋转部分有参数

$w$：生成的随机数的二进制长度；

$n$：参与旋转的随机数个数(旋转的深度)；

$m$：参与旋转的中间项；

$r$：$vec x_{k}^{(u)}$ 和 $vec x_{k + 1}^{(l)}$ 的切分位置；

$vec a$：矩阵 $mathbf A$ 的最后一行。

提取部分有参数

$u$, $s$, $t$, $l$：整数参数，移位运算的移动距离；

$vec b$, $vec c$：比特遮罩。

于是，我们得到 MT 算法的完整描述如下。

常量初始化。

begin{align*}

vec{lower} &{}gets 0 \

vec{lower} &{}gets (vec{lower} << 1) mathop{mathrm{AND}} 1quad text{for $r$ times.} \

vec{upper} &{}gets mathop{mathrm{COMPL}}vec{lower} \

vec a &{}gets a_{w - 1}a_{w - 2}cdots a_{1}a_{0} \

i &{}gets 0

end{align*}

工作区非零初始化。

$$vec x[0], vec x[1], ldots, vec x[n - 1]$$

旋转

begin{align*}

vec t &{}gets (vec x[i]mathop{mathrm{AND}} vec{upper}) mathop{mathrm{OR}} (vec x[(i + 1) mod n]mathop{mathrm{AND}} vec{lower}) \

vec x[i] &{}gets vec x[(i + m) mod n]mathop{mathrm{XOR}} (vec t >> 1) mathop{mathrm{XOR}} begin{cases}vec 0& text{if $t_0 = 0$}\ vec a& text{otherwise}end{cases}

end{align*}

提取输出。

begin{align*}

vec y &{}gets vec x[i] \

vec y &{}gets vec y mathop{mathbf{XOR}} (vec y >> u) \

vec y &{}gets vec y mathop{mathbf{XOR}} ((vec y << s) mathop{mathbf{AND}} vec b) \

vec y &{}gets vec y mathop{mathbf{XOR}} ((vec y << t) mathop{mathbf{AND}} vec c) \

vec y &{}gets vec y mathop{mathbf{XOR}} (vec y >> l) \

&{} text{output $vec y$}

end{align*}

更新循环变量 $igets (i + 1)mod n$，而后跳转至步骤 3 继续进行。

再探梅森旋转

至此，我们已经探索了梅森旋转算法表面上的全部内容：我们已经知道梅森旋转算法是什么，也知道梅森旋转算法为什么这样起名，也有了完整的算法描述。但是，关于梅森旋转算法还有很多深层的问题我们未曾探索。比如说，对于 $n$, $w$ 和 $r$ 的组合，我们是否有必要追求最长周期 $P$ 使得 $P = w^{nw - r} - 1$？又比如说，我们提到 LFSR 取得最长周期的充要条件是反馈函数是 $mathbb F_{2}$ 上的素多项式，那么怎样验证反馈函数是否是素的？

这一节，我们来讨论这些问题。

关于周期

前文提到，梅森旋转的过程，实际上是对长度为 $nw - r$ 的二进制串做了一个 LFSR 的变体。这里，我们将它记作 $mathbf B$。

我们已经知道，这个 LFSR 的变体，其周期的上限是 $2^{nw - r} - 1$。这样一来，整个序列的周期达到这一上限就意味着除去全零的状态，整个序列每一种可能状态都被遍历了一遍；而全零的状态则被遍历了 1 遍。考虑在这 $nw - r$ 比特的序列中，$\{vec x_n\}$ 有 $n - 1$ 个完整的 $w$-比特向量；因此，$\{vec x_n\}$ 显然是 $(n - 1)$-维的。这也就是说，选择不同的随机数种子，至多只能改变 ${vec x_n}$ 序列的起始相位。

这样一来，我们有：当梅森旋转达到最大周期时，若 $n$ 确定，$n - 1$ 就确定了，进而整个序列同分布的维数 $n - 1$ 也就确定了。因此，对于梅森旋转而言，提升维数是很容易的事情。

这即是努力使梅森旋转达到最大周期的意义。

多项式素检测与参数调优

在梅森旋转算法中，反馈函数($mathbf B$ 的特征多项式)的素检测是很容易的。这是因为，对于 $p = nw - r$(其中 $p$ 是梅森素数的幂)级的 LFSR 来说，其反馈函数在 $mathbb F_2$ 上的素检测的复杂度是 $O(p^2)$。这一方面得益于梅森素数的性质，另一方面得益于 MT 是工作在 $mathbb F_2$ 上的算法。(Matsumoto and Nishimura, 1997) 这一特性的证明，牵扯到很多抽象代数和数论方面的知识；此处我们按下不表，留待后续用专门的文章来证明。

梅森旋转算法中，要实现 PRNG 的最佳性能，需要对旋转和提取两部分参数做细致的调整。调整这部分参数，寻得最优参数组合，是有特定算法可寻的。这部分内容十分繁琐，此处也不表。有兴趣的用户可阅读梅森旋转算法原始论文第四节、第五节。

梅森旋转算法的 Python 实现

此处给出一个 Python 实现的梅森旋转算法(mt19937)，为后续对算法的「爆破」提供素材。

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

56

57

58

59

60

61

62class :

__n = 624

__m = 397

__a = 0x9908b0df

__b = 0x9d2c5680

__c = 0xefc60000

__kInitOperand = 0x6c078965

__kMaxBits = 0xffffffff

__kUpperBits = 0x80000000

__kLowerBits = 0x7fffffff

def __init__(self, seed = 0):

self.__register = [0] * self.__n

self.__state = 0

self.__register[0] = seed

for i in range(1, self.__n):

prev = self.__register[i - 1]

temp = self.__kInitOperand * (prev ^ (prev >> 30)) + i

self.__register[i] = temp & self.__kMaxBits

def __twister(self):

for i in range(self.__n):

y = (self.__register[i] & self.__kUpperBits) +

(self.__register[(i + 1) % self.__n] & self.__kLowerBits)

self.__register[i] = self.__register[(i + self.__m) % self.__n] ^ (y >> 1)

if y % 2:

self.__register[i] ^= self.__a

return None

def __temper(self):

if self.__state == 0:

self.__twister()

y = self.__register[self.__state]

y = y ^ (y >> 11)

y = y ^ (y << 7) & self.__b

y = y ^ (y << 15) & self.__c

y = y ^ (y >> 18)

self.__state = (self.__state + 1) % self.__n

return y

def __call__(self):

return self.__temper()

def load_register(self, register):

self.__state = 0

self.__register = register

if __name__ == "__main__":

mt = MersenneTwister(0)

tank = set()

kLen = 100

for i in range(kLen):

t = mt()

tank.add(t)

print(t)

print(len(tank) == kLen)

爆破梅森旋转算法

梅森旋转算法的设计目的是优秀的伪随机数发生算法，而不是产生密码学上安全的随机数。从梅森旋转算法的结构上说，其提取算法 __temper 完全基于二进制的按位异或；而二进制按位异或是可逆的，故而 __temper 是可逆的。这就意味着，攻击者可以从梅森旋转算法的输出，逆推出产生该输出的内部寄存器状态 __register[__state]。若攻击者能够获得连续的至少 __n 个寄存器状态，那么攻击者就能预测出接下来的随机数序列。

现在我们遵循这个思路，爆破梅森旋转算法。

逆向 __temper

我们以向右移位后异或为例，首先观察原函数。

1

2

3

4def right_shift_xor(value, shift):

result = value

result ^= (result >> shift)

return result

简单起见，我们观察一个 8 位二进制数，右移 3 位后异或的过程。

1

2

3value: 1101 0010

shifted: 0001 1010 # 010 (>> 3)

result: 1100 1000

首先，观察到 result 的最高 shift 位与 value 的最高 shift 位是一样的。因此，在 result 的基础上，我们可以将其与一个二进制遮罩取与，得到 value 的最高 shift 位。这个遮罩应该是：1111 1111 << (8 - 3) = 1110 0000。于是我们得到 1100 0000。

其次，注意到对于异或运算有如下事实：a ^ b ^ b = a。依靠二进制遮罩，我们已经获得了 value 的最高 shift 位。因此，我们也就能得到 shifted 的最高 2 * shift 位。它应该是 1100 0000 >> 3 = 0001 1000。将其与 result 取异或，则能得到 value 的最高 2 * shift 位。于是我们得到 1101 0000。

如此往复，即可复原 value。据此有代码

1

2

3

4

5

6

7

8

9def inverse_right_shift_xor(value, shift):

i, result = 0, 0

while i * shift < 32:

part_mask = ((0xffffffff << (32 - shift)) & 0xffffffff) >> (i * shift)

part = value & part_mask

value ^= part >> shift

result |= part

i += 1

return result

对左移后取异或，也有类似分析。于是，得到对 __temper 的完整求逆代码。

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

35class TemperInverser:

__b = 0x9d2c5680

__c = 0xefc60000

__kMaxBits = 0xffffffff

def __inverse_right_shift_xor(self, value, shift):

i, result = 0, 0

while i * shift < 32:

part_mask = ((self.__kMaxBits << (32 - shift)) & self.__kMaxBits) >> (i * shift)

part = value & part_mask

value ^= part >> shift

result |= part

i += 1

return result

def __inverse_left_shift_xor(self, value, shift, mask):

i, result = 0, 0

while i * shift < 32:

part_mask = (self.__kMaxBits >> (32 - shift)) << (i * shift)

part = value & part_mask

value ^= (part << shift) & mask

result |= part

i += 1

return result

def __inverse_temper(self, tempered):

value = tempered

value = self.__inverse_right_shift_xor(value, 18)

value = self.__inverse_left_shift_xor(value, 15, self.__c)

value = self.__inverse_left_shift_xor(value, 7, self.__b)

value = self.__inverse_right_shift_xor(value, 11)

return value

def __call__(self, tempered):

return self.__inverse_temper(tempered)

爆破

逆向 __temper() 之后，只要获得足够的状态，即可构建出梅森旋转内部的寄存器状态。因此有如下验证代码。

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19class MersenneTwisterCracker:

__n = 624

def __init__(self, mt_obj):

inverser = TemperInverser()

register = [inverser(mt_obj()) for i in range(self.__n)]

self.__mt = MersenneTwister(0)

self.__mt.load_register(register)

def __call__(self):

return self.__mt()

if __name__ == "__main__":

mt = MersenneTwister(0)

for i in range(100):

mt()

mtc = MersenneTwisterCracker(mt)

for i in range(100):

assert(mt() == mtc())

运行后，Python 没有抛出异常，顺利推出。这说明 mtc 已能够成功预测 mt 之后的任意顺序输出。

总结

梅森旋转算法，是一个优秀的伪随机数发生算法。在伪随机数的评价体系中，它是一个相当优秀的算法：周期长、均匀性好、速度快(基本都是位运算)。在条件允可的情形下，若有使用随机数的需求，应首先考虑梅森旋转算法。

同时也应该注意到，梅森旋转算法不是为了密码学随机而设计的——在获得足够连续输出的情况下，梅森旋转算法接下来的输出值是可以准确预测的。梅森旋转算法容易被爆破的根源在于，其提取输出函数是可逆的，因此暴露了其内部状态。若要产生密码学上的随机数，可考虑在梅森旋转算法之后，拼接一值得信赖的单向杂凑函数(如 sha256)。否则，若直接用梅森旋转算法的输出值作密码学用途，则有信息泄露的风险，应引起注意。

错误应用梅森旋转算法，导致高危漏洞的一个典型是 Discuz! 的密码重置漏洞。