1. Introduction

The extra-accurate fused multiply–add (FMA) operation of the RS/6000* and PowerPC* family of RISC microprocessors offers many opportunities to use mathematical innovation to produce fast algorithms for numerically intensive computation (NIC). In this paper we illustrate this assertion by giving several examples and by demonstrating a 50-fold increase in performance (over generic algorithms [1]) for pseudorandom-number generation. The results obtained are bit-wise identical to the results of the generic algorithms.

For multiplicative congruential generators [2] of the type

si+1 = asi mod 2k, xi+1 = 2–ksi+1 or 2–k+1si+1 – 1,

(1)

with k 52, we show that the cost per random number for the two distribution intervals (0, 1) and (–1, 1) is respectively 3 and 3.125 multiply–adds on RS/6000 processors.

We also specify a linear congruential generator of the form

s_i+1 = (as_i + c) mod 2^k,

x_i+1 = 2^–ks_i+1,

related to the generator in [2], which has a full period of 2^k. The cost per random number [in the range (0, 1)] for this generator is four multiply–adds on RS/6000 processors.

Finally, we present a multiplicative congruential generator for the interval (0, 1) of the form

s_i+1 = as_i mod (2^k – 1),

xi+1 = si+1  , 2k – 1

(2)

(discussed in [3]), for which the modulus is not a power of 2. Such a generator, as well as a generator of type (1), is selectable as the generator used in the RANDOM_NUMBER intrinsic function of IBM XL Fortran (XLF) [4] and XL High Performance Fortran (XLHPF) [5]. [By “a generator of type (n),” or “generator (n),” we mean a generator whose description is given by equation(s) (n).]

We first introduce some notation. We use q to represent the modulus, either 2^k or 2^k – 1, depending on the generator. Arithmetic operators in equations are exact. Operands are IEEE normalized numbers. We use the operators , , , ø for IEEE double-precision (64-bit) floating-point arithmetic, and fl(x) for the correctly rounded double-precision value corresponding to x. In most cases, x = ab + c, where a, b, and c are IEEE numbers. Unless otherwise stated, the rounding mode we use is round-to-zero (chop). This leads to the best performance.

Let a, b, and c be arbitrary IEEE 64-bit floating-point numbers. The fused multiply–add operation on RS/6000 and PowerPC computes the correctly rounded d = fl(ab + c) for any of the four IEEE rounding modes. On RS/6000 machines, all 106 bits of the product a * b are added to c in order to guarantee that d will be correctly rounded for all possible values of a, b, and c. POWER and POWER2 RS/6000 machines can on a continuous basis compute one and two FMAs every machine cycle if the results of the FMAs do not cause any pipeline delays. The pipeline delay for these machines is two or three cycles. Loop unrolling can be used to avoid any pipeline delays for the pseudorandom-number calculation.

We first consider a multiplicative congruential pseudorandom-number generator for which the modulus is a power of 2. See Knuth [2] for a thorough discussion of pseudorandom-number generators. Let 0 < s₀ < 2^k, s₀ odd, 1 < a < 2^k, k 52, and for i 0, let

si+1 = asi mod 2k, xi+1 = 2–ksi+1.

(3)

For k = 46 and a = 5¹³ we have a specific instance of such a generator. This generator has period 2^k–2, and it is extensively used by the Numerical Aerodynamic Simulation (NAS) suite of paper and pencil benchmarks [1]. We mention here that the random-number generator (3) is embarrassingly parallel, or EP. In [1, p. 31, bullet 3], this point is made for the EP benchmark. Also, on page 29 of [1], Bailey describes the binary algorithm for exponentiation that allows one to compute aⁿ mod 2^k in log₂(n) steps. The fact that aⁿ mod 2^k is computable in log₂(n) steps is crucial to making the random-number generator (3) embarrassingly parallel. Bailey implicitly points this out in [1, p. 31, bullet 2]. On the IBM POWER2 Model 590, our bit-wise identical algorithms corresponding to Equation (3) compute more than 40 million random numbers per second. On the IBM SP2 machine, using the Model 590 for its nodes, and using p such nodes, we can compute more than 40p million random numbers per second because of the EP nature of these algorithms. Thus, using 25 nodes our algorithms can compute more than a billion random numbers in a second.

