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 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, or (at your option)
  14 ;; 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; see the file COPYING.  If not, write to the
  23 ;; Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
  24 ;; Boston, MA 02110-1301, USA.
  25
  26 ;;; Commentary:
  27
  28 ;;; Code:
  29
  30 ;; This file is autoloaded from calc-ext.el.
  31
  32 (require 'calc-ext)
  33 (require 'calc-macs)
  34
  35 (defun calc-mdet (arg)
  36   (interactive "P")
  37   (calc-slow-wrapper
  38    (calc-unary-op "mdet" 'calcFunc-det arg)))
  39
  40 (defun calc-mtrace (arg)
  41   (interactive "P")
  42   (calc-slow-wrapper
  43    (calc-unary-op "mtr" 'calcFunc-tr arg)))
  44
  45 (defun calc-mlud (arg)
  46   (interactive "P")
  47   (calc-slow-wrapper
  48    (calc-unary-op "mlud" 'calcFunc-lud arg)))
  49
  50
  51 ;;; Coerce row vector A to be a matrix.  [V V]
  52 (defun math-row-matrix (a)
  53   (if (and (Math-vectorp a)
  54            (not (math-matrixp a)))
  55       (list 'vec a)
  56     a))
  57
  58 ;;; Coerce column vector A to be a matrix.  [V V]
  59 (defun math-col-matrix (a)
  60   (if (and (Math-vectorp a)
  61            (not (math-matrixp a)))
  62       (cons 'vec (mapcar (function (lambda (x) (list 'vec x))) (cdr a)))
  63     a))
  64
  65
  66
  67 ;;; Multiply matrices A and B.  [V V V]
  68 (defun math-mul-mats (a b)
  69   (let ((mat nil)
  70         (cols (length (nth 1 b)))
  71         row col ap bp accum)
  72     (while (setq a (cdr a))
  73       (setq col cols
  74             row nil)
  75       (while (> (setq col (1- col)) 0)
  76         (setq ap (cdr (car a))
  77               bp (cdr b)
  78               accum (math-mul (car ap) (nth col (car bp))))
  79         (while (setq ap (cdr ap) bp (cdr bp))
  80           (setq accum (math-add accum (math-mul (car ap) (nth col (car bp))))))
  81         (setq row (cons accum row)))
  82       (setq mat (cons (cons 'vec row) mat)))
  83     (cons 'vec (nreverse mat))))
  84
  85 (defun math-mul-mat-vec (a b)
  86   (cons 'vec (mapcar (function (lambda (row)
  87                                  (math-dot-product row b)))
  88                      (cdr a))))
  89
  90
  91
  92 (defun calcFunc-tr (mat)   ; [Public]
  93   (if (math-square-matrixp mat)
  94       (math-matrix-trace-step 2 (1- (length mat)) mat (nth 1 (nth 1 mat)))
  95     (math-reject-arg mat 'square-matrixp)))
  96
  97 (defun math-matrix-trace-step (n size mat sum)
  98   (if (<= n size)
  99       (math-matrix-trace-step (1+ n) size mat
 100                               (math-add sum (nth n (nth n mat))))
 101     sum))
 102
 103
 104 ;;; Matrix inverse and determinant.
 105 (defun math-matrix-inv-raw (m)
 106   (let ((n (1- (length m))))
 107     (if (<= n 3)
 108         (let ((det (math-det-raw m)))
 109           (and (not (math-zerop det))
 110                (math-div
 111                 (cond ((= n 1) 1)
 112                       ((= n 2)
 113                        (list 'vec
 114                              (list 'vec
 115                                    (nth 2 (nth 2 m))
 116                                    (math-neg (nth 2 (nth 1 m))))
 117                              (list 'vec
 118                                    (math-neg (nth 1 (nth 2 m)))
 119                                    (nth 1 (nth 1 m)))))
 120                       ((= n 3)
 121                        (list 'vec
 122                              (list 'vec
 123                                    (math-sub (math-mul (nth 3 (nth 3 m))
 124                                                        (nth 2 (nth 2 m)))
 125                                              (math-mul (nth 3 (nth 2 m))
 126                                                        (nth 2 (nth 3 m))))
 127                                    (math-sub (math-mul (nth 3 (nth 1 m))
 128                                                        (nth 2 (nth 3 m)))
 129                                              (math-mul (nth 3 (nth 3 m))
 130                                                        (nth 2 (nth 1 m))))
 131                                    (math-sub (math-mul (nth 3 (nth 2 m))
 132                                                        (nth 2 (nth 1 m)))
 133                                              (math-mul (nth 3 (nth 1 m))
 134                                                        (nth 2 (nth 2 m)))))
 135                              (list 'vec
 136                                    (math-sub (math-mul (nth 3 (nth 2 m))
 137                                                        (nth 1 (nth 3 m)))
 138                                              (math-mul (nth 3 (nth 3 m))
 139                                                        (nth 1 (nth 2 m))))
 140                                    (math-sub (math-mul (nth 3 (nth 3 m))
 141                                                        (nth 1 (nth 1 m)))
 142                                              (math-mul (nth 3 (nth 1 m))
 143                                                        (nth 1 (nth 3 m))))
 144                                    (math-sub (math-mul (nth 3 (nth 1 m))
 145                                                        (nth 1 (nth 2 m)))
 146                                              (math-mul (nth 3 (nth 2 m))
 147                                                        (nth 1 (nth 1 m)))))
 148                              (list 'vec
 149                                    (math-sub (math-mul (nth 2 (nth 3 m))
 150                                                        (nth 1 (nth 2 m)))
 151                                              (math-mul (nth 2 (nth 2 m))
 152                                                        (nth 1 (nth 3 m))))
 153                                    (math-sub (math-mul (nth 2 (nth 1 m))
 154                                                        (nth 1 (nth 3 m)))
 155                                              (math-mul (nth 2 (nth 3 m))
 156                                                        (nth 1 (nth 1 m))))
 157                                    (math-sub (math-mul (nth 2 (nth 2 m))
 158                                                        (nth 1 (nth 1 m)))
 159                                              (math-mul (nth 2 (nth 1 m))
 160                                                        (nth 1 (nth 2 m))))))))
 161                 det)))
 162       (let ((lud (math-matrix-lud m)))
 163         (and lud
 164              (math-lud-solve lud (calcFunc-idn 1 n)))))))
 165
 166 (defun calcFunc-det (m)
 167   (if (math-square-matrixp m)
 168       (math-with-extra-prec 2 (math-det-raw m))
 169     (if (and (eq (car-safe m) 'calcFunc-idn)
 170              (or (math-zerop (nth 1 m))
 171                  (math-equal-int (nth 1 m) 1)))
 172         (nth 1 m)
 173       (math-reject-arg m 'square-matrixp))))
 174
 175 ;; The variable math-det-lu is local to math-det-raw, but is
 176 ;; used by math-det-step, which is called by math-det-raw.
 177 (defvar math-det-lu)
 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 ((math-det-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 math-det-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 (provide 'calc-mtx)
 369
 370 ;;; arch-tag: fc0947b1-90e1-4a23-8950-d8ead9c3a306
 371 ;;; calc-mtx.el ends here