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

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

;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc.

;; Author: David Gillespie <daveg@synaptics.com>
;; Maintainer: Jay Belanger <belanger@truman.edu>

;; This file is part of GNU Emacs.

;; GNU Emacs is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY.  No author or distributor
;; accepts responsibility to anyone for the consequences of using it
;; or for whether it serves any particular purpose or works at all,
;; unless he says so in writing.  Refer to the GNU Emacs General Public
;; License for full details.

;; Everyone is granted permission to copy, modify and redistribute
;; GNU Emacs, but only under the conditions described in the
;; GNU Emacs General Public License.   A copy of this license is
;; supposed to have been given to you along with GNU Emacs so you
;; can know your rights and responsibilities.  It should be in a
;; file named COPYING.  Among other things, the copyright notice
;; and this notice must be preserved on all copies.

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