diff --git a/basis/random/sfmt/authors.txt b/basis/random/sfmt/authors.txt new file mode 100644 index 0000000000..b4bd0e7b35 --- /dev/null +++ b/basis/random/sfmt/authors.txt @@ -0,0 +1 @@ +Doug Coleman \ No newline at end of file diff --git a/basis/random/sfmt/sfmt-tests.factor b/basis/random/sfmt/sfmt-tests.factor new file mode 100644 index 0000000000..53a134dbb8 --- /dev/null +++ b/basis/random/sfmt/sfmt-tests.factor @@ -0,0 +1,6 @@ +! Copyright (C) 2009 Doug Coleman. +! See http://factorcode.org/license.txt for BSD license. +USING: tools.test random.sfmt ; +IN: random.sfmt.tests + +[ ] [ 100 drop ] unit-test diff --git a/basis/random/sfmt/sfmt.factor b/basis/random/sfmt/sfmt.factor new file mode 100644 index 0000000000..a8fab9f3ab --- /dev/null +++ b/basis/random/sfmt/sfmt.factor @@ -0,0 +1,118 @@ +! Copyright (C) 2009 Doug Coleman. +! See http://factorcode.org/license.txt for BSD license. +USING: accessors alien.c-types kernel locals math math.ranges +math.bitwise math.vectors math.vectors.simd random +sequences specialized-arrays ; +IN: random.sfmt + +SIMD: uint +SPECIALIZED-ARRAY: uint +SPECIALIZED-ARRAY: uint-4 + +CONSTANT: SFMT_N 156 +CONSTANT: SFMT_M 122 + +CONSTANT: state-multiplier 1812433253 + +TUPLE: sfmt +sl1 sl2 sr1 sr2 mask parity +{ seed integer } { n fixnum } { m fixnum } +{ ix fixnum } { state uint-4-array } ; + +: init-state ( sfmt -- sfmt' ) + dup [ n>> 4 * iota >uint-array ] [ seed>> ] bi + [ + [ + [ + [ -30 shift ] [ ] bi bitxor + state-multiplier * 32 bits + ] dip + + ] unless-zero 32 bits + ] uint-array{ } accumulate-as nip underlying>> byte-array>uint-4-array + >>state ; + +: wA ( w -- wA ) + dup 1 hlshift vbitxor ; inline + +: wB ( w mask -- wB ) + [ 11 vrshift ] dip vbitand ; inline + +: wC ( w -- wC ) + 1 hrshift ; inline + +: wD ( w -- wD ) + 18 vlshift ; inline + +: formula ( a b mask c d -- r ) + [ wC ] dip wD vbitxor + [ wB ] dip vbitxor + [ wA ] dip vbitxor ; inline + +GENERIC: generate ( sfmt -- sfmt' ) + +M:: sfmt generate ( sfmt -- sfmt' ) + sfmt n>> 2 - sfmt state>> nth :> r1! + sfmt n>> 1 - sfmt state>> nth :> r2! + 0 :> r! + 0 :> i! + sfmt n>> :> n + sfmt m>> :> m + sfmt mask>> :> mask + sfmt state>> :> state + + n m - iota [ + i! + i state nth + i m + state nth + mask r1 r2 formula r! + + r i state set-nth + r2 r1! + r r2! + ] each + + i 1 + n [a,b) [ + i! + i state nth + m n - i + state nth + mask r1 r2 formula r! + + r i state set-nth + r2 r1! + r r2! + ] each + + sfmt 0 >>ix ; + +: ( seed n m sl1 sl2 sr1 sr2 mask parity -- sfmt ) + sfmt new + swap >>parity + swap >>mask + swap >>sr2 + swap >>sr1 + swap >>sl2 + swap >>sl1 + swap >>m + swap >>n + swap 32 bits >>seed + 0 >>ix + init-state + generate ; + +: ( seed -- sfmt ) + 348 330 5 3 9 3 + uint-4{ HEX: BFFFFFF6 HEX: BFFAFFFF HEX: DDFECB7F HEX: DFFFFFEF } + uint-4{ HEX: ecc1327a HEX: a3ac4000 HEX: 0 HEX: 1 } + ; + +: refill-sfmt? ( sfmt -- ? ) + [ ix>> ] [ n>> 4 * ] bi >= ; + +: nth-sfmt ( sfmt -- n ) + [ ix>> 4 /mod swap ] + [ state>> nth nth ] + [ [ 1 + ] change-ix drop ] tri ; inline + +M: sfmt random-32* ( sfmt -- n ) + dup refill-sfmt? [ generate ] when + nth-sfmt ;