Mersenne Twister Algorithm Analysis

发表于 更新于 分类于 本文字数： 5.8k 阅读时长 ≈ 10 分钟

概述

梅森旋转算法（Mersenne Twister Algorithm） 是一个伪随机发生算法。由松本真和西村拓士在 1997 年开发，基于有限二进制字段上的矩阵线性递归。虽然在实现上，MT算法使用了反馈移位寄存器（FSR），但注意其不能作为密码学意义上的伪随机数发生器。

最为广泛使用 Mersenne Twister 的一种算法变体是 MT19937，可以产生 32 位整数序列。

预备知识

梅森素数

如果 p 是素数，我们将 $ 2^p - 1 $ 的一类数记为 $ M_p $，如果 $ M_p $ 也是素数，那么它就被称为 梅森素数（Mersenne Prime Number）。

首先，用因式分解法可以证明，若 $ 2^p - 1 $ 是素数，则 p 也一定是素数，证明如下：
首先，易知 k方差的因式分解如下：
$$ a^k - 1 = (a - 1)(a^{k-1} + a^{k-2} + … + a + 1) $$

我们假设 p 不是素数，其因式分解为 $ p = sk (s > 1, k > 1)$，则：
$$ 2^p - 1 = 2^{sk} - 1 = ({2^s})^k - 1 = (2^s - 1)[({2^s})^{k-1} + ({2^s})^{k-2} + … + 2^s + 1] $$

即 $ 2^p - 1 $ 可以因式分解为大于 1 的整数，这与其为素数的前提条件矛盾，故原假设不成立，所以 p 是素数。

但反过来，如果 p 为素数，则 $ 2^p - 1 $ 不一定为素数。目前已经发现 52 个梅森素数，最大的是 $ M_{136279841} $（即 $ 2^{136279841} - 1 $）。

算法分析

MT19937 算法的周期为 $ 2^{19937} - 1 $，该算法的优点如下：
a. 周期非常长，该周期为梅森素数
b. 在 1 <= k <= 623 都满足 k-分布
c. 产生速度快。

算法分为三个阶段：

  • a. 初始化，获得基础的梅森旋转链
  • b. 对旋转链进行旋转算法
  • c. 对旋转算法所得结果进行处理，提取伪随机数

初始化

初始化 state = 0, 初始化向量 iv 含有 n 个元素，每一个元素为一个 word (32位或者64位) 长。其中将 seed 设置为 iv 中的第一个向量，剩下的 n - 1 个 word 按照如下的方式初始化：
$$ x_i = f \times (x_{i-1} \bigoplus (x_{i-1} >> (w - 2))) + i $$
其中 f 为常数 1812433253。

从而得到一条初始化的梅森链。

梅森旋转

对于 k 从 0 到 n - 1:
$$ \vec x_k = \vec x_{k+m} \bigoplus ((\vec x_k ^{(u)}||\vec x_{k+1} ^{(l)})A) $$
其中，$ \vec x_k ^{(u)} $ 代表第 k 个向量的高 $ w - r $ 位，$ \vec x_{k+1} ^{(l)} $ 代表第 k + 1 个向量的低 r 位。A 矩阵定义如下：

$$A = \begin{pmatrix}
0&I_{w-1}\\
a_{w-1}&(a_{w-2},…,a_0)\\
\end{pmatrix}$$

因此，$ \vec xA $ 可以直接表示为：
$$ xA = \begin{cases}
x >> 1 & x_0 = 0 \\
(x >> 1) \bigoplus a & x_0 = 1
\end{cases}$$

其中 $x_0$ 是向量 $\vec x$的最低位；a 为常量。

获取随机数

首先，在梅森链中获取 state 位置的 word，赋值给 x，然后通过下面的方式获取最终的随机数：
$$\begin{align}
y \equiv x \bigoplus ((x >> u) & d) \\
y \equiv y \bigoplus ((y << s) & b) \\
y \equiv y \bigoplus ((y << t) & c) \\
z \equiv y \bigoplus ((y >> l)) \\
\end{align}$$

获取完随机数后，state + 1。

代码实现

参考 梅森旋转算法

