Crypto
May 8, 2025

PRNG-MT19937分析

随机数生成过程

  1. 从种子生成初始状态数组
  2. 更新状态以获得新的内部状态
  3. 当需要生成一个数字时，返回一个经过回火的状态值(tempered state value)
  4. 当所有状态值都被使用完后，返回步骤 2

数字生成过程

给出底层数字生成的代码实现genrand_uint32 函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
/* generates a random number on [0,0xffffffff]-interval */
static uint32_t
genrand_uint32(RandomObject *self)
{
    uint32_t y;
    //...
    uint32_t *mt;

    mt = self->state;    
    if (self->index >= N) {
        // state update
    }

    y = mt[self->index++];
    y ^= (y >> 11);
    y ^= (y << 7) & 0x9d2c5680U;
    y ^= (y << 15) & 0xefc60000U;
    y ^= (y >> 18);
    return y;
}

变量y是最终返回的随机数结果，mt 是指向状态数组的指针，保存了当前生成器的状态（通常是长度为 624 的数组）

我们看一下其中的一些代码

1
mt = self->state;

这里将指向状态数组的指针赋给 mt，这是梅森旋转算法的核心状态数组

代码中if判断的那一步先跳过，后续会介绍状态使用完时应该如何进行状态数组的更新，这里只要先知道624个随机数结果输出后需要对其进行更新

1
y = mt[self->index++]

这里从状态数组取出下一个状态值，并自增index

1
2
3
4
y ^= (y >> 11);
y ^= (y << 7) & 0x9d2c5680U;
y ^= (y << 15) & 0xefc60000U;
y ^= (y >> 18);

这几行代码进行了一系列按位异或和位移操作，也就是上面所提到的tempering回火操作，用来提升输出的随机性，MT19937可被攻击的关键就是回火操作可逆，可利用以下python代码实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
def unshiftRight(x, shift):
    res = x
    for i in range(32):
        res = x ^ res >> shift
    return res

def unshiftLeft(x, shift, mask):
    res = x
    for i in range(32):
        res = x ^ (res << shift & mask)
    return res

def untemper(v):
    v = unshiftRight(v, 18)
    v = unshiftLeft(v, 15, 0xefc60000)
    v = unshiftLeft(v, 7, 0x9d2c5680)
    v = unshiftRight(v, 11)
    return v

我们详细的讲一下恢复的过程，这里我们举一个例子：

1
2
3
x = 0b10101010  # 原始值
x >> 3 = 0b00010101
x ^ (x >> 3) = y = 0b10111111

从上面的右移操作可知，高位的值被传递到了位移后的低位，在传递的过程中空出的高位用0来填充，那么我们根据异或的性质：

可知原始和移位的异或结果的高位(位数与移位数量一致)与原始一致，现在我们知道tempering操作的右移式子如下所示：

那么 的高位和原始的高位一致，我们现在只有，目的是恢复原始的，还以上面的移位三位为例，我们有如下表格

位索引 y[i] x[i] x[i+3] y[i] = x[i] ^ x[i+3] 说明
0–2 y[i] x[i] ? y[i] = x[i] x和y高位一致
3 y[3] x[3] x[0] x[3] = y[3] ^ x[0] 因为x右移了三位，最高位落到了移位后的第4位
4 y[4] x[4] x[1] x[4] = y[4] ^ x[1] x第二位落到了移位后的第5位
5 y[5] x[5] x[2] x[5] = y[5] ^ x[2] x第三位落到了移位后的第6位
6 y[6] x[6] x[3] x[6] = y[6] ^ x[3] 依赖之前恢复的 x[3]
7 y[7] x[7] x[4] x[7] = y[7] ^ x[4] 依赖 x[4]

由上面的图表可知，每次恢复的每一位均依赖已知的 高位和已知的 相异或得到，因此对于整个tempering操作，无论是左移还是右移均可被我们恢复，因此它是可逆的

利用获取到的 624 个连续的 32 位输出我们可对PRNG进行预测，因为我们已经还原了整个PRNG的内部状态

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
import random

# 逆向 tempering过程
def unshiftRight(x, shift):
    res = x
    for i in range(32):
        res = x ^ res >> shift
    return res

def unshiftLeft(x, shift, mask):
    res = x
    for i in range(32):
        res = x ^ (res << shift & mask)
    return res

def untemper(v):
    v = unshiftRight(v, 18)
    v = unshiftLeft(v, 15, 0xefc60000)
    v = unshiftLeft(v, 7, 0x9d2c5680)
    v = unshiftRight(v, 11)
    return v


myseed = int(b'Chromos2me'.hex(), 16)
print(myseed)
random.seed(myseed)

# 一些未知的输出
for _ in range(myseed):
    random.getrandbits(32)

# 捕获 624 个连续的 32 位输出
state = [untemper(random.getrandbits(32)) for _ in range(624)]

print("正常情况的输出 :")

