source: branches/f4grobner/grobner.lisp@ 69

;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;                                                                              
;;;  Copyright (C) 1999, 2002, 2009 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.  
;;;                                                                              
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(in-package :maxima)

(macsyma-module cgb-maxima)

(eval-when
    #+gcl (load eval)
    #-gcl (:load-toplevel :execute)
    (format t "~&Loading maxima-grobner ~a ~a~%"
            "$Revision: 2.0 $" "$Date: 2015/06/02 0:34:17 $"))

;;FUNCTS is loaded because it contains the definition of LCM
($load "functs")






;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Global switches
;; (Can be used in Maxima just fine)
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defmvar $poly_monomial_order '$lex
  "This switch controls which monomial order is used in polynomial
and Grobner basis calculations. If not set, LEX will be used")

(defmvar $poly_coefficient_ring '$expression_ring
  "This switch indicates the coefficient ring of the polynomials
that will be used in grobner calculations. If not set, Maxima's
general expression ring will be used. This variable may be set
to RING_OF_INTEGERS if desired.")

(defmvar $poly_primary_elimination_order nil
  "Name of the default order for eliminated variables in elimination-based functions.
If not set, LEX will be used.")

(defmvar $poly_secondary_elimination_order nil
  "Name of the default order for kept variables in elimination-based functions.
If not set, LEX will be used.")

(defmvar $poly_elimination_order nil
  "Name of the default elimination order used in elimination calculations.
If set, it overrides the settings in variables POLY_PRIMARY_ELIMINATION_ORDER
and SECONDARY_ELIMINATION_ORDER. The user must ensure that this is a true
elimination order valid for the number of eliminated variables.")

(defmvar $poly_return_term_list nil
  "If set to T, all functions in this package will return each polynomial as a
list of terms in the current monomial order rather than a Maxima general expression.")

(defmvar $poly_grobner_debug nil
  "If set to TRUE, produce debugging and tracing output.")

(defmvar $poly_grobner_algorithm '$buchberger
  "The name of the algorithm used to find grobner bases.")

(defmvar $poly_top_reduction_only nil
  "If not FALSE, use top reduction only whenever possible.
Top reduction means that division algorithm stops after the first reduction.")



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Coefficient ring operations
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; These are ALL operations that are performed on the coefficients by
;; the package, and thus the coefficient ring can be changed by merely
;; redefining these operations.
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defstruct (ring)
  (parse #'identity :type function)
  (unit #'identity :type function)
  (zerop #'identity :type function)
  (add #'identity :type function)
  (sub #'identity :type function)
  (uminus #'identity :type function)
  (mul #'identity :type function)
  (div #'identity :type function)
  (lcm #'identity :type function)
  (ezgcd #'identity :type function)
  (gcd #'identity :type function))

(defparameter *ring-of-integers*
    (make-ring
     :parse #'identity
     :unit #'(lambda () 1)
     :zerop #'zerop
     :add #'+
     :sub #'-
     :uminus #'-
     :mul #'*
     :div #'/
     :lcm #'lcm
     :ezgcd #'(lambda (x y &aux (c (gcd x y))) (values c (/ x c) (/ y c)))
     :gcd #'gcd)
  "The ring of integers.")



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; This is how we perform operations on coefficients
;; using Maxima functions. 
;;
;; Functions and macros dealing with internal representation structure
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;




;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Low-level polynomial arithmetic done on 
;; lists of terms
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defmacro termlist-lt (p) `(car ,p))
(defun termlist-lm (p) (term-monom (termlist-lt p)))
(defun termlist-lc (p) (term-coeff (termlist-lt p)))

(define-modify-macro scalar-mul (c) coeff-mul)

(defun scalar-times-termlist (ring c p)
  "Multiply scalar C by a polynomial P. This function works
even if there are divisors of 0."
  (mapcan
   #'(lambda (term)
       (let ((c1 (funcall (ring-mul ring) c (term-coeff term))))
         (unless (funcall (ring-zerop ring) c1)
           (list (make-term (term-monom term) c1)))))
   p))


(defun term-mul (ring term1 term2)
  "Returns (LIST TERM) wheter TERM is the product of the terms TERM1 TERM2,
or NIL when the product is 0. This definition takes care of divisors of 0
in the coefficient ring."
  (let ((c (funcall (ring-mul ring) (term-coeff term1) (term-coeff term2))))
    (unless (funcall (ring-zerop ring) c)
      (list (make-term (monom-mul (term-monom term1) (term-monom term2)) c)))))

(defun term-times-termlist (ring term f)
  (declare (type ring ring))
  (mapcan #'(lambda (term-f) (term-mul ring term term-f)) f))

(defun termlist-times-term (ring f term)
  (mapcan #'(lambda (term-f) (term-mul ring term-f term)) f))

(defun monom-times-term (m term)
  (make-term (monom-mul m (term-monom term)) (term-coeff term)))

(defun monom-times-termlist (m f)
  (cond
   ((null f) nil)
   (t
    (mapcar #'(lambda (x) (monom-times-term m x)) f))))

(defun termlist-uminus (ring f)
  (mapcar #'(lambda (x)
              (make-term (term-monom x) (funcall (ring-uminus ring) (term-coeff x))))
          f))

(defun termlist-add (ring p q)
  (declare (type list p q))
  (do (r)
      ((cond
        ((endp p)
         (setf r (revappend r q)) t)
        ((endp q)
         (setf r (revappend r p)) t)
        (t
         (multiple-value-bind
             (lm-greater lm-equal)
             (monomial-order (termlist-lm p) (termlist-lm q))
           (cond
            (lm-equal
             (let ((s (funcall (ring-add ring) (termlist-lc p) (termlist-lc q))))
               (unless (funcall (ring-zerop ring) s)    ;check for cancellation
                 (setf r (cons (make-term (termlist-lm p) s) r)))
               (setf p (cdr p) q (cdr q))))
            (lm-greater
             (setf r (cons (car p) r)
                   p (cdr p)))
            (t (setf r (cons (car q) r)
                     q (cdr q)))))
         nil))
       r)))

(defun termlist-sub (ring p q)
  (declare (type list p q))
  (do (r)
      ((cond
        ((endp p)
         (setf r (revappend r (termlist-uminus ring q)))
         t)
        ((endp q)
         (setf r (revappend r p))
         t)
        (t
         (multiple-value-bind
             (mgreater mequal)
             (monomial-order (termlist-lm p) (termlist-lm q))
           (cond
            (mequal
             (let ((s (funcall (ring-sub ring) (termlist-lc p) (termlist-lc q))))
               (unless (funcall (ring-zerop ring) s)    ;check for cancellation
                 (setf r (cons (make-term (termlist-lm p) s) r)))
               (setf p (cdr p) q (cdr q))))
            (mgreater
             (setf r (cons (car p) r)
                   p (cdr p)))
            (t (setf r (cons (make-term (termlist-lm q) (funcall (ring-uminus ring) (termlist-lc q))) r)
                     q (cdr q)))))
         nil))
       r)))

;; Multiplication of polynomials
;; Non-destructive version
(defun termlist-mul (ring p q)
  (cond ((or (endp p) (endp q)) nil)    ;p or q is 0 (represented by NIL)
        ;; If p=p0+p1 and q=q0+q1 then pq=p0q0+p0q1+p1q
        ((endp (cdr p))
         (term-times-termlist ring (car p) q))
        ((endp (cdr q))
         (termlist-times-term ring p (car q)))
        (t
         (let ((head (term-mul ring (termlist-lt p) (termlist-lt q)))
               (tail (termlist-add ring (term-times-termlist ring (car p) (cdr q))
                                   (termlist-mul ring (cdr p) q))))
           (cond ((null head) tail)
                 ((null tail) head)
                 (t (nconc head tail)))))))
                    
(defun termlist-unit (ring dimension)
  (declare (fixnum dimension))
  (list (make-term (make-monom dimension :initial-element 0)
                   (funcall (ring-unit ring)))))

(defun termlist-expt (ring poly n &aux (dim (monom-dimension (termlist-lm poly))))
  (declare (type fixnum n dim))
  (cond
   ((minusp n) (error "termlist-expt: Negative exponent."))
   ((endp poly) (if (zerop n) (termlist-unit ring dim) nil))
   (t
    (do ((k 1 (ash k 1))
         (q poly (termlist-mul ring q q))       ;keep squaring
         (p (termlist-unit ring dim) (if (not (zerop (logand k n))) (termlist-mul ring p q) p)))
        ((> k n) p)
      (declare (fixnum k))))))






;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Debugging/tracing
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defmacro debug-cgb (&rest args)
  `(when $poly_grobner_debug (format *terminal-io* ,@args)))





;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; These are provided mostly for debugging purposes To enable
;; verification of grobner bases with BUCHBERGER-CRITERION, do
;; (pushnew :grobner-check *features*) and compile/load this file.
;; With this feature, the calculations will slow down CONSIDERABLY.
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defun grobner-test (ring g f)
  "Test whether G is a Grobner basis and F is contained in G. Return T
upon success and NIL otherwise."
  (debug-cgb "~&GROBNER CHECK: ")
  (let (($poly_grobner_debug nil)
        (stat1 (buchberger-criterion ring g))
        (stat2
          (every #'poly-zerop
                 (makelist (normal-form ring (copy-tree (elt f i)) g nil)
                           (i 0 (1- (length f)))))))
    (unless stat1 (error "~&Buchberger criterion failed."))
    (unless stat2
      (error "~&Original polys not in ideal spanned by Grobner.")))
  (debug-cgb "~&GROBNER CHECK END")
  t)


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Selection of algorithm and pair heuristic
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defun find-grobner-function (algorithm)
  "Return a function which calculates Grobner basis, based on its
names. Names currently used are either Lisp symbols, Maxima symbols or
keywords."
  (ecase algorithm
    ((buchberger :buchberger $buchberger) #'buchberger)
    ((parallel-buchberger :parallel-buchberger $parallel_buchberger) #'parallel-buchberger)
    ((gebauer-moeller :gebauer_moeller $gebauer_moeller) #'gebauer-moeller)))

(defun grobner (ring f &optional (start 0) (top-reduction-only nil))
  ;;(setf F (sort F #'< :key #'sugar))
  (funcall
   (find-grobner-function $poly_grobner_algorithm)
   ring f start top-reduction-only))

(defun reduced-grobner (ring f &optional (start 0) (top-reduction-only $poly_top_reduction_only))
  (reduction ring (grobner ring f start top-reduction-only)))

(defun set-pair-heuristic (method)
  "Sets up variables *PAIR-KEY-FUNCTION* and *PAIR-ORDER* used
to determine the priority of critical pairs in the priority queue."
  (ecase method
    ((sugar :sugar $sugar)
     (setf *pair-key-function* #'sugar-pair-key
           *pair-order* #'sugar-order))
;     ((minimal-mock-spoly :minimal-mock-spoly $minimal_mock_spoly)
;      (setf *pair-key-function* #'mock-spoly
;          *pair-order* #'mock-spoly-order))
    ((minimal-lcm :minimal-lcm $minimal_lcm)
     (setf *pair-key-function* #'(lambda (p q)
                                   (monom-lcm (poly-lm p) (poly-lm q)))
           *pair-order* #'reverse-monomial-order))
    ((minimal-total-degree :minimal-total-degree $minimal_total_degree)
     (setf *pair-key-function* #'(lambda (p q)
                                   (monom-total-degree
                                    (monom-lcm (poly-lm p) (poly-lm q))))
           *pair-order* #'<))
    ((minimal-length :minimal-length $minimal_length)
     (setf *pair-key-function* #'(lambda (p q)
                                   (+ (poly-length p) (poly-length q)))
           *pair-order* #'<))))




;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Set up the coefficients to be polynomials
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; (defun poly-ring (ring vars)
;;   (make-ring 
;;    :parse #'(lambda (expr) (poly-eval ring expr vars))
;;    :unit #'(lambda () (poly-unit ring (length vars)))
;;    :zerop #'poly-zerop
;;    :add #'(lambda (x y) (poly-add ring x y))
;;    :sub #'(lambda (x y) (poly-sub ring x y))
;;    :uminus #'(lambda (x) (poly-uminus ring x))
;;    :mul #'(lambda (x y) (poly-mul ring x y))
;;    :div #'(lambda (x y) (poly-exact-divide ring x y))
;;    :lcm #'(lambda (x y) (poly-lcm ring x y))
;;    :ezgcd #'(lambda (x y &aux (gcd (poly-gcd ring x y)))
;;            (values gcd
;;                    (poly-exact-divide ring x gcd)
;;                    (poly-exact-divide ring y gcd)))
;;    :gcd #'(lambda (x y) (poly-gcd x y))))



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Conversion from internal to infix form
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defun coerce-to-infix (poly-type object vars)
  (case poly-type
    (:termlist
     `(+ ,@(mapcar #'(lambda (term) (coerce-to-infix :term term vars)) object)))
    (:polynomial
     (coerce-to-infix :termlist (poly-termlist object) vars))
    (:poly-list
     `([ ,@(mapcar #'(lambda (p) (coerce-to-infix :polynomial p vars)) object)))
    (:term
     `(* ,(term-coeff object)
         ,@(mapcar #'(lambda (var power) `(expt ,var ,power))
                   vars (monom-exponents (term-monom object)))))
    (otherwise
     object)))



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Maxima expression ring
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defparameter *expression-ring*
    (make-ring 
     ;;(defun coeff-zerop (expr) (meval1 `(($is) (($equal) ,expr 0))))
     :parse #'(lambda (expr)
                (when modulus (setf expr ($rat expr)))
                expr)
     :unit #'(lambda () (if modulus ($rat 1) 1))
     :zerop #'(lambda (expr)
                ;;When is exactly a maxima expression equal to 0?
                (cond ((numberp expr)
                       (= expr 0))
                      ((atom expr) nil)
                      (t
                       (case (caar expr)
                         (mrat (eql ($ratdisrep expr) 0))
                         (otherwise (eql ($totaldisrep expr) 0))))))
     :add #'(lambda (x y) (m+ x y))
     :sub #'(lambda (x y) (m- x y))
     :uminus #'(lambda (x) (m- x))
     :mul #'(lambda (x y) (m* x y))
     ;;(defun coeff-div (x y) (cadr ($divide x y)))
     :div #'(lambda (x y) (m// x y))
     :lcm #'(lambda (x y) (meval1 `((|$LCM|) ,x ,y)))
     :ezgcd #'(lambda (x y) (apply #'values (cdr ($ezgcd ($totaldisrep x) ($totaldisrep y)))))
     ;; :gcd #'(lambda (x y) (second ($ezgcd x y)))))
     :gcd #'(lambda (x y) ($gcd x y))))

(defvar *maxima-ring* *expression-ring*
  "The ring of coefficients, over which all polynomials 
are assumed to be defined.")


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Order utilities
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun find-order (order)
  "This function returns the order function bases on its name."
  (cond
   ((null order) nil)
   ((symbolp order)
    (case order
      ((lex :lex $lex) #'lex>) 
      ((grlex :grlex $grlex) #'grlex>)
      ((grevlex :grevlex $grevlex) #'grevlex>)
      ((invlex :invlex $invlex) #'invlex>)
      ((elimination-order-1 :elimination-order-1 elimination_order_1) #'elimination-order-1)
      (otherwise
       (mtell "~%Warning: Order ~M not found. Using default.~%" order))))
   (t
    (mtell "~%Order specification ~M is not recognized. Using default.~%" order)
    nil)))

(defun find-ring (ring)
  "This function returns the ring structure bases on input symbol."
  (cond
   ((null ring) nil)
   ((symbolp ring)
    (case ring
      ((expression-ring :expression-ring $expression_ring) *expression-ring*) 
      ((ring-of-integers :ring-of-integers $ring_of_integers) *ring-of-integers*) 
      (otherwise
       (mtell "~%Warning: Ring ~M not found. Using default.~%" ring))))
   (t
    (mtell "~%Ring specification ~M is not recognized. Using default.~%" ring)
    nil)))

(defmacro with-monomial-order ((order) &body body)
  "Evaluate BODY with monomial order set to ORDER."
  `(let ((*monomial-order* (or (find-order ,order) *monomial-order*)))
     . ,body))

(defmacro with-coefficient-ring ((ring) &body body)
  "Evaluate BODY with coefficient ring set to RING."
  `(let ((*maxima-ring* (or (find-ring ,ring) *maxima-ring*)))
     . ,body))

(defmacro with-elimination-orders ((primary secondary elimination-order)
                                   &body body)
  "Evaluate BODY with primary and secondary elimination orders set to PRIMARY and SECONDARY."
  `(let ((*primary-elimination-order* (or (find-order ,primary)  *primary-elimination-order*))
         (*secondary-elimination-order* (or (find-order ,secondary) *secondary-elimination-order*))
         (*elimination-order* (or (find-order ,elimination-order) *elimination-order*)))
     . ,body))



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Conversion from internal form to Maxima general form
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defun maxima-head ()
  (if $poly_return_term_list
      '(mlist)
    '(mplus)))

(defun coerce-to-maxima (poly-type object vars)
  (case poly-type
    (:polynomial 
     `(,(maxima-head) ,@(mapcar #'(lambda (term) (coerce-to-maxima :term term vars)) (poly-termlist object))))
    (:poly-list
     `((mlist) ,@(mapcar #'(lambda (p) ($ratdisrep (coerce-to-maxima :polynomial p vars))) object)))
    (:term
     `((mtimes) ,($ratdisrep (term-coeff object))
       ,@(mapcar #'(lambda (var power) `((mexpt) ,var ,power))
                 vars (monom-exponents (term-monom object)))))
    ;; Assumes that Lisp and Maxima logicals coincide
    (:logical object)
    (otherwise
     object)))



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Macro facility for writing Maxima-level wrappers for
;; functions operating on internal representation
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defmacro with-parsed-polynomials (((maxima-vars &optional (maxima-new-vars nil new-vars-supplied-p))
                                    &key (polynomials nil)
                                         (poly-lists nil)
                                         (poly-list-lists nil)
                                         (value-type nil))
                                   &body body
                                   &aux (vars (gensym))
                                        (new-vars (gensym)))
  `(let ((,vars (coerce-maxima-list ,maxima-vars))
         ,@(when new-vars-supplied-p
             (list `(,new-vars (coerce-maxima-list ,maxima-new-vars)))))
     (coerce-to-maxima
      ,value-type
      (with-coefficient-ring ($poly_coefficient_ring)
        (with-monomial-order ($poly_monomial_order)
          (with-elimination-orders ($poly_primary_elimination_order
                                    $poly_secondary_elimination_order
                                    $poly_elimination_order)
            (let ,(let ((args nil))
                    (dolist (p polynomials args)
                      (setf args (cons `(,p (parse-poly ,p ,vars)) args)))
                    (dolist (p poly-lists args)
                      (setf args (cons `(,p (parse-poly-list ,p ,vars)) args)))
                    (dolist (p poly-list-lists args)
                      (setf args (cons `(,p (parse-poly-list-list ,p ,vars)) args))))
              . ,body))))
      ,(if new-vars-supplied-p
           `(append ,vars ,new-vars)
         vars))))

(defmacro define-unop (maxima-name fun-name
                       &optional (documentation nil documentation-supplied-p))
  "Define a MAXIMA-level unary operator MAXIMA-NAME corresponding to unary function FUN-NAME."
  `(defun ,maxima-name (p vars
                             &aux
                             (vars (coerce-maxima-list vars))
                             (p (parse-poly p vars)))
     ,@(when documentation-supplied-p (list documentation))
     (coerce-to-maxima :polynomial (,fun-name *maxima-ring* p) vars)))

(defmacro define-binop (maxima-name fun-name
                        &optional (documentation nil documentation-supplied-p))
  "Define a MAXIMA-level binary operator MAXIMA-NAME corresponding to binary function FUN-NAME."
  `(defmfun ,maxima-name (p q vars
                             &aux
                             (vars (coerce-maxima-list vars))
                             (p (parse-poly p vars))
                             (q (parse-poly q vars)))
     ,@(when documentation-supplied-p (list documentation))
     (coerce-to-maxima :polynomial (,fun-name *maxima-ring* p q) vars)))



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Maxima-level interface functions
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; Auxillary function for removing zero polynomial
(defun remzero (plist) (remove #'poly-zerop plist))

;;Simple operators

(define-binop $poly_add poly-add
  "Adds two polynomials P and Q")

(define-binop $poly_subtract poly-sub
  "Subtracts a polynomial Q from P.")

(define-binop $poly_multiply poly-mul
  "Returns the product of polynomials P and Q.")

(define-binop $poly_s_polynomial spoly
  "Returns the syzygy polynomial (S-polynomial) of two polynomials P and Q.")

(define-unop $poly_primitive_part poly-primitive-part
  "Returns the polynomial P divided by GCD of its coefficients.")

(define-unop $poly_normalize poly-normalize
  "Returns the polynomial P divided by the leading coefficient.")

;;Functions

(defmfun $poly_expand (p vars)
  "This function is equivalent to EXPAND(P) if P parses correctly to a polynomial.
If the representation is not compatible with a polynomial in variables VARS,
the result is an error."
  (with-parsed-polynomials ((vars) :polynomials (p)
                            :value-type :polynomial)
                           p))

(defmfun $poly_expt (p n vars)
  (with-parsed-polynomials ((vars) :polynomials (p) :value-type :polynomial)
    (poly-expt *maxima-ring* p n)))

(defmfun $poly_content (p vars)
  (with-parsed-polynomials ((vars) :polynomials (p))
    (poly-content *maxima-ring* p)))

(defmfun $poly_pseudo_divide (f fl vars
                            &aux (vars (coerce-maxima-list vars))
                                 (f (parse-poly f vars))
                                 (fl (parse-poly-list fl vars)))
  (multiple-value-bind (quot rem c division-count)
      (poly-pseudo-divide *maxima-ring* f fl)
    `((mlist)
      ,(coerce-to-maxima :poly-list quot vars)
      ,(coerce-to-maxima :polynomial rem vars)
      ,c
      ,division-count)))

(defmfun $poly_exact_divide (f g vars)
  (with-parsed-polynomials ((vars) :polynomials (f g) :value-type :polynomial)
    (poly-exact-divide *maxima-ring* f g)))

(defmfun $poly_normal_form (f fl vars)
  (with-parsed-polynomials ((vars) :polynomials (f)
                                   :poly-lists (fl)
                                   :value-type :polynomial)
    (normal-form *maxima-ring* f (remzero fl) nil)))

(defmfun $poly_buchberger_criterion (g vars)
  (with-parsed-polynomials ((vars) :poly-lists (g) :value-type :logical)
    (buchberger-criterion *maxima-ring* g)))

(defmfun $poly_buchberger (fl vars)
  (with-parsed-polynomials ((vars) :poly-lists (fl) :value-type :poly-list)
    (buchberger *maxima-ring*  (remzero fl) 0 nil)))

(defmfun $poly_reduction (plist vars)
  (with-parsed-polynomials ((vars) :poly-lists (plist)
                                   :value-type :poly-list)
    (reduction *maxima-ring* plist)))

(defmfun $poly_minimization (plist vars)
  (with-parsed-polynomials ((vars) :poly-lists (plist)
                                   :value-type :poly-list)
    (minimization plist)))

(defmfun $poly_normalize_list (plist vars)
  (with-parsed-polynomials ((vars) :poly-lists (plist)
                                   :value-type :poly-list)
    (poly-normalize-list *maxima-ring* plist)))

(defmfun $poly_grobner (f vars)
  (with-parsed-polynomials ((vars) :poly-lists (f)
                                   :value-type :poly-list)
    (grobner *maxima-ring* (remzero f))))

(defmfun $poly_reduced_grobner (f vars)
  (with-parsed-polynomials ((vars) :poly-lists (f)
                                   :value-type :poly-list)
    (reduced-grobner *maxima-ring* (remzero f))))

(defmfun $poly_depends_p (p var mvars
                        &aux (vars (coerce-maxima-list mvars))
                             (pos (position var vars)))
  (if (null pos)
      (merror "~%Variable ~M not in the list of variables ~M." var mvars)
    (poly-depends-p (parse-poly p vars) pos)))

(defmfun $poly_elimination_ideal (flist k vars)
  (with-parsed-polynomials ((vars) :poly-lists (flist)
                                   :value-type :poly-list)
    (elimination-ideal *maxima-ring* flist k nil 0)))

(defmfun $poly_colon_ideal (f g vars)
  (with-parsed-polynomials ((vars) :poly-lists (f g) :value-type :poly-list)
    (colon-ideal *maxima-ring* f g nil)))

(defmfun $poly_ideal_intersection (f g vars)
  (with-parsed-polynomials ((vars) :poly-lists (f g) :value-type :poly-list)  
    (ideal-intersection *maxima-ring* f g nil)))

(defmfun $poly_lcm (f g vars)
  (with-parsed-polynomials ((vars) :polynomials (f g) :value-type :polynomial)
    (poly-lcm *maxima-ring* f g)))

(defmfun $poly_gcd (f g vars)
  ($first ($divide (m* f g) ($poly_lcm f g vars))))

(defmfun $poly_grobner_equal (g1 g2 vars)
  (with-parsed-polynomials ((vars) :poly-lists (g1 g2))
    (grobner-equal *maxima-ring* g1 g2)))

(defmfun $poly_grobner_subsetp (g1 g2 vars)
  (with-parsed-polynomials ((vars) :poly-lists (g1 g2))
    (grobner-subsetp *maxima-ring* g1 g2)))

(defmfun $poly_grobner_member (p g vars)
  (with-parsed-polynomials ((vars) :polynomials (p) :poly-lists (g))
    (grobner-member *maxima-ring* p g)))

(defmfun $poly_ideal_saturation1 (f p vars)
  (with-parsed-polynomials ((vars) :poly-lists (f) :polynomials (p)
                                   :value-type :poly-list)
    (ideal-saturation-1 *maxima-ring* f p 0)))

(defmfun $poly_saturation_extension (f plist vars new-vars)
  (with-parsed-polynomials ((vars new-vars)
                            :poly-lists (f plist)
                            :value-type :poly-list)
    (saturation-extension *maxima-ring* f plist)))

(defmfun $poly_polysaturation_extension (f plist vars new-vars)
  (with-parsed-polynomials ((vars new-vars)
                            :poly-lists (f plist)
                            :value-type :poly-list)
    (polysaturation-extension *maxima-ring* f plist)))

(defmfun $poly_ideal_polysaturation1 (f plist vars)
  (with-parsed-polynomials ((vars) :poly-lists (f plist)
                                   :value-type :poly-list)
    (ideal-polysaturation-1 *maxima-ring* f plist 0 nil)))

(defmfun $poly_ideal_saturation (f g vars)
  (with-parsed-polynomials ((vars) :poly-lists (f g)
                                   :value-type  :poly-list)
    (ideal-saturation *maxima-ring* f g 0 nil)))

(defmfun $poly_ideal_polysaturation (f ideal-list vars)
  (with-parsed-polynomials ((vars) :poly-lists (f)
                                   :poly-list-lists (ideal-list)
                                   :value-type :poly-list)
    (ideal-polysaturation *maxima-ring* f ideal-list 0 nil)))

(defmfun $poly_lt (f vars)
  (with-parsed-polynomials ((vars) :polynomials (f) :value-type :polynomial)
    (make-poly-from-termlist (list (poly-lt f)))))

(defmfun $poly_lm (f vars)
  (with-parsed-polynomials ((vars) :polynomials (f) :value-type :polynomial)
    (make-poly-from-termlist (list (make-term (poly-lm f) (funcall (ring-unit *maxima-ring*)))))))