close Warning: Can't synchronize with repository "(default)" (The repository directory has changed, you should resynchronize the repository with: trac-admin $ENV repository resync '(default)'). Look in the Trac log for more information.

source: branches/f4grobner/ngrobner.lisp@ 220

Last change on this file since 220 was 220, checked in by Marek Rychlik, 9 years ago

* empty log message *
File size: 11.2 KB
Line 
1;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*-
2;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3;;;
4;;; Copyright (C) 1999, 2002, 2009, 2015 Marek Rychlik <rychlik@u.arizona.edu>
5;;;
6;;; This program is free software; you can redistribute it and/or modify
7;;; it under the terms of the GNU General Public License as published by
8;;; the Free Software Foundation; either version 2 of the License, or
9;;; (at your option) any later version.
10;;;
11;;; This program is distributed in the hope that it will be useful,
12;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
13;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14;;; GNU General Public License for more details.
15;;;
16;;; You should have received a copy of the GNU General Public License
17;;; along with this program; if not, write to the Free Software
18;;; Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
19;;;
20;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
21
22(in-package :ngrobner)
23
24;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
25;;
26;; Global switches
27;;
28;; Can be used in Maxima just fine, as they observe the
29;; Maxima naming convention, i.e. all names visible at the
30;; Maxima toplevel begin with a '$'.
31;;
32;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
33
34(defvar $poly_monomial_order '$lex
35 "This switch controls which monomial order is used in polynomial
36and Grobner basis calculations. If not set, LEX will be used")
37
38(defvar $poly_coefficient_ring '$expression_ring
39 "This switch indicates the coefficient ring of the polynomials
40that will be used in grobner calculations. If not set, Maxima's
41general expression ring will be used. This variable may be set
42to RING_OF_INTEGERS if desired.")
43
44(defvar $poly_primary_elimination_order nil
45 "Name of the default order for eliminated variables in elimination-based functions.
46If not set, LEX will be used.")
47
48(defvar $poly_secondary_elimination_order nil
49 "Name of the default order for kept variables in elimination-based functions.
50If not set, LEX will be used.")
51
52(defvar $poly_elimination_order nil
53 "Name of the default elimination order used in elimination calculations.
54If set, it overrides the settings in variables POLY_PRIMARY_ELIMINATION_ORDER
55and SECONDARY_ELIMINATION_ORDER. The user must ensure that this is a true
56elimination order valid for the number of eliminated variables.")
57
58(defvar $poly_return_term_list nil
59 "If set to T, all functions in this package will return each polynomial as a
60list of terms in the current monomial order rather than a Maxima general expression.")
61
62(defvar $poly_grobner_debug nil
63 "If set to TRUE, produce debugging and tracing output.")
64
65(defvar $poly_grobner_algorithm '$buchberger
66 "The name of the algorithm used to find grobner bases.")
67
68(defvar $poly_top_reduction_only nil
69 "If not FALSE, use top reduction only whenever possible.
70Top reduction means that division algorithm stops after the first reduction.")
71
72
73
74;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
75;;
76;; Coefficient ring operations
77;;
78;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
79;;
80;; These are ALL operations that are performed on the coefficients by
81;; the package, and thus the coefficient ring can be changed by merely
82;; redefining these operations.
83;;
84;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
85
86(defstruct (ring)
87 (parse #'identity :type function)
88 (unit #'identity :type function)
89 (zerop #'identity :type function)
90 (add #'identity :type function)
91 (sub #'identity :type function)
92 (uminus #'identity :type function)
93 (mul #'identity :type function)
94 (div #'identity :type function)
95 (lcm #'identity :type function)
96 (ezgcd #'identity :type function)
97 (gcd #'identity :type function))
98
99(defparameter *ring-of-integers*
100 (make-ring
101 :parse #'identity
102 :unit #'(lambda () 1)
103 :zerop #'zerop
104 :add #'+
105 :sub #'-
106 :uminus #'-
107 :mul #'*
108 :div #'/
109 :lcm #'lcm
110 :ezgcd #'(lambda (x y &aux (c (gcd x y))) (values c (/ x c) (/ y c)))
111 :gcd #'gcd)
112 "The ring of integers.")
113
114(defvar *expression-ring* *ring-of-integers*
115 "The ring of coefficients, over which all polynomials are assumed to
116 be defined.")
117
118(defvar *ratdisrep-fun* #'identity
119 "A function applied to polynomials after conversion to Maxima representation.")
120
121
122
123;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
124;;
125;; This is how we perform operations on coefficients
126;; using Maxima functions.
127;;
128;; Functions and macros dealing with internal representation structure
129;;
130;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
131
132
133;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
134;;
135;; Debugging/tracing
136;;
137;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
138(defmacro debug-cgb (&rest args)
139 `(when $poly_grobner_debug (format *terminal-io* ,@args)))
140
141
142
143;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
144;;
145;; These are provided mostly for debugging purposes To enable
146;; verification of grobner bases with BUCHBERGER-CRITERION, do
147;; (pushnew :grobner-check *features*) and compile/load this file.
148;; With this feature, the calculations will slow down CONSIDERABLY.
149;;
150;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
151
152(defun grobner-test (ring g f)
153 "Test whether G is a Grobner basis and F is contained in G. Return T
154upon success and NIL otherwise."
155 (debug-cgb "~&GROBNER CHECK: ")
156 (let (($poly_grobner_debug nil)
157 (stat1 (buchberger-criterion ring g))
158 (stat2
159 (every #'poly-zerop
160 (makelist (normal-form ring (copy-tree (elt f i)) g nil)
161 (i 0 (1- (length f)))))))
162 (unless stat1 (error "~&Buchberger criterion failed."))
163 (unless stat2
164 (error "~&Original polys not in ideal spanned by Grobner.")))
165 (debug-cgb "~&GROBNER CHECK END")
166 t)
167
168
169;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
170;;
171;; Selection of algorithm and pair heuristic
172;;
173;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
174
175(defun find-grobner-function (algorithm)
176 "Return a function which calculates Grobner basis, based on its
177names. Names currently used are either Lisp symbols, Maxima symbols or
178keywords."
179 (ecase algorithm
180 ((buchberger :buchberger $buchberger) #'buchberger)
181 ((parallel-buchberger :parallel-buchberger $parallel_buchberger) #'parallel-buchberger)
182 ((gebauer-moeller :gebauer_moeller $gebauer_moeller) #'gebauer-moeller)))
183
184(defun grobner (ring f &optional (start 0) (top-reduction-only nil))
185 ;;(setf F (sort F #'< :key #'sugar))
186 (funcall
187 (find-grobner-function $poly_grobner_algorithm)
188 ring f start top-reduction-only))
189
190(defun reduced-grobner (ring f &optional (start 0) (top-reduction-only $poly_top_reduction_only))
191 (reduction ring (grobner ring f start top-reduction-only)))
192
193
194
195
196;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
197;;
198;; Conversion from internal to infix form
199;;
200;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
201
202(defun coerce-to-infix (poly-type object vars)
203 (case poly-type
204 (:termlist
205 `(+ ,@(mapcar #'(lambda (term) (coerce-to-infix :term term vars)) object)))
206 (:polynomial
207 (coerce-to-infix :termlist (poly-termlist object) vars))
208 (:poly-list
209 `([ ,@(mapcar #'(lambda (p) (coerce-to-infix :polynomial p vars)) object)))
210 (:term
211 `(* ,(term-coeff object)
212 ,@(mapcar #'(lambda (var power) `(expt ,var ,power))
213 vars (monom-exponents (term-monom object)))))
214 (otherwise
215 object)))
216
217
218;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
219;;
220;; Order utilities
221;;
222;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
223(defun find-order (order)
224 "This function returns the order function bases on its name."
225 (cond
226 ((null order) nil)
227 ((symbolp order)
228 (case order
229 ((lex :lex $lex) #'lex>)
230 ((grlex :grlex $grlex) #'grlex>)
231 ((grevlex :grevlex $grevlex) #'grevlex>)
232 ((invlex :invlex $invlex) #'invlex>)
233 ((elimination-order-1 :elimination-order-1 elimination_order_1) #'elimination-order-1)
234 (otherwise
235 (warn "~%Warning: Order ~A not found. Using default.~%" order))))
236 (t
237 (warn "~%Order specification ~A is not recognized. Using default.~%" order)
238 nil)))
239
240(defun find-ring (ring)
241 "This function returns the ring structure bases on input symbol."
242 (cond
243 ((null ring) nil)
244 ((symbolp ring)
245 (case ring
246 ((expression-ring :expression-ring $expression_ring) *expression-ring*)
247 ((ring-of-integers :ring-of-integers $ring_of_integers) *ring-of-integers*)
248 (otherwise
249 (warn "~%Warning: Ring ~A not found. Using default.~%" ring))))
250 (t
251 (warn "~%Ring specification ~A is not recognized. Using default.~%" ring)
252 nil)))
253
254
255
256;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
257;;
258;; Conversion from internal form to Maxima general form
259;;
260;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
261
262(defun maxima-head ()
263 (if $poly_return_term_list
264 '(mlist)
265 '(mplus)))
266
267(defun coerce-to-maxima (poly-type object vars)
268 (case poly-type
269 (:polynomial
270 `(,(maxima-head) ,@(mapcar #'(lambda (term) (coerce-to-maxima :term term vars)) (poly-termlist object))))
271 (:poly-list
272 `((mlist) ,@(mapcar #'(lambda (p) (funcall *ratdisrep-fun* (coerce-to-maxima :polynomial p vars))) object)))
273 (:term
274 `((mtimes) ,(funcall *ratdisrep-fun* (term-coeff object))
275 ,@(mapcar #'(lambda (var power) `((mexpt) ,var ,power))
276 vars (monom-exponents (term-monom object)))))
277 ;; Assumes that Lisp and Maxima logicals coincide
278 (:logical object)
279 (otherwise
280 object)))
281
282
283
284;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
285;;
286;; Macro facility for writing Maxima-level wrappers for
287;; functions operating on internal representation
288;;
289;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
290
291(defmacro with-parsed-polynomials (((maxima-vars &optional (maxima-new-vars nil new-vars-supplied-p))
292 &key (polynomials nil)
293 (poly-lists nil)
294 (poly-list-lists nil)
295 (value-type nil))
296 &body body
297 &aux (vars (gensym))
298 (new-vars (gensym)))
299 `(let ((,vars (coerce-maxima-list ,maxima-vars))
300 ,@(when new-vars-supplied-p
301 (list `(,new-vars (coerce-maxima-list ,maxima-new-vars)))))
302 (coerce-to-maxima
303 ,value-type
304 (with-coefficient-ring ($poly_coefficient_ring)
305 (with-monomial-order ($poly_monomial_order)
306 (with-elimination-orders ($poly_primary_elimination_order
307 $poly_secondary_elimination_order
308 $poly_elimination_order)
309 (let ,(let ((args nil))
310 (dolist (p polynomials args)
311 (setf args (cons `(,p (parse-poly ,p ,vars)) args)))
312 (dolist (p poly-lists args)
313 (setf args (cons `(,p (parse-poly-list ,p ,vars)) args)))
314 (dolist (p poly-list-lists args)
315 (setf args (cons `(,p (parse-poly-list-list ,p ,vars)) args))))
316 . ,body))))
317 ,(if new-vars-supplied-p
318 `(append ,vars ,new-vars)
319 vars))))
320
321
Note: See TracBrowser for help on using the repository browser.