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