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  ***********************