patternMinor
Why is the period of Mersenne Twister prng $2^{19937}-1$ and not $2^{19968}$?
Viewed 0 times
mersenneprngwhythe19968period19937andnottwister
Problem
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?
$32\cdot 624 = 19968$.
How come the exponent is 19937 instead?
Solution
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.
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
How does
Note that all bits in the array MT at the start of
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
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 > (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 > 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])
}Code Snippets
// 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
}// 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])
}Context
StackExchange Computer Science Q#109059, answer score: 4
Revisions (0)
No revisions yet.