print(random.getrandbits(32))
print(random.random())
print(random.randbytes(4).hex())
print(random.randrange(1, 100000))

print("\n预测输出 :")

# 利用逆向 tempering过程的结果设置random内部状态
random.setstate((3, tuple(state + [624]), None))

print(random.getrandbits(32))
print(random.random())
print(random.randbytes(4).hex())
print(random.randrange(1, 100000))

微信截图_20250503181701

内部状态更新过程

状态更新过程也在genrand_uint32 函数中，每624个输出执行一次更新，还有一种情况是在种子初始化之后第一次生成随机数时，上面的代码中我们并未给出具体执行更新操作的代码，详细代码如下所示：

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
#define N 624
#define M 397    
#define MATRIX_A 0x9908b0dfU    /* constant vector a */
#define UPPER_MASK 0x80000000U  /* most significant w-r bits */
#define LOWER_MASK 0x7fffffffU  /* least significant r bits */
/* generates a random number on [0,0xffffffff]-interval */
static uint32_t
genrand_uint32(RandomObject *self)
{
    uint32_t y;
    static const uint32_t mag01[2] = {0x0U, MATRIX_A};
    /* mag01[x] = x * MATRIX_A  for x=0,1 */
    uint32_t *mt;

    mt = self->state;
    if (self->index >= N) { /* generate N words at one time */
        int kk;

        for (kk=0;kk<N-M;kk++) {
            y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
            mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1U];
        }
        for (;kk<N-1;kk++) {
            y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
            mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1U];
        }
        y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
        mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1U];

        self->index = 0;
    }

    // tempering
}

看起来在处理内部状态更新的时候比较复杂，先简单介绍一下twist主要过程：

我们将(y >> 1) ^ mag01[y & 0x1U]这一部分称为 twist 过程，我们将其变为双变量函数的形式以便于我们理解：

注：x表示mt[kk]，y表示mt[kk+1]

接下来令初始状态为，当前状态为 ，我们可从初始状态推出完整的内部更新过程：

注意到中

很明显 主要依赖于 ，如果我们知道了第 228 次生成的输出伪随机数，我们实际上可以：

我们这里讨论python中的特别之处， 中的条件异或依赖于的 LSB

下面给出进行twist逆转的代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
def invertStep(si, si227):
    # S[i] ^ S[i-227] == (((I[i] & 0x80000000) | (I[i+1] & 0x7FFFFFFF)) >> 1) ^ (0x9908b0df if I[i+1] & 1 else 0)
    X = si ^ si227
    # we know the LSB of I[i+1] because MSB of 0x9908b0df is set, we can see if the XOR has been applied
    mti1 = (X & 0x80000000) >> 31
    if mti1:
        X ^= 0x9908b0df
    # undo shift right
    X <<= 1
    # now recover MSB of state I[i]
    mti = X & 0x80000000 
    # recover the rest of state I[i+1]
    mti1 += X & 0x7FFFFFFF
    return mti, mti1

random.getstate() 返回一个内部状态元组，其中 [1] 是 Mersenne Twister 的内部状态数组（共 624 + 1 个元素，第一个是 index）

调用一次 random.random()，这会：

  • 从当前状态数组中生成一个伪随机数（浮点数）
  • 并推进 Mersenne Twister 的状态
  • 如果这次调用跨越了 624 个数字的边界，会触发一次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
import random

def invertStep(si, si227):
    # S[i] ^ S[i-227] == (((I[i] & 0x80000000) | (I[i+1] & 0x7FFFFFFF)) >> 1) ^ (0x9908b0df if I[i+1] & 1 else 0)
    X = si ^ si227
    # we know the LSB of I[i+1] because MSB of 0x9908b0df is set, we can see if the XOR has been applied
    mti1 = (X & 0x80000000) >> 31
    if mti1:
        X ^= 0x9908b0df
    # undo shift right
    X <<= 1
    # now recover MSB of state I[i]
    mti = X & 0x80000000
    # recover the rest of state I[i+1]
    mti1 += X & 0x7FFFFFFF
    return mti, mti1

I = random.getstate()[1]
# this will force a state update
random.random()
S = random.getstate()[1]

for i in range(227, 240):
    Ii, Ii1 = invertStep(S[i], S[i-227])
    print(f"{Ii} == {I[i]&0x80000000}")
    print(f"{Ii1} == {I[i+1]&0x7FFFFFFF}")

Python 的内部种子生成机制

