Mersenne Twister - a verfy fast random number generator

2019/03/09 |

Programming/ Java

| Programming Java Randomization

Mersenne Twister

Mersenne Twister (mt19937) is a verfy fast random number generator of period 2^19937-1, and basically for Monte-Carlo simulations - it is not cryptographically secure “as is”.

What is Mersenne Twister (MT)

Mersenne Twister(MT) is a pseudorandom number generating algorithm developped by Makoto Matsumoto and Takuji Nishimura (alphabetical order) in 1996/1997. An improvement on initialization was given on 2002 Jan.

MT has the following merits:

• It is designed with consideration on the flaws of various existing generators.
• Far longer period and far higher order of equidistribution than any other implemented generators. (It is proved that the period is 2^19937-1, and 623-dimensional equidistribution property is assured.)
• Fast generation. (Although it depends on the system, it is reported that there are no much difference in speed between MT and the standard ANSI-C library function rand().

Boost C++ Libraries

#include <boost/random/mersenne_twister.hpp>
#include <boost/random/discrete_distribution.hpp>
#include <boost/format.hpp>

int main()
{
    boost::mt19937 rng(5489);
    for(int i = 0; i < 10; ++i) {
        std::cout << (boost::format("0x%08X") % rng()) << std::endl;
    }
    std::cout << std::endl;
}

/*
0xD091BB5C
0x22AE9EF6
0xE7E1FAEE
0xD5C31F79
0x2082352C
0xF807B7DF
0xE9D30005
0x3895AFE1
0xA1E24BBA
0x4EE4092B
 */

Java implementation

public class MT19937 {
    private final static int UPPER_MASK = 0x80000000;
    private final static int LOWER_MASK = 0x7fffffff;
    private final static int N = 624;
    private final static int M = 397;
    private final static int[] MAGIC = {0x0, 0x9908b0df};
    private final static int MAGIC_FACTOR1 = 1812433253;
    private final static int MAGIC_FACTOR2 = 1664525;
    private final static int MAGIC_FACTOR3 = 1566083941;
    private final static int MAGIC_MASK1 = 0x9d2c5680;
    private final static int MAGIC_MASK2 = 0xefc60000;
    private final static int MAGIC_SEED = 0x12BD6AA;

    // Internal state
    private final transient int[] mt = new int[N];
    private transient int mti;

    public MT19937() {
        this(Double.hashCode(System.currentTimeMillis() * (double) System.nanoTime()));
    }

    public MT19937(int seed) {
        setSeed(seed);
    }

    public MT19937(int[] seed) {
        setSeed(seed);
    }

    public MT19937(byte[] seed) {
        setSeed(pack(seed));
    }

    /*
        byte[] { 0x01, 0x02, 0x03, 0x04, 0x05, 0x06 }
        ->
        int[]  { 0x04030201, 0x00000605 }
     */
    private final static int[] pack(byte[] buf) {
        int up = ((buf.length + 3) >>> 2);
        int[] ibuf = new int[up];
        for (int n = 0; n < up; n++) {
            int k, m = (n + 1) << 2;
            if (m > buf.length)
                m = buf.length;
            for (k = buf[--m] & 0xff; (m & 0x3) != 0; k = (k << 8) | buf[--m] & 0xff)
                ;
            ibuf[n] = k;
        }
        return ibuf;
    }

    // 1 <= buf.length << 624
    public final void setSeed(int[] buf) {
        setSeed(MAGIC_SEED);

        int length = buf.length;
        if (length == 0)
            return;

        int i = 1, j = 0, k = (N > length ? N : length);
        for (; k > 0; k--) {
            mt[i] = (mt[i] ^ ((mt[i - 1] ^ (mt[i - 1] >>> 30)) * MAGIC_FACTOR2)) + buf[j] + j;
            i++;
            j++;
            if (i >= N) {
                mt[0] = mt[N - 1];
                i = 1;
            }
            if (j >= length)
                j = 0;
        }
        for (k = N - 1; k > 0; k--) {
            mt[i] = (mt[i] ^ ((mt[i - 1] ^ (mt[i - 1] >>> 30)) * MAGIC_FACTOR3)) - i;
            i++;
            if (i >= N) {
                mt[0] = mt[N - 1];
                i = 1;
            }
        }
        mt[0] |= UPPER_MASK; // MSB is 1; assuring non-zero initial array
    }

    public final void setSeed(int seed) {
        mt[0] = seed;
        for (mti = 1; mti < N; mti++) {
            mt[mti] = (MAGIC_FACTOR1 * (mt[mti - 1] ^ (mt[mti - 1] >>> 30)) + mti);
        }
    }

    public final int next() {
        int y;
        if (mti >= N) {
            int kk;
            for (kk = 0; kk < N - M; kk++) {
                y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
                mt[kk] = mt[kk + M] ^ (y >>> 1) ^ MAGIC[y & 0x1];
            }
            for (; kk < N - 1; kk++) {
                y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
                mt[kk] = mt[kk + (M - N)] ^ (y >>> 1) ^ MAGIC[y & 0x1];
            }
            y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
            mt[N - 1] = mt[M - 1] ^ (y >>> 1) ^ MAGIC[y & 0x1];

            mti = 0;
        }

        y = mt[mti++];

        // Tempering
        y ^= (y >>> 11);
        y ^= (y << 7) & MAGIC_MASK1;
        y ^= (y << 15) & MAGIC_MASK2;
        y ^= (y >>> 18);

        return y;
    }

    public final void nextBytes(byte[] bytes) {
        int i, value, n = bytes.length & 0x7FFF_FFFC;
        for (i = 0; i < n; i += 4) {
            value = next();
            bytes[i] = (byte) (value >> 24);
            bytes[i + 1] = (byte) (value >> 16);
            bytes[i + 2] = (byte) (value >> 8);
            bytes[i + 3] = (byte) (value);
        }

        n = bytes.length & 0x03;
        if (n > 0) {
            value = next();
            if (n >= 3) {
                bytes[i++] = (byte) (value >> 16);
            }
            if (n >= 2) {
                bytes[i++] = (byte) (value >> 8);
            }
            if (n >= 1) {
                bytes[i++] = (byte) (value >> 8);
            }
        }
    }
}