;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;                                                                              
;;;  Copyright (C) 1999, 2002, 2009, 2015 Marek Rychlik <rychlik@u.arizona.edu>		 
;;;  		       								 
;;;  This program is free software; you can redistribute it and/or modify	 
;;;  it under the terms of the GNU General Public License as published by	 
;;;  the Free Software Foundation; either version 2 of the License, or		 
;;;  (at your option) any later version.					 
;;; 		       								 
;;;  This program is distributed in the hope that it will be useful,		 
;;;  but WITHOUT ANY WARRANTY; without even the implied warranty of		 
;;;  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the		 
;;;  GNU General Public License for more details.				 
;;; 		       								 
;;;  You should have received a copy of the GNU General Public License		 
;;;  along with this program; if not, write to the Free Software 		 
;;;  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.	 
;;;										 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


(defpackage "POLYNOMIAL"
  (:use :cl :ring :ring-and-order :monomial :order :term :termlist :infix)
  (:export "POLY" 
	   "POLY-TERMLIST" 
	   "POLY-SUGAR" 
	   "POLY-LT"
	   "MAKE-POLY-FROM-TERMLIST"
	   "MAKE-POLY-ZERO"
	   "MAKE-VARIABLE"
	   "POLY-UNIT"
	   "POLY-LM"
	   "POLY-SECOND-LM"
	   "POLY-SECOND-LT"
	   "POLY-LC"
	   "POLY-SECOND-LC"
	   "POLY-ZEROP"
	   "POLY-LENGTH"
	   "SCALAR-TIMES-POLY"
	   "SCALAR-TIMES-POLY-1"
	   "MONOM-TIMES-POLY"
	   "TERM-TIMES-POLY"
	   "POLY-ADD"
	   "POLY-SUB"
	   "POLY-UMINUS"
	   "POLY-MUL"
	   "POLY-EXPT"
	   "POLY-APPEND"
	   "POLY-NREVERSE"
	   "POLY-CONTRACT"
	   "POLY-EXTEND"
	   "POLY-ADD-VARIABLES"
	   "POLY-LIST-ADD-VARIABLES"
	   "POLY-STANDARD-EXTENSION"
	   "SATURATION-EXTENSION"
	   "POLYSATURATION-EXTENSION"
	   "SATURATION-EXTENSION-1"
	   "COERCE-COEFF"
	   "POLY-EVAL"
	   "POLY-EVAL-SCALAR"
	   "SPOLY"
	   "POLY-PRIMITIVE-PART"
	   "POLY-CONTENT"
	   "READ-INFIX-FORM"
	   "READ-POLY"
	   "STRING->POLY"
	   ))

(in-package :polynomial)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Polynomials
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defstruct (poly
	     ;;
	     ;; BOA constructor, by default constructs zero polynomial
	     (:constructor make-poly-from-termlist (termlist &optional (sugar (termlist-sugar termlist))))
	     (:constructor make-poly-zero (&aux (termlist nil) (sugar -1)))
	     ;; Constructor of polynomials representing a variable
	     (:constructor make-variable (ring nvars pos &optional (power 1)
					       &aux
					       (termlist (list
							  (make-term-variable ring nvars pos power)))
					       (sugar power)))
	     (:constructor poly-unit (ring dimension
					   &aux
					   (termlist (termlist-unit ring dimension))
					   (sugar 0))))
  (termlist nil :type list)
  (sugar -1 :type fixnum))

;; Leading term
(defmacro poly-lt (p) `(car (poly-termlist ,p)))

;; Second term
(defmacro poly-second-lt (p) `(cadar (poly-termlist ,p)))

;; Leading monomial
(defun poly-lm (p) (term-monom (poly-lt p)))

;; Second monomial
(defun poly-second-lm (p) (term-monom (poly-second-lt p)))

;; Leading coefficient
(defun poly-lc (p) (term-coeff (poly-lt p)))

;; Second coefficient
(defun poly-second-lc (p) (term-coeff (poly-second-lt p)))

;; Testing for a zero polynomial
(defun poly-zerop (p) (null (poly-termlist p)))

;; The number of terms
(defun poly-length (p) (length (poly-termlist p)))

