Lehmer random number generator

In 1988, Park and Miller[3] suggested a Lehmer RNG with particular parameters m = 231 − 1 = 2,147,483,647 (a Mersenne prime M31) and a = 75 = 16,807 (a primitive root modulo M31), now known as MINSTD.

Park, Miller and Stockmeyer responded to the criticism (1993),[6] saying: Given the dynamic nature of the area, it is difficult for nonspecialists to make decisions about what generator to use.

"Give me something I can understand, implement and port... it needn't be state-of-the-art, just make sure it's reasonably good and efficient."

The Sinclair ZX81 and its successors use the Lehmer RNG with parameters m = 216 + 1 = 65,537 (a Fermat prime F4) and a = 75 (a primitive root modulo F4).

[7][8] The CRAY random number generator RANF is a Lehmer RNG with the power-of-two modulus m = 248 and a = 44,485,709,377,909.

[9] Most commonly, the modulus is chosen as a prime number, making the choice of a coprime seed trivial (any 0 < X0 < m will do).

Therefore, the application using these random numbers must use the most significant bits; reducing to a smaller range using a modulo operation with an even modulus will produce disastrous results.

Using a composite modulus is possible, but the generator must be seeded with a value coprime to m, or the period will be greatly reduced.

In particular, for the Lehmer RNG, the initial seed X0 must be coprime to the modulus m, which is not required for LCGs in general.

A prime modulus requires the computation of a double-width product and an explicit reduction step.

If a modulus just less than a power of 2 is used (the Mersenne primes 231 − 1 and 261 − 1 are popular, as are 232 − 5 and 264 − 59), reduction modulo m = 2e − d can be implemented more cheaply than a general double-width division using the identity 2e ≡ d (mod m).

This avoids the need to consider equivalent e-bit representations of the state; only values where the high bits are non-zero need reduction.

Even in high-level languages, if the multiplier a is limited to √m, then the double-width product ax can be computed using two single-width multiplications, and reduced using the techniques described above.

However, each step also divides x by an ever-increasing quotient q = ⌊m/a⌋, and quickly a point is reached where the argument is 0 and the recursion may be terminated.

That test has been designed to catch exactly the defect of this type of generator: since the modulus is a power of 2, the period of the lowest bit in the output is only 262, rather than 2126.

The following core routine improves upon the speed of the above code for integer workloads (if the constant declaration is allowed to be optimized out of a calculation loop by the compiler): However, because the multiplication is deferred, it is not suitable for hashing, since the first call simply returns the upper 64 bits of the seed state.