-! Copyright (C) 2006, 2008 Slava Pestov.
+! Copyright (C) 2006, 2010 Slava Pestov.
! See http://factorcode.org/license.txt for BSD license.
USING: kernel locals math math.vectors math.matrices
namespaces sequences fry sorting ;
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 ;
[ first-col ] keep
dup 1 + rows-from clear-col ;
+: find-row ( row# quot -- i elt )
+ [ rows-from ] dip find ; inline
+
: pivot-row ( col# row# -- n )
- rows-from swap '[ [ _ ] dip nth-row nth abs ] sort-with last ;
+ [ dupd nth-row nth zero? not ] find-row 2nip ;
: (echelon) ( col# row# -- )
over cols < over rows < and [
[ 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 ;