生成和为0的随机数组_浅析MT19937伪随机数生成算法

 

引言

近年来MT19937在各大CTF赛事中出现的频率越来越高，本文主要讨论CTF中的常见的MT相关问题，并利用几道典型的赛题进行实战分析，尤其是今年二月份pwnhub的一道零解题(详见扩展题型)。感谢zzh师傅的细心指点，赛后才勉强解出了这道CoinFlip2。  

前置知识

MT19937是一种周期很长的的伪随机数生成算法，可以快速的产生高质量的伪随机数，主要分为三部分。 如果读者对该算法不了解，可以先参考 wiki 1.利用seed初始化624的状态
2.对状态进行旋转
3.根据状态提取伪随机数 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

python中内置的Random类就是采用了MT19937算法，getrandbits(32)方法可以获得一个32位随机数  

题型1 逆向 extract_number

首先分析extract_number函数，可以发现输出的伪随机数是对state[i]进行了异或，位运算后的结果。 逐步分析extract_number

y = y ^ (y >> 18)

我们可以发现这一步的运算对操作数y的高18位并没有产生影响，也就是说运算结束后 y1 = y ^ ( y>> 18) ，得到 y1 高18位其实就是 y 的高18位，那么我们就可以得到 y>>18 的高36位，进而得到 y1 的高36位，也就是说我们可以在有限步内，获得 y 的所有信息，即我们可以根据 y1 逆向 y 。代码如下：

o = 2080737669y = o^o>>18# 控制位移的次数for i in range(32//18):    y = y^(y>>18)print(y==o)

我们可以发现逆向的过程和正向的过程惊人的相似，不得不感叹数学的奇妙。 继续分析

y = y ^ y << 15 & 4022730752

可以看到这一步位移的方向发生改变，而且增加了掩码，根据运算符的优先级，可以得到 y1 = y^((y<<15)&4022730752) ，实际上增加的掩码并没有太大的作用，因为 y1 的低15位实际上就是 y 的低15位和 4022730752 的低15位异或运算的结果，我们只需要 y1^4022730752 便可以得到 y 的低15位，于是得到 y<<15 的低30位，同理可以得到 y1 的低30位，经过有限步，最终可以得到 y 的全部信息。代码如下：