为了理解python下的种子机制，我们看一下Random类中的seed方法是如何实现的

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
class Random(_random.Random):
    def seed(self, a=None, version=2):
        """Initialize internal state from a seed.

        The only supported seed types are None, int, float,
        str, bytes, and bytearray.

        None or no argument seeds from current time or from an operating
        system specific randomness source if available.

        If *a* is an int, all bits are used.

        For version 2 (the default), all of the bits are used if *a* is a str,
        bytes, or bytearray.  For version 1 (provided for reproducing random
        sequences from older versions of Python), the algorithm for str and
        bytes generates a narrower range of seeds.

        """

        if version == 1 and isinstance(a, (str, bytes)):
            a = a.decode('latin-1') if isinstance(a, bytes) else a
            x = ord(a[0]) << 7 if a else 0
            for c in map(ord, a):
                x = ((1000003 * x) ^ c) & 0xFFFFFFFFFFFFFFFF
            x ^= len(a)
            a = -2 if x == -1 else x

        elif version == 2 and isinstance(a, (str, bytes, bytearray)):
            if isinstance(a, str):
                a = a.encode()
            a = int.from_bytes(a + _sha512(a).digest(), 'big')

        elif not isinstance(a, (type(None), int, float, str, bytes, bytearray)):
            _warn('Seeding based on hashing is deprecated\n'
                  'since Python 3.9 and will be removed in a subsequent '
                  'version. The only \n'
                  'supported seed types are: None, '
                  'int, float, str, bytes, and bytearray.',
                  DeprecationWarning, 2)

        super().seed(a)
        self.gauss_next = None

我们可以利用random.seed 函数手动为随机数生成器设置种子，同时这个函数在导入random 模块时会被自动调用，用于初始化PRNG。自动设置种子的情况下，参数a，即我们所说的种子为None。

seed函数内部在实现时存在两个版本的初始化过程，版本1兼容老版python，会将str 或 bytes哈希成一个64位整数，生成方式类似python的哈希算法，用于将字符串或者字节串形式的种子变成一个整数用于初始化。

版本2支持str，bytes，bytearray，首先将字符串编码成字节，拼接字节对象和自己的sha512哈希，然后转为一个大端整数作为种子。

注：当使用 int、float 或 None 类型进行种子初始化时，种子会被直接传递给内部的种子设置过程。

下面就可以举一个利用等价种子的例子来验证seed内部的运行过程和我们上面说到的一致（VERSION == 2）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import random
import hashlib

random.seed('my seed')
print(random.randbytes(32).hex())

sha_512 = hashlib.sha512('my seed'.encode()).digest().hex()
hex_str = 'my seed'.encode().hex()
# print(sha_512)  c046d02141256723ba893a44f1a54538aba62218287bc56767433568ca7548571fe1760ba924e2f8b1e9e8887814c985f3827c68fdb09e639f313ec1c83476be
# print(hex_str)  6d792073656564

# 等价种子 0x6d792073656564c046d02141256723ba893a44f1a54538aba62218287bc56767433568ca7548571fe1760ba924e2f8b1e9e8887814c985f3827c68fdb09e639f313ec1c83476be
equivalent_seed = 0x6d792073656564c046d02141256723ba893a44f1a54538aba62218287bc56767433568ca7548571fe1760ba924e2f8b1e9e8887814c985f3827c68fdb09e639f313ec1c83476be
random.seed(equivalent_seed)
print(random.randbytes(32).hex())

微信截图_20250504170914

实际上在Random类底层调用_randommodule.c的 random_seed 函数实现Random.seed

默认情况使用None进行种子生成，将会执行下面的代码：

1
2
3
4
5
6
7
8
9
10
11
12
if (arg == NULL || arg == Py_None) {
    if (random_seed_urandom(self) < 0) {
        PyErr_Clear();

        /* Reading system entropy failed, fall back on the worst entropy:
            use the current time and process identifier. */
        if (random_seed_time_pid(self) < 0) {
            return -1;
        }
    }
    return 0;
}

PRNG优先使用操作系统提供的加密安全伪随机数生成器（CSPRNG）进行种子初始化，这些来源的随机数足够不可预测，可以用于密码学用途，比如密钥、Token生成等；如果系统不支持这些加密安全接口，或者接口调用失败（例如沙盒限制、权限问题），Python 会降级使用较弱的熵源，如当前时间和PID

利用CSPRNG

系统特定的种子初始化过程由以下函数负责处理：

1
2
3
4
5
6
7
8
9
10
11
static int
random_seed_urandom(RandomObject *self)
{
    uint32_t key[N];

    if (_PyOS_URandomNonblock(key, sizeof(key)) < 0) {
        return -1;
    }
    init_by_array(self, key, Py_ARRAY_LENGTH(key));
    return 0;
}

关于熵源的详细信息如下摘自 pyurandom 函数的文档字符串：

Used sources of entropy ordered by preference, preferred source first:

  • BCryptGenRandom() on Windows
  • getrandom() function (ex: Linux and Solaris): call py_getrandom()
  • getentropy() function (ex: OpenBSD): call py_getentropy()
  • /dev/urandom device

系统 CSPRNG 用于生成 624 个 32 位整数，足以填充 MT PRNG 的整个状态，但这只是初始化工作的一部分，剩余工作由init_by_array完成

