Wil Baden 1999-10-19

If you are like most of us, and need a few random numbers
for a game or such, then a single *linear congruential
sequence* (LCS) should be good enough. With a little more
effort we can get millions and millions of "random"
numbers.

According to Knuth, Chapter 3 of *The Art of Computer
Programming*, LCS was introduced by D. H. Lehmer in 1949.

Quoting Knuth,

We choose four magic integers,_

m, the modulus; 0 <m.a, the multiplier; 0 <=a<m.c, the increment; 0 <=c<m.X[0], the starting value; 0 <=X[0]<m.

The desired sequence of random numbersX[n]is then obtained by settingX[n+1]= (a*X[n]+c) modm

Knuth gives the following principles (paraphrased) for selecting those numbers.

i. The "seed" number

X[0]may be chosen arbitrarily.

ii. The numbermshould be large, say at least 2^30. Conveniently it may be the computer's word size, since that makes the computation quite efficient.

16-bit word size cannot satisfy this principle.

In Forth we can write a LCS with *m* as the word size:

VARIABLE RAND-X : RAND-NEXT ( -- u ) RAND-X @ a * c + DUP RAND-X ! ;

Another approach is to use for *m* an easily referenced large
prime within the word size. For 32-bit arithmetic,
2^31-1, which is 2147483647, is a popular choice. The
value of *c* should be 0.

: RAND-NEXT ( -- u ) RAND-X @ a m */MOD DROP DUP RAND-X ! ;

iii. If

mis a power of 2, pickaso thatamod 8 is 5. Ifmis a power of 10, chooseaso thatamod 200 is 21.

This with *c* as chosen below ensures a cycle of *m*
values that pass a test of "potency".

In *Starting Forth*, *a* is 31421.

iv. The multiplier

ashould preferably be chosen between .01*mand .99*m, and its binary or decimal digits shouldnothave a simple, regular pattern.

31421 for a 32-bit word size fails this principle -- it's too small.

Knuth recommends a "haphazard" constant like 3141592621. I call this "Pi21". It's enough digits of pi, with 21 tacked on, to fill a 32-bit word. I think Starting Forth's 31421 was chosen the same way for 16-bit words.

I remember Pi21 by "Now I want a large container of coffee-21".

v. The value of

cis immaterial whenais a good multiplier, except thatcmust have no factor in common withm.

So 1 and *a* look like good values for *c*.

This gives:

: RAND-NEXT ( -- u ) RAND-X @ 3141562621 * 1+ DUP RAND-X ! ;

or

: RAND-NEXT ( -- u ) RAND-X @ 1+ 3141562621 * DUP RAND-X ! ;

vi. When

mis the word size the least significant digits of random numbers are not very random, so decisions based on the number should always be influenced primarily by the most significant digits.

In other words, don't use ** MOD** to select a value. The

: RAND-UNIF ( u -- n ) RAND-NEXT UM* NIP ;

*Starting Forth* calls this ** CHOOSE** .

vii. The randomness in

tdimensions is only one part in thet-th root ofm.

Don't use a LCS for simulations requiring high resolution.

viii. At most

m/1000 numbers should be generated; otherwise the future will behave more and more like the past.

For a 32-bit word size a new scheme or a new multiplier should chosen at least every two million random numbers.

For a 16-bit word size a new scheme or a new multiplier should chosen every few 64 or 65 random numbers.

From this, you can see that you can't get a good single LCS for 16-bit arithmetic. This can be fixed by using multiple LCSs or other methods. Later we'll give a multiple LCS definition for 16-bit Forth.

Here is an exhibition of LCSs that have been popular. I have assigned names for ease of reference here.

Definitions that do not return 0 should not be initialized with 0.

Some of the tests that have been made on the sequences are named, but details are postponed till later.

** PI-RAND** is LCS with the multiplier based on PI. This is
given in the summary of LCSs in the last section of
chapter 3 of Knuth's

: PI-RAND-NEXT ( -- 0..4294967295 ) RAND-X @ 3141592621 * 1+ DUP RAND-X ! ;

