Mercurial > emacs
annotate lisp/calc/calc-mtx.el @ 74425:7bbd3732cfed
(HAVE_SYNC): Define.
author | Eli Zaretskii <eliz@gnu.org> |
---|---|
date | Mon, 04 Dec 2006 20:34:53 +0000 |
parents | 6bf177f8065b |
children | 7a3f13e2dd57 c5406394f567 |
rev | line source |
---|---|
41271
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
parents:
41047
diff
changeset
|
1 ;;; calc-mtx.el --- matrix functions for Calc |
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
parents:
41047
diff
changeset
|
2 |
64325
1db49616ce05
Update copyright information.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
62442
diff
changeset
|
3 ;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004, |
68636
6bf177f8065b
Update copyright year.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
64325
diff
changeset
|
4 ;; 2005, 2006 Free Software Foundation, Inc. |
41271
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
parents:
41047
diff
changeset
|
5 |
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
parents:
41047
diff
changeset
|
6 ;; Author: David Gillespie <daveg@synaptics.com> |
58548
f2bcc073e61e
(calc-det-lu): New variable.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
52401
diff
changeset
|
7 ;; Maintainer: Jay Belanger <belanger@truman.edu> |
40785 | 8 |
9 ;; This file is part of GNU Emacs. | |
10 | |
11 ;; GNU Emacs is distributed in the hope that it will be useful, | |
12 ;; but WITHOUT ANY WARRANTY. No author or distributor | |
13 ;; accepts responsibility to anyone for the consequences of using it | |
14 ;; or for whether it serves any particular purpose or works at all, | |
15 ;; unless he says so in writing. Refer to the GNU Emacs General Public | |
16 ;; License for full details. | |
17 | |
18 ;; Everyone is granted permission to copy, modify and redistribute | |
19 ;; GNU Emacs, but only under the conditions described in the | |
20 ;; GNU Emacs General Public License. A copy of this license is | |
21 ;; supposed to have been given to you along with GNU Emacs so you | |
22 ;; can know your rights and responsibilities. It should be in a | |
23 ;; file named COPYING. Among other things, the copyright notice | |
24 ;; and this notice must be preserved on all copies. | |
25 | |
41271
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
parents:
41047
diff
changeset
|
26 ;;; Commentary: |
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
parents:
41047
diff
changeset
|
27 |
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
parents:
41047
diff
changeset
|
28 ;;; Code: |
40785 | 29 |
30 ;; This file is autoloaded from calc-ext.el. | |
31 | |
58666
b21f320635a7
Add a provide statement.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
58548
diff
changeset
|
32 (require 'calc-ext) |
40785 | 33 (require 'calc-macs) |
34 | |
35 (defun calc-mdet (arg) | |
36 (interactive "P") | |
37 (calc-slow-wrapper | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
38 (calc-unary-op "mdet" 'calcFunc-det arg))) |
40785 | 39 |
40 (defun calc-mtrace (arg) | |
41 (interactive "P") | |
42 (calc-slow-wrapper | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
43 (calc-unary-op "mtr" 'calcFunc-tr arg))) |
40785 | 44 |
45 (defun calc-mlud (arg) | |
46 (interactive "P") | |
47 (calc-slow-wrapper | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
48 (calc-unary-op "mlud" 'calcFunc-lud arg))) |
40785 | 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) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
56 a)) |
40785 | 57 |
58 ;;; Coerce column vector A to be a matrix. [V V] | |
59 (defun math-col-matrix (a) | |
60 (if (and (Math-vectorp a) | |
61 (not (math-matrixp a))) | |
62 (cons 'vec (mapcar (function (lambda (x) (list 'vec x))) (cdr a))) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
63 a)) |
40785 | 64 |
65 | |
66 | |
67 ;;; Multiply matrices A and B. [V V V] | |
68 (defun math-mul-mats (a b) | |
69 (let ((mat nil) | |
70 (cols (length (nth 1 b))) | |
71 row col ap bp accum) | |
72 (while (setq a (cdr a)) | |
73 (setq col cols | |
74 row nil) | |
75 (while (> (setq col (1- col)) 0) | |
76 (setq ap (cdr (car a)) | |
77 bp (cdr b) | |
78 accum (math-mul (car ap) (nth col (car bp)))) | |
79 (while (setq ap (cdr ap) bp (cdr bp)) | |
80 (setq accum (math-add accum (math-mul (car ap) (nth col (car bp)))))) | |
81 (setq row (cons accum row))) | |
82 (setq mat (cons (cons 'vec row) mat))) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
83 (cons 'vec (nreverse mat)))) |
40785 | 84 |
85 (defun math-mul-mat-vec (a b) | |
86 (cons 'vec (mapcar (function (lambda (row) | |
87 (math-dot-product row b))) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
88 (cdr a)))) |
40785 | 89 |
90 | |
91 | |
92 (defun calcFunc-tr (mat) ; [Public] | |
93 (if (math-square-matrixp mat) | |
94 (math-matrix-trace-step 2 (1- (length mat)) mat (nth 1 (nth 1 mat))) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
95 (math-reject-arg mat 'square-matrixp))) |
40785 | 96 |
97 (defun math-matrix-trace-step (n size mat sum) | |
98 (if (<= n size) | |
99 (math-matrix-trace-step (1+ n) size mat | |
100 (math-add sum (nth n (nth n mat)))) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
101 sum)) |
40785 | 102 |
103 | |
104 ;;; Matrix inverse and determinant. | |
105 (defun math-matrix-inv-raw (m) | |
106 (let ((n (1- (length m)))) | |
107 (if (<= n 3) | |
108 (let ((det (math-det-raw m))) | |
109 (and (not (math-zerop det)) | |
110 (math-div | |
111 (cond ((= n 1) 1) | |
112 ((= n 2) | |
113 (list 'vec | |
114 (list 'vec | |
115 (nth 2 (nth 2 m)) | |
116 (math-neg (nth 2 (nth 1 m)))) | |
117 (list 'vec | |
118 (math-neg (nth 1 (nth 2 m))) | |
119 (nth 1 (nth 1 m))))) | |
120 ((= n 3) | |
121 (list 'vec | |
122 (list 'vec | |
123 (math-sub (math-mul (nth 3 (nth 3 m)) | |
124 (nth 2 (nth 2 m))) | |
125 (math-mul (nth 3 (nth 2 m)) | |
126 (nth 2 (nth 3 m)))) | |
127 (math-sub (math-mul (nth 3 (nth 1 m)) | |
128 (nth 2 (nth 3 m))) | |
129 (math-mul (nth 3 (nth 3 m)) | |
130 (nth 2 (nth 1 m)))) | |
131 (math-sub (math-mul (nth 3 (nth 2 m)) | |
132 (nth 2 (nth 1 m))) | |
133 (math-mul (nth 3 (nth 1 m)) | |
134 (nth 2 (nth 2 m))))) | |
135 (list 'vec | |
136 (math-sub (math-mul (nth 3 (nth 2 m)) | |
137 (nth 1 (nth 3 m))) | |
138 (math-mul (nth 3 (nth 3 m)) | |
139 (nth 1 (nth 2 m)))) | |
140 (math-sub (math-mul (nth 3 (nth 3 m)) | |
141 (nth 1 (nth 1 m))) | |
142 (math-mul (nth 3 (nth 1 m)) | |
143 (nth 1 (nth 3 m)))) | |
144 (math-sub (math-mul (nth 3 (nth 1 m)) | |
145 (nth 1 (nth 2 m))) | |
146 (math-mul (nth 3 (nth 2 m)) | |
147 (nth 1 (nth 1 m))))) | |
148 (list 'vec | |
149 (math-sub (math-mul (nth 2 (nth 3 m)) | |
150 (nth 1 (nth 2 m))) | |
151 (math-mul (nth 2 (nth 2 m)) | |
152 (nth 1 (nth 3 m)))) | |
153 (math-sub (math-mul (nth 2 (nth 1 m)) | |
154 (nth 1 (nth 3 m))) | |
155 (math-mul (nth 2 (nth 3 m)) | |
156 (nth 1 (nth 1 m)))) | |
157 (math-sub (math-mul (nth 2 (nth 2 m)) | |
158 (nth 1 (nth 1 m))) | |
159 (math-mul (nth 2 (nth 1 m)) | |
160 (nth 1 (nth 2 m)))))))) | |
161 det))) | |
162 (let ((lud (math-matrix-lud m))) | |
163 (and lud | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
164 (math-lud-solve lud (calcFunc-idn 1 n))))))) |
40785 | 165 |
166 (defun calcFunc-det (m) | |
167 (if (math-square-matrixp m) | |
168 (math-with-extra-prec 2 (math-det-raw m)) | |
169 (if (and (eq (car-safe m) 'calcFunc-idn) | |
170 (or (math-zerop (nth 1 m)) | |
171 (math-equal-int (nth 1 m) 1))) | |
172 (nth 1 m) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
173 (math-reject-arg m 'square-matrixp)))) |
40785 | 174 |
58548
f2bcc073e61e
(calc-det-lu): New variable.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
52401
diff
changeset
|
175 ;; The variable math-det-lu is local to math-det-raw, but is |
f2bcc073e61e
(calc-det-lu): New variable.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
52401
diff
changeset
|
176 ;; used by math-det-step, which is called by math-det-raw. |
f2bcc073e61e
(calc-det-lu): New variable.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
52401
diff
changeset
|
177 (defvar math-det-lu) |
f2bcc073e61e
(calc-det-lu): New variable.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
52401
diff
changeset
|
178 |
40785 | 179 (defun math-det-raw (m) |
180 (let ((n (1- (length m)))) | |
181 (cond ((= n 1) | |
182 (nth 1 (nth 1 m))) | |
183 ((= n 2) | |
184 (math-sub (math-mul (nth 1 (nth 1 m)) | |
185 (nth 2 (nth 2 m))) | |
186 (math-mul (nth 2 (nth 1 m)) | |
187 (nth 1 (nth 2 m))))) | |
188 ((= n 3) | |
189 (math-sub | |
190 (math-sub | |
191 (math-sub | |
192 (math-add | |
193 (math-add | |
194 (math-mul (nth 1 (nth 1 m)) | |
195 (math-mul (nth 2 (nth 2 m)) | |
196 (nth 3 (nth 3 m)))) | |
197 (math-mul (nth 2 (nth 1 m)) | |
198 (math-mul (nth 3 (nth 2 m)) | |
199 (nth 1 (nth 3 m))))) | |
200 (math-mul (nth 3 (nth 1 m)) | |
201 (math-mul (nth 1 (nth 2 m)) | |
202 (nth 2 (nth 3 m))))) | |
203 (math-mul (nth 3 (nth 1 m)) | |
204 (math-mul (nth 2 (nth 2 m)) | |
205 (nth 1 (nth 3 m))))) | |
206 (math-mul (nth 1 (nth 1 m)) | |
207 (math-mul (nth 3 (nth 2 m)) | |
208 (nth 2 (nth 3 m))))) | |
209 (math-mul (nth 2 (nth 1 m)) | |
210 (math-mul (nth 1 (nth 2 m)) | |
211 (nth 3 (nth 3 m)))))) | |
212 (t (let ((lud (math-matrix-lud m))) | |
213 (if lud | |
58548
f2bcc073e61e
(calc-det-lu): New variable.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
52401
diff
changeset
|
214 (let ((math-det-lu (car lud))) |
40785 | 215 (math-det-step n (nth 2 lud))) |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
216 0)))))) |
40785 | 217 |
218 (defun math-det-step (n prod) | |
219 (if (> n 0) | |
58548
f2bcc073e61e
(calc-det-lu): New variable.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
52401
diff
changeset
|
220 (math-det-step (1- n) (math-mul prod (nth n (nth n math-det-lu)))) |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
221 prod)) |
40785 | 222 |
42206 | 223 ;;; This returns a list (LU index d), or nil if not possible. |
40785 | 224 ;;; Argument M must be a square matrix. |
41271
fcd507927105
Change all toplevel `setq' forms to `defvar' forms, and move them
Colin Walters <walters@gnu.org>
parents:
41047
diff
changeset
|
225 (defvar math-lud-cache nil) |
40785 | 226 (defun math-matrix-lud (m) |
227 (let ((old (assoc m math-lud-cache)) | |
228 (context (list calc-internal-prec calc-prefer-frac))) | |
229 (if (and old (equal (nth 1 old) context)) | |
230 (cdr (cdr old)) | |
231 (let* ((lud (catch 'singular (math-do-matrix-lud m))) | |
232 (entry (cons context lud))) | |
233 (if old | |
234 (setcdr old entry) | |
235 (setq math-lud-cache (cons (cons m entry) math-lud-cache))) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
236 lud)))) |
40785 | 237 |
238 ;;; Numerical Recipes section 2.3; implicit pivoting omitted. | |
239 (defun math-do-matrix-lud (m) | |
240 (let* ((lu (math-copy-matrix m)) | |
241 (n (1- (length lu))) | |
242 i (j 1) k imax sum big | |
243 (d 1) (index nil)) | |
244 (while (<= j n) | |
245 (setq i 1 | |
246 big 0 | |
247 imax j) | |
248 (while (< i j) | |
249 (math-working "LUD step" (format "%d/%d" j i)) | |
250 (setq sum (nth j (nth i lu)) | |
251 k 1) | |
252 (while (< k i) | |
253 (setq sum (math-sub sum (math-mul (nth k (nth i lu)) | |
254 (nth j (nth k lu)))) | |
255 k (1+ k))) | |
256 (setcar (nthcdr j (nth i lu)) sum) | |
257 (setq i (1+ i))) | |
258 (while (<= i n) | |
259 (math-working "LUD step" (format "%d/%d" j i)) | |
260 (setq sum (nth j (nth i lu)) | |
261 k 1) | |
262 (while (< k j) | |
263 (setq sum (math-sub sum (math-mul (nth k (nth i lu)) | |
264 (nth j (nth k lu)))) | |
265 k (1+ k))) | |
266 (setcar (nthcdr j (nth i lu)) sum) | |
267 (let ((dum (math-abs-approx sum))) | |
268 (if (Math-lessp big dum) | |
269 (setq big dum | |
270 imax i))) | |
271 (setq i (1+ i))) | |
272 (if (> imax j) | |
273 (setq lu (math-swap-rows lu j imax) | |
274 d (- d))) | |
275 (setq index (cons imax index)) | |
276 (let ((pivot (nth j (nth j lu)))) | |
277 (if (math-zerop pivot) | |
278 (throw 'singular nil) | |
279 (setq i j) | |
280 (while (<= (setq i (1+ i)) n) | |
281 (setcar (nthcdr j (nth i lu)) | |
282 (math-div (nth j (nth i lu)) pivot))))) | |
283 (setq j (1+ j))) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
284 (list lu (nreverse index) d))) |
40785 | 285 |
286 (defun math-swap-rows (m r1 r2) | |
287 (or (= r1 r2) | |
288 (let* ((r1prev (nthcdr (1- r1) m)) | |
289 (row1 (cdr r1prev)) | |
290 (r2prev (nthcdr (1- r2) m)) | |
291 (row2 (cdr r2prev)) | |
292 (r2next (cdr row2))) | |
293 (setcdr r2prev row1) | |
294 (setcdr r1prev row2) | |
295 (setcdr row2 (cdr row1)) | |
296 (setcdr row1 r2next))) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
297 m) |
40785 | 298 |
299 | |
300 (defun math-lud-solve (lud b &optional need) | |
301 (if lud | |
302 (let* ((x (math-copy-matrix b)) | |
303 (n (1- (length x))) | |
304 (m (1- (length (nth 1 x)))) | |
305 (lu (car lud)) | |
306 (col 1) | |
307 i j ip ii index sum) | |
308 (while (<= col m) | |
309 (math-working "LUD solver step" col) | |
310 (setq i 1 | |
311 ii nil | |
312 index (nth 1 lud)) | |
313 (while (<= i n) | |
314 (setq ip (car index) | |
315 index (cdr index) | |
316 sum (nth col (nth ip x))) | |
317 (setcar (nthcdr col (nth ip x)) (nth col (nth i x))) | |
318 (if (null ii) | |
319 (or (math-zerop sum) | |
320 (setq ii i)) | |
321 (setq j ii) | |
322 (while (< j i) | |
323 (setq sum (math-sub sum (math-mul (nth j (nth i lu)) | |
324 (nth col (nth j x)))) | |
325 j (1+ j)))) | |
326 (setcar (nthcdr col (nth i x)) sum) | |
327 (setq i (1+ i))) | |
328 (while (>= (setq i (1- i)) 1) | |
329 (setq sum (nth col (nth i x)) | |
330 j i) | |
331 (while (<= (setq j (1+ j)) n) | |
332 (setq sum (math-sub sum (math-mul (nth j (nth i lu)) | |
333 (nth col (nth j x)))))) | |
334 (setcar (nthcdr col (nth i x)) | |
335 (math-div sum (nth i (nth i lu))))) | |
336 (setq col (1+ col))) | |
337 x) | |
338 (and need | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
339 (math-reject-arg need "*Singular matrix")))) |
40785 | 340 |
341 (defun calcFunc-lud (m) | |
342 (if (math-square-matrixp m) | |
343 (or (math-with-extra-prec 2 | |
344 (let ((lud (math-matrix-lud m))) | |
345 (and lud | |
346 (let* ((lmat (math-copy-matrix (car lud))) | |
347 (umat (math-copy-matrix (car lud))) | |
348 (n (1- (length (car lud)))) | |
349 (perm (calcFunc-idn 1 n)) | |
350 i (j 1)) | |
351 (while (<= j n) | |
352 (setq i 1) | |
353 (while (< i j) | |
354 (setcar (nthcdr j (nth i lmat)) 0) | |
355 (setq i (1+ i))) | |
356 (setcar (nthcdr j (nth j lmat)) 1) | |
357 (while (<= (setq i (1+ i)) n) | |
358 (setcar (nthcdr j (nth i umat)) 0)) | |
359 (setq j (1+ j))) | |
360 (while (>= (setq j (1- j)) 1) | |
361 (let ((pos (nth (1- j) (nth 1 lud)))) | |
362 (or (= pos j) | |
363 (setq perm (math-swap-rows perm j pos))))) | |
364 (list 'vec perm lmat umat))))) | |
365 (math-reject-arg m "*Singular matrix")) | |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
366 (math-reject-arg m 'square-matrixp))) |
40785 | 367 |
58666
b21f320635a7
Add a provide statement.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
58548
diff
changeset
|
368 (provide 'calc-mtx) |
b21f320635a7
Add a provide statement.
Jay Belanger <jay.p.belanger@gmail.com>
parents:
58548
diff
changeset
|
369 |
52401 | 370 ;;; arch-tag: fc0947b1-90e1-4a23-8950-d8ead9c3a306 |
41047
73f364fd8aaa
Style cleanup; don't put closing parens on their
Colin Walters <walters@gnu.org>
parents:
40785
diff
changeset
|
371 ;;; calc-mtx.el ends here |