Mersenne twister

From Wikipedia, the free encyclopedia

The Mersenne twister is a pseudorandom number generator developed in 1997 by Makoto Matsumoto (松本眞) and Takuji Nishimura (西村拓士)^[1] that is based on a matrix linear recurrence over a finite binary field $\mathbb{F}_2$ . It provides for fast generation of very high quality pseudorandom numbers, having been designed specifically to rectify many of the flaws found in older algorithms.

Its name derives from the fact that period length is chosen to be a Mersenne prime. There are at least two common variants of the algorithm, differing only in the size of the Mersenne primes used. The newer and more commonly used one is the Mersenne Twister MT19937, with 32-bit word length. There is also a variant with 64-bit word length, MT19937-64, which generates a different sequence.

1 Advantages
2 Algorithmic detail
3 Application
4 Pseudocode
5 References
6 External links

[edit] Advantages

As the commonly used variant of Mersenne Twister, MT19937 has e.g. the following desirable properties:

It was designed to have a colossal period of 2¹⁹⁹³⁷ − 1 (the creators of the algorithm proved this property). In practice, there is little reason to use larger ones, as most applications do not require 2¹⁹⁹³⁷ unique combinations.
It has a very high order of dimensional equidistribution (see linear congruential generator). Note that this means, by default, that there is negligible serial correlation between successive values in the output sequence.
It is faster than all but the most statistically unsound generators.
It passes numerous tests for statistical randomness, including the stringent Diehard tests.

[edit] Algorithmic detail

The Mersenne Twister algorithm is a twisted generalised feedback shift register^[2] (twisted GFSR, or TGFSR) of rational normal form (TGFSR(R)), with state bit reflection and tempering. It is characterized by the following quantities:

w: word size (in number of bits)
n: degree of recurrence
m: middle word, or the number of parallel sequences, $1 \le m \le n$
r: separation point of one word, or the number of bits of the lower bitmask, $0 \le r \le 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, l: additional Mersenne Twister tempering bit shifts

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.

For a word x with w bit width, it is expresses as the recurrence relation

x_k + n := x_k + m $\oplus$ (x_k^u | x_k + 1^l) A, k = 0, 1, ...

with $|$ as the bitwise or and $\oplus$ as the bitwise exclusive or (XOR), x^u, x^l being x with upper and lower bitmasks applied. The twist transformation A is defined in rational normal form

$A = R = \begin{pmatrix} 0 & I_{n - 1} \\ a_n & (a_{n - 1}, \ldots a_0) \end{pmatrix}$

with I_n − 1 as the (n − 1) × (n − 1) identity matrix (and in contrast to normal matrix multiplication, bitwise XOR replaces addition). The rational normal form has the benefit that it can be efficiently expressed as

$\boldsymbol{x}A = \begin{cases}\boldsymbol{x} \gg 1 & x_0 = 0\\(\boldsymbol{x} \gg 1) \oplus \boldsymbol{a} & x_0 = 1\end{cases}$

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

$B = \begin{pmatrix} 0 & I_{w} & \cdots & 0 & 0 \\ \vdots & & & & \\ I_{w} & \vdots & \ddots & \vdots & \vdots \\ \vdots & & & & \\ 0 & 0 & \cdots & I_{w} & 0 \\ 0 & 0 & \cdots & 0 & I_{w - r} \\ S & 0 & \cdots & 0 & 0 \end{pmatrix} \begin{matrix} \\ \\ \leftarrow m\hbox{-th row} \\ \\ \\ \\ \end{matrix}$
$S = \begin{pmatrix} 0 & I_{r} \\ I_{w - r} & 0 \end{pmatrix} A$

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

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

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), which is equivalent to the transformation A = R → A = T⁻¹RT, T invertible. The tempering is defined in the case of Mersenne Twister as

y := x $\oplus$ (x >> u)
y := y $\oplus$ ((y << s) & b)
y := y $\oplus$ ((y << t) & c)
z := y $\oplus$ (y >> l)

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, $s + t \ge \lfloor w/2 \rfloor - 1$ 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 = 11
(s, b) = (7, 9D2C5680₁₆)
(t, c) = (15, EFC60000₁₆)
l = 18

[edit] Application

Unlike Blum Blum Shub, the algorithm in its native form is not suitable for cryptography. Observing a sufficient number of iterates (624 in the case of MT19937) allows one to predict all future iterates. Combining the Mersenne twister with a hash function solves this problem, but slows down generation.

For many other applications, however, the Mersenne twister is fast becoming the random number generator of choice.

[edit] Pseudocode

The following generates uniformly 32 bit integers in the range [0, 2³² − 1] with the MT19937 algorithm:

 // Create a length 624 array to store the state of the generator
 var int[0..623] MT
 var int y
 // Initialise the generator from a seed
 function initialiseGenerator ( 32-bit int seed ) {
     MT[0] := seed
     for i from 1 to 623 { // loop over each other element
         MT[i] := last_32bits_of((69069 * MT[i-1]) + 1)
     }
 }

 // Generate an array of 624 untempered numbers
 function generateNumbers() {
     for i from 0 to 622 {
         y := 32nd_bit_of(MT[i]) + last_31bits_of(MT[i+1])
         if y even {
             MT[i] := MT[(i + 397) % 624] bitwise xor (right_shift_by_1_bit(y))
         } else if y odd {
             MT[i] := MT[(i + 397) % 624] bitwise_xor (right_shift_by_1_bit(y)) bitwise_xor (2567483615)
         }
     }
     y := 32nd_bit_of(MT[623]) + last_31bits_of(MT[0])
     if y even {
         MT[623] := MT[396] bitwise_xor (right_shift_by_1_bit(y))
     } else if y odd {
         MT[623] := MT[396] bitwise_xor (right_shift_by_1_bit(y)) bitwise_xor (2567483615)
     }
 }
 
 // Extract a tempered pseudorandom number based on the i-th value
 function extractNumber(int i) {
     y := MT[i]
     y := y bitwise_xor (right_shift_by_11_bits(y))
     y := y bitwise_xor (left_shift_by_7_bits(y) bitwise_and (2636928640))
     y := y bitwise_xor (left_shift_by_15_bits(y) bitwise_and (4022730752))
     y := y bitwise_xor (right_shift_by_18_bits(y))
     return y
 }