** SJ-RAND** was proposed by Lewis, Goodman, and Miller in the

It was used in APL and IMLS subroutine library. It was also an
option in SwiftForth 1. The main reason for continued use is
that the square of *a* is less than modulus *m* and it can be
implemented without arithmetic overflow. However such small
multipliers have known defects. (16807 is 7^5.)

: SJ-RAND-NEXT ( -- 1..2147483646 ) RAND-X @ 16807 2147483647 */MOD DROP DUP RAND-NEXT ! ;

** EASY-RAND** was nominated by George Marsaglia, 1972, as a
candidate for the best multiplier, perhaps because
69069 is easy to remember.

: EASY-RAND-NEXT ( -- 0..4294967295 ) RAND-X @ 69069 * 1+ DUP RAND-NEXT ! ;

** BS-RAND** uses the best spectral primitive root for modulus
2147483647. G. S. Fishman found it by brute force in 1986.

: BS-RAND-NEXT ( -- 1..2147483646 ) RAND-X @ 62089911 2147483647 */MOD DROP DUP RAND-X ! ;

** EP-RAND** is an efficiently portable multiplier found by
Fishman in 1988.

: EP-RAND-NEXT ( -- 1..2147483646 ) RAND-X @ 48271 2147483647 */MOD DROP DUP RAND-X ! ;

** EP2-RAND** is an efficiently portable multiplier found by
L'Ecuyer in 1988.

: EP-RAND-NEXT ( -- 1..2147483338 ) RAND-X @ 40692 2147483647 248 - */MOD DROP DUP RAND-X ! ;

** SF-RAND** is the random number generator from Brodie's

`65535 AND`

: SF-RAND-NEXT ( -- 0..32767 ) RAND-X @ 31421 * 6927 + 65535 AND DUP RAND-NEXT ! ;

** C-RAND** is the default random number generator for the
Standard C Library.

** C-RAND-NEXT** is C's

`rand`

: C-RAND-NEXT ( -- 0..32767 ) RAND-X @ 1103515245 * 12345 + DUP RAND-X ! 16 RSHIFT 32767 AND ;

** RANDU** is the egregious RANDU of the '60s and '70s. It
must be initialized to odd values only. Note that the
multiplier in hex is 10003. For any three successive
values,

: RANDU-NEXT ( -- 1..2147483647 ) RAND-X @ 65539 * 2147483647 AND DUP RAND-X ! ; Implementations

In the *m*-is-word-size definitions, the *n* low-order bits
cycle in a 2^*n* period, as mentioned in (6) above.

( n is 1 ) MARKER ONCE : GO CR 17 0 DO PI-RAND-NEXT 1 AND . LOOP ; GO ONCE 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 ( n is 2 ) MARKER ONCE : GO CR 17 0 DO PI-RAND-NEXT 3 AND . LOOP ; GO ONCE 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 ( n is 3 ) MARKER ONCE : GO CR 17 0 DO PI-RAND-NEXT 7 AND . LOOP ; GO ONCE 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 Etcetera.

In *m*-is-word-size definitions, the second half of the
period is the same as the first half with sign-bit inverted.
The quarter-cyles are also the same except for the two
top bits.

Thanks to the first edition of *The Art of Computer
Programming*, I've used the function ** PI-RAND** (with various
names) in different languages for 30 years. But things
are bigger and faster than they used to be. Typically I
now use random numbers to test something a million times.
So I suggest for a simple, powerful, random number
generator for today, to combine

`EP-RAND`

`EP2-RAND`

VARIABLE RAND-X VARIABLE RAND-Y : RAND-NEXT ( -- 1..2147483647 ) RAND-X @ 48271 2147483647 */MOD DROP DUP RAND-X ! RAND-Y @ 40692 2147483399 */MOD DROP DUP RAND-Y ! - DUP 0> NOT IF 2147483647 + THEN ;

According to Knuth the period is about 74 quadrillion.

