A linear congruential generator (LCG) is an algorithm that yields a sequence of pseudo-randomized numbers calculated with a discontinuous piecewise linear equation. The method represents one of the oldest and best-known pseudorandom number generator algorithms. The theory behind them is relatively easy to understand, and they are easily implemented and fast, especially on computer hardware which can provide modular arithmetic by storage-bit truncation.
The generator is defined by the recurrence relation:
where is the sequence of pseudo-random values, and
— the "modulus" — the "multiplier" — the "increment" — the "seed" or "start value"
are integer constants that specify the generator. If c = 0, the generator is often called a multiplicative congruential generator (MCG), or Lehmer RNG. If c ≠ 0, the method is called a mixed congruential generator.
When c ≠ 0, a mathematician would call the recurrence an affine transformation, not a linear one, but the misnomer is well-established in computer science. Associated software and data at https://github.com/vigna/CPRNG.
History The Lehmer generator was published in 1951 and the Linear congruential generator was published in 1958 by W. E. Thomson and A. Rotenberg.
Period length A benefit of LCGs is that an appropriate choice of parameters results in a period which is both known and long. Although not the only criterion, too short a period is a fatal flaw in a pseudorandom number generator.
While LCGs are capable of producing pseudorandom numbers which can pass formal tests for randomness, the quality of the output is extremely sensitive to the choice of the parameters m and a. For example, a = 1 and c = 1 produces a simple modulo-m counter, which has a long period, but is obviously non-random. Other values of c coprime to m produce a Weyl sequence, which is better distributed but still obviously non-random.
Historically, poor choices for a have led to ineffective implementations of LCGs. A particularly illustrative example of this is RANDU, which was widely used in the early 1970s and led to many results which are currently being questioned because of the use of this poor LCG.
There are three common families of parameter choice:
m prime, c = 0
This is the original Lehmer RNG construction. The period is m−1 if the multiplier a is chosen to be a primitive element of the integers modulo m. The initial state must be chosen between 1 and m−1.
One disadvantage of a prime modulus is that the modular reduction requires a double-width product and an explicit reduction step. Often a prime just less than a power of 2 is used (the Mersenne primes 231−1 and 261−1 are popular), so that the reduction modulo m = 2e − d can be computed as (ax mod 2e) + d . This must be followed by a conditional subtraction of m if the result is too large, but the number of subtractions is limited to ad/m, which can be easily limited to one if d is small.
If a double-width product is unavailable, and the multiplier is chosen carefully, Schrage's method may be used. To do this, factor m = qa+r, i.e. and . Then compute ax mod m = . Since x mod q < q ≤ m/a, the first term is strictly less than am/a = m. If a is chosen so that r ≤ q (and thus r/q ≤ 1), then the second term is also less than m: r ≤ rx/q = x(r/q) ≤ x < m. Thus, both products can be computed with a single-width product, and the difference between them lies in the range [1−m, m−1], so can be reduced to [0, m−1] with a single conditional add.
The most expensive operation in Schrage's method is the division (with remainder) of x by q; fast algorithms for division by a constant are not available since they also rely on double-width products.
A second disadvantage of a prime modulus is that it is awkward to convert the value 1 ≤ x < m to uniform random bits. If a prime just less than a power of 2 is used, sometimes the missing values are simply ignored.
m a power of 2, c = 0 Choosing m to be a power of two, most often m = 232 or m = 264, produces a particularly efficient LCG, because this allows the modulus operation to be computed by simply truncating the binary representation. In fact, the most significant bits are usually not computed at all. There are, however, disadvantages.
This form has maximal period m/4, achieved if a ≡ ±3 (mod 8) and the initial state X0 is odd. Even in this best case, the low three bits of X alternate between two values and thus only contribute one bit to the state. X is always odd (the lowest-order bit never changes), and only one of the next two bits ever changes. If a ≡ +3, X alternates ±1↔±3, while if a ≡ −3, X alternates ±1↔∓3 (all modulo 8).
It can be shown that this form is equivalent to a generator with modulus m/4 and c ≠ 0.
A more serious issue with the use of a power-of-two modulus is that the low bits have a shorter period than the high bits. Its simplicity of implementation comes from the fact that bits are never affected by higher-order bits, so the low b bits of such a generator form a modulo-2b LCG by themselves, repeating with a period of 2b−2. Only the most significant bit of X achieves the full period.
m a power of 2, c ≠ 0 When c ≠ 0, correctly chosen parameters allow a period equal to m, for all seed values. This will occur if and only if: and are coprime, is divisible by all prime factors of , is divisible by 4 if is divisible by 4.
These three requirements are referred to as the Hull–Dobell Theorem.
This form may be used with any m, but only works well for m with many repeated prime factors, such as a power of 2; using a computer's word size is the most common choice. If m were a square-free integer, this would only allow a ≡ 1 (mod m), which makes a very poor PRNG; a selection of possible full-period multipliers is only available when m has repeated prime factors.
Although the Hull–Dobell theorem provides maximum period, it is not sufficient to guarantee a good generator. For example, it is desirable for a − 1 to not be any more divisible by prime factors of m than necessary. If m is a power of 2, then a − 1 should be divisible by 4 but not divisible by 8, i.e. a ≡ 5 (mod 8).
Indeed, most multipliers produce a sequence which fails one test for non-randomness or another, and finding a multiplier which is satisfactory to all applicable criteria is quite challenging. The spectral test is one of the most important tests.
Note that a power-of-2 modulus shares the problem as described above for c = 0: the low k bits form a generator with modulus 2k and thus repeat with a period of 2k; only the most significant bit achieves the full period. If a pseudorandom number less than r is desired, is a much higher-quality result than X mod r. Unfortunately, most programming languages make the latter much easier to write (X % r), so it is very commonly used.
The generator is not sensitive to the choice of c, as long as it is relatively prime to the modulus (e.g. if m is a power of 2, then c must be odd), so the value c=1 is commonly chosen.
The sequence produced by other choices of c can be written as a simple function of the sequence when c=1. Specifically, if Y is the prototypical sequence defined by Y0 = 0 and Yn+1 = aYn + 1 mod m, then a general sequence Xn+1 = aXn + c mod m can be written as an affine function of Y:
More generally, any two sequences X and Z with the same multiplier and modulus are related by
In the common case where m is a power of 2 and a ≡ 5 (mod 8) (a desirable property for other reasons), it is always possible to find an initial value X0 so that the denominator X1 − X0 ≡ ±1 (mod m), producing an even simpler relationship. With this choice of X0, Xn = X0 ± Yn will remain true for all n. The sign is determined by c ≡ ±1 (mod 4), and the constant X0 is determined by 1 ∓ c ≡ (1 − a)X0 (mod m).
As a simple example, consider the generators Xn+1 = 157Xn + 3 mod 256 and Yn+1 = 157Yn + 1 mod 256; i.e. m = 256, a = 157, and c = 3. Because 3 ≡ −1 (mod 4), we are searching for a solution to 1 + 3 ≡ (1 − 157)X0 (mod 256). This is satisfied by X0 ≡ 41 (mod 64), so if we start with that, then Xn ≡ X0 − Yn (mod 256) for all n.
For example, using X0 = 233 = 3×64 + 41: X = 233, 232, 75, 2, 61, 108, ... Y = 0, 1, 158, 231, 172, 125, ... X + Y mod 256 = 233, 233, 233, 233, 233, 233, ...
Parameters in common use The following table lists the parameters of LCGs in common use, including built-in rand() functions in runtime libraries of various compilers. This table is to show popularity, not examples to emulate; many of these parameters are poor. Tables of good parameters are available.
Source modulusm multipliera incrementc output bits of seed in rand() or Random(L) ZX81, ZX Spectrum 216 + 1 75 0 (xn − 1) / 216 Numerical Recipes ranqd1, Chapter 7.1, §An Even Quicker Generator, Eq. 7.1.6 parameters from Knuth and H. W. Lewis 232 1664525 1013904223 Borland C/C++ 231 22695477 1 bits 30..16 in rand(), 30..0 in lrand() glibc (used by GCC)Implementation in glibc-2.26 release. See the code after the test for "TYPE_0"; the GNU C library's rand() in stdlib.h uses a simple (single state) linear congruential generator only in case that the state is declared as 8 bytes. If the state is larger (an array), the generator becomes an additive feedback generator (initialized using minstd_rand0) and the period increases. See the simplified code that reproduces the random sequence from this library. 231 1103515245 12345 bits 30..0 ANSI C: Watcom, Digital Mars, CodeWarrior, IBM VisualAge C/C++C90, C99, C11: Suggestion in the ISO/IEC 9899, C17 231 1103515245 12345 bits 30..16 Borland Delphi, Virtual Pascal 232 134775813 1 bits 63..32 of (seed × L) Turbo Pascal 4.0 et seq. 232 134775813 (808840516) 1 Microsoft Visual/Quick C/C++ 231 214013 (343FD16) 2531011 (269EC316) bits 30..16 Microsoft Visual Basic (6 and earlier) 224 1140671485 (43FD43FD16) 12820163 (C39EC316) RtlUniform from Native APIIn spite of documentation on MSDN, RtlUniform uses LCG, and not Lehmer's algorithm, implementations before Windows Vista are flawed, because the result of multiplication is cut to 32 bits, before modulo is applied 231 − 1 −18 (7FFFFFED16) −60 (7FFFFFC316) Apple CarbonLib, C++11's minstd_rand0, MATLAB's v4 legacy generator mcg16807 231 − 1 16807 0 see MINSTD C++11's minstd_rand 231 − 1 48271 0 see MINSTD rsixfour.c, from the MMIX Supplement 264 6364136223846793005 9754186451795953191 Newlib 263 6364136223846793005 1 bits 62..32 (46..32 for 16-bit int) Musl 264 6364136223846793005 1 bits 63..33 VMS's MTH$RANDOM, old versions of glibc 232 69069 (10DCD16) 1 Java's java.utilRandom, POSIX [ln]rand48, glibc [ln]rand48[_r] 248 25214903917 (5DEECE66D16) 11 bits 47..17 POSIXThe Open Group Base Specifications Issue 7 IEEE Std 1003.1, 2013 Edition [dejm]rand48, glibc [dejm]rand48[_r] 248 25214903917 (5DEECE66D16) 11 bits 47..0 or bits 47..16, as required random0Stephen J. Chapman. "Example 6.4 – Random Number Generator". "MATLAB Programming for Engineers". 2015. pp. 253–256.Stephen J. Chapman. "Example 6.4 – Random Number Generator". "MATLAB Programming with Applications for Engineers". 2012. pp. 292–295.S. J. Chapman. random0. 2004.Stephen J. Chapman. "Introduction to Fortran 90/95". 1998. pp. 322–324.Wu-ting Tsai. "'Module': A Major Feature of the Modern Fortran" . pp. 6–7. 134456 = 2375 8121 28411 cc65's rand 232 16843009 (101010116) 826366247 (3141592716) [orig] bits 22..8 {effective modulus 224} bits 22..16,31..24 {change #323} cc65's rand 232 16843009 (101010116) 3014898611 (B3B3B3B316) bits 22..16,31..24 {change #951} bits 22..16,31..24 xor bits 6..0,15..8 {change #1107} Formerly common: 231 65539 0