// A straightforward C++ implementation of MT19937
// Haruhiko Okumura <okumura@matsusaka-u.ac.jp>
// Based on C source by Takuji Nishimura and Makoto Matsumoto
// (initialization improved 2002-01-26)

#include <iostream>
#include <iomanip>

using namespace std;

typedef unsigned long uint32;

class MersenneTwister {
    static const int N = 624;
    static const int M = 397;
    static const uint32 MAT_A = 0x9908b0dfUL;
    static const uint32 UMASK = 0x80000000UL;
    static const uint32 LMASK = 0x7fffffffUL;
    int mti;
    uint32 mt[N];               // array for the state vector

    inline uint32 twist(uint32 u, uint32 v) {
        return (((u & UMASK) | (v & LMASK)) >> 1) ^ ((v & 1) * MAT_A);
    }

    void init(uint32 s) {
        mti = N;
        mt[0] = s & 0xffffffffUL;
        for (int i = 1; i < N; i++) {
            mt[i] = 1812433253UL * (mt[i-1] ^ (mt[i-1] >> 30)) + i;
            mt[i] &= 0xffffffffUL;
        }
    }

public:
    MersenneTwister(uint32 s = 5489) {
        init(s);
    }

    MersenneTwister(uint32 init_key[], int key_length) {
        init(19650218UL);
        int i = 1;
        int j = 0;
        for (int k = (N > key_length) ? N : key_length; k != 0; k--) {
            mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) + init_key[j] + j;
            mt[i] &= 0xffffffffUL;
            if (++i >= N) { mt[0] = mt[N-1]; i = 1; }
            if (++j >= key_length) j = 0;
        }
        for (int k = N - 1; k != 0; k--) {
            mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) - i;
            mt[i] &= 0xffffffffUL;
            if (++i >= N) { mt[0] = mt[N-1]; i = 1; }
        }
        mt[0] = 0x80000000UL;
    }

    uint32 genrand_int32() {
        if (mti == N) {
            for (int k = 0; k < N - M; k++)
                mt[k] = mt[k+M] ^ twist(mt[k], mt[k+1]);
            for (int k = N - M; k < N - 1; k++)
                mt[k] = mt[k-(N-M)] ^ twist(mt[k], mt[k+1]);
            mt[N-1] = mt[M-1] ^ twist(mt[N-1], mt[0]);
            mti = 0;
        }
        uint32 y = mt[mti++];
        y ^= (y >> 11);
        y ^= (y <<  7) & 0x9d2c5680UL;
        y ^= (y << 15) & 0xefc60000UL;
        y ^= (y >> 18);
        return y;
    }

    // generates a random number on [0,0x7fffffff]-interval
    long genrand_int31() {
        return long(genrand_int32() >> 1);
    }

    // What follows are by Isaku Wada
    // generates a random number on [0,1]-real-interval
    double genrand_real1() {
        return genrand_int32() * (1.0 / 4294967295.0);
    }

    // generates a random number on [0,1)-real-interval
    double genrand_real2() {
        return genrand_int32() * (1.0/4294967296.0);
    }

    // generates a random number on (0,1)-real-interval
    double genrand_real3() {
        return (((double)genrand_int32()) + 0.5) * (1.0/4294967296.0);
    }

    // generates a random number on [0,1) with 53-bit resolution
    double genrand_res53() { 
        uint32 a = genrand_int32() >> 5;
        uint32 b = genrand_int32() >> 6;
        return(a * 67108864.0 + b) * (1.0 / 9007199254740992.0);
    }
};

#if 1
int main()
{
    uint32 init[4] = {0x123, 0x234, 0x345, 0x456};
    MersenneTwister MT(init, 4);
    cout << "1000 outputs of genrand_int32()\n";
    for (int i = 0; i < 1000; i++) {
        cout << setw(10) << MT.genrand_int32() << ' ';
        if (i % 5 == 4) cout << '\n';
    }
    cout << "\n1000 outputs of genrand_real2()\n";
    cout << fixed << setprecision(8);
    for (int i = 0; i < 1000; i++) {
        cout << setw(10) << MT.genrand_real2() << ' ';
        if (i % 5 == 4) cout << '\n';
    }
}
#else
int main()
{
    uint32 init[4] = {0x123, 0x234, 0x345, 0x456};
    MersenneTwister MT(init, 4);
    cout << "Start\n";
    uint32 x;
    for (int i = 0; i < 100000000; i++)
        x = MT.genrand_int32();
    cout << x << "\n";
}
#endif