  -*- Mode : Common - Lisp ; Package : RAT ;  -*-

;;
  Written by :    ,   
;; File: rat.lisp
;; Contains: rational function (in partially-factored form) arithmetic 
(provide 'rat1)

  ( c ) Copyright 1990 , 1991    ,   

  A fpe is a ( possibly ) factored polynomial expression implemented
;;; as a list of pairs of polynomials and exponents. For example,
;;;
                              2   2          4
                    12 ( x + 1 ) ( y   + y - 5 )
;;;
  is represented as a list of three pairs , the first pair being the
  the number 12 and the exponent 1 , the second being the
  polynomial ( x + 1 ) and the exponent 2 , and the third being the
  polynomial ( y^2 + y - 5 ) and the exponent 4 .
  Currently , the first pair is always   an integer with exponent 1 .

;;; Although the polynomials could have coefficients in any ring
  for some operations , the idea of gcd as used here assumes coefs form a
;;; unique factorization domain, and therefore it is principally useful
  to have integer coefficients .   ( Rationals are not a UFD since
   1/2 * 2 = 1/3 * 3 = 1 etc . )

;;; Getting back to our representation --
  There is a first pair that always looks like ( integer . 1 ) .
;;; Subsequent pairs consist of
;;; distinct polynomials with integer coefficients, and
  exponents which are integers > = 1 .
;;; It is not required that polys  be squarefree, monic, irreducible
;;; or restricted in some other ways.  Examples of possible pairs:

   ( ( x^2 ) . 2 )   ;    ( ( 3 x^2 - 3 x ) . 1 )
;;;  actually, the polynomials will be in some better form, so it is
;;; more accurate to depict these examples as 
   ( # ( 1 0 0 1 ) . 2 )   ;   ( # ( 1 0 -3 3 ) . 1 )    where x < -- > 1

  The polynomials in an fpe , except for intermediate forms produced
;;; during gcd computation, are maintained in lexicographically
  increasing order . Note that two polynomials in fpe form
;;; may not be IDENTICAL even though they are mathematically
;;; equivalent, because of different factorings.

  For example ,   ( ( 1 . 1 ) ( ( x^2 + 2x+1 ) . 1 ) )   and ( ( 1 . 1 ) ( ( x+1 ) . 2 ) )
;;; are equivalent but not identical.  If they were added together,
;;; the result would be  (if we refrained from computing
;;; gcds)
   ( ( 1 . 1 ) ( ( 2x^2 + 4*x+2 ) . 1 ) )

;;; Another equivalent answer (more effort to compute .. we don't do
;;; this...) would get
   ( ( 2 . 1 ) ( ( x+1 ) . 2 ) )

  On the other hand , if we were adding ( ( 1 . 1 ) ( ( x+1 ) . 2 ) ) and
  ( ( 1 . 1 ) ( ( x+1 ) . 3 ) ) , the common factor of x+1 would be obvious
  and we would get   ( ( 1 . 1 ) ( ( x+1 ) . 2 ) ( ( x+2 ) . 1 ) )

(proclaim '(optimize (speed 3) (safety 0) (space 0) (compilation-speed 0)))

(require 'poly) ;; need macro defs
(in-package :mma)
;; export what other modules may need
#+ignore (export '(make-fpe fpe-expand fpe-coef-p fpe-coefone-p fpe-negativep
		   fpe-insert fpe-coefzero-p
		   make-rat rat-p rat make-good-rat rat-numerator 
		   rat-denominator  rat+ rat* rat^ rat/ poly2rat
		   *rat-gcd-level*))

;; default is to compute the gcd, but only between factors of
;; the numerator with factors of the denominator.

(defvar *rat-gcd-level* 'gcd) 
(defvar *expand* nil)

;; other possible settings:
   cdi : just find the obvious common factors by identity . That is
;;  u and v have a common divisor iff u = v.

  cdd : common divisor by division .   That is , u and v have a
  common divisor if one divides the other exactly .

  current implementation uses gcd as default .


  NOT IMPLEMENTED techniques :


;; supergcd: all gcds of factors in numerator OR denominator are
  considered for pair - wise gcds .   That is , ( x^2 - 1)*(x+1 ) -- > ( x+1)^2*(x-1 )

;; square-free: all factors are run through a square-free factorization
;; program. (they are tagged as square-free to save redundant calculations).

;; factors: all are factored over the integers. (also tagged)

  make - fpe : make a " normal " fpe p^e from a polynomial p , exponent e>=0
;; examples:

   ( make - fpe poly 4 ) -- >   ( ( 1 . 1 ) ( poly 4 ) ) ; ; usual case
   ( make - fpe 0 1 )     -- >   ( ( 0    . 1 ) )
    by virtue of using fpe - insert with monomial flag = t , if
;;   make-fpe is given something like p = 3x^2y^3, e=2, it will
    produce ( ( 9 . 1 ) ( # ( x 0 1 ) . 4 ) ( # ( y 0 1 ) . 6 ) ) .
;;   {actually, you won't see x and y, but some number-encoding of them.}

(defun make-fpe (p e) 
  (cond ((coefp p) 
	 (list (cons (coef^ p e) 1)))
	(t (fpe-insert p e '((1 . 1)) t))))

(defun poly2rat (p e)
  (cond ((>= e 0)(make-rat :numerator (make-fpe p e)
			  :denominator (list '(1 . 1))))
	(t (make-good-rat (list '(1 . 1))
		      (make-fpe p (- e))))))
	
  fpe - copy : return a copy of the fpe u

(defun fpe-copy (u) (copy-tree u))

  fpe - expand : expand the fpe into a fully - expanded form ; the result
;;             is returned as a polynomial

(defun fpe-expand (u)
  (let ((*expand* t))
  (do ((ul u (cdr ul)) (p 1))
      ((null ul) p)
      (setq p (p* p (p^ (caar ul) (cdar ul)))))))

  A fpe is " normal " if its first factor is an integer and it has
;;; no other integer factors. Also, the coefficient of all
  monomial factors should be 1 . ( reason : they look odd otherwise ...
  consider 3*x^2 * 4*y^2   instead of 12*x^2*y^2 )
;;; The factors are unique, and
;;; sorted. All coefficients in the domain
  are represented as c^1 , including coefzero and coefone .

  fpe - norm : normalize a fpe

;;; this program could be reprogrammed to do much less consing.

(defun fpe-norm (u &aux ctest (const 1))
  (labels
   ((fn1(ul) ; aux function to sort terms.
	(cond ((null ul) nil)
  stick coefs up front
	       (setq const (coef* const
				  (coef^ (caar ul)(cdar ul))))
	       (fn1 (cdr ul)))
	      (t 
	       (fn1-insert (car ul) (fn1 (cdr ul))))))
    (fn1-insert (h ul)
		(cond ((null ul)
		       ;; Since factor (car h) is nowhere to be seen further in
  this fpe , make the pair of this factor and its power
		       ;; appear at the end of the list.
		       (cons h ul))
		      ((eq 'e (setq ctest (pcompare (car h) (caar ul))))
		       ;; We found this factor. Increment the exponent
		       (cons (cons (car h)(+ (cdr h)(cdar ul))) (cdr ul)))
		      ;; Otherwise we keep looking for the factor.
		      ((eq 'l ctest)
  this poly is less in ordering that (   ) ,
		       ;; so place it right here.
		       (cons h ul))
		      ;; This keep searching 
		      (t 
		       (cons  (car ul)(fn1-insert h (cdr ul)))
			 ul)))
  end of local fns
			  
   (setq u (fn1 u))
   (cons  (cons const 1) u)))
			 
  
(defun fpe-coef-p (u) (null(cdr u)))

(defun fpe-coefone-p (u) (and (null (cdr u)) (coefonep (caar u))))

(defun fpe-coefzero-p (u) (coefzerop (caar u)))



  fpe - insert : insert a pair of polynomial and exponent into the fpe u.
;;             This checks to see whether the polynomial already exists in
              the fpe ; if yes , increment the exponent ; otherwise
              adds the pair to the fpe . This will not change 

  in effect , this multiplies u in fpe form , by poly^u .
  compare to fpe- * , which multiplies 2 polys in fpe form .
  if monom is t , then maybe poly is a monomial . Check it
;; and insert it in pieces if appropriate.

(defun fpe-insert (poly exp u &optional (monom nil) &aux ctest)
  (labels
   ((fi1(ul)
	     ;; this routine recurses down the list of factors
	     (cond ((null ul)
		    ;; Since this factor is nowhere to be seen in
  this fpe , insert the new factor and its power
		    ;; at the end of the list.
		    (list (cons poly exp)))
		   
		   ((eq 'e (setq ctest (pcompare poly (caar ul))))
		    ;; We found this factor. Increment the exponent.
		    (cons (cons (caar ul) (+ exp (cdar ul)))(cdr ul)))
		   ;; Otherwise we keep looking for the factor.
		   ((eq 'l ctest)
  this poly is less in ordering than (   ) ,
		    ;; so place it right here.
		    (cons (cons poly exp) ul))
		   ;; This keep searching 
		   (t (cons (car ul)(fi1 (cdr ul)))))))
	(cond ((coefp poly)
  inserting an integer factor : modify the first element
	       (cond ((coefonep poly) u)
		     ((coefzerop poly) (make-fpe 0 1))
		     (t (cons (cons (coef* (caar u)
						 (coef^ poly exp))
					  1)
			      (cdr u)))))
multiplying u by z^0 gives u
	      ((and monom (monomialp poly)(or (> (degree poly) 1)
					      (not (coefonep (lc poly)))))
  the poly is something like ( 3*x^2 ) ^4
  insert 81   ( i.e. 3 ^ 4 ) into   1*(x)^8 .
  should work recursively for ( 3*y^2*x^2)^4
	       
	       (fpe-insert (lc poly) exp (fpe-insert 
					  (vector (svref poly 0) 0 1)
					  (* (degree poly) exp) u )
			   monom))
	      
	      ;; we can, more generally, not insist on monomial,
	      ;; but only 0 x^1 term.  That is,
	      ;; (a*x^5+b*x^2) --> x^2*(a*x^3+b).  or
	      ;; (x^5+x^2)^3 --> x^6*(x^3+1)^3. etc.
	      ;; note that 0 constant term is not useful because
	      ;;  x is encoded as (vector n 0 1) and has 0 const term

	      ((and monom 
		    (not *expand*)
		    (coefzerop(constc poly))
		    (> (length (the simple-vector poly)) 3)
		    (coefzerop (svref poly 2)))
	       (do ((i 2 (1+ i)))
		   ((null (coefzerop (svref poly i)));there's a non-0
		    (fpe-insert
		     (vector (svref poly 0) 0 1)
		     (* exp (1- i))  ;; the degree of the factor
		     (fpe-insert (polyshift poly (1- i)) exp u)))
		   (declare (fixnum i))
		   ;; no do-body
		   ))
		    		    
	      (t (fi1 u)))))

(defun polyshift(p i)  ;; this divides a poly in x by x^i
  (let ((z(subseq p i)))
    (setf (svref z 0) (svref p 0)) z))


  fpe- *   multiplies two fpe 's u and v. To be canonical ,
;;identical polynomials will appear only once in the result.
;; Both u and v are preserved.  No wasted conses.

(defun fpe-* (u v &aux ctest)
  (labels ((fm1(ul vl)
	    (cond ((null ul) vl)
		  ((null vl) ul)
		  ((eq 'e (setq ctest (pcompare (caar ul) (caar vl))))
		   (cons (cons (caar ul) (+ (cdar ul)(cdar vl)))
			 (fm1 (cdr ul)(cdr vl))))
		   ((eq 'l ctest)
  (   ) is less in ordering that (   ) ,
		    ;; so place it right here.
		    (cons (car ul) (fm1 (cdr ul) vl)))
		   ;; Otherwise (car vl) goes here
		   (t (cons (car vl)(fm1 ul (cdr vl)))))))

  u is coefone
  v is coefone
	      
	      (t (cons (cons (coef* (caar u)(caar v)) 1) ; mult consts.
		       (fm1 (cdr u)(cdr v)))))))

powering an fpe is " easy "
  ;; presumably n is a positive integer, although if it is some
  ;; other number, this program won't break.
  (labels((fp^1 (a)
	      (cond ((null a)nil)
		    (t (cons (cons (caar a) (* n (cdar a)))
			     (fp^1 (cdr a)))))))
       (cons (cons (expt (caar a) n) 1)
	     (fp^1 (cdr a)))))

(defun fpe-negative-p (u) (coef-negative-p (caar u)))

(defun fpe-negate (u)
  (cons (cons (coefneg (caar u))(cdar u))(cdr u)))

  These procedures are concerned with finding the gcd and cd of fpe 's .
  They are independent of the implementation of  's .

  fpe - gcd - cofac :
  returns 3 values : g= the gcd of u and v ;
;;                   ubyg = u/g;
;;                   vbyg = v/g.

;; u and v are unchanged, though any of the outputs may share structure.
  all values are fpe in normal form .

(defun fpe-gcd-cofac (u v)
  (setq u (append u nil) v (append v nil));; copy only top levels of inputs
   g= gcd = 1 , initially

    (do ((i u (cdr i))) 
	;; for each polynomial in u do the following
	;; until we run off the end of u. Then return the triple 
  < gcd , u / gcd , v / gcd > .
	
	((null i) (values g (fpe-norm u)(fpe-norm v)))
  do this only if ( caar i ) is not 1

	(do ((j v (cdr j)));; for each polynomial in v do the following
	    ((null j))
	    (if (not (coefonep (caar j))) ;otherwise next j
	    (let* ((up (caar i))	;  next poly in u
		   (vp (caar j))	;  next poly in v
		   (ue (cdar i))
		   (ve (cdar j)))
	      (multiple-value-bind 
gcd , up / gp , vp / gp
	       (pgcdswitch-cofac up vp)
	       (cond ((coefonep gp))	; no non-trivial divisor! just advance i,j.
		     (t
		      ;; ok, we've found a common factor gp
		      (let ((ge (min ue ve)));; gcd is really gp^ge
			;; add the factor to the ones already found.
			;; incidentally, this factor may already be
  in g because of some earlier discovered gcd .
			(setq g (fpe-insert gp ge g t));; g = g*gp^ge
			
  blot out the two factors in u and v.
			;; we may have to put some remnants at the end, though.
			;; splice in what's left of powers of up and powers
			;; of vp, if any, and then consider up/gp and vp/gp
			
			;; possible missed factorization here.. if we compute
			;; gcd of x^2-1 and (x^2*y^2-y^2)^3 we get a 
			;; co-factor of y^2.  This is a monomial. It would be
			;; better to put y^6 on the factor list, rather than
  ( y^2)^3 .   Could check in fpe - norm for monomials ,
			;; but this would add to expense.

			(cond ((= ge ue)
			       (setf (car i) (cons uq ue))
			       (setq up uq))
			      (t ; that is, ge < ue
			       ;; first decrement the exponent on this factor
			       ;; to make it reflect the number of times
			       ;; we've divided out gp
			       (setf (car i)(cons gp (- ue ge)))
			       (setq up gp)
			       (if (coefonep uq) nil
				 (nconc i (list (cons uq ue))))))
			(cond ((= ge ve)
			       (setf (car j) (cons vq ve))
			       (setq vp vq))
  that is , ge < ve
			       (setf (car j)(cons gp (- ve ge)))
			       (setq vp gp)
			       (if (coefonep vq) nil
				 (nconc j (list (cons vq ve))))))
			)))))))))))

(defun pgcdswitch-cofac(u v)
  (case *rat-gcd-level*
	(gcd 
compute the gcd . Which algorithm depends perhaps
	 ;; on other switches in polynomial package.
	 (pgcd-cofac u v))
	(cdi
	 ;; look for identical common factors only
	 ;; equalp checks element by element in arrays
	 (if (equalp u v)(values u 1 1)(values 1 u v)))
	(cdd
	 ;; look to see if u divides v or vice versa
	 (let (r)
	   (cond ((setq r (p/-test u v)) (values v r 1))
		 ((setq r (p/-test v u)) (values u 1 r))
		 (t (values 1 u v)))))
	   
	;; put other choices here
  use gcd as default
	 (pgcd-cofac u v))))
	 

;;; A rat is a rational polynomial implemented as a structure
  consisting of two fpe 's : the numerator and the denominator .

(defstruct rat numerator denominator)

;;; A "good" rat is one whose denominator is always positive.
;;; If the rat as a whole is negative, the negative sign
;;; is in the numerator.

;; make-good-rat

(defun make-good-rat (n d)
  (cond ((fpe-negative-p d)
	 (setq d (fpe-negate d))
	 (setq n (fpe-negate n))))
  (make-rat :numerator n :denominator d))

;;; The following procedures perform simple arithmetic on rats.
;;; The returned rats from these procedures will be "normal,"
;;; provided that the arguments are "normal." A rat is said to
;;; be "normal" if
;;;       1) it is in reduced form (no common factors between numerator
;;;          and denominator).
        2 ) each polynomial in a fpe that makes up the numerator
;;;          or the denominator appears only once within that rat.
        3 ) the fpe 's are themselves normal : leading coefficient followed
;;;          by terms sorted etc.

  rat * : non - destructively multiply two rats r1 and r2
        reference -      paper


(defun rat* (r1 r2) ;;  r1 * r2 = a/b * c/d
  (let ((a (rat-numerator r1))    
	(b (rat-denominator r1))  
	(c (rat-numerator r2))    
	(d (rat-denominator r2))  
	g num den)
  we assume that gcd(a , b)=1 , and also that gcd(c , d)=1 .
    ;; This may not be, strictly speaking, true, depending on *rat-gcd-level*.
    (if (or (fpe-coefzero-p a)
	    (fpe-coefzero-p c))
	    (return-from rat* (poly2rat 0 1)))
  First , if a and c have any factors in common ,
    ;; set them aside, because if g is a factor of a, then gcd(g,b)
  is predictably 1 , so no need to compute it .  
    ;; We could set aside all common factors regardless of multiplicity
  as an additional optimization . (  -- this is what you were
    ;; doing, I guess).  An alternative would be to save all poly gcd
    ;; results in an eq-hash table..what do you think?  We'd have
    ;; to be careful not to make equal but not-eq polynomials... RJF
    
    (multiple-value-setq (num  a c) (fpe-gcd-cofac a c))
    ;; Ditto for b and d
    (multiple-value-setq (den  b d) (fpe-gcd-cofac b d))
  we can conclude that gcd(num , den)=1 by construction , and
    ;; that with these new values, r1*r2 = (num^2*a*c)/(den^2*b*d),
  subject , however , to more gcd removal . Next ,
    ;; remove the common factors from a and d,  then from b and c.
    (multiple-value-setq (g a d)(fpe-gcd-cofac a d))
    (multiple-value-setq (g b c)(fpe-gcd-cofac b c))
    ;; put back the factors that were removed
    (setq a (fpe-* a (fpe-^ num 2)))
    (setq b (fpe-* b (fpe-^ den 2)))
    ;; finally put the pieces together
    (make-good-rat  (fpe-* a c) (fpe-* b d))))

;; division. almost same as *.  
(defun rat/ (r1 r2) ;;  r1 * r2 = a/b  / d/c
  (let ((a (rat-numerator r1))    
	(b (rat-denominator r1))  
	(d (rat-numerator r2))      ;; note change from rat*
	(c (rat-denominator r2))    ;;  ditto
	g num den)
    (if (fpe-coefzero-p a) (return-from rat/ (poly2rat 0 1)))
    (if (fpe-coefzero-p d) (error "Rational division by zero"))

    (multiple-value-setq (num  a c) (fpe-gcd-cofac a c))
    (multiple-value-setq (den  b d) (fpe-gcd-cofac b d))

    (multiple-value-setq (g a d)(fpe-gcd-cofac a d))
    (multiple-value-setq (g b c)(fpe-gcd-cofac b c))

    (setq a (fpe-* a (fpe-^ num 2)))
    (setq b (fpe-* b (fpe-^ den 2)))

    (make-good-rat  (fpe-* a c) (fpe-* b d))))


  rat+ : non - destructively add two rats r1 and r2

(defun rat+ (r1 r2)
  (let ((a (rat-numerator r1))
	(b (rat-denominator r1))
	(c (rat-numerator r2))   ;; want to perform
	(d (rat-denominator r2)) ;;  a/b + c/d
	(n (make-fpe (coefone) 1))
	num den	g)
    (if (fpe-coefzero-p a) (return-from rat+ r2))
    (if (fpe-coefzero-p c) (return-from rat+ r1))
    
  extract common factors from a and b , also b and d.
    (multiple-value-setq
     (num a c)(fpe-gcd-cofac a c))
    (multiple-value-setq
     (den b d)(fpe-gcd-cofac b d))
    
  now   ( num / den ) * ( a / b+c / d ) and gcd(b , d)=1 .
    ;; n=ad+bc
    (setq n (pnorm(p+ (fpe-expand (fpe-* a d));  n is a polynomial
		(fpe-expand (fpe-* b c)))))
    (if (coefzerop n) (return-from rat+ (make-rat :numerator '((0 . 1))
						  :denominator '((1 . 1)))))
  now n is an fpe .

  set denom to bd * den
        
    ;; remove common factors between  ad+bc and bd*den
    (multiple-value-setq (g n den)(fpe-gcd-cofac n den))
    (setq num (fpe-* num n)); set num to num *(ad + bc)
    (make-good-rat num den)))

;; rat^: raise r to the power e

(defun rat^ (r e)
	 (cond ((= e 0) (make-rat :numerator '((1 . 1))
				  :denominator '((1 . 1))))
	       ((> e 0)
		(make-rat :numerator (fpe-^ (rat-numerator r) e)
			  :denominator (fpe-^ (rat-denominator r) e)))
	       (t (if (fpe-coefzero-p (rat-numerator r))
		      (error "Rational division by zero"))
check sign of  ..
		    (fpe-^ (rat-denominator r) (- e))
		    (fpe-^ (rat-numerator r) (- e))))))



 Package : RAT ;  -*-

 File: rat.lisp
 Contains: rational function (in partially-factored form) arithmetic 
 as a list of pairs of polynomials and exponents. For example,


 Although the polynomials could have coefficients in any ring
 unique factorization domain, and therefore it is principally useful
 Getting back to our representation --
 Subsequent pairs consist of
 distinct polynomials with integer coefficients, and
 It is not required that polys  be squarefree, monic, irreducible
 or restricted in some other ways.  Examples of possible pairs:
    ( ( 3 x^2 - 3 x ) . 1 )
  actually, the polynomials will be in some better form, so it is
 more accurate to depict these examples as 
   ( # ( 1 0 -3 3 ) . 1 )    where x < -- > 1
 during gcd computation, are maintained in lexicographically
 may not be IDENTICAL even though they are mathematically
 equivalent, because of different factorings.
 are equivalent but not identical.  If they were added together,
 the result would be  (if we refrained from computing
 gcds)
 Another equivalent answer (more effort to compute .. we don't do
 this...) would get
 need macro defs
 export what other modules may need
 default is to compute the gcd, but only between factors of
 the numerator with factors of the denominator.
 other possible settings:
  u and v have a common divisor iff u = v.
 supergcd: all gcds of factors in numerator OR denominator are
 square-free: all factors are run through a square-free factorization
 program. (they are tagged as square-free to save redundant calculations).
 factors: all are factored over the integers. (also tagged)
 examples:
 ; usual case
   make-fpe is given something like p = 3x^2y^3, e=2, it will
   {actually, you won't see x and y, but some number-encoding of them.}
 the result
             is returned as a polynomial
 no other integer factors. Also, the coefficient of all
 The factors are unique, and
 sorted. All coefficients in the domain
 this program could be reprogrammed to do much less consing.
 aux function to sort terms.
 Since factor (car h) is nowhere to be seen further in
 appear at the end of the list.
 We found this factor. Increment the exponent
 Otherwise we keep looking for the factor.
 so place it right here.
 This keep searching 
             This checks to see whether the polynomial already exists in
 if yes , increment the exponent ; otherwise
 and insert it in pieces if appropriate.
 this routine recurses down the list of factors
 Since this factor is nowhere to be seen in
 at the end of the list.
 We found this factor. Increment the exponent.
 Otherwise we keep looking for the factor.
 so place it right here.
 This keep searching 
 we can, more generally, not insist on monomial,
 but only 0 x^1 term.  That is,
 (a*x^5+b*x^2) --> x^2*(a*x^3+b).  or
 (x^5+x^2)^3 --> x^6*(x^3+1)^3. etc.
 note that 0 constant term is not useful because
  x is encoded as (vector n 0 1) and has 0 const term
there's a non-0
 the degree of the factor
 no do-body
 this divides a poly in x by x^i
identical polynomials will appear only once in the result.
 Both u and v are preserved.  No wasted conses.
 so place it right here.
 Otherwise (car vl) goes here
 mult consts.
 presumably n is a positive integer, although if it is some
 other number, this program won't break.

                   ubyg = u/g;
                   vbyg = v/g.
 u and v are unchanged, though any of the outputs may share structure.
 copy only top levels of inputs
 for each polynomial in u do the following
 until we run off the end of u. Then return the triple 
 for each polynomial in v do the following
otherwise next j
  next poly in u
  next poly in v
 no non-trivial divisor! just advance i,j.
 ok, we've found a common factor gp
 gcd is really gp^ge
 add the factor to the ones already found.
 incidentally, this factor may already be
 g = g*gp^ge
 we may have to put some remnants at the end, though.
 splice in what's left of powers of up and powers
 of vp, if any, and then consider up/gp and vp/gp
 possible missed factorization here.. if we compute
 gcd of x^2-1 and (x^2*y^2-y^2)^3 we get a 
 co-factor of y^2.  This is a monomial. It would be
 better to put y^6 on the factor list, rather than
 but this would add to expense.
 that is, ge < ue
 first decrement the exponent on this factor
 to make it reflect the number of times
 we've divided out gp
 on other switches in polynomial package.
 look for identical common factors only
 equalp checks element by element in arrays
 look to see if u divides v or vice versa
 put other choices here
 A rat is a rational polynomial implemented as a structure
 A "good" rat is one whose denominator is always positive.
 If the rat as a whole is negative, the negative sign
 is in the numerator.
 make-good-rat
 The following procedures perform simple arithmetic on rats.
 The returned rats from these procedures will be "normal,"
 provided that the arguments are "normal." A rat is said to
 be "normal" if
       1) it is in reduced form (no common factors between numerator
          and denominator).
          or the denominator appears only once within that rat.
          by terms sorted etc.
  r1 * r2 = a/b * c/d
 This may not be, strictly speaking, true, depending on *rat-gcd-level*.
 set them aside, because if g is a factor of a, then gcd(g,b)
 We could set aside all common factors regardless of multiplicity
 doing, I guess).  An alternative would be to save all poly gcd
 results in an eq-hash table..what do you think?  We'd have
 to be careful not to make equal but not-eq polynomials... RJF
 Ditto for b and d
 that with these new values, r1*r2 = (num^2*a*c)/(den^2*b*d),
 remove the common factors from a and d,  then from b and c.
 put back the factors that were removed
 finally put the pieces together
 division. almost same as *.  
  r1 * r2 = a/b  / d/c
 note change from rat*
  ditto
 want to perform
  a/b + c/d
 n=ad+bc
  n is a polynomial
 remove common factors between  ad+bc and bd*den
 set num to num *(ad + bc)
 rat^: raise r to the power e

  Written by :    ,   
(provide 'rat1)

  ( c ) Copyright 1990 , 1991    ,   

  A fpe is a ( possibly ) factored polynomial expression implemented
                              2   2          4
                    12 ( x + 1 ) ( y   + y - 5 )
  is represented as a list of three pairs , the first pair being the
  the number 12 and the exponent 1 , the second being the
  polynomial ( x + 1 ) and the exponent 2 , and the third being the
  polynomial ( y^2 + y - 5 ) and the exponent 4 .
  Currently , the first pair is always   an integer with exponent 1 .

  for some operations , the idea of gcd as used here assumes coefs form a
  to have integer coefficients .   ( Rationals are not a UFD since
   1/2 * 2 = 1/3 * 3 = 1 etc . )

  There is a first pair that always looks like ( integer . 1 ) .
  exponents which are integers > = 1 .


  The polynomials in an fpe , except for intermediate forms produced
  increasing order . Note that two polynomials in fpe form

  For example ,   ( ( 1 . 1 ) ( ( x^2 + 2x+1 ) . 1 ) )   and ( ( 1 . 1 ) ( ( x+1 ) . 2 ) )
   ( ( 1 . 1 ) ( ( 2x^2 + 4*x+2 ) . 1 ) )

   ( ( 2 . 1 ) ( ( x+1 ) . 2 ) )

  On the other hand , if we were adding ( ( 1 . 1 ) ( ( x+1 ) . 2 ) ) and
  ( ( 1 . 1 ) ( ( x+1 ) . 3 ) ) , the common factor of x+1 would be obvious
  and we would get   ( ( 1 . 1 ) ( ( x+1 ) . 2 ) ( ( x+2 ) . 1 ) )

(proclaim '(optimize (speed 3) (safety 0) (space 0) (compilation-speed 0)))

(in-package :mma)
#+ignore (export '(make-fpe fpe-expand fpe-coef-p fpe-coefone-p fpe-negativep
		   fpe-insert fpe-coefzero-p
		   make-rat rat-p rat make-good-rat rat-numerator 
		   rat-denominator  rat+ rat* rat^ rat/ poly2rat
		   *rat-gcd-level*))


(defvar *rat-gcd-level* 'gcd) 
(defvar *expand* nil)

   cdi : just find the obvious common factors by identity . That is

  cdd : common divisor by division .   That is , u and v have a
  common divisor if one divides the other exactly .

  current implementation uses gcd as default .


  NOT IMPLEMENTED techniques :


  considered for pair - wise gcds .   That is , ( x^2 - 1)*(x+1 ) -- > ( x+1)^2*(x-1 )



  make - fpe : make a " normal " fpe p^e from a polynomial p , exponent e>=0

   ( make - fpe 0 1 )     -- >   ( ( 0    . 1 ) )
    by virtue of using fpe - insert with monomial flag = t , if
    produce ( ( 9 . 1 ) ( # ( x 0 1 ) . 4 ) ( # ( y 0 1 ) . 6 ) ) .

(defun make-fpe (p e) 
  (cond ((coefp p) 
	 (list (cons (coef^ p e) 1)))
	(t (fpe-insert p e '((1 . 1)) t))))

(defun poly2rat (p e)
  (cond ((>= e 0)(make-rat :numerator (make-fpe p e)
			  :denominator (list '(1 . 1))))
	(t (make-good-rat (list '(1 . 1))
		      (make-fpe p (- e))))))
	
  fpe - copy : return a copy of the fpe u

(defun fpe-copy (u) (copy-tree u))


(defun fpe-expand (u)
  (let ((*expand* t))
  (do ((ul u (cdr ul)) (p 1))
      ((null ul) p)
      (setq p (p* p (p^ (caar ul) (cdar ul)))))))

  A fpe is " normal " if its first factor is an integer and it has
  monomial factors should be 1 . ( reason : they look odd otherwise ...
  consider 3*x^2 * 4*y^2   instead of 12*x^2*y^2 )
  are represented as c^1 , including coefzero and coefone .

  fpe - norm : normalize a fpe


(defun fpe-norm (u &aux ctest (const 1))
  (labels
	(cond ((null ul) nil)
  stick coefs up front
	       (setq const (coef* const
				  (coef^ (caar ul)(cdar ul))))
	       (fn1 (cdr ul)))
	      (t 
	       (fn1-insert (car ul) (fn1 (cdr ul))))))
    (fn1-insert (h ul)
		(cond ((null ul)
  this fpe , make the pair of this factor and its power
		       (cons h ul))
		      ((eq 'e (setq ctest (pcompare (car h) (caar ul))))
		       (cons (cons (car h)(+ (cdr h)(cdar ul))) (cdr ul)))
		      ((eq 'l ctest)
  this poly is less in ordering that (   ) ,
		       (cons h ul))
		      (t 
		       (cons  (car ul)(fn1-insert h (cdr ul)))
			 ul)))
  end of local fns
			  
   (setq u (fn1 u))
   (cons  (cons const 1) u)))
			 
  
(defun fpe-coef-p (u) (null(cdr u)))

(defun fpe-coefone-p (u) (and (null (cdr u)) (coefonep (caar u))))

(defun fpe-coefzero-p (u) (coefzerop (caar u)))



  fpe - insert : insert a pair of polynomial and exponent into the fpe u.
              adds the pair to the fpe . This will not change 

  in effect , this multiplies u in fpe form , by poly^u .
  compare to fpe- * , which multiplies 2 polys in fpe form .
  if monom is t , then maybe poly is a monomial . Check it

(defun fpe-insert (poly exp u &optional (monom nil) &aux ctest)
  (labels
   ((fi1(ul)
	     (cond ((null ul)
  this fpe , insert the new factor and its power
		    (list (cons poly exp)))
		   
		   ((eq 'e (setq ctest (pcompare poly (caar ul))))
		    (cons (cons (caar ul) (+ exp (cdar ul)))(cdr ul)))
		   ((eq 'l ctest)
  this poly is less in ordering than (   ) ,
		    (cons (cons poly exp) ul))
		   (t (cons (car ul)(fi1 (cdr ul)))))))
	(cond ((coefp poly)
  inserting an integer factor : modify the first element
	       (cond ((coefonep poly) u)
		     ((coefzerop poly) (make-fpe 0 1))
		     (t (cons (cons (coef* (caar u)
						 (coef^ poly exp))
					  1)
			      (cdr u)))))
multiplying u by z^0 gives u
	      ((and monom (monomialp poly)(or (> (degree poly) 1)
					      (not (coefonep (lc poly)))))
  the poly is something like ( 3*x^2 ) ^4
  insert 81   ( i.e. 3 ^ 4 ) into   1*(x)^8 .
  should work recursively for ( 3*y^2*x^2)^4
	       
	       (fpe-insert (lc poly) exp (fpe-insert 
					  (vector (svref poly 0) 0 1)
					  (* (degree poly) exp) u )
			   monom))
	      

	      ((and monom 
		    (not *expand*)
		    (coefzerop(constc poly))
		    (> (length (the simple-vector poly)) 3)
		    (coefzerop (svref poly 2)))
	       (do ((i 2 (1+ i)))
		    (fpe-insert
		     (vector (svref poly 0) 0 1)
		     (fpe-insert (polyshift poly (1- i)) exp u)))
		   (declare (fixnum i))
		   ))
		    		    
	      (t (fi1 u)))))

  (let ((z(subseq p i)))
    (setf (svref z 0) (svref p 0)) z))


  fpe- *   multiplies two fpe 's u and v. To be canonical ,

(defun fpe-* (u v &aux ctest)
  (labels ((fm1(ul vl)
	    (cond ((null ul) vl)
		  ((null vl) ul)
		  ((eq 'e (setq ctest (pcompare (caar ul) (caar vl))))
		   (cons (cons (caar ul) (+ (cdar ul)(cdar vl)))
			 (fm1 (cdr ul)(cdr vl))))
		   ((eq 'l ctest)
  (   ) is less in ordering that (   ) ,
		    (cons (car ul) (fm1 (cdr ul) vl)))
		   (t (cons (car vl)(fm1 ul (cdr vl)))))))

  u is coefone
  v is coefone
	      
		       (fm1 (cdr u)(cdr v)))))))

powering an fpe is " easy "
  (labels((fp^1 (a)
	      (cond ((null a)nil)
		    (t (cons (cons (caar a) (* n (cdar a)))
			     (fp^1 (cdr a)))))))
       (cons (cons (expt (caar a) n) 1)
	     (fp^1 (cdr a)))))

(defun fpe-negative-p (u) (coef-negative-p (caar u)))

(defun fpe-negate (u)
  (cons (cons (coefneg (caar u))(cdar u))(cdr u)))

  These procedures are concerned with finding the gcd and cd of fpe 's .
  They are independent of the implementation of  's .

  fpe - gcd - cofac :

  all values are fpe in normal form .

(defun fpe-gcd-cofac (u v)
   g= gcd = 1 , initially

    (do ((i u (cdr i))) 
  < gcd , u / gcd , v / gcd > .
	
	((null i) (values g (fpe-norm u)(fpe-norm v)))
  do this only if ( caar i ) is not 1

	    ((null j))
		   (ue (cdar i))
		   (ve (cdar j)))
	      (multiple-value-bind 
gcd , up / gp , vp / gp
	       (pgcdswitch-cofac up vp)
		     (t
  in g because of some earlier discovered gcd .
			
  blot out the two factors in u and v.
			
  ( y^2)^3 .   Could check in fpe - norm for monomials ,

			(cond ((= ge ue)
			       (setf (car i) (cons uq ue))
			       (setq up uq))
			       (setf (car i)(cons gp (- ue ge)))
			       (setq up gp)
			       (if (coefonep uq) nil
				 (nconc i (list (cons uq ue))))))
			(cond ((= ge ve)
			       (setf (car j) (cons vq ve))
			       (setq vp vq))
  that is , ge < ve
			       (setf (car j)(cons gp (- ve ge)))
			       (setq vp gp)
			       (if (coefonep vq) nil
				 (nconc j (list (cons vq ve))))))
			)))))))))))

(defun pgcdswitch-cofac(u v)
  (case *rat-gcd-level*
	(gcd 
compute the gcd . Which algorithm depends perhaps
	 (pgcd-cofac u v))
	(cdi
	 (if (equalp u v)(values u 1 1)(values 1 u v)))
	(cdd
	 (let (r)
	   (cond ((setq r (p/-test u v)) (values v r 1))
		 ((setq r (p/-test v u)) (values u 1 r))
		 (t (values 1 u v)))))
	   
  use gcd as default
	 (pgcd-cofac u v))))
	 

  consisting of two fpe 's : the numerator and the denominator .

(defstruct rat numerator denominator)



(defun make-good-rat (n d)
  (cond ((fpe-negative-p d)
	 (setq d (fpe-negate d))
	 (setq n (fpe-negate n))))
  (make-rat :numerator n :denominator d))

        2 ) each polynomial in a fpe that makes up the numerator
        3 ) the fpe 's are themselves normal : leading coefficient followed

  rat * : non - destructively multiply two rats r1 and r2
        reference -      paper


  (let ((a (rat-numerator r1))    
	(b (rat-denominator r1))  
	(c (rat-numerator r2))    
	(d (rat-denominator r2))  
	g num den)
  we assume that gcd(a , b)=1 , and also that gcd(c , d)=1 .
    (if (or (fpe-coefzero-p a)
	    (fpe-coefzero-p c))
	    (return-from rat* (poly2rat 0 1)))
  First , if a and c have any factors in common ,
  is predictably 1 , so no need to compute it .  
  as an additional optimization . (  -- this is what you were
    
    (multiple-value-setq (num  a c) (fpe-gcd-cofac a c))
    (multiple-value-setq (den  b d) (fpe-gcd-cofac b d))
  we can conclude that gcd(num , den)=1 by construction , and
  subject , however , to more gcd removal . Next ,
    (multiple-value-setq (g a d)(fpe-gcd-cofac a d))
    (multiple-value-setq (g b c)(fpe-gcd-cofac b c))
    (setq a (fpe-* a (fpe-^ num 2)))
    (setq b (fpe-* b (fpe-^ den 2)))
    (make-good-rat  (fpe-* a c) (fpe-* b d))))

  (let ((a (rat-numerator r1))    
	(b (rat-denominator r1))  
	g num den)
    (if (fpe-coefzero-p a) (return-from rat/ (poly2rat 0 1)))
    (if (fpe-coefzero-p d) (error "Rational division by zero"))

    (multiple-value-setq (num  a c) (fpe-gcd-cofac a c))
    (multiple-value-setq (den  b d) (fpe-gcd-cofac b d))

    (multiple-value-setq (g a d)(fpe-gcd-cofac a d))
    (multiple-value-setq (g b c)(fpe-gcd-cofac b c))

    (setq a (fpe-* a (fpe-^ num 2)))
    (setq b (fpe-* b (fpe-^ den 2)))

    (make-good-rat  (fpe-* a c) (fpe-* b d))))


  rat+ : non - destructively add two rats r1 and r2

(defun rat+ (r1 r2)
  (let ((a (rat-numerator r1))
	(b (rat-denominator r1))
	(n (make-fpe (coefone) 1))
	num den	g)
    (if (fpe-coefzero-p a) (return-from rat+ r2))
    (if (fpe-coefzero-p c) (return-from rat+ r1))
    
  extract common factors from a and b , also b and d.
    (multiple-value-setq
     (num a c)(fpe-gcd-cofac a c))
    (multiple-value-setq
     (den b d)(fpe-gcd-cofac b d))
    
  now   ( num / den ) * ( a / b+c / d ) and gcd(b , d)=1 .
		(fpe-expand (fpe-* b c)))))
    (if (coefzerop n) (return-from rat+ (make-rat :numerator '((0 . 1))
						  :denominator '((1 . 1)))))
  now n is an fpe .

  set denom to bd * den
        
    (multiple-value-setq (g n den)(fpe-gcd-cofac n den))
    (make-good-rat num den)))


(defun rat^ (r e)
	 (cond ((= e 0) (make-rat :numerator '((1 . 1))
				  :denominator '((1 . 1))))
	       ((> e 0)
		(make-rat :numerator (fpe-^ (rat-numerator r) e)
			  :denominator (fpe-^ (rat-denominator r) e)))
	       (t (if (fpe-coefzero-p (rat-numerator r))
		      (error "Rational division by zero"))
check sign of  ..
		    (fpe-^ (rat-denominator r) (- e))
		    (fpe-^ (rat-numerator r) (- e))))))


