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

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

;; Copyright (C) 1990-1993, 2001-2013 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))))


(defun math-lud-pivot-check (a)
  "Determine a useful value for checking the size of potential pivots
in LUD decomposition."
  (cond ((eq (car-safe a) 'mod)
         (if (and (math-integerp (nth 1 a))
                  (math-integerp (nth 2 a))
                  (eq (math-gcd (nth 1 a) (nth 2 a)) 1))
             1
           0))
        (t
         (math-abs-approx a))))


;;; 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-lud-pivot-check 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)

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