# Multiplicative, congruential random-number generators with multiplier +/- 2(super k1) +/- 2(super k2) and modulus 2(super p) - 1.

1. INTRODUCTIONThe demand for random numbers in many scientific applications is increasing due to the use of high-performance computer systems for large-scale scientific problems. Even a low-end workstation can consume [10.sup.7] random numbers in a few minutes. However, the most widely used multiplicative, congruential random-number generators with modulus [2.sup.31] - 1 have a cycle length of about 2.1 x [10.sup.9]. Kirkpatrick and Stoll [1981] presented a lagged-Fibonacci generator (cf. Anderson [1990]), called R250, which is very fast and has a cycle length of [2.sup.250] - 1. Although R250 passed a number of statistical tests, a recent experience [Selke 1993] showed that R250 gives wrong results for a clustered Monte Carlo simulation, while the multiplicative congruential generator

[x.sub.i] = 16807 [multiplied by] [x.sub.i-1] mod ([2.sup.31] - 1),

the minimal standard generator suggested by Park and Miller [1988], worked very well.

Multiplicative congruential generators have been widely used and tested. It would be better if this kind of generator provides a longer cycle length to meet the current needs of large-scale simulations. Payne et al. [1969] predicted that, due to increases in computer speed and the next Mersenne prime of [2.sup.31] - 1 being [2.sup.61] - 1, multiplicative congruential generators with modulus [2.sup.61] - 1 would be needed. The time has now arrived. The cycle length of these generators, up to [2.sup.61] - 2 [approximately equal to] 2.3 x [10.sup.18], is long enough for most current scientific applications. In the future, generators with modulus [2.sup.127] - 1 may be required to provide a longer cycle length of up to 1.7 x [10.sup.38]. However, developing portable and efficient multiplicative congruential generators is not straightforward [Park and Miller 1988; Schrage 1979]. Moreover, developing generators with larger moduli such as [2.sup.61] - 1 is more difficult than developing those with modulus [2.sup.31] - 1.

This article presents the development of multiplicative congruential generators with prime modulus m = [2.sup.p] - 1 and four forms of multipliers: [2.sup.k1] - [2.sup.k2], [2.sup.k1] + [2.sup.k2], m - [2.sup.k2] + [2.sup.k2], and m - [2.sup.k1] - [2.sup.k2], k1 [greater than] k2. The multipliers for modulus [2.sup.31] - 1 and [2.sup.61] - 1 are measured by spectral tests, and the best ones are presented. The generators with these multipliers are portable and very fast. The performance results compared with other implementation techniques are also presented. These generators have also passed several empirical tests, including the frequency test, the run test, and the maximum-of-t test [Knuth 1981].

2. RELATED WORK

A linear congruential generator can be defined by the following equation:

[x.sub.i] = a [multiplied by] [x.sub.i-1] + c mod m. (1)

Computing the above equation requires one integer addition, one multiplication, and one division. Usually m = [2.sup.e] is chosen to avoid the division operation and the portability problem in getting product a [multiplied by] [x.sub.i-1] in high-level languages. Knuth [1981, pp. 22-24] analyzed multiplier a = [2.sup.k] + 1 for generators with modulus [2.sup.e], k [less than] e, and c = 1. The equation is shown below:

[x.sub.i] = (([2.sup.k] + 1) [multiplied by] [x.sub.i-1]) + 1 mod [2.sup.e] (2)

The above equation can be computed by merely shifting and adding, and it avoids the multiplication and division operation. However, Knuth [1981, p. 22] showed that this kind of multiplier should be avoided based on the concept of potency. Knuth [1981, p. 24] suggested that (for multiplier a = [2.sup.k] + 1) multiplicative congruential generators with a prime modulus should be used instead.

A multiplicative congruential generator is a special case of Eq. (1) taking c = 0:

[x.sub.i] = a [multiplied by] [x.sub.i-1] mod m. (3)

To get a period of maximal length m - 1, m must be a prime; a is a primitive root of m; and [x.sub.0] [Epsilon][1, m - 1] (cf. Knuth [1981, p.19]). Developing portable and efficient multiplicative congruential generators is not straightforward, because simply taking m = [2.sup.e] does not get an "almost correct" implementation. Park and Miller [1988] gave a sampling of simple but "bad" multiplicative congruential generators, many of which use modulus [2.sup.e], e.g., e = 16, e = 31, and e = 32.

