ANEW --GB_FLIP-- DECIMAL \ Wil Baden 1993-09-11
\ GB_FLIP
\ 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 style with Forth style.
\ INIT-RAND ( seed -- )
\ To use the routines in this file, first call the function
\ `_seed_ INIT-RAND`.
\ NEXT-RAND ( -- random )
\ Subsequent use of `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 [2nd edition], 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.
\ UNIF-RAND ( m -- u )
\ Uniform integers. Here is a simple routine that produces a
\ uniform integer between 0 and _m_-1 inclusive, when _m_ is any
\ positive integer less than 2^31. It avoids the bias toward
\ small values that would occur if we simply calculated
\ `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.)
\
\ This routine consumes fewer than two random numbers, on the
\ average, for any fixed _m_.
\ Flip-Cycle ( -- random )
\ Compute 55 more pseudo-random numbers.
\ Standard Forth CORE and CORE EXT with `NOT`.
\ A validation program is provided so installers can tell if
\ it's working properly.
: Mod-Diff S" - 2147483647 AND " EVALUATE ; IMMEDIATE
CREATE Flip-A 56 CELLS ALLOT \ Pseudo-random values
VARIABLE Flip-PTR \ The next A value to be exported
: Flip-Cycle ( -- random )
Flip-A CELL+ Flip-A 32 CELLS + ( ii jj)
BEGIN
over >R ( R: ii)
over @ over @ Mod-Diff ( . . diff)
R> ! ( ii jj)( R: )
>R CELL+ R> CELL+
dup Flip-A 55 CELLS + U>
UNTIL
DROP ( ii)
Flip-A CELL+ ( ii jj)
BEGIN
over >R ( R: ii)
over @ over @ Mod-Diff ( . . diff)
R> ! ( ii jj)( R: )
>R CELL+ R> CELL+
over Flip-A 55 CELLS + U>
UNTIL 2DROP
Flip-A 54 CELLS + Flip-PTR !
Flip-A 55 CELLS + @ ( random)
;
: Init-Rand ( seed -- )
-1 Flip-A ! \ Initialize Flip-A.
0 Mod-Diff \ Strip off the sign.
dup Flip-A 55 CELLS + !
dup 1 ( seed prev next)
\ Compute a new _next_ value.
1 21 DO
dup Flip-A I CELLS + !
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
Flip-A I CELLS + @ SWAP ( seed prev next)
I 21 + 55 MOD I - +LOOP 2DROP DROP
\ Get the array values "warmed up".
Flip-Cycle DROP
Flip-Cycle DROP
Flip-Cycle DROP
Flip-Cycle DROP
Flip-Cycle DROP
;
: Next-Rand ( -- +random )
Flip-PTR @ @ ( +random)
dup 0< IF
DROP Flip-Cycle
ELSE
-1 CELLS Flip-PTR +!
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 ( -- flag )
-314159 Init-Rand
Next-Rand
119318998 = NOT IF CR ." Failure on the first try. "
-1 EXIT THEN
133 0 DO Next-Rand DROP LOOP
[ HEX ] 55555555 [ DECIMAL ] Unif-Rand
748103812 = NOT IF CR ." Failure on the second try. "
-2 EXIT THEN
CR ." OK, the gb_flip routines seem to work. "
0 ;
Test-Flip [IF] CR .( Test-Flip failed. ) [THEN] VALIDATE
( With one's complement arithmetic, try --
: MOD-DIFF
- dup 0< IF [ HEX ] 40000000 + 40000000 + [ DECIMAL ] THEN ;
)
\\ *********************** End of GB_FLIP ***********************