code.delx.au - gnu-emacs/blob - lisp/calc/calc-mtx.el

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