From 1239984625846787e6b57ab4a3e42f44b8d651bf Mon Sep 17 00:00:00 2001
From: Doug Coleman <doug.coleman@gmail.com>
Date: Fri, 2 Oct 2009 01:18:18 -0500
Subject: [PATCH] optimizing random.sfmt

---
 basis/random/sfmt/sfmt-tests.factor |   9 +-
 basis/random/sfmt/sfmt.factor       | 145 +++++++++++++++-------------
 2 files changed, 85 insertions(+), 69 deletions(-)

diff --git a/basis/random/sfmt/sfmt-tests.factor b/basis/random/sfmt/sfmt-tests.factor
index 674693805d..7e6a996278 100644
--- a/basis/random/sfmt/sfmt-tests.factor
+++ b/basis/random/sfmt/sfmt-tests.factor
@@ -1,8 +1,13 @@
 ! Copyright (C) 2009 Doug Coleman.
 ! See http://factorcode.org/license.txt for BSD license.
-USING: accessors kernel random.sfmt sequences tools.test ;
+USING: accessors kernel random random.sfmt random.sfmt.private
+sequences tools.test ;
 IN: random.sfmt.tests
 
 [ ] [ 100 <sfmt-19937> drop ] unit-test
+
 [ 1096298955 ]
-[ 100 <sfmt-19937> generate generate state>> first first ] unit-test
+[ 100 <sfmt-19937> dup generate dup generate uint-array>> first ] unit-test
+
+[ 2556114782 ]
+[ 100 <sfmt-19937> random-32* ] unit-test
diff --git a/basis/random/sfmt/sfmt.factor b/basis/random/sfmt/sfmt.factor
index 15dc23df5c..23f4b8ad33 100644
--- a/basis/random/sfmt/sfmt.factor
+++ b/basis/random/sfmt/sfmt.factor
@@ -2,35 +2,30 @@
 ! 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 sequences.private ;
-IN: random.sfmt
-
+sequences specialized-arrays sequences.private classes.struct ;
 SIMD: uint
 SPECIALIZED-ARRAY: uint
 SPECIALIZED-ARRAY: uint-4
+IN: random.sfmt
 
-CONSTANT: SFMT_N 156
-CONSTANT: SFMT_M 122
+<PRIVATE
 
 CONSTANT: state-multiplier 1812433253
 
-TUPLE: sfmt
-sl1 sl2 sr1 sr2 mask parity
-{ r1 uint-4 } { r2 uint-4 }
-{ seed fixnum } { n fixnum } { m fixnum }
-{ m-n fixnum } { ix fixnum } { state uint-4-array } ;
+STRUCT: sfmt-state
+    { seed uint }
+    { n uint }
+    { m uint }
+    { m-n int }
+    { ix uint }
+    { mask uint-4 }
+    { r1 uint-4 }
+    { r2 uint-4 } ;
 
-: 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 ;
+TUPLE: sfmt
+{ state sfmt-state }
+{ uint-array uint-array }
+{ uint-4-array uint-4-array } ;
 
 : wA ( w -- wA )
    dup 1 hlshift vbitxor ; inline
@@ -49,69 +44,85 @@ sl1 sl2 sr1 sr2 mask parity
     [ wB ] dip vbitxor
     [ wA ] dip vbitxor ; inline
 
-GENERIC: generate ( sfmt -- sfmt' )
+GENERIC: generate ( sfmt -- )
 
-M:: sfmt generate ( sfmt -- sfmt' )
+M:: sfmt generate ( sfmt -- )
     sfmt state>> :> state
-    sfmt n>> 2 - state nth-unsafe sfmt (>>r1)
-    sfmt n>> 1 - state nth-unsafe sfmt (>>r2)
-    sfmt m>> :> m 
-    sfmt n>> :> n 
-    sfmt m-n>> :> m-n
-    sfmt mask>> :> mask
+    sfmt uint-4-array>> :> array
+    state n>> 2 - array nth state (>>r1)
+    state n>> 1 - array nth state (>>r2)
+    state m>> :> m 
+    state n>> :> n 
+    state m-n>> :> m-n
+    state mask>> :> mask
 
-    n m - iota [| i |
-        i state nth-unsafe 
-        i m + state nth-unsafe
-        mask sfmt r1>> sfmt r2>> formula :> r
+    n m - >fixnum iota [| i |
+        i array nth-unsafe 
+        i m + array nth-unsafe
+        mask state r1>> state r2>> formula :> r
 
-        r i state set-nth-unsafe
-        sfmt r2>> sfmt (>>r1)
-        r sfmt (>>r2)
+        r i array set-nth-unsafe
+        state r2>> state (>>r1)
+        r state (>>r2)
     ] each
 
-    n m - 1 + n [a,b) [| i |
-        i state nth-unsafe
-        m-n i + state nth-unsafe
-        mask sfmt r1>> sfmt r2>> formula :> r
+    ! n m - 1 + n [a,b) [
+    m 1 - iota [
+        n m - 1 + + >fixnum :> i
+        i array nth-unsafe
+        m-n i + array nth-unsafe
+        mask state r1>> state r2>> formula :> r
 
-        r i state set-nth-unsafe
-        sfmt r2>> sfmt (>>r1)
-        r sfmt (>>r2)
+        r i array set-nth-unsafe
+        state r2>> state (>>r1)
+        r state (>>r2)
     ] each
     
-    sfmt 0 >>ix ;
+    0 state (>>ix) ;
 
-: <sfmt> ( seed n m sl1 sl2 sr1 sr2 mask parity -- sfmt )
-    sfmt new
-        swap >>parity
+: <sfmt-array> ( sfmt -- uint-array uint-4-array )
+    state>> 
+    [ 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
+    dup underlying>> byte-array>uint-4-array ;
+
+: <sfmt-state> ( seed n m mask -- sfmt )
+    sfmt-state <struct>
         swap >>mask
-        swap >>sr2
-        swap >>sr1
-        swap >>sl2
-        swap >>sl1
         swap >>m
         swap >>n
-        swap 32 bits >>seed
+        swap >>seed
         dup [ m>> ] [ n>> ] bi - >>m-n
-        0 >>ix
-        init-state
-        generate ;
+        0 >>ix ;
 
-: <sfmt-19937> ( 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 }
-    <sfmt> ;
+: <sfmt> ( seed n m mask -- sfmt )
+    <sfmt-state>
+    sfmt new
+        swap >>state
+        dup <sfmt-array> [ >>uint-array ] [ >>uint-4-array ] bi*
+        [ generate ] keep ;
 
 : refill-sfmt? ( sfmt -- ? )
-    [ ix>> ] [ n>> 4 * ] bi >= ;
+    state>> [ ix>> ] [ n>> 4 * ] bi >= ;
 
-: nth-sfmt ( sfmt -- n )
-    [ ix>> 4 /mod swap ]
-    [ state>> nth nth ]
-    [ [ 1 + ] change-ix drop ] tri ; inline
+: next-ix ( sfmt -- ix )
+    state>> [ dup 1 + ] change-ix drop ; inline
+
+: next ( sfmt -- n )
+    [ next-ix ] [ uint-array>> ] bi nth-unsafe ; inline
+
+PRIVATE>
 
 M: sfmt random-32* ( sfmt -- n )
-    dup refill-sfmt? [ generate ] when
-    nth-sfmt ;
+    dup refill-sfmt? [ dup generate ] when next ;
+
+: <sfmt-19937> ( seed -- sfmt )
+    348 330 uint-4{ HEX: BFFFFFF6 HEX: BFFAFFFF HEX: DDFECB7F HEX: DFFFFFEF }
+    <sfmt> ;