In [1], Bailey gives a generic algorithm that simulates base 2²³ multiple-precision arithmetic to compute the s_i's in (3). This algorithm requires 18 floating-point operations and four convert-to-integer operations per random number. We implement the same generator using three multiply–adds. Any pseudorandom s_i from (3) can then be placed in the range 0 < x_i < 1 by the scaling x_i = 2^–ks_i or be put in the range –1 < x_i < 1 by computing x_i = 2^–k+1s_i – 1. These computations respectively require one and two additional floating-point operations. In this paper, we redefine (3) to compute each x_i directly, without first computing s_i, and thus compute

x_i+1 = ax_i mod 1.

This change avoids the actual scalings done above. We also show how these computations can be done by three and 3.125 multiply–adds per random number on RS/6000 machines.

When doing modular arithmetic on integers, it is natural to use the greatest integer function. In floating-point arithmetic this is achieved by using the IEEE round-to-zero (chop) rounding mode. Throughout this paper, except for one place in Section 2, we use the chop rounding mode.

Many NIC algorithms are rated by their megaflop rate, although this popular measure is often misleading. We think that pseudorandom-number generation is such a case. An NIC computation related to pseudorandom generation is the EP NAS benchmark [1, 6]. By using mathematical innovation and the fast random-number generation described here, we were able to demonstrate that the RS/6000 POWER2 Model 590 could currently outperform a Cray YMP for EP [6]. In this paper we compare the ratio of the computing time for the generic algorithm [1] to the time of our algorithm. In both cases we use the same model of RS/6000. The new algorithm is written entirely in Fortran and is compiled using the XL Fortran (XLF) compiler [4].

We also implemented another generator of the type (3) in which a = 44485709377909 and k = 48. This generator is the RANF( ) pseudorandom function used on the CDC Cyber 174 computer, and is also generator 2 in the RANDOM_NUMBER intrinsic function provided with the IBM XLF and XLHPF compilers. It is described in the book Stochastic Simulation by Brian Ripley [7, p. 216], who presents statistical test results for several generators; he finds this generator “quite acceptable.” Gordon Slishman has implemented this generator¹ using three multiply–adds for single-processor execution of RANDOM_NUMBER. The parallel implementation in XLF and XLHPF is described in Section 4.

Slishman also implemented the generator (2), for which the modulus is not a power of 2, for generator 1 of single-processor RANDOM_NUMBER. Parallel versions of both generators 1 and 2 in RANDOM_NUMBER were implemented by Robert Enenkel for SP (distributed-memory) machines, and by Enenkel and Xinmin Tian² for SMP (shared-memory) machines. Generator (2) is discussed in detail in Section 5.

In Section 2 we describe some elementary mathematical ideas that can be used to compute (3) rapidly. These ideas are used to derive an algorithm to compute pseudorandom numbers in the range (0, 1) in three multiply–adds and pseudorandom numbers in the range (–1, 1) in 3.125 multiply–adds. Also, using the IEEE round-to-nearest mode, we show how the latter computation can be done in three multiply–adds.

We now discuss POWER and PowerPC models. Before the advent of vector processors, integer arithmetic was generally faster than floating-point arithmetic. It was then customary to produce pseudorandom numbers using fixed-point arithmetic and then convert these integers to floating point with scaling to get floating-point pseudorandom numbers. When fast floating-point processors became available, it became more economical to produce the random numbers directly in floating point. For example, the first equation of generator (1) could be computed on an RS/6000 by setting 0 < s₀ < 2^k, s₀ odd, and letting, for i 0,

u = fl(as_i + 2^52+k),

v = u 2^52+k,

s_i+1 = fl(as_i – v).

The above three statements define pseudorandom integers in the range 0 < s_i < 2^k that are bit-wise identical to the generator (1). However, one usually wants random floating-point numbers in the range (0, 1) or (–1, 1). The formulas above would then require modification; i.e.,

x_i = 2^–ks_i or x_i = 2^–k+1s_i – 1.

However, these computations require an extra multiply–add in addition to the four needed to compute s_i. One intention of this paper is to demonstrate that this additional multiply–add can be removed. This results in improvement factors of 4/3 and 4/3.125, which come to 33% and 28% improvements per iteration over computation based on the above four statements.

