Algorithmic detail Mersenne Twister
visualisation of generation of pseudo-random 32-bit integers using mersenne twister. extract number section shows example integer 0 has been output , index @ integer 1. generate numbers run when integers have been output.
for w-bit word length, mersenne twister generates integers in range [0, 2−1].
the mersenne twister algorithm based on matrix linear recurrence on finite binary field f2. algorithm twisted generalised feedback shift register (twisted gfsr, or tgfsr) of rational normal form (tgfsr(r)), state bit reflection , tempering. basic idea define series
x
i
{\displaystyle x_{i}}
through simple recurrence relation, , output numbers of form
x
i
t
{\displaystyle x_{i}t}
,
t
{\displaystyle t}
invertible f2 matrix called tempering matrix.
the general algorithm characterized following quantities (some of these explanations make sense after reading rest of algorithm):
w: word size (in number of bits)
n: degree of recurrence
m: middle word, offset used in recurrence relation defining series x, 1 ≤ m < n
r: separation point of 1 word, or number of bits of lower bitmask, 0 ≤ r ≤ w - 1
a: coefficients of 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 restriction 2 − 1 mersenne prime. choice simplifies primitivity test , k-distribution test needed in parameter search.
the series x defined series of w-bit quantities recurrence relation:
x
k
+
n
:=
x
k
+
m
⊕
(
(
x
k
u
∣∣
x
k
+
1
l
)
a
)
k
=
0
,
1
,
…
{\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 left),
⊕
{\displaystyle \oplus }
bitwise exclusive or (xor),
x
k
u
{\displaystyle {x_{k}}^{u}}
means upper
w
−
r
{\displaystyle w-r}
bits of
x
k
{\displaystyle x_{k}}
, ,
x
k
+
1
l
{\displaystyle x_{k+1}^{l}}
means lower
r
{\displaystyle r}
bits of
x
k
+
1
{\displaystyle x_{k+1}}
. twist transformation defined in rational normal form as:
a
=
(
0
i
w
−
1
a
w
−
1
(
a
w
−
2
,
…
,
a
0
)
)
{\displaystyle a={\begin{pmatrix}0&i_{w-1}\\a_{w-1}&(a_{w-2},\ldots ,a_{0})\end{pmatrix}}}
with in − 1 (n − 1) × (n − 1) identity matrix. rational normal form has benefit multiplication can efficiently expressed as: (remember here matrix multiplication being done in f2, , therefore bitwise xor takes place of addition)
x
a
=
{
x
≫
1
x
0
=
0
(
x
≫
1
)
⊕
a
x
0
=
1
{\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 x0 lowest order bit of x.
as tgfsr(r), mersenne twister cascaded tempering transform compensate reduced dimensionality of equidistribution (because of choice of being in rational normal form). note equivalent using matrix a = tat t invertible matrix, , therefore analysis of characteristic polynomial mentioned below still holds.
as a, choose tempering transform computable, , not construct t itself. tempering defined in 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 next value series, y temporary intermediate value, z value returned algorithm, <<, >> bitwise left , right shifts, , & bitwise and. first , last transforms added in order improve lower-bit equidistribution. property of tgfsr,
s
+
t
≥
⌊
w
/
2
⌋
−
1
{\displaystyle s+t\geq \lfloor w/2\rfloor -1}
required reach upper bound of equidistribution upper bits.
the coefficients mt19937 are:
(w, n, m, r) = (32, 624, 397, 31)
a = 9908b0df16
(u, d) = (11, ffffffff16)
(s, b) = (7, 9d2c568016)
(t, c) = (15, efc6000016)
l = 18
note 32-bit implementations of mersenne twister have d = ffffffff16. result, d omitted algorithm description, since bitwise , d in case has no effect.
the coefficients mt19937-64 are:
(w, n, m, r) = (64, 312, 156, 31)
a = b5026f5aa96619e916
(u, d) = (29, 555555555555555516)
(s, b) = (17, 71d67fffeda6000016)
(t, c) = (37, fff7eee00000000016)
l = 43
initialization
as should apparent above description, state needed mersenne twister implementation array of n values of w bits each. initialize array, w-bit seed value used supply x0 through xn − 1 setting x0 seed value , thereafter setting
xi = f × (xi-1 ⊕ (xi-1 >> (w-2))) + i
for 1 n-1. first value algorithm generates based on xn, not based on x0. constant f forms parameter generator, though not part of algorithm proper. value f mt19937 1812433253 , mt19937-64 6364136223846793005.
comparison classical gfsr
in order achieve 2 − 1 theoretical upper limit of period in tgfsr, φb(t) must primitive polynomial, φb(t) being characteristic polynomial of
b
=
(
0
i
w
⋯
0
0
⋮
i
w
⋮
⋱
⋮
⋮
⋮
0
0
⋯
i
w
0
0
0
⋯
0
i
w
−
r
s
0
⋯
0
0
)
←
m
-th row
{\displaystyle 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{\hbox{-th row}}\\\\\\\\\end{matrix}}}
s
=
(
0
i
r
i
w
−
r
0
)
a
{\displaystyle s={\begin{pmatrix}0&i_{r}\\i_{w-r}&0\end{pmatrix}}a}
the twist transformation improves classical gfsr following key properties:
the period reaches theoretical upper limit 2 − 1 (except if initialized 0)
equidistribution in n dimensions (e.g. linear congruential generators can @ best manage reasonable distribution in 5 dimensions)
^ matsumoto, m.; kurita, y. (1992). twisted gfsr generators . acm transactions on modeling , computer simulation. 2 (3): 179–194. doi:10.1145/146382.146383.
^ std::mersenne_twister_engine . pseudo random number generation. retrieved 2015-07-20.
^ std::mersenne_twister_engine . pseudo random number generation. retrieved 2015-07-20.
Comments
Post a Comment