# 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:

${\displaystyle X_{k+1}=a\cdot X_{k}~~{\bmod {~}}~m}$

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

## Parameters in common use

In 1988, Park and Miller[2] 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. Although MINSTD was later criticized by Marsaglia and Sullivan (1993),[3][4] 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),[5] 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 m = 216 + 1 = 65,537 (a Fermat prime F4) and a = 75 (a primitive root modulo F4). [6] The CRAY random number generator RANF is a Lehmer RNG with m = 248 − 1 and a = 44,485,709,377,909.[7] The GNU Scientific Library includes several random number generators of the Lehmer form, including MINSTD, RANF, and the infamous IBM random number generator RANDU.[7]

## Power of two modulus

Using a modulus m which is a power of two makes for a particularly convenient computer implementation, but comes at a cost: the period is at most m/4, and the low bits have periods shorter than that. This is because the low k bits form a modulo-2k generator all by themselves; the higher-order bits never affect lower-order bits. To achieve this period, the multiplier must satisfy a ≡ ±3 (mod 8), and X0 must be odd. The Xi are always odd (bit 0 never changes), they alternate between two values mod 8 (bits 2 and 1 alternate), bit 3 repeats with a period of 4, but 4 has a period of 8, and so on.

## 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 X0 must be coprime to the modulus m that is not required for LCGs in general. The choice of the modulus m and the multiplier a is also more restrictive for the Lehmer RNG. In contrast to LCG, the maximum period of the Lehmer RNG equals m−1 and it is such when m is prime and a is a primitive root modulo m.

On the other hand, the discrete logarithms (to base a or any primitive root modulo m) of Xk in ${\displaystyle \mathbb {Z} _{n}}$ represent linear congruential sequence modulo Euler totient ${\displaystyle \varphi (n)}$.

## Sample C99 code

Using C code, the Park-Miller RNG can be written as follows:

uint32_t lcg_parkmiller(uint32_t *state)
{
return *state = ((uint64_t)*state * 48271u) % 0x7fffffff;
}


This function can be called repeatedly to generate pseudorandom numbers, as long as the caller is careful to initialize the state to any number greater than zero and less than the modulus. In this implementation, 64-bit arithmetic is required; otherwise, the product of two 32-bit integers may overflow.

Though less straightforward, it is also possible to implement the Park-Miller RNG using only 32-bit arithmetic. Alternative equations must be derived which cannot result in overflow. An example C implementation is shown below:

uint32_t lcg_parkmiller(uint32_t *state)
{
const uint32_t N = 0x7fffffff;
const uint32_t G = 48271u;

/*
Indirectly compute state*G%N.

Let:
div = state/(N/G)
rem = state%(N/G)

Then:
rem + div*(N/G) == state
rem*G + div*(N/G)*G == state*G

Now:
div*(N/G)*G == div*(N - N%G) === -div*(N%G)  (mod N)

Therefore:
rem*G - div*(N%G) === state*G  (mod N)

Add N if necessary so that the result is between 1 and N-1.
*/
uint32_t div = *state / (N / G);  /* max : 2,147,483,646 / 44,488 = 48,271 */
uint32_t rem = *state % (N / G);  /* max : 2,147,483,646 % 44,488 = 44,487 */

uint32_t a = rem * G;        /* max : 44,487 * 48,271 = 2,147,431,977 */
uint32_t b = div * (N % G);  /* max : 48,271 * 3,399 = 164,073,129 */

return *state = (a > b) ? (a - b) : (a + (N - b));
}


Another popular pair of Lehmer generator parameters uses the prime modulus 232−5:

uint32_t lcg_rand(uint32_t *state)
{
return *state = ((uint64_t)*state * 279470273u) % 0xfffffffb;
}


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.[8] It has a period of 2126:

const uint64_t mult_hi = 0x12e15e35b500f16e; // Removing 'const' declarations from 'mult_hi' and 'mult_lo' may
const uint64_t mult_lo = 0x2e714eb2b37916a5; // increase speed under GCC on some modern CPU architectures.

static union {
__int128 x;
uint64_t s[2];
} state;

/* The state must be seeded with two 64-bit values, among which s[0] MUST be odd. */
void seed(__int128 seed) {
state.x = seed;
state.s[0] |= 1;
}

uint64_t next(void) {
state.x *= ((__int128)mult_hi << 64 | mult_lo);
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 BigCrush from TestU01, but fails the TMFn test from PractRand. The 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 262, rather than 2126. Linear congruential generators have a similar behavior.

## References

1. ^ W.H. Payne; J.R. Rabung; T.P. Bogyo (1969). "Coding the Lehmer pseudo-random number generator" (PDF). Communications of the ACM. 12 (2): 85–86. doi:10.1145/362848.362860.
2. ^ Stephen K. Park; Keith W. Miller (1988). "Random Number Generators: Good Ones Are Hard To Find" (PDF). Communications of the ACM. 31 (10): 1192–1201. doi:10.1145/63039.63042.
3. ^
4. ^ Sullivan, Stephen (1993). "Technical correspondence: Another test for randomness" (PDF). Communications of the ACM. 36 (7): 108. doi:10.1145/159544.376068.
5. ^ Stephen K. Park; Keith W. Miller; Paul K. Stockmeyer (1988). "Technical Correspondence: Response" (PDF). Communications of the ACM. 36 (7): 108–110. doi:10.1145/159544.376068.
6. ^ Vickers, Steve. "Chapter 5 - Functions". ZX81 Basic Programming. Sinclair Research Ltd. The ZX81 uses p=65537 & a=75 [...] (Note the ZX81 manual incorrectly states that 65537 is a Mersenne prime that equals 216-1. The ZX Spectrum manual fixed that and correctly states that it is a Fermat prime that equals 216+1.)
7. ^ a b GNU Scientific Library: Other random number generators
8. ^ Pierre L’Ecuyer (January 1999). "Tables of linear congruential generators of different sizes and good lattice structure" (PDF). Mathematics of Computation. 68 (225): 249–260. doi:10.1090/s0025-5718-99-00996-5.