Editing Linear congruential generator (section)

== 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.<ref name=History>{{cite conference
 |title=History of Uniform Random Number Generation
 |first=Pierre |last=L'Ecuyer
 |date=13 July 2017<!--Conference is in December 2017-->
 |conference=Proceedings of the 2017 Winter Simulation Conference (to appear)
 |editor1-first=W. K. V. |editor1-last=Chan   |editor2-first=A. |editor2-last=D'Ambrogio
 |editor3-first=G. |editor3-last=Zacharewicz  |editor4-first=N. |editor4-last=Mustafee
 |editor5-first=G. |editor5-last=Wainer       |editor6-first=E. |editor6-last=Page
 |url=https://www.iro.umontreal.ca/~lecuyer/myftp/papers/wsc17rng-history.pdf
 |location=Las Vegas, United States
 |id=[https://hal.inria.fr/hal-01561551 hal-01561551]
}}</ref>

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''.{{r|KnuthV2|Steele20|Marsaglia68|Park88|Hörmann92|LEcuyer99}}  For example, ''a''&nbsp;=&nbsp;1 and ''c''&nbsp;=&nbsp;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.<ref name="Press">{{cite book |last=Press |first=William H. |year=1992 |title=Numerical Recipes in Fortran 77: The Art of Scientific Computing |edition=2nd |isbn=978-0-521-43064-7 |display-authors=etal |title-link=Numerical Recipes |page=268}}</ref>{{r|Park88|p=1198–9}}

There are three common families of parameter choice:

=== ''m'' prime, ''c'' = 0 ===
{{main|Lehmer random number generator}}
This is the original Lehmer RNG construction. The period is ''m''−1 if the multiplier ''a'' is chosen to be a [[Primitive element (finite field)|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 prime]]s 2<sup>31</sup>−1 and 2<sup>61</sup>−1 are popular), so that the reduction modulo ''m''&nbsp;=&nbsp;2<sup>''e''</sup>&nbsp;−&nbsp;''d'' can be computed as (''ax''&nbsp;mod&nbsp;2<sup>''e''</sup>)&nbsp;+&nbsp;''d''&nbsp;{{floor|''ax''/2<sup>''e''</sup>}}. 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'''<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 |pages=132–138 |date=June 1979
 |doi=10.1145/355826.355828
 |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–20
 |url=http://www.cse.wustl.edu/~jain/iucee/ftp/k_26rng.pdf#page=19
 |access-date=2017-10-31
}}</ref> may be used. To do this, factor ''m''&nbsp;=&nbsp;''qa''+''r'', i.e. {{nobr|1=''q'' = {{floor|''m''/''a''}}}} and {{nobr|1=''r'' = ''m'' mod ''a''}}. Then compute ''ax''&nbsp;mod&nbsp;''m'' = {{nobr|''a''(''x'' mod ''q'') − ''r''{{floor|''x''/''q''}}}}. Since ''x''&nbsp;mod&nbsp;''q'' < ''q'' ≤ ''m''/''a'', the first term is strictly less than ''am''/''a''&nbsp;=&nbsp;''m''. 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''{{floor|''x''/''q''}} ≤ ''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'',&nbsp;''m''−1], so 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
 |archive-url=https://web.archive.org/web/20181230061306/http://home.earthlink.net/~pfenerty/pi/schrages_method.html
 |archive-date=2018-12-30 |url-status=dead
}}</ref>

The most expensive operation in Schrage's method is the division (with remainder) of ''x'' by ''q''; fast [[Division algorithm#Division by a constant|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&nbsp;≤&nbsp;''x''&nbsp;<&nbsp;''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'' = 2<sup>32</sup> or ''m'' = 2<sup>64</sup>, 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''&nbsp;≡&nbsp;±3 (mod 8) and the initial state ''X''<sub>0</sub> 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''&nbsp;≡&nbsp;+3, ''X'' alternates ±1↔±3, while if ''a''&nbsp;≡&nbsp;−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.{{r|KnuthV2}}

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-2<sup>''b''</sup> LCG by themselves, repeating with a period of 2<sup>''b''−2</sup>. 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]]:{{r|KnuthV2|p=17–19}}
# <math>m</math> and <math>c</math> are coprime,
# <math>a - 1</math> is divisible by all [[prime factor]]s of <math>m</math>,
# <math>a - 1</math> is divisible by 4 if <math>m</math> is divisible by 4.

These three requirements are referred to as the Hull–Dobell Theorem.<ref>{{Cite journal
 |last1=Hull |first1=T. E. |last2=Dobell |first2=A. R.
 |date=July 1962
 |title=Random Number Generators
 |url=http://chagall.med.cornell.edu/BioinfoCourse/PDFs/Lecture4/random_number_generator.pdf
 |journal=SIAM Review |volume=4 |issue=3 |pages=230–254 
 |access-date=2016-06-26 |doi=10.1137/1004061
|bibcode=1962SIAMR...4..230H |hdl=1828/3142 |hdl-access=free }}</ref><ref name="HullDobell">{{cite book
 |title=System Modeling and Simulation
 |publisher=John Wiley & Sons, Ltd.
 |last=Severance |first=Frank |year=2001 |page=86
 |isbn=978-0-471-49694-6
}}</ref>

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''&nbsp;≡&nbsp;1 (mod&nbsp;''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.{{r|Park88|p=1199}}  For example, it is desirable for ''a''&nbsp;−&nbsp;1 to not be any more divisible by prime factors of ''m'' than necessary.  If ''m'' is a power of 2, then ''a''&nbsp;−&nbsp;1 should be divisible by 4 but not divisible by 8, i.e.&nbsp;''a''&nbsp;≡&nbsp;5&nbsp;(mod&nbsp;8).{{r|KnuthV2|at=§3.2.1.3}}

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{{r|KnuthV2|at=§3.3.3}} is quite challenging.{{r|Park88}}  The [[spectral test]] is one of the most important tests.<ref name=Austin2008>{{cite web |first=David |last=Austin |title=Random Numbers: Nothing Left to Chance |date=March 2008 |publisher=American Mathematical Society |work=Feature Column |url=https://www.ams.org/publicoutreach/feature-column/fcarc-random}}</ref>

Note that a power-of-2 modulus shares the problem as described above for ''c''&nbsp;=&nbsp;0: the low ''k'' bits form a generator with modulus 2<sup>''k''</sup> and thus repeat with a period of 2<sup>''k''</sup>; only the most significant bit achieves the full period.  If a pseudorandom number less than ''r'' is desired, {{floor|''rX''/''m''}} is a much higher-quality result than ''X'' mod ''r''.  Unfortunately, most programming languages make the latter much easier to write (<code>X % r</code>), 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.{{r|KnuthV2|p=11}}  Specifically, if ''Y'' is the prototypical sequence defined by ''Y''<sub>0</sub> =&nbsp;0 and ''Y''<sub>''n''+1</sub> =&nbsp;''aY<sub>n</sub>''&nbsp;+&nbsp;1 mod&nbsp;m, then a general sequence ''X''<sub>''n''+1</sub> =&nbsp;''aX<sub>n</sub>''&nbsp;+&nbsp;''c'' mod&nbsp;''m'' can be written as an affine function of ''Y'':
:<math>X_n = (X_0(a-1)+c)Y_n + X_0 = (X_1 - X_0)Y_n + X_0 \pmod m.</math>
More generally, any two sequences ''X'' and ''Z'' with the same multiplier and modulus are related by
:<math>{ X_n - X_0 \over X_1 - X_0 } = Y_n = {a^n - 1 \over a - 1} = { Z_n - Z_0 \over Z_1 - Z_0 } \pmod m.</math>

In the common case where ''m'' is a power of 2 and ''a''&nbsp;≡&nbsp;5&nbsp;(mod&nbsp;8) (a desirable property for other reasons), it is always possible to find an initial value ''X''<sub>0</sub> so that the denominator ''X''<sub>1</sub>&nbsp;−&nbsp;''X''<sub>0</sub> ≡&nbsp;±1&nbsp;(mod&nbsp;''m''), producing an even simpler relationship.  With this choice of ''X''<sub>0</sub>, ''X<sub>n</sub>'' =&nbsp;''X''<sub>0</sub>&nbsp;±&nbsp;''Y<sub>n</sub>'' will remain true for all ''n''.{{r|Steele20|p=10-11}} The sign is determined by ''c''&nbsp;≡&nbsp;±1&nbsp;(mod&nbsp;4), and the constant ''X''<sub>0</sub> is determined by 1&nbsp;∓&nbsp;''c''&nbsp;≡&nbsp;(1&nbsp;−&nbsp;''a'')''X''<sub>0</sub>&nbsp;(mod&nbsp;''m'').<!--1∓c is a multiple of 4, and 1−a ≡ 4 (mod 8), so this has solutions mod m/4-->

As a simple example, consider the generators ''X''<sub>''n''+1</sub> =&nbsp;157''X<sub>n</sub>''&nbsp;+&nbsp;3 mod&nbsp;256 and ''Y''<sub>''n''+1</sub> =&nbsp;157''Y<sub>n</sub>''&nbsp;+&nbsp;1 mod&nbsp;256; i.e. ''m''&nbsp;=&nbsp;256, ''a''&nbsp;=&nbsp;157, and ''c''&nbsp;=&nbsp;3.  Because 3&nbsp;≡&nbsp;−1&nbsp;(mod&nbsp;4), we are searching for a solution to 1&nbsp;+&nbsp;3&nbsp;≡&nbsp;(1&nbsp;−&nbsp;157)''X''<sub>0</sub>&nbsp;(mod&nbsp;256).  This is satisfied by ''X''<sub>0</sub>&nbsp;≡&nbsp;41&nbsp;(mod&nbsp;64), so if we start with that, then ''X<sub>n</sub>'' ≡&nbsp;''X''<sub>0</sub>&nbsp;−&nbsp;''Y<sub>n</sub>''&nbsp;(mod&nbsp;256) for all ''n''.

For example, using ''X''<sub>0</sub> = 233 = 3&times;64 + 41:
* ''X'' = 233, 232, 75, 2, 61, 108, ...<!--63, 166, 209, 48, 115, 138, 165, 52, 231-->
* ''Y'' = 0, 1, 158, 231, 172, 125, ...<!--170, 67, 24, 185, 118, 95, 68, 181, 2-->
* ''X''&nbsp;+&nbsp;''Y'' mod 256 = 233, 233, 233, 233, 233, 233, ...