3 (defun tuple-multiple-nsubtrn (u k v
)
4 "Vector /u/ minus scalar /k/ times vector /v/"
5 (do ((u-rem u
(cdr u-rem
))
7 ((not (or u-rem v-rem
)) u
)
8 (decf (car u-rem
) (* k
(car v-rem
)))))
10 (defun row-clear-leading-with (u rows
)
11 (declare (inline tuple-multiple-nsubtrn
))
12 "Clear the first column with the leading-1 row."
13 (loop :for to-clear
:in rows
:collect
14 (cdr (tuple-multiple-nsubtrn
15 to-clear
(car to-clear
) u
))))
17 (defun row-reduction (a incheck-fn nolead-fn
)
21 :with toprow
= (car a
)
23 :unless
(zerop (caar rows
)) :do
25 (rotatef (caar rows
) (car toprow
))
26 (rotatef (cdar rows
) (cdr toprow
))
27 (funcall incheck-fn
(cdr nat-zero
)
30 :collect
(cdar rows
) :into nat-zero
31 :finally
(funcall nolead-fn nat-zero
))
32 (funcall incheck-fn nil
(car a
) (cdr a
)))))
34 (defun det (a &optional
(destroy nil
)
40 (setq curr-det
(* curr-det
(car row
)))
41 (row-clear-leading-with
42 (tuple-scalar-ndivisn row
(car row
))
47 (setq a tmp-a
) #'incheck-fn
#'nolead-fn
))
49 (nat-zero lead toclear
)
50 (when nat-zero
(setq curr-det
(- curr-det
)))
53 (fix-in-lead lead toclear
))))
56 (declare (ignore nat-zero
))
58 (det-recurse (if destroy a
(copy-tree a
)))
59 ;^ First call the recursion...^
60 (if a
0 curr-det
)));->... then return the result.
63 (a &optional
(destroy nil
)
64 &aux
(orig-mtrix (if destroy a
(copy-tree a
))))
70 (setq a tmp-a
) #'incheck-fn
#'nolead-fn
))
72 (nat-zero lead toclear
)
76 (row-clear-leading-with
77 (tuple-scalar-ndivisn lead
(car lead
))
79 (nolead-fn (nat-zero) (r-ef-recurse nat-zero
)))
80 (r-ef-recurse orig-mtrix
)
84 (a &optional
(destroy nil
)
86 (orig-mtrix (if destroy a
(copy-tree a
))))
87 "Reduced row-echelon form"
92 (setq a tmp-a
) #'incheck-fn
#'nolead-fn
))
94 (nat-zero lead toclear
)
95 (tuple-scalar-ndivisn lead
(car lead
))
98 (row-clear-leading-with
102 (row-clear-leading-with
106 ;V No leading 1 can be obtained, go to V
107 ;V next column in /also-clear/ without V
108 ;V clearing the current one. V
109 (do ((toclear also-clear
(cdr toclear
)))
111 (rplaca toclear
(cdar toclear
)))
112 (rr-ef-recurse nat-zero
)))
113 (rr-ef-recurse orig-mtrix
)
117 (a &key
(destroy nil
)
118 ((dimdom dimdom-left
) (length (car a
)))
119 &aux
(also-clear nil
)
120 ;V Collect basis for nullspace.V
121 (null-basis nil
) (nullity-sofar 0))
122 "A basis for the nullspace of a matrix."
126 (do ((toclear also-clear
(cdr toclear
))
129 (make-list nullity-sofar
131 (cons -
1 (make-list (decf dimdom-left
)
132 :initial-element
0)))))
135 (push newcol null-basis
))
136 (push (pop (car toclear
)) newcol
)))
140 (setq a tmp-a
) #'incheck-fn
#'nolead-fn
))
142 (nat-zero lead toclear
)
143 (tuple-scalar-ndivisn lead
(car lead
))
146 (row-clear-leading-with
151 (row-clear-leading-with
156 (nullspace-recurse nat-zero
)))
157 (nullspace-recurse (if destroy a
(copy-tree a
)))
158 ;V When done with rr-ef, add all V
159 ;V remaining columns to the basis.V
160 (dotimes (iters dimdom-left
) (add-to-basis))
163 (defun mtrix-square-identity (n)
164 "Create a square identity matrix with /n/ rows and columns."
167 :until
(= x n
) :collect
168 (nconc (make-list x
:initial-element
0)
171 :initial-element
0)))))
173 (defun mtrix-augment (a b
)
174 "Augment matrix /a/ with /b/"
178 (rplacd (last arow
) brow
))
180 ;V Return an "undo" procedure.V
182 (loop :for atail
:in tails
183 :do
(rplacd atail nil
)))))
185 ; todo: make row swaps more efficient,
186 ; perhaps a full-fledged row-reduction function is needed
188 (a &optional
(row-count (length a
)))
189 "Inverse of a given matrix."
190 (setq a
(copy-tree a
))
191 (mapc #'nconc a
(mtrix-square-identity row-count
))
192 (mapcar (lambda (row) (nthcdr row-count row
))