;;; -*-  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))))

