Mersenne Twister

The Mersenne Twister is a pseudorandom number generator (PRNG). It is by far the most widely used general-purpose PRNG.^[1] Its name derives from the fact that its period length is chosen to be a Mersenne prime.

The Mersenne Twister was developed in 1997 by Makoto Matsumoto (松本 眞) and Takuji Nishimura (西村 拓士).^[2] It was designed specifically to rectify most of the flaws found in older PRNGs. It was the first PRNG to provide fast generation of high-quality pseudorandom integers.

The most commonly used version of the Mersenne Twister algorithm is based on the Mersenne prime 2¹⁹⁹³⁷−1. The standard implementation of that, MT19937, uses a 32-bit word length. There is another implementation that uses a 64-bit word length, MT19937-64; it generates a different sequence.

Adoption in software systems

The Mersenne Twister is the default PRNG for the following software systems:

Microsoft Visual C++,^[3] Microsoft Excel,^[4] GAUSS,^[5] GLib,^[6] GNU Multiple Precision Arithmetic Library,^[7] GNU Octave,^[8] GNU Scientific Library,^[9] gretl,^[10] IDL,^[11] Julia,^[12] CMU Common Lisp,^[13] Embeddable Common Lisp,^[14] Steel Bank Common Lisp,^[15] Maple,^[16] MATLAB,^[17] Free Pascal,^[18] PHP,^[19] Python,^[20][21] R,^[22] Ruby,^[23] SageMath,^[24] Scilab,^[25] Stata.^[26] It is also available in Apache Commons,^[27] in standard C++ (since C++11),^[28][29] and in Mathematica.^[30] Add-on implementations are provided in many program libraries, including the Boost C++ Libraries,^[31] the CUDA Library,^[32] and the NAG Numerical Library.^[33]

The Mersenne Twister is one of two PRNGs in SPSS: the other generator is kept only for compatibility with older programs, and the Mersenne Twister is stated to be "more reliable".^[34] The Mersenne Twister is similarly one of the PRNGs in SAS: the other generators are older and deprecated.^[35]

Advantages

The commonly used version of Mersenne Twister, MT19937, which produces a sequence of 32-bit integers, has the following desirable properties:

1. It has a very long period of 2¹⁹⁹³⁷ − 1. While a long period is not a guarantee of quality in a random number generator, short periods (such as the 2³² common in many older software packages) can be problematic.^[36]
2. It is k-distributed to 32-bit accuracy for every 1 ≤ k ≤ 623 (see definition below).
3. It passes numerous tests for statistical randomness, including the Diehard tests.

Disadvantages

The large state space comes with a performance cost: the 2.5 KiB state buffer will place a load on the memory caches. In 2011, Saito & Matsumoto proposed a version of the Mersenne Twister to address this issue. The tiny version, TinyMT, uses just 127 bits of state space.^[37]

By today's standards, the Mersenne Twister is somewhat slow, ^[38] unless the SFMT implementation is used (see section below).

It passes most, but not all, of the stringent TestU01 randomness tests.^[39] Because it is based on simple linear (xor) operations, it fails tests based on linear complexity after relatively few bits of output, despite its extremely large state. Passing the output through a simple hash function can remedy this weakness.

Multiple Mersenne Twister instances that differ only in seed value (but not other parameters) are not generally appropriate for Monte-Carlo simulations that require independent random number generators, though there exists a method for choosing multiple sets of parameter values.^[40][41]

It can take a long time to start generating output that passes randomness tests, if the initial state is highly non-random—particularly if the initial state has many zeros. A consequence of this is that two instances of the generator, started with initial states that are almost the same, will usually output nearly the same sequence for many iterations, before eventually diverging. The 2002 update to the MT algorithm has improved initialization, so that beginning with such a state is very unlikely.^[42]

k-distribution

A pseudorandom sequence x_i of w-bit integers of period P is said to be k-distributed to v-bit accuracy if the following holds.

Let trunc_v(x) denote the number formed by the leading v bits of x, and consider P of the kv-bit vectors

.

