-! Copyright (C) 2006, 2008 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 ;
+USING: kernel locals math math.vectors math.matrices
+namespaces sequences fry sorting ;
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: ( 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 ;
[ [ 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 -- )
[ 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 )
[ rows-from ] dip find ; inline
: (echelon) ( col# row# -- )
over cols < over rows < and [
- 2dup pivot-row [ over do-row 1+ ] when*
- [ 1+ ] dip (echelon)
+ 2dup pivot-row [ over do-row 1 + ] when*
+ [ 1 + ] dip (echelon)
] [
2drop
] if ;
[ 0 0 (echelon) ] with-matrix ;
: nonzero-rows ( matrix -- matrix' )
- [ [ zero? ] all? not ] filter ;
+ [ [ zero? ] all? ] reject ;
: null/rank ( matrix -- null rank )
echelon dup length swap nonzero-rows length [ - ] keep ;
: reduced ( matrix' -- matrix'' )
[
- rows <reversed> [
+ rows iota <reversed> [
dup nth-row leading drop
- dup [ swap dup clear-col ] [ 2drop ] if
+ dup [ swap dup iota clear-col ] [ 2drop ] if
] each
] with-matrix ;
-: basis-vector ( row col# -- )
- [ clone ] dip
- [ swap nth neg recip ] 2keep
- [ 0 spin set-nth ] 2keep
- [ n*v ] dip
- 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? [