Some research has been devoted to portable implementations of multiplicative congruential generators with m = [2.sup.31] - 1. Schrage [1979] presented a Fortran program that uses two integer multiplication operations and several bit operations to get the product 16807 [multiplied by] x. This technique is more costly to apply for generators with a larger modulus. Park and Miller [1988] presented a portable program based on simple integer arithmetic. The technique also works for generators with modulus [2.sup.61] - 1, if the compiler used supports 64-bit integer arithmetic. All these portable implementations allow only small multipliers and are slower than those using assembly codes directly.

3. THE [+ or -] [2.sup.k1] [+ or -] [2.sup.k2]. MULTIPLIERS

The division operation (rood m) in multiplicative congruential generators with (prime) modulus m = [2.sup.p] - 1 can be performed by shifting and addition [Payne et al. 1969]. To further replace multiplication with shifting and addition, multiplier a must be in the form of simple expressions of [2.sup.k]. To guarantee that the resulting generators are good, we need to meet the following criteria: A multiplier a is adequate if

(1) The number a is a primitive root of m (or primitive element modulo m): a is a primitive root of m if [a.sup.n] mod m [not equal to] 1 for n = 1, 2, ..., m - 2 [Knuth 1981, p. 10]. Knuth [1981, p. 20] presented a method to test whether a number is a primitive root of m:

The number a is a primitive root of m if and only if

a [not equal to] 0(mod m), and [a.sup.(m-1/q)] [not equal to] 1(mod m),

for any prime divisor q of m - 1.

(2) The generator with multiplier a and modulus m should pass a number of statistical tests. Knuth [1981, p. 89] and Fishman and Moore [1986] all suggested the spectral test.

The multipliers with the simplest form are [+ or -] [2.sup.k] mod m, i.e., multipliers [2.sup.k] and m - [2.sup.k], [Mathematical Expression Omitted], where lgm denotes [log.sub.2] m. Arithmetic modulo [2.sup.p] - 1, for any p [greater than] 1, is well known as "one's complement" arithmetic. For this arithmetic, multiplication by [2.sup.k] results in a k-place left cyclic shift of the binary digits. That is, if integer x has binary digits

[b.sub.p-1] . . . [b.sub.0]

then the integer [2.sup.k] [multiplied by] x mod ([2.sup.p] - 1) has digits

[b.sub.p-k-1] . . . [b.sub.0] [b.sub.p-1] . . . [b.sub.p-k].

Also, -x mod ([2.sup.p] - 1) = [2.sup.p] - 1 - x results in bitwise negation. Thus, with m = [2.sup.p] - 1, a multiplier of form a = [2.sup.k] must cause the sequence to cycle after p steps, as p cyclic k-place shift restores the original digits. If a = -[2.sup.k], each step is a k-place cyclic shift followed by negation, and the sequence will cycle after p steps if p is even, or 2p steps otherwise. Hence, no such multiplier can be a primitive root.

Now, we pay attention to multipliers with the form of [+ or -] [2.sup.k1] [+ or -] [2.sup.k2] mod m, k1 [greater than] k2, i.e., four kinds of multipliers: [2.sup.k1] - [2.sup.k2], [2.sup.k1] + [2.sup.k2], m - [2.sup.k1] + [2.sup.k2], and m - [2.sup.k1] - [2.sup.k2]. There is an efficient algorithm for generators with these multipliers.

Let a = [+ or -][2.sup.k1] [+ or -] [2.sup.k2], m = [2.sup.p] - 1.

[x.sub.i] = a [multiplied by] [x.sub.i-1] mod m.

= [x.sub.i-1]([+ or -] [2.sup.k1] [+ or -] [2.sup.k2]) mod m.

= ([+ or -] [2.sup.k1] [multiplied by] [x.sub.i-1] [+ or -] [2.sup.k2] [multiplied by] [x.sub.i-1]) mod m.

= ([+ or -] ([2.sup.k1] [multiplied by] [x.sub.i-1] mod m) [+ or -] ([2.sup.k2] [multiplied by] [x.sub.i-1] mod m)) mod m.

Multiplication by [2.sup.k] is a k-place left cyclic shift and can be rewritten as follows:

Let x be a p-bit number, k [less than] p.

[2.sup.k] [multiplied by] x mod ([2.sup.p] - 1) = x (high-order k bits) + x (low-order p - k bits) [multiplied by] [2.sup.k].

Thus, the multiplication and division can be reduced to shifting and addition [ILLUSTRATION FOR FIGURE 1 OMITTED].

