1 ;;; calc-mtx.el --- matrix functions for Calc
3 ;; Copyright (C) 1990-1993, 2001-2016 Free Software Foundation, Inc.
5 ;; Author: David Gillespie <daveg@synaptics.com>
7 ;; This file is part of GNU Emacs.
9 ;; GNU Emacs is free software: you can redistribute it and/or modify
10 ;; it under the terms of the GNU General Public License as published by
11 ;; the Free Software Foundation, either version 3 of the License, or
12 ;; (at your option) any later version.
14 ;; GNU Emacs is distributed in the hope that it will be useful,
15 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
16 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 ;; GNU General Public License for more details.
19 ;; You should have received a copy of the GNU General Public License
20 ;; along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>.
26 ;; This file is autoloaded from calc-ext.el.
31 (defun calc-mdet (arg)
34 (calc-unary-op "mdet" 'calcFunc-det arg)))
36 (defun calc-mtrace (arg)
39 (calc-unary-op "mtr" 'calcFunc-tr arg)))
41 (defun calc-mlud (arg)
44 (calc-unary-op "mlud" 'calcFunc-lud arg)))
47 ;;; Coerce row vector A to be a matrix. [V V]
48 (defun math-row-matrix (a)
49 (if (and (Math-vectorp a)
50 (not (math-matrixp a)))
54 ;;; Coerce column vector A to be a matrix. [V V]
55 (defun math-col-matrix (a)
56 (if (and (Math-vectorp a)
57 (not (math-matrixp a)))
58 (cons 'vec (mapcar (function (lambda (x) (list 'vec x))) (cdr a)))
63 ;;; Multiply matrices A and B. [V V V]
64 (defun math-mul-mats (a b)
66 (cols (length (nth 1 b)))
68 (while (setq a (cdr a))
71 (while (> (setq col (1- col)) 0)
72 (setq ap (cdr (car a))
74 accum (math-mul (car ap) (nth col (car bp))))
75 (while (setq ap (cdr ap) bp (cdr bp))
76 (setq accum (math-add accum (math-mul (car ap) (nth col (car bp))))))
77 (setq row (cons accum row)))
78 (setq mat (cons (cons 'vec row) mat)))
79 (cons 'vec (nreverse mat))))
81 (defun math-mul-mat-vec (a b)
82 (cons 'vec (mapcar (function (lambda (row)
83 (math-dot-product row b)))
88 (defun calcFunc-tr (mat) ; [Public]
89 (if (math-square-matrixp mat)
90 (math-matrix-trace-step 2 (1- (length mat)) mat (nth 1 (nth 1 mat)))
91 (math-reject-arg mat 'square-matrixp)))
93 (defun math-matrix-trace-step (n size mat sum)
95 (math-matrix-trace-step (1+ n) size mat
96 (math-add sum (nth n (nth n mat))))
100 ;;; Matrix inverse and determinant.
101 (defun math-matrix-inv-raw (m)
102 (let ((n (1- (length m))))
104 (let ((det (math-det-raw m)))
105 (and (not (math-zerop det))
112 (math-neg (nth 2 (nth 1 m))))
114 (math-neg (nth 1 (nth 2 m)))
119 (math-sub (math-mul (nth 3 (nth 3 m))
121 (math-mul (nth 3 (nth 2 m))
123 (math-sub (math-mul (nth 3 (nth 1 m))
125 (math-mul (nth 3 (nth 3 m))
127 (math-sub (math-mul (nth 3 (nth 2 m))
129 (math-mul (nth 3 (nth 1 m))
132 (math-sub (math-mul (nth 3 (nth 2 m))
134 (math-mul (nth 3 (nth 3 m))
136 (math-sub (math-mul (nth 3 (nth 3 m))
138 (math-mul (nth 3 (nth 1 m))
140 (math-sub (math-mul (nth 3 (nth 1 m))
142 (math-mul (nth 3 (nth 2 m))
145 (math-sub (math-mul (nth 2 (nth 3 m))
147 (math-mul (nth 2 (nth 2 m))
149 (math-sub (math-mul (nth 2 (nth 1 m))
151 (math-mul (nth 2 (nth 3 m))
153 (math-sub (math-mul (nth 2 (nth 2 m))
155 (math-mul (nth 2 (nth 1 m))
156 (nth 1 (nth 2 m))))))))
158 (let ((lud (math-matrix-lud m)))
160 (math-lud-solve lud (calcFunc-idn 1 n)))))))
162 (defun calcFunc-det (m)
163 (if (math-square-matrixp m)
164 (math-with-extra-prec 2 (math-det-raw m))
165 (if (and (eq (car-safe m) 'calcFunc-idn)
166 (or (math-zerop (nth 1 m))
167 (math-equal-int (nth 1 m) 1)))
169 (math-reject-arg m 'square-matrixp))))
171 ;; The variable math-det-lu is local to math-det-raw, but is
172 ;; used by math-det-step, which is called by math-det-raw.
175 (defun math-det-raw (m)
176 (let ((n (1- (length m))))
180 (math-sub (math-mul (nth 1 (nth 1 m))
182 (math-mul (nth 2 (nth 1 m))
190 (math-mul (nth 1 (nth 1 m))
191 (math-mul (nth 2 (nth 2 m))
193 (math-mul (nth 2 (nth 1 m))
194 (math-mul (nth 3 (nth 2 m))
196 (math-mul (nth 3 (nth 1 m))
197 (math-mul (nth 1 (nth 2 m))
199 (math-mul (nth 3 (nth 1 m))
200 (math-mul (nth 2 (nth 2 m))
202 (math-mul (nth 1 (nth 1 m))
203 (math-mul (nth 3 (nth 2 m))
205 (math-mul (nth 2 (nth 1 m))
206 (math-mul (nth 1 (nth 2 m))
207 (nth 3 (nth 3 m))))))
208 (t (let ((lud (math-matrix-lud m)))
210 (let ((math-det-lu (car lud)))
211 (math-det-step n (nth 2 lud)))
214 (defun math-det-step (n prod)
216 (math-det-step (1- n) (math-mul prod (nth n (nth n math-det-lu))))
219 ;;; This returns a list (LU index d), or nil if not possible.
220 ;;; Argument M must be a square matrix.
221 (defvar math-lud-cache nil)
222 (defun math-matrix-lud (m)
223 (let ((old (assoc m math-lud-cache))
224 (context (list calc-internal-prec calc-prefer-frac)))
225 (if (and old (equal (nth 1 old) context))
227 (let* ((lud (catch 'singular (math-do-matrix-lud m)))
228 (entry (cons context lud)))
231 (setq math-lud-cache (cons (cons m entry) math-lud-cache)))
235 (defun math-lud-pivot-check (a)
236 "Determine a useful value for checking the size of potential pivots
237 in LUD decomposition."
238 (cond ((eq (car-safe a) 'mod)
239 (if (and (math-integerp (nth 1 a))
240 (math-integerp (nth 2 a))
241 (eq (math-gcd (nth 1 a) (nth 2 a)) 1))
245 (math-abs-approx a))))
248 ;;; Numerical Recipes section 2.3; implicit pivoting omitted.
249 (defun math-do-matrix-lud (m)
250 (let* ((lu (math-copy-matrix m))
252 i (j 1) k imax sum big
259 (math-working "LUD step" (format "%d/%d" j i))
260 (setq sum (nth j (nth i lu))
263 (setq sum (math-sub sum (math-mul (nth k (nth i lu))
266 (setcar (nthcdr j (nth i lu)) sum)
269 (math-working "LUD step" (format "%d/%d" j i))
270 (setq sum (nth j (nth i lu))
273 (setq sum (math-sub sum (math-mul (nth k (nth i lu))
276 (setcar (nthcdr j (nth i lu)) sum)
277 (let ((dum (math-lud-pivot-check sum)))
278 (if (Math-lessp big dum)
283 (setq lu (math-swap-rows lu j imax)
285 (setq index (cons imax index))
286 (let ((pivot (nth j (nth j lu))))
287 (if (math-zerop pivot)
288 (throw 'singular nil)
290 (while (<= (setq i (1+ i)) n)
291 (setcar (nthcdr j (nth i lu))
292 (math-div (nth j (nth i lu)) pivot)))))
294 (list lu (nreverse index) d)))
296 (defun math-swap-rows (m r1 r2)
298 (let* ((r1prev (nthcdr (1- r1) m))
300 (r2prev (nthcdr (1- r2) m))
305 (setcdr row2 (cdr row1))
306 (setcdr row1 r2next)))
310 (defun math-lud-solve (lud b &optional need)
312 (let* ((x (math-copy-matrix b))
314 (m (1- (length (nth 1 x))))
319 (math-working "LUD solver step" col)
326 sum (nth col (nth ip x)))
327 (setcar (nthcdr col (nth ip x)) (nth col (nth i x)))
333 (setq sum (math-sub sum (math-mul (nth j (nth i lu))
334 (nth col (nth j x))))
336 (setcar (nthcdr col (nth i x)) sum)
338 (while (>= (setq i (1- i)) 1)
339 (setq sum (nth col (nth i x))
341 (while (<= (setq j (1+ j)) n)
342 (setq sum (math-sub sum (math-mul (nth j (nth i lu))
343 (nth col (nth j x))))))
344 (setcar (nthcdr col (nth i x))
345 (math-div sum (nth i (nth i lu)))))
349 (math-reject-arg need "*Singular matrix"))))
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)))
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))
364 (setcar (nthcdr j (nth i lmat)) 0)
366 (setcar (nthcdr j (nth j lmat)) 1)
367 (while (<= (setq i (1+ i)) n)
368 (setcar (nthcdr j (nth i umat)) 0))
370 (while (>= (setq j (1- j)) 1)
371 (let ((pos (nth (1- j) (nth 1 lud))))
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)))
380 ;;; calc-mtx.el ends here