1 ! Copyright (C) 2006, 2008 Slava Pestov.
2 ! See http://factorcode.org/license.txt for BSD license.
3 USING: kernel locals math math.vectors math.matrices
4 namespaces sequences fry sorting ;
5 IN: math.matrices.elimination
9 : with-matrix ( matrix quot -- )
10 [ swap matrix set call matrix get ] with-scope ; inline
12 : nth-row ( row# -- seq ) matrix get nth ;
14 : change-row ( row# quot: ( seq -- seq ) -- )
15 matrix get swap change-nth ; inline
17 : exchange-rows ( row# row# -- ) matrix get exchange ;
19 : rows ( -- n ) matrix get length ;
21 : cols ( -- n ) 0 nth-row length ;
23 : skip ( i seq quot -- n )
24 over [ find-from drop ] dip length or ; inline
26 : first-col ( row# -- n )
27 #! First non-zero column
28 0 swap nth-row [ zero? not ] skip ;
30 : clear-scale ( col# pivot-row i-row -- n )
31 [ over ] dip nth dup zero? [
34 [ nth dup zero? ] dip swap [
41 : (clear-col) ( col# pivot-row i -- )
42 [ [ clear-scale ] 2keep [ n*v ] dip v+ ] change-row ;
44 : rows-from ( row# -- slice )
47 : clear-col ( col# row# rows -- )
48 [ nth-row ] dip [ [ 2dup ] dip (clear-col) ] each 2drop ;
50 : do-row ( exchange-with row# -- )
51 [ exchange-rows ] keep
53 dup 1 + rows-from clear-col ;
55 : pivot-row ( col# row# -- n )
56 rows-from swap '[ [ _ ] dip nth-row nth abs ] sort-with last ;
58 : (echelon) ( col# row# -- )
59 over cols < over rows < and [
60 2dup pivot-row [ over do-row 1 + ] when*
66 : echelon ( matrix -- matrix' )
67 [ 0 0 (echelon) ] with-matrix ;
69 : nonzero-rows ( matrix -- matrix' )
70 [ [ zero? ] all? not ] filter ;
72 : null/rank ( matrix -- null rank )
73 echelon dup length swap nonzero-rows length [ - ] keep ;
75 : leading ( seq -- n elt ) [ zero? not ] find ;
77 : reduced ( matrix' -- matrix'' )
80 dup nth-row leading drop
81 dup [ swap dup clear-col ] [ 2drop ] if
85 :: basis-vector ( row col# -- )
87 col# row' nth neg recip :> a
89 a row n*v col# matrix get set-nth ;
91 : nullspace ( matrix -- seq )
92 echelon reduced dup empty? [
93 dup first length identity-matrix [
96 dup [ basis-vector ] [ 2drop ] if
98 ] with-matrix flip nonzero-rows
101 : 1-pivots ( matrix -- matrix )
102 [ dup leading nip [ recip v*n ] when* ] map ;
104 : solution ( matrix -- matrix )
105 echelon nonzero-rows reduced 1-pivots ;
107 : inverse ( matrix -- matrix ) ! Assumes an invertible matrix
109 [ identity-matrix [ append ] 2map solution ] keep