source: CGBLisp/src/RCS/grobner.lisp,v@ 1

Last change on this file since 1 was 1, checked in by Marek Rychlik, 15 years ago

First import of a version circa 1997.

File size: 69.1 KB
Line 
1head 1.8;
2access;
3symbols;
4locks; strict;
5comment @;;; @;
6
7
81.8
9date 2009.01.23.10.40.26; author marek; state Exp;
10branches;
11next 1.7;
12
131.7
14date 2009.01.23.10.28.38; author marek; state Exp;
15branches;
16next 1.6;
17
181.6
19date 2009.01.23.10.28.05; author marek; state Exp;
20branches;
21next 1.5;
22
231.5
24date 2009.01.22.04.02.22; author marek; state Exp;
25branches;
26next 1.4;
27
281.4
29date 2009.01.19.09.25.46; author marek; state Exp;
30branches;
31next 1.3;
32
331.3
34date 2009.01.19.07.50.10; author marek; state Exp;
35branches;
36next 1.2;
37
381.2
39date 2009.01.19.06.49.18; author marek; state Exp;
40branches;
41next 1.1;
42
431.1
44date 2009.01.19.06.48.15; author marek; state Exp;
45branches;
46next ;
47
48
49desc
50@@
51
52
531.8
54log
55@*** empty log message ***
56@
57text
58@;;; -*- Mode: Lisp; Syntax: Common-Lisp; Package: Grobner; Base: 10 -*-
59#|
60 $Id$
61 *--------------------------------------------------------------------------*
62 | Copyright (C) 1994, Marek Rychlik (e-mail: rychlik@@math.arizona.edu) |
63 | Department of Mathematics, University of Arizona, Tucson, AZ 85721 |
64 | |
65 | Everyone is permitted to copy, distribute and modify the code in this |
66 | directory, as long as this copyright note is preserved verbatim. |
67 *--------------------------------------------------------------------------*
68|#
69(defpackage "GROBNER"
70 (:use "ORDER" "PARSE" "PRINTER" "DIVISION" "MONOM"
71 "MAKELIST" "COEFFICIENT-RING" "TERM" "POLY"
72 "POLY-WITH-SUGAR" "COMMON-LISP")
73 (:export
74 spoly
75 grobner
76 reduction
77 normal-form
78 reduced-grobner
79 normalize-poly
80 normalize-basis
81 elimination-ideal
82 ring-intersection
83 poly-contract
84 ideal-intersection
85 poly-lcm
86 grobner-gcd
87 pseudo-divide
88 buchberger-criterion
89 grobner-test
90 grobner-equal
91 grobner-subsetp
92 grobner-member
93 ideal-equal
94 ideal-subsetp
95 ideal-member
96 ideal-saturation
97 ideal-polysaturation
98 ideal-saturation-1
99 ideal-polysaturation-1
100 grobner-primitive-part
101 grobner-content
102 colon-ideal
103 *grobner-debug*
104 buchberger-set-pair-heuristic
105 gebauer-moeller-set-pair-heuristic
106 select-grobner-algorithm
107 spoly-with-sugar
108 normal-form-with-sugar
109 grobner-with-sugar
110 ))
111
112(in-package "GROBNER")
113
114#+debug(proclaim '(optimize (speed 0) (debug 3)))
115#-debug(proclaim '(optimize (speed 3) (debug 0)))
116
117#+debug
118(defvar *grobner-debug* nil
119 "If T, produce debugging and tracing output.")
120
121(defvar *buchberger-merge-pairs* 'buchberger-merge-pairs-smallest-lcm
122 "A variable holding the predicate used to sort critical pairs
123in the order of decreasing priority for the Buchberger algorithm.")
124
125(defvar *gebauer-moeller-merge-pairs* 'gebauer-moeller-merge-pairs-use-mock-spoly
126 "A variable holding the predicate used to sort critical pairs
127in the order of decreasing priority for the Gebauer-Moeller algorithm")
128
129(defvar *grobner-function* 'buchberger
130 "The name of the function used to calculate Grobner basis")
131
132(defun select-grobner-algorithm (algorithm)
133 "Select one of the several implementations of the Grobner basis
134algorithm."
135 (setf *grobner-function*
136 (ecase algorithm
137 (:buchberger 'buchberger)
138 (:cached-buchberger 'cached-buchberger)
139 (:gebauer-moeller 'gebauer-moeller)
140 (:sugar 'buchberger-with-sugar)
141 (:cached-sugar 'cached-buchberger-with-sugar)
142 (:gebauer-moeller-with-sugar 'gebauer-moeller-with-sugar))))
143
144(defun grobner (F &optional (pred #'lex>)
145 (start 0)
146 (top-reduction-only nil)
147 (ring *coefficient-ring*))
148 "Return Grobner basis of an ideal generated by polynomial list F.
149Assume that the monomials of each polynomial are ordered according to
150the admissible monomial order PRED. The result will be a list of
151polynomials ordered as well according to PRED. The polynomials from 0
152to START-1 in the list F are assumed to be a Grobner basis, and thus
153certain critical pairs do not have to be calculated. If
154TOP-REDUCTION-ONLY is not NIL then only top reductions will be
155performed, instead of full division, and thus the returned Grobner
156basis will have its polynomials only top-reduced with respect to the
157preceding polynomials. RING should be set to the RING structure (see
158COEFFICIENT-RING package) corresponding to the coefficients of the
159polynomials in the list F. This function is currently just an
160interface to several algorithms which fine Grobner bases. Use
161SELECT-GROBNER-ALGORITHM in order to change the algorithm."
162 (funcall *grobner-function* F pred start top-reduction-only ring))
163
164
165#+debug
166(defmacro debug-cgb (&rest args)
167 "This macro is used in printing debugging and tracing messages
168and its arguments are analogous to the FORMAT function, except
169that the stream argument is omitted. By default, the output is
170sent to *trace-output*, which can be set to produce output
171in a file."
172 `(when *grobner-debug* (format *trace-output* ,@@args)))
173
174(defun spoly (f g pred ring)
175 "Return the S-polynomial of F and G, which should
176be two polynomials with terms ordered by PRED and
177coefficients in a ring whose operations are the slots
178of the RING structure."
179 (declare (optimize (speed 3) (safety 0)))
180 (let* ((lcm (monom-lcm (lm f) (lm g)))
181 (m1 (monom/ lcm (lm f)))
182 (m2 (monom/ lcm (lm g))))
183 (let* ((c (funcall (ring-lcm ring) (lc f) (lc g)))
184 (cf (funcall (ring-/ ring) c (lc f)))
185 (cg (funcall (ring-/ ring) c (lc g))))
186 (poly-
187 (term-times-poly (cons m1 cf) (rest f) ring)
188 (term-times-poly (cons m2 cg) (rest g) ring)
189 pred
190 ring))))
191
192
193(defun grobner-primitive-part (p ring)
194 "Divide polynomial P with integer coefficients by gcd of its
195coefficients and return the result. Use the RING structure from the
196COEFFICIENT-RING package to perform the operations on the
197coefficients."
198 (if (poly-zerop p)
199 (values p 1)
200 (let ((c (grobner-content p ring)))
201 (values (mapcar
202 #'(lambda (x)
203 (cons (car x)
204 (funcall
205 (ring-/ ring)
206 (cdr x) c)))
207 p)
208 c))))
209
210
211(defun grobner-content (p ring)
212 "Greatest common divisor of the coefficients of the polynomial P. Use the RING structure
213to compute the greatest common divisor."
214 (let ((c (apply (ring-gcd ring) (mapcar #'cdr p))))
215 (when (/= (funcall (ring-signum ring) (cdar p))
216 (funcall (ring-signum ring) c))
217 (setf c (funcall (ring-- ring) c)))
218 c))
219
220(defun normal-form (f fl pred top-reduction-only ring)
221 "Calculate the normal form of the polynomial F with respect to
222the polynomial list FL. Destructive to F; thus, copy-tree should be used
223where F needs to be preserved. It returns multiple values.
224The first value is the polynomial R which is reduced (or top-reduced
225if TOP-REDUCTION-ONLY is not NIL). The second value is an integer
226which is the number of reductions actually performed. The third value
227is the coefficient by which F needs to be multiplied in order to be
228able to perform the division without actually having to divide in the
229coefficient ring (a form of pseudo-division is used). All operations
230on the coefficients are performed using the RING structure from the
231COEFFICIENT-RING package."
232 (declare (optimize (speed 3) (safety 0)))
233 ;; Loop invariant: c*f=sum ai*fi+r+p
234 (do (r (c (ring-unit ring))
235 (division-count 0)
236 (p f))
237 ((or (endp p)
238 (endp fl)
239 (and top-reduction-only r))
240 #+debug(progn
241 (debug-cgb "~&~3T~d reduction~:p" division-count)
242 (when (endp r)
243 (debug-cgb " ---> 0")))
244 (values (append (nreverse r) p)
245 division-count c))
246 (declare (fixnum division-count) (integer c))
247 (let ((g (find (lm p) fl :test #'monom-divisible-by-p :key #'lm)))
248 (cond
249 (g ;division possible
250 (incf division-count)
251 (let* ((gcd (funcall (ring-gcd ring) (lc g) (lc p)))
252 (c1 (funcall (ring-/ ring) (lc g) gcd))
253 (c2 (funcall (ring-/ ring) (lc p) gcd))
254 (m (monom/ (lm p) (lm g))))
255 ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
256 (setf r (scalar-times-poly c1 r ring)
257 c (funcall (ring-* ring) c c1)
258 p (grobner-op c2 c1 m (rest p) (rest g) pred ring))))
259 (t ;no division possible
260 (setf r (cons (lt p) r) ;move lt(p) to remainder
261 p (rest p))))))) ;remove lt(p) from p
262
263(defun buchberger (F pred start top-reduction-only ring
264 &aux (s (1- (length F))) B M)
265 "An implementation of the Buchberger algorithm. Return Grobner
266basis of the ideal generated by the polynomial list F.
267The terms of each polynomial are ordered by the admissible
268monomial order PRED. Polynomials 0 to START-1 are assumed to
269be a Grobner basis already, so that certain critical pairs
270will not be examined. If TOP-REDUCTION-ONLY set, top reduction
271will be preformed. The RING structure is used to perform operations
272on coefficents (see COEFFICIENT-RING package)."
273 (declare (fixnum s))
274 (unless (plusp s) (return-from buchberger F)) ;cut startup costs
275 #+debug (debug-cgb "~&GROBNER BASIS - BUCHBERGER ALGORITHM")
276 #+debug (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
277 #+grobner-check (when (plusp start)
278 (grobner-test (subseq F 0 start) (subseq F 0 start) pred ring))
279 (setf B (nconc (makelist (list i j) (i 0 (1- start)) (j start s))
280 (makelist (list i j) (i start (1- s)) (j (1+ i) s)))
281 M (make-hash-table :test #'equal)
282 B (buchberger-sort-pairs B F pred ring))
283 ;;Initialize treated pairs
284 (dotimes (i (1- start))
285 (do ((j (1+ i) (1+ j))) ((>= j start))
286 (setf (gethash (list i j) M) t)))
287 (do ((G F))
288 ((endp B)
289 #+grobner-check(grobner-test G F pred ring)
290 #+debug(debug-cgb "~&GROBNER END")
291 G)
292 (let ((pair (pop B)))
293 (unless (or (Criterion-1 pair G)
294 (Criterion-2 pair G M))
295 (let ((SP (normal-form (spoly (elt G (first pair))
296 (elt G (second pair))
297 pred
298 ring)
299 G pred top-reduction-only ring)))
300 (unless (poly-zerop SP) ;S-polynomial ---> 0;
301 (setf s (1+ s)
302 SP (grobner-primitive-part SP ring)
303 G (nconc G (list SP))
304 B (funcall *buchberger-merge-pairs* B (makelist (list i s)
305 (i 0 (1- s)))
306 G pred ring))
307 #+debug (debug-cgb "~&Polynomials: ~d; Pairs left: ~d; Treated pairs: ~d"
308 (length G) (length B)
309 (hash-table-count M)))))
310 (setf (gethash pair M) t))))
311
312(defun grobner-op (c1 c2 m A B pred ring &aux n A1 A2)
313 "Perform a calculation equivalent to
314(POLY- (SCALAR-TIMES-POLY C2 A)
315 (SCALAR-TIMES-POLY C1 (MONOM-TIMES-POLY M B))
316 PRED)
317but more efficiently; it destructively modifies A and stores the
318result in it; A is returned."
319 (dolist (x A) (setf (cdr x) (funcall (ring-* ring) c2 (cdr x))))
320 (unless B
321 (return-from grobner-op A))
322 (unless A
323 (return-from grobner-op
324 (mapcar
325 #'(lambda (x) (cons (monom* m (car x))
326 (funcall (ring-- ring)
327 (funcall (ring-* ring)
328 c1 (cdr x)))))
329 B)))
330 (setf A1 (cons nil A) A2 A1)
331 (tagbody
332 begin-outer-loop
333 (unless B
334 (return-from grobner-op (cdr A2)))
335 (setf n (monom* m (caar B)))
336 (tagbody
337 begin-inner-loop
338 (unless (cdr A1)
339 (setf (cdr A1)
340 (mapcar #'(lambda (x) (cons (monom* m (car x))
341 (funcall (ring-- ring)
342 (funcall (ring-* ring)
343 c1 (cdr x))))) B))
344 (return-from grobner-op (cdr A2)))
345 (multiple-value-bind (lm-greater lm-equal)
346 (funcall pred (caadr A1) n)
347 (cond
348 (lm-equal
349 (let ((s (funcall (ring-- ring)
350 (cdadr A1)
351 (funcall (ring-* ring) c1 (cdar B)))))
352 (if (funcall (ring-zerop ring) s) ;check for cancellation
353 (setf (cdr A1) (cddr A1))
354 (setf (cdadr A1) s))
355 (setf B (cdr B))
356 (go begin-outer-loop)))
357 (lm-greater
358 (setf A1 (cdr A1))
359 (go begin-inner-loop))
360 (t ;(lm B) is greater than (lm A)
361 (setf (cdr A1)
362 (cons
363 (cons
364 n
365 (funcall (ring-- ring)
366 (funcall (ring-* ring) c1 (cdar B))))
367 (cdr A1))
368 B (cdr B)
369 A1 (cdr A1))
370 (go begin-outer-loop)))))))
371
372;;----------------------------------------------------------------
373;; Pairs
374;;----------------------------------------------------------------
375
376(defun buchberger-sort-pairs (B G pred ring)
377 "Sort critical pairs B by the function which is
378the value of the vairable *buchberger-merge-pairs*.
379G is the list of polynomials which the elemnts
380of B point into. PRED is the admissible monomial order.
381The RING structure holds operations on the coefficients
382as described in the COEFFICIENT-RING package."
383 (funcall *buchberger-merge-pairs* nil B G pred ring))
384
385
386(defun mock-spoly (f g pred)
387 "Returns an upper-bound on the leading
388monomial of the S-polynomial of F and G, which are two polynomials
389sorted by the admissible monomial order PRED."
390 (let* ((lcm (monom-lcm (lm f) (lm g)))
391 (m1 (monom/ lcm (lm f)))
392 (m2 (monom/ lcm (lm g))))
393 (cond
394 ((endp (second f))
395 (monom* m2 (car (second g))))
396 ((endp (second g))
397 (monom* m1 (car (second f))))
398 (t (let ((n1 (monom* m2 (car (second g))))
399 (n2 (monom* m1 (car (second f)))))
400 (if (funcall pred n1 n2)
401 n1
402 n2))))))
403
404;;----------------------------------------------------------------
405;; A stragegy based on a predicted approximate S-polynomials
406;;----------------------------------------------------------------
407(defun buchberger-merge-pairs-use-mock-spoly (B C G pred ring)
408 "Merges lists of critical pairs B and C into one list.
409The pairs are assumed to be ordered according to the
410increasing value of MOCK-SPOLY, as determined by the admissible
411monomial order PRED. G is the list of polynomials which the pairs
412of integers in B and C point into, and RING is the ring structure
413defined in the COEFFICIENT-RING package."
414 (declare (ignore ring))
415 ;; Delete pairs consisting of two single monomials
416 (setf C (delete-if
417 #'(lambda (p) (and (endp (cdr (elt G (first p))))
418 (endp (cdr (elt G (second p))))))
419 C))
420 (labels ((pair-key (p) (mock-spoly (elt G (first p))
421 (elt G (second p)) pred))
422 (pair-order (x y) (funcall pred y x)))
423 (merge 'list B (sort C #'pair-order :key #'pair-key)
424 #'pair-order :key #'pair-key)))
425
426
427;;----------------------------------------------------------------
428;; A strategy based on the smallest lcm of the leading
429;; monomials of the two polynomials - so called normal strategy
430;;----------------------------------------------------------------
431(defun buchberger-merge-pairs-smallest-lcm (B C G pred ring)
432 "Merges lists of critical pairs B and C into a single ordered
433list. Implements a strategy of ordering pairs according to the
434smallest lcm of the leading monomials of the two polynomials - so
435called normal strategy."
436 (declare (ignore ring))
437 (labels ((pair-key (p) (monom-lcm (lm (elt G (first p)))
438 (lm (elt G (second p)))))
439 (pair-order (x y) (funcall pred y x)))
440 (merge 'list B (sort C #'pair-order :key #'pair-key)
441 #'pair-order :key #'pair-key)))
442
443;;----------------------------------------------------------------
444;; A strategy based on the smallest degree of the lcm of the leading
445;; monomials of the two polynomials
446;;----------------------------------------------------------------
447(defun buchberger-merge-pairs-use-smallest-degree (B C G pred ring)
448 "Mergest lists B and C of critical pairs. This strategy is based on
449 ordering the pairs according to the smallest degree of the lcm of the
450 leading monomials of the two polynomials."
451 (declare (ignore pred ring))
452 (labels ((pair-key (p) (total-degree
453 (monom-lcm (lm (elt G (first p)))
454 (lm (elt G (second p)))))))
455 (merge 'list B (sort C #'< :key #'pair-key) #'< :key #'pair-key)))
456
457(defun buchberger-merge-pairs-use-smallest-length (B C G pred ring)
458 "Mergest lists B and C of critical pairs. This strategy is based on
459 ordering the pairs according to the smallest total length of the
460 two polynomials, where length is the number of terms."
461 (declare (ignore pred ring))
462 (labels ((pair-key (p) (+
463 (length (elt G (first p)))
464 (length (elt G (second p))))))
465 (merge 'list B (sort C #'< :key #'pair-key) #'< :key #'pair-key)))
466
467(defun buchberger-merge-pairs-use-smallest-coefficient-length (B C G pred ring)
468 "Mergest lists B and C of critical pairs. This strategy is based on
469 ordering the pairs according to the smallest combined length of the
470 coefficients of the two polynomials."
471 (declare (ignore pred))
472 (labels ((coefficient-length (poly)
473 (apply #'+ (mapcar (ring-length ring)
474 (mapcar #'cdr poly))))
475 (pair-key (p) (+
476 (coefficient-length (elt G (first p)))
477 (coefficient-length (elt G (second p))))))
478 (merge 'list B (sort C #'< :key #'pair-key) #'< :key #'pair-key)))
479
480(defun buchberger-set-pair-heuristic
481 (method
482 &aux
483 (strategy-fn
484 (ecase method
485 (:minimal-mock-spoly #'buchberger-merge-pairs-use-mock-spoly)
486 (:minimal-lcm #'buchberger-merge-pairs-smallest-lcm)
487 (:minimal-total-degree #'buchberger-merge-pairs-use-smallest-degree)
488 (:minimal-length #'buchberger-merge-pairs-use-smallest-length)
489 (:minimal-coefficient-length #'buchberger-merge-pairs-use-smallest-coefficient-length))))
490 "Simply sets the variable *buchberger-merge-pairs* to the heuristic
491function STRATEGY-FN implementing one of several strategies introduces
492before of selecting the most promising. METHOD should be one of the
493listed keywords."
494 (setf *buchberger-merge-pairs* strategy-fn))
495
496;;----------------------------------------------------------------
497;; Grobner Criteria
498;;----------------------------------------------------------------
499(defun Criterion-1 (pair G &aux (i (first pair)) (j (second pair)))
500 "Returns T if the leading monomials of the two polynomials
501in G pointed to by the integers in PAIR have disjoint (relatively prime)
502monomials. This test is known as the first Buchberger criterion."
503 (monom-rel-prime (lm (elt G i))
504 (lm (elt G j))))
505
506
507(defun Criterion-2 (pair G M &aux (i (first pair)) (j (second pair)))
508 "Returns T if the leading monomial of some element P of G divides
509the LCM of the leading monomials of the two polynomials in the
510polynomial list G, pointed to by the two integers in PAIR, and P
511paired with each of the polynomials pointed to by the the PAIR has
512already been treated, as indicated by the hash table M, which stores
513value T for every treated pair so far."
514 (labels ((make-pair (u v) (if (< u v) (list u v) (list v u))))
515 (dotimes (k (length G) nil)
516 (when (and (/= k i)
517 (/= k j)
518 (monom-divides-p
519 (lm (elt G k))
520 (monom-lcm (lm (elt G i)) (lm (elt G j))))
521 (gethash (make-pair i k) M)
522 (gethash (make-pair k j) M))
523 #+debug(when *grobner-debug* (format t ":B2"))
524 (return-from Criterion-2 t)))))
525
526;;----------------------------------------------------------------
527;; End of GROBNER implementation
528;;----------------------------------------------------------------
529
530(defun normalize-poly (p ring)
531 "Divide a polynomial by its leading coefficient. It assumes
532that the division is possible, which may not always be the
533case in rings which are not fields. The exact division operator
534is assumed to be provided by the RING structure of the
535COEFFICIENT-RING package."
536 (mapcar #'(lambda (term)
537 (cons
538 (car term)
539 (funcall (ring-/ ring) (cdr term) (lc p)))) p))
540
541(defun normalize-basis (plist ring)
542 "Divide every polynomial in a list PLIST by its leading coefficient.
543Use RING structure to perform the division of the coefficients."
544 (mapcar #'(lambda (x) (normalize-poly x ring)) plist))
545
546(defun reduction (P pred ring)
547 "Reduce a list of polynomials P, so that non of the terms in any of the
548polynomials is divisible by a leading monomial of another polynomial.
549Return the reduced list."
550 (do ((Q P)
551 (found t))
552 ((not found)
553 (mapcar #'(lambda (x) (grobner-primitive-part x ring)) Q))
554 ;;Find p in Q such that p is reducible mod Q\{p}
555 (setf found nil)
556 (dolist (x Q)
557 (let ((Q1 (remove x Q)))
558 (multiple-value-bind (h div-count)
559 (normal-form x Q1
560 pred
561 nil ;not a top reduction
562 ring)
563 (unless (zerop div-count)
564 (setf found t Q Q1)
565 (when h
566 (setf Q (nconc Q1 (list h))))
567 (return)))))))
568
569
570(defun reduced-grobner (F &optional (pred #'lex>)
571 (start 0)
572 (top-reduction-only nil)
573 (ring *coefficient-ring*))
574 "Return the reduced Grobner basis of the ideal generated by a polynomial
575list F. This combines calls to two functions: GROBNER and REDUCTION.
576The parameters have the same meaning as in GROBNER or BUCHBERGER."
577 (reduction (grobner F pred start top-reduction-only ring) pred ring))
578
579#+kcl
580(defvar *vars* nil "Keep track of variables for tracing SPOLY.")
581
582;; Does the monomial P depend on variable K
583(defun monom-depends-p (m k)
584 "Return T if the monomial M depends on variable number K."
585 (plusp (elt m k)))
586
587;; Does the term depend on variable K
588(defun term-depends-p (term k)
589 "Return T if the term TERM depends on variable number K."
590 (monom-depends-p (car term) k))
591
592;; Does the polynomial P depend on variable K
593(defun poly-depends-p (p k)
594 "Return T if the term polynomial P depends on variable number K."
595 (some #'(lambda (term) (term-depends-p term k)) p))
596
597(defun ring-intersection (plist k &key (key #'identity))
598 "This function assumes that polynomial list PLIST is a Grobner basis
599and it calculates the intersection with the ring R[x[k+1],...,x[n]], i.e.
600it discards polynomials which depend on variables x[0], x[1], ..., x[k]."
601 (dotimes (i k plist)
602 (setf plist
603 (remove-if #'(lambda (p)
604 (poly-depends-p p (funcall key i)))
605 plist))))
606
607(defun elimination-ideal (flist k
608 &key
609 (primary-order #'lex>)
610 (secondary-order #'lex>)
611 (start 0)
612 (order (elimination-order
613 k
614 :primary-order primary-order
615 :secondary-order secondary-order
616 ))
617 (top-reduction-only nil)
618 (ring *coefficient-ring*))
619 "Returns the K-th elimination ideal of the ideal generated by polynomial list FLIST.
620Thus, a Grobner basis of the ideal generated by FLIST is calculated and
621polynomials depending on variables 0-K are discarded. The monomial order in the
622calculation is an elimination order formed by combining two orders: PRIMARY-ORDER
623used on variables 0-K and SECONDARY-ORDER used on variables starting from variable K+1.
624Thus, if both orders are set to the lexicographic order #'LEX> then the resulting
625order is simply equivalent to #'LEX>. But one could use arbitrary two orders in this
626function, for instance, two copies of #'GREVLEX> (see ORDER package). When doing so,
627the caller must ensure that the same order has been used to sort the terms of
628the polynomials FLIST."
629 (ring-intersection
630 (reduced-grobner flist order start top-reduction-only ring)
631 k))
632
633
634(defun ideal-intersection (F G pred top-reduction-only ring)
635 "Return the Grobner basis of the intersection of the ideals
636generated by the polynomial lists F and G. PRED is an admissible
637monomial order with respect to which the terms of F and G have been
638sorted. This order is going to be used in the Grobner basis
639calculation. If TOP-REDUCTION-ONLY is not NIL, internally the Grobner
640basis algorithm will perform top reduction only. The RING parameter,
641as usual, should be set to a structure defined in the package
642COEFFICIENT-RING, and it will be used to perform all operations on
643coefficients."
644 (mapcar
645 #'poly-contract
646 (ring-intersection
647 (reduced-grobner
648 (append
649 (mapcar #'(lambda (p) (poly-extend p (list 1))) F)
650 (mapcar #'(lambda (p)
651 (append (poly-extend (minus-poly p ring) (list 1))
652 (poly-extend p)))
653 G))
654 (elimination-order-1 pred)
655 0
656 top-reduction-only
657 ring)
658 1)))
659
660
661(defun poly-contract (f &optional (k 1))
662 "Return a polynomial obtained from polynomial F by dropping the
663first K variables, if F does not depend on them. Note that this
664operation makes the polynomial incompatible for certain arithmetical
665operations with the original polynomial F. Calling this function on a polynomial F
666which does depend on variables 0-K will result in error."
667 (when (some #'(lambda (term) (find-if-not #'zerop (car term) :end k)) f)
668 (error "Contracting a polynomial which depends on contracted variables"))
669 (mapcar #'(lambda (term) (cons (nthcdr k (car term)) (cdr term))) f))
670
671
672(defun poly-lcm (f g &optional (pred #'lex>) (ring *coefficient-ring*))
673 "Return LCM (least common multiple) of two polynomials F and G.
674The polynomials must be ordered according to monomial order PRED
675and their coefficients must be compatible with the RING structure
676defined in the COEFFICIENT-RING package."
677 (cond
678 ((endp f) nil)
679 ((endp g) nil)
680 ((and (endp (cdr f)) (endp (cdr g)))
681 (list (cons (monom-lcm (caar f) (caar g))
682 (lcm (cdar f) (cdar g)))))
683 (t
684 (multiple-value-bind (f f-cont)
685 (grobner-primitive-part f ring)
686 (multiple-value-bind (g g-cont)
687 (grobner-primitive-part g ring)
688 (scalar-times-poly
689 (funcall (ring-lcm ring) f-cont g-cont)
690 (grobner-primitive-part (car (ideal-intersection (list f) (list g) pred nil ring))
691 ring)
692 ring))))))
693
694;; GCD of two polynomials, using Grobner bases
695(defun grobner-gcd (f g &optional (pred #'lex>) (ring *coefficient-ring*))
696 "Return GCD (greatest common divisor) of two polynomials F and G.
697The polynomials must be ordered according to monomial order PRED
698and their coefficients must be compatible with the RING structure
699defined in the COEFFICIENT-RING package."
700 (poly-exact-divide (poly* f g pred ring)
701 (poly-lcm f g pred ring)
702 pred ring))
703
704;; Do two Grobner bases yield the same ideal?
705(defun grobner-equal (g1 g2 &optional (pred #'lex>) (ring *coefficient-ring*))
706 "Returns T if two lists of polynomials G1 and G2, assumed to be Grobner bases,
707generate the same ideal, and NIL otherwise."
708 (and (grobner-subsetp g1 g2 pred ring)
709 (grobner-subsetp g2 g1 pred ring)))
710
711(defun grobner-subsetp (g1 g2 &optional (pred #'lex>) (ring *coefficient-ring*))
712 "Returns T if a list of polynomials G1 generates
713an ideal contained in the ideal generated by a polynomial list G2,
714both G1 and G2 assumed to be Grobner bases. Returns NIL otherwise."
715 (every #'(lambda (p) (grobner-member p g2 pred ring)) g1))
716
717(defun grobner-member (p g &optional (pred #'lex>) (ring *coefficient-ring*))
718 "Returns T if a polynomial P belongs to the ideal generated by the
719polynomial list G, which is assumed to be a Grobner basis. Returns NIL otherwise."
720 (endp (normal-form (copy-tree p) g pred nil ring)))
721
722(defun ideal-equal (f1 f2 pred top-reduction-only ring)
723 "Returns T if two ideals generated by polynomial lists F1 and F2 are identical.
724Returns NIL otherwise."
725 (grobner-equal (grobner f1 pred 0 top-reduction-only ring)
726 (grobner f2 pred 0 top-reduction-only ring)
727 pred ring))
728
729(defun ideal-subsetp (f1 f2 &optional (pred #'lex>) (ring *coefficient-ring*))
730 "Returns T if the ideal spanned by the polynomial list F1 is contained
731in the ideal spanned by the polynomial list F2. Returns NIL otherwise."
732 (grobner-subsetp
733 (grobner f1 pred 0 nil ring)
734 (grobner f2 pred 0 nil ring)
735 pred ring))
736
737(defun ideal-member (p plist pred top-reduction-only ring)
738 "Returns T if the polynomial P belongs to the ideal spanned by the polynomial
739list PLIST, and NIL otherwise."
740 (endp (normal-form (copy-tree p) (grobner plist pred 0 top-reduction-only ring)
741 pred nil ring)))
742
743;; Calculate F : p^inf
744(defun ideal-saturation-1 (F p pred start top-reduction-only ring
745 &aux (pred (elimination-order-1 pred)))
746 "Returns the reduced Grobner basis of the saturation of the ideal
747generated by a polynomial list F in the ideal generated by a single
748polynomial P. The saturation ideal is defined as the set of
749polynomials H such for some natural number n (* (EXPT P N) H) is in the ideal
750F. Geometrically, over an algebraically closed field, this is the set
751of polynomials in the ideal generated by F which do not identically
752vanish on the variety of P."
753 (mapcar
754 #'poly-contract
755 (ring-intersection
756 (reduced-grobner
757 (saturation-extension-1 F p ring)
758 pred
759 start top-reduction-only ring)
760 1)))
761
762
763(defun add-variables (plist n)
764 "Add N new variables, adn the first N variables, to every polynomial in polynomial list PLIST.
765The resulting polynomial will not depend on these N variables."
766 (mapcar
767 #'(lambda (g)
768 (mapcar #'(lambda (q) (cons (append (make-list n :initial-element 0)
769 (car q))
770 (cdr q))) g))
771 plist))
772
773;; Calculate <u1*p1,u2*p2,...,uk*pk>
774(defun extend-polynomials (plist &aux (k (length plist)))
775 "Returns polynomial list {U1*P1', U2*P2', ... , UK*PK'}
776where Ui are new variables and PLIST={P1, P2, ... , PK} is a polynomial
777list. PI' is obtained from PI by adding new variables U1, U2, ..., UN
778UK as the first K variables. Thus, the resulting list is consists of polynomials whose
779terms are ordered by the lexicographic order on variables UI, with
780ties resolved by the monomial order which was used to order the terms
781of the polynomials in PLIST. We note that the monomial order does not
782explicitly participate in this calculation, and neither do the
783variable names."
784 (labels ((incf-power (g i)
785 (dolist (x g g)
786 (incf (nth i (car x))))))
787 (setf plist (add-variables plist k))
788 (dotimes (i k plist)
789 (setf (nth i plist) (incf-power (nth i plist) i)))))
790
791;; Calculate <F, u1*p1-1,u2*p2-1,...,uk*pk-1>
792(defun saturation-extension (F plist ring &aux (k (length plist))
793 (d (+ k (length (caaar plist)))))
794 "Returns F' union {U1*P1-1,U2*P2-1,...,UK*PK-1} where Ui are new variables
795and F' is a polynomial list obtained from F by adding variables U1, U2, ..., Uk
796as the first K variables to each polynomial in F. Thus, the resulting list is
797consists of polynomials whose terms are ordered by the lexicographic order on
798variables Ui, with ties resolved by the monomial order which was used to order the
799terms of the polynomials in F. We note that the monomial order does not explicitly
800participate in this calculation, and neither do the variable names."
801 (setf F (add-variables F k)
802 plist (mapcar #'(lambda (x)
803 (nconc x (list (cons (make-list d :initial-element 0)
804 (funcall (ring-- ring)
805 (ring-unit ring))))))
806 (extend-polynomials plist)))
807 (append F plist))
808
809;; Calculate <F, u1*p1+u2*p2+...+uk*pk-1>
810(defun polysaturation-extension (F plist ring &aux (k (length plist))
811 (d (+ k (length (caaar plist)))))
812 "Returns F' union {U1*P1+U2*P2+UK*PK-1} where Ui are new variables
813and F' is a polynomial list obtained from F by adding variables U1, U2, ..., Uk
814as the first K variables to each polynomial in F. Thus, the resulting list is
815consists of polynomials whose terms are ordered by the lexicographic order on
816variables Ui, with ties resolved by the monomial order which was used to order the
817terms of the polynomials in F. We note that the monomial order does not explicitly
818participate in this calculation, and neither do the variable names."
819 (setf F (add-variables F k)
820 plist (apply #'append (extend-polynomials plist))
821 (cdr (last plist)) (list (cons (make-list d :initial-element 0)
822 (funcall (ring-- ring)
823 (ring-unit ring)))))
824 (append F (list plist)))
825
826(defun saturation-extension-1 (F p ring)
827 "Return the list F' union {U*P-1} where U is a new variable,
828where F' is obtained from the polynomial list F by adding U as the first variable.
829The variable U will be added as the first variable, so that it can be easily
830eliminated."
831 (saturation-extension F (list p) ring))
832
833;; Calculate F : p1^inf : p2^inf : ... : ps^inf
834(defun ideal-polysaturation-1 (F plist pred start top-reduction-only ring)
835 "Returns the reduced Grobner basis of the ideal obtained by a
836sequence of successive saturations in the polynomials
837of the polynomial list PLIST of the ideal generated by the
838polynomial list F."
839 (cond
840 ((endp plist) (reduced-grobner F pred start top-reduction-only ring))
841 (t (let ((G (ideal-saturation-1 F (car plist) pred start top-reduction-only ring)))
842 (ideal-polysaturation-1 G (rest plist) pred (length G) top-reduction-only ring)))))
843
844(defun ideal-saturation (F G pred start top-reduction-only ring
845 &aux (k (length G))
846 (pred (elimination-order k
847 :primary-order #'lex>
848 :secondary-order pred)))
849 "Returns the reduced Grobner basis of the saturation of the ideal
850generated by a polynomial list F in the ideal generated a polynomial
851list G The saturation ideal is defined as the set of polynomials H
852such for some natural number n and some P in the ideal generated by G
853the polynomial P**N * H is in the ideal spanned by F. Geometrically,
854over an algebraically closed field, this is the set of polynomials in
855the ideal generated by F which do not identically vanish on the
856variety of G."
857 (mapcar
858 #'(lambda (q) (poly-contract q k))
859 (ring-intersection
860 (reduced-grobner (polysaturation-extension F G ring)
861 pred start
862 top-reduction-only ring)
863 k)))
864
865(defun ideal-polysaturation (F ideal-list pred start top-reduction-only ring)
866 "Returns the reduced Grobner basis of the ideal obtained by a
867successive applications of IDEAL-SATURATIONS to F and
868lists of polynomials in the list IDEAL-LIST."
869 (cond
870 ((endp ideal-list) F)
871 (t (let ((H (ideal-saturation F (car ideal-list) pred start top-reduction-only ring)))
872 (ideal-polysaturation H (rest ideal-list) pred (length H) top-reduction-only ring)))))
873
874(defun buchberger-criterion (G pred ring)
875 "Returns T if G is a Grobner basis, by using the Buchberger
876criterion: for every two polynomials h1 and h2 in G the S-polynomial
877S(h1,h2) reduces to 0 modulo G."
878 (every
879 #'endp
880 (makelist (normal-form (spoly (elt G i) (elt G j) pred ring) G pred nil ring)
881 (i 0 (- (length G) 2))
882 (j (1+ i) (1- (length G))))))
883
884
885(defun grobner-test (G F pred ring)
886 "Test whether G is a Grobner basis and F is contained in G. Return T
887upon success and NIL otherwise."
888 #+debug(debug-cgb "~&GROBNER CHECK: ")
889 (let ((*grobner-debug* nil)
890 (stat1 (buchberger-criterion G pred ring))
891 (stat2
892 (every #'endp
893 (makelist (normal-form (copy-tree (elt F i)) G pred nil ring)
894 (i 0 (1- (length F)))))))
895 (unless stat1 (error "~&Buchberger criterion failed."))
896 (unless stat2
897 (error "~&Original polys not in ideal spanned by Grobner.")))
898 #+debug(debug-cgb "~&GROBNER CHECK END")
899 T)
900
901
902(defun minimization (P pred)
903 "Returns a sublist of the polynomial list P spanning the same
904monomial ideal as P but minimal, i.e. no leading monomial
905of a polynomial in the sublist divides the leading monomial
906of another polynomial."
907 (declare (ignore pred))
908 (do ((Q P)
909 (found t))
910 ((not found) Q)
911 ;;Find p in Q such that lm(p) is in LM(Q\{p})
912 (setf found nil
913 Q (dolist (x Q Q)
914 (let ((Q1 (remove x Q)))
915 (when (member-if #'(lambda (p) (monom-divides-p (lm x) (lm p))) Q1)
916 (setf found t)
917 (return Q1)))))))
918
919
920(defun add-minimized (f Gred pred)
921 "Adds a polynomial f to GRED, reduced Grobner basis, preserving the
922property described in the documentation for MINIMIZATION."
923 (if (member-if #'(lambda (p) (monom-divides-p (lm p) (lm f))) Gred)
924 Gred
925 (merge 'list (list f)
926 (delete-if #'(lambda (p) (monom-divides-p (lm f) (lm p))) Gred)
927 pred :key #'lm)))
928
929(defun colon-ideal (F G pred top-reduction-only ring)
930 "Returns the reduced Grobner basis of the colon ideal Id(F):Id(G), where
931F and G are two lists of polynomials. The colon ideal I:J is defined as the
932set of polynomials H such that there is a polynomial W in J for which
933W*H belongs to I."
934 (cond
935 ((endp G)
936 (if (every #'endp F) nil
937 (list (list
938 (cons
939 (make-list (length (caar (find-if-not #'endp F)))
940 :initial-element 0)
941 1)))))
942 ((endp (cdr G))
943 (colon-ideal-1 F (car G) pred top-reduction-only ring))
944 (t
945 (ideal-intersection (colon-ideal-1 F (car G) pred top-reduction-only ring)
946 (colon-ideal F (rest G) pred top-reduction-only ring)
947 pred
948 top-reduction-only
949 ring))))
950
951(defun colon-ideal-1 (F g pred top-reduction-only ring)
952 "Returns the reduced Grobner basis of the colon ideal Id(F):Id({G}), where
953F is a list of polynomials and G is a polynomial."
954 (mapcar #'(lambda (x) (poly-exact-divide x g pred ring))
955 (ideal-intersection F (list g) pred top-reduction-only ring)))
956
957
958(defun pseudo-divide (f fl pred ring)
959 "Pseudo-divide a polynomial F by the list of polynomials FL. Return
960multiple values. The first value is a list of quotients A.
961The second value is the remainder R. The third value is an integer
962count of the number of reductions performed. Finally, the fourth
963argument is a scalar coefficient C, such that C*F can be divided by FL
964within the ring of coefficients, which is not necessarily a field.
965The resulting objects satisfy the equation:
966 C*F= sum A[i]*FL[i] + R
967"
968 (declare (optimize (speed 3) (safety 0)))
969 (do (r
970 (c (ring-unit ring))
971 (a (make-list (length fl)))
972 (division-count 0)
973 (p f))
974 ((endp p)
975 #+debug(debug-cgb "~&~3T~d reduction~:p" division-count)
976 #+debug(when (endp r) (debug-cgb " ---> 0"))
977 (values a (nreverse r) division-count c))
978 (declare (fixnum division-count) (integer c))
979 (do ((fl fl (rest fl)) ;scan list of divisors
980 (b a (rest b)))
981 ((cond
982 ((endp fl) ;no division occurred
983 (setf r (cons (lt p) r) ;move lt(p) to remainder
984 p (rest p)) ;remove lt(p) from p
985 t)
986 ((monom-divides-p (lm (car fl)) (lm p)) ;division occurred
987 (incf division-count)
988 (let* ((gcd (funcall (ring-gcd ring) (lc (car fl)) (lc p)))
989 (c1 (funcall (ring-/ ring) (lc (car fl)) gcd))
990 (c2 (funcall (ring-/ ring) (lc p) gcd))
991 (quot (cons (monom/ (lm p) (lm (car fl))) c2)))
992 ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
993 (mapl #'(lambda (x)
994 (setf (car x) (scalar-times-poly c1 (car x) ring)))
995 a)
996 (setf p (scalar-times-poly c1 (rest p) ring) ;drop first term
997 r (scalar-times-poly c1 r ring)
998 c (funcall (ring-* ring) c c1)
999 p (poly-op p quot (rest (car fl)) pred ring)
1000 (car b) (nconc (car b) (list quot))))
1001 t))))))
1002
1003;;----------------------------------------------------------------
1004;; An implementation of the algorithm of Gebauer and Moeller,
1005;; as described in the book of Becker-Weispfenning, p. 232
1006;;----------------------------------------------------------------
1007(defun gebauer-moeller (F pred start top-reduction-only ring
1008 &aux B G F1)
1009 "Compute Grobner basis by using the algorithm of Gebauer and Moeller.
1010This algorithm is described as BUCHBERGERNEW2 in the book by
1011Becker-Weispfenning entitled ``Grobner Bases''"
1012 ;;cut startup costs if (length F)<=1
1013 (cond
1014 ((endp F) (return-from gebauer-moeller nil))
1015 ((endp (cdr F))
1016 (return-from gebauer-moeller (list (grobner-primitive-part (car F) ring)))))
1017 #+debug (debug-cgb "~&GROBNER BASIS - GEBAUER MOELLER ALGORITHM")
1018 #+debug (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
1019 #+grobner-check (when (plusp start)
1020 (grobner-test (subseq F 0 start) (subseq F 0 start)
1021 pred ring))
1022 (setf B nil
1023 G (subseq F 0 start)
1024 F1 (subseq F start))
1025 (do () ((endp F1))
1026 (multiple-value-setq (G B)
1027 (update G B (grobner-primitive-part (pop F1) ring) pred ring)))
1028 (do () ((endp B))
1029 (let* ((pair (pop B))
1030 (g1 (first pair))
1031 (g2 (second pair))
1032 (h (normal-form (spoly g1 g2 pred ring)
1033 G pred top-reduction-only ring)))
1034 (unless (poly-zerop h)
1035 (setf h (grobner-primitive-part h ring))
1036 (multiple-value-setq (G B)
1037 (update G B h pred ring))
1038 #+debug (debug-cgb "~&Polynomials: ~d; Pairs left: ~d~%"
1039 (length G) (length B))
1040 )))
1041 #+grobner-check(grobner-test G F pred ring)
1042 #+debug(debug-cgb "~&GROBNER END")
1043 G)
1044
1045(defun update (G B h pred ring &aux C D E B-new G-new pair)
1046 "An implementation of the auxillary UPDATE algorithm used by the
1047Gebauer-Moeller algorithm. G is a list of polynomials, B is a list of
1048critical pairs and H is a new polynomial which possibly will be added
1049to G. The naming conventions used are very close to the one used in
1050the book of Becker-Weispfenning."
1051 (setf C G D nil)
1052 (do (g1) ((endp C))
1053 (setf g1 (pop C))
1054 (when (or (monom-rel-prime (lm h) (lm g1))
1055 (and
1056 (not (some
1057 #'(lambda (g2) (monom-divides-p (monom-lcm (lm h) (lm g2))
1058 (monom-lcm (lm h) (lm g1))))
1059 C))
1060 (not (some
1061 #'(lambda (g2) (monom-divides-p (monom-lcm (lm h) (lm g2))
1062 (monom-lcm (lm h) (lm g1))))
1063 D))))
1064 (push g1 D)))
1065 (setf E nil)
1066 (do (g1) ((endp D))
1067 (setf g1 (pop D))
1068 (unless (monom-rel-prime (lm h) (lm g1))
1069 (push g1 E)))
1070 (setf B-new nil)
1071 (do (g1 g2) ((endp B))
1072 (setf pair (pop B)
1073 g1 (first pair)
1074 g2 (second pair))
1075 (when (or (not (monom-divides-p (lm h) (monom-lcm (lm g1) (lm g2))))
1076 (equal (monom-lcm (lm g1) (lm h))
1077 (monom-lcm (lm g1) (lm g2)))
1078 (equal (monom-lcm (lm h) (lm g2))
1079 (monom-lcm (lm g1) (lm g2))))
1080 (setf B-new (funcall *gebauer-moeller-merge-pairs* B-new (list (list g1 g2)) pred ring))))
1081 (dolist (g3 E)
1082 (setf B-new (funcall *gebauer-moeller-merge-pairs* B-new (list (list h g3)) pred ring)))
1083 (setf G-new nil)
1084 (do (g1) ((endp G))
1085 (setf g1 (pop G))
1086 (unless (monom-divides-p (lm h) (lm g1))
1087 (push g1 G-new)))
1088 (push h G-new)
1089 (values G-new B-new))
1090
1091(defun gebauer-moeller-merge-pairs-use-mock-spoly (B C pred ring)
1092 "The merging strategy used by the Gebauer-Moeller algorithm, based
1093on MOCK-SPOLY; see the documentation of
1094BUCHBERGER-MERGE-PAIRS-USE-MOCK-SPOLY."
1095 (declare (ignore ring))
1096 ;; Delete pairs consisting of two single monomials
1097 (setf C (delete-if
1098 #'(lambda (p) (and (endp (cdr (first p)))
1099 (endp (cdr (second p)))))
1100 C))
1101 (labels ((pair-key (p) (mock-spoly (first p) (second p) pred))
1102 (pair-order (x y) (funcall pred y x)))
1103 (merge 'list B C #'pair-order :key #'pair-key)))
1104
1105(defun gebauer-moeller-merge-pairs-smallest-lcm (B C pred ring)
1106 "The merging strategy based on the smallest-lcm (normal strategy); see
1107the documentation of BUCHBERGER-MERGE-PAIRS-SMALLEST-LCM."
1108 (declare (ignore ring))
1109 ;; Delete pairs consisting of two single monomials
1110 (labels ((pair-key (p) (monom-lcm (lm (first p))
1111 (lm (second p))))
1112 (pair-order (x y) (funcall pred y x)))
1113 (merge 'list B C #'pair-order :key #'pair-key)))
1114
1115
1116(defun gebauer-moeller-merge-pairs-use-smallest-degree (B C pred ring)
1117 "The merging strategy based on the smallest-lcm (normal strategy); see
1118the documentation of BUCHBERGER-MERGE-PAIRS-USE-SMALLEST-DEGREE."
1119 (declare (ignore ring))
1120 (declare (ignore pred))
1121 (labels ((pair-key (p) (total-degree
1122 (monom-lcm (lm (first p))
1123 (lm (second p)))))
1124 (pair-order (x y) (funcall #'< x y)))
1125 (merge 'list B C #'pair-order :key #'pair-key)))
1126
1127
1128(defun gebauer-moeller-merge-pairs-use-smallest-length (B C pred ring)
1129 (declare (ignore ring))
1130 "The merging strategy based on the smallest-lcm (normal strategy); see
1131the documentation of BUCHBERGER-MERGE-PAIRS-USE-SMALLEST-LENGTH."
1132 (declare (ignore pred))
1133 (labels ((pair-key (p) (+
1134 (length (first p))
1135 (length (second p))))
1136 (pair-order (x y) (funcall #'< x y)))
1137 (if B
1138 (merge 'list B C #'pair-order :key #'pair-key)
1139 (sort C #'pair-order :key #'pair-key))))
1140
1141
1142(defun gebauer-moeller-merge-pairs-use-smallest-coefficient-length (B C pred ring)
1143 "The merging strategy based on the smallest-lcm (normal strategy); see
1144the documentation of BUCHBERGER-MERGE-PAIRS-USE-SMALLEST-COEFFICIENT-LENGTH."
1145 (declare (ignore pred))
1146 (labels ((coefficient-length (poly)
1147 (apply #'+ (mapcar (ring-length ring) (mapcar #'cdr poly))))
1148 (pair-key (p) (+
1149 (coefficient-length (first p))
1150 (coefficient-length (second p))))
1151 (pair-order (x y) (funcall #'< x y)))
1152 (merge 'list B C #'pair-order :key #'pair-key)))
1153
1154(defun gebauer-moeller-set-pair-heuristic
1155 (method
1156 &aux (strategy-fn
1157 (ecase method
1158 (:minimal-mock-spoly #'gebauer-moeller-merge-pairs-use-mock-spoly)
1159 (:minimal-lcm #'gebauer-moeller-merge-pairs-smallest-lcm)
1160 (:minimal-total-degree #'gebauer-moeller-merge-pairs-use-smallest-degree)
1161 (:minimal-length #'gebauer-moeller-merge-pairs-use-smallest-length)
1162 (:minimal-coefficient-length #'gebauer-moeller-merge-pairs-use-smallest-coefficient-length))))
1163 (setf *gebauer-moeller-merge-pairs* strategy-fn))
1164
1165;;------------------------------------------------------------------------------------------------
1166;;
1167;; An implementation of the sugar strategy
1168;;
1169;;------------------------------------------------------------------------------------------------
1170
1171;; Calculate sugar of the S-polynomial of f and g
1172(defun spoly-sugar (f-with-sugar g-with-sugar ring)
1173 "Calculate the ``sugar'' of the S-polynomial of two polynomials with
1174sugar. A polynomial with sugar is simply a pair (P . SUGAR) where
1175SUGAR is an integer constant defined according to the following
1176algorighm: the sugar of sum or difference of two polynomials with
1177sugar is the MAX of the sugars of those two polynomials. The sugar of
1178a product of a term and a polynomial is the sum of the degree of the
1179term and the sugar of the polynomial. The idea is to accumulate sugar
1180as we perform the arithmetic operations, and that polynomials with
1181small (little) sugar should be given priority in the calculations.
1182Thus, the ``sugar strategy'' of the critical pair selection is to
1183select a pair with the smallest value of sugar of the resulting
1184S-polynomial."
1185 (let* ((f (poly-with-sugar-poly f-with-sugar))
1186 (g (poly-with-sugar-poly g-with-sugar))
1187 (lcm (monom-lcm (lm f) (lm g)))
1188 (m1 (monom/ lcm (lm f)))
1189 (m2 (monom/ lcm (lm g))))
1190 (let* ((c (funcall (ring-lcm ring) (lc f) (lc g)))
1191 (cf (funcall (ring-/ ring) c (lc f)))
1192 (cg (funcall (ring-/ ring) c (lc g))))
1193 (max
1194 (+ (term-sugar (cons m1 cf) ring) (poly-with-sugar-sugar f-with-sugar))
1195 (+ (term-sugar (cons m2 cg) ring) (poly-with-sugar-sugar g-with-sugar))))))
1196
1197;; Calculate the S-polynomial of f and g with sugar
1198(defun spoly-with-sugar (f-with-sugar g-with-sugar pred ring)
1199 "The S-polynomials of two polynomials with SUGAR strategy."
1200 (let* ((f (poly-with-sugar-poly f-with-sugar))
1201 (g (poly-with-sugar-poly g-with-sugar))
1202 (lcm (monom-lcm (lm f) (lm g)))
1203 (m1 (monom/ lcm (lm f)))
1204 (m2 (monom/ lcm (lm g))))
1205 (let* ((c (funcall (ring-lcm ring) (lc f) (lc g)))
1206 (cf (funcall (ring-/ ring) c (lc f)))
1207 (cg (funcall (ring-/ ring) c (lc g))))
1208 (poly-with-sugar-
1209 (term-times-poly-with-sugar (cons m1 cf) (poly-with-sugar-tail f-with-sugar) ring)
1210 (term-times-poly-with-sugar (cons m2 cg) (poly-with-sugar-tail g-with-sugar) ring)
1211 pred
1212 ring))))
1213
1214
1215(defun normal-form-with-sugar (f fl pred top-reduction-only ring)
1216 "Normal form of the polynomial with sugar F with respect to
1217a list of polynomials with sugar FL. Assumes that FL is not empty.
1218Parameters PRED, TOP-REDUCTION-ONLY and RING have the same meaning as in NORMAL-FORM.
1219The parameter SUGAR-LIMIT should be set to a positive integer. If the sugar limit of
1220the partial remainder exceeds SUGAR-LIMIT then the calculation is stopped and the partial
1221remainder is returned, although it is not fully reduced with respect to FL."
1222 (declare (optimize (speed 3) (safety 0)))
1223 ;; Loop invariant: c*f=sum ai*fi+r+p
1224 (do ((r (cons nil 0))
1225 (c (ring-unit ring))
1226 (division-count 0)
1227 (p f))
1228 ((or (poly-with-sugar-zerop p)
1229 (and top-reduction-only (not (poly-with-sugar-zerop r))))
1230 #+debug
1231 (progn
1232 (debug-cgb "~&~3T~d reduction~:p" division-count)
1233 (cond
1234 ((poly-with-sugar-zerop r)
1235 (debug-cgb " ---> 0"))))
1236 (values (poly-with-sugar-append (poly-with-sugar-nreverse r) p)
1237 division-count c))
1238 (declare (fixnum division-count) (integer c))
1239 (let ((g (find (poly-with-sugar-lm p) fl
1240 :test #'monom-divisible-by-p
1241 :key #'poly-with-sugar-lm)))
1242 (cond
1243 (g ;division possible
1244 (incf division-count)
1245 (let* ((gcd (funcall (ring-gcd ring) (poly-with-sugar-lc g) (poly-with-sugar-lc p)))
1246 (c1 (funcall (ring-/ ring) (poly-with-sugar-lc g) gcd))
1247 (c2 (funcall (ring-/ ring) (poly-with-sugar-lc p) gcd))
1248 (m (monom/ (poly-with-sugar-lm p) (poly-with-sugar-lm g))))
1249 ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
1250 (setf r (scalar-times-poly-with-sugar c1 r ring)
1251 c (funcall (ring-* ring) c c1)
1252 p (scalar-times-poly-with-sugar c1 p ring)
1253 p (poly-with-sugar-op p (cons m c2) g pred ring))))
1254 (t ;no division possible
1255 ;;move lt(p) to remainder
1256 (setf r (cons (cons (poly-with-sugar-lt p) (poly-with-sugar-poly r))
1257 (poly-with-sugar-sugar r))
1258 p (poly-with-sugar-tail p)))))))
1259
1260;; Parameter START is used when Grobner basis is calculated incrementally
1261;; and it signals that the elements of F from 0 to START-1 form a Grobner
1262;; basis already, so their pairs are not formed.
1263(defun buchberger-with-sugar (F-no-sugar pred start top-reduction-only ring
1264 &aux (s (1- (length F-no-sugar))) B M F)
1265 "Returns a Grobner basis of an ideal generated by a list of ordinary
1266polynomials F-NO-SUGAR. Adds initial sugar to the polynomials and
1267performs the Buchberger algorithm with ``sugar strategy''. It returns
1268an ordinary list of polynomials with no sugar. One of the most
1269effective algorithms for Grobner basis calculation."
1270 (declare (fixnum s))
1271 (unless (plusp s) (return-from buchberger-with-sugar F-no-sugar)) ;cut startup costs
1272 #+debug (debug-cgb "~&GROBNER BASIS - BUCHBERGER WITH SUGER STRATEGY ALGORITHM")
1273 #+debug (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
1274 #+grobner-check (when (plusp start)
1275 (grobner-test (subseq F-no-sugar 0 start)
1276 (subseq F-no-sugar 0 start) pred ring))
1277 (setf F (mapcar #'(lambda (x) (poly-add-sugar x ring)) F-no-sugar)
1278 B (nconc (makelist (list i j) (i 0 (1- start)) (j start s))
1279 (makelist (list i j) (i start (1- s)) (j (1+ i) s)))
1280 M (make-hash-table :test #'equal))
1281 ;;Initialize treated pairs
1282 (dotimes (i (1- start))
1283 (do ((j (1+ i) (1+ j))) ((>= j start))
1284 (setf (gethash (list i j) M) t)))
1285 (setf B (buchberger-with-sugar-sort-pairs B F M pred ring))
1286 (do ((G F))
1287 ((endp B)
1288 (setf G (mapcar #'poly-with-sugar-poly G))
1289 ;;#+grobner-check(grobner-test G F-no-sugar pred ring)
1290 #+debug(debug-cgb "~&GROBNER END")
1291 G)
1292 (let ((pair (pop B)))
1293 (cond
1294 ((Criterion-1-with-sugar pair G))
1295 ((Criterion-2-with-sugar pair G M))
1296 (t
1297 (let ((SP (normal-form-with-sugar
1298 (cons (spoly (poly-with-sugar-poly (elt G (first pair)))
1299 (poly-with-sugar-poly (elt G (second pair))) pred ring)
1300 (third pair))
1301 G pred top-reduction-only ring)))
1302 (cond
1303 ((poly-with-sugar-zerop SP))
1304 (t
1305 (setf s (1+ s)
1306 (poly-with-sugar-poly SP) (grobner-primitive-part
1307 (poly-with-sugar-poly SP) ring)
1308 (cdr (last G)) (list SP)
1309 B (buchberger-with-sugar-merge-pairs B (makelist (list i s)
1310 (i 0 (1- s)))
1311 G M pred ring))
1312 #+debug (debug-cgb "~&[~d] Polynomials: ~d; Pairs left: ~d; Treated pairs: ~d"
1313 (third pair)
1314 (length G) (length B)
1315 (hash-table-count M)))))))
1316 (setf (cddr pair) nil
1317 (gethash pair M) t))))
1318
1319
1320(defun buchberger-with-sugar-merge-pairs (B C G M pred ring)
1321 "Merges lists of critical pairs. It orders pairs according to
1322increasing sugar, with ties broken by smaller MOCK-SPOLY value of the
1323two polynomials. In this function B must already be sorted."
1324 ;;Remove pairs consisting of two monomials
1325 ;;These pairs produce zero S-polynomial
1326 (setf C (delete-if
1327 #'(lambda (p)
1328 (when (and (endp (cdr (poly-with-sugar-poly (elt G (first p)))))
1329 (endp (cdr (poly-with-sugar-poly (elt G (second p))))))
1330 ;;Record the pair as treated
1331 (setf (gethash p M) t)
1332 t))
1333 C))
1334 (labels ((pair-sugar (p)
1335 (spoly-sugar (elt G (first p)) (elt G (second p)) ring))
1336 (pair-mock-spoly (p)
1337 (mock-spoly (poly-with-sugar-poly (elt G (first p)))
1338 (poly-with-sugar-poly (elt G (second p))) pred))
1339 (pair-order (p q)
1340 (let ((sugar-p (third p))
1341 (sugar-q (third q)))
1342 (cond
1343 ((< sugar-p sugar-q) t)
1344 ((> sugar-p sugar-q) nil)
1345 (t
1346 (funcall pred (fourth q) (fourth p)))))))
1347 ;;Add sugar and mock S-polynomial as the third and fourth element
1348 ;;of new pair in C
1349 (dolist (p C)
1350 (setf (cdr (last p))
1351 (list (pair-sugar p) (pair-mock-spoly p))))
1352 (merge 'list B (sort C #'pair-order) #'pair-order)))
1353
1354(defun buchberger-with-sugar-sort-pairs (C G M pred ring)
1355 "Sorts critical pairs C according to sugar strategy."
1356 (buchberger-with-sugar-merge-pairs nil C G M pred ring))
1357
1358(defun Criterion-1-with-sugar (pair G &aux (i (first pair)) (j (second pair)))
1359 "An implementation of Buchberger's first criterion for polynomials
1360with sugar. See the documentation of BUCHBERGER-CRITERION-1."
1361 (monom-rel-prime (poly-with-sugar-lm (elt G i))
1362 (poly-with-sugar-lm (elt G j))))
1363
1364(defun Criterion-2-with-sugar (pair G M &aux (i (first pair)) (j (second pair)))
1365 "An implementation of Buchberger's first criterion for polynomials
1366with sugar. See the documentation of BUCHBERGER-CRITERION-2."
1367 (labels ((make-pair (u v) (if (< u v) (list u v) (list v u))))
1368 (dotimes (k (length G) nil)
1369 (when (and (/= k i)
1370 (/= k j)
1371 (monom-divides-p
1372 (poly-with-sugar-lm (elt G k))
1373 (monom-lcm (poly-with-sugar-lm (elt G i))
1374 (poly-with-sugar-lm (elt G j))))
1375 (gethash (make-pair i k) M)
1376 (gethash (make-pair k j) M))
1377 #+debug(when *grobner-debug* (format t ":B2"))
1378 (return-from Criterion-2-with-sugar t)))))
1379
1380
1381;;----------------------------------------------------------------
1382;; An implementation of the algorithm of Gebauer and Moeller,
1383;; as described in the book of Becker-Weispfenning, p. 232
1384;; with sugar strategy
1385;;----------------------------------------------------------------
1386(defun gebauer-moeller-with-sugar (F pred start top-reduction-only ring
1387 &aux B G F1)
1388 "Compute Grobner basis by using the algorithm of Gebauer and Moeller.
1389This algorithm is described as BUCHBERGERNEW2 in the book by
1390Becker-Weispfenning entitled ``Grobner Bases''"
1391 ;;cut startup costs if (length F)<=1
1392 (cond
1393 ((endp F) (return-from gebauer-moeller-with-sugar nil))
1394 ((endp (cdr F))
1395 (return-from gebauer-moeller-with-sugar (list (grobner-primitive-part (car F) ring)))))
1396 #+debug (debug-cgb "~&GROBNER BASIS - GEBAUER MOELLER ALGORITHM WITH SUGAR STRATEGY")
1397 #+debug (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
1398 #+grobner-check (when (plusp start)
1399 (grobner-test (subseq F 0 start) (subseq F 0 start)
1400 pred ring))
1401 (setf B nil
1402 G (subseq F 0 start)
1403 F1 (subseq F start))
1404 ;;Initialize sugar
1405 (setf G (mapcar #'(lambda (x) (poly-add-sugar x ring)) G)
1406 F1 (mapcar #'(lambda (x) (poly-add-sugar x ring)) F1))
1407 (do () ((endp F1))
1408 (multiple-value-setq (G B)
1409 (update-with-sugar G B (grobner-primitive-part-with-sugar (pop F1) ring) pred ring)))
1410 (do ()
1411 ((endp B))
1412 (let* ((pair (pop B))
1413 (g1 (first pair))
1414 (g2 (second pair))
1415 (h (normal-form-with-sugar (spoly-with-sugar g1 g2 pred ring)
1416 G pred top-reduction-only ring)))
1417 (cond
1418 ((poly-with-sugar-zerop h))
1419 (t
1420 (setf h (grobner-primitive-part-with-sugar h ring))
1421 (multiple-value-setq (G B)
1422 (update-with-sugar G B h pred ring))
1423 #+debug (debug-cgb "~&[[~d]] Polynomials: ~d; Pairs left: ~d~%"
1424 (poly-with-sugar-sugar h)
1425 (length G) (length B))
1426 ))))
1427 (setf G (mapcar #'poly-with-sugar-poly G))
1428 ;;#+grobner-check(grobner-test G F pred ring)
1429 #+debug(debug-cgb "~&GROBNER END")
1430 G)
1431
1432(defun update-with-sugar (G B h pred ring &aux C D E B-new G-new pair)
1433 "An implementation of the auxillary UPDATE algorithm used by the
1434Gebauer-Moeller algorithm. G is a list of polynomials, B is a list of
1435critical pairs and H is a new polynomial which possibly will be added
1436to G. The naming conventions used are very close to the one used in
1437the book of Becker-Weispfenning. Operates on polynomials with sugar."
1438 (setf C G D nil)
1439 (do (g1) ((endp C))
1440 (setf g1 (pop C))
1441 (when (or (monom-rel-prime (poly-with-sugar-lm h) (poly-with-sugar-lm g1))
1442 (and
1443 (not (some
1444 #'(lambda (g2) (monom-divides-p (monom-lcm (poly-with-sugar-lm h)
1445 (poly-with-sugar-lm g2))
1446 (monom-lcm (poly-with-sugar-lm h)
1447 (poly-with-sugar-lm g1))))
1448 C))
1449 (not (some
1450 #'(lambda (g2) (monom-divides-p (monom-lcm (poly-with-sugar-lm h)
1451 (poly-with-sugar-lm g2))
1452 (monom-lcm (poly-with-sugar-lm h)
1453 (poly-with-sugar-lm g1))))
1454 D))))
1455 (push g1 D)))
1456 (setf E nil)
1457 (do (g1) ((endp D))
1458 (setf g1 (pop D))
1459 (unless (monom-rel-prime (poly-with-sugar-lm h) (poly-with-sugar-lm g1))
1460 (push g1 E)))
1461 (setf B-new nil)
1462 (do (g1 g2) ((endp B))
1463 (setf pair (pop B)
1464 g1 (first pair)
1465 g2 (second pair))
1466 (when (or (not (monom-divides-p (poly-with-sugar-lm h) (monom-lcm (poly-with-sugar-lm g1)
1467 (poly-with-sugar-lm g2))))
1468 (equal (monom-lcm (poly-with-sugar-lm g1) (poly-with-sugar-lm h))
1469 (monom-lcm (poly-with-sugar-lm g1) (poly-with-sugar-lm g2)))
1470 (equal (monom-lcm (poly-with-sugar-lm h) (poly-with-sugar-lm g2))
1471 (monom-lcm (poly-with-sugar-lm g1) (poly-with-sugar-lm g2))))
1472 (setf B-new (gebauer-moeller-with-sugar-merge-pairs B-new (list (list g1 g2)) pred ring))))
1473 (dolist (g3 E)
1474 (setf B-new (gebauer-moeller-with-sugar-merge-pairs B-new (list (list h g3)) pred ring)))
1475 (setf G-new nil)
1476 (do (g1) ((endp G))
1477 (setf g1 (pop G))
1478 (unless (monom-divides-p (poly-with-sugar-lm h) (poly-with-sugar-lm g1))
1479 (push g1 G-new)))
1480 (push h G-new)
1481 (values G-new B-new))
1482
1483(defun gebauer-moeller-with-sugar-merge-pairs (B C pred ring)
1484 "Merges lists of critical pairs. It orders pairs according to
1485increasing sugar, with ties broken by smaller lcm of head monomials In
1486this function B must already be sorted. Operates on polynomials with sugar."
1487 (setf C (delete-if
1488 #'(lambda (p) (and (endp (cdr (first p)))
1489 (endp (cdr (second p)))))
1490 C))
1491 (labels ((pair-sugar (p)
1492 (spoly-sugar (first p) (second p) ring))
1493 (pair-mock-spoly (p)
1494 (mock-spoly (poly-with-sugar-poly (first p))
1495 (poly-with-sugar-poly (second p)) pred))
1496 (pair-order (p q)
1497 (let ((sugar-p (pair-sugar p))
1498 (sugar-q (pair-sugar q)))
1499 (cond
1500 ((< sugar-p sugar-q) t)
1501 ((> sugar-p sugar-q) nil)
1502 (t
1503 (funcall pred (pair-mock-spoly q) (pair-mock-spoly p)))))))
1504 (merge 'list B (sort C #'pair-order) #'pair-order)))
1505
1506(defun grobner-primitive-part-with-sugar (h ring)
1507 (cons (grobner-primitive-part (car h) ring) (cdr h)))
1508
1509;;----------------------------------------------------------------
1510;;
1511;; An implementation of the algorithms using product caching
1512;;
1513;;----------------------------------------------------------------
1514(defun cached-normal-form (f fl pred top-reduction-only ring cache)
1515 "This version of NORMAL-FORM uses a cache which stores products of monomials
1516and polynomials found in FL. The cache is a hash table storing pairs (G . G-HASH-TABLE)
1517where G is a polynomial in FL serving as a key. The hash table G-HASH-TABLE is a list
1518of pairs (MONOMIAL . PRODUCT-POLY) where PRODUCT-POLY is (POLY* MONOMIAL G)."
1519
1520 (declare (optimize (speed 3) (safety 0)))
1521 ;; Loop invariant: c*f=sum ai*fi+r+p
1522 (do (r (c (ring-unit ring))
1523 (division-count 0)
1524 (p f))
1525 ((or (endp p)
1526 (endp fl)
1527 (and top-reduction-only r))
1528 #+debug(progn
1529 (debug-cgb "~&~3T~d reduction~:p" division-count)
1530 (when (endp r)
1531 (debug-cgb " ---> 0")))
1532 (values (append (nreverse r) p)
1533 division-count c))
1534 (declare (fixnum division-count) (integer c))
1535 (let ((i (position (lm p) fl :test #'monom-divisible-by-p :key #'lm)))
1536 (cond
1537 (i ;division possible
1538 (incf division-count)
1539 (let* ((g (nth i fl))
1540 (hash-tab (nth i cache))
1541 (gcd (funcall (ring-gcd ring) (lc g) (lc p)))
1542 (c1 (funcall (ring-/ ring) (lc g) gcd))
1543 (c2 (funcall (ring-/ ring) (lc p) gcd))
1544 (m (monom/ (lm p) (lm g))))
1545 ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
1546 (setf r (scalar-times-poly c1 r ring)
1547 c (funcall (ring-* ring) c c1)
1548 p (cached-grobner-op c2 c1 m p g pred ring hash-tab))))
1549 (t ;no division possible
1550 (setf r (cons (lt p) r) ;move lt(p) to remainder
1551 p (rest p)))))))
1552
1553(defun cached-grobner-op (c1 c2 m A B pred ring cache)
1554 "Returns the same result as (GROBNER-OP C1 C2 M A B PRED RING) but
1555it used variable CACHE which is described in CACHED-NORMAL-FORM for
1556efficient calculation of products of monomials with polynomials."
1557 (poly- (scalar-times-poly c2 A)
1558 (scalar-times-poly c1 (cached-monom-times-poly m B cache) ring)
1559 pred ring))
1560
1561(defun cached-buchberger (F pred start top-reduction-only ring
1562 &aux (s (1- (length F))) B M CACHE)
1563 "It returns the same result as BUCHBERGER but it uses caching
1564of products of monomials with polynomials."
1565 (declare (fixnum s))
1566 (unless (plusp s) (return-from cached-buchberger F)) ;cut startup costs
1567 #+debug (debug-cgb "~&GROBNER BASIS - CACHED BUCHBERGER ALGORITHM")
1568 #+debug (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
1569 #+grobner-check (when (plusp start)
1570 (grobner-test (subseq F 0 start) (subseq F 0 start) pred ring))
1571 (setf B (nconc (makelist (list i j) (i 0 (1- start)) (j start s))
1572 (makelist (list i j) (i start (1- s)) (j (1+ i) s)))
1573 M (make-hash-table :test #'equal)
1574 B (buchberger-sort-pairs B F pred ring))
1575 ;;Initialize cache
1576 (dotimes (k (1+ s)) (push (make-hash-table :test #'equal) CACHE))
1577 ;;Initialize treated pairs
1578 (dotimes (i (1- start))
1579 (do ((j (1+ i) (1+ j))) ((>= j start))
1580 (setf (gethash (list i j) M) t)))
1581 (do ((G F))
1582 ((endp B)
1583 #+grobner-check(grobner-test G F pred ring)
1584 #+debug(debug-cgb "~&GROBNER END")
1585 G)
1586 (let ((pair (pop B)))
1587 (unless (or (Criterion-1 pair G)
1588 (Criterion-2 pair G M))
1589 (let ((SP (cached-normal-form (cached-spoly (elt G (first pair))
1590 (elt G (second pair))
1591 pred
1592 ring
1593 (elt CACHE (first pair))
1594 (elt CACHE (second pair)))
1595 G pred top-reduction-only ring CACHE)))
1596 (unless (poly-zerop SP) ;S-polynomial ---> 0;
1597 (setf s (1+ s)
1598 SP (grobner-primitive-part SP ring)
1599 (cdr (last G)) (list SP)
1600 (cdr (last Cache)) (list (make-hash-table :test #'equal))
1601 B (funcall *buchberger-merge-pairs* B (makelist (list i s)
1602 (i 0 (1- s)))
1603 G pred ring))
1604 #+debug (debug-cgb "~&Polynomials: ~d; Pairs left: ~d; Treated pairs: ~d"
1605 (length G) (length B)
1606 (hash-table-count M)))))
1607 (setf (gethash pair M) t))))
1608
1609(defun cached-spoly (f g pred ring cache-f cache-g)
1610 "Uses CACHE to fetch products of monomials and polynomials."
1611 (declare (optimize (speed 3) (safety 0)))
1612 (let* ((lcm (monom-lcm (lm f) (lm g)))
1613 (m1 (monom/ lcm (lm f)))
1614 (m2 (monom/ lcm (lm g))))
1615 (let* ((c (funcall (ring-lcm ring) (lc f) (lc g)))
1616 (cf (funcall (ring-/ ring) c (lc f)))
1617 (cg (funcall (ring-/ ring) c (lc g))))
1618 (poly-
1619 (scalar-times-poly cf (cached-monom-times-poly m1 f cache-f) ring)
1620 (scalar-times-poly cg (cached-monom-times-poly m2 g cache-g) ring)
1621 pred
1622 ring))))
1623
1624
1625(defun cached-monom-times-poly (m f cache)
1626 (let ((prod (gethash m cache)))
1627 (cond
1628 (prod #+debug(debug-cgb "c") prod)
1629 (t
1630 #+debug(debug-cgb ".")
1631 (setf (gethash m cache) (monom-times-poly m f))))))
1632
1633;;------------------------------------------------------------------------------------------------
1634;;
1635;; An implementation of the sugar strategy with product caching
1636;;
1637;;------------------------------------------------------------------------------------------------
1638(defun cached-normal-form-with-sugar (f fl pred top-reduction-only ring cache)
1639 "Normal form of the polynomial with sugar F with respect to
1640a list of polynomials with sugar FL. Assumes that FL is not empty.
1641Parameters PRED, TOP-REDUCTION-ONLY and RING have the same meaning as in NORMAL-FORM.
1642The parameter SUGAR-LIMIT should be set to a positive integer. If the sugar limit of
1643the partial remainder exceeds SUGAR-LIMIT then the calculation is stopped and the partial
1644remainder is returned, although it is not fully reduced with respect to FL."
1645 (declare (optimize (speed 3) (safety 0)))
1646 ;; Loop invariant: c*f=sum ai*fi+r+p
1647 (do ((r (cons nil 0))
1648 (c (ring-unit ring))
1649 (division-count 0)
1650 (p f))
1651 ((or (poly-with-sugar-zerop p)
1652 (and top-reduction-only (not (poly-with-sugar-zerop r))))
1653 #+debug
1654 (progn
1655 (debug-cgb "~&~3T~d reduction~:p" division-count)
1656 (cond
1657 ((poly-with-sugar-zerop r)
1658 (debug-cgb " ---> 0"))))
1659 (values (poly-with-sugar-append (poly-with-sugar-nreverse r) p)
1660 division-count c))
1661 (declare (fixnum division-count) (integer c))
1662 (let ((i (position (poly-with-sugar-lm p) fl
1663 :test #'monom-divisible-by-p
1664 :key #'poly-with-sugar-lm)))
1665 (cond
1666 (i ;division possible
1667 (incf division-count)
1668 (let* ((g (nth i fl))
1669 (hash-tab (nth i cache))
1670 (gcd (funcall (ring-gcd ring) (poly-with-sugar-lc g) (poly-with-sugar-lc p)))
1671 (c1 (funcall (ring-/ ring) (poly-with-sugar-lc g) gcd))
1672 (c2 (funcall (ring-/ ring) (poly-with-sugar-lc p) gcd))
1673 (m (monom/ (poly-with-sugar-lm p) (poly-with-sugar-lm g))))
1674 ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
1675 (setf r (scalar-times-poly-with-sugar c1 r ring)
1676 c (funcall (ring-* ring) c c1)
1677 p (scalar-times-poly-with-sugar c1 p ring)
1678 p (cached-poly-with-sugar-op p c2 m g pred ring hash-tab))))
1679 (t ;no division possible
1680 ;;move lt(p) to remainder
1681 (setf r (cons (cons (poly-with-sugar-lt p) (poly-with-sugar-poly r))
1682 (poly-with-sugar-sugar r))
1683 p (poly-with-sugar-tail p)))))))
1684
1685;; Parameter START is used when Grobner basis is calculated incrementally
1686;; and it signals that the elements of F from 0 to START-1 form a Grobner
1687;; basis already, so their pairs are not formed.
1688(defun cached-buchberger-with-sugar (F-no-sugar pred start top-reduction-only ring
1689 &aux (s (1- (length F-no-sugar))) B M F Cache)
1690 "Returns a Grobner basis of an ideal generated by a list of ordinary
1691polynomials F-NO-SUGAR. Adds initial sugar to the polynomials and
1692performs the Buchberger algorithm with ``sugar strategy''. It returns
1693an ordinary list of polynomials with no sugar. One of the most
1694effective algorithms for Grobner basis calculation."
1695 (declare (fixnum s))
1696 (unless (plusp s) (return-from cached-buchberger-with-sugar F-no-sugar)) ;cut startup costs
1697 #+debug (debug-cgb "~&GROBNER BASIS - BUCHBERGER WITH SUGER STRATEGY ALGORITHM")
1698 #+debug (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
1699 #+grobner-check (when (plusp start)
1700 (grobner-test (subseq F-no-sugar 0 start)
1701 (subseq F-no-sugar 0 start) pred ring))
1702 (setf F (mapcar #'(lambda (x) (poly-add-sugar x ring)) F-no-sugar)
1703 B (nconc (makelist (list i j) (i 0 (1- start)) (j start s))
1704 (makelist (list i j) (i start (1- s)) (j (1+ i) s)))
1705 M (make-hash-table :test #'equal))
1706 ;;Initialize cache
1707 (dotimes (k (1+ s)) (push (make-hash-table :test #'equal) Cache))
1708 ;;Initialize treated pairs
1709 (dotimes (i (1- start))
1710 (do ((j (1+ i) (1+ j))) ((>= j start))
1711 (setf (gethash (list i j) M) t)))
1712 (setf B (buchberger-with-sugar-sort-pairs B F M pred ring))
1713 (do ((G F))
1714 ((endp B)
1715 (setf G (mapcar #'poly-with-sugar-poly G))
1716 ;;#+grobner-check(grobner-test G F-no-sugar pred ring)
1717 #+debug(debug-cgb "~&GROBNER END")
1718 G)
1719 (let ((pair (pop B)))
1720 (cond
1721 ((Criterion-1-with-sugar pair G))
1722 ((Criterion-2-with-sugar pair G M))
1723 (t
1724 (let ((SP (cached-normal-form-with-sugar
1725 (cons (cached-spoly (poly-with-sugar-poly (elt G (first pair)))
1726 (poly-with-sugar-poly (elt G (second pair)))
1727 pred ring
1728 (elt Cache (first pair))
1729 (elt Cache (second pair)))
1730 (third pair))
1731 G pred top-reduction-only ring Cache)))
1732 (cond
1733 ((poly-with-sugar-zerop SP))
1734 (t
1735 (setf s (1+ s)
1736 (poly-with-sugar-poly SP) (grobner-primitive-part
1737 (poly-with-sugar-poly SP) ring)
1738 (cdr (last G)) (list SP)
1739 (cdr (last Cache)) (list (make-hash-table :test #'equal))
1740 B (buchberger-with-sugar-merge-pairs B (makelist (list i s)
1741 (i 0 (1- s)))
1742 G M pred ring))
1743 #+debug (debug-cgb "~&[~d] Polynomials: ~d; Pairs left: ~d; Treated pairs: ~d"
1744 (third pair)
1745 (length G) (length B)
1746 (hash-table-count M)))))))
1747 (setf (cddr pair) nil
1748 (gethash pair M) t))))
1749
1750;; Implement f-term*g
1751(defun cached-poly-with-sugar-op (f c m g pred ring Cache)
1752 (poly-with-sugar- f (scalar-times-poly-with-sugar c (cached-monom-times-poly-with-sugar m g Cache) ring) pred ring))
1753
1754(defun cached-monom-times-poly-with-sugar (m f Cache)
1755 (cons
1756 (cached-monom-times-poly m (poly-with-sugar-poly f) Cache)
1757 (+ (poly-with-sugar-sugar f) (monom-sugar m))))
1758@
1759
1760
17611.7
1762log
1763@*** empty log message ***
1764@
1765text
1766@a86 3
1767(eval-when (compile)
1768 (proclaim '(inline grobner)))
1769
1770a134 3
1771(eval-when (compile)
1772 (proclaim '(inline grobner-primitive-part grobner-content)))
1773
1774a317 2
1775(eval-when (compile)
1776 (proclaim '(inline buchberger-sort-pairs)))
1777a346 6
1778(eval-when (compile)
1779 (proclaim '(inline
1780 buchberger-merge-pairs-use-mock-spoly
1781 buchberger-merge-pairs-smallest-lcm
1782 buchberger-merge-pairs-use-smallest-degree)))
1783
1784a441 3
1785(eval-when (compile)
1786 (proclaim '(inline Criterion-1 Criterion-2)))
1787
1788a524 3
1789(eval-when (compile)
1790 (proclaim '(inline monom-depends-p term-depends-p poly-depends-p)))
1791
1792a1300 3
1793(eval-when (compile)
1794 (proclaim '(inline Criterion-1-with-sugar Criterion-2-with-sugar)))
1795
1796@
1797
1798
17991.6
1800log
1801@*** empty log message ***
1802@
1803text
1804@a974 1
1805 (declare (ignore sorted))
1806@
1807
1808
18091.5
1810log
1811@*** empty log message ***
1812@
1813text
1814@a1357 1
1815 (declare (ignore sorted))
1816@
1817
1818
18191.4
1820log
1821@*** empty log message ***
1822@
1823text
1824@d57 2
1825a58 2
1826;;(proclaim '(optimize (speed 0) (debug 3)))
1827(proclaim '(optimize (speed 3) (debug 0)))
1828@
1829
1830
18311.3
1832log
1833@*** empty log message ***
1834@
1835text
1836@d57 2
1837a58 1
1838(proclaim '(optimize (speed 0) (debug 3)))
1839@
1840
1841
18421.2
1843log
1844@*** empty log message ***
1845@
1846text
1847@d57 1
1848a57 2
1849(eval-when (compile)
1850 (proclaim '(optimize (safety 0) (speed 3))))
1851@
1852
1853
18541.1
1855log
1856@Initial revision
1857@
1858text
1859@d3 1
1860a3 1
1861 $Id: grobner.lisp,v 1.288 1998/01/18 09:06:58 marek Exp $
1862d15 1
1863a15 1
1864 "POLY-WITH-SUGAR" "LISP")
1865@
Note: See TracBrowser for help on using the repository browser.