/[emacs]/emacs/lisp/calc/calc-mtx.el
ViewVC logotype

Diff of /emacs/lisp/calc/calc-mtx.el

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by eliz, Tue Nov 6 18:59:06 2001 UTC revision 1.2 by walters, Wed Nov 14 09:06:29 2001 UTC
# Line 1  Line 1 
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.
# Line 32  Line 32 
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]
# Line 53  Line 50 
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    
# Line 82  Line 77 
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.
# Line 167  Line 158 
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)
# Line 177  Line 167 
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))))
# Line 217  Line 206 
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.
# Line 238  Line 225 
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.
# Line 288  Line 274 
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)
# Line 302  Line 287 
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)
# Line 345  Line 329 
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)
# Line 373  Line 356 
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

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26