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 SPECIALIZED-ARRAY: uint
8 SPECIALIZED-ARRAY: uint-4
13 CONSTANT: state-multiplier 1812433253
27 { uint-array uint-array }
28 { uint-4-array uint-4-array } ;
30 : endian-shuffle ( v -- w )
32 uchar-16{ 3 2 1 0 7 6 5 4 11 10 9 8 15 14 13 12 } vshuffle
35 : hlshift* ( v n -- w )
36 [ endian-shuffle ] dip hlshift endian-shuffle ; inline
38 : hrshift* ( v n -- w )
39 [ endian-shuffle ] dip hrshift endian-shuffle ; inline
42 dup 1 hlshift* vbitxor ; inline
45 [ 11 vrshift ] dip vbitand ; inline
53 : formula ( a b mask c d -- r )
56 [ wA ] dip vbitxor ; inline
58 GENERIC: generate ( sfmt -- )
60 M:: sfmt generate ( sfmt -- )
62 sfmt uint-4-array>> :> array
63 state n>> 2 - array nth state r1<<
64 state n>> 1 - array nth state r2<<
69 n m - >fixnum <iota> [| i |
71 i m + array nth-unsafe
72 mask state r1>> state r2>> formula :> r
74 r i array set-nth-unsafe
81 n m - 1 + + >fixnum :> i
83 m n - i + array nth-unsafe
84 mask state r1>> state r2>> formula :> r
86 r i array set-nth-unsafe
93 : period-certified? ( sfmt -- ? )
94 [ uint-4-array>> first ]
95 [ state>> parity>> ] bi vbitand odd-parity? ;
97 : first-set-bit ( x -- n )
99 dup { [ 0 > ] [ 1 bitand 0 = ] } 1&&
101 [ 1 + ] [ -1 shift ] bi*
104 : correct-period ( sfmt -- )
106 [ state>> parity>> first first-set-bit ]
107 [ uint-array>> swap '[ _ toggle-bit ] change-nth ] tri ;
109 : certify-period ( sfmt -- sfmt )
110 dup period-certified? [ dup correct-period ] unless ;
112 : <sfmt-array> ( sfmt -- uint-array uint-4-array )
114 [ n>> 4 * [1,b] uint >c-array ] [ seed>> ] bi
117 [ -30 shift ] [ ] bi bitxor
120 ] uint-array{ } accumulate-as nip
121 dup uint-4 cast-array ;
123 : <sfmt-state> ( seed n m mask parity -- sfmt )
132 : init-sfmt ( sfmt -- sfmt' )
133 dup <sfmt-array> [ >>uint-array ] [ >>uint-4-array ] bi*
134 certify-period [ generate ] keep ; inline
136 : <sfmt> ( seed n m mask parity -- sfmt )
142 : refill-sfmt? ( sfmt -- ? )
143 state>> [ index>> ] [ n>> 4 * ] bi >= ; inline
145 : next-index ( sfmt -- index )
146 state>> [ dup 1 + ] change-index drop ; inline
149 [ next-index ] [ uint-array>> ] bi nth-unsafe ; inline
154 dup refill-sfmt? [ dup generate ] when next ; inline
157 [ [ state>> ] dip >>seed drop ]
158 [ drop init-sfmt ] 2bi ;
160 : <sfmt-19937> ( seed -- sfmt )
162 uint-4{ 0xdfffffef 0xddfecb7f 0xbffaffff 0xbffffff6 }
163 uint-4{ 0x1 0x0 0x0 0x13c9e684 }
166 : default-sfmt ( -- sfmt )
167 [ random-32 ] with-secure-random <sfmt-19937> ;