** RAND-X** and

`RAND-Y`

** RAND-X** does not have the problem with low-order bits.

1 RAND-X ! 1 RAND-Y ! ( n is 1 ) MARKER ONCE : GO CR 17 0 DO RAND-NEXT 1 AND . LOOP ; GO ONCE 1 1 0 1 0 1 1 1 0 1 1 0 0 0 1 1 0 ( n is 2 ) MARKER ONCE : GO CR 17 0 DO RAND-NEXT 3 AND . LOOP ; GO ONCE 3 2 0 3 3 3 0 3 0 0 2 1 2 0 3 2 2 ( n is 3 ) MARKER ONCE : GO CR 17 0 DO RAND-NEXT 7 AND . LOOP ; GO ONCE 5 7 5 3 5 6 0 3 5 2 5 5 3 4 3 2 4

With 16-bit arithmetic Forth, **three** LCSs should be
combined. The following is adapted from a suggestion by
L'Ecuyer in 1988.

VARIABLE RAND-X VARIABLE RAND-Y VARIABLE RAND-Z : RAND-NEXT ( -- 1..32363 ) RAND-X @ 157 32363 */MOD DROP DUP RAND-X ! RAND-Y @ 146 31727 */MOD DROP DUP RAND-Y ! - DUP 0> NOT IF 32363 + THEN RAND-Z @ 142 31657 */MOD DROP DUP RAND-Z ! - DUP 0> NOT IF 32363 + THEN ;

The exhibits here were originally presented in profane
languages -- outside of the temple of Forth. Some of those
languages do not have a double accumulator for
multiplication, and may fail with integer overflow. When
*a* times *a* is less than *m* and the remainder *m*
divided by *a* is less than the quotient then ** a m */MOD
DROP** can be replaced by the possibly more efficient:

[ m a / ] LITERAL /MOD [ m a MOD ] LITERAL * >R a * R> - DUP 0< IF m + THEN

All intermediate results will be in the official range of integers.

Many conversions from Fortran or another profane language to
Forth use a precomputed quotient and remainder. For a Forth
definition, ** */MOD DROP** can be used to simplify.

You have been warned not to use ** MOD** to select random numbers.
With a small range, the problem is neglible, but the following
shows the danger.

1 RAND-X ! 1 RAND-Y ! CREATE BINS 10 CELLS ALLOT : MIL 1000000 * ; MARKER ONCE : GO CR CR ." Distribution Using MOD " CR BINS 10 CELLS ERASE 20000 0 DO \ 9 Digit Number RAND-NEXT 1000 MIL MOD ( u) \ Into 10 Bins 100 MIL / CELLS BINS + 1 SWAP +! ( ) LOOP ; GO ONCE MARKER ONCE : GO 10 0 DO CR I CELLS BINS + @ DUP 4 .R SPACE 50 / [CHAR] * EMITS LOOP ; GO ONCE Distribution Using MOD 2805 ******************************************************** 2288 ********************************************* 1816 ************************************ 1871 ************************************* 1828 ************************************ 1807 ************************************ 1908 ************************************** 1940 ************************************** 1884 ************************************* 1853 ************************************* : RAND-UNIF ( +n -- u ) RAND-NEXT 2* UM* NIP ; MARKER ONCE : GO CR CR ." Distribution Using Multiply " CR BINS 10 CELLS ERASE 20000 0 DO \ 9 Digit Number 1000 MIL RAND-UNIF ( u) \ Into 10 Bins 100 MIL / CELLS BINS + 1 SWAP +! ( ) LOOP ; GO ONCE MARKER ONCE : GO 10 0 DO CR I CELLS BINS + @ DUP 4 .R SPACE 50 / [CHAR] * EMITS LOOP ; GO ONCE Distribution Using Multiply 1977 *************************************** 2042 **************************************** 1974 *************************************** 2024 **************************************** 2021 **************************************** 2007 **************************************** 1981 *************************************** 1991 *************************************** 2001 **************************************** 1982 ***************************************