factor: trim using lists

[factor.git] / extra / math / matrices / elimination / elimination.factor

diff --git a/extra/math/matrices/elimination/elimination.factor b/extra/math/matrices/elimination/elimination.factor
old mode 100755 (executable)
new mode 100644 (file)
index 73f6dd7..970e316
--- a/extra/math/matrices/elimination/elimination.factor
+++ b/extra/math/matrices/elimination/elimination.factor
@@ -1,17 +1,17 @@
-! Copyright (C) 2006, 2007 Slava Pestov.
+! Copyright (C) 2006, 2010 Slava Pestov.
 ! See http://factorcode.org/license.txt for BSD license.
-USING: kernel math math.vectors math.matrices namespaces
-sequences parser ;
+USING: kernel math math.matrices math.vectors namespaces
+sequences ;
 IN: math.matrices.elimination
 
 SYMBOL: matrix
 
 : with-matrix ( matrix quot -- )
-    [ swap matrix set call matrix get ] with-scope ; inline
+    matrix swap [ matrix get ] compose with-variable ; inline
 
 : nth-row ( row# -- seq ) matrix get nth ;
 
-: change-row ( row# quot -- | quot: seq -- seq )
+: change-row ( ..a row# quot: ( ..a seq -- ..b seq ) -- ..b )
     matrix get swap change-nth ; inline
 
 : exchange-rows ( row# row# -- ) matrix get exchange ;
@@ -20,45 +20,48 @@ SYMBOL: matrix
 
 : cols ( -- n ) 0 nth-row length ;
 
+: skip ( i seq quot -- n )
+    over [ find-from drop ] dip swap [ ] [ length ] ?if ; inline
+
 : first-col ( row# -- n )
-    #! First non-zero column
+    ! First non-zero column
     0 swap nth-row [ zero? not ] skip ;
 
 : clear-scale ( col# pivot-row i-row -- n )
-    >r over r> nth dup zero? [
+    overd nth dup zero? [
         3drop 0
     ] [
-        >r nth dup zero? [
-            r> 2drop 0
+        [ nth dup zero? ] dip swap [
+            2drop 0
         ] [
-            r> swap / neg
+            swap / neg
         ] if
     ] if ;
 
 : (clear-col) ( col# pivot-row i -- )
-    [ [ clear-scale ] 2keep >r n*v r> v+ ] change-row ;
+    [ [ clear-scale ] 2keep [ n*v ] dip v+ ] change-row ;
 
 : rows-from ( row# -- slice )
-    rows dup <slice> ;
+    rows dup <iota> <slice> ;
 
 : clear-col ( col# row# rows -- )
-    >r nth-row r> [ >r 2dup r> (clear-col) ] each 2drop ;
+    [ nth-row ] dip [ [ 2dup ] dip (clear-col) ] each 2drop ;
 
 : do-row ( exchange-with row# -- )
     [ exchange-rows ] keep
     [ first-col ] keep
-    dup 1+ rows-from clear-col ;
+    dup 1 + rows-from clear-col ;
 
 : find-row ( row# quot -- i elt )
-    >r rows-from r> find ; inline
+    [ rows-from ] dip find ; inline
 
 : pivot-row ( col# row# -- n )
     [ dupd nth-row nth zero? not ] find-row 2nip ;
 
 : (echelon) ( col# row# -- )
     over cols < over rows < and [
-        2dup pivot-row [ over do-row 1+ ] when*
-        >r 1+ r> (echelon)
+        2dup pivot-row [ over do-row 1 + ] when*
+        [ 1 + ] dip (echelon)
     ] [
         2drop
     ] if ;
@@ -66,7 +69,8 @@ SYMBOL: matrix
 : echelon ( matrix -- matrix' )
     [ 0 0 (echelon) ] with-matrix ;
 
-: nonzero-rows [ [ zero? ] all? not ] subset ;
+: nonzero-rows ( matrix -- matrix' )
+    [ [ zero? ] all? ] reject ;
 
 : null/rank ( matrix -- null rank )
     echelon dup length swap nonzero-rows length [ - ] keep ;
@@ -75,25 +79,24 @@ SYMBOL: matrix
 
 : reduced ( matrix' -- matrix'' )
     [
-        rows <reversed> [
+        rows <iota> <reversed> [
             dup nth-row leading drop
-            dup [ swap dup clear-col ] [ 2drop ] if
+            [ swap dup <iota> clear-col ] [ drop ] if*
         ] each
     ] with-matrix ;
 
-: basis-vector ( row col# -- )
-    >r clone r>
-    [ swap nth neg recip ] 2keep
-    [ 0 spin set-nth ] 2keep
-    >r n*v r>
-    matrix get set-nth ;
+:: basis-vector ( row col# -- )
+    row clone :> row'
+    col# row' nth neg recip :> a
+    0 col# row' set-nth
+    a row n*v col# matrix get set-nth ;
 
 : nullspace ( matrix -- seq )
     echelon reduced dup empty? [
-        dup first length identity-matrix [
+        dup first length <identity-matrix> [
             [
                 dup leading drop
-                dup [ basis-vector ] [ 2drop ] if
+                [ basis-vector ] [ drop ] if*
             ] each
         ] with-matrix flip nonzero-rows
     ] unless ;
@@ -106,5 +109,5 @@ SYMBOL: matrix
 
 : inverse ( matrix -- matrix ) ! Assumes an invertible matrix
     dup length
-    [ identity-matrix [ append ] 2map solution ] keep
+    [ <identity-matrix> [ append ] 2map solution ] keep
     [ tail ] curry map ;