lisp/calc/calc-mtx.el

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