Pseudorandom number generators

In the 1960s IBM’s scientific library had a random number generator called RANDU: multiply the previous number by 65539 and take the remainder modulo 2³¹. It was used for simulations for many years. Below are 6000 triples of consecutive RANDU numbers, placed as points in a cube. It looks like a cloud. Press “Turn to the planes”.

All the triples lie on 15 parallel planes, with emptiness between them. Any experiment that needs three random numbers at once (a point in a volume, three steps of a walk, three parameters of a model) gets 15 sheets instead of random points. With splitmix64 the cloud looks the same from any angle.

The idea

A pseudorandom number generator is a deterministic program: a state and a step rule. From the same seed it produces the same sequence, and sooner or later the sequence cycles, since there are finitely many states. A good generator cannot be told apart from true randomness by any reasonable test, and its repetition is out of reach. A bad one gives itself away by structure: planes, periodic bits, predictability.

What a seed is and why Python runs the Mersenne Twister is told in the course chapter. Here: what is inside generators, how they are tested and how they are broken.

The linear congruential generator

The oldest and simplest generator was proposed by Derrick Lehmer in 1949: x_n+1 = (a·x_n + c) mod m. One multiplication and one remainder per number. The period is at most m, and it equals m exactly when c is coprime to m, a − 1 is divisible by all prime factors of m, and, if m is divisible by 4, a − 1 is divisible by 4 as well; Thomas Hull and Alan Dobell proved this in 1962. A long period guarantees nothing by itself: RANDU and its planes show that.

The listing checks this on 65,534 consecutive triples: all lie on the planes, and all 15 are occupied. RANDU is an extreme case but not an exception. George Marsaglia proved in 1968 that for any linear generator the t-tuples of consecutive numbers lie on at most (t!·m)^1/t hyperplanes. For m = 2³¹ and triples that is about 2344 planes for the best multipliers, and 15 for RANDU.

The spectral test

In the plane the same can be seen on pairs. Take a generator with the small modulus 4096 and mark all the pairs (x_n, x_n+1) over a full period. Move the multiplier: for some a the points lie on three or four lines, for others there are many lines and they are almost invisible.

Why there are always lines: if s₁ + a·s₂ is divisible by m, then s₁x_n + s₂x_n+1 modulo m is the same number for all n, and all the pairs divided by m lie on the lines s₁u + s₂v = an integer plus a constant offset. Neighbouring lines of such a family are 1/√(s₁² + s₂²) apart. The worst is the family with the sparsest lines, that is, the shortest vector (s₁, s₂). Its length is what the spectral test measures, proposed by Robert Coveyou and Robert MacPherson in 1967. The shortest vector is never longer than (4/3)^1/4·√m, so the ratio of its length to this bound is a score from 0 to 1.

In the plane RANDU looks good: 0.93. Its trouble is in the volume, where the shortest vector is (9, −6, 1), shorter than 11 out of a possible 1450 or so. That is why the test is run in several dimensions at once. The MINSTD of Stephen Park and Keith Miller with the multiplier 16807 is weak in the plane: a small number x goes to 16807·x, also relatively small. In 1993 the authors themselves proposed the multiplier 48271, and the score rose from 0.34 to 0.90.

The shortest vector is found quickly. All solutions of s₁ + a·s₂ ≡ 0 form a lattice generated by the vectors (m, 0) and (−a, 1). Lagrange–Gauss reduction is a two-dimensional Euclidean algorithm: subtract from the long vector the nearest multiple of the short one, and repeat until the long one stops getting shorter. There are O(log m) steps, as with Euclid. This is the second problem of the article. In three or more dimensions the same is done by lattice reduction, for example the Lenstra–Lenstra–Lovász algorithm.

The low bits

For generators with modulus 2³² or 2⁴⁸ the remainder costs almost nothing: it is just the high bits being dropped on overflow. The price is the low bits. Below is a picture of 16,384 consecutive numbers: each point is one bit of the next number.

With odd a and c the lowest bit simply alternates: 0, 1, 0, 1. The low byte repeats every 256 numbers, and in the listing’s test it is “too even”: each value exactly 256 times in 65,536 numbers, which never happens with honest random numbers. Once this already spoiled a measurement on this site: in the article on empirical complexity long numbers assembled from the low bytes of such a generator came out periodic, and Karatsuba multiplication in Go ran 19 times faster than on real data. That is why java.util.Random keeps 48 bits of state and hands out only the upper 32, and MINSTD takes the prime modulus 2³¹ − 1, which does not have this trouble.

Better generators

xorshift. In 2003 George Marsaglia proposed generators without multiplication: the state is mixed by shifts and exclusive or, x ^= x << 13, x ^= x >> 17, x ^= x << 5. This is very fast, and it passes simple tests. But every step is linear over the bits: any output bit is the exclusive or of some bits of the state. The Berlekamp–Massey algorithm finds in such a bit sequence a linear recurrence of length at most 32: xorshift32 has a linear complexity of 32 instead of the expected 500 over a thousand bits. François Panneton and Pierre L’Ecuyer showed the weaknesses of xorshift on large test batteries in 2005 as well. The cure is a nonlinear step at the output: multiply the result by a constant or add two states.

PCG. In 2014 Melissa O’Neill took a linear generator with modulus 2⁶⁴, whose high bits are good, and put a permutation at the output: the high bits give the size of a shift or rotation for the rest. The lattice stays inside but is not visible from outside. PCG is one of the generators in Go’s math/rand/v2 package; Go’s global generator since version 1.22 is ChaCha8.

