|
Wil Baden 1993-09-11
Knuth's recommended Subtractive RNG from The Stanford GraphBase.
This is a port of GB_FLIP, Donald E. Knuth's random number
generator package from his The Stanford GraphBase, 1993,
ISBN 0-201-54275-7.
Words have been renamed to agree with Forth naming conventions.
These routines give identical results whether working in C or Forth.
From the book, replacing C jargon with Forth jargon.
INIT-RAND ( n -- )
seed INIT-RAND.
NEXT-RAND ( -- +n )
NEXT-RAND will return pseudo-random integers
between 0 and 2^31 - 1, inclusive.
GraphBase programs are designed to produce identical results on almost all existing computers and operating systems. An improved version of the portable subtractive method recommended in Seminumerical Algorithms, Section 3.6, is used to generate random numbers in the routines below. The period length is at least 2^55 - 1, and is in fact plausibly conjectured to be 2^85 - 2^30 for all but at most one choice of the seed value. The low-order bits of the generated numbers are just as random as the high-order bits.
A validation program is provided so installers can tell if it's working properly.
UNIF-RAND ( m -- u )
NEXT-RAND m MOD. (The bias is insignificant when m is
small, but can be serious when m is large. For example, if
m is approximately 2^32/3, the simple remainder algorithm
would give an answer less than m / 2 about 2/3 of the time.)
( Standard Forth CORE and CORE EXT with NOT )
: MOD-DIFF S" - 2147483647 AND " EVALUATE ; IMMEDIATE
CREATE A 56 CELLS ALLOT
VARIABLE FPTR
: FLIP-CYCLE ( -- +random )
A CELL+ 32 CELLS A + ( ii jj)
BEGIN
OVER >R ( R: ii)
OVER @ OVER @ MOD-DIFF ( . . diff)
R> ! ( ii jj)( R: )
SWAP CELL+ SWAP CELL+
DUP 55 CELLS A + U>
UNTIL
DROP ( ii)
A CELL+ ( ii jj)
BEGIN
OVER >R ( R: ii)
OVER @ OVER @ MOD-DIFF ( . . diff)
R> ! ( ii jj)( R: )
SWAP CELL+ SWAP CELL+
OVER 55 CELLS A + U>
UNTIL 2DROP
54 CELLS A + FPTR !
55 CELLS A + @ ( +random)
;
: INIT-RAND ( seed -- )
-1 A !
0 MOD-DIFF
DUP 55 CELLS A + !
DUP 1 ( seed prev next)
1 21 DO
DUP I CELLS A + !
MOD-DIFF ( seed next)
>R ( seed)( R: next)
DUP 1 AND IF
2/
[ HEX ] 40000000 [ DECIMAL ] +
ELSE
2/
THEN
R> ( seed next)( R: )
OVER MOD-DIFF
I CELLS A + @ SWAP ( seed prev next)
I 21 + 55 MOD I - +LOOP 2DROP DROP
FLIP-CYCLE DROP
FLIP-CYCLE DROP
FLIP-CYCLE DROP
FLIP-CYCLE DROP
FLIP-CYCLE DROP
;
: NEXT-RAND ( -- +random )
FPTR @ @ ( +random)
DUP 0< IF
DROP FLIP-CYCLE
ELSE
-1 CELLS FPTR +!
THEN ;
: UNIF-RAND ( +range -- +random )
>R ( )( R: +range)
[ HEX ] 80000000 [ DECIMAL ] ( t)
DUP 0 R@ UM/MOD DROP -
BEGIN
NEXT-RAND ( t +random)
2DUP U> NOT
WHILE
DROP ( t)
REPEAT ( t +random)
NIP ( +random)
R> MOD ( +random)( R: )
;
MARKER VALIDATE
: TEST-FLIP
-314159 INIT-RAND
NEXT-RAND
119318998 <> ABORT" Failure on the first try. "
133 0 DO NEXT-RAND DROP LOOP
[ HEX ] 55555555 [ DECIMAL ] UNIF-RAND
748103812 <> ABORT" Failure on the second try. "
;
TEST-FLIP VALIDATE
( With one's complement arithmetic, try --
: MOD-DIFF
- DUP 0< IF [ HEX ] 40000000 + 40000000 + [ DECIMAL ] THEN
;
)