【MT19937】学习分析

MT19937伪随机数生成算法题型学习

32位的MT19937的python代码:

def _int32(x):return int(0xFFFFFFFF & x)class MT19937:# 根据seed初始化624的statedef __init__(self, seed):self.mt = [0] * 624self.mt[0] = seedself.mti = 0for 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 >> 11y = y ^ y << 7 & 2636928640y = y ^ y << 15 & 4022730752y = y ^ y >> 18self.mti = (self.mti + 1) % 624return _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

其中_int32()函数固定了输出位数为32位。

__init__构造函数()通过初始种子seed来生成后面的623位mt[i],有一种密钥扩展的感觉,得到了一个state块。

twist()函数对状态进行旋转，该旋转与第i块，第i+1块,第i+397块有关(皆为state块中的循环数，比如624块的下一块为第0块)

extract_number()函数通过状态量mti选出state中的第mti块，然后将该块进行一些逻辑运算以后输出32位，当624块用完时，即当mti=0或回到0时，进行twist（），state的刷新。

参考文章:https://www.anquanke.com/post/id/205861#h2-2

本文中的函数分析参考都来自于[badmonkey]https://www.anquanke.com/member/146200)

Case1 reverse_extract_number()

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

倒着分析一下

y = y ^ y >> 18

y ^ (y>>18)，由于y是32位的，向右移18位后与原来异或，不改变y的高18位，而低14位也可以通过异或得到，方法可以通过将y的高18位中的高14位取出后，与得到的y异或即可。

大概是这样

y=_int32(1232132132132120)
temp = y
y = y^(y>>18)
x = ((y>>14))>>4
y = y^x
print(bool(y==temp))
#True

第一步分析是这样子的，中间的y两次位移是为了把思路理清楚，整理一下可以变成x=(y>>18)

所以逆推变成了

y=_int32(1232132132132120)
temp = y
y = y^(y>>18)
y = y^(y>>18)
print(bool(y==temp))
#True

和正向的时候是一样的!

y1 = y ^ y << 15 & 4022730752

y1 = y^ [(y<<15)&4022730752]

4022730752=1110 1111 1100 0110 0000 0000 0000 0000

可以看到该掩码的低17位全为0，而在其中y向左移动了15位，则最后获得的y的低15位是不变的。

然后我们就可以像第一轮正向的时候，把y1与低15位左移15位后&4022730752后异或,得到高30位。

但总长度是32位，我们同理可以再来一次就行。

