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.

About these ads

Pages: 1 2 3

3 Responses to “Mersenne Twister”

  1. Mike said

    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))
    
  2. Graham said

    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.

  3. Jamie Hope said

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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 615 other followers

%d bloggers like this: