1;;; calc-mtx.el --- matrix functions for Calc 2 3;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004, 4;; 2005, 2006, 2007 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 2, 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 372