跳至主要內容

MT19937分析

原创Xenny约 3309 字大约 14 分钟密码学PRNGMT19937

MT19937分析

前言

  • 本内容是2022发在cryptography-wiki上的文章了,现在转回到blog上

正文

  • MT19937即梅森旋转算法(Mersenne twister)由松本眞(日语:松本真)和西村拓士在1997年开发,基于二进制有限域F2\mathbb{F}_2上的矩阵线性递归,可以快速产生高质量的伪随机数。

    该算法的周期为21993712^{19937}-1,故名为MT19937。该算法具有以下优点

    1. 周期非常长,为21993712^{19937}-1
    2. 1k6231\le k \le 623都满足kk-分布
    3. 能比硬件实现的方法更快产生随机数
  • kk-分布:一个周期为PPww位整数的伪随机数序列xix_i,如果下列成立则称其vv-比特精度的kk-分布成立。

    truncv(x)\mathrm{trunc_v(x)}表示由xx的前vv位形成的数,并考虑PPkkvv位向量。

    (truncv(xi),truncv(xi+1),,truncv(xi+k1)) (0i<P) (\mathrm{trunc_v(x_i)},\mathrm{trunc_v(x_{i+1})},\dots,\mathrm{trunc_v(x_{i+k-1})})\ (0\le i< P)

    然后,2kv2^{kv}个组合中每一个都在一个周期内出现次数相同(全0组合出现次数较少除外)。

代码实现

  • MT19937算法可分为三个部分

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

    其中32位的MT19937 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
    

安全分析

