Commit | Line | Data |
---|---|---|
3132f345 CW |
1 | ;;; calc-mtx.el --- matrix functions for Calc |
2 | ||
bf77c646 | 3 | ;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc. |
3132f345 CW |
4 | |
5 | ;; Author: David Gillespie <daveg@synaptics.com> | |
a1506d29 | 6 | ;; Maintainers: D. Goel <deego@gnufans.org> |
6e1c888a | 7 | ;; Colin Walters <walters@debian.org> |
136211a9 EZ |
8 | |
9 | ;; This file is part of GNU Emacs. | |
10 | ||
11 | ;; GNU Emacs is distributed in the hope that it will be useful, | |
12 | ;; but WITHOUT ANY WARRANTY. No author or distributor | |
13 | ;; accepts responsibility to anyone for the consequences of using it | |
14 | ;; or for whether it serves any particular purpose or works at all, | |
15 | ;; unless he says so in writing. Refer to the GNU Emacs General Public | |
16 | ;; License for full details. | |
17 | ||
18 | ;; Everyone is granted permission to copy, modify and redistribute | |
19 | ;; GNU Emacs, but only under the conditions described in the | |
20 | ;; GNU Emacs General Public License. A copy of this license is | |
21 | ;; supposed to have been given to you along with GNU Emacs so you | |
22 | ;; can know your rights and responsibilities. It should be in a | |
23 | ;; file named COPYING. Among other things, the copyright notice | |
24 | ;; and this notice must be preserved on all copies. | |
25 | ||
3132f345 CW |
26 | ;;; Commentary: |
27 | ||
28 | ;;; Code: | |
136211a9 EZ |
29 | |
30 | ||
31 | ;; This file is autoloaded from calc-ext.el. | |
32 | (require 'calc-ext) | |
33 | ||
34 | (require 'calc-macs) | |
35 | ||
36 | (defun calc-Need-calc-mat () nil) | |
37 | ||
38 | ||
39 | (defun calc-mdet (arg) | |
40 | (interactive "P") | |
41 | (calc-slow-wrapper | |
bf77c646 | 42 | (calc-unary-op "mdet" 'calcFunc-det arg))) |
136211a9 EZ |
43 | |
44 | (defun calc-mtrace (arg) | |
45 | (interactive "P") | |
46 | (calc-slow-wrapper | |
bf77c646 | 47 | (calc-unary-op "mtr" 'calcFunc-tr arg))) |
136211a9 EZ |
48 | |
49 | (defun calc-mlud (arg) | |
50 | (interactive "P") | |
51 | (calc-slow-wrapper | |
bf77c646 | 52 | (calc-unary-op "mlud" 'calcFunc-lud arg))) |
136211a9 EZ |
53 | |
54 | ||
55 | ;;; Coerce row vector A to be a matrix. [V V] | |
56 | (defun math-row-matrix (a) | |
57 | (if (and (Math-vectorp a) | |
58 | (not (math-matrixp a))) | |
59 | (list 'vec a) | |
bf77c646 | 60 | a)) |
136211a9 EZ |
61 | |
62 | ;;; Coerce column vector A to be a matrix. [V V] | |
63 | (defun math-col-matrix (a) | |
64 | (if (and (Math-vectorp a) | |
65 | (not (math-matrixp a))) | |
66 | (cons 'vec (mapcar (function (lambda (x) (list 'vec x))) (cdr a))) | |
bf77c646 | 67 | a)) |
136211a9 EZ |
68 | |
69 | ||
70 | ||
71 | ;;; Multiply matrices A and B. [V V V] | |
72 | (defun math-mul-mats (a b) | |
73 | (let ((mat nil) | |
74 | (cols (length (nth 1 b))) | |
75 | row col ap bp accum) | |
76 | (while (setq a (cdr a)) | |
77 | (setq col cols | |
78 | row nil) | |
79 | (while (> (setq col (1- col)) 0) | |
80 | (setq ap (cdr (car a)) | |
81 | bp (cdr b) | |
82 | accum (math-mul (car ap) (nth col (car bp)))) | |
83 | (while (setq ap (cdr ap) bp (cdr bp)) | |
84 | (setq accum (math-add accum (math-mul (car ap) (nth col (car bp)))))) | |
85 | (setq row (cons accum row))) | |
86 | (setq mat (cons (cons 'vec row) mat))) | |
bf77c646 | 87 | (cons 'vec (nreverse mat)))) |
136211a9 EZ |
88 | |
89 | (defun math-mul-mat-vec (a b) | |
90 | (cons 'vec (mapcar (function (lambda (row) | |
91 | (math-dot-product row b))) | |
bf77c646 | 92 | (cdr a)))) |
136211a9 EZ |
93 | |
94 | ||
95 | ||
96 | (defun calcFunc-tr (mat) ; [Public] | |
97 | (if (math-square-matrixp mat) | |
98 | (math-matrix-trace-step 2 (1- (length mat)) mat (nth 1 (nth 1 mat))) | |
bf77c646 | 99 | (math-reject-arg mat 'square-matrixp))) |
136211a9 EZ |
100 | |
101 | (defun math-matrix-trace-step (n size mat sum) | |
102 | (if (<= n size) | |
103 | (math-matrix-trace-step (1+ n) size mat | |
104 | (math-add sum (nth n (nth n mat)))) | |
bf77c646 | 105 | sum)) |
136211a9 EZ |
106 | |
107 | ||
108 | ;;; Matrix inverse and determinant. | |
109 | (defun math-matrix-inv-raw (m) | |
110 | (let ((n (1- (length m)))) | |
111 | (if (<= n 3) | |
112 | (let ((det (math-det-raw m))) | |
113 | (and (not (math-zerop det)) | |
114 | (math-div | |
115 | (cond ((= n 1) 1) | |
116 | ((= n 2) | |
117 | (list 'vec | |
118 | (list 'vec | |
119 | (nth 2 (nth 2 m)) | |
120 | (math-neg (nth 2 (nth 1 m)))) | |
121 | (list 'vec | |
122 | (math-neg (nth 1 (nth 2 m))) | |
123 | (nth 1 (nth 1 m))))) | |
124 | ((= n 3) | |
125 | (list 'vec | |
126 | (list 'vec | |
127 | (math-sub (math-mul (nth 3 (nth 3 m)) | |
128 | (nth 2 (nth 2 m))) | |
129 | (math-mul (nth 3 (nth 2 m)) | |
130 | (nth 2 (nth 3 m)))) | |
131 | (math-sub (math-mul (nth 3 (nth 1 m)) | |
132 | (nth 2 (nth 3 m))) | |
133 | (math-mul (nth 3 (nth 3 m)) | |
134 | (nth 2 (nth 1 m)))) | |
135 | (math-sub (math-mul (nth 3 (nth 2 m)) | |
136 | (nth 2 (nth 1 m))) | |
137 | (math-mul (nth 3 (nth 1 m)) | |
138 | (nth 2 (nth 2 m))))) | |
139 | (list 'vec | |
140 | (math-sub (math-mul (nth 3 (nth 2 m)) | |
141 | (nth 1 (nth 3 m))) | |
142 | (math-mul (nth 3 (nth 3 m)) | |
143 | (nth 1 (nth 2 m)))) | |
144 | (math-sub (math-mul (nth 3 (nth 3 m)) | |
145 | (nth 1 (nth 1 m))) | |
146 | (math-mul (nth 3 (nth 1 m)) | |
147 | (nth 1 (nth 3 m)))) | |
148 | (math-sub (math-mul (nth 3 (nth 1 m)) | |
149 | (nth 1 (nth 2 m))) | |
150 | (math-mul (nth 3 (nth 2 m)) | |
151 | (nth 1 (nth 1 m))))) | |
152 | (list 'vec | |
153 | (math-sub (math-mul (nth 2 (nth 3 m)) | |
154 | (nth 1 (nth 2 m))) | |
155 | (math-mul (nth 2 (nth 2 m)) | |
156 | (nth 1 (nth 3 m)))) | |
157 | (math-sub (math-mul (nth 2 (nth 1 m)) | |
158 | (nth 1 (nth 3 m))) | |
159 | (math-mul (nth 2 (nth 3 m)) | |
160 | (nth 1 (nth 1 m)))) | |
161 | (math-sub (math-mul (nth 2 (nth 2 m)) | |
162 | (nth 1 (nth 1 m))) | |
163 | (math-mul (nth 2 (nth 1 m)) | |
164 | (nth 1 (nth 2 m)))))))) | |
165 | det))) | |
166 | (let ((lud (math-matrix-lud m))) | |
167 | (and lud | |
bf77c646 | 168 | (math-lud-solve lud (calcFunc-idn 1 n))))))) |
136211a9 EZ |
169 | |
170 | (defun calcFunc-det (m) | |
171 | (if (math-square-matrixp m) | |
172 | (math-with-extra-prec 2 (math-det-raw m)) | |
173 | (if (and (eq (car-safe m) 'calcFunc-idn) | |
174 | (or (math-zerop (nth 1 m)) | |
175 | (math-equal-int (nth 1 m) 1))) | |
176 | (nth 1 m) | |
bf77c646 | 177 | (math-reject-arg m 'square-matrixp)))) |
136211a9 EZ |
178 | |
179 | (defun math-det-raw (m) | |
180 | (let ((n (1- (length m)))) | |
181 | (cond ((= n 1) | |
182 | (nth 1 (nth 1 m))) | |
183 | ((= n 2) | |
184 | (math-sub (math-mul (nth 1 (nth 1 m)) | |
185 | (nth 2 (nth 2 m))) | |
186 | (math-mul (nth 2 (nth 1 m)) | |
187 | (nth 1 (nth 2 m))))) | |
188 | ((= n 3) | |
189 | (math-sub | |
190 | (math-sub | |
191 | (math-sub | |
192 | (math-add | |
193 | (math-add | |
194 | (math-mul (nth 1 (nth 1 m)) | |
195 | (math-mul (nth 2 (nth 2 m)) | |
196 | (nth 3 (nth 3 m)))) | |
197 | (math-mul (nth 2 (nth 1 m)) | |
198 | (math-mul (nth 3 (nth 2 m)) | |
199 | (nth 1 (nth 3 m))))) | |
200 | (math-mul (nth 3 (nth 1 m)) | |
201 | (math-mul (nth 1 (nth 2 m)) | |
202 | (nth 2 (nth 3 m))))) | |
203 | (math-mul (nth 3 (nth 1 m)) | |
204 | (math-mul (nth 2 (nth 2 m)) | |
205 | (nth 1 (nth 3 m))))) | |
206 | (math-mul (nth 1 (nth 1 m)) | |
207 | (math-mul (nth 3 (nth 2 m)) | |
208 | (nth 2 (nth 3 m))))) | |
209 | (math-mul (nth 2 (nth 1 m)) | |
210 | (math-mul (nth 1 (nth 2 m)) | |
211 | (nth 3 (nth 3 m)))))) | |
212 | (t (let ((lud (math-matrix-lud m))) | |
213 | (if lud | |
214 | (let ((lu (car lud))) | |
215 | (math-det-step n (nth 2 lud))) | |
bf77c646 | 216 | 0)))))) |
136211a9 EZ |
217 | |
218 | (defun math-det-step (n prod) | |
219 | (if (> n 0) | |
220 | (math-det-step (1- n) (math-mul prod (nth n (nth n lu)))) | |
bf77c646 | 221 | prod)) |
136211a9 | 222 | |
f0529b5b | 223 | ;;; This returns a list (LU index d), or nil if not possible. |
136211a9 | 224 | ;;; Argument M must be a square matrix. |
3132f345 | 225 | (defvar math-lud-cache nil) |
136211a9 EZ |
226 | (defun math-matrix-lud (m) |
227 | (let ((old (assoc m math-lud-cache)) | |
228 | (context (list calc-internal-prec calc-prefer-frac))) | |
229 | (if (and old (equal (nth 1 old) context)) | |
230 | (cdr (cdr old)) | |
231 | (let* ((lud (catch 'singular (math-do-matrix-lud m))) | |
232 | (entry (cons context lud))) | |
233 | (if old | |
234 | (setcdr old entry) | |
235 | (setq math-lud-cache (cons (cons m entry) math-lud-cache))) | |
bf77c646 | 236 | lud)))) |
136211a9 EZ |
237 | |
238 | ;;; Numerical Recipes section 2.3; implicit pivoting omitted. | |
239 | (defun math-do-matrix-lud (m) | |
240 | (let* ((lu (math-copy-matrix m)) | |
241 | (n (1- (length lu))) | |
242 | i (j 1) k imax sum big | |
243 | (d 1) (index nil)) | |
244 | (while (<= j n) | |
245 | (setq i 1 | |
246 | big 0 | |
247 | imax j) | |
248 | (while (< i j) | |
249 | (math-working "LUD step" (format "%d/%d" j i)) | |
250 | (setq sum (nth j (nth i lu)) | |
251 | k 1) | |
252 | (while (< k i) | |
253 | (setq sum (math-sub sum (math-mul (nth k (nth i lu)) | |
254 | (nth j (nth k lu)))) | |
255 | k (1+ k))) | |
256 | (setcar (nthcdr j (nth i lu)) sum) | |
257 | (setq i (1+ i))) | |
258 | (while (<= i n) | |
259 | (math-working "LUD step" (format "%d/%d" j i)) | |
260 | (setq sum (nth j (nth i lu)) | |
261 | k 1) | |
262 | (while (< k j) | |
263 | (setq sum (math-sub sum (math-mul (nth k (nth i lu)) | |
264 | (nth j (nth k lu)))) | |
265 | k (1+ k))) | |
266 | (setcar (nthcdr j (nth i lu)) sum) | |
267 | (let ((dum (math-abs-approx sum))) | |
268 | (if (Math-lessp big dum) | |
269 | (setq big dum | |
270 | imax i))) | |
271 | (setq i (1+ i))) | |
272 | (if (> imax j) | |
273 | (setq lu (math-swap-rows lu j imax) | |
274 | d (- d))) | |
275 | (setq index (cons imax index)) | |
276 | (let ((pivot (nth j (nth j lu)))) | |
277 | (if (math-zerop pivot) | |
278 | (throw 'singular nil) | |
279 | (setq i j) | |
280 | (while (<= (setq i (1+ i)) n) | |
281 | (setcar (nthcdr j (nth i lu)) | |
282 | (math-div (nth j (nth i lu)) pivot))))) | |
283 | (setq j (1+ j))) | |
bf77c646 | 284 | (list lu (nreverse index) d))) |
136211a9 EZ |
285 | |
286 | (defun math-swap-rows (m r1 r2) | |
287 | (or (= r1 r2) | |
288 | (let* ((r1prev (nthcdr (1- r1) m)) | |
289 | (row1 (cdr r1prev)) | |
290 | (r2prev (nthcdr (1- r2) m)) | |
291 | (row2 (cdr r2prev)) | |
292 | (r2next (cdr row2))) | |
293 | (setcdr r2prev row1) | |
294 | (setcdr r1prev row2) | |
295 | (setcdr row2 (cdr row1)) | |
296 | (setcdr row1 r2next))) | |
bf77c646 | 297 | m) |
136211a9 EZ |
298 | |
299 | ||
300 | (defun math-lud-solve (lud b &optional need) | |
301 | (if lud | |
302 | (let* ((x (math-copy-matrix b)) | |
303 | (n (1- (length x))) | |
304 | (m (1- (length (nth 1 x)))) | |
305 | (lu (car lud)) | |
306 | (col 1) | |
307 | i j ip ii index sum) | |
308 | (while (<= col m) | |
309 | (math-working "LUD solver step" col) | |
310 | (setq i 1 | |
311 | ii nil | |
312 | index (nth 1 lud)) | |
313 | (while (<= i n) | |
314 | (setq ip (car index) | |
315 | index (cdr index) | |
316 | sum (nth col (nth ip x))) | |
317 | (setcar (nthcdr col (nth ip x)) (nth col (nth i x))) | |
318 | (if (null ii) | |
319 | (or (math-zerop sum) | |
320 | (setq ii i)) | |
321 | (setq j ii) | |
322 | (while (< j i) | |
323 | (setq sum (math-sub sum (math-mul (nth j (nth i lu)) | |
324 | (nth col (nth j x)))) | |
325 | j (1+ j)))) | |
326 | (setcar (nthcdr col (nth i x)) sum) | |
327 | (setq i (1+ i))) | |
328 | (while (>= (setq i (1- i)) 1) | |
329 | (setq sum (nth col (nth i x)) | |
330 | j i) | |
331 | (while (<= (setq j (1+ j)) n) | |
332 | (setq sum (math-sub sum (math-mul (nth j (nth i lu)) | |
333 | (nth col (nth j x)))))) | |
334 | (setcar (nthcdr col (nth i x)) | |
335 | (math-div sum (nth i (nth i lu))))) | |
336 | (setq col (1+ col))) | |
337 | x) | |
338 | (and need | |
bf77c646 | 339 | (math-reject-arg need "*Singular matrix")))) |
136211a9 EZ |
340 | |
341 | (defun calcFunc-lud (m) | |
342 | (if (math-square-matrixp m) | |
343 | (or (math-with-extra-prec 2 | |
344 | (let ((lud (math-matrix-lud m))) | |
345 | (and lud | |
346 | (let* ((lmat (math-copy-matrix (car lud))) | |
347 | (umat (math-copy-matrix (car lud))) | |
348 | (n (1- (length (car lud)))) | |
349 | (perm (calcFunc-idn 1 n)) | |
350 | i (j 1)) | |
351 | (while (<= j n) | |
352 | (setq i 1) | |
353 | (while (< i j) | |
354 | (setcar (nthcdr j (nth i lmat)) 0) | |
355 | (setq i (1+ i))) | |
356 | (setcar (nthcdr j (nth j lmat)) 1) | |
357 | (while (<= (setq i (1+ i)) n) | |
358 | (setcar (nthcdr j (nth i umat)) 0)) | |
359 | (setq j (1+ j))) | |
360 | (while (>= (setq j (1- j)) 1) | |
361 | (let ((pos (nth (1- j) (nth 1 lud)))) | |
362 | (or (= pos j) | |
363 | (setq perm (math-swap-rows perm j pos))))) | |
364 | (list 'vec perm lmat umat))))) | |
365 | (math-reject-arg m "*Singular matrix")) | |
bf77c646 | 366 | (math-reject-arg m 'square-matrixp))) |
136211a9 | 367 | |
bf77c646 | 368 | ;;; calc-mtx.el ends here |