MT19937伪随机数生成算法题型学习
32位的MT19937的python代码:
def _int32(x):
return int(0xFFFFFFFF & x)
class MT19937:
# 根据seed初始化624的state
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
其中_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 >> 11
y = y ^ y << 7 & 2636928640
y = y ^ y << 15 & 4022730752
y = y ^ y >> 18
self.mti = (self.mti + 1) % 624
return _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 = res
for i in range(bits // shift):
tmp = res ^ tmp >> shift
return tmp
右移带掩码的:
def inverse_right_mask(res, shift, mask, bits=32):
tmp = res
for i in range(bits // shift):
tmp = res ^ tmp >> shift & mask
return tmp
左移不带掩码的:
def inverse_left(res, shift, bits=32):
tmp = res
for i in range(bits // shift):
tmp = res ^ tmp << shift
return tmp
左移带掩码的:
def inverse_left_mask(res, shift, mask, bits=32):
tmp = res
for i in range(bits // shift):
tmp = res ^ tmp << shift & mask
return 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 md5
def 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 Random
def invert_right(m,l,val=''):
length = 32
mx = 0xffffffff
if val == '':
val = mx
i,res = 0,0
while i*l<length:
mask = (mx<<(length-l)&mx)>>i*l
tmp = m & mask
m = m^tmp>>l&val
res += tmp
i += 1
return res
def invert_left(m,l,val):
length = 32
mx = 0xffffffff
i,res = 0,0
while i*l < length:
mask = (mx>>(length-l)&mx)<<i*l
tmp = m & mask
m ^= tmp<<l&val
res |= tmp
i += 1
return res
def 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 m
def clone_mt(record):
state = [invert_temper(i) for i in record]
gen = Random()
gen.setstate((3,tuple(state+[0]),None))
return gen
f = 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+=1
g = 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 >> 11
y = y ^ y << 7 & 2636928640
y = y ^ y << 15 & 4022730752
y = y ^ y >> 18
self.mti = (self.mti + 1) % 624
return _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 = 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
当一步一步理解了分析以后，看这个代码就会很轻松了。
Case4 Reverse_init()
我们再来回顾一下MT19937的代码组成:
def _int32(x):
return int(0xFFFFFFFF & x)
class MT19937:
# 根据seed初始化624的state
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
extract_number() 我们在Case1中已经分析完了，twist()我们在刚刚的Case3中也分析完了，现在整个MT19937还剩下__init__函数我们没有分析。
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)
此中的关键在于最后一行代码
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 = res
for i in range(bits // shift):
tmp = res ^ tmp >> shift
return tmp
相比于之前的情况就是多了三部分,_int32,1812433253*,+i这三部分操作。
_int32相当于
(
m
o
d
0
x
F
F
F
F
F
F
F
F
)
\pmod{0xFFFFFFFF}
(mod0xFFFFFFFF)所以相当于整个运算是在
(
m
o
d
0
x
F
F
F
F
F
F
F
F
)
\pmod{0xFFFFFFFF}
(mod0xFFFFFFFF)下进行的。
由于:
print(gmpy2.gcd(1812433253,0xffffffff))
#1
所以在模运算下，1812433253存在逆元，那整个过程这相当于解一个简单的一次同余式，然后再进行inverse_right()。
1812433253 x + i = s t a t e [ i ] ( m o d 0 x F F F F F F F F ) 1812433253x+i=state[i]\pmod{0xFFFFFFFF} 1812433253x+i=state[i](mod0xFFFFFFFF)
x = ( s t a t e [ i ] − i ) ∗ 1812433253 − 1 ( m o d 0 x F F F F F F F F ) x=(state[i]-i)*{1812433253}^{-1}\pmod{0xFFFFFFFF} x=(state[i]−i)∗1812433253−1(mod0xFFFFFFFF)
def Reverse_init(last):
n = 1<<32
inv = invert(1812433253,n)
for i in range(623,0,-1):
last = ((last-i)*inv)%n
last = invert_right(last,30)
return last
通过该函数我们就可以获取最初的seed了,达到了逆向init()函数的目的。
def init(seed):
mt = [0] * 624
mt[0] = seed
for 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
372
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
