[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-Need-calc-mat () nil)


(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)))

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