random memes }

Mersenne twister - random number generator

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