lisp/calc/calc-mtx.el

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