Changeset 65 for branches/f4grobner/grobner.lisp

Timestamp:

2015-06-05T11:27:08-07:00 (9 years ago)

Author:

Marek Rychlik

Message:

* empty log message *

File:

• branches/f4grobner/grobner.lisp (modified) (3 diffs)

branches/f4grobner/grobner.lisp

@@ -293 +293 @@
   `(when $poly_grobner_debug (format *terminal-io* ,@args)))
 
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;
-;; An implementation of Grobner basis
-;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-
-
-
-
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;
-;; An implementation of the division algorithm
-;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-
-(defun grobner-op (ring c1 c2 m f g)
-  "Returns C2*F-C1*M*G, where F and G are polynomials M is a monomial.
-Assume that the leading terms will cancel."
-  #+grobner-check(funcall (ring-zerop ring)
-                          (funcall (ring-sub ring)
-                                   (funcall (ring-mul ring) c2 (poly-lc f))
-                                   (funcall (ring-mul ring) c1 (poly-lc g))))
-  #+grobner-check(monom-equal-p (poly-lm f) (monom-mul m (poly-lm g)))
-  ;; Note that we can drop the leading terms of f ang g
-  (poly-sub ring
-            (scalar-times-poly-1 ring c2 f)
-            (scalar-times-poly-1 ring c1 (monom-times-poly m g))))
-
-(defun poly-pseudo-divide (ring f fl)
-  "Pseudo-divide a polynomial F by the list of polynomials FL. Return
-multiple values. The first value is a list of quotients A.  The second
-value is the remainder R. The third argument is a scalar coefficient
-C, such that C*F can be divided by FL within the ring of coefficients,
-which is not necessarily a field. Finally, the fourth value is an
-integer count of the number of reductions performed.  The resulting
-objects satisfy the equation: C*F= sum A[i]*FL[i] + R."
-  (declare (type poly f) (list fl))
-  (do ((r (make-poly-zero))
-       (c (funcall (ring-unit ring)))
-       (a (make-list (length fl) :initial-element (make-poly-zero)))
-       (division-count 0)
-       (p f))
-      ((poly-zerop p)
-       (debug-cgb "~&~3T~d reduction~:p" division-count)
-       (when (poly-zerop r) (debug-cgb " ---> 0"))
-       (values (mapcar #'poly-nreverse a) (poly-nreverse r) c division-count))
-    (declare (fixnum division-count))
-    (do ((fl fl (rest fl))                              ;scan list of divisors
-         (b a (rest b)))
-        ((cond
-          ((endp fl)                                    ;no division occurred
-           (push (poly-lt p) (poly-termlist r))         ;move lt(p) to remainder
-           (setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
-           (pop (poly-termlist p))                      ;remove lt(p) from p
-           t)
-          ((monom-divides-p (poly-lm (car fl)) (poly-lm p)) ;division occurred
-           (incf division-count)
-           (multiple-value-bind (gcd c1 c2)
-               (funcall (ring-ezgcd ring) (poly-lc (car fl)) (poly-lc p))
-             (declare (ignore gcd))
-             (let ((m (monom-div (poly-lm p) (poly-lm (car fl)))))
-               ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
-               (mapl #'(lambda (x)
-                         (setf (car x) (scalar-times-poly ring c1 (car x))))
-                     a)
-               (setf r (scalar-times-poly ring c1 r)
-                     c (funcall (ring-mul ring) c c1)
-                     p (grobner-op ring c2 c1 m p (car fl)))
-               (push (make-term m c2) (poly-termlist (car b))))
-             t)))))))
-
-(defun poly-exact-divide (ring f g)
-  "Divide a polynomial F by another polynomial G. Assume that exact division
-with no remainder is possible. Returns the quotient."
-  (declare (type poly f g))
-  (multiple-value-bind (quot rem coeff division-count)
-      (poly-pseudo-divide ring f (list g))
-    (declare (ignore division-count coeff)
-             (list quot)
-             (type poly rem)
-             (type fixnum division-count))
-    (unless (poly-zerop rem) (error "Exact division failed."))
-    (car quot)))
-
-
-
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;
-;; An implementation of the normal form
-;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-
-(defun normal-form-step (ring fl p r c division-count
-                         &aux (g (find (poly-lm p) fl
-                                       :test #'monom-divisible-by-p
-                                       :key #'poly-lm)))
-  (cond
-   (g                                   ;division possible
-    (incf division-count)
-    (multiple-value-bind (gcd cg cp)
-        (funcall (ring-ezgcd ring) (poly-lc g) (poly-lc p))
-      (declare (ignore gcd))
-      (let ((m (monom-div (poly-lm p) (poly-lm g))))
-        ;; Multiply the equation c*f=sum ai*fi+r+p by cg.
-        (setf r (scalar-times-poly ring cg r)
-              c (funcall (ring-mul ring) c cg)
-              ;; p := cg*p-cp*m*g
-              p (grobner-op ring cp cg m p g))))
-    (debug-cgb "/"))
-   (t                                                   ;no division possible
-    (push (poly-lt p) (poly-termlist r))                ;move lt(p) to remainder
-    (setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
-    (pop (poly-termlist p))                             ;remove lt(p) from p
-    (debug-cgb "+")))
-  (values p r c division-count))
-
-;; Merge it sometime with poly-pseudo-divide
-(defun normal-form (ring f fl &optional (top-reduction-only $poly_top_reduction_only))
-  ;; Loop invariant: c*f0=sum ai*fi+r+f, where f0 is the initial value of f
-  #+grobner-check(when (null fl) (warn "normal-form: empty divisor list."))
-  (do ((r (make-poly-zero))
-       (c (funcall (ring-unit ring)))
-       (division-count 0))
-      ((or (poly-zerop f)
-           ;;(endp fl)
-           (and top-reduction-only (not (poly-zerop r))))
-       (progn
-         (debug-cgb "~&~3T~d reduction~:p" division-count)
-         (when (poly-zerop r)
-           (debug-cgb " ---> 0")))
-       (setf (poly-termlist f) (nreconc (poly-termlist r) (poly-termlist f)))
-       (values f c division-count))
-    (declare (fixnum division-count)
-             (type poly r))
-    (multiple-value-setq (f r c division-count)
-      (normal-form-step ring fl f r c division-count))))
+
 
 
@@ -440 +305 @@
 ;;
 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-
-(defun buchberger-criterion (ring g)
-  "Returns T if G is a Grobner basis, by using the Buchberger
-criterion: for every two polynomials h1 and h2 in G the S-polynomial
-S(h1,h2) reduces to 0 modulo G."
-  (every
-   #'poly-zerop
-   (makelist (normal-form ring (spoly ring (elt g i) (elt g j)) g nil)
-             (i 0 (- (length g) 2))
-             (j (1+ i) (1- (length g))))))
 
 (defun grobner-test (ring g f)
@@ -467 +322 @@
   t)
 
-
-
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;
-;; Pair queue implementation
-;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-
-(defun sugar-pair-key (p q &aux (lcm (monom-lcm (poly-lm p) (poly-lm q)))
-                                (d (monom-sugar lcm)))
-  "Returns list (S LCM-TOTAL-DEGREE) where S is the sugar of the S-polynomial of
-polynomials P and Q, and LCM-TOTAL-DEGREE is the degree of is LCM(LM(P),LM(Q))."
-  (declare (type poly p q) (type monom lcm) (type fixnum d))
-  (cons (max 
-         (+  (- d (monom-sugar (poly-lm p))) (poly-sugar p))
-         (+  (- d (monom-sugar (poly-lm q))) (poly-sugar q)))
-        lcm))
-
-(defstruct (pair
-            (:constructor make-pair (first second
-                                           &aux
-                                           (sugar (car (sugar-pair-key first second)))
-                                           (division-data nil))))
-  (first nil :type poly)
-  (second nil :type poly)
-  (sugar 0 :type fixnum)
-  (division-data nil :type list))
-  
-;;(defun pair-sugar (pair &aux (p (pair-first pair)) (q (pair-second pair)))
-;;  (car (sugar-pair-key p q)))
-
-(defun sugar-order (x y)
-  "Pair order based on sugar, ties broken by normal strategy."
-  (declare (type cons x y))
-  (or (< (car x) (car y))
-      (and (= (car x) (car y))
-           (< (monom-total-degree (cdr x))
-              (monom-total-degree (cdr y))))))
-
-(defvar *pair-key-function* #'sugar-pair-key
-  "Function that, given two polynomials as argument, computed the key
-in the pair queue.")
-
-(defvar *pair-order* #'sugar-order
-  "Function that orders the keys of pairs.")
-
-(defun make-pair-queue ()
-  "Constructs a priority queue for critical pairs."
-  (make-priority-queue
-   :element-type 'pair
-   :element-key #'(lambda (pair) (funcall *pair-key-function* (pair-first pair) (pair-second pair)))
-   :test *pair-order*))
-
-(defun pair-queue-initialize (pq f start
-                              &aux
-                              (s (1- (length f)))
-                              (b (nconc (makelist (make-pair (elt f i) (elt f j))
-                                                 (i 0 (1- start)) (j start s))
-                                        (makelist (make-pair (elt f i) (elt f j))
-                                                 (i start (1- s)) (j (1+ i) s)))))
-  "Initializes the priority for critical pairs. F is the initial list of polynomials.
-START is the first position beyond the elements which form a partial
-grobner basis, i.e. satisfy the Buchberger criterion."
-  (declare (type priority-queue pq) (type fixnum start))
-  (dolist (pair b pq)
-    (priority-queue-insert pq pair)))
-
-(defun pair-queue-insert (b pair)
-  (priority-queue-insert b pair))
-
-(defun pair-queue-remove (b)
-  (priority-queue-remove b))
-
-(defun pair-queue-size (b)
-  (priority-queue-size b))
-
-(defun pair-queue-empty-p (b)
-  (priority-queue-empty-p b))
-
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;
-;; Buchberger Algorithm Implementation
-;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-
-(defun buchberger (ring f start &optional (top-reduction-only $poly_top_reduction_only))
-  "An implementation of the Buchberger algorithm. Return Grobner basis
-of the ideal generated by the polynomial list F.  Polynomials 0 to
-START-1 are assumed to be a Grobner basis already, so that certain
-critical pairs will not be examined. If TOP-REDUCTION-ONLY set, top
-reduction will be preformed. This function assumes that all polynomials
-in F are non-zero."
-  (declare (type fixnum start))
-  (when (endp f) (return-from buchberger f)) ;cut startup costs
-  (debug-cgb "~&GROBNER BASIS - BUCHBERGER ALGORITHM")
-  (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
-  #+grobner-check  (when (plusp start)
-                     (grobner-test ring (subseq f 0 start) (subseq f 0 start)))
-  ;;Initialize critical pairs
-  (let ((b (pair-queue-initialize (make-pair-queue)
-                                  f start))
-        (b-done (make-hash-table :test #'equal)))
-    (declare (type priority-queue b) (type hash-table b-done))
-    (dotimes (i (1- start))
-      (do ((j (1+ i) (1+ j))) ((>= j start))
-        (setf (gethash (list (elt f i) (elt f j)) b-done) t)))
-    (do ()
-        ((pair-queue-empty-p b)
-         #+grobner-check(grobner-test ring f f)
-         (debug-cgb "~&GROBNER END")
-         f)
-      (let ((pair (pair-queue-remove b)))
-        (declare (type pair pair))
-        (cond
-          ((criterion-1 pair) nil)
-          ((criterion-2 pair b-done f) nil)
-          (t 
-           (let ((sp (normal-form ring (spoly ring (pair-first pair)
-                                              (pair-second pair))
-                                  f top-reduction-only)))
-             (declare (type poly sp))
-             (cond
-               ((poly-zerop sp)
-                nil)
-               (t
-                (setf sp (poly-primitive-part ring sp)
-                      f (nconc f (list sp)))
-                ;; Add new critical pairs
-                (dolist (h f)
-                  (pair-queue-insert b (make-pair h sp)))
-                (debug-cgb "~&Sugar: ~d Polynomials: ~d; Pairs left: ~d; Pairs done: ~d;"
-                           (pair-sugar pair) (length f) (pair-queue-size b)
-                           (hash-table-count b-done)))))))
-        (setf (gethash (list (pair-first pair) (pair-second pair)) b-done)
-              t)))))
-
-(defun parallel-buchberger (ring f start &optional (top-reduction-only $poly_top_reduction_only))
-  "An implementation of the Buchberger algorithm. Return Grobner basis
-of the ideal generated by the polynomial list F.  Polynomials 0 to
-START-1 are assumed to be a Grobner basis already, so that certain
-critical pairs will not be examined. If TOP-REDUCTION-ONLY set, top
-reduction will be preformed."
-  (declare (ignore top-reduction-only)
-           (type fixnum start))
-  (when (endp f) (return-from parallel-buchberger f)) ;cut startup costs
-  (debug-cgb "~&GROBNER BASIS - PARALLEL-BUCHBERGER ALGORITHM")
-  (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
-  #+grobner-check  (when (plusp start)
-                     (grobner-test ring (subseq f 0 start) (subseq f 0 start)))
-  ;;Initialize critical pairs
-  (let ((b (pair-queue-initialize (make-pair-queue) f start))
-        (b-done (make-hash-table :test #'equal)))
-    (declare (type priority-queue b)
-             (type hash-table b-done))
-    (dotimes (i (1- start))
-      (do ((j (1+ i) (1+ j))) ((>= j start))
-        (declare (type fixnum j))
-        (setf (gethash (list (elt f i) (elt f j)) b-done) t)))
-    (do ()
-        ((pair-queue-empty-p b)
-         #+grobner-check(grobner-test ring f f)
-         (debug-cgb "~&GROBNER END")
-         f)
-      (let ((pair (pair-queue-remove b)))
-        (when (null (pair-division-data pair))
-          (setf (pair-division-data pair) (list (spoly ring
-                                                       (pair-first pair)
-                                                       (pair-second pair))
-                                                (make-poly-zero)
-                                                (funcall (ring-unit ring))
-                                                0)))
-        (cond
-          ((criterion-1 pair) nil)
-          ((criterion-2 pair b-done f) nil)
-          (t
-           (let* ((dd (pair-division-data pair))
-                  (p (first dd))
-                  (sp (second dd))
-                  (c (third dd))
-                  (division-count (fourth dd)))
-             (cond
-               ((poly-zerop p)          ;normal form completed
-                (debug-cgb "~&~3T~d reduction~:p" division-count)
-                (cond 
-                  ((poly-zerop sp)
-                   (debug-cgb " ---> 0")
-                   nil)
-                  (t
-                   (setf sp (poly-nreverse sp)
-                         sp (poly-primitive-part ring sp)
-                         f (nconc f (list sp)))
-                   ;; Add new critical pairs
-                   (dolist (h f)
-                     (pair-queue-insert b (make-pair h sp)))
-                   (debug-cgb "~&Sugar: ~d Polynomials: ~d; Pairs left: ~d; Pairs done: ~d;"
-                              (pair-sugar pair) (length f) (pair-queue-size b)
-                              (hash-table-count b-done))))
-                (setf (gethash (list (pair-first pair) (pair-second pair))
-                               b-done) t))
-               (t                               ;normal form not complete
-                (do ()
-                    ((cond
-                       ((> (poly-sugar sp) (pair-sugar pair))
-                        (debug-cgb "(~a)?" (poly-sugar sp))
-                        t)
-                       ((poly-zerop p)
-                        (debug-cgb ".")
-                        t)
-                       (t nil))
-                     (setf (first dd) p
-                           (second dd) sp
-                           (third dd) c
-                           (fourth dd) division-count
-                           (pair-sugar pair) (poly-sugar sp))
-                     (pair-queue-insert b pair))
-                  (multiple-value-setq (p sp c division-count)
-                    (normal-form-step ring f p sp c division-count))))))))))))
-
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;
-;; Grobner Criteria
-;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-
-(defun criterion-1 (pair)
-  "Returns T if the leading monomials of the two polynomials
-in G pointed to by the integers in PAIR have disjoint (relatively prime)
-monomials. This test is known as the first Buchberger criterion."
-  (declare (type pair pair))
-  (let ((f (pair-first pair))
-        (g (pair-second pair)))
-    (when (monom-rel-prime-p (poly-lm f) (poly-lm g))
-      (debug-cgb ":1")
-      (return-from criterion-1 t))))
-
-(defun criterion-2 (pair b-done partial-basis
-                    &aux (f (pair-first pair)) (g (pair-second pair))
-                         (place :before))
-  "Returns T if the leading monomial of some element P of
-PARTIAL-BASIS divides the LCM of the leading monomials of the two
-polynomials in the polynomial list PARTIAL-BASIS, and P paired with
-each of the polynomials pointed to by the the PAIR has already been
-treated, as indicated by the absence in the hash table B-done."
-  (declare (type pair pair) (type hash-table b-done)
-           (type poly f g))
-  ;; In the code below we assume that pairs are ordered as follows: 
-  ;; if PAIR is (I J) then I appears before J in the PARTIAL-BASIS.
-  ;; We traverse the list PARTIAL-BASIS and keep track of where we
-  ;; are, so that we can produce the pairs in the correct order
-  ;; when we check whether they have been processed, i.e they
-  ;; appear in the hash table B-done
-  (dolist (h partial-basis nil)
-    (cond
-     ((eq h f)
-      #+grobner-check(assert (eq place :before))
-      (setf place :in-the-middle))
-     ((eq h g)
-      #+grobner-check(assert (eq place :in-the-middle))
-      (setf place :after))
-     ((and (monom-divides-monom-lcm-p (poly-lm h) (poly-lm f) (poly-lm g))
-           (gethash (case place
-                      (:before (list h f))
-                      ((:in-the-middle :after) (list f h)))
-                    b-done)
-           (gethash (case place
-                      ((:before :in-the-middle) (list h g))
-                      (:after (list g h)))
-                    b-done))
-      (debug-cgb ":2")
-      (return-from criterion-2 t)))))
-
-
-
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;
-;; An implementation of the algorithm of Gebauer and Moeller, as
-;; described in the book of Becker-Weispfenning, p. 232
-;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-
-(defun gebauer-moeller (ring f start &optional (top-reduction-only $poly_top_reduction_only))
-  "Compute Grobner basis by using the algorithm of Gebauer and
-Moeller.  This algorithm is described as BUCHBERGERNEW2 in the book by
-Becker-Weispfenning entitled ``Grobner Bases''. This function assumes
-that all polynomials in F are non-zero."
-  (declare (ignore top-reduction-only)
-           (type fixnum start))
-  (cond
-   ((endp f) (return-from gebauer-moeller nil))
-   ((endp (cdr f))
-    (return-from gebauer-moeller (list (poly-primitive-part ring (car f))))))
-   (debug-cgb "~&GROBNER BASIS - GEBAUER MOELLER ALGORITHM")
-   (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
-  #+grobner-check  (when (plusp start)
-                     (grobner-test ring (subseq f 0 start) (subseq f 0 start)))
-  (let ((b (make-pair-queue))
-        (g (subseq f 0 start))
-        (f1 (subseq f start)))
-    (do () ((endp f1))
-      (multiple-value-setq (g b)
-        (gebauer-moeller-update g b (poly-primitive-part ring (pop f1)))))
-    (do () ((pair-queue-empty-p b))
-      (let* ((pair (pair-queue-remove b))
-             (g1 (pair-first pair))
-             (g2 (pair-second pair))
-             (h (normal-form ring (spoly ring g1 g2)
-                             g
-                             nil #| Always fully reduce! |#
-                             )))
-        (unless (poly-zerop h)
-          (setf h (poly-primitive-part ring h))
-          (multiple-value-setq (g b)
-            (gebauer-moeller-update g b h))
-          (debug-cgb "~&Sugar: ~d Polynomials: ~d; Pairs left: ~d~%"
-                     (pair-sugar pair) (length g) (pair-queue-size b))
-          )))
-    #+grobner-check(grobner-test ring g f)
-    (debug-cgb "~&GROBNER END")
-    g))
-
-(defun gebauer-moeller-update (g b h
-                 &aux
-                 c d e
-                 (b-new (make-pair-queue))
-                 g-new)
-  "An implementation of the auxillary UPDATE algorithm used by the
-Gebauer-Moeller algorithm. G is a list of polynomials, B is a list of
-critical pairs and H is a new polynomial which possibly will be added
-to G. The naming conventions used are very close to the one used in
-the book of Becker-Weispfenning."
-  (declare
-   #+allegro (dynamic-extent b)
-   (type poly h)
-   (type priority-queue b))
-  (setf c g d nil) 
-  (do () ((endp c))
-    (let ((g1 (pop c)))
-      (declare (type poly g1))
-      (when (or (monom-rel-prime-p (poly-lm h) (poly-lm g1))
-                (and
-                 (notany #'(lambda (g2) (monom-lcm-divides-monom-lcm-p
-                                         (poly-lm h) (poly-lm g2)
-                                         (poly-lm h) (poly-lm g1)))
-                         c)
-                 (notany #'(lambda (g2) (monom-lcm-divides-monom-lcm-p
-                                         (poly-lm h) (poly-lm g2)
-                                         (poly-lm h) (poly-lm g1)))
-                         d)))
-        (push g1 d))))
-  (setf e nil)
-  (do () ((endp d))
-    (let ((g1 (pop d)))
-      (declare (type poly g1))
-      (unless (monom-rel-prime-p (poly-lm h) (poly-lm g1))
-        (push g1 e))))
-  (do () ((pair-queue-empty-p b))
-    (let* ((pair (pair-queue-remove b))
-           (g1 (pair-first pair))
-           (g2 (pair-second pair)))
-      (declare (type pair pair)
-               (type poly g1 g2))
-      (when (or (not (monom-divides-monom-lcm-p
-                      (poly-lm h)
-                      (poly-lm g1) (poly-lm g2)))
-                (monom-lcm-equal-monom-lcm-p
-                 (poly-lm g1) (poly-lm h)
-                 (poly-lm g1) (poly-lm g2))
-                (monom-lcm-equal-monom-lcm-p
-                 (poly-lm h) (poly-lm g2)
-                 (poly-lm g1) (poly-lm g2)))
-        (pair-queue-insert b-new (make-pair g1 g2)))))
-  (dolist (g3 e)
-    (pair-queue-insert b-new (make-pair h g3)))
-  (setf g-new nil)
-  (do () ((endp g))
-    (let ((g1 (pop g)))
-      (declare (type poly g1))
-      (unless (monom-divides-p (poly-lm h) (poly-lm g1))
-        (push g1 g-new))))
-  (push h g-new)
-  (values g-new b-new))
-
-
-
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;
-;; Standard postprocessing of Grobner bases
-;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-
-(defun reduction (ring plist)
-  "Reduce a list of polynomials PLIST, so that non of the terms in any of
-the polynomials is divisible by a leading monomial of another
-polynomial.  Return the reduced list."
-  (do ((q plist)
-       (found t))
-      ((not found)
-       (mapcar #'(lambda (x) (poly-primitive-part ring x)) q))
-    ;;Find p in Q such that p is reducible mod Q\{p}
-    (setf found nil)
-    (dolist (x q)
-      (let ((q1 (remove x q)))
-        (multiple-value-bind (h c div-count)
-            (normal-form ring x q1 nil #| not a top reduction! |# )
-          (declare (ignore c))
-          (unless (zerop div-count)
-            (setf found t q q1)
-            (unless (poly-zerop h)
-              (setf q (nconc q1 (list h))))
-            (return)))))))
-
-(defun minimization (p)
-  "Returns a sublist of the polynomial list P spanning the same
-monomial ideal as P but minimal, i.e. no leading monomial
-of a polynomial in the sublist divides the leading monomial
-of another polynomial."
-  (do ((q p)
-       (found t))
-      ((not found) q)
-    ;;Find p in Q such that lm(p) is in LM(Q\{p})
-    (setf found nil
-          q (dolist (x q q)
-              (let ((q1 (remove x q)))
-                (when (member-if #'(lambda (p) (monom-divides-p (poly-lm x) (poly-lm p))) q1)
-                  (setf found t)
-                  (return q1)))))))
-
-(defun poly-normalize (ring p &aux (c (poly-lc p)))
-  "Divide a polynomial by its leading coefficient. It assumes
-that the division is possible, which may not always be the
-case in rings which are not fields. The exact division operator
-is assumed to be provided by the RING structure of the
-COEFFICIENT-RING package."
-  (mapc #'(lambda (term)
-            (setf (term-coeff term) (funcall (ring-div ring) (term-coeff term) c)))
-        (poly-termlist p))
-  p)
-
-(defun poly-normalize-list (ring plist)
-  "Divide every polynomial in a list PLIST by its leading coefficient. "
-  (mapcar #'(lambda (x) (poly-normalize ring x)) plist))
-
-
-
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-;;
-;; Algorithm and Pair heuristic selection
-;;
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
 (defun find-grobner-function (algorithm)