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(macsyma-module ufact)
13
14(load-macsyma-macros ratmac rzmac)
15
16;; Dense Polynomial Representation
17
18(defun dprep (p)
19  (do ((n (car p))
20       (e (car p) (f1- e))
21       (l))
22      ((< e 0) (cons n (nreverse l)))
23    (cond ((equal e (car p))
24	   (push (cadr p) l)
25	   (setq p (cddr p)))
26	  (t (push 0 l)))))
27
28(defun dpdisrep (l)
29  (cond ((zerop (car l)) (cadr l))
30	((do ((l (nreverse (cdr l)) (cdr l))
31	      (n 0 (f1+ n))
32	      (ll))
33	     ((null l) ll)
34	   (or (= (car l) 0)
35	       (setq ll (cons n (cons (car l) ll))))))))
36
37;; not currently called
38;;(DEFUN PGCDU* (P Q)
39;;       (COND ((OR (PCOEFP P) (PCOEFP Q))  1)
40;;	     ((NULL MODULUS)
41;;	      (merror  "Illegal CALL TO PGCDU"))
42;;	     ((> (CADR P) (CADR Q))
43;;	      (PSIMP (CAR P) (DPDISREP (DPGCD (DPREP (CDR P)) (DPREP (CDR Q))))))
44;;	     ((PSIMP (CAR P) (DPDISREP (DPGCD (DPREP (CDR Q)) (DPREP (CDR P))))))))
45;;
46;;(DEFUN PMODSQFRU (P)
47;;       (DO ((DPL (DPSQFR (DPREP (CDR P))) (CDR DPL))
48;;	    (PL NIL (CONS (PSIMP (CAR P) (DPDISREP (CDAR DPL))) (CONS (CAAR DPL) PL))))
49;;	   ((NULL DPL) PL)))
50
51(defun dpgcd (p q)
52  (if (< (car p) (car q)) (rotatef p q))
53  (do ((p (copy-list p) q)
54       (q (copy-list q) (dpremquo p q nil)))
55      ((= (car q) 0)
56       (if (= (cadr q) 0) p '(0 1)))))
57
58(defun dpdif (p q)
59  (cond ((> (car p) (car q))
60	 (do ((i (car p) (f1- i))
61	      (pl (cdr p) (cdr pl))
62	      (l nil (cons (car pl) l)))
63	     ((= i (car q)) (dpdif1 pl (cdr q) l)) ))
64	((< (car p) (car q))
65	 (do ((i (car q) (f1- i))
66	      (ql (cdr q) (cdr ql))
67	      (l nil (cons (cminus (car ql)) l)))
68	     ((= i (car p)) (dpdif1 (cdr p) ql l))))
69	(t (dpdif1 (cdr p) (cdr q) nil))))
70
71(defun dpdif1 (p1 q1 l)
72  (do ((pl p1 (cdr pl))
73       (ql q1 (cdr ql))
74       (ll l (cons (cdifference (car pl) (car ql)) ll)))
75      ((null pl) (dpsimp (nreverse ll)))))
76
77(defun dpsimp (pl) (setq pl (ufact-strip-zeroes pl))
78       (cond ((null pl) '(0 0))
79	     (t (cons (f1- (length pl)) pl))))
80
81(defun dpderiv (p)
82  (cond ((= 0 (car p)) '(0 0))
83	(t (do ((l (cdr p) (cdr l))
84		(i (car p) (f1- i))
85		(dp nil (cons (ctimes i (car l)) dp)))
86	       ((= i 0) (cons (f1- (car p)) (nreverse dp)))))))
87
88(defun dpsqfr (q)			;ASSUMES MOD > DEGREE
89  (do ((c q (dpmodquo c p))
90       (d (dpderiv q) (dpmodquo d p))
91       (i 0 (f1+ i))
92       (p)
93       (pl))
94      ((= 0 (car c)) pl)
95    (cond (p (setq d (dpdif d (dpderiv c))
96		   p (dpgcd c d))
97	     (and (> (car p) 0)
98		  (setq pl (cons (cons i p) pl))))
99	  (t (setq p (dpgcd c d))
100	     (cond ((= (car p) 0) (return (ncons (cons 1 c)))))))))
101
102
103
104(defun dpmodrem (p q)
105  (cond ((< (car p) (car q)) p)
106	((= (car q) 0) '(0 0))
107	((dpremquo (copy-list p) (copy-list q) nil))))
108
109(defun dpmodquo (p q)
110  (cond ((< (car p) (car q)) '(0 0))
111	((= (car q) 0)
112	 (cond ((equal (cadr q) 1) p)
113	       (t (cons (car p)
114			(mapcar #'(lambda (c) (cquotient c (cadr q))) (cdr p))
115			))))
116	((dpremquo (copy-list p) (copy-list q) t))))
117
118;; If FLAG is T, return quotient.  Otherwise return remainder.
119
120(defun dpremquo (p q flag)
121  (prog (lp lq l alpha)
122     (cond ((= (cadr q) 1)
123	    (setq alpha 1))
124	   (t (setq alpha (crecip (cadr q)))
125	      (do ((l (cddr q) (cdr l)))
126		  ((null l)
127		   (rplaca (cdr q) 1))
128		(rplaca l (ctimes (car l) alpha)))))
129     a    (and flag (setq l (cons (ctimes (cadr p) alpha) l)))
130     (setq lp (cddr p) lq (cddr q))
131     b    (rplaca lp (cdifference (car lp) (ctimes (car lq) (cadr p))))
132     (cond ((null (setq lq (cdr lq)))
133	    (do ((e (f1- (car p)) (f1- e))
134		 (pp (cddr p) (cdr pp)))
135		((null pp) (setq p '(0 0)))
136	      (cond ((signp e (car pp))
137		     (and flag (not (< e (car q)))
138			  (setq l (cons 0 l))))
139		    ((return (setq p (cons e pp))))))
140	    (cond ((< (car p) (car q))
141		   (return (cond (flag (dpsimp (nreverse l)));GET EXP?
142				 (p))))
143		  ((go a))))
144	   (t (setq lp (cdr lp))
145	      (go b)))))
146
147(defun ufact-strip-zeroes (l)
148  (do ((l l (cdr l)))
149      ((null (pzerop (car l))) l)))
150
151(defun cpres1 (a b)
152  (prog (res (v 0) a3) (declare (fixnum v))
153	(setq  a (dprep a) b (dprep b))
154	(setq res 1)
155	again (setq a3 (dpmodrem a b))
156	(setq v (boole boole-xor v (logand 1 (car a) (car b) )))
157	(setq res (ctimes res (cexpt (cadr b)
158				     (f- (car a) (car a3)))))
159	(cond ((= 0 (car a3))
160	       (setq res (ctimes res (cexpt (cadr a3) (car b))))
161	       (return (cond ((oddp v) (cminus res))
162			     (t res))) ))
163	(setq a b)
164	(setq b a3)
165	(go again) ))
166