src/matrix.lisp

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