32位版本

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
62
63
64
65
class MT19937:  
    # w: word 大小  
    # n: 参与梅森旋转的随机数个数  
    # m: 中间词, 在[1,n)  
    # r: [0,w-1)    w, n, m, r = 32, 624, 397, 31  
  
    a = 0x9908b0df  
    u, d = 11, 0xffffffff  
    s, b = 7, 0x9d2c5680  
    t, c = 15, 0xefc60000  
    l = 18  
    f = 1812433253  
  
    def __init__(self):  
        seed = 0  
        self.state = 0  
        # 初始化，生成n行w列的向量  
        self.iv = [0] * self.n  
        self.iv[0] = seed  
        for i in range(1, self.n):  
            prev = self.iv[i - 1]  
            temp = self.f * (prev ^ (prev >> (self.w - 2))) + i  
            self.iv[i] = temp & 0xffffffff  # 转换为int32  
  
    def twister(self):  
        # 梅森旋转  
        for i in range(self.n):  
            y = (self.iv[i] & 0x80000000) + (self.iv[(i + 1) % self.n] & 0x7fffffff)  
            if y % 2:  
                y >>= 1  
                y ^= self.a  
            else:  
                y >>= 1  
            self.iv[i] = self.iv[(i + self.m) % self.n] ^ y  
  
    def generate(self):  
        if self.state == 0:  # 此处为效率优化，可以不进行判断，每次生成随机数都进行旋转
            self.twister()  
  
        y = self.iv[self.state]  
        y = y ^ (y >> self.u) & self.d  
        y = y ^ (y << self.s) & self.b  
        y = y ^ (y << self.t) & self.c  
        y = y ^ (y >> self.l)  
  
        self.state = (self.state + 1) % self.n  
        return y  
  
    def __call__(self):  
        return self.generate()  
  
  
if __name__ == "__main__":  
    mt = MT19937()  
    tank = set()  
    kLen = 100000  
    odd = 0  
    for _ in range(kLen):  
        t = mt()  
        odd += 1 if t % 2 else 0  
        tank.add(t)  
        print(t)  
    print(odd / kLen)  
    print(f"set: {len(tank)}")  
    print(f"kLen: {kLen}")

64位版本：

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
62
63
64
65
66
class MT19937x64:  
    w, n, m, r = 64, 312, 156, 31  
  
    a = 0xb5026f5aa96619e9  
    u, d = 29, 0x5555555555555555  
    s, b = 17, 0x71d67fffeda60000  
    t, c = 37, 0xfff7eee000000000  
    l = 43  
    f = 6364136223846793005  
    upperMask = 0xFFFFFFFF80000000  
    lowerMask = 0x7FFFFFFF  
  
    def __init__(self):  
        seed = 0  
        self.state = 0  
        # 初始化，生成n行w列的向量  
        self.iv = [0] * self.n  
        self.iv[0] = seed  
        for i in range(1, self.n):  
            prev = self.iv[i - 1]  
            temp = self.f * (prev ^ (prev >> (self.w - 2))) + i  
            self.iv[i] = temp & 0xffffffffffffffff  # 转换为int64  
  
    def twister(self):  
        # 梅森旋转  
        for i in range(self.n):  
            y = (self.iv[i] & 0x80000000) + (self.iv[(i + 1) % self.n] & 0x7fffffff)  
            if y % 2:  
                y >>= 1  
                y ^= self.a  
            else:  
                y >>= 1  
            self.iv[i] = self.iv[(i + self.m) % self.n] ^ y  
  
    def generate(self):  
        if self.state == 0:   # 此处为效率优化，可以不进行判断，每次生成随机数都进行旋转
            self.twister()  
  
        y = self.iv[self.state]  
        y = y ^ (y >> self.u) & self.d  
        y = y ^ (y << self.s) & self.b  
        y = y ^ (y << self.t) & self.c  
        y = y ^ (y >> self.l)  
  
        self.state = (self.state + 1) % self.n  
  
        return y  
  
    def __call__(self):  
        return self.generate()  
  
  
if __name__ == "__main__":  
    mt = MT19937x64()  
    tank = set()  
    kLen = 100000  
    odd = 0  
    for _ in range(kLen):  
        t = mt()  
        odd += 1 if t % 2 else 0  
        tank.add(t)  
        print(t)  
    print(odd / kLen)  
    print(f"set: {len(tank)}")  
    print(f"kLen: {kLen}")