[bpt/emacs.git] / lisp / calc / calc-mtx.el

;;; calc-mtx.el --- matrix functions for Calc

;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004,
;;   2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.

;; Author: David Gillespie <daveg@synaptics.com>
;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>

;; This file is part of GNU Emacs.

;; GNU Emacs is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.

;; GNU Emacs is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.

;; You should have received a copy of the GNU General Public License
;; along with GNU Emacs.  If not, see <http://www.gnu.org/licenses/>.

;;; Commentary:

;;; Code:

;; This file is autoloaded from calc-ext.el.

(require 'calc-ext)
(require 'calc-macs)

(defun calc-mdet (arg)
  (interactive "P")
  (calc-slow-wrapper
   (calc-unary-op "mdet" 'calcFunc-det arg)))

(defun calc-mtrace (arg)
  (interactive "P")
  (calc-slow-wrapper
   (calc-unary-op "mtr" 'calcFunc-tr arg)))

(defun calc-mlud (arg)
  (interactive "P")
  (calc-slow-wrapper
   (calc-unary-op "mlud" 'calcFunc-lud arg)))


;;; Coerce row vector A to be a matrix.  [V V]
(defun math-row-matrix (a)
  (if (and (Math-vectorp a)
	   (not (math-matrixp a)))
      (list 'vec a)
    a))

;;; Coerce column vector A to be a matrix.  [V V]
(defun math-col-matrix (a)
  (if (and (Math-vectorp a)
	   (not (math-matrixp a)))
      (cons 'vec (mapcar (function (lambda (x) (list 'vec x))) (cdr a)))
    a))


;;; Multiply matrices A and B.  [V V V]
(defun math-mul-mats (a b)
  (let ((mat nil)
	(cols (length (nth 1 b)))
	row col ap bp accum)
    (while (setq a (cdr a))
      (setq col cols
	    row nil)
      (while (> (setq col (1- col)) 0)
	(setq ap (cdr (car a))
	      bp (cdr b)
	      accum (math-mul (car ap) (nth col (car bp))))
	(while (setq ap (cdr ap) bp (cdr bp))
	  (setq accum (math-add accum (math-mul (car ap) (nth col (car bp))))))
	(setq row (cons accum row)))
      (setq mat (cons (cons 'vec row) mat)))
    (cons 'vec (nreverse mat))))

(defun math-mul-mat-vec (a b)
  (cons 'vec (mapcar (function (lambda (row)
				 (math-dot-product row b)))
		     (cdr a))))


(defun calcFunc-tr (mat)   ; [Public]
  (if (math-square-matrixp mat)
      (math-matrix-trace-step 2 (1- (length mat)) mat (nth 1 (nth 1 mat)))
    (math-reject-arg mat 'square-matrixp)))

(defun math-matrix-trace-step (n size mat sum)
  (if (<= n size)
      (math-matrix-trace-step (1+ n) size mat
			      (math-add sum (nth n (nth n mat))))
    sum))


;;; Matrix inverse and determinant.
(defun math-matrix-inv-raw (m)
  (let ((n (1- (length m))))
    (if (<= n 3)
	(let ((det (math-det-raw m)))
	  (and (not (math-zerop det))
	       (math-div
		(cond ((= n 1) 1)
		      ((= n 2)
		       (list 'vec
			     (list 'vec
				   (nth 2 (nth 2 m))
				   (math-neg (nth 2 (nth 1 m))))
			     (list 'vec
				   (math-neg (nth 1 (nth 2 m)))
				   (nth 1 (nth 1 m)))))
		      ((= n 3)
		       (list 'vec
			     (list 'vec
				   (math-sub (math-mul (nth 3 (nth 3 m))
						       (nth 2 (nth 2 m)))
					     (math-mul (nth 3 (nth 2 m))
						       (nth 2 (nth 3 m))))
				   (math-sub (math-mul (nth 3 (nth 1 m))
						       (nth 2 (nth 3 m)))
					     (math-mul (nth 3 (nth 3 m))
						       (nth 2 (nth 1 m))))
				   (math-sub (math-mul (nth 3 (nth 2 m))
						       (nth 2 (nth 1 m)))
					     (math-mul (nth 3 (nth 1 m))
						       (nth 2 (nth 2 m)))))
			     (list 'vec
				   (math-sub (math-mul (nth 3 (nth 2 m))
						       (nth 1 (nth 3 m)))
					     (math-mul (nth 3 (nth 3 m))
						       (nth 1 (nth 2 m))))
				   (math-sub (math-mul (nth 3 (nth 3 m))
						       (nth 1 (nth 1 m)))
					     (math-mul (nth 3 (nth 1 m))
						       (nth 1 (nth 3 m))))
				   (math-sub (math-mul (nth 3 (nth 1 m))
						       (nth 1 (nth 2 m)))
					     (math-mul (nth 3 (nth 2 m))
						       (nth 1 (nth 1 m)))))
			     (list 'vec
				   (math-sub (math-mul (nth 2 (nth 3 m))
						       (nth 1 (nth 2 m)))
					     (math-mul (nth 2 (nth 2 m))
						       (nth 1 (nth 3 m))))
				   (math-sub (math-mul (nth 2 (nth 1 m))
						       (nth 1 (nth 3 m)))
					     (math-mul (nth 2 (nth 3 m))
						       (nth 1 (nth 1 m))))
				   (math-sub (math-mul (nth 2 (nth 2 m))
						       (nth 1 (nth 1 m)))
					     (math-mul (nth 2 (nth 1 m))
						       (nth 1 (nth 2 m))))))))
		det)))
      (let ((lud (math-matrix-lud m)))
	(and lud
	     (math-lud-solve lud (calcFunc-idn 1 n)))))))

(defun calcFunc-det (m)
  (if (math-square-matrixp m)
      (math-with-extra-prec 2 (math-det-raw m))
    (if (and (eq (car-safe m) 'calcFunc-idn)
	     (or (math-zerop (nth 1 m))
		 (math-equal-int (nth 1 m) 1)))
	(nth 1 m)
      (math-reject-arg m 'square-matrixp))))

;; The variable math-det-lu is local to math-det-raw, but is
;; used by math-det-step, which is called by math-det-raw.
(defvar math-det-lu)

(defun math-det-raw (m)
  (let ((n (1- (length m))))
    (cond ((= n 1)
	   (nth 1 (nth 1 m)))
	  ((= n 2)
	   (math-sub (math-mul (nth 1 (nth 1 m))
			       (nth 2 (nth 2 m)))
		     (math-mul (nth 2 (nth 1 m))
			       (nth 1 (nth 2 m)))))
	  ((= n 3)
	   (math-sub
	    (math-sub
	     (math-sub
	      (math-add
	       (math-add
		(math-mul (nth 1 (nth 1 m))
			  (math-mul (nth 2 (nth 2 m))
				    (nth 3 (nth 3 m))))
		(math-mul (nth 2 (nth 1 m))
			  (math-mul (nth 3 (nth 2 m))
				    (nth 1 (nth 3 m)))))
	       (math-mul (nth 3 (nth 1 m))
			 (math-mul (nth 1 (nth 2 m))
				   (nth 2 (nth 3 m)))))
	      (math-mul (nth 3 (nth 1 m))
			(math-mul (nth 2 (nth 2 m))
				  (nth 1 (nth 3 m)))))
	     (math-mul (nth 1 (nth 1 m))
		       (math-mul (nth 3 (nth 2 m))
				 (nth 2 (nth 3 m)))))
	    (math-mul (nth 2 (nth 1 m))
		      (math-mul (nth 1 (nth 2 m))
				(nth 3 (nth 3 m))))))
	  (t (let ((lud (math-matrix-lud m)))
	       (if lud
		   (let ((math-det-lu (car lud)))
		     (math-det-step n (nth 2 lud)))
		 0))))))

(defun math-det-step (n prod)
  (if (> n 0)
      (math-det-step (1- n) (math-mul prod (nth n (nth n math-det-lu))))
    prod))

;;; This returns a list (LU index d), or nil if not possible.
;;; Argument M must be a square matrix.
(defvar math-lud-cache nil)
(defun math-matrix-lud (m)
  (let ((old (assoc m math-lud-cache))
	(context (list calc-internal-prec calc-prefer-frac)))
    (if (and old (equal (nth 1 old) context))
	(cdr (cdr old))
      (let* ((lud (catch 'singular (math-do-matrix-lud m)))
	     (entry (cons context lud)))
	(if old
	    (setcdr old entry)
	  (setq math-lud-cache (cons (cons m entry) math-lud-cache)))
	lud))))

;;; Numerical Recipes section 2.3; implicit pivoting omitted.
(defun math-do-matrix-lud (m)
  (let* ((lu (math-copy-matrix m))
	 (n (1- (length lu)))
	 i (j 1) k imax sum big
	 (d 1) (index nil))
    (while (<= j n)
      (setq i 1
	    big 0
	    imax j)
      (while (< i j)
	(math-working "LUD step" (format "%d/%d" j i))
	(setq sum (nth j (nth i lu))
	      k 1)
	(while (< k i)
	  (setq sum (math-sub sum (math-mul (nth k (nth i lu))
					    (nth j (nth k lu))))
		k (1+ k)))
	(setcar (nthcdr j (nth i lu)) sum)
	(setq i (1+ i)))
      (while (<= i n)
	(math-working "LUD step" (format "%d/%d" j i))
	(setq sum (nth j (nth i lu))
	      k 1)
	(while (< k j)
	  (setq sum (math-sub sum (math-mul (nth k (nth i lu))
					    (nth j (nth k lu))))
		k (1+ k)))
	(setcar (nthcdr j (nth i lu)) sum)
	(let ((dum (math-abs-approx sum)))
	  (if (Math-lessp big dum)
	      (setq big dum
		    imax i)))
	(setq i (1+ i)))
      (if (> imax j)
	  (setq lu (math-swap-rows lu j imax)
		d (- d)))
      (setq index (cons imax index))
      (let ((pivot (nth j (nth j lu))))
	(if (math-zerop pivot)
	    (throw 'singular nil)
	  (setq i j)
	  (while (<= (setq i (1+ i)) n)
	    (setcar (nthcdr j (nth i lu))
		    (math-div (nth j (nth i lu)) pivot)))))
      (setq j (1+ j)))
    (list lu (nreverse index) d)))

(defun math-swap-rows (m r1 r2)
  (or (= r1 r2)
      (let* ((r1prev (nthcdr (1- r1) m))
	     (row1 (cdr r1prev))
	     (r2prev (nthcdr (1- r2) m))
	     (row2 (cdr r2prev))
	     (r2next (cdr row2)))
	(setcdr r2prev row1)
	(setcdr r1prev row2)
	(setcdr row2 (cdr row1))
	(setcdr row1 r2next)))
  m)


(defun math-lud-solve (lud b &optional need)
  (if lud
      (let* ((x (math-copy-matrix b))
	     (n (1- (length x)))
	     (m (1- (length (nth 1 x))))
	     (lu (car lud))
	     (col 1)
	     i j ip ii index sum)
	(while (<= col m)
	  (math-working "LUD solver step" col)
	  (setq i 1
		ii nil
		index (nth 1 lud))
	  (while (<= i n)
	    (setq ip (car index)
		  index (cdr index)
		  sum (nth col (nth ip x)))
	    (setcar (nthcdr col (nth ip x)) (nth col (nth i x)))
	    (if (null ii)
		(or (math-zerop sum)
		    (setq ii i))
	      (setq j ii)
	      (while (< j i)
		(setq sum (math-sub sum (math-mul (nth j (nth i lu))
						  (nth col (nth j x))))
		      j (1+ j))))
	    (setcar (nthcdr col (nth i x)) sum)
	    (setq i (1+ i)))
	  (while (>= (setq i (1- i)) 1)
	    (setq sum (nth col (nth i x))
		  j i)
	    (while (<= (setq j (1+ j)) n)
	      (setq sum (math-sub sum (math-mul (nth j (nth i lu))
						(nth col (nth j x))))))
	    (setcar (nthcdr col (nth i x))
		    (math-div sum (nth i (nth i lu)))))
	  (setq col (1+ col)))
	x)
    (and need
	 (math-reject-arg need "*Singular matrix"))))

(defun calcFunc-lud (m)
  (if (math-square-matrixp m)
      (or (math-with-extra-prec 2
	    (let ((lud (math-matrix-lud m)))
	      (and lud
		   (let* ((lmat (math-copy-matrix (car lud)))
			  (umat (math-copy-matrix (car lud)))
			  (n (1- (length (car lud))))
			  (perm (calcFunc-idn 1 n))
			  i (j 1))
		     (while (<= j n)
		       (setq i 1)
		       (while (< i j)
			 (setcar (nthcdr j (nth i lmat)) 0)
			 (setq i (1+ i)))
		       (setcar (nthcdr j (nth j lmat)) 1)
		       (while (<= (setq i (1+ i)) n)
			 (setcar (nthcdr j (nth i umat)) 0))
		       (setq j (1+ j)))
		     (while (>= (setq j (1- j)) 1)
		       (let ((pos (nth (1- j) (nth 1 lud))))
			 (or (= pos j)
			     (setq perm (math-swap-rows perm j pos)))))
		     (list 'vec perm lmat umat)))))
	  (math-reject-arg m "*Singular matrix"))
    (math-reject-arg m 'square-matrixp)))

(provide 'calc-mtx)

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