;;; -*-  Mode: Lisp -*- 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;                                                                              
;;;  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.	 
;;;										 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Load this file into Maxima to bootstrap the Grobner package. 
;; NOTE: This file does use symbols defined by Maxima, so it
;; will not work when loaded in Common Lisp.
;;
;; DETAILS: This file implements an interface between the Grobner
;; basis package NGROBNER, which is a pure Common Lisp package, and
;; Maxima. NGROBNER for efficiency uses its own representation of
;; polynomials. Thus, it is necessary to convert Maxima representation
;; to the internal representation and back. The facilities to do so
;; are implemented in this file.
;;
;; Also, since the NGROBNER package consists of many Lisp files, it is
;; necessary to load the files. It is possible and preferrable to use
;; ASDF for this purpose. The default is ASDF. It is also possible to
;; simply used LOAD and COMPILE-FILE to accomplish this task.
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(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")
#+sbcl(progn (require 'asdf) (load "ngrobner.asd")(asdf:load-system :ngrobner))

(use-package :ngrobner)


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Maxima expression ring
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; This is how we perform operations on coefficients
;; using Maxima functions. 
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defparameter +maxima-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))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Maxima expression parsing
;;
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Functions and macros dealing with internal representation
;; structure.
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defun equal-test-p (expr1 expr2)
  (alike1 expr1 expr2))

(defun coerce-maxima-list (expr)
  "Convert a Maxima list to Lisp list."
  (cond
   ((and (consp (car expr)) (eql (caar expr) 'mlist)) (cdr expr))
   (t expr)))

(defun free-of-vars (expr vars) (apply #'$freeof `(,@vars ,expr)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Order utilities
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defun find-ring-by-name (ring)
  "This function returns the ring structure bases on input symbol."
  (cond
   ((null ring) nil)
   ((symbolp ring)
    (case ring
      ((maxima-ring  :maxima-ring #:maxima-ring $expression_ring #:expression_ring)
       +maxima-ring+) 
      ((ring-of-integers :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)))

(defun find-order-by-name (order)
  "This function returns the order function bases on its name."
  (cond
   ((null order) nil)
   ((symbolp order)
    (case order
      ((lex :lex $lex #:lex)    
       #'lex>) 
      ((grlex :grlex $grlex #:grlex)
       #'grlex>)
      ((grevlex :grevlex $grevlex #:grevlex)
       #'grevlex>)
      ((invlex :invlex $invlex #:invlex) 
       #'invlex>)
      (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-and-order-by-name (&optional
				      (ring (find-ring-by-name $poly_coefficient_ring))
				      (order (find-order-by-name $poly_monomial_order))
				      (primary-elimination-order (find-order-by-name $poly_primary_elimination_order))
				      (secondary-elimination-order (find-order-by-name $poly_secondary_elimination_order))
				    &aux 
				      (ring-and-order (make-ring-and-order 
						       :ring  ring
						       :order order
						       :primary-elimination-order primary-elimination-order
						       :secondary-elimination-order secondary-elimination-order)))
  ring-and-order)

(defun maxima->poly (expr vars 
		     &optional 
		       (ring-and-order (find-ring-and-order-by-name))
		     &aux
		       (vars (coerce-maxima-list vars))
		       (ring (ro-ring ring-and-order)))
  "Convert a maxima polynomial expression EXPR in variables VARS to
internal form. This works by first converting the expression to Lisp,
and then evaluating the expression using polynomial arithmetic
implemented by the POLYNOMIAL package."
  (labels ((parse (arg) (maxima->poly arg vars ring-and-order))
	   (parse-list (args) (mapcar #'parse args)))
    (cond
     ((eql expr 0) (make-poly-zero))
     ((member expr vars :test #'equal-test-p)
      (let ((pos (position expr vars :test #'equal-test-p)))
	(make-poly-variable ring (length vars) pos)))
     ((free-of-vars expr vars)
      ;;This means that variable-free CRE and Poisson forms will be converted
      ;;to coefficients intact
      (coerce-coeff ring expr vars))
     (t
      (case (caar expr)
	(mplus (reduce #'(lambda (x y) (poly-add ring-and-order x y)) (parse-list (cdr expr))))
	(mminus (poly-uminus ring (parse (cadr expr))))
	(mtimes
	 (if (endp (cddr expr))		;unary
	     (parse (cdr expr))
	     (reduce #'(lambda (p q) (poly-mul ring-and-order p q)) (parse-list (cdr expr)))))
	(mexpt
	 (cond
	  ((member (cadr expr) vars :test #'equal-test-p)
	   ;;Special handling of (expt var pow)
	   (let ((pos (position (cadr expr) vars :test #'equal-test-p)))
	     (make-poly-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
	   (mtell "~%Warning: Expression ~%~M~%contains power which is not a positive integer. Parsing as coefficient.~%"
		  expr)
	   (coerce-coeff ring expr vars))
	  (t (poly-expt ring (parse (cadr expr)) (caddr expr)))))
	(mrat (parse ($ratdisrep expr)))
	(mpois (parse ($outofpois expr)))
	(otherwise
	 (coerce-coeff ring expr vars)))))))

(defun maxima->poly-list (expr vars 
			  &optional 
			    (ring-and-order (find-ring-and-order-by-name)))
  "Convert a Maxima representation of a list of polynomials to the internal form."
  (case (caar expr)
    (mlist (mapcar #'(lambda (p) 
		       (maxima->poly p vars ring-and-order))
		   (cdr expr)))
    (otherwise (merror "Expression ~M is not a list of polynomials in variables ~M."
		       expr vars))))

(defun maxima->poly-list-of-lists (poly-list-of-lists vars
				   &optional 
				     (ring-and-order (find-ring-and-order-by-name)))
  "Parse a Maxima representation of a list of lists of polynomials."
  (mapcar #'(lambda (g) (maxima->poly-list g vars ring-and-order))
	  (coerce-maxima-list poly-list-of-lists)))



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

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

(defun poly->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) (funcall *ratdisrep-fun* (coerce-to-maxima :polynomial p vars))) object)))
    (:term
     `((mtimes) ,(funcall *ratdisrep-fun* (term-coeff object))
       ,@(mapcar #'(lambda (var power) `((mexpt) ,var ,power))
		 vars (monom->list object))))
    ;; Assumes that Lisp and Maxima logicals coincide
    (:logical object)
    (otherwise
     object)))


#|


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Unary and binary operation definition facility
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(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)))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Facilities for evaluating Grobner package expressions
;; within a prepared environment
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(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-ring-and-order ((ring order) &body body)
  "Evaluate BODY with monomial order set to ORDER and coefficient ring set to RING."
  `(let ((*monomial-order* (or (find-order ,order) *monomial-order*))
	 (+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))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; 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.")

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; 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))))


;;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+)))))))

|#