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