(defun scalar-times-poly (ring c p)
  (declare (type ring ring) (poly p))
  (make-poly-from-termlist (scalar-times-termlist ring c (poly-termlist p)) (poly-sugar p)))

;; The scalar product omitting the head term
(defun scalar-times-poly-1 (ring c p)
  (declare (type ring ring) (poly p))
  (make-poly-from-termlist (scalar-times-termlist ring c (cdr (poly-termlist p))) (poly-sugar p)))

(defun monom-times-poly (m p)
  (declare (poly p))
  (make-poly-from-termlist 
   (monom-times-termlist m (poly-termlist p))
   (+ (poly-sugar p) (monom-sugar m))))

(defun term-times-poly (ring term p)
  (declare (type ring ring) (type term term) (type poly p))
  (make-poly-from-termlist 
   (term-times-termlist ring term (poly-termlist p))
   (+ (poly-sugar p) (term-sugar term))))

(defun poly-add (ring-and-order p q)
  (declare (type ring-and-order ring-and-order) (type poly p q))
  (make-poly-from-termlist 
   (termlist-add ring-and-order 
		 (poly-termlist p) 
		 (poly-termlist q))
   (max (poly-sugar p) (poly-sugar q))))

(defun poly-sub (ring-and-order p q)
  (declare (type ring-and-order ring-and-order) (type poly p q))
  (make-poly-from-termlist 
   (termlist-sub ring-and-order (poly-termlist p) (poly-termlist q))
   (max (poly-sugar p) (poly-sugar q))))

(defun poly-uminus (ring p)
  (declare (type ring ring) (type poly p))
  (make-poly-from-termlist 
   (termlist-uminus ring (poly-termlist p))
   (poly-sugar p)))

(defun poly-mul (ring-and-order p q)
  (declare (type ring-and-order ring-and-order) (type poly p q))
  (make-poly-from-termlist
   (termlist-mul ring-and-order (poly-termlist p) (poly-termlist q))
   (+ (poly-sugar p) (poly-sugar q))))

(defun poly-expt (ring-and-order p n)
  (declare (type ring-and-order ring-and-order) (type poly p))
  (make-poly-from-termlist (termlist-expt ring-and-order (poly-termlist p) n) (* n (poly-sugar p))))

