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
14 CONSTANT: state-multiplier 1812433253
28 { uint-array uint-array }
29 { uint-4-array uint-4-array } ;
31 : endian-shuffle ( v -- w )
33 uchar-16{ 3 2 1 0 7 6 5 4 11 10 9 8 15 14 13 12 } vshuffle
36 : hlshift* ( v n -- w )
37 [ endian-shuffle ] dip hlshift endian-shuffle ; inline
39 : hrshift* ( v n -- w )
40 [ endian-shuffle ] dip hrshift endian-shuffle ; inline
43 dup 1 hlshift* vbitxor ; inline
46 [ 11 vrshift ] dip vbitand ; inline
54 : formula ( a b mask c d -- r )
57 [ wA ] dip vbitxor ; inline
59 GENERIC: generate ( sfmt -- )
61 M:: sfmt generate ( sfmt -- )
63 sfmt uint-4-array>> :> array
64 state n>> 2 - array nth state r1<<
65 state n>> 1 - array nth state r2<<
70 n m - >fixnum iota [| i |
72 i m + array nth-unsafe
73 mask state r1>> state r2>> formula :> r
75 r i array set-nth-unsafe
82 n m - 1 + + >fixnum :> i
84 m n - i + array nth-unsafe
85 mask state r1>> state r2>> formula :> r
87 r i array set-nth-unsafe
94 : period-certified? ( sfmt -- ? )
95 [ uint-4-array>> first ]
96 [ state>> parity>> ] bi vbitand odd-parity? ;
98 : first-set-bit ( x -- n )
100 dup { [ 0 > ] [ 1 bitand 0 = ] } 1&&
102 [ 1 + ] [ -1 shift ] bi*
105 : correct-period ( sfmt -- )
107 [ state>> parity>> first first-set-bit ]
108 [ uint-array>> swap '[ _ toggle-bit ] change-nth ] tri ;
110 : certify-period ( sfmt -- sfmt )
111 dup period-certified? [ dup correct-period ] unless ;
113 : <sfmt-array> ( sfmt -- uint-array uint-4-array )
115 [ n>> 4 * [1,b] uint >c-array ] [ seed>> ] bi
118 [ -30 shift ] [ ] bi bitxor
119 state-multiplier * 32 bits
121 ] uint-array{ } accumulate-as nip
122 dup uint-4 cast-array ;
124 : <sfmt-state> ( seed n m mask parity -- sfmt )
133 : init-sfmt ( sfmt -- sfmt' )
134 dup <sfmt-array> [ >>uint-array ] [ >>uint-4-array ] bi*
135 certify-period [ generate ] keep ; inline
137 : <sfmt> ( seed n m mask parity -- sfmt )
143 : refill-sfmt? ( sfmt -- ? )
144 state>> [ index>> ] [ n>> 4 * ] bi >= ; inline
146 : next-index ( sfmt -- index )
147 state>> [ dup 1 + ] change-index drop ; inline
150 [ next-index ] [ uint-array>> ] bi nth-unsafe ; inline
154 M: sfmt random-32* ( sfmt -- n )
155 dup refill-sfmt? [ dup generate ] when next ; inline
157 M: sfmt seed-random ( sfmt seed -- sfmt )
158 [ [ state>> ] dip >>seed drop ]
159 [ drop init-sfmt ] 2bi ;
161 : <sfmt-19937> ( seed -- sfmt )
163 uint-4{ HEX: dfffffef HEX: ddfecb7f HEX: bffaffff HEX: bffffff6 }
164 uint-4{ HEX: 1 HEX: 0 HEX: 0 HEX: 13c9e684 }
167 : default-sfmt ( -- sfmt )
168 [ random-32 ] with-secure-random <sfmt-19937> ;