1
;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*-
2
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4
;;; $Id: grobner.lisp,v 1.1 2006/02/07 04:49:49 robert_dodier Exp $
5
;;; Copyright (C) 1999, 2002 Marek Rychlik <rychlik@u.arizona.edu>
7
;;; This program is free software; you can redistribute it and/or modify
8
;;; it under the terms of the GNU General Public License as published by
9
;;; the Free Software Foundation; either version 2 of the License, or
10
;;; (at your option) any later version.
12
;;; This program is distributed in the hope that it will be useful,
13
;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
14
;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15
;;; GNU General Public License for more details.
17
;;; You should have received a copy of the GNU General Public License
18
;;; along with this program; if not, write to the Free Software
19
;;; Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
21
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
24
(macsyma-module cgb-maxima)
26
(eval-when (load eval)
27
(format t "~&Loading maxima-grobner ~a ~a~%"
28
"$Revision: 1.1 $" "$Date: 2006/02/07 04:49:49 $"))
30
;;FUNCTS is loaded because it contains the definition of LCM
33
;; Macros for making lists with iterators - an exammple of GENSYM
34
;; MAKELIST-1 makes a list with one iterator, while MAKELIST accepts an
35
;; arbitrary number of iterators
39
;; >(makelist-1 (* 2 i) i 0 10)
40
;; (0 2 4 6 8 10 12 14 16 18 20)
42
;; >(makelist-1 (* 2 i) i 0 10 3)
45
;; Generate sums of squares of numbers between 1 and 4:
46
;; >(makelist (+ (* i i) (* j j)) (i 1 4) (j 1 i))
47
;; (2 5 8 10 13 18 17 20 25 32)
48
;; >(makelist (list i j '---> (+ (* i i) (* j j))) (i 1 4) (j 1 i))
49
;; ((1 1 ---> 2) (2 1 ---> 5) (2 2 ---> 8) (3 1 ---> 10) (3 2 ---> 13)
50
;; (3 3 ---> 18) (4 1 ---> 17) (4 2 ---> 20) (4 3 ---> 25) (4 4 ---> 32))
52
;; Evaluate expression expr with variable set to lo, lo+1,... ,hi
53
;; and put the results in a list.
54
(defmacro makelist-1 (expr var lo hi &optional (step 1))
56
`(do ((,var ,lo (+ ,var ,step))
57
(,l nil (cons ,expr ,l)))
58
((> ,var ,hi) (reverse ,l))
59
(declare (fixnum ,var)))))
61
(defmacro makelist (expr (var lo hi &optional (step 1)) &rest more)
63
`(makelist-1 ,expr ,var ,lo ,hi ,step)
65
`(do ((,var ,lo (+ ,var ,step))
66
(,l nil (nconc ,l `,(makelist ,expr ,@more))))
68
(declare (fixnum ,var))))))
70
;;----------------------------------------------------------------
71
;; This package implements BASIC OPERATIONS ON MONOMIALS
72
;;----------------------------------------------------------------
73
;; DATA STRUCTURES: Monomials are represented as lists:
75
;; monom: (n1 n2 ... nk) where ni are non-negative integers
77
;; However, lists may be implemented as other sequence types,
78
;; so the flexibility to change the representation should be
79
;; maintained in the code to use general operations on sequences
80
;; whenever possible. The optimization for the actual representation
81
;; should be left to declarations and the compiler.
82
;;----------------------------------------------------------------
83
;; EXAMPLES: Suppose that variables are x and y. Then
85
;; Monom x*y^2 ---> (1 2)
87
;;----------------------------------------------------------------
90
"Type of exponent in a monomial."
93
(deftype monom (&optional dim)
95
`(simple-array exponent (,dim)))
97
(declaim (optimize (speed 3) (safety 0)))
99
(declaim (ftype (function (monom) fixnum) monom-dimension monom-sugar)
100
(ftype (function (monom &optional fixnum fixnum) fixnum) monom-total-degree)
101
(ftype (function (monom monom) monom) monom-div monom-mul monom-lcm monom-gcd)
102
(ftype (function (monom monom) (member t nil)) monom-divides-p monom-divisible-by-p monom-rel-prime-p)
103
(ftype (function (monom monom monom) (member t nil)) monom-divides-monom-lcm-p)
104
(ftype (function (monom monom monom monom) (member t nil)) monom-lcm-divides-monom-lcm-p)
105
(ftype (function (monom fixnum) (member t nil)) monom-depends-p)
106
;;(ftype (function (t monom &optional monom) monom) monom-map)
107
;;(ftype (function (monom monom) monom) monom-append)
110
(declaim (inline monom-mul monom-div
111
monom-total-degree monom-divides-p
112
monom-divisible-by-p monom-rel-prime monom-lcm))
114
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
116
;; Construction of monomials
118
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
120
(defmacro make-monom (dim &key (initial-contents nil initial-contents-supplied-p)
121
(initial-element 0 initial-element-supplied-p))
122
"Make a monomial with DIM variables. Additional argument
123
INITIAL-CONTENTS specifies the list of powers of the consecutive
124
variables. The alternative additional argument INITIAL-ELEMENT
125
specifies the common power for all variables."
126
(declare (fixnum dim))
128
:element-type 'exponent
129
,@(when initial-contents-supplied-p `(:initial-contents ,initial-contents))
130
,@(when initial-element-supplied-p `(:initial-element ,initial-element))))
133
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
135
;; Operations on monomials
137
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
139
(defmacro monom-elt (m index)
140
"Return the power in the monomial M of variable number INDEX."
143
(defun monom-dimension (m)
144
"Return the number of variables in the monomial M."
147
(defun monom-total-degree (m &optional (start 0) (end (length m)))
148
"Return the todal degree of a monomoal M. Optinally, a range
149
of variables may be specified with arguments START and END."
150
(declare (type monom m) (fixnum start end))
151
(reduce #'+ m :start start :end end))
153
(defun monom-sugar (m &aux (start 0) (end (length m)))
154
"Return the sugar of a monomial M. Optinally, a range
155
of variables may be specified with arguments START and END."
156
(declare (type monom m) (fixnum start end))
157
(monom-total-degree m start end))
159
(defun monom-div (m1 m2 &aux (result (copy-seq m1)))
160
"Divide monomial M1 by monomial M2."
161
(declare (type monom m1 m2 result))
162
(map-into result #'- m1 m2))
164
(defun monom-mul (m1 m2 &aux (result (copy-seq m1)))
165
"Multiply monomial M1 by monomial M2."
166
(declare (type monom m1 m2 result))
167
(map-into result #'+ m1 m2))
169
(defun monom-divides-p (m1 m2)
170
"Returns T if monomial M1 divides monomial M2, NIL otherwise."
171
(declare (type monom m1 m2))
174
(defun monom-divides-monom-lcm-p (m1 m2 m3)
175
"Returns T if monomial M1 divides MONOM-LCM(M2,M3), NIL otherwise."
176
(declare (type monom m1 m2 m3))
177
(every #'(lambda (x y z) (declare (type exponent x y z)) (<= x (max y z))) m1 m2 m3))
179
(defun monom-lcm-divides-monom-lcm-p (m1 m2 m3 m4)
180
"Returns T if monomial MONOM-LCM(M1,M2) divides MONOM-LCM(M3,M4), NIL otherwise."
181
(declare (type monom m1 m2 m3 m4))
182
(every #'(lambda (x y z w) (declare (type exponent x y z w)) (<= (max x y) (max z w))) m1 m2 m3 m4))
184
(defun monom-lcm-equal-monom-lcm-p (m1 m2 m3 m4)
185
"Returns T if monomial MONOM-LCM(M1,M2) equals MONOM-LCM(M3,M4), NIL otherwise."
186
(declare (type monom m1 m2 m3 m4))
187
(every #'(lambda (x y z w) (declare (type exponent x y z w)) (= (max x y) (max z w))) m1 m2 m3 m4))
189
(defun monom-divisible-by-p (m1 m2)
190
"Returns T if monomial M1 is divisible by monomial M2, NIL otherwise."
191
(declare (type monom m1 m2))
194
(defun monom-rel-prime-p (m1 m2)
195
"Returns T if two monomials M1 and M2 are relatively prime (disjoint)."
196
(declare (type monom m1 m2))
197
(every #'(lambda (x y) (declare (type exponent x y)) (zerop (min x y))) m1 m2))
199
(defun monom-equal-p (m1 m2)
200
"Returns T if two monomials M1 and M2 are equal."
201
(declare (type monom m1 m2))
204
(defun monom-lcm (m1 m2 &aux (result (copy-seq m1)))
205
"Returns least common multiple of monomials M1 and M2."
206
(declare (type monom m1 m2))
207
(map-into result #'max m1 m2))
209
(defun monom-gcd (m1 m2 &aux (result (copy-seq m1)))
210
"Returns greatest common divisor of monomials M1 and M2."
211
(declare (type monom m1 m2))
212
(map-into result #'min m1 m2))
214
(defun monom-depends-p (m k)
215
"Return T if the monomial M depends on variable number K."
216
(declare (type monom m) (fixnum k))
219
(defmacro monom-map (fun m &rest ml &aux (result `(copy-seq ,m)))
220
`(map-into ,result ,fun ,m ,@ml))
222
(defmacro monom-append (m1 m2)
223
`(concatenate 'monom ,m1 ,m2))
225
(defmacro monom-contract (k m)
228
(defun monom-exponents (m)
229
(declare (type monom m))
233
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
235
;; Implementations of various admissible monomial orders
237
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
239
;; pure lexicographic
240
(defun lex> (p q &optional (start 0) (end (monom-dimension p)))
241
"Return T if P>Q with respect to lexicographic order, otherwise NIL.
242
The second returned value is T if P=Q, otherwise it is NIL."
243
(declare (type monom p q) (type fixnum start end))
244
(do ((i start (1+ i)))
245
((>= i end) (values NIL T))
246
(declare (type fixnum i))
248
((> (monom-elt p i) (monom-elt q i))
249
(return-from lex> (values t nil)))
250
((< (monom-elt p i) (monom-elt q i))
251
(return-from lex> (values nil nil))))))
253
;; total degree order , ties broken by lexicographic
254
(defun grlex> (p q &optional (start 0) (end (monom-dimension p)))
255
"Return T if P>Q with respect to graded lexicographic order, otherwise NIL.
256
The second returned value is T if P=Q, otherwise it is NIL."
257
(declare (type monom p q) (type fixnum start end))
258
(let ((d1 (monom-total-degree p start end))
259
(d2 (monom-total-degree q start end)))
261
((> d1 d2) (values t nil))
262
((< d1 d2) (values nil nil))
264
(lex> p q start end)))))
267
;; total degree, ties broken by reverse lexicographic
268
(defun grevlex> (p q &optional (start 0) (end (monom-dimension p)))
269
"Return T if P>Q with respect to graded reverse lexicographic order,
270
NIL otherwise. The second returned value is T if P=Q, otherwise it is NIL."
271
(declare (type monom p q) (type fixnum start end))
272
(let ((d1 (monom-total-degree p start end))
273
(d2 (monom-total-degree q start end)))
275
((> d1 d2) (values t nil))
276
((< d1 d2) (values nil nil))
278
(revlex> p q start end)))))
281
;; reverse lexicographic
282
(defun revlex> (p q &optional (start 0) (end (monom-dimension p)))
283
"Return T if P>Q with respect to reverse lexicographic order, NIL
284
otherwise. The second returned value is T if P=Q, otherwise it is
285
NIL. This is not and admissible monomial order because some sets do
286
not have a minimal element. This order is useful in constructing other
288
(declare (type monom p q) (type fixnum start end))
289
(do ((i (1- end) (1- i)))
290
((< i start) (values NIL T))
291
(declare (type fixnum i))
293
((< (monom-elt p i) (monom-elt q i))
294
(return-from revlex> (values t nil)))
295
((> (monom-elt p i) (monom-elt q i))
296
(return-from revlex> (values nil nil))))))
299
(defun invlex> (p q &optional (start 0) (end (monom-dimension p)))
300
"Return T if P>Q with respect to inverse lexicographic order, NIL otherwise
301
The second returned value is T if P=Q, otherwise it is NIL."
302
(declare (type monom p q) (type fixnum start end))
303
(do ((i (1- end) (1- i)))
304
((< i start) (values NIL T))
305
(declare (type fixnum i))
307
((> (monom-elt p i) (monom-elt q i))
308
(return-from invlex> (values t nil)))
309
((< (monom-elt p i) (monom-elt q i))
310
(return-from invlex> (values nil nil))))))
313
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
315
;; Order making functions
317
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
319
(declaim (type function *monomial-order* *primary-elimination-order* *secondary-elimination-order*))
321
(defvar *monomial-order* #'lex>
322
"Default order for monomial comparisons")
324
(defmacro monomial-order (x y)
325
`(funcall *monomial-order* ,x ,y))
327
(defun reverse-monomial-order (x y)
328
(monomial-order y x))
330
(defvar *primary-elimination-order* #'lex>)
332
(defvar *secondary-elimination-order* #'lex>)
334
(defvar *elimination-order* nil
335
"Default elimination order used in elimination-based functions.
336
If not NIL, it is assumed to be a proper elimination order. If NIL,
337
we will construct an elimination order using the values of
338
*PRIMARY-ELIMINATION-ORDER* and *SECONDARY-ELIMINATION-ORDER*.")
340
(defun elimination-order (k)
341
"Return a predicate which compares monomials according to the
342
K-th elimination order. Two variables *PRIMARY-ELIMINATION-ORDER*
343
and *SECONDARY-ELIMINATION-ORDER* control the behavior on the first K
344
and the remaining variables, respectively."
345
(declare (type fixnum k))
346
#'(lambda (p q &optional (start 0) (end (monom-dimension p)))
347
(declare (type monom p q) (type fixnum start end))
348
(multiple-value-bind (primary equal)
349
(funcall *primary-elimination-order* p q start k)
351
(funcall *secondary-elimination-order* p q k end)
352
(values primary nil)))))
354
(defun elimination-order-1 (p q &optional (start 0) (end (monom-dimension p)))
355
"Equivalent to the function returned by the call to (ELIMINATION-ORDER 1)."
356
(declare (type monom p q) (type fixnum start end))
358
((> (monom-elt p start) (monom-elt q start)) (values t nil))
359
((< (monom-elt p start) (monom-elt q start)) (values nil nil))
360
(t (funcall *secondary-elimination-order* p q (1+ start) end))))
363
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
365
;; Priority queue stuff
367
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
369
(declaim (integer *priority-queue-allocation-size*))
371
(defparameter *priority-queue-allocation-size* 16)
373
(defun priority-queue-make-heap (&key (element-type 'fixnum))
374
(make-array *priority-queue-allocation-size* :element-type element-type :fill-pointer 1
377
(defstruct (priority-queue (:constructor priority-queue-construct))
378
(heap (priority-queue-make-heap))
381
(defun make-priority-queue (&key (element-type 'fixnum)
383
(element-key #'identity))
384
(priority-queue-construct
385
:heap (priority-queue-make-heap :element-type element-type)
386
:test #'(lambda (x y) (funcall test (funcall element-key y) (funcall element-key x)))))
388
(defun priority-queue-insert (pq item)
389
(priority-queue-heap-insert (priority-queue-heap pq) item (priority-queue-test pq)))
391
(defun priority-queue-remove (pq)
392
(priority-queue-heap-remove (priority-queue-heap pq) (priority-queue-test pq)))
394
(defun priority-queue-empty-p (pq)
395
(priority-queue-heap-empty-p (priority-queue-heap pq)))
397
(defun priority-queue-size (pq)
398
(fill-pointer (priority-queue-heap pq)))
400
(defun priority-queue-upheap (a k
405
(assert (< 0 k (fill-pointer a)))
407
(let ((parent (ash k -1)))
408
(when (zerop parent) (return))
409
(unless (funcall test (aref a parent) v)
411
(setf (aref a k) (aref a parent)
417
(defun priority-queue-heap-insert (a item &optional (test #'<=))
418
(vector-push-extend item a)
419
(priority-queue-upheap a (1- (fill-pointer a)) test))
421
(defun priority-queue-downheap (a k
424
&aux (v (aref a k)) (j 0) (n (fill-pointer a)))
425
(declare (fixnum k n j))
427
(unless (<= k (ash n -1))
430
(if (and (< j n) (not (funcall test (aref a (1+ j)) (aref a j))))
432
(when (funcall test (aref a j) v)
434
(setf (aref a k) (aref a j)
439
(defun priority-queue-heap-remove (a &optional (test #'<=) &aux (v (aref a 1)))
440
(when (<= (fill-pointer a) 1) (error "Empty queue."))
441
(setf (aref a 1) (vector-pop a))
442
(priority-queue-downheap a 1 test)
445
(defun priority-queue-heap-empty-p (a)
446
(<= (fill-pointer a) 1))
449
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
452
;; (Can be used in Maxima just fine)
454
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
456
(defmvar $poly_monomial_order '$lex
457
"This switch controls which monomial order is used in polynomial
458
and Grobner basis calculations. If not set, LEX will be used")
460
(defmvar $poly_coefficient_ring '$expression_ring
461
"This switch indicates the coefficient ring of the polynomials
462
that will be used in grobner calculations. If not set, Maxima's
463
general expression ring will be used. This variable may be set
464
to RING_OF_INTEGERS if desired.")
466
(defmvar $poly_primary_elimination_order NIL
467
"Name of the default order for eliminated variables in elimination-based functions.
468
If not set, LEX will be used.")
470
(defmvar $poly_secondary_elimination_order NIL
471
"Name of the default order for kept variables in elimination-based functions.
472
If not set, LEX will be used.")
474
(defmvar $poly_elimination_order NIL
475
"Name of the default elimination order used in elimination calculations.
476
If set, it overrides the settings in variables POLY_PRIMARY_ELIMINATION_ORDER
477
and SECONDARY_ELIMINATION_ORDER. The user must ensure that this is a true
478
elimination order valid for the number of eliminated variables.")
480
(defmvar $poly_return_term_list NIL
481
"If set to T, all functions in this package will return each polynomial as a
482
list of terms in the current monomial order rather than a Maxima general expression.")
484
(defmvar $poly_grobner_debug nil
485
"If set to TRUE, produce debugging and tracing output.")
487
(defmvar $poly_grobner_algorithm '$buchberger
488
"The name of the algorithm used to find grobner bases.")
490
(defmvar $poly_top_reduction_only NIL
491
"If not FALSE, use top reduction only whenever possible.
492
Top reduction means that division algorithm stops after the first reduction.")
495
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
497
;; Coefficient ring operations
499
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
501
;; These are ALL operations that are performed on the coefficients by
502
;; the package, and thus the coefficient ring can be changed by merely
503
;; redefining these operations.
505
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
508
(parse #'identity :type function)
509
(unit #'identity :type function)
510
(zerop #'identity :type function)
511
(add #'identity :type function)
512
(sub #'identity :type function)
513
(uminus #'identity :type function)
514
(mul #'identity :type function)
515
(div #'identity :type function)
516
(lcm #'identity :type function)
517
(ezgcd #'identity :type function)
518
(gcd #'identity :type function))
520
(declaim (type ring *RingOfIntegers* *FieldOfRationals*))
522
(defparameter *RingOfIntegers*
525
:unit #'(lambda () 1)
533
:ezgcd #'(lambda (x y &aux (c (gcd x y))) (values c (/ x c) (/ y c)))
535
"The ring of integers.")
538
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
540
;; This is how we perform operations on coefficients
541
;; using Maxima functions.
543
;; Functions and macros dealing with internal representation structure
545
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
548
(:constructor make-term (monom coeff))
549
(:constructor make-term-variable)
552
(monom (make-monom 0) :type monom)
555
(defun make-term-variable (ring nvars pos
558
(coeff (funcall (ring-unit ring)))
560
(monom (make-monom nvars :initial-element 0)))
561
(declare (fixnum nvars pos power))
562
(incf (monom-elt monom pos) power)
563
(make-term monom coeff))
565
(defun term-sugar (term)
566
(monom-sugar (term-monom term)))
568
(defun termlist-sugar (p &aux (sugar -1))
569
(declare (fixnum sugar))
570
(dolist (term p sugar)
571
(setf sugar (max sugar (term-sugar term)))))
575
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
577
;; Low-level polynomial arithmetic done on
580
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
582
(defmacro termlist-lt (p) `(car ,p))
583
(defun termlist-lm (p) (term-monom (termlist-lt p)))
584
(defun termlist-lc (p) (term-coeff (termlist-lt p)))
586
(define-modify-macro scalar-mul (c) coeff-mul)
588
(declaim (ftype (function (ring t t) t) scalar-times-termlist))
590
(defun scalar-times-termlist (ring c p)
591
"Multiply scalar C by a polynomial P. This function works
592
even if there are divisors of 0."
595
(let ((c1 (funcall (ring-mul ring) c (term-coeff term))))
596
(unless (funcall (ring-zerop ring) c1)
597
(list (make-term (term-monom term) c1)))))
601
(declaim (ftype (function (ring term term) list) term-mul))
603
(defun term-mul (ring term1 term2)
604
"Returns (LIST TERM) wheter TERM is the product of the terms TERM1 TERM2,
605
or NIL when the product is 0. This definition takes care of divisors of 0
606
in the coefficient ring."
607
(let ((c (funcall (ring-mul ring) (term-coeff term1) (term-coeff term2))))
608
(unless (funcall (ring-zerop ring) c)
609
(list (make-term (monom-mul (term-monom term1) (term-monom term2)) c)))))
611
(declaim (ftype (function (ring term list) list) term-times-termlist))
613
(defun term-times-termlist (ring term f)
614
(declare (type ring ring))
615
(mapcan #'(lambda (term-f) (term-mul ring term term-f)) f))
617
(declaim (ftype (function (ring list term) list) termlist-times-term))
619
(defun termlist-times-term (ring f term)
620
(mapcan #'(lambda (term-f) (term-mul ring term-f term)) f))
622
(declaim (ftype (function (monom term) term) monom-times-term))
624
(defun monom-times-term (m term)
625
(make-term (monom-mul m (term-monom term)) (term-coeff term)))
627
(declaim (ftype (function (monom list) list) monom-times-poly))
629
(defun monom-times-termlist (m f)
633
(mapcar #'(lambda (x) (monom-times-term m x)) f))))
635
(declaim (ftype (function (ring list) list) termlist-uminus))
637
(defun termlist-uminus (ring f)
638
(mapcar #'(lambda (x)
639
(make-term (term-monom x) (funcall (ring-uminus ring) (term-coeff x))))
642
(declaim (ftype (function (ring list list) list) termlist-add termlist-sub termlist-mul))
644
(defun termlist-add (ring p q)
645
(declare (type list p q))
649
(setf r (revappend r q)) t)
651
(setf r (revappend r p)) t)
654
(lm-greater lm-equal)
655
(monomial-order (termlist-lm p) (termlist-lm q))
658
(let ((s (funcall (ring-add ring) (termlist-lc p) (termlist-lc q))))
659
(unless (funcall (ring-zerop ring) s) ;check for cancellation
660
(setf r (cons (make-term (termlist-lm p) s) r)))
661
(setf p (cdr p) q (cdr q))))
663
(setf r (cons (car p) r)
665
(t (setf r (cons (car q) r)
670
(defun termlist-sub (ring p q)
671
(declare (type list p q))
675
(setf r (revappend r (termlist-uminus ring q)))
678
(setf r (revappend r p))
683
(monomial-order (termlist-lm p) (termlist-lm q))
686
(let ((s (funcall (ring-sub ring) (termlist-lc p) (termlist-lc q))))
687
(unless (funcall (ring-zerop ring) s) ;check for cancellation
688
(setf r (cons (make-term (termlist-lm p) s) r)))
689
(setf p (cdr p) q (cdr q))))
691
(setf r (cons (car p) r)
693
(t (setf r (cons (make-term (termlist-lm q) (funcall (ring-uminus ring) (termlist-lc q))) r)
698
;; Multiplication of polynomials
699
;; Non-destructive version
700
(defun termlist-mul (ring p q)
701
(cond ((or (endp p) (endp q)) nil) ;p or q is 0 (represented by NIL)
702
;; If p=p0+p1 and q=q0+q1 then pq=p0q0+p0q1+p1q
704
(term-times-termlist ring (car p) q))
706
(termlist-times-term ring p (car q)))
708
(let ((head (term-mul ring (termlist-lt p) (termlist-lt q)))
709
(tail (termlist-add ring (term-times-termlist ring (car p) (cdr q))
710
(termlist-mul ring (cdr p) q))))
711
(cond ((null head) tail)
713
(t (nconc head tail)))))))
715
(defun termlist-unit (ring dimension)
716
(declare (fixnum dimension))
717
(list (make-term (make-monom dimension :initial-element 0)
718
(funcall (ring-unit ring)))))
720
(defun termlist-expt (ring poly n &aux (dim (monom-dimension (termlist-lm poly))))
721
(declare (type fixnum n dim))
723
((minusp n) (error "termlist-expt: Negative exponent."))
724
((endp poly) (if (zerop n) (termlist-unit ring dim) nil))
727
(q poly (termlist-mul ring q q)) ;keep squaring
728
(p (termlist-unit ring dim) (if (not (zerop (logand k n))) (termlist-mul ring p q) p)))
730
(declare (fixnum k))))))
733
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
735
;; Additional structure operations on a list of terms
737
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
739
(defun termlist-contract (p &optional (k 1))
740
"Eliminate first K variables from a polynomial P."
741
(mapcar #'(lambda (term) (make-term (monom-contract k (term-monom term))
745
(defun termlist-extend (p &optional (m (list 0)))
746
"Extend every monomial in a polynomial P by inserting at the
747
beginning of every monomial the list of powers M."
748
(mapcar #'(lambda (term) (make-term (monom-append m (term-monom term))
752
(defun termlist-add-variables (p n)
753
"Add N variables to a polynomial P by inserting zero powers
754
at the beginning of each monomial."
756
(mapcar #'(lambda (term)
757
(make-term (monom-append (make-monom n :initial-element 0)
763
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
765
;; Arithmetic on polynomials
767
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
770
;;BOA constructor, by default constructs zero polynomial
771
(:constructor make-poly-from-termlist (termlist &optional (sugar (termlist-sugar termlist))))
772
(:constructor make-poly-zero (&aux (termlist nil) (sugar -1)))
773
;;Constructor of polynomials representing a variable
774
(:constructor make-variable (ring nvars pos &optional (power 1)
777
(make-term-variable ring nvars pos power)))
779
(:constructor poly-unit (ring dimension
781
(termlist (termlist-unit ring dimension))
783
(termlist nil :type list)
784
(sugar -1 :type fixnum))
787
(defmacro poly-lt (p) `(car (poly-termlist ,p)))
790
(defmacro poly-second-lt (p) `(cadar (poly-termlist ,p)))
793
(defun poly-lm (p) (term-monom (poly-lt p)))
796
(defun poly-second-lm (p) (term-monom (poly-second-lt p)))
798
;; Leading coefficient
799
(defun poly-lc (p) (term-coeff (poly-lt p)))
801
;; Second coefficient
802
(defun poly-second-lc (p) (term-coeff (poly-second-lt p)))
804
;; Testing for a zero polynomial
805
(defun poly-zerop (p) (null (poly-termlist p)))
807
;; The number of terms
808
(defun poly-length (p) (length (poly-termlist p)))
810
(declaim (ftype (function (ring t poly) poly) scalar-times-poly))
812
(defun scalar-times-poly (ring c p)
813
(make-poly-from-termlist (scalar-times-termlist ring c (poly-termlist p)) (poly-sugar p)))
815
(declaim (ftype (function (monom poly) poly) monom-times-poly))
817
(defun monom-times-poly (m p)
818
(make-poly-from-termlist (monom-times-termlist m (poly-termlist p)) (+ (poly-sugar p) (monom-sugar m))))
820
(declaim (ftype (function (ring term poly) poly) term-times-poly))
822
(defun term-times-poly (ring term p)
823
(make-poly-from-termlist (term-times-termlist ring term (poly-termlist p)) (+ (poly-sugar p) (term-sugar term))))
825
(declaim (ftype (function (ring poly poly) poly) poly-add poly-sub poly-mul))
827
(defun poly-add (ring p q)
828
(make-poly-from-termlist (termlist-add ring (poly-termlist p) (poly-termlist q)) (max (poly-sugar p) (poly-sugar q))))
830
(defun poly-sub (ring p q)
831
(make-poly-from-termlist (termlist-sub ring (poly-termlist p) (poly-termlist q)) (max (poly-sugar p) (poly-sugar q))))
833
(declaim (ftype (function (ring poly) poly) poly-uminus))
835
(defun poly-uminus (ring p)
836
(make-poly-from-termlist (termlist-uminus ring (poly-termlist p)) (poly-sugar p)))
838
(defun poly-mul (ring p q)
839
(make-poly-from-termlist (termlist-mul ring (poly-termlist p) (poly-termlist q)) (+ (poly-sugar p) (poly-sugar q))))
841
(declaim (ftype (function (ring poly fixnum) poly) poly-expt))
843
(defun poly-expt (ring p n)
844
(make-poly-from-termlist (termlist-expt ring (poly-termlist p) n) (* n (poly-sugar p))))
846
(defun poly-append (&rest plist)
847
(make-poly-from-termlist (apply #'append (mapcar #'poly-termlist plist))
848
(apply #'max (mapcar #'poly-sugar plist))))
850
(declaim (ftype (function (poly) poly) poly-nreverse))
852
(defun poly-nreverse (p)
853
(setf (poly-termlist p) (nreverse (poly-termlist p)))
856
(declaim (ftype (function (poly &optional fixnum) poly) poly-contract))
858
(defun poly-contract (p &optional (k 1))
859
(make-poly-from-termlist (termlist-contract (poly-termlist p) k)
862
(declaim (ftype (function (poly &optional sequence)) poly-extend))
864
(defun poly-extend (p &optional (m (list 0)))
865
(make-poly-from-termlist
866
(termlist-extend (poly-termlist p) m)
867
(+ (poly-sugar p) (monom-sugar m))))
869
(declaim (ftype (function (poly fixnum)) poly-add-variables))
871
(defun poly-add-variables (p k)
872
(setf (poly-termlist p) (termlist-add-variables (poly-termlist p) k))
875
(defun poly-list-add-variables (plist k)
876
(mapcar #'(lambda (p) (poly-add-variables p k)) plist))
878
(defun poly-standard-extension (plist &aux (k (length plist)))
879
"Calculate [U1*P1,U2*P2,...,UK*PK], where PLIST=[P1,P2,...,PK]."
880
(declare (list plist) (fixnum k))
881
(labels ((incf-power (g i)
882
(dolist (x (poly-termlist g))
883
(incf (monom-elt (term-monom x) i)))
884
(incf (poly-sugar g))))
885
(setf plist (poly-list-add-variables plist k))
887
(incf-power (nth i plist) i))))
889
(defun saturation-extension (ring F plist &aux (k (length plist)) (d (monom-dimension (poly-lm (car plist)))))
890
"Calculate [F, U1*P1-1,U2*P2-1,...,UK*PK-1], where PLIST=[P1,P2,...,PK]."
891
(setf F (poly-list-add-variables F k)
892
plist (mapcar #'(lambda (x)
893
(setf (poly-termlist x) (nconc (poly-termlist x)
894
(list (make-term (make-monom d :initial-element 0)
895
(funcall (ring-uminus ring) (funcall (ring-unit ring)))))))
897
(poly-standard-extension plist)))
901
(defun polysaturation-extension (ring F plist &aux (k (length plist))
902
(d (+ k (length (poly-lm (car plist))))))
903
"Calculate [F, U1*P1+U2*P2+...+UK*PK-1], where PLIST=[P1,P2,...,PK]."
904
(setf F (poly-list-add-variables F k)
905
plist (apply #'poly-append (poly-standard-extension plist))
906
(cdr (last (poly-termlist plist))) (list (make-term (make-monom d :initial-element 0)
907
(funcall (ring-uminus ring) (funcall (ring-unit ring))))))
908
(append F (list plist)))
910
(defun saturation-extension-1 (ring F p) (polysaturation-extension ring F (list p)))
914
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
916
;; Evaluation of polynomial (prefix) expressions
918
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
920
(defun coerce-coeff (ring expr vars)
921
"Coerce an element of the coefficient ring to a constant polynomial."
922
;; Modular arithmetic handler by rat
923
(make-poly-from-termlist (list (make-term (make-monom (length vars) :initial-element 0)
924
(funcall (ring-parse ring) expr)))
927
(defun poly-eval (ring expr vars &optional (list-marker '[))
928
(labels ((p-eval (arg) (poly-eval ring arg vars))
929
(p-eval-list (args) (mapcar #'p-eval args))
930
(p-add (x y) (poly-add ring x y)))
932
((eql expr 0) (make-poly-zero))
933
((member expr vars :test #'equalp)
934
(let ((pos (position expr vars :test #'equalp)))
935
(make-variable ring (length vars) pos)))
937
(coerce-coeff ring expr vars))
938
((eq (car expr) list-marker)
939
(cons list-marker (p-eval-list (cdr expr))))
942
(+ (reduce #'p-add (p-eval-list (cdr expr))))
943
(- (case (length expr)
945
(2 (poly-uminus ring (p-eval (cadr expr))))
946
(3 (poly-sub ring (p-eval (cadr expr)) (p-eval (caddr expr))))
947
(otherwise (poly-sub ring (p-eval (cadr expr))
948
(reduce #'p-add (p-eval-list (cddr expr)))))))
950
(if (endp (cddr expr)) ;unary
952
(reduce #'(lambda (p q) (poly-mul ring p q)) (p-eval-list (cdr expr)))))
955
((member (cadr expr) vars :test #'equalp)
956
;;Special handling of (expt var pow)
957
(let ((pos (position (cadr expr) vars :test #'equalp)))
958
(make-variable ring (length vars) pos (caddr expr))))
959
((not (and (integerp (caddr expr)) (plusp (caddr expr))))
960
;; Negative power means division in coefficient ring
961
;; Non-integer power means non-polynomial coefficient
962
(coerce-coeff ring expr vars))
963
(t (poly-expt ring (p-eval (cadr expr)) (caddr expr)))))
965
(coerce-coeff ring expr vars)))))))
968
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
970
;; Global optimization/debugging options
972
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
975
;;All inline functions of this module
976
(declaim (inline free-of-vars make-pair-queue pair-queue-insert
977
pair-queue-remove pair-queue-empty-p
978
pair-queue-remove pair-queue-size Criterion-1
979
Criterion-2 grobner reduced-grobner sugar-pair-key
980
sugar-order normal-form normal-form-step grobner-op spoly
984
;;Optimization options
985
(declaim (optimize (speed 3) (safety 0)))
988
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
992
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
996
(defmacro debug-cgb (&rest args)
997
`(when $poly_grobner_debug (format *terminal-io* ,@args)))
999
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1001
;; An implementation of Grobner basis
1003
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1005
(defun spoly (ring f g)
1006
"It yields the S-polynomial of polynomials F and G."
1007
(declare (type poly f g))
1008
(let* ((lcm (monom-lcm (poly-lm f) (poly-lm g)))
1009
(mf (monom-div lcm (poly-lm f)))
1010
(mg (monom-div lcm (poly-lm g))))
1011
(declare (type monom mf mg))
1012
(multiple-value-bind (c cf cg)
1013
(funcall (ring-ezgcd ring) (poly-lc f) (poly-lc g))
1014
(declare (ignore c))
1017
(scalar-times-poly ring cg (monom-times-poly mf f))
1018
(scalar-times-poly ring cf (monom-times-poly mg g))))))
1021
(defun poly-primitive-part (ring p)
1022
"Divide polynomial P with integer coefficients by gcd of its
1023
coefficients and return the result."
1024
(declare (type poly p))
1027
(let ((c (poly-content ring p)))
1028
(values (make-poly-from-termlist (mapcar
1030
(make-term (term-monom x)
1031
(funcall (ring-div ring) (term-coeff x) c)))
1036
(defun poly-content (ring p)
1037
"Greatest common divisor of the coefficients of the polynomial P. Use the RING structure
1038
to compute the greatest common divisor."
1039
(declare (type poly p))
1040
(reduce (ring-gcd ring) (mapcar #'term-coeff (rest (poly-termlist p))) :initial-value (poly-lc p)))
1043
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1045
;; An implementation of the division algorithm
1047
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1049
(declaim (ftype (function (ring t t monom poly poly) poly) grobner-op))
1051
(defun grobner-op (ring c1 c2 m f g)
1052
"Returns C2*F-C1*M*G, where F and G are polynomials M is a monomial.
1053
Assume that the leading terms will cancel."
1054
#+grobner-check(funcall (ring-zerop ring)
1055
(funcall (ring-sub ring)
1056
(funcall (ring-mul ring) c2 (poly-lc f))
1057
(funcall (ring-mul ring) c1 (poly-lc g))))
1058
#+grobner-check(monom-equal-p (poly-lm f) (monom-mul m (poly-lm g)))
1060
(scalar-times-poly ring c2 f)
1061
(scalar-times-poly ring c1 (monom-times-poly m g))))
1063
(defun poly-pseudo-divide (ring f fl)
1064
"Pseudo-divide a polynomial F by the list of polynomials FL. Return
1065
multiple values. The first value is a list of quotients A. The second
1066
value is the remainder R. The third argument is a scalar coefficient
1067
C, such that C*F can be divided by FL within the ring of coefficients,
1068
which is not necessarily a field. Finally, the fourth value is an
1069
integer count of the number of reductions performed. The resulting
1070
objects satisfy the equation: C*F= sum A[i]*FL[i] + R."
1071
(declare (type poly f) (list fl))
1072
(do ((r (make-poly-zero))
1073
(c (funcall (ring-unit ring)))
1074
(a (make-list (length fl) :initial-element (make-poly-zero)))
1078
(debug-cgb "~&~3T~d reduction~:p" division-count)
1079
(when (poly-zerop r) (debug-cgb " ---> 0"))
1080
(values (mapcar #'poly-nreverse a) (poly-nreverse r) c division-count))
1081
(declare (fixnum division-count))
1082
(do ((fl fl (rest fl)) ;scan list of divisors
1085
((endp fl) ;no division occurred
1086
(push (poly-lt p) (poly-termlist r)) ;move lt(p) to remainder
1087
(setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
1088
(pop (poly-termlist p)) ;remove lt(p) from p
1090
((monom-divides-p (poly-lm (car fl)) (poly-lm p)) ;division occurred
1091
(incf division-count)
1092
(multiple-value-bind (gcd c1 c2)
1093
(funcall (ring-ezgcd ring) (poly-lc (car fl)) (poly-lc p))
1094
(declare (ignore gcd))
1095
(let ((m (monom-div (poly-lm p) (poly-lm (car fl)))))
1096
;; Multiply the equation c*f=sum ai*fi+r+p by c1.
1098
(setf (car x) (scalar-times-poly ring c1 (car x))))
1100
(setf r (scalar-times-poly ring c1 r)
1101
c (funcall (ring-mul ring) c c1)
1102
p (grobner-op ring c2 c1 m p (car fl)))
1103
(push (make-term m c2) (poly-termlist (car b))))
1106
(defun poly-exact-divide (ring f g)
1107
"Divide a polynomial F by another polynomial G. Assume that exact division
1108
with no remainder is possible. Returns the quotient."
1109
(declare (type poly f g))
1110
(multiple-value-bind (quot rem coeff division-count)
1111
(poly-pseudo-divide ring f (list g))
1112
(declare (ignore division-count coeff)
1113
(type poly quot rem) (type fixnum division-count))
1114
(unless (poly-zerop rem) (error "Exact division failed."))
1118
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1120
;; An implementation of the normal form
1122
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1124
(declaim (ftype (function (ring t poly poly t fixnum)
1125
(values poly poly t fixnum))
1128
(defun normal-form-step (ring fl p r c division-count
1129
&aux (g (find (poly-lm p) fl
1130
:test #'monom-divisible-by-p
1133
(g ;division possible
1134
(incf division-count)
1135
(multiple-value-bind (gcd cg cp)
1136
(funcall (ring-ezgcd ring) (poly-lc g) (poly-lc p))
1137
(declare (ignore gcd))
1138
(let ((m (monom-div (poly-lm p) (poly-lm g))))
1139
;; Multiply the equation c*f=sum ai*fi+r+p by cg.
1140
(setf r (scalar-times-poly ring cg r)
1141
c (funcall (ring-mul ring) c cg)
1142
p (grobner-op ring cp cg m p g))))
1144
(t ;no division possible
1145
(push (poly-lt p) (poly-termlist r)) ;move lt(p) to remainder
1146
(setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
1147
(pop (poly-termlist p)) ;remove lt(p) from p
1149
(values p r c division-count))
1151
(declaim (ftype (function (ring poly t &optional t) (values poly t fixnum)) normal-form))
1153
;; Merge it sometime with poly-pseudo-divide
1154
(defun normal-form (ring f fl &optional (top-reduction-only $poly_top_reduction_only))
1155
;; Loop invariant: c*f0=sum ai*fi+r+f, where f0 is the initial value of f
1156
#+grobner-check(when (null fl) (warn "normal-form: empty divisor list."))
1157
(do ((r (make-poly-zero))
1158
(c (funcall (ring-unit ring)))
1162
(and top-reduction-only (not (poly-zerop r))))
1164
(debug-cgb "~&~3T~d reduction~:p" division-count)
1165
(when (poly-zerop r)
1166
(debug-cgb " ---> 0")))
1167
(setf (poly-termlist f) (nreconc (poly-termlist r) (poly-termlist f)))
1168
(values f c division-count))
1169
(declare (fixnum division-count)
1171
(multiple-value-setq (f r c division-count)
1172
(normal-form-step ring fl f r c division-count))))
1175
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1177
;; These are provided mostly for debugging purposes To enable
1178
;; verification of grobner bases with BUCHBERGER-CRITERION, do
1179
;; (pushnew :grobner-check *features*) and compile/load this file.
1180
;; With this feature, the calculations will slow down CONSIDERABLY.
1182
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1184
(defun buchberger-criterion (ring G)
1185
"Returns T if G is a Grobner basis, by using the Buchberger
1186
criterion: for every two polynomials h1 and h2 in G the S-polynomial
1187
S(h1,h2) reduces to 0 modulo G."
1190
(makelist (normal-form ring (spoly ring (elt G i) (elt G j)) G nil)
1191
(i 0 (- (length G) 2))
1192
(j (1+ i) (1- (length G))))))
1194
(defun grobner-test (ring G F)
1195
"Test whether G is a Grobner basis and F is contained in G. Return T
1196
upon success and NIL otherwise."
1197
(debug-cgb "~&GROBNER CHECK: ")
1198
(let (($poly_grobner_debug nil)
1199
(stat1 (buchberger-criterion ring G))
1202
(makelist (normal-form ring (copy-tree (elt F i)) G nil)
1203
(i 0 (1- (length F)))))))
1204
(unless stat1 (error "~&Buchberger criterion failed."))
1206
(error "~&Original polys not in ideal spanned by Grobner.")))
1207
(debug-cgb "~&GROBNER CHECK END")
1211
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1213
;; Pair queue implementation
1215
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1217
(defun sugar-pair-key (p q &aux (lcm (monom-lcm (poly-lm p) (poly-lm q)))
1218
(d (monom-sugar lcm)))
1219
"Returns list (S LCM-TOTAL-DEGREE) where S is the sugar of the S-polynomial of
1220
polynomials P and Q, and LCM-TOTAL-DEGREE is the degree of is LCM(LM(P),LM(Q))."
1221
(declare (type poly p q) (type monom lcm) (type fixnum d))
1223
(+ (- d (monom-sugar (poly-lm p))) (poly-sugar p))
1224
(+ (- d (monom-sugar (poly-lm q))) (poly-sugar q)))
1228
(:constructor make-pair (first second
1230
(sugar (car (sugar-pair-key first second)))
1231
(division-data nil))))
1232
(first nil :type poly)
1233
(second nil :type poly)
1234
(sugar 0 :type fixnum)
1235
(division-data nil :type list))
1237
;;(defun pair-sugar (pair &aux (p (pair-first pair)) (q (pair-second pair)))
1238
;; (car (sugar-pair-key p q)))
1240
(defun sugar-order (x y)
1241
"Pair order based on sugar, ties broken by normal strategy."
1242
(declare (type cons x y))
1243
(or (< (car x) (car y))
1244
(and (= (car x) (car y))
1245
(< (monom-total-degree (cdr x))
1246
(monom-total-degree (cdr y))))))
1248
(defvar *pair-key-function* #'sugar-pair-key
1249
"Function that, given two polynomials as argument, computed the key
1250
in the pair queue.")
1252
(defvar *pair-order* #'sugar-order
1253
"Function that orders the keys of pairs.")
1255
(defun make-pair-queue ()
1256
"Constructs a priority queue for critical pairs."
1257
(make-priority-queue
1259
:element-key #'(lambda (pair) (funcall *pair-key-function* (pair-first pair) (pair-second pair)))
1260
:test *pair-order*))
1262
(defun pair-queue-initialize (pq F start
1265
(B (nconc (makelist (make-pair (elt F i) (elt F j))
1266
(i 0 (1- start)) (j start s))
1267
(makelist (make-pair (elt F i) (elt F j))
1268
(i start (1- s)) (j (1+ i) s)))))
1269
"Initializes the priority for critical pairs. F is the initial list of polynomials.
1270
START is the first position beyond the elements which form a partial
1271
grobner basis, i.e. satisfy the Buchberger criterion."
1272
(declare (type priority-queue pq) (type fixnum start))
1274
(priority-queue-insert pq pair)))
1276
(defun pair-queue-insert (B pair)
1277
(priority-queue-insert B pair))
1279
(defun pair-queue-remove (B)
1280
(priority-queue-remove B))
1282
(defun pair-queue-size (B)
1283
(priority-queue-size B))
1285
(defun pair-queue-empty-p (B)
1286
(priority-queue-empty-p B))
1288
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1290
;; Buchberger Algorithm Implementation
1292
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1294
(defun buchberger (ring F start &optional (top-reduction-only $poly_top_reduction_only)
1296
"An implementation of the Buchberger algorithm. Return Grobner basis
1297
of the ideal generated by the polynomial list F. Polynomials 0 to
1298
START-1 are assumed to be a Grobner basis already, so that certain
1299
critical pairs will not be examined. If TOP-REDUCTION-ONLY set, top
1300
reduction will be preformed. This function assumes that all polynomials
1302
(declare (type fixnum start) (type priority-queue B) (type hash-table B-done))
1303
(when (endp F) (return-from buchberger F)) ;cut startup costs
1304
(debug-cgb "~&GROBNER BASIS - BUCHBERGER ALGORITHM")
1305
(when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
1306
#+grobner-check (when (plusp start)
1307
(grobner-test ring (subseq F 0 start) (subseq F 0 start)))
1308
;;Initialize critical pairs
1309
(setf B (pair-queue-initialize (make-pair-queue)
1311
B-done (make-hash-table :test #'equal))
1312
(dotimes (i (1- start))
1313
(do ((j (1+ i) (1+ j))) ((>= j start))
1314
(setf (gethash (list (elt F i) (elt F j)) B-done) t)))
1316
((pair-queue-empty-p B)
1317
#+grobner-check(grobner-test ring F F)
1318
(debug-cgb "~&GROBNER END")
1320
(let ((pair (pair-queue-remove B)))
1321
(declare (type pair pair))
1323
((Criterion-1 pair) nil)
1324
((Criterion-2 pair B-done F) nil)
1326
(let ((SP (normal-form ring
1327
(spoly ring (pair-first pair) (pair-second pair))
1328
F top-reduction-only)))
1329
(declare (type poly SP))
1334
(setf SP (poly-primitive-part ring SP)
1335
F (nconc F (list SP)))
1336
;; Add new critical pairs
1338
(pair-queue-insert B (make-pair h SP)))
1339
(debug-cgb "~&Sugar: ~d Polynomials: ~d; Pairs left: ~d; Pairs done: ~d;"
1340
(pair-sugar pair) (length F) (pair-queue-size B)
1341
(hash-table-count B-done)))))))
1342
(setf (gethash (list (pair-first pair) (pair-second pair)) B-done) t))))
1344
(defun parallel-buchberger (ring F start &optional (top-reduction-only $poly_top_reduction_only)
1346
"An implementation of the Buchberger algorithm. Return Grobner basis
1347
of the ideal generated by the polynomial list F. Polynomials 0 to
1348
START-1 are assumed to be a Grobner basis already, so that certain
1349
critical pairs will not be examined. If TOP-REDUCTION-ONLY set, top
1350
reduction will be preformed."
1351
(declare (ignore top-reduction-only)
1353
(type priority-queue B)
1354
(type hash-table B-done))
1355
(when (endp F) (return-from parallel-buchberger F)) ;cut startup costs
1356
(debug-cgb "~&GROBNER BASIS - PARALLEL-BUCHBERGER ALGORITHM")
1357
(when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
1358
#+grobner-check (when (plusp start)
1359
(grobner-test ring (subseq F 0 start) (subseq F 0 start)))
1360
;;Initialize critical pairs
1361
(setf B (pair-queue-initialize (make-pair-queue)
1363
B-done (make-hash-table :test #'equal))
1364
(dotimes (i (1- start))
1365
(do ((j (1+ i) (1+ j))) ((>= j start))
1366
(declare (type fixnum j))
1367
(setf (gethash (list (elt F i) (elt F j)) B-done) t)))
1369
((pair-queue-empty-p B)
1370
#+grobner-check(grobner-test ring F F)
1371
(debug-cgb "~&GROBNER END")
1373
(let ((pair (pair-queue-remove B)))
1374
(when (null (pair-division-data pair))
1375
(setf (pair-division-data pair) (list (spoly ring
1379
(funcall (ring-unit ring))
1382
((Criterion-1 pair) nil)
1383
((Criterion-2 pair B-done F) nil)
1385
(let* ((dd (pair-division-data pair))
1389
(division-count (fourth dd)))
1391
((poly-zerop p) ;normal form completed
1392
(debug-cgb "~&~3T~d reduction~:p" division-count)
1395
(debug-cgb " ---> 0")
1398
(setf SP (poly-nreverse SP)
1399
SP (poly-primitive-part ring SP)
1400
F (nconc F (list SP)))
1401
;; Add new critical pairs
1403
(pair-queue-insert B (make-pair h SP)))
1404
(debug-cgb "~&Sugar: ~d Polynomials: ~d; Pairs left: ~d; Pairs done: ~d;"
1405
(pair-sugar pair) (length F) (pair-queue-size B)
1406
(hash-table-count B-done))))
1407
(setf (gethash (list (pair-first pair) (pair-second pair)) B-done) t))
1408
(t ;normal form not complete
1411
((> (poly-sugar SP) (pair-sugar pair))
1412
(debug-cgb "(~a)?" (poly-sugar SP))
1421
(fourth dd) division-count
1422
(pair-sugar pair) (poly-sugar SP))
1423
(pair-queue-insert B pair))
1424
(multiple-value-setq (p SP c division-count)
1425
(normal-form-step ring F p SP c division-count)))))))))))
1428
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1432
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1434
(defun Criterion-1 (pair &aux (f (pair-first pair)) (g (pair-second pair)))
1435
"Returns T if the leading monomials of the two polynomials
1436
in G pointed to by the integers in PAIR have disjoint (relatively prime)
1437
monomials. This test is known as the first Buchberger criterion."
1438
(declare (type pair pair))
1439
(when (monom-rel-prime-p (poly-lm f) (poly-lm g))
1441
(return-from Criterion-1 t)))
1443
(defun Criterion-2 (pair B-done partial-basis
1444
&aux (f (pair-first pair)) (g (pair-second pair))
1446
"Returns T if the leading monomial of some element P of
1447
PARTIAL-BASIS divides the LCM of the leading monomials of the two
1448
polynomials in the polynomial list PARTIAL-BASIS, and P paired with
1449
each of the polynomials pointed to by the the PAIR has already been
1450
treated, as indicated by the absence in the hash table B-done."
1451
(declare (type pair pair) (type hash-table B-done)
1453
;; In the code below we assume that pairs are ordered as follows:
1454
;; if PAIR is (I J) then I appears before J in the PARTIAL-BASIS.
1455
;; We traverse the list PARTIAL-BASIS and keep track of where we
1456
;; are, so that we can produce the pairs in the correct order
1457
;; when we check whether they have been processed, i.e they
1458
;; appear in the hash table B-done
1459
(dolist (h partial-basis nil)
1462
#+grobner-check(assert (eq place :before))
1463
(setf place :in-the-middle))
1465
#+grobner-check(assert (eq place :in-the-middle))
1466
(setf place :after))
1467
((and (monom-divides-monom-lcm-p (poly-lm h) (poly-lm f) (poly-lm g))
1468
(gethash (case place
1469
(:before (list h f))
1470
((:in-the-middle :after) (list f h)))
1472
(gethash (case place
1473
((:before :in-the-middle) (list h g))
1474
(:after (list g h)))
1477
(return-from Criterion-2 t)))))
1480
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1482
;; An implementation of the algorithm of Gebauer and Moeller, as
1483
;; described in the book of Becker-Weispfenning, p. 232
1485
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1487
(defun gebauer-moeller (ring F start &optional (top-reduction-only $poly_top_reduction_only)
1489
"Compute Grobner basis by using the algorithm of Gebauer and
1490
Moeller. This algorithm is described as BUCHBERGERNEW2 in the book by
1491
Becker-Weispfenning entitled ``Grobner Bases''. This function assumes
1492
that all polynomials in F are non-zero."
1493
(declare (ignore top-reduction-only)
1495
(type priority-queue B))
1497
((endp F) (return-from gebauer-moeller nil))
1499
(return-from gebauer-moeller (list (poly-primitive-part ring (car F))))))
1500
(debug-cgb "~&GROBNER BASIS - GEBAUER MOELLER ALGORITHM")
1501
(when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
1502
#+grobner-check (when (plusp start)
1503
(grobner-test ring (subseq F 0 start) (subseq F 0 start)))
1504
(setf B (make-pair-queue)
1505
G (subseq F 0 start)
1506
F1 (subseq F start))
1508
(multiple-value-setq (G B)
1509
(gebauer-moeller-update G B (poly-primitive-part ring (pop F1)))))
1510
(do () ((pair-queue-empty-p B))
1511
(let* ((pair (pair-queue-remove B))
1512
(g1 (pair-first pair))
1513
(g2 (pair-second pair))
1514
(h (normal-form ring (spoly ring g1 g2)
1516
nil #| Always fully reduce! |#
1518
(unless (poly-zerop h)
1519
(setf h (poly-primitive-part ring h))
1520
(multiple-value-setq (G B)
1521
(gebauer-moeller-update G B h))
1522
(debug-cgb "~&Sugar: ~d Polynomials: ~d; Pairs left: ~d~%"
1523
(pair-sugar pair) (length G) (pair-queue-size B))
1525
#+grobner-check(grobner-test ring G F)
1526
(debug-cgb "~&GROBNER END")
1529
(defun gebauer-moeller-update (G B h
1532
(B-new (make-pair-queue))
1535
"An implementation of the auxillary UPDATE algorithm used by the
1536
Gebauer-Moeller algorithm. G is a list of polynomials, B is a list of
1537
critical pairs and H is a new polynomial which possibly will be added
1538
to G. The naming conventions used are very close to the one used in
1539
the book of Becker-Weispfenning."
1541
#+allegro (dynamic-extent B)
1543
(type priority-queue B)
1547
(declare (type poly g1))
1549
(when (or (monom-rel-prime-p (poly-lm h) (poly-lm g1))
1551
(notany #'(lambda (g2) (monom-lcm-divides-monom-lcm-p
1552
(poly-lm h) (poly-lm g2)
1553
(poly-lm h) (poly-lm g1)))
1555
(notany #'(lambda (g2) (monom-lcm-divides-monom-lcm-p
1556
(poly-lm h) (poly-lm g2)
1557
(poly-lm h) (poly-lm g1)))
1562
(declare (type poly g1))
1564
(unless (monom-rel-prime-p (poly-lm h) (poly-lm g1))
1566
(do (g1 g2) ((pair-queue-empty-p B))
1567
(declare (type poly g1 g2))
1568
(setf pair (pair-queue-remove B)
1569
g1 (pair-first pair)
1570
g2 (pair-second pair))
1571
(when (or (not (monom-divides-monom-lcm-p
1573
(poly-lm g1) (poly-lm g2)))
1574
(monom-lcm-equal-monom-lcm-p
1575
(poly-lm g1) (poly-lm h)
1576
(poly-lm g1) (poly-lm g2))
1577
(monom-lcm-equal-monom-lcm-p
1578
(poly-lm h) (poly-lm g2)
1579
(poly-lm g1) (poly-lm g2)))
1580
(pair-queue-insert B-new (make-pair g1 g2))))
1582
(pair-queue-insert B-new (make-pair h g3)))
1585
(declare (type poly g1))
1587
(unless (monom-divides-p (poly-lm h) (poly-lm g1))
1590
(values G-new B-new))
1593
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1595
;; Standard postprocessing of Grobner bases
1597
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1599
(defun reduction (ring plist)
1600
"Reduce a list of polynomials PLIST, so that non of the terms in any of
1601
the polynomials is divisible by a leading monomial of another
1602
polynomial. Return the reduced list."
1606
(mapcar #'(lambda (x) (poly-primitive-part ring x)) Q))
1607
;;Find p in Q such that p is reducible mod Q\{p}
1610
(let ((Q1 (remove x Q)))
1611
(multiple-value-bind (h c div-count)
1612
(normal-form ring x Q1 nil #| not a top reduction! |# )
1613
(declare (ignore c))
1614
(unless (zerop div-count)
1616
(unless (poly-zerop h)
1617
(setf Q (nconc Q1 (list h))))
1620
(defun minimization (P)
1621
"Returns a sublist of the polynomial list P spanning the same
1622
monomial ideal as P but minimal, i.e. no leading monomial
1623
of a polynomial in the sublist divides the leading monomial
1624
of another polynomial."
1628
;;Find p in Q such that lm(p) is in LM(Q\{p})
1631
(let ((Q1 (remove x Q)))
1632
(when (member-if #'(lambda (p) (monom-divides-p (poly-lm x) (poly-lm p))) Q1)
1636
(defun poly-normalize (ring p &aux (c (poly-lc p)))
1637
"Divide a polynomial by its leading coefficient. It assumes
1638
that the division is possible, which may not always be the
1639
case in rings which are not fields. The exact division operator
1640
is assumed to be provided by the RING structure of the
1641
COEFFICIENT-RING package."
1642
(mapc #'(lambda (term)
1643
(setf (term-coeff term) (funcall (ring-div ring) (term-coeff term) c)))
1647
(defun poly-normalize-list (ring plist)
1648
"Divide every polynomial in a list PLIST by its leading coefficient. "
1649
(mapcar #'(lambda (x) (poly-normalize ring x)) plist))
1652
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1654
;; Algorithm and Pair heuristic selection
1656
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1658
(defun find-grobner-function (algorithm)
1659
"Return a function which calculates Grobner basis, based on its
1660
names. Names currently used are either Lisp symbols, Maxima symbols or
1663
((buchberger :buchberger $buchberger) #'buchberger)
1664
((parallel-buchberger :parallel-buchberger $parallel_buchberger) #'parallel-buchberger)
1665
((gebauer-moeller :gebauer_moeller $gebauer_moeller) #'gebauer-moeller)))
1667
(defun grobner (ring F &optional (start 0) (top-reduction-only nil))
1668
;;(setf F (sort F #'< :key #'sugar))
1670
(find-grobner-function $poly_grobner_algorithm)
1671
ring F start top-reduction-only))
1673
(defun reduced-grobner (ring F &optional (start 0) (top-reduction-only $poly_top_reduction_only))
1674
(reduction ring (grobner ring F start top-reduction-only)))
1676
(defun set-pair-heuristic (method)
1677
"Sets up variables *PAIR-KEY-FUNCTION* and *PAIR-ORDER* used
1678
to determine the priority of critical pairs in the priority queue."
1680
((sugar :sugar $sugar)
1681
(setf *pair-key-function* #'sugar-pair-key
1682
*pair-order* #'sugar-order))
1683
; ((minimal-mock-spoly :minimal-mock-spoly $minimal_mock_spoly)
1684
; (setf *pair-key-function* #'mock-spoly
1685
; *pair-order* #'mock-spoly-order))
1686
((minimal-lcm :minimal-lcm $minimal_lcm)
1687
(setf *pair-key-function* #'(lambda (p q)
1688
(monom-lcm (poly-lm p) (poly-lm q)))
1689
*pair-order* #'reverse-monomial-order))
1690
((minimal-total-degree :minimal-total-degree $minimal_total_degree)
1691
(setf *pair-key-function* #'(lambda (p q)
1693
(monom-lcm (poly-lm p) (poly-lm q))))
1695
((minimal-length :minimal-length $minimal_length)
1696
(setf *pair-key-function* #'(lambda (p q)
1697
(+ (poly-length p) (poly-length q)))
1698
*pair-order* #'<))))
1701
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1703
;; Operations in ideal theory
1705
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1707
;; Does the term depend on variable K?
1708
(defun term-depends-p (term k)
1709
"Return T if the term TERM depends on variable number K."
1710
(monom-depends-p (term-monom term) k))
1712
;; Does the polynomial P depend on variable K?
1713
(defun poly-depends-p (p k)
1714
"Return T if the term polynomial P depends on variable number K."
1715
(some #'(lambda (term) (term-depends-p term k)) (poly-termlist p)))
1717
(defun ring-intersection (plist k)
1718
"This function assumes that polynomial list PLIST is a Grobner basis
1719
and it calculates the intersection with the ring R[x[k+1],...,x[n]], i.e.
1720
it discards polynomials which depend on variables x[0], x[1], ..., x[k]."
1721
(dotimes (i k plist)
1723
(remove-if #'(lambda (p)
1724
(poly-depends-p p i))
1727
(defun elimination-ideal (ring flist k
1728
&optional (top-reduction-only $poly_top_reduction_only) (start 0)
1729
&aux (*monomial-order*
1730
(or *elimination-order*
1731
(elimination-order k))))
1732
(ring-intersection (reduced-grobner ring flist start top-reduction-only) k))
1734
(defun colon-ideal (ring F G &optional (top-reduction-only $poly_top_reduction_only))
1735
"Returns the reduced Grobner basis of the colon ideal Id(F):Id(G),
1736
where F and G are two lists of polynomials. The colon ideal I:J is
1737
defined as the set of polynomials H such that for all polynomials W in
1738
J the polynomial W*H belongs to I."
1741
;;Id(G) consists of 0 only so W*0=0 belongs to Id(F)
1742
(if (every #'poly-zerop F)
1743
(error "First ideal must be non-zero.")
1746
(make-monom (monom-dimension (poly-lm (find-if-not #'poly-zerop F)))
1748
(funcall (ring-unit ring))))))))
1750
(colon-ideal-1 ring F (car G) top-reduction-only))
1752
(ideal-intersection ring
1753
(colon-ideal-1 ring F (car G) top-reduction-only)
1754
(colon-ideal ring F (rest G) top-reduction-only)
1755
top-reduction-only))))
1757
(defun colon-ideal-1 (ring F g &optional (top-reduction-only $poly_top_reduction_only))
1758
"Returns the reduced Grobner basis of the colon ideal Id(F):Id({G}), where
1759
F is a list of polynomials and G is a polynomial."
1760
(mapcar #'(lambda (x) (poly-exact-divide ring x g)) (ideal-intersection ring F (list g) top-reduction-only)))
1763
(defun ideal-intersection (ring F G &optional (top-reduction-only $poly_top_reduction_only)
1764
&aux (*monomial-order* (or *elimination-order*
1765
#'elimination-order-1)))
1766
(mapcar #'poly-contract
1770
(append (mapcar #'(lambda (p) (poly-extend p (list 1))) F)
1771
(mapcar #'(lambda (p)
1772
(poly-append (poly-extend (poly-uminus ring p) (list 1))
1779
(defun poly-lcm (ring f g)
1780
"Return LCM (least common multiple) of two polynomials F and G.
1781
The polynomials must be ordered according to monomial order PRED
1782
and their coefficients must be compatible with the RING structure
1783
defined in the COEFFICIENT-RING package."
1787
((and (endp (cdr (poly-termlist f))) (endp (cdr (poly-termlist g))))
1788
(let ((m (monom-lcm (poly-lm f) (poly-lm g))))
1789
(make-poly-from-termlist (list (make-term m (funcall (ring-lcm ring) (poly-lc f) (poly-lc g)))))))
1791
(multiple-value-bind (f f-cont)
1792
(poly-primitive-part ring f)
1793
(multiple-value-bind (g g-cont)
1794
(poly-primitive-part ring g)
1797
(funcall (ring-lcm ring) f-cont g-cont)
1798
(poly-primitive-part ring (car (ideal-intersection ring (list f) (list g) nil)))))))))
1800
;; Do two Grobner bases yield the same ideal?
1801
(defun grobner-equal (ring G1 G2)
1802
"Returns T if two lists of polynomials G1 and G2, assumed to be Grobner bases,
1803
generate the same ideal, and NIL otherwise."
1804
(and (grobner-subsetp ring G1 G2) (grobner-subsetp ring G2 G1)))
1806
(defun grobner-subsetp (ring G1 G2)
1807
"Returns T if a list of polynomials G1 generates
1808
an ideal contained in the ideal generated by a polynomial list G2,
1809
both G1 and G2 assumed to be Grobner bases. Returns NIL otherwise."
1810
(every #'(lambda (p) (grobner-member ring p G2)) G1))
1812
(defun grobner-member (ring p G)
1813
"Returns T if a polynomial P belongs to the ideal generated by the
1814
polynomial list G, which is assumed to be a Grobner basis. Returns NIL otherwise."
1815
(poly-zerop (normal-form ring p G nil)))
1817
;; Calculate F : p^inf
1818
(defun ideal-saturation-1 (ring F p start &optional (top-reduction-only $poly_top_reduction_only)
1819
&aux (*monomial-order* (or *elimination-order*
1820
#'elimination-order-1)))
1821
"Returns the reduced Grobner basis of the saturation of the ideal
1822
generated by a polynomial list F in the ideal generated by a single
1823
polynomial P. The saturation ideal is defined as the set of
1824
polynomials H such for some natural number n (* (EXPT P N) H) is in the ideal
1825
F. Geometrically, over an algebraically closed field, this is the set
1826
of polynomials in the ideal generated by F which do not identically
1827
vanish on the variety of P."
1833
(saturation-extension-1 ring F p)
1834
start top-reduction-only)
1839
;; Calculate F : p1^inf : p2^inf : ... : ps^inf
1840
(defun ideal-polysaturation-1 (ring F plist start &optional (top-reduction-only $poly_top_reduction_only))
1841
"Returns the reduced Grobner basis of the ideal obtained by a
1842
sequence of successive saturations in the polynomials
1843
of the polynomial list PLIST of the ideal generated by the
1846
((endp plist) (reduced-grobner ring F start top-reduction-only))
1847
(t (let ((G (ideal-saturation-1 ring F (car plist) start top-reduction-only)))
1848
(ideal-polysaturation-1 ring G (rest plist) (length G) top-reduction-only)))))
1850
(defun ideal-saturation (ring F G start &optional (top-reduction-only $poly_top_reduction_only)
1853
(*monomial-order* (or *elimination-order*
1854
(elimination-order k))))
1855
"Returns the reduced Grobner basis of the saturation of the ideal
1856
generated by a polynomial list F in the ideal generated a polynomial
1857
list G. The saturation ideal is defined as the set of polynomials H
1858
such for some natural number n and some P in the ideal generated by G
1859
the polynomial P**N * H is in the ideal spanned by F. Geometrically,
1860
over an algebraically closed field, this is the set of polynomials in
1861
the ideal generated by F which do not identically vanish on the
1864
#'(lambda (q) (poly-contract q k))
1866
(reduced-grobner ring
1867
(polysaturation-extension ring F G)
1872
(defun ideal-polysaturation (ring F ideal-list start &optional (top-reduction-only $poly_top_reduction_only))
1873
"Returns the reduced Grobner basis of the ideal obtained by a
1874
successive applications of IDEAL-SATURATION to F and lists of
1875
polynomials in the list IDEAL-LIST."
1877
((endp ideal-list) F)
1878
(t (let ((H (ideal-saturation ring F (car ideal-list) start top-reduction-only)))
1879
(ideal-polysaturation ring H (rest ideal-list) (length H) top-reduction-only)))))
1882
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1884
;; Set up the coefficients to be polynomials
1886
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1888
;; (defun poly-ring (ring vars)
1890
;; :parse #'(lambda (expr) (poly-eval ring expr vars))
1891
;; :unit #'(lambda () (poly-unit ring (length vars)))
1892
;; :zerop #'poly-zerop
1893
;; :add #'(lambda (x y) (poly-add ring x y))
1894
;; :sub #'(lambda (x y) (poly-sub ring x y))
1895
;; :uminus #'(lambda (x) (poly-uminus ring x))
1896
;; :mul #'(lambda (x y) (poly-mul ring x y))
1897
;; :div #'(lambda (x y) (poly-exact-divide ring x y))
1898
;; :lcm #'(lambda (x y) (poly-lcm ring x y))
1899
;; :ezgcd #'(lambda (x y &aux (gcd (poly-gcd ring x y)))
1901
;; (poly-exact-divide ring x gcd)
1902
;; (poly-exact-divide ring y gcd)))
1903
;; :gcd #'(lambda (x y) (poly-gcd x y))))
1906
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1908
;; Conversion from internal to infix form
1910
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1912
(defun coerce-to-infix (poly-type object vars)
1915
`(+ ,@(mapcar #'(lambda (term) (coerce-to-infix :term term vars)) object)))
1917
(coerce-to-infix :termlist (poly-termlist object) vars))
1919
`([ ,@(mapcar #'(lambda (p) (coerce-to-infix :polynomial p vars)) object)))
1921
`(* ,(term-coeff object)
1922
,@(mapcar #'(lambda (var power) `(expt ,var ,power))
1923
vars (monom-exponents (term-monom object)))))
1928
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1930
;; Maxima expression ring
1932
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1934
(defparameter *ExpressionRing*
1936
;;(defun coeff-zerop (expr) (meval1 `(($is) (($equal) ,expr 0))))
1937
:parse #'(lambda (expr)
1938
(when modulus (setf expr ($rat expr)))
1940
:unit #'(lambda () (if modulus ($rat 1) 1))
1941
:zerop #'(lambda (expr)
1942
;;When is exactly a maxima expression equal to 0?
1943
(cond ((numberp expr)
1948
(mrat (eql ($ratdisrep expr) 0))
1949
(otherwise (eql ($totaldisrep expr) 0))))))
1950
:add #'(lambda (x y) (m+ x y))
1951
:sub #'(lambda (x y) (m- x y))
1952
:uminus #'(lambda (x) (m- x))
1953
:mul #'(lambda (x y) (m* x y))
1954
;;(defun coeff-div (x y) (cadr ($divide x y)))
1955
:div #'(lambda (x y) (m// x y))
1956
:lcm #'(lambda (x y) (meval1 `((|$lcm|) ,x ,y)))
1957
:ezgcd #'(lambda (x y) (apply #'values (cdr ($ezgcd x y))))
1958
:gcd #'(lambda (x y) (second ($ezgcd x y)))))
1960
(defvar *MaximaRing* *ExpressionRing*
1961
"The ring of coefficients, over which all polynomials
1962
are assumed to be defined.")
1965
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1967
;; Maxima expression parsing
1969
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1971
(defun equal-test-p (expr1 expr2)
1972
(alike1 expr1 expr2))
1974
(defun coerce-maxima-list (expr)
1975
"Convert a Maxima list to Lisp list."
1977
((and (consp (car expr)) (eql (caar expr) 'mlist)) (cdr expr))
1980
(defun free-of-vars (expr vars) (apply #'$freeof `(,@vars ,expr)))
1982
(defun parse-poly (expr vars &aux (vars (coerce-maxima-list vars)))
1983
"Convert a maxima polynomial expression EXPR in variables VARS to internal form."
1984
(labels ((parse (arg) (parse-poly arg vars))
1985
(parse-list (args) (mapcar #'parse args)))
1987
((eql expr 0) (make-poly-zero))
1988
((member expr vars :test #'equal-test-p)
1989
(let ((pos (position expr vars :test #'equal-test-p)))
1990
(make-variable *MaximaRing* (length vars) pos)))
1991
((free-of-vars expr vars)
1992
;;This means that variable-free CRE and Poisson forms will be converted
1993
;;to coefficients intact
1994
(coerce-coeff *MaximaRing* expr vars))
1997
(mplus (reduce #'(lambda (x y) (poly-add *MaximaRing* x y)) (parse-list (cdr expr))))
1998
(mminus (poly-uminus *MaximaRing* (parse (cadr expr))))
2000
(if (endp (cddr expr)) ;unary
2002
(reduce #'(lambda (p q) (poly-mul *MaximaRing* p q)) (parse-list (cdr expr)))))
2005
((member (cadr expr) vars :test #'equal-test-p)
2006
;;Special handling of (expt var pow)
2007
(let ((pos (position (cadr expr) vars :test #'equal-test-p)))
2008
(make-variable *MaximaRing* (length vars) pos (caddr expr))))
2009
((not (and (integerp (caddr expr)) (plusp (caddr expr))))
2010
;; Negative power means division in coefficient ring
2011
;; Non-integer power means non-polynomial coefficient
2012
(mtell "~%Warning: Expression ~%~M~%contains power which is not a positive integer. Parsing as coefficient.~%"
2014
(coerce-coeff *MaximaRing* expr vars))
2015
(t (poly-expt *MaximaRing* (parse (cadr expr)) (caddr expr)))))
2016
(mrat (parse ($ratdisrep expr)))
2017
(mpois (parse ($outofpois expr)))
2019
(coerce-coeff *MaximaRing* expr vars)))))))
2021
(defun parse-poly-list (expr vars)
2023
(mlist (mapcar #'(lambda (p) (parse-poly p vars)) (cdr expr)))
2024
(t (merror "Expression ~M is not a list of polynomials in variables ~M."
2026
(defun parse-poly-list-list (poly-list-list vars)
2027
(mapcar #'(lambda (G) (parse-poly-list G vars)) (coerce-maxima-list poly-list-list)))
2030
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2034
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2035
(defun find-order (order)
2036
"This function returns the order function bases on its name."
2041
((lex :lex $lex) #'lex>)
2042
((grlex :grlex $grlex) #'grlex>)
2043
((grevlex :grevlex $grevlex) #'grevlex>)
2044
((invlex :invlex $invlex) #'invlex>)
2045
((elimination-order-1 :elimination-order-1 elimination_order_1) #'elimination-order-1)
2047
(mtell "~%Warning: Order ~M not found. Using default.~%" order))))
2049
(mtell "~%Order specification ~M is not recognized. Using default.~%" order)
2052
(defun find-ring (ring)
2053
"This function returns the ring structure bases on input symbol."
2058
((expression-ring :expression-ring $expression_ring) *ExpressionRing*)
2059
((ring-of-integers :ring-of-integers $ring_of_integers) *RingOfIntegers*)
2061
(mtell "~%Warning: Ring ~M not found. Using default.~%" ring))))
2063
(mtell "~%Ring specification ~M is not recognized. Using default.~%" ring)
2066
(defmacro with-monomial-order ((order) &body body)
2067
"Evaluate BODY with monomial order set to ORDER."
2068
`(let ((*monomial-order* (or (find-order ,order) *monomial-order*)))
2071
(defmacro with-coefficient-ring ((ring) &body body)
2072
"Evaluate BODY with coefficient ring set to RING."
2073
`(let ((*MaximaRing* (or (find-ring ,ring) *MaximaRing*)))
2076
(defmacro with-elimination-orders ((primary secondary elimination-order)
2078
"Evaluate BODY with primary and secondary elimination orders set to PRIMARY and SECONDARY."
2079
`(let ((*primary-elimination-order* (or (find-order ,primary) *primary-elimination-order*))
2080
(*secondary-elimination-order* (or (find-order ,secondary) *secondary-elimination-order*))
2081
(*elimination-order* (or (find-order ,elimination-order) *elimination-order*)))
2085
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2087
;; Conversion from internal form to Maxima general form
2089
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2091
(defun maxima-head ()
2092
(if $poly_return_term_list
2096
(defun coerce-to-maxima (poly-type object vars)
2099
`(,(maxima-head) ,@(mapcar #'(lambda (term) (coerce-to-maxima :term term vars)) (poly-termlist object))))
2101
`((mlist) ,@(mapcar #'(lambda (p) (coerce-to-maxima :polynomial p vars)) object)))
2103
`((mtimes) ,(term-coeff object)
2104
,@(mapcar #'(lambda (var power) `((mexpt) ,var ,power))
2105
vars (monom-exponents (term-monom object)))))
2110
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2112
;; Macro facility for writing Maxima-level wrappers for
2113
;; functions operating on internal representation
2115
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2117
(defmacro with-parsed-polynomials (((maxima-vars &optional (maxima-new-vars nil new-vars-supplied-p))
2118
&key (polynomials nil)
2120
(poly-list-lists nil)
2123
&aux (vars (gensym))
2124
(new-vars (gensym)))
2125
`(let ((,vars (coerce-maxima-list ,maxima-vars))
2126
,@(when new-vars-supplied-p
2127
(list `(,new-vars (coerce-maxima-list ,maxima-new-vars)))))
2130
(with-coefficient-ring ($poly_coefficient_ring)
2131
(with-monomial-order ($poly_monomial_order)
2132
(with-elimination-orders ($poly_primary_elimination_order
2133
$poly_secondary_elimination_order
2134
$poly_elimination_order)
2135
(let ,(let ((args nil))
2136
(dolist (p polynomials args)
2137
(setf args (cons `(,p (parse-poly ,p ,vars)) args)))
2138
(dolist (p poly-lists args)
2139
(setf args (cons `(,p (parse-poly-list ,p ,vars)) args)))
2140
(dolist (p poly-list-lists args)
2141
(setf args (cons `(,p (parse-poly-list-list ,p ,vars)) args))))
2143
,(if new-vars-supplied-p
2144
`(append ,vars ,new-vars)
2147
(defmacro define-unop (maxima-name fun-name
2148
&optional (documentation nil documentation-supplied-p))
2149
"Define a MAXIMA-level unary operator MAXIMA-NAME corresponding to unary function FUN-NAME."
2150
`(defun ,maxima-name (p vars
2152
(vars (coerce-maxima-list vars))
2153
(p (parse-poly p vars)))
2154
,@(when documentation-supplied-p (list documentation))
2155
(coerce-to-maxima :polynomial (,fun-name *MaximaRing* p) vars)))
2157
(defmacro define-binop (maxima-name fun-name
2158
&optional (documentation nil documentation-supplied-p))
2159
"Define a MAXIMA-level binary operator MAXIMA-NAME corresponding to binary function FUN-NAME."
2160
`(defmfun ,maxima-name (p q vars
2162
(vars (coerce-maxima-list vars))
2163
(p (parse-poly p vars))
2164
(q (parse-poly q vars)))
2165
,@(when documentation-supplied-p (list documentation))
2166
(coerce-to-maxima :polynomial (,fun-name *MaximaRing* p q) vars)))
2169
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2171
;; Maxima-level interface functions
2173
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2175
;; Auxillary function for removing zero polynomial
2176
(defun remzero (plist) (remove #'poly-zerop plist))
2180
(define-binop $poly_add poly-add
2181
"Adds two polynomials P and Q")
2183
(define-binop $poly_subtract poly-sub
2184
"Subtracts a polynomial Q from P.")
2186
(define-binop $poly_multiply poly-mul
2187
"Returns the product of polynomials P and Q.")
2189
(define-binop $poly_s_polynomial spoly
2190
"Returns the syzygy polynomial (S-polynomial) of two polynomials P and Q.")
2192
(define-unop $poly_primitive_part poly-primitive-part
2193
"Returns the polynomial P divided by GCD of its coefficients.")
2195
(define-unop $poly_normalize poly-normalize
2196
"Returns the polynomial P divided by the leading coefficient.")
2200
(defmfun $poly_expand (p vars)
2201
"This function is equivalent to EXPAND(P) if P parses correctly to a polynomial.
2202
If the representation is not compatible with a polynomial in variables VARS,
2203
the result is an error."
2204
(with-parsed-polynomials ((vars) :polynomials (p)
2205
:value-type :polynomial)
2208
(defmfun $poly_expt (p n vars)
2209
(with-parsed-polynomials ((vars) :polynomials (p) :value-type :polynomial)
2210
(poly-expt *MaximaRing* p n)))
2212
(defmfun $poly_content (p vars)
2213
(with-parsed-polynomials ((vars) :polynomials (p))
2214
(poly-content *MaximaRing* p)))
2216
(defmfun $poly_pseudo_divide (f fl vars
2217
&aux (vars (coerce-maxima-list vars))
2218
(f (parse-poly f vars))
2219
(fl (parse-poly-list fl vars)))
2220
(multiple-value-bind (quot rem c division-count)
2221
(poly-pseudo-divide *MaximaRing* f fl)
2223
,(coerce-to-maxima :poly-list quot vars)
2224
,(coerce-to-maxima :polynomial rem vars)
2228
(defmfun $poly_exact_divide (f g vars)
2229
(with-parsed-polynomials ((vars) :polynomials (f g) :value-type :polynomial)
2230
(poly-exact-divide *MaximaRing* f g)))
2232
(defmfun $poly_normal_form (f fl vars)
2233
(with-parsed-polynomials ((vars) :polynomials (f)
2235
:value-type :polynomial)
2236
(normal-form *MaximaRing* f (remzero fl) nil)))
2238
(defmfun $poly_buchberger_criterion (G vars)
2239
(with-parsed-polynomials ((vars) :poly-lists (G))
2240
(buchberger-criterion *MaximaRing* G)))
2242
(defmfun $poly_buchberger (fl vars)
2243
(with-parsed-polynomials ((vars) :poly-lists (fl) :value-type :poly-list)
2244
(buchberger *MaximaRing* (remzero fl) 0 nil)))
2246
(defmfun $poly_reduction (plist vars)
2247
(with-parsed-polynomials ((vars) :poly-lists (plist)
2248
:value-type :poly-list)
2249
(reduction *MaximaRing* plist)))
2251
(defmfun $poly_minimization (plist vars)
2252
(with-parsed-polynomials ((vars) :poly-lists (plist)
2253
:value-type :poly-list)
2254
(minimization plist)))
2256
(defmfun $poly_normalize_list (plist vars)
2257
(with-parsed-polynomials ((vars) :poly-lists (plist)
2258
:value-type :poly-list)
2259
(poly-normalize-list *MaximaRing* plist)))
2261
(defmfun $poly_grobner (F vars)
2262
(with-parsed-polynomials ((vars) :poly-lists (F)
2263
:value-type :poly-list)
2264
(grobner *MaximaRing* (remzero F))))
2266
(defmfun $poly_reduced_grobner (F vars)
2267
(with-parsed-polynomials ((vars) :poly-lists (F)
2268
:value-type :poly-list)
2269
(reduced-grobner *MaximaRing* (remzero F))))
2271
(defmfun $poly_depends_p (p var mvars
2272
&aux (vars (coerce-maxima-list mvars))
2273
(pos (position var vars)))
2275
(merror "~%Variable ~M not in the list of variables ~M." var mvars)
2276
(poly-depends-p (parse-poly p vars) pos)))
2278
(defmfun $poly_elimination_ideal (flist k vars)
2279
(with-parsed-polynomials ((vars) :poly-lists (flist)
2280
:value-type :poly-list)
2281
(elimination-ideal *MaximaRing* flist k nil 0)))
2283
(defmfun $poly_colon_ideal (F G vars)
2284
(with-parsed-polynomials ((vars) :poly-lists (F G) :value-type :poly-list)
2285
(colon-ideal *MaximaRing* F G nil)))
2287
(defmfun $poly_ideal_intersection (F G vars)
2288
(with-parsed-polynomials ((vars) :poly-lists (F G) :value-type :poly-list)
2289
(ideal-intersection *MaximaRing* F G nil)))
2291
(defmfun $poly_lcm (f g vars)
2292
(with-parsed-polynomials ((vars) :polynomials (f g) :value-type :polynomial)
2293
(poly-lcm *MaximaRing* f g)))
2295
(defmfun $poly_gcd (f g vars)
2296
($first ($divide (m* f g) ($poly_lcm f g vars))))
2298
(defmfun $poly_grobner_equal (G1 G2 vars)
2299
(with-parsed-polynomials ((vars) :poly-lists (G1 G2))
2300
(grobner-equal *MaximaRing* G1 G2)))
2302
(defmfun $poly_grobner_subsetp (G1 G2 vars)
2303
(with-parsed-polynomials ((vars) :poly-lists (G1 G2))
2304
(grobner-subsetp *MaximaRing* G1 G2)))
2306
(defmfun $poly_grobner_member (p G vars)
2307
(with-parsed-polynomials ((vars) :polynomials (p) :poly-lists (G))
2308
(grobner-member *MaximaRing* p G)))
2310
(defmfun $poly_ideal_saturation1 (F p vars)
2311
(with-parsed-polynomials ((vars) :poly-lists (F) :polynomials (p)
2312
:value-type :poly-list)
2313
(ideal-saturation-1 *MaximaRing* F p)))
2315
(defmfun $poly_saturation_extension (F plist vars new-vars)
2316
(with-parsed-polynomials ((vars new-vars)
2317
:poly-lists (F plist)
2318
:value-type :poly-list)
2319
(saturation-extension *MaximaRing* F plist)))
2321
(defmfun $poly_polysaturation_extension (F plist vars new-vars)
2322
(with-parsed-polynomials ((vars new-vars)
2323
:poly-lists (F plist)
2324
:value-type :poly-list)
2325
(polysaturation-extension *MaximaRing* F plist)))
2327
(defmfun $poly_ideal_polysaturation1 (F plist vars)
2328
(with-parsed-polynomials ((vars) :poly-lists (F plist)
2329
:value-type :poly-list)
2330
(ideal-polysaturation-1 *MaximaRing* F plist 0 nil)))
2332
(defmfun $poly_ideal_saturation (F G vars)
2333
(with-parsed-polynomials ((vars) :poly-lists (F G)
2334
:value-type :poly-list)
2335
(ideal-saturation *MaximaRing* F G 0 nil)))
2337
(defmfun $poly_ideal_polysaturation (F ideal-list vars)
2338
(with-parsed-polynomials ((vars) :poly-lists (F)
2339
:poly-list-lists (ideal-list)
2340
:value-type :poly-list)
2341
(ideal-polysaturation *MaximaRing* F ideal-list 0 nil)))
added
removed
share/contrib/Grobner/grobner.lisp
