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

Last change on this file since 269 was 269, checked in by Marek Rychlik, 10 years ago

* empty log message *
File size: 17.0 KB
RevLine 
[98]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
[133]22;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
23;;
[268]24;; Load this file into Maxima to bootstrap the Grobner package.
[133]25;;
[268]26;; DETAILS: This file implements an interface between the Grobner
[269]27;; basis package NGROBNER and Maxima. NGROBNER for efficiency uses its
28;; own representation of polynomials. Thus, it is necessary to convert
29;; Maxima representation to the internal representation and back. The
30;; facilities to do so are implemented in this file.
[268]31;;
[133]32;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
33
[98]34(in-package :maxima)
35
36(macsyma-module cgb-maxima)
37
38(eval-when
39 #+gcl (load eval)
40 #-gcl (:load-toplevel :execute)
41 (format t "~&Loading maxima-grobner ~a ~a~%"
42 "$Revision: 2.0 $" "$Date: 2015/06/02 0:34:17 $"))
43
44;;FUNCTS is loaded because it contains the definition of LCM
45($load "functs")
[152]46
[258]47
48(defvar *ngrobner-files* '("ngrobner-package" "utils" "ngrobner" "monomial"
49 "order" "order-mk" "term" "termlist" "polynomial" "priority-queue"
50 "pair-queue" "division" "criterion" "buchberger" "gebauer-moeller"
51 "gb-postprocessing" "ideal")
52 "List of files in the NGROBNER package")
53
[254]54;; Compile/load NGROBNER package files
[252]55(eval-when
[262]56 #+gcl (load)
[252]57 #-gcl (:load-toplevel :execute)
[258]58 (dolist (file *ngrobner-files*)
[265]59 (load file :verbose t)))
[157]60
[229]61(use-package :ngrobner)
62
[98]63;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
64;;
65;; Maxima expression ring
66;;
67;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
68
[230]69(defparameter *maxima-ring*
70 (make-ring
[98]71 ;;(defun coeff-zerop (expr) (meval1 `(($is) (($equal) ,expr 0))))
72 :parse #'(lambda (expr)
73 (when modulus (setf expr ($rat expr)))
74 expr)
75 :unit #'(lambda () (if modulus ($rat 1) 1))
76 :zerop #'(lambda (expr)
77 ;;When is exactly a maxima expression equal to 0?
78 (cond ((numberp expr)
79 (= expr 0))
80 ((atom expr) nil)
81 (t
82 (case (caar expr)
83 (mrat (eql ($ratdisrep expr) 0))
84 (otherwise (eql ($totaldisrep expr) 0))))))
85 :add #'(lambda (x y) (m+ x y))
86 :sub #'(lambda (x y) (m- x y))
87 :uminus #'(lambda (x) (m- x))
88 :mul #'(lambda (x y) (m* x y))
89 ;;(defun coeff-div (x y) (cadr ($divide x y)))
90 :div #'(lambda (x y) (m// x y))
91 :lcm #'(lambda (x y) (meval1 `((|$LCM|) ,x ,y)))
92 :ezgcd #'(lambda (x y) (apply #'values (cdr ($ezgcd ($totaldisrep x) ($totaldisrep y)))))
93 ;; :gcd #'(lambda (x y) (second ($ezgcd x y)))))
94 :gcd #'(lambda (x y) ($gcd x y))))
95
[237]96(setf *expression-ring* *maxima-ring*)
97
[114]98;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
99;;
100;; Maxima expression parsing
101;;
102;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
103
104(defun equal-test-p (expr1 expr2)
105 (alike1 expr1 expr2))
106
107(defun coerce-maxima-list (expr)
108 "convert a maxima list to lisp list."
109 (cond
110 ((and (consp (car expr)) (eql (caar expr) 'mlist)) (cdr expr))
111 (t expr)))
112
113(defun free-of-vars (expr vars) (apply #'$freeof `(,@vars ,expr)))
114
115(defun parse-poly (expr vars &aux (vars (coerce-maxima-list vars)))
116 "Convert a maxima polynomial expression EXPR in variables VARS to internal form."
117 (labels ((parse (arg) (parse-poly arg vars))
118 (parse-list (args) (mapcar #'parse args)))
119 (cond
120 ((eql expr 0) (make-poly-zero))
121 ((member expr vars :test #'equal-test-p)
122 (let ((pos (position expr vars :test #'equal-test-p)))
[233]123 (make-variable *expression-ring* (length vars) pos)))
[114]124 ((free-of-vars expr vars)
125 ;;This means that variable-free CRE and Poisson forms will be converted
126 ;;to coefficients intact
[233]127 (coerce-coeff *expression-ring* expr vars))
[114]128 (t
129 (case (caar expr)
[233]130 (mplus (reduce #'(lambda (x y) (poly-add *expression-ring* x y)) (parse-list (cdr expr))))
131 (mminus (poly-uminus *expression-ring* (parse (cadr expr))))
[114]132 (mtimes
133 (if (endp (cddr expr)) ;unary
134 (parse (cdr expr))
[233]135 (reduce #'(lambda (p q) (poly-mul *expression-ring* p q)) (parse-list (cdr expr)))))
[114]136 (mexpt
137 (cond
138 ((member (cadr expr) vars :test #'equal-test-p)
139 ;;Special handling of (expt var pow)
140 (let ((pos (position (cadr expr) vars :test #'equal-test-p)))
[233]141 (make-variable *expression-ring* (length vars) pos (caddr expr))))
[114]142 ((not (and (integerp (caddr expr)) (plusp (caddr expr))))
143 ;; Negative power means division in coefficient ring
144 ;; Non-integer power means non-polynomial coefficient
145 (mtell "~%Warning: Expression ~%~M~%contains power which is not a positive integer. Parsing as coefficient.~%"
146 expr)
[233]147 (coerce-coeff *expression-ring* expr vars))
148 (t (poly-expt *expression-ring* (parse (cadr expr)) (caddr expr)))))
[114]149 (mrat (parse ($ratdisrep expr)))
150 (mpois (parse ($outofpois expr)))
151 (otherwise
[233]152 (coerce-coeff *expression-ring* expr vars)))))))
[114]153
154(defun parse-poly-list (expr vars)
155 (case (caar expr)
156 (mlist (mapcar #'(lambda (p) (parse-poly p vars)) (cdr expr)))
157 (t (merror "Expression ~M is not a list of polynomials in variables ~M."
158 expr vars))))
159(defun parse-poly-list-list (poly-list-list vars)
160 (mapcar #'(lambda (g) (parse-poly-list g vars)) (coerce-maxima-list poly-list-list)))
161
162
[111]163;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
164;;
[241]165;; Conversion from internal form to Maxima general form
166;;
167;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
168
169(defun maxima-head ()
170 (if $poly_return_term_list
171 '(mlist)
172 '(mplus)))
173
174(defun coerce-to-maxima (poly-type object vars)
175 (case poly-type
176 (:polynomial
177 `(,(maxima-head) ,@(mapcar #'(lambda (term) (coerce-to-maxima :term term vars)) (poly-termlist object))))
178 (:poly-list
179 `((mlist) ,@(mapcar #'(lambda (p) (funcall *ratdisrep-fun* (coerce-to-maxima :polynomial p vars))) object)))
180 (:term
181 `((mtimes) ,(funcall *ratdisrep-fun* (term-coeff object))
182 ,@(mapcar #'(lambda (var power) `((mexpt) ,var ,power))
183 vars (monom-exponents (term-monom object)))))
184 ;; Assumes that Lisp and Maxima logicals coincide
185 (:logical object)
186 (otherwise
187 object)))
188
189
190;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
191;;
[111]192;; Unary and binary operation definition facility
193;;
194;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
[98]195
[111]196(defmacro define-unop (maxima-name fun-name
197 &optional (documentation nil documentation-supplied-p))
198 "Define a MAXIMA-level unary operator MAXIMA-NAME corresponding to unary function FUN-NAME."
199 `(defun ,maxima-name (p vars
200 &aux
201 (vars (coerce-maxima-list vars))
202 (p (parse-poly p vars)))
203 ,@(when documentation-supplied-p (list documentation))
[233]204 (coerce-to-maxima :polynomial (,fun-name *expression-ring* p) vars)))
[111]205
206(defmacro define-binop (maxima-name fun-name
207 &optional (documentation nil documentation-supplied-p))
208 "Define a MAXIMA-level binary operator MAXIMA-NAME corresponding to binary function FUN-NAME."
209 `(defmfun ,maxima-name (p q vars
210 &aux
211 (vars (coerce-maxima-list vars))
212 (p (parse-poly p vars))
213 (q (parse-poly q vars)))
214 ,@(when documentation-supplied-p (list documentation))
[233]215 (coerce-to-maxima :polynomial (,fun-name *expression-ring* p q) vars)))
[111]216
217
[219]218;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
219;;
220;; Facilities for evaluating Grobner package expressions
221;; within a prepared environment
222;;
223;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
224
225(defmacro with-monomial-order ((order) &body body)
226 "Evaluate BODY with monomial order set to ORDER."
227 `(let ((*monomial-order* (or (find-order ,order) *monomial-order*)))
228 . ,body))
229
230(defmacro with-coefficient-ring ((ring) &body body)
231 "Evaluate BODY with coefficient ring set to RING."
[233]232 `(let ((*expression-ring* (or (find-ring ,ring) *expression-ring*)))
[219]233 . ,body))
234
235(defmacro with-elimination-orders ((primary secondary elimination-order)
236 &body body)
237 "Evaluate BODY with primary and secondary elimination orders set to PRIMARY and SECONDARY."
238 `(let ((*primary-elimination-order* (or (find-order ,primary) *primary-elimination-order*))
239 (*secondary-elimination-order* (or (find-order ,secondary) *secondary-elimination-order*))
240 (*elimination-order* (or (find-order ,elimination-order) *elimination-order*)))
241 . ,body))
242
243
[98]244;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
245;;
246;; Maxima-level interface functions
247;;
248;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
249
250;; Auxillary function for removing zero polynomial
251(defun remzero (plist) (remove #'poly-zerop plist))
252
253;;Simple operators
254
255(define-binop $poly_add poly-add
256 "Adds two polynomials P and Q")
257
258(define-binop $poly_subtract poly-sub
259 "Subtracts a polynomial Q from P.")
260
261(define-binop $poly_multiply poly-mul
262 "Returns the product of polynomials P and Q.")
263
264(define-binop $poly_s_polynomial spoly
265 "Returns the syzygy polynomial (S-polynomial) of two polynomials P and Q.")
266
267(define-unop $poly_primitive_part poly-primitive-part
268 "Returns the polynomial P divided by GCD of its coefficients.")
269
270(define-unop $poly_normalize poly-normalize
271 "Returns the polynomial P divided by the leading coefficient.")
272
[222]273;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
274;;
275;; Macro facility for writing Maxima-level wrappers for
276;; functions operating on internal representation
277;;
278;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
279
280(defmacro with-parsed-polynomials (((maxima-vars &optional (maxima-new-vars nil new-vars-supplied-p))
281 &key (polynomials nil)
282 (poly-lists nil)
283 (poly-list-lists nil)
284 (value-type nil))
285 &body body
286 &aux (vars (gensym))
287 (new-vars (gensym)))
288 `(let ((,vars (coerce-maxima-list ,maxima-vars))
289 ,@(when new-vars-supplied-p
290 (list `(,new-vars (coerce-maxima-list ,maxima-new-vars)))))
291 (coerce-to-maxima
292 ,value-type
293 (with-coefficient-ring ($poly_coefficient_ring)
294 (with-monomial-order ($poly_monomial_order)
295 (with-elimination-orders ($poly_primary_elimination_order
296 $poly_secondary_elimination_order
297 $poly_elimination_order)
298 (let ,(let ((args nil))
299 (dolist (p polynomials args)
300 (setf args (cons `(,p (parse-poly ,p ,vars)) args)))
301 (dolist (p poly-lists args)
302 (setf args (cons `(,p (parse-poly-list ,p ,vars)) args)))
303 (dolist (p poly-list-lists args)
304 (setf args (cons `(,p (parse-poly-list-list ,p ,vars)) args))))
305 . ,body))))
306 ,(if new-vars-supplied-p
307 `(append ,vars ,new-vars)
308 vars))))
309
310
[98]311;;Functions
312
313(defmfun $poly_expand (p vars)
314 "This function is equivalent to EXPAND(P) if P parses correctly to a polynomial.
315If the representation is not compatible with a polynomial in variables VARS,
316the result is an error."
317 (with-parsed-polynomials ((vars) :polynomials (p)
318 :value-type :polynomial)
319 p))
320
321(defmfun $poly_expt (p n vars)
322 (with-parsed-polynomials ((vars) :polynomials (p) :value-type :polynomial)
[233]323 (poly-expt *expression-ring* p n)))
[98]324
325(defmfun $poly_content (p vars)
326 (with-parsed-polynomials ((vars) :polynomials (p))
[233]327 (poly-content *expression-ring* p)))
[98]328
329(defmfun $poly_pseudo_divide (f fl vars
330 &aux (vars (coerce-maxima-list vars))
331 (f (parse-poly f vars))
332 (fl (parse-poly-list fl vars)))
333 (multiple-value-bind (quot rem c division-count)
[233]334 (poly-pseudo-divide *expression-ring* f fl)
[98]335 `((mlist)
336 ,(coerce-to-maxima :poly-list quot vars)
337 ,(coerce-to-maxima :polynomial rem vars)
338 ,c
339 ,division-count)))
340
341(defmfun $poly_exact_divide (f g vars)
342 (with-parsed-polynomials ((vars) :polynomials (f g) :value-type :polynomial)
[233]343 (poly-exact-divide *expression-ring* f g)))
[98]344
345(defmfun $poly_normal_form (f fl vars)
346 (with-parsed-polynomials ((vars) :polynomials (f)
347 :poly-lists (fl)
348 :value-type :polynomial)
[233]349 (normal-form *expression-ring* f (remzero fl) nil)))
[98]350
351(defmfun $poly_buchberger_criterion (g vars)
352 (with-parsed-polynomials ((vars) :poly-lists (g) :value-type :logical)
[233]353 (buchberger-criterion *expression-ring* g)))
[98]354
355(defmfun $poly_buchberger (fl vars)
356 (with-parsed-polynomials ((vars) :poly-lists (fl) :value-type :poly-list)
[233]357 (buchberger *expression-ring* (remzero fl) 0 nil)))
[98]358
359(defmfun $poly_reduction (plist vars)
360 (with-parsed-polynomials ((vars) :poly-lists (plist)
361 :value-type :poly-list)
[233]362 (reduction *expression-ring* plist)))
[98]363
364(defmfun $poly_minimization (plist vars)
365 (with-parsed-polynomials ((vars) :poly-lists (plist)
366 :value-type :poly-list)
367 (minimization plist)))
368
369(defmfun $poly_normalize_list (plist vars)
370 (with-parsed-polynomials ((vars) :poly-lists (plist)
371 :value-type :poly-list)
[233]372 (poly-normalize-list *expression-ring* plist)))
[98]373
374(defmfun $poly_grobner (f vars)
375 (with-parsed-polynomials ((vars) :poly-lists (f)
376 :value-type :poly-list)
[233]377 (grobner *expression-ring* (remzero f))))
[98]378
379(defmfun $poly_reduced_grobner (f vars)
380 (with-parsed-polynomials ((vars) :poly-lists (f)
381 :value-type :poly-list)
[233]382 (reduced-grobner *expression-ring* (remzero f))))
[98]383
384(defmfun $poly_depends_p (p var mvars
385 &aux (vars (coerce-maxima-list mvars))
386 (pos (position var vars)))
387 (if (null pos)
388 (merror "~%Variable ~M not in the list of variables ~M." var mvars)
389 (poly-depends-p (parse-poly p vars) pos)))
390
391(defmfun $poly_elimination_ideal (flist k vars)
392 (with-parsed-polynomials ((vars) :poly-lists (flist)
393 :value-type :poly-list)
[233]394 (elimination-ideal *expression-ring* flist k nil 0)))
[98]395
396(defmfun $poly_colon_ideal (f g vars)
397 (with-parsed-polynomials ((vars) :poly-lists (f g) :value-type :poly-list)
[233]398 (colon-ideal *expression-ring* f g nil)))
[98]399
400(defmfun $poly_ideal_intersection (f g vars)
401 (with-parsed-polynomials ((vars) :poly-lists (f g) :value-type :poly-list)
[233]402 (ideal-intersection *expression-ring* f g nil)))
[98]403
404(defmfun $poly_lcm (f g vars)
405 (with-parsed-polynomials ((vars) :polynomials (f g) :value-type :polynomial)
[233]406 (poly-lcm *expression-ring* f g)))
[98]407
408(defmfun $poly_gcd (f g vars)
409 ($first ($divide (m* f g) ($poly_lcm f g vars))))
410
411(defmfun $poly_grobner_equal (g1 g2 vars)
412 (with-parsed-polynomials ((vars) :poly-lists (g1 g2))
[233]413 (grobner-equal *expression-ring* g1 g2)))
[98]414
415(defmfun $poly_grobner_subsetp (g1 g2 vars)
416 (with-parsed-polynomials ((vars) :poly-lists (g1 g2))
[233]417 (grobner-subsetp *expression-ring* g1 g2)))
[98]418
419(defmfun $poly_grobner_member (p g vars)
420 (with-parsed-polynomials ((vars) :polynomials (p) :poly-lists (g))
[233]421 (grobner-member *expression-ring* p g)))
[98]422
423(defmfun $poly_ideal_saturation1 (f p vars)
424 (with-parsed-polynomials ((vars) :poly-lists (f) :polynomials (p)
425 :value-type :poly-list)
[233]426 (ideal-saturation-1 *expression-ring* f p 0)))
[98]427
428(defmfun $poly_saturation_extension (f plist vars new-vars)
429 (with-parsed-polynomials ((vars new-vars)
430 :poly-lists (f plist)
431 :value-type :poly-list)
[233]432 (saturation-extension *expression-ring* f plist)))
[98]433
434(defmfun $poly_polysaturation_extension (f plist vars new-vars)
435 (with-parsed-polynomials ((vars new-vars)
436 :poly-lists (f plist)
437 :value-type :poly-list)
[233]438 (polysaturation-extension *expression-ring* f plist)))
[98]439
440(defmfun $poly_ideal_polysaturation1 (f plist vars)
441 (with-parsed-polynomials ((vars) :poly-lists (f plist)
442 :value-type :poly-list)
[233]443 (ideal-polysaturation-1 *expression-ring* f plist 0 nil)))
[98]444
445(defmfun $poly_ideal_saturation (f g vars)
446 (with-parsed-polynomials ((vars) :poly-lists (f g)
447 :value-type :poly-list)
[233]448 (ideal-saturation *expression-ring* f g 0 nil)))
[98]449
450(defmfun $poly_ideal_polysaturation (f ideal-list vars)
451 (with-parsed-polynomials ((vars) :poly-lists (f)
452 :poly-list-lists (ideal-list)
453 :value-type :poly-list)
[233]454 (ideal-polysaturation *expression-ring* f ideal-list 0 nil)))
[98]455
456(defmfun $poly_lt (f vars)
457 (with-parsed-polynomials ((vars) :polynomials (f) :value-type :polynomial)
458 (make-poly-from-termlist (list (poly-lt f)))))
459
460(defmfun $poly_lm (f vars)
461 (with-parsed-polynomials ((vars) :polynomials (f) :value-type :polynomial)
[233]462 (make-poly-from-termlist (list (make-term (poly-lm f) (funcall (ring-unit *expression-ring*)))))))
[98]463
Note: See TracBrowser for help on using the repository browser.