Programming Praxis

4 Responses to “Mersenne Twister”

1. Mike said

   September 9, 2011 at 11:20 PM

   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

   September 13, 2011 at 12:29 PM

   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

   September 13, 2011 at 7:02 PM

   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();
     }
   }
   
4. Random Number Test | Programming Praxis said

   July 7, 2017 at 9:01 AM

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