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.

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))))

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

    public class TestMersenne {
      public static void main(String[] args) {
        MersenneTwister mt = new MersenneTwister();
        for (int j = 0; j < 1000; ++j) {
          System.out.printf("%10d ", mt.genrand());
          if (j % 8 == 7) System.out.println();

Leave a Reply

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

You are commenting using your 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


Get every new post delivered to your Inbox.

Join 709 other followers

%d bloggers like this: