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?
Computer Science Stack Exchange is a question and answer site for students, researchers and practitioners of computer science. It only takes a minute to sign up.
Sign up to join this communityFor 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?
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.
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$.
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?
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.
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])
}