1;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3;;;     The data in this file contains enhancments.                    ;;;;;
4;;;                                                                    ;;;;;
5;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
6;;;     All rights reserved                                            ;;;;;
7;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
8;;;     (c) Copyright 1982 Massachusetts Institute of Technology         ;;;
9;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
10
11(in-package :maxima)
12
13(macsyma-module matrix)
14
15(declare-top (special *ech* *tri* *inv*
16		      mdl $detout vlist mul* top* *det* genvar $ratfac
17		      varlist header $scalarmatrixp $sparse
18		      $algebraic *rank* *mat*))
19
20(defmvar $detout nil)
21(defmvar top* nil)
22(defmvar $ratmx nil)
23(defmvar $matrix_element_mult "*")  ;;; Else, most useful when "."
24(defmvar $matrix_element_add "+")
25(defmvar $matrix_element_transpose nil)
26
27(defvar *mat*)
28
29;;I believe that all the code now stores arrays in the value cell
30(defun get-array-pointer (symbol)
31  "There may be nesting of functions and we may well need to apply
32   this twice in a row"
33  (if (arrayp symbol) symbol (symbol-value symbol)))
34
35(defun mxc (x)
36  (mapcar #'(lambda (y) (cons '(mlist) y)) x)) ; Matrix to MACSYMA conversion
37
38(defun mcx (x)
39  (mapcar #'cdr x))			; MACSYMA to Matrix conversion
40
41;; Transpose a list of lists ll. Example: ((1 2) (3 4)) --> ((1 3) (2 4)).
42
43(defun transpose (ll)
44  (let ((acc nil))
45    (while (car ll)
46      (push (mapcar #'car ll) acc)
47      (setq ll (mapcar #'cdr ll)))
48    (nreverse acc)))
49
50(defun nthcol (x nn)
51  (if (or (null x) (> nn (length (car x))))
52      nil
53      (nthcol1 x nn)))
54
55(defun nthcol1 (x nn)
56  (if (or (null x) (= nn 0))
57      nil
58      (cons (ith (car x) nn) (nthcol1 (cdr x) nn))))
59
60;; MAYBE THIS FUNCTION SHOULD HAVE AN ARGUMENT TO INDICATE WHO CALLED IT (TO SMASH INTO ERROR MESSAGES)
61(defun check (x)
62  (cond ((atom x) (merror (intl:gettext "not a matrix: ~M") x))
63	((eq (caar x) '$matrix) x)
64	((eq (caar x) 'mlist) (list '($matrix) x))
65	(t (merror (intl:gettext "not a matrix: ~M") x))))
66
67(defun check1 (x)
68  (cond ((atom x) nil)
69	((eq (caar x) '$matrix) x)
70	((eq (caar x) 'mlist) (list '($matrix) x))))
71
72(defmfun $matrixp (x)
73  (and (not (atom x)) (eq (caar x) '$matrix)))
74
75(defmfun $charpoly (mat var)
76  (setq mat (check mat))
77  (unless (= (length mat) (length (cadr mat)))
78    (merror (intl:gettext "charpoly: matrix must be square; found ~M rows, ~M columns.") (length mat) (length (cadr mat))))
79  (cond ((not $ratmx)
80	 (det1 (addmatrix1
81		(setq mat (mcx (cdr mat)))
82		(diagmatrix (length mat) (list '(mtimes) -1 var) '$charpoly))))
83	(t (newvar var) (newvarmat1 mat)
84	   (setq mat (mcx (cdr mat)))
85	   (determinant1 (addmatrix mat (diagmatrix (length mat)
86						    (list '(mtimes) -1 var)
87						    '$charpoly))))))
88
89(defun disreplist1 (a)
90  (setq header (list 'mrat 'simp varlist genvar))
91  (mapcar #'disreplist a))
92
93(defun disreplist (a)
94  (mapcar #'(lambda (e) (cons header e)) a))
95
96(defun replist1 (a)
97  (mapcar #'replist a))
98
99(defun replist (a)
100  (mapcar #'(lambda (e) (cdr (ratrep* e))) a))
101
102(defun timex (mat1 mat2)
103  (cond ((equal mat1 1) mat2)
104	((and ($matrixp mat1) ($matrixp mat2) (null (cdr mat1)))
105	 (ncons '($matrix simp)))
106	(t (newvarmat mat1 mat2)
107	   (let (($scalarmatrixp
108		  (if (and ($listp mat1) ($listp mat2)) t $scalarmatrixp)))
109	     (simplifya (timex0 mat1 mat2) nil)))))
110
111(defun lnewvar (a)
112  (let ((vlist nil))
113    (lnewvar1 a)
114    (setq varlist (nconc (sortgreat vlist) varlist))))
115
116(defun lnewvar1 (a)
117  (cond ((atom a) (newvar1 a))
118	((mbagp a) (mapc #'lnewvar1 (cdr a)))
119	(t (newvar1 a))))
120
121(defun newvarmat (mat1 mat2)
122  (when $ratmx
123    (let ((vlist nil))
124      (lnewvar1 mat1) (lnewvar1 mat2)
125      (setq varlist (nconc (sortgreat vlist) varlist)))))
126
127(defun newvarmat1 (a)
128  (cond ($ratmx (lnewvar a))))
129
130(defun addmatrix (x y)
131  (setq x (replist1 x) y (replist1 y))
132  (disreplist1 (addmatrix1 x y)))
133
134(defun addmatrix1 (b c)
135  (unless (and (= (length b) (length c))
136	       (= (length (car b)) (length (car c))))
137    (merror (intl:gettext "ADDMATRIX1: attempt to add nonconformable matrices.")))
138  (mapcar #'addrows b c))
139
140(defun addrows (a b)
141  (if (not $ratmx)
142      (mapcar #'(lambda (i j) (simplus (list '(mplus) i j) 1 nil)) a b)
143      (mapcar #'ratplus a b)))
144
145(defmfun $determinant (mat)
146  (cond ((not (or (mbagp mat) ($matrixp mat)))
147         (if ($scalarp mat) mat (list '(%determinant) mat)))
148	(t (setq mat (check mat))
149	   (unless  (= (length mat) (length (cadr mat)))
150             (merror
151               (intl:gettext
152                 "determinant: matrix must be square; found ~M rows, ~M columns.")
153               (length (cdr mat))
154               (length (cdadr mat))))
155           (cond ((not $ratmx) (det1 (mcx (cdr mat))))
156	         (t (newvarmat1 mat) (determinant1 (mcx (cdr mat))))))))
157
158(defun det (m)
159  (if (= (length m) 1)
160      (caar m)
161      (let (*det* mul*)
162	(mtoa '*mat* (setq *det* (length m)) *det* m)
163	(setq *det* (tfgeli0 '*mat* *det* *det*))
164	(ratreduce *det* mul*))))
165
166(defun determinant1 (x)
167  (catch 'dz (rdis (det (replist1 x)))))
168
169(defun treedet (mat)
170  (prog (row mdl lindex tuplel n id md lt)
171     (setq mat (reverse mat))
172     (setq n (length mat) md (car mat))
173     (setq mat (cdr mat)) (setq lindex (nreverse (index* n)) tuplel (mapcar #'list lindex))
174     loop1 (when (null mat) (return (car md)))
175     (setq mdl nil)
176     (mapcar #'(lambda(a b) (setq mdl (nconc mdl (list a b)))) tuplel md)
177     (setq md nil)
178     (setq row (car mat)
179	   mat (cdr mat))
180     (setq lt (setq tuplel (nextlevel tuplel lindex)))
181     loop2 (when (null lt)
182	     (setq md (nreverse md))
183	     (go loop1))
184     (setq id (car lt) lt (cdr lt))
185     (setq md (cons (compumd id row) md))
186     (go loop2)))
187
188(defun assoo (e l)
189  (prog ()
190   loop (cond ((null l) (return nil))
191	      ((equal e (car l)) (return (cadr l))))
192   (setq l (cddr l))
193   (go loop)))
194
195(defun compumd (id row)
196  (prog (e minor i d sign ans)
197     (setq ans 0 sign -1 i id)
198     loop (when (null i) (return ans))
199     (setq d (car i) i (cdr i) sign (* -1 sign))
200     (cond ((equal (setq e (ith row d)) 0)
201	    (go loop))
202	   ((equal (setq minor (assoo (delete d (copy-tree id) :test #'equal) mdl)) 0)
203	    (go loop)))
204     (setq ans
205	   (if (and (equal $matrix_element_mult "*")
206		    (equal $matrix_element_add "+"))
207	       (add ans (mul sign e minor)) ;fast common case
208	     (mapply $matrix_element_add
209		     (list ans
210			   (mapply $matrix_element_mult
211				   (list sign e minor)
212				   $matrix_element_mult))
213		     $matrix_element_add)))
214     (go loop)))
215
216(defun apdl (l1 l2)
217  (mapcar #'(lambda (j) (append l1 (list j))) l2))
218
219(defun nextlevel (tuplel lindex)
220  (prog (ans l li)
221   loop (when (null tuplel) (return ans))
222   (setq l (car tuplel)
223	 tuplel (cdr tuplel)
224	 li (cdr (nthcdr (1- (car (last l))) lindex)))
225   (when (null li) (go loop))
226   (setq ans (nconc ans (apdl l li)))
227   (go loop)))
228
229(defun det1 (x)
230  (cond ($sparse (mtoa '*mat* (length x) (length x)
231		       (mapcar #'(lambda (x) (mapcar #'(lambda (y) (ncons y)) x))x))
232		 (sprdet '*mat* (length x)))
233	(t (treedet x))))
234
235(defmfun $ident (n)
236  (cons '($matrix) (mxc (diagmatrix n 1 '$ident))))
237
238(defmfun $diagmatrix (n var)
239  (cons '($matrix) (mxc (diagmatrix n var '$diagmatrix))))
240
241(defun diagmatrix (n var fn)
242  (prog (i ans)
243     (when (or (not (fixnump n)) (minusp n))
244       (improper-arg-err n fn))
245     (setq i n)
246     loop (if (zerop i) (return ans))
247     (setq ans (cons (onen i n var 0) ans) i (1- i))
248     (go loop)))
249
250;; ATOMAT GENERATES A MATRIX FROM A MXN ARRAY BY TAKING COLUMNS S TO N
251
252(defun atomat (name m n s)
253  (setq name (get-array-pointer name))
254  (prog (j d row mat)
255     (incf m)
256     (incf n)
257     loop1 (when (= m 1) (return mat))
258     (decf m)
259     (setq j n)
260     loop2 (when (= j s)
261	     (push row mat)
262	     (setq row nil)
263	     (go loop1))
264     (decf j)
265     (setq d (if top*
266		 (meval (list (list name 'array) m j))
267		 (aref name m j)))
268     (push (or d '(0 . 1)) row)
269     (go loop2)))
270
271(defmfun $invert_by_gausselim (k)
272  (let ((*inv* t) *det* top* mul* ($ratmx t) (ratmx $ratmx) $ratfac $sparse)
273    (cond ((atom k) ($nounify '$inverx) (list '(%inverx) k))
274	  (t (newvarmat1 (setq k (check k)))
275	     (setq k (invert1 (replist1 (mcx (cdr k)))))
276	     (setq k (cond ($detout `((mtimes)
277				      ((mexpt) ,(rdis (or *det* '(1 . 1))) -1)
278				      (($matrix) ,@(mxc (disreplist1 k)))))
279			   (t (cons '($matrix) (mxc (disreplist1 k))))))
280	     (cond ((and ratmx (not $detout))
281		    (fmapl1 #'(lambda (x) x) k))
282		   ((not ratmx) ($totaldisrep k))
283		   (t k))))))
284
285(defun diaginv (ax m)
286  (setq ax (get-array-pointer ax))
287  (cond ($detout (setq *det* 1)
288		 (do ((i 1 (1+ i)))
289		     ((> i m))
290		   (setq *det* (plcm *det* (car (aref ax i i)))))
291		 (setq *det* (cons *det* 1))))
292  (do ((i 1 (1+ i))
293       (elm))
294      ((> i m))
295    (setq elm (aref ax i i))
296    (setf (aref ax i (+ m i))
297	   (cond ($detout (cons (ptimes (cdr elm)
298					(pquotient (car *det*) (car elm))) 1))
299		 (t (ratinvert elm))))))
300
301(defun invert1 (k)
302  (prog (l r g i m n)
303     (setq l (length k) i 1)
304     (cond ((= l (length (car k))) nil)
305	   (t(merror (intl:gettext "invert: matrix must be square; found ~M rows, ~M columns.") l (length (car k)))))
306     loop (cond ((null k) (go l1)))
307     (setq r (car k))
308     (setq g (nconc g (list (nconc r (onen i l '(1 . 1) '(0 . 1))))))
309     (setq k (cdr k) i (1+ i))
310     (go loop)
311     l1   (setq k g)
312     (mtoa '*mat* (setq m (length k)) (setq n (length (car k))) k)
313     (setq k nil)
314     (cond ((diagp '*mat* m) (diaginv '*mat* m)) (t (tfgeli0 '*mat* m n)))
315     (setq k (atomat '*mat* m n (1+ m)))
316     (return k)))
317
318(defun diagp (ax m)
319  (declare (fixnum m))
320  (prog ((i 0) (j 0))
321     (declare (fixnum i j))
322     (setq ax (get-array-pointer ax))
323     loop1 (setq i (1+ i) j 0)
324     (cond ((> i m) (return t)))
325     loop2 (incf j)
326     (cond ((> j m) (go loop1))
327	   ((and (not (= i j)) (equal (aref ax i j) '(0 . 1))) nil)
328	   ((and (= i j) (not (equal (aref ax i j) '(0 . 1)))) nil)
329	   (t (return nil)))
330     (go loop2)))
331
332(defun tfgeli0 (x m n)
333  (cond ((or $sparse *det*) (tfgeli x m n))
334	(t (tfgeli x m n) (diaglize1 x m n))))
335
336;;  TWO-STEP FRACTION-FREE GAUSSIAN ELIMINATION ROUTINE
337
338(defun ritediv (x m n a)
339  (declare (fixnum m n))
340  (setq x (get-array-pointer x))
341  (prog ((j 0) (i 0) d)
342     (declare (fixnum i j))
343     (setq i m)
344     loop1 (when (zerop i) (return nil))
345     (setf (aref x i i) nil)
346     (setq j m)
347     loop (cond ((= j n) (decf i) (go loop1)))
348     (incf j)
349     (cond ((equal a 1)
350	    (setf (aref x i j) (cons (aref x i j) 1))
351	    (go loop)))
352     (setq d (ignore-rat-err (pquotient (aref x i j) a)))
353     (setq d (cond (d (cons d 1))
354		   (t (ratreduce (aref x i j) a))))
355     (setf (aref x i j) d)
356     (go loop)))
357
358(defun diaglize1 (x m n)
359  (setq x (get-array-pointer x))
360  (prog nil
361     (cond (*det* (return (ptimes *det* (aref x m m)))))
362     (setq *det* (cons (aref x m m) 1))
363     (cond ((not $detout) (return (ritediv x m n (aref x m m))))
364	   (t (return (ritediv x m n 1))))))
365
366;; Takes an M by N matrix and creates an array containing the elements
367;; of the matrix.  The array is associated "functionally" with the
368;; symbol NAME.
369;; For CL we have put it in the value cell-WFS.  Things still work.
370
371(defun mtoa (name m n mat)
372  (declare (fixnum m n))
373  (proclaim (list 'special name))
374  (setf (symbol-value name) (make-array (list (1+ m) (1+ n))))
375  (setq name (get-array-pointer name))
376  (do ((i 1 (1+ i))
377       (mat mat (cdr mat)))
378      ((> i m) nil)
379    (declare (fixnum i))
380    (do ((j 1 (1+ j))
381	 (row (car mat) (cdr row)))
382	((> j n))
383      (declare (fixnum j))
384      (setf (aref name i j) (car row)))))
385
386
387(defmfun $echelon (x)
388  (let (($ratmx t))
389    (newvarmat1 (setq x (check x))))
390  (let ((*ech* t) ($algebraic $algebraic))
391    (and (not $algebraic) (some #'algp varlist) (setq $algebraic t))
392    (setq x (cons '($matrix) (mxc (disreplist1 (echelon1 (replist1 (mcx (cdr x)))))))))
393  (if $ratmx x ($totaldisrep x)))
394
395(defun echelon1 (x)
396  (let ((m (length x))
397	(n (length (car x))))
398    (mtoa '*mat* m n x)
399    (setq x (catch 'rank (tfgeli '*mat* m n)))
400    (cond ((and *rank* x)
401	   (throw 'rnk x))
402	  (t (echelon2 '*mat* m n)))))
403
404(defun echelon2 (name m n)
405  (declare (fixnum m n))
406  (setq name (symbol-value name))
407  (prog ((j 0) row mat a)
408     (declare (fixnum j))
409     (incf m)
410     loop1 (when (= m 1) (return mat))
411     (setq m (1- m) j 0 a nil)
412     loop2 (when (= j n)
413	     (setq mat (cons row mat) row nil)
414	     (go loop1))
415     (incf j)
416     (setq row (nconc row (ncons (cond ((or (> m j) (equal (aref name m j) 0))
417					'(0 . 1))
418				       (a (ratreduce (aref name m j)a))
419				       (t (setq a (aref name m j)) '(1 . 1))))))
420     (go loop2)))
421
422(defun triang (x)
423  (let ((m (length x))
424	(n (length (car x)))
425	(*tri* t))
426    (mtoa '*mat* m n x)
427    (tfgeli '*mat* m n)
428    (triang2 '*mat* m n)))
429
430(defun triang2 (nam m n)
431  (declare (fixnum m n))
432  (setq nam (get-array-pointer nam))
433  (prog ((j 0) row mat)
434     (declare (fixnum j))
435     (setf (aref nam 0 0) 1)
436     (incf m)
437     loop1 (when (= m 1) (return mat))
438     (decf m)
439     (setq j 0)
440     loop2 (when (= j n)
441	     (setq mat (cons row mat) row nil)
442	     (go loop1))
443     (incf j)
444     (setq row (nconc row (ncons (if (> m j) '(0 . 1) (cons (aref nam m j) 1)))))
445     (go loop2)))
446
447(defun onen (n i var fill)
448  (prog (g)
449   loop (cond ((= i n) (setq g (cons var g)))
450	      ((zerop i) (return g))
451	      (t (setq g (cons fill g))))
452   (decf i)
453   (go loop)))
454
455(defun timex0 (x y)
456  (let ((u (check1 x))
457	(v (check1 y)))
458    (cond ((and (null u) (null v)) (list '(mtimes) x y))
459	  ((null u) (timex1 x (cons '($matrix) (mcx (cdr v)))))
460	  ((null v) (timex1 y (cons '($matrix) (mcx (cdr u)))))
461	  (t (cons '($matrix mult) (mxc (multiplymatrices (mcx (cdr u)) (mcx (cdr v)))))))))
462
463(defun timex1 (x y)
464  (setq y (check y))
465  (cond ((not $ratmx) (setq y (cdr y)))
466	(t (setq x (cdr (ratf x)) y (replist1 (cdr y)))))
467  (ctimesx x y))
468
469(defun ctimesx (x y)
470  (prog (c)
471   loop (cond ((null y)
472	       (return (cons '($matrix mult)
473			     (mxc (cond ((not $ratmx) c) (t (disreplist1 c))))))))
474   (setq c (nconc c (list (timesrow x (car y)))) y (cdr y))
475   (go loop)))
476
477(defun multiplymatrices (x y)
478  (cond ((and (null (cdr y)) (null (cdr x)))
479	 (and (cdar x) (setq y (transpose y))))
480	((and (null (cdar x)) (null (cdar y)))
481	 (and (cdr y) (setq x (transpose x)))))
482  (cond ((not (= (length (car x)) (length y)))
483	 (cond ((and (null (cdr y)) (= (length (car x)) (length (car y))))
484		(setq y (transpose y)))
485	       (t (merror (intl:gettext "MULTIPLYMATRICES: attempt to multiply nonconformable matrices."))))))
486  (cond ((not $ratmx) (multmat x y))
487	(t (setq x (replist1 x) y (replist1 y))
488	   (disreplist1 (multmat x y)))))
489
490(defun multmat (x y)
491  (prog (mat row yt rowx)
492     (setq yt (transpose y))
493     loop1 (when (null x) (return mat))
494     (setq rowx (car x) y yt)
495     loop2 (when (null y)
496	     (setq mat (nconc mat (ncons row)) x (cdr x) row nil)
497	     (go loop1))
498     (setq row (nconc row (ncons (multl rowx (car y)))) y (cdr y))
499     (go loop2)))
500
501;;; This actually takes the inner product of the two vectors.
502;;; I check for the most common cases for speed. "*" is a slight
503;;; violation of data abstraction here. The parser should turn "*" into
504;;; MTIMES, well, it may someday, which will break this code. Don't
505;;; hold your breath.
506
507(defun multl (a b)
508  (cond ((equal $matrix_element_add "+")
509	 (do ((ans (if (not $ratmx) 0 '(0 . 1))
510		   (cond ((not $ratmx)
511			  (cond ((equal $matrix_element_mult "*")
512				 (add ans (mul (car a) (car b))))
513				((equal $matrix_element_mult ".")
514				 (add ans (ncmul (car a) (car b))))
515				(t
516				 (add ans
517				      (meval `((,(getopr $matrix_element_mult))
518					       ((mquote simp) ,(car a))
519					       ((mquote simp) ,(car b))))))))
520			 (t
521			  (ratplus ans (rattimes (car a) (car b) t)))))
522	      (a a (cdr a))
523	      (b b (cdr b)))
524	     ((null a) ans)))
525	(t
526	 (mapply (getopr $matrix_element_add)
527		 (mapcar #'(lambda (u v)
528			     (meval `((,(getopr $matrix_element_mult))
529				      ((mquote simp) ,u)
530				      ((mquote simp) ,v))))
531			 a b)
532		 (getopr $matrix_element_add)))))
533
534(defun bbsort (l fn)
535  (nreverse (sort (copy-list l) fn)))
536
537(defun powerx (mat x)
538  (prog (n y)
539     (cond ((not (fixnump x))
540	    (return (list '(mncexpt simp) mat x)))
541	   ((= x 1) (return mat))
542	   ((minusp x)
543	    (setq x (- x) mat ($invert mat))
544	    (cond ($detout
545		   (return (let ((*inv* '$detout))
546			     (mul2* (power* (cadr mat) x)
547				    (fmapl1 #'(lambda (x) x)
548					    (powerx (caddr mat) x)))))))))
549     (newvarmat1 (setq mat (check mat)))
550     (setq n 1 mat (mcx (cdr mat)) y mat)
551     loop (if (= n x)
552	      (let (($scalarmatrixp (if (eq $scalarmatrixp '$all) '$all)))
553		(return (simplify (cons '($matrix mult) (mxc y))))))
554     (setq y (multiplymatrices y mat) n (1+ n))
555     (go loop)))
556
557;; The following $ALGEBRAIC code is so that
558;; RANK(MATRIX([1-SQRT(5),2],[-2,1+SQRT(5)])); will give 1.
559;; - JPG and BMT
560
561(defmfun $rank (x)
562  (let ((*rank* t) ($ratmx t) ($algebraic $algebraic))
563    (newvarmat1 (setq x (check x)))
564    (and (not $algebraic) (some #'algp varlist) (setq $algebraic t))
565    (setq x (replist1 (mcx (cdr x))))
566    (mtoa '*mat* (length x) (length (car x)) x)
567    (tfgeli '*mat* (length x) (length (car x)))))
568
569(defun replacerow (i y x)
570  (if (= i 1)
571      (cons y (cdr x))
572      (cons (car x) (replacerow (1- i) y (cdr x)))))
573
574(defun timesrow (y row)
575  (prog (ans)
576     (when (and $ratmx (atom y) y)
577       (setq y (cdr (ratf y))))
578     loop (when (null row) (return ans))
579     (setq ans (nconc ans (list (if (not $ratmx)
580				    (simptimes (list '(mtimes) y (car row)) 1 nil)
581				    (rattimes y (car row) t)))))
582     (setq row (cdr row))
583     (go loop)))
584
585(defmfun $triangularize (x)
586  (let (($ratmx t))
587    (newvarmat1 (setq x (check x))))
588  (let (($algebraic $algebraic))
589    (and (not $algebraic) (some #'algp varlist) (setq $algebraic t))
590    (setq x (cons '($matrix) (mxc (disreplist1 (triang (replist1 (mcx (cdr x)))))))))
591  (if $ratmx x ($totaldisrep x)))
592
593(defmfun $col (mat n)
594  (cons '($matrix) (mxc (transpose (list (nthcol (mcx (cdr (check mat))) n))))))
595
596(defun deletecol (n x)
597  (prog (m g)
598     (setq m x)
599     loop (when (null m) (return g))
600     (setq g (nconc g (ncons (deleterow n (car m)))) m (cdr m))
601     (go loop)))
602
603(defun deleterow (i m)
604  (cond ((or (null m) (< i 0)) (merror (intl:gettext "DELETEROW: matrix is null, or index is negative.")))
605	((= i 1) (cdr m))
606	(t (cons (car m) (deleterow (1- i) (cdr m))))))
607
608(defmfun $minor (mat m n)
609  (cons '($matrix) (mxc (minor m n (mcx (cdr (check mat)))))))
610
611(defun minor (i j m)
612  (deletecol j (deleterow i m)))
613
614(defmfun $row (mat m)
615  (cons '($matrix) (mxc (list (ith (mcx (cdr (check mat))) m)))))
616
617(defmfun $setelmx (elm m n mat)
618  (cond
619    ((not (and (integerp m) (integerp n)))
620	 (merror (intl:gettext "setelmx: indices must be integers; found: ~M, ~M") m n))
621    ((not ($matrixp mat))
622     (merror (intl:gettext "setelmx: last argument must be a matrix; found: ~M") mat))
623	((not (and (> m 0) (> n 0) (> (length mat) m) (> (length (cadr mat)) n)))
624	 (merror (intl:gettext "setelmx: no such element [~M, ~M]") m n)))
625  (rplaca (nthcdr n (car (nthcdr m mat))) elm) mat)
626
627;;; Here the function transpose can actually do simplification of
628;;; its argument. TRANSPOSE(TRANSPOSE(FOO)) => FOO.
629;;; If you think this is a hack, well, realize that the hack is
630;;; actually the fact that TRANSPOSE can return a noun form.
631
632
633
634(defmfun $transpose (mat)
635  (cond ((not (mxorlistp mat))
636	 (cond ((and (not (atom mat)) (member (mop mat) '($transpose %transpose) :test #'eq))
637		(cadr mat))
638	       (($scalarp mat) mat)
639	       ((mplusp mat)
640		`((mplus) .,(mapcar #'$transpose (cdr mat))))
641	       ((mtimesp mat)
642		`((mtimes) .,(mapcar #'$transpose (cdr mat))))
643	       ((mnctimesp mat)
644		`((mnctimes) .,(nreverse (mapcar #'$transpose (cdr mat)))))
645	       ((mncexptp mat)
646		(destructuring-let (((mat pow) (cdr mat)))
647		  `((mncexpt) ,($transpose mat) ,pow)))
648	       (t ($nounify '$transpose) (list '(%transpose) mat))))
649	(t
650	 (let ((ans (transpose (mcx (cdr (check mat))))))
651	   (cond ($matrix_element_transpose
652		  (setq ans (mapcar #'(lambda (u) (mapcar #'transpose-els u))
653				    ans))))
654	   `(($matrix) . ,(mxc ans))))))
655
656;;; THIS IS FOR TRANSPOSING THE ELEMENTS OF A MATRIX
657;;; A hack for Block matricies and tensors.
658
659(defun transpose-els (elem)
660  (cond ((eq $matrix_element_transpose '$transpose)
661	 ($transpose elem))
662	((eq $matrix_element_transpose '$nonscalars)
663	 (if ($nonscalarp elem)
664	     ($transpose elem)
665	     elem))
666	(t
667	 (meval `((,(getopr $matrix_element_transpose)) ((mquote simp) ,elem))))))
668
669
670(defmfun $submatrix (&rest x)
671  (prog (r c)
672     l1   (when (numberp (car x))
673	    (push (car x) r)
674	    (setq x (cdr x))
675	    (go l1))
676     (setq c (nreverse (bbsort (cdr x) #'>))
677	   r (nreverse (bbsort r #'>)))
678     (setq x (mcx (cdar x)))
679     l2   (if (null r)
680	      (go b)
681	      (setq x (deleterow (car r) x)))
682     (setq r (cdr r))
683     (go l2)
684     b    (when (null c) (return (cons '($matrix) (mxc x))))
685     (setq x (deletecol (car c) x) c (cdr c))
686     (go b)))
687
688
689(defmfun $list_matrix_entries (m)
690  (unless ($matrixp m)
691    (merror (intl:gettext "list_matrix_entries: argument must be a matrix; found: ~M") m))
692  (cons (if (null (cdr m)) '(mlist) (caadr m))
693	(loop for row in (cdr m) append (cdr row))))
694