o = 2080737669y = o ^ o << 15 & 4022730752tmp = yfor i in range(32 // 15):    # (y<<15)&40022730752 每次可以恢复y的15位    y = tmp ^ y << 15 & 4022730752print(y==o)

剩下的两步操作，可以采用同样的分析方法进行逆行。最终完整的代码如下：

o = 2080737669# right shift inversedef inverse_right(res, shift, bits=32):    tmp = res    for i in range(bits // shift):        tmp = res ^ tmp >> shift    return tmp# right shift with mask inversedef 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 inversedef inverse_left(res, shift, bits=32):    tmp = res    for i in range(bits // shift):        tmp = res ^ tmp << shift    return tmp# left shift with mask inversedef inverse_left_mask(res, shift, mask, bits=32):    tmp = res    for i in range(bits // shift):        tmp = res ^ tmp << shift & mask    return tmpdef extract_number(y):    y = y ^ y >> 11    y = y ^ y << 7 & 2636928640    y = y ^ y << 15 & 4022730752    y = y ^ y >> 18    return y&0xffffffffdef 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&0xffffffffy = extract_number(o)print(recover(y) == o)

上述讨论的部分是，基于对extract_number函数的观察得到的逆向方法，除此之外我们还可以从extract_number的运算本质上进行逆向。 设state[i]的二进制表示形式为：

输出的随机数二进制形式为：

而z和x具有如下线性关系： 也就是说

其中X,Z是GF(2)上的1 x 32的向量，T是GF(2)上的 32 x 32的矩阵。我们只需要在GF(2)上求解X即可。已知Z，如果T也已知，可以快速的求解出X。那么如何计算T呢？ 实际上我们可以采用黑盒测试的方法，猜解T。例如当X为(1,0,0,0,…..0)时，经过T得到的Z其实就是T中第一行。采用这种类似选择明文攻击的方法，我们可以得到T矩阵的每一行，进而还原T。最后再利用T和Z得到原始的X。代码如下：

# sagemath 9.0from sage.all import *from random import Randomdef buildT():    rng = Random()    T = matrix(GF(2),32,32)    for i in range(32):        s = [0]*624        # 构造特殊的state        s[0] = 1<31-i)        rng.setstate((3,tuple(s+[0]),None))        tmp = rng.getrandbits(32)        # 获取T矩阵的每一行        row = vector(GF(2),[int(x) for x in bin(tmp)[2:].zfill(32)])        T[i] = row    return Tdef reverse(T,leak):    Z = vector(GF(2),[int(x) for x in bin(leak)[2:].zfill(32)])    X = T.solve_left(Z)    state = int(''.join([str(i) for i in X]),2)    return statedef test():    rng = Random()    # 泄露信息    leak = [rng.getrandbits(32) for i in range(32)]    originState = [i for i in rng.getstate()[1][:32]]    # 构造矩阵T    T = buildT()    recoverState = [reverse(T,i) for i in leak]    print(recoverState==originState)test()

这里的黑盒测试思想，会在后面的题型x，得到应用。希望读者理解后再继续阅读  

题型2 预测随机数

这一类题型是基于对extract_number 函数的逆向，而发展来的。这里以2020年网鼎杯白虎组的random为例

2020网鼎杯白虎组random

定位到关键函数

def generate():    fw = open("random", "w")    for i in range(700):        fw.write(str(random.getrandbits(32))+"n")    fw.close()generate()f = open("flag.txt", "w")key = str(random.getrandbits(32))ciphertext = encryption(flag, key)f.write(ciphertext)f.close()

可以发现加密的key时第701个随机数，而且我们已知前700个随机数。于是可以根据输出的随机数逆向extract_number得到对应的700个state，实际上只需要前624个随机数恢复前624个state，就可以预测此后生成的随机数。代码如下：

# -*- coding: utf-8 -*-from random import Randomdef invert_right(m,l,val=''):    length = 32    mx = 0xffffffff    if val == '':        val = mx    i,res = 0,0    while i*l        mask = (mx<>i*l        tmp = m & mask        m = m^tmp>>l&val        res += tmp        i += 1    return resdef invert_left(m,l,val):    length = 32    mx = 0xffffffff    i,res = 0,0    while i*l < length:        mask = (mx>>(length-l)&mx)<        tmp = m & mask        m ^= tmp<        res |= tmp        i += 1    return 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",'r').readlines()prng = []for i in f:    i = i.strip('n')    prng.append(int(i))g = clone_mt(prng[:624])for i in range(700):    g.getrandbits(32)key = g.getrandbits(32)print(key)#2990136190

代码中逆向extract_number 部分的invert_right函数和题型1中的代码不同，但是效果是相同的(很久之前写的，所以直接粘过来了。)得到key之后就是常规的des解密，最后可以得到flag。  

题型3 逆向twist

在已知连续624个随机数时(state没有进行twist)，可以还原state，预测后续的随机数。那么如何获得624个随机数之前的随机数呢？此时可以考虑，逆向twist获得前一组的state，进而获得前一组的624个随机数。 首先看一下twist函数

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

考虑下面的例子：

1. 11100110110101000100101111000001 // state[i]2. 10101110111101011001001001011111 // state[i + 1]3. 11101010010001001010000001001001 // state[i + 397]// y = state[i] & 0x80000000 | state[i + 1] & 0x7fffffff4. 10101110111101011001001001011111 // y5. 01010111011110101100100100101111 // next = y >>> 16. 11001110011100100111100111110000 // next ^= 0x9908b0df0x9908b0df => 100110010000100010110000110111117. 00100100001101101101100110111001 // next ^= state[i + 397]

我们知道生成新的state[i]只与原来的state[i],state[i+1],state[i+397]有关。第7步是必须进行的一步(异或的次序不影响结果，所以异或state[i+397]可以看成最后一步),第6步是根据第4步结果的奇偶性确定的,不一定有第6步,但是因为第7步是第5步或者第6步异或state[i+397]的结果,我们可以考察新生成的state[i]异或state[i+397]的结果，来判断是否进行了第六步的操作。 由于0x9908b0df => 10011001000010001011000011011111,而第5步的最高位必定是0,但是如果执行了第6步那么执行后的结果首位则会变成1,于是我们可以根据第7步逆向出的结果的首位判断是否进行了第6步.进而推出第5步,第5步的后31位包含了state[i]的第1位和state[i+1]的第2位至第31位,根据第6步是否进行可以得到state[i+1]的最后1位,所以根据现在的state[i]和以前的state[i+397],可以获得原来state[i]的1位信息和state[i+1]的31位信息,要获得state[i]剩下的31位信息,需要对现在的state[i-1]进行同样的运算.当需要计算第一位state时,剩下的state都已经恢复了,可以利用恢复了的最后一位state获得还未恢复的state[0]的后31位,非常巧妙! 实现代码如下：

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

2020 V&N 招新赛 Backtrace

# !/usr/bin/env/python3import randomflag = "flag{" + ''.join(str(random.getrandbits(32)) for _ in range(4)) + "}"with open('output.txt', 'w') as f:    for i in range(1000):        f.write(str(random.getrandbits(32)) + "n")print(flag)

题目信息很简练，我们需要恢复前四个随机数。我们需要连续的624个随机数才能获得上一组随机数的状态，这里只有1000个随机数，但是之后前4位丢失，而前4位产生的新状态对应第624，625，626，627位随机数的状态，而他们的状态是可逆的。所以不需要完整的连续624个随机数，也可以求解完整的state。代码如下：

#!/usr/bin/python3from random import Random# right shift inversedef inverse_right(res,shift,bits=32):    tmp = res    for i in range(bits//shift):        tmp = res ^ tmp >> shift    return tmp# right shift with mask inversedef inverse_right_values(res,shift,mask,bits=32):    tmp = res    for i in range(bits//shift):        tmp = res ^ tmp>>shift & mask    return tmp# left shift inversedef inverse_left(res,shift,bits=32):    tmp = res    for i in range(bits//shift):        tmp = res ^ tmp << shift    return tmp# left shift with mask inversedef inverse_left_values(res,shift,mask,bits=32):    tmp = res    for i in range(bits//shift):        tmp = res ^ tmp << shift & mask    return tmpdef backtrace(cur):    high = 0x80000000    low = 0x7fffffff    mask = 0x9908b0df    state = cur    for i in range(3,-1,-1):        tmp = state[i+624]^state[i+397]        # 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+624]^state[i+396]        # recover Y,tmp = Y        if tmp & high == high:            tmp ^= mask            tmp <<= 1            tmp |= 1        else:            tmp <<=1        res |= (tmp)&low        state[i] = res    return statedef recover_state(out):    state = []    for i in out:        i = inverse_right(i,18)        i = inverse_left_values(i,15,0xefc60000)        i = inverse_left_values(i,7,0x9d2c5680)        i = inverse_right(i,11)        state.append(i)    return statef = open("output.txt","r").readlines()c = []for i in range(1000):    c.append(int(f[i].strip()))partS = recover_state(c)state = backtrace([0]*4+partS)[:624]# print(state)prng = Random()prng.setstate((3,tuple(state+[0]),None))flag = "flag{" + ''.join(str(prng.getrandbits(32)) for _ in range(4)) + "}"print(flag)

 

题型4 逆向init

其实这一部分还没有出现过类似的题目，这里就我的理解，简单的叙述一下。 前文分别叙述了，根据输出逆向对应的state，以及根据输出逆向上一组的state。此部分则是根据第一次的state，逆向seed。 首先看一下，init函数

def _int32(x):    return int(0xFFFFFFFF & x)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

定位到关键函数 mt[i] = _int32(1812433253 * (mt[i - 1] ^ mt[i - 1] >> 30) + i) ，可以发现其中的 (mt[i - 1] ^ mt[i - 1] >> 30) 是可逆的(详见题型1)，注意到 _int32 其实相当于取低32位，相当于mod 2^32。而 gcd(1812433253,2**32) == 1 于是 1812433253 存在逆元。那么逆向的过程其实就是模逆操作再组合上题型1的操作，便可以得到seed。代码如下：

from gmpy2 import invertdef _int32(x):    return int(0xFFFFFFFF & x)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 mtseed = 2080737669def invert_right(res,shift):    tmp = res    for i in range(32//shift):        res = tmp^res>>shift    return _int32(res)def recover(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 laststate = init(seed)print(recover(state[-1]) == seed)

 

扩展题型

在题型2中我们采用的黑盒测试的方法，找到了状态和输出的转移矩阵T。这是因为每个状态的二进制位与其输出的随机数的二进制位，存在线性关系。而每次twist后新生成的状态都是前一组状态的线性组合，也就是说任意一组状态都是第一组状态的线性组合。那么由任意一组状态生成的伪随机数的二进制位和第一组状态的二进制位存在线性关系。 不妨设 初始状态的二进制形式为： 设19968个伪随机数的MSB (注意这里不要求一定是MSB，可以是随机数的任意二进制位，也不一定要求是19968个随机数，只要能构成19968个二进制位即可)构成的向量为： 于是有在GF(2)上有如下等式,其中X,Z为1 x 19968 的向量，T为19968 x 19968 的矩阵：

采用题型2中同样的方法，对T进行黑盒测试。即设置X(1,0,0,0,…0)得到Z为T的第一行，遍历所有情况得到T。最后还原出X。

2020 pwnhub CoinFlip2

#!/usr/bin/env -S python3 -uimport randomcount = 0while True:    print("Progress: %s/50" % count)    coin = random.getrandbits(1)    if int(input("Your guess: ")) == coin:        print("Correct!")        count += 1        if count == 50:            print(open("flag").read())            exit()    else:        print("Wrong!")        count = 0

python中的 getrandbits(1) 其实截取的是 getrandbits(32) 的最高位即 MSB 。于是我们需要收集19968个 MSB ，并构造 T ，才能还原state。需要注意的是这里构造的 T 并不是满秩的，但是sage仍然可以求解，不过求解得到的 X 的前32位，并不正确。但是我们仍然能够得到623位的state，我们可以利用题型3中的方法还原第一位state。进而利用state预测接下来的随机数。不过需要注意的是，构造矩阵 T 是相当耗费时间的，笔者在2h2g的机器上跑了大约一个多小时，才算出结果。为了读者的方便，这里直接给出矩阵的数据 坚果云(需要注册，但速度很快访问密码 : badmonkey) 腾讯云(密码 monkey) 具体代码如下，可以根据需要更改：

#! /bin/bash/env python3from sage.all import *from random import Randomfrom tqdm import tqdmprng = Random()length = 19968def myState():    state = [0]*624    i = 0    while i        ind = i//32        expont = i%32        state[ind] = 1<31-expont)        s = (3,tuple(state+[0]),None)        yield s        state[ind] = 0        i += 1def getRow():    rng = Random()    gs = myState()    for i in range(length):        s = next(gs)        rng.setstate(s)#         print(s[1][0])        row = vector(GF(2),[rng.getrandbits(1) for j in range(length)])        yield rowdef buildBox():    b = matrix(GF(2),length,length)    rg = getRow()    for i in tqdm(range(length)):        b[i] = next(rg)    return bdef test():    prng = Random()    originState = prng.getstate()    # 这里都是用的MSB,如果采用不同的二进制位(如LSB)最后的矩阵T 也会不同    leak = vector(GF(2),[prng.getrandbits(1) for i in range(length)])    b = buildBox()    f = open("Matrix","w")    for i in range(b.nrows()):        for j in range(b.ncols()):            f.write(str(b[i,j])+"n")    f.close()    x = b.solve_left(leak)    x = ''.join([str(i) for i in x])    state = []    for i in range(624):        tmp = int(x[i*32:(i+1)*32],2)        state.append(tmp)    prng.setstate(originState)    prng.getrandbits(1)    originState = [x for x in prng.getstate()[1][:-1]]    print(originState[1:] == state[1:])#     print(state)    return state,btest()

利用脚本如下：

from sage.all import *from random import Randomfrom tqdm import tqdm# 根据文件中的信息，构造矩阵def buildMatrix():    length = 19968    cnt = 0    m = matrix(GF(2), length, length)    for line in tqdm(open("Matrix", "r")):        row = cnt // 19968        col = cnt % 19968        m[row, col] = int(line.strip('n'))        cnt += 1    return mm = buildMatrix()# X = Z*(T^-1)def recoverState(leak):    x = m.solve_left(leak)    x = ''.join([str(i) for i in x])    state = []    for i in range(624):        tmp = int(x[i * 32:(i + 1) * 32], 2)        state.append(tmp)    return state# 根据题型2,还原state,有两种可能,这时候可以用暴破def backfirst(state):    high = 0x80000000    low = 0x7fffffff    mask = 0x9908b0df    tmp = state[623] ^ state[396]    if tmp & high == high:        tmp ^= mask        tmp <<= 1        tmp |= 1    else:        tmp <<= 1    return (1 << 32 - 1) | tmp & low, tmp & lowdef pwn(leak):    state = recoverState(leak)    L = [leak[i] for i in range(100)]    prng = Random()    guess1, guess2 = backfirst(state)    print(guess1, guess2)    state[0] = guess1    s = state    prng.setstate((3, tuple(s + [0]), None))    g1 = [prng.getrandbits(1) for i in range(100)]    if g1 == L:        print("first")        prng.setstate((3, tuple(s + [0]), None))        return prng    state[0] = guess2    s = state    prng.setstate((3, tuple(s + [0]), None))    g2 = [prng.getrandbits(1) for i in range(100)]    if g2 == L:        print("second")        prng.setstate((3, tuple(s + [0]), None))        return prngdef test():    length = 19968    prng = Random()    originState = prng.getstate()    leak = vector(GF(2), [prng.getrandbits(1) for i in range(length)])    # 恢复state    state = recoverState(leak)    prng.setstate(originState)    prng.getrandbits(1)    originState = [x for x in prng.getstate()[1][:-1]]    # 成功恢复623个state    print(originState[1:] == state[1:])    # 获取泄露信息    L = [leak[i] for i in range(100)]    # 两种可能    guess1, guess2 = backfirst(state)    print(guess1, guess2)    state[0] = guess1    s = state    prng.setstate((3, tuple(s + [0]), None))    g1 = [prng.getrandbits(1) for i in range(100)]    if g1 == L:        print("first")        prng.setstate((3, tuple(s + [0]), None))        now = vector(GF(2), [prng.getrandbits(1) for i in range(length)])        if now == leak:            print("true")            return    state[0] = guess2    s = state    prng.setstate((3, tuple(s + [0]), None))    g2 = [prng.getrandbits(1) for i in range(100)]    if g2 == L:        print("second")        prng.setstate((3, tuple(s + [0]), None))        now = vector(GF(2), [prng.getrandbits(1) for i in range(length)])        if now == leak:            print("true")            returntest()

大概需要6-7分钟可以构建完矩阵。  

结语

本人水平有限，如有错误欢迎指出，希望大家能从中学到一些东西。

点击此处“阅读全文”查看更多内容