ANEW --MT19937--                              \  Wil Baden  2003-01-31

\  *******************************************************************
\  *                                                                 *
\  *  Makoto Matsumoto and Takuji Nishimura  2002-01-09              *
\  *                                                                 *
\  *          Mersenne Twister 2002 Update                           *
\  *                                                                 *
\  *  http://www.math.keio.ac.jp/matumoto/MT2002/emt19937ar.html     *
\  *                                                                 *
\  *  Generate sequence of "random" numbers with a cycle of          *
\  *  2^19937-1.  That is 6002 decimal digits.                       *
\  *                                                                 *
\  *  Each block of 19937 bits is different.  Every possible block   *
\  *  (except all zeros) is generated.                               *
\  *                                                                 *
\  *  This is a tempered linear congruential sequence of binary      *
\  *  vectors of length 19937 bits.                                  *
\  *                                                                 *
\  *    A C-program for MT19937, with initialization improved        *
\  *    2002/1/26. Coded by Takuji Nishimura and Makoto Matsumoto.   *
\  *                                                                 *
\  *    Before using, initialize the state by using                  *
\  *    init_genrand(seed) or init_by_array(init_key, key_length).   *
\  *                                                                 *
\  *    C version copyright (C) 1997 - 2002, Makoto Matsumoto and    *
\  *    Takuji Nishimura, All rights reserved.                       *
\  *                                                                 *
\  *  Converted to Standard Forth by Wil Baden, 2003-01-30.          *
\  *                                                                 *
\  *   genrand_int31  genrand_real1  genrand_real3   init_by_array   * 
\  *   genrand_int32  genrand_real2  genrand_real53  init_genrand    *
\  *                                                                 *
\  *******************************************************************

\          Billions and Billions of Random Numbers

\  In 1997, Makoto Matsumoto and Takuji Nishimura developed the 
\  Mersenne Twister. The Mersenne Twister is "designed with 
\  consideration of the flaws of various existing generators," has a 
\  super-astronomical period of 2^19937 - 1, gives a sequence that is 
\  623-dimensionally equidistributed, and "has passed many stringent 
\  tests, including the die-hard test of G. Marsaglia and the load 
\  test of P. Hellekalek and S. Wegenkittl."  It is efficient in 
\  memory usage (typically using 2506 bytes of static data, and the 
\  code is quite short as well).  It generates random numbers in 
\  batches of 624 at a time, so the caching and pipelining of modern 
\  systems is exploited. It is also divide- and mod-free. 

\  REFERENCE

\    M. Matsumoto and T. Nishimura,

\    "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform
\    Pseudo-Random Number Generator",

\    ACM Transactions on Modeling and Computer Simulation,
\    Vol. 8, No. 1, January 1998, pp 3--30.

\  *******************************************************************
\  *                        Common Functions                         *
\  *******************************************************************

TRUE 0<> [IF]  \  Comment out what you already have.

    : H#  ( "hexnumber" -- u )  \  Simplified for easy porting.
        0 0 BL WORD COUNT             ( 0 0 str len)
        BASE @ >R  HEX  >NUMBER  R> BASE !
            ABORT" Not Hex " 2DROP      ( u)
        STATE @ IF  postpone LITERAL  THEN
        ; IMMEDIATE

    : 'th                             ( n "addr" -- &addr[n] )
        S" CELLS "    EVALUATE
        BL WORD COUNT EVALUATE
        S" + "        EVALUATE
        ; IMMEDIATE
        
    : THIRD  ( x y z -- x y z x )      2 PICK ;  \  Should be in CODE.
    
    : FOURTH ( w x y z -- w x y z w )  3 PICK ;  \  Should be in CODE.

\  `NOT` is equivalent to `0=` or `INVERT`.

[THEN]

\  *******************************************************************
\  *                        Period Parameters                        *
\  *******************************************************************

624 CONSTANT  N
397 CONSTANT  M
H# 9908B0DF CONSTANT  MATRIX_A    \  constant vector a
H# 80000000 CONSTANT  UPPER_MASK  \  most significant w-r bits
H# 7FFFFFFF CONSTANT  LOWER_MASK  \  least significant r bits

CREATE  MT  N CELLS ALLOT  \  the array for the state vector
CREATE  MTI  N 1+ ,  \  mti==N+1 means mt[N] is not initialized

\  *******************************************************************
\  *                  init_genrand  init_by_array                    *
\  *******************************************************************

\  init_genrand      ( s -- )
\     initialize mt[N] with a seed

\  init_by_array     ( &init_key key_length -- )
\     initialize by an array with array-length
\     init_key is the array for initializing keys
\     key_length is its length

