# Mersenne twister

The Mersenne twister is a pseudo random 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 $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.

For a k-bit word length, the Mersenne Twister generates numbers with an almost uniform distribution in the range $[0, 2^k-1]$.

## Applications

The Mersenne Twister has been optimized for use with Monte Carlo simulations in a number of fields, including simulating complex biochemical pathways,[2] photon migration,[3] genome coalescence,[4] cellular biology,[5] and computational finance.[6]

The Mersenne twister is the default random number generator for Python,[7][8] Ruby,[9] R,[10] PHP,[11] MATLAB and also available in C++[12] since C++11.

It is also used as an optional Random Number Generator in SPSS set through "Transform" menu. Add-on implementations are provided by the Boost C++ Libraries,[13] NAG Numerical Library,[14] and Glib.[15]

The commonly used variant 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 219937 − 1. While a long period is not a guarantee of quality in a random number generator, short periods (such as the 232 common in many software packages) can be problematic.[16]
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. It passes most, but not all, of the even more stringent TestU01 Crush randomness tests.[17]

The algorithm in its native form is not suitable for cryptography (unlike Blum Blum Shub). 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 Mersenne twister has been proposed by Makoto Matsumoto et al. The authors claim speeds 1.5 to 2 times faster than Advanced Encryption Standard in counter mode.[18]

Another issue is that it can take a long time to turn a non-random initial state (notably the presence of many zeros) into output that passes randomness tests. A small lagged Fibonacci generator or linear congruential generator gets started much more quickly and usually is used to seed the Mersenne Twister with random initial values.[citation needed]

## k-distribution

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

Let truncv(x) denote the number formed by the leading v bits of x, and consider P of the kv-bit vectors
$(\text{trunc}_v(x_i), \, \text{trunc}_v(x_{i+1}), \, ..., \, \text{trunc}_v(x_{i+k-1})) \quad (0\leq i< P)$.
Then each of the 2kv 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 Mersenne Twister is sensitive to poor initialization and can take a long time to recover from a zero-excess initial state. An alternative, WELL ("Well Equidistributed Long-period Linear"), has quicker recovery, the same or better performance and equal randomness.[19]

## Algorithmic detail

The Mersenne Twister algorithm is a twisted generalised feedback shift register[20] (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 ≤ mn
• r: separation point of one word, or the number of bits of the lower bitmask, 0 ≤ rw - 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 2nw − 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 expressed as the recurrence relation

$x_{k+n} := x_{k+m} \oplus ({x_k}^u \mid {x_{k+1}}^l) A \qquad \qquad k=0,1,\ldots$

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

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

with In − 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}$

where

$\boldsymbol{x} := ({x_k}^u \mid {x_{k+1}}^l) \qquad \qquad k=0,1,\ldots$

In order to achieve the 2nw − 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 2nw − 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 = RA = T−1RT, T invertible. The tempering is defined in the case of Mersenne Twister as