In Section 3 we describe two full-period linear congruential generators,

xi+1 = (axi + c) mod 1,

(4)

that compute pseudorandom numbers in the range (0, 1) in four multiply–adds for c = 2^–k and c = a2^–k. We show a proof that these generators are EP. The proof involves showing that (1 + a + · · · + a^n–1) mod q can be computed in log₂(n) steps. Our four-multiply–add implementation of these generators requires us to avoid conditional operations in the unrolled inner loop. It turns out that ordinary unrolling via vectorization fails. To overcome this failure, we introduce “generalized unrolling,” which becomes possible because the generator is EP.

In [2, Section 3.6, pp. 170–173], Knuth provides a summary on “how to generate random numbers.” He recommends using a linear congruential generator,

X (aX + c) mod q,

(5)

that satisfies seven properties. However, he states that this class of generator applies primarily to machine-level coding and hence is not portable. For IEEE arithmetic, it makes sense to choose q = 2^k. The fused multiply–add instruction is provided by many computer architectures, including the IBM RS/6000 POWER and PowerPC, IBM S/390* G5, Intel/HP IA-64 Itanium**, Apple Power Mac**, HP Precision Architecture RISC 2.0 (e.g., HP900/800), and SGI MIPS** R-8000.³ Thus, within this more restrictive domain, we can propose our four-multiply–add generator as being “portable.”

In Section 5 we extend the ideas of the previous sections to the generator (2). The modulus of this generator is 2³¹ – 1, and not a power of 2 as for the previous generators. To achieve high performance, nontrivial changes to the algorithm are required.

In Section 6 we describe the performance of these algorithms relative to the generic algorithm of [1].

2. Three multiplicative congruential generators

Let a and y be positive integers less than 2^k, where k 52. Clearly a and y are IEEE numbers. Let x = 2^–ky so that 0 < x < 1. Also, x is an IEEE number. Let p = ax. The base-2 (bit) representation of p is p₁p₂ · · · p_k.p_k+1p_k+2 · · · p_2k, where each p_i is 0 or 1. Represent p = I + F, where I = p₁ · · · p_k and F = .p_k+1 · · · p_2k; i.e., I and F are the integer and fractional parts of p. Note that ax mod 1 = F, that I and F are IEEE numbers, and that 2⁵² p + 2⁵² 2⁵³ – 1. Let

u = fl(ax + 252).

(6)

All 2⁵² positive integers in the range 2⁵² to 2⁵³ – 1 are all of the IEEE numbers in that range. Thus, u in (6) is an integer, and since we are in chop mode, u is the largest integer not exceeding p + 2⁵²; i.e., u = p + I. (See Lemma 1 for a detailed constructive proof.) Let

v = u 252.

(7)

The computation in (7) is done exactly, since v is computed using IEEE arithmetic, in which u, 2⁵², and v = I are all IEEE numbers. Now let

w = fl(ax – v).

(8)

On RS/6000, w is the fractional part F of ax because the fused multiply–add instruction always delivers the correct answer when its operands and result are IEEE numbers. Note that ax – v = F. Thus,

w = ax mod 1.

(9)

We have just proved the following.

Theorem 1
The computation in (6), (7), and (8) above produces the result (9). The value computed by (8) is bit-wise identical to the infinitely precise result (9).

When unrolled, Equations (6), (7), and (8) constitute a three-multiply–add implementation of the random-number generator (3). Let a = 5¹³ and choose 1 < s₀ < 2^k with s₀ an odd integer and k = 46. Set x₀ = s₀2^–k. Then 0 < x₀ < 1. Note that the seven lower-order bits of x₀ are zero. In fact, each x_i, i 0, given by

xi+1 = axi mod 1

(10)

has this property.

