[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>
;; Maintainers: D. Goel <deego@gnufans.org>
;;              Colin Walters <walters@debian.org>

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

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