View source for Mersenne twister

[[:en:Mersenne twister]]
{{TOCleft}}
== My Draft ==

A version of this adapted as an example toy '''ai-integration.biz''' application is the purpose of the draft here.

== Overview ==

The '''Mersenne twister''' is a [[pseudorandom number generator]] developed in 1997 by {{nihongo|[[Makoto Matsumoto]]|松本 眞}} and {{nihongo|[[Takuji Nishimura]]|西村 拓士}}<ref>{{cite doi|10.1145/272991.272995}}</ref> that is based on a [[matrix linear recurrence]] over a finite [[binary numeral system|binary]] [[field (mathematics)|field]] <math>F_{2}</math>. 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.

== Application ==

The algorithm in its native form is not suitable for [[cryptography]] (unlike [[Blum Blum Shub]]). Observing a sufficient number of iterates (624 in the case of MT19937) allows one to predict all future iterates. A  pair of cryptographic stream ciphers based on output from Mersenne twister has been proposed by Makoto Matsumoto et&nbsp;al. The authors claim speeds 1.5 to 2 times faster than [[AES]] in [[counter mode]].<ref>{{citation|title=Cryptographic Mersenne Twister and Fubuki Stream/Block Cipher|last1=Matsumoto|first1=Makoto|last2=Nishimura|first2=Takuji|last3=Hagita|first3=Mariko|last4=Saito|first4=Mutsuo|year=2005|url=http://eprint.iacr.org/2005/165.pdf}}</ref>

Another issue is that it can take a long time to turn a non-random initial state into output that passes [[randomness tests]]. A small [[lagged Fibonacci generator]] or [[linear congruential generator]] gets started much quicker and usually is used to seed the Mersenne Twister. If only a few numbers are required and standards aren't high, it is simpler to use the seed generator. But the Mersenne Twister will still work.

For many applications, however, the Mersenne twister is quickly becoming the pseudorandom number generator of choice; for example, it is the default in [[R (programming language)|R]], [[Maple (programming language)|Maple]],  [[MATLAB]], and the two popular scripting languages [[Python (programming language)|Python]]<ref>{{cite web|title=10.6. random — Generate pseudo-random numbers|work=Python v2.6.2 documentation|url=http://docs.python.org/library/random.html#module-random|accessdate=2009-04-27}}</ref> and [[Ruby (programming language|Ruby]]<ref>{{cite web|title=Module: Kernel|work=Ruby documentation|url=http://www.ruby-doc.org/core/classes/Kernel.html#M005977|accessdate=2009-04-27}}</ref>. Since the library is portable, freely available, and quickly generates high quality pseudorandom numbers, it is rarely a bad choice.

It is designed with [[Monte Carlo simulation]]s and other statistical simulations in mind. Researchers primarily want high quality numbers but also benefit from its speed and portability.

==Advantages==

The commonly used variant of Mersenne Twister, MT19937 has the following desirable properties:

# It was designed to have a period of 2<sup>19937</sup>&nbsp;&minus;&nbsp;1 (the creators of the algorithm proved this property). In practice, there is little reason to use a larger period, as most applications do not require 2<sup>19937</sup> unique combinations (2<sup>19937</sup> is approximately 4.3&nbsp;×&nbsp;10<sup>6001</sup>; this is many orders of magnitude larger than the estimated number of particles in the [[observable universe#Matter content|observable universe]], which is 10<sup>87</sup>).
# It has a very high order of dimensional [[equidistributed|equidistribution]] (see [[linear congruential generator]]). This implies that there is negligible serial correlation between successive values in the output sequence.
# 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.

== Alternatives ==

The Mersenne Twister algorithm has received some criticism in the computer science field, notably by [[George Marsaglia]]. These critics claim that while it is good at generating random numbers, it is not very elegant and is overly complex to implement. Marsaglia has provided several examples of random number generators that are less complex yet which he claims provide significantly larger periods. For example, a simple complementary [[multiply-with-carry]] generator can have a period 10<sup>33000</sup> times as long, be significantly faster, and maintain better or equal randomness.<ref>[http://groups.google.com/group/comp.lang.c/browse_thread/thread/a9915080a4424068/ Marsaglia on Mersenne Twister 2003]</ref><ref>[http://groups.google.com/group/sci.crypt/browse_thread/thread/305c507efbe85be4 Marsaglia on Mersenne Twister 2005]</ref>

== Algorithmic detail ==

The Mersenne Twister algorithm is a twisted [[generalised feedback shift register]]<ref>{{cite doi|10.1145/146382.146383}}</ref> (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 ≤ ''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'', ''l'': additional Mersenne Twister tempering bit shifts

with the restriction that 2<sup>''nw''&nbsp;−&nbsp;''r''</sup>&nbsp;−&nbsp;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

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

with | as the bitwise [[Logical disjunction|or]] and ⊕ as the bitwise [[exclusive or]] (XOR), '''''x'''''<sup>''u''</sup>, '''''x'''''<sup>''l''</sup> being '''''x''''' with upper and lower bitmasks applied. The twist transformation ''A'' is defined in rational normal form

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

with ''I''<sub>''n''&nbsp;&minus;&nbsp;1</sub> as the (''n''&nbsp;&minus;&nbsp;1)&nbsp;&times;&nbsp;(''n''&nbsp;&minus;&nbsp;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

<math>
\boldsymbol{x}A = \begin{cases}\boldsymbol{x} \gg 1 & x_0 = 0\\(\boldsymbol{x} \gg 1) \oplus \boldsymbol{a} & x_0 = 1\end{cases}
</math>

where
:<math>\boldsymbol{x} := ({x_k}^u \mid {x_{k+1}}^l) \qquad \qquad k=0,1,\ldots</math>

In order to achieve the 2<sup>''nw''&nbsp;−&nbsp;''r''</sup>&nbsp;−&nbsp;1 theoretical upper limit of the period in a TGFSR, ''φ''<sub>''B''</sub>(''t'') must be a [[primitive polynomial]], ''φ''<sub>''B''</sub>(''t'') being the [[characteristic polynomial]] of

<math>
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}
</math>

<math>
S = \begin{pmatrix} 0 & I_{r} \\ I_{w - r} & 0 \end{pmatrix} A
</math>

The twist transformation improves the classical GFSR with the following key properties:
* Period reaches the theoretical upper limit 2<sup>''nw''&nbsp;−&nbsp;''r''</sup>&nbsp;−&nbsp;1 (except if initialized with 0)
* Equidistribution in ''n'' dimensions (e.g. [[linear congruential generator]]s 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''&nbsp;=&nbsp;''R'' → ''A''&nbsp;=&nbsp;''T''<sup>−1</sup>''RT'', ''T'' invertible. The tempering is defined in the case of Mersenne Twister as

:'''''y''''' := '''''x'''''&nbsp;&oplus;&nbsp;('''''x'''''&nbsp;>>&nbsp;''u'')
:'''''y''''' := :'''''y'''''&nbsp;&oplus;&nbsp;(('''''y'''''&nbsp;<<&nbsp;''s'')&nbsp;&amp;&nbsp;'''''b''''')
:'''''y''''' := :'''''y'''''&nbsp;&oplus;&nbsp;(('''''y'''''&nbsp;<<&nbsp;''t'')&nbsp;&amp;&nbsp;'''''c''''')
:'''''z''''' := '''''y'''''&nbsp;&oplus;&nbsp;('''''y'''''&nbsp;>>&nbsp;''l'')

with <<, >> as the bitwise left and right shifts, and & as the bitwise [[Logical conjunction|and]]. The first and last transforms are added in order to improve lower bit equidistribution. From the property of TGFSR, <math>s + t \ge \lfloor w/2 \rfloor - 1</math> is required to reach the upper bound of equidistribution for the upper bits.

The coefficients for MT19937 are:
* (''w'', ''n'', ''m'', ''r'')&nbsp;=&nbsp;(32, 624, 397, 31)
* '''''a'''''&nbsp;=&nbsp;9908B0DF<sub>16</sub>
* ''u''&nbsp;=&nbsp;11
* (''s'', '''''b''''')&nbsp;=&nbsp;(7, 9D2C5680<sub>16</sub>)
* (''t'', '''''c''''')&nbsp;=&nbsp;(15, EFC60000<sub>16</sub>)
* ''l''&nbsp;=&nbsp;18

==Pseudocode==

The following generates uniformly 32-bit integers in the range [0, 2<sup>32</sup>&nbsp;&minus;&nbsp;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''' initializeGenerator('''int''' seed) {
      MT[0] := seed
      '''for''' i '''from''' 1 '''to''' 623 { ''// loop over each other element''
          MT[i] := '''last 32 bits of'''(1812433253 * (MT[i-1] '''[[Bitwise operation#XOR|xor]]''' ('''right shift by 30 bits'''(MT[i-1]))) + i) ''// 0x6c078965''
      }
  }
  
  ''// Extract a tempered pseudorandom number based on the index-th value,''
  ''// calling generateNumbers() every 624 numbers''
  '''function''' extractNumber() {
      '''if''' index == 0 {
          generateNumbers()
      }
      
      '''int''' y := MT[index]
      y := y '''xor''' ('''right shift by 11 bits'''(y))
      y := y '''xor''' ('''left shift by 7 bits'''(y) '''[[Bitwise operation#AND|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) '''[[modulo operation|mod]]''' 624
      '''return''' y
  }
  
  ''// Generate an array of 624 untempered numbers''
  '''function''' generateNumbers() {
      '''for''' i '''from''' 0 '''to''' 623 {
          '''int''' y := '''32nd bit of'''(MT[i]) + '''last 31 bits of'''(MT[(i+1) '''mod''' 624])
          MT[i] := MT[(i + 397) '''mod''' 624] '''xor''' ('''right shift by 1 bit'''(y))
          '''if''' (y '''mod''' 2) == 1 { ''// y is odd''
              MT[i] := MT[i] '''xor''' (2567483615) ''// 0x9908b0df''
          }
      }
  }

== SFMT ==
{{Expand-section|date=June 2007}}

SFMT, the [[SIMD]]-oriented Fast Mersenne Twister, is a variant of Mersenne Twister, introduced in 2006<ref>[http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/index.html SIMD-oriented Fast Mersenne Twister (SFMT)<!-- Bot generated title -->]</ref>, designed to be fast when it runs on 128-bit SIMD. 
*It is roughly twice as fast as Mersenne Twister.<ref>[http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/speed.html SFMT:Comparison of speed<!-- Bot generated title -->]</ref>
*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<sup>607</sup>-1 to 2<sup>216091</sup>-1.

Intel [[SSE2]] and [[PowerPC]] AltiVec are supported by SFMT. It is also used for games with the [[Cell (microprocessor)|Cell BE]] in the [[Playstation 3]].<ref>[http://www.scei.co.jp/ps3-license/index.html PLAYSTATION 3 License<!-- Bot generated title -->]</ref>

=== Implementations ===
* [http://www.pimpmyexcel.com/ XLL Excel Addin Implementation of mersenne twister random number generator]
* [http://www.cs.gmu.edu/~sean/research/ Two implementations of Mersenne Twister in Java: one is the fastest known, and the other is a drop-in replacement for java.util.Random]
* [http://www.gnu.org/software/gsl/ The GNU Scientific Library (GSL), containing an implementation of the Mersenne Twister]
* [http://www.agner.org/random/ C++ and binary function libraries for several platforms. Multithreaded. Includes Mersenne Twister and SFMT]
* [http://www.cs.hmc.edu/~geoff/mtwist.html Implementations of the Mersenne Twister in C and C++]
* [http://www-personal.engin.umich.edu/~wagnerr/MersenneTwister.html Implementation of the Mersenne Twister in C++]
* [http://babel.isa.uma.es/mrpt/index.php/Example:Random_number_generation Implementation of the Mersenne Twister in C++] within the [[Mobile Robot Programming Toolkit]]
* [http://www.numtech.com/NtRand/ Implementation of Mersenne Twister as an add-in for Microsoft Excel]
* [http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/BASIC/basic.html Implementation of Mersenne Twister as a free module for Visual Basic (Microsoft Excel, Microsoft Access and VB compilers) and for other Basic versions in the official site of the Mersenne Twister]
* [http://www.aaronballman.com/programming/REALbasic/Rand.php Implementation of Mersenne Twister for REALbasic (requires REALbasic 2006r1 or greater)]
* [http://cybertiggyr.com/gene/jmt/ Implementation of Mersenne Twister for Lisp]
* [http://www.rapideuphoria.com/mt.zip Implementation of Mersenne Twister in Euphoria]
* [http://code.msdn.microsoft.com/MersenneTwister Implementation of Mersenne Twister for C# (newer, System.Random drop-in replacement)] ([http://www.c-sharpcorner.com/Code/2003/April/MersenneTwisterAlgo.asp Older implementation])
* [http://adrianhoe.com/adrianhoe/projects/adamt19937/ Implementation of Mersenne Twister for Ada]
* [http://www.coyotegulch.com/products/libcoyotl/twisted_road/index.html Implementation of Mersenne Twister for Fortran 95]
* [http://modp.com/release/mma_mersenne_twister/ Implementation of Mersenne Twister for Mathematica]
* [http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=6614&objectType=file Implementation of Mersenne Twister for MATLAB]
* [http://www.mitrionics.com/?page=page49058d1d215e8 Implementation of Mersenne Twister for Mitrion-C]
* [http://clean.cs.ru.nl/Download/Download_Libraries/mt/body_mt.html Implementation of Mersenne Twister for Clean]
* [http://herbert.gandraxa.com/herbert/rng.asp High-speed Implementation of Mersenne Twister] in [[Linoleum programming language|Linoleum]] (a cross-platform [[Assembly language|Assembler]]), by Herbert Glarner
* [http://search.cpan.org/search?module=Math%3A%3ARandom%3A%3AMT%3A%3AAuto CPAN module implementing the Mersenne Twister for use with Perl]
* [http://www.augustsson.net/Darcs/MT/ Implementation of Mersenne Twister for Haskell]
* [http://mlton.org/cgi-bin/viewsvn.cgi/mltonlib/trunk/org/mlton/ville/mersenne-twister/unstable/ Implementation of Mersenne Twister for Standard ML]
* [http://mbishop.esoteriq.org/code/Fsharp/mt.fs Implementation of Mersenne Twister in F#]
* It also is implemented in [[gLib]] and the standard libraries of at least [[PHP]], [[Python (programming language)|Python]] and [[Ruby programming language|Ruby]].
* [http://charles.karney.info/random C++ class implementing Mersenne Twister and SFMT]
* [http://adam.ierymenko.name/files/MersenneTwister32_spe.cpp C++ implementation of Mersenne Twister for the IBM/Sony Cell Broadband Engine (Cell BE) specialized processing units]
* [http://flashexperiments.insh-allah.com/#Mersenne_Twister_ported_to_ActionScript Mersenne Twister ported to ActionScript]

==References==
{{Reflist|2}}

== External links ==
* [http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/earticles.html  The academic paper for MT, and related articles by Makoto Matsumoto]
* [http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html Mersenne Twister home page, with codes in C, Fortran, Java, Lisp and some other languages]
* [http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/index.html SIMD-oriented Fast Mersenne Twister (SFMT)]

[[Category:Pseudorandom number generators]]
[[Category:Articles with example pseudocode]]

[[de:Mersenne-Twister]]
[[fr:Mersenne Twister]]
[[ko:메르센 트위스터]]
[[it:Mersenne Twister]]
[[nl:Mersenne-twister]]
[[ja:メルセンヌ・ツイスタ]]
[[ru:Вихрь Мерсенна]]
[[pl:Mersenne Twister]]