We now discuss some aspects of the Algorithm and Architecture approach [8] and how it relates to unrolling. In Section 1.3 of [1], Bailey et al. describe sample codes for the NAS benchmarks. There were two codes distributed for random-number generation, namely RANDLC and VRANLC; the latter is of interest to this paper. Subroutine VRANLC generates n REAL*8 uniform pseudorandom numbers in the range (0, 1) by using Equation (3). The documentation states that VRANLC is the standard version designed for scalar or RISC systems. A comment in this code states that the DO loop below in (11) which generates the n uniform pseudorandom numbers is not vectorizable.

       T1 =  R23 * A 

       A1 =  AINT (T1) 

       A2 =  A – T23 * A1 

  C 

  C Generate N results.  This loop is not

  C                      vectorizable. 

  C 

       DO 120 I = 1, N 

  C 

  C Break X into two parts such that

  C X = 223 * X1 + X2, compute 

  C Z = A1 * X2 + A2 * X1 (mod 223), and then 

  C X = 223 * Z + A2 * X2 (mod 246). 

  C 

       T1 =  R23 * X 

       X1 =  AINT (T1) 

       X2 =  X – T23 * X1 

       T1 =  A1 * X2 + A2 * X1 

       T2 =  AINT (R23 * T1) 

       Z  =  T1 – T23 * T2 

       T3 =  T23 * Z + A2 * X2 

       T4 =  AINT (R46 * T3) 

       X  =  T3 – T46 * T4 

       Y(I)= R46 * X 

We believe the authors' statement about the code in (11). Today's compilers are not able to vectorize this loop. On the other hand, Equation (3) is vectorizable. The Algorithm and Architecture approach to handling Equation (3) would be to rewrite the code so that it performs well when run on some actual machine, e.g., an RS/6000 RISC-type machine. In the present case we would also need to unroll Equations (6), (7), and (8). Please note that the code in (11) can be replaced by the following code, in which TP52 = 2⁵²:

DO i = 0, n – 1 

       u = a * x(i) + TP52 

       v = u – TP52 

       x(i+1) = a * x(i) – v

However, the above code will not execute in three cycles because of pipeline delays. This code is also not vectorizable by a compiler. In producing the code in (12), we have used the extra accuracy feature of the fused multiply–add instruction that is present on RS/6000 and PowerPC machines. Also, the code in (12) is equivalent to using Equation (9). A compiler cannot recognize this fact. This is where we use the Architecture feature of our approach. We now show how to “unroll” the code in (12). In doing so, we vectorize the code in (12).

Let a_j = a^j mod 2^k for 1 j m. Note that, from (3),

si+j	= ajsi mod q
	= ajsi mod q,

so that

xi+j  =  ajsi mod q q = ajxi mod 1, 1 j m.

(13)

Thus, given x_i and (13), we have defined the next m iterations of (10). This unrolling of (10) by m constitutes using a vector of length m to compute these next m iterations of (10). In effect, RS/6000 machines can be viewed as vector machines with a short vector length. Because RS/6000 machines possess functional parallelism (see [9]), RS/6000 machines have very low vector start-up costs. This fact implies that using a short vector length is not a performance drawback. However, to achieve this vectorization it must be programmed. Now let m = 4. Thus, the code in (12) becomes

DO    i = 0, n – 4, 4

      u = a * x(i) + TP52

      u2 = a2 * x(i) + TP52

      u3 = a3 * x(i) + TP52

      u4 = a4 * x(i) + TP52

      v = u – TP52

      v2 = u2 – TP52

      v3 = u3 – TP52

      v4 = u4 – TP52

      x(i+1) = a * x(i) – v

      x(i+2) = a2 * x(i) – v2

      x(i+3) = a3 * x(i) – v3

      x(i+4) = a4 * x(i) – v4

ENDDO

The above code eliminates most of the pipeline delays. Every target of an FMA in that code is separated by three independent FMAs. However, the last FMA of the loop is followed by the first FMA of the loop with i replaced by i + 4; thus, there is no separation between the target x(i + 4) and the target u. In order to eliminate all pipeline delays, we need a more complicated “unrolling” of the code in (12). Consider

xi = x(4) 

xi3 = x(3) 

xi2 = x(2) 

xi1 = x(1) 

