;;; -*- 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) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; ;; Global switches ;; ;; Can be used in Maxima just fine, as they observe the ;; Maxima naming convention, i.e. all names visible at the ;; Maxima toplevel begin with a '$'. ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defvar $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") (defvar $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.") (defvar $poly_primary_elimination_order nil "Name of the default order for eliminated variables in elimination-based functions. If not set, LEX will be used.") (defvar $poly_secondary_elimination_order nil "Name of the default order for kept variables in elimination-based functions. If not set, LEX will be used.") (defvar $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.") (defvar $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.") ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; ;; 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." (declare (type ring-and-order ring-and-order)) (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." (declare (type ring-and-order ring-and-order)) (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))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; ;; 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)))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; ;; N-ary (unary and binary) operation definition facility ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defmacro define-op (maxima-name ;Name of maxima level function (fun-name env &rest args) ;Lisp level form to evaluate &optional (documentation nil documentation-supplied-p) &aux ;; The argument passed as first arg (env-arg (ecase env (:ring-and-order 'ring-and-order) (:ring 'ring)))) "Define a MAXIMA-level unary operator MAXIMA-NAME corresponding to unary function FUN-NAME. The second argument should be :RING or :RING-AND-ORDER, and it signals the type of the first argument that should be passed to function FUN-NAME. ARGS is a list of formal parameters passed to the function, i.e. symbols used as arguments. The macro expands to a Maxima-level function definition with name MAXIMA-NAME, which wraps FUN-NAME." `(defmfun ,maxima-name (,@args vars) ,@(when documentation-supplied-p (list documentation)) (with-ring-and-order ((vars) :polynomials (,@args) :value-type :polynomial) (,fun-name ,env-arg ,@args)))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; ;; Maxima-level interface functions ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Auxillary function for removing zero polynomial (defun remzero (plist) (remove #'poly-zerop plist)) ;;Simple operators (define-op $poly_add (poly-add :ring-and-order p q) "Adds two polynomials P and Q") (define-op $poly_subtract (poly-sub :ring-and-order p q) "Subtracts a polynomial Q from P.") (define-op $poly_multiply (poly-mul :ring-and-order p q) "Returns the product of polynomials P and Q.") (define-op $poly_s_polynomial (spoly :ring-and-order p q) "Returns the syzygy polynomial (S-polynomial) of two polynomials P and Q.") (define-op $poly_primitive_part (poly-primitive-part :ring p) "Returns the polynomial P divided by GCD of its coefficients.") (define-op $poly_normalize (poly-normalize :ring p) "Returns the polynomial P divided by the leading coefficient.") ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; ;; More complex 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 mvars &aux (vars (coerce-maxima-list mvars))) (with-ring-and-order ((mvars) :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-and-order 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-and-order 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-and-order 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 :monom (poly-lm f) :coeff (funcall (ring-unit ring)))))))