basis/random/sfmt/sfmt.factor

   1 ! Copyright (C) 2009 Doug Coleman.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: accessors alien.c-types kernel locals math math.ranges
   4 math.bitwise math.vectors math.vectors.simd random
   5 sequences specialized-arrays sequences.private classes.struct ;
   6 SIMD: uint
   7 SPECIALIZED-ARRAY: uint
   8 SPECIALIZED-ARRAY: uint-4
   9 IN: random.sfmt
  10
  11 <PRIVATE
  12
  13 CONSTANT: state-multiplier 1812433253
  14
  15 STRUCT: sfmt-state
  16     { seed uint }
  17     { n uint }
  18     { m uint }
  19     { ix uint }
  20     { mask uint-4 }
  21     { r1 uint-4 }
  22     { r2 uint-4 } ;
  23
  24 TUPLE: sfmt
  25     { state sfmt-state }
  26     { uint-array uint-array }
  27     { uint-4-array uint-4-array } ;
  28
  29 : wA ( w -- wA )
  30    dup 1 hlshift vbitxor ; inline
  31
  32 : wB ( w mask -- wB )
  33    [ 11 vrshift ] dip vbitand ; inline
  34
  35 : wC ( w -- wC )
  36    1 hrshift ; inline
  37
  38 : wD ( w -- wD )
  39    18 vlshift ; inline
  40
  41 : formula ( a b mask c d -- r )
  42     [ wC ] dip wD vbitxor
  43     [ wB ] dip vbitxor
  44     [ wA ] dip vbitxor ; inline
  45
  46 GENERIC: generate ( sfmt -- )
  47
  48 M:: sfmt generate ( sfmt -- )
  49     sfmt state>> :> state
  50     sfmt uint-4-array>> :> array
  51     state n>> 2 - array nth state (>>r1)
  52     state n>> 1 - array nth state (>>r2)
  53     state m>> :> m
  54     state n>> :> n
  55     state mask>> :> mask
  56
  57     n m - >fixnum iota [| i |
  58         i array nth-unsafe
  59         i m + array nth-unsafe
  60         mask state r1>> state r2>> formula :> r
  61
  62         r i array set-nth-unsafe
  63         state r2>> state (>>r1)
  64         r state (>>r2)
  65     ] each
  66
  67     ! n m - 1 + n [a,b) [
  68     m 1 - iota [
  69         n m - 1 + + >fixnum :> i
  70         i array nth-unsafe
  71         m n - i + array nth-unsafe
  72         mask state r1>> state r2>> formula :> r
  73
  74         r i array set-nth-unsafe
  75         state r2>> state (>>r1)
  76         r state (>>r2)
  77     ] each
  78
  79     0 state (>>ix) ;
  80
  81 : <sfmt-array> ( sfmt -- uint-array uint-4-array )
  82     state>>
  83     [ n>> 4 * iota >uint-array ] [ seed>> ] bi
  84     [
  85         [
  86             [
  87                 [ -30 shift ] [ ] bi bitxor
  88                 state-multiplier * 32 bits
  89             ] dip +
  90         ] unless-zero 32 bits
  91     ] uint-array{ } accumulate-as nip
  92     dup underlying>> byte-array>uint-4-array ;
  93
  94 : <sfmt-state> ( seed n m mask -- sfmt )
  95     sfmt-state <struct>
  96         swap >>mask
  97         swap >>m
  98         swap >>n
  99         swap >>seed
 100         0 >>ix ;
 101
 102 : init-sfmt ( sfmt -- sfmt' )
 103     dup <sfmt-array> [ >>uint-array ] [ >>uint-4-array ] bi*
 104     [ generate ] keep ; inline
 105
 106 : <sfmt> ( seed n m mask -- sfmt )
 107     <sfmt-state>
 108     sfmt new
 109         swap >>state
 110         init-sfmt ; inline
 111
 112 : refill-sfmt? ( sfmt -- ? )
 113     state>> [ ix>> ] [ n>> 4 * ] bi >= ;
 114
 115 : next-ix ( sfmt -- ix )
 116     state>> [ dup 1 + ] change-ix drop ; inline
 117
 118 : next ( sfmt -- n )
 119     [ next-ix ] [ uint-array>> ] bi nth-unsafe ; inline
 120
 121 PRIVATE>
 122
 123 M: sfmt random-32* ( sfmt -- n )
 124     dup refill-sfmt? [ dup generate ] when next ; inline
 125
 126 M: sfmt seed-random ( sfmt seed -- sfmt )
 127     [ [ state>> ] dip >>seed drop ]
 128     [ drop init-sfmt ] 2bi ;
 129
 130 : <sfmt-19937> ( seed -- sfmt )
 131     348 330 uint-4{ HEX: BFFFFFF6 HEX: BFFAFFFF HEX: DDFECB7F HEX: DFFFFFEF }
 132     <sfmt> ; inline