DO i = 0, n – 8, 8 

       u4 = a4 * xi + TP52 

       u3 = a3 * xi + TP52 

       u2 = a2 * xi + TP52 

       u = a * xi + TP52 

       x(i+3) = xi3 

       x(i+4) = xi 

       v4 = u4 – TP52 

       v3 = u3 – TP52 

       v2 = u2 – TP52 

       v = u – TP52 

       xi0 = a4 * xi – v4 

       xi3 = a3 * xi – v3 

       x(i+1) = xi1 

       x(i+2) = xi2 

       xi2 = a2 * xi – v2 

       xi1 = a * xi – v 

       u4 = a4 * xi0 + TP52 

       u3 = a3 * xi0 + TP52 

       u2 = a2 * xi0 + TP52 

       u = a * xi0 + TP52 

       x(i+7) = xi3 

       x(i+8) = xi0 

       v4 = u4 – TP52 

       v3 = u3 – TP52 

       v2 = u2 – TP52 

       v = u – TP52 

       xi = a4 * xi0 – v4 

       xi3 = a3 * xi0 – v3 

       x(i+5) = xi1 

       x(i+6) = xi2 

       xi2 = a2 * xi0 – v2 

       xi1 = a * xi0 – v

The code in (14) eliminates all pipeline delays and hence should execute at the rate of one pseudorandom number every three multiply–adds. To start the pipeline, we need to precompute outside of the loop the first four pseudorandom numbers, x(i), i = 1, · · · , 4. The code in (14) is unrolled by 8.

Now we demonstrate another feature of the Algorithm and Architecture approach. Suppose we decide to use (13) to unroll by m = 8. We have seen that a compiler cannot unroll the present loop (12). We were forced to program the unrolling and to introduce the concept of a “vector instruction” into RS/6000 coding. We now claim that this necessity of having to program a “vector instruction” gives rise to a feature of superscalar machines that vector machines do not possess. This feature allows us to produce “customized vector instructions.” We remark that this is a part of the Algorithm and Architecture approach, and we now illustrate this concept with the example of generating uniform pseudorandom numbers in the range (–1, 1). A standard implementation requires an extra operation to convert the range from (0, 1) to (–1, 1).

Let c1 = 2⁵³ and c2 = 2⁵³ – 1. Let 0 < s < 2^k with s odd and k = 46. Let x = 2^–k+1s. Then 0 < x < 2. Let, for i 1,

u = fl(ajx + c1),

(15)

v = u c2,

(16)

and

xi+j = fl(ajx – v),

(17)

where 1 j < 8. For j = 8 we let

u = fl(a7x + c1),

(18)

v = u c1,

(19)

x = fl(a7x – v),

(20)

and

xi+8 = x 1.

(21)

Equations (15) to (21) define a loop, unrolled by 8, in which the next eight random iterates are computed and also the “seed” x is updated in Equation (20).

The first seven iterations cost three multiply–adds, while the last iteration costs four multiply–adds. The functional parallelism of RS/6000 indicates that this loop will execute in 25 multiply–adds, or 25/8 = 3.125 multiply–adds per iteration.

We close Section 2 by demonstrating a use of the IEEE “round-to-nearest” rounding mode. Let 0 < u_i < 1 be the x_i computed by Equation (10). For each u_i, let x_i = 2u_i – 1. Then –1 < x_i < 1. Let c1 = 2⁵³ + 2^k. Let SEED be the seed x₀ for the u_i computation given by Equation (10):

C    use round-to-nearest mode 

     x(0) = TWO * SEED – ONE 

     DO i = 0, n – 1 

            uc = a * x(i) + c1 

            v = uc – c1 

            x(i+1) = a * x(i) – v 

     ENDDO 

     SEED = HALF * x(n) + HALF

We now prove the following.

Theorem 2
The above code produces results that are bit-wise identical to x_i = 2u_i – 1, where u_i is the x_i given by (10).

Proof    Let au_i = I + F, where 0 < F < 1, so u_i+1 = F. We want to show x_i+1 = 2F – 1. Let u = ax_i + c1. Since a is odd, a = 2a₁ + 1 for some a₁, and r = 2(I – a₁) + c1 + 2F – 1. Note that
–2^k < ax_i < 2^k, so that 2⁵³ < u < 2⁵³ + 2^k+1. For IEEE numbers x with exponent 53, namely those x that satisfy 2⁵³ x 2⁵⁴ – 2, x is an even integer. For round-to-nearest, the computed value of u, namely uc, is the IEEE number closest to u. We claim uc = c1 + 2(I – a₁). Since uc is an even number with exponent 53, it is an IEEE number. Also, u – uc = 2F – 1 or –1 < u – uc < 1 as 0 < F < 1. Thus, uc is the closest IEEE number to u. The fused multiply–add instruction computes v and x_i+1 exactly. Thus, v = 2(I – a₁) and x_i+1 = 2F – 1.

