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 alien.data kernel locals math
   4 math.ranges math.bitwise math.vectors math.vectors.simd random
   5 sequences specialized-arrays sequences.private classes.struct
   6 combinators.short-circuit fry ;
   7 FROM: sequences => change-nth ;
   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 : endian-shuffle ( v -- w )
  32     little-endian? [
  33         uchar-16{ 3 2 1 0 7 6 5 4 11 10 9 8 15 14 13 12 } vshuffle
  34     ] unless ; inline
  35
  36 : hlshift* ( v n -- w )
  37     [ endian-shuffle ] dip hlshift endian-shuffle ; inline
  38
  39 : hrshift* ( v n -- w )
  40     [ endian-shuffle ] dip hrshift endian-shuffle ; inline
  41
  42 : wA ( w -- wA )
  43    dup 1 hlshift* vbitxor ; inline
  44
  45 : wB ( w mask -- wB )
  46    [ 11 vrshift ] dip vbitand ; inline
  47
  48 : wC ( w -- wC )
  49    1 hrshift* ; inline
  50
  51 : wD ( w -- wD )
  52    18 vlshift ; inline
  53
  54 : formula ( a b mask c d -- r )
  55     [ wC ] dip wD vbitxor
  56     [ wB ] dip vbitxor
  57     [ wA ] dip vbitxor ; inline
  58
  59 GENERIC: generate ( sfmt -- )
  60
  61 M:: sfmt generate ( sfmt -- )
  62     sfmt state>> :> state
  63     sfmt uint-4-array>> :> array
  64     state n>> 2 - array nth state r1<<
  65     state n>> 1 - array nth state r2<<
  66     state m>> :> m
  67     state n>> :> n
  68     state mask>> :> mask
  69
  70     n m - >fixnum iota [| i |
  71         i array nth-unsafe
  72         i m + array nth-unsafe
  73         mask state r1>> state r2>> formula :> r
  74
  75         r i array set-nth-unsafe
  76         state r2>> state r1<<
  77         r state r2<<
  78     ] each
  79
  80     ! n m - 1 + n [a,b) [
  81     m 1 - iota [
  82         n m - 1 + + >fixnum :> i
  83         i array nth-unsafe
  84         m n - i + array nth-unsafe
  85         mask state r1>> state r2>> formula :> r
  86
  87         r i array set-nth-unsafe
  88         state r2>> state r1<<
  89         r state r2<<
  90     ] each
  91
  92     0 state index<< ;
  93
  94 : period-certified? ( sfmt -- ? )
  95     [ uint-4-array>> first ]
  96     [ state>> parity>> ] bi vbitand odd-parity? ;
  97
  98 : first-set-bit ( x -- n )
  99     0 swap [
 100         dup { [ 0 > ] [ 1 bitand 0 = ] } 1&&
 101     ] [
 102         [ 1 + ] [ -1 shift ] bi*
 103     ] while drop ;
 104
 105 : correct-period ( sfmt -- )
 106     [ drop 0 ]
 107     [ state>> parity>> first first-set-bit ]
 108     [ uint-array>> swap '[ _ toggle-bit ] change-nth ] tri ;
 109
 110 : certify-period ( sfmt -- sfmt )
 111     dup period-certified? [ dup correct-period ] unless ;
 112
 113 : <sfmt-array> ( sfmt -- uint-array uint-4-array )
 114     state>>
 115     [ n>> 4 * [1,b] uint >c-array ] [ seed>> ] bi
 116     [
 117         [
 118             [ -30 shift ] [ ] bi bitxor
 119             state-multiplier * 32 bits
 120         ] dip + 32 bits
 121     ] uint-array{ } accumulate-as nip
 122     dup uint-4 cast-array ;
 123
 124 : <sfmt-state> ( seed n m mask parity -- sfmt )
 125     sfmt-state <struct>
 126         swap >>parity
 127         swap >>mask
 128         swap >>m
 129         swap >>n
 130         swap >>seed
 131         0 >>index ;
 132
 133 : init-sfmt ( sfmt -- sfmt' )
 134     dup <sfmt-array> [ >>uint-array ] [ >>uint-4-array ] bi*
 135     certify-period [ generate ] keep ; inline
 136
 137 : <sfmt> ( seed n m mask parity -- sfmt )
 138     <sfmt-state>
 139     sfmt new
 140         swap >>state
 141         init-sfmt ; inline
 142
 143 : refill-sfmt? ( sfmt -- ? )
 144     state>> [ index>> ] [ n>> 4 * ] bi >= ; inline
 145
 146 : next-index ( sfmt -- index )
 147     state>> [ dup 1 + ] change-index drop ; inline
 148
 149 : next ( sfmt -- n )
 150     [ next-index ] [ uint-array>> ] bi nth-unsafe ; inline
 151
 152 PRIVATE>
 153
 154 M: sfmt random-32* ( sfmt -- n )
 155     dup refill-sfmt? [ dup generate ] when next ; inline
 156
 157 M: sfmt seed-random ( sfmt seed -- sfmt )
 158     [ [ state>> ] dip >>seed drop ]
 159     [ drop init-sfmt ] 2bi ;
 160
 161 : <sfmt-19937> ( seed -- sfmt )
 162     156 122
 163     uint-4{ 0xdfffffef 0xddfecb7f 0xbffaffff 0xbffffff6 }
 164     uint-4{ 0x1 0x0 0x0 0x13c9e684 }
 165     <sfmt> ; inline
 166
 167 : default-sfmt ( -- sfmt )
 168     [ random-32 ] with-secure-random <sfmt-19937> ;