前一篇文章分析了Python中随机算法的实现细节，本文就来对其进行逆向。

由前文所述，MT19937提取随机数可分为两部分：twist 、tempering

def extract_number(self):
    if self._mti >= self.N:
        for i in range(self.N):
            self.twist(i)
        self._mti = 0

    y = self._mt[self._mti]
    y = self.tempering(y)

    self._mti += 1

    return _int32(y)

那么，其逆向过程就先从termpering操作开始。

逆向 tempering

def tempering(y):
    y ^= (y >> 11)
    y ^= (y << 7) & 0x9d2c5680
    y ^= (y << 15) & 0xefc60000
    y ^= (y >> 18)
    return y

我们倒着一步一步分析，约定记号如下：

  • 「异或」运算记为$\oplus$，「与」运算记为$\wedge$
  • 每一步运算前的变量为$y$，得到的结果为$z$​
  • 记变量最高位的下标为0，第二高位的下标为1，以此类推
  • 变量从高位到低位的连续一段切片，以上下标标记，下标为起点，上标为终点。例如$y$的高18字节记为$y_0^{17}$

先看最后一步：y ^= (y >> 18)

我们知道$z$是32位整数，根据这个公式，显而易见的结论有：

于是：

注意到这个$z\to y$的公式与前面$y\to z$的在形式上一模一样，故这一步的逆向我们只需照抄正向：

y ^= (y >> 18)

再看倒数第二步：y ^= (y << 15) & 0xefc60000

记0xefc60000为$c$。

注意到bin(c) == 0b11101111110001100000000000000000，这个二进制数的低17位全为0。

故我们可以写出这一步的正向公式：

同理，容易写出逆公式：

发现形式也相同，故这一步也直接抄正向：

y ^= (y << 15) & 0xefc60000

接下来分析倒数第三步：y ^= (y << 7) & 0x9d2c5680

记0x9d2c5680为$d_1$。

这里不容易像前面那样直接写出逆公式，不过我们可以用类似于递归的方法来求解。

首先我们有：

因此：

将此表达式直接代入右边的$y$，得到：

记上式为

我们来计算$X$：

这里$d_2=(d_1\ll7)\wedge d_1=\text{0x94284000}$​

同理，我们可以不断将下式

代入到等号右侧的$y$并展开，我们会得到：

我们记右侧的异或项序列为${X_i}$，即

其中，

计算得：

由此可知，我们在展开到第五项时，彻底消去了等号右侧的$y$，因此：

至此，我们已经可以写出这一步的逆向代码：

y ^= ((y << 7) & 0x9d2c5680) ^ ((y << 14) & 0x94284000) ^ ((y << 21) & 0x14200000) ^ ((y << 28) & 0x10000000)

最后，逆向第一步：y ^= (y >> 11)

类似于上一步，我们可以不断右移再异或，直到右侧的$y$变成0：

注意，$y$是32位整数，右移33位就归零了，因此，第一步的逆向如下：

y ^= (y >> 11) ^ (y >> 22)

综合上述内容，我们可以写出untempering：

def untempering(y):
    y ^= (y >> 18)
    y ^= (y << 15) & 0xefc60000
    y ^= ((y << 7) & 0x9d2c5680) ^ ((y << 14) & 0x94284000) ^ ((y << 21) & 0x14200000) ^ ((y << 28) & 0x10000000)
    y ^= (y >> 11) ^ (y >> 22)
    return y

逆向 twist

def twist(self, i):
    y = (self._mt[i] & self.UPPER_MASK) | (self._mt[(i + 1) % self.N] & self.LOWER_MASK)
    self._mt[i] = y >> 1
    if y & 0x1 == 1:
        self._mt[i] ^= self.MATRIX_A
    self._mt[i] ^= self._mt[(i + self.M) % self.N]

逆向twist其实相当于恢复_mt[i]。

我们首先写出最后一步的逆向：

tmp = self._mt[i] ^ self._mt[(i + self.M) % self.N]

这里tmp的值应为上面twist函数中经过3、4、5三行之后的_mt[i]的值，那么如何判断在正向过程中是否进入了这个if分支？其实我们只要关注tmp的最高位（32位整数的意义下）即可。

如果未曾进入if分支，那么tmp的值为y>>1，是某个32位整数右移得来的，故此时tmp最高位必为0；反之，若进入了if分支，其还会异或一个MATRIX_A，而我们知道MATRIX_A的最高位为1，因此这时tmp最高位也一定会变成1。

再考虑到进入if语句的条件是y的最低位为1，因此我们根据tmp的最高位的值，其实已经可以推导出twist函数里的变量y的值了。

接下来的几行代码：

if tmp & self.UPPER_MASK == self.UPPER_MASK:
    tmp ^= self.MATRIX_A
    tmp <<= 1
    tmp |= 1
else:
    tmp <<= 1

此时tmp的值相当于twist函数里的变量y，我们看看它包含了哪些信息：

y = (self._mt[i] & self.UPPER_MASK) | (self._mt[(i + 1) % self.N] & self.LOWER_MASK)

这说明y的最高位是_mt[i]的最高位，y的后31位则是_mt[i+1]的后31位。

因此我们已经恢复出_mt[i]的最高位了，接下来只要恢复其后31位即可。

显然，想要恢复_mt[i]的后31位，只需将前面所有操作的下标减去1，即可在tmp的后31位得到_mt[i]的后31位。这一步非常巧妙。

于是，untwist的完整代码就呼之欲出了：

def untwist(self, i):
    tmp = self._mt[i] ^ self._mt[(i + self.M) % self.N]
    if tmp & self.UPPER_MASK == self.UPPER_MASK:
        tmp ^= self.MATRIX_A
        tmp <<= 1
        tmp |= 1
    else:
        tmp <<= 1
    res = tmp & self.UPPER_MASK
	
    # 进行与前面一模一样的操作，不过将下标减去了1
    tmp = self._mt[i - 1] ^ self._mt[(i + self.M - 1) % self.N]
    if tmp & self.UPPER_MASK == self.UPPER_MASK:
        tmp ^= self.MATRIX_A
        tmp <<= 1
        tmp |= 1
    else:
        tmp <<= 1
    res |= tmp & self.LOWER_MASK
    self._mt[i] = res

逆向 extract_number

def extract_number(self):
    if self._mti >= self.N:
        for i in range(self.N):
            self.twist(i)
        self._mti = 0

    y = self._mt[self._mti]
    y = self.tempering(y)

    self._mti += 1

    return _int32(y)

这就很容易了，可以直接写出：

def unextract_number(self):
    self._mti = (self._mti - 1) % self.N
    if self._mti == 0:
        for i in range(self.N - 1, -1, -1):
            self.untwist(i)
        self._mti = self.N

调用此函数即可将当前的随机数内部状态倒回去一次迭代。

以上，我们已基本实现了对MT19937算法的逆向，以此为基础，我们便有了预测Python随机数的能力，具体内容见下一篇文章。