Linear congruential method

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

(a1)

for all . The fractions are the derived pseudo-random numbers in the interval (cf. also Random and pseudo-random numbers; Pseudo-random numbers). The constants , the modulus, , the multiplicator, , the increment, and , the starting number, are suitably chosen non-negative integers. Three choices of , and are common on most computers:

1) , , , and . All are generated.

2) , , prime, a primitive root modulo . All are generated.

3) , , . All integers are generated.

For selecting "good" random number generators one has to study the distribution of the -tuplets . Geometrically these may be considered as points of a lattice in the -dimensional hypercube . The lattice points can also be seen as intersection points of sets of parallel hyperplanes. Consequently, the following questions may be raised:

i) Determine the minimal number of parallel hyperplanes on which all points lie.

ii) Determine the maximal distance of parallel hyperplanes on which all points lie.

Question i) was asked by G. Marsaglia [a12], who derived upper bounds for using Minkowski's convex body theorem (cf. also Minkowski theorem). The "wave numbers" were introduced by R.R. Coveyou and R.D. MacPherson [a4]. Their algorithm for calculating 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 the -norm is used, and for the Euclidean norm is appropriate. The algorithm of U. Dieter (1973) gives exact values for both quantities; no exact values for 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 of -uniformly distributed random numbers, the local deviations

and their largest value, the (global) discrepancy

(a2)

are of great importance. For example, for the calculation of -dimensional integrals by Monte-Carlo methods, the difference of the integral and its approximation by a Riemann sum is bounded by the discrepancy multiplied by the variation of the function (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 and . is equal to if and and in the two other cases. Define a sequence by

(a3)

where

(a4)

and is the integer function (or floor function). Associated is the sequence

(a5)

Then

(a6)

and

(a7)

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

(a8)

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 lie decreases considerably. Therefore a different procedure was proposed by Knuth [a10]: A sequence of integers is initialized to and updated by

(a9)

for . Here the factors are given integers, and for the modulus only prime numbers are considered. Again, the fractions are taken as random samples from the interval .

Since there are only possible -tuplets and must not occur, the period length of (a9) is at most . It can be shown that this maximum period length of may in fact be achieved for suitable choices of the factors . This means that all -tuplets must occur, resulting in perfectly uniform distributions of the , the pairs , the triplets (if ), etc. In other words, for all dimensions the hypercube is now evenly filled. For dimensions the quantities and can be calculated exactly, since the points are again points of a lattice. The numerical values show that a sequence of type (a9) behaves as good in dimension as a linear congruential generator behaves in dimension . 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 is of special interest; it is called the Tausworth generator. Here, bits are produced and a word of, say, 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