Linear congruential method

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