random模块内置伪随机数生成—MT19937算法

random模块内置伪随机数生成—MT19937算法

所谓MT19937算法也就是梅森旋转算法中的一种，是伪随机数发生算法

有关梅森旋转算法可以看：梅森旋转算法 - 维基百科，自由的百科全书 (wikipedia.org)

$SimplyanalyzeMT19937$​

Python应用MT19937算法生成范围在$[0,2^{32}-1]$​的均匀分布的32位整数（来自wiki）

def _int32(x):
    return int(0xFFFFFFFF & x)

class MT19937:
    def __init__(self, seed):
        self.mt = [0] * 624
        self.mt[0] = seed
        self.mti = 0
        for i in range(1, 624):
            self.mt[i] = _int32(1812433253 * (self.mt[i - 1] ^ self.mt[i - 1] >> 30) + i)


    def extract_number(self):
        if self.mti == 0:
            self.twist()
        y = self.mt[self.mti]
        y = y ^ y >> 11
        y = y ^ y << 7 & 2636928640
        y = y ^ y << 15 & 4022730752
        y = y ^ y >> 18
        self.mti = (self.mti + 1) % 624
        return _int32(y)


    def twist(self):
        for i in range(0, 624):
            y = _int32((self.mt[i] & 0x80000000) + (self.mt[(i + 1) % 624] & 0x7fffffff))
            self.mt[i] = (y >> 1) ^ self.mt[(i + 397) % 624]

            if y % 2 != 0:
                self.mt[i] = self.mt[i] ^ 0x9908b0df

可以看到代码分为三段

• 导入seed，初始化伪随机数发生器

可以看到每次生成mt[i]是由前一项进行运算得到的，而第一项也就是mt[0] == seed

总的生成了有$624$​​项整数的列表（以及其他参数）作为state；

state是伪随机数发生器的重要变量，它决定了输出怎样的随机数，也就是说知道了state你可以预测出它将要输出的数

Python中的random模块就是内置的MT19937算法

举个例子验证state在伪随机数发生器中的重要性

import random
random_object = random.Random() # 生成伪随机数对象,实际上不必手动调用此函数，在直接调用random的其他生成随机数的函数时，会自动调用此函数生成对象
ori_state = random_object.getstate() # 得到该伪随机数发生器的state,也就是状态量
# 输出随机数
random_object.getrandbits(32) # 4221216727
random_object.getrandbits(32) # 1519334418
random_object2 = random.Random()
random_object2.setstate(ori_state)
# 验证随机数是否相等
random_object2.getrandbits(32) # 4221216727
random_object2.getrandbits(32) # 1519334418

观察state的组成

print(ori_state)
# (3, (2147483648,...,624), None)

中间的元组类型有$625$​项，前$624$​项即是之前提到的整数，而第$625$​项（此时等于$624$）是该state此时对应的第几个整数，也就是说，再下一次执行extract_number()函数时的state里面应该调用第几个整数，换句话说，mti = 624

其他参数似乎是默认的，在之后我们求出$624$​个整数之后，需要按这个格式进行构造真正的state；为了方便叙述，之后的state都暂时指代这$624$个整数

• 进行twist()函数

注意是在extract_number()函数内可能会执行twist()，进行twist()函数需要满足的条件，也就是每$624$项作为一轮直接执行extract_number()函数里后面的位运算以及抽取随机数的操作，而下一轮$624$项需要进行twist()函数替换state

这样一来，就意味着每$624$项的列表作为一次state，当输出完$624$​个伪随机数之后，需要通过twist()函数变换state，变换完成之后再进行输出伪随机数，以此类推

设新生成的state列表为mt'，而原来的state列表为mt（注意在MT19937原代码中并没有赋值新的列表，只是不断地更新原来地列表；这里这样讲只是为了方便理解）

在twist()函数中，每生成的新的一个state列表中的mt'[i]，只和在原来的state中的mt[i+1]以及mt[i+397]有关（具体的相关计算在逆转部分分析）

而当i+1或者i+397大于$624$时，也就是超过了一个state列表总长度时，那么mt[i+1]或者mt[i+397]实际上是mt'[(i+1) % 624]或者mt'[(i+397) % 624]）

• 继续extract_number()函数

在判断执行twist()函数与否之后，从目前的state列表中单独抽取每个数进行运算，运算结果即是输出的随机数。也就是说，state列表里的每个数都单独的进行与常量的位运算。

$ReverseMT19937$

分别按对应的函数进行逆转（$reverse$​​部分涉及的思路以及脚本均来自badmonkey）

