40785
|
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
|