;;; -*- Mode:Common-Lisp; Package:ma; Base:10 -*-
;;
;; Written by: Richard J. Fateman
;; File: poly.lisp
;; Contains: vector-based polynomial recursive representation arithmetic

;;; (c) Copyright 1990 Richard J. Fateman, Peter Klier, Peter Norvig, Tak Yan

;;;(provide 'poly)
;;; A polynomial is either a number (integer -- for now)
;;; or a vector of length d+2 where d>0 is the degree of S
;;; S=#(mainvar coeff coeff ... coeff).
;;;              ^                ^
;;;              |                |-coeff of mainvar^d : non-zero
;;;              |-coeff of mainvar^0  
;;; e.g.
;;; 4 + 5*x + x^2  = #(mainvar 4 5 1), where mainvar is an integer
;;; representing x in the system.

;;; Coeffs can be polynomials in other variables.

;;; Note: this is not the shortest code, but it should not be much slower
;;; than the fastest code absolutely necessary to get the job done.

;;(provide 'vector-math)
(declaim (optimize (speed 3) (safety 0) (space 0) (compilation-speed 0)))
(declaim (inline svref = eql make-array))
(eval-when (compile load) (provide 'poly))
(in-package :ma)
;; coefp: is defined to return true (non-nil, anyway)
;;        for an element of the coefficient domain. By default, we
;;        use integers.  This could be changed in concert with other
;;        programs to support other domains (see coef+).
;;        Coefp is usually called to see if its argument is
;;        (recursively) a polynomial in yet more variables, or is the
;;        ``base case'' of a polynomial which is, in fact, a coefficient.
;;        By defining it as "numberp" rather than "integerp" it works more
;;        smoothly (no change needed for rationals or floats).

(defmacro coefp (x) `(numberp ,x))

;; mainvar:  extract the main variable from a polynomial

(defmacro mainvar (x) `(svref (the simple-vector ,x) 0))

;; samevar:  see if two polynomials have the same main variable

(defmacro samevar (x y) `(eq (mainvar ,x) (mainvar ,y)))

;; degree: extract the degree from a polynomial (degree in main variable)

(defmacro degree (x) `(- (length (the simple-vector ,x)) 2))

;; lc:  extract the leading coefficient of a polynomial

(defmacro lc (x) `(svref ,x (1- (length ,x))))

;; constc: extract the constant coefficient of a polynomial
;; with respect to its main variable

(defmacro constc(x) `(svref ,x 1))

;; monomialp: check if v is a monomial (like x, or 3x^2)

(defun monomialp (v)
  (declare (simple-vector v))
  (= (length v) (1+ (position-if-not #'(lambda(x) (equal x 0)) v :start 1))))

;; var>:  an ordering predicate on variables to determine which
;;        is more "main". We could use strings or symbols for the variables
;;        instead of numbers, and use alphabetic ordering for var>.
;;        This seems simpler overall, since the variables might be more
;;        complicated in applications (e.g. (cos z)).
;;        Some easy encoding of names to numbers can be used (e.g. the
;;        first symbol mentioned is 1, and then count up.)

(defmacro var> (x y) `(> ,x ,y))

;;  The next set of macro definitions supplies the arithmetic
;;  for coefficients.  In general, any domain that can support
;;  the coefficient operations below is fair game for the package.
;;  More advanced operations may place additional requirements
;;  on the coefficient domain (E.g. they form an algebraic Ring or Field).

;;  Add/subtract/multiply two coefficients, etc.

(defmacro coef+ (x y) `(+ ,x ,y))
(defmacro coef- (x y) `(- ,x ,y))
(defmacro coef* (x y) `(* ,x ,y))
(defmacro coef> (x y) `(> ,x ,y))

;; Dividing two coefficients is not so obvious.  If the coefficient
;; domain is the common-lisp rational numbers (an algebraic field) and
;; if y is non-zero, then (/ x y) is also in the same field.  If the
;; coefficient domain is the integers, then division may throw the
;; computation out of that domain.  (/ 1 2) is not an integer.
;; Here we assume the person or program using coef/ is knowledgeable
;; about its uses.  Actually, we go further than that here. We
;; don't use coef/ in this file at all.  In the places we use
;; division (in p/ ) we make the (perhaps unwarranted) assumption
;; that exact division over the integers is required. See p/.

(defmacro coef/ (x y) `(/ ,x ,y))

;;  coef^ computer a coefficient to an integer power n

(defmacro coef^ (x n) `(expt ,x ,n))

;; some extra pieces just in case they're needed: remainder, negation, and
;; absolute value

(defmacro coefrem (x y) `(mod ,x ,y))
(defun coefneg (x) (- x))
(defun coefabs (x) (abs x))

;; Tests for zero and unity are important, as are constructors

(defmacro coefzero () 0)               ;;  zero in the coefficient domain
(defmacro coefone  () 1)               ;; unity in the coefficient domain
(defmacro coefzerop (x) `(eql 0 ,x))
(defmacro coefonep (x) `(eql 1 ,x))
(defmacro coef-negative-p (x) `(and (coefp ,x) (< ,x 0)))
(defun coefzerofunp (x) (eql 0 x))    ;;  to funcall, we need a function

;;; This product function preserves both arguments and returns the product
;;; of two polynomials.

;; p*: return the product of v1*v2

(defun p* (v1 v2)
  (declare (inline make-array))
  (cond ((coefp v1) (p*cv v1 v2)) ;; call function to multiply coef*poly
	((coefp v2) (p*cv v2 v1))
	((samevar v1 v2)
	 ;; two polynomials in the same main variable
	 ;; do  n X m multiplications
	 (let (ilim jlim index ival res)
	   (declare (fixnum ilim jlim index i j))
	   (setq ilim (1- (length v1)) jlim (1- (length v2)))
	   (setq res (make-array (+ ilim jlim) :initial-element 0))
	   (setf (svref res 0) (mainvar v1))
	   (do ((i 1 (1+ i)))
	       ((> i ilim) res)
	       (setq ival (svref v1 i))
	       (unless (coefzerop ival)	;; shortcut 0 * anything = 0
		       (do ((j 1 (1+ j)))
			   ((> j jlim))
			   (declare (fixnum i j))
			   (setq index (+ i j -1))
			   (setf (svref res index)
				 (p+vv-zero-check (svref res index)
						(p* ival (svref v2 j)))))))))
	((var> (mainvar v1) (mainvar v2)) (p*cv v2 v1))
	(t (p*cv v1 v2))))

;; p*cv: coefficient times polynomial vector;
;;       preserves both inputs and although the result can
;;       share substructure with the vector v, the top-level
;;       vector is new. (true recursively)

(defun p*cv (c v)
  (declare (inline make-array))
  (cond
   ((coefp v) (coef* c v))
   ((coefzerop c) (coefzero)) ;; 0 * anything is 0	
   ((coefonep c) v) ;; 1 * v = v. 
   ;; run down the length of the vector, multiplying.
   ;; p* is not destructive of its arguments either.
   (t (let* ((v v) (len (length (the simple-vector v))) u)
      (declare (simple-vector v u) (fixnum len) (inline svref setf mainvar))
      (setq u (make-array len))
      (setf (mainvar u) (mainvar v))
      (do
       ((i (1- len) (1- i)))
       ((= i 0) u)
       (declare (fixnum i))
       (setf (svref u i) (p* c (svref v i))))))))

;; pnegate: negate p
;; there are trivial other ways of doing this
;; like multiplying by -1. This takes less time
;; and space.

(defun pnegate(v)
  (declare (simple-vector v))
  (if (coefp v)
      (coefneg v)
      (let ((u (make-array (length v))))
        (declare (simple-vector u))
        (setf (svref u 0)(svref v 0))
        (do ((i (1-(length v)) (1- i)))
            ((= i 0) u)  ;; no need to call pnorm since -u is normalized
            (declare (fixnum i))
            (setf (svref u i)(pnegate (svref v i)))))))

;; p+: sum two polynomials as vectors, non-destructively

(defun p+ (v1 v2)
  (cond ((coefp v1) (p+cv v1 v2))
        ((coefp v2) (p+cv v2 v1)) ;; reverse args
        ((samevar v1 v2) ;; same main var
	 (let
	   ((lv1 (length (the simple-vector v1)))
	    (lv2 (length (the simple-vector v2)))
	    (v1 v1) (v2 v2)) ;; not redundant
	   (declare (simple-vector v1 v2) (fixnum lv1 lv2))
	   (cond ((> lv1 lv2)
		  (p+into v2 v1 lv2 lv1)) ;; v1 is longer
		 (t (p+into v1 v2 lv1 lv2)))))
	((var> (mainvar v1) (mainvar v2)) (p+cv v2 v1))
	(t (p+cv v1 v2))))

;; p-: compute difference of two polynomials as vectors, non-destructively
;; this could be done trivially by u+(-v) but this is
;; faster and uses less space.

(defun p- (v1 v2)
  (declare (simple-vector v1 v2))
  (cond ((coefp v1) (p+cv v1 (pnegate v2)))
        ((coefp v2) (p+cv (pnegate v2) v1))
        ((samevar v1 v2);; same main var
         ;;; aw, could do this in a pickier fashion...
         (p+ v1 (pnegate v2)))
        ((var> (mainvar v1) (mainvar v2))(p+cv (pnegate v2) v1))
        (t (p+cv v1 (pnegate v2)))))

;; p+cv: add coeff to vector, non-destructively

(defun p+cv(c v)
  (if (coefp v)
      (coef+ c v)
      (let ((v (copy-seq v)))
	(declare (simple-vector v))
	(setf (svref v 1) (p+ c (svref v 1)))
	v)))

;; p+into: add terms from one polynomial v1 into another

(defun p+into (v1 v2 shorter longer)
  (let (res)
    (declare (fixnum shorter longer)
	     (simple-vector v1 v2 res)
	     (inline p+vv-zero-check))
    (setq res (make-array longer))
    (setf (svref res 0) (mainvar v2))
    (do ((i 1 (1+ i)))
	((= i shorter))	
	(declare (fixnum i)) 
	(setf (svref res i) (p+vv-zero-check (svref v1 i) (svref v2 i))))
    (do ((i shorter (1+ i)))
	((= i longer) (pnorm res))
	(declare (fixnum i))
	(setf (svref res i) (svref v2 i)))))

;; p+vv-zero-check: helper for p+into

(defun p+vv-zero-check (place form) 
  (if (coefzerop place) form (p+ place form)))

;; pnorm converts a polynomial into a normal form in case it is
;; really zero or a constant or has trailing (high degree) zero
;; coeffs.  pnorm is destructive. pnorm is Not recursive except
;; in converting constant terms to manifest non-vectors.
;; Assume x is an arbitrary main-variable index:
;; #(x 5) -> 5.  #(x 5 4 3 0 0) -> #(x 5 4 3).  #(x 0 0) -> 0. 
;; #(x 0 1) -> #(x 0 1) [no change]

;; pnorm: return the normal form of x

(defun pnorm (x)
  (if (coefp x)
      x
      (let ((x x) pos)
	(declare (simple-vector x) (fixnum pos)
		 (inline position-if-not coefzerofunp delete-if))
	(setq pos (position-if-not #'coefzerofunp x :from-end t))
	(cond ((= pos 0)
	       (coefzero)) ;; nothing left but the header: zero polynomial
	      ((= pos 1) ;; just the constant itself
	       (pnorm (svref x 1))) ;; constant polynomial
	      ((= pos (1- (length x))) x)
	      (t #-kcl (delete-if #'coefzerofunp x :start pos)
		 ;fix bug in kcl delete-if
		 #+kcl (setq x  (delete-if #'coefzerofunp x :start pos))
		 )))))

;; p^v: this may seem like a dumb way to compute powers, but it's not
;; Repeated multiplication is generally faster than squaring.  In this
;; representation, binomial expansion, a good bet for sparse representation,
;; is only sometimes advantageous, and then not by very much, usually.

(defun p^ (x n)             ;; x^n -  n integer, x polynomial
  (cond ((integerp n)
	 (cond ((minusp n) (error "negative powers not allowed"))
	       ((zerop n) 1) ;; x^0 = 1 (even if x = 0)
	       ((eql n 1) x) ;; x^1 = x
	       ((coefp x) (expt x n))
	       (t (p* x (p^ x (1- n))))))
	(t (error "only integer powers allowed"))))

;; p/: divide p by s; only exact divisions are performed, and
;;            the result is q such that q*s = p; exact means that p, q,
;;            and s all have integer coefficients.

;; This does NOT work over random coefficient domains that you might
;; construct.

;; Ordinarily a single value is returned unless there is an
;; error condition, in which case the first value is nil and the second
;; is a string which describes the problem.

;; Call p/-test if you want a "semi-predicate."
;; Look at the first (or only value). If it is non-nil, you've got
;; the exact quotient.

;; Call p/-and-barf if you want to signal an error in case
;; an exact division is not possible.

(defun p/-and-barf (p s)
  (multiple-value-bind (v err)
		       (catch 'p/ (p/ p s))
		       (cond ((null v) (error err))
			     (t v))))

(defun p/-test(p s)
  (catch 'p/ (p/ p s)))

(defun p/ (p s)	
  (cond ((coefzerop s) (throw 'p/ (values nil "Division by zero")))
	((coefonep s) p)
        ((coefzerop p) p)
        ((coefp p)
         (cond ((coefp s)
                (cond ((eql 0 (mod p s))(/ p s))
                      (t (throw 'p/
				(values nil "Inexact integer division")))))
               (t (throw	
                   'p/
                   (values nil "Division by a polynomial of higher degree")))))
        ((or (coefp s) (var> (mainvar p) (mainvar s))) (p/vc p s))
        ((var> (mainvar s) (mainvar p))
         (throw 'p/
                (values nil "Division by a polynomial of higher degree")))
        (t (p/helper p s))))

(defun p/vc (p s)
  (let ((q (copy-seq p)))
    (do ((i (1- (length p)) (1- i)))
        ((= i 0) q)
        (setf (svref q i) (p/ (svref q i) s)))))

(defun p/helper (p s)
  (declare (simple-vector p s))
  (let ((l1 (length p)) (l2 (length s)))
    (cond ((> l2 l1)
           (throw 'p/
                  (values nil "Division by a polynomial of higher degree")))
          (t (let ((q (make-array (+ 2 l1 (- l2)) :initial-element 0))
		   (sneg (pnegate s))
		   (slc (lc s)))
	       (setf (mainvar q) (mainvar p))
	       (do* ((i l1 (length p)))
		    ((< i l2) (throw 'p/
				     (values nil "Quotient not exact")))
		    (setf (svref q (+ 1 i (- l2))) (p/ (lc p) slc))
		    (setq p (pnorm (p+ p (p* (subseq q 0 (+ 2 i (- l2)))
						 sneg))))
		    (if (coefzerop p) (return (pnorm q)))
		    (if (coefp p)
			(throw 'p/
			       (values nil "Quotient not exact")))))))))

;; prem: return the pseudo remainder when p is divided by s
;; (You didn't learn this in high school.) Ref. Knuth, Art of Comptr. Prog.
;; vol. 2.  This turns out to be important for GCD computations.

(defun prem (p s)
  (cond ((coefzerop s) (error "Division by zero"))
	((coefzerop p) p)
        ((coefp p)
         (cond ((coefp s) (coefrem p s))
               (t p)))
        ((or (coefp s) (var> (mainvar p) (mainvar s))) (coefzero))
        ((var> (mainvar s) (mainvar p)) p)
	(t (prem2 p s))))

;; prem2: return the pseudo remainder when p is divided by s; p and s
;;        must have the same variable. This is the usual situation.
(defun prem2 (p s)
  (if (> (length s) (length p))
      p
    (let* ((slc (lc s))
             (l (- (length p) (length s)))
             (temp (make-array (+ 2 l) :initial-element 0)))
        (setf (mainvar temp) (mainvar p))       
        (do ((k l m)
             (m))
            (nil)
          #-kcl(delete-if #' identity temp :start (+ 2 k))
          ;; fix bug in kcl delete-if, which is non-destructive
          #+kcl (setq temp (delete-if #' identity temp :start (+ 2 k)))
            (setf (lc temp) (lc p))
            (setq p (pnorm (p+ (p* p slc) (p*cv -1 (p* temp s)))))
            (cond ((or (coefp p)
                       (var> (mainvar s) (mainvar p))) 
                   (return (cond ((= k 0) p)                      ; PATCH 12/8/94
                                 (t (p* p (p^ slc k))))))         ; philliao
                  ((< (setq m (- (length p) (length s))) 0) 
                   (return (cond ((= k 0) p) 
                                 (t (p* p (p^ slc k)))))))  
            (if (> (1- k) m)
                (setq p (p* p (p^ slc (- (1- k) m)))))))))



;; Subresultant polynomial gcd.
;; See Knuth vol 2 2nd ed p 410. or TOMS Sept. 1978.
;; pgcd2sr: return the gcd of u and v;

#+ignore  ;; code taken from maxima. reprogrammed without loop, later
(defun pgcd2sr (p q)					
  ;; set up so q is of same or lower degree
  ;; p and q have the same main variable
  (cond ((> (length q) (length p)) (rotatef p q)))

    (let*((pcont (pcontentsr p))
	  (qcont (pcontentsr q))
	  (p (p/ p pcont))
	  (q (p/ q qcont))
	  (content (pgcdsr pcont qcont)))
      (loop for g = 1 then (lc p)
	  for h = 1 then (p/ (p^ g d) h^1-d)
	  for d = (- (length p) (length q))
	  for h^1-d = (if (eql h 1) 1 (p^ h (1- d)))
	  do (psetq p q
		    q (p/ (prem p q) (p* g (p* h h^1-d))))
	  if (or (coefp q)(var> (mainvar p)(mainvar q)))
	  return (if (coefzerop q) (p* content (ppartsr p))content))))



(defun pgcdsr (u v)
  (cond ((coefzerop u) v)
	((coefzerop v) u)
	((coefp u) (cond ((coefp v) (gcd u v))
			 (t (pcontentxsr v u))))
	((coefp v) (pcontentxsr u v))
	((var> (mainvar v) (mainvar u)) (pcontentxsr v u))
	((var> (mainvar u) (mainvar v)) (pcontentxsr u v))
	(t (pgcd2sr u v))))

;; pgcd2sr: return the gcd of p and q, which have the same main variable

;#+ignore
(defun pgcd2sr (p q)
  ;; set up so q is of same or lower degree
  (cond ((> (length q) (length p)) (rotatef p q)))
  ;; next, remove the content
  (let*((pcont (pcontentsr p))
	(qcont (pcontentsr q))
	(p (p/ p pcont))
	(q (p/ q qcont))
	(content (pgcdsr pcont qcont)))
    
    (do* ((g 1 (lc p)) 
	 (h 1 (p/ (p^ g d) hpow))
	 (d (- (length p)(length q)) (- (length p)(length q)))
	 (hpow 1 (if (eql h 1) 1 (p^ h (1- d)))))
	(nil)
	(psetq p q
	       q (p/ (prem p q)(p* g (p* h hpow))))
		
;	(format  t "~%p=~s~%q=~s~% d=~s, h=~s" p q d h)
	(if (coefzerop q) ;; poly remainder seq ended with zero
;;	      (return (p* content p)) ;bug in version prior to 11/94
	    (return (p* content (p/ p (pcontentsr p))))
	  )
	(if (or (coefp q) (var> (mainvar p) (mainvar q)))
	    ;; the remainder sequence ended with non-zero, just
	    ;;return the gcd of the contents.
	    (return content)))))

#+ignore  ;; we used this when the real pgcd2sr had a bug
(defun pgcd2sr (p q)(pgcd2prim p q))

;; this is a primitive prs, just like pgcd2sr but not as fast sometimes

(defun pgcd2prim (u v)
  (cond ((> (length v) (length u)) (rotatef u v)))
  (let ((ucont (pcontentsr u))
	(vcont (pcontentsr v)))
    (do ((c (p/ u ucont) d)
	 (d (p/ v vcont) (ppartsr (prem c d))))
	(nil)
	(if (coefzerop d)
	    (return (p* (pgcdsr ucont vcont) c)))
	(if (or (coefp d) (var> (mainvar c) (mainvar d)))
	    (return (pgcdsr ucont vcont))))))

;; pcontentsr: return the content of p (see Knuth, op. cit)
;; the content is the GCD of the coefficients of the top-level in the
;; polynomial p.

(defun pcontentsr (p)
  (cond ((coefzerop p) p)
        ((coefp p) (coefone))
	(t (let ((len (length p)))
	     (do ((i 2 (1+ i)) (g (svref p 1) (pgcdsr g (svref p i))))
		 ((or (= i len) (coefonep g)) g))))))


; pcontentxsr: return the gcd of u and the content of p
;; this is a handy time-saver

(defun pcontentxsr (p u)
  (cond 
   ;;((coefzerop p) (format t "~%pcontentxsr line 1 p=~s u=~s" p u)u)
   ;;((coefp p) (format t "~%pcontentxsr line 2 p=~s u=~s" p u)(gcd p u))
	(t (do ((i (1- (length p)) (1- i)) (g u (pgcdsr g (svref p i))))
	       ((or (coefonep g) (= i 0)) g)))))

; ppartsr: return the primitive part of p

(defun ppartsr (p)
  (cond ((or (coefzerop p) (coefp p)) p)
        (t (p/ p (pcontentsr p)))))

;; in some of the algorithms, the cofactors will be "free".  Here we're
;; just dividing the inputs by the gcd.

(defun pgcd-cofac (u v)(let ((g (pgcdsr u v)))
			    (values  g (p/ u g)(p/ v g))))

;; this program provides an ordering on polynomials, lexicographically.
;; it returns 'l (less than) 'e (equal)  or 'g (greater than).

;; this ordering is recursive by main variable.  There are other orderings,
;; like total degree ordering.

(defun pcompare (u v) 
  (cond ((eq u v) 'e) ;; quick check for equality
	;; if u and v are coefficients, then compare their numerical values.
	((coefp u)(cond ((coefp v)(cond((coef> v u) 'l)
				       ((coef> u v) 'g)
				       (t           'e)))
			;;u is coefficient, v is not, so  u << v
			(t 'l)))
	;; v is coefficient, u is not, so u >> v
	((coefp v) 'g)
	;; neither is a coefficient
	(t (let ((um (mainvar u))
		 (vm (mainvar v))
		 (ud (degree u))
		 (vd (degree v))
		 (ctest nil))
	     (cond ((var> um vm) 'g)  ; u has a more-main variable than v
		   ((var> vm um) 'l)  ; the reverse case
		   ((> ud vd)  'g) ; u is higher degree  u >> v
		   ((> vd ud)  'l) ; v is higher degree, so  v >> u
		   ;; Same main variable, same degree.
		   ;; Compare the coefficients, starting from the highest
		   ;; degree. The first one that differs determines
		   ;; which direction the comparison goes.
		   
		   (t (do ((i (1+ ud) (1- i)))
			  ((= i 0) 'e); no differences found.  u = v
			  (setq ctest (pcompare (svref u i)(svref v i)))
			  (if (eq ctest 'e) nil (return ctest))
			  )))))))


;; compute the derivative wrt v, the encoding of a variable.
;; This assumes that there are no other dependencies on kernels, etc.

(defun pderiv (p v)
  (labels
   ((pd1 (p &aux r)
	 (cond ((coefp p) 0)
	       ((var> v (svref p 0)) 0)
	       ((eql (svref p 0) v)
		(setq r (make-array (1- (length p))))
		(setf (svref r 0) v)
		(do ((i (1- (length p)) (1- i)))
		    ((= i 1) (pnorm r))
		    (setf (svref r (1- i))(p* (- i 1) (svref p i)))))
	       (t (setq r(copy-seq p))
		  (do ((i (1- (length r))(1- i)))
		      ((= i 0) (pnorm r))
		      (setf (svref r i)(pd1 (svref r i))))))))
   (pd1 p)))