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

(COMMENT THIS IS THE ALG FACTOR PACKAGE) 
;;; Toplevel functions are: CPBGZASS CPTOM

(load-macsyma-macros ratmac)

(DECLARE-TOP(*LEXPR $FACTOR)
	 (SPECIAL TRA* TRL* *XN VAR INTBS* PLIM MANY* SPLIT* ALC IND P L)
	 (GENPREFIX AFC))


;; (DEFUN FIXMINPOLY NIL
;;        (PROG (ANS)
;; 	     (COND ((NOT (EQUAL (CADDR MINPOLY*) 1))
;; 		    (SETQ MPLC* (CADDR MINPOLY*))
;; 		    (SETQ MINPOLY* (PMONZ MINPOLY*))))
;; 	     (SETQ ANS (CDR MINPOLY*))
;; 	     (COND ((AND (EQUAL (CAR ANS) 2.)
;; 			 (EQUAL (CADR ANS) 1.)
;; 			 (EQUAL (CADDR ANS) 0.)
;; 			 (NOT (EQUAL 1.
;; 				     ((LAMBDA(MODULUS)
;; 					     (CMOD(MINUS (CADDDR ANS))))
;; 				      4.)))
;; 			 (CSQFRP (CADDDR ANS)))
;; 		    (SETQ INTBS* 1.))
;; 		   ((EQUAL (CDR ANS) '(1. 0. -1.)) (SETQ INTBS* 1.)))))


(DEFUN ZIREDUP (P)
       ((LAMBDA (MODULUS ALPHA MINPOLY* ALGFAC* GAUSS TELLRATLIST MANY* MM* $GCD)
		(NULL (CDDR(PFACTOR P))))
	NIL NIL NIL NIL NIL NIL NIL 1. '$EZ))

(DEFUN INTBASEHK (P)
       (PROG (MODULUS)
	     (SETQ MODULUS PLIM)
	     (SETQ P (PCTIMES INTBS* P))
	     (SETQ MODULUS NIL)
	     (RETURN (CAR (RATREDUCE P INTBS*))))) 

(DEFUN FINDIBASE (P) 
       (PROG (MAINVAR) 
	     (SETQ MAINVAR (CAR P))
	     (SETQ P (REDRESULT P (PDERIVATIVE P MAINVAR)))
	     (SETQ P (CFACTORW P))
	     (SETQ MAINVAR 1.)
	LOOP (COND ((NULL P) (RETURN MAINVAR)))
	     (SETQ MAINVAR (TIMES MAINVAR (EXPT (CAR P) (QUOTIENT (CADR P) 2.))))
	     (SETQ P (CDDR P))
	     (GO LOOP)))


(DEFUN CPBGZASS (QLIST V M) 
       (PROG (F Y VJ FACTORS U W LC J P2 FNJ FNQ OLDFAC) 
	     (COND ((EQUAL M 1.) (RETURN (LIST V)))
((EQUAL M (CADR V))  (RETURN((LAMBDA (VAR) (GFSPLIT V))
		     (LIST VAR 1. 1.)))))
	     (SETQ F (PMOD V))
	     (SETQ LC (CADDR F))
	     (SETQ F (MONIZE F))
	     (SETQ P2 1. QLIST (CDR (NREVERSE QLIST)))
	     (SETQ OLDFAC (LIST NIL F))
	NEXTQ(SETQ V (CAR QLIST))
	     (SETQ QLIST (CDR QLIST))
	     (SETQ J (FINDSES V F))
	     (SETQ OLDFAC (NCONC OLDFAC FNQ))
	     (SETQ FNQ NIL)
	INCRJ(SETQ FACTORS (NCONC OLDFAC FNJ))
	     (SETQ FNJ NIL)
	     (SETQ VJ (PPLUS V (CAR J)) J (CDR J))
	TAG2 (SETQ U (CADR FACTORS))
	     (SETQ W (PGCDU VJ U))
	     (COND ((OR (NUMBERP W) (AND ALPHA (ALG W))(= (CADR U) (CADR W))) (GO NEXTFAC)))
	     (SETQ Y (CAR (PMODQUO U W)))
	     (SETQ FNQ (CONS W FNQ))
	     (SETQ FNJ (CONS Y FNJ))
	     (SETQ P2 (f1+ P2))
	     (RPLACD FACTORS (CDDR FACTORS))
	     (COND ((EQUAL P2 M) (GO OUT)) (T (GO TAG1)))
	NEXTFAC
	     (SETQ FACTORS (CDR FACTORS))
	TAG1 (COND ((CDR FACTORS) (GO TAG2)) (J (GO INCRJ)) (QLIST (GO NEXTQ)))
	OUT  (SETQ FNQ (NCONC FNQ FNJ (CDR OLDFAC)))
	     (RETURN (CONS (PTIMES LC (CAR FNQ)) (CDR FNQ))))) 


;; The function PMONZ used to be defined here.  It is also defined in
;; RAT;RAT3A and BMT claims the definitions are equivalent.

(DEFUN FINDSES (G F) 
       (PROG (VAR TRA* TRL*) 
	     (SETQ G (ZASSG (CDR G) (CDR F) (CAR G)))
	     (SETQ VAR (LIST (CAR F) 1. 1.))
	     (SETQ F (GFSPLIT G))
	     (RETURN (MAPCAR (FUNCTION (LAMBDA (A) (CAR (LAST A)))) F)))) 

(DEFUN COEFVEC (P N VEC) 
       (PROG NIL 
	LOOP (COND ((ZEROP N) (RETURN VEC)))
	     (SETQ N (f1- N))
	     (SETQ VEC (CONS (PTERM P N) VEC))
	     (GO LOOP))) 

(DEFUN ZASSG (G F VAR) 
       (PROG (I MAT GN ANS N) 
	     (SETQ N (CAR F))
	     (SETQ GN G)
	     (SETQ I 1. MAT (LIST (COEFVEC '(0. 1.) N (LIST 1.))))
	     (GO ON)
	LOOP (SETQ I (ADD1 I))
	     (SETQ GN (PGCD1 (PTIMES1 GN G) F))
	ON   (SETQ ANS (LINDEP MAT (COEFVEC GN N (LIST (LIST VAR I 1.)))))
	     (COND (ANS (RETURN ANS)))
	     (GO LOOP))) 

(DEFUN DIVL (J A) (MAPCAR (FUNCTION (LAMBDA (L) (CAR (PMODQUO L A)))) J)) 

;; (DEFUN PADDROWS (A B) (MAPCAR (FUNCTION PPLUS) A B))

(DEFUN PDIFROWS (A B) (MAPCAR (FUNCTION PDIFFERENCE) A B)) 

(DEFUN PTIMESROW (VAR ROW) (MAPCAR (FUNCTION (LAMBDA (A) (PTIMES VAR A))) ROW)) 

(DEFUN DDIV (J) 
       (PROG (A B) 
	     (SETQ B J)
	AG   (SETQ A (CAR B))
	     (COND ((ZEROP A) (SETQ B (CDR B)) (GO AG)))
	     (RETURN (DIVL J A)))) 

(DEFUN LINDEP (MAT VEC) 
       (PROG (E D M ROW ROWD VECD) 
	     (SETQ M MAT)
	     (COND ((EQUAL 0. (CAR VEC)) (SETQ VEC (CDR VEC)))
		   (T (SETQ VEC (PDIFROWS (CDR VEC)
					  (PTIMESROW (CAR VEC) (CDAR MAT))))))
	LOOP (COND ((NULL (CDR M))
		    (COND ((ZEROLP (CDR (REVERSE VEC)))
			   (RETURN (CAR (LAST VEC))))
			  (T (RPLACD M (CONS (DDIV VEC) (CDR M)))
			     (RETURN NIL)))))
	     (SETQ ROW (CADR M))
	     (SETQ ROWD ROW VECD VEC)
	LOOP1(SETQ D (CAR ROWD))
	     (SETQ E (CAR VECD))
	     (COND ((EQUAL 0. E)
		    (COND ((EQUAL 0. D)
			   (SETQ VECD (CDR VECD) ROWD (CDR ROWD))
			   (GO LOOP1))
			  (T (SETQ VEC (CDR VEC)) (SETQ M (CDR M)) (GO LOOP))))
		   ((EQUAL 0. D)
		    (RPLACD M
			    (CONS (DIVL VEC E) (MAPCAR (FUNCTION CDR) (CDR M))))
		    (RETURN NIL)))
	     (SETQ VEC (PDIFROWS (CDR VEC) (PTIMESROW E (CDR ROW))))
	     (SETQ M (CDR M))
	     (GO LOOP)))

(DEFUN GFSPLIT (F) 

       (PROG (TR FL (N 0) ANS TRA* (I 0) NFL)
	       (declare(FIXNUM N I))
	     (SETQ FL (LIST F) N (CADR F))
	LOOP (COND ((NULL FL)(COND((NULL NFL)  
(COND ((= N (LENGTH ANS))(SETQ TRL* NIL)(RETURN ANS))(T (merror "GFSPLIT"))))(T(SETQ FL NFL NFL NIL I (f1+ I))))))
	     (SETQ F (CAR FL) FL (CDR FL))
(COND((> I MM*) (merror "GFSPLIT")))

(COMMENT	     (COND ((AND (NULL MODULUS) (NOT ALPHA) (EQUAL 2. (CADR F)))
		    (SETQ ANS (NCONC ANS (PFACTORQUAD F)))
		    (GO LOOP))))

	     (SETQ TR (TRACEMOD0 F I))
   (COND((OR (PCOEFP TR)(AND ALGFAC* (ALG TR)))(SETQ NFL (CONS F NFL))(GO LOOP)))
	     (SETQ F (CPBG0 TR F))
	     (SETQ ANS (NCONC ANS (CAR F)))
	     (COND ((NULL (CDR F)) (GO LOOP)))
	     (SETQ NFL (NCONC NFL (CDR F)))
	     (GO LOOP)))

(DEFUN CPBG0 (TR F) 
       (PROG (M F1 F2 G ALC TRM) 
	     (SETQ M 0.)
	     (COND ((AND (NOT (NUMBERP (CADDR TR))) (ALG (CADDR TR)))
		    (SETQ ALC (PAINVMOD (CADDR TR)) TR (PTIMES ALC TR)))
		   (T (SETQ ALC 1.)))
	BK   (COND ((PCOEFP F) (RETURN
(COND ((AND(NULL F1)(NULL F2))(merror "WRONG TRACE"))(T (CONS F1 F2)))))
		   ((EQUAL (CADR F) 1.) (RETURN (CONS (CONS F F1) F2)))
((EQUAL M MODULUS)(RETURN (CONS F1 (CONS F F2)))))
	     (SETQ TRM (PDIFFERENCE TR (PTIMES M ALC)))
	     (SETQ G (PGCDU TRM F))
	     (COND ((OR (NUMBERP G) (AND ALPHA (ALG G)))
		    (SETQ M (f1+ M))
		    (GO BK)))
	     (SETQ F (CAR (PMODQUO F G)))
	     (COND ((EQUAL (CADR G) 1.) (SETQ F1 (CONS G F1)))
		   (T (SETQ F2 (CONS G F2))))
	     (GO BK)))

(DEFUN CPOL2P (P VAR)
  (PROG((I 0) ANS)
         (declare(FIXNUM I))
       (SETQ  P (NREVERSE P))
    LOOP (COND((NULL P) (RETURN (CONS VAR ANS)))
	      ((EQUAL 0.  (CAR P)) NIL)
	      (T (SETQ ANS (CONS I(CONS (CAR P) ANS)))))
       (SETQ P (CDR P)I(f1+ I))(GO LOOP)
       ))

(DEFUN TRACEMOD (V) 
       (PROG (ANS TR QLARGE TERM) 
	     (SETQ ANS 0 TR (NREVERSE TRL*) TRL* NIL)
(COND ((AND (ATOM (CAAR TR))(NOT (NUMBERP (CAAR TR))))(SETQ QLARGE T)))
	LOOP (COND((NULL TR)(RETURN ANS)))
(SETQ TERM(COND (QLARGE (CAR TR))(T (CPOL2P(CAR TR)V))) TR (CDR TR))
	     (SETQ ANS (PPLUS ANS TERM))
(SETQ TRL* (CONS TERM TRL*))
	     (GO LOOP)))

(DEFUN OTRACEMOD (TERM Q M PRIME) 
       (PROG (ANS I) 
	     (SETQ ANS TERM I 1. TRL* (LIST TERM))
	LOOP (COND ((EQUAL I M) (RETURN ANS)))
	     (SETQ ANS (PPLUS ANS (SETQ TERM (PEXPTMOD TERM PRIME Q))))
(SETQ TRL* (CONS TERM TRL*))
	     (SETQ I (f1+ I))
	     (GO LOOP)))

(DEFUN TRACEMOD0 (Q I) 
  (declare(FIXNUM I))
  (PROG ( L ANS A DL) 
	(COND ((= I 0.)(RETURN (COND (TRL*(TRACEMOD (CAR Q)))
				     (T (OTRACEMOD VAR Q MM* MODULUS)))))
	      (TRL* (SETQ DL TRL* TRL* (MAPCAR (FUNCTION (LAMBDA(X)
							   (CONS (CAR X) (PGCD1 (CDR X) (CDR Q))))) TRL*))))
	(COND (TRA* (GO TAG))(T(SETQ L (CDR TRL*) TRA* (LIST ALPHA) A ALPHA)))
     LOOP(COND ((NULL L)(GO TAG)))
	(SETQ L (CDR L) A (PEXPT A MODULUS) TRA* (CONS A TRA*))
	(GO LOOP)
     TAG
	(SETQ ANS(TRACEMOD1 I TRA* TRL*))
	(COND (DL (SETQ TRL* DL)))
	(RETURN ANS)
	))

(DEFUN TRACEMOD1 (N A L)
       (PROG(ANS)
	    (SETQ ANS 0)
   LOOP     (COND ((NULL L)(RETURN ANS)))
            (SETQ ANS (PPLUS ANS (PTIMES(PEXPT (CAR A) N) (CAR L))))
	    (SETQ L (CDR L) A (CDR A))(GO LOOP) ))


;; The way arrays are manipulated has been changed to make this code reentrant.
;; Previously, arrays were kept on the array properties of symbols.  Now, the
;; symbols are bound to the arrays, so they can be rebound on re-entry.
;; The ANOTYPE, INVC, and FCTC arrays are set up in RAT;FACTOR.
;; 
;; (DECLARE (ARRAY* (NOTYPE A 2 INVC 1 FCTC 1)))

(declare-top (special anotype invc fctc))
(defmacro a (row col) `(arraycall t anotype ,row ,col))
(defmacro invc (ind) `(arraycall t invc ,ind))
(defmacro fctc (ind) `(arraycall t fctc ,ind))

(DEFUN CPTOMEXP (P M U N)
(PROG(B1 B2 J N-1 I L)
(SETQ B1 (X**Q1(LIST (CAR U) 1 1) U M P))
(COND ((EQUAL (CDR B1) '(1 1)) (SETQ SPLIT* T)(RETURN NIL)))
(SETQ B2 B1 J 1. N-1 (f1- N))
(GO TAG1)
TAG(SETQ J(f1+ J))
  (COND ((= J N)(RETURN NIL)))
(SETQ B1 (PMODREM(PTIMES B1 B2) U))
TAG1 (SETQ L (P2CPOL B1 N-1) I N-1)
	sharp2   (COND ((NULL L) (GO ON)))
	     (STORE (A J I) (CAR L))
	     (SETQ L (CDR L))
	     (SETQ I (f1- I))
	     (GO sharp2)
	ON   (STORE (A J J) (PDIFFERENCE (A J J) 1))
	     (GO TAG)))

(DEFVAR THR* 100.)

(DEFUN CPTOM (P M U N) 
  #-Multics (DECLARE (FIXNUM N M  P))
  #+Multics (DECLARE (FIXNUM N M Q I J))
  (PROG (( Q (EXPT P M)) L S *XN (J 0) (I 0) IND N-1)
	(declare (special Q I J))
	(SETQ  N-1 (f1- N))
	(COND ((> Q THR*) (RETURN (CPTOMEXP P M U N))))
   LOOP (SETQ J (f1+ J))
	(COND ((= J N) (RETURN NIL))
	      (IND (GO SA))
	      (T
	       (SETQ *XN (MAPCAR #'PMINUS (P2CPOL (CDDR U) N-1))
		     S (X**Q (P2CPOL(LIST VAR 1 1) N-1) P M)
		     IND T)))
	(GO ST)
     SA (CPTIMESXA S Q)
     ST (COND ((AND (= J 1)
		    (EQUAL '(1 0) (NCDR S (f1- (LENGTH S) )))
		    (= 1 (LENGTH (zl-DELETE 0 (COPY S)))))
	       (RETURN (SETQ SPLIT* T))))
        (SETQ L S)
	(SETQ I N-1)
	(COMMENT(PRINT L))
	sharp2   (COND ((NULL L) (GO ON)))
	     (STORE (A J I) (CAR L))
	     (SETQ L (CDR L))
	     (SETQ I (f1- I))
	     (GO sharp2)
	ON   (STORE (A J J) (PDIFFERENCE (A J J) 1))
	     (GO LOOP)))

(DEFUN CPTIMESXA (P I) 
  (declare(FIXNUM I))
  (PROG (XN Q LC) 
     AG    (COND ((= I 0.) (RETURN P)))
	(SETQ XN *XN Q P LC (CAR P))
     LOOP (COND ((CDR Q)
		 (RPLACA Q (PPLUS (CADR Q) (PTIMES LC (CAR XN))))
		 (SETQ Q (CDR Q) XN (CDR XN)))
		(T (RPLACA Q (PTIMES LC (CAR XN))) 
		   (SETQ I(f1- I))(GO AG)))
	(GO LOOP)))

(DEFUN X**Q (X P M) 
  #-Multics (DECLARE (FIXNUM   M P))
  #+Multics (DECLARE (FIXNUM D PP M I))
  (PROG ((I 1)  (PP 1) (D 0))
	(declare (fixnum i pp d))
	(SETQ I 1. TRL* (LIST X) PP 1.)
   LOOP (COND ((= I M) (RETURN (CPTIMESXA X (f- (f* PP P) PP)))))
	(SETQ D PP)
	(CPTIMESXA X (f- (SETQ PP(f* PP P)) D))
	(SETQ TRL* (CONS(COPY X) TRL*))
	(SETQ I (f1+ I))
	(GO LOOP)))

(DEFUN CMNULL (N) 
       (DECLARE (FIXNUM N ))
       (PROG (NULLSP 	     (SUB1N (f1- N)) MONE (K 1) (J 0) (S 0) NULLV (I 0) VJ M AKS)
	     (declare (fixnum sub1n i j k s))
	     (SETQ MONE (CMOD -1))
	SHARP    (COND ((> I SUB1N) (GO ON)))
	     (STORE (FCTC I) -1)
	     (STORE (INVC I) -1)
	     (SETQ I (f1+ I))
	     (GO SHARP)
	ON   (SETQ K 1 NULLSP (LIST 1))
	N2   (COND ((> K SUB1N) (RETURN NULLSP)))
	     (SETQ J 0)
	N3A  (COND ((> J SUB1N) (GO NULL))
		   ((OR (EQUAL (A K J) 0) (GREATERP (FCTC J) -1))
		    (SETQ J (f1+ J))
		    (GO N3A)))
	     (STORE (INVC K) J)
	     (STORE (FCTC J) K)
	     (SETQ M (A K J))
	     (SETQ M (PTIMES MONE (PAINVMOD M)))
	     (SETQ S K)
	sharp1   (COND ((> S SUB1N) (GO ON1)))
	     (STORE (A S J) (PTIMES M (A S J)))
	     (SETQ S (f1+ S))
	     (GO sharp1)
	ON1  (COMMENT (GO THROUGH COLUMNS))
	     (SETQ S 0)
	sharp2   (COND ((> S SUB1N) (GO NEXTK)))
	     (COMMENT (GO THROUGH ROWS IN EACH COLUMN))
	     (COND ((= S J) NIL)
		   (T (PROG (I) 
			    (SETQ AKS (A K S))
			    (SETQ I K)
		       sharp3   (COND ((> I SUB1N) (RETURN NIL)))
			    (STORE (A I S)
				   (PPLUS (A I S)
					  (PTIMES (A I J) AKS)))
			    (SETQ I (f1+ I))
			    (GO sharp3))))
	     (SETQ S (f1+ S))
	     (GO sharp2)
	NULL (SETQ NULLV NIL)
	     (SETQ S 0)
	sharp4   (COND ((> S SUB1N) (GO ON4))
		   ((= S K) (SETQ NULLV (CONS S (CONS 1 NULLV))))
		   ((> (INVC S) -1)
		    (SETQ VJ (A K (INVC S)))
		    (COND ((EQUAL VJ 0) NIL)
			  (T (SETQ NULLV (CONS S (CONS VJ NULLV))))))
)	     (SETQ S (f1+ S))
	     (GO sharp4)
	ON4  (COND ((EQUAL (CAR NULLV) 0) (SETQ NULLV (CADR NULLV)))
		   ((SETQ NULLV (CONS VAR NULLV))))
	     (SETQ NULLSP (CONS NULLV NULLSP))
	NEXTK(SETQ K (f1+ K))
	     (GO N2)))