splitmix64. In 2014 Guy Steele, Doug Lea and Christine Flood described a generator whose state is just a counter growing by an odd constant γ, with all the quality in the mixing function: two multiplications and three shifts with exclusive or. The k-th number equals mix(seed + k·γ); it can be obtained without all the previous ones, and independent streams for parallel tasks are easy to set up. This catalogue’s listings use it wherever Python and Go must produce the same random numbers: five lines give the same stream in Python, Go and the scenes on the page, so the listing and the scene show the same experiment.

The Mersenne Twister. The 1998 generator of Makoto Matsumoto and Takuji Nishimura, with period 2¹⁹⁹³⁷ − 1, sits in Python, R and many other systems. Its state is 624 32-bit words, about 2.5 KB. It is linear over the bits, like xorshift, and each number it outputs is a reversibly scrambled word of the state, so 624 consecutive numbers recover the state completely, as the course chapter explains.

How a generator is tested

It is impossible to prove that a generator is good: a test can only catch a bad one. The listing runs four simple tests on 65,536 numbers of each generator. The first three are chi-squared: the low byte should take all 256 values roughly equally; pairs of low four-bit groups, all 256 combinations; triples of the high four bits, all 4096 cells of a cube. The deviation is converted to z: for an honest generator |z| is almost always below 3.3; a large positive value means too uneven, a large negative one suspiciously even. The fourth test is the linear complexity of the lowest bit over a thousand numbers, about 500 for random bits.

All of RANDU’s numbers are odd (the multiplier and the seed are odd, and the modulus is a power of two), so the lowest bit is always 1, and the triples in the cube see the planes. The generator from Numerical Recipes fails on the low bits: the byte is “too even”, the pairs repeat. xorshift32 passes three tests and fails the fourth. MINSTD, java.util.Random and splitmix64 pass all four, which only means that these four tests did not catch them.

Real test batteries are larger. TestU01 by Pierre L’Ecuyer and Richard Simard, released in 2007, runs over a hundred tests on hundreds of billions of numbers in its BigCrush set. But even a battery is not the last word. In 1992 Alan Ferrenberg, David Landau and Joanna Wong simulated the two-dimensional Ising model with Wolff’s cluster algorithm using the generator R250, which was considered good, and got a systematically wrong energy: the algorithm turned out to be sensitive to subtle correlations between bits that the standard tests did not notice. Since then a proven practice is to repeat an important computation with a generator of a different design.

Prediction

A generator is deterministic: whoever knows the state knows everything it will output. java.util.Random hides the low 16 of its 48 bits, but from two consecutive numbers they can be brute-forced: 65,536 variants, one multiplication per check. As a rule exactly one fits, and from then on every number is known in advance. This is the first problem of the article.

A poor seed is more dangerous still. In 2008 it came out that after an edit to the OpenSSL package of the Debian distribution, its key generator had for almost two years taken only the process number from all its sources of randomness: about 32 thousand variants. All possible keys of each type could be enumerated in advance. That is why keys, passwords and tokens need cryptographic generators that take randomness from the operating system: secrets in Python and crypto/rand in Go; more in the course chapter on ciphers.

Code

The listing sets up six generators and runs the four tests on them, checks RANDU’s planes, computes the spectral test for five real generators, recovers the state of java.util.Random from two numbers, and shows how splitmix64 produces the hundred-thousandth number directly.

What it costs

Each number is O(1): for a linear generator one multiplication and a remainder, for xorshift three shifts and three exclusive ors, for splitmix64 an addition, two multiplications and three shifts. Memory ranges from 4 bytes for xorshift32 and 8 for splitmix64 to 2.5 KB for the Mersenne Twister; the period from 2³¹ − 2 for MINSTD to 2¹⁹⁹³⁷ − 1. The spectral test in the plane takes O(log m) reduction steps; in higher dimensions lattice reduction is more expensive, but for a generator it is done once.

Where it works

Simulation. The Monte Carlo method eats billions of numbers, and quality matters here: RANDU’s planes turn straight into wrong answers.

Algorithms. Shuffling, sampling from a stream, the random pivot of quicksort, random hash bases.

Tests. The test generators for this catalogue’s problems create inputs from a fixed seed, so every run gives the same tests and the answers can be checked.

Games and generation. A world, a level or a layout is built from a single seed: pass on the seed and you have passed on the whole world.

History

In the 1940s John von Neumann proposed the middle-square method: square the number and take the middle digits. The method quickly falls into a cycle, and von Neumann knew it. In 1949 Derrick Lehmer described the linear congruential generator, and in 1962 Hull and Dobell found the conditions for a full period. In the 1960s RANDU spread with IBM’s scientific library; in 1967 Coveyou and MacPherson proposed the spectral test, in 1968 Marsaglia proved the theorem on planes, and in 1969 Donald Knuth gathered the theory of generators in the second volume of The Art of Computer Programming. Then came Park and Miller’s MINSTD in 1988, the Mersenne Twister in 1998, xorshift in 2003, the TestU01 battery in 2007, and PCG and splitmix in 2014.

Problems

The first is to predict java.util.Random from two numbers. The second is the spectral test in the plane for thousands of generators: Lagrange–Gauss reduction.