1 |
;; Calculator for GNU Emacs, part II [calc-mat.el] |
;; Calculator for GNU Emacs, part II [calc-mat.el] |
2 |
;; Copyright (C) 1990, 1991, 1992, 1993 Free Software Foundation, Inc. |
;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc. |
3 |
;; Written by Dave Gillespie, daveg@synaptics.com. |
;; Written by Dave Gillespie, daveg@synaptics.com. |
4 |
|
|
5 |
;; This file is part of GNU Emacs. |
;; This file is part of GNU Emacs. |
32 |
(defun calc-mdet (arg) |
(defun calc-mdet (arg) |
33 |
(interactive "P") |
(interactive "P") |
34 |
(calc-slow-wrapper |
(calc-slow-wrapper |
35 |
(calc-unary-op "mdet" 'calcFunc-det arg)) |
(calc-unary-op "mdet" 'calcFunc-det arg))) |
|
) |
|
36 |
|
|
37 |
(defun calc-mtrace (arg) |
(defun calc-mtrace (arg) |
38 |
(interactive "P") |
(interactive "P") |
39 |
(calc-slow-wrapper |
(calc-slow-wrapper |
40 |
(calc-unary-op "mtr" 'calcFunc-tr arg)) |
(calc-unary-op "mtr" 'calcFunc-tr arg))) |
|
) |
|
41 |
|
|
42 |
(defun calc-mlud (arg) |
(defun calc-mlud (arg) |
43 |
(interactive "P") |
(interactive "P") |
44 |
(calc-slow-wrapper |
(calc-slow-wrapper |
45 |
(calc-unary-op "mlud" 'calcFunc-lud arg)) |
(calc-unary-op "mlud" 'calcFunc-lud arg))) |
|
) |
|
46 |
|
|
47 |
|
|
48 |
;;; Coerce row vector A to be a matrix. [V V] |
;;; Coerce row vector A to be a matrix. [V V] |
50 |
(if (and (Math-vectorp a) |
(if (and (Math-vectorp a) |
51 |
(not (math-matrixp a))) |
(not (math-matrixp a))) |
52 |
(list 'vec a) |
(list 'vec a) |
53 |
a) |
a)) |
|
) |
|
54 |
|
|
55 |
;;; Coerce column vector A to be a matrix. [V V] |
;;; Coerce column vector A to be a matrix. [V V] |
56 |
(defun math-col-matrix (a) |
(defun math-col-matrix (a) |
57 |
(if (and (Math-vectorp a) |
(if (and (Math-vectorp a) |
58 |
(not (math-matrixp a))) |
(not (math-matrixp a))) |
59 |
(cons 'vec (mapcar (function (lambda (x) (list 'vec x))) (cdr a))) |
(cons 'vec (mapcar (function (lambda (x) (list 'vec x))) (cdr a))) |
60 |
a) |
a)) |
|
) |
|
61 |
|
|
62 |
|
|
63 |
|
|
77 |
(setq accum (math-add accum (math-mul (car ap) (nth col (car bp)))))) |
(setq accum (math-add accum (math-mul (car ap) (nth col (car bp)))))) |
78 |
(setq row (cons accum row))) |
(setq row (cons accum row))) |
79 |
(setq mat (cons (cons 'vec row) mat))) |
(setq mat (cons (cons 'vec row) mat))) |
80 |
(cons 'vec (nreverse mat))) |
(cons 'vec (nreverse mat)))) |
|
) |
|
81 |
|
|
82 |
(defun math-mul-mat-vec (a b) |
(defun math-mul-mat-vec (a b) |
83 |
(cons 'vec (mapcar (function (lambda (row) |
(cons 'vec (mapcar (function (lambda (row) |
84 |
(math-dot-product row b))) |
(math-dot-product row b))) |
85 |
(cdr a))) |
(cdr a)))) |
|
) |
|
86 |
|
|
87 |
|
|
88 |
|
|
89 |
(defun calcFunc-tr (mat) ; [Public] |
(defun calcFunc-tr (mat) ; [Public] |
90 |
(if (math-square-matrixp mat) |
(if (math-square-matrixp mat) |
91 |
(math-matrix-trace-step 2 (1- (length mat)) mat (nth 1 (nth 1 mat))) |
(math-matrix-trace-step 2 (1- (length mat)) mat (nth 1 (nth 1 mat))) |
92 |
(math-reject-arg mat 'square-matrixp)) |
(math-reject-arg mat 'square-matrixp))) |
|
) |
|
93 |
|
|
94 |
(defun math-matrix-trace-step (n size mat sum) |
(defun math-matrix-trace-step (n size mat sum) |
95 |
(if (<= n size) |
(if (<= n size) |
96 |
(math-matrix-trace-step (1+ n) size mat |
(math-matrix-trace-step (1+ n) size mat |
97 |
(math-add sum (nth n (nth n mat)))) |
(math-add sum (nth n (nth n mat)))) |
98 |
sum) |
sum)) |
|
) |
|
99 |
|
|
100 |
|
|
101 |
;;; Matrix inverse and determinant. |
;;; Matrix inverse and determinant. |
158 |
det))) |
det))) |
159 |
(let ((lud (math-matrix-lud m))) |
(let ((lud (math-matrix-lud m))) |
160 |
(and lud |
(and lud |
161 |
(math-lud-solve lud (calcFunc-idn 1 n)))))) |
(math-lud-solve lud (calcFunc-idn 1 n))))))) |
|
) |
|
162 |
|
|
163 |
(defun calcFunc-det (m) |
(defun calcFunc-det (m) |
164 |
(if (math-square-matrixp m) |
(if (math-square-matrixp m) |
167 |
(or (math-zerop (nth 1 m)) |
(or (math-zerop (nth 1 m)) |
168 |
(math-equal-int (nth 1 m) 1))) |
(math-equal-int (nth 1 m) 1))) |
169 |
(nth 1 m) |
(nth 1 m) |
170 |
(math-reject-arg m 'square-matrixp))) |
(math-reject-arg m 'square-matrixp)))) |
|
) |
|
171 |
|
|
172 |
(defun math-det-raw (m) |
(defun math-det-raw (m) |
173 |
(let ((n (1- (length m)))) |
(let ((n (1- (length m)))) |
206 |
(if lud |
(if lud |
207 |
(let ((lu (car lud))) |
(let ((lu (car lud))) |
208 |
(math-det-step n (nth 2 lud))) |
(math-det-step n (nth 2 lud))) |
209 |
0))))) |
0)))))) |
|
) |
|
210 |
|
|
211 |
(defun math-det-step (n prod) |
(defun math-det-step (n prod) |
212 |
(if (> n 0) |
(if (> n 0) |
213 |
(math-det-step (1- n) (math-mul prod (nth n (nth n lu)))) |
(math-det-step (1- n) (math-mul prod (nth n (nth n lu)))) |
214 |
prod) |
prod)) |
|
) |
|
215 |
|
|
216 |
;;; This returns a list (LU index d), or NIL if not possible. |
;;; This returns a list (LU index d), or NIL if not possible. |
217 |
;;; Argument M must be a square matrix. |
;;; Argument M must be a square matrix. |
225 |
(if old |
(if old |
226 |
(setcdr old entry) |
(setcdr old entry) |
227 |
(setq math-lud-cache (cons (cons m entry) math-lud-cache))) |
(setq math-lud-cache (cons (cons m entry) math-lud-cache))) |
228 |
lud))) |
lud)))) |
|
) |
|
229 |
(defvar math-lud-cache nil) |
(defvar math-lud-cache nil) |
230 |
|
|
231 |
;;; Numerical Recipes section 2.3; implicit pivoting omitted. |
;;; Numerical Recipes section 2.3; implicit pivoting omitted. |
274 |
(setcar (nthcdr j (nth i lu)) |
(setcar (nthcdr j (nth i lu)) |
275 |
(math-div (nth j (nth i lu)) pivot))))) |
(math-div (nth j (nth i lu)) pivot))))) |
276 |
(setq j (1+ j))) |
(setq j (1+ j))) |
277 |
(list lu (nreverse index) d)) |
(list lu (nreverse index) d))) |
|
) |
|
278 |
|
|
279 |
(defun math-swap-rows (m r1 r2) |
(defun math-swap-rows (m r1 r2) |
280 |
(or (= r1 r2) |
(or (= r1 r2) |
287 |
(setcdr r1prev row2) |
(setcdr r1prev row2) |
288 |
(setcdr row2 (cdr row1)) |
(setcdr row2 (cdr row1)) |
289 |
(setcdr row1 r2next))) |
(setcdr row1 r2next))) |
290 |
m |
m) |
|
) |
|
291 |
|
|
292 |
|
|
293 |
(defun math-lud-solve (lud b &optional need) |
(defun math-lud-solve (lud b &optional need) |
329 |
(setq col (1+ col))) |
(setq col (1+ col))) |
330 |
x) |
x) |
331 |
(and need |
(and need |
332 |
(math-reject-arg need "*Singular matrix"))) |
(math-reject-arg need "*Singular matrix")))) |
|
) |
|
333 |
|
|
334 |
(defun calcFunc-lud (m) |
(defun calcFunc-lud (m) |
335 |
(if (math-square-matrixp m) |
(if (math-square-matrixp m) |
356 |
(setq perm (math-swap-rows perm j pos))))) |
(setq perm (math-swap-rows perm j pos))))) |
357 |
(list 'vec perm lmat umat))))) |
(list 'vec perm lmat umat))))) |
358 |
(math-reject-arg m "*Singular matrix")) |
(math-reject-arg m "*Singular matrix")) |
359 |
(math-reject-arg m 'square-matrixp)) |
(math-reject-arg m 'square-matrixp))) |
|
) |
|
360 |
|
|
361 |
|
;;; calc-mtx.el ends here |