Some people collect shells. Others collect stamps. I stockpile random-number generators. Who knows? Maybe someday, they'll be as valuable as baseball cards.
I still remember the first random- number generator I acquired. It gave me a case of the "chip-is-bad" fever, a disease that causes hardware engineers to repair erratic circuits by replacing every chip in sight. The generator was implemented as an interrupt service routine in the operating system for an RCA VIP computer with an 1802 microprocessor. At each clock tick, the subroutine incremented register R9--one of 16 16-bit registers in the mighty 1802. To obtain a number at random, you simply used the current register value. One day, forgetting about that subroutine and noticing R9 values changing willfully, I concluded the chip had gone bananas. Fortunately, just before pulling the processor, I remembered the interrupt, which I easily disabled by modifying a copy of the operating-system code.
Of course, it's hard to imagine a worse random-number generator than an interrupt-driven clock subroutine! Because most VIP programs were games, however, user-response times determined when programs accessed values in R9, and the resulting sequences were as random as alphabet soup. The mix of two nonrandom processes--a fast cycling register and human interaction--created apparently random sequences: an excellent example of how combination generators can produce "more random" results than their individual parts.
Last month, I examined the chi-square test for randomness. This time, I list three random-number generators you can implement in most programming languages. I'll also show you how to mix two methods to create a combination generator. If you like scrambled eggs, this column's for you.
Most first-year computer-science students learn the middle-square random-number algorithm. The method was heavily used in computing's early days, but its output tends to break down into nonrandom, repeating patterns or, worse, a stream of zeros. Better random-number algorithms are known, and the middle-square method is now considered obsolete. It provides a useful platform, however, for investigating related issues that apply to all generators.
The middle-square method is simply explained. Square an N-digit seed value, producing a product with 2N digits, then use the middle N digits of that product as a random number. Save the result for the next seed. Example 1, Algorithm #15, shows the middle-square method in pseudocode for eight-digit decimal random numbers. Because of the potential for overflow in the multiplication of Seed*Seed, you probably can't implement Algorithm #15 directly on most computers. Squaring the seed value 87654321, for example, gives 7683279989971041, which is way too large to store in a 16- or 32-bit register. Because only the middle eight digits of that value (27998997) are needed, it's possible to avoid overflow by breaking the multiplication into pieces, as shown in Listing One, MIDSQR.PAS (page 136). The program uses a technique explained in Robert Sedgewick's Algorithms in C++ (Addison-Wesley, 1992), where p*q is expressed as (104p1+p0)(104q1+q0), assuming p1=p/104, q1=q/104, p0=p%104, and q0=q%104 (with % meaning modulo).
That math and the program's logic are easier to follow if you think of 104v and v/104 as shift operations that move decimal digits left and right. If v=1234, for instance, then 104v=12340000. If v=87654321, then v/104=00008765. If v=12345678, then v%104=00005678. With these formulas, it's possible to multiply eight-digit values using 32-bit integers without concern for overflow--an important consideration because integer wraparound may create a weak link in the random-number chain. Despite that fact, even top-rated C compilers typically ignore multiplication overflows in their stock random-number generators. For more trustworthy sequences, use the expressions in Listing Two, LONGMULT.PAS (page 136), to multiply eight-digit values and extract the high, middle, or low eight digits (whichever the algorithm specifies) from the 16-digit product.
Some programs require floating-point (aka real) numbers selected at random so that 0<=r<1.0. The easy way to create random reals is to imagine a decimal place to the left of an integer value k. In MIDSQR.PAS, for example, divide RandomInt by m8 (equal to 108). If k=52478293, k/108 gives 0.52478293. You may use a similar technique with any integer random function.
The hands-down, most popular random-number generator goes by the name "linear-congruential method." You've probably seen it expressed as the formula a =ab+c mod m where a is a starting seed for each successive value, b, a constant such as 31415621, c, 0 or 1, and m, the computer's word size. With careful programming, it's also possible for m to be a power of ten, as in Example 2, Algorithm #16.
Here again, implementations of Algorithm #16 may fail to produce reliable random sequences unless you avoid integer wraparound in the multiplication of Seed*b. In this case, we need only the eight least-significant digits of the 16-digit product, obtained by breaking the multiplication into pieces, as shown in Listing Three, LINEAR.PAS (page 136), using the same formulas explained earlier.
The pseudocode in Example 2 also solves a subtle problem ignored by stock generators that simply return Seed%m. That may be a mistake because, as LINEAR.PAS demonstrates, the least-significant digits of successive seed values can be nonrandom. Run the program and examine the sequence produced by function RandomInt: 88971108, 8878069, 50915850, 46492851, 86225472, 48898113, 85623174_.The rightmost digit in each value cycles from 0 to 9--an obviously nonrandom, but expected, characteristic of the linear-congruential method. Instead of truncating the product as is commonly done, a properly written generator should extract the higher-order digits as shown in function RandomRange.
Another important consideration is the choice of values for constant b. In Seminumerical Algorithms, Donald Knuth lists several criteria for selecting a proper constant, and most generators are written to conform to his proposals. Just for fun, I examined the source code for several random-number generators based on Algorithm #16. Without naming sources, some of the constants used included 22695477, 7654321, 31415821, 3141592621, and 134775813. Not wanting to be outdone, I derived my own constant, 31415621, for LINEAR.PAS. You might want to replace that value with others and compare the results using last month's chi-square analysis.
A relatively obscure method for generating random sequences is based on the Fibonacci series, 1, 2, 3, 5, 8, 13_ in which each successive value, starting with the third, equals the sum of the preceding two. The technique uses a small list of unsigned words initialized in reverse to the first 17 Fibonacci values so that List[1]=2584, List[2]=1597_List[16]=2, List[17]=1. You also need two index variables, I=17 and J=5. After initialization, use the method shown in Example 3, Algorithm #17, to produce the next number in a random sequence. First, assign List[I]+List[J] to a temporary variable, K, and save that value in List[I]. Then, decrease I and J by one, and if either index equals 0, reset it to 17.
Unlike the classic Fibonacci sequence, which sums successive numbers, Algorithm #17 adds the nonadjacent values List[I] and List[J]. The method produces sequences with extremely long periods--according to one source, of 16,777,088 on 8-bit systems; even longer using 16- or 32-bit registers--and is well suited for implementation in assembly language. (Values are stored backwards in List so you can use a fast assembly-language instruction to test when an index is decremented to 0.) Because it uses no multiplications, the algorithm cannot suffer from overflow. Listing Four, FIBO.PAS (page 136), implements Algorithm #17 in Pascal.
There is no such thing as the perfect random-number generator--any algorithm can produce a sequence that fails one or more tests for randomness. For better results, you might try combining two or more techniques, which would produce a sequence that often seems "more random" than values that are generated by each method alone. The VIP's interrupt subroutine coupled with human interaction is a good example of a combination generator. Another, more-general technique works with any two integer random functions. Listing Five, COMBINE.PAS (page 137), demonstrates the basic idea. The program defines an array called Values that holds 100 32-bit integers (the size of a LongInt type). Each value in the array is initialized using the linear-congruential method. In order to obtain a random number, the program uses the Fibonacci method to produce an index from 1 to 100, used as an index V into the Values array. That value is returned as the next random number, and Values[V] is assigned a new number using the linear-congruential technique. The downside of combination generators is that they run more slowly and consume more memory than do other methods. For Monte Carlo testing and other critical applications, however, many experts recommend a mix of generators rather than relying on just one.
I've gotten over my case of the chip-is-bad fever, but I haven't lost my fondness for random-number generators. If you know a good one, or if you have another algorithm to share, I'd like to hear from you. Write to me in care of DDJ, or send CompuServe mail to my ID, 73627,3241.
Copyright © 1994, Dr. Dobb's JournalMiddle-square Method
Linear-congruential Method
Fibonacci Method
Combining the Classics
Your Turn
Example 1: Pseudocode for Algorithm #15 (middle-square method, obsolete).
const
m4=10000;
m8=100000000;
var
Seed: Integer
function NextRandom: Integer;
begin
Seed<-((Seed*Seed)/m4) mod m8;
NextRandom<-Seed;
end;
Example 2: Pseudocode for Algorithm #16 (linear-congruential method).
const
m4=10000;
m8=100000000;
b=31415621;
r=65536;
var
Seed: Integer
function NextRandom: Integer;
begin
Seed<-((Seed*b)+1) mod m8;
NextRandom<-((Seed/m4)*r) div m4;
end;
Example 3: Pseudocode for Algorithm #17 (Fibonacci method).
var
List: array[1..17] of WORD;
I, J: Integer;
function NextRandom: Integer;
var
K: Integer;
begin
K<-List[I]+List[J];
List[I]<-K;
I<-I--1;
J<-J--1;
if I=0 then I<-17;
if J=0 then J<-17;
NextRandom<-K;
end;
[LISTING ONE]
(* ------------------------------------------------------------------------ *(
** midsqr.pas -- Middle-Square Random Number Generator. Note: This method **
** method is considered obsolete and is listed for reference only. **
** Copyright (c) 1993 by Tom Swan. All rights reserved. **
)* ------------------------------------------------------------------------ *)
program MidSqr;
const
m4 = 10000; { 10 ^^ 4 }
m8 = 100000000; { 10 ^^ 8 }
var
Seed: LongInt;
procedure RandomInit(StartingSeed: LongInt);
begin
Seed := StartingSeed
end;
{ Square the eight-digit seed; return middle eight digits of 16-digit result. }
function RandomInt: LongInt;
var
p0, p1: LongInt;
begin
p0 := seed mod m4;
p1 := seed div m4;
Seed :=
(((p1 * p0 + p1 * p0) + (p0 * p0 div m4)) +
((p1 * p1 mod m4) * m4) ) mod m8;
RandomInt := Seed
end;
var
N: Integer;
begin
Writeln;
Writeln('Middle-Square Random Number Generator (obsolete)');
Writeln;
RandomInit(45086273);
for N := 1 to 100 do
Write(RandomInt:10);
Writeln
end.
[LISTING TWO]
(* ------------------------------------------------------------------------- *(
** longmult.pas -- Demonstrate 32-bit multiplication **
** Copyright (c) 1993 by Tom Swan. All rights reserved. **
)* ------------------------------------------------------------------------- *)
program LongMultTest;
const
m4 = 10000; { 10 ^^ 4 }
m8 = 100000000; { 10 ^^ 8 }
var
p, q, p1, p0, q1, q0: LongInt;
low, mid, high: LongInt;
begin
p := 12345678; { First value }
q := 87654321; { Second value }
p1 := p div m4;
p0 := p mod m4;
q1 := q div m4;
q0 := q mod m4;
low :=
(((p0 * q1 + p1 * q0) mod m4) * m4 + p0 * q0) mod m8;
mid :=
(((p1 * q0 + p0 * q1) +
((p0 * q0) div m4)) +
((((p1 * q1)) mod m4) * m4)) mod m8;
high :=
((((((p1 * q0 + p0 * q1) +
((p0 * q0) div m4)) +
((((p1 * q1)) mod m4) * m4)) div m4) +
((((p1 * q1)) div m4) * m4))) mod m8;
Writeln('low = ', low);
Writeln('mid = ', mid);
Writeln('high = ', high);
Writeln;
Writeln(p, ' * ', q, ' = ', high, low)
end.
[LISTING THREE]
(* ------------------------------------------------------------------------- *(
** linear.pas -- Linear Congruential Random Number Generator **
** Copyright (c) 1993 by Tom Swan. All rights reserved. **
)* ------------------------------------------------------------------------- *)
program Linear;
const
m4 = 10000; { 10 ^^ 4 }
m8 = 100000000; { 10 ^^ 8 }
b = 31415621; { constant multiplier }
var
Seed: LongInt;
I: Integer;
function LongMult(p, q: LongInt): LongInt;
var
p1, p0, q1, q0: LongInt;
begin
p1 := p div m4;
p0 := p mod m4;
q1 := q div m4;
q0 := q mod m4;
LongMult :=
(((p0 * q1 + p1 * q0) mod m4) * m4 + p0 * q0) mod m8
end;
procedure RandomInit(StartingSeed: LongInt);
begin
Seed := StartingSeed
end;
function RandomInt: LongInt;
begin
Seed := (LongMult(Seed, b) + 1) mod m8;
RandomInt := Seed
end;
function RandomRange(R: LongInt): LongInt;
begin
RandomRange := ((RandomInt div m4) * R) div m4
end;
begin
Writeln;
Writeln('Linear Congruential Random Number Generator');
Writeln;
RandomInit(1234567);
for I := 1 to 32 do
Write(RandomInt:10);
Writeln;
for I := 1 to 32 do
Write(RandomRange(32768):10);
end.
[LISTING FOUR]
(* ------------------------------------------------------------------------- *(
** fibo.pas -- Fibonacci Random Number Generator **
** Copyright (c) 1993 by Tom Swan. All rights reserved. **
)* ------------------------------------------------------------------------- *)
program Fibo;
var
List: array[1 .. 17] of WORD;
I, J: Integer;
procedure RandomInit;
var
N: Integer;
begin
List[17] := 1;
List[16] := 2;
for N := 15 downto 1 do
List[N] := List[N+1] + List[N+2];
I := 17;
J := 5
end;
function RandomInt: WORD;
var
K: WORD;
begin
K := List[I] + List[J];
List[I] := K;
I := I - 1;
J := J - 1;
if I = 0 then I := 17;
if J = 0 then J := 17;
RandomInt := K
end;
var
N: Integer;
begin
Writeln;
Writeln('Fibonacci Random Number Generator');
Writeln;
RandomInit;
for N := 1 to 100 do
Write(RandomInt:10);
Writeln
end.
[LISTING FIVE]
(* ------------------------------------------------------------------------- *(
** combine.pas - Combination Random Number Generator. Combines Algorithm #16 **
** (Linear Congruential Random Number Generator) with #17 (Fibonaccii Random **
** Number Generator) to create a table-driven combination generator. **
** Copyright (c) 1993 by Tom Swan. All rights reserved. **
)* ------------------------------------------------------------------------- *)
program Combine;
const
m4 = 10000; { 10 ^^ 4 }
m8 = 100000000; { 10 ^^ 8 }
b = 31415621; { constant multiplier }
maxIndex = 100; { array of Values }
var
List: array[1 .. 17] of WORD;
Values: array[1 .. maxIndex] of LongInt;
I, J: Integer;
Seed: LongInt;
function LongMult(p, q: LongInt): LongInt;
var
p1, p0, q1, q0: LongInt;
begin
p1 := p div m4;
p0 := p mod m4;
q1 := q div m4;
q0 := q mod m4;
LongMult :=
(((p0 * q1 + p1 * q0) mod m4) * m4 + p0 * q0) mod m8
end;
function LinearRandomInt: LongInt;
begin
Seed := (LongMult(Seed, b) + 1) mod m8;
LinearRandomInt := Seed
end;
procedure RandomInit(StartingSeed: LongInt);
var
N: Integer;
begin
Seed := StartingSeed;
List[17] := 1;
List[16] := 2;
for N := 15 downto 1 do
List[N] := List[N+1] + List[N+2];
I := 17;
J := 5;
for N := 1 to maxIndex do
Values[N] := LinearRandomInt
end;
function FiboRandomInt: WORD;
var
K: WORD;
begin
K := List[I] + List[J];
List[I] := K;
I := I - 1;
J := J - 1;
if I = 0 then I := 17;
if J = 0 then J := 17;
FiboRandomInt := K
end;
function RandomRange(R: LongInt): LongInt;
var
V: Integer;
begin
V := 1 + (FiboRandomInt mod 100);
RandomRange := ((Values[V] div m4) * R) div m4;
Values[V] := LinearRandomInt
end;
var
N: Integer;
begin
Writeln;
Writeln('Combination Random Number Generator');
Writeln;
RandomInit(7654321);
for N := 1 to 1000 do
Write(RandomRange(65536):10)
end.
End Listings