Mersenne Twisters

The Mersenne Twister1 is so far the most widely used PRNG. Mersenne Twisters are taken as default random number generators of a great many of software systems, including the Julia language until the current 1.0 version. The most commonly used version of Mersenne Twisters is MT19937, which has a very long period of $2^{19937}-1$ and passes numerous tests including the Diehard tests.

However, it also has many flaws by today's standards. For the large period, MT19937 has to use a 2.5 KiB state buffer, which place a load on the memory caches. More severely, it cannot pass all the TestU01 statistical tests2, and the speed is not so fast. So it is not recommended in most situations.

The MersenneTwisters in this package currently only provides one RNG: MT19937. MT19937 can only produce UInt32 numbers as output. Its state is an array of 624 UInt32 numbers, so it takes 624 UInt32s as seed. A default function is also provided to deal with one UInt32 as seed.

Examples

To use the Mersenne Twisters, firstly import the module:

julia> using RandomNumbers.MersenneTwisters


A certain sequence can be used to initialize an instance of MT19937:

julia> seed = Tuple(UInt32(i) for i in 1:624);  # This is a Tuple of 1..624

julia> r = MT19937(seed);


Since MT19937 is a RNG based on linear-feedback shift-register techniques, this approach is not recommended for an obivous reason:

julia> rand(r, UInt32, 10)
10-element Array{UInt32,1}:
0x23864fee
0xdd51a593
0x20c049a2
0xde17a3df
0x20804931
0xde57a34c
0x20c049a0
0xde17a3dd
0x20804933
0xde57a34e


The firstly generated numbers are so poorly random. This is because the most bits of states are zeros. So it is better to create a MT19937 in this way:

julia> r = MT19937();


In this case, all the 624 states will be filled with truly random numbers produced by RandomDevice. If someone needs the reproducibility, just save the state r.mt and use it for next time.

An initialization function described in the original paper1 is also implemented here, so the seed can also be just one UInt32 number (or an Integer whose least 32 bits will be truncated):

julia> import Random

julia> Random.seed!(r, 0xabcdef12);

julia> rand(r, UInt32, 10)
10-element Array{UInt32,1}:
0x986c5150
0xbb9e20f5
0xfb59a25a
0x189eb49c
0xbcca5395
0x48d9fdf5
0x3193f581
0xe0b6d080
0x63154ca2
0x72b28f0d


Note that if you use one UInt32 number as seed, you will always get in a bias way.3

1. Matsumoto M, Nishimura T. Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator[J]. ACM Transactions on Modeling and Computer Simulation (TOMACS), 1998, 8(1): 3-30. doi:10.1145/272991.272995

2. L'Ecuyer P, Simard R. TestU01: AC library for empirical testing of random number generators[J]. ACM Transactions on Mathematical Software (TOMS), 2007, 33(4): 22. doi:10.1145/1268776.1268777