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 ;
   6 IN: random.sfmt
   7
   8 SIMD: uint
   9 SPECIALIZED-ARRAY: uint
  10 SPECIALIZED-ARRAY: uint-4
  11
  12 CONSTANT: SFMT_N 156
  13 CONSTANT: SFMT_M 122
  14
  15 CONSTANT: state-multiplier 1812433253
  16
  17 TUPLE: sfmt
  18 sl1 sl2 sr1 sr2 mask parity
  19 { seed integer } { n fixnum } { m fixnum }
  20 { m-n fixnum } { ix fixnum } { state uint-4-array } ;
  21
  22 : init-state ( sfmt -- sfmt' )
  23     dup [ n>> 4 * iota >uint-array ] [ seed>> ] bi
  24     [
  25         [
  26             [
  27                 [ -30 shift ] [ ] bi bitxor
  28                 state-multiplier * 32 bits
  29             ] dip +
  30         ] unless-zero 32 bits
  31     ] uint-array{ } accumulate-as nip underlying>> byte-array>uint-4-array
  32     >>state ;
  33
  34 : wA ( w -- wA )
  35    dup 1 hlshift vbitxor ; inline
  36
  37 : wB ( w mask -- wB )
  38    [ 11 vrshift ] dip vbitand ; inline
  39
  40 : wC ( w -- wC )
  41    1 hrshift ; inline
  42
  43 : wD ( w -- wD )
  44    18 vlshift ; inline
  45
  46 : formula ( a b mask c d -- r )
  47     [ wC ] dip wD vbitxor
  48     [ wB ] dip vbitxor
  49     [ wA ] dip vbitxor ; inline
  50
  51 GENERIC: generate ( sfmt -- sfmt' )
  52
  53 M:: sfmt generate ( sfmt -- sfmt' )
  54     sfmt n>> 2 - sfmt state>> nth :> r1!
  55     sfmt n>> 1 - sfmt state>> nth :> r2!
  56     0 :> r!
  57     0 :> i!
  58     sfmt m>> :> m
  59     sfmt n>> :> n
  60     sfmt m-n>> :> m-n
  61     sfmt mask>> :> mask
  62     sfmt state>> :> state
  63
  64     n m - iota [
  65         i!
  66         i state nth
  67         i m + state nth
  68         mask r1 r2 formula r!
  69
  70         r i state set-nth
  71         r2 r1!
  72         r r2!
  73     ] each
  74
  75     n m - 1 + n [a,b) [
  76         i!
  77         i state nth
  78         m-n i + state nth
  79         mask r1 r2 formula r!
  80
  81         r i state set-nth
  82         r2 r1!
  83         r r2!
  84     ] each
  85
  86     sfmt 0 >>ix ;
  87
  88 : <sfmt> ( seed n m sl1 sl2 sr1 sr2 mask parity -- sfmt )
  89     sfmt new
  90         swap >>parity
  91         swap >>mask
  92         swap >>sr2
  93         swap >>sr1
  94         swap >>sl2
  95         swap >>sl1
  96         swap >>m
  97         swap >>n
  98         swap 32 bits >>seed
  99         dup [ m>> ] [ n>> ] bi - >>m-n
 100         0 >>ix
 101         init-state
 102         generate ;
 103
 104 : <sfmt-19937> ( seed -- sfmt )
 105     348 330 5 3 9 3
 106     uint-4{ HEX: BFFFFFF6 HEX: BFFAFFFF HEX: DDFECB7F HEX: DFFFFFEF }
 107     uint-4{ HEX: ecc1327a HEX: a3ac4000 HEX: 0 HEX: 1 }
 108     <sfmt> ;
 109
 110 : refill-sfmt? ( sfmt -- ? )
 111     [ ix>> ] [ n>> 4 * ] bi >= ;
 112
 113 : nth-sfmt ( sfmt -- n )
 114     [ ix>> 4 /mod swap ]
 115     [ state>> nth nth ]
 116     [ [ 1 + ] change-ix drop ] tri ; inline
 117
 118 M: sfmt random-32* ( sfmt -- n )
 119     dup refill-sfmt? [ generate ] when
 120     nth-sfmt ;