1;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;; 2;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3;;; The data in this file contains enhancments. ;;;;; 4;;; ;;;;; 5;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;; 6;;; All rights reserved ;;;;; 7;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 8;;; (c) Copyright 1980 Massachusetts Institute of Technology ;;; 9;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 10 11(in-package :maxima) 12 13(macsyma-module sprdet) 14 15;; THIS IS THE NEW DETERMINANT PACKAGE 16 17(declare-top (special x *ptr* *ptc* *blk* $ratmx ml* *detsign* rzl*)) 18 19(defun sprdet (ax n) 20 (declare (fixnum n)) 21 (setq ax (get-array-pointer ax)) 22 (prog ((j 0) rodr codr bl det (dm 0) (r 0) (i 0)) 23 (declare (fixnum i j dm r)) 24 (setq det 1) 25 (setq *ptr* (make-array (1+ n))) 26 (setq *ptc* (make-array (1+ n))) 27 (setq bl (tmlattice ax '*ptr* '*ptc* n)) ;tmlattice isn't defined anywhere -- are_muc 28 (when (null bl) (return 0)) 29 (setq rodr (apply #'append bl)) 30 (setq codr (mapcar #'cadr rodr)) 31 (setq rodr (mapcar #'car rodr)) 32 (setq det (* (prmusign rodr) (prmusign codr))) 33 (setq bl (mapcar #'length bl)) 34 loop1 (when (null bl) (return det)) 35 (setq i (car bl)) 36 (setq dm i) 37 (setq *blk* (make-array (list (1+ dm) (1+ dm)))) 38 (cond ((= dm 1) 39 (setq det (gptimes det (car (aref ax (aref *ptr* (1+ r)) (aref *ptc* (1+ r)))))) 40 (go next)) 41 ((= dm 2) 42 (setq det (gptimes det 43 (gpdifference 44 (gptimes (car (aref ax (aref *ptr* (1+ r)) (aref *ptc* (1+ r)))) 45 (car (aref ax (aref *ptr* (+ 2 r)) (aref *ptc* (+ 2 r))))) 46 (gptimes (car (aref ax (aref *ptr* (1+ r)) (aref *ptc* (+ 2 r)))) 47 (car (aref ax (aref *ptr* (+ 2 r)) (aref *ptc* (1+ r)))))))) 48 (go next))) 49 loop2 (when (= i 0) (go cmp)) 50 (setq j dm) 51 loop3 (when (= j 0) (decf i) (go loop2)) 52 (setf (aref *blk* i j) (car (aref ax (aref *ptr* (+ r i)) (aref *ptc*(+ r j))))) 53 (decf j) (go loop3) 54 cmp 55 (setq det (gptimes det (tdbu '*blk* dm))) 56 next 57 (incf r dm) 58 (setq bl (cdr bl)) 59 (go loop1))) 60 61(defun minorl (x n l nz) 62 (declare (fixnum n)) 63 (prog (ans s rzl* (col 1) (n2 (truncate n 2)) d dl z a elm rule) 64 (declare (fixnum n2 col )) 65 (decf n2) 66 (setq dl l l nil nz (cons nil nz)) 67 l1 (when (null nz) (return ans)) 68 l3 (setq z (car nz)) 69 (cond ((null l) (if dl (push dl ans) (return nil)) 70 (setq nz (cdr nz) col (1+ col) l dl dl nil) 71 (go l1))) 72 (setq a (caar l) ) 73 l2 (cond ((null z) 74 (cond (rule (rplaca (car l) (list a rule)) 75 (setq rule nil) (setq l (cdr l))) 76 ((null (cdr l)) 77 (rplaca (car l) (list a 0)) 78 (setq l (cdr l))) 79 (t (rplaca l (cadr l)) 80 (rplacd l (cddr l)))) 81 (go l3))) 82 (setq elm (car z) z (cdr z)) 83 (setq s (signnp elm a)) 84 (cond (s (setq d (delete elm (copy-list a) :test #'equal)) 85 (cond ((membercar d dl) 86 (go on)) 87 (t 88 (when (or (< col n2) (not (singp x d col n))) 89 (push (cons d 1) dl) 90 (go on)))))) 91 (go l2) 92 on (setq rule (cons (list d s elm (1- col)) rule)) 93 (go l2))) 94 95(declare-top (special j)) 96 97(defun singp (x ml col n) 98 (declare (fixnum col n)) 99 (prog (i (j col) l) 100 (declare (fixnum j)) 101 (setq l ml) 102 (if (null ml) 103 (go loop) 104 (setq i (car ml) 105 ml (cdr ml))) 106 (cond ((member i rzl* :test #'equal) (return t)) 107 ((zrow x i col n) (return (setq rzl*(cons i rzl*))))) 108 loop (cond ((> j n) (return nil)) 109 ((every #'(lambda (i) (equal (aref x i j) 0)) l) 110 (return t))) 111 (incf j) 112 (go loop))) 113 114(declare-top (unspecial j)) 115 116(defun tdbu (x n) 117 (declare (fixnum n)) 118 (prog (a ml* nl nml dd) 119 (setq *detsign* 1) 120 (setq x (get-array-pointer x)) 121 (detpivot x n) 122 (setq x (get-array-pointer 'x*)) 123 (setq nl (nzl x n)) 124 (when (member nil nl :test #'eq) (return 0)) 125 (setq a (minorl x n (list (cons (nreverse(index* n)) 1)) nl)) 126 (setq nl nil) 127 (when (null a) (return 0)) 128 (tb2 x (car a) n) 129 tag2 130 (setq ml* (cons (cons nil nil) (car a))) 131 (setq a (cdr a)) 132 (when (null a) (return (if (= *detsign* 1) 133 (cadadr ml*) 134 (gpctimes -1 (cadadr ml*))))) 135 (setq nml (car a)) 136 tag1 (when (null nml) (go tag2)) 137 (setq dd (car nml)) 138 (setq nml (cdr nml)) 139 (nbn dd) 140 (go tag1))) 141 142(defun nbn (rule) 143 (declare (special x)) 144 (prog (ans r a) 145 (setq ans 0 r (cadar rule)) 146 (when (equal r 0) (return 0)) 147 (rplaca rule (caar rule)) 148 loop (when (null r) (return (rplacd rule (cons ans (cdr rule))))) 149 (setq a (car r) r(cdr r)) 150 (setq ans (gpplus ans (gptimes (if (= (cadr a) 1) 151 (aref x (caddr a) (cadddr a)) 152 (gpctimes (cadr a) (aref x (caddr a) (cadddr a)))) 153 (getminor (car a))))) 154 (go loop))) 155 156(defun getminor (index) 157 (cond ((null (setq index (assoc index ml* :test #'equal))) 0) 158 (t (rplacd (cdr index) (1- (cddr index))) 159 (when (= (cddr index) 0) 160 (setf index (delete index ml* :test #'equal))) 161 (cadr index)))) 162 163(defun tb2 (x l n) 164 (declare (fixnum n )) 165 (prog ((n-1 (1- n)) b a) 166 (declare (fixnum n-1)) 167 loop (when (null l) (return nil)) 168 (setq a (car l) l (cdr l) b (car a)) 169 (rplacd a (cons (gpdifference (gptimes (aref x (car b) n-1) (aref x (cadr b) n)) 170 (gptimes (aref x (car b) n) (aref x (cadr b) n-1))) 171 (cdr a))) 172 (go loop))) 173 174(defun zrow (x i col n) 175 (declare (fixnum i col n)) 176 (prog ((j col)) 177 (declare (fixnum j)) 178 loop (cond ((> j n) 179 (return t)) 180 ((equal (aref x i j) 0) 181 (incf j) 182 (go loop))))) 183 184(defun nzl (a n) 185 (declare (fixnum n)) 186 (prog ((i 0) (j (- n 2)) d l) 187 (declare (fixnum i j)) 188 loop0 (when (= j 0) (return l)) 189 (setq i n) 190 loop1 (when (= i 0) 191 (push d l) 192 (setq d nil) 193 (decf j) 194 (go loop0)) 195 (when (not (equal (aref a i j) 0)) 196 (push i d)) 197 (decf i) 198 (go loop1))) 199 200(defun signnp (e l) 201 (prog (i) 202 (setq i 1) 203 loop (cond ((null l) (return nil)) 204 ((equal e (car l)) (return i))) 205 (setq l (cdr l) i (- i)) 206 (go loop))) 207 208(defun membercar (e l) 209 (prog() 210 loop (cond ((null l) 211 (return nil)) 212 ((equal e (caar l)) 213 (return (rplacd (car l) (1+ (cdar l)))))) 214 (setq l (cdr l)) 215 (go loop))) 216 217(declare-top (unspecial x ml* rzl*)) 218 219(defun atranspose (a n) 220 (prog (i j d) 221 (setq i 0) 222 loop1 (setq i (1+ i) j i) 223 (when (> i n) (return nil)) 224 loop2 (incf j) 225 (when (> j n) (go loop1)) 226 (setq d (aref a i j)) 227 (setf (aref a i j) (aref a j i)) 228 (setf (aref a j i) d) 229 (go loop2))) 230 231(defun mxcomp (l1 l2) 232 (prog() 233 loop (cond ((null l1) (return t)) 234 ((car> (car l1) (car l2)) (return t)) 235 ((car> (car l2) (car l1)) (return nil))) 236 (setq l1 (cdr l1) l2 (cdr l2)) 237 (go loop))) 238 239(defun prmusign (l) 240 (prog ((b 0) a d) 241 (declare (fixnum b)) 242 loop (when (null l) (return (if (even b) 1 -1))) 243 (setq a (car l) l (cdr l) d l) 244 loop1 (cond ((null d) (go loop)) 245 ((> a (car d)) (incf b))) 246 (setq d (cdr d)) 247 (go loop1))) 248 249(defun detpivot (x n) 250 (prog (r0 c0) 251 (setq c0 (colrow0 x n nil) 252 r0 (colrow0 x n t)) 253 (setq c0 (nreverse (bbsort c0 #'car>))) 254 (setq r0 (nreverse (bbsort r0 #'car>))) 255 (when (not (mxcomp c0 r0)) 256 (atranspose x n) 257 (setq c0 r0)) 258 (setq *detsign* (prmusign (mapcar #'car c0))) 259 (newmat 'x* x n c0))) 260 261(defun newmat(x y n l) 262 (prog (i j jl) 263 (setf (symbol-value x) (make-array (list (1+ n) (1+ n)))) 264 (setq x (get-array-pointer x)) 265 (setq j 0) 266 loop (setq i 0 j (1+ j)) 267 (when (null l) (return nil)) 268 (setq jl (cdar l) l (cdr l)) 269 tag (incf i) 270 (when (> i n) (go loop)) 271 (setf (aref x i j) (aref y i jl)) 272 (go tag))) 273 274(defun car> (a b) 275 (> (car a) (car b))) 276 277;; ind=t for row ortherwise col 278 279(defun colrow0 (a n ind) 280 (declare (fixnum n)) 281 (prog ((i 0) (j n) l (c 0)) 282 (declare (fixnum i c j)) 283 loop0 (cond ((= j 0) (return l))) 284 (setq i n) 285 loop1 (when (= i 0) 286 (push (cons c j) l) 287 (setq c 0) 288 (decf j) 289 (go loop0)) 290 (when (equal (if ind (aref a j i) (aref a i j)) 0) 291 (incf c)) 292 (decf i) 293 (go loop1))) 294 295(defun gpdifference (a b) 296 (if $ratmx 297 (pdifference a b) 298 (simplus(list '(mplus) a (list '(mtimes) -1 b)) 1 nil))) 299 300(defun gpctimes(a b) 301 (if $ratmx 302 (pctimes a b) 303 (simptimes(list '(mtimes) a b) 1 nil))) 304 305(defun gptimes(a b) 306 (if $ratmx 307 (ptimes a b) 308 (simptimes (list '(mtimes) a b) 1 nil))) 309 310(defun gpplus(a b) 311 (if $ratmx 312 (pplus a b) 313 (simplus(list '(mplus) a b) 1 nil))) 314