A Pseudorandom Number Generator for the Parallel NAS Kernels

D. Bailey

Suppose that n uniform pseudorandom numbers are to be generated. Set and let be a specified initial ``seed,'' i.e., an integer in the range . Generate the integers for using the linear congruential recursion

and return as the results. Thus , and the are very nearly uniformly distributed on the unit interval. See [2], beginning on page 9 for further discussion of this type of pseudorandom number generator.

Note that any particular value of the sequence can be computed directly from the initial seed s by using the binary algorithm for exponentiation, taking remainders modulo after each multiplication. To be specific, let m be the smallest integer such that , set b = s and t = a. Then repeat the following for i from 1 to m:

The final value of b is . See [2] for further discussion of the binary algorithm for exponentiation.

The operation of multiplying two large integers modulo can be implemented using 64 bit floating point arithmetic by splitting the arguments into two words with 23 bits each. To be specific, suppose one wishes to compute . Then perform the following steps, where int denotes the greatest integer:

An implementation of the complete pseudorandom number generator algorithm using this scheme produces the same sequence of results on any system that satisfies the following requirements:

• The input multiplier a and the initial seed s, as well as the constants and , can be represented exactly as 64 bit floating point constants.

• The truncation of a nonnegative 64 bit floating point value less than is exact.

• The addition, subtraction and multiplication of 64 bit floating point values, where the arguments and results are nonnegative whole numbers less than , produce exact results.

• The multiplication of a 64 bit floating point value, which is a nonnegative whole number less than , by the 64 bit floating point value , produces an exact result.

These requirements are met by virtually all scientific computers in use today. Any system based on the IEEE-754 floating point arithmetic standard [1] easily meets these requirements using double precision. However, it should be noted that obtaining an exact power of two constant on some systems requires a loop rather than merely an assignment statement with **.

Other Features

• The period of this pseudorandom number generator is , and it passes all reasonable statistical tests.

• This calculation can be vectorized on vector computers by generating results in batches of size equal to the hardware vector length.

• By using the scheme described above for computing directly, the starting seed of a particular segment of the sequence can be quickly and independently determined. Thus numerous separate segments can be generated on separate processors of a multiprocessor system.

• Once the IEEE-754 floating point arithmetic standard gains universal acceptance among scientific computers, the radix can be safely increased to , although the scheme described above for multiplying two such numbers must be correspondingly changed. This will increase the period of the pseudorandom sequence by a factor of 64 to approximately .