15 : MIN-VAL-adsoda ( -- x ) 0.00000001
23 ! [ number>string string>number ] map
25 : with-matrix ( matrix quot -- )
26 [ swap matrix set call matrix get ] with-scope ; inline
28 : nth-row ( row# -- seq ) matrix get nth ;
30 : change-row ( row# quot -- seq ) ! row# quot -- | quot: seq -- seq )
31 matrix get swap change-nth ; inline
33 : exchange-rows ( row# row# -- ) matrix get exchange ;
35 : rows ( -- n ) matrix get length ;
37 : cols ( -- n ) 0 nth-row length ;
39 : skip ( i seq quot -- n )
40 over [ find-from drop ] dip length or ; inline
42 : first-col ( row# -- n )
43 #! First non-zero column
44 0 swap nth-row [ zero? not ] skip ;
46 : clear-scale ( col# pivot-row i-row -- n )
47 [ over ] dip nth dup zero? [
50 [ nth dup zero? ] dip swap [
57 : (clear-col) ( col# pivot-row i -- )
58 [ [ clear-scale ] 2keep [ n*v ] dip v+ ] change-row ;
60 : rows-from ( row# -- slice )
63 : clear-col ( col# row# rows -- )
64 [ nth-row ] dip [ [ 2dup ] dip (clear-col) ] each 2drop ;
66 : do-row ( exchange-with row# -- )
67 [ exchange-rows ] keep
69 dup 1 + rows-from clear-col ;
71 : find-row ( row# quot -- i elt )
72 [ rows-from ] dip find ; inline
74 : pivot-row ( col# row# -- n )
75 [ dupd nth-row nth zero? not ] find-row 2nip ;
77 : (echelon) ( col# row# -- )
78 over cols < over rows < and [
79 2dup pivot-row [ over do-row 1 + ] when*
85 : echelon ( matrix -- matrix' )
86 [ 0 0 (echelon) ] with-matrix ;
88 : nonzero-rows ( matrix -- matrix' )
89 [ [ zero? ] all? ] reject ;
91 : null/rank ( matrix -- null rank )
92 echelon dup length swap nonzero-rows length [ - ] keep ;
94 : leading ( seq -- n elt ) [ zero? not ] find ;
96 : reduced ( matrix' -- matrix'' )
99 dup nth-row leading drop
100 dup [ swap dup clear-col ] [ 2drop ] if
104 : basis-vector ( row col# -- )
106 [ swap nth neg recip ] 2keep
107 [ 0 spin set-nth ] 2keep
111 : nullspace ( matrix -- seq )
112 echelon reduced dup empty? [
113 dup first length identity-matrix [
116 dup [ basis-vector ] [ 2drop ] if
118 ] with-matrix flip nonzero-rows
121 : 1-pivots ( matrix -- matrix )
122 [ dup leading nip [ recip v*n ] when* ] map ;
124 : solution ( matrix -- matrix )
125 echelon nonzero-rows reduced 1-pivots ;