梅森旋转算法

梅森旋转算法(Mersense Twister, MT)是一种通用的伪随机数生成器 (pseudo-random number generator, PRNG),由松本刚明和西村在 1997 年提出,算法基于扭曲广义反馈移位寄存器TGFSR。事实上,MT 名字的由来是因为其循环周期为 ,而这正是梅森素数。

MT 算法专为纠正旧版 PRNG 算法的缺陷而提出,其中最常用的版本是 MT19937 和 MT 19937-64 分别是 32bit 和 64bit。在现代语言中,C++,Python 等众多语言内部随机数的实现有很多都是梅森旋转算法。

下面我们来详细学习 MT 算法的流程。

算法可以分为 3 部分:

  • 初始化,生成一个梅森链
  • 梅森旋转
  • 获取随机数

下面以 32 bit 的 MT 19937 为例讲解这三个过程。

初始化

初始化 state=0,初始向量 iv 含有 n 个元素,每一个元素为一个 word (32 bit 长)。

设:

1
iv[0] = seed

也就是将种子设置为 iv 中的第一个向量。剩下的 n-1 个 word 采用如下的方式初始化。

其中 f 为常数 1812433253

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

梅森旋转

对于 k 从 0 到 n-1:

其中, 代表拼接操作,即 concatenation。 代表第 个向量的高 位, 代表第 个向量的低 位。A 矩阵定义如下:

因此, 可以表示为 (如果有点混乱就直接记忆下面的结果):

其中 是向量 的最低位。

代码如下:

1
2
3
4
5
6
7
8
9
10
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

获取随机数

在获取随机数时用到了大量的常数项:

首先在梅森链中选取 state 位置的 word,赋值给 x,然后通过下面的方式获取最终的随机数。

在获取完随机数后,state+1。

代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
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

以上就是梅森随机算法的 3 个阶段,里面固然存在很多不能理解的部分,但是我们可以先对其中的算法流程进行掌握,在今后的学习中,不断深入学习。

完整代码

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
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}")

一些思考

在上面的代码中,我们可以看到,梅森旋转貌似是"成批旋转"的。其实这只是为了提高代码效率,事实上我们也可以旋转一次生成一次随机数。

例如:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
def twister(self):  
# 梅森旋转
i = self.state
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):
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

参考

https://en.wikipedia.org/wiki/Mersenne_Twister#Algorithmic_detail https://liam.page/2018/01/12/Mersenne-twister https://blog.csdn.net/shshss64/article/details/127326357 https://blog.csdn.net/tick_tock97/article/details/78657851 https://www.cnblogs.com/CHNmuxii/p/12232475.html https://blog.csdn.net/ACdreamers/article/details/44656743 https://blog.csdn.net/wlmnzf/article/details/83109567 张琳琳. 一种基于混合映射策略的改进梅森旋转演算法[J]. 信息技术,2022,46(10):12-17,23. DOI:10.13274/j.cnki.hdzj.2022.10.003.


梅森旋转算法
https://d4wnnn.github.io/2024/02/08/Academic/随机数相关/梅森旋转算法/
作者
D4wn
发布于
2024年2月8日
许可协议