关于脚本中每个函数的具体分析请见原文：浅析MT19937伪随机数生成算法 - 安全客，安全资讯平台 (anquanke.com)

$Reverseextract\_number()$​

假设已知624个从开始连续输出的随机数，需要预测之后每一个输出的随机数；

我们只需要逆转输出的随机数组成state列表，然后设置伪随机数发生器的state即可，然后输出需要预测的随机数

# Author:badmonkey

# right shift inverse
def inverse_right(res, shift, bits=32):
    tmp = res
    for i in range(bits // shift):
        tmp = res ^ tmp >> shift
    return tmp


# right shift with mask inverse
def inverse_right_mask(res, shift, mask, bits=32):
    tmp = res
    for i in range(bits // shift):
        tmp = res ^ tmp >> shift & mask
    return tmp

# left shift inverse
def inverse_left(res, shift, bits=32):
    tmp = res
    for i in range(bits // shift):
        tmp = res ^ tmp << shift
    return tmp


# left shift with mask inverse
def inverse_left_mask(res, shift, mask, bits=32):
    tmp = res
    for i in range(bits // shift):
        tmp = res ^ tmp << shift & mask
    return tmp

def recover(y):
    y = inverse_right(y,18)
    y = inverse_left_mask(y,15,4022730752)
    y = inverse_left_mask(y,7,2636928640)
    y = inverse_right(y,11)
    return y&0xffffffff

执行逆转时直接调用recover()函数即可；

random_number = []
state = [recover(i) for i in random_number]

$Reversetwist$​

$Case1$

假设已知该state中的$624$​个整数，想要知道之前输出的随机数（而不是之后输出的随机数），那么我们需要逆转twist()得到上一个state或者更之前的state

def backtrace(cur):
    high = 0x80000000
    low = 0x7fffffff
    mask = 0x9908b0df
    state = cur
    for i in range(623,-1,-1):
        tmp = state[i]^state[(i+397)%624]
        # recover Y,tmp = Y
        if tmp & high == high:
            tmp ^= mask
            tmp <<= 1
            tmp |= 1
        else:
            tmp <<=1
        # recover highest bit
        res = tmp&high
        # recover other 31 bits,when i =0,it just use the method again it so beautiful!!!!
        tmp = state[i-1]^state[(i+396)%624]
        # recover Y,tmp = Y
        if tmp & high == high:
            tmp ^= mask
            tmp <<= 1
            tmp |= 1
        else:
            tmp <<=1
        res |= (tmp)&low
        state[i] = res    
    return state

使用方法：直接向该函数传入逆转extract_number()得到的state即可得到上一个state，然后再执行extract_number()输出对应的随机数即可

$Case2$

$Question:$​假设已知$1000$​个从开始第$5$​个输出的连续随机数，求前$4$​个输出的随机数的大小（来自V&N 2020题目backtrace）

$Analysis:$​那么这$1000$个随机数可以分成两段（对应不同的state生成的随机数），也就是前$620$个数为第一个state生成的，后$380$​​个数为第二个state生成的​；需要求第一个state生成的前四个随机数

这与之前的$Case1$不同，这里不知道一个完整state的$624$个整数

仔细观察twist()函数，与之前提到的一样，第二个state列表中的mt'[i]（注意在这里的前提条件下mt'[i] == mt[i+624]），只和在第一个的state中的mt[i+1]以及mt[i+397]有关

又因为此时相当于求在mt[i]中，$i$​等于$0$​到$3$​​​​的情况，那么我们只需要知道在已知$1000$个的随机数中第$397$到$400$个数以及第$624$到$628$​个数，还有第$1$个数；接下来就是相同的逆转，只是由于前提条件不同，部分参数改变了

这里与$Case1$相同，也从大到小进行逆转

def backtrace(cur):
    high = 0x80000000
    low = 0x7fffffff
    mask = 0x9908b0df
    state = cur
    for i in range(3,-1,-1):
        tmp = state[i+624]^state[i+397] # mt'[i] == mt [i+624]
        if tmp & high == high:
            tmp ^= mask
            tmp <<= 1
            tmp |= 1
        else:
            tmp <<=1
        res = tmp&high
        tmp = state[i-1+624]^state[i+396] # 这里也与Case 1不同，但是实际上是一个意思
        if tmp & high == high:
            tmp ^= mask
            tmp <<= 1
            tmp |= 1
        else:
            tmp <<=1
        res |= (tmp)&low
        state[i] = res
    return state

$Peference$

浅析MT19937伪随机数生成算法 - 安全客，安全资讯平台 (anquanke.com)