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

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