1 ;;; calc-mtx.el --- matrix functions for Calc
3 ;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc.
5 ;; Author: David Gillespie <daveg@synaptics.com>
6 ;; Maintainer: Jay Belanger <belanger@truman.edu>
8 ;; This file is part of GNU Emacs.
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.
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.
29 ;; This file is autoloaded from calc-ext.el.
34 (defun calc-mdet (arg)
37 (calc-unary-op "mdet" 'calcFunc-det arg
)))
39 (defun calc-mtrace (arg)
42 (calc-unary-op "mtr" 'calcFunc-tr arg
)))
44 (defun calc-mlud (arg)
47 (calc-unary-op "mlud" 'calcFunc-lud arg
)))
50 ;;; Coerce row vector A to be a matrix. [V V]
51 (defun math-row-matrix (a)
52 (if (and (Math-vectorp a
)
53 (not (math-matrixp a
)))
57 ;;; Coerce column vector A to be a matrix. [V V]
58 (defun math-col-matrix (a)
59 (if (and (Math-vectorp a
)
60 (not (math-matrixp a
)))
61 (cons 'vec
(mapcar (function (lambda (x) (list 'vec x
))) (cdr a
)))
66 ;;; Multiply matrices A and B. [V V V]
67 (defun math-mul-mats (a b
)
69 (cols (length (nth 1 b
)))
71 (while (setq a
(cdr a
))
74 (while (> (setq col
(1- col
)) 0)
75 (setq ap
(cdr (car a
))
77 accum
(math-mul (car ap
) (nth col
(car bp
))))
78 (while (setq ap
(cdr ap
) bp
(cdr bp
))
79 (setq accum
(math-add accum
(math-mul (car ap
) (nth col
(car bp
))))))
80 (setq row
(cons accum row
)))
81 (setq mat
(cons (cons 'vec row
) mat
)))
82 (cons 'vec
(nreverse mat
))))
84 (defun math-mul-mat-vec (a b
)
85 (cons 'vec
(mapcar (function (lambda (row)
86 (math-dot-product row b
)))
91 (defun calcFunc-tr (mat) ; [Public]
92 (if (math-square-matrixp mat
)
93 (math-matrix-trace-step 2 (1- (length mat
)) mat
(nth 1 (nth 1 mat
)))
94 (math-reject-arg mat
'square-matrixp
)))
96 (defun math-matrix-trace-step (n size mat sum
)
98 (math-matrix-trace-step (1+ n
) size mat
99 (math-add sum
(nth n
(nth n mat
))))
103 ;;; Matrix inverse and determinant.
104 (defun math-matrix-inv-raw (m)
105 (let ((n (1- (length m
))))
107 (let ((det (math-det-raw m
)))
108 (and (not (math-zerop det
))
115 (math-neg (nth 2 (nth 1 m
))))
117 (math-neg (nth 1 (nth 2 m
)))
122 (math-sub (math-mul (nth 3 (nth 3 m
))
124 (math-mul (nth 3 (nth 2 m
))
126 (math-sub (math-mul (nth 3 (nth 1 m
))
128 (math-mul (nth 3 (nth 3 m
))
130 (math-sub (math-mul (nth 3 (nth 2 m
))
132 (math-mul (nth 3 (nth 1 m
))
135 (math-sub (math-mul (nth 3 (nth 2 m
))
137 (math-mul (nth 3 (nth 3 m
))
139 (math-sub (math-mul (nth 3 (nth 3 m
))
141 (math-mul (nth 3 (nth 1 m
))
143 (math-sub (math-mul (nth 3 (nth 1 m
))
145 (math-mul (nth 3 (nth 2 m
))
148 (math-sub (math-mul (nth 2 (nth 3 m
))
150 (math-mul (nth 2 (nth 2 m
))
152 (math-sub (math-mul (nth 2 (nth 1 m
))
154 (math-mul (nth 2 (nth 3 m
))
156 (math-sub (math-mul (nth 2 (nth 2 m
))
158 (math-mul (nth 2 (nth 1 m
))
159 (nth 1 (nth 2 m
))))))))
161 (let ((lud (math-matrix-lud m
)))
163 (math-lud-solve lud
(calcFunc-idn 1 n
)))))))
165 (defun calcFunc-det (m)
166 (if (math-square-matrixp m
)
167 (math-with-extra-prec 2 (math-det-raw m
))
168 (if (and (eq (car-safe m
) 'calcFunc-idn
)
169 (or (math-zerop (nth 1 m
))
170 (math-equal-int (nth 1 m
) 1)))
172 (math-reject-arg m
'square-matrixp
))))
174 ;; The variable math-det-lu is local to math-det-raw, but is
175 ;; used by math-det-step, which is called by math-det-raw.
178 (defun math-det-raw (m)
179 (let ((n (1- (length m
))))
183 (math-sub (math-mul (nth 1 (nth 1 m
))
185 (math-mul (nth 2 (nth 1 m
))
193 (math-mul (nth 1 (nth 1 m
))
194 (math-mul (nth 2 (nth 2 m
))
196 (math-mul (nth 2 (nth 1 m
))
197 (math-mul (nth 3 (nth 2 m
))
199 (math-mul (nth 3 (nth 1 m
))
200 (math-mul (nth 1 (nth 2 m
))
202 (math-mul (nth 3 (nth 1 m
))
203 (math-mul (nth 2 (nth 2 m
))
205 (math-mul (nth 1 (nth 1 m
))
206 (math-mul (nth 3 (nth 2 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
)))
213 (let ((math-det-lu (car lud
)))
214 (math-det-step n
(nth 2 lud
)))
217 (defun math-det-step (n prod
)
219 (math-det-step (1- n
) (math-mul prod
(nth n
(nth n math-det-lu
))))
222 ;;; This returns a list (LU index d), or nil if not possible.
223 ;;; Argument M must be a square matrix.
224 (defvar math-lud-cache nil
)
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
))
230 (let* ((lud (catch 'singular
(math-do-matrix-lud m
)))
231 (entry (cons context lud
)))
234 (setq math-lud-cache
(cons (cons m entry
) math-lud-cache
)))
237 ;;; Numerical Recipes section 2.3; implicit pivoting omitted.
238 (defun math-do-matrix-lud (m)
239 (let* ((lu (math-copy-matrix m
))
241 i
(j 1) k imax sum big
248 (math-working "LUD step" (format "%d/%d" j i
))
249 (setq sum
(nth j
(nth i lu
))
252 (setq sum
(math-sub sum
(math-mul (nth k
(nth i lu
))
255 (setcar (nthcdr j
(nth i lu
)) sum
)
258 (math-working "LUD step" (format "%d/%d" j i
))
259 (setq sum
(nth j
(nth i lu
))
262 (setq sum
(math-sub sum
(math-mul (nth k
(nth i lu
))
265 (setcar (nthcdr j
(nth i lu
)) sum
)
266 (let ((dum (math-abs-approx sum
)))
267 (if (Math-lessp big dum
)
272 (setq lu
(math-swap-rows lu j imax
)
274 (setq index
(cons imax index
))
275 (let ((pivot (nth j
(nth j lu
))))
276 (if (math-zerop pivot
)
277 (throw 'singular nil
)
279 (while (<= (setq i
(1+ i
)) n
)
280 (setcar (nthcdr j
(nth i lu
))
281 (math-div (nth j
(nth i lu
)) pivot
)))))
283 (list lu
(nreverse index
) d
)))
285 (defun math-swap-rows (m r1 r2
)
287 (let* ((r1prev (nthcdr (1- r1
) m
))
289 (r2prev (nthcdr (1- r2
) m
))
294 (setcdr row2
(cdr row1
))
295 (setcdr row1 r2next
)))
299 (defun math-lud-solve (lud b
&optional need
)
301 (let* ((x (math-copy-matrix b
))
303 (m (1- (length (nth 1 x
))))
308 (math-working "LUD solver step" col
)
315 sum
(nth col
(nth ip x
)))
316 (setcar (nthcdr col
(nth ip x
)) (nth col
(nth i x
)))
322 (setq sum
(math-sub sum
(math-mul (nth j
(nth i lu
))
323 (nth col
(nth j x
))))
325 (setcar (nthcdr col
(nth i x
)) sum
)
327 (while (>= (setq i
(1- i
)) 1)
328 (setq sum
(nth col
(nth i x
))
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
)))))
338 (math-reject-arg need
"*Singular matrix"))))
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
)))
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
))
353 (setcar (nthcdr j
(nth i lmat
)) 0)
355 (setcar (nthcdr j
(nth j lmat
)) 1)
356 (while (<= (setq i
(1+ i
)) n
)
357 (setcar (nthcdr j
(nth i umat
)) 0))
359 (while (>= (setq j
(1- j
)) 1)
360 (let ((pos (nth (1- j
) (nth 1 lud
))))
362 (setq perm
(math-swap-rows perm j pos
)))))
363 (list 'vec perm lmat umat
)))))
364 (math-reject-arg m
"*Singular matrix"))
365 (math-reject-arg m
'square-matrixp
)))
369 ;;; arch-tag: fc0947b1-90e1-4a23-8950-d8ead9c3a306
370 ;;; calc-mtx.el ends here