The Mersenne Twister

I stumbled upon the Mersenne Twister after investigating Pseudorandom Number Generators (PRNGs) stimulated by an article about a hacker who had made a speciality of exploiting weaknesses in the PRNGs incorporated into certain types of slot machines.

Here is an excerpt from the article:

In the course of reverse engineering Novomatic’s software, Alex encountered his first PRNG. He was instantly fascinated by the elegance of this sort of algorithm, which is designed to spew forth an endless series of results that appear impossible to forecast. It does this by taking an initial number, known as a seed, and then mashing it together with various hidden and shifting inputs—the time from a machine’s internal clock, for example. Writing such algorithms requires tremendous mathematical skill, since they’re supposed to produce an output that defies human comprehension; ideally, a PRNG should approximate the utter unpredictability of radioactive decay.

The intriguingly named Mersenne Twister is described thus by Wikipedia:

The Mersenne Twister is a pseudorandom number generator (PRNG). It is by far the most widely used general-purpose PRNG. Its name derives from the fact that its period length is chosen to be a Mersenne prime. 

The Mersenne Twister was developed in 1997 by Makoto Matsumoto and Takuji Nishimura. It was designed specifically to rectify most of the flaws found in older PRNGs. 

The most commonly used version of the Mersenne Twister algorithm is based on the Mersenne prime \(2^{19937}−1\). The standard implementation of that, MT19937, uses a 32-bit word length. 

The Mersenne Twister is the default PRNG for the following software systems: Microsoft Excel, GAUSS,  GLib, GNU Multiple Precision Arithmetic Library, GNU Octave, GNU Scientific Library, gretl, IDL, Julia,  CMU Common Lisp, Embeddable Common Lisp, Steel Bank Common Lisp, Maple, MATLAB, Free Pascal, PHP, Python, R, Ruby, SageMath, Scilab, Stata.

I don't want to delve too deeply into the how the Mersenne Twister actually generates its random numbers. However, the following article provides a simple overview of how it works.

How does the Mersenne Twister work?

Posted February 2016

Someone asked that question on reddit, and so I replied with a high level answer that should provide a clear enough view of the algorithm:

From a high level, here's what a PRNG is supposed to look like (DIAGRAM 1):

DIAGRAM 1

You start with a seed (if you re-use the same seed you will obtain the same random numbers), you initialise it into a state. Then, every time you want to obtain a random number, you transform that state with a one-way function \(g\). This is because you don't want people to find out the state out of the random output.

You want another random number? You first transform the state with a one way function
\(f\): this is because you don't want people who found out the state to be able to retrieve past states (forward secrecy). And then you use your function \(g\) again to output a random number.

Mersenne Twister (MT) is like that, except:

• your first state is not used to output any random numbers
• a state allows you to output not only one, but 624 random numbers (although this could be thought as one big random number)
• the \(g\) function is reversible, it's not a one-way function, so MT it is not a cryptographically secure PRNG.

With more details, here's what MT looks like (DIAGRAM 2):

DIAGRAM 2

The \(f\) function is called "twist", the \(g\) function is called "temper". You can find out how each functions work by looking at the working code on the wikipedia page of MT.

*********************************

The slot machine hacking via knowledge of the PRNG within got me thinking about my own gambling on Lotto, a vice I indulge in weekly. How are the lotto numbers selected? Are there physical, numbered balls involved or is a PRNG used? There's no shortage of sites that will generate random numbers for you to enter in any lottery of your choosing but how are the actual winning numbers generated? Well it seems that the numbers are physically selected in Australian lotteries as shown in the following video:

However, Lucky Lotteries is different.

To quote from the website:

Lucky Lotteries is a raffle style jackpot game that guarantees over 10,000 prizes in every draw! Unlike other lottery games, each number is unique so there is no sharing of prizes. 

There are two Lucky Lotteries games – Super Jackpot and Mega Jackpot. 

Each game has a set amount of numbers per draw and once all available numbers are sold, the winning numbers are drawn using a random number generator.

So clearly a PRNG is used in certain lotteries. Of course, I use SageMath extensively and so it's of interest to know about what PRNG it uses. SageMath documentation explains that:

It is possible to select which random number generator drives the sampling as well as the seed. The default is the Mersenne twister. Also available are the RANDLXS algorithm and the Tausworthe generator (see the gsl reference manual for more details). These are all supposed to be simulation quality generators. For RANDLXS use rng = 'luxury' and for tausworth use rng = 'taus': 

sage: T = RealDistribution('gaussian', 1, rng = 'luxury', seed = 10)

Lotto Simulations

I've written about Lotto before in this post (Oz Lotto) and another (Losing at Lotto). The current post is not about probabilities but simulation. How can we use SageMath to simulate a lottery draw? Well as it turns out it's not too difficult.

The key is to use the ZZ.random_element(lowerbound, upperbound+1) function as in the following code. The trick is not to allow duplicate numbers and that is the reason for the inclusion of Set. This function takes our lotto list which might look like this the list [20, 37, 6, 37, 1, 19, 16] and turns it, after the sort, into the ordered set {1, 6, 16, 19, 20, 37}.

Here is the code I developed. I'm really just a tyro at both Python and SageMath so there may well be much better ways to achieve my goals but this is what I came up with:

# Generate a set of numbers for lotto games
# let N represent how many numbers are required in the game
# let G represent the number of games
# let U represent the upper bound of numbers
G=8
N=7
U=40
for i in range(G):
    lotto=[]
    while len(Set(lotto))<N:
        lotto.append(ZZ.random_element(1,U+1))
    print sorted(Set(lotto))

The following is an example of the output:

[11, 23, 27, 30, 37, 38, 40]
[4, 7, 8, 15, 20, 29, 38]
[8, 9, 13, 14, 21, 28, 36]
[7, 11, 15, 21, 28, 33, 40]
[1, 15, 17, 18, 28, 33, 39]
[4, 7, 11, 14, 18, 29, 32]
[1, 10, 12, 16, 21, 31, 33]
[2, 7, 8, 12, 13, 19, 27]

How really random this output is I don't know. I'm sure SageMath has other ways to generate more truly random numbers but this is sufficient for the modest purposes of the exercise.

The code has the added advantage of being able to accommodate "fixed" numbers. For example, let's say we always want the numbers 23 and 32 included in the output. All that needs to be done is to start with lotto = [23, 32] instead of it being an empty list.

For example:

# Generate a set of numbers for lotto games
# let N represent how many numbers are required in the game
# let G represent the number of games
# let U represent the upper bound of numbers
G=8
N=7
U=40
for i in range(G):
    lotto=[23, 32]
    while len(Set(lotto))<N:
        lotto.append(ZZ.random_element(1,U+1))
    print sorted(Set(lotto))

Here an example of the output:

[10, 11, 21, 23, 32, 35, 37]
[5, 6, 17, 23, 31, 32, 39]
[6, 11, 17, 23, 31, 32, 39]
[3, 15, 18, 23, 26, 30, 32]
[2, 12, 13, 23, 27, 28, 32]
[7, 18, 23, 26, 32, 35, 38]
[2, 11, 17, 19, 23, 32, 33]
[3, 16, 23, 26, 27, 31, 32] 

Not that I expect it will make any difference to my chances of winning but I'll test this out in some coming Lotto games. Oddly enough what got me thinking about generating some random lotto numbers was this article titled JERRY AND MARGE GO LARGE in the Huffington Post. It's about an American couple who exploited a feature in their state's lottery to make millions of dollars. It's quite an interesting story.

As I've been doing I'll include a SageMathCell with the code so that anyone curious can experiment: