Mersenne twister - Wikipedia, the free encyclopedia The Mersenne twister is a pseudorandom number generator developed in 1997 by Makoto Matsumoto (松本 眞) and Takuji Nishimura (西村 拓士)[1] that is based on a matrix linear recurrence over a finite binary field \mathbb{F}_2. It provides for fast generation of very high quality pseudorandom numbers, having been designed specifically to rectify many of the flaws found in older algorithms.

Once again, ran across a reference to the “Mersenne twister” as an especially good random number generator. The various implementations I found on the web had rather too many bells and whistles for my taste. Transliterated the algorithm in the Wikipedia article into C++, then applied a bit of optimization. The result generates 32-bit unsigned integers, and is simple, short … and fast (41.1 million random values per second on my AMD Athlon 2200+ laptop).

File: ZRandom.h

#ifndef __ZRANDOM_H__
#define __ZRANDOM_H__

//
//  Mersenne twister - random number generator.
//  Generate uniform distribution of 32 bit integers with the MT19937 algorithm.
//

class ZRandom
{

private:
    enum { N = 624, M = 397 };
    unsigned MT[N+1];
    unsigned* map[N];
    int nValues;

public:
    ZRandom(unsigned iSeed = 20070102);
    unsigned getValue();

};

#endif

File: ZRandom.cpp

#include "ZRandom.h"

ZRandom::ZRandom(unsigned iSeed) : nValues(0)
{
    // Seed the array used in random number generation.
    MT[0] = iSeed;
    for (int i=1; i<N ; ++i) {
        MT[i] = 1 + (69069 * MT[i-1]);
    }
    // Compute map once to avoid % in inner loop.
    for (i=0; i<N; ++i) {
        map[i] = MT + ((i + M) % N);
    }
}

unsigned ZRandom::getValue()
{
    if (0 == nValues) {
        MT[N] = MT[0];
        for (int i=0; i<N ; ++i) {
            register unsigned y = (0x80000000 & MT[i]) | (0x7FFFFFFF & MT[i+1]);
            register unsigned v = *(map[i]) ^ (y >> 1);
            if (1 & y) v ^= 2567483615;
            MT[i] = v;
        }
        nValues = N;
    }
    register unsigned y = MT[N - nValues--];
    y ^= y >> 11;
    y ^= (y << 7) & 2636928640;
    y ^= (y << 15) & 4022730752;
    y ^= y >> 18;
    return y;
}

… and a simple test program main.cpp

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include "ZRandom.h"

class ZTimer
{
    clock_t t1;
    clock_t t2;
public:
    ZTimer() : t1(::clock()), t2(t1) {}
    int elapsed() { return (int)(((t2 - t1) * 1000.0) / CLOCKS_PER_SEC); }
    int split() { t2 = ::clock(); return elapsed(); }
};

int main(int ac,char** av)
{
    enum { dtSample = 30000 };
    {
        ZRandom o;
        ZTimer t;
        int n = 0;
        do {
            for (int i=0; i<1000000; ++i) {
                unsigned v = o.getValue();
            }
            ++n;
        } while (t.split() < dtSample);
        int dt = t.elapsed();
        printf("Computed %d million primes in %d MS - %0.1f m/s \n",n,dt,((n * 1000.0)/dt));
    }

    return 0;
}