📘

Mersenne Twister (MT19937) で未来と過去の乱数列を予測してみる【Python】

2023/09/04に公開

はじめに

Mersenne Twister (MT, メルセンヌ・ツイスタ) は擬似乱数を生成するアルゴリズムの一種です。その中でも、MT19937 という手法には次のような長所があります。

周期が $2^{19937}-1$ と非常に長い。
高次元（ $623$ 次元）でも均等に分布する。
擬似乱数の生成が高速。

このような特徴から、様々なプログラミング言語の標準的な擬似乱数生成器に Mersenne Twister が採用されています。たとえば、Python の random には MT19937 が採用されています。

しかし、Mersenne Twister は暗号論的擬似乱数ではありません^[1]。暗号論的擬似乱数とは、暗号論的に安全な擬似乱数のことを指し、具体的には予測不可能であることが求められます。つまり、Mersenne Twister は一定の条件下で予測可能な擬似乱数です。

この記事では、MT19937 によって出力された乱数列から、未来と過去の乱数列を予測します。より詳細には次の 2 ステップを取ります。

現在の内部状態を復元する。
過去の内部状態を復元する。

MT19937

数学的な背景はともかく、MT19937 のアルゴリズムを知るためには上記の C 言語実装を読むのが最速だと思います。ここでは、ソースコードを適宜抜粋しながら簡単に解説します^[2]。

ライブラリと定数

stdio.h を include していますが、これは main 関数でのサンプル用です。そのため、mt19937ar.c は単体で動作します^[3]。

#include <stdio.h>

次に、乱数生成に必要な定数が定義されています。N は内部状態を表す状態ベクトルのサイズ、他は内部状態の更新で使われる定数です。

#define N 624
#define M 397
#define MATRIX_A 0x9908b0dfUL   /* constant vector a */
#define UPPER_MASK 0x80000000UL /* most significant w-r bits */
#define LOWER_MASK 0x7fffffffUL /* least significant r bits */

状態ベクトル

内部状態を表す状態ベクトル mt は次のように定義されます。mti は次に生成される乱数を表す index です。mti が N+1 の場合、内部状態が初期化されていない、つまり後述の init_genrand が一度も呼び出されていないことを示します。

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 */

内部状態の初期化関数

このプログラムで定義されている初期化の関数は以下の通りです。

init_genrand
シード値を受け取り、内部状態を初期化する。初期化後の mti は N になります。
init_by_array
シード値を配列とその長さとして与えて、内部状態を初期化する。

乱数生成関数

このプログラムで定義されている乱数生成の関数は以下の通りです。この記事では genrand_int32 を解説します。

genrand_int32
32 ビットの符号なし整数を生成する。
genrand_int31
31 ビットの符号付き整数を生成する。
genrand_real1
$[0,1]$ の浮動小数点数を生成する。
genrand_real2
$[0,1)$ の浮動小数点数を生成する。
genrand_real3
$(0,1)$ の浮動小数点数を生成する。
genrand_res53
$[0,1)$ の浮動小数点数（53-bit resolution）を生成する。

genrand_int32

genrand_int32 は 32 ビットの符号なし整数を生成する関数です。

変数の定義

y は後で生成される乱数を格納、mag01 は内部状態を更新するために使います。

unsigned long y;
static unsigned long mag01[2]={0x0UL, MATRIX_A};
/* mag01[x] = x * MATRIX_A  for x=0,1 */

内部状態の更新

内部状態 mt は、初期化または N 個の 32 ビット整数が出力されて枯渇したときに更新されます。init_genrand が呼ばれていない場合（mti が N+1）には、デフォルトのシード値 5489 を使って init_genrand を呼び出します。更新後、mti は 0 にリセットされます。

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;
}

乱数の生成

内部状態 mt から値を取り出し、mti を更新します。

y = mt[mti++];

