Editing Mersenne Twister (section)

== Algorithmic detail ==
[[File:Mersenne_Twister_visualisation.svg|thumb|288x288px|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 <math>[0, 2^w-1]</math>.

The Mersenne Twister algorithm is based on a [[Recurrence relation|matrix linear recurrence]] over a finite [[Binary numeral system|binary]] [[Field (mathematics)|field]] <math>\textbf{F}_2</math>. The algorithm is a twisted [[generalised feedback shift register]]<ref>{{Cite journal|last1=Matsumoto|first1=M.|last2=Kurita|first2=Y.|year=1992|title=Twisted GFSR generators|url=http://ir.lib.hiroshima-u.ac.jp/00015037|journal=ACM Transactions on Modeling and Computer Simulation|volume=2|issue=3|pages=179–194|doi=10.1145/146382.146383|s2cid=15246234}}</ref> (twisted GFSR, or TGFSR) of [[rational normal form]] (TGFSR(R)), with state bit reflection and tempering. The basic idea is to define a series <math>x_i</math> through a simple recurrence relation, and then output numbers of the form <math>x_i^T</math>, where ''T'' is an invertible <math>\textbf{F}_2</math>-matrix called a [[Tempered representation|tempering matrix]].

The general algorithm is characterized by the following quantities:

* ''w'': word size (in number of bits)
* ''n'': degree of recurrence
* ''m'': middle word, an offset used in the recurrence relation defining the series <math>x</math>, <math>1 \le m < n</math>
* ''r'': separation point of one word, or the number of bits of the lower bitmask, <math>0 \le r \le w - 1</math>
* ''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 <math>2^{nw-r}-1</math> is a Mersenne prime. This choice simplifies the primitivity test and [[K-distribution|''k''-distribution]] test needed in the parameter search.

The series <math>x</math> is defined as a series of ''w''-bit quantities with the recurrence relation:

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

where <math>\mid</math> denotes [[concatenation]] of bit vectors (with upper bits on the left), <math>
\oplus
</math> the bitwise [[exclusive or]] (XOR), <math>
x_{k}^{u}
</math> means the upper {{nowrap|''w'' − ''r''}} bits of <math>
x_k
</math>, and <math>
x_{k+1}^{l}
</math> means the lower ''r'' bits of <math>
x_{k+1}
</math>.

The subscripts may all be offset by ''-n''

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

where now the LHS, <math>
x_k
</math>, is the next generated value in the series in terms of values generated in the past, which are on the RHS.

The twist transformation ''A'' is defined in rational normal form as:<math display="block">
A = \begin{pmatrix} 0 & I_{w - 1} \\ a_{w-1} & (a_{w - 2}, \ldots , a_0) \end{pmatrix} </math>
with <math>
I_{w-1}
</math> as the <math>
(w-1)(w-1)
</math> 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 <math>
\textbf{F}_{2}
</math>, and therefore bitwise XOR takes the place of addition)<math display="block">
\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>
x_0
</math> is the lowest order bit of <math>
x
</math>.

As like TGFSR(R), the Mersenne Twister is cascaded with a [[Tempered representation|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 <math>
A = T^{-1}*AT
</math> 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. This tempering is defined in the case of Mersenne Twister as

: <math>
\begin{aligned}
y &\equiv x \oplus ((x\gg u)~\And~d)\\
y &\equiv y \oplus ((y\ll s)~\And~b)\\
y &\equiv y \oplus ((y\ll t)~\And~c)\\
z &\equiv y \oplus (y\gg l)
\end{aligned}
</math>

where <math>x</math> is the next value from the series, <math>y</math> is a temporary intermediate value, and <math>z</math> is the value returned from the algorithm, with <math>\ll</math>and <math>\gg</math> as the [[Bit-shift|bitwise left and right shifts]], and <math>\&</math> 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 \left\lfloor{\frac{w}{2}}\right\rfloor - 1</math> is required to reach the upper bound of equidistribution for the upper bits.

The coefficients for MT19937 are:

<math>
\begin{aligned}
(w, n, m, r) &= (32, 624, 397, 31)\\
a &= \textrm{9908B0DF}_{16}\\
(u, d) &= (11, \textrm{FFFFFFFF}_{16})\\
(s, b) &= (7, \textrm{9D2C5680}_{16})\\
(t, c) &= (15, \textrm{EFC60000}_{16})\\
l &= 18\\
\end{aligned}
</math>

Note that 32-bit implementations of the Mersenne Twister generally have ''d''&nbsp;=&nbsp;FFFFFFFF<sub>16</sub>. As a result, the ''d'' is occasionally omitted from the algorithm description, since the bitwise [[Logical conjunction|and]] with ''d'' in that case has no effect.

The coefficients for MT19937-64 are:<ref name="std::mersenne_twister_engine">{{cite web|title=std::mersenne_twister_engine|url=http://en.cppreference.com/w/cpp/numeric/random/mersenne_twister_engine|access-date=2015-07-20|work=Pseudo Random Number Generation}}</ref>

<math>
\begin{aligned}
(w, n, m, r) = (64, 312, 156, 31)\\
a = \textrm{B5026F5AA96619E9}_{16}\\
(u, d) = (29, \textrm{5555555555555555}_{16})\\
(s, b) = (17, \textrm{71D67FFFEDA60000}_{16})\\
(t, c) = (37, \textrm{FFF7EEE000000000}_{16})\\
l = 43\\
\end{aligned}

</math>

=== 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 <math>x_0</math> through <math>x_{n-1}</math> by setting <math>x_0</math> to the seed value and thereafter setting

: <math>
x_i = f \times (x_{i-1} \oplus (x_{i-1} \gg (w-2))) + i
</math>

for <math>i</math> from <math>1</math> to <math>n-1</math>.

* The first value the algorithm then generates is based on <math>x_n</math>, not on <math>x_0</math>.
* The constant ''f'' forms another parameter to the generator, though not part of the algorithm proper.
* The value for ''f'' for MT19937 is 1812433253.
* The value for ''f'' for MT19937-64 is 6364136223846793005.<ref name="std::mersenne_twister_engine" />

=== C code ===

<syntaxhighlight lang="c">#include <stdint.h>

#define n 624
#define m 397
#define w 32
#define r 31
#define UMASK (0xffffffffUL << r)
#define LMASK (0xffffffffUL >> (w-r))
#define a 0x9908b0dfUL
#define u 11
#define s 7
#define t 15
#define l 18
#define b 0x9d2c5680UL
#define c 0xefc60000UL
#define f 1812433253UL

typedef struct
{
    uint32_t state_array[n];         // the array for the state vector 
    int state_index;                 // index into state vector array, 0 <= state_index <= n-1   always
} mt_state;


void initialize_state(mt_state* state, uint32_t seed) 
{
    uint32_t* state_array = &(state->state_array[0]);
    
    state_array[0] = seed;                          // suggested initial seed = 19650218UL
    
    for (int i=1; i<n; i++)
    {
        seed = f * (seed ^ (seed >> (w-2))) + i;    // Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
        state_array[i] = seed; 
    }
    
    state->state_index = 0;
}


uint32_t random_uint32(mt_state* state)
{
    uint32_t* state_array = &(state->state_array[0]);
    
    int k = state->state_index;      // point to current state location
                                     // 0 <= state_index <= n-1   always
    
//  int k = k - n;                   // point to state n iterations before
//  if (k < 0) k += n;               // modulo n circular indexing
                                     // the previous 2 lines actually do nothing
                                     //  for illustration only
    
    int j = k - (n-1);               // point to state n-1 iterations before
    if (j < 0) j += n;               // modulo n circular indexing

    uint32_t x = (state_array[k] & UMASK) | (state_array[j] & LMASK);
    
    uint32_t xA = x >> 1;
    if (x & 0x00000001UL) xA ^= a;
    
    j = k - (n-m);                   // point to state n-m iterations before
    if (j < 0) j += n;               // modulo n circular indexing
    
    x = state_array[j] ^ xA;         // compute next value in the state
    state_array[k++] = x;            // update new state value
    
    if (k >= n) k = 0;               // modulo n circular indexing
    state->state_index = k;
    
    uint32_t y = x ^ (x >> u);       // tempering 
             y = y ^ ((y << s) & b);
             y = y ^ ((y << t) & c);
    uint32_t z = y ^ (y >> l);
    
    return z; 
}</syntaxhighlight>

=== Comparison with classical GFSR ===

In order to achieve the <math>2^{nw-r}-1</math> theoretical upper limit of the period in a T[[GFSR]], <math>\phi_{B}(t)</math> must be a [[Primitive polynomial (field theory)|primitive polynomial]], <math>\phi_{B}(t)</math> 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\text{-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:

* The period reaches the theoretical upper limit <math>2^{nw-r}-1</math> (except if initialized with 0)
* Equidistribution in ''n'' dimensions (e.g. [[Linear congruential generator|linear congruential generators]] can at best manage reasonable distribution in five dimensions)