利用时间+PID

降级熵源

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
static int
random_seed_time_pid(RandomObject *self)
{
    PyTime_t now;
    if (PyTime_Time(&now) < 0) {
        return -1;
    }

    uint32_t key[5];
    key[0] = (uint32_t)(now & 0xffffffffU);
    key[1] = (uint32_t)(now >> 32);

#if defined(MS_WINDOWS) && !defined(MS_WINDOWS_DESKTOP) && !defined(MS_WINDOWS_SYSTEM)
    key[2] = (uint32_t)GetCurrentProcessId();
#elif defined(HAVE_GETPID)
    key[2] = (uint32_t)getpid();
#else
    key[2] = 0;
#endif

    if (PyTime_Monotonic(&now) < 0) {
        return -1;
    }
    key[3] = (uint32_t)(now & 0xffffffffU);
    key[4] = (uint32_t)(now >> 32);

    init_by_array(self, key, Py_ARRAY_LENGTH(key));
    return 0;
}

此过程生成四个 32 位整数，这些整数再次用于使用 init_by_array 函数初始化 MT PRNG 的内部状态。（生成数组较小，理想情况下只有128位熵，较差情况取决于当前时间和PID是否被成功猜测出来）

直接利用数字

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
/* This algorithm relies on the number being unsigned.
 * So: if the arg is a PyLong, use its absolute value.
 * Otherwise use its hash value, cast to unsigned.
 */
if (PyLong_CheckExact(arg)) {
    n = PyNumber_Absolute(arg);
} else if (PyLong_Check(arg)) {
    /* Calling int.__abs__() prevents calling arg.__abs__(), which might
        return an invalid value. See issue #31478. */
    _randomstate *state = _randomstate_type(Py_TYPE(self));
    n = PyObject_CallOneArg(state->Long___abs__, arg);
}
else {
    Py_hash_t hash = PyObject_Hash(arg);
    if (hash == -1)
        goto Done;
    n = PyLong_FromSize_t((size_t)hash);
}
if (n == NULL)
    goto Done;

/* Now split n into 32-bit chunks, from the right. */
bits = _PyLong_NumBits(n);
if (bits == (size_t)-1 && PyErr_Occurred())
    goto Done;

/* Figure out how many 32-bit chunks this gives us. */
keyused = bits == 0 ? 1 : (bits - 1) / 32 + 1;

/* Convert seed to byte sequence. */
key = (uint32_t *)PyMem_Malloc((size_t)4 * keyused);
if (key == NULL) {
    PyErr_NoMemory();
    goto Done;
}
res = _PyLong_AsByteArray((PyLongObject *)n,
                            (unsigned char *)key, keyused * 4,
                            PY_LITTLE_ENDIAN,
                            0, /* unsigned */
                            1); /* with exceptions */
if (res == -1) {
    goto Done;
}

init_by_array(self, key, keyused);

代码的大致逻辑是如果种子是一个正整数（无论长度多长），它会被分割为若干个32位的块以构造关键数组；如果是负整数，则使用其绝对值；如果是浮点数，则使用其hash值

利用原文章的例子介绍一下每段代码中出现的key数组的构造过程

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
import random

def split_into_32bit_chunks(n):
    chunks = []
    while n > 0:
        chunks.append(n & 0xFFFFFFFF)
        n >>= 32
    if not chunks:
        chunks = [0]
    return chunks

# random.seed(-0x909192939495969798)
key1 = [hex(chunk) for chunk in  split_into_32bit_chunks(abs(-0x909192939495969798))]
print(key1)
# key1 = [0x95969798, 0x91929394, 0x90]

# random.seed(0x909192939495969798)
key2 = [hex(chunk) for chunk in split_into_32bit_chunks(0x909192939495969798)]
print(key2)
# key2 = [0x95969798, 0x91929394, 0x90]

random.seed(hash(871628728198271972.8621))
key3 = [chunk for chunk in split_into_32bit_chunks(hash(871628728198271972.8621))]
print(key3)
# key4 = [3494417408, 202941877]

random.seed("my seed")
equivalent_seed = 0x6d792073656564c046d02141256723ba893a44f1a54538aba62218287bc56767433568ca7548571fe1760ba924e2f8b1e9e8887814c985f3827c68fdb09e639f313ec1c83476be
key4 = [chunk for chunk in split_into_32bit_chunks(equivalent_seed)]
print(key4)
# key = [3358881470, 2670804673, 4256210531, 4085415016, 2014628229, 2984896648, 2837766904, 534869515, 3396683863, 1732457832, 679200103, 2879791640, 4054140216, 3129555524, 1092970275, 3225866273, 1936024932, 7174432]

之后key数组（即种子数组）仍然和之前一样传入 init_by_array 函数

现在我们已经了解了不同类型的种子是如何去构造key数组的，接下来介绍如何通过这个可变长数组计算mt的初始状态，我们记key数组为，长度为

