Random Numbers Part 1: Modulus and Primes
I’ve wondered a few times how we might generate a one-to-one mapping between a range of integers, and a “shuffled” counterpart in that range. For example, imagine a website where we expose primary keys as part of the URL path:
http://www.example.com/pageId=123
This invites someone to try plugging in different pageId values into the URL, as it would be a fair guess that it just represents a database sequence. A scraper could enumerate all pages on the website by incrementing pageId by 1 on each request.
Now this isn’t the right way to solve this problem at all, but just for fun let’s imagine we wanted to solve this problem by mapping each database key to an “encrypted” version of itself, which is some other integer:
http://www.example.com/pageId=7834190
When the server receives the request, it “decrypts” 7834190 back to 123 and the correct page can then be retrieved from the database.
Is this possible? Can we come up with a way to map every integer in a certain range to exactly one other integer? There are ways to trivially achieve this of course: simply add 1 - we’ve achieved the stated aim but this clearly isn’t what we’re looking for. We’re looking for something more “random” than that - a series of incrementing integers maps to outputs which jump all over the place, for example:
1 -> 72911
2 -> 412
3 -> 9910844
...
There’s no way for a third party to figure out how to retrieve the record for a particular database key.
What we need is some kind of function which is reversible, so one input always gives a single output. This way, by definition we can get back to the input. If two inputs mapped to one output, then there would be no way to get back to the input - which one would it be?
Hash functions like MD5 and SHA-256 are known for being lossy and irreversible. Once you’ve put a value into the hash function and got the result, you can’t recover the original value just from its output hash.
What we’re looking for is something more like an encryption cipher - given a value we can encrypt it to some output ciphertext, and we can then decrypt that ciphertext back to the plain text input.
It’s a bit simpler in our case as we’re just dealing with numbers. Let’s say we’re only dealing with 32-bit numbers - we’re going to map every 32-bit number to some other 32-bit number, pseudorandomly. This is called a bjijection - there’s a set of inputs (all the possible numbers we might want to encrypt) and outputs (all the numbers we get as output). Every possible input number points to a single output number, and every number in the output set maps perfectly back to its input, and without any obvious pattern. Can we do this?
Of course one easy way would be to come up with a mapping table. Just decide which value each possible input will map to. Possible, but there are several problems: (i) how do we decide what the output numbers are, cleanly, and without duplicating any one? (ii) how much storage would it take?
The answer to question (ii) is as follows. We need to store 32 bit values for each value we’re mapping. But we’re mapping all the 32 bit values. So there are 2^32 mappings (for each of the possible values we can hold in a 32 bit value), and each mapping costs us 32 * 2 bits (2 32-bit numbers: the input and the output). This adds up to 2^32 * 64 bits, or 2^6, or 2^(32+6), or 2^38 bits which is 2^(38-3) bytes, which is 2^35 bytes, or 34GB! To be fair, we don’t need to store the input values, as this is just the index in the list. So this cuts storage in half: 17GB which is still a pretty big ask.
It’s actually much easier than this - we can use a simple calculation with almost no storage, and very few instructions.
This is how pseudo-random number generators work, and the principle hasn’t changed much from 1949 when the approach was proposed by D.H. Lehmer, and now when the approach is still used for generating random numbers in places like the Java language, and the standard C libraries, at least when we don’t need strong cryptographic guarantees. The approach is called the Linear Congruential Generator (or LCG).
An LCG uses several mathematical parts: prime numbers, and the modular arithemetic. Let’s understand how all this works step-by-step. We’ll use a simpler type of LCG - one which only performs a multiplicative update, and with no additive part. This is known as the Lehmer generator, or a Park-Miller generator or an MLCG, and is simpler than a full LCG which has an additive component, but still has been used in real systems including the 1982 ZX Spectrum, and systems with limited processing power such as microcontrollers.
Modular Arithmetic
Modular arithmetic is like clock arithmetic. If it’s now 1 ‘o’ clock then what time will it be in 14 hours time? It’ll be 3 ‘o’ clock. We work this out by adding 1 and 14 (15), and then dividing by 12 and taking the remainder. The remainder of 15/12 is 3, and this is the answer. We write this opeartion as:
\(1 + 14 \equiv 3 \pmod{12}\)
(1 + 14, or 15, is congruent to 3 modulo 12).
Random number generators are often based on modular arithmetic. How can this work? Let’s try a few ideas.
We’ll define our random number generator as a recurrence relation like this:
\(x_{n+1} = a x_n \pmod{m}\)
So the next term is equal to the previous term, multiplied by \(a\) and we then find the remainder to some modulus \(m\).
Let’s try this with our “hours of the clock” example. Our goal is to find how to ‘scramble’ the hours of the clock so we get a sequence of the hours in the wrong order.
Let’s take a seed value of 1, and a multiplier of 2:
\(x_{n+1} = 2 . x_n \pmod{12}\)
This gives the sequence:
1 (seed value)
2
4
8
4
8
4
...
Not a very good random number generator! In fact we get a cycle forever between 4 and 8. Let’s think about what we really want. We want a sequence of numbers with no particular order, which covers the entire set of possible numbers under that modulus system. So in \(\pmod{12}\) we expect to generate nearly all 12 numbers, 1 to 11, in some ‘shuffled’ order. Incidentally we need to avoid a 0 output with the multiplicative LCGs we’re playing with here, as otherwise the next term will also be 0 and that pattern will continue!
So what’s happening in our example above? Firstly, our multiplier is 2. This means every “step” we take around the clock will be a multiple of 2. So once we land on 2 (or any other even hour) we’ll never land on an odd number again. That’s not going to work.
Another problem is that cycle between 4 and 8 - what’s going on there, it’s like the clock is getting stuck. \(4 * 2 = 8\), and \(8 * 2 = 16\) and when we take the modulus 12 we get 4 again. It’s a cycle of length 2. These short cycles are particularly pathological for random number generators. How do these occur?
12 (our modulus) is a composite number; in other words it’s not a prime number, but it’s comprised from prime numbers. 12 is made up of two prime factors 2 and 3:
\(2^2 + 3 = 12\)
The factors 4 and 3 are co-prime factors of 12, which means they share no common divisor other than 1. Another way to think about this in terms of a clock face is that 12 hours are made up of 3 lots of 4 hours, or 4 lots of 3 hours.
Why are we breaking this up like this? It’s quite hard to think about what’s happening in a modulus 12 system, but it’s easier to see what’s happening when you consider there are two separate systems happening, superimposed on each other, and this will explain the cycles we were seeing earlier!
Let’s look at what happens in these different modulus systems as we run our first random number generator.
As a reminder, we’re working in \(\pmod{12}\) with a seed value of 1, and a multiplier \(a\) of 2. We’re also looking at the values of \(\pmod{3}\) and \(\pmod{4}\) at each step to see how these relate to the true value we’re interested in (\(\pmod{12}\)).
| Index | \(\pmod{3}\) | \(\pmod{4}\) | \(\pmod{12}\) |
|---|---|---|---|
| 0 (seed) | 1 | 1 | 1 |
| 1 | 2 | 2 | 2 |
| 2 | 1 | 0 | 4 |
| 3 | 2 | 0 | 8 |
| 4 | 1 | 0 | 4 |
| … |
We see our repeating pattern again in the \(\pmod{12}\) column once again, but this time accompanied by what’s happening to the values of the co-prime factors 3 and 4. Fairly quickly, the value in the \(\pmod{4}\) space becomes 0, and this is significant. What’s happened is that in the world of multiples of 4, the result has become 0. We started on 1, multiplied by 2 giving 2, and the next multiplication with 2 gives 4, which is \(0 \pmod{4}\). What does this tell us? The result in the world of \(\pmod{12}\) is now “stuck” in the world of the \(\pmod{4}\) space. If this seems unclear, instead consider each result always being \(0 \pmod{2}\). This means every result will be even as the remainder when dividing by 2 is 0. This means every possible value would be even and our random number generator is broken. In the case we’re looking at, every result is \(0 \pmod{4}\), which means every result from now on is a multiple of 4 (or 0).
And why are we bouncing between 4 and 8? The reason is the \(\pmod{3}\) space is still contributing information. Since the \(\pmod{3}\) result is always changing, the overall result is always changing. \(2 * 2 \pmod{3} \equiv 1\) and \(1 * 2 \pmod{3} \equiv 2\) and so the cycle repeats in the world of \(\pmod{3}\).
In fact, what’s happened is that our generator has reduced down to one purely over \(\pmod{3}\) and we’re just scaling up into the world of \(\pmod{12}\).
Chinese Remainder Theorem
If it’s still not obvious how the two parallel worlds of \(\pmod{3}\) and \(\pmod{4}\) play a role in the overall world of \(\pmod{12}\) then this will hopefully help.
Since we’ve broken our overall modulus system of 12 down into co-prime factors 3 and 4, we can exploit a useful property. The Chinese Remainder Theorem (from now on, CRT) tells us that for any pair of co-prime factors (in this case 3 and 4), if we know the remainder in each of these moduli, then there exists a single result in the combined product of the two moduli - \(\pmod{12}\).
Moreover, because we’re considering results in the two co-prime moduli and their product (the target modulus) we can treat the results in the table above as a linear system of equations to solve! Ignoring the seed row (which contributes no information):
\(x_0.2 + x_1.2 = 2\)
\(x_0.1 + x_1.0 = 4\)
This is not much of a system of equations as there is only one datapoint for which we get a contribution from \(\pmod{4}\). Solving it anyway, we can trivially find the co-efficient of \(\pmod{3}\) as: \(x_0 = 4\) by considering the second equation alone. We can then find the co-efficient for \(\pmod{4}\) by plugging in our result for \(x_0\) into the first equation:
\(4.2 + x_1.2 = 2\)
\(8 + x_1.2 = 2\)
\(x_1.2 = -6\)
\(x_1 = -3\)
There’s hardly any point, but testing these results for index 3:
\(4.2 + (-3).0 = 8\)
And we get the correct answer as observed in the table above.
The fact we only have one equation contributing information for the \(\pmod{4}\) component gives a geomtric perspective to all this. This is a linear system of two dimensions, and due to the lack of information on the \(\pmod{4}\) dimension, instead of the results lying on a plane they describe a one-dimensional line (pretty obvious since we’re only obtaining powers of 2 in the first place).
Fixing the short loops
Jumping abruptly to the big conclusion, we can get a much better result by choosing a prime number for the modulus. For example:
- modulus = 11
- multiplier = 3
- seed = 1
So our new generator is defined as:
\(x_{n+1} = 3 . x_n \pmod{11}\)
We’ll tabulate the results of the generator. There are no factors of 11, other than 1 and 11 itself, since 11 is prime.
| Index | \(\pmod{11}\) |
|---|---|
| 0 (seed) | 1 |
| 1 | 3 |
| 2 | 9 |
| 3 | 5 |
| 4 | 4 |
| 5 | 1 |
| … |
This is slightly better than our first attempt as we now get a period of 5 rather than 4. Not a big improvement, and in fact we don’t even get full cycle of values 1 to 10. Can we achieve this? Let’s try some other multipliers:
- modulus = 11
- multiplier = 2
- seed = 1
So our new generator is defined as:
\(x_{n+1} = 2 . x_n \pmod{11}\)
| Index | \(\pmod{11}\) |
|---|---|
| 0 (seed) | 1 |
| 1 | 2 |
| 2 | 4 |
| 3 | 8 |
| 4 | 5 |
| 5 | 10 |
| 6 | 9 |
| 7 | 7 |
| 8 | 3 |
| 9 | 6 |
| 10 | 1 |
| … |
That’s it! We designed a random number generator with cycle 10. It’s not the greatest ever, since it has such a short period and obviously starts with powers of 2. But this illustrates the idea in a simple way.
Why Primes?
Why did the prime number modulus work so well?
A big clue is that there are no factors to consider. We can think of the factors of a compound modulus (say 12) forming sub-groups which we can get stuck within. This all sounds vague, so let’s try to explain a bit better.
With our earlier example in \(\pmod{12}\) we can easily get “stuck” in a factor, or sub-group. Say if we use a multiplier of 2, we only get results which are powers of 2: 1, 2, 4, 8, 4, 8, … Here’s what happens:
| Index | \(\pmod{12}\) | What’s happening |
|---|---|---|
| 0 (seed) | 1 | This is our seed. From here we could get anywhere in \(\pmod{12}\) if we could choose our multipler. But it’s fixed at 2. |
| 1 | 2 | We now can only reach multiples of 2 - reducing future values down to even numbers only. |
| 2 | 4 | We now can only reach multiples of 4 - we can now only reach 4 or 8 in the world of \(\pmod{12}\)! |
| 3 | 8 | Sure enough - this is the only other value we could get to. |
| 4 | 4 | No surprise again - we’re stuck with only 4 or 8 being reachable. |
| … |
In \(\pmod{12}\) Other multiplers don’t do much better. We either get stuck in these “factor traps”, or something else happens. Say with a multiplier of 5:
| Index | \(\pmod{12}\) | What’s happening |
|---|---|---|
| 0 (seed) | 1 | This is our seed. From here we could get anywhere in \(\pmod{12}\) if we could choose our multipler. But it’s fixed at 5. |
| 1 | 5 | We now can only reach multiples of 5 - reducing future values down to only 1 or 5 in \(\pmod{12}\).. |
| 2 | 1 | Sure enough … |
| 3 | 5 | |
| … |
This is a bit of a special case - a modulus system has the concept of a multiplicative inverse. We’re used to the idea of a reciprocal:
\(x * 1/x = 1\)
Some x has a reciprocal which when multiplied together give 1. In a modulus system the same idea applies:
\(x * r \equiv 1 \pmod{M}\)
In the table above we see the output oscillating between 1 and 5. This is because 5 is its own multiplicative inverse in \(\pmod{12}\):
\(5 * 5 \equiv 1 \pmod{12}\)
When we see a 1 in the output of a purely multiplicative LCG its a sign we’re now stuck in a loop, sometimes for this reason and sometimes just because we’ve reached the end of the cycle!
Do primes always work?
With a prime modulus, we don’t have to worry about the factor traps described above, nor can we hit the multiplicative inverse case. But we still have to be careful which multiplier we pick. For example in \(\pmod{11}\) if we pick \(a = 3\) we only get a cycle of 5, unlike when using multiplier 2 where we saw a cycle of 10.
Even worse, if we pick \(a = 10\) in \(\pmod{11}\) we get a cycle of length 2. We start at 1, then we get 10, and the next result is 1 and the cycle repeats. This is always the case when the multiplier is one less than the modulus, prime or not.
So, good parameters have been provided for use in purely multiplicative LCGs. For example, in 1982 the ZX Spectrum came out with its random number generator based on the following MLCG:
\(x_{n+1} = 75 . x_n \pmod{65537}\)
The choice of modulus is \(2^{16} + 1\) - a Fermat Prime, which is a particularly good choice for a computer which naturally works in binary (base-2)!
Ignoring the implementation details and purely looking at the results of the Lehmer generator used in the ZX Spectrum we get:
1, 75, 5,625, 28,653, 51,791, 17,642, 11,013, ...
Another good choice for the modulus is \(2^{31} - 1\), a Mersenne Prime, and also a good choice for a computational target. This modulus was proposed in a 1988 paper by Park and Miller (Communications of the ACM article “Random Number Generators: Good Ones Are Hard To Find”), leading to the name Park-Miller Generator for this type of generator.