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

;;
;; Written by: Tak W. Yan, Richard J. Fateman
;; File: rat.lisp
;; Contains: rational function (in partially-factored form) arithmetic 
(provide 'rat1)

;;; (c) Copyright 1990, 1991 Richard J. Fateman, Tak W. Yan

;;; 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)
	      ((coefp (caar ul)) ; 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 (caar ul),
		       ;; 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 u.

;; 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 (caar ul),
		    ;; 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)))))
	      ((= exp 0) 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)
		    ;; (caar ul) is less in ordering that (caar vl),
		    ;; 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)))))))

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

(defun fpe-^ (a n) ;;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 fpe'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
  (let((g '((1 . 1))))	;  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)))
	(if (not(coefonep (caar i)))  ;; 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 
	       (gp uq vq)		;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))
			      (t ; 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
	(t ; 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 - W. S.  Brown's 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. (Tak -- 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 r1+r2 = (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)))))
    (setq n (make-fpe n 1)) ; now n is an fpe.

    (setq den(fpe-* den (fpe-* b d))) ;; 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"))
		  (make-good-rat ;check sign of denom..
		    (fpe-^ (rat-denominator r) (- e))
		    (fpe-^ (rat-numerator r) (- e))))))