初始状态的构建

init_by_array函数如下所示：

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
/* initialize by an array with array-length */
/* init_key is the array for initializing keys */
/* key_length is its length */
static void
init_by_array(RandomObject *self, uint32_t init_key[], size_t key_length)
{
    size_t i, j, k, l;       /* was signed in the original code. RDH 12/16/2002 */
    uint32_t *mt;

    mt = self->state;
    init_genrand(self, 19650218U);
    i=1; j=0;
    k = (N>key_length ? N : key_length);

    for (; k; k--) {
        mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525U))
                 + init_key[j] + (uint32_t)j; /* non linear */
        i++; j++;
        if (i>=N) { mt[0] = mt[N-1]; i=1; }
        if (j>=key_length) j=0;
    }

    for (k=N-1; k; k--) {
        mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941U))
                 - (uint32_t)i; /* non linear */
        i++;
        if (i>=N) { mt[0] = mt[N-1]; i=1; }
    }

    mt[0] = 0x80000000U; /* MSB is 1; assuring non-zero initial array */
}

它用于通过用户提供的种子数组 init_key[] 来初始化状态数组 mt，从而确定随机数序列的起点，函数需要传入三个参数：

接下来介绍init_by_array函数的逻辑：

init_genrand 函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
/* initializes mt[N] with a seed */
static void
init_genrand(RandomObject *self, uint32_t s)
{
int mti;
uint32_t *mt;

mt = self->state;
mt[0]= s;
for (mti=1; mti<N; mti++) {
    mt[mti] =
    (1812433253U * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
    /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
    /* In the previous versions, MSBs of the seed affect   */
    /* only MSBs of the array mt[].                                */
    /* 2002/01/09 modified by Makoto Matsumoto                     */
}
self->index = mti;
return;
}

我们不去关注函数内部的具体实现，我们只记初始化完成的mt为

设A=1664525，函数定义如下：

第一个循环执行如下的操作：

==注：新中间状态 要通过预初始化后的mt得到，两式出现的原因是下面的代码==

1
2
3
4
5
6
7
for (; k; k--) {
        mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525U))
                 + init_key[j] + (uint32_t)j; /* non linear */
        i++; j++;
        if (i>=N) { mt[0] = mt[N-1]; i=1; }    // the code for (2)(3)
        if (j>=key_length) j=0;				   
    }

首先，i 和 j 同时自增，当i>=N时，将mt[0]的值修改为mt[623]的值；当j>=key_length时，将j置为0，因为循环的条件是max(N,key_length)，因此if的两个判断均能被执行，当执行完 if (i>=N) { mt[0] = mt[N-1]; i=1; }之后，状态被改变， 被重置到1，也即下一次更新的状态是（这里的更新是在已被更新过的基础上），那么在这一次循环中被改变了两次，在一般情况下循环到这里就结束了。

那么这里i和j的作用就很明显了，i 控制状态数组的索引，如果 i >= N，回绕到头部；j 控制 key 的索引，如果 j >= key_length，回绕 key 数组

这里设B = 1566083941，第二个循环执行了下面的操作：

这里的第(7)个式子从i=2开始，因为i=1在(6)中被使用，剩下的操作基本上从代码中都能看懂

从少量输出中恢复种子

假如我们现在知道了初始状态的两个连续值 $Ji和J{i-1}$可从(4)式写出

注：因为的原因 要从3开始，状态的值固定

同时假设我们知道种子数组 的长度 ，那么我们就能恢复索引值 （个人感觉能够恢复 的原因是，在if (j>=key_length) j=0;这行代码里， 和数组的长度存在关联），之后就能恢复 的实际值。这样的推导是合理的，我们在上面说过 key 数组的 取决于种子的类型或其长度，即使我们不知道 的实际长度，猜测它对于我们来说也是非常容易的。

注：通常情况下

实现上面操作的代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
def init_genrand(seed):
        MT = [0] * 624
        MT[0] = seed & 0xffffffff
        for i in range(1, 623+1): # loop over each element
            MT[i] = ((0x6c078965 * (MT[i-1] ^ (MT[i-1] >> 30))) + i) & 0xffffffff
        return MT

def recover_kj_from_Ji(ji, ji1, i):
    # ji => J[i]
    # ji1 => J[i-1]
    const = init_genrand(19650218)
    key = ji - (const[i] ^ ((ji1 ^ (ji1 >> 30))*1664525))
    key &= 0xffffffff
    # return K[j] + j
    return key

到目前为止，我们实现了中间状态 到种子 的逆过程，接下来寻找初始状态 到中间状态 的过程

若已知两个连续的初始状态值 $Ii和I{i-1}，由(7)$可得

注：由于，从3开始上式才成立

上述操作的代码如下所示:

