cgit.sxemacs.org Git - packages/blob - xemacs-packages/calc/calc-mat.el

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