Random Numbers
Volume Number:   8

Issue Number:   3

Column Tag:   Coding Efficiently

Fast Random Numbers
A random number generator that is 10 times faster
By Jon Bell, Clinton South Carolina
Note: Source code files accompanying article are located on MacTech CDROM or source code disks.
Socalled “random numbers” are an important ingredient in many applications. This article discusses some of the theory behind random number generators, and shows how to replace the Toolbox random number generators with one that produces the same results but is more than ten times faster on a Mac SE/30.
[NOTE: The material in the following two sections is summarized from the article by Park and Miller which is listed in the references at the end of this article.]
Random Numbers
All random number generators produce sequences of numbers by using a fixed set of rules to produce each number in the sequence from its predecessors. This means that the numbers are not really “random” at all, in a technical sense, because they are completely predictable, provided that you know the rules. The only way to get true “random” numbers would be to use some physical process which is inherently random, such as rolling dice or counting the number of decays per second from a radioactive source. Nevertheless, a computer can produce numbers which “look” random to a person who doesn’t know the rule by which they were generated, and which satisfy many of the common statistical tests of randomness. Such numbers are often called pseudorandom numbers.
Depending on the application, there are many kinds of random number sequences. For example, we might want real (floatingpoint) numbers which lie between some given minimum and maximum value; or we might want integers, instead. We might want “uniformly” random numbers, in which the probability of getting any single number is the same as the probability of getting any other number; or we might want random numbers which are distributed according to some pattern, such as the familiar Gaussian (bellshaped) curve. Any of these kinds of random numbers can be produced by applying a suitable transformation to a sequence of integers which is uniformly distributed between 1 and some maximum value, and so most random number generators have a random integer generator as their core.
Before actually looking at specific random number generators, it’s worth summarizing the desirable properties which they should have. First, we want to get as many different random numbers as possible, within the limits we specify. The problem is that any random number generator will begin to repeat itself eventually. If we calculate each number in the sequence based only on the previous number in the sequence, as soon as we generate a number we’ve generated before, the sequence will repeat itself from that point on. The length of the sequence before it begins to repeat is called its period. Clearly we would like the period to be as long as possible.
Second, the sequence of numbers, although completely deterministic, should “look” random enough for our purposes. At the very least, the numbers should be uniformly distributed, so that we are just as likely to get 1236 as we are to get 8490792, or 52804, or whatever. In addition, the numbers should not be correlated: if we take two (or three, or four...) numbers from the sequence according to some specified rule (say, two consecutive numbers), those numbers should not be “related” to each other in any “obvious” way. All random number generators can be made to fail this test, but some fail more “obviously” than others.
Finally, it should be possible to implement the generator correctly and (hopefully) efficiently on a wide variety of machines and in a wide variety of languages, without having to depend on “tricks” which are peculiar to a particular machine or language. No matter what the computing environment, we should be able to get repeatable results.
Many random number generators which have been implemented on various systems, or recommended in various programming textbooks, have failed to meet one or more of these criteria. Below I will describe a socalled “minimal standard” random number generator which meets all three to some extent. It is almost certainly not the “best” random number generator available, but it is “good enough” for most purposes, and can serve as a benchmark for evaluating other random number generators.
The “Minimal Standard” Generator
There are many different schemes for producing uniformly random integers on a computer. One of the simplest was devised by D.H. Lehmer in 1951. After choosing two parameters, the modulus (m) and the multiplier (a), and the first integer, z1, in the sequence, we generates successive numbers in the sequence as follows:
z2 = az1 mod m
z3 = az2 mod m
...
zn = azn1 mod m
Here “mod” means “take the remainder after dividing by,” as in Pascal and other programming languages. By the nature of the mod operation, the result always lies between 0 and m1, inclusive.
It turns out that if we let m = 231  1, there are 534,600,000 different multipliers, a, which give us a generator with the maximum possible period, namely 231  1. This is the longest period which we can attain using 32bit integers.
It looks as if we have an embarassingly rich choice of multipliers. Unfortunately, most of them cannot be used correctly on a computer with 32bit integers, because they would cause the intermediate product anz to overflow. In 1979, G. L. Schrage devised an implementation of the linear congruential generator which is less susceptible to integer overflow than the simpleminded version described above. Even with Schrage’s method, nevertheless, it turns out that only 23,093 of the 534,600,000 fullperiod multipliers can be used in a correct implementation. To decide which of these multipliers are “best”, we must study in detail the sequences they produce.
As far as randomness is concerned, linear congruential random number generators have a wellknown flaw. If we take consecutive pairs of numbers, say (z1, z2), (z3, z4), (z5, z6), etc., from a linear congruential sequence, and plot them on a graph, we find that they invariably lie on a series of parallel diagonal lines. Similarly, consecutive triplets lie on a series of parallel planes in threedimensional space; consecutive quadruplets lie on parallel hyperplanes in fourdimensional space; and so on. Depending on how we use the numbers, these correlations can occasionally produce strange nonrandom results. The spacing between the planes varies according to the choice of a and m, and should be made as small as possible.
In 1969, Lewis, Goodman and Miller designed a 32bit random number generator for the IBM System/360. They decided to use a linear congruential generator with m = 231  1 = 2147483647 and a = 75 = 16807. Since then, this generator has been implemented on a variety of systems, and has been dubbed a “minimal standard” random number generator.
The “minimal standard” generator does not give the smallest hyperplane spacing, though. That honor appears to be shared by the two generators which use the multipliers a = 48271 and a = 69621. Nevertheless, I will use the “minimal standard” in my implementation. The program listings indicate which parameters you need to change if you want to experiment with the “better” multipliers.
Before looking at how to implement the minimal standard generator, though, let’s look at what the Macintosh Toolbox already provides us.
The QuickDraw Random Number Generator
The Macintosh Toolbox actually contains two random number generators. The first one is part of QuickDraw. Its seed is a QuickDraw global variable named randSeed, a longint which is initialized to 1 by InitGraf at the start of a program. The function Random updates the seed and returns an integer in the range from 32767 to 32767. I know nothing about the properties of this generator, but clearly its period must fall far short of the period of a decent 32bit generator, and so I will not consider it further.
The SANE Random Number Generator
The Toolbox’s second random number generator is part of the SANE floatingpoint arithmetic package:
function RandomX (var x : extended) : extended;
The parameter x is a variable which keeps track of the sequence of random integers (which I called z in the discussion above). It is commonly called the seed for the random number generator. At the beginning of a program, we must initialize this variable to a suitable integer value in the range [1, 231  2], and make sure that we preserve its value between calls to RandomX. Whenever we call RandomX, it updates the seed to a new value, and, in addition, returns that value as its result. Note that the seed and result are always integer numbers, even though they are stored in extendedprecision floatingpoint variables.
According to the Apple Numerics Manual, RandomX is in fact a “minimal standard” random number generator, as described above.
Having to preserve the value of the seed between calls to RandomX is a nuisance if you use RandomX in several different procedures whose scopes do not all overlap. One quick and easy solution is to make the seed a global variable, which, however, offends structuredprogramming purists like me. A better solution (in MPW or Think Pascal, anyway) is to set up a unit which declares the seed as a global variable in the implementation section. This effectively “hides” the seed from the rest of your program, and prevents you from trashing it accidentally in some other routine.
Listing 2 shows an MPW Pascal unit, SANERandomNumbers.p, which carries out this idea. The procedure InitRandomSeed initializes the seed to a specified value, and the parameterless function RandomSeed returns its value for inspection. These are the only two ways that another routine can gain access to the seed. The parameterless function RandomReal updates the seed and returns an extendedprecision real number, uniformly distributed in the range (0,1). The function RandomInteger updates the seed and returns a longint, uniformly distributed in the range [1, max], where max is specified as a parameter.
Listing 1 shows a simple MPW tool which I used to test this generator, and the other ones to be described in this article. It initializes the seed to 1, displays the first ten random numbers in the sequence, displays the value of the seed after 10,000 iterations, and finally counts how many random real numbers are generated per second. (Actually, it counts for one minute, then divides by 60, to get a reliable average.) The 10,000th seed tests whether the generator has been implemented correctly; the minimal standard generator must produce the value 1043618065. If this value does not appear, either the generator is not the minimal standard one, or else it was not implemented correctly.
Likning this program with SANERandomNumbers gave the following results:
The first ten random numbers are:
0.000007826369259426
0.131537788143166242
0.755605322195033227
0.458650131923449287
0.532767237412169221
0.218959186328090348
0.047044616214486126
0.678864716868318951
0.679296405836612175
0.934692895940827623
After 10000 iterations,
the random seed is 1043618065.
In one second, 1381 random numbers were generated.
The “Minimal Standard” Generator in Pascal
Having to go through the Toolbox trap dispatcher slows down the SANE random number generator significantly, which is a drawback if we need to generate lots of random numbers quickly. Since we know the algorithm, why not simply implement it directly in Pascal and bypass the trap dispatcher completely? This would seem to be trivial: simply replace calls to RandomX with the statement
seed := (a * seed) mod m;
The problem is that the product a * seed can easily overflow the limits of a 32bit integer variable, causing incorrect results from that point on. G. L. Schrage came to the rescue in 1979 with an algorithm which splits the seed into its high and low 16bit halves. Schrage’s algorithm is implemented in the unit PasRandomNumbers.p, given in Listing 3. It is gratifyingly quicker than the SANE implementation, generating 3049 random reals per second versus 1381, a speedup of 2.2x. (These timings, and all those that follow, were measured on an unmodified Mac SE/30, using the program given in Listing 1.)
Schrage’s method has two hidden “bottlenecks” when implemented on a 68000 processor (Mac Plus and SE). The first one involves the operations seed div q and seed mod q, which require dividing a 32bit integer by another 32bit integer. On the 68000, the DIVS instruction can only divide a 32bit quantity by a 16bit quantity! Therefore a compiler must perform these operations in software. MPW Pascal inserts calls to the routines %I_DIV4 and %I_MOD4, which are located in the library Paslib.o. Similarly, the products a * lo and r * hi must be carried out as 32bit multiplications, via the routine %I_MUL4 in Paslib.o.
If we tell the compiler to generate 68020/68030 code (in MPW Pascal, by using the mc68020 option), these operations can each be performed with a single machine instruction. DIVS.L divides two 32bit numbers to produce a 32bit quotient, DIVSL.L divides two 32bit numbers to produce both a 32bit quotient and a 32bit remainder, and MULS.L multiplies two 32bit numbers to produce a 32bit result. Recompiling PasRandomNumbers.p with the mc68020 option allows us to produce 4491 random numbers per second, a speed increase of 1.5x over the preceding version.
In generating random real numbers (function RandomReal), there is another bottleneck in the final step, which is a floatingpoint division operation. If we do not instruct the compiler otherwise, it generates calls to the SANE package for all floatingpoint operations, which allows the code to run on all Macs, but is rather slow. If we are running on a Mac with a 68020 or 68030 CPU, we probably have a 68881 or 68882 floatingpoint coprocessor (FPU) available as well. (Exceptions include the Mac LC, and some accelerated SEs and Pluses.) Therefore, if we’re using the mc68020 option, we may as well use mc68881 as well, to tell the compiler to use the FPU for floatingpoint arithmetic. Recompiling PasRandomNumbers.p yet again, with both options enabled, allows us to produce 15067 random numbers per second, a speed of increase of 3.6x over the previous version. [Note: we must also recompile the main program with the mc68881 option, because the FPU expects extended variables to be stored in 10 bytes, whereas SANE expects extended variables to be stored in 8 bytes.]
The “Minimal Standard” Generatorin
Assembly Language
Examining the compiler’s output for the version with 68020/30 and FPU code, we can (as usual) spot things which could be improved in handcoded assembly language. In procedure UpdateSeed, there are two things in particular. First, the compiler uses separate operations for seed div q and seed mod q (DIVS.L and TDIVS.L respectively), even though TDIVS.L produces both the quotient and remainder, in separate registers. Second, there is quite a bit of “unnecessary” shuffling of data between registers.
The obvious solution here is to recode the entire unit in assembly language, which is not a very difficult task. The results are shown in Listing 4. To use it with the existing test program, we follow these steps:
• Assemble RandomNumbers.a to produce an object file RandomNumbers.a.o;
• Create an interface file, RandomNumbers.p, by removing the implementation section from PasRandomNumbers.p;
• Compile Test.p, using the interface file RandomNumbers.p to define the interfaces to the randomnumber routines. (Actually, I simply recycled Test.p.o from the previous test, since the interface to the random numberroutines hadn’t changed);
• Link Test.p.o with RandomNumbers.a.o.
Results: 18758 random numbers per second, a speed increase of 1.2x over the previous version, or a net increase of 13.6x over the original SANE implementation!
For Fortran Freaks Only
Since many scientific applications are still written in Fortran, I ought to point out that any of the units described above can be easily used in Language Systems Fortran under MPW. Listing 5 shows a Fortran version of the test program given in Listing 1. LS Fortran compiles this to a standalone application rather than a MPW tool, but otherwise the two programs work the same way.
The most important thing to note here is that standard Fortran passes all subroutine arguments by reference, like Pascal’s var parameters. However, two of our randomnumber routines (InitRandomSeed and RandomInteger) expect their arguments to be passed by value. LS Fortran has an extension which enables this, in the form of the %VAL function:
CALL INITRANDOMSEED (%VAL(1))
This passes the value 1 directly to the subroutine, rather than creating a temporary variable, storing 1 in it, and passing the address of this temporary variable.
References
Apple Computer, Inc. Inside Macintosh, Volume I. AddisonWesley, 1985.
Apple Computer, Inc. Apple Numerics Manual, 2nd ed. AddisonWesley, 1988.
Stephen K. Park and Keith W. Miller. “Random Number Generators: Good Ones Are Hard to Find”, Communications of the ACM, vol. 31, p. 1192 (October 1988).
Steve Williams. 68030 Assembly Language Reference. AddisonWesley, 1989.
Listing 1: Test.p
program TestRandomNumbers (input, output);
{ This program tests the random number generator imple
mented by the unit RandomNumbers. If the generator is
the "minimal standard generator" (multiplicative linear
congruential algorithm with modulus 2147483647 and mul
tiplier 16807), and the initial seed is 1, the value of
the seed after 10000 iterations should be 1043618065.
If this value is not obtained, then the generator either
is not a "minimal standard" generator, or has failed due
to integer overflow at some point.
Under the Macintosh Programmer's Workshop, this program
can be compiled and run as an MPW Tool. Under other
development systems, it may need some modification.
Jon Bell
Dept. of Physics & Comp. Sci.
Presbyterian College
Clinton SC 29325
CIS #70441,353 }
USES RandomNumbers, Events;
var
k : longint;
x : extended;
stopTime : longint;
begin
{ First verify that the generator is working correctly. }
InitRandomSeed (1);
writeln;
writeln ('The first ten random numbers are:');
writeln;
for k := 1 to 10 do
writeln (RandomReal:20:18);
for k := 11 to 10000 do
x := RandomReal;
writeln;
writeln ('After 10000 iterations,');
writeln ('the random seed is ', RandomSeed:1, '.');
writeln;
{ Now let's find out how fast the generator is. }
k := 0;
stopTime := TickCount + 3600;
while TickCount <= stopTime do
begin
x := RandomReal;
k := k + 1
end;
writeln ('In one second, ', trunc(k/60):1,
' random numbers were generated.');
writeln
end.
Listing 2: SANERandomNumbers.p
UNIT RandomNumbers;
{ This unit provides a convenient interface to the SANE
random number generator. }
INTERFACE
USES Sane;
procedure InitRandomSeed (newSeed : longint);
{ Initializes the random number seed to "newSeed". You
must call this procedure once, at the beginning of your
program, before you use any of the following functions.
As far as randomness is concerned, it makes no difference
what value you use for "newSeed", so long as it isn't
zero. Using different seeds merely gives you different
sequences of "random" numbers. Using the same seed each
time you run the program gives you the same sequence of
"random" numbers each time, which may be useful for
debugging. }
function RandomSeed : longint;
{ Returns the current value of the random number seed. }
function RandomReal : extended;
{ Returns a real number, "randomly" selected from the
interval [0,1]. }
function RandomInteger (max : longint) : longint;
{ Returns an integer, "randomly" selected from the
interval [1,max]. }
IMPLEMENTATION
const
m = 2147483647;
var
seed : extended;
procedure InitRandomSeed (newSeed : longint);
begin
seed := newSeed
end;
function RandomSeed : longint;
begin
RandomSeed := round(seed)
end;
function RandomReal : extended;
begin
RandomReal := RandomX (seed) / m;
end;
function RandomInteger (max : longint) : longint;
var
ignore : extended;
begin
ignore := RandomX (seed);
RandomInteger := (round(seed) mod max) + 1
end;
END.
Listing 3: PasRandomNumbers.p
UNIT RandomNumbers;
{ This unit implements a "minimal standard" random number
generator. }
INTERFACE
procedure InitRandomSeed (newSeed : longint);
function RandomSeed : longint;
function RandomReal : extended;
function RandomInteger (max : longint) : longint;
IMPLEMENTATION
const
a = 16807;
m = 2147483647;
q = 127773; { m div a }
r = 2836; { m mod a }
{ NOTE: Other possible values for these constants, which
may give "better" results, are:
a = 48271, q = 44488, r = 3399 or
a = 69621, q = 30845, r = 23902,
keeping m = 2147483647. }
var
seed : longint;
procedure InitRandomSeed (newSeed : longint);
begin
seed := newSeed
end;
function RandomSeed : longint;
begin
RandomSeed := seed
end;
{private} procedure UpdateSeed;
var
lo, hi, test : longint;
begin
hi := seed div q;
lo := seed mod q;
test := a * lo  r * hi;
if test > 0 then
seed := test
else
seed := test + m
end;
function RandomReal : extended;
begin
UpdateSeed;
RandomReal := seed / m
end;
function RandomInteger (max : longint) : longint;
begin
UpdateSeed;
RandomInteger := (seed mod max) + 1
end;
END.
Listing 4: RandomNumbers.a
MACHINE MC68020
MC68881
;
; WARNING: These routines require a 68020 or 68030 CPU,
; and a 68881 or 68882 floatingpoint coprocessor!!!
;
; Numerical constants used in the seedupdating algorithm.
a EQU 16807 ; multiplier
m EQU 2147483647 ; modulus
q EQU 127773 ; m div a
r EQU 2836 ; m mod a
;
; The random number seed (a 32bit integer) is a global
; variable whose value is preserved between calls to
; the random number functions.
Seed DS.L 1
InitRandomSeed PROC EXPORT
;                             
; Initializes the random number seed to a specified value.
;
; Pascal interface:
; procedure InitRandomSeed (newSeed : longint);
;                             
newSeed EQU 8
ParamSize EQU 4
LocalSize EQU 0
LINK A6, #LocalSize
MOVE.L newSeed(A6), seed(A5)
UNLK A6
RTD #ParamSize
DC.B 'INITRANDOMSEED'
RandomSeed FUNC EXPORT
;                             
; Returns the current value of the random number seed.
;
; Pascal interface:
; function RandomSeed : longint;
;                             
result EQU 8
LocalSize EQU 0
LINK A6, #LocalSize
MOVE.L Seed(A5), result(A6)
UNLK A6
RTS
DC.B 'RANDOMSEED'
UpdateSeed PROC
;                             
; Applies the Lehmer / LewisGoodmanMiller / Schrage
; algorithm to update the random number seed. This
; procedure is not to be used outside this unit, therefore
; it is not EXPORTed.
;
; It stores the new seed in the global variable "Seed",
; and leaves a copy in register D1 for use by the calling
; routine.
;                             
hi EQU D0
lo EQU D1
MOVE.L Seed(A5), hi
TDIVS.L #q, lo:hi
MULS.L #a, lo
MULS.L #r, hi
SUB.L hi, lo
BGT.S StoreNewSeed
ADD.L #m, lo ; if lo <= 0
StoreNewSeed
MOVE.L lo, Seed(A5)
RTS
DC.B 'UPDATESEED'
RandomReal FUNC EXPORT
;                             
; Updates the random number seed and returns a real number
; in the range (0,1).
;
; Pascal interface:
; function RandomReal : extended;
;                             
quotient EQU FP0
newSeed EQU D1
result EQU 8
LocalSize EQU 0
LINK A6, #LocalSize
JSR UpdateSeed
FMOVE.L newSeed, quotient
FDIV.L #m, quotient
FMOVE.X quotient, ([result,A6])
UNLK A6
RTS
DC.B 'RANDOMREAL'
RandomInteger FUNC EXPORT
;                             
; Updates the random number seed and returns a longint
; in the range [1,max].
;
; Pascal interface:
; function RandomInteger (max : longint) : longint;
;                             
dividend EQU D1
divisor EQU D0
result EQU 8
max EQU 12
LocalSize EQU 0
ParamSize EQU 4
LINK A6, #LocalSize
JSR UpdateSeed
; The new seed is now in "dividend."
TDIVS.L max(A6), divisor:dividend
; "Divisor" now contains the remainder.
ADDQ.L #1, divisor
MOVE.L divisor, result(A6)
UNLK A6
RTD #ParamSize
DC.B 'RANDOMINTEGER'
END
Listing 5: TestF77.f
C
C This program demonstrates the use of a handcoded
C randomnumber generator in a Language Systems Fortran
C program. It can be linked with any of the following
C object files:
C
C SANERandomNumbers.p.o (SANE's random number generator)
C PasRandomNumbers.p.o (handcoded in Pascal)
C RandomNumbers.a.o (handcoded in assembly language)
C
!!M Inlines.f
IMPLICIT NONE
INTEGER K, STOPTIME, RANDOMSEED
EXTENDED X, RANDOMREAL
EXTERNAL RANDOMREAL, RANDOMSEED
CALL INITRANDOMSEED (%VAL(1))
PRINT *
PRINT *, 'The first ten random numbers are:'
PRINT *
DO K = 1, 10
PRINT '(1X, F20.18)', RANDOMREAL()
END DO
DO K = 11, 10000
X = RANDOMREAL()
END DO
PRINT *
PRINT *, 'After 10000 iterations,'
PRINT *, 'the random seed is ', RANDOMSEED(), '.'
PRINT *
K = 0
STOPTIME = TICKCOUNT() + 3600
DO WHILE (TICKCOUNT .LE. STOPTIME)
X = RANDOMREAL()
K = K + 1
END DO
PRINT *, 'In one second, ', K/60,
& ' random numbers were generated.'
PRINT *
END