Mersenne Twister

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]

• 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 2¹⁹⁹³⁷ − 1. Note that while a long period is not a guarantee of quality in a random number generator, short periods, such as the 2³² 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]

• 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 F₂ algebra.[39] There are a number of other generators that pass all the tests (and numerous generators that fail badly).
• 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.

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 F₂-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.[50]

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

Let trunc_v(x) denote the number formed by the leading v bits of x, and consider P of the k v-bit vectors

${\displaystyle ({\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 2^kv possible combinations of bits occurs the same number of times in a period, except for the all-zero combination that occurs once less often.

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 [0, 2^w−1].

The Mersenne Twister algorithm is based on a matrix linear recurrence over a finite binary field F₂. The algorithm is a twisted generalised feedback shift register[51] (twisted GFSR, or TGFSR) of rational normal form (TGFSR(R)), with state bit reflection and tempering. The basic idea is to define a series ${\displaystyle x_{i}}$ through a simple recurrence relation, and then output numbers of the form ${\displaystyle x_{i}T}$, where ${\displaystyle T}$ is an invertible F₂ 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 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.

The series x is defined as a series of w-bit quantities with the recurrence relation:

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

where ${\displaystyle \mid \mid }$ denotes concatenation of bit vectors (with upper bits on the left), ${\displaystyle \oplus }$ the bitwise exclusive or (XOR), ${\displaystyle {x_{k}}^{u}}$ means the upper ${\displaystyle w-r}$ bits of ${\displaystyle x_{k}}$, and ${\displaystyle x_{k+1}^{l}}$ means the lower ${\displaystyle r}$ bits of ${\displaystyle x_{k+1}}$. The twist transformation A is defined in rational normal form as:

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

with I_n − 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 F₂, and therefore bitwise XOR takes the place of addition)

${\displaystyle {\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 x₀ 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 ${\displaystyle A}$ where ${\displaystyle A=T^{-1}AT}$ for ${\displaystyle 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. 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, ${\displaystyle s+t\geq \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, d) = (11, FFFFFFFF₁₆)
• (s, b) = (7, 9D2C5680₁₆)
• (t, c) = (15, EFC60000₁₆)
• l = 18

Note that 32-bit implementations of the Mersenne Twister generally have d = FFFFFFFF₁₆. 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:[52]

• (w, n, m, r) = (64, 312, 156, 31)
• a = B5026F5AA96619E9₁₆
• (u, d) = (29, 5555555555555555₁₆)
• (s, b) = (17, 71D67FFFEDA60000₁₆)
• (t, c) = (37, FFF7EEE000000000₁₆)
• l = 43

 // 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[53]
         }
         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
 }

from warnings import warn


class Mersenne:
    """Pseudorandom number generator"""

    def __init__(self, seed=1234):
        """
        Initialize pseudorandom number generator. Accepts an
        integer or floating-point seed, which is used in
        conjunction with an integer multiplier, k, and the
        Mersenne prime, j, to "twist" pseudorandom numbers out
        of the latter. This member also initializes the order
        of the generator's period, so that members floating and
        integer can emit a warning when generation is about to
        cycle and thus become not so pseudorandom.
        """
        self.seed = seed
        self.j = 2 ** 31 - 1
        self.k = 16807
        self.period = 2 ** 30

    def floating(self, interval=None, count=1):
        """
        Return a pseudorandom float. Default is one floating-
        point number between zero and one. Pass in a tuple or
        list, (a,b), to return a floating-point number on
        [a,b]. If count is 1, a single number is returned,
        otherwise a list of numbers is returned.
        """
        results = []
        for i in range(count):
            self.seed = (self.k * self.seed) % self.j
            if interval is not None:
                results.append(
                    (interval[1] - interval[0]) * (self.seed / self.j) + interval[0]
                )
            else:
                results.append(self.seed / self.j)
            self.period -= 1
            if self.period == 0:
                warn("Pseudorandom period nearing!!", category=ResourceWarning)
                self.period = 2 ** 30 # reset period
        if count == 1:
            return results.pop()
        else:
            return results

    def integer(self, interval=None, count=1):
        """
        Return a pseudorandom integer. Default is one integer
        number in {0,1}. Pass in a tuple or list, (a,b), to
        return an integer number on [a,b]. If count is 1, a
        single number is returned, otherwise a list of numbers
        is returned.
        """
        results = []
        for i in range(count):
            self.seed = (self.k * self.seed) % self.j
            if interval is not None:
                results.append(
                    int(
                        (interval[1] - interval[0] + 1) * (self.seed / self.j)
                        + interval[0]
                    )
                )
            else:
                result = self.seed / self.j
                if result < 0.50:
                    results.append(0)
                else:
                    results.append(1)
            self.period -= 1
            if self.period == 0:
                warn("Pseudorandom period nearing!!", category=ResourceWarning)
                self.period = 2 ** 30 # reset period
        if count == 1:
            return results.pop()
        else:
            return results

• Harase, S. (2014), "On the 𝔽₂-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.

• 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