(defun poly-append (&rest plist)
  (make-poly-from-termlist (apply #'append (mapcar #'poly-termlist plist))
			   (apply #'max (mapcar #'poly-sugar plist))))

(defun poly-nreverse (p)
  (declare (type poly p))
  (setf (poly-termlist p) (nreverse (poly-termlist p)))
  p)

(defun poly-contract (p &optional (k 1))
  (declare (type poly p))
  (make-poly-from-termlist (termlist-contract (poly-termlist p) k)
			   (poly-sugar p)))

(defun poly-extend (p &optional (m (make-monom :dimension 1)))
  (declare (type poly p))
  (make-poly-from-termlist
   (termlist-extend (poly-termlist p) m)
   (+ (poly-sugar p) (monom-sugar m))))

(defun poly-add-variables (p k)
  (declare (type poly p))
  (setf (poly-termlist p) (termlist-add-variables (poly-termlist p) k))
  p)

(defun poly-list-add-variables (plist k)
  (mapcar #'(lambda (p) (poly-add-variables p k)) plist))

(defun poly-standard-extension (plist &aux (k (length plist)))
  "Calculate [U1*P1,U2*P2,...,UK*PK], where PLIST=[P1,P2,...,PK]."
  (declare (list plist) (fixnum k))
  (labels ((incf-power (g i)
	     (dolist (x (poly-termlist g))
	       (incf (monom-elt (term-monom x) i)))
	     (incf (poly-sugar g))))
    (setf plist (poly-list-add-variables plist k))
    (dotimes (i k plist)
      (incf-power (nth i plist) i))))

(defun saturation-extension (ring f plist &aux (k (length plist)) (d (monom-dimension (poly-lm (car plist)))))
  "Calculate [F, U1*P1-1,U2*P2-1,...,UK*PK-1], where PLIST=[P1,P2,...,PK]."
  (setf f (poly-list-add-variables f k)
	plist (mapcar #'(lambda (x)
			  (setf (poly-termlist x) (nconc (poly-termlist x)
							 (list (make-term (make-monom :dimension d)
									  (funcall (ring-uminus ring) (funcall (ring-unit ring)))))))
			  x)
		      (poly-standard-extension plist)))
  (append f plist))


(defun polysaturation-extension (ring f plist &aux (k (length plist))
						(d (+ k (monom-dimension (poly-lm (car plist))))))
  "Calculate [F, U1*P1+U2*P2+...+UK*PK-1], where PLIST=[P1,P2,...,PK]."
  (setf f (poly-list-add-variables f k)
	plist (apply #'poly-append (poly-standard-extension plist))
	(cdr (last (poly-termlist plist))) (list (make-term (make-monom :dimension d)
							    (funcall (ring-uminus ring) (funcall (ring-unit ring))))))
  (append f (list plist)))

(defun saturation-extension-1 (ring f p) (polysaturation-extension ring f (list p)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Evaluation of polynomial (prefix) expressions
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defun coerce-coeff (ring expr vars)
  "Coerce an element of the coefficient ring to a constant polynomial."
  ;; Modular arithmetic handler by rat
  (make-poly-from-termlist (list (make-term (make-monom :dimension (length vars))
					    (funcall (ring-parse ring) expr)))
			   0))

(defun poly-eval (expr vars 
		  &optional 
		    (ring *ring-of-integers*) 
		    (order #'lex>)
		    (list-marker '[)
		  &aux 
		    (ring-and-order (make-ring-and-order :ring ring :order order)))
  (labels ((p-eval (arg) (poly-eval arg vars ring order))
	   (p-eval-scalar (arg) (poly-eval-scalar arg))
	   (p-eval-list (args) (mapcar #'p-eval args))
	   (p-add (x y) (poly-add ring-and-order x y)))
    (cond
      ((null expr) (error "Empty expression"))
      ((eql expr 0) (make-poly-zero))
      ((member expr vars :test #'equalp)
       (let ((pos (position expr vars :test #'equalp)))
	 (make-variable ring (length vars) pos)))
      ((atom expr)
       (coerce-coeff ring expr vars))
      ((eq (car expr) list-marker)
       (cons list-marker (p-eval-list (cdr expr))))
      (t
       (case (car expr)
	 (+ (reduce #'p-add (p-eval-list (cdr expr))))
	 (- (case (length expr)
	      (1 (make-poly-zero))
	      (2 (poly-uminus ring (p-eval (cadr expr))))
	      (3 (poly-sub ring-and-order (p-eval (cadr expr)) (p-eval (caddr expr))))
	      (otherwise (poly-sub ring-and-order (p-eval (cadr expr))
				   (reduce #'p-add (p-eval-list (cddr expr)))))))
	 (*
	  (if (endp (cddr expr))		;unary
	      (p-eval (cdr expr))
	      (reduce #'(lambda (p q) (poly-mul ring-and-order p q)) (p-eval-list (cdr expr)))))
	 (/ 
	  ;; A polynomial can be divided by a scalar
	  (cond 
	    ((endp (cddr expr))
	     ;; A special case (/ ?), the inverse
	     (coerce-coeff ring (apply (ring-div ring) (cdr expr)) vars))
	    (t
	     (let ((num (p-eval (cadr expr)))
		   (denom-inverse (apply (ring-div ring)
					 (cons (funcall (ring-unit ring)) 
					       (mapcar #'p-eval-scalar (cddr expr))))))
	       (scalar-times-poly ring denom-inverse num)))))
	 (expt
	  (cond
	    ((member (cadr expr) vars :test #'equalp)
	     ;;Special handling of (expt var pow)
	     (let ((pos (position (cadr expr) vars :test #'equalp)))
	       (make-variable ring (length vars) pos (caddr expr))))
	    ((not (and (integerp (caddr expr)) (plusp (caddr expr))))
	     ;; Negative power means division in coefficient ring
	     ;; Non-integer power means non-polynomial coefficient
	     (coerce-coeff ring expr vars))
	    (t (poly-expt ring-and-order (p-eval (cadr expr)) (caddr expr)))))
	 (otherwise
	  (coerce-coeff ring expr vars)))))))

(defun poly-eval-scalar (expr 
			 &optional
			   (ring *ring-of-integers*)
			 &aux 
			   (order #'lex>))
  "Evaluate a scalar expression EXPR in ring RING."
  (poly-lc (poly-eval expr nil ring order)))

(defun spoly (ring f g)
  "It yields the S-polynomial of polynomials F and G."
  (declare (type poly f g))
  (let* ((lcm (monom-lcm (poly-lm f) (poly-lm g)))
	  (mf (monom-div lcm (poly-lm f)))
	  (mg (monom-div lcm (poly-lm g))))
    (declare (type monom mf mg))
    (multiple-value-bind (c cf cg)
	(funcall (ring-ezgcd ring) (poly-lc f) (poly-lc g))
      (declare (ignore c))
      (poly-sub 
       ring
       (scalar-times-poly ring cg (monom-times-poly mf f))
       (scalar-times-poly ring cf (monom-times-poly mg g))))))


(defun poly-primitive-part (ring p)
  "Divide polynomial P with integer coefficients by gcd of its
coefficients and return the result."
  (declare (type poly p))
  (if (poly-zerop p)
      (values p 1)
    (let ((c (poly-content ring p)))
      (values (make-poly-from-termlist (mapcar
			  #'(lambda (x)
			      (make-term (term-monom x)
					 (funcall (ring-div ring) (term-coeff x) c)))
			  (poly-termlist p))
			 (poly-sugar p))
	       c))))

(defun poly-content (ring p)
  "Greatest common divisor of the coefficients of the polynomial P. Use the RING structure
to compute the greatest common divisor."
  (declare (type poly p))
  (reduce (ring-gcd ring) (mapcar #'term-coeff (rest (poly-termlist p))) :initial-value (poly-lc p)))

(defun read-infix-form (&key (stream t))
  "Parser of infix expressions with integer/rational coefficients
The parser will recognize two kinds of polynomial expressions:

- polynomials in fully expanded forms with coefficients
  written in front of symbolic expressions; constants can be optionally
  enclosed in (); for example, the infix form
  	X^2-Y^2+(-4/3)*U^2*W^3-5
  parses to
	(+ (- (EXPT X 2) (EXPT Y 2)) (* (- (/ 4 3)) (EXPT U 2) (EXPT W 3)) (- 5))

- lists of polynomials; for example
        [X-Y, X^2+3*Z]          
  parses to
	  (:[ (- X Y) (+ (EXPT X 2) (* 3 Z)))
  where the first symbol [ marks a list of polynomials.

-other infix expressions, for example
        [(X-Y)*(X+Y)/Z,(X+1)^2]
parses to:
	(:[ (/ (* (- X Y) (+ X Y)) Z) (EXPT (+ X 1) 2))
Currently this function is implemented using M. Kantrowitz's INFIX package."
  (read-from-string
   (concatenate 'string
     "#I(" 
     (with-output-to-string (s)
       (loop
	 (multiple-value-bind (line eof)
	     (read-line stream t)
	   (format s "~A" line)
	   (when eof (return)))))
     ")")))
	
(defun read-poly (vars &key (stream t))
  "Reads an expression in prefix form from a stream STREAM.
If the expression represents a polynomial or a list of polynomials in variables VARS then
the polynomial or list of polynomials is returned."
  (poly-eval (read-infix-form :stream stream) vars))

(defun string->poly (str vars)
  "Converts a string STR to a polynomial in variables VARS"
  (with-input-from-string (s str)
    (read-poly vars :stream s)))
  

(defun poly->alist (p)
  "Convert a polynomial P to an association list. Thus, the format of the
returned value is  ((MONOM[0] . COEFF[0]) (MONOM[1] . COEFF[1]) ...), where
MONOM[I] is a list of exponents in the monomial and COEFF[I] is the
corresponding coefficient in the ring."
  (mapcar term->cons (poly-termlist p)))

(defun string->alist (str vars)
  "Convert a string STR representing a polynomial or polynomial list to
an association list
  (poly->alist (str->poly str vars))