src/matrix.lisp

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