Then each of the 2^kv possible combinations of bits occurs the same number of times in a period, except for the all-zero combination that occurs once less often.

Alternatives

The algorithm in its native form is not cryptographically secure. The reason is that observing a sufficient number of iterations (624 in the case of MT19937, since this is the size of the state vector from which future iterations are produced) allows one to predict all future iterations.

A pair of cryptographic stream ciphers based on output from the Mersenne Twister has been proposed by Matsumoto, Nishimura, and co-authors. The authors claim speeds 1.5 to 2 times faster than Advanced Encryption Standard in counter mode.^[43]

An alternative generator, WELL ("Well Equidistributed Long-period Linear"), offers quicker recovery, and equal randomness, and nearly equal speed.^[44] Marsaglia's xorshift generators and variants are the fastest in this class.^[45]

Algorithmic detail

Visualisation of generation of pseudo-random 32-bit integers using a Mersenne Twister. The 'Extract number' section shows an example where integer 0 has already been output and the index is at integer 1. 'Generate numbers' is run when all integers have been output.

For a w-bit word length, the Mersenne Twister generates integers in the range [0, 2^w−1].

The Mersenne Twister algorithm is based on a matrix linear recurrence over a finite binary field F₂. The algorithm is a twisted generalised feedback shift register^[46] (twisted GFSR, or TGFSR) of rational normal form (TGFSR(R)), with state bit reflection and tempering. The basic idea is to define a series through a simple recurrence relation, and then output numbers of the form , where is an invertible F₂ matrix called a tempering matrix.

The general algorithm is characterized by the following quantities (some of these explanations make sense only after reading the rest of the algorithm):

• w: word size (in number of bits)
• n: degree of recurrence
• m: middle word, an offset used in the recurrence relation defining the series x, 1 ≤ m < n
• r: separation point of one word, or the number of bits of the lower bitmask, 0 ≤ r ≤ w - 1
• a: coefficients of the rational normal form twist matrix
• b, c: TGFSR(R) tempering bitmasks
• s, t: TGFSR(R) tempering bit shifts
• u, d, l: additional Mersenne Twister tempering bit shifts/masks

with the restriction that 2^nw − r − 1 is a Mersenne prime. This choice simplifies the primitivity test and k-distribution test that are needed in the parameter search.

The series x is defined as a series of w-bit quantities with the recurrence relation:

where denotes concatenation of bit vectors (with upper bits on the left), the bitwise exclusive or (XOR), means the upper bits of , and means the lower bits of . The twist transformation A is defined in rational normal form as:

with I_n − 1 as the (n − 1) × (n − 1) identity matrix. The rational normal form has the benefit that multiplication by A can be efficiently expressed as: (remember that here matrix multiplication is being done in F₂, and therefore bitwise XOR takes the place of addition)

where x₀ is the lowest order bit of x.

As like TGFSR(R), the Mersenne Twister is cascaded with a tempering transform to compensate for the reduced dimensionality of equidistribution (because of the choice of A being in the rational normal form). Note that this is equivalent to using the matrix A where A = T⁻¹AT for T an invertible matrix, and therefore the analysis of characteristic polynomial mentioned below still holds.

As with A, we choose a tempering transform to be easily computable, and so do not actually construct T itself. The tempering is defined in the case of Mersenne Twister as

y := x ⊕ ((x >> u) & d)

y := y ⊕ ((y << s) & b)

y := y ⊕ ((y << t) & c)

z := y ⊕ (y >> l)

where x is the next value from the series, y a temporary intermediate value, z the value returned from the algorithm, with <<, >> as the bitwise left and right shifts, and & as the bitwise and. The first and last transforms are added in order to improve lower-bit equidistribution. From the property of TGFSR, is required to reach the upper bound of equidistribution for the upper bits.

The coefficients for MT19937 are:

• (w, n, m, r) = (32, 624, 397, 31)
• a = 9908B0DF₁₆
• (u, d) = (11, FFFFFFFF₁₆)
• (s, b) = (7, 9D2C5680₁₆)
• (t, c) = (15, EFC60000₁₆)
• l = 18

