## Mersenne Twister

### September 9, 2011

We begin with some definitions:

`(define n 624)`

(define m 397)

```
```(define upper-mask #x80000000)

(define lower-mask #x7FFFFFFF)

`(define mag01 (vector #x0 #x9908B0DF))`

The current state of the generator is stored in two variables: the `mt`

vector stores the state, and the `mti`

variable points to the current position within the vector:

`(define mt (make-vector n 0))`

(define mti (+ n 1))

The `sgenrand`

function uses a simple linear-congruential generator to create the full state from a single seed:

`(define (sgenrand seed)`

(vector-set! mt 0 (logand seed #xFFFFFFFF))

(do ((mti 1 (+ mti 1))) ((= mti n))

(vector-set! mt mti

(logand (* 69069 (vector-ref mt (- mti 1))) #xFFFFFFFF)))

(set! mti n))

Function `genrand`

gets the next random number, resetting the current state of the generator as a side-effect. There are two clauses. The outer `when`

recalculates the current state every *n* calls to the generator. The `let*`

tempers the current state to produce a random number:

`(define (genrand)`

(when (<= n mti)

(when (= mti (+ n 1)) (sgenrand 4357))

(do ((kk 0 (+ kk 1))) ((= kk (- n m)))

(let ((y (logior

(logand (vector-ref mt kk) upper-mask)

(logand (vector-ref mt (+ kk 1)) lower-mask))))

(vector-set! mt kk

(logxor

(logxor (vector-ref mt (+ kk m)) (ash y -1))

(vector-ref mag01 (logand y #x1))))))

(do ((kk (- n m) (+ kk 1))) ((= kk (- n 1)))

(let ((y (logior

(logand (vector-ref mt kk) upper-mask)

(logand (vector-ref mt (+ kk 1)) lower-mask))))

(vector-set! mt kk

(logxor

(logxor (vector-ref mt (+ kk m (- n))) (ash y -1))

(vector-ref mag01 (logand y #x1))))))

(let ((y (logior

(logand (vector-ref mt (- n 1)) upper-mask)

(logand (vector-ref mt 0) lower-mask))))

(vector-set! mt (- n 1)

(logxor

(logxor (vector-ref mt (- m 1)) (ash y -1))

(vector-ref mag01 (logand y #x1)))))

(set! mti 0))

(let* ((y (vector-ref mt mti))

(y (logxor y (ash y -11)))

(y (logxor y (logand (ash y 7) #x9D2C5680)))

(y (logxor y (logand (ash y 15) #xEFC60000)))

(y (logxor y (ash y -18))))

(set! mti (+ mti 1))

y))

The `test-rand`

function produces the same output at the original C program:

`(define (test-rand)`

(sgenrand 4357)

(do ((j 0 (+ j 1))) ((= j 1000))

(printf "~10d " (genrand))

(if (= (modulo j 8) 7) (newline) (display " "))))

You can run the program at http://programmingpraxis.codepad.org/xDgJGH16, where you also find the bit operators `logand`

, `logior`

, `logxor`

and `ash`

from the Standard Prelude, which you will need if your Scheme system doesn’t provide them natively.

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.

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.

The compiled MersenneTwister can then be used from Java, too:

[…] 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 […]