3. Two full-period generators

We now specify a high-level description of the full-period generators. Recall that the full-cycle generator is a linear congruential generator of the form (5) with q = 2^k. We follow the recommendation of Knuth (see [2, p. 171]) and choose the value of c to be 1 or a. These two choices of c give rise to the two versions of our full-period generator. These generators have implementations that are nearly the same as the three-multiply–add generator given by Equations (6), (7), and (8).

Let x be the current iterate, in which 0 < x < 1. Then the following four multiply–add statements in (22) to (25) below describe the full-period generator:

u = fl(ax + 252),

(22)

v = u 252,

(23)

w = fl(ax – v),

(24)

x = w .

(25)

In (25), the last operation, x = w , computes the next random number. In (25), = 2^–k when c = 1, and = a2^–k when c = a. We first consider the case = 2^–k. We were not able to reduce the number of multiply–adds to three as we did in Equations (15) to (21). We would need to define c2 = 2⁵² + 2^–46, and this value cannot be represented in a double-precision word. In Equation (25), note that w < 1. This fact follows from (9). Since is the smallest random number, we are guaranteed that the new x 1. In fact, x will equal 1 only once in 2^k iterations of the generator. For k = 46, this is a very rare event. Suppose x = 1. Then (22), (23), and (24) compute w = 0 because x is an integer. However, the new x = , and we reach the conclusion that the generator is self-correcting! This feature is very important (see [3] for another example), because it is not necessary to “check” x during each iteration of the calculation. Suppose we chose
2^–46 <  < 1. Then, x in (25) could be greater than 1. In that case we would have to test and possibly correct x during every iteration:

IF (x .gt. 1) x = x – 1

(26)

A conditional statement of this sort, when present in a pipelined inner loop, can significantly degrade performance. Thus, the choice of = 2^–k is essential to making the performance of the full-cycle generator fast.

To obtain a new random number every four multiply–adds, it is necessary to unroll (22), (23), (24), and (25) by at least a factor of 2. Our numerical experiments on POWER2 showed that unrolling by 8 was necessary. The code for a linear congruential generator without unrolling is

DO i = 0, n – 1 

       u = a * x(i) + TP52 

       v = u – TP52 

       w = a * x(i) – v 

       x(i+1) = w + 

Because of pipeline delays, this code will not execute at the rate of one random number produced every four multiply–adds. It turns out that unrolling via “program vectorization” described in Section 2 does not work well. To see this, note that

xj+i = ajxi + j,

(28)

where a_j = a^j mod 2^k and _j = [(1 + a + · · · + a^j–1) mod 2^k]. The _j's are now variable! We have previously seen that the choice of = 2^–k was essential in order to maintain good performance. Any other value of gave rise to the use of conditional statements in the code, such as the statement in (26). Thus, the type of unrolling that comes from ordinary vectorization fails to provide peak performance. However, another form of unrolling works. In the present context, we call this “generalized unrolling.”

Let N = 2n be an arbitrary even integer. Suppose that we can compute a_n and _n of Equation (28) in log₂(n) steps. Then we can unroll the code by 2 in (27) as follows:

DO i = 0, n – 1 

       u = a * x(i) + TP52 

       u1 = a * x(i+n) + TP52 

       v = u – TP52 

       v1 = u1 – TP52 

       w = a * x(i) – v 

       w1 = a * x(i+n) – v1 

       x(i+1) = w +  

       x(n+i+1) = w1 +  

The code in (29) constitutes “generalized unrolling.” We have divided the vector x(0:N – 1) into two vectors x(0:n – 1) and y(0:n – 1) = x(n:N – 1) and worked on them independently. This type of unrolling gives rise to a form of single-instruction, multiple-data (SIMD) parallelism on a single processor. In [9] we also exploit this form of SIMD parallelism on POWER2 machines. The crucial point of the code in (29) for the present application is that the two vectors x and y both use the same .

We show how to compute _n in log₂(n) steps. First note that