Note that 32-bit implementations of the Mersenne Twister generally have d = FFFFFFFF₁₆. As a result, the d is occasionally omitted from the algorithm description, since the bitwise and with d in that case has no effect.

The coefficients for MT19937-64 are:^[47]

• (w, n, m, r) = (64, 312, 156, 31)
• a = B5026F5AA96619E9₁₆
• (u, d) = (29, 5555555555555555₁₆)
• (s, b) = (17, 71D67FFFEDA60000₁₆)
• (t, c) = (37, FFF7EEE000000000₁₆)
• l = 43

Initialization

As should be apparent from the above description, the state needed for a Mersenne Twister implementation is an array of n values of w bits each. To initialize the array, a w-bit seed value is used to supply x₀ through x_n − 1 by setting x₀ to the seed value and thereafter setting

x_i = f × (x_i-1 ⊕ (x_i-1 >> (w-2))) + i

for i from 1 to n-1. The first value the algorithm then generates is based on x_n, not based on x₀. The constant f forms another parameter to the generator, though not part of the algorithm proper. The value for f for MT19937 is 1812433253 and for MT19937-64 is 6364136223846793005.^[48]

Comparison with classical GFSR

In order to achieve the 2^nw − r − 1 theoretical upper limit of the period in a TGFSR, φ_B(t) must be a primitive polynomial, φ_B(t) being the characteristic polynomial of

The twist transformation improves the classical GFSR with the following key properties:

• The period reaches the theoretical upper limit 2^nw − r − 1 (except if initialized with 0)
• Equidistribution in n dimensions (e.g. linear congruential generators can at best manage reasonable distribution in five dimensions)

Pseudocode

The following piece of pseudocode implements the general Mersenne Twister algorithm. The constants w, n, m, r, a, u, d, s, b, t, c, l, and f are as in the algorithm description above. It is assumed that int represents a type sufficient to hold values with w bits:

 // Create a length n array to store the state of the generator
 int[0..n-1] MT
 int index := n+1
 const int lower_mask = (1 << r) - 1 // That is, the binary number of r 1's
 const int upper_mask = lowest w bits of (not lower_mask)
 
 // Initialize the generator from a seed
 function seed_mt(int seed) {
     index := n
     MT[0] := seed
     for i from 1 to (n - 1) { // loop over each element
         MT[i] := lowest w bits of (f * (MT[i-1] xor (MT[i-1] >> (w-2))) + i)
     }
 }
 
 // Extract a tempered value based on MT[index]
 // calling twist() every n numbers
 function extract_number() {
     if index >= n {
         if index > n {
           error "Generator was never seeded"
           // Alternatively, seed with constant value; 5489 is used in reference C code[49]
         }
         twist()
     }
 
     int y := MT[index]
     y := y xor ((y >> u) and d)
     y := y xor ((y << s) and b)
     y := y xor ((y << t) and c)
     y := y xor (y >> l)
 
     index := index + 1
     return lowest w bits of (y)
 }
 
 // Generate the next n values from the series x_i 
 function twist() {
     for i from 0 to (n-1) {
         int x := (MT[i] and upper_mask)
                   + (MT[(i+1) mod n] and lower_mask)
         int xA := x >> 1
         if (x mod 2) != 0 { // lowest bit of x is 1
             xA := xA xor a
         }
         MT[i] := MT[(i + m) mod n] xor xA
     }
     index := 0
 }

Python implementation

This python implementation hard-codes the constants for MT19937:

def _int32(x):
    # Get the 32 least significant bits.
    return int(0xFFFFFFFF & x)