y=_int32(1232132132132120)
temp = y
y = y ^ y << 15 & 4022730752
tmp = y
for i in range(32 // 15):# (y<<15)&40022730752 每次可以恢复y的15位y = tmp ^ y << 15 & 4022730752
print(bool(y==temp))
#True

y = y ^ y << 7 & 2636928640

y = y^ [(y<<7)&2636928640]

2636928640= 1001 1101 0010 1100 0101 0110 1000 0000

这个刚好是低7位为0，所以可以得到最初的低7位的值，然后得到14位，21位，28位，32位。

y=_int32(1232132132132120)
temp = y
y = y ^ y << 7 & 2636928640
tmp = y
for i in range(32 // 7):#(y<<15)&40022730752 每次可以恢复y的15位y = tmp ^ y << 7 & 2636928640 #这种是每次只能恢复7位
print(bool(y==temp))
#True

y = y ^ y >> 11

y = y ^ (y>>11) 高11位不变

与一开始的有点不一样，这里要通过高11位得到高22位，然后得到32位

y=_int32(1232132132132120)
temp = y
y = y^y>>18
# 控制位移的次数
for i in range(32//18):y = y^(y>>18)
print(bool(y==temp))
#True

所以可以总结出四种情况:

右移不带掩码的:

    def inverse_right(res, shift, bits=32):tmp = resfor i in range(bits // shift):tmp = res ^ tmp >> shiftreturn tmp

右移带掩码的:

def inverse_right_mask(res, shift, mask, bits=32):tmp = resfor i in range(bits // shift):tmp = res ^ tmp >> shift & maskreturn tmp

左移不带掩码的:

def inverse_left(res, shift, bits=32):tmp = resfor i in range(bits // shift):tmp = res ^ tmp << shiftreturn tmp

左移带掩码的:

def inverse_left_mask(res, shift, mask, bits=32):tmp = resfor i in range(bits // shift):tmp = res ^ tmp << shift & maskreturn tmp

所以可以把reverse_extract_number写成这样:

def reverse_extract_number(m):m = inverse_right(m,18)m = inverse_left_mask(m,15,4022730752)m = inverse_left_mask(m,7,2636928640)m = inverse_right(m,11)return _int32(m)

基于以上分析发展而来了第二类问题

Case2 预测随机数

相关的例子在之前的文章Random里面也提到，也是从那里了解到了MT19937，用到了这里的解密脚本，现在回过头来分析一下细节上的东西。

题目代码:

import random
from hashlib import md5def get_mask():file = open("random.txt","w")for i in range(104):file.write(str(random.getrandbits(32))+"\n")file.write(str(random.getrandbits(64))+"\n")file.write(str(random.getrandbits(96))+"\n")file.close()
get_mask()flag = md5(str(random.getrandbits(32)).encode()).hexdigest()
print(flag)

解密脚本:

from random import Randomdef invert_right(m,l,val=''):length = 32mx = 0xffffffffif val == '':val = mxi,res = 0,0while i*l<length:mask = (mx<<(length-l)&mx)>>i*ltmp = m & maskm = m^tmp>>l&valres += tmpi += 1return resdef invert_left(m,l,val):length = 32mx = 0xffffffffi,res = 0,0while i*l < length:mask = (mx>>(length-l)&mx)<<i*ltmp = m & maskm ^= tmp<<l&valres |= tmpi += 1return resdef invert_temper(m):m = invert_right(m,18)m = invert_left(m,15,4022730752)m = invert_left(m,7,2636928640)m = invert_right(m,11)return mdef clone_mt(record):state = [invert_temper(i) for i in record]gen = Random()gen.setstate((3,tuple(state+[0]),None))return genf = open("random.txt",'r').readlines()
prng = []
j=0
for i in f:i = i.strip('\n')print(int(i).bit_length())if(j%3==0):prng.append(int(i))elif(j%3==1):#将生成两次随机数的两个随机数分离出来prng.append(int(i)& (2 ** 32 - 1))prng.append(int(i)>> 32)else:#将生成三次随机数的三个随机数分离出来prng.append(int(i)& (2 ** 32 - 1))prng.append(int(i)& (2 ** 64 - 1) >> 32)prng.append(int(i)>>64)j+=1g = clone_mt(prng[:624])
for i in range(624):g.getrandbits(32)#产生前624个随机数，让state状态到生成flag前key = g.getrandbits(32)
print(key)
from hashlib import md5
flag = md5(str(key).encode()).hexdigest()
print(flag)
#14c71fec812b754b2061a35a4f6d8421

现在我们来一部分一部分的看一看

invert_right()和invert_left()这里是合并写法，相当于把我们上面分析的四个函数合并成了两个，效果都是一样的

可以直接把上面的代码替换下来也可。

clone_mt()函数

def clone_mt(record):state = [invert_temper(i) for i in record]gen = Random()gen.setstate((3,tuple(state+[0]),None))return(gen)

Random().setstate()，用来设置当前随机状态，与之对应的还有Random().getstate(),获取当前随机状态

比如:

a=[1,2,3,4,5,6]
b=[1,2,3,4,5,6]
state=random.getstate()
random.shuffle(a)
print(a)
random.shuffle(b)
print(b)
b=[1,2,3,4,5,6]
random.setstate(state)
random.shuffle(b)
print(b)[5, 3, 4, 1, 6, 2]
[1, 6, 4, 3, 5, 2]
[5, 3, 4, 1, 6, 2]

看一看state是怎样构成的

print(state)
#(3, (2147483648, 1853904046, ......,624),None)所以我们按照这个格式把state传进去就可以了
gen.setstate((3,tuple(state+[0]),None))
最后加一个0表示当前所在位置，即种子seed所在的0位

然后当state一样的时候，输出的随机数就会按照我们设定好的种子，或者是已知随机数序列来逆推出的种子来预测随机数啦。

比如例题中的

g = clone_mt(prng[:624])
for i in range(624):g.getrandbits(32)#产生前624个随机数，让state状态到生成flag前key = g.getrandbits(32)

生成的key就是我们"预测"出来的随机数。

Case3 reverse_twist()

由case2 预测随机数部分分析我们可以得知，如果已知state状态的话（已知连续的624个随机数），我们就可以把之后产生的随机数都预测出来，但是该怎么预测出之前的随机数呢？

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

由extract_number()函数可知，当mti==0时，状态量会发生一个旋转，即调用twist()函数。每当产生一个随机数，mti值就加1，当产生624次以后，mti归0,调用twist()函数。

所以预测出之前产生随机数的关键在于，怎么将twist()函数逆回去，获取前一轮的state状态

    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

我们在开头简单分析过，twist()函数对状态进行旋转，该旋转与第i块，第i+1块,第i+397块有关(皆为state块中的循环数，比如624块的下一块为第0块)

我们跟着文章中的例子一起来分析一下

1. 11100110110101000100101111000001 // state[i]
2. 10101110111101011001001001011111 // state[i + 1]
3. 11101010010001001010000001001001 // state[i + 397]// y = state[i] & 0x80000000 | state[i + 1] & 0x7fffffff
4. 10101110111101011001001001011111 // y
5. 01010111011110101100100100101111 // next = y >> 1
6. 11001110011100100111100111110000 // next ^= 0x9908b0df
7. 00100100001101101101100110111001 // next ^= state[i + 397]

第1、2、3步中就是与twist()相关的块(当i=i时)

在第4步中得到的y就可以进行逻辑运算得到最后的mt[i]，且之后都没有改变过y的值

根据这一点以及异或的顺序不改变，在此例中原作者将第6、7步与twist()函数中互换了位置，便于分析

第6步是否运行取决于第5步得出来的y：

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

而第7步是每一次都需要运行的

我们可以通过将得到的新state[i]与state[i+397]异或回到第七步之前的状态。(这里需要注意的是，twist()中的每一步只改变state[i],不改变state[i+1]和state[i+397])

第七步之前我们进行了两步关键的操作,一步是y>>1，必须要进行的操作和异或0x9908b0df这一步可能进行的操作

怎么判断有没有进行第六步呢

分析一下：

在判断有没有进行第六步的时候，此时y的高位一定是0，因为刚刚才进行了右移1位的操作。

我们分两种情况：

第一种：进行了第六步异或操作

0x9908b0df =10011001000010001011000011011111

异或以后，得到的值的高位会是0^1=1，则高位为1。

第二种：没有进行第六步异或操作

则高位保持第五步时的0，则高位为0。

所以我们可以通过异或state[i+397]之后的结果，查看该结果的高位来判断是否进行了第六步异或操作，进而判断我们逆向的时候是否应该再异或0x9908b0df，从而推出第五步移位之后的值，得到该值相当于我们得到了y的高31位，但低位的1位我们不知道。

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

此时就要通过这一步的判断条件来恢复最低位，如果是异或了0x9908b0df，则我们的低位为0，若没有，则是1。

到此我们就逆向得到了y。

再来看看y的生成

y = state[i] & 0x80000000 | state[i + 1] & 0x7fffffff
y = _int32((self.mt[i] & 0x80000000) + (self.mt[(i + 1) % 624] & 0x7fffffff))
0x80000000=10000000000000000000000000000000
0x7fffffff=01111111111111111111111111111111

很显然可以看出，y的最高位是state[i]的最高位，y的低31位是state[i+1]的低31位。

我们的目的是获得32位的state[i]，现在还差31位该怎么获得呢，看到31这个数字，我们刚刚才提到过，y的低31位是state[i+1]的低31位

那如果想要state[i]的低31位，是不是在逆state[i-1]的时候就可以得到。这一巧妙的处理让我们获得了高1位和低32位。

当需要计算第1位state[0]时,剩下的state都已经恢复了,可以利用恢复了的最后一位state获得还未恢复的state[0]的后31位,即state[0-624]=state[624]

分析得出的代码如下：

def reverse_twist(cur):high = 0x80000000low = 0x7fffffffmask = 0x9908b0dfstate = curfor i in range(623,-1,-1):tmp = state[i]^state[(i+397)%624]# recover Y,tmp = Yif tmp & high == high:tmp ^= masktmp <<= 1tmp |= 1else:tmp <<=1# recover highest bitres = 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 = Yif tmp & high == high:tmp ^= masktmp <<= 1tmp |= 1else:tmp <<=1res |= (tmp)&lowstate[i] = res    return state

当一步一步理解了分析以后，看这个代码就会很轻松了。

Case4 Reverse_init()

我们再来回顾一下MT19937的代码组成:

def _int32(x):return int(0xFFFFFFFF & x)class MT19937:# 根据seed初始化624的statedef __init__(self, seed):self.mt = [0] * 624self.mt[0] = seedself.mti = 0for 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 >> 11y = y ^ y << 7 & 2636928640y = y ^ y << 15 & 4022730752y = y ^ y >> 18self.mti = (self.mti + 1) % 624return _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

extract_number() 我们在Case1中已经分析完了，twist()我们在刚刚的Case3中也分析完了，现在整个MT19937还剩下__init__函数我们没有分析。

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

此中的关键在于最后一行代码

self.mt[i] = _int32(1812433253 * (self.mt[i - 1] ^ self.mt[i - 1] >> 30) + i)

基于Case1、2中的分析，self.mt[i - 1] ^ self.mt[i - 1] >> 30可逆，函数可以调用

def inverse_right(res, shift, bits=32):tmp = resfor i in range(bits // shift):tmp = res ^ tmp >> shiftreturn tmp

相比于之前的情况就是多了三部分,_int32,1812433253*,+i这三部分操作。

_int32相当于(mod0xFFFFFFFF)\pmod{0xFFFFFFFF}(mod0xFFFFFFFF)所以相当于整个运算是在(mod0xFFFFFFFF)\pmod{0xFFFFFFFF}(mod0xFFFFFFFF)下进行的。

由于:

print(gmpy2.gcd(1812433253,0xffffffff))
#1

所以在模运算下，1812433253存在逆元，那整个过程这相当于解一个简单的一次同余式，然后再进行inverse_right()。

1812433253x+i=state[i](mod0xFFFFFFFF)1812433253x+i=state[i]\pmod{0xFFFFFFFF}1812433253x+i=state[i](mod0xFFFFFFFF)

x=(state[i]−i)∗1812433253−1(mod0xFFFFFFFF)x=(state[i]-i)*{1812433253}^{-1}\pmod{0xFFFFFFFF}x=(state[i]−i)∗1812433253−1(mod0xFFFFFFFF)

def Reverse_init(last):n = 1<<32inv = invert(1812433253,n)for i in range(623,0,-1):last = ((last-i)*inv)%nlast = invert_right(last,30)return last

通过该函数我们就可以获取最初的seed了,达到了逆向init()函数的目的。

def init(seed):mt = [0] * 624mt[0] = seedfor i in range(1, 624):mt[i] = _int32(1812433253 * (mt[i - 1] ^ mt[i - 1] >> 30) + i)return mt
seed = 234567890
state = init(seed)
print(reverse_init(state[-1]) == seed)
#True