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/grobner.lisp@ 131

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

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