梅森旋转生成随机数

梅森旋转生成随机数的逆向预测和恢复

  • 预测下一个随机数，恢复被隐藏的随机数据

首发于先知社区 ：子集和问题的两种解决方式-先知社区

算法过程

  • wiki上有算法实现的源代码
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
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 extract_number(self,x=None):
        if self.mti == 0 and x is None:
            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

根据代码分析过程。因为社区已经有人讲过很详细的过程，我们这里就只说一下大致分为三个过程

  • 初始化

根据种子seed，生成初始状态mt(共有624个数据。然后mti==0.

  • 处理mt生成随机数

对mt进行一系列线性操作然后生成随机数。

  • 进行旋转

当已经生成了624个随机数时，相当于耗光了一轮的状态。就需要利用旋转生成新一轮数据，观察这个twist代码，我们可以利用这个预测624位之后出现的下一轮的mt状态然后预测下一个随机数。

1
2
3
4
5
6
7
8
9
10
11
M1=[]
M2=[]
rng=MT19937(seed=1123)
for i in range(1248):
    a=rng.extract_number()
    if i<624:
        M1.append(a)
    else:
        M2.append(a)
print(M1)
print(M2)

如图可以调用函数生成随机数我们可以利用M1预测M2的数据。

我们可以用RandCrack库预测生成的随机数和生成之前生成的随机数，也可以自己写一个预测的代码。

  • 库函数
1
2
3
4
5
rc=RandCrack()
for i in M1:
    rc.submit(i)
next=rc.predict_randrange(0, 0xFFFFFFFF)#(这里确定输出的起始和末尾bit,，这里在后面讲解为何这样)
print(next)
  • 逆向函数进行预测
  • 逆向extract_number()函数

​ y = y ^ y >> 11
​ y = y ^ y << 7 & 2636928640
​ y = y ^ y << 15 & 4022730752
​ y = y ^ y >> 18

这是这个函数最关键的部分，我们在一轮中可以利用生成一轮的随机数恢复mt，然后预测下一轮的mt恢复数据，从下往上。

我们把操作后的y写成bit形式y=$y_1..y_{32}$ |同时记操作前$x=x_1..x_{32}$。。。

​ $$y=x\oplus x>>n$$ n>=32-n

那么有$y_1..y_{32}=x_1..x_{32}\oplus 0..0x_1..x_{32-n}$

那么$y_1..y_n=x_1..x_n$ $x_{n+1}..x_{32}\oplus x_1..x_{32-n}=y_{n+1}..y_{32}$

可以得到$x_{n+1}..x_{32}=y_{n+1}..y_{32}\oplus y_1..y_{32-n}$ 因为n>32-n,所以$y_1..y_{32-n}$是已知的

所以$x=y>>n$

同理，当n<32-n ,情况有所转变

最后这里无法$y_{n+1}..y_{32}\oplus y_1..y_{32-n}$ 无法一次性完全恢复多轮几次就好了。

我们可以$x_{n+1}..x_{2n+1}=y_{n+1}..y_{32}\oplus y_1..y_{n+1}$

我们多进行几次这样的变化就行了。

逆向代码

1
2
3
4
5
6
7
8
9
10
 def re_right(self,x,bit,mask=0xffffffff):
     tmp=x
     for _ in range(32//bit):
         tmp=x^tmp>>bit&mask
     return tmp
def re_left(self,x,bit,mask=0xffffffff):
     tmp=x
     for _ in range(32//bit):
         tmp=x^tmp<<bit&mask
     return tmp

这里如果我们只进行正向的预测，就暂时不用反向写一个twist代码了

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
    def re_extract(self,m):
        m=self.re_right(m,18,0xffffffff)
        m=self.re_left(m,15,4022730752)
        m=self.re_left(m,7,2636928640)
        m=self.re_right(m,11)
        return m&0xffffffff
    def re_state(self,outputs):
        if len(outputs)!=624:
            raise ValueError("Invalid number of outputs")

        self.mt=[self.re_extract(m)for m in outputs]
        self.mti=0
        return self.mt
M1=[]
M2=[]
rng=MT19937(seed=1123)
for i in range(1248):
    a=rng.extract_number()
    if i<624:
        M1.append(a)
    else:
        M2.append(a)
print(M1)
print(M2)
pre=MT19937()
pre.re_state(M1)
print("预测 RNG:", [pre.extract_number() for _ in range(10)])

结合上述代码可以很好的预测下一个。不过我们要如果进行恢复前面的数和更精简的操作。便需要逆向twist

1
2
3
4
5
6
7
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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
def re_twist(self,mt):
    re_tw=[0]*624#生成列表
    a=deepcopy(mt[623])
    c=deepcopy(mt[396])
    for i in range(623,-1,-1):#从大到小遍历，以便twist[(i+397)%624]是符合条件的
        k=mt[i]^mt[(i+397) % 624]
        if (k&0x80000000)>>31==1:#判断y>>1的第一位
            k=k^0x9908b0df
            low=1
            high=(k&0x40000000)>>30
            re_tw[i]=high<<31
            re_tw[(i+1)%624]=re_tw[(i+1)%624]+((k&0x3fffffff)<<1)+low
            if i !=623:
                mt[(i+1)%624]=re_tw[(i+1)%624]#还原正确的依赖
        elif (k & 0x80000000) >> 31 == 0:
            low = 0
            high = (k & 0x40000000) >> 30
            re_tw[i] = high << 31
            re_tw[(i + 1) % 624] = re_tw[(i + 1) % 624] + ((k & 0x3fffffff) << 1) + low
            if i != 623:
                mt[(i + 1) % 624] = re_tw[(i + 1) % 624]

    return re_tw

这个反转函数，可读性不太好，我们从大到小遍历。

先将异或mt[(i+397) % 624] (从大到小的很明显的是能够直接利用这个反解，并且每一次循环后都会更改依赖)

不过这关有个缺陷就是不能还原出mt[0]，

因为最后一个计算mt[0]的时候会自动把re_tw[0]变成re_tw[i]=high<<31。

我恢复出来结果恢复出了它下一个的mt1[0]所以最后我换了一个方法。

比如我们要求mt[i],已知它旋转后mt1[]

我们利用mt1[i]恢复出mt[i]的最高的1位，然后再利用mt1[i-1]恢复出第的31位

具体恢复过程

1
2
0x7fffffff=0b01111111111111111111111111111111
0x80000000=0b10000000000000000000000000000000

如果这种操作应用给一个规定的32bit的数那么

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

在y等于mt[i]第32Bit 和.mt[(i + 1) % 624]的后31bit之和

由于都是线性操作，我们先

k=mt[i]^mt[(i+397) % 624]，

  • 如果k与0x9908b0df异或过，那么由于(y<<1)^0x9908b0df=k，

k<<31==1。

  • 如果k没与0x9908b0df异或

k<<31==0

便能通过y得到我们想要的信息。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
def re_tw1(mt):		
    high=0x80000000
    low=0x7fffffff
    mask=0x9908b0df
    for i in range(623,-1,-1):
        t=mt[i]^mt[(i+397) % 624]
        if t&high==high:
            t=t^mask
            t=t<<1
            t|=1#确定位奇数
        else:
            t=t<<1
        res=t&high #取得高位
        t=mt[i-1]^mt[(i+396) % 624]
        if t&high==high:
            t=t^mask
            t=t<<1
            t|=1
        else:
            t=t<<1
        res=res+(t&low)
        mt[i]=res
    return mt

验证代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
if __name__ == '__main__':
    M1 = []
    M2 = []
    rng = MT19937(seed=1123)
    for i in range(1248):
        a = rng.extract_number()
        if i < 624:
            M1.append(a)
        else:
            M2.append(a)
    print(f'M1={M1}')
    print(f'M2={M2}')
    pre = MT19937()
    m = pre.re_state(M2)
    pre.mt=re_tw1(m)
    print(pre.mt)
    print("预测 RNG:", [pre.extract_number(x=1) for _ in range(10)])

可以跑一下是一样的

这里的x=1是控制是否旋转的，具体翻阅上述源代码。

[TM19973/MT.py at main · rockfox0/TM19973]—这是MT19973的源码和逆向的代码的github链接

如果是不连续的输出信息

在一些比较难的梅森旋转的题目中，我们得到的随机数是间断，来自于不同的bit值。在正常情况下，我们如果知道了连续的19937个bit的信息(也就是一轮)，我们恢复之前的和预测之后的是很轻易的。

因为梅森旋转是线性操作，我们可以利用线性代数的知识来做这个。

假设已知的数据是b，b可以转化为一个向量。y表示初始的state。y经过M这个矩阵的线性运算可以得到b。

$ yM=b—y=b(M)^{-1}$

我们只需要构造出M表示线性运算的过程。

举个例子。。在python 代码 

1
a=getrandombits(t)

如果t<32，就会先生成正常的一组32bit数，然后把32bit右移16bit。就能得到16bit以内的a。如果我们得到的是一些t=20bit的数据想要恢复它，显然用前面直接逆向的方式是不太可行的。我们这是就可以构造矩阵M，利用线性代数的运算表示线性的运算过程。

  • 我们要先自定义一个随机数生成器的函数，来表示我们生成随机数的方式。
1
2
3
4
5
def getRows(rng):
    row=[]
    for i in range(n):
        row+=list(map(int, (bin(rng.getrandbits(20))[2:].zfill(20))))
    return row

这里表示就是生成20bit内的随机数然后累加1000个到row上得到的row就是一个有20000个0/1数据的列表

其实我们要求19968Bit的已知数据就够了，不过这里多点也行。

  • 然后就是构建M线性运算的矩阵,M将初始矩阵变成
1
2
3
4
5
6
7
8
for i in tqdm(range(19968)):
    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)) 
    M.append(getRows(rng))
M=Matrix(GF(2),M)

这个生成M线性运算的矩阵代码是不变的可以直接。无论生成随机数是隐藏了哪些部分，这些代码都是不变的。

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
from Crypto.Util.number import *
from random import *
from tqdm import *
n=1000
D=[]#1000个得到20
rng=Random()

def getRows(rng):
    row=[]
    for i in range(n):
        row+=list(map(int, (bin(rng.getrandbits(20))[2:].zfill(20))))
    return row
M=[]

for i in tqdm(range(19968)):
    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)) 
    M.append(getRows(rng))
