Mersenne Twisterの出力を推測してみる

「Z3Pyでglibc rand(3)の出力を推測してみる」および「Z3Pyでxorshift128+の出力を推測してみる」では、Z3Pyを使って疑似乱数生成アルゴリズムのいくつかの出力を観測することで、後続の出力を推測した。 ここでは、Pythonのrandomモジュールなどで利用されている疑似乱数生成アルゴリズムMersenne Twisterの出力を推測してみる。

環境

Ubuntu 14.04.3 LTS 64bit版

$ uname -a
Linux vm-ubuntu64 3.19.0-25-generic #26~14.04.1-Ubuntu SMP Fri Jul 24 21:16:20 UTC 2015 x86_64 x86_64 x86_64 GNU/Linux

$ lsb_release -a
No LSB modules are available.
Distributor ID: Ubuntu
Description:    Ubuntu 14.04.3 LTS
Release:        14.04
Codename:       trusty

$ gcc --version
gcc (Ubuntu 4.8.4-2ubuntu1~14.04) 4.8.4

Mersenne Twisterソースコードを読んでみる

Mersenne Twisterソースコードから、乱数生成部分を抜き出したものを次に示す。

#define N 624

(snip)

static unsigned long mt[N]; /* the array for the state vector  */
static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */

(snip)

/* generates a random number on [0,0xffffffff]-interval */
unsigned long genrand_int32(void)
{
    unsigned long y;
    static unsigned long mag01[2]={0x0UL, MATRIX_A};
    /* mag01[x] = x * MATRIX_A  for x=0,1 */

    if (mti >= N) { /* generate N words at one time */
        int kk;

        if (mti == N+1)   /* if init_genrand() has not been called, */
            init_genrand(5489UL); /* a default initial seed is used */

        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 & 0x1UL];
        }
        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 & 0x1UL];
        }
        y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
        mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];

        mti = 0;
    }

    y = mt[mti++];

    /* Tempering */
    y ^= (y >> 11);
    y ^= (y << 7) & 0x9d2c5680UL;
    y ^= (y << 15) & 0xefc60000UL;
    y ^= (y >> 18);

    return y;
}

上のコードから、Mersenne Twisterは内部状態として624個の32bit整数(unsigned long型)を持つことがわかる。 また、624回出力するたびに内部状態が全更新され、出力はそれぞれの内部状態に対してTemperingと呼ばれる処理を加えたものであることがわかる。

Pythonのrandomモジュールの出力を推測してみる

Pythonのrandomモジュールは、Mersenne Twisterを用いて疑似乱数の生成を行う。

上の結果をもとに、624個の出力からMersenne Twisterの内部状態を復元し、出力を推測するスクリプトを書いてみる。

# predict_mt.py
import random

def untemper(x):
    x = unBitshiftRightXor(x, 18)
    x = unBitshiftLeftXor(x, 15, 0xefc60000)
    x = unBitshiftLeftXor(x, 7, 0x9d2c5680)
    x = unBitshiftRightXor(x, 11)
    return x

def unBitshiftRightXor(x, shift):
    i = 1
    y = x
    while i * shift < 32:
        z = y >> shift
        y = x ^ z
        i += 1
    return y

def unBitshiftLeftXor(x, shift, mask):
    i = 1
    y = x
    while i * shift < 32:
        z = y << shift
        y = x ^ (z & mask)
        i += 1
    return y

values1 = [random.getrandbits(32) for i in xrange(624)]
values2 = [random.getrandbits(32) for i in xrange(624)]

mt_state = tuple([untemper(x) for x in values1] + [624])
random.setstate((3, mt_state, None))

predicted = [random.getrandbits(32) for i in xrange(624)]
print predicted == values2

ここで、untemper関数はTempering処理の逆演算を行う関数である。 また、randomモジュールが用いるMersenne Twisterの内部状態はrandom.setstate()関数でセットできる。 セットすべき値の詳細は、実際のrandom.getstate()関数の出力およびソースコードの次の箇所から調べられる。

    def getstate(self):
        """Return internal state; can be passed to setstate() later."""
        return self.VERSION, super(Random, self).getstate(), self.gauss_next
static PyObject *
random_getstate(RandomObject *self)
{
    PyObject *state;
    PyObject *element;
    int i;

    state = PyTuple_New(N+1);
    if (state == NULL)
        return NULL;
    for (i=0; i<N ; i++) {
        element = PyLong_FromUnsignedLong(self->state[i]);
        if (element == NULL)
            goto Fail;
        PyTuple_SET_ITEM(state, i, element);
    }
    element = PyLong_FromLong((long)(self->index));
    if (element == NULL)
        goto Fail;
    PyTuple_SET_ITEM(state, i, element);
    return state;

Fail:
    Py_DECREF(state);
    return NULL;
}

実際にスクリプトを実行してみると次のようになる。

$ python predict_mt.py
True

実行結果より、Mersenne Twisterの出力を624個観測することにより、後続の疑似乱数列が推測できていることが確認できる。

関連リンク