Why is the period of Mersenne Twister prng $2^{19937}-1$ and not $2^{19968}$?

Question

For the Mersenne Twister prng, wikipedia, I feel the period exponent should be the size of the MT array, $32\cdot 624 = 19968$.

How come the exponent is 19937 instead?

Answer 1

Here is the one-line explanation.

The number of effective bits in MT19937 is 624*32-31=19937. They cannot be all 0. So the period of MT19937 is $2^{19937}-1$.

---

All notations if not introduced below come from this Wikipedia article. The parameters for MT19937, (w, n, m, r) = (32, 624, 397, 31) will be used from time to time.

Explanation from the recurrence

The series $x$ is defined as a series of w-bit quantities with the recurrence relation: $$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$$

Note the lower part of the lowest term $x_k$, ${x_k}^l$ is not used in the recurrence relation. Hence the number of effective bits used by the recurrence is $(k+n-k)w-r=nw-r$. Since 0 is excluded, the theoretical limit of period of Mersernne Twister is $2^{nw-r}-1=2^{19937}-1$.

Detailed explanation.

The above explanation is not rigorous enough because of a lack of proper reference. To convince myself, I followed the established pseudocode implementation.

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

The code above shows that the internal state of a Mersenne Twister is kept in MT[0], MT[1], ..., MT[n-1], which is an array of n values of w bits each. It seems that the period of a Mersenne Twister should be $2^{nw}=2^{19968}$.

Let us check the situation closely.

Every call to twist() will update the internal state. After each call to twist(), function extract_number() will be called $n$ times, each time extracting from MT[index] one word of w bits to output where index is one of 0, 1, ..., n-1.

How does twist() update the internal state?

• the upper w-r bits of MT[0], the lower r bits of MT[1], and all bits of MT[317] are used to update MT[0].
• the upper w-r bits of MT[1], the lower r bits of MT[2], and all bits of MT[318] are used to update MT[1].
• the upper w-r bits of MT[2], the lower r bits of MT[3], and all bits of MT[319] are used to update MT[2].
• ...
• the upper w-r bits of MT[306], the lower r bits of MT[307], and all bits of MT[623] are used to update MT[306].
• the upper w-r bits of MT[307], the lower r bits of MT[308], and all bits of the updated MT[0] are used to update MT[307].
• the upper w-r bits of MT[308], the lower r bits of MT[309], and all bits of the updated MT[1] are used to update MT[308].
• ...
• the upper w-r bits of MT[622], the lower r bits of MT[623], and all bits of the updated MT[315] are used to update MT[622].
• the upper w-r bits of MT[623], the lower r bits of updated MT[0], and all bits of the updated MT[316] are used to update MT[623].

Note that all bits in the array MT at the start of twist() are used except the lower r bits of MT[0]. So an effective internal state of MT is r bits less than all bits in array MT. Since an effective internal state can never be all 0, the number of effective internal states is at most $2^{nw-r}-1=2^{19937}-1$, which is the theoretical limit of the period for a linear-feedback shift register like Mersenne Twister.

It is another level of interesting mathematical reasoning how MT19937 attains that theoretical limit of its period.

Another explanation by pseudocode

The following pseudocode is another implementation of Mersenne Twister.

Note the code MT[0] := MT[0] and upper_mask and MT[0] := MT[1] and upper_mask, which limits the number of effective bits in MT[0] to w-r. There are (n-1)w bits in MT[1], MT[2], ..., MT[n-1]. In total, there are (w-r)+(n-1)w = nw-r bits in the (effective) internal state of a Mersenne Twister.

// Adapted from https://en.wikipedia.org/w/index.php?title=Mersenne_Twister&oldid=896072298#Pseudocode by Apass.Jack

// 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, which must be done once before repeated calls to next(). 
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)
    }
    MT[0] := MT[0] and upper_mask
}

function temper(int y) {
    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)

    return lowest w bits of (y)
}

function twist() {
    int x := MT[0]
              + (MT[1] and lower_mask)
    int xA := x >> 1
    if (x mod 2) != 0 { // lowest bit of x is 1
        xA := xA xor a
    }
    x := MT[m] xor xA

    MT[0] := MT[1] and upper_mask
    for i from 1 to (n-2) {
      MT[i] := MT[i+1]
    }
    MT[n-1] := x
}

// Generate the next value from the series MT[]  
function next() {
    twist()
    return temper(MT[n-1])
}