M=Matrix(GF(2),M)
print("行数:", M.nrows())
print("列数:", M.ncols())
r=[]
for i in range(n):
    r+=list(map(int, (bin(D[i])[2:].zfill(20))))
r=vector(GF(2),r)
y=M.solve_left(r)
G=[]
for i in range(624):
    C=0
    for j in range(32):
        C<<=1
        C|=int(y[32*i+j])
    G.append(C)#将恢复的state转化为每个元素32bit的列表
import random
RNG1 = random.Random()
for i in range(624):
    G[i]=int(G[i])
RNG1.setstate((int(3),tuple(G+[int(624)]),None))#将state作为随机数生成器的初始状态，测试生成的state是否正确
P=[RNG1.getrandbits(20) for _ in range(75)]
print(P)
print(D[:75])

知道这两个就比较容易了剩下可以看注释

  • ps:

多数人的sagemath应该是装在wsl等其他虚拟环境平台上的，因为这个题目比较消耗内存，wsl的内存默认只给6.5G。这个已知20bit的还好，如果未知bit太多消耗内存肯定是不够的。建议把wsl内存上限拉到16G。网上有很多更改教程这里就不说了

tpctf 2025

接下来可以看看一道tp的ctf题目

1
2
3
4
5
6
7
import random
with open("flag.txt","rb") as f:
    flag=f.read()
