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的基本代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
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文章的分析）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
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$个随机数处理成随机数产生器内部状态，然后进行一次推测就行了

1
2
3
4
5
6
7
8
9
10
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()函数：

1
2
3
4
5
6
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
2
3
4
5
6
7
8
9
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位。

逆向代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
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的工具包供使用

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
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()

0x03 给出任意19937个bit-upd:2025.01.22

xs，3 年半了，这篇博客竟然还更新了

经过测试，可以得到下面的结论：

getrandbits(64) 中，后面 32 bit 为高位，前面 32 bit 为低位（96，128同理，越靠后的32bit组位置越高），可以直接按32bit原方法逆向。

getrandbits(16) 中，只返回高 16 bit，低 16 bit 被丢弃，而 getrandbits(33) 中，低 32 bit 为第一组比特，而高 1 位取自第二组的最高位bit，第二组低31bit被丢弃。

那，如果给出了任意不连续的 19937 个比特呢？我们可以用这个代码构造矩阵进行求解（此处数据包括1250个 getrandbits(16) 的结果，如果是其他奇奇怪怪的取值方式，则 getRows 函数的取样方法和取样顺序保持和题目中完全一致就OK了）

提示：在自己电脑上跑，耗时 $\mathsf{8 min}$ （比想象得快），内存消耗峰值为 $\mathsf{7GB540MB}$（也比想象的低，估计Sage对 $\mathbb F_2$ 上计算做了优化，但此代码切勿在内存小于 8G 的服务器上跑）。

REF:2024-同济大学第二届网络安全新生赛CatCTF-wp-crypto | 糖醋小鸡块的blog (tangcuxiaojikuai.xyz)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
Dall=list(map(int,open('data3.txt','r').readlines()))
from Crypto.Util.number import *
from random import *
from tqdm import *
n=1250
D=Dall[:n]
rng=Random()
def getRows(rng):
    #这一部分根据题目实际编写，必须和题目实际比特获取顺序和方式完全一致，且确保比特数大于19937，并且请注意zfill。
    row=[]
    for i in range(n):
        row+=list(map(int, (bin(rng.getrandbits(16))[2:].zfill(16))))
    return row
M=[]
for i in tqdm_notebook(range(19968)):#这一部分为固定套路，具体原因已经写在注释中了
    """
    referennce:
    糖醋小鸡块 2025/1/21 20:26:51
    这部分代码相当于取了一组线性基

    糖醋小鸡块 2025/1/21 20:26:56
    因为mt19937是线性的
    """
    state = [0]*624
    temp = "0"*i + "1"*1 + "0"*(19968-1-i)
    for j in range(624):
        state[j] = int(temp[32*j:32*j+32],2)
    rng.setstate((3,tuple(state+[624]),None)) #这个setstate也是固定格式，已于2025.1.21测试
    M.append(getRows(rng))
M=Matrix(GF(2),M)
y=[]
for i in range(n):
    y+=list(map(int, (bin(D[i])[2:].zfill(16))))
y=vector(GF(2),y)
s=M.solve_left(y)
#print(s)
G=[]
for i in range(624):
    C=0
    for j in range(32):
        C<<=1
        C|=int(s[32*i+j])
    G.append(C)
import random
RNG1 = random.Random()
for i in range(624):
    G[i]=int(G[i])
RNG1.setstate((int(3),tuple(G+[int(624)]),None))

print([RNG1.getrandbits(16) for _ in range(75)])
print(D[:75])
#[24611, 61540, 27000, 22927, 22627, 55288, 63953, 49489, 34260, 42784, 11039, 23642, 10110, 16227, 39041, 62962, 56808, 36780, 3159, 51528, 32460, 35716, 47277, 16978, 11858, 52799, 13693, 32000, 60366, 40153, 65005, 61189, 33210, 43238, 25824, 22460, 8863, 57878, 36381, 56343, 55634, 1060, 8091, 61140, 8566, 27917, 60023, 7332, 46069, 28344, 1102, 14223, 24904, 56001, 35032, 60197, 25690, 38921, 38024, 52458, 33339, 16684, 55115, 48897, 27067, 43509, 1101, 3168, 35237, 27651, 2668, 5881, 12436, 24064, 30439]
#[24611, 61540, 27000, 22927, 22627, 55288, 63953, 49489, 34260, 42784, 11039, 23642, 10110, 16227, 39041, 62962, 56808, 36780, 3159, 51528, 32460, 35716, 47277, 16978, 11858, 52799, 13693, 32000, 60366, 40153, 65005, 61189, 33210, 43238, 25824, 22460, 8863, 57878, 36381, 56343, 55634, 1060, 8091, 61140, 8566, 27917, 60023, 7332, 46069, 28344, 1102, 14223, 24904, 56001, 35032, 60197, 25690, 38921, 38024, 52458, 33339, 16684, 55115, 48897, 27067, 43509, 1101, 3168, 35237, 27651, 2668, 5881, 12436, 24064, 30439]
CTFSeason3
#Notes
Explore MT19937
http://example.com/2021/07/10/Explore-MT19937/
作者
huangx607087
发布于
2021年7月10日
许可协议