;;; -*- Mode: Lisp -*- ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; ;;; Copyright (C) 1999, 2002, 2009, 2015 Marek Rychlik ;;; ;;; 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))) "Build RING-AND-ORDER structure. The defaults are determined by various Maxima-level switches, which are names of ring and orders." 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-and-order (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-list (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 (:custom object) ;Bypass processing (:polynomial `(,(maxima-head) ,@(mapcar #'(lambda (term) (poly->maxima :term term vars)) (poly-termlist object)))) (:poly-list `((mlist) ,@(mapcar #'(lambda (p) ($ratdisrep (poly->maxima :polynomial p vars))) object))) (:term `((mtimes) ,($ratdisrep (term-coeff object)) ,@(mapcar #'(lambda (var power) `((mexpt) ,var ,power)) vars (monom->list (term-monom 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-ring-and-order (((maxima-vars &optional (maxima-new-vars nil new-vars-supplied-p)) &key (polynomials nil) (poly-lists nil) (poly-list-lists nil) (value-type nil) (ring-and-order-var 'ring-and-order) (ring-var 'ring)) &body body &aux (vars (gensym)) (new-vars (gensym))) "Evaluate a polynomial expression BODY in an environment constructred from Maxima switches. The supplied arguments POLYNOMIALS, POLY-LISTS and POLY-LIST-LISTS should be polynomials, polynomial lists an lists of lists of polynomials, in Maxima general form. These are translated to NGROBNER package internal form and evaluated using operations in the NGROBNER package. The BODY should be defined in terms of those operations. MAXIMA-VARS is set to the list of variable names used at the Maxima level. The evaluation is performed by the NGROBNER package which ignores variable names, thus MAXIMA-VARS is used only to translate the polynomial expression to NGROBNER internal form. After evaluation, the value of BODY is translated back to the Maxima general form. When MAXIMA-NEW-VARS is present, it is appended to MAXIMA-VARS upon translation from the internal form back to Maxima general form, thus allowing extra variables which may have been created by the evaluation process. The value type can be either :POLYNOMIAL, :POLY-LIST or :TERM, depending on the form of the result returned by the top NGROBNER operation. During evaluation, symbols supplied by RING-AND-ORDER-VAR (defaul value 'RING-AND-ORDER), and RING-VAR (default value 'RING) are bound to RING-AND-ORDER and RING instances." `(let ((,vars (coerce-maxima-list ,maxima-vars)) ,@(when new-vars-supplied-p (list `(,new-vars (coerce-maxima-list ,maxima-new-vars))))) (poly->maxima ,value-type (let ((,ring-and-order-var ,(find-ring-and-order-by-name))) ;; Define a shorthand to RING (symbol-macrolet ((,ring-var (ro-ring ring-and-order))) (let ,(let ((args nil)) (dolist (p polynomials args) (setf args (cons `(,p (maxima->poly ,p ,vars ,ring-and-order-var)) args))) (dolist (p poly-lists args) (setf args (cons `(,p (maxima->poly-list ,p ,vars ,ring-and-order-var)) args))) (dolist (p poly-list-lists args) (setf args (cons `(,p (maxima->poly-list-list ,p ,vars ,ring-and-order-var)) 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-ring-and-order ((vars) :polynomials (p) :value-type :polynomial) p)) (defmfun $poly_expt (p n vars) (with-ring-and-order ((vars) :polynomials (p) :value-type :polynomial) (poly-expt ring-and-order p n))) (defmfun $poly_content (p vars) (with-ring-and-order ((vars) :polynomials (p)) (poly-content ring p))) #| (defmfun $poly_pseudo_divide (f fl vars) (with-ring-and-order ((vars) :polynomials (f) :poly-lists (fl) :value-type :custom) (multiple-value-bind (quot rem c division-count) (poly-pseudo-divide ring-and-order f fl) `((mlist) ,(poly->maxima :poly-list quot vars) ,(poly->maxima :polynomial rem vars) ,c ,division-count)))) |# (defmfun $poly_exact_divide (f g vars) (with-ring-and-order ((vars) :polynomials (f g) :value-type :polynomial) (poly-exact-divide ring-and-order f g))) (defmfun $poly_normal_form (f fl vars) (with-ring-and-order ((vars) :polynomials (f) :poly-lists (fl) :value-type :polynomial) (normal-form ring-and-order f (remzero fl) nil))) (defmfun $poly_buchberger_criterion (g vars) (with-ring-and-order ((vars) :poly-lists (g) :value-type :logical) (buchberger-criterion ring-and-order g))) (defmfun $poly_buchberger (fl vars) (with-ring-and-order ((vars) :poly-lists (fl) :value-type :poly-list) (buchberger ring-and-order (remzero fl) 0 nil))) (defmfun $poly_reduction (plist vars) (with-ring-and-order ((vars) :poly-lists (plist) :value-type :poly-list) (reduction ring-and-order plist))) (defmfun $poly_minimization (plist vars) (with-ring-and-order ((vars) :poly-lists (plist) :value-type :poly-list) (minimization plist))) (defmfun $poly_normalize_list (plist vars) (with-ring-and-order ((vars) :poly-lists (plist) :value-type :poly-list) (poly-normalize-list ring plist))) (defmfun $poly_grobner (f vars) (with-ring-and-order ((vars) :poly-lists (f) :value-type :poly-list) (grobner ring-and-order (remzero f)))) (defmfun $poly_reduced_grobner (f vars) (with-ring-and-order ((vars) :poly-lists (f) :value-type :poly-list) (reduced-grobner ring-and-order (remzero f)))) (defmfun $poly_depends_p (p var mvars &aux (vars (coerce-maxima-list mvars)) (pos (position var vars))) (with-ring-and-order ((mvars) :polynomials (p) :value-type :custom) (if (null pos) (merror "~%Variable ~M not in the list of variables ~M." var mvars) (poly-depends-p p pos)))) (defmfun $poly_elimination_ideal (flist k vars) (with-ring-and-order ((vars) :poly-lists (flist) :value-type :poly-list) (elimination-ideal ring-and-order flist k nil 0))) (defmfun $poly_colon_ideal (f g vars) (with-ring-and-order ((vars) :poly-lists (f g) :value-type :poly-list) (colon-ideal ring-and-order f g nil))) (defmfun $poly_ideal_intersection (f g vars) (with-ring-and-order ((vars) :poly-lists (f g) :value-type :poly-list) (ideal-intersection ring-and-order f g nil))) (defmfun $poly_lcm (f g vars) (with-ring-and-order ((vars) :polynomials (f g) :value-type :polynomial) (poly-lcm ring-and-order 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-ring-and-order ((vars) :poly-lists (g1 g2)) (grobner-equal ring-and-order g1 g2))) (defmfun $poly_grobner_subsetp (g1 g2 vars) (with-ring-and-order ((vars) :poly-lists (g1 g2)) (grobner-subsetp ring-and-order g1 g2))) (defmfun $poly_grobner_member (p g vars) (with-ring-and-order ((vars) :polynomials (p) :poly-lists (g)) (grobner-member ring-and-order p g))) (defmfun $poly_ideal_saturation1 (f p vars) (with-ring-and-order ((vars) :poly-lists (f) :polynomials (p) :value-type :poly-list) (ideal-saturation-1 ring f p 0))) (defmfun $poly_saturation_extension (f plist vars new-vars) (with-ring-and-order ((vars new-vars) :poly-lists (f plist) :value-type :poly-list) (saturation-extension ring f plist))) (defmfun $poly_polysaturation_extension (f plist vars new-vars) (with-ring-and-order ((vars new-vars) :poly-lists (f plist) :value-type :poly-list) (polysaturation-extension ring f plist))) (defmfun $poly_ideal_polysaturation1 (f plist vars) (with-ring-and-order ((vars) :poly-lists (f plist) :value-type :poly-list) (ideal-polysaturation-1 ring f plist 0 nil))) (defmfun $poly_ideal_saturation (f g vars) (with-ring-and-order ((vars) :poly-lists (f g) :value-type :poly-list) (ideal-saturation ring f g 0 nil))) (defmfun $poly_ideal_polysaturation (f ideal-list vars) (with-ring-and-order ((vars) :poly-lists (f) :poly-list-lists (ideal-list) :value-type :poly-list) (ideal-polysaturation ring-and-order f ideal-list 0 nil))) (defmfun $poly_lt (f vars) (with-ring-and-order ((vars) :polynomials (f) :value-type :polynomial) (make-poly-from-termlist (list (poly-lt f))))) (defmfun $poly_lm (f vars) (with-ring-and-order ((vars) :polynomials (f) :value-type :polynomial) (make-poly-from-termlist (list (make-term (poly-lm f) (funcall (ring-unit ring)))))))