Evaluating Parallel Random Number Generators^{∗}

W.E. Mahoney^{†}  P.D. Coddington^{‡}

1993 NPAC REU Project

Abstract Parallel random number generators were studied and tested on a variety of parallel computers. The performance of each random number generator was evaluated using parallel Monte Carlo simulations. A model with known results was simulated to enable a more precise analysis of the random number generator that was used. The results of these parallel simulations were similar to those achieved with sequential simulations.

1 Introduction

The study of Random Number Generators (RNGs) has become (and will continue to be) very important because of the expanding magnitude of scientific research, and the increasing demand for more accurate results. RNGs are utilized in the studies of probability and statistics and cryptography, and are an integral part of Monte Carlo simulations done in physics and chemistry [2, 8]. For example, they can be used to simulate probabilistic events such as nuclear reactions. These Monte Carlo simulations are very dependent upon the RNGs that they use because any correlations in the RNGs sequence can lead to biases, and incorrect results will be obtained [7].

To conceive of how a RNG functions, it is important to understand that the numbers being generated are not truly random; they are pseudo-random. This means that the behavior of the RNG algorithm is inherently predeterministic (i.e., each successive number is known before it is generated). In fact, many RNGs prompt the user for an initial value called the seed. By knowing the seed, it is possible to calculate the sequence of random numbers that a particular algorithm will generate.

There is a finite quantity of random numbers that an RNG algorithm can generate before the same exact sequence begins to repeat. This quantity is known as the RNG’s period. The period should always be longer than the number of calls to the RNG in a given application/program, or else the results of the process will be correlated in some manner. A poor RNG algorithm or badly chosen parameters can have an adverse effect on the period by shortening it, giving the user fewer random numbers to work with.

Without going into further detail about how random numbers are generated, it is known where the generated numbers start, stop, their sequence, and how many to expect. This

∗Project conducted as part of the 1993 Research Experiences for Undergraduates (REU) Program in High-Performance Computing conducted by the Northeast Parallel Architectures Center (NPAC) at Syracuse University. The National Science Foundation provides primary funding for the 1993 NPAC REU Program, through NSF Grant CDA-9200577. Additional funding is provided by the GE Foundation Faculty for the Future Program; the Vice President for Research and Computing, the Dean of the College of Engineering and Computer Science, the Dean of the College of Arts and Sciences, and NPAC, Syracuse University.

†Research Apprentice, 1993 NPAC REU Program; Computer Engineering Major, Syracuse University. ‡Project Leader, Northeast Parallel Architectures Center, Syracuse University.

predictability is considered good, not bad, because the results of many simulations that utilize RNGs must be exactly reproduced in order to effectively study them. Also, this characteristic is a necessity when debugging programs that use RNGs.

Now, just because RNGs can mysteriously supply a seemingly random sequence of integers does not mean that it is a “good” RNG. To be considered “good” an RNG should have the following characteristics [4]:

• Large period

• Non-correlated numbers

• Uniform distribution

• Adjustable seed

• Fast algorithm

• Uses little memory

There are statistical tests that check the randomness properties of RNGs, however a Monte Carlo simulation provides a more stringent test for RNGs. A Monte Carlo simulation is a probabilistic method of approximating a multi-dimensional integral [2, 8]. This type of simulation is responsible for exposing “bad” RNG algorithms that were previously thought to be “good” [6].

2 Kinds of RNGs 2.1 Multiplicative Linear Congruential

The Multiplicative Linear Congruential (MLC) generator was first proposed by D. H. Lehmer in 1948, and is the oldest RNG in use today [1, 9]. The sequence is expressed by:

x_{n} = (zx_{n−1} + c) mod m

Where z is the multiplier, x_{n−1} is the previous number in the sequence, c is an added constant, and m is the modulus. The initial value of x_{n−1} is the seed value (x_{0}).

Parameters for the MLC RNG must be chosen carefully in order to get the maximum performance. For example, one implementation of this generator will give a sequence of alternating odd and even integers, while another will show a correlation in the low order bits [1]. The most popular MLC [12] does not achieve a full period, but the numbers generated have very good randomness in their low order bits. This sequence is achieved by setting z = 16807, c = 0, and m= 2^{31}-1, with the seed being in the range 1 ≤ m.

