MT19937 是一种使用梅森旋转算法的伪随机数生成器。基于二进制有限域 $\mathbb{F}_2$ 上的矩阵线性递归，可以快速产生高质量且均匀分布的伪随机数。MT19937 的名称来源于它的周期长度为 $2^{19937}-1$

原理

（此处仅讨论 32 位的 MT19937）

MT19937 的状态数组是一个长度为 624 的整数数组，称为 MT。每一个元素都是一个 32 位的整数。它使用一个整数作为种子来初始化状态数组。每次调用随机数生成函数时，它会从状态数组中提取一个整数，并使用梅森旋转算法来更新状态数组。

MT19937 主要有三个部分组成

1. 初始化
2. 旋转状态
3. 提取伪随机数

初始化。MT19937 初始化的时候会使用给定的数字作为种子，或者使用系统时间戳等随机数作为种子。然后会初始化一个长度为 624 的状态数组 `MT`，并从种子开始通过特定算法来填充这个数组，作为初始状态。

旋转。旋转是指将状态数组中的每个元素进行线性变换，生成新的状态数组。这个过程是通过对当前状态数组中的每个元素进行位运算和异或运算来实现的。旋转的目的是为了增加随机数的均匀性和不可预测性。

提取伪随机数。这一步会从状态数组中依次提取一个整数，并对其进行位运算和异或运算，生成的数即为输出的伪随机数。当所有数组全被遍历过之后，就会对状态数组再次进行一次旋转，重新生成新的状态数组。

一个简易的 Python 代码实现如下

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):
        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

⚠ 需要注意的是，这里只是实现了简易的 MT19937 算法，虽然 Python 的 random 库中使用的也是同样的 MT19937 算法来产生伪随机数，但是它对 seed 的传入经过了两步处理 `init_genrand`和 `init_by_array`，因此和上面的实现是有区别的，也就是说两者产生的状态矩阵和伪随机数是不一样的。后面如果有机会，可以详细分析一下之间的差别。

<!– https://github.com/python/cpython/blob/3.11/Modules/_randommodule.c –>

可能会注意到，经过 MT19937 算法产生的伪随机数都是一个 32 位的整数，但事实上可以使用 randm.getrandbits()来获取任意位数的随机数。

于是翻一翻 Random 库的调用链，最后发现 getrandbits()通过调用 _randommodule.c文件中的 _random_Random_getrandbits_imp函数来获取随机数。

random 的实现可以查看源码https://github.com/python/cpython/blob/3.11/Modules/\_randommodule.c

static PyObject *
_random_Random_getrandbits_impl(RandomObject *self, int k)
/*[clinic end generated code: output=b402f82a2158887f input=8c0e6396dd176fc0]*/
{
    int i, words;
    uint32_t r;
    uint32_t *wordarray;
    PyObject *result;

    if (k < 0) {
        PyErr_SetString(PyExc_ValueError,
                        "number of bits must be non-negative");
        return NULL;
    }

    if (k == 0)
        return PyLong_FromLong(0);

    if (k <= 32)  /* Fast path */
        return PyLong_FromUnsignedLong(genrand_uint32(self) >> (32 - k));

    words = (k - 1) / 32 + 1;
    wordarray = (uint32_t *)PyMem_Malloc(words * 4);
    if (wordarray == NULL) {
        PyErr_NoMemory();
        return NULL;
    }

    /* Fill-out bits of long integer, by 32-bit words, from least significant
       to most significant. */
#if PY_LITTLE_ENDIAN
    for (i = 0; i < words; i++, k -= 32)
#else
    for (i = words - 1; i >= 0; i--, k -= 32)
#endif
    {
        r = genrand_uint32(self);
        if (k < 32)
            r >>= (32 - k);  /* Drop least significant bits */
        wordarray[i] = r;
    }

    result = _PyLong_FromByteArray((unsigned char *)wordarray, words * 4,
                                   PY_LITTLE_ENDIAN, 0 /* unsigned */);
    PyMem_Free(wordarray);
    return result;
}

这里实际上是当位数小于 32 时，会截去低位部分，直接取高位作为输出，例如 0x12345678 的 32 位整数，经过 getrandbits(16)后会得到 0x1234，经过 getrandbits(8)后会得到 0x12。

当位数大于 32 时，会使用多个伪随机数进行拼接，例如生成 64 位随机数时，会先按顺序产生两个 32 位整数，0x12345678 和 0x87654321，然后按先生成的排在低位的顺序，拼接成 0x8765432112345678。

连续的 624 个 32 位伪随机数

这里要求这 624 个数是在同一个状态下产生的，以便可以直接对上述的三部分操作进行逆向

1. extract_number 逆向

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)

发现对于 y 的四次操作逻辑上很类似，我们以 y = y ^ ((y << 7) & 2636928640) 为例来进行逆向处理。

首先能注意到 2636928640的二进制表示为 0b10011101001011000101011010000000，那么就意味着最终输出的结果低 7 位和原来 y 的低 7 位是相同的。于是就能从低位倒推得到完整的 y

