Mersenne Twister
September 9, 2011
The original article included versions for both integer and floating-point random numbers; we give here the integer version, on which the floating-point version is based:
/* A C-program for MT19937: Integer version (1998/4/6) */
/* genrand() generates one pseudorandom unsigned integer (32bit) */
/* which is uniformly distributed among 0 to 2^32-1 for each */
/* call. sgenrand(seed) set initial values to the working area */
/* of 624 words. Before genrand(), sgenrand(seed) must be */
/* called once. (seed is any 32-bit integer except for 0). */
/* Coded by Takuji Nishimura, considering the suggestions by */
/* Topher Cooper and Marc Rieffel in July-Aug. 1997. Comments */
/* should be addressed to: matumoto@math.keio.ac.jp */
#include<stdio.h>
/* Period parameters */
#define N 624
#define M 397
#define MATRIX_A 0x9908b0df /* constant vector a */
#define UPPER_MASK 0x80000000 /* most significant w-r bits */
#define LOWER_MASK 0x7fffffff /* least significant r bits */
/* Tempering parameters */
#define TEMPERING_MASK_B 0x9d2c5680
#define TEMPERING_MASK_C 0xefc60000
#define TEMPERING_SHIFT_U(y) (y >> 11)
#define TEMPERING_SHIFT_S(y) (y << 7)
#define TEMPERING_SHIFT_T(y) (y << 15)
#define TEMPERING_SHIFT_L(y) (y >> 18)
static unsigned long mt[N]; /* the array for the state vector */
static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */
/* initializing the array with a NONZERO seed */
void
sgenrand(seed)
unsigned long seed;
{
/* setting initial seeds to mt[N] using */
/* the generator Line 25 of Table 1 in */
/* [KNUTH 1981, The Art of Computer Programming */
/* Vol. 2 (2nd Ed.), pp102] */
mt[0]= seed & 0xffffffff;
for (mti=1; mti<N; mti++)
mt[mti] = (69069 * mt[mti-1]) & 0xffffffff;
}
unsigned long
genrand()
{
unsigned long y;
static unsigned long mag01[2]={0x0, MATRIX_A};
/* mag01[x] = x * MATRIX_A for x=0,1 */
if (mti >= N) { /* generate N words at one time */
int kk;
if (mti == N+1) /* if sgenrand() has not been called, */
sgenrand(4357); /* a default initial seed is used */
for (kk=0;kk<N-M;kk++) {
y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[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) ^ mag01[y & 0x1];
}
y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];
mti = 0;
}
y = mt[mti++];
y ^= TEMPERING_SHIFT_U(y);
y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
y ^= TEMPERING_SHIFT_L(y);
return y;
}
/* this main() outputs first 1000 generated numbers */
main()
{
int j;
sgenrand(4357); /* any nonzero integer can be used as a seed */
for (j=0; j<1000; j++) {
printf("%10lu ", genrand());
if (j%8==7) printf("\n");
}
printf("\n");
}
You can run this program at http://codepad.org/cwWsQUbZ.
Programming Praxis 
Python 3 version, written as a generator. There are no global variables, so multiple generators run independently.
First, the state vector is intialized from the seed and then “stirred”. On each call, a random number is generated from the next number in the state vector (mt[mti]). When all numbers in the state vector have been used (mti == N), the state vector is stirred again.
I reorganized the c-code so I could understand what was happening. The state vector is initialized from the seed. For each random number, the index in to the state vector is advanced and only the indexed number in the state vector is stirred. A random number is then derived from the indexed number.
N = 624 M = 397 MATRIX_A = 0x9908b0df UPPER_MASK = 0x80000000 LOWER_MASK = 0x7fffffff TEMPERING_MASK_B = 0x9d2c5680 TEMPERING_MASK_C = 0xefc60000 TEMPERING_SHIFT_U = lambda y: (y >> 11) TEMPERING_SHIFT_S = lambda y: (y << 7) TEMPERING_SHIFT_T = lambda y: (y << 15) TEMPERING_SHIFT_L = lambda y: (y >> 18) MASK32 = 0xffffffff def MTwister(seed=4357): if not seed: raise ValueError("Seed must be non-zero") # initialize the state vector fro mthe seed sv = [seed & MASK32] for n in range(1,N): sv.append(69069 * sv[-1] & MASK32) ndx = -1 while True: # advance the index into the state vector ndx = (ndx + 1)%N # stir the indexed item y = (sv[ndx]&UPPER_MASK)|(sv[(ndx+1)%N]&LOWER_MASK) sv[ndx] = sv[(ndx+M)%N] ^ (y >> 1) if y & 0x1: sv[ndx] ^= MATRIX_A # derive random number from the indexed item rn = sv[ndx] rn ^= TEMPERING_SHIFT_U(rn) rn ^= TEMPERING_SHIFT_S(rn) & TEMPERING_MASK_B rn ^= TEMPERING_SHIFT_T(rn) & TEMPERING_MASK_C rn ^= TEMPERING_SHIFT_L(rn) yield rn # test from itertools import islice res = list(islice(MTwister(),1000))I’ve wanted to try my hand at this exercise, since I’m interested in PRNGs; however, I don’t think I’ll be able to come up with more elegant solutions than the two here! Nice work.
Here’s a Kawa Scheme version, based upon a combination of the posted Scheme and Python versions.
(define-simple-class MersenneTwister () (N ::int allocation: 'static 624) (M ::int allocation: 'static 397) (mag01 ::long[] allocation: 'static [#x0 #x9908b0df]) (mt ::long[] (long[] length: MersenneTwister:N)) (mti ::int) ((*init*) (sgenrand 4357)) ((*init* (seed ::long)) (sgenrand seed)) ((sgenrand (seed ::long)) ::void (set! (mt 0) (bitwise-and seed #xffffffff)) (do ((i 1 (+ i 1))) ((= i N)) (set! (mt i) (bitwise-and (* 69069 (mt (- i 1))) #xffffffff))) (set! mti -1)) ((genrand) ::long (set! mti (remainder (+ mti 1) N)) (let ((y (bitwise-ior (bitwise-and (mt mti) #x80000000) (bitwise-and (mt (remainder (+ mti 1) N)) #x7fffffff)))) (set! (mt mti) (bitwise-xor (mt (remainder (+ mti M) N)) (ash y -1) (mag01 (bitwise-and y #x1))))) (let* ((y (mt mti)) (y (bitwise-xor y (ash y -11))) (y (bitwise-xor y (bitwise-and (ash y 7) #x9d2c5680))) (y (bitwise-xor y (bitwise-and (ash y 15) #xefc60000)))) (bitwise-xor y (ash y -18))))) (define (test-rand) (let ((mt (MersenneTwister))) (do ((j 0 (+ j 1))) ((= j 1000)) (format #t "~10@A " (mt:genrand)) (if (= (remainder j 8) 7) (newline)))) (newline)) (test-rand)The compiled MersenneTwister can then be used from Java, too:
public class TestMersenne { public static void main(String[] args) { MersenneTwister mt = new MersenneTwister(); mt.sgenrand(4357); for (int j = 0; j < 1000; ++j) { System.out.printf("%10d ", mt.genrand()); if (j % 8 == 7) System.out.println(); } System.out.println(); } }[…] built several random number generators: [1], [2], [3], [4], [5], [6], [7], [8], [9] (I didn’t realize it was so many until I went back and looked). In […]