2.2 Lagged Fibonacci

Like the MLC, the Lagged Fibonacci (LF) generator uses a modulo operator, but the similarities end there. An LF generator utilizes a table and lags into it twice and performs a specified operation on the two numbers found. The new number is then placed at the beginning of the table. The LF RNG has the form:

x_{n} = (x_{n−p} OP x_{n−q}) mod m

Where p and q are the lags into the table from n (current number), and OP = +, −, or ∗.

A problem for LF generators is that for small lags the low order bits show some correlation, but this seems to disappear when very long lags are used. This kind of generator is fairly new and currently has no preferred implementation.

2.3 Shift Register

The Shift Register (SR) generator is, in fact, an LF generator that uses XOR as its OP operator, and does not utilize the modulo operator at all. An SR’s sequence is generated by:

x_{n} = (x_{n−p} XOR x_{n−q})

RNGs of this form are known to be some of the worst generators of all [10], but amazingly enough are still found in use. Most of the reason for the SR’s bad reputation stems from the fact that the XOR operator does not mix neighboring bits. Therefore the numbers remain highly correlated throughout their run.

2.4 Cellular Automata

The Cellular Automata (CA) generator performs logical operations on the bits within a number to obtain the next number of the sequence. The sequence is generated by:

a_{n+1}^{i} = a_{n}^{i−1} XOR (a_{n}^{i} OR a_{n}^{i+1})

Where i represents bits being operated on in the current set of n, and n+1 is a newly generated bit that will contribute to forming the next number.

The CA RNG mixes neighboring bits and shows no correlations in the low order bits of the numbers it generates. In fact, despite its simple structure the CA is capable of complex and apparently chaotic behavior [17]. Although all the implementations of CA generators have not been studied and tested, they seem to perform quite well.

3 Method of Testing

To test the various Parallel Random Number Generators (PRNGs), the 2-D Ising model of magnetism on a 128 x 128 lattice was used. The 2-D Ising model is a very simple model that, given the relationship between specified parameters, will change a cell on the lattice to a 1 or a 0 (on or off etc.). More importantly, this model has been solved and its exact result was compared to the results of the simulation.

This model was implemented using two of the three most widely used Monte Carlo simulations [4]. The Metropolis algorithm [2, 11] updates a single site of the lattice at a time, and the Swendson-Wang algorithm [13] forms clusters of sites that are updated collectively. In the near future the Wolff algorithm [16] will be used, and it updates a single cluster of sites at a time. The PRNGs are used differently by each algorithm, therefore it cannot be assumed that if a PRNG does well with one of these algorithms that it will do well with another. Each of these algorithms approximates the energy and specific heat of the Ising model. Therefore, a “good” PRNG will give better approximations of the exact results.

The University of Maryland’s CM2 was used to test a cellular automata PRNG called CMF-RANDOM [14] which was utilized by both the Metropolis and Swendson-Wang algorithms. The instruction CMF-RANDOMIZE initializes CMF-RANDOM with a user specified seed value. Since this is a PRNG, each bit of the seed value is placed on an individual processor, so each processor performs the appropriate cellular automaton operations when CMF-RANDOM is called.

Syracuse University’s CM5 was used to implement a lagged Fibonacci PRNG using addition as its operator. This PRNG is found in the CMSSL library and is called FAST-RNG [15]. FAST-RNG was tested via the Metropolis algorithm. Again CMF-RANDOMIZE is used to initialize CMF-RANDOM with a seed value, but then INITIALIZE-FAST-RNG is called to set up the lag table using CMF-RANDOM. INITIALIZE-FAST-RNG allows the user to specify the size of the lags to be used. Lags of 17 and 5 were used for these tests. Finally, FAST-RNG is called to generate a random number.

4 Results

All Figures and Tables are located at the end of this paper.

The results were as expected for each of the simulations. CMF-Random (CM2) was expected to do well based upon the research of Wolfram [17], and it did perform up to expectations when utilized by the Swendson-Wang and Metropolis Monte Carlo simulations. Also, results from the use of FAST-RNG (CM5) were as expected, since its results concurred with that of its sequential counterpart [4]. FAST-RNG performed well despite the small lags used (17 and 5). As previously discussed, correlations in LF generators with small lags can lead to biased or erroneous results. The flawed use of this LF generators sequential counterpart is not discovered until tested with the Wolff algorithm, but did perform well using the Metropolis algorithm [4] as did this parallel version.