class MT19937:

    def __init__(self, seed):
        # Initialize the index to 0
        self.index = 624
        self.mt = [0] * 624
        self.mt[0] = seed  # Initialize the initial state to the seed
        for i in range(1, 624):
            self.mt[i] = _int32(
                1812433253 * (self.mt[i - 1] ^ self.mt[i - 1] >> 30) + i)

    def extract_number(self):
        if self.index >= 624:
            self.twist()

        y = self.mt[self.index]

        # Right shift by 11 bits
        y = y ^ y >> 11
        # Shift y left by 7 and take the bitwise and of 2636928640
        y = y ^ y << 7 & 2636928640
        # Shift y left by 15 and take the bitwise and of y and 4022730752
        y = y ^ y << 15 & 4022730752
        # Right shift by 18 bits
        y = y ^ y >> 18

        self.index = self.index + 1

        return _int32(y)

    def twist(self):
        for i in range(624):
            # Get the most significant bit and add it to the less significant
            # bits of the next number
            y = _int32((self.mt[i] & 0x80000000) +
                       (self.mt[(i + 1) % 624] & 0x7fffffff))
            self.mt[i] = self.mt[(i + 397) % 624] ^ y >> 1

            if y % 2 != 0:
                self.mt[i] = self.mt[i] ^ 0x9908b0df
        self.index = 0
        #test

Then MT19937(seed).extract_number() returns the random number, where seed is the initial seed.

C# implementation

This C# implementation hard-codes the constants for MT19937-64:

using System;

public class MersenneTwister
{
    public const int w = 64;
    public const ulong n = 312;
    public const ulong m = 156;
    public const ulong r = 31;
    public const ulong a = 0xB5026F5AA96619E9;
    public const int u = 29;
    public const ulong d = 0x5555555555555555;
    public const int s = 17;
    public const ulong b = 0x71D67FFFEDA60000;
    public const int t = 37;
    public const ulong c = 0xFFF7EEE000000000;
    public const int l = 43;
    public const ulong f = 6364136223846793005;

    public const ulong lower_mask = 0x7FFFFFFF;
    public const ulong upper_mask = ~lower_mask;

    private ulong[] MT = new ulong[n];
    private ulong index = n + 1;

    public MersenneTwister(ulong seed)
    {
        seed_mt(seed);
    }

    private void seed_mt(ulong seed)
    {
        index = n;
        MT[0] = seed;

        for (ulong i = 1; i < n; ++i)
        {
            MT[i] = (f * (MT[i - 1] ^ (MT[i - 1] >> (w - 2))) + i);
        }
    }

    public ulong extract_number()
    {
        if (index >= n)
        {
            if (index > n)
            {
                throw new Exception("Generator was never seeded");
            }
            twist();
        }

        ulong y = MT[index];
        y = y ^ ((y >> u) & d);
        y = y ^ ((y << s) & b);
        y = y ^ ((y << t) & c);
        y = y ^ (y >> l);

        ++index;

        return y;
    }

    private void twist()
    {
        for (ulong i = 0; i < n; ++i)
        {
            ulong x = (MT[i] & upper_mask) + (MT[(i + 1) % n] & lower_mask);
            ulong xA = x >> 1;

            if (x % 2 != 0)
            {
                xA = xA ^ a;
            }

            MT[i] = MT[(i + m) % n] ^ xA;
        }

        index = 0;
    }
}

C/C++ implementation

Simple 32-bit C/C++ implementation (tested using GCC for ARM):

#include <stdint.h>

// Define MT19937 constants (32-bit RNG)
enum
{
    // Assumes W = 32 (omitting this)
    N = 624,
    M = 397,
    R = 31,
    A = 0x9908B0DF,

    F = 1812433253,

    U = 11,
    // Assumes D = 0xFFFFFFFF (omitting this)

    S = 7,
    B = 0x9D2C5680,

    T = 15,
    C = 0xEFC60000,

    L = 18,

    MASK_LOWER = (1ull << R) - 1,
    MASK_UPPER = (1ull << R)
};

static uint32_t  mt[N];
static uint16_t  index;

// Re-init with a given seed
void Initialize(const uint32_t  seed)
{
    uint32_t  i;

    mt[0] = seed;

    for ( i = 1; i < N; i++ )
    {
        mt[i] = (F * (mt[i - 1] ^ (mt[i - 1] >> 30)) + i);
    }

    index = N;
}

