Namespaces
Variants
Actions

Difference between revisions of "Linear congruential method"

From Encyclopedia of Mathematics
Jump to: navigation, search
(Importing text file)
 
m (AUTOMATIC EDIT (latexlist): Replaced 102 formulas out of 102 by TEX code with an average confidence of 2.0 and a minimal confidence of 2.0.)
Line 1: Line 1:
A method widely used for generating random numbers from the [[Uniform distribution|uniform distribution]]: A sequence of integers is initialized with a value <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l1300601.png" /> and continued as
+
<!--This article has been texified automatically. Since there was no Nroff source code for this article,
 +
the semi-automatic procedure described at https://encyclopediaofmath.org/wiki/User:Maximilian_Janisch/latexlist
 +
was used.
 +
If the TeX and formula formatting is correct, please remove this message and the {{TEX|semi-auto}} category.
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l1300602.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a1)</td></tr></table>
+
Out of 102 formulas, 102 were replaced by TEX code.-->
  
for all <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l1300603.png" />. The fractions <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l1300604.png" /> are the derived pseudo-random numbers in the interval <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l1300605.png" /> (cf. also [[Random and pseudo-random numbers|Random and pseudo-random numbers]]; [[Pseudo-random numbers|Pseudo-random numbers]]). The constants <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l1300606.png" />, the modulus, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l1300607.png" />, the multiplicator, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l1300608.png" />, the increment, and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l1300609.png" />, the starting number, are suitably chosen non-negative integers. Three choices of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006010.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006011.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006012.png" /> are common on most computers:
+
{{TEX|semi-auto}}{{TEX|done}}
 +
A method widely used for generating random numbers from the [[Uniform distribution|uniform distribution]]: A sequence of integers is initialized with a value $z_0$ and continued as
  
1) <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006013.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006014.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006015.png" />, and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006016.png" />. All <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006017.png" /> are generated.
+
\begin{equation} \tag{a1} z _ { i  + 1} \equiv a z _ { i } + r ( \operatorname { mod } m ) ,\, 0 \leq z _ { i } &lt; m, \end{equation}
  
2) <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006018.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006019.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006020.png" /> prime, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006021.png" /> a [[Primitive root|primitive root]] modulo <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006022.png" />. All <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006023.png" /> are generated.
+
for all $i$. The fractions $u _ { i } = z _ { i } / m$ are the derived pseudo-random numbers in the interval $[ 0,1 )$ (cf. also [[Random and pseudo-random numbers|Random and pseudo-random numbers]]; [[Pseudo-random numbers|Pseudo-random numbers]]). The constants $m$, the modulus, $a$, the multiplicator, $r$, the increment, and $z_0$, the starting number, are suitably chosen non-negative integers. Three choices of $m$, $a$ and $r$ are common on most computers:
  
3) <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006024.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006025.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006026.png" />. All integers <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006027.png" /> are generated.
+
1) $r = 0$, $m = 2 ^ { E }$, $a \equiv 5 ( \operatorname { mod } 8 )$, and $z _ { 0 } \equiv 1 ( \operatorname { mod } 4 )$. All $z _ { i } \equiv 1 ( \operatorname { mod } 4 )$ are generated.
  
For selecting  "good"  random number generators one has to study the distribution of the <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006028.png" />-tuplets <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006029.png" />. Geometrically these <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006030.png" /> may be considered as points of a lattice <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006031.png" /> in the <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006032.png" />-dimensional hypercube <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006033.png" />. The lattice points can also be seen as intersection points of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006034.png" /> sets of parallel hyperplanes. Consequently, the following questions may be raised:
+
2) $r = 0$, $m = p$, $p$ prime, $a$ a [[Primitive root|primitive root]] modulo $p$. All $z _ { i } = 1 , \dots , p - 1$ are generated.
  
i) Determine the minimal number <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006035.png" /> of parallel hyperplanes on which all points <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006036.png" /> lie.
+
3) $r \equiv 1 ( \operatorname { mod } 2 )$, $m = 2 ^ { E }$, $a \equiv 1 ( \operatorname { mod } 4 )$. All integers $0 , \ldots , 2 ^ { E } - 1$ are generated.
  
ii) Determine the maximal distance <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006037.png" /> of parallel hyperplanes on which all points <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006038.png" /> lie.
+
For selecting  "good" random number generators one has to study the distribution of the $k$-tuplets $P _ { k } = ( u _ { i  + 1} , \dots , u _ { i  + k})$. Geometrically these $P _ { k }$ may be considered as points of a lattice $G$ in the $k$-dimensional hypercube $[ 0,1 ) ^ { k }$. The lattice points can also be seen as intersection points of $k$ sets of parallel hyperplanes. Consequently, the following questions may be raised:
  