Both of the Metropolis simulations were run 20 times and the Swendson-Wang 10 times, each run utilizing different seed values. The standard deviation from the exact result, of both the energy and specific heat, was calculated for each run (see Figures ?? and ??), and the error in the mean was calculated for each of these data sets. Finally, these calculations were used to calculate values for chi squared per degree of freedom (χ^{2}) and the deviation (Δ) between the mean results and exact value, as a multiple error in the mean (σ). The first test checks for differences in the statistical fluctuations between each run, the second checks for any bias in the average of all runs. A PRNG was considered “good” if χ^{2} and Δ met the following requirements:

0.34 < χ^{2} < 2.0

and Δ < 3.3σ

Any value outside of each range should occur with the probability of less than 0.001 for a truly “good” PRNG [3].

χ^{2} was calculated by summing the squares of the standard deviation of each run and dividing by the number of runs. Δ was calculated by subtracting the exact value from the mean results and dividing that difference by the error in the mean. The error in the mean was determined by taking the square root of the average of the square of the results minus the square of the average result divided by the number of runs. The results of these calculations are depicted in Table 1. As you can see all of the values fall into the required range and can be considered “good” at this point in time.

5 Conclusion

There does not seem to be any correlations produced by parallelizing RNG algorithms, including inter-processor correlations for the two PRNGs tested. By the standards described herein we can say that the tested PRNGs are “good”, but further testing is expected to reveal that at least FAST-RNG is “bad” when used by certain algorithms (such as Wolff). Currently, more testing is being done on the PRNGs presented here, and other PRNGs, in order to properly evaluate their quality.

Studying and Testing Parallel Random Number Generators 5

Acknowledgments

The first author would like to thank the second author for all of his answers to my many questions concerning this topic. Both authors would like to thank John Apostolakis for his help with the CM2.

References

[1] Anderson, S.L., SIAM Rev. 32 (1990) 221.

[2] Binder, K., ed., Monte Carlo Methods in Statistical Physics, (Springer-Verlag, Berlin, 1986); Binder, K., Heermann, D.W., Monte Carlo Simulations in Statistical Physics, (Springer-Verlag, Berlin, 1988)

[3] Bevington, P.R., Data Reduction and Error Analysis for the Physical Sciences, (McGraw-Hill, New York, 1969).

[4] Coddington, P.D., Analysis of Random Number Generators Using Monte Carlo Simulation, Northeast Parallel Architectures Center Technical Report, Syracuse University, June 1993.

[5] Ferdinand, A.E., Fisher, M.E., Phys. Rev. 185 (1969) 832.

[6] Ferrenberg, A.M., Landau, D.P., and Wong, Y.J., Phys. Rev. Lett. 69 (1992) 3382.

[7] Fox, G., et al, Solving Problems On Concurrent Processors, Vol. 1, (Prentice Hall, Englewood Cliffs, New Jersey, 1988).

[8] Gould, H., Tobochnik, J., An Introduction to Computer Simulation Methods, Vol. 2, (Addison-Wesley, Reading, Mass., 1988)

[9] Knuth, D.E., The Art of Computer Programming Vol. 2: Seminumerical Methods (Addison-Wesley, Reading, Mass., 1981).

[10] Marsaglia, G.A., A Current View of Random Number Generators, in Computer Science and Statistics: The Interface, L. Billard, ed., (Elsevier, Amsterdam, 1985)

[11] Metropolis, N., et al, J. Chem. Phys. 21 (1953) 1087.

[12] Park, S.K., Miller, K.W., Comm. ACM 33:10 (1988) 1192.

[13] Swendson, R.H., Wang, J.-S., Phys. Rev. Lett. 58 (1987) 86.

[14] Thinking Machines Corp., CM Fortran User’s Guide, Ver. 1.0, (Cambridge, Mass., 1991), p. 66.

[15] Thinking Machines Corp., CMSSL for CM Fortran, Ver. 2.2, (Cambridge, Mass., 1991), p. 298.

[16] Wolff, U., Phys. Rev. Lett. 62 (1989) 361.

[17] Wolfram, S., Advances in Applied Mathematics 7 (1986) 123.