1
2
3
4
5
6
7
def recover_Ji_from_Ii(Ii, Ii1, i):
    # Ii => I[i]
    # Ii1 => I[i-1]
    ji = (Ii + i) ^ ((Ii1 ^ (Ii1 >> 30)) * 1566083941)
    ji &= 0xffffffff
    # return J[i]
    return ji

到目前为止我们证明了两个连续的 值可以恢复单个 值，又证明了两个连续的 值足以恢复 值，如原文配图所示

微信截图_20250507221432

上图的过程用代码写出即

1
2
3
4
5
6
7
8
9
def recover_Kj_from_Ii(Ii, Ii1, Ii2, i):
    # Ii => I[i]
    # Ii1 => I[i-1]
    # Ii2 => I[i-2]
    # Ji => J[i]
    # Ji1 => J[i-1]
    Ji = recover_Ji_from_Ii(Ii, Ii1, i)
    Ji1 = recover_Ji_from_Ii(Ii1, Ii2, i-1)
    return recover_kj_from_Ji(Ji, Ji1, i)

上述的恢复过程在下才成立

完整POC

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
import random

def init_genrand(seed):
    MT = [0] * 624
    MT[0] = seed & 0xffffffff
    for i in range(1, 623 + 1):  # loop over each element
        MT[i] = ((0x6c078965 * (MT[i - 1] ^ (MT[i - 1] >> 30))) + i) & 0xffffffff
    return MT


def recover_kj_from_Ji(ji, ji1, i):
    # ji => J[i]
    # ji1 => J[i-1]
    const = init_genrand(19650218)
    key = ji - (const[i] ^ ((ji1 ^ (ji1 >> 30)) * 1664525))
    key &= 0xffffffff
    # return K[j] + j
    return key

def recover_Ji_from_Ii(Ii, Ii1, i):
    # Ii => I[i]
    # Ii1 => I[i-1]
    ji = (Ii + i) ^ ((Ii1 ^ (Ii1 >> 30)) * 1566083941)
    ji &= 0xffffffff
    # return J[i]
    return ji

def recover_Kj_from_Ii(Ii, Ii1, Ii2, i):
    # Ii => I[i]
    # Ii1 => I[i-1]
    # Ii2 => I[i-2]
    # Ji => J[i]
    # Ji1 => J[i-1]
    Ji = recover_Ji_from_Ii(Ii, Ii1, i)
    Ji1 = recover_Ji_from_Ii(Ii1, Ii2, i-1)
    return recover_kj_from_Ji(Ji, Ji1, i)

random.seed(12345)
# K = [12345]
# k = 1

I = random.getstate()[1]

for i in range(4, 624):
    print(i, recover_Kj_from_Ii(I[i], I[i-1], I[i-2], i))

微信截图_20250507221936

在上面我们提到知道当前状态的一对值 $(Si,\ S{i+227})可以恢复I{i+228}的低位以及I{i+227}的最高位，那么现在如果知道了另一对状态的值(S{i+1},\ S{i+228})，我们能够恢复I{i+229}的低位以及I{i+228}的最高位，那么此时我们就能获得I_{i+228}$ 的精确数值

那么现在继续添加一组会发生什么？添加 $(S{i+2},\ S{i+229})，我们恢复I{i+230}的低位以及I{i+229}的最高位，我们又知道了一个当前状态I{i+229}的精确值，那么现在我们有状态I{i+228}，I{i+229}，I{i+230}的低位，我们就可以说我们拥有了个初始状态连续值，我们又可以利用这种方式恢复两个Kj的潜在值了（因为I{i+230}$最高位值不确定）

==注： 的MSB只影响 的MSB==

如果Python的PRNG使用一个 32 位整数进行手动播种，那么只需要前 624 个输出中的任意 6 个，就可以恢复出该种子。

下面是验证代码：

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
import random

def unshiftRight(x, shift):
    res = x
    for i in range(32):
        res = x ^ res >> shift
    return res

def unshiftLeft(x, shift, mask):
    res = x
    for i in range(32):
        res = x ^ (res << shift & mask)
    return res

def untemper(v):
    v = unshiftRight(v, 18)
    v = unshiftLeft(v, 15, 0xefc60000)
    v = unshiftLeft(v, 7, 0x9d2c5680)
    v = unshiftRight(v, 11)
    return v

def invertStep(si, si227):
    # S[i] ^ S[i-227] == (((I[i] & 0x80000000) | (I[i+1] & 0x7FFFFFFF)) >> 1) ^ (0x9908b0df if I[i+1] & 1 else 0)
    X = si ^ si227
    # we know the LSB of I[i+1] because MSB of 0x9908b0df is set, we can see if the XOR has been applied
    mti1 = (X & 0x80000000) >> 31
    if mti1:
        X ^= 0x9908b0df
    # undo shift right
    X <<= 1
    # now recover MSB of state I[i]
    mti = X & 0x80000000
    # recover the rest of state I[i+1]
    mti1 += X & 0x7FFFFFFF
    return mti, mti1

