疑似乱数メルセンヌツイスタは連続 624 項が分かれば完全に再現できる

標準ライブラリのメルセンヌツイスタは実装と定数が公開されている。すると、出力された疑似乱数列から内部状態をコピーすることができる。

メルセンヌツイスタの実装

定数

標準ライブラリは周期が 2^19937 - 1 になるような定数があらかじめ設定されている。

関数

配列 unsigned mt[624] は長い周期で変化するカウンタの役割をする。mt[624] は実質リングバッファで、index でずれていく。
wikipedia の実装では「初めて乱数を生成するまで初期化しない」「twistは 624 回に 1 回だけまとめて行う」ような効率化がされていた。

乱数を1つ生成する関数の中で tempering （調律）をしている

unsigned y = mt[index];
y ^= (y >> U) & D;
y ^= (y << S) & B;
y ^= (y << T) & C;
y ^= (y >> L);

最終的にこの y を出力するので、tempering の逆関数を作ると配列 mt をコピーできるはず。

逆関数

定数に依存すれば結構簡単な逆関数ができる。

static unsigned inverse_tempering(unsigned y) {
    // inverse y ^= y >> L;
    y ^= y >> L;

    // inverse y ^= (y << T) & C; 
    y ^= (y << T) & C;

    // inverse y ^= (y << S) & B;
    y ^= (y << S) & 0x00001680;
    y ^= (y << S) & 0x000C4000;
    y ^= (y << S) & 0x0D200000;
    y ^= (y << S) & 0x90000000;

    // inverse y ^= (y >> U) & D;
    y ^= (y >> U) & 0x001FFC00;
    y ^= (y >> U) & 0x000003FF;

    return y;
}

コード全体

struct Random {
private:
    static const int W = 32; // word size;
    static const int N = 624; // degree of recurrence;
    static const int M = 397; // shift size;
    static const int R = 31; // mask bits;
    static const unsigned A = 0x9908B0DF; // XOR mask;
    
    // tempering;
    static const int U = 11;
    static const unsigned D = 0xFFFFFFFF;
    static const int S = 7;
    static const unsigned B = 0x9D2C5680;
    static const int T = 15;
    static const unsigned C = 0xEFC60000;
    static const int L = 18;

    static const unsigned F = 1812433253; // initialization;

    static const unsigned LOWER_MASK = (1u << R) - 1u;
    static const unsigned UPPER_MASK = ~0u ^ LOWER_MASK;

    static const unsigned DEFAULT_SEED = 5489;

    unsigned mt[N];
    int index;

    void twist(int i) {
	unsigned x = (mt[i] & UPPER_MASK) | (mt[(i + 1) % N] & LOWER_MASK);
	unsigned xA = x >> 1;
	if (x % 2 != 0) {
	    xA ^= A;
	}
	mt[i] = mt[(i + M) % N] ^ xA;
    }

    void initialize(unsigned seed) {
	mt[0] = seed;
	for (int i=1; i<N; i++) {
	    mt[i] = F * (mt[i-1] ^ (mt[i-1] >> (W - 2))) + i;
	}
	for (int i=0; i<N; i++) {
	    twist(i);
	}
	index = 0;
    }

    unsigned extract_number() {
	if (index >= N) {
	    index = 0;
	}

	unsigned y = mt[index];
	y ^= (y >> U) & D;
	y ^= (y << S) & B;
	y ^= (y << T) & C;
	y ^= (y >> L);

	twist(index);

	index++;
	return y;
    }

    static unsigned inverse_tempering(unsigned y) {
	// inverse y ^= y >> L;
	y ^= y >> L;

	// inverse y ^= (y << T) & C; 
	y ^= (y << T) & C;

	// inverse y ^= (y << S) & B;
	y ^= (y << S) & 0x00001680;
	y ^= (y << S) & 0x000C4000;
	y ^= (y << S) & 0x0D200000;
	y ^= (y << S) & 0x90000000;

	// inverse y ^= (y >> U) & D;
	y ^= (y >> U) & 0x001FFC00;
	y ^= (y >> U) & 0x000003FF;

	return y;
    }

public:
    Random(int seed=DEFAULT_SEED) {
	initialize(seed);
    }

    Random(const vector<unsigned> &prv) {
	assert(prv.size() == N);
	for (int i=0; i<N; i++) {
	    mt[i] = inverse_tempering(prv[i]);
	}
	for (int i=0; i<N; i++) {
	    twist(i);
	}
	index = 0;
    }

    unsigned operator()() {
	return extract_number();
    }
};

void MAIN() {
    mt19937 engine(100);

    vector<unsigned> prv;
    REP (i, 624) {
	prv.push_back(engine());
    }
    Random random(prv);

    assert(engine() == random());
}