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