code.delx.au - gnu-emacs/blob - 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: Jay Belanger <belanger@truman.edu>
   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 ;; The variable math-det-lu is local to math-det-raw, but is
 179 ;; used by math-det-step, which is called by math-det-raw.
 180 (defvar math-det-lu)
 181
 182 (defun math-det-raw (m)
 183   (let ((n (1- (length m))))
 184     (cond ((= n 1)
 185            (nth 1 (nth 1 m)))
 186           ((= n 2)
 187            (math-sub (math-mul (nth 1 (nth 1 m))
 188                                (nth 2 (nth 2 m)))
 189                      (math-mul (nth 2 (nth 1 m))
 190                                (nth 1 (nth 2 m)))))
 191           ((= n 3)
 192            (math-sub
 193             (math-sub
 194              (math-sub
 195               (math-add
 196                (math-add
 197                 (math-mul (nth 1 (nth 1 m))
 198                           (math-mul (nth 2 (nth 2 m))
 199                                     (nth 3 (nth 3 m))))
 200                 (math-mul (nth 2 (nth 1 m))
 201                           (math-mul (nth 3 (nth 2 m))
 202                                     (nth 1 (nth 3 m)))))
 203                (math-mul (nth 3 (nth 1 m))
 204                          (math-mul (nth 1 (nth 2 m))
 205                                    (nth 2 (nth 3 m)))))
 206               (math-mul (nth 3 (nth 1 m))
 207                         (math-mul (nth 2 (nth 2 m))
 208                                   (nth 1 (nth 3 m)))))
 209              (math-mul (nth 1 (nth 1 m))
 210                        (math-mul (nth 3 (nth 2 m))
 211                                  (nth 2 (nth 3 m)))))
 212             (math-mul (nth 2 (nth 1 m))
 213                       (math-mul (nth 1 (nth 2 m))
 214                                 (nth 3 (nth 3 m))))))
 215           (t (let ((lud (math-matrix-lud m)))
 216                (if lud
 217                    (let ((math-det-lu (car lud)))
 218                      (math-det-step n (nth 2 lud)))
 219                  0))))))
 220
 221 (defun math-det-step (n prod)
 222   (if (> n 0)
 223       (math-det-step (1- n) (math-mul prod (nth n (nth n math-det-lu))))
 224     prod))
 225
 226 ;;; This returns a list (LU index d), or nil if not possible.
 227 ;;; Argument M must be a square matrix.
 228 (defvar math-lud-cache nil)
 229 (defun math-matrix-lud (m)
 230   (let ((old (assoc m math-lud-cache))
 231         (context (list calc-internal-prec calc-prefer-frac)))
 232     (if (and old (equal (nth 1 old) context))
 233         (cdr (cdr old))
 234       (let* ((lud (catch 'singular (math-do-matrix-lud m)))
 235              (entry (cons context lud)))
 236         (if old
 237             (setcdr old entry)
 238           (setq math-lud-cache (cons (cons m entry) math-lud-cache)))
 239         lud))))
 240
 241 ;;; Numerical Recipes section 2.3; implicit pivoting omitted.
 242 (defun math-do-matrix-lud (m)
 243   (let* ((lu (math-copy-matrix m))
 244          (n (1- (length lu)))
 245          i (j 1) k imax sum big
 246          (d 1) (index nil))
 247     (while (<= j n)
 248       (setq i 1
 249             big 0
 250             imax j)
 251       (while (< i j)
 252         (math-working "LUD step" (format "%d/%d" j i))
 253         (setq sum (nth j (nth i lu))
 254               k 1)
 255         (while (< k i)
 256           (setq sum (math-sub sum (math-mul (nth k (nth i lu))
 257                                             (nth j (nth k lu))))
 258                 k (1+ k)))
 259         (setcar (nthcdr j (nth i lu)) sum)
 260         (setq i (1+ i)))
 261       (while (<= i n)
 262         (math-working "LUD step" (format "%d/%d" j i))
 263         (setq sum (nth j (nth i lu))
 264               k 1)
 265         (while (< k j)
 266           (setq sum (math-sub sum (math-mul (nth k (nth i lu))
 267                                             (nth j (nth k lu))))
 268                 k (1+ k)))
 269         (setcar (nthcdr j (nth i lu)) sum)
 270         (let ((dum (math-abs-approx sum)))
 271           (if (Math-lessp big dum)
 272               (setq big dum
 273                     imax i)))
 274         (setq i (1+ i)))
 275       (if (> imax j)
 276           (setq lu (math-swap-rows lu j imax)
 277                 d (- d)))
 278       (setq index (cons imax index))
 279       (let ((pivot (nth j (nth j lu))))
 280         (if (math-zerop pivot)
 281             (throw 'singular nil)
 282           (setq i j)
 283           (while (<= (setq i (1+ i)) n)
 284             (setcar (nthcdr j (nth i lu))
 285                     (math-div (nth j (nth i lu)) pivot)))))
 286       (setq j (1+ j)))
 287     (list lu (nreverse index) d)))
 288
 289 (defun math-swap-rows (m r1 r2)
 290   (or (= r1 r2)
 291       (let* ((r1prev (nthcdr (1- r1) m))
 292              (row1 (cdr r1prev))
 293              (r2prev (nthcdr (1- r2) m))
 294              (row2 (cdr r2prev))
 295              (r2next (cdr row2)))
 296         (setcdr r2prev row1)
 297         (setcdr r1prev row2)
 298         (setcdr row2 (cdr row1))
 299         (setcdr row1 r2next)))
 300   m)
 301
 302
 303 (defun math-lud-solve (lud b &optional need)
 304   (if lud
 305       (let* ((x (math-copy-matrix b))
 306              (n (1- (length x)))
 307              (m (1- (length (nth 1 x))))
 308              (lu (car lud))
 309              (col 1)
 310              i j ip ii index sum)
 311         (while (<= col m)
 312           (math-working "LUD solver step" col)
 313           (setq i 1
 314                 ii nil
 315                 index (nth 1 lud))
 316           (while (<= i n)
 317             (setq ip (car index)
 318                   index (cdr index)
 319                   sum (nth col (nth ip x)))
 320             (setcar (nthcdr col (nth ip x)) (nth col (nth i x)))
 321             (if (null ii)
 322                 (or (math-zerop sum)
 323                     (setq ii i))
 324               (setq j ii)
 325               (while (< j i)
 326                 (setq sum (math-sub sum (math-mul (nth j (nth i lu))
 327                                                   (nth col (nth j x))))
 328                       j (1+ j))))
 329             (setcar (nthcdr col (nth i x)) sum)
 330             (setq i (1+ i)))
 331           (while (>= (setq i (1- i)) 1)
 332             (setq sum (nth col (nth i x))
 333                   j i)
 334             (while (<= (setq j (1+ j)) n)
 335               (setq sum (math-sub sum (math-mul (nth j (nth i lu))
 336                                                 (nth col (nth j x))))))
 337             (setcar (nthcdr col (nth i x))
 338                     (math-div sum (nth i (nth i lu)))))
 339           (setq col (1+ col)))
 340         x)
 341     (and need
 342          (math-reject-arg need "*Singular matrix"))))
 343
 344 (defun calcFunc-lud (m)
 345   (if (math-square-matrixp m)
 346       (or (math-with-extra-prec 2
 347             (let ((lud (math-matrix-lud m)))
 348               (and lud
 349                    (let* ((lmat (math-copy-matrix (car lud)))
 350                           (umat (math-copy-matrix (car lud)))
 351                           (n (1- (length (car lud))))
 352                           (perm (calcFunc-idn 1 n))
 353                           i (j 1))
 354                      (while (<= j n)
 355                        (setq i 1)
 356                        (while (< i j)
 357                          (setcar (nthcdr j (nth i lmat)) 0)
 358                          (setq i (1+ i)))
 359                        (setcar (nthcdr j (nth j lmat)) 1)
 360                        (while (<= (setq i (1+ i)) n)
 361                          (setcar (nthcdr j (nth i umat)) 0))
 362                        (setq j (1+ j)))
 363                      (while (>= (setq j (1- j)) 1)
 364                        (let ((pos (nth (1- j) (nth 1 lud))))
 365                          (or (= pos j)
 366                              (setq perm (math-swap-rows perm j pos)))))
 367                      (list 'vec perm lmat umat)))))
 368           (math-reject-arg m "*Singular matrix"))
 369     (math-reject-arg m 'square-matrixp)))
 370
 371 ;;; arch-tag: fc0947b1-90e1-4a23-8950-d8ead9c3a306
 372 ;;; calc-mtx.el ends here