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 CD-ROM orsource code disks.
So-called “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.]
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 (floating-point) 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 (bell-shaped) 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 so-called “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 = azn-1 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 m-1, 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 32-bit 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 32-bit
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 simple-minded version described above.
Even with Schrage’s method, nevertheless, it turns out that only 23,093 of the
534,600,000 full-period 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 well-known 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 three-dimensional space;
consecutive quadruplets lie on parallel hyperplanes in four-dimensional space; and so
on. Depending on how we use the numbers, these correlations can occasionally produce
strange non-random 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 32-bit 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 32-bit 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
floating-point 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
extended-precision floating-point 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
structured-programming 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 pr events 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 extended-precision 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
32-bit 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 16-bit 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 32-bit integer by another 32-bit integer. On the
68000, the DIVS instruction can only divide a 32-bit quantity by a 16-bit 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 32-bit 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 32-bit numbers to produce a 32-bit quotient,
DIVSL.L divides two 32-bit numbers to produce both a 32-bit quotient and a 32-bit
remainder, and MULS.L multiplies two 32-bit numbers to produce a 32-bit 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 floating-point division operation. If we do not
instruct the compiler otherwise, it generates calls to the SANE package for all
floating-point 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 floating-point 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
floating-point 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 hand-coded 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 random-number routines. (Actually, I simply recycled Test.p.o
from the previous test, since the interface to the random number-routines 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
random-number 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. Addison-Wesley, 1985.
Apple Computer, Inc. Apple Numerics Manual, 2nd ed. Addison-Wesley, 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. Addison-Wesley, 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 floating-point coprocessor!!!
;-----------------------------------------------------------
; Numerical constants used in the seed-updating 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 32-bit 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 / Lewis-Goodman-Miller / 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 hand-coded
C random-number generator in a Language Systems Fortran
C program. It can be linked with any of the following
C object files:
C SANERandomNumbers.p.o (SANE's random number generator)
C PasRandomNumbers.p.o (hand-coded in Pascal)
C RandomNumbers.a.o (hand-coded 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