y := x ⊕ (x >> u)
y := :y ⊕ ((y << s) & b)
y := :y ⊕ ((y << t) & c)
z := y ⊕ (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 = 9908B0DF16
• u = 11
• (s, b) = (7, 9D2C568016)
• (t, c) = (15, EFC6000016)
• l = 18

## Pseudocode

The following piece of pseudocode generates uniformly distributed 32-bit integers in the range [0, 232 − 1] with the MT19937 algorithm:

 // Create a length 624 array to store the state of the generator
int[0..623] MT
int index = 0

// Initialize the generator from a seed
function initialize_generator(int seed) {
i := 0
MT[0] := seed
for i from 1 to 623 { // loop over each other element
MT[i] := last 32 bits of(1812433253 * (MT[i-1] xor (right shift by 30 bits(MT[i-1]))) + i) // 0x6c078965
}
}

// Extract a tempered pseudorandom number based on the index-th value,
// calling generate_numbers() every 624 numbers
function extract_number() {
if index == 0 {
generate_numbers()
}

int y := MT[index]
y := y xor (right shift by 11 bits(y))
y := y xor (left shift by 7 bits(y) and (2636928640)) // 0x9d2c5680
y := y xor (left shift by 15 bits(y) and (4022730752)) // 0xefc60000
y := y xor (right shift by 18 bits(y))

index := (index + 1) mod 624
return y
}

// Generate an array of 624 untempered numbers
function generate_numbers() {
for i from 0 to 623 {
int y := (MT[i] & 0x80000000)                       // bit 31 (32nd bit) of MT[i]
+ (MT[(i+1) mod 624] & 0x7fffffff)   // bits 0-30 (first 31 bits) of MT[...]
MT[i] := MT[(i + 397) mod 624] xor (right shift by 1 bit(y))
if (y mod 2) != 0 { // y is odd
MT[i] := MT[i] xor (2567483615) // 0x9908b0df
}
}
}


## SFMT

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

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

## MTGP

MTGP is a variant of Mersenne Twister optimised for GPUs published by Mutsuo Saito and Makoto Matsumoto.[24] 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. Sample code [1] for CUDA includes parameter sets suitable for 256, 512 and 1024 parallel threads per block, and up to 200 blocks generating independent random streams. The paper claims improved equidistribution over MT and performance on a high specification GPU (Nvidia GTX260 with 192 cores) of 4.7ms for 5x107 random 32-bit integers.

## References

1. ^ Matsumoto, M.; Nishimura, T. (1998). "Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator". ACM Transactions on Modeling and Computer Simulation 8 (1): 3–30. doi:10.1145/272991.272995. edit
2. ^ Hoops, Stefan; Sven Sahle, Ralph Gauges, Christine Lee, Jürgen Pahle, Natalia Simus, Mudita Singhal, Liang Xu, Pedro Mendes,Ursula Kummer (10 October 2006). "COPASI—a COmplex PAthway SImulator". Bioinformatics 22 (24): 3067–3074.
3. ^ Fang, Qianqian; David A. Boas (26 October 2009). "Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units". Optics Express 17 (22): 20178–20190. doi:10.1364/OE.17.020178.
4. ^ Liang, Liming; Sebastian Zöllner, Gonçalo R. Abecasis (25 April 2007). "GENOME: a rapid coalescent-based whole genome simulator". Bioinformatics 23 (12): 1565–1567.
5. ^ Takahashi, Kouichi; Kazunari Kaizu, Bin Hu, Masaru Tomita (22 January 2004). "A multi-algorithm, multi-timescale method for cell simulation". Bioinformatics 20 (4): 538–546.
6. ^ Singla, N; M. Hall, B. Shands, R. D. Chamberlain (16 November 2008). "Financial Monte Carlo simulation on architecturally diverse systems". Workshop on High Performance Computational Finance (WHPCF): 1–7.
7. ^ "9.6 random — Generate pseudo-random numbers". Python v2.6.8 documentation. Retrieved 2012-05-29.
8. ^ "8.6 random — Generate pseudo-random numbers". Python v3.2 documentation. Retrieved 2012-05-29.
9. ^ ""Random" class documentation". Ruby 1.9.3 documentation. Retrieved 2012-05-29.
10. ^ "Random Number Generators". CRAN Task View: Probability Distributions. Retrieved 2012-05-29.
11. ^ "mt_srand". php documentation. Retrieved 2012-05-29.
12. ^ "std::mersenne_twister_engine". Pseudo Random Number Generation. Retrieved 2012-09-25.
13. ^ "boost/random/mersenne_twister.hpp". Boost C++ Libraries. Retrieved 2012-05-29.
14. ^ "G05 – Random Number Generators". NAG Library Chapter Introduction. Retrieved 2012-05-29.
15. ^ "Changes to GLib". GLib Reference Manual. Retrieved 2012-05-29.
16. ^ Note: 219937 is approximately 4.3 × 106001; this is many orders of magnitude larger than the estimated number of particles in the observable universe, which is 1087.
17. ^ P. L'Ecuyer and R. Simard, TestU01: "A C Library for Empirical Testing of Random Number Generators", ACM Transactions on Mathematical Software, 33, 4, Article 22, August 2007.
18. ^ Matsumoto, Makoto; Nishimura, Takuji; Hagita, Mariko; Saito, Mutsuo (2005). "Cryptographic Mersenne Twister and Fubuki Stream/Block Cipher"
19. ^ P. L'Ecuyer, Uniform Random Number Generators, in International Encyclopedia of Statistical Science, Lovric, Miodrag (Ed.), Springer-Verlag, 2010.
20. ^ Matsumoto, M.; Kurita, Y. (1992). "Twisted GFSR generators". ACM Transactions on Modeling and Computer Simulation 2 (3): 179–194. doi:10.1145/146382.146383. edit
21. ^ SIMD-oriented Fast Mersenne Twister (SFMT)
22. ^ SFMT:Comparison of speed