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

;; Calculator for GNU Emacs, part II [calc-mat.el]
;; Copyright (C) 1990, 1991, 1992, 1993 Free Software Foundation, Inc.
;; Written by Dave Gillespie, daveg@synaptics.com.

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


;; 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.
(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)))
)
(defvar math-lud-cache nil)

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