Explore MT19937

huangx607087对PRNG-MT19937的深一步学习

0x01 对MT19937更深的理解

MT19937，即梅森旋转算法，是一种高质量的伪随机数生成器，在2020年10月26日的PRNG MT Notes文章中，我们已经对这个算法的概念性的内容进行了初步讲述，现在对于MT19937的一些其他性质进行了解。

CTF中很多会让你预测随机数的情形，都是基于MT19937的。对于MT19937的预测，根据理论可知，我们只需要知道连续的$624$个$32$位二进制数，就可以获取该算法的后面的所有状态或者前面的所有状态。

而在实际应用中，并不是所有的伪随机数都是32位的，有的时候也会出现8位、16位、64位、96位等32的倍数或约数。当然也有可能出现9位、33位，34位、70位等不规则的情形。因此，上面提到的“连续$624$个$32$位二进制数”可以被替换成“连续$19968$个比特位”。（这里有个很神奇的地方，就是$623×32$的值恰好是$19936$）。

下面上一下MT19937的基本代码

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 getstate(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
for i in range(5):
    print(D.getstate())        

可以看出，在MT19937中，由于getstate函数并不是直接把内部数组的内容输出，而是经过一个可逆并且一一对应的处理后输出的。这就使得输出的数字随机性较高。

而如果我们想要生成的随机数的比特位数是$32$的倍数的话，那么这个随机数也是可以预测出来的。这个时候我们只需要将每$32$个比特位看作一个整体，然后按照最低$32$位先生成的原则，就可以预测出整个随机数序列了。也就是说，如果我们想要预测$64$位的随机数序列，那么我们只需要将每个$64$位的随机数拆成两个$32$位的数字，然后把代表着$2^{31}$到$2^{0}$的比特位放前面，$2^{63}$到$2^{32}$的比特位放后面，这样就可以恢复所有的$32$位序列进行预测了。同理可以预测$96,128$等$32$的倍数的随机数的序列。

0x02 CTF中MT19937常见题型

2o01 题型一：逆向随机数产生器内部的state中的内容

根据上面的代码分析内容可知，最后输出的内容并不是随机数产生时原本的state。根据去年10月26日的那篇文章的分析。我们利用下面的代码中的recover()函数和extract_number()函数，就可以实现内部state和实际输出之间数值的互相转换。（具体原理见去年10.26文章的分析）

def inverse_right(res, shift, mask=0xffffffff, bits=32):
    tmp = res
    for i in range(bits // shift):
        tmp = res ^ tmp >> shift & mask
    return tmp
def inverse_left(res, shift, mask=0xffffffff, bits=32):
    tmp = res
    for i in range(bits // shift):
        tmp = res ^ tmp << shift & mask
    return tmp
def extract_number(y):
    y = y ^ y >> 11
    y = y ^ y << 7 & 2636928640
    y = y ^ y << 15 & 4022730752
    y = y ^ y >> 18
    return y&0xffffffff
def recover(y):
    y = inverse_right(y,18)
    y = inverse_left(y,15,4022730752)
    y = inverse_left(y,7,2636928640)
    y = inverse_right(y,11)
    return y&0xffffffff

2o02 题型二：预测后续的随机数

我们可以直接用以下代码来进行预测随机数。也就是直接对已知的$624$个随机数处理成随机数产生器内部状态，然后进行一次推测就行了

def setstate(self,s):
    if(len(s)!=624):
        raise ValueError("The length of prediction must be 624!")
    for i in range(624):
        self.mt[i]=self.recover(s[i])
    self.mti=0
def predict(self,s):
    self.setstate(s)
    self.twist()
    return self.getstate(True)

2o03 题型三：求之前的随机数，即逆向twist()函数

我们先来观察一下这个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

很显然，对于内部的每一状态，该函数先继承了原状态的最高位，加上下一状态的最低$31$位，然后又是一次异或的过程，最后有根据临时值$y$的奇偶性又会对内部状态进行一次异或。由于异或过程都是可逆的，如果我们举一个特殊值进行简单的分析，可以得到：

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
0x9908b0df => 10011001000010001011000011011111
7. 00100100001101101101100110111001 // next ^= state[i + 397]

可以看出，处理后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 invtwist(self):
    high = 0x80000000
    low = 0x7fffffff
    mask = 0x9908b0df
    for i in range(623,-1,-1):
        tmp = self.mt[i]^self.mt[(i+397)%624]
        if tmp & high == high:
            tmp ^= mask
            tmp <<= 1
            tmp |= 1
        else:
            tmp <<=1
        res = tmp&high
        tmp = self.mt[i-1]^self.mt[(i+396)%624]
        if tmp & high == high:
            tmp ^= mask
            tmp <<= 1
            tmp |= 1
        else:
            tmp <<=1
        res |= (tmp)&low
        self.mt[i] = res

由此，我们可以上一下整个MT19937的工具包供使用

from Crypto.Util.number import *
from hashlib import md5
import random
def _int32(x):
    return int(0xFFFFFFFF & x)
class MT19937:
    def __init__(self, seed=0):
        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 getstate(self,op=False):
        if self.mti == 0 and op==False:
            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
    def inverse_right(self,res, shift, mask=0xffffffff, bits=32):
        tmp = res
        for i in range(bits // shift):
            tmp = res ^ tmp >> shift & mask
        return tmp
    def inverse_left(self,res, shift, mask=0xffffffff, bits=32):
        tmp = res
        for i in range(bits // shift):
            tmp = res ^ tmp << shift & mask
        return tmp
    def extract_number(self,y):
        y = y ^ y >> 11
        y = y ^ y << 7 & 2636928640
        y = y ^ y << 15 & 4022730752
        y = y ^ y >> 18
        return y&0xffffffff
    def recover(self,y):
        y = self.inverse_right(y,18)
        y = self.inverse_left(y,15,4022730752)
        y = self.inverse_left(y,7,2636928640)
        y = self.inverse_right(y,11)
        return y&0xffffffff
    def setstate(self,s):
        if(len(s)!=624):
            raise ValueError("The length of prediction must be 624!")
        for i in range(624):
            self.mt[i]=self.recover(s[i])
        #self.mt=s
        self.mti=0
    def predict(self,s):
        self.setstate(s)
        self.twist()
        return self.getstate(True)
    def invtwist(self):
        high = 0x80000000
        low = 0x7fffffff
        mask = 0x9908b0df
        for i in range(623,-1,-1):
            tmp = self.mt[i]^self.mt[(i+397)%624]
            if tmp & high == high:
                tmp ^= mask
                tmp <<= 1
                tmp |= 1
            else:
                tmp <<=1
            res = tmp&high
            tmp = self.mt[i-1]^self.mt[(i+396)%624]
            if tmp & high == high:
                tmp ^= mask
                tmp <<= 1
                tmp |= 1
            else:
                tmp <<=1
            res |= (tmp)&low
            self.mt[i] = res
def example():
    D=MT19937(48)
    print(D.getstate())
    print(D.mt[:5])
    print(D.recover(90324435))
    print(D.extract_number(90324435))
    D.twist()
    print(D.mt[:5])
    D.invtwist()
    print(D.mt[:5])
#Main Below#
example()