Question i) was asked by G. Marsaglia [[#References|[a12]]], who derived upper bounds for <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006039.png" /> using Minkowski's convex body theorem (cf. also [[Minkowski theorem|Minkowski theorem]]). The  "wave numbers"  <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006040.png" /> were introduced by R.R. Coveyou and R.D. MacPherson [[#References|[a4]]]. Their algorithm for calculating <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006041.png" /> was simplified by D.E. Knuth [[#References|[a10]]].
+
i) Determine the minimal number $N _ { k } ^ { * }$ of parallel hyperplanes on which all points $P _ { k }$ lie.
  
The calculation of both quantities is based on a general procedure to determine non-zero vectors of shortest length in the dual lattice of covering hyperplanes. For the determination of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006042.png" /> the <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006043.png" />-norm is used, and for <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006044.png" /> the Euclidean norm is appropriate. The algorithm of U. Dieter (1973) gives exact values for both quantities; no exact values for <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006045.png" /> were known before. Knuth included a variant of this algorithm in [[#References|[a10]]]. A completely different approach was proposed by L. Lovász, J.K. Lenstra and H.W. Lenstra, Jr., called the LLL-algorithm (cf. also [[LLL basis reduction method|LLL basis reduction method]]). In the case of the Euclidean norm, the final search can be shortened by an idea of U. Finke and M. Pohst [[#References|[a8]]].
+
ii) Determine the maximal distance $D _ { k } ^ { * }$ of parallel hyperplanes on which all points $P _ { k }$ lie.
  
For any sequence <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006046.png" /> of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006047.png" />-uniformly distributed random numbers, the local deviations
+
Question i) was asked by G. Marsaglia [[#References|[a12]]], who derived upper bounds for $N _ { k } ^ { * }$ using Minkowski's convex body theorem (cf. also [[Minkowski theorem|Minkowski theorem]]). The  "wave numbers" $W _ { k } ^ { * } = 1 / D _ { k } ^ { * }$ were introduced by R.R. Coveyou and R.D. MacPherson [[#References|[a4]]]. Their algorithm for calculating $W _ { k } ^ { * }$ was simplified by D.E. Knuth [[#References|[a10]]].
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006048.png" /></td> </tr></table>
+
The calculation of both quantities is based on a general procedure to determine non-zero vectors of shortest length in the dual lattice of covering hyperplanes. For the determination of $N _ { k } ^ { * }$ the $\mathbf{l}_{1}$-norm is used, and for $D _ { k } ^ { * }$ the Euclidean norm is appropriate. The algorithm of U. Dieter (1973) gives exact values for both quantities; no exact values for $N _ { k } ^ { * }$ were known before. Knuth included a variant of this algorithm in [[#References|[a10]]]. A completely different approach was proposed by L. Lovász, J.K. Lenstra and H.W. Lenstra, Jr., called the LLL-algorithm (cf. also [[LLL basis reduction method|LLL basis reduction method]]). In the case of the Euclidean norm, the final search can be shortened by an idea of U. Finke and M. Pohst [[#References|[a8]]].
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006049.png" /></td> </tr></table>
+
For any sequence $\{ u _ { i } \}$ of $[ 0,1 )$-uniformly distributed random numbers, the local deviations
 +
 
 +
\begin{equation*} \Delta _ { k } ( \mathbf{s} , \mathbf{t} ) = - \prod _ { j = 1 } ^ { k } ( t _ { j } - s _ { j } ) + \end{equation*}
 +
 
 +
\begin{equation*} + \frac {\#  \{ U _ { i } = ( u _ { i + 1} , \ldots , u _ { i  + k} ) : s _ { j } &lt; u_{i + j} \leq t _ { j } , 1 \leq j \leq k \} } { \# \{ U _ { i } = ( u _ { i  + 1} , \ldots , u _ { i  + k}) \} } \end{equation*}
  
 
and their largest value, the (global) discrepancy
 
and their largest value, the (global) discrepancy
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006050.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a2)</td></tr></table>
+
\begin{equation} \tag{a2} \Delta _ { k } = \operatorname { sup } \{ | \Delta _ { k } ( \mathbf{s} , \mathbf{t} ) | : 0 \leq s _ { j } \leq t _ { j } &lt; 1,1 \leq j \leq k \}, \end{equation}
  
are of great importance. For example, for the calculation of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006051.png" />-dimensional integrals by Monte-Carlo methods, the difference of the integral and its approximation by a Riemann sum is bounded by the discrepancy <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006052.png" /> multiplied by the variation of the function <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006053.png" /> (in the sense of Hardy–Krause, cf. also [[Variation of a function|Variation of a function]]). Since the variation of the function is fixed, the discrepancy has to be as small as possible. See [[#References|[a9]]].
+
are of great importance. For example, for the calculation of $k$-dimensional integrals by Monte-Carlo methods, the difference of the integral and its approximation by a Riemann sum is bounded by the discrepancy $\Delta _ { k }$ multiplied by the variation of the function $V ( f )$ (in the sense of Hardy–Krause, cf. also [[Variation of a function|Variation of a function]]). Since the variation of the function is fixed, the discrepancy has to be as small as possible. See [[#References|[a9]]].
  
 
No methods are known for calculating the discrepancy in dimension greater than two. Dieter derived a lower bound in 1973, and H. Niederreiter found an upper bound in 1978 [[#References|[a13]]]. Even these bounds are difficult to calculate.
 
No methods are known for calculating the discrepancy in dimension greater than two. Dieter derived a lower bound in 1973, and H. Niederreiter found an upper bound in 1978 [[#References|[a13]]]. Even these bounds are difficult to calculate.
  
In dimension two, i.e. in the case of pairs, all three quantities can be calculated by the [[Euclidean algorithm|Euclidean algorithm]] for the period length <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006054.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006055.png" />. <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006056.png" /> is equal to <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006057.png" /> if <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006058.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006059.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006060.png" /> in the two other cases. Define a sequence <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006061.png" /> by
+
In dimension two, i.e. in the case of pairs, all three quantities can be calculated by the [[Euclidean algorithm|Euclidean algorithm]] for the period length $n$ and $a$. $n$ is equal to $m /4$ if $m = 2 ^ { E }$ and $r = 0$ and $n = m$ in the two other cases. Define a sequence $\{ m_i \}$ by
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006062.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a3)</td></tr></table>
+
\begin{equation} \tag{a3} m_0 = n , m_1 = a, \end{equation}
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006063.png" /></td> </tr></table>
+
\begin{equation*} m _ { i - 1 } = a _ { i - 1 } m _ { i } + m _ { i + 1 } , i = 1,2 , \dots , \end{equation*}
  
 
where
 
where
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006064.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a4)</td></tr></table>
+
\begin{equation} \tag{a4} a _ { i - 1 } = \lfloor \frac { m _ { i - 1 } } { m _ { i } } \rfloor , \end{equation}
  
and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006065.png" /> is the integer function (or [[Floor function|floor function]]). Associated is the sequence
+
and $\lfloor x \rfloor$ is the integer function (or [[Floor function|floor function]]). Associated is the sequence
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006066.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a5)</td></tr></table>
+
\begin{equation} \tag{a5} p _ { 0 } = 0 , p _ { 1 } = 1, \end{equation}
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006067.png" /></td> </tr></table>
+
\begin{equation*} p _ { i + 1 } = a _ { i - 1 } p _ { i } + p _ { i - 1 } ,\,  i = 1,2, \dots . \end{equation*}
  
 
Then
 
Then
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006068.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a6)</td></tr></table>
+
\begin{equation} \tag{a6} N _ { 2 } ^ { * } = \operatorname { min } _ { i } \{ m _ { i } + p _ { i } \} \end{equation}
  
 
and
 
and
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006069.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a7)</td></tr></table>
+
\begin{equation} \tag{a7} W _ { 2 } ^ { * } = \frac { 1 } { D _ { 2 } ^ { * } } = \operatorname { min } _ { i } \sqrt { m _ { i } ^ { 2 } + p _ { i } ^ { 2 } }. \end{equation}
  
 
Finally, for the discrepancy the following rather sharp bounds hold [[#References|[a6]]]:
 
Finally, for the discrepancy the following rather sharp bounds hold [[#References|[a6]]]:
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006070.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a8)</td></tr></table>
+
\begin{equation} \tag{a8} \frac { 1 } { 4 n } \operatorname { max } \{ a _ { i } : 0 \leq i \leq t \} \leq \Delta _ { 2 } \leq \frac { 1 } { 4 n } \left( \sum _ { i = 0 } ^ { t } a _ { i } + 2 \right). \end{equation}
  
 
See [[#References|[a2]]], [[#References|[a7]]] for methods for calculating exact values of the discrepancy. The numerical values differ only slightly from the upper bound in (a8).
 
See [[#References|[a2]]], [[#References|[a7]]] for methods for calculating exact values of the discrepancy. The numerical values differ only slightly from the upper bound in (a8).
  
If the dimensions become larger, the number of hyperplanes on which the lattice points <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006071.png" /> lie decreases considerably. Therefore a different procedure was proposed by Knuth [[#References|[a10]]]: A sequence of integers <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006072.png" /> is initialized to <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006073.png" /> and updated by
+
If the dimensions become larger, the number of hyperplanes on which the lattice points $P _ { k }$ lie decreases considerably. Therefore a different procedure was proposed by Knuth [[#References|[a10]]]: A sequence of integers $z_i$ is initialized to $( z_0 , \dots , z _ { r  - 1} ) \neq ( 0 , \dots , 0 )$ and updated by
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006074.png" /></td> <td valign="top" style="width:5%;text-align:right;">(a9)</td></tr></table>
+
\begin{equation} \tag{a9} z _ { i } \equiv a _ { i } z _ { i - 1 } + \ldots + a _ { i } z _ { i - r } ( \operatorname { mod } p ) \end{equation}
  
for <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006075.png" />. Here the factors <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006076.png" /> are given integers, and for the modulus <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006077.png" /> only prime numbers are considered. Again, the fractions <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006078.png" /> are taken as random samples from the interval <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006079.png" />.
+
for $0 \leq z _ { i } &lt; p$. Here the factors $a_i$ are given integers, and for the modulus $p$ only prime numbers are considered. Again, the fractions $u _ { i } = z _ { i } / p$ are taken as random samples from the interval $[ 0,1 )$.
  
Since there are only <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006080.png" /> possible <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006081.png" />-tuplets <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006082.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006083.png" /> must not occur, the period length of (a9) is at most <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006084.png" />. It can be shown that this maximum period length of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006085.png" /> may in fact be achieved for suitable choices of the factors <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006086.png" />. This means that all <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006087.png" />-tuplets <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006088.png" /> must occur, resulting in perfectly uniform distributions of the <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006089.png" />, the pairs <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006090.png" />, the triplets <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006091.png" /> (if <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006092.png" />), etc. In other words, for all dimensions <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006093.png" /> the hypercube <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006094.png" /> is now evenly filled. For dimensions <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006095.png" /> the quantities <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006096.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006097.png" /> can be calculated exactly, since the points <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006098.png" /> are again points of a lattice. The numerical values show that a sequence of type (a9) behaves as good in dimension <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l13006099.png" /> as a linear congruential generator behaves in dimension <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l130060100.png" />. Furthermore, reduced bases (in the sense of H. Minkowski) can be determined which show how  "good"  the specific generator behaves. See [[#References|[a3]]] and recent publications of L. Afflerbach and H. Grothe, starting with [[#References|[a1]]].
+
Since there are only $p ^ { r}$ possible $r$-tuplets $( z _ { k } , \ldots , z _ { k  + r - 1})$ and $( 0 , \ldots , 0 )$ must not occur, the period length of (a9) is at most $p ^ { r } - 1$. It can be shown that this maximum period length of $p ^ { r } - 1$ may in fact be achieved for suitable choices of the factors $a _ { 1 } , \dots , a _ { r }$. This means that all $r$-tuplets $( z _ { k } , \ldots , z _ { k + r - 1} ) \neq ( 0 , \dots , 0 )$ must occur, resulting in perfectly uniform distributions of the $u _ { i } = z _ { i } / p$, the pairs $( u _ { i } , u _ { i  + 1} )$, the triplets $( u _ { i } , u _ { i + 1} , u _ { i + 2} )$ (if $r \geq 3$), etc. In other words, for all dimensions $k \leq r$ the hypercube $[ 0,1 ) ^ { k }$ is now evenly filled. For dimensions $k &gt; r$ the quantities $N _ { k } ^ { * }$ and $D _ { k } ^ { * }$ can be calculated exactly, since the points $P _ { k }$ are again points of a lattice. The numerical values show that a sequence of type (a9) behaves as good in dimension $k \times r$ as a linear congruential generator behaves in dimension $k$. Furthermore, reduced bases (in the sense of H. Minkowski) can be determined which show how  "good"  the specific generator behaves. See [[#References|[a3]]] and recent publications of L. Afflerbach and H. Grothe, starting with [[#References|[a1]]].
  
The case <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l130060101.png" /> is of special interest; it is called the Tausworth generator. Here, bits are produced and a word of, say, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/l/l130/l130060/l130060102.png" /> bits is taken as a sample from the uniform distribution. See [[#References|[a11]]] for further information.
+
The case $p = 2$ is of special interest; it is called the Tausworth generator. Here, bits are produced and a word of, say, $8$ bits is taken as a sample from the uniform distribution. See [[#References|[a11]]] for further information.
  
 
====References====
 
====References====
<table><TR><TD valign="top">[a1]</TD> <TD valign="top">  L. Afflerbach,  H. Grothe,  "Calculation of Minkowski-reduced lattice bases"  ''Computing'' , '''35'''  (1985)  pp. 269–276</TD></TR><TR><TD valign="top">[a2]</TD> <TD valign="top">  L. Afflerbach,  R. Weilbächer,  "The exact determination of rectangle discrepancy for linear congruential pseudorandom generators"  ''Math. Comput.'' , '''53'''  (1989)  pp. 343–354</TD></TR><TR><TD valign="top">[a3]</TD> <TD valign="top">  W.A. Beyer,  R.B. Roof,  D. Williamson,  "The lattice structure of multiplicative congruential pseudo-random vectors"  ''Math. Comput.'' , '''25'''  (1971)  pp. 345–360</TD></TR><TR><TD valign="top">[a4]</TD> <TD valign="top">  R.R. Coveyou,  R.D. MacPherson,  "Fourier analysis of uniform random number generators"  ''J. ACM'' , '''14'''  (1967)  pp. 100–119</TD></TR><TR><TD valign="top">[a5]</TD> <TD valign="top">  U. Dieter,  "Pseudo-random numbers: the exact distribution of pairs"  ''Math. Comput.'' , '''25'''  (1971)  pp. 855–883</TD></TR><TR><TD valign="top">[a6]</TD> <TD valign="top">  U. Dieter,  "How to calculate shortest vectors in a lattice"  ''Math. Comput.'' , '''29'''  (1975)  pp. 827–833</TD></TR><TR><TD valign="top">[a7]</TD> <TD valign="top">  U. Dieter,  "Pseudo-random numbers: The discrepancy in two dimensions"  ''to appear''  (2002)</TD></TR><TR><TD valign="top">[a8]</TD> <TD valign="top">  U. Finke,  M. Pohst,  "Improved methods for calculating vectors of short length in a lattice, including a complexity analysis"  ''Math. Comput.'' , '''44'''  (1985)  pp. 463–471</TD></TR><TR><TD valign="top">[a9]</TD> <TD valign="top">  G.S. Fishman,  "Monte Carlo: Concepts, algorithms, and applications" , Springer  (1996)</TD></TR><TR><TD valign="top">[a10]</TD> <TD valign="top">  D.E. Knuth,  "The art of computer programming" , '''II: Seminumerical algorithms''' , Addison-Wesley  (1969)  (Edition: First)</TD></TR><TR><TD valign="top">[a11]</TD> <TD valign="top">  T.G. Lewis,  W.H. Payne,  "Generalized feedback shift register pseudorandom number algorithms"  ''J. ACM'' , '''20'''  (1973)  pp. 456–468</TD></TR><TR><TD valign="top">[a12]</TD> <TD valign="top">  G. Marsaglia,  "Random numbers fall mainly in the planes"  ''Proc. Nat. Acad. Sci.'' , '''61'''  (1968)  pp. 25–28</TD></TR><TR><TD valign="top">[a13]</TD> <TD valign="top">  H. Niederreiter,  "Quasi-Monte-Carlo methods and pseudo-random numbers"  ''Bull. Amer. Math. Soc.'' , '''84'''  (1978)  pp. 957–1041</TD></TR></table>
+
<table><tr><td valign="top">[a1]</td> <td valign="top">  L. Afflerbach,  H. Grothe,  "Calculation of Minkowski-reduced lattice bases"  ''Computing'' , '''35'''  (1985)  pp. 269–276</td></tr><tr><td valign="top">[a2]</td> <td valign="top">  L. Afflerbach,  R. Weilbächer,  "The exact determination of rectangle discrepancy for linear congruential pseudorandom generators"  ''Math. Comput.'' , '''53'''  (1989)  pp. 343–354</td></tr><tr><td valign="top">[a3]</td> <td valign="top">  W.A. Beyer,  R.B. Roof,  D. Williamson,  "The lattice structure of multiplicative congruential pseudo-random vectors"  ''Math. Comput.'' , '''25'''  (1971)  pp. 345–360</td></tr><tr><td valign="top">[a4]</td> <td valign="top">  R.R. Coveyou,  R.D. MacPherson,  "Fourier analysis of uniform random number generators"  ''J. ACM'' , '''14'''  (1967)  pp. 100–119</td></tr><tr><td valign="top">[a5]</td> <td valign="top">  U. Dieter,  "Pseudo-random numbers: the exact distribution of pairs"  ''Math. Comput.'' , '''25'''  (1971)  pp. 855–883</td></tr><tr><td valign="top">[a6]</td> <td valign="top">  U. Dieter,  "How to calculate shortest vectors in a lattice"  ''Math. Comput.'' , '''29'''  (1975)  pp. 827–833</td></tr><tr><td valign="top">[a7]</td> <td valign="top">  U. Dieter,  "Pseudo-random numbers: The discrepancy in two dimensions"  ''to appear''  (2002)</td></tr><tr><td valign="top">[a8]</td> <td valign="top">  U. Finke,  M. Pohst,  "Improved methods for calculating vectors of short length in a lattice, including a complexity analysis"  ''Math. Comput.'' , '''44'''  (1985)  pp. 463–471</td></tr><tr><td valign="top">[a9]</td> <td valign="top">  G.S. Fishman,  "Monte Carlo: Concepts, algorithms, and applications" , Springer  (1996)</td></tr><tr><td valign="top">[a10]</td> <td valign="top">  D.E. Knuth,  "The art of computer programming" , '''II: Seminumerical algorithms''' , Addison-Wesley  (1969)  (Edition: First)</td></tr><tr><td valign="top">[a11]</td> <td valign="top">  T.G. Lewis,  W.H. Payne,  "Generalized feedback shift register pseudorandom number algorithms"  ''J. ACM'' , '''20'''  (1973)  pp. 456–468</td></tr><tr><td valign="top">[a12]</td> <td valign="top">  G. Marsaglia,  "Random numbers fall mainly in the planes"  ''Proc. Nat. Acad. Sci.'' , '''61'''  (1968)  pp. 25–28</td></tr><tr><td valign="top">[a13]</td> <td valign="top">  H. Niederreiter,  "Quasi-Monte-Carlo methods and pseudo-random numbers"  ''Bull. Amer. Math. Soc.'' , '''84'''  (1978)  pp. 957–1041</td></tr></table>

Revision as of 17:00, 1 July 2020

A method widely used for generating random numbers from the uniform distribution: A sequence of integers is initialized with a value $z_0$ and continued as

\begin{equation} \tag{a1} z _ { i + 1} \equiv a z _ { i } + r ( \operatorname { mod } m ) ,\, 0 \leq z _ { i } < m, \end{equation}

for all $i$. The fractions $u _ { i } = z _ { i } / m$ are the derived pseudo-random numbers in the interval $[ 0,1 )$ (cf. also Random and pseudo-random numbers; Pseudo-random numbers). The constants $m$, the modulus, $a$, the multiplicator, $r$, the increment, and $z_0$, the starting number, are suitably chosen non-negative integers. Three choices of $m$, $a$ and $r$ are common on most computers:

1) $r = 0$, $m = 2 ^ { E }$, $a \equiv 5 ( \operatorname { mod } 8 )$, and $z _ { 0 } \equiv 1 ( \operatorname { mod } 4 )$. All $z _ { i } \equiv 1 ( \operatorname { mod } 4 )$ are generated.

2) $r = 0$, $m = p$, $p$ prime, $a$ a primitive root modulo $p$. All $z _ { i } = 1 , \dots , p - 1$ are generated.

3) $r \equiv 1 ( \operatorname { mod } 2 )$, $m = 2 ^ { E }$, $a \equiv 1 ( \operatorname { mod } 4 )$. All integers $0 , \ldots , 2 ^ { E } - 1$ are generated.

For selecting "good" random number generators one has to study the distribution of the $k$-tuplets $P _ { k } = ( u _ { i + 1} , \dots , u _ { i + k})$. Geometrically these $P _ { k }$ may be considered as points of a lattice $G$ in the $k$-dimensional hypercube $[ 0,1 ) ^ { k }$. The lattice points can also be seen as intersection points of $k$ sets of parallel hyperplanes. Consequently, the following questions may be raised:

i) Determine the minimal number $N _ { k } ^ { * }$ of parallel hyperplanes on which all points $P _ { k }$ lie.

ii) Determine the maximal distance $D _ { k } ^ { * }$ of parallel hyperplanes on which all points $P _ { k }$ lie.

Question i) was asked by G. Marsaglia [a12], who derived upper bounds for $N _ { k } ^ { * }$ using Minkowski's convex body theorem (cf. also Minkowski theorem). The "wave numbers" $W _ { k } ^ { * } = 1 / D _ { k } ^ { * }$ were introduced by R.R. Coveyou and R.D. MacPherson [a4]. Their algorithm for calculating $W _ { k } ^ { * }$ was simplified by D.E. Knuth [a10].

The calculation of both quantities is based on a general procedure to determine non-zero vectors of shortest length in the dual lattice of covering hyperplanes. For the determination of $N _ { k } ^ { * }$ the $\mathbf{l}_{1}$-norm is used, and for $D _ { k } ^ { * }$ the Euclidean norm is appropriate. The algorithm of U. Dieter (1973) gives exact values for both quantities; no exact values for $N _ { k } ^ { * }$ were known before. Knuth included a variant of this algorithm in [a10]. A completely different approach was proposed by L. Lovász, J.K. Lenstra and H.W. Lenstra, Jr., called the LLL-algorithm (cf. also LLL basis reduction method). In the case of the Euclidean norm, the final search can be shortened by an idea of U. Finke and M. Pohst [a8].

For any sequence $\{ u _ { i } \}$ of $[ 0,1 )$-uniformly distributed random numbers, the local deviations

\begin{equation*} \Delta _ { k } ( \mathbf{s} , \mathbf{t} ) = - \prod _ { j = 1 } ^ { k } ( t _ { j } - s _ { j } ) + \end{equation*}

\begin{equation*} + \frac {\# \{ U _ { i } = ( u _ { i + 1} , \ldots , u _ { i + k} ) : s _ { j } < u_{i + j} \leq t _ { j } , 1 \leq j \leq k \} } { \# \{ U _ { i } = ( u _ { i + 1} , \ldots , u _ { i + k}) \} } \end{equation*}

and their largest value, the (global) discrepancy

\begin{equation} \tag{a2} \Delta _ { k } = \operatorname { sup } \{ | \Delta _ { k } ( \mathbf{s} , \mathbf{t} ) | : 0 \leq s _ { j } \leq t _ { j } < 1,1 \leq j \leq k \}, \end{equation}

are of great importance. For example, for the calculation of $k$-dimensional integrals by Monte-Carlo methods, the difference of the integral and its approximation by a Riemann sum is bounded by the discrepancy $\Delta _ { k }$ multiplied by the variation of the function $V ( f )$ (in the sense of Hardy–Krause, cf. also Variation of a function). Since the variation of the function is fixed, the discrepancy has to be as small as possible. See [a9].

No methods are known for calculating the discrepancy in dimension greater than two. Dieter derived a lower bound in 1973, and H. Niederreiter found an upper bound in 1978 [a13]. Even these bounds are difficult to calculate.

In dimension two, i.e. in the case of pairs, all three quantities can be calculated by the Euclidean algorithm for the period length $n$ and $a$. $n$ is equal to $m /4$ if $m = 2 ^ { E }$ and $r = 0$ and $n = m$ in the two other cases. Define a sequence $\{ m_i \}$ by

\begin{equation} \tag{a3} m_0 = n , m_1 = a, \end{equation}

\begin{equation*} m _ { i - 1 } = a _ { i - 1 } m _ { i } + m _ { i + 1 } , i = 1,2 , \dots , \end{equation*}

where

\begin{equation} \tag{a4} a _ { i - 1 } = \lfloor \frac { m _ { i - 1 } } { m _ { i } } \rfloor , \end{equation}

and $\lfloor x \rfloor$ is the integer function (or floor function). Associated is the sequence

\begin{equation} \tag{a5} p _ { 0 } = 0 , p _ { 1 } = 1, \end{equation}

\begin{equation*} p _ { i + 1 } = a _ { i - 1 } p _ { i } + p _ { i - 1 } ,\, i = 1,2, \dots . \end{equation*}

Then

\begin{equation} \tag{a6} N _ { 2 } ^ { * } = \operatorname { min } _ { i } \{ m _ { i } + p _ { i } \} \end{equation}

and

\begin{equation} \tag{a7} W _ { 2 } ^ { * } = \frac { 1 } { D _ { 2 } ^ { * } } = \operatorname { min } _ { i } \sqrt { m _ { i } ^ { 2 } + p _ { i } ^ { 2 } }. \end{equation}

Finally, for the discrepancy the following rather sharp bounds hold [a6]:

\begin{equation} \tag{a8} \frac { 1 } { 4 n } \operatorname { max } \{ a _ { i } : 0 \leq i \leq t \} \leq \Delta _ { 2 } \leq \frac { 1 } { 4 n } \left( \sum _ { i = 0 } ^ { t } a _ { i } + 2 \right). \end{equation}

See [a2], [a7] for methods for calculating exact values of the discrepancy. The numerical values differ only slightly from the upper bound in (a8).

If the dimensions become larger, the number of hyperplanes on which the lattice points $P _ { k }$ lie decreases considerably. Therefore a different procedure was proposed by Knuth [a10]: A sequence of integers $z_i$ is initialized to $( z_0 , \dots , z _ { r - 1} ) \neq ( 0 , \dots , 0 )$ and updated by

\begin{equation} \tag{a9} z _ { i } \equiv a _ { i } z _ { i - 1 } + \ldots + a _ { i } z _ { i - r } ( \operatorname { mod } p ) \end{equation}

for $0 \leq z _ { i } < p$. Here the factors $a_i$ are given integers, and for the modulus $p$ only prime numbers are considered. Again, the fractions $u _ { i } = z _ { i } / p$ are taken as random samples from the interval $[ 0,1 )$.

Since there are only $p ^ { r}$ possible $r$-tuplets $( z _ { k } , \ldots , z _ { k + r - 1})$ and $( 0 , \ldots , 0 )$ must not occur, the period length of (a9) is at most $p ^ { r } - 1$. It can be shown that this maximum period length of $p ^ { r } - 1$ may in fact be achieved for suitable choices of the factors $a _ { 1 } , \dots , a _ { r }$. This means that all $r$-tuplets $( z _ { k } , \ldots , z _ { k + r - 1} ) \neq ( 0 , \dots , 0 )$ must occur, resulting in perfectly uniform distributions of the $u _ { i } = z _ { i } / p$, the pairs $( u _ { i } , u _ { i + 1} )$, the triplets $( u _ { i } , u _ { i + 1} , u _ { i + 2} )$ (if $r \geq 3$), etc. In other words, for all dimensions $k \leq r$ the hypercube $[ 0,1 ) ^ { k }$ is now evenly filled. For dimensions $k > r$ the quantities $N _ { k } ^ { * }$ and $D _ { k } ^ { * }$ can be calculated exactly, since the points $P _ { k }$ are again points of a lattice. The numerical values show that a sequence of type (a9) behaves as good in dimension $k \times r$ as a linear congruential generator behaves in dimension $k$. Furthermore, reduced bases (in the sense of H. Minkowski) can be determined which show how "good" the specific generator behaves. See [a3] and recent publications of L. Afflerbach and H. Grothe, starting with [a1].

The case $p = 2$ is of special interest; it is called the Tausworth generator. Here, bits are produced and a word of, say, $8$ bits is taken as a sample from the uniform distribution. See [a11] for further information.

References

[a1] L. Afflerbach, H. Grothe, "Calculation of Minkowski-reduced lattice bases" Computing , 35 (1985) pp. 269–276
[a2] L. Afflerbach, R. Weilbächer, "The exact determination of rectangle discrepancy for linear congruential pseudorandom generators" Math. Comput. , 53 (1989) pp. 343–354
[a3] W.A. Beyer, R.B. Roof, D. Williamson, "The lattice structure of multiplicative congruential pseudo-random vectors" Math. Comput. , 25 (1971) pp. 345–360
[a4] R.R. Coveyou, R.D. MacPherson, "Fourier analysis of uniform random number generators" J. ACM , 14 (1967) pp. 100–119
[a5] U. Dieter, "Pseudo-random numbers: the exact distribution of pairs" Math. Comput. , 25 (1971) pp. 855–883
[a6] U. Dieter, "How to calculate shortest vectors in a lattice" Math. Comput. , 29 (1975) pp. 827–833
[a7] U. Dieter, "Pseudo-random numbers: The discrepancy in two dimensions" to appear (2002)
[a8] U. Finke, M. Pohst, "Improved methods for calculating vectors of short length in a lattice, including a complexity analysis" Math. Comput. , 44 (1985) pp. 463–471
[a9] G.S. Fishman, "Monte Carlo: Concepts, algorithms, and applications" , Springer (1996)
[a10] D.E. Knuth, "The art of computer programming" , II: Seminumerical algorithms , Addison-Wesley (1969) (Edition: First)
[a11] T.G. Lewis, W.H. Payne, "Generalized feedback shift register pseudorandom number algorithms" J. ACM , 20 (1973) pp. 456–468
[a12] G. Marsaglia, "Random numbers fall mainly in the planes" Proc. Nat. Acad. Sci. , 61 (1968) pp. 25–28
[a13] H. Niederreiter, "Quasi-Monte-Carlo methods and pseudo-random numbers" Bull. Amer. Math. Soc. , 84 (1978) pp. 957–1041
How to Cite This Entry:
Linear congruential method. Encyclopedia of Mathematics. URL: http://encyclopediaofmath.org/index.php?title=Linear_congruential_method&oldid=14102
This article was adapted from an original article by U. Dieter (originator), which appeared in Encyclopedia of Mathematics - ISBN 1402006098. See original article