-! 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
+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 ;
: cols ( -- n ) 0 nth-row length ;
: skip ( i seq quot -- n )
- over >r find* drop r> length or ; inline
+ 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 ;
: 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 ;
: 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 ;
: inverse ( matrix -- matrix ) ! Assumes an invertible matrix
dup length
- [ identity-matrix [ append ] 2map solution ] keep
+ [ <identity-matrix> [ append ] 2map solution ] keep
[ tail ] curry map ;