Mercurial > emacs
comparison lisp/calc/calc-mtx.el @ 40785:2fb9d407ae73
Initial import of Calc 2.02f.
author | Eli Zaretskii <eliz@gnu.org> |
---|---|
date | Tue, 06 Nov 2001 18:59:06 +0000 |
parents | |
children | 73f364fd8aaa |
comparison
equal
deleted
inserted
replaced
40784:d57f74c55909 | 40785:2fb9d407ae73 |
---|---|
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 |