def invert(res, shift, right=True, mask=0xffffffff, bits=32):
    tmp = res
    if right:
        for i in range(bits // shift):
            tmp = res ^ tmp >> shift & mask
        return tmp
    else:
        for i in range(bits // shift):
            tmp = res ^ tmp << shift & mask
        return tmp

def inv_extract_number(y):
    y = invert(y,18,True)
    y = invert(y,15,False,4022730752)
    y = invert(y,7,False,2636928640)
    y = invert(y,11,True)
    return _int32(y)

于是我们就能得到整个完整的状态数组

2. 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

先关注旋转的后几步操作，能注意到由于 y>>1的最高位一定为 0，所以最终 self.mt[i]的最高位一定由 self.mt[(i + 397) % 624]或 self.mt[(i + 397) % 624] ^ 0x9908b0df控制，所以可以判断出是否经历了异或 0x9908b0df 操作。然后由于是否异或操作同时受最低位控制，那么逆向的时候即可，通过是否异或来恢复因为就右移而丢失的最低位。于是我们就得到了 y

然后分析旋转的第一步，y 是由 self.mt[i]的最高位和 self.mt[(i + 1) % 624]的除最高位部分组合得到的。所以我们只要计算 self.mt[i]和 self.mt[(i - 1) % 624]两个位置的 y 就能得到 self.mt[i]的值了

def inv_twist(state):
    high = 0x80000000
    low = 0x7fffffff
    mask = 0x9908b0df

    def _recover(i):
        y = state[i] ^ state[(i + 397) % 624]
        if y & high == high:
            y ^= mask
            y <<= 1
            y |= 1
        else:
            y <<= 1
        return y

    for i in range(len(state)-625, -1, -1):
        state[i] = _recover(i) & high
        state[i] |= _recover(i-1) & low
    return state

3. init 逆向

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)

这个部分的逆向就和 extract_number部分的比较类似了，加法乘法都有对应的逆运算，只需要逆向 self.mt[i - 1] ^ self.mt[i - 1] >> 30即可，前面提过如何逆向

def inv_init(last):
    n = 1<<32
    inv = pow(1812433253,-1,n)
    for i in range(623,0,-1):
        last = ((last-i)*inv)%n
        last = invert(last,30)
    return last

ps. 虽然这里 init 的逆向有点鸡肋，因为很显然在最开始有 seed 初始化状态数组的时候，第一个就是 seed

不连续的 624*32 位伪随机数

这种情况是只要能获得总计 624*32 位（19968 位）的伪随机数就能进行逆向了，不需要是连续的 32 位即可进行随机数预测。

首先已经提到过 MT19937 是线性的，线性便能想到矩阵。此外描述一个 MT19937 的状态需要 19968 位来确定，于是可以思考有这样一个列向量$s_{19968\times 1}$来表示初始状态和一个递推矩阵$T_{19968\times 19968}$。

$$T_{19968\times 19968}\cdot s_{19968\times 1} = b_{19968\times 1}$$

用$b$表示对初始状态经过线性变化后得到的伪随机数序列。也就是说，不管$b$中的伪随机数是怎么排列的，一定是由$s$经过$T$的线性变化得到的，那么我们就需要找到这个递推矩阵$T$，然后就能通过线性代数的方法来求解出$s$了。从而恢复出 MT19937 的状态，也就是能实现伪随机数的预测。

获取 T

可以确定的一件事是，$T$表示的是一次的线性操作变换，这个变换只和想要的伪随机数的位置有关，和初始状态是什么无关。

那么是否可以通过构造初始状态来得到$T$呢？显然是可以的，这种方法称为黑盒攻击。将$s$设置为全 0，然后令第 i 个位置为 1，经过相同的递推矩阵$T$线性变换后得到的$b$就是$T$的第 i 行的行向量。于是遍历$s$所有位置就可以构造出$T$了。

例如我能得到 2496 个连续的 random.getrandbits(8)，可以这样预测随机数

# Sage
from random import Random

def construct_a_row(RNG):
    row = []
    for _ in range(19968//8):
        bits = RNG.getrandbits(8)
        row += list(map(int, bin(bits)[2:].zfill(2)))
    return row

T = []
for i in trange(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 = Random()
    RNG.setstate((3,tuple(state+[624]),None))
    T.append(construct_a_row(RNG))

T = Matrix(GF(2),T)
b = [int(j) for i in ouput for j in bin(i)[2:].zfill(8)]
B = vector(GF(2),b)
s = L.solve_left(R)

state = []
for i in range(624):
    state.append(int("".join(list(map(str,s)))[32*i:32*i+32],2))

RNG = Random()
RNG.setstate((3,tuple(state+[624]),None))

于是就得到了状态相同的 RNG，接下来就可以使用 RNG.getrandbits()来预测伪随机数了。

* random 库分析

如果你去尝试用上述给出的 MT19937 的实现来预测随机数，你会发现它和 Python 的 random 库产生的伪随机数是不同的。原因在于 random 库在初始化状态数组时会经过两步处理 init_genrand()和 init_by_array()，而我们上面给出的实现只是简单地使用了种子来初始化状态数组。

去翻查 random 库的调用链，就能发现 init_by_array()这个函数用来初始化状态数组，详细实现如下

static void
init_by_array(RandomObject *self, uint32_t init_key[], size_t key_length)
{
    size_t i, j, k;       /* 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_genrand(self, 19650218U);来预初始化。这个函数的实现如下

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);
    }
    self->index = mti;
    return;
}

这下发现 init_genrand()函数和上述的 __init__()函数是一样的，在 random 库中是使用固定的常量先预初始化了状态数组，然后再用 seed 通过 init_by_array()函数来进行初始化，使得状态数组分布更加均匀。

这里在之前 MT19937 的基础上用 python 实现一下 random 库的主要函数 getrandbits()

class Py_MT19937:
    def __init__(self, seed):
        seed = self._int32(seed)
        self.mt = [0] * 624
        self.mti = 0
        self._init_by_array([seed], 1)

    def _int32(self, x):
        return int(0xFFFFFFFF & x)

    def _init_genrand(self, seed):
        self.mt[0] = seed
        for i in range(1, 624):
            self.mt[i] = self._int32(1812433253 * (self.mt[i - 1] ^ self.mt[i - 1] >> 30) + i)

    def _init_by_array(self, init_key, key_length):
        self._init_genrand(19650218)
        i = 1
        j = 0
        k = max(624, key_length)
        for _ in range(k):
            self.mt[i] = self._int32((self.mt[i] ^ ((self._int32(self.mt[i - 1]) ^ (self._int32(self.mt[i - 1]) >> 30)) * 1664525)) + init_key[j] + j)
            i += 1
            j = (j + 1) % key_length
            if i >= 624:
                self.mt[0] = self.mt[623]
                i = 1

        for _ in range(623):
            self.mt[i] = self._int32((self.mt[i] ^ ((self._int32(self.mt[i - 1]) ^ (self._int32(self.mt[i - 1]) >> 30)) * 1566083941)) - i)
            i += 1
            if i >= 624:
                self.mt[0] = self.mt[623]
                i = 1
        self.mt[0] = 0x80000000

    def _twist(self):
        for i in range(624):
            y = self._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] ^= 0x9908B0DF

    def _extract_number(self, n=32):
        if self.mti%624 == 0:
            self._twist()
        y = self.mt[self.mti]
        y ^= (y >> 11)
        y ^= (y << 7) & 0x9D2C5680
        y ^= (y << 15) & 0xEFC60000
        y ^= (y >> 18)

        self.mti = (self.mti + 1) % 624
        return self._int32(y) >> (32-n)

    def getrandbits(self, n):
        res = 0
        for i in range(n//32):
            res += self._extract_number(32) << (32 * i)
        res += self._extract_number(n % 32)
        return res

<!– // TODO: 有空研究一下该初始化的逆向，如果可行的话就能求出seed了–>

* 尝试逆向

_init_by_array()函数主要分为两个部分，这两个部分基本都是对使用固定种子预初始化过的状态数组进行类似的变化操作，唯一区别在于第一部分使用了 seed 进行额外操作

仔细分析能发现所谓的变化操作是对当前值和上个值操作，上个值是已知的，于是就能逆向得到变化操作之前的原值。但注意一点，第二部分的 i 遍历是从 2-623-1 的顺序，而且最后还把状态数组的首位定为常量。那么必然状态数组经过第二步操作后会有前三个值无法还原。不过这不影响获取 seed。

然后尝试对第一步操作逆向，会发现原值其实是已知的，这就是预处理的状态数组，而预处理过程的 seed 是固定的，那么结合原值，当前值和上个值就能倒退出 seed 了。

def _int32(x):
    return int(0xFFFFFFFF & x)

def _re_init_by_array_part(index, mt, multiplier):
    return _int32((mt[index]+index) ^ (mt[index-1] ^ mt[index-1] >> 30) * multiplier)

def _init_genrand(seed,mt):
    mt[0] = seed
    for i in range(1, 624):
        mt[i] = _int32(1812433253 * (mt[i - 1] ^ mt[i - 1] >> 30) + i)

def re_init_by_array(mt=None):
    if mt is None:
        seed = random.randint(0, 2**32-1)
        RNG = random.Random(seed)
        _, mt, _ = RNG.getstate()

    tmp = [_re_init_by_array_part(i, mt[:-1], 1566083941) for i in [622,623]]

    original_mt = [0] * 624
    _init_genrand(19650218, original_mt)

    predict_seed = _int32(tmp[-1] - _int32((tmp[-2] ^ (tmp[-2] >> 30)) * 1664525 ^ original_mt[-1]))
    if "seed" in locals():
        print(predict_seed, seed)
        print(predict_seed == seed % 2**32)
    else:
        return predict_seed

⚠ 显然这里输出的预测种子一定是 32 位以内的数，所以只能保证得到的是种子的低 32 位，至于高位部分就无法恢复了。因为对于 32 位的 MT19937 来说 seed 真正起作用的是低 32 位部分，高位是无效的