j+i mod 1 = (aj mod 2k)(i mod 1) + j mod 1

(30)

and

aj+i mod 2k = (aj mod 2k)(ai mod 2k)

(31)

hold for all i and j. We construct a table = TBL(2, 0:46) whose ith entries TBL(1:2, i) are a²ⁱ and _2ⁱ, for 0 i 46. Note that

TBL(1, i + 1) = TBL(1, i) * TBL(1, i) mod 2k,

(32)

TBL(2, i + 1) = [TBL(1, i) + 1] * TBL(2, i) mod< 1,

(33)

along with TBL(1, 0) = a and TBL(2, 0) = 2^–46, generates the table in 46 steps. We can now, given x₀ and the precomputed table, use Equations (28), (30), and (31) to compute x_n in log₂(n) steps as follows:

   t = x(0) 

   n0 = n 

   lb = 0 

1  nn = n0/2 

   IF (n0 .gt. 2 * nn) THEN 

       u = t * TBL (1, lb) + TP52 

       v = u – TP52 

       w = t * TBL (1, lb) – v 

       t = t + TBL (2, lb) 

       IF (t .gt. 1) t = t – 1 

   ENDIF 

   lb = lb + 1 

   n0 = nn 

   IF (n0 .gt. 0) goto 1 

We now consider the second generator; this generator has = a2^–k. Consider the following code, in which TPMK = 2^–k:

DO i = 0, n – 1 

       u = x(i) + TPMK 

       v = a * u + TP52 

       w = v – TP52 

       x(i+1) = a * u – w 

Here we have 0 x_i < 1. When x_i = 1 – 2^–k, x_i+1 = 0 as u = 1. Again, this generator is self-correcting. The unrolled-by-2 version of (35) becomes

DO i = 0, n – 1 

       u = x(i) + TPMK 

       u1 = x(i+n) + TPMK 

       v = a * u + TP52 

       v1 = a * u1 + TP52 

       w = v – TP52 

       w1 = v1 – TP52 

       x(i+1) = a * u – w 

       x(n+i+1) = a * u1 – w1 

Equations (30), (31), (32), and (34) remain unchanged. Equation (33) is modified to

TBL(2, i + 1) = [TBL(1, i) + a] * TBL(2, i) mod 1.

(37)

4. Parallel implementation for HPF

The IBM XL Fortran (XLF) [4] and XL High Performance Fortran (HPF) [5] languages include a RANDOM_NUMBER intrinsic function that implements two multiplicative congruential generators. These generators are also part of PESSL, the IBM Parallel Engineering and Scientific Subroutine Library for AIX [10]. Generator 1 is of type (2), with a = 7⁵, k = 31, and modulus q = 2^k – 1. Generator 2 is of type (3), with a = 44 485 709 377 909, k = 48, and modulus q = 2^k. The generator that is to be used is selected by setting the generator argument of RANDOM_NUMBER to 1 or 2, respectively. In this section, we discuss implementation issues arising from the parallel directives in HPF for the modulus 2⁴⁸ generator. The modulus 2³¹ – 1 generator is discussed in Section 5.

HPF parallel directives
The HPF statement

CALL RANDOM_NUMBER (A)

fills the scalar, vector, or array A with random numbers. HPF provides directives (ALIGN, DISTRIBUTE, BLOCK, CYCLIC, etc.) that allow the programmer to specify what elements of A are to reside on what processor. To achieve high performance, the parallel algorithm works by computing on each processor only those random numbers resident on that processor, ultimately achieving linear speedups for large quantities of random numbers.

The N random numbers required on a particular processor j {0, · · · , P – 1} are sub-sequences of x₀, x₁, x₂, · · · , x_N from (3) of the form

xij, xij + k, xij + 2k, · · · , xij + (N – 1)k,

(38)

where k is called the stride of the sequence. If A is BLOCK-distributed, k = 1 and i_j = jN, j = 0, · · · , P – 1. If A is CYCLIC-distributed, k = P and i_j = j, j = 0, · · · , P – 1. A contiguous sequence is one with stride k = 1, and one with k > 1 is called strided. (There is also a BLOCK-CYCLIC distribution and other variations which are not discussed here for brevity, but which are handled by applying techniques similar to those presented here.)