def init_genrand(seed):
    MT = [0] * 624
    MT[0] = seed & 0xffffffff
    for i in range(1, 623 + 1):  # loop over each element
        MT[i] = ((0x6c078965 * (MT[i - 1] ^ (MT[i - 1] >> 30))) + i) & 0xffffffff
    return MT


def recover_kj_from_Ji(ji, ji1, i):
    # ji => J[i]
    # ji1 => J[i-1]
    const = init_genrand(19650218)
    key = ji - (const[i] ^ ((ji1 ^ (ji1 >> 30)) * 1664525))
    key &= 0xffffffff
    # return K[j] + j
    return key

def recover_Ji_from_Ii(Ii, Ii1, i):
    # Ii => I[i]
    # Ii1 => I[i-1]
    ji = (Ii + i) ^ ((Ii1 ^ (Ii1 >> 30)) * 1566083941)
    ji &= 0xffffffff
    # return J[i]
    return ji

def recover_Kj_from_Ii(Ii, Ii1, Ii2, i):
    # Ii => I[i]
    # Ii1 => I[i-1]
    # Ii2 => I[i-2]
    # Ji => J[i]
    # Ji1 => J[i-1]
    Ji = recover_Ji_from_Ii(Ii, Ii1, i)
    Ji1 = recover_Ji_from_Ii(Ii1, Ii2, i-1)
    return recover_kj_from_Ji(Ji, Ji1, i)

for i in range(10):
    random.seed(i)
    # k = 1
    # K = [i]

    S = [untemper(random.getrandbits(32)) for _ in range(624)]

    I_227_, I_228 = invertStep(S[0], S[227])
    I_228_, I_229 = invertStep(S[1], S[228])
    I_229_, I_230 = invertStep(S[2], S[229])

    I_228 += I_228_
    I_229 += I_229_

    # two possibilities for I_230
    seed1 = recover_Kj_from_Ii(I_230, I_229, I_228, 230)
    seed2 = recover_Kj_from_Ii(I_230+0x80000000, I_229, I_228, 230)
    # only the MSB differs
    print(hex(seed1), hex(seed2))

微信截图_20250508170455

上面我们所讨论的是恢复单个32位整数种子的推导过程，那么对于这种情况来说种子数组的长度是 ，如果现在种子数组的长度是2呢？回想上面我们已经获得3个连续的初始状态，如果现在再利用一对当前状态反推出一组初始状态的 MSB 和低31位，我们就能继续恢复整个数组，无论恢复多少个初始状态只有的MSB会缺失，因此这个状态都会存在两个值，这也就导致存在两个可能的值

微信截图_20250508172831

原文章给出了在每种情形下恢复种子所需的最小输出数量，如下所示：

微信截图_20250508173112

下面我们将探究一个特殊的情况，，这种情况下我们只需要624个输出即可完成种子的恢复，因为随着关系对的不断加入，我们可以在他们之间建立更多的关系，从而降低恢复种子所需要的输出。

从完整状态恢复种子

状态更新操作可逆，知道完整的状态就足以从任何状态倒回 初始状态

状态更新过程中仅使用了的MSB，其他位均未被使用，因此我们无法去恢复这些确切的位，因为这个原因回溯后的状态在第一个值上将会有所不同

下面代码基于实现了倒带过程

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
def rewindState(state):
    prev = [0]*624
    # copy to not modify input array
    s = state[:]
    I, I0 = invertStep(s[623], s[396])
    prev[623] += I
    # update state 0
    # this does nothing when working with a known full state, but is important we rewinding more than 1 time
    s[0] = (s[0]&0x80000000) + I0
    for i in range(227, 623):
        I, I1 = invertStep(s[i], s[i-227])
        prev[i] += I
        prev[i+1] += I1
    for i in range(227):
        I, I1 = invertStep(s[i], prev[i+397])
        prev[i] += I
        prev[i+1] += I1
    # The LSBs of prev[0] do not matter, they are 0 here
    return prev

上面的代码实际上利用所有生成的当前状态去反推初始状态，逻辑还是很简单的

下面是POC

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
import random

def unshiftRight(x, shift):
    res = x
    for i in range(32):
        res = x ^ res >> shift
    return res

def unshiftLeft(x, shift, mask):
    res = x
    for i in range(32):
        res = x ^ (res << shift & mask)
    return res

def untemper(v):
    v = unshiftRight(v, 18)
    v = unshiftLeft(v, 15, 0xefc60000)
    v = unshiftLeft(v, 7, 0x9d2c5680)
    v = unshiftRight(v, 11)
    return v

