View source for Lehmer random number generator

The '''Lehmer random number generator'''<ref>{{cite journal |author1=W. H. Payne |author2=J. R. Rabung |author3=T. P. Bogyo |title=Coding the Lehmer pseudo-random number generator |journal=[[Communications of the ACM]] |year=1969 |volume=12 |issue=2 |pages=85–86 |url=http://www.firstpr.com.au/dsp/rand31/p85-payne.pdf |doi=10.1145/362848.362860|s2cid=2749316 }}[https://spinthewheel.io]</ref> (named after [[D. H. Lehmer]]), sometimes also referred to as the '''Park–Miller random number generator''' (after Stephen K.&nbsp;Park and Keith W.&nbsp;Miller), is a type of [[linear congruential generator]] (LCG) that operates in [[multiplicative group of integers modulo n]]. The general formula is

: <math>X_{k+1} = a \cdot X_k \bmod m,</math>

where the modulus ''m'' is a [[prime number]] or a [[Prime power|power of a prime number]], the multiplier ''a'' is an element of high [[multiplicative order]] modulo ''m'' (e.g., a [[primitive root modulo n|primitive root modulo ''n'']]), and the seed ''X''{{sub|0}} is [[coprime]] to ''m''.

Other names are '''multiplicative linear congruential generator (MLCG)'''<ref name=cacm88>{{cite journal |title=Efficient and Portable Combined Random Number Generators |first=Pierre |last=L'Ecuyer |journal=[[Communications of the ACM]] |date=June 1988 |volume=31 |issue=6 |pages=742–774 |doi=10.1145/62959.62969 |s2cid=9593394 |url=http://www.iro.umontreal.ca/~lecuyer/myftp/papers/cacm88.pdf}}</ref> and '''multiplicative congruential generator (MCG)'''.

== Parameters in common use ==
In 1988, Park and Miller<ref>{{cite journal |first1=Stephen K. |last1=Park |first2=Keith W. |last2=Miller |title=Random Number Generators: Good Ones Are Hard To Find |journal=[[Communications of the ACM]] |year=1988 |volume=31 |issue=10 |pages=1192–1201 |url=http://www.firstpr.com.au/dsp/rand31/p1192-park.pdf |doi=10.1145/63039.63042|s2cid=207575300 }}</ref> suggested a Lehmer RNG with particular parameters ''m'' = 2{{sup|31}} − 1 = 2,147,483,647 (a [[Mersenne prime]] ''M''{{sub|31}}) and ''a'' = 7{{sup|5}} = 16,807 (a primitive root modulo ''M''{{sub|31}}), now known as '''MINSTD'''. Although MINSTD was later criticized by Marsaglia and Sullivan (1993),<ref>{{cite journal |title=Technical correspondence: Remarks on Choosing and Implementing Random Number Generators |journal=[[Communications of the ACM]] |year=1993 |volume=36 |issue=7 |pages=105–108 |url=http://www.firstpr.com.au/dsp/rand31/p105-crawford.pdf |doi=10.1145/159544.376068 |last1=Marsaglia |first1=George |s2cid=26156905 |author-link=George Marsaglia}}</ref><ref>{{cite journal |title=Technical correspondence: Another test for randomness |journal=[[Communications of the ACM]] |year=1993 |volume=36 |issue=7 |page=108 |url=http://www.firstpr.com.au/dsp/rand31/p105-crawford.pdf#page=4 |doi=10.1145/159544.376068 |author=Sullivan, Stephen|s2cid=26156905 }}</ref> it is still in use today (in particular, in [[CarbonLib]] and [[C++11]]'s <code>minstd_rand0</code>). Park, Miller and Stockmeyer responded to the criticism (1993),<ref>{{cite journal |first1=Stephen K. |last1=Park |first2=Keith W. |last2=Miller |first3=Paul K. |last3=Stockmeyer |title=Technical Correspondence: Response |journal=[[Communications of the ACM]] |year=1988 |volume=36 |issue=7 |pages=108–110 |url=http://www.firstpr.com.au/dsp/rand31/p105-crawford.pdf#page=4 |doi=10.1145/159544.376068|s2cid=26156905 }}</ref> saying:
<blockquote>
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''&nbsp;=&nbsp;48271 in place of 16807.
</blockquote>
This revised constant is used in  [[C++11]]'s <code>minstd_rand</code> random number generator.

The [[ZX81|Sinclair ZX81]] and its successors use the Lehmer RNG with parameters ''m''&nbsp;=&nbsp;2{{sup|16}}&nbsp;+&nbsp;1 = 65,537 (a [[Fermat prime]] ''F''{{sub|4}}) and ''a''&nbsp;=&nbsp;75 (a primitive root modulo ''F''{{sub|4}}).<ref>{{cite book
 |last=Vickers |first=Steve |author-link=Steve Vickers (computer scientist)
 |title=ZX81 Basic Programming
 |url=https://www.worldofspectrum.net/ZX81BasicProgramming/
 |chapter=Chapter 5. Functions
 |chapter-url=https://www.worldofspectrum.net/ZX81BasicProgramming/chap05.html
 |publisher=Sinclair Research Ltd.
 |year=1981 |edition=2nd
 |quote=The ZX81 uses p=65537 &amp; a=75 [...]
 |access-date=2024-04-21
}}
(Note that the ZX81 manual incorrectly states that 65537 is a Mersenne prime that equals 2<sup>16</sup>&nbsp;−&nbsp;1. The ZX Spectrum manual fixed that and correctly states that it is a Fermat prime that equals 2<sup>16</sup>&nbsp;+&nbsp;1.)
</ref><ref>{{cite book
 |last=Vickers |first=Steve |author-link=Steve Vickers (computer scientist)
 |title=Sinclair ZX Spectrum Basic Programming
 |url=https://worldofspectrum.org/ZXBasicManual/
 |chapter=Chapter 11. Random numbers
 |pages=73–75
 |chapter-url=https://worldofspectrum.org/ZXBasicManual/zxmanchap11.html
 |publisher=Sinclair Research Ltd.
 |year=1983 |edition=2nd
 |quote=The ZX Spectrum uses p=65537 and a=75, and stores some bi-1 in memory.
 |access-date=2022-05-26
}}</ref>
The [[CRAY]] random number generator '''RANF''' is a Lehmer RNG with the power-of-two modulus ''m''&nbsp;=&nbsp;2{{sup|48}} and ''a''&nbsp;=&nbsp;44,485,709,377,909.<ref name="gnu">[https://www.gnu.org/software/gsl/doc/html/rng.html#other-random-number-generators GNU Scientific Library: Other random number generators].</ref> The [[GNU Scientific Library]] includes several random number generators of the Lehmer form, including MINSTD, RANF, and the infamous [[IBM]] random number generator [[RANDU]].<ref name="gnu" />

== Choice of modulus ==

Most commonly, the modulus is chosen as a prime number, making the choice of a coprime seed trivial (any 0 &lt; ''X''{{sub|0}} &lt; ''m'' will do). This produces the best-quality output, but introduces some implementation complexity, and the range of the output is unlikely to match the desired application; converting to the desired range requires an additional multiplication.

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 lower bits have periods shorter than that. This is because the lowest ''k'' bits form a modulo-2{{sup|''k''}} generator all by themselves; the higher-order bits never affect lower-order bits.<ref>{{TAOCP|volume=2|edition=2|pages=12–14}}</ref> The values ''X{{sub|i}}'' are always odd (bit 0 never changes), bits 2 and 1 alternate (the lower 3 bits repeat with a period of 2), the lower 4 bits repeat with a period of 4, and so on. 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.<ref>{{cite conference |title=Fast Generation of High-Quality Pseudorandom Numbers and Permutations Using MPI and OpenMP on the Cray XD1 |first1=Stephen |last1=Bique |first2=Robert |last2=Rosenberg |date=May 2009 |conference=Cray User Group 2009 |url=https://cug.org/5-publications/proceedings_attendee_lists/CUG09CD/S09_Proceedings/pages/authors/06-10Tuesday/8A-Bique/Bique.htm |quote=The die is determined using modular arithmetic, e.g., <code>lrand48() % 6 + 1</code>,&nbsp;... The CRAY RANF function only rolls three of the six possible outcomes (which three sides depends on the seed)!}}</ref>

To achieve this period, the multiplier must satisfy ''a''&nbsp;≡&nbsp;±3 (mod&nbsp;8),<ref>{{cite journal
 |first=Martin |last=Greenberger
 |title=Notes on a New Pseudo-Random Number Generator
 |journal=Journal of the ACM
 |volume=8 |issue=2 |date=April 1961 |pages=163–167
 |doi=10.1145/321062.321065
|s2cid=17196519
 |doi-access=free}}</ref> and the seed ''X''{{sub|0}} must be odd.

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. For example, a modulus of [[Fermat number|F{{sub|5}}]] = 2{{sup|32}}&nbsp;+&nbsp;1 might seem attractive, as the outputs can be easily mapped to a 32-bit word 0 ≤ ''X{{sub|i}}''&nbsp;−&nbsp;1 &lt; 2{{sup|32}}. However, a seed of ''X''{{sub|0}}&nbsp;=&nbsp;6700417 (which divides 2{{sup|32}}&nbsp;+&nbsp;1) or any multiple would lead to an output with a period of only 640.

A more popular implementation for large periods is a [[combined linear congruential generator]]; combining (e.g. by summing their outputs) several generators is equivalent to the output of a single generator whose modulus is the product of the component generators' moduli.<ref name=mcomp91>{{cite journal
 |title=Structural Properties for Two Classes of Combined Random Number Generators
 |first1=Pierre |last1=L'Ecuyer |first2=Shu |last2=Tezuka
 |journal=Mathematics of Computation |volume=57 |issue=196 |date=October 1991 |pages=735–746
 |doi=10.2307/2938714
 |jstor=2938714
}}</ref> and whose period is the [[least common multiple]] of the component periods. Although the periods will share a [[common divisor]] of 2, the moduli can be chosen so that is the only common divisor and the resultant period is (''m''{{sub|1}}&nbsp;−&nbsp;1)(''m''{{sub|2}}&nbsp;−&nbsp;1)···(''m{{sub|k}}''&nbsp;−&nbsp;1)/2{{sup|''k''−1}}.{{r|cacm88|p=744}} One example of this is the [[Wichmann–Hill]] generator.

== Relation to LCG ==

While the Lehmer RNG can be viewed as a particular case of the [[linear congruential generator]] with {{gaps|''c''|{{=}}|0}}, it is a special case that implies certain restrictions and properties. In particular, for the Lehmer RNG, the initial seed ''X''{{sub|0}} must be [[coprime]] to the modulus ''m'', which 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''&nbsp;−&nbsp;1, and it is such when ''m'' is prime and ''a'' is a primitive root modulo ''m''.

On the other hand, the [[discrete logarithm]]s (to base ''a'' or any primitive root modulo ''m'') of ''X{{sub|k}}'' in <math>\mathbb{Z}_m</math> represent a linear congruential sequence modulo the [[Euler totient]] <math>\varphi(m)</math>.

== Implementation ==

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 prime]]s 2<sup>31</sup>&nbsp;−&nbsp;1 and 2<sup>61</sup>&nbsp;−&nbsp;1 are popular, as are 2<sup>32</sup>&nbsp;−&nbsp;5 and 2<sup>64</sup>&nbsp;−&nbsp;59), reduction modulo {{nobr|1=''m'' = 2<sup>''e''</sup> − ''d''}} can be implemented more cheaply than a general double-width division using the identity {{nobr|2<sup>''e''</sup> ≡ ''d'' (mod ''m'')}}.

The basic reduction step divides the product into two ''e''-bit parts, multiplies the high part by ''d'', and adds them: {{nobr|(''ax'' mod 2<sup>''e''</sup>) + ''d''&thinsp;{{floor|''ax''/2<sup>''e''</sup>}}}}. This can be followed by subtracting ''m'' until the result is in range. The number of subtractions is limited to ''ad''/''m'', which can be easily limited to one if ''d'' is small and {{nobr|''a'' &lt; ''m''/''d''}} is chosen. (This condition also ensures that ''d''&thinsp;{{floor|''ax''/2<sup>''e''</sup>}} is a single-width product; if it is violated, a double-width product must be computed.)

When the modulus is a Mersenne prime (''d''&nbsp;=&nbsp;1), the procedure is particularly simple. Not only is multiplication by ''d'' trivial, but the [[Conditional (computer programming)|conditional]] subtraction can be replaced by an unconditional shift and addition. To see this, note that the algorithm guarantees that {{nobr|''x'' ≢ 0 (mod ''m'')}}, meaning that ''x''&nbsp;=&nbsp;0 and ''x''&nbsp;=&nbsp;''m'' are both impossible. This avoids the need to consider equivalent ''e''-bit representations of the state; only values where the high bits are non-zero need reduction.

The low ''e'' bits of the product ''ax'' cannot represent a value larger than ''m'', and the high bits will never hold a value greater than ''a''&nbsp;−&nbsp;1&nbsp;≤&nbsp;m&nbsp;−&nbsp;2. Thus the first reduction step produces a value at most ''m''&nbsp;+&nbsp;''a''&nbsp;−&nbsp;1 ≤ 2''m''&nbsp;−&nbsp;2 = 2<sup>''e''+1</sup>&nbsp;−&nbsp;4. This is an (''e''&nbsp;+&nbsp;1)-bit number, which can be greater than ''m'' (i.e. might have bit ''e'' set), but the high half is at most 1, and if it is, the low ''e'' bits will be strictly less than ''m''. Thus whether the high bit is 1 or 0, a second reduction step (addition of the halves) will never overflow ''e'' bits, and the sum will be the desired value.

If ''d''&nbsp;&gt;&nbsp;1, conditional subtraction can also be avoided, but the procedure is more intricate. The fundamental challenge of a modulus like 2<sup>32</sup>&nbsp;−&nbsp;5 lies in ensuring that we produce only one representation for values such as 1&nbsp;≡&nbsp;2<sup>32</sup>&nbsp;−&nbsp;4. The solution is to temporarily add ''d'', so that the range of possible values is ''d'' through 2<sup>''e''</sup>&nbsp;−&nbsp;1, and reduce values larger than ''e'' bits in a way that never generates representations less than ''d''. Finally subtracting the temporary offset produces the desired value.

Begin by assuming that we have a partially reduced value ''y'' bounded so that 0&nbsp;≤&nbsp;''y''&nbsp;&lt; 2''m'' =&nbsp;2<sup>''e''+1</sup>&nbsp;−&nbsp;2''d''. In this case, a single offset subtraction step will produce 0&nbsp;≤&nbsp;''{{prime|y}}''  {{nobr|{{=}} ((''y''&nbsp;+&nbsp;''d'') mod 2<sup>''e''</sup>) + ''d''&thinsp;{{floor|(''y''&nbsp;+&nbsp;''d'')/2<sup>''e''</sup>}} − ''d''}} &lt;&nbsp;''m''. To see this, consider two cases:
; 0 ≤ ''y'' &lt; ''m'' = 2<sup>''e''</sup> − ''d''
: In this case, ''y'' + ''d'' &lt; 2<sup>''e''</sup> and ''{{prime|y}}''&nbsp;=&nbsp;''y''&nbsp;&lt;&nbsp;''m'', as desired.
; ''m'' ≤ ''y'' &lt; 2''m''
: In this case, 2<sup>''e''</sup> ≤&nbsp;''y''&nbsp;+&nbsp;''d'' &lt;&nbsp;2<sup>''e''+1</sup> is an (''e''&nbsp;+&nbsp;1)-bit number, and {{floor|(''y''&nbsp;+&nbsp;''d'')/2<sup>''e''</sup>}}&nbsp;=&nbsp;1. Thus, ''{{prime|y}}''&nbsp;=&nbsp;(''y''&nbsp;+&nbsp;''d'')&nbsp;−&nbsp;2<sup>''e''</sup>&nbsp;+&nbsp;''d''&nbsp;−&nbsp;''d'' =&nbsp;''y''&nbsp;−&nbsp;2<sup>''e''</sup>&nbsp;+&nbsp;''d'' =&nbsp;''y''&nbsp;−&nbsp;''m''&nbsp;&lt;&nbsp;''m'', as desired. Because the multiplied high part is ''d'', the sum is at least ''d'', and subtracting the offset never causes underflow.

(For the case of a Lehmer generator specifically, a zero state or its image ''y''&nbsp;=&nbsp;''m'' will never occur, so an offset of ''d''&nbsp;−&nbsp;1 will work just the same, if that is more convenient. This reduces the offset to 0 in the Mersenne prime case, when ''d''&nbsp;=&nbsp;1.)

Reducing a larger product ''ax'' to less than 2''m''&nbsp;=&nbsp;2<sup>''e''+1</sup>&nbsp;−&nbsp;2''d'' can be done by one or more reduction steps without an offset.

If ''ad''&nbsp;≤&nbsp;''m'', then one additional reduction step suffices. Since ''x''&nbsp;&lt;&nbsp;''m'', ''ax''&nbsp;&lt;&nbsp;''am'' ≤&nbsp;(''a''&nbsp;−&nbsp;1)2<sup>''e''</sup>, and one reduction step converts this to at most 2<sup>''e''</sup>&nbsp;−&nbsp;1&nbsp;+&nbsp;(''a''&nbsp;−&nbsp;1)''d''&nbsp;=&nbsp;''m''&nbsp;+&nbsp;''ad''&nbsp;−&nbsp;1. This is within the limit of 2''m'' if ''ad''&nbsp;−&nbsp;1&nbsp;&lt;&nbsp;''m'', which is the initial assumption.

If ''ad''&nbsp;&gt;&nbsp;''m'', then it is possible for the first reduction step to produce a sum greater than 2''m''&nbsp;=&nbsp;2<sup>''e''+1</sup>&nbsp;−&nbsp;2''d'', which is too large for the final reduction step. (It also requires the multiplication by ''d'' to produce a product larger than ''e'' bits, as mentioned above.) However, as long as ''d''<sup>2</sup>&nbsp;&lt;&nbsp;2<sup>''e''</sup>,{{^|This bound is not exactly tight; it's more like ''ad''<sup>2</sup> < 2<sup>2''e''</sup>.}} the first reduction will produce a value in the range required for the preceding case of two reduction steps to apply.

== Schrage's method ==

If a double-width product is not available, '''Schrage's method''',<ref>{{cite journal
 |title=A More Portable Fortran Random Number Generator
 |first=Linus |last=Schrage
 |journal=ACM Transactions on Mathematical Software
 |volume=5 |issue=2 |date=June 1979 |pages=132–138 
 |doi=10.1145/355826.355828 |citeseerx=10.1.1.470.6958
 |s2cid=14090729 |url=http://degiorgi.math.hr/aaa_sem/Rand_Gen/p132-schrage.pdf
}}</ref><ref>{{cite web
 |title=Computer Systems Performance Analysis Chapter 26: Random-Number Generation
 |first=Raj |last=Jain |date=9 July 2010
 |pages=19–22
 |url=http://www.cse.wustl.edu/~jain/iucee/ftp/k_26rng.pdf#page=19
 |access-date=2017-10-31
}}</ref> also called the approximate factoring method,<ref name=toms91>{{cite journal
 |title=Implementing a random number package with splitting facilities
 |first=Pierre |last1=L'Ecuyer |first2=Serge |last2=Côté
 |journal=ACM Transactions on Mathematical Software |volume=17 |issue=1 |date=March 1991 |pages=98–111
 |doi=10.1145/103147.103158
|s2cid=14385204 |doi-access=free }}  This explores several different implementations of modular multiplication by a constant.</ref> may be used to compute {{nobr|''ax'' mod ''m''}}, but this comes at the cost:
* The modulus must be representable in a [[signed integer]]; arithmetic operations must allow a range of ±''m''.
* The choice of multiplier ''a'' is restricted.  We require that {{nobr|''m'' mod ''a'' ≤ {{floor|''m''/''a''}}}}, commonly achieved by choosing ''a''&nbsp;≤&nbsp;{{sqrt|''m''}}.
* One division (with remainder) per iteration is required.

While this technique is popular for portable implementations in [[high-level language]]s which lack double-width operations,{{r|cacm88|p=744}} on modern computers [[Division algorithm#Division by a constant|division by a constant]] is usually implemented using double-width multiplication, so this technique should be avoided if efficiency is a concern.  Even in high-level languages, if the multiplier ''a'' is limited to {{sqrt|''m''}}, then the double-width product ''ax'' can be computed using two single-width multiplications, and reduced using the techniques described above.

To use Schrage's method, first factor {{nobr|1=''m'' = ''qa'' + ''r''}}, i.e. precompute the auxiliary constants {{nobr|1=''r'' = ''m'' mod ''a''}} and {{nobr|1=''q'' = {{floor|''m''/''a''}}}} = (''m''−''r'')/''a''. Then, each iteration, compute {{nobr|''ax'' ≡ ''a''(''x'' mod ''q'') − ''r''&thinsp;{{floor|''x''/''q''}} (mod ''m'')}}.

This equality holds because
:<math>
\begin{align}
aq &= a \cdot (m - r) / a \\
   &= (a / a) \cdot (m - r) \\
   &= (m - r) \\
   &\equiv -r \pmod m \\
\end{align}
</math>
so if we factor ''x'' = (''x'' mod ''q'') + ''q''{{floor|''x''/''q''}}, we get:
:<math>
\begin{align}
ax &= a(x \bmod q) + aq\lfloor x/q \rfloor \\
   &= a(x \bmod q) + (m-r)\lfloor x/q \rfloor \\
   &\equiv a(x \bmod q) - r\lfloor x/q \rfloor \pmod m \\
\end{align}
</math>

The reason it does not overflow is that both terms are less than ''m''.  Since ''x''&nbsp;mod&nbsp;''q'' < ''q'' ≤ ''m''/''a'', the first term is strictly less than ''am''/''a''&nbsp;=&nbsp;''m'' and may be computed with a single-width product.

If ''a'' is chosen so that ''r''&nbsp;≤&nbsp;''q'' (and thus ''r''/''q''&nbsp;≤&nbsp;1), then the second term is also less than ''m'': ''r''&thinsp;{{floor|''x''/''q''}} ≤ ''rx''/''q'' = ''x''(''r''/''q'') ≤ ''x''(1) = ''x'' < ''m''. Thus, the difference lies in the range [1−''m'',&nbsp;''m''−1] and can be reduced to [0,&nbsp;''m''−1] with a single conditional add.<ref>{{cite web
 |title=Schrage's Method
 |url=http://home.earthlink.net/~pfenerty/pi/schrages_method.html
 |first=Paul |last=Fenerty |date=11 September 2006
 |access-date=2017-10-31
}}</ref>

This technique may be extended to allow a negative ''r'' (−''q''&nbsp;≤&nbsp;''r''&nbsp;&lt;&nbsp;0), changing the final reduction to a conditional subtract.

The technique may also be extended to allow larger ''a'' by applying it recursively.{{r|toms91|p=102}}  Of the two terms subtracted to produce the final result, only the second (''r''&thinsp;{{floor|''x''/''q''}}) risks overflow.  But this is itself a modular multiplication by a [[Constant folding|compile-time constant]] ''r'', and may be implemented by the same technique.  Because each step, on average, halves the size of the multiplier (0&nbsp;≤&nbsp;''r''&nbsp;&lt;&nbsp;''a'', average value (''a''−1)/2), this would appear to require one step per bit and be spectacularly inefficient.  However, each step also divides ''x'' by an ever-increasing quotient {{nobr|1=''q'' = {{floor|''m''/''a''}}}}, and quickly a point is reached where the argument is 0 and the recursion may be terminated.

== Sample C99 code ==
Using [[ISO C|C]] code, the Park-Miller RNG can be written as follows:
<syntaxhighlight lang="c">
uint32_t lcg_parkmiller(uint32_t *state)
{
	return *state = (uint64_t)*state * 48271 % 0x7fffffff;
}
</syntaxhighlight>
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.

To avoid the 64-bit division, do the reduction by hand:
<syntaxhighlight lang="c">
uint32_t lcg_parkmiller(uint32_t *state)
{
	uint64_t product = (uint64_t)*state * 48271;
	uint32_t x = (product & 0x7fffffff) + (product >> 31);

	x = (x & 0x7fffffff) + (x >> 31);
	return *state = x;
}
</syntaxhighlight>

To use only 32-bit arithmetic, use Schrage's method:
<syntaxhighlight lang="c">
uint32_t lcg_parkmiller(uint32_t *state)
{
	// Precomputed parameters for Schrage's method
	const uint32_t M = 0x7fffffff;
	const uint32_t A = 48271;
	const uint32_t Q = M / A;    // 44488
	const uint32_t R = M % A;    //  3399

	uint32_t div = *state / Q;	// max: M / Q = A = 48,271
	uint32_t rem = *state % Q;	// max: Q - 1     = 44,487

	int32_t s = rem * A;	// max: 44,487 * 48,271 = 2,147,431,977 = 0x7fff3629
	int32_t t = div * R;	// max: 48,271 *  3,399 =   164,073,129
	int32_t result = s - t;

	if (result < 0)
		result += M;

	return *state = result;
}
</syntaxhighlight>
or use two 16×16-bit multiplies:
<syntaxhighlight lang="c">
uint32_t lcg_parkmiller(uint32_t *state)
{
	const uint32_t A = 48271;

	uint32_t low  = (*state & 0x7fff) * A;			// max: 32,767 * 48,271 = 1,581,695,857 = 0x5e46c371
	uint32_t high = (*state >> 15)    * A;			// max: 65,535 * 48,271 = 3,163,439,985 = 0xbc8e4371
	uint32_t x = low + ((high & 0xffff) << 15) + (high >> 16);	// max: 0x5e46c371 + 0x7fff8000 + 0xbc8e = 0xde46ffff

	x = (x & 0x7fffffff) + (x >> 31);
	return *state = x;
}
</syntaxhighlight>

Another popular Lehmer generator uses the prime modulus 2{{sup|32}}−5:
<syntaxhighlight lang="c">
uint32_t lcg_rand(uint32_t *state)
{
    return *state = (uint64_t)*state * 279470273u % 0xfffffffb;
}
</syntaxhighlight>

This can also be written without a 64-bit division:
<syntaxhighlight lang="c">
uint32_t lcg_rand(uint32_t *state)
{
	uint64_t product = (uint64_t)*state * 279470273u;
	uint32_t x;

	// Not required because 5 * 279470273 = 0x5349e3c5 fits in 32 bits.
	// product = (product & 0xffffffff) + 5 * (product >> 32);
	// A multiplier larger than 0x33333333 = 858,993,459 would need it.

	// The multiply result fits into 32 bits, but the sum might be 33 bits.
	product = (product & 0xffffffff) + 5 * (uint32_t)(product >> 32);

	product += 4;
	// This sum is guaranteed to be 32 bits.
	x = (uint32_t)product + 5 * (uint32_t)(product >> 32);
	return *state = x - 4;
}
</syntaxhighlight>

Many other Lehmer generators have good properties. The following modulo-2<sup>128</sup> Lehmer generator requires 128-bit support from the compiler and uses a multiplier computed by L'Ecuyer.<ref>{{cite journal
 |first=Pierre |last=L’Ecuyer
 |title=Tables of linear congruential generators of different sizes and good lattice structure
 |journal=[[Mathematics of Computation]] |date=January 1999 |volume=68 |number=225 |pages=249–260
 |url=https://www.ams.org/journals/mcom/1999-68-225/S0025-5718-99-00996-5/S0025-5718-99-00996-5.pdf
 |doi=10.1090/s0025-5718-99-00996-5 |citeseerx=10.1.1.34.1024
}}</ref> It has a period of 2<sup>126</sup>:
<syntaxhighlight lang="c">
static unsigned __int128 state;

/* The state must be seeded with an odd value. */
void seed(unsigned __int128 seed)
{
	state = seed << 1 | 1;
}

uint64_t next(void)
{
	// GCC cannot write 128-bit literals, so we use an expression
	const unsigned __int128 mult =
		(unsigned __int128)0x12e15e35b500f16e << 64 |
		0x2e714eb2b37916a5;
	state *= mult;
	return state >> 64;
}
</syntaxhighlight>
The generator computes an odd 128-bit value and returns its upper 64 bits.

This generator passes BigCrush from [[TestU01]], but fails the TMFn test from [http://pracrand.sourceforge.net/ PractRand]. 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 2<sup>62</sup>, rather than 2<sup>126</sup>. [[Linear congruential generator]]s with a power-of-2 modulus have a similar behavior.

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):
<syntaxhighlight lang="c">
uint64_t next(void)
{
	uint64_t result = state >> 64;
	// GCC cannot write 128-bit literals, so we use an expression
	const unsigned __int128 mult =
		(unsigned __int128)0x12e15e35b500f16e << 64 |
		0x2e714eb2b37916a5;
	state *= mult;
	return result;
}
</syntaxhighlight>
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. 
== References ==
{{reflist|2}}

* {{cite book |last=Lehmer |first=D. H. |contribution=Mathematical methods in large-scale computing units |url=https://archive.org/details/proceedings_of_a_second_symposium_on_large-scale_ |title=Proceedings of a Second Symposium on Large-Scale Digital Calculating Machinery |year=1949 |pages=[https://archive.org/details/proceedings_of_a_second_symposium_on_large-scale_/page/n178 141]–146 |mr=0044899}} (journal version: ''Annals of the Computation Laboratory of Harvard University'', Vol. 26 (1951)).
* Steve Park, [http://www.cs.wm.edu/~va/software/park/ Random Number Generators]

== External links ==
* [https://primes.utm.edu/lists/2small/0bit.html Primes just less than a power of two] may be useful for choosing moduli.  Part of [[Prime Pages]].

[[Category:Pseudorandom number generators]]
[[Category:Modular arithmetic]]