src/sprdet.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;     The data in this file contains enhancements.                   ;;;;;
   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
  18               *ptr* *ptc* *blk* ml* *detsign* rzl*))
  19
  20 (defun sprdet (ax n)
  21   (declare (fixnum n))
  22   (setq ax (get-array-pointer ax))
  23   (prog ((j 0) rodr codr bl det (dm 0) (r 0) (i 0))
  24      (declare (fixnum i j dm r))
  25      (setq det 1)
  26      (setq  *ptr* (make-array (1+ n)))
  27      (setq  *ptc* (make-array (1+ n)))
  28      (setq bl (tmlattice ax '*ptr* '*ptc* n)) ;tmlattice isn't defined anywhere -- are_muc
  29      (when (null bl) (return 0))
  30      (setq rodr (apply #'append bl))
  31      (setq codr (mapcar #'cadr rodr))
  32      (setq rodr (mapcar #'car rodr))
  33      (setq det (* (prmusign rodr) (prmusign codr)))
  34      (setq bl (mapcar #'length bl))
  35      loop1 (when (null bl) (return det))
  36      (setq i (car bl))
  37      (setq dm i)
  38      (setq *blk* (make-array (list (1+ dm) (1+ dm))))
  39      (cond ((= dm 1)
  40             (setq det (gptimes det (car (aref ax (aref *ptr* (1+ r)) (aref *ptc* (1+ r))))))
  41             (go next))
  42            ((= dm 2)
  43             (setq det (gptimes det
  44                                (gpdifference
  45                                 (gptimes (car (aref ax (aref *ptr* (1+ r)) (aref *ptc* (1+ r))))
  46                                          (car (aref ax (aref *ptr* (+ 2 r)) (aref *ptc* (+ 2 r)))))
  47                                 (gptimes (car (aref ax (aref *ptr* (1+ r)) (aref *ptc* (+ 2 r))))
  48                                          (car (aref ax (aref *ptr* (+ 2 r)) (aref *ptc* (1+ r))))))))
  49             (go next)))
  50      loop2 (when (= i 0) (go cmp))
  51      (setq j dm)
  52      loop3 (when (= j 0) (decf i) (go loop2))
  53      (setf (aref *blk* i j) (car (aref ax (aref *ptr* (+ r i)) (aref *ptc*(+ r j)))))
  54      (decf j) (go loop3)
  55      cmp
  56      (setq det (gptimes det (tdbu '*blk* dm)))
  57      next
  58      (incf r dm)
  59      (setq bl (cdr bl))
  60      (go loop1)))
  61
  62 (defun minorl (x n l nz)
  63   (declare (fixnum n))
  64   (prog (ans s rzl* (col 1) (n2 (truncate n 2)) d dl z a elm rule)
  65      (declare (fixnum n2  col ))
  66      (decf n2)
  67      (setq dl l l nil nz (cons nil nz))
  68      l1 (when (null nz) (return ans))
  69      l3 (setq z (car nz))
  70      (cond ((null l) (if dl (push dl ans) (return nil))
  71             (setq nz (cdr nz) col (1+ col) l dl dl nil)
  72             (go l1)))
  73      (setq a (caar l) )
  74      l2 (cond ((null z)
  75                (cond (rule (rplaca (car l) (list a rule))
  76                            (setq rule nil) (setq l (cdr l)))
  77                      ((null (cdr l))
  78                       (rplaca (car l) (list a 0))
  79                       (setq l (cdr l)))
  80                      (t (rplaca l (cadr l))
  81                         (rplacd l (cddr l))))
  82                (go l3)))
  83      (setq elm (car z) z (cdr z))
  84      (setq s (signnp elm a))
  85      (cond (s (setq d (delete elm (copy-list a) :test #'equal))
  86               (cond ((membercar d dl)
  87                      (go on))
  88                     (t
  89                      (when (or (< col n2) (not (singp x d col n)))
  90                        (push (cons d 1) dl)
  91                        (go on))))))
  92      (go l2)
  93      on (setq rule (cons (list d s elm (1- col)) rule))
  94      (go l2)))
  95
  96 (declare-top (special j))
  97
  98 (defun singp (x ml col n)
  99   (declare (fixnum col n))
 100   (prog (i (j col) l)
 101      (declare (fixnum  j))
 102      (setq l ml)
 103      (if (null ml)
 104          (go loop)
 105          (setq i (car ml)
 106                ml (cdr ml)))
 107      (cond ((member i rzl* :test #'equal) (return t))
 108            ((zrow x i col n) (return (setq rzl*(cons i rzl*)))))
 109      loop (cond ((> j n) (return nil))
 110                 ((every #'(lambda (i) (equal (aref x i j) 0)) l)
 111                  (return t)))
 112      (incf j)
 113      (go loop)))
 114
 115 (declare-top (unspecial j))
 116
 117 (defun tdbu (x n)
 118   (declare (fixnum n))
 119   (prog (a ml* nl nml dd)
 120      (setq *detsign* 1)
 121      (setq x (get-array-pointer x))
 122      (detpivot x n)
 123      (setq x (get-array-pointer 'x*))
 124      (setq nl (nzl x n))
 125      (when (member nil nl :test #'eq) (return 0))
 126      (setq a (minorl x n (list (cons (nreverse(index* n)) 1)) nl))
 127      (setq nl nil)
 128      (when (null a) (return 0))
 129      (tb2 x (car a) n)
 130      tag2
 131      (setq ml* (cons (cons nil nil) (car a)))
 132      (setq a (cdr a))
 133      (when (null a) (return (if (= *detsign* 1)
 134                                 (cadadr ml*)
 135                                 (gpctimes -1  (cadadr ml*)))))
 136      (setq nml (car a))
 137      tag1 (when (null nml) (go tag2))
 138      (setq dd  (car nml))
 139      (setq nml (cdr nml))
 140      (nbn dd x)
 141      (go tag1)))
 142
 143 (defun nbn (rule 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 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)))