: init_genrand      ( s -- )
    H# FFFFFFFF AND  MT !   ( )
    N 1 DO
        1812433253  \  Borosch-Niederreiter
            I 1- 'th MT @  dup 30 RSHIFT  XOR  *  I +      
        I 'th MT !
        \  See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. In the
        \  previous versions, MSBs of the seed affect only MSBs of
        \  the array mt[]. 2002/01/09 modified by Makoto Matsumoto
        I 'th MT dup @  H# FFFFFFFF AND  SWAP !
        \  for >32 bit machines
    LOOP 
    N MTI ! ;

: init_by_array     ( &init_key key_length -- )
    19650218 init_genrand
    N over > IF  N  ELSE  dup  THEN  >R  ( R: k)
    0 1                                    ( &init_key key_length j i)
    1 R> ( 1 k) DO
        dup 1- 'th MT @  
        dup 30 RSHIFT XOR  
        1664525 *  
        over 'th MT @ XOR
        >R  FOURTH THIRD CELLS + @  R> +  
        THIRD +  
        over 'th MT !  \  non-linear
        
        dup 'th MT  dup @  H# FFFFFFFF AND  SWAP !
                                     \  for WORDSIZE > 32 machines
        1+ >R  1+ R>
        dup N < NOT IF  N 1- 'th MT @  MT !  DROP 1  THEN
        >R  2dup > NOT IF  DROP 0  THEN  R>
    -1 +LOOP
    1 N 1- DO
        dup 1- 'th MT @  
        dup 30 RSHIFT XOR  
        1566083941 * 
        over 'th MT @ XOR  
        over -  
        over 'th MT !
        dup 'th MT  dup @  H# FFFFFFFF AND  SWAP !
                                     \  for WORDSIZE > 32 machines
        1+
        dup N < NOT IF  N 1- 'th MT @  MT !  DROP 1  THEN
    -1 +LOOP
    2drop 2drop                                                    ( )
    H# 80000000  MT !  \  MSB is 1; assuring non-zero initial array
    ;

\  *******************************************************************
\  *                  genrand_int32  genrand_int31                   *
\  *******************************************************************

\  genrand_int32            ( -- u )
\     generate a random number on [0,0xffffffff]-interval

\  genrand_int31             ( -- u )
\     generate a random number on [0,0x7fffffff]-interval 

    CREATE  MAG01  0 , MATRIX_A ,
    \  mag01[x] = x * MATRIX_A  for x=0,1

: genrand_int32             ( -- u )
    MTI @ N < NOT IF  \  generate N words at one time
        MTI @ N 1+ = IF  \  if init_genrand() has not been called,
            5489 init_genrand  \  a default initial seed is used
        THEN
        N M -  0  DO
            I 'th MT @ UPPER_MASK AND  I 1+ 'th MT @ LOWER_MASK AND  OR
	            dup 1 RSHIFT  SWAP  1 AND  'th MAG01 @ XOR
	            I M + 'th MT @ XOR  
	            I 'th MT !
        LOOP
        N 1-  N M -  DO
            I 'th MT @ UPPER_MASK AND  I 1+ 'th MT @ LOWER_MASK AND  OR
	            dup 1 RSHIFT  SWAP  1 AND  'th MAG01 @ XOR
	            I M + N - 'th MT @ XOR  
	            I 'th MT !
        LOOP
        N 1- 'th MT @ UPPER_MASK AND  MT @ LOWER_MASK AND  OR
        dup 1 RSHIFT  SWAP  1 AND  'th MAG01 @ XOR
            N 1- 'th MT @ XOR  N 1- 'th MT !
        0 MTI !
    THEN
    MTI @ 'th MT @  1 MTI +!
    \  Tempering
	    dup 11 RSHIFT XOR
	    dup 7 LSHIFT  H# 9D2C5680 AND XOR
	    dup 15 LSHIFT  H# EFC60000 AND XOR
	    dup 18 RSHIFT XOR 
    ;

: genrand_int31              ( -- u )
    genrand_int32  1 rshift ;

\  *******************************************************************
\  *           genrand_real1  genrand_real2  genrand_real3           *
\  *******************************************************************

\  genrand_real1             ( F: -- r )
\     generate a random number on [0,1]-real-interval 

\  genrand_real2             ( F: -- r )
\     generate a random number on [0,1)-real-interval 

\  genrand_real3             ( F: -- r )
\     generate a random number on (0,1)-real-interval 

    1.0E 4294967295.0E F/  FCONSTANT  (1.0/4294967295.0)

: genrand_real1             ( F: -- r )
    genrand_int32 0 D>F  (1.0/4294967295.0) F* ;
    \  divided by 2^32-1 

    1.0E 4294967296.0E F/  FCONSTANT  (1.0/4294967296.0)

:  genrand_real2             ( F: -- r )
     genrand_int32 0 D>F (1.0/4294967296.0) F* ;
    \  divided by 2^32 
    
