梅森旋转算法

梅森旋转算法(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 长)。

设：

iv[0] = seed

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

其中 f 为常数 1812433253。

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

梅森旋转

对于 k 从 0 到 n-1：

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

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

其中 是向量 的最低位。

代码如下：

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。

代码如下：

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 位版本：

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 位版本：

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

一些思考

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

例如：

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.