\ mersenne twister 05.3.21 NAB
\ Original C-code (c) 1997, 1999
\ Makoto Matsumoto,
\ Takuji Nishimura
\ Ported to ANS-Forth
\ (c) 2000 Neil Bawd
\ Ported to Quartus Forth
\ (c) 2001 Robert Ryan

needs dblmath-ext

decimal

module mersenne

\ period parameters
624 constant MTN
397 constant MTM

(hex) 9908.B0DF 2constant Matrix-A inline
(hex) 8000.0000 2constant Upper-Mask inline
(hex) 7FFF.FFFF 2constant Lower-Mask inline

\ tempering masks
(hex) 9D2C.5680 2constant TMB inline
(hex) EFC6.0000 2constant TMC inline

\ tempering shifts
: TSU   11 drshift ; inline
: TSS   7 dlshift ; inline
: TST   15 dlshift ; inline
: TSL   18 drshift ; inline

: 'th   ( n "addr" -- &addr[n] )
   s" 2 lshift " evaluate
   bl word count evaluate
   s" + " evaluate
   ; immediate

: double-cells
   s" 2 lshift " evaluate
   ; immediate

create MT 624 double-cells allot inline
create MTI -1 , inline

public:

: SGENRAND   ( seed. -- )
   MTN 0
   DO
      2dup (hex) FFFF.0000 dand
      i 'th MT 2!
      69069. d* 1. d+
      2dup (hex) FFFF.0000 dand
      16 drshift  i 'th MT 2@ dor
      i 'th MT 2!
      69069. d* 1. d+
   LOOP 2drop
   MTN MTI ! ;

: LSGENRAND   ( &seed-array -- )
   \ length of seed-array must be
   \ at least mtn double-cells
   MTN 0
   DO
      i 'th over 2@  i 'th MT 2!
   LOOP drop
   MTN MTI ! ;

: GENRAND   ( -- u. )
   MTI @ MTN u< 0=
   IF
      MTI @ MTN = 0=
      IF
         4357. SGENRAND
      THEN

      \ 0 ... n-m-1
      MTN MTM -  0
      DO
         i 'th MT 2@ Upper-Mask dand
         i 1+ 'th MT 2@ Lower-Mask
         dand dor  ( y. )
         2dup 1 drshift 2swap ( x. y.)
         drop 1 and
         IF Matrix-A dxor THEN  ( x. )
         i MTM + 'th MT 2@ dxor
         i 'th MT 2!
      LOOP

      \ n-m ... n-2
      MTN 1-  MTN MTM -
      DO
         i 'th MT 2@ Upper-Mask dand
         i 1+ 'th MT 2@ Lower-Mask
         dand dor  ( y. )
         2dup 1 drshift 2swap ( x. y.)
         drop 1 and
         IF Matrix-A dxor THEN  ( x. )
         i MTM + MTN - 'th MT 2@ dxor
         i 'th MT 2!
      LOOP

     \ n-1, 0
     MTN 1- 'th MT 2@ Upper-Mask dand
     MT 2@ Lower-Mask dand dor  ( y. )
     2dup 1 drshift 2swap ( x. y.)
     drop 1 and
     IF Matrix-A dxor THEN  ( x. )
     MTM 1- 'th MT 2@ dxor
     MTN 1- 'th MT 2!

      0 MTI !
   THEN

   MTI @ 'th MT 2@  1 MTI +!
   2dup TSU dxor
   2dup TSS TMB dand dxor
   2dup TST TMC dand dxor
   2dup TSL dxor ;

end-module