The rest of this section shows how it is possible to compute (38) in O(N + log i_j) time, resulting in asymptotically linear speedup. The algorithm used depends on whether the required sequence is contiguous or strided.

Contiguous random-number sequences
A contiguous random-number sequence has the form

xi, xi+1, xi+2, · · · , xi+N–1.

Contiguous sequences are required in HPF on each processor when the random numbers are written to a BLOCK-distributed vector or array. Once the starting number x_i is known for a particular processor, the rest of the sequence can be directly computed by the code (14). The fast calculation of the starting numbers on each processor is dealt with next.

Strided random-number sequences
If the CYCLIC distribution directive is used in HPF with more than one processor, the sequence of random numbers resident on each processor consists of strided (k > 1) sub-sequences of the form (38). In addition, for a BLOCK distribution, the starting values x_ij, j = 0, · · · , P – 1 in (38) must be computed from some available previous x_i, for example, x₀. Both of these situations require the efficient computation of multiple steps of the recurrence (13). That is, given x_i and n, we must compute x_i+n. The code in (14) does this for several fixed n to achieve the unrolling. However, since n is not known a priori in the current context, we cannot simply precompute the value of

an = an mod q,

(39)

and must therefore devise an efficient means to compute it for general n.

We now show how to compute (13) in O(log₂n) steps, which is crucial to the asymptotically linear speedup of the parallel algorithm. Binary algorithms for exponentiation are described in [2], and one is used in Bailey's random-number generator implementation [1]. We use a binary exponentiation algorithm here, but achieve high performance through the use of the FMA instruction.

Let the binary representation of n be

n = b_k · · · b₀, b_k 0, k = log₂n, n 0.

Then

n =	k	2jbj
	j = 0

and

anxi mod 1  = [(a kj = 02jbj) mod q] xi mod 1 =  k a2j mod q xi mod 1 j=0 bj0 =  k (a2j mod q) xi mod 1 j=0 bj0

(40)

The values T_i = a²ⁱ mod q, i = 0, · · · , log₂q – 1, are precomputed and stored in a table, allowing the computation of the product (40) in k log₂n steps. The product is accumulated in a loop, as in the following pseudocode:

t = X(i) 

DO i = 0, k 

       IF (bi  0) THEN 

C         Compute t = tTi mod 1. 

          u = t * T(i) + TP52 

          v = u – TP52 

          t = t * T(i) – v 

       ENDIF 

ENDDO 

5. A modulus 2³¹ – 1 generator

We now consider efficient implementation of the generator of type (2) used in RANDOM_NUMBER. This is an extension of the work in the previous section, but is significantly more complicated since the modulus is no longer a power of 2. For the modulus 2³¹ – 1 generator, high performance is achieved by scaling the integer recurrence by 2^–31, finding a fast implementation of the scaled recurrence, and finally transforming the resulting numbers to the range (0, 1). Both the mathematical issues and the implementation issues arising from the parallel directives in HPF are discussed next.

Generator recurrence
The random-number sequence

xi (0, 1), i = 0, 1, · · ·

(42)

produced by the generator is defined by the recurrence

si+1 = asi mod q,

(43)

xi+1 =  si+1 q

(44)

where s₀ is an integer, 0 < s₀ < q, and

a = 7⁵ = 16807,

q = 2³¹ – 1.

It has a period of q – 1.

This value of q is well suited to an implementation on a machine with 32-bit floating-point arithmetic, such as the IBM S/360 [3]. Although the target architecture of XLF and XLHPF is the IBM RS/6000, which uses IEEE floating-point, this generator has been provided for compatibility with existing applications. The mathematical properties of this generator are discussed in [3].

As for the modulus 2^k generator, the algorithm used depends on whether the required sequence is contiguous or strided.

Contiguous random-number sequences
Before discussing the main contribution of this section, parallel algorithms for strided random-number sequences, it is useful to describe the sequential algorithms (due to G. Slishman⁴) on which they are based. Although the ideas here are based on earlier sections of this paper, the implementation for this generator is more involved, since q is not a power of 2.

RANDOM_NUMBER uses a sequential implementation of (43) and (44) to produce contiguous random-number sequences of the form

x_i, x_i+1, x_i+2, · · · , x_i+N–1,