Let [w.sub.1] = [2.sup.k1] [multiplied by] [x.sub.i-1] mod m, [w.sub.2] = [2.sup.k2] [multiplied by] [x.sub.i-1] mod m. Table I lists the equations for all these multipliers. Because 0 [less than] [w.sub.1], [w.sub.2] [less than] m, the range of x[prime] is - m [less than] x[prime] [less than] m. When x[prime] [greater than] 0, [x.sub.i] = x[prime]; otherwise, [x.sub.i] = m + x[prime]. Algorithm RAND shows the code.

Algorithm RAND

Input: multiplier a (p bits) and seed [x.sub.i-1] (p bits), a is one of the forms of Table I.

Output: a random number [x.sub.i] = a [multiplied by] [x.sub.i-1] mod [2.sup.p] - 1.

Method:

Compute x[prime] using the corresponding equation in Table I.

IF (x[prime] [less than] 0)

THEN Output x[prime] + m;

ELSE Output x[prime];

End of algorithm.

Table I. Equations of Four Forms of Multipliers and Their Equations Form Equation [2.sup.k1] - [2.sup.k2] x[prime] = [w.sub.1]-[w.sub.2] [2.sup.k1] + [2.sup.k2] x[prime] = [w.sub.1] + [w.sub.2] - m m - [2.sup.k1] + [2.sup.k2] x[prime] = -[w.sub.1] + [w.sub.2] m - [2.sup.k1] - [2.sup.k2] x[prime] = m - [w.sub.1] - [w.sub.2]

4. SPECTRAL TEST AND EMPIRICAL TEST RESULTS

This section first presents spectral test results for generators with multipliers of [+ or -] [2.sup.k1] [+ or -] [2.sup.k2] and modulus [2.sup.31] - 1 and [2.sup.61] - 1. The notation used here follows that of Anderson [1990]. The relative quality measure is

[q.sub.k] = [v.sub.k]/([[Gamma].sub.k] [multiplied by] [m.sup.1/k]),

where [[Nu].sub.k] is the spectral test result on k-dimensional space, and m is the prime modulus. The constants [[Gamma].sub.k], 2 [less than or equal to] k [less than or equal to] 8 were shown in Knuth [1981, p. 105] and Anderson [1990, Table 4.1]. The measure [q.sub.k] is a real number in [0, 1]. Large values of [q.sub.k] are best. We applied spectral test to all multipliers [+ or -] [2.sup.k1] [+ or -] [2.sup.k2] for m = [2.sup.31] - 1 and m = [2.sup.61] - 1. These multipliers are rated by the minimum of [q.sub.k], 2 [less than or equal to] k [less than or equal to] 8. based on this measurement, the best multipliers for m = [2.sup.31] - 1 are m - [2.sup.16] - [2.sup.11] and [2.sup.15] - [2.sup.10] (see Table II), and the best multipliers for [2.sup.61] - 1 are [2.sup.42] - [2.sup.31] and [2.sup.30] - [2.sup.19] (see Table III). The measure [Mathematical Expression Omitted] is interpreted [TABULAR DATA FOR TABLE II OMITTED] [TABULAR DATA FOR TABLE III OMITTED] [TABULAR DATA FOR TABLE IV OMITTED] as the number of bits that are "random" when k-tuples are considered [Anderson 1990].

How good are these generators? Fishman and Moore [1986] presented an exhaustive spectral test of multipliers for modulus [2.sup.31] - 1 (see Table IV). The 410 multipliers selected as best (as corrected by Park and Miller [1988]) all meet the quality measure [q.sub.k] [greater than or equal to] 0.8, for 2 [less than or equal to] k [less than or equal to] 6. By extending k to 8 in spectral tests on all these multipliers, we find that the largest min [TABULAR DATA FOR TABLE V OMITTED] [q.sub.k] = 0.7229 when a = 1754050460. The multiplier 742938285, which was considered the best, has min [q.sub.k] = 0.6211. On the other hand, in Table II, the multiplier m - [2.sup.16] - [2.sup.11] also has min qk = 0.6211. Table V summarizes the quality measure of these multipliers. based on the spectral test results, multipliers m - [2.sup.16] - [2.sup.11] and [2.sup.15] - [2.sup.10] are better than the three widely used multipliers 16807, 397204094, and 630360016.

Because there are no well-known multipliers for modulus m = [2.sup.61] - 1, we compare the multipliers [2.sup.42] - [2.sup.31] and [2.sup.30] - [2.sup.19] with two multipliers found by the following methods:

(1) Search for the best multiplier a = [37.sup.b] mod m, where b is relatively prime to m - 1. Because 37 is the minimal primitive root of m, all these multipliers are primitive roots of m [Fishman and Moore 1986]. This search is done in the range b[Epsilon] [1,[10.sup.6]].

(2) Search for the best multiplier a = [2.sup.k] [+ or -] 1. The number a must also be a primitive root.

Method (1) performs a time-consuming search to find a rather good multiplier. Method (2) searches for multipliers of a more restricted form, which is easier to compute than those of the general form [+ or -] [2.sup.k1] [+ or -] [2.sup.k2] mod m. The results of these searches are shown in Table VI. The best multiplier in method (1), a = [37.sup.458191] mod m, has min [q.sub.k] = 0.7129, and the best multiplier in method (2), a = [2.sup.38] - 1, has min [q.sub.k] = 0.0073. This result shows that [2.sup.42] - [2.sup.31] and [2.sup.30] - [2.sup.19] are far better than [2.sup.38] - 1, and using the form a = [2.sup.k] [+ or -] 1 is not likely to yield a good multiplier.

Table VII summarizes spectral test results on multipliers for [2.sup.61] - 1. Although the min [q.sub.k] values of [2.sup.42] - [2.sup.31] and [2.sup.30] - [2.sup.19] are not very high, they are still higher than that of 16807, a commonly used multiplier. Moreover, because [q.sub.k] is not an absolute measure, multipliers with different m cannot be compared using [q.sub.k]. A large modulus m is better based on the measure [[Beta].sub.k]. Generators with multipliers [2.sup.42] - [2.sup.31] and [2.sup.30] - [2.sup.19] and modulus [2.sup.61] - 1 are better than any generators with modulus [2.sup.31] - 1.

We apply three empirical tests on the 61-bit generators with multipliers [2.sup.30] - [2.sup.19] and [2.sup.42] - [2.sup.31] and the 31-bit generators with m - [2.sup.16] - [2.sup.11], [TABULAR DATA FOR TABLE VI OMITTED] [TABULAR DATA FOR TABLE VII OMITTED] [2.sup.15] - [2.sup.10]. The frequency test takes x rood 12. The run test examines the length of "run up" sequences, categorized into [1, 6] and [greater than] 6. The maximum-of-t test takes the maximum of 5 consecutive numbers and examines whether max (x) [less than] 7/8m. The degrees of freedom of these tests are 11, 6, and 1, respectively. To analyze the test results, we use the chi-square test [Knuth 1981, p. 44] in 6 rounds. Each test consumes 2 million random numbers. In total, our experiment consumes 36 million consecutive random numbers generated by seed 1. The chi-square result V is rejected if V is outside [1%, 99%]. The result is "suspect" if V is in [1%, 5%] or [95%, 99%]. The result is "almost suspect" if V is in [5%, 10%] or [90%, 95%]. The following are the results for the 61-bit generators with [2.sup.30] - [2.sup.19] and [2.sup.42] - [2.sup.31] and the 31-bit generators with m - [2.sup.16] - [2.sup.11], [2.sup.42] - [2.sup.31], 16807, and 1754050460. The latter two are included merely for comparison. All these generators are satisfactory.

5. PORTABILITY AND PERFORMANCE RESULTS

Figure 2 shows the C program for the generator [x.sub.i] ([2.sup.30] - [2.sup.19]) [x.sub.i-1] mod ([2.sup.61] - 1). This program can easily be adapted for any generator of the form [+ or -] [2.sup.k1] [+ or -] [2.sup.k2] mod m. LOG_W denotes the number of bits in integer types, and LOG_M denotes [Mathematical Expression Omitted], i.e., LOG_M = p, letting m = [2.sup.p] - 1. K1 and K2 represent k1 and k2. The equation from Table I should be plugged into the code. The "unsigned" variable x0 makes the right shift operator ([much greater than]) fill zeros on the left. The program assumes that the size of "long long int" type in C is 64 bits wide. For C compilers not supporting such an integer type, [TABULAR DATA FOR TABLE VIII OMITTED] the programmer needs to handle double-word integers by hand. Because the program contains only shifting and addition operations, coding these operations by hand is not difficult.

The code in Figure 2 is general for any k1 and k2. For a specific generator, the code can be simplified. For example, the recurrence for a = [2.sup.30] - [2.sup.19] can be rewritten in the C fragment:

#define M 0x1fffffffffffffffll

long long int x;

x = (x [much greater than] 31) + ((x [much less than] 30) & M)

- (x [much greater than] 42) - ((x [much less than] 19) & M);

if (x [less than] 0)

X + = M;

The bitwise-and operator (&) is used to get the high-order bits. This gets simpler code; however, compilers may generate the code that loads constant M from memory, which is usually less efficient than bit shift.

We compared the performance of our generators with other development techniques. A test was conducted that generated 2 million random integers for all following programs. Program m61-port is the generator by Park and Miller [1988] with modulus [2.sup.61] - 1. The codes on Sun SPARC and HP PA-RISC were compiled by GNU CC (gcc -O); the codes in Alpha and RS/6000 were compiled by the vendor's C compiler (cc -O). Program m61-asm denotes generators with modulus [2.sup.61] - 1 developed using assembly languages. We had two implementations: one on the Sun SPARC and one on the DEC Alpha. The SPARC assembly code allows only 32-bit multipliers, while the Alpha assembly code allows any 61-bit multipliers. Program m61-p3019 denotes the 61-bit generator with a = [2.sup.30] - [2.sup.19]. All these 61-bit generators generated the same sequence of random numbers. Program m31-asm is an assembly version of a 31-bit generator with any multiplier a and modulus m = [2.sup.31] - 1. The time was measured with a = 16,807. Program m31-p1611 denotes a 31-bit generator with a - m - [2.sup.16] - [2.sup.11]. Program m31-p1510 denotes a 31-bit generator with a = [2.sup.15] - [2.sup.10]. Finally, the test on the generator with modulus [2.sup.32] was also performed. The results show that our generators (m61-p3019, m31-p1611, and m31-p1510) are very fast.

[TABULAR DATA FOR TABLE IX OMITTED]

6. CONCLUSION AND FUTURE WORK

This article has presented the development of multiplicative congruential generators with prime modulus [2.sup.p] - 1 and four forms of multipliers: [2.sup.k1] - [2.sup.k2], [2.sup.k1] + [2.sup.k2], m - [2.sup.k1] + [2.sup.k2], and m - [2.sup.k1] - [2.sup.k2], k1 [greater than] k2. The multipliers for modulus [2.sup.31] - 1 and [2.sup.61] - 1 have been measured by spectral tests. The generators with these multipliers are portable and very fast. They have also passed several empirical tests, including the frequency test, the run test, and the maximum-of-t test. In the future, more tests and more experience with these generators are needed.

ACKNOWLEDGMENTS

The author would like to thank the referees and the editor, whose comments helped to improve the overall presentation.

REFERENCES

ANDERSON, S.L. 1990. Random number generations generators on Vector supercomputers and other advanced architectures. SIAM Rev. 32, 2 (June), 221-251.

FISHMAN, G. A. AND MOORE, L.R. 1986. An exhaustive analysis of multiplicative congruential random number generators with modulus [2.sup.31]-1. SIAM J. Sci. Stat. Comput. 7, 1, 24-45.

KIRKPATRICK, S. AND STOLL, E. 1981. A very fast shift-register sequence random number generator. J. Comput. Phys. 40, 517-526.

KNUTH, D.E. 1981. The Art of Computer Programming. Vol. 2, Seminumerical Algorithms. Addison-Wesley, Reading, Mass.

PARK, S. K. AND MILLER, K.W. 1988. Random number generators: Good ones are hard to find. Commun. ACM 31, 10 (Oct.), 1192-1201.

PAYNE, W. H., RABUNG, J. R., AND BOGYO, T.P. 1969. Coding the Lehmer pseudo-random number generator. Commun. ACM 12, 2, 85-86.

SCHRAGE, L. 1979. A more portable Fortran random number generator. ACM Trans. Math. Softw. 5, 2, 132-138.

SELKE, W. 1993. Cluster-flipping Monte Carlo algorithm and correlations in "good" random number generators. JETP Lett. 58, 8, 665-668.

Printer friendly Cite/link Email Feedback | |

Author: | Wu, Pei-Chi |
---|---|

Publication: | ACM Transactions on Mathematical Software |

Date: | Jun 1, 1997 |

Words: | 3600 |

Previous Article: | Algorithm 770: BVSPIS - a package for computing boundary-valued shape-preserving interpolating splines. |

Next Article: | Computational investigations of low-discrepancy sequences. |

Topics: |