: genrand_real3              ( F: -- r )
    genrand_int32 0 D>F 0.5E F+ (1.0/4294967296.0) F* ;
    \  divided by 2^32 

\  *******************************************************************
\  *                         genrand_real53                          *
\  *******************************************************************

\  genrand_real53             ( F: -- r )
\     generate a random number on [0,1) with 53-bit resolution

    1.0E 9007199254740992.0E F/  FCONSTANT  (1.0/9007199254740992.0)

: genrand_real53              ( F: -- r )
    genrand_int32 5 RSHIFT  0 D>F  
    67108864.0E  F*  genrand_int32 6 RSHIFT  0 D>F  F+
    (1.0/9007199254740992.0) F* ;

\  These real versions are due to Isaku Wada, 2002/01/09 added 

\  *******************************************************************
\  *                            MAIN DEMO                            *
\  *******************************************************************

MARKER DEMO

    CREATE INIT  H# 123 , H# 234 , H# 345 , H# 456 ,

: MAIN  CR
    INIT 4 init_by_array 
    ." 1000 outputs of genrand_int32 " CR
    1000 0 DO
        genrand_int32 12 U.R
        I 5 MOD 4 = IF  CR  THEN
    LOOP CR
    ." 1000 outputs of genrand_real2 " CR
    1000 0 DO
        genrand_real2 FE.
        I 5 MOD 4 = IF  CR  THEN
    LOOP CR ;

MAIN DEMO

\  *******************************************************************
\  *                                                                 *
\  *                    How Mersenne Twister Works                   *
\  *                                                                 *
\  *******************************************************************

\  A _Mersenne number_ is a number of the form 2^_w_-1.  A _Mersenne
\  prime_ is a Mersenne number that is prime.  A necessary but not
\  sufficient condition for a Mersenne prime is that _w_ is prime.
\  Another condition will make it sufficient.

\  As of May 2000, there are 38 known Mersenne primes.  2^19937-1
\  is the 24th and has 6002 decimal digits.  The 38th is
\  2^6972593-1 and I don't know how many decimal digits.

\  One way to define a linear congruential series is to pick numbers 
\  _m_ and _p_, where _p_ is a prime, and _m_ is a "generator" for 
\  it.  A generator for a prime _p_ is a number such that its powers 
\  modulo _p_ yield all positive numbers less than _p_.  Thus you can 
\  start with any positive number less than _p_, and continue 
\  multiplying by _m_ to get all positive numbers less than _p_. In 
\  the Mersenne Twister, instead of numbers, we are working with 
\  vectors of 19937 bits.  In this, the equivalent of a generator is 
\  a 19937 by 19937 matrix that can successively multiply any of the 
\  vectors that is not all 0 to obtain every vector that is not all 
\  0.  The arithmetic is done modulo 2. Addition is _x + y mod 2_ and 
\  multiplication is _x * y mod 2_. These are the same as _x xor y_ 
\  and _x and y_. 

\  In the program, `Matrix-A` is a value that can be used to form 
\  such a matrix.  The once-every-624-times part of the program takes 
\  this value and does calculations that are equivalent to 
\  multiplying by the full matrix. 

\  This gives 2^19937-1 combinations of 19937 bits. 

\  The matrix has the following form, but 19937 by 19937. 

\      0 1 0 0 ... 0 0 0
\      0 0 1 0 ... 0 0 0
\      0 0 0 1 ... 0 0 0
\              ...
\      0 0 0 0 ... 1 0 0
\      0 0 0 0 ... 0 1 0
\      a b c d ... x y z

\  The done-every-time part twists the output to obscure the
\  algebraic connection between successive elements.  It's
\  equivalent to multiplying by an invertible 19937 by 19937
\  matrix.

\\   //   \\   //   \\   //   \\   //   \\   //   \\   //   \\   //   

genrand_int31             ( -- u )                        MT19937.2002
   generate a random number on [0,0x7fffffff]-interval

genrand_int32             ( -- u )                        MT19937.2002
   generate a random number on [0,0xffffffff]-interval

genrand_real1             ( F: -- r )                     MT19937.2002
   generate a random number on [0,1]-real-interval

genrand_real2             ( F: -- r )                     MT19937.2002
   generate a random number on [0,1)-real-interval

genrand_real3             ( F: -- r )                     MT19937.2002
   generate a random number on (0,1)-real-interval

genrand_real53             ( F: -- r )                    MT19937.2002
   generate a random number on [0,1) with 53-bit resolution

init_by_array              ( &init_key key_length -- )    MT19937.2002
   initialize by an array with array-length
   init_key is the array for initializing keys
   key_length is its length

init_genrand               ( s -- )                       MT19937.2002
   initializes mt[N] with a seed