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@ 123

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

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