lisp/calc/calc-mtx.el

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