逆向extract_number函数

  • 对于extract_number函数,可以发现就是将状态中的数据进行位运算后给出。

    显然我们可以发现这里的四次位运算类似,我们以y = y ^ y << 7 & 2636928640为例来进行逆向处理。

    为了方便描述,我们将代码改为y = x ^ ((x << 7) & 2636928640),其中yy可以看作已知,xx未知。

    xaaaabbbbccccddddeeeeffffgggghhhhx<<7bccccddddeeeeffffgggghhhh0000000&263692864010011101001011000101011010000000=ynaannnbnccncnnddenenfnnfnggghhhh \begin{aligned} x &\rightarrow \mathtt{aaaabbbbccccddddee{\color{blue}eeffffg}\color{red}{ggghhhh}}\\ &\oplus\\ x<<7 &\rightarrow \mathtt{bccccddddee{\color{blue}eeffffg}{\color{red}{ggghhhh}}0000000}\\ &\&\\ 2636928640 &\rightarrow \mathtt{10011101001011000101011010000000}\\ &=\\ y &\rightarrow \mathtt{naannnbnccncnndden{\color{purple}enfnnfn}\color{red}{ggghhhh}} \end{aligned}

    这里为了方便编写便使用相同的字母作为一段,实际上其值可能不同。

    显然我们可以发现yy的低7位就是xx的低7位和2636928640的低7位异或的结果。且2636928640低7位为0,所以yy的低7位就是xx的低7位。即我们能够得到ggghhhh\mathtt{\color{red}{ggghhhh}}的结果,基于此,我们便可以向前恢复得到eeffffg\mathtt{\color{blue}eeffffg}

    eeffffg=(ggghhhh&0101101)enfnnfn \mathtt{\color{blue}eeffffg} = (\mathtt{\color{red}{ggghhhh}} \&\mathtt{0101101}) \oplus \mathtt{\color{purple}enfnnfn}

    同理可以向前恢复得到xx的所有值。完整Python代码实现为

    o = 3989032602
    
    def inverse_right_mask(res, shift, mask=0xffffffff, bits=32):
        tmp = res
        for i in range(bits // shift):
            tmp = res ^ tmp >> shift & mask
        return tmp
    
    def inverse_left_mask(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_mask(y,18)
        y = inverse_left_mask(y,15,4022730752)
        y = inverse_left_mask(y,7,2636928640)
        y = inverse_right_mask(y,11)
        return y&0xffffffff
    
    print(recover(extract_number(o)) == o)
    

预测随机数

  • 当我们能够逆向extract_number中的数据时,这也意味着我们能够提取得到state中的原始数据,同时可以发现后续的随机数完全依赖于上一轮的state,所以如果我们能够得到某一轮的全部state数据,便可以向后调用twist来预测随机数。

    此类题型在CTF中经常出现,部分题型没有明确的给出随机数信息,但可以泄漏的其他信息得到随机数,解题时可以考虑加以注意是否能够直接或间接的得到题目中的随机数信息。

逆向twist函数

  • 在上文中我们提到如果得到了某一轮state的全部信息便可以向后预测随机数,那么如果我们需要向前恢复随机数,则需要对twist函数进行逆向。

    yaaaaaaaaaaaaaaaaaaaaaaaaaaaaaauxy10aaaaaaaaaaaaaaaaaaaaaaaaaaaaaau0x9908b0df10011001000010001011000011011111=bbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbsi+397=cccccccccccccccccccccccccccccccv=out=dddddddddddddddddddddddddddddddd \begin{aligned} y \rightarrow& \mathtt{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaa{\color{blue}{u}}\color{red}{x}}\\ y \gg 1 \rightarrow& \mathtt{0aaaaaaaaaaaaaaaaaaaaaaaaaaaaaa{\color{blue}{u}}}\\ \oplus&\\ \mathrm{0x9908b0df} \rightarrow& \mathtt{10011001000010001011000011011111}\\ =& \mathtt{bbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbb}\\ \oplus&\\ s_{i+397} =& \mathtt{ccccccccccccccccccccccccccccccc{\color{blue}{v}}}\\ =&\\ out=& \mathtt{ddddddddddddddddddddddddddddddd{\color{blue}{d}}} \end{aligned}

    其中y = state[i] & 0x80000000 | state[i + 1] & 0x7fffffff,我们可以发现新的值只和si,si+1,si+397s_i,s_{i+1},s_{i+397}有关,在上面的两次异或中,第一次异或和yy的奇偶性来确定

    由于y1y\gg 1的最高位必定为0,此时我们可以列出下列关系表,

    01
    0NY
    1YN

    其中行为最终输出out的最高位值,列为si+397s_{i+397}的最高位值,单元值代表是否发生第一次异或。由此通过判断是否发生异或,可以得到yy的最低位值。

    此时按照extract_number中一样的方法,我们有

    u=dvx ,or u=dvx1 \mathtt{\color{blue}{u}} = {\color{blue}{d}}\oplus {\color{blue}{v}} \oplus {\color{red}{x}}\ ,\mathrm{or}\ \mathtt{\color{blue}{u}} = {\color{blue}{d}}\oplus {\color{blue}{v}} \oplus {\color{red}{x}}\oplus 1

    同理,我们便可以向前恢复得到yy的所有值。

    同时yy的值即为sis_i的最高位和si+1s_{i+1}的第[2, 32]位。同时我们可以对si+1s_{i+1}si+2s_{i+2}生成的yy值做同样操作得到同样si+1s_{i+1}的最高位,对si1s_{i-1}sis_i生成的yy值做同样操作得到sis_i的第[2, 32]位,至此原state恢复完成。

    Python代码实现为:

    def inv_twist(state):
        high = 0x80000000
        low = 0x7fffffff
        mask = 0x9908b0df
        
        def recover(i):
            y = state[i + 624] ^ state[i + 397]
            if y & high == high: # 异或了常数
                y ^= mask
                y <<= 1
                y |= 1
            else: # 没有异或常数
                y <<= 1
            return y
        
        for i in range(len(state)-625, -1, -1):
            # 得到s_i的最高位
            state[i] = recover(i) & high
            # 对s_{i-1}做同样操作得到2-32位
            state[i] |= recover(i-1) & low
        return state
    

逆向init函数

  • 我们可以发现第一轮的初始状态是通过seed生成的。

    关键操作为

    self.mt[i] = _int32(1812433253 * (self.mt[i - 1] ^ self.mt[i - 1] >> 30) + i)
    

    显然这里的加法和乘法都存在逆运算,而中间的mti1(mti130)mt_{i-1}\oplus (mt_{i-1} \gg 30)可以通过上文所述方法进行逆运算。

    from gmpy2 import invert
    
    def _int32(x):
        return int(0xFFFFFFFF & x)
    
    def init(seed):
        mt = [0] * 624
        mt[0] = seed
        for i in range(1, 624):
            mt[i] = _int32(1812433253 * (mt[i - 1] ^ mt[i - 1] >> 30) + i)
        return mt
    
    seed = 3989032602
    
    def invert_right(res,shift):
        tmp = res
        for i in range(32//shift):
            res = tmp^res>>shift
        return _int32(res)
    
    def recover(last):
        n = 1<<32
        inv = invert(1812433253,n)
        for i in range(623,0,-1):
            last = ((last-i)*inv)%n
            last = invert_right(last,30)
        return last
    
    state = init(seed)
    
    print(recover(state[-1]) == seed)
    

扩展

  • 在上文中我们介绍了对于MT19937各个函数的逆向分析以及算法实现,但其中我们针对的都是以32bit为一组,并且能够获取到足够多连续随机数值的情况,那么如果我们得到的并不是连续的随机数此时我们该做何分析。

    以Python中的random实现的getrandbits为例。

    // getrandbits(k) -> x.  Generates an int with k random bits.
    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;
    }
    

    可以发现,当参数为0时,返回0,当参数小于32时,生成一个32位的随机数取其高位,当参数大于32时,生成多个随机数进行拼接。我们不妨考虑一些最极端的情况

情况一

  • 考虑我们不能获取连续的32bit随机数,只能隔1个取1个,考虑此时恢复随机数。

    显然这对我们来说没有太多影响,因为si+624s_{i+624}只和si,si+1,si+397s_i,s_{i+1},s_{i+397}有关。假设此时ii为偶数,且我们只取第奇数位数,则我们可以得到si+1,si+397s_{i+1}, s_{i+397}

    yaaaaaaaaaaaaaaaaaaaaaaaaaaaaaauxy10aaaaaaaaaaaaaaaaaaaaaaaaaaaaaau0x9908b0df10011001000010001011000011011111=bbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbsi+397=vcccccccccccccccccccccccccccccc=out=ddddddddddddddddddddddddddddddd \begin{aligned} y \rightarrow& \mathtt{{\color{red}{a}}aaaaaaaaaaaaaaaaaaaaaaaaaaaaa{\color{blue}{u}}x}\\ y \gg 1 \rightarrow& \mathtt{0{\color{red}{a}}aaaaaaaaaaaaaaaaaaaaaaaaaaaaa{\color{blue}{u}}}\\ \oplus&\\ \mathrm{0x9908b0df} \rightarrow& \mathtt{10011001000010001011000011011111}\\ =& \mathtt{bbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbb}\\ \oplus&\\ s_{i+397} =& \mathtt{{\color{blue}{v}}{\color{blue}{c}}ccccccccccccccccccccccccccccc}\\ =&\\ out=& \mathtt{{\color{blue}{d}}{\color{purple}{d}}ddddddddddddddddddddddddddddd} \end{aligned}

    此时我们可以得到yy的[2-32]位,则可以由此得到si+624s_{i+624}的第1和第[3-32]位。同时又有第2位的约束

    {ac=d \begin{cases} \mathtt{{\color{red}{a}}}\oplus\mathtt{{\color{blue}{c}}} = \mathtt{{\color{purple}{d}}}\\ \end{cases}

    显然这是一个较强的约束,因为si+624s_{i+624}中的错误会扩散至si,si+624397,si+6241,...s_{i},s_{i+624-397},s_{i+624-1},...

    而此时我们能够所有第奇数个数的值和第偶数个数的1,332{1,3-32}位值。故我们可以得到si+624s_{i+624}以及sis_{i}的全部信息。

情况二

  • 考虑我们现在只能获取所有第偶数个数的情况,此时显然我们只能够获取si+397s_{i+397}的第2位。但是si+397s_{i+397}的第二位配合si+397+1s_{i+397+1},我们便能知道

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

    的全部信息,此时再异或si+397+397s_{i+397+397}便可以得到si+397+624s_{i+397+624}的所有信息。最终通过此信息我们便能向前回推出sj (1j)s_j\ (1 \le j)的所有信息。

情况三

  • 现在我们考虑一种更加极端的情况,假设我们只能获取每个随机数的第kk位比特信息。我们是否还能利用扩散得到足够多的约束信息。

    这里我们可以回顾之前的内容,我们的本质是通过已知信息推出与之直接关联的未知信息,再利用已知信息和未知信息之间的约束条件进行扩散,因为我们有足够多的信息,我们总是能够利用足够多的约束来求解所有信息。

    对于这种情况,我们来考虑MT19937在F2\mathbb{F}_2上的表现形式。

    设state[i]的二进制表示为:

    x={x0,x1,,x30,x31} \mathbf{x} = \{x_0,x_1,\dots,x_{30},x_{31}\}

    输出output的二进制表示位:

    z={z0,z1,,z30,z31} \mathbf{z} = \{z_0,z_1,\dots,z_{30},z_{31}\}

    我们有

    z0=x0x4x7x15z1=x1x5x16z2=x2x6x13x17x24z3=x3x10z4=x0x4x8x11x15x19x26z5=x1x5x9x12x20z6=x6x10x17x21x28z7=x3x7x11x14x18x22x29z8=x8x12x23z9=x9x13x20x24x31z10=x6x10x17z11=x0x11z12=x1x8x12x19z13=x2x9x13x17x20x28z14=x3x14x18x29z15=x4x15z16=x5x16z17=x6x13x17x24z18=x0x4x15x18z19=x1x5x8x15x16x19x26z20=x2x6x9x13x17x20x24z21=x3x17x21x28z22=x0x4x8x15x18x19x22x26x29z23=x1x5x9x20x23z24=x6x10x13x17x20x21x24x28x31z25=x3x7x11x18x22x25x29z26=x8x12x15x23x26z27=x9x13x16x20x24x27x31z28=x6x10x28z29=x0x11x18x29z30=x1x8x12x30z31=x2x9x13x17x28x31 \begin{aligned} z_0=&x_0 \oplus x_4 \oplus x_7 \oplus x_{15}\\ z_1=&x_1 \oplus x_5 \oplus x_{16}\\ z_2=&x_2 \oplus x_6 \oplus x_{13} \oplus x_{17} \oplus x_{24}\\ z_3=&x_3 \oplus x_{10}\\ z_4=&x_0 \oplus x_4 \oplus x_8 \oplus x_{11} \oplus x_{15} \oplus x_{19} \oplus x_{26}\\ z_5=&x_1 \oplus x_5 \oplus x_9 \oplus x_{12} \oplus x_{20}\\ z_6=&x_6 \oplus x_{10} \oplus x_{17} \oplus x_{21} \oplus x_{28}\\ z_7=&x_3 \oplus x_7 \oplus x_{11} \oplus x_{14} \oplus x_{18} \oplus x_{22} \oplus x_{29}\\ z_8=&x_8 \oplus x_{12} \oplus x_{23}\\ z_9=&x_9 \oplus x_{13} \oplus x_{20} \oplus x_{24} \oplus x_{31}\\ z_{10}=&x_6 \oplus x_{10} \oplus x_{17}\\ z_{11}=&x_0 \oplus x_{11}\\ z_{12}=&x_1 \oplus x_8 \oplus x_{12} \oplus x_{19}\\ z_{13}=&x_2 \oplus x_9 \oplus x_{13} \oplus x_{17} \oplus x_{20} \oplus x_{28}\\ z_{14}=&x_3 \oplus x_{14} \oplus x_{18} \oplus x_{29}\\ z_{15}=&x_4 \oplus x_{15}\\ z_{16}=&x_5 \oplus x_{16}\\ z_{17}=&x_6 \oplus x_{13} \oplus x_{17} \oplus x_{24}\\ z_{18}=&x_0 \oplus x_4 \oplus x_{15} \oplus x_{18}\\ z_{19}=&x_1 \oplus x_5 \oplus x_8 \oplus x_{15} \oplus x_{16} \oplus x_{19} \oplus x_{26}\\ z_{20}=&x_2 \oplus x_6 \oplus x_9 \oplus x_{13} \oplus x_{17} \oplus x_{20} \oplus x_{24}\\ z_{21}=&x_3 \oplus x_{17} \oplus x_{21} \oplus x_{28}\\ z_{22}=&x_0 \oplus x_4 \oplus x_8 \oplus x_{15} \oplus x_{18} \oplus x_{19} \oplus x_{22} \oplus x_{26} \oplus x_{29}\\ z_{23}=&x_1 \oplus x_5 \oplus x_9 \oplus x_{20} \oplus x_{23}\\ z_{24}=&x_6 \oplus x_{10} \oplus x_{13} \oplus x_{17} \oplus x_{20} \oplus x_{21} \oplus x_{24} \oplus x_{28} \oplus x_{31}\\ z_{25}=&x_3 \oplus x_7 \oplus x_{11} \oplus x_{18} \oplus x_{22} \oplus x_{25} \oplus x_{29}\\ z_{26}=&x_8 \oplus x_{12} \oplus x_{15} \oplus x_{23} \oplus x_{26}\\ z_{27}=&x_9 \oplus x_{13} \oplus x_{16} \oplus x_{20} \oplus x_{24} \oplus x_{27} \oplus x_{31}\\ z_{28}=&x_6 \oplus x_{10} \oplus x_{28}\\ z_{29}=&x_0 \oplus x_{11} \oplus x_{18} \oplus x_{29}\\ z_{30}=&x_1 \oplus x_8 \oplus x_{12} \oplus x_{30}\\ z_{31}=&x_2 \oplus x_9 \oplus x_{13} \oplus x_{17} \oplus x_{28} \oplus x_{31} \end{aligned}

    显然我们可以构造一个F232\mathbb{F}_2^{32}方阵TT满足

    xT=z\mathbf{x}\cdot T = \mathbf{z}

    对于如何构造这个方阵有很多方式,在[1]中使用了一种选择乘数的方式得到TT的每一行。

    现在我们把这种表示扩展到更高位空间,对于MT19937中的某一轮,其state可以表示为

    x={x0,x1,,x19966,x19967} \mathbf{x} = \{x_0,x_1,\dots,x_{19966},x_{19967}\}

    此时如果我们能获取每个随机数的第kk位比特,则我们可以得到

    s={zi+k}s = \{z_{i+k}\}

    若我们能够获取足够多的数据,则可以构造

    s={s0,s1,s19966,s19967} \mathbf{s} = \{s_0,s_1,s_{19966},s_{19967}\}

    同时构造

    xT=s \mathbf{x}\cdot T = \mathbf{s}

    最后在GF2上求解这个代数方程即可。

参考

上次编辑于:
贡献者: X3NNY