Lehmer random number generator

The Lehmer random number generator^[1] (named after D. H. Lehmer), sometimes also referred to as the Park–Miller random number generator (after Stephen K. Park and Keith W. Miller), is a type of linear congruential generator (LCG) that operates in multiplicative group of integers modulo n. The general formula is:

where the modulus n is a prime number or a power of a prime number, the multiplier g is an element of high multiplicative order modulo n (e.g., a primitive root modulo n), and the seed X₀ is coprime to n.

Parameters in common use

In 1988, Park and Miller^[2] suggested a Lehmer RNG with particular parameters n = 2³¹ − 1 = 2,147,483,647 (a Mersenne prime M₃₁) and g = 7⁵ = 16,807 (a primitive root modulo M₃₁), now known as MINSTD. Although MINSTD was later criticized by Marsaglia and Sullivan (1993),^[3] it is still in use today (in particular, in CarbonLib and C++11's minstd_rand0). Park, Miller and Stockmeyer responded to the criticism (1993),^[4] 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." Our article and the associated minimal standard generator was an attempt to respond to this request. Five years later, we see no need to alter our response other than to suggest the use of the multiplier a = 48271 in place of 16807.

This revised constant is used in C++11's minstd_rand random number generator.

The Sinclair ZX81 and its successors use the Lehmer RNG with parameters n = 2¹⁶ + 1 = 65,537 (a Fermat prime F₄) and g = 75 (a primitive root modulo F₄). ^[5] The CRAY random number generator RANF is a Lehmer RNG with n = 2⁴⁸ − 1 and g = 44,485,709,377,909.^[6] The GNU Scientific Library includes several random number generators of the Lehmer form, including MINSTD, RANF, and the infamous IBM random number generator RANDU.^[6]

Relation to LCG

While the Lehmer RNG can be viewed as a particular case of the linear congruential generator with c=0, it is a special case that implies certain restrictions and properties. In particular, for the Lehmer RNG, the initial seed X₀ must be coprime to the modulus n that is not required for LCGs in general. The choice of the modulus n and the multiplier g is also more restrictive for the Lehmer RNG. In contrast to LCG, the maximum period of the Lehmer RNG equals n−1 and it is such when n is prime and g is a primitive root modulo n.

On the other hand, the discrete logarithms (to base g or any primitive root modulo n) of X_k in represent linear congruential sequence modulo Euler totient .

Sample C99 code

Using C code, the Lehmer RNG can be written as follows.

#define M 2147483647 /* 2^31 - 1 (A large prime number) */
#define A 16807      /* Prime root of M, passes statistical tests and produces a full cycle */
#define Q 127773 /* M / A (To avoid overflow on A * seed) */
#define R 2836   /* M % A (To avoid overflow on A * seed) */

uint32_t lcg_parkmiller(uint32_t seed)
{
    uint32_t hi = seed / Q;
    uint32_t lo = seed % Q;
    int32_t test = A * lo - R * hi;
    if (test <= 0)
        test += M;
    return test;
}

This function's output can be fed back to its input to generate pseudorandom numbers, as long as the caller is careful to begin with any number except zero. As the product of two 32 bit integers may overflow, the cast to uint64_t is necessary.

Another popular pair of Lehmer generator parameters uses the prime modulus 2³²-5:

uint32_t lcg_rand(uint32_t a)
{
    return ((uint64_t)a * 279470273UL) % 4294967291UL;
}

Many other Lehmer generators and primes have good properties. The following 128-bit Lehmer generator requires 128-bit support from the compiler and uses a multiplier computed by L'Ecuyer.^[7] It has a period of 2¹²⁶:

/* The state must be seeded with two 64-bit values, among which s[0] MUST be odd. */
static union {
	__int128 x;
	uint64_t s[2];
} state;

uint64_t next(void) {
	state.x *= ((__int128)0x12e15e35b500f16e << 64 | 0x2e714eb2b37916a5);
	return state.s[1];
}

The generator computes an odd 128-bit value and returns its upper 64 bits. Note that the role of s[0] and s[1] must be exchanged in a big-endian architecture.

This generator passes the most stringent statistical tests such as BigCrush from TestU01 and PractRand.

References

• Lehmer, D. H. (1949). "Mathematical methods in large-scale computing units". Proceedings of a Second Symposium on Large-Scale Digital Calculating Machinery: 141–146. MR 0044899.  (journal version: Annals of the Computation Laboratory of Harvard University, Vol. 26 (1951)).
• Martin Greenberger (1961). "Notes on a New Pseudo-Random Number Generator". Journal of the ACM. 8 (2): 163–167. doi:10.1145/321062.321065. 
• Steve Park, Random Number Generators
• Park–Miller–Carta Pseudo-Random Number Generator