;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;     The data in this file contains enhancments.                    ;;;;;
;;;                                                                    ;;;;;
;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
;;;     All rights reserved                                            ;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;     (c) Copyright 1980 Massachusetts Institute of Technology         ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(in-package "MAXIMA")
(macsyma-module sprdet)

;; THIS IS THE NEW DETERMINANT PACKAGE

(DECLARE-TOP(SPECIAL X *ptr* *ptc* *blk* $SPARSE $RATMX ML* *DETSIGN* RZL*) (GENPREFIX ND))

(DEFUN SPRDET(AX N)
  (DECLARE(FIXNUM N))
  (setq ax (get-array-pointer ax))
  (PROG ((J 0) RODR CODR BL DET (DM 0) (R 0) (I 0))
	(DECLARE(FIXNUM I J DM R))
       (SETQ DET 1.)
       (setq  *PTR* (*ARRAY nil T  (f1+ N)))
       (setq  *PTC* (*ARRAY nil T  (f1+ N)))
       (SETQ BL (TMLATTICE AX '*PTR* '*PTC* N))
       (COND ((NULL BL)(RETURN 0)))
       (SETQ RODR(APPLY(FUNCTION APPEND) BL))
       (SETQ CODR(MAPCAR (FUNCTION CADR) RODR))
       (SETQ RODR(MAPCAR (FUNCTION CAR) RODR))
       (SETQ DET(*(PRMUSIGN RODR)(PRMUSIGN CODR)))
       (SETQ BL (MAPCAR (FUNCTION LENGTH) BL ))
 LOOP1 (COND ((NULL BL) (RETURN DET)))
       (SETQ I (CAR BL) )(SETQ DM I)
       (setq  *BLK* (*ARRAY nil T (f1+ DM)(f1+ DM)))
       (COND ((= DM 1.)
	      (SETQ DET (GPTIMES DET (CAR(AREF AX (AREF *PTR* (f1+ R))(AREF *PTC*(f1+ R))))))
	      (GO NEXT))
	    ((= DM 2.)
	     (SETQ DET (GPTIMES DET
				(GPDIFFERENCE
				 (GPTIMES (CAR (AREF AX (AREF *PTR* (f1+ R))(AREF *PTC* (f1+ R))))
					  (CAR (AREF AX (AREF *PTR* (f+ 2. R))(AREF *PTC* (f+ 2. R)))))
				 (GPTIMES (CAR (AREF AX (AREF *PTR* (f1+ R))(AREF *PTC* (f+ 2. R))))
					  (CAR (AREF AX (AREF *PTR* (f+ 2. R))(AREF *PTC* (f1+  R))))))))
	     (GO NEXT)))
       LOOP2 (COND ((= I 0)(GO CMP)))
       (SETQ J DM)
       LOOP3 (COND ((= J 0) (SETQ I (f1- I)) (GO LOOP2)))
       (STORE (aref *BLK* I J)(CAR  (AREF AX  (AREF *PTR* (f+ R I))  (AREF *PTC*(f+ R J)))))
       (SETQ J (f1- J)) (GO LOOP3)
       CMP
       (SETQ DET (GPTIMES DET (TDBU '*BLK* DM)))
       NEXT   
       (SETQ R(f+ R DM))
       (SETQ BL (CDR BL))
       (GO LOOP1)
       ))

(DEFUN MINORL (X N L NZ)
  (DECLARE(FIXNUM  N ))
  (PROG (ANS S RZL* (COL 1) ( N2 (// N 2.))  D DL Z A ELM RULE)
	  (DECLARE(FIXNUM N2  COL ))
       (SETQ N2(f1- N2))
       (SETQ DL L L NIL NZ (CONS NIL NZ))
       L1(COND((NULL NZ)(RETURN ANS)))
       L3(SETQ Z (CAR NZ))
       (COND ((NULL L) (COND (DL (SETQ ANS (CONS DL ANS)))
			     (T (RETURN NIL)))
		       (SETQ NZ (CDR NZ) COL (f1+ COL) L DL DL NIL)
		       (GO L1)))
       (SETQ A (CAAR L) )
       L2(COND((NULL Z)
	       (COND (RULE (RPLACA (CAR L) (LIST A RULE))
			   (SETQ RULE NIL) (SETQ L (CDR L)))
		     ((NULL (CDR L))
		      (RPLACA (CAR L) (LIST A 0))
		      (SETQ L (CDR L)))
		     (T (RPLACA L (CADR L))
			(RPLACD L (CDDR L))))
	       (GO L3)))
       (SETQ ELM (CAR Z) Z (CDR Z))
       (SETQ S(SIGNNP ELM A))
       (COND(S(SETQ D (zl-DELETE ELM (COPY1 A)))
	      (COND((MEMBERCAR D DL) (GO ON))
		   (T
		    (COND((OR(< COL N2)(NOT(SINGP X D COL N)))(SETQ DL (CONS (CONS D 1) DL))(GO ON)))
		    ))))
       (GO L2)
       ON(SETQ RULE(CONS (LIST D S ELM (f1- COL)) RULE))
       (GO L2)))

#-NIL
(DECLARE-TOP(SPECIAL J))

(DEFUN SINGP (X ML COL N)
    #+cl (DECLARE (FIXNUM COL N))
  ;#-Multics (DECLARE (FIXNUM COL N I J))
  (PROG (I (J col) L) 
	(DECLARE (FIXNUM  J))
	(SETQ L ML)
	(COND((NULL ML)(GO LOOP))
	     (T (SETQ I (CAR ML) ML (CDR ML))))
	(COND((zl-MEMBER I RZL*)(RETURN T))
	     ((ZROW X I COL N)(RETURN (SETQ RZL*(CONS I RZL*)))))
	LOOP(COND((> J N)(RETURN NIL))
		 ((EVERY #'(LAMBDA (I) (EQUAL (AREF X I J) 0)) L)
		  (RETURN T)))
	(SETQ J(f1+ J))(GO LOOP)
	))
#-NIL
(DECLARE-TOP(UNSPECIAL J))

(DEFUN TDBU (X N)
  (DECLARE(FIXNUM N))
  (PROG(A ML* NL NML DD)
       (SETQ *DETSIGN* 1)
      (setq x ( get-array-pointer x))
       (DETPIVOT X N)
       (SETQ X (get-array-pointer 'X*))
;       (setq x ( get-array-pointer x))
       (SETQ NL (NZL X N))
       (COND ((MEMQ NIL NL)(RETURN 0)))
       (SETQ A (MINORL X N (LIST (CONS (NREVERSE(INDEX* N)) 1)) NL))
       (SETQ NL NIL)
       (COND ((NULL A)(RETURN 0)))
       (TB2 X  (CAR A)N)
       TAG2
       (SETQ ML*(CONS (CONS NIL NIL)(CAR A)))
       (SETQ A (CDR A))
       (COND ((NULL A) (RETURN (COND ((= *DETSIGN* 1) (CADADR ML*))
				     (T (GPCTIMES -1  (CADADR ML*)))))))
       (SETQ NML (CAR A))
       TAG1(COND((NULL NML)(GO TAG2)))
       (SETQ DD  (CAR NML))
       (SETQ NML (CDR NML))
       (NBN DD)
       (GO TAG1)
       ))

(DEFUN NBN (RULE)
  (declare (special x))
  (PROG (ANS R A)
       (SETQ ANS 0 R (CADAR RULE))
       (COND ((EQUAL R 0) (RETURN 0)))
       (RPLACA RULE (CAAR RULE))
       LOOP(COND((NULL R) (RETURN(RPLACD RULE(CONS ANS (CDR RULE))))))
       (SETQ A (CAR R) R(CDR R))
       (SETQ ANS(GPPLUS ANS
			(GPTIMES
			 (COND ((= (CADR A) 1)
				(AREF X (CADDR A) (CADDDR A)))
			       (T (GPCTIMES (CADR A) (AREF X (CADDR A) (CADDDR A)))))
			 (GETMINOR (CAR A)))))
       (GO LOOP)))

(DEFUN GETMINOR (INDEX)
       (COND((NULL(SETQ INDEX(zl-ASSOC INDEX ML*)))0)
	    (T(RPLACD (CDR INDEX)(f1- (CDDR INDEX)))
	      (COND((= (CDDR INDEX )0)
		    (zl-DELETE INDEX ML*)))
	      (CADR INDEX)))
       )

(DEFUN TB2 (X L N)
  (DECLARE(FIXNUM N ))
;  (setq x (get-array-pointer x))
  (PROG( ( N-1(f1- N)) B A)
  (DECLARE(FIXNUM  N-1))
       LOOP(COND((NULL L) (RETURN NIL)))
       (SETQ A (CAR L) L (CDR L)B (CAR A))
       (RPLACD A (CONS (GPDIFFERENCE(GPTIMES (AREF X (CAR B) N-1) (AREF X (CADR B) N))
				    (GPTIMES (AREF X (CAR B) N) (AREF X (CADR B) N-1)))
		       (CDR A)))
       (GO LOOP)
       ))

(DEFUN ZROW (X I COL N)
       (DECLARE(FIXNUM I COL N ))
;    (setq x (get-array-pointer x))
       (PROG((J COL))
	    (DECLARE(FIXNUM  J))
	    LOOP(COND((> J N)(RETURN T))
		     ((EQUAL (AREF X I J) 0)(SETQ J(f1+ J))(GO LOOP)))
	    ))

(DEFUN NZL (A N)
  (DECLARE(FIXNUM N ))
  
;  (setq a (get-array-pointer a))
       (PROG((I 0)( J (f- N 2)) D L)
	    (DECLARE(FIXNUM  I J))
	    LOOP0(COND((= J 0) (RETURN L)))
	    (SETQ I N)
	    LOOP1(COND((= I 0) (SETQ L (CONS D L)) (SETQ D NIL)(SETQ J (f1- J))(GO LOOP0)))
	    (COND((NOT(EQUAL (AREF A I J) 0))(SETQ D (CONS I D))))
	    (SETQ I (f1- I))(GO LOOP1)
	    ))

(DEFUN SIGNNP (E L)
       (PROG(I)
	    (SETQ I 1)
	    LOOP (COND ((NULL L)(RETURN NIL))
		       ((EQUAL E (CAR L)) (RETURN I)))
	    (SETQ L(CDR L) I (f- I))
	    (GO LOOP)
	    ))

(DEFUN MEMBERCAR (E L)
       (PROG()
	    LOOP(COND((NULL L)(RETURN NIL))
		     ((EQUAL E (CAAR L))(RETURN(RPLACD (CAR L) (f1+ (CDAR L))))))
	    (SETQ L (CDR L))(GO LOOP)
	    ))

(DECLARE-TOP (UNSPECIAL X ML* RZL*))

(DEFUN ATRANSPOSE (A N)
       (PROG(I J D) (SETQ I 0)
	    LOOP1(SETQ I (f1+ I) J I)
	    (COND ((> I N) (RETURN NIL)))
	    LOOP2 (SETQ J (f1+ J))
	    (COND ((> J N) (GO LOOP1)))
	    (SETQ D (AREF A I J))
	    (STORE (AREF A I J) (AREF A J I))
	    (STORE (AREF A J I) D)
	    (GO LOOP2)
	    ))

(DEFUN MXCOMP (L1 L2)
       (PROG()
	    LOOP(COND((NULL L1)(RETURN T))
		     ((CAR> (CAR L1) (CAR L2))(RETURN T))
		     ((CAR> (CAR L2) (CAR L1))(RETURN NIL)))
	    (SETQ L1 (CDR L1) L2 (CDR L2))(GO LOOP)
	    ))

(DEFUN PRMUSIGN (L)
       (PROG((B 0) A D)
	    (DECLARE(FIXNUM B))
	    LOOP (COND((NULL L)(RETURN (COND((EVEN B) 1)(T -1)))))
	    (SETQ A (CAR L) L (CDR L) D L )
	    LOOP1 (COND ((NULL D) (GO LOOP))
			((> A (CAR D)) (SETQ B (f1+ B))))
	    (SETQ D (CDR D))(GO LOOP1)
	    ))

(DEFUN DETPIVOT (X N)
       (PROG(R0 C0)
	    (SETQ C0 (COLROW0 X N NIL) R0(COLROW0 X N T))
	    (SETQ C0 (NREVERSE(BBSORT C0 (FUNCTION CAR>))))
	    (SETQ  R0 (NREVERSE(BBSORT R0 (FUNCTION CAR>))))
	    (COND ((NOT(MXCOMP C0 R0))(ATRANSPOSE X N)(SETQ C0 R0)))
	    (SETQ *DETSIGN* (PRMUSIGN (MAPCAR (FUNCTION CAR) C0)))
	    (NEWMAT 'X* X N C0)
	    (*REARRAY X)))

(DEFUN NEWMAT(X Y N L)
;  (setq y (get-array-pointer y))
  (PROG (I J JL)
	;(set x   (*ARRAY  nil T (f1+ N) (f1+ N)))
	(set x   (*ARRAY  nil T (f1+ N) (f1+ N)))
	(setq x (get-array-pointer x))
	(SETQ J 0.)
   LOOP (SETQ I 0 J (f1+ J))
	(COND ((NULL L) (RETURN NIL)))
	(SETQ JL (CDAR L) L (CDR L))
    TAG (SETQ I (f1+ I))
	(COND ((> I N)(GO LOOP)))
	(STORE (AREF X I J) (AREF Y I JL))
	(GO TAG)))

(DEFUN CAR> (A B) (> (CAR A) (CAR B)))

(COMMENT IND=T FOR ROW ORTHERWISE COL)

(DEFUN COLROW0 (A N IND)
  (DECLARE(FIXNUM N ))
;  (setq a (get-array-pointer a))
  (PROG ((I 0) (J n)  L (C 0))
  (DECLARE(FIXNUM i  C J))
  LOOP0 (COND((= J 0) (RETURN L)))
        (SETQ I N)
  LOOP1 (COND ((= I 0)
	       (SETQ L (CONS (CONS C J) L))
	       (SETQ C 0.)
	       (SETQ J (f1- J))
	       (GO LOOP0)))
        (COND ((EQUAL (COND (IND (AREF A J I))
			    (T (AREF A I J))) 0)
	       (SETQ C (f1+ C))))
       (SETQ I (f1- I))(GO LOOP1)
       ))

(DEFUN GPDIFFERENCE (A B)
       (COND ($RATMX (PDIFFERENCE A B))
	     (T (SIMPLUS(LIST '(MPLUS) A (LIST '(MTIMES) -1 B)) 1 NIL))))

(DEFUN GPCTIMES(A B) (COND ($RATMX (PCTIMES A B)) (T (SIMPTIMES(LIST '(MTIMES) A B) 1 NIL))))

(DEFUN GPTIMES(A B) (COND ($RATMX (PTIMES A B)) (T(SIMPTIMES (LIST '(MTIMES) A B) 1 NIL))))

(DEFUN GPPLUS(A B) (COND ($RATMX (PPLUS A B)) (T (SIMPLUS(LIST '(MPLUS) A B) 1 NIL))))