Mercurial > emacs
annotate lisp/calc/calc-mtx.el @ 89690:dde1f8021ed6
*** empty log message ***
author | Kenichi Handa <handa@m17n.org> |
---|---|
date | Mon, 29 Dec 2003 07:53:55 +0000 |
parents | 375f2633d815 |
children | 68c22ea6027c |
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> |
49598
0d8b17d428b5
Trailing whitepace deleted.
Juanma Barranquero <lekktu@gmail.com>
parents:
49263
diff
changeset
|
6 ;; Maintainers: D. Goel <deego@gnufans.org> |
49263
f4d68f97221e
Add new maintainer (deego).
Deepak Goel <deego@gnufans.org>
parents:
42206
diff
changeset
|
7 ;; Colin Walters <walters@debian.org> |
40785 | 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 | |
41271
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
parents:
41047
diff
changeset
|
26 ;;; Commentary: |
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
parents:
41047
diff
changeset
|
27 |
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
parents:
41047
diff
changeset
|
28 ;;; Code: |
40785 | 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 | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
42 (calc-unary-op "mdet" 'calcFunc-det arg))) |
40785 | 43 |
44 (defun calc-mtrace (arg) | |
45 (interactive "P") | |
46 (calc-slow-wrapper | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
47 (calc-unary-op "mtr" 'calcFunc-tr arg))) |
40785 | 48 |
49 (defun calc-mlud (arg) | |
50 (interactive "P") | |
51 (calc-slow-wrapper | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
52 (calc-unary-op "mlud" 'calcFunc-lud arg))) |
40785 | 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) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
60 a)) |
40785 | 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))) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
67 a)) |
40785 | 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))) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
87 (cons 'vec (nreverse mat)))) |
40785 | 88 |
89 (defun math-mul-mat-vec (a b) | |
90 (cons 'vec (mapcar (function (lambda (row) | |
91 (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
|
92 (cdr a)))) |
40785 | 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))) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
99 (math-reject-arg mat 'square-matrixp))) |
40785 | 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)))) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
105 sum)) |
40785 | 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 | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
168 (math-lud-solve lud (calcFunc-idn 1 n))))))) |
40785 | 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) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
177 (math-reject-arg m 'square-matrixp)))) |
40785 | 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))) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
216 0)))))) |
40785 | 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)))) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
221 prod)) |
40785 | 222 |
42206 | 223 ;;; This returns a list (LU index d), or nil if not possible. |
40785 | 224 ;;; 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
|
225 (defvar math-lud-cache nil) |
40785 | 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))) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
236 lud)))) |
40785 | 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))) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
284 (list lu (nreverse index) d))) |
40785 | 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))) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
297 m) |
40785 | 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 | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
339 (math-reject-arg need "*Singular matrix")))) |
40785 | 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")) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
366 (math-reject-arg m 'square-matrixp))) |
40785 | 367 |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
368 ;;; calc-mtx.el ends here |