static void Twist()
{
    uint32_t  i, x, xA;

    for ( i = 0; i < N; i++ )
    {
        x = (mt[i] & MASK_UPPER) + (mt[(i + 1) % N] & MASK_LOWER);

        xA = x >> 1;

        if ( x & 0x1 )
            xA ^= A;

        mt[i] = mt[(i + M) % N] ^ xA;
    }

    index = 0;
}

// Obtain a 32-bit random number
uint32_t ExtractU32()
{
    uint32_t  y;
    int       i = index;

    if ( index >= N )
    {
        Twist();
        i = index;
    }

    y = mt[i];
    index = i + 1;

    y ^= (mt[i] >> U);
    y ^= (y << S) & B;
    y ^= (y << T) & C;
    y ^= (y >> L);

    return y;
}

Haskell implementation

This Haskell implementation is optimized for readability and uses an inefficient list representation of the state array :

{-# LANGUAGE ParallelListComp #-}

import Data.Bits
import Data.Word

(w, n, m, r) = (32, 624, 397, 31)
a            = 0x9908B0DF
(u, d)       = (11, 0xFFFFFFFF)
(s, b)       = (7, 0x9D2C5680)
(t, c)       = (15, 0xEFC60000)
l            = 18

-- | Initialize the MT state of 624 32-bit values from a 32-bit seed.
--
-- Usage example:
--
-- > st0 = mtInit 1
-- > z0  = mtNextValue st0
-- > st1 = mtNextState st0
-- > z1  = mtNextValue st1
-- > st2 = mtNextState st1
--
-- That initializes the state from the value 1, and then computes two
-- random values @z0@ and @z1@; the next random value would be
-- @mtNextValue st2@.
--
-- Note that the first 624 state values are not used to produce random
-- numbers, so the initial state we produce comprises the @x@ values
-- @[x_1,...,x_n]@.
mtInit :: Word32 -> [Word32]
mtInit seed = mtNextState st0
  where
    st0 = [seed] ++
          [ f * (x `xor` (x `shiftR` (w - 2))) + i
          | x <- st0
          | i <- [1 .. n - 1] ]
    f = 1812433253

-- | Compute the next MT state from the current state.
--
-- I.e., if the input state is @[x_i,...,x_{i+n-1}]@ then the output
-- state is @[x_{i+1},...,x_{i+n}]@.
mtNextState :: [Word32] -> [Word32]
mtNextState st = tail st ++ [nextStateValue]
  where
    nextStateValue = (st !! m) `xor`
                     multA (upper (st !! 0) .|. lower (st !! 1))
    multA x = if (x .&. 1) == 0
              then (x `shiftR` 1)
              else (x `shiftR` 1) `xor` a
    upper x = (x `shiftR` r)       `shiftL` r
    lower x = (x `shiftL` (w - r)) `shiftR` (w - r)

-- | Compute the next MT random value @z@ from the last value
-- @x_{i+n-1}@ of the current state @[x_i,...,x_{i+n-1}]@.
mtNextValue :: [Word32] -> Word32
mtNextValue st = z
  where
    x  = last st
    y0 = x  `xor` ((x  `shiftR` u) .&. d)
    y1 = y0 `xor` ((y0 `shiftL` s) .&. b)
    y2 = y1 `xor` ((y1 `shiftL` t) .&. c)
    z  = y2 `xor` (y2  `shiftR` l) 

-- | Compute an infinite stream of MT random values from a seed.
mtStream :: Word32 -> [Word32]
mtStream = map mtNextValue . iterate mtNextState . mtInit

Freepascal implementation

Output checked against G++ C++11 mt19937 and FPC's default random.
This is NOT the Mersenne twister implementation that Freepascal uses as its default prng, but the algorithm is the same.
Usage:

program testmer32;
uses mersenne;
begin
  writeln(mer32.random);
end.

unit mersenne;
{$mode objfpc}{$A8}{$R-}{$Q-}
interface
type
   TMersenne32 = class sealed
   strict private
   const 
     // Define MT19937 constants (32-bit RNG)
     N = 624;M = 397;R = 31;A = $9908B0DF;F = 1812433253;
     U = 11;S = 7;B = $9D2C5680;T = 15;C = $EFC60000;L = 18;
     MASK_LOWER = (QWORD(1) shl R) - 1;
     MASK_UPPER = QWORD(1) shl R;
     class var mt:array[0..N-1] of dword;
     class var index:word;
     class procedure twist;inline;static;
   public
     class constructor create;
     class procedure initialize(const seed:dword);inline;static;
     class function Random:dword;inline;static;
   end;
   
   mer32 = type TMersenne32;

implementation 

class constructor TMersenne32.Create;
begin
  initialize(5489);
end;

class procedure TMersenne32.Initialize(const seed:dword);inline;static;
var 
  i:dword;
begin
  mt[0] := seed;
 for  i := 1 to pred(N) do 
   mt[i] := F * (mt[i - 1] xor (mt[i - 1] shr 30)) + i;
 index := N;
end;

// works around the use of modulo, which is slow on some platforms   
class procedure TMersenne32.Twist;inline;static;
var 
  i:integer;
begin
  for i:=0 to N-M-1 do
    mt[i]:=mt[i+M] xor {twist} (((mt[i] and MASK_UPPER) or 
    (mt[i+1] and MASK_LOWER)) shr 1)xor(dword(-(mt[i+1] and 1)) and A);
  for i:=N-M to N-2 do
    mt[i]:=mt[i+(M-N)]xor{twist}(((mt[i] and MASK_UPPER) or 
    (mt[i+1] and MASK_LOWER)) shr 1)xor(dword(-(mt[i+1] and 1)) and A);
    mt[N-1]:=mt[M-1] xor {twist} (((mt[n-1] and MASK_UPPER) or (mt[0] and 
    MASK_LOWER)) shr 1)xor(dword(-(mt[0] and 1)) and A);
  index:=0;
end;

class function TMersenne32.Random:dword;inline;static;
var
  i:integer;
begin
  i := index;
  if  index >= N then
  begin
    Twist;
    i := index;
  end;
  Result := mt[i];
  index := i + 1;
  Result := Result xor (mt[i] shr U);
  Result := Result xor (Result shl S) and B;
  Result := Result xor (Result shl T) and C;
  Result := Result xor (Result shr L);
end;

end.

SFMT

The SFMT (SIMD-oriented Fast Mersenne Twister) is a variant of Mersenne Twister, introduced in 2006,^[50] designed to be fast when it runs on 128-bit SIMD.

• It is roughly twice as fast as Mersenne Twister.^[51]
• It has a better equidistribution property of v-bit accuracy than MT but worse than WELL ("Well Equidistributed Long-period Linear").
• It has quicker recovery from zero-excess initial state than MT, but slower than WELL.
• It supports various periods from 2⁶⁰⁷−1 to 2²¹⁶⁰⁹¹−1.

Intel SSE2 and PowerPC AltiVec are supported by SFMT. It is also used for games with the Cell BE in the PlayStation 3.^[52]

MTGP

MTGP is a variant of Mersenne Twister optimised for graphics processing units published by Mutsuo Saito and Makoto Matsumoto.^[53] The basic linear recurrence operations are extended from MT and parameters are chosen to allow many threads to compute the recursion in parallel, while sharing their state space to reduce memory load. The paper claims improved equidistribution over MT and performance on a high specification GPU (Nvidia GTX260 with 192 cores) of 4.7 ms for 5×10⁷ random 32-bit integers.

References

External links

• The academic paper for MT, and related articles by Makoto Matsumoto
• Mersenne Twister home page, with codes in C, Fortran, Java, Lisp and some other languages
• SFMT in Action: Part I – Generating a DLL Including SSE2 Support – on Code Project
• Harase, S. (2014), "On the F₂-linear relations of Mersenne Twister pseudorandom number generators", Mathematics and Computers in Simulation, 100: 103–113, doi:10.1016/j.matcom.2014.02.002