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