最後に、Tempering と呼ばれる操作が実行されます。これにより、生成された乱数が均等に分布するように調整されます。

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

return y;

必要な乱数列の生成

この記事では Python の random を用いて、MT19937 の乱数予測を行います。必要な乱数の長さは、内部状態に対応して $32 \times 624 = 19968$ ビットです。ここでは、32 ビットずつ生成しています^[4]。

import random

N = 624 # 状態ベクトルのサイズ
xs0 = [random.getrandbits(32) for _ in range(N)] # 検証用の乱数列（過去）
xs1 = [random.getrandbits(32) for _ in range(N)] # 復元に使用する乱数列
xs2 = [random.getrandbits(32) for _ in range(N)] # 検証用の乱数列（未来）

現在の内部状態の復元

現在の内部状態を復元するためには、各乱数に Tempering の逆操作をする必要があります。実装は上記のサイトから引用しています。

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

これらの関数を用いて Tempering の逆操作をし、random.setstate で内部状態をセットします。

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

random.setstate の引数として与えられるタプル（またはリスト）の要素は次の通りです。ドキュメントにほとんど記載されていないので、CPython の実装を参考にしました。

version
バージョンごとに処理を分岐させる。3 または 2。通常用途では 3 でよさそう。

In version 2, the state was saved as signed ints, which causes inconsistencies between 32/64-bit systems. The state is really unsigned 32-bit ints, so we convert negative ints from version 2 to positive longs for version 3.
internalstate
内部状態。mt19937ar.c における mt と mti のタプル。
gauss_next
ガウス分布の計算で使われるパラメータ。通常用途では None でよさそう。 $\mu + z \sigma$ の $z$ に対応する。

最後に、復元した内部状態から生成される乱数列 predicted が、検証用の乱数列（未来） xs2 と一致するか検証します。結果は無事一致しました。ちなみに、mti を 0 にすると predicted は xs1 と一致します。

predicted = [random.getrandbits(32) for _ in range(N)]
print(predicted == xs2)

True

まとめるとプログラムは次の通りです。

プログラム

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


N = 624 # 状態ベクトルのサイズ
xs0 = [random.getrandbits(32) for _ in range(N)] # 検証用の乱数列（過去）
xs1 = [random.getrandbits(32) for _ in range(N)] # 復元に使用する乱数列
xs2 = [random.getrandbits(32) for _ in range(N)] # 検証用の乱数列（未来）

mt_state = [untemper(x) for x in xs1]
random.setstate((3, tuple(mt_state + [N]), None))
# random.setstate((3, tuple(mt_state + [0]), None))

predicted = [random.getrandbits(32) for _ in range(N)]
print(predicted == xs2)
# print(predicted == xs1)

過去の内部状態の復元

過去の内部状態を復元するためには、内部状態の更新の逆操作をする必要があります。実装は上記のサイトから引用しています。

def get_prev_state(state):
    for i in range(623, -1, -1):
        result = 0
        tmp = state[i]
        tmp ^= state[(i + 397) % 624]
        if ((tmp & 0x80000000) == 0x80000000):
            tmp ^= 0x9908b0df
        result = (tmp << 1) & 0x80000000
        tmp = state[(i - 1 + 624) % 624]
        tmp ^= state[(i + 396) % 624]
        if ((tmp & 0x80000000) == 0x80000000):
            tmp ^= 0x9908b0df
            result |= 1
        result |= (tmp << 1) & 0x7fffffff
        state[i] = result
    return state

まずはさきほどの方法で現在の内部状態を復元します。その後、get_prev_state で内部状態を遡り、検証用の乱数列（過去） xs0 と比較します。今回は mti を 0 にします。

mt_state = [untemper(x) for x in xs1]
prev_mt_state = get_prev_state(mt_state)
random.setstate((3, tuple(prev_mt_state + [0]), None))

predicted = [random.getrandbits(32) for _ in range(N)]
print(predicted == xs0)

True