first draft of threads documentation
[emacs.git] / lisp / calc / calc-mtx.el
blob6fc2d9463d4c87d4b6d21feb9b5051272dd5e24c
1 ;;; calc-mtx.el --- matrix functions for Calc
3 ;; Copyright (C) 1990-1993, 2001-2012 Free Software Foundation, Inc.
5 ;; Author: David Gillespie <daveg@synaptics.com>
6 ;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>
8 ;; This file is part of GNU Emacs.
10 ;; GNU Emacs is free software: you can redistribute it and/or modify
11 ;; it under the terms of the GNU General Public License as published by
12 ;; the Free Software Foundation, either version 3 of the License, or
13 ;; (at your option) any later version.
15 ;; GNU Emacs is distributed in the hope that it will be useful,
16 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
17 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 ;; GNU General Public License for more details.
20 ;; You should have received a copy of the GNU General Public License
21 ;; along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>.
23 ;;; Commentary:
25 ;;; Code:
27 ;; This file is autoloaded from calc-ext.el.
29 (require 'calc-ext)
30 (require 'calc-macs)
32 (defun calc-mdet (arg)
33 (interactive "P")
34 (calc-slow-wrapper
35 (calc-unary-op "mdet" 'calcFunc-det arg)))
37 (defun calc-mtrace (arg)
38 (interactive "P")
39 (calc-slow-wrapper
40 (calc-unary-op "mtr" 'calcFunc-tr arg)))
42 (defun calc-mlud (arg)
43 (interactive "P")
44 (calc-slow-wrapper
45 (calc-unary-op "mlud" 'calcFunc-lud arg)))
48 ;;; Coerce row vector A to be a matrix. [V V]
49 (defun math-row-matrix (a)
50 (if (and (Math-vectorp a)
51 (not (math-matrixp a)))
52 (list 'vec a)
53 a))
55 ;;; Coerce column vector A to be a matrix. [V V]
56 (defun math-col-matrix (a)
57 (if (and (Math-vectorp a)
58 (not (math-matrixp a)))
59 (cons 'vec (mapcar (function (lambda (x) (list 'vec x))) (cdr a)))
60 a))
64 ;;; Multiply matrices A and B. [V V V]
65 (defun math-mul-mats (a b)
66 (let ((mat nil)
67 (cols (length (nth 1 b)))
68 row col ap bp accum)
69 (while (setq a (cdr a))
70 (setq col cols
71 row nil)
72 (while (> (setq col (1- col)) 0)
73 (setq ap (cdr (car a))
74 bp (cdr b)
75 accum (math-mul (car ap) (nth col (car bp))))
76 (while (setq ap (cdr ap) bp (cdr bp))
77 (setq accum (math-add accum (math-mul (car ap) (nth col (car bp))))))
78 (setq row (cons accum row)))
79 (setq mat (cons (cons 'vec row) mat)))
80 (cons 'vec (nreverse mat))))
82 (defun math-mul-mat-vec (a b)
83 (cons 'vec (mapcar (function (lambda (row)
84 (math-dot-product row b)))
85 (cdr a))))
89 (defun calcFunc-tr (mat) ; [Public]
90 (if (math-square-matrixp mat)
91 (math-matrix-trace-step 2 (1- (length mat)) mat (nth 1 (nth 1 mat)))
92 (math-reject-arg mat 'square-matrixp)))
94 (defun math-matrix-trace-step (n size mat sum)
95 (if (<= n size)
96 (math-matrix-trace-step (1+ n) size mat
97 (math-add sum (nth n (nth n mat))))
98 sum))
101 ;;; Matrix inverse and determinant.
102 (defun math-matrix-inv-raw (m)
103 (let ((n (1- (length m))))
104 (if (<= n 3)
105 (let ((det (math-det-raw m)))
106 (and (not (math-zerop det))
107 (math-div
108 (cond ((= n 1) 1)
109 ((= n 2)
110 (list 'vec
111 (list 'vec
112 (nth 2 (nth 2 m))
113 (math-neg (nth 2 (nth 1 m))))
114 (list 'vec
115 (math-neg (nth 1 (nth 2 m)))
116 (nth 1 (nth 1 m)))))
117 ((= n 3)
118 (list 'vec
119 (list 'vec
120 (math-sub (math-mul (nth 3 (nth 3 m))
121 (nth 2 (nth 2 m)))
122 (math-mul (nth 3 (nth 2 m))
123 (nth 2 (nth 3 m))))
124 (math-sub (math-mul (nth 3 (nth 1 m))
125 (nth 2 (nth 3 m)))
126 (math-mul (nth 3 (nth 3 m))
127 (nth 2 (nth 1 m))))
128 (math-sub (math-mul (nth 3 (nth 2 m))
129 (nth 2 (nth 1 m)))
130 (math-mul (nth 3 (nth 1 m))
131 (nth 2 (nth 2 m)))))
132 (list 'vec
133 (math-sub (math-mul (nth 3 (nth 2 m))
134 (nth 1 (nth 3 m)))
135 (math-mul (nth 3 (nth 3 m))
136 (nth 1 (nth 2 m))))
137 (math-sub (math-mul (nth 3 (nth 3 m))
138 (nth 1 (nth 1 m)))
139 (math-mul (nth 3 (nth 1 m))
140 (nth 1 (nth 3 m))))
141 (math-sub (math-mul (nth 3 (nth 1 m))
142 (nth 1 (nth 2 m)))
143 (math-mul (nth 3 (nth 2 m))
144 (nth 1 (nth 1 m)))))
145 (list 'vec
146 (math-sub (math-mul (nth 2 (nth 3 m))
147 (nth 1 (nth 2 m)))
148 (math-mul (nth 2 (nth 2 m))
149 (nth 1 (nth 3 m))))
150 (math-sub (math-mul (nth 2 (nth 1 m))
151 (nth 1 (nth 3 m)))
152 (math-mul (nth 2 (nth 3 m))
153 (nth 1 (nth 1 m))))
154 (math-sub (math-mul (nth 2 (nth 2 m))
155 (nth 1 (nth 1 m)))
156 (math-mul (nth 2 (nth 1 m))
157 (nth 1 (nth 2 m))))))))
158 det)))
159 (let ((lud (math-matrix-lud m)))
160 (and lud
161 (math-lud-solve lud (calcFunc-idn 1 n)))))))
163 (defun calcFunc-det (m)
164 (if (math-square-matrixp m)
165 (math-with-extra-prec 2 (math-det-raw m))
166 (if (and (eq (car-safe m) 'calcFunc-idn)
167 (or (math-zerop (nth 1 m))
168 (math-equal-int (nth 1 m) 1)))
169 (nth 1 m)
170 (math-reject-arg m 'square-matrixp))))
172 ;; The variable math-det-lu is local to math-det-raw, but is
173 ;; used by math-det-step, which is called by math-det-raw.
174 (defvar math-det-lu)
176 (defun math-det-raw (m)
177 (let ((n (1- (length m))))
178 (cond ((= n 1)
179 (nth 1 (nth 1 m)))
180 ((= n 2)
181 (math-sub (math-mul (nth 1 (nth 1 m))
182 (nth 2 (nth 2 m)))
183 (math-mul (nth 2 (nth 1 m))
184 (nth 1 (nth 2 m)))))
185 ((= n 3)
186 (math-sub
187 (math-sub
188 (math-sub
189 (math-add
190 (math-add
191 (math-mul (nth 1 (nth 1 m))
192 (math-mul (nth 2 (nth 2 m))
193 (nth 3 (nth 3 m))))
194 (math-mul (nth 2 (nth 1 m))
195 (math-mul (nth 3 (nth 2 m))
196 (nth 1 (nth 3 m)))))
197 (math-mul (nth 3 (nth 1 m))
198 (math-mul (nth 1 (nth 2 m))
199 (nth 2 (nth 3 m)))))
200 (math-mul (nth 3 (nth 1 m))
201 (math-mul (nth 2 (nth 2 m))
202 (nth 1 (nth 3 m)))))
203 (math-mul (nth 1 (nth 1 m))
204 (math-mul (nth 3 (nth 2 m))
205 (nth 2 (nth 3 m)))))
206 (math-mul (nth 2 (nth 1 m))
207 (math-mul (nth 1 (nth 2 m))
208 (nth 3 (nth 3 m))))))
209 (t (let ((lud (math-matrix-lud m)))
210 (if lud
211 (let ((math-det-lu (car lud)))
212 (math-det-step n (nth 2 lud)))
213 0))))))
215 (defun math-det-step (n prod)
216 (if (> n 0)
217 (math-det-step (1- n) (math-mul prod (nth n (nth n math-det-lu))))
218 prod))
220 ;;; This returns a list (LU index d), or nil if not possible.
221 ;;; Argument M must be a square matrix.
222 (defvar math-lud-cache nil)
223 (defun math-matrix-lud (m)
224 (let ((old (assoc m math-lud-cache))
225 (context (list calc-internal-prec calc-prefer-frac)))
226 (if (and old (equal (nth 1 old) context))
227 (cdr (cdr old))
228 (let* ((lud (catch 'singular (math-do-matrix-lud m)))
229 (entry (cons context lud)))
230 (if old
231 (setcdr old entry)
232 (setq math-lud-cache (cons (cons m entry) math-lud-cache)))
233 lud))))
236 (defun math-lud-pivot-check (a)
237 "Determine a useful value for checking the size of potential pivots
238 in LUD decomposition."
239 (cond ((eq (car-safe a) 'mod)
240 (if (and (math-integerp (nth 1 a))
241 (math-integerp (nth 2 a))
242 (eq (math-gcd (nth 1 a) (nth 2 a)) 1))
246 (math-abs-approx a))))
249 ;;; Numerical Recipes section 2.3; implicit pivoting omitted.
250 (defun math-do-matrix-lud (m)
251 (let* ((lu (math-copy-matrix m))
252 (n (1- (length lu)))
253 i (j 1) k imax sum big
254 (d 1) (index nil))
255 (while (<= j n)
256 (setq i 1
257 big 0
258 imax j)
259 (while (< i j)
260 (math-working "LUD step" (format "%d/%d" j i))
261 (setq sum (nth j (nth i lu))
262 k 1)
263 (while (< k i)
264 (setq sum (math-sub sum (math-mul (nth k (nth i lu))
265 (nth j (nth k lu))))
266 k (1+ k)))
267 (setcar (nthcdr j (nth i lu)) sum)
268 (setq i (1+ i)))
269 (while (<= i n)
270 (math-working "LUD step" (format "%d/%d" j i))
271 (setq sum (nth j (nth i lu))
272 k 1)
273 (while (< k j)
274 (setq sum (math-sub sum (math-mul (nth k (nth i lu))
275 (nth j (nth k lu))))
276 k (1+ k)))
277 (setcar (nthcdr j (nth i lu)) sum)
278 (let ((dum (math-lud-pivot-check sum)))
279 (if (Math-lessp big dum)
280 (setq big dum
281 imax i)))
282 (setq i (1+ i)))
283 (if (> imax j)
284 (setq lu (math-swap-rows lu j imax)
285 d (- d)))
286 (setq index (cons imax index))
287 (let ((pivot (nth j (nth j lu))))
288 (if (math-zerop pivot)
289 (throw 'singular nil)
290 (setq i j)
291 (while (<= (setq i (1+ i)) n)
292 (setcar (nthcdr j (nth i lu))
293 (math-div (nth j (nth i lu)) pivot)))))
294 (setq j (1+ j)))
295 (list lu (nreverse index) d)))
297 (defun math-swap-rows (m r1 r2)
298 (or (= r1 r2)
299 (let* ((r1prev (nthcdr (1- r1) m))
300 (row1 (cdr r1prev))
301 (r2prev (nthcdr (1- r2) m))
302 (row2 (cdr r2prev))
303 (r2next (cdr row2)))
304 (setcdr r2prev row1)
305 (setcdr r1prev row2)
306 (setcdr row2 (cdr row1))
307 (setcdr row1 r2next)))
311 (defun math-lud-solve (lud b &optional need)
312 (if lud
313 (let* ((x (math-copy-matrix b))
314 (n (1- (length x)))
315 (m (1- (length (nth 1 x))))
316 (lu (car lud))
317 (col 1)
318 i j ip ii index sum)
319 (while (<= col m)
320 (math-working "LUD solver step" col)
321 (setq i 1
322 ii nil
323 index (nth 1 lud))
324 (while (<= i n)
325 (setq ip (car index)
326 index (cdr index)
327 sum (nth col (nth ip x)))
328 (setcar (nthcdr col (nth ip x)) (nth col (nth i x)))
329 (if (null ii)
330 (or (math-zerop sum)
331 (setq ii i))
332 (setq j ii)
333 (while (< j i)
334 (setq sum (math-sub sum (math-mul (nth j (nth i lu))
335 (nth col (nth j x))))
336 j (1+ j))))
337 (setcar (nthcdr col (nth i x)) sum)
338 (setq i (1+ i)))
339 (while (>= (setq i (1- i)) 1)
340 (setq sum (nth col (nth i x))
341 j i)
342 (while (<= (setq j (1+ j)) n)
343 (setq sum (math-sub sum (math-mul (nth j (nth i lu))
344 (nth col (nth j x))))))
345 (setcar (nthcdr col (nth i x))
346 (math-div sum (nth i (nth i lu)))))
347 (setq col (1+ col)))
349 (and need
350 (math-reject-arg need "*Singular matrix"))))
352 (defun calcFunc-lud (m)
353 (if (math-square-matrixp m)
354 (or (math-with-extra-prec 2
355 (let ((lud (math-matrix-lud m)))
356 (and lud
357 (let* ((lmat (math-copy-matrix (car lud)))
358 (umat (math-copy-matrix (car lud)))
359 (n (1- (length (car lud))))
360 (perm (calcFunc-idn 1 n))
361 i (j 1))
362 (while (<= j n)
363 (setq i 1)
364 (while (< i j)
365 (setcar (nthcdr j (nth i lmat)) 0)
366 (setq i (1+ i)))
367 (setcar (nthcdr j (nth j lmat)) 1)
368 (while (<= (setq i (1+ i)) n)
369 (setcar (nthcdr j (nth i umat)) 0))
370 (setq j (1+ j)))
371 (while (>= (setq j (1- j)) 1)
372 (let ((pos (nth (1- j) (nth 1 lud))))
373 (or (= pos j)
374 (setq perm (math-swap-rows perm j pos)))))
375 (list 'vec perm lmat umat)))))
376 (math-reject-arg m "*Singular matrix"))
377 (math-reject-arg m 'square-matrixp)))
379 (provide 'calc-mtx)
381 ;;; calc-mtx.el ends here