Mersenne Twister
September 9, 2011
We have written random number generators in several previous exercises. Today we look at the Mersenne Twister of Makoto Matsumoto and Takuji Nishimura. Here is their description of their algorithm:
We propose a new random number generator Mersenne Twister. An implemented C-code MT19937 has the period 219937−1 and 523-dimensional equidistribution property, which seems to be the best among all generators ever implemented. There are two new ideas added to the previous twisted GFSR to attain these records. One is an incomplete array to realize a Mersenne-prime period. The other is a fast algorithm to test the primitivity of the characteristic polynomial of a linear recurrence, named inversive-decimation method. This algorithm does not require even the explicit form of the characteristic polynomial. It needs only (1) the defining recurrence, and (2) some fast algorithm that obtains the present state vector from its 1-bit output stream. The computational complexity of the inversive-decimation method is the order of the algorithm in (2) multiplied by the degree of the characteristic polynomial. To attain higher order equidistribution properties, we used the resolution-wise lattice method, with Lenstra’s algorithm for successive minima.
If that’s not quite clear, their reference implementation is given on the next page. Note that this implementation is the original Mersenne Twister MT19937; there are several variants.
Though George Marsaglia has strongly criticized it as needlessly complex, the Mersenne Twister is quite a popular random number generator. It takes a 32-bit integer other than 0 as a seed, returns a 32-bit integer each time it is called, and has a period of 219937−1 ≅ 4·106001. The Mersenne twister is suitable for simulation but not for cryptography.
Your task is to implement the Mersenne twister in your favorite language. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
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(); } }