for i in range(2**64):
    print(random.getrandbits(32)+flag[random.getrandbits(32)%len(flag)])
    input()

这里把最后8bit加上了另一个随机数进行了混淆，不过经过测试，加法可能会影响20bit以上甚至更高。不过在10000bit数据以内大多出现20bit以上数据概率比较小。我们这里设置混淆的是末尾24bit并且是间隔获得信息，获取第一个信息后，后一个随机数就被混淆掉了。我们这里可以获得2700组随机的数据

构建恢复state的代码如下

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
from Crypto.Util.number import *
from random import *
from tqdm import *
n=2500
D=[]
print(len(D))
rng=Random()

def getRows(rng):
    row=[]
    for i in range(n):
        row+=list(map(int, (bin(rng.getrandbits(32)>>24)[2:].zfill(8))))
        _=rng.getrandbits(32)
    return row
M=[]

for i in tqdm(range(19968)):
    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)) 
    M.append(getRows(rng))
M=Matrix(GF(2),M)
print("行数:", M.nrows())
print("列数:", M.ncols())
r=[]
for i in range(n):
    r+=list(map(int, (bin(D[i]>>24)[2:].zfill(8))))
r=vector(GF(2),r)
y=M.solve_left(r)
G=[]
for i in range(624):
    C=0
    for j in range(32):
        C<<=1
        C|=int(y[32*i+j])
    G.append(C)#将恢复的state转化为每个元素32bit的列表
import random
RNG1 = random.Random()
for i in range(624):
    G[i]=int(G[i])
RNG1.setstate((int(3),tuple(G+[int(624)]),None))#将state作为随机数生成器的初始状态，测试生成的state是否正确

恢复出状态后我们可以得到未被混淆的随机数。我们先可以利用爆破求出flag的长度和flag

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
G=[]
for i in range(624):
    G[i]=int(G[i])
RNG1 = Random()
RNG1.setstate((int(3),tuple(G+[int(624)]),None))

flag=[]
ind=[]
for i in range(len(X)):
    a = RNG1.getrandbits(32)
    b = RNG1.getrandbits(32)
    f = X[i]-a
    if 0 < f < 255 :
        flag.append(f)
        ind.append(b)
for fl in range(len(set(flag)),50):#fl是flag的长度，利用爆破判断长度是否符合条件
    TN = [i%fl for i in ind]
    TF = []
    for i in range(fl):
        t = TN.index(i)
        TF.append(flag[t])
    TF =  bytes(TF)
    if TF[:6]==(b'TPCTF{'):
       print(TF)
crypto
#crypto
梅森旋转生成随机数
http://example.com/2025/03/30/逆向梅森算法/
Author
fox
Posted on
March 30, 2025
Licensed under