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.
The most commonly used version of the Mersenne Twister algorithm is based on the Mersenne prime 219937−1. The standard implementation of that, MT19937, uses a 32-bit word length. There is another implementation (with five variants[3]) 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: Dyalog APL,[4] Microsoft Excel,[5] GAUSS,[6] GLib,[7] GNU Multiple Precision Arithmetic Library,[8] GNU Octave,[9] GNU Scientific Library,[10] gretl,[11] IDL,[12] Julia,[13] CMU Common Lisp,[14] Embeddable Common Lisp,[15] Steel Bank Common Lisp,[16] Maple,[17] MATLAB,[18] Free Pascal,[19] PHP,[20] Python,[21][22][23] R,[24] Ruby,[25] SageMath,[26] Scilab,[27] Stata.[28]
It is also available in Apache Commons,[29] in standard C++ (since C++11),[30][31] and in Mathematica.[32] Add-on implementations are provided in many program libraries, including the Boost C++ Libraries,[33] the CUDA Library,[34] and the NAG Numerical Library.[35]
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".[36] The Mersenne Twister is similarly one of the PRNGs in SAS: the other generators are older and deprecated.[37] The Mersenne Twister is the default PRNG in Stata, the other one is KISS, for compatibility with older versions of Stata.[38]
Advantages
- Permissively-licensed and patent-free for all variants except CryptMT.
- Passes numerous tests for statistical randomness, including the Diehard tests and most, but not all of the TestU01 tests.[39]
- A very long period of 219937 − 1. Note that while a long period is not a guarantee of quality in a random number generator, short periods, such as the 232 common in many older software packages, can be problematic.[40]
- k-distributed to 32-bit accuracy for every 1 ≤ k ≤ 623 (for a definition of k-distributed, see below)
- Implementations generally create random numbers faster than true random methods. A study found that the Mersenne Twister creates 64-bit floating point random numbers approximately twenty times faster than the hardware-implemented, processor-based RDRAND instruction set.[41]
Disadvantages
- Relatively large state buffer, of 2.5 KiB, unless the TinyMT variant (discussed below) is used.
- Mediocre throughput by modern standards, unless the SFMT variant (discussed below) is used.[42]
- Exhibits two clear failures (linear complexity) in both Crush and BigCrush in the TestU01 suite. The test, like Mersenne Twister, is based on F2 algebra.[39] There are a number of other generators that pass all the tests (and numerous generators that fail badly)[clarification needed].
- Multiple 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.[43][44]
- Poor diffusion: 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.[45] The GPU version (MTGP) is said to be even better.[46]
- Contains subsequences with more 0's than 1's. This adds to the poor diffusion property to make recovery from many-zero states difficult.
- Is not cryptographically secure, unless the CryptMT variant (discussed below) is used. 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.
Alternatives
An alternative generator, WELL ("Well Equidistributed Long-period Linear"), offers quicker recovery, and equal randomness, and nearly equal speed.[47]
Marsaglia's xorshift generators and variants are the fastest in the class of LFSRs.[48]
64-bit MELGs ("64-bit Maximally Equidistributed F2-Linear Generators with Mersenne Prime Period") are completely optimized in terms of the k-distribution properties.[49]
The ACORN family (published 1989) is another k-distributed PRNG, which shows similar computational speed to MT, and better statistical properties as it satisfies all the current (2019) TestU01 criteria; when used with appropriate choices of parameters, ACORN can have arbitrarily long period and precision.
The PCG family is a more modern long-period generator, with better cache locality, and less detectable bias using modern analysis methods.Cite error: A <ref>
tag is missing the closing </ref>
(see the help page). (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 F2 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 2nw − 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 In − 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 F2, and therefore bitwise XOR takes the place of addition)
where x0 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 where for 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 = 9908B0DF16
- (u, d) = (11, FFFFFFFF16)
- (s, b) = (7, 9D2C568016)
- (t, c) = (15, EFC6000016)
- l = 18
Note that 32-bit implementations of the Mersenne Twister generally have d = FFFFFFFF16. 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:[50]
- (w, n, m, r) = (64, 312, 156, 31)
- a = B5026F5AA96619E916
- (u, d) = (29, 555555555555555516)
- (s, b) = (17, 71D67FFFEDA6000016)
- (t, c) = (37, FFF7EEE00000000016)
- l = 43
Initialization
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 x0 through xn − 1 by setting x0 to the seed value and thereafter setting
- xi = f × (xi−1 ⊕ (xi−1 >> (w−2))) + i
for i from 1 to n−1. The first value the algorithm then generates is based on xn, not on x0. 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.[50]
Comparison with classical GFSR
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
The twist transformation improves the classical GFSR with the following key properties:
- The 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 five dimensions)
Pseudocode
The following 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[51] } 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 }
Variants
CryptMT
CryptMT is a stream cipher and cryptographically secure pseudorandom number generator which uses Mersenne Twister internally.[52][53] It was developed by Matsumoto and Nishimura alongside Mariko Hagita and Mutsuo Saito. It has been submitted to the eSTREAM project of the eCRYPT network.[52] Unlike Mersenne Twister or its other derivatives, CryptMT is patented.
MTGP
MTGP is a variant of Mersenne Twister optimised for graphics processing units published by Mutsuo Saito and Makoto Matsumoto.[54] 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 very old GPU (Nvidia GTX260 with 192 cores) of 4.7 ms for 5×107 random 32-bit integers.
SFMT
![]() | This section needs expansion. You can help by adding to it. (June 2007) |
The SFMT (SIMD-oriented Fast Mersenne Twister) is a variant of Mersenne Twister, introduced in 2006,[55] designed to be fast when it runs on 128-bit SIMD.
- It is roughly twice as fast as Mersenne Twister.[56]
- 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 2607−1 to 2216091−1.
Intel SSE2 and PowerPC AltiVec are supported by SFMT. It is also used for games with the Cell BE in the PlayStation 3.[57]
TinyMT
TinyMT is a variant of Mersenne Twister, proposed by Saito and Matsumoto in 2011.[58] TinyMT uses just 127 bits of state space, a significant decrease compared to the original's 2.5 KiB of state. However, it has a period of 2127−1, far shorter than the original, so it is only recommended by the authors in cases where memory is at a premium.
References
- ^ E.g. Marsland S. (2011) Machine Learning (CRC Press), §4.1.1. Also see the section "Adoption in software systems".
- ^ Matsumoto, M.; Nishimura, T. (1998). "Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator" (PDF). ACM Transactions on Modeling and Computer Simulation. 8 (1): 3–30. CiteSeerX 10.1.1.215.1141. doi:10.1145/272991.272995. S2CID 3332028.
- ^ John Savard. "The Mersenne Twister".
A subsequent paper, published in the year 2000, gave five additional forms of the Mersenne Twister with period 2^19937-1. All five were designed to be implemented with 64-bit arithmetic instead of 32-bit arithmetic.
- ^ "Random link". Dyalog Language Reference Guide. Retrieved 2020-06-04.
- ^ Mélard, G. (2014), "On the accuracy of statistical procedures in Microsoft Excel 2010", Computational Statistics, 29 (5): 1095–1128, CiteSeerX 10.1.1.455.5508, doi:10.1007/s00180-014-0482-5, S2CID 54032450.
- ^ GAUSS 14 Language Reference
- ^ Random Numbers: GLib Reference Manual
- ^ "Random Number Algorithms". GNU MP. Retrieved 2013-11-21.
- ^ "16.3 Special Utility Matrices". GNU Octave.
Built-in Function: rand
- ^ "Random number environment variables". GNU Scientific Library. Retrieved 2013-11-24.
- ^ "uniform". Gretl Function Reference.
- ^ "RANDOMU (IDL Reference)". Exelis VIS Docs Center. Retrieved 2013-08-23.
- ^ "Random Numbers". Julia Language Documentation—The Standard Library.
- ^ "Design choices and extensions". CMUCL User's Manual. Retrieved 2014-02-03.
- ^ "Random states". The ECL manual. Retrieved 2015-09-20.
- ^ "Random Number Generation". SBCL User's Manual.
- ^ "random number generator". Maple Online Help. Retrieved 2013-11-21.
- ^ "Random number generator algorithms". Documentation Center, MathWorks.
- ^ "random". free pascal documentation. Retrieved 2013-11-28.
- ^ "mt_rand — Generate a better random value". PHP Manual. Retrieved 2016-03-02.
- ^ "9.6 random — Generate pseudo-random numbers". Python v2.6.8 documentation. Retrieved 2012-05-29.
- ^ "8.6 random — Generate pseudo-random numbers". Python v3.2 documentation. Retrieved 2012-05-29.
- ^ "random — Generate pseudo-random numbers — Python 3.8.3 documentation". Python 3.8.3 documentation. Retrieved 2020-06-23.
- ^ "Random Number Generators". CRAN Task View: Probability Distributions. Retrieved 2012-05-29.
- ^ ""Random" class documentation". Ruby 1.9.3 documentation. Retrieved 2012-05-29.
- ^ Probability Distributions — Sage Reference Manual v7.2: Probablity
- ^ "grand - Random numbers". Scilab Help.
- ^ New random-number generator—64-bit Mersenne Twister
- ^ "Data Generation". Apache Commons Math User Guide.
- ^ "Random Number Generation in C++11" (PDF). Standard C++ Foundation.
- ^ "std::mersenne_twister_engine". Pseudo Random Number Generation. Retrieved 2012-09-25.
- ^ [1] Mathematica Documentation
- ^ "boost/random/mersenne_twister.hpp". Boost C++ Libraries. Retrieved 2012-05-29.
- ^ "Host API Overview". CUDA Toolkit Documentation. Retrieved 2016-08-02.
- ^ "G05 – Random Number Generators". NAG Library Chapter Introduction. Retrieved 2012-05-29.
- ^ "Random Number Generators". IBM SPSS Statistics. Retrieved 2013-11-21.
- ^ "Using Random-Number Functions". SAS Language Reference. Retrieved 2013-11-21.
- ^ Stata help: set rng -- Set which random-number generator (RNG) to use
- ^ a b 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).
- ^ 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.
- ^ Route, Matthew (August 10, 2017). "Radio-flaring Ultracool Dwarf Population Synthesis". The Astrophysical Journal. 845 (1): 66. arXiv:1707.02212. Bibcode:2017ApJ...845...66R. doi:10.3847/1538-4357/aa7ede. S2CID 118895524.
{{cite journal}}
: CS1 maint: unflagged free DOI (link) - ^ "SIMD-oriented Fast Mersenne Twister (SFMT): twice faster than Mersenne Twister". Japan Society for the Promotion of Science. Retrieved 27 March 2017.
- ^ Makoto Matsumoto; Takuji Nishimura. "Dynamic Creation of Pseudorandom Number Generators" (PDF). Retrieved 19 July 2015.
- ^ Hiroshi Haramoto; Makoto Matsumoto; Takuji Nishimura; François Panneton; Pierre L’Ecuyer. "Efficient Jump Ahead for F2-Linear Random Number Generators" (PDF). Retrieved 12 Nov 2015.
- ^ "mt19937ar: Mersenne Twister with improved initialization". hiroshima-u.ac.jp. Retrieved 4 October 2015.
- ^ Fog, Agner (1 May 2015). "Pseudo-Random Number Generators for Vector Processors and Multicore Processors". Journal of Modern Applied Statistical Methods. 14 (1): 308–334. doi:10.22237/jmasm/1430454120.
- ^ P. L'Ecuyer, "Uniform Random Number Generators", International Encyclopedia of Statistical Science, Lovric, Miodrag (Ed.), Springer-Verlag, 2010.
- ^ "xorshift*/xorshift+ generators and the PRNG shootout".
- ^ Harase, S.; Kimoto, T. (2018). "Implementing 64-bit Maximally Equidistributed F2-Linear Generators with Mersenne Prime Period". ACM Transactions on Mathematical Software. 44 (3): 30:1–30:11. arXiv:1505.06582. doi:10.1145/3159444. S2CID 14923086.
- ^ a b "std::mersenne_twister_engine". Pseudo Random Number Generation. Retrieved 2015-07-20.
- ^ Takuji Nishimura; Makoto Matsumoto. "A C-program for MT19937, with initialization improved 2002/1/26". Retrieved 20 July 2015.
- ^ a b "CryptMt and Fubuki". eCRYPT. Retrieved 2017-11-12.
- ^ Matsumoto, Makoto; Nishimura, Takuji; Hagita, Mariko; Saito, Mutsuo (2005). "Cryptographic Mersenne Twister and Fubuki Stream/Block Cipher" (PDF).
- ^ Mutsuo Saito; Makoto Matsumoto (2010). "Variants of Mersenne Twister Suitable for Graphic Processors". arXiv:1005.4973v3 [cs.MS].
{{cite arXiv}}
: Cite has empty unknown parameter:|accessdate=
(help) - ^ "SIMD-oriented Fast Mersenne Twister (SFMT)". hiroshima-u.ac.jp. Retrieved 4 October 2015.
- ^ "SFMT:Comparison of speed". hiroshima-u.ac.jp. Retrieved 4 October 2015.
- ^ "PlayStation®3 License". scei.co.jp. Retrieved 4 October 2015.
- ^ "Tiny Mersenne Twister (TinyMT)". hiroshima-u.ac.jp. Retrieved 4 October 2015.
Further reading
- Harase, S. (2014), "On the 𝔽2-linear relations of Mersenne Twister pseudorandom number generators", Mathematics and Computers in Simulation, 100: 103–113, arXiv:1301.5435, doi:10.1016/j.matcom.2014.02.002, S2CID 6984431.
- Harase, S. (2019), "Conversion of Mersenne Twister to double-precision floating-point numbers", Mathematics and Computers in Simulation, 161: 76–83, arXiv:1708.06018, doi:10.1016/j.matcom.2018.08.006, S2CID 19777310.
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
- Mersenne Twister examples —a collection of Mersenne Twister implementations, in several programming languages - at GitHub
- SFMT in Action: Part I – Generating a DLL Including SSE2 Support – at Code Project