def invertStep(si, si227):
    # S[i] ^ S[i-227] == (((I[i] & 0x80000000) | (I[i+1] & 0x7FFFFFFF)) >> 1) ^ (0x9908b0df if I[i+1] & 1 else 0)
    X = si ^ si227
    # we know the LSB of I[i+1] because MSB of 0x9908b0df is set, we can see if the XOR has been applied
    mti1 = (X & 0x80000000) >> 31
    if mti1:
        X ^= 0x9908b0df
    # undo shift right
    X <<= 1
    # now recover MSB of state I[i]
    mti = X & 0x80000000
    # recover the rest of state I[i+1]
    mti1 += X & 0x7FFFFFFF
    return mti, mti1

def rewindState(state):
    prev = [0]*624
    # copy to not modify input array
    s = state[:]
    I, I0 = invertStep(s[623], s[396])
    prev[623] += I
    # update state 0
    # this does nothing when working with a known full state, but is important we rewinding more than 1 time
    s[0] = (s[0]&0x80000000) + I0
    for i in range(227, 623):
        I, I1 = invertStep(s[i], s[i-227])
        prev[i] += I
        prev[i+1] += I1
    for i in range(227):
        I, I1 = invertStep(s[i], prev[i+397])
        prev[i] += I
        prev[i+1] += I1
    # The LSBs of prev[0] do not matter, they are 0 here
    return prev

I = list(random.getstate()[1][:-1])

S1 = [untemper(random.getrandbits(32)) for _ in range(624)]
S2 = [untemper(random.getrandbits(32)) for _ in range(624)]
S3 = [untemper(random.getrandbits(32)) for _ in range(624)]

# rewind once
I_ = rewindState(S1)
S2_ = rewindState(S3)
S1_ = rewindState(S2)

print(I_ == I)
print(S1_[1:] == S1[1:])
print(S2_[1:] == S2[1:])

# rewind multiple times
I_ = rewindState(rewindState(rewindState(S3)))
print(I_ == I)
print(I_[:5])
print(I[:5])

微信截图_20250508182451

执行完上面的POC之后我们恢复了完整的初始状态，接下来就是恢复种子，我们将利用下面的代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
def seedArrayFromState(s, subtractIndices=True):
    s_ = [0]*624
    for i in range(623, 2, -1):
        s_[i] = recover_Ji_from_Ii(s[i], s[i-1], i)
    s_[0]=s_[623]
    s_[1]=recover_Ji_from_Ii(s[1], s[623],  1)
    s_[2]=recover_Ji_from_Ii(s[2], s_[1], 2)
    seed = [0]*624
    for i in range(623, 2, -1):
        seed[i-1] = recover_kj_from_Ji(s_[i], s_[i-1], i)
    # system overdefined for seed[0,1,623]
    seed[0] = 0
    # thus s1 = (const[1] ^ ((const[0] ^ (const[0] >> 30))*1664525))
    s1_old = ((2194844435 ^ ((19650218 ^ (19650218 >> 30))*1664525))) & 0xffffffff
    seed[1] = recover_kj_from_Ji(s_[2], s1_old, 2)
    seed[623] = (s_[1] - (s1_old ^ ((s_[0] ^ (s_[0] >> 30))*1664525))) & 0xffffffff
    # subtract the j indices
    if subtractIndices:
        seed = [(2**32+e-i)%2**32 for i,e in enumerate(seed)]
    return seed

但是 init_by_array 函数执行种子初始化操作的方式导致了不同的init_key数组可能会生成出相同的初始化状态，即种子与初始化状态不是一一对应的，这就造成了种子冲突的可能性

这里原文章中提出了一种解决方案，我们可以固定为一个数，这样使得我们反推init_key时简化解空间，实现种子恢复的唯一性，这里其实也好理解就相当于一个不定方程，我们存在非常多种可能解，我们这里固定，那么 ，即我们从大量的可能中选择一个作为我们的结果，这里原文章选择固定

实现这一过程的POC如下：

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
import random

def seedArrayToInt(s):
    seed = 0
    for e in s[::-1]:
        seed += e
        seed <<= 32
    return seed >> 32

# random is not manually seeded so it uses 624 random values

S = [untemper(random.getrandbits(32)) for _ in range(624)]
I = rewindState(S)

print("Normal run :")

print(random.getrandbits(32))
print(random.random())
print(random.randbytes(4).hex())
print(random.randrange(1, 100000))

print("\nReseeded run :")

seed_array = seedArrayFromState(I)
seed = seedArrayToInt(seed_array)

# the recovered seed is very big. Too big to be printed in decimal
# print(hex(seed))
random.seed(seed)

S_ = [untemper(random.getrandbits(32)) for _ in range(624)]

assert(S_ == S)

print(random.getrandbits(32))
print(random.random())
print(random.randbytes(4).hex())
print(random.randrange(1, 100000))

上述的假设建立在已知的 624 个连续输出要么是直接在种子初始化操作之后生成的，要么与其之间相隔了 624 的倍数个输出

微信截图_20250508194110

About this Post

This post is written by Chromos2me, licensed under CC BY-NC 4.0.

#MT19937#PRNG