A high-speed normal random number generator using table look-up

Abstract

Random number generators (RNGs) output uniformly distributed numbers. To get other distributions, RNG output may be transformed, which usually requires much more computation than the RNG itself. A highly efficient alternative is the random selection of numbers from a core array containing numbers with the desired distribution (Lehman & Bailey, 1968). Since most distributions of interest are continuous with infinite bounds, the speed of look-up is bought at the cost of discontinuities and finite bounds. Increasing the dimension, K, of the array reduces these; but beyond moderate K, efficiency is lost to core costs and/or array loading costs. For a given K, the discontinuities can be minimized by cutting the distribution into successive vertical slices, each with area l/K, and loading the array with the middle abscissa value under each slice. This has the further advantage of reducing discontinuities in denser areas of the distribution at the cost of larger ones among the less frequent numbers. For the unit normal curve (P= 0, a = 1), K = (4,096) appears to be a good array size. With it, the difference between adjacent numbers is .0000610 at /1, increasing to .0061108 at ±3a. The largest magnitude array number is 3.48738. Doubling K halves the step sizes (and halving K doubles them), but there is little gain in maximum magnitude number. For the normal curve, summing two or three randomly selected array numbers is more efficient than increasing K. It is most efficient to compose the table only once and place it in peripheral storage, from which it can be rapidly read into core whenever the RNG is needed. The following program does this for the unit normal curve, using an inverse integral approximation from the Handbook of Mathematical Functions (1964). To use it, one must first choose K, which must be a power of in order to preserve the uniform distribution of the random subscripts used to access the array. The user must also: (1) dimension TABLE to the selected K, (2) set KHALF to K/2, and (3) store the calculated values. C ADD STATEMENT HERE TO GIVE TABLE DIMENSION K C SET KHALf TO K/2 CALCK = KHALF J = 2 KHALF DENOM = ~. 0 CALCK RNUM = 1.0 DO 20 I = I, KHALF AREA = RNUM/DENDM RNUM = RNUM + 2.0 TSQR = ALDGCI.O/(AREAAREA)) T = SQRT(TSQR) . TABLE(I) = T (2.515517 + 0.802853T + 0.010328TSQR)1 (1.0 + (1.~32788 + 0.189269T + 0.001308TSQR)T) TABLEeJ) = TABLE(l)(-l.O) 20 J = J I C INSERT STATEMENTS HERE TO STORE TABLE END