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 combinators.short-circuit fry ;
   7 SIMD: uint
   8 SPECIALIZED-ARRAY: uint
   9 SPECIALIZED-ARRAY: uint-4
  10 IN: random.sfmt
  11
  12 <PRIVATE
  13
  14 CONSTANT: state-multiplier 1812433253
  15
  16 STRUCT: sfmt-state
  17     { seed uint }
  18     { n uint }
  19     { m uint }
  20     { index uint }
  21     { mask uint-4 }
  22     { parity uint-4 }
  23     { r1 uint-4 }
  24     { r2 uint-4 } ;
  25
  26 TUPLE: sfmt
  27     { state sfmt-state }
  28     { uint-array uint-array }
  29     { uint-4-array uint-4-array } ;
  30
  31 : wA ( w -- wA )
  32    dup 1 hlshift vbitxor ; inline
  33
  34 : wB ( w mask -- wB )
  35    [ 11 vrshift ] dip vbitand ; inline
  36
  37 : wC ( w -- wC )
  38    1 hrshift ; inline
  39
  40 : wD ( w -- wD )
  41    18 vlshift ; inline
  42
  43 : formula ( a b mask c d -- r )
  44     [ wC ] dip wD vbitxor
  45     [ wB ] dip vbitxor
  46     [ wA ] dip vbitxor ; inline
  47
  48 GENERIC: generate ( sfmt -- )
  49
  50 M:: sfmt generate ( sfmt -- )
  51     sfmt state>> :> state
  52     sfmt uint-4-array>> :> array
  53     state n>> 2 - array nth state (>>r1)
  54     state n>> 1 - array nth state (>>r2)
  55     state m>> :> m
  56     state n>> :> n
  57     state mask>> :> mask
  58
  59     n m - >fixnum iota [| i |
  60         i array nth-unsafe
  61         i m + array nth-unsafe
  62         mask state r1>> state r2>> formula :> r
  63
  64         r i array set-nth-unsafe
  65         state r2>> state (>>r1)
  66         r state (>>r2)
  67     ] each
  68
  69     ! n m - 1 + n [a,b) [
  70     m 1 - iota [
  71         n m - 1 + + >fixnum :> i
  72         i array nth-unsafe
  73         m n - i + array nth-unsafe
  74         mask state r1>> state r2>> formula :> r
  75
  76         r i array set-nth-unsafe
  77         state r2>> state (>>r1)
  78         r state (>>r2)
  79     ] each
  80
  81     0 state (>>index) ;
  82
  83 : period-certified? ( sfmt -- ? )
  84     [ uint-4-array>> first ]
  85     [ state>> parity>> ] bi vbitand odd-parity? ;
  86
  87 : first-set-bit ( x -- n )
  88     0 swap [
  89         dup { [ 0 > ] [ 1 bitand 0 = ] } 1&&
  90     ] [
  91         [ 1 + ] [ -1 shift ] bi*
  92     ] while drop ;
  93
  94 : correct-period ( sfmt -- )
  95     [ drop 0 ]
  96     [ state>> parity>> first first-set-bit ]
  97     [ uint-array>> swap '[ _ toggle-bit ] change-nth ] tri ;
  98
  99 : certify-period ( sfmt -- sfmt )
 100     dup period-certified? [ dup correct-period ] unless ;
 101
 102 : <sfmt-array> ( sfmt -- uint-array uint-4-array )
 103     state>>
 104     [ n>> 4 * 1 swap [a,b] >uint-array ] [ seed>> ] bi
 105     [
 106         [
 107             [ -30 shift ] [ ] bi bitxor
 108             state-multiplier * 32 bits
 109         ] dip + 32 bits
 110     ] uint-array{ } accumulate-as nip
 111     dup underlying>> byte-array>uint-4-array ;
 112
 113 : <sfmt-state> ( seed n m mask parity -- sfmt )
 114     sfmt-state <struct>
 115         swap >>parity
 116         swap >>mask
 117         swap >>m
 118         swap >>n
 119         swap >>seed
 120         0 >>index ;
 121
 122 : init-sfmt ( sfmt -- sfmt' )
 123     dup <sfmt-array> [ >>uint-array ] [ >>uint-4-array ] bi*
 124     certify-period [ generate ] keep ; inline
 125
 126 : <sfmt> ( seed n m mask parity -- sfmt )
 127     <sfmt-state>
 128     sfmt new
 129         swap >>state
 130         init-sfmt ; inline
 131
 132 : refill-sfmt? ( sfmt -- ? )
 133     state>> [ index>> ] [ n>> 4 * ] bi >= ; inline
 134
 135 : next-index ( sfmt -- index )
 136     state>> [ dup 1 + ] change-index drop ; inline
 137
 138 : next ( sfmt -- n )
 139     [ next-index ] [ uint-array>> ] bi nth-unsafe ; inline
 140
 141 PRIVATE>
 142
 143 M: sfmt random-32* ( sfmt -- n )
 144     dup refill-sfmt? [ dup generate ] when next ; inline
 145
 146 M: sfmt seed-random ( sfmt seed -- sfmt )
 147     [ [ state>> ] dip >>seed drop ]
 148     [ drop init-sfmt ] 2bi ;
 149
 150 : <sfmt-19937> ( seed -- sfmt )
 151     156 122
 152     uint-4{ HEX: dfffffef HEX: ddfecb7f HEX: bffaffff HEX: bffffff6 }
 153     uint-4{ HEX: 1 HEX: 0 HEX: 0 HEX: 13c9e684 }
 154     <sfmt> ; inline