head	1.8;
access;
symbols;
locks; strict;
comment	@;;; @;


1.8
date	2009.01.23.10.40.26;	author marek;	state Exp;
branches;
next	1.7;

1.7
date	2009.01.23.10.28.38;	author marek;	state Exp;
branches;
next	1.6;

1.6
date	2009.01.23.10.28.05;	author marek;	state Exp;
branches;
next	1.5;

1.5
date	2009.01.22.04.02.22;	author marek;	state Exp;
branches;
next	1.4;

1.4
date	2009.01.19.09.25.46;	author marek;	state Exp;
branches;
next	1.3;

1.3
date	2009.01.19.07.50.10;	author marek;	state Exp;
branches;
next	1.2;

1.2
date	2009.01.19.06.49.18;	author marek;	state Exp;
branches;
next	1.1;

1.1
date	2009.01.19.06.48.15;	author marek;	state Exp;
branches;
next	;


desc
@@


1.8
log
@*** empty log message ***
@
text
@;;; -*- Mode: Lisp; Syntax: Common-Lisp; Package: Grobner; Base: 10 -*-
#|
	$Id$
  *--------------------------------------------------------------------------*
  |  Copyright (C) 1994, Marek Rychlik (e-mail: rychlik@@math.arizona.edu)    |
  |    Department of Mathematics, University of Arizona, Tucson, AZ 85721    |
  |                                                                          |
  | Everyone is permitted to copy, distribute and modify the code in this    |
  | directory, as long as this copyright note is preserved verbatim.         |
  *--------------------------------------------------------------------------*
|#
(defpackage "GROBNER"
  (:use "ORDER" "PARSE" "PRINTER" "DIVISION" "MONOM"
	"MAKELIST" "COEFFICIENT-RING" "TERM" "POLY"
	"POLY-WITH-SUGAR" "COMMON-LISP")
  (:export
   spoly
   grobner
   reduction
   normal-form
   reduced-grobner
   normalize-poly
   normalize-basis
   elimination-ideal
   ring-intersection
   poly-contract
   ideal-intersection
   poly-lcm
   grobner-gcd
   pseudo-divide
   buchberger-criterion
   grobner-test
   grobner-equal
   grobner-subsetp
   grobner-member
   ideal-equal
   ideal-subsetp
   ideal-member
   ideal-saturation
   ideal-polysaturation
   ideal-saturation-1
   ideal-polysaturation-1
   grobner-primitive-part
   grobner-content
   colon-ideal
   *grobner-debug*
   buchberger-set-pair-heuristic
   gebauer-moeller-set-pair-heuristic
   select-grobner-algorithm
   spoly-with-sugar
   normal-form-with-sugar
   grobner-with-sugar
   ))

(in-package "GROBNER")

#+debug(proclaim '(optimize (speed 0) (debug 3)))
#-debug(proclaim '(optimize (speed 3) (debug 0)))

#+debug
(defvar *grobner-debug* nil
  "If T, produce debugging and tracing output.")

(defvar *buchberger-merge-pairs* 'buchberger-merge-pairs-smallest-lcm
  "A variable holding the predicate used to sort critical pairs
in the order of decreasing priority for the Buchberger algorithm.")

(defvar *gebauer-moeller-merge-pairs* 'gebauer-moeller-merge-pairs-use-mock-spoly
  "A variable holding the predicate used to sort critical pairs
in the order of decreasing priority for the Gebauer-Moeller algorithm")

(defvar *grobner-function* 'buchberger
  "The name of the function used to calculate Grobner basis")

(defun select-grobner-algorithm (algorithm)
  "Select one of the several implementations of the Grobner basis
algorithm."
  (setf *grobner-function*
    (ecase algorithm
      (:buchberger 'buchberger)
      (:cached-buchberger 'cached-buchberger)
      (:gebauer-moeller 'gebauer-moeller)
      (:sugar 'buchberger-with-sugar)
      (:cached-sugar 'cached-buchberger-with-sugar)
      (:gebauer-moeller-with-sugar 'gebauer-moeller-with-sugar))))

(defun grobner (F &optional (pred #'lex>)
			    (start 0)
			    (top-reduction-only nil)
			    (ring *coefficient-ring*))
  "Return Grobner basis of an ideal generated by polynomial list F.
Assume that the monomials of each polynomial are ordered according to
the admissible monomial order PRED. The result will be a list of
polynomials ordered as well according to PRED. The polynomials from 0
to START-1 in the list F are assumed to be a Grobner basis, and thus
certain critical pairs do not have to be calculated. If
TOP-REDUCTION-ONLY is not NIL then only top reductions will be
performed, instead of full division, and thus the returned Grobner
basis will have its polynomials only top-reduced with respect to the
preceding polynomials. RING should be set to the RING structure (see
COEFFICIENT-RING package) corresponding to the coefficients of the
polynomials in the list F. This function is currently just an
interface to several algorithms which fine Grobner bases. Use
SELECT-GROBNER-ALGORITHM in order to change the algorithm."
  (funcall *grobner-function* F pred start top-reduction-only ring))


#+debug
(defmacro debug-cgb (&rest args)
  "This macro is used in printing debugging and tracing messages
and its arguments are analogous to the FORMAT function, except
that the stream argument is omitted. By default, the output is
sent to *trace-output*, which can be set to produce output 
in a file."
  `(when *grobner-debug* (format *trace-output* ,@@args)))

(defun spoly (f g pred ring)
  "Return the S-polynomial of F and G, which should
be two polynomials with terms ordered by PRED and
coefficients in a ring whose operations are the slots
of the RING structure."
  (declare (optimize (speed 3) (safety 0)))
  (let* ((lcm (monom-lcm (lm f) (lm g)))
	 (m1 (monom/ lcm (lm f)))
	 (m2 (monom/ lcm (lm g))))
    (let* ((c (funcall (ring-lcm ring) (lc f) (lc g)))
	   (cf (funcall (ring-/ ring) c (lc f)))
	   (cg (funcall (ring-/ ring) c (lc g))))
      (poly-
       (term-times-poly (cons m1 cf) (rest f) ring)
       (term-times-poly (cons m2 cg) (rest g) ring)
       pred
       ring))))


(defun grobner-primitive-part (p ring)
  "Divide polynomial P with integer coefficients by gcd of its
coefficients and return the result. Use the RING structure from the
COEFFICIENT-RING package to perform the operations on the
coefficients."
  (if (poly-zerop p)
      (values p 1)
    (let ((c (grobner-content p ring)))
      (values (mapcar
	       #'(lambda (x)
		   (cons (car x)
			 (funcall
			  (ring-/ ring)
			  (cdr x) c)))
	       p)
	      c))))


(defun grobner-content (p ring)
  "Greatest common divisor of the coefficients of the polynomial P. Use the RING structure
to compute the greatest common divisor."
  (let ((c (apply (ring-gcd ring) (mapcar #'cdr p))))
    (when (/= (funcall (ring-signum ring) (cdar p))
	      (funcall (ring-signum ring) c))
      (setf c (funcall (ring-- ring) c)))
    c))

(defun normal-form (f fl pred  top-reduction-only ring)
  "Calculate the normal form of the polynomial F with respect to
the polynomial list FL. Destructive to F; thus, copy-tree should be used
where F needs to be preserved. It returns multiple values.
The first value is the polynomial R which is reduced (or top-reduced
if TOP-REDUCTION-ONLY is not NIL). The second value is an integer
which is the number of reductions actually performed. The third value
is the coefficient by which F needs to be multiplied in order to be
able to perform the division without actually having to divide in the
coefficient ring (a form of pseudo-division is used). All operations
on the coefficients are performed using the RING structure from the
COEFFICIENT-RING package."
  (declare (optimize (speed 3) (safety 0)))
  ;; Loop invariant: c*f=sum ai*fi+r+p
  (do (r (c (ring-unit ring))
       (division-count 0)
       (p f))
      ((or (endp p)
	   (endp fl)
	   (and top-reduction-only r))
       #+debug(progn
		(debug-cgb "~&~3T~d reduction~:p" division-count)
		(when (endp r)
		  (debug-cgb " ---> 0")))
       (values (append (nreverse r) p)
	       division-count c))
    (declare (fixnum division-count) (integer c))
    (let ((g (find (lm p) fl :test #'monom-divisible-by-p :key #'lm)))
      (cond 
       (g				;division possible
	(incf division-count)
	(let* ((gcd (funcall (ring-gcd ring) (lc g) (lc p)))
	       (c1 (funcall (ring-/ ring) (lc g) gcd))
	       (c2 (funcall (ring-/ ring) (lc p) gcd))
	       (m (monom/ (lm p) (lm g))))
	  ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
	  (setf r (scalar-times-poly c1 r ring)
		c (funcall (ring-* ring) c c1)
		p (grobner-op c2 c1 m (rest p) (rest g) pred ring))))
       (t				;no division possible
	(setf r (cons (lt p) r)		;move lt(p) to remainder
	      p (rest p)))))))		;remove lt(p) from p	

(defun buchberger (F pred start top-reduction-only ring 
		   &aux (s (1- (length F))) B M)
  "An implementation of the Buchberger algorithm. Return Grobner
basis of the ideal generated by the polynomial list F.
The terms of each polynomial are ordered by the admissible
monomial order PRED. 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. The RING structure is used to perform operations
on coefficents (see COEFFICIENT-RING package)."
  (declare (fixnum s))
  (unless (plusp s) (return-from buchberger F)) ;cut startup costs
  #+debug (debug-cgb "~&GROBNER BASIS - BUCHBERGER ALGORITHM")
  #+debug (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
  #+grobner-check  (when (plusp start)
		     (grobner-test (subseq F 0 start) (subseq F 0 start) pred ring))
  (setf B (nconc (makelist (list i j) (i 0 (1- start)) (j start s))
		 (makelist (list i j) (i start (1- s)) (j (1+ i) s)))
	M (make-hash-table :test #'equal)
	B (buchberger-sort-pairs B F pred ring))
  ;;Initialize treated pairs
  (dotimes (i (1- start))
    (do ((j (1+ i) (1+ j))) ((>= j start))
      (setf (gethash (list i j) M) t)))
  (do ((G F))
      ((endp B)
       #+grobner-check(grobner-test G F pred ring)
       #+debug(debug-cgb "~&GROBNER END")
       G)
    (let ((pair (pop B)))
      (unless (or (Criterion-1 pair G)
		  (Criterion-2 pair G M))
	(let ((SP (normal-form (spoly (elt G (first pair))
				      (elt G (second pair))
				      pred
				      ring)
			       G pred top-reduction-only ring)))
	  (unless (poly-zerop SP)	;S-polynomial ---> 0;
	    (setf s (1+ s)
		  SP (grobner-primitive-part SP ring)
		  G (nconc G (list SP))
		  B (funcall *buchberger-merge-pairs* B (makelist (list i s)
							       (i 0 (1- s)))
			     G pred ring))
	    #+debug (debug-cgb "~&Polynomials: ~d; Pairs left: ~d; Treated pairs: ~d"
			       (length G) (length B)
			       (hash-table-count M)))))
      (setf (gethash pair M) t))))

(defun grobner-op (c1 c2 m A B pred ring &aux n A1 A2)
  "Perform a calculation equivalent to
(POLY- (SCALAR-TIMES-POLY C2 A)
	 (SCALAR-TIMES-POLY C1 (MONOM-TIMES-POLY M B))
	 PRED)
but more efficiently; it destructively modifies A and stores the 
result in it; A is returned."
  (dolist (x A) (setf (cdr x) (funcall (ring-* ring) c2 (cdr x))))
  (unless B
    (return-from grobner-op A))
  (unless A
    (return-from grobner-op
      (mapcar 
       #'(lambda (x) (cons (monom* m (car x))
			   (funcall (ring-- ring)
				    (funcall (ring-* ring)
					     c1 (cdr x)))))
       B)))
  (setf A1 (cons nil A)	A2 A1)
  (tagbody
   begin-outer-loop
    (unless B
      (return-from grobner-op (cdr A2)))
    (setf n (monom* m (caar B)))
    (tagbody
     begin-inner-loop
      (unless (cdr A1)
	(setf (cdr A1) 
	  (mapcar #'(lambda (x) (cons (monom* m (car x))
				      (funcall (ring-- ring)
					       (funcall (ring-* ring)
							c1 (cdr x))))) B))
	(return-from grobner-op (cdr A2)))
      (multiple-value-bind (lm-greater lm-equal)
	  (funcall pred (caadr A1) n)
	(cond
	 (lm-equal
	  (let ((s (funcall (ring-- ring)
			    (cdadr A1)
			    (funcall (ring-* ring) c1 (cdar B)))))
	    (if (funcall (ring-zerop ring) s) ;check for cancellation
		(setf (cdr A1) (cddr A1))
	      (setf (cdadr A1) s))
	    (setf B (cdr B))
	    (go begin-outer-loop)))
	 (lm-greater
	  (setf A1 (cdr A1))
	  (go begin-inner-loop))
	 (t				;(lm B) is greater than (lm A)
	  (setf (cdr A1)
		(cons
		 (cons
		  n
		  (funcall (ring-- ring)
			   (funcall (ring-* ring) c1 (cdar B))))
		 (cdr A1))
		B (cdr B)
		A1 (cdr A1))
	  (go begin-outer-loop)))))))

;;----------------------------------------------------------------
;; Pairs
;;----------------------------------------------------------------

(defun buchberger-sort-pairs (B G pred ring)
  "Sort critical pairs B by the function which is
the value of the vairable *buchberger-merge-pairs*.
G is the list of polynomials which the elemnts
of B point into. PRED is the admissible monomial order.
The RING structure holds operations on the coefficients
as described in the COEFFICIENT-RING package."
  (funcall *buchberger-merge-pairs* nil B G pred ring))


(defun mock-spoly (f g pred)
  "Returns an upper-bound on the leading
monomial of the S-polynomial of F and G, which are two polynomials
sorted by the admissible monomial order PRED."
  (let* ((lcm (monom-lcm (lm f) (lm g)))
	 (m1 (monom/ lcm (lm f)))
	 (m2 (monom/ lcm (lm g))))
    (cond 
     ((endp (second f))
      (monom* m2 (car (second g))))
     ((endp (second g))
      (monom* m1 (car (second f))))
     (t (let ((n1 (monom* m2 (car (second g))))
	      (n2 (monom* m1 (car (second f)))))
	  (if (funcall pred n1 n2)
	      n1
	    n2))))))

;;----------------------------------------------------------------
;; A stragegy based on a predicted approximate S-polynomials
;;----------------------------------------------------------------
(defun buchberger-merge-pairs-use-mock-spoly (B C G pred ring)
  "Merges lists of critical pairs B and C into one list. 
The pairs are assumed to be ordered according to the
increasing value of MOCK-SPOLY, as determined by the admissible
monomial order PRED. G is the list of polynomials which the pairs
of integers in B and C point into, and RING is the ring structure
defined in the COEFFICIENT-RING package."
  (declare (ignore ring))
  ;; Delete pairs consisting of two single monomials
  (setf C (delete-if
	   #'(lambda (p) (and (endp (cdr (elt G (first p))))
			      (endp (cdr (elt G (second p))))))
	   C))
  (labels ((pair-key (p) (mock-spoly (elt G (first p))
				     (elt G (second p)) pred))
	   (pair-order (x y) (funcall pred y x)))
    (merge 'list B (sort C #'pair-order :key #'pair-key)
	   #'pair-order :key #'pair-key)))


;;----------------------------------------------------------------
;; A strategy based on the smallest lcm of the leading
;; monomials of the two polynomials - so called normal strategy
;;----------------------------------------------------------------
(defun buchberger-merge-pairs-smallest-lcm (B C G pred ring)
  "Merges lists of critical pairs B and C into a single ordered
list. Implements a strategy of ordering pairs according to the
smallest lcm of the leading monomials of the two polynomials - so
called normal strategy."
  (declare (ignore ring))
  (labels ((pair-key (p) (monom-lcm (lm (elt G (first p)))
				    (lm (elt G (second p)))))
	   (pair-order (x y) (funcall pred y x)))
    (merge 'list B (sort C #'pair-order :key #'pair-key)
	   #'pair-order :key #'pair-key)))

;;----------------------------------------------------------------
;; A strategy based on the smallest degree of the lcm of the leading
;; monomials of the two polynomials
;;----------------------------------------------------------------
(defun buchberger-merge-pairs-use-smallest-degree (B C G pred ring)
  "Mergest lists B and C of critical pairs. This strategy is based on
 ordering the pairs according to the smallest degree of the lcm of the
 leading monomials of the two polynomials."
  (declare (ignore pred ring))
  (labels ((pair-key (p) (total-degree
			  (monom-lcm (lm (elt G (first p)))
				     (lm (elt G (second p)))))))
    (merge 'list B (sort C #'< :key #'pair-key) #'< :key #'pair-key)))

(defun buchberger-merge-pairs-use-smallest-length (B C G pred ring)
  "Mergest lists B and C of critical pairs. This strategy is based on
 ordering the pairs according to the smallest total length of the 
 two polynomials, where length is the number of terms."
  (declare (ignore pred ring))
  (labels ((pair-key (p) (+
			  (length (elt G (first p)))
			  (length (elt G (second p))))))
    (merge 'list B (sort C #'< :key #'pair-key) #'< :key #'pair-key)))

(defun buchberger-merge-pairs-use-smallest-coefficient-length (B C G pred ring)
  "Mergest lists B and C of critical pairs. This strategy is based on
 ordering the pairs according to the smallest combined length of the 
 coefficients of the two polynomials."
  (declare (ignore pred))
  (labels ((coefficient-length (poly)
	     (apply #'+ (mapcar (ring-length ring)
				(mapcar #'cdr poly))))
	   (pair-key (p) (+
			  (coefficient-length (elt G (first p)))
			  (coefficient-length (elt G (second p))))))
    (merge 'list B (sort C #'< :key #'pair-key) #'< :key #'pair-key)))

(defun buchberger-set-pair-heuristic
    (method
     &aux
     (strategy-fn
      (ecase method
	(:minimal-mock-spoly #'buchberger-merge-pairs-use-mock-spoly)
	(:minimal-lcm #'buchberger-merge-pairs-smallest-lcm)
	(:minimal-total-degree #'buchberger-merge-pairs-use-smallest-degree)
	(:minimal-length #'buchberger-merge-pairs-use-smallest-length)
	(:minimal-coefficient-length #'buchberger-merge-pairs-use-smallest-coefficient-length))))
  "Simply sets the variable *buchberger-merge-pairs* to the heuristic
function STRATEGY-FN implementing one of several strategies introduces
before of selecting the most promising. METHOD should be one of the
listed keywords."
  (setf *buchberger-merge-pairs* strategy-fn))

;;----------------------------------------------------------------
;; Grobner Criteria
;;----------------------------------------------------------------
(defun Criterion-1 (pair G &aux (i (first pair)) (j (second 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."
  (monom-rel-prime (lm (elt G i))
		   (lm (elt G j))))


(defun Criterion-2 (pair G M &aux (i (first pair)) (j (second pair)))
  "Returns T if the leading monomial of some element P of G divides
the LCM of the leading monomials of the two polynomials in the
polynomial list G, pointed to by the two integers in PAIR, and P
paired with each of the polynomials pointed to by the the PAIR has
already been treated, as indicated by the hash table M, which stores
value T for every treated pair so far."
  (labels ((make-pair (u v) (if (< u v) (list u v) (list v u))))
    (dotimes (k (length G) nil)
      (when (and (/= k i)
		 (/= k j)
		 (monom-divides-p
		  (lm (elt G k))
		  (monom-lcm (lm (elt G i)) (lm (elt G j))))
		 (gethash (make-pair i k) M)
		 (gethash (make-pair k j) M))
	#+debug(when *grobner-debug* (format t ":B2"))
	(return-from Criterion-2 t)))))

;;----------------------------------------------------------------
;; End of GROBNER implementation
;;----------------------------------------------------------------

(defun normalize-poly (p ring)
  "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."
  (mapcar #'(lambda (term)
	      (cons
	       (car term)
	       (funcall (ring-/ ring) (cdr term) (lc p)))) p))

(defun normalize-basis (plist ring)
  "Divide every polynomial in a list PLIST by its leading coefficient.
Use RING structure to perform the division of the coefficients."
  (mapcar #'(lambda (x) (normalize-poly x ring)) plist))

(defun reduction (P pred ring)
  "Reduce a list of polynomials P, 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 P)
       (found t))
      ((not found)
       (mapcar #'(lambda (x) (grobner-primitive-part x ring)) 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 div-count)
	    (normal-form x Q1
			 pred
			 nil		;not a top reduction
			 ring)
	  (unless (zerop div-count)
	    (setf found t Q Q1)
	    (when h
	      (setf Q (nconc Q1 (list  h))))
	    (return)))))))


(defun reduced-grobner (F &optional (pred #'lex>)
				    (start 0)
				    (top-reduction-only nil)
				    (ring *coefficient-ring*))
  "Return the reduced Grobner basis of the ideal generated by a polynomial
list F. This combines calls to two functions: GROBNER and REDUCTION.
The parameters have the same meaning as in GROBNER or BUCHBERGER."
  (reduction (grobner F pred start top-reduction-only ring) pred ring))

#+kcl
(defvar *vars* nil "Keep track of variables for tracing SPOLY.")

;; Does the monomial P depend on variable K
(defun monom-depends-p (m k)
  "Return T if the monomial M depends on variable number K."
  (plusp (elt m k)))

;; Does the term depend on variable K
(defun term-depends-p (term k)
  "Return T if the term TERM depends on variable number K."
  (monom-depends-p (car term) k))

;; Does the polynomial P depend on variable K
(defun poly-depends-p (p k)
  "Return T if the term polynomial P depends on variable number K."
  (some #'(lambda (term) (term-depends-p term k)) p))

(defun ring-intersection (plist k &key (key #'identity))
  "This function assumes that polynomial list PLIST is a Grobner basis
and it calculates the intersection with the ring R[x[k+1],...,x[n]], i.e.
it discards polynomials which depend on variables x[0], x[1], ..., x[k]."
  (dotimes (i k plist)
    (setf plist
	  (remove-if #'(lambda (p)
			 (poly-depends-p p (funcall key i)))
		     plist))))

(defun elimination-ideal (flist k
				&key
				(primary-order #'lex>)
				(secondary-order #'lex>)
				(start 0)
				(order (elimination-order
					k
					:primary-order primary-order
					:secondary-order secondary-order
					))
				(top-reduction-only nil)
				(ring *coefficient-ring*))
  "Returns the K-th elimination ideal of the ideal generated by polynomial list FLIST.
Thus, a Grobner basis of the ideal generated by FLIST is calculated and
polynomials depending on variables 0-K are discarded. The monomial order in the
calculation is an elimination order formed by combining two orders: PRIMARY-ORDER
used on variables 0-K and SECONDARY-ORDER used on variables starting from variable K+1.
Thus, if both orders are set to the lexicographic order #'LEX> then the resulting
order is simply equivalent to #'LEX>. But one could use arbitrary two orders in this
function, for instance, two copies of #'GREVLEX> (see ORDER package). When doing so,
the caller must ensure that the same order has been used to sort the terms of
the polynomials FLIST."
  (ring-intersection
   (reduced-grobner flist order start top-reduction-only ring)
   k))


(defun ideal-intersection (F G pred top-reduction-only ring)
  "Return the Grobner basis of the intersection of the ideals
generated by the polynomial lists F and G. PRED is an admissible
monomial order with respect to which the terms of F and G have been
sorted. This order is going to be used in the Grobner basis
calculation. If TOP-REDUCTION-ONLY is not NIL, internally the Grobner
basis algorithm will perform top reduction only.  The RING parameter,
as usual, should be set to a structure defined in the package
COEFFICIENT-RING, and it will be used to perform all operations on
coefficients."
  (mapcar
   #'poly-contract
   (ring-intersection
    (reduced-grobner
     (append
      (mapcar #'(lambda (p) (poly-extend p (list 1))) F)
      (mapcar #'(lambda (p)
		  (append (poly-extend (minus-poly p ring) (list 1))
			  (poly-extend p)))
	      G))
     (elimination-order-1 pred)
     0
     top-reduction-only
     ring)
    1)))


(defun poly-contract (f &optional (k 1))
  "Return a polynomial obtained from polynomial F by dropping the
first K variables, if F does not depend on them.  Note that this
operation makes the polynomial incompatible for certain arithmetical
operations with the original polynomial F. Calling this function on a polynomial F
which does depend on variables 0-K will result in error."
  (when (some #'(lambda (term) (find-if-not #'zerop (car term) :end k)) f)
    (error "Contracting a polynomial which depends on contracted variables"))
  (mapcar #'(lambda (term) (cons (nthcdr k (car term)) (cdr term))) f))


(defun poly-lcm (f g &optional (pred #'lex>) (ring *coefficient-ring*))
  "Return LCM (least common multiple) of two polynomials F and G.
The polynomials must be ordered according to monomial order PRED
and their coefficients must be compatible with the RING structure
defined in the COEFFICIENT-RING package."
  (cond
    ((endp f) nil)
    ((endp g) nil)
    ((and (endp (cdr f)) (endp (cdr g)))
     (list (cons (monom-lcm (caar f) (caar g))
		 (lcm (cdar f) (cdar g)))))
    (t
     (multiple-value-bind (f f-cont)
	 (grobner-primitive-part f ring)
       (multiple-value-bind (g g-cont)
	   (grobner-primitive-part g ring)
	 (scalar-times-poly
	  (funcall (ring-lcm ring) f-cont g-cont)
	  (grobner-primitive-part (car (ideal-intersection (list f) (list g) pred nil ring))
				  ring)
	  ring))))))

;; GCD of two polynomials, using Grobner bases
(defun grobner-gcd (f g &optional (pred #'lex>) (ring *coefficient-ring*))
  "Return GCD (greatest common divisor) of two polynomials F and G.
The polynomials must be ordered according to monomial order PRED
and their coefficients must be compatible with the RING structure
defined in the COEFFICIENT-RING package."
  (poly-exact-divide (poly* f g pred ring)
		     (poly-lcm f g pred ring)
		     pred ring))

;; Do two Grobner bases yield the same ideal?
(defun grobner-equal (g1 g2 &optional (pred #'lex>) (ring *coefficient-ring*))
  "Returns T if two lists of polynomials G1 and G2, assumed to be Grobner bases,
generate  the same ideal, and NIL otherwise."
  (and (grobner-subsetp g1 g2 pred ring)
       (grobner-subsetp g2 g1 pred ring)))

(defun grobner-subsetp (g1 g2 &optional (pred #'lex>) (ring *coefficient-ring*))
  "Returns T if a list of polynomials G1 generates
an ideal contained in the ideal generated by a polynomial list G2,
both G1 and G2 assumed to be Grobner bases. Returns NIL otherwise."
  (every #'(lambda (p) (grobner-member p g2 pred ring)) g1))

(defun grobner-member (p g &optional (pred #'lex>) (ring *coefficient-ring*))
  "Returns T if a polynomial P belongs to the ideal generated by the
polynomial list G, which is assumed to be a Grobner basis. Returns NIL otherwise."
  (endp (normal-form (copy-tree p) g pred nil ring)))

(defun ideal-equal (f1 f2 pred top-reduction-only ring)
  "Returns T if two ideals generated by polynomial lists F1 and F2 are identical.
Returns NIL otherwise."
  (grobner-equal (grobner f1 pred 0 top-reduction-only ring)
		 (grobner f2 pred 0 top-reduction-only ring)
		 pred ring))

(defun ideal-subsetp (f1 f2 &optional (pred #'lex>) (ring *coefficient-ring*))
  "Returns T if the ideal spanned by the polynomial list F1 is contained
in the ideal spanned by the polynomial list F2. Returns NIL otherwise."
  (grobner-subsetp
   (grobner f1 pred 0 nil ring)
   (grobner f2 pred 0 nil ring)
   pred ring))

(defun ideal-member (p plist pred top-reduction-only ring)
  "Returns T if the polynomial P belongs to the ideal spanned by the polynomial
list PLIST, and NIL otherwise."
  (endp (normal-form (copy-tree p) (grobner plist pred 0 top-reduction-only ring)
		     pred nil ring)))

;; Calculate F : p^inf
(defun ideal-saturation-1 (F p pred start top-reduction-only ring
			   &aux (pred (elimination-order-1 pred)))
  "Returns the reduced Grobner basis of the saturation of the ideal
generated by a polynomial list F in the ideal generated by a single
polynomial P. The saturation ideal is defined as the set of
polynomials H such for some natural number n (* (EXPT P N) H) is in the ideal
F. Geometrically, over an algebraically closed field, this is the set
of polynomials in the ideal generated by F which do not identically
vanish on the variety of P."
  (mapcar
   #'poly-contract
   (ring-intersection
    (reduced-grobner
     (saturation-extension-1 F p ring)
     pred
     start top-reduction-only ring)
    1)))


(defun add-variables (plist n)
  "Add N new variables, adn the first N variables, to every polynomial in polynomial list PLIST.
The resulting polynomial will not depend on these N variables."
  (mapcar
   #'(lambda (g)
       (mapcar #'(lambda (q) (cons (append (make-list n :initial-element 0)
					   (car q))
				   (cdr q))) g))
   plist))

;; Calculate <u1*p1,u2*p2,...,uk*pk>
(defun extend-polynomials (plist &aux (k (length plist)))
  "Returns polynomial list {U1*P1', U2*P2', ... , UK*PK'}
where Ui are new variables and PLIST={P1, P2, ... , PK} is a polynomial
list. PI' is obtained from PI by adding new variables U1, U2, ..., UN
UK as the first K variables. Thus, the resulting list is consists of polynomials whose
terms are ordered by the lexicographic order on variables UI, with
ties resolved by the monomial order which was used to order the terms
of the polynomials in PLIST. We note that the monomial order does not
explicitly participate in this calculation, and neither do the
variable names."
  (labels ((incf-power (g i)
	     (dolist (x g g)
	       (incf (nth i (car x))))))
    (setf plist (add-variables plist k))
    (dotimes (i k plist)
      (setf (nth i plist) (incf-power (nth i plist) i)))))

;; Calculate <F, u1*p1-1,u2*p2-1,...,uk*pk-1>
(defun saturation-extension (F plist ring &aux (k (length plist))
					  (d (+ k (length (caaar plist)))))
  "Returns F' union {U1*P1-1,U2*P2-1,...,UK*PK-1} where Ui are new variables
and F' is a polynomial list  obtained from F by adding variables U1, U2, ..., Uk
as the first K variables to each polynomial in F. Thus, the resulting list is
consists of polynomials whose terms are ordered by the lexicographic order on
variables Ui, with ties resolved by the monomial order which was used to order the
terms of the polynomials in F. We note that the monomial order does not explicitly
participate in this calculation, and neither do the variable names."
  (setf F (add-variables F k)
	plist (mapcar #'(lambda (x)
			  (nconc x (list (cons (make-list d :initial-element 0)
					       (funcall (ring-- ring)
							(ring-unit ring))))))
		      (extend-polynomials plist)))
  (append F plist))

;; Calculate <F, u1*p1+u2*p2+...+uk*pk-1>
(defun polysaturation-extension (F plist ring &aux (k (length plist))
					      (d (+ k (length (caaar plist)))))
  "Returns F' union {U1*P1+U2*P2+UK*PK-1} where Ui are new variables
and F' is a polynomial list  obtained from F by adding variables U1, U2, ..., Uk
as the first K variables to each polynomial in F. Thus, the resulting list is
consists of polynomials whose terms are ordered by the lexicographic order on
variables Ui, with ties resolved by the monomial order which was used to order the
terms of the polynomials in F. We note that the monomial order does not explicitly
participate in this calculation, and neither do the variable names."
    (setf F (add-variables F k)
	  plist (apply #'append (extend-polynomials plist))
	  (cdr (last plist)) (list (cons (make-list d :initial-element 0)
					 (funcall (ring-- ring)
						  (ring-unit ring)))))
    (append F (list plist)))

(defun saturation-extension-1 (F p ring)
  "Return the list F' union {U*P-1} where U is a new variable,
where F' is obtained from the polynomial list F by adding U as the first variable.
The variable U will be added as the first variable, so that it can be easily
eliminated."
  (saturation-extension F (list p) ring))

;; Calculate F : p1^inf : p2^inf : ... : ps^inf
(defun ideal-polysaturation-1 (F plist pred start top-reduction-only ring)
  "Returns the reduced Grobner basis of the ideal obtained by a
sequence of successive saturations in the polynomials
of the polynomial list PLIST of the ideal generated by the
polynomial list F."
  (cond
   ((endp plist) (reduced-grobner F pred start top-reduction-only ring))
   (t (let ((G (ideal-saturation-1 F (car plist) pred start top-reduction-only ring)))
	(ideal-polysaturation-1 G (rest plist) pred (length G) top-reduction-only ring)))))

(defun ideal-saturation (F G pred start top-reduction-only ring
			 &aux (k (length G))
			      (pred (elimination-order k
						       :primary-order #'lex>
						       :secondary-order pred)))
  "Returns the reduced Grobner basis of the saturation of the ideal
generated by a polynomial list F in the ideal generated a polynomial
list G The saturation ideal is defined as the set of polynomials H
such for some natural number n and some P in the ideal generated by G
the polynomial P**N * H is in the ideal spanned by F.  Geometrically,
over an algebraically closed field, this is the set of polynomials in
the ideal generated by F which do not identically vanish on the
variety of G."
  (mapcar
   #'(lambda (q) (poly-contract q k))
   (ring-intersection
    (reduced-grobner  (polysaturation-extension  F G ring)
		      pred start
		      top-reduction-only ring)
    k)))

(defun ideal-polysaturation (F ideal-list pred start top-reduction-only ring)
    "Returns the reduced Grobner basis of the ideal obtained by a
successive applications of IDEAL-SATURATIONS to F and
lists of polynomials in the list IDEAL-LIST."
  (cond
   ((endp ideal-list) F)
   (t (let ((H (ideal-saturation F (car ideal-list) pred start top-reduction-only ring)))
	(ideal-polysaturation H (rest ideal-list) pred (length H) top-reduction-only ring)))))

(defun buchberger-criterion (G pred ring)
  "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
   #'endp
   (makelist (normal-form (spoly (elt G i) (elt G j) pred ring) G pred nil ring)
	     (i 0 (- (length G) 2))
	     (j (1+ i) (1- (length G))))))


(defun grobner-test (G F pred ring)
  "Test whether G is a Grobner basis and F is contained in G. Return T
upon success and NIL otherwise."
  #+debug(debug-cgb "~&GROBNER CHECK: ")
  (let ((*grobner-debug* nil)
	(stat1 (buchberger-criterion G pred ring))
	(stat2
	  (every #'endp
		 (makelist (normal-form (copy-tree (elt F i)) G pred nil ring)
			   (i 0 (1- (length F)))))))
    (unless stat1 (error "~&Buchberger criterion failed."))
    (unless stat2
      (error "~&Original polys not in ideal spanned by Grobner.")))
  #+debug(debug-cgb "~&GROBNER CHECK END")
  T)


(defun minimization (P pred)
  "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."
  (declare (ignore pred))
  (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 (lm x) (lm p))) Q1)
		  (setf found t)
		  (return Q1)))))))


(defun add-minimized (f Gred pred)
  "Adds a polynomial f to GRED, reduced Grobner basis, preserving the
property described in the documentation for MINIMIZATION."
  (if (member-if #'(lambda (p) (monom-divides-p (lm p) (lm f))) Gred)
      Gred
    (merge 'list (list f)
	   (delete-if #'(lambda (p) (monom-divides-p (lm f) (lm p))) Gred)
	   pred :key #'lm)))

(defun colon-ideal (F G pred top-reduction-only ring)
  "Returns the reduced Grobner basis of the colon ideal Id(F):Id(G), where
F and G are two lists of polynomials. The colon ideal I:J is defined as the
set of polynomials H such that there is a polynomial W in J for which
W*H belongs to I."
  (cond
    ((endp G)
     (if (every #'endp F) nil
       (list (list
	      (cons
	       (make-list (length (caar (find-if-not #'endp F)))
			  :initial-element 0)
	       1)))))
    ((endp (cdr G))
     (colon-ideal-1 F (car G) pred top-reduction-only ring))
    (t
     (ideal-intersection (colon-ideal-1 F (car G) pred top-reduction-only ring)
			 (colon-ideal F (rest G) pred top-reduction-only ring)
			 pred
			 top-reduction-only
			 ring))))

(defun colon-ideal-1 (F g pred top-reduction-only ring)
  "Returns the reduced Grobner basis of the colon ideal Id(F):Id({G}), where
F is a list of polynomials and G is a polynomial."
  (mapcar #'(lambda (x) (poly-exact-divide x g pred ring))
	  (ideal-intersection F (list g) pred top-reduction-only ring)))


(defun pseudo-divide (f fl pred ring)
  "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 value is an integer
count of the number of reductions performed. Finally, the fourth
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.
The resulting objects satisfy the equation:
  C*F= sum A[i]*FL[i] + R
"
  (declare (optimize (speed 3) (safety 0)))
  (do (r
       (c (ring-unit ring))
       (a (make-list (length fl)))
       (division-count 0)
       (p f))
      ((endp p)
       #+debug(debug-cgb "~&~3T~d reduction~:p" division-count)
       #+debug(when (endp r) (debug-cgb " ---> 0"))
       (values a (nreverse r) division-count c))
    (declare (fixnum division-count) (integer c))
    (do ((fl fl (rest fl))		;scan list of divisors
	 (b a (rest b)))
	((cond
	   ((endp fl)			;no division occurred
	    (setf r (cons (lt p) r) ;move lt(p) to remainder
		  p (rest p))		;remove lt(p) from p
	    t)
	   ((monom-divides-p (lm (car fl)) (lm p)) ;division occurred
	    (incf division-count)
	    (let* ((gcd (funcall (ring-gcd ring) (lc (car fl)) (lc p)))
		   (c1 (funcall (ring-/ ring) (lc (car fl)) gcd))
		   (c2 (funcall (ring-/ ring) (lc p) gcd))
		   (quot (cons (monom/ (lm p) (lm (car fl))) c2)))
	      ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
	      (mapl #'(lambda (x)
			(setf (car x) (scalar-times-poly c1 (car x) ring)))
		    a)
	      (setf p (scalar-times-poly c1 (rest p) ring) ;drop first term
		    r (scalar-times-poly c1 r ring)
		    c (funcall (ring-* ring) c c1)
		    p (poly-op p quot (rest (car fl)) pred ring)
		    (car b) (nconc (car b) (list quot))))
	    t))))))

;;----------------------------------------------------------------
;; An implementation of the algorithm of Gebauer and Moeller,
;; as described in the book of Becker-Weispfenning, p. 232
;;----------------------------------------------------------------
(defun gebauer-moeller (F pred start top-reduction-only ring
			&aux B G F1)
  "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''"
  ;;cut startup costs if (length F)<=1
  (cond
   ((endp F) (return-from gebauer-moeller nil))
   ((endp (cdr F))
    (return-from gebauer-moeller (list (grobner-primitive-part (car F) ring)))))
  #+debug (debug-cgb "~&GROBNER BASIS - GEBAUER MOELLER ALGORITHM")
  #+debug (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
  #+grobner-check  (when (plusp start)
		     (grobner-test (subseq F 0 start) (subseq F 0 start)
				   pred ring))
  (setf B nil
	G (subseq F 0 start)
	F1 (subseq F start))
  (do () ((endp F1))
    (multiple-value-setq (G B)
      (update G B (grobner-primitive-part (pop F1) ring) pred ring)))
  (do () ((endp B))
    (let* ((pair (pop B))
	   (g1 (first pair))
	   (g2 (second pair))
	   (h (normal-form (spoly g1 g2 pred ring)
			   G pred top-reduction-only ring)))
      (unless (poly-zerop h)
	(setf h (grobner-primitive-part h ring))
	(multiple-value-setq (G B)
	  (update G B h pred ring))
	#+debug (debug-cgb "~&Polynomials: ~d; Pairs left: ~d~%"
			   (length G) (length B))
	)))
  #+grobner-check(grobner-test G F pred ring)
  #+debug(debug-cgb "~&GROBNER END")
  G)

(defun update (G B h pred ring &aux C D E B-new G-new pair)
  "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."
  (setf C G D nil) 
  (do (g1) ((endp C))
    (setf g1 (pop C))
    (when (or (monom-rel-prime (lm h) (lm g1))
	      (and
	       (not (some
		     #'(lambda (g2) (monom-divides-p (monom-lcm (lm h) (lm g2))
						     (monom-lcm (lm h) (lm g1))))
		     C))
	       (not (some
		     #'(lambda (g2) (monom-divides-p (monom-lcm (lm h) (lm g2))
						     (monom-lcm (lm h) (lm g1))))
		     D))))
      (push g1 D)))
  (setf E nil)
  (do (g1) ((endp D))
    (setf g1 (pop D))
    (unless (monom-rel-prime (lm h) (lm g1))
      (push g1 E)))
  (setf B-new nil)
  (do (g1 g2) ((endp B))
    (setf pair (pop B)
	  g1 (first pair)
	  g2 (second pair))
    (when (or (not (monom-divides-p (lm h) (monom-lcm (lm g1) (lm g2))))
	      (equal (monom-lcm (lm g1) (lm h))
		     (monom-lcm (lm g1) (lm g2)))
	      (equal (monom-lcm (lm h) (lm g2))
		     (monom-lcm (lm g1) (lm g2))))
      (setf B-new (funcall *gebauer-moeller-merge-pairs* B-new (list (list g1 g2)) pred ring))))
  (dolist (g3 E)
    (setf B-new (funcall *gebauer-moeller-merge-pairs* B-new (list (list h g3)) pred ring)))
  (setf G-new nil)
  (do (g1) ((endp G))
    (setf g1 (pop G))
    (unless (monom-divides-p (lm h) (lm g1))
      (push g1 G-new)))
  (push h G-new)
  (values G-new B-new))
      
(defun gebauer-moeller-merge-pairs-use-mock-spoly (B C pred ring)
  "The merging strategy used by the Gebauer-Moeller algorithm, based
on MOCK-SPOLY; see the documentation of
BUCHBERGER-MERGE-PAIRS-USE-MOCK-SPOLY."
  (declare (ignore ring))
  ;; Delete pairs consisting of two single monomials
  (setf C (delete-if
	   #'(lambda (p) (and (endp (cdr (first p)))
			      (endp (cdr (second p)))))
	   C))
  (labels ((pair-key (p) (mock-spoly (first p) (second p) pred))
	   (pair-order (x y) (funcall pred y x)))
    (merge 'list B C #'pair-order :key #'pair-key)))

(defun gebauer-moeller-merge-pairs-smallest-lcm (B C pred ring)
  "The merging strategy based on the smallest-lcm (normal strategy); see 
the documentation of BUCHBERGER-MERGE-PAIRS-SMALLEST-LCM."
  (declare (ignore ring))
  ;; Delete pairs consisting of two single monomials
  (labels ((pair-key (p) (monom-lcm (lm (first p))
				    (lm (second p))))
	   (pair-order (x y) (funcall pred y x)))
    (merge 'list B C #'pair-order :key #'pair-key)))


(defun gebauer-moeller-merge-pairs-use-smallest-degree (B C pred ring)
  "The merging strategy based on the smallest-lcm (normal strategy); see 
the documentation of BUCHBERGER-MERGE-PAIRS-USE-SMALLEST-DEGREE."
  (declare (ignore ring))
  (declare (ignore pred))
  (labels ((pair-key (p) (total-degree
			  (monom-lcm (lm (first p))
				     (lm (second p)))))
	   (pair-order (x y) (funcall #'< x y)))
    (merge 'list B C #'pair-order :key #'pair-key)))


(defun gebauer-moeller-merge-pairs-use-smallest-length (B C pred ring)
  (declare (ignore ring))
  "The merging strategy based on the smallest-lcm (normal strategy); see 
the documentation of BUCHBERGER-MERGE-PAIRS-USE-SMALLEST-LENGTH."
  (declare (ignore pred))
  (labels ((pair-key (p) (+
			  (length (first p))
			  (length (second p))))
	   (pair-order (x y) (funcall #'< x y)))
    (if B
	(merge 'list B C #'pair-order :key #'pair-key)
      (sort C #'pair-order :key #'pair-key))))


(defun gebauer-moeller-merge-pairs-use-smallest-coefficient-length (B C pred ring)
  "The merging strategy based on the smallest-lcm (normal strategy); see 
the documentation of BUCHBERGER-MERGE-PAIRS-USE-SMALLEST-COEFFICIENT-LENGTH."
  (declare (ignore pred))
  (labels ((coefficient-length (poly)
	     (apply #'+ (mapcar (ring-length ring) (mapcar #'cdr poly))))
	   (pair-key (p) (+
			  (coefficient-length (first p))
			  (coefficient-length (second p))))
	   (pair-order (x y) (funcall #'< x y)))
    (merge 'list B C #'pair-order :key #'pair-key)))

(defun gebauer-moeller-set-pair-heuristic
    (method
     &aux (strategy-fn
	(ecase method
	      (:minimal-mock-spoly #'gebauer-moeller-merge-pairs-use-mock-spoly)
	      (:minimal-lcm #'gebauer-moeller-merge-pairs-smallest-lcm)
	      (:minimal-total-degree #'gebauer-moeller-merge-pairs-use-smallest-degree)
	      (:minimal-length #'gebauer-moeller-merge-pairs-use-smallest-length)
	      (:minimal-coefficient-length #'gebauer-moeller-merge-pairs-use-smallest-coefficient-length))))
  (setf *gebauer-moeller-merge-pairs* strategy-fn))

;;------------------------------------------------------------------------------------------------
;;
;; An implementation of the sugar strategy
;;
;;------------------------------------------------------------------------------------------------

;; Calculate sugar of the S-polynomial of f and g
(defun spoly-sugar (f-with-sugar g-with-sugar ring)
  "Calculate the ``sugar'' of the S-polynomial of two polynomials with
sugar. A polynomial with sugar is simply a pair (P . SUGAR) where
SUGAR is an integer constant defined according to the following
algorighm: the sugar of sum or difference of two polynomials with
sugar is the MAX of the sugars of those two polynomials. The sugar of
a product of a term and a polynomial is the sum of the degree of the
term and the sugar of the polynomial.  The idea is to accumulate sugar
as we perform the arithmetic operations, and that polynomials with
small (little) sugar should be given priority in the calculations.
Thus, the ``sugar strategy'' of the critical pair selection is to
select a pair with the smallest value of sugar of the resulting
S-polynomial."
  (let* ((f (poly-with-sugar-poly f-with-sugar))
	 (g (poly-with-sugar-poly g-with-sugar))
	 (lcm (monom-lcm (lm f) (lm g)))
	 (m1 (monom/ lcm (lm f)))
	 (m2 (monom/ lcm (lm g))))
    (let* ((c (funcall (ring-lcm ring) (lc f) (lc g)))
	   (cf (funcall (ring-/ ring) c (lc f)))
	   (cg (funcall (ring-/ ring) c (lc g))))
      (max
       (+ (term-sugar (cons m1 cf) ring) (poly-with-sugar-sugar f-with-sugar))
       (+ (term-sugar (cons m2 cg) ring) (poly-with-sugar-sugar g-with-sugar))))))

;; Calculate the S-polynomial of f and g with sugar
(defun spoly-with-sugar (f-with-sugar g-with-sugar pred ring)
  "The S-polynomials of two polynomials with SUGAR strategy."
  (let* ((f (poly-with-sugar-poly f-with-sugar))
	 (g (poly-with-sugar-poly g-with-sugar))
	 (lcm (monom-lcm (lm f) (lm g)))
	 (m1 (monom/ lcm (lm f)))
	 (m2 (monom/ lcm (lm g))))
    (let* ((c (funcall (ring-lcm ring) (lc f) (lc g)))
	   (cf (funcall (ring-/ ring) c (lc f)))
	   (cg (funcall (ring-/ ring) c (lc g))))
      (poly-with-sugar-
       (term-times-poly-with-sugar (cons m1 cf) (poly-with-sugar-tail f-with-sugar) ring)
       (term-times-poly-with-sugar (cons m2 cg) (poly-with-sugar-tail g-with-sugar) ring)
       pred
       ring))))


(defun normal-form-with-sugar (f fl pred top-reduction-only ring)
  "Normal form of the polynomial with sugar F with respect to
a list of polynomials with sugar FL. Assumes that FL is not empty.
Parameters PRED, TOP-REDUCTION-ONLY and RING have the same meaning as in NORMAL-FORM.
The parameter SUGAR-LIMIT should be set to a positive integer. If the sugar limit of
the partial remainder exceeds SUGAR-LIMIT then the calculation is stopped and the partial
remainder is returned, although it is not fully reduced with respect to FL."
  (declare (optimize (speed 3) (safety 0)))
  ;; Loop invariant: c*f=sum ai*fi+r+p
  (do ((r (cons nil 0))
       (c (ring-unit ring))
       (division-count 0)
       (p f))
      ((or (poly-with-sugar-zerop p)
	   (and top-reduction-only (not (poly-with-sugar-zerop r))))
       #+debug
       (progn
	 (debug-cgb "~&~3T~d reduction~:p" division-count)
	 (cond
	  ((poly-with-sugar-zerop r)
	   (debug-cgb " ---> 0"))))
       (values (poly-with-sugar-append (poly-with-sugar-nreverse r) p)
	       division-count c))
    (declare (fixnum division-count) (integer c))
    (let ((g (find (poly-with-sugar-lm p) fl
		   :test #'monom-divisible-by-p
		   :key #'poly-with-sugar-lm)))
      (cond 
       (g				;division possible
	(incf division-count)
	(let* ((gcd (funcall (ring-gcd ring) (poly-with-sugar-lc g) (poly-with-sugar-lc p)))
	       (c1 (funcall (ring-/ ring) (poly-with-sugar-lc g) gcd))
	       (c2 (funcall (ring-/ ring) (poly-with-sugar-lc p) gcd))
	       (m (monom/ (poly-with-sugar-lm p) (poly-with-sugar-lm g))))
	  ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
	  (setf r (scalar-times-poly-with-sugar c1 r ring)
		c (funcall (ring-* ring) c c1)
		p (scalar-times-poly-with-sugar c1 p ring)
		p (poly-with-sugar-op p (cons m c2) g pred ring))))
       (t				;no division possible
	;;move lt(p) to remainder
	(setf r (cons (cons (poly-with-sugar-lt p) (poly-with-sugar-poly r))
		      (poly-with-sugar-sugar r))
	      p (poly-with-sugar-tail p)))))))

;; Parameter START is used when Grobner basis is calculated incrementally
;; and it signals that the elements of F from 0 to START-1 form a Grobner
;; basis already, so their pairs are not formed.
(defun buchberger-with-sugar (F-no-sugar pred start top-reduction-only ring
			      &aux (s (1- (length F-no-sugar))) B M F)
  "Returns a Grobner basis of an ideal generated by a list of ordinary
polynomials F-NO-SUGAR. Adds initial sugar to the polynomials and
performs the Buchberger algorithm with ``sugar strategy''. It returns
an ordinary list of polynomials with no sugar. One of the most
effective algorithms for Grobner basis calculation."
  (declare (fixnum s))
  (unless (plusp s) (return-from buchberger-with-sugar F-no-sugar)) ;cut startup costs
  #+debug (debug-cgb "~&GROBNER BASIS - BUCHBERGER WITH SUGER STRATEGY ALGORITHM")
  #+debug (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
  #+grobner-check  (when (plusp start)
		     (grobner-test (subseq F-no-sugar 0 start)
				   (subseq F-no-sugar 0 start) pred ring))
  (setf F (mapcar #'(lambda (x) (poly-add-sugar x ring)) F-no-sugar)
	B (nconc (makelist (list i j) (i 0 (1- start)) (j start s))
		 (makelist (list i j) (i start (1- s)) (j (1+ i) s)))
	M (make-hash-table :test #'equal))
  ;;Initialize treated pairs
  (dotimes (i (1- start))
    (do ((j (1+ i) (1+ j))) ((>= j start))
      (setf (gethash (list i j) M) t)))
  (setf B (buchberger-with-sugar-sort-pairs B F M pred ring))
  (do ((G F))
      ((endp B)
       (setf G (mapcar #'poly-with-sugar-poly G))
       ;;#+grobner-check(grobner-test G F-no-sugar pred ring)
       #+debug(debug-cgb "~&GROBNER END")
       G)
    (let ((pair (pop B)))
      (cond
       ((Criterion-1-with-sugar pair G))
       ((Criterion-2-with-sugar pair G M))
       (t
	(let ((SP (normal-form-with-sugar
		   (cons (spoly (poly-with-sugar-poly (elt G (first pair)))
				(poly-with-sugar-poly (elt G (second pair))) pred ring)
			 (third pair))
		   G pred top-reduction-only ring)))
	  (cond
	   ((poly-with-sugar-zerop SP))
	   (t
	    (setf s (1+ s)
		  (poly-with-sugar-poly SP) (grobner-primitive-part
					     (poly-with-sugar-poly SP) ring)
		  (cdr (last G)) (list SP)
		  B (buchberger-with-sugar-merge-pairs B (makelist (list i s)
								   (i 0 (1- s)))
						       G M pred ring))
	    #+debug (debug-cgb "~&[~d] Polynomials: ~d; Pairs left: ~d; Treated pairs: ~d"
			       (third pair)
			       (length G) (length B)
			       (hash-table-count M)))))))
      (setf (cddr pair) nil
	    (gethash pair M) t))))


(defun buchberger-with-sugar-merge-pairs (B C G M pred ring)
  "Merges lists of critical pairs. It orders pairs according to
increasing sugar, with ties broken by smaller MOCK-SPOLY value of the
two polynomials. In this function B must already be sorted."
  ;;Remove pairs consisting of two monomials
  ;;These pairs produce zero S-polynomial
  (setf C (delete-if
	   #'(lambda (p)
	       (when (and (endp (cdr (poly-with-sugar-poly (elt G (first p)))))
			  (endp (cdr (poly-with-sugar-poly (elt G (second p))))))
		 ;;Record the pair as treated
		 (setf (gethash p M) t)
		 t))
	   C))
  (labels ((pair-sugar (p)
	     (spoly-sugar (elt G (first p)) (elt G (second p)) ring))
	   (pair-mock-spoly (p)
	     (mock-spoly (poly-with-sugar-poly (elt G (first p)))
			 (poly-with-sugar-poly (elt G (second p))) pred))
	   (pair-order (p q)
	     (let ((sugar-p (third p))
		   (sugar-q (third q)))
	       (cond
		((< sugar-p sugar-q) t)
		((> sugar-p sugar-q) nil)
		(t
		 (funcall pred (fourth q) (fourth p)))))))
    ;;Add sugar and mock S-polynomial as the third and fourth element
    ;;of new pair in C
    (dolist (p C)
      (setf (cdr (last p))
	    (list (pair-sugar p) (pair-mock-spoly p))))
    (merge 'list B (sort C #'pair-order) #'pair-order)))

(defun buchberger-with-sugar-sort-pairs (C G M pred ring) 
  "Sorts critical pairs C according to sugar strategy."
  (buchberger-with-sugar-merge-pairs nil C G M pred ring))

(defun Criterion-1-with-sugar (pair G &aux (i (first pair)) (j (second pair)))
  "An implementation of Buchberger's first criterion for polynomials 
with sugar. See the documentation of BUCHBERGER-CRITERION-1."
  (monom-rel-prime (poly-with-sugar-lm (elt G i))
		   (poly-with-sugar-lm (elt G j))))

(defun Criterion-2-with-sugar (pair G M &aux (i (first pair)) (j (second pair)))
  "An implementation of Buchberger's first criterion for polynomials 
with sugar. See the documentation of BUCHBERGER-CRITERION-2."
  (labels ((make-pair (u v) (if (< u v) (list u v) (list v u))))
    (dotimes (k (length G) nil)
      (when (and (/= k i)
		 (/= k j)
		 (monom-divides-p
		  (poly-with-sugar-lm (elt G k))
		  (monom-lcm (poly-with-sugar-lm (elt G i))
			     (poly-with-sugar-lm (elt G j))))
		 (gethash (make-pair i k) M)
		 (gethash (make-pair k j) M))
	#+debug(when *grobner-debug* (format t ":B2"))
	(return-from Criterion-2-with-sugar t)))))


;;----------------------------------------------------------------
;; An implementation of the algorithm of Gebauer and Moeller,
;; as described in the book of Becker-Weispfenning, p. 232
;; with sugar strategy
;;----------------------------------------------------------------
(defun gebauer-moeller-with-sugar (F pred start top-reduction-only ring
			&aux B G F1)
  "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''"
  ;;cut startup costs if (length F)<=1
  (cond
   ((endp F) (return-from gebauer-moeller-with-sugar nil))
   ((endp (cdr F))
    (return-from gebauer-moeller-with-sugar (list (grobner-primitive-part (car F) ring)))))
  #+debug (debug-cgb "~&GROBNER BASIS - GEBAUER MOELLER ALGORITHM WITH SUGAR STRATEGY")
  #+debug (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
  #+grobner-check  (when (plusp start)
		     (grobner-test (subseq F 0 start) (subseq F 0 start)
				   pred ring))
  (setf B nil
	G (subseq F 0 start)
	F1 (subseq F start))
  ;;Initialize sugar
  (setf G (mapcar #'(lambda (x) (poly-add-sugar x ring)) G)
	F1 (mapcar #'(lambda (x) (poly-add-sugar x ring)) F1))
  (do () ((endp F1))
    (multiple-value-setq (G B)
      (update-with-sugar G B (grobner-primitive-part-with-sugar (pop F1) ring) pred ring)))
  (do ()
      ((endp B))
    (let* ((pair (pop B))
	   (g1 (first pair))
	   (g2 (second pair))
	   (h (normal-form-with-sugar (spoly-with-sugar g1 g2 pred ring)
				      G pred top-reduction-only ring)))
      (cond
       ((poly-with-sugar-zerop h))
       (t
	(setf h (grobner-primitive-part-with-sugar h ring))
	(multiple-value-setq (G B)
	  (update-with-sugar G B h pred ring))
	#+debug (debug-cgb "~&[[~d]] Polynomials: ~d; Pairs left: ~d~%"
			   (poly-with-sugar-sugar h)
			   (length G) (length B))
	))))
  (setf G (mapcar #'poly-with-sugar-poly G))
  ;;#+grobner-check(grobner-test G F pred ring)
  #+debug(debug-cgb "~&GROBNER END")
  G)

(defun update-with-sugar (G B h pred ring &aux C D E B-new G-new pair)
  "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. Operates on polynomials with sugar."
  (setf C G D nil) 
  (do (g1) ((endp C))
    (setf g1 (pop C))
    (when (or (monom-rel-prime (poly-with-sugar-lm h) (poly-with-sugar-lm g1))
	      (and
	       (not (some
		     #'(lambda (g2) (monom-divides-p (monom-lcm (poly-with-sugar-lm h)
								(poly-with-sugar-lm g2))
						     (monom-lcm (poly-with-sugar-lm h)
								(poly-with-sugar-lm g1))))
		     C))
	       (not (some
		     #'(lambda (g2) (monom-divides-p (monom-lcm (poly-with-sugar-lm h)
								(poly-with-sugar-lm g2))
						     (monom-lcm (poly-with-sugar-lm h)
								(poly-with-sugar-lm g1))))
		     D))))
      (push g1 D)))
  (setf E nil)
  (do (g1) ((endp D))
    (setf g1 (pop D))
    (unless (monom-rel-prime (poly-with-sugar-lm h) (poly-with-sugar-lm g1))
      (push g1 E)))
  (setf B-new nil)
  (do (g1 g2) ((endp B))
    (setf pair (pop B)
	  g1 (first pair)
	  g2 (second pair))
    (when (or (not (monom-divides-p (poly-with-sugar-lm h) (monom-lcm (poly-with-sugar-lm g1)
								      (poly-with-sugar-lm g2))))
	      (equal (monom-lcm (poly-with-sugar-lm g1) (poly-with-sugar-lm h))
		     (monom-lcm (poly-with-sugar-lm g1) (poly-with-sugar-lm g2)))
	      (equal (monom-lcm (poly-with-sugar-lm h) (poly-with-sugar-lm g2))
		     (monom-lcm (poly-with-sugar-lm g1) (poly-with-sugar-lm g2))))
      (setf B-new (gebauer-moeller-with-sugar-merge-pairs B-new (list (list g1 g2)) pred ring))))
  (dolist (g3 E)
    (setf B-new (gebauer-moeller-with-sugar-merge-pairs B-new (list (list h g3)) pred ring)))
  (setf G-new nil)
  (do (g1) ((endp G))
    (setf g1 (pop G))
    (unless (monom-divides-p (poly-with-sugar-lm h) (poly-with-sugar-lm g1))
      (push g1 G-new)))
  (push h G-new)
  (values G-new B-new))

(defun gebauer-moeller-with-sugar-merge-pairs (B C pred ring)
  "Merges lists of critical pairs. It orders pairs according to
increasing sugar, with ties broken by smaller lcm of head monomials In
this function B must already be sorted. Operates on polynomials with sugar."
  (setf C (delete-if
	   #'(lambda (p) (and (endp (cdr (first p)))
			      (endp (cdr (second p)))))
	   C))
  (labels ((pair-sugar (p)
	     (spoly-sugar (first p) (second p) ring))
	   (pair-mock-spoly (p)
	     (mock-spoly (poly-with-sugar-poly (first p))
			 (poly-with-sugar-poly (second p)) pred))
	   (pair-order (p q)
	     (let ((sugar-p (pair-sugar p))
		   (sugar-q (pair-sugar q)))
	       (cond
		((< sugar-p sugar-q) t)
		((> sugar-p sugar-q) nil)
		(t
		 (funcall pred (pair-mock-spoly q) (pair-mock-spoly p)))))))
    (merge 'list B (sort C #'pair-order) #'pair-order)))

(defun grobner-primitive-part-with-sugar (h ring)
  (cons (grobner-primitive-part (car h) ring) (cdr h)))

;;----------------------------------------------------------------
;;
;; An implementation of the algorithms using product caching
;;
;;----------------------------------------------------------------
(defun cached-normal-form (f fl pred  top-reduction-only ring cache)
  "This version of NORMAL-FORM uses a cache which stores products of monomials
and polynomials found in FL. The cache is a hash table storing pairs (G . G-HASH-TABLE)
where G is a polynomial in FL serving as a key. The hash table G-HASH-TABLE is a list
of pairs (MONOMIAL . PRODUCT-POLY) where PRODUCT-POLY is (POLY* MONOMIAL G)."

  (declare (optimize (speed 3) (safety 0)))
  ;; Loop invariant: c*f=sum ai*fi+r+p
  (do (r (c (ring-unit ring))
       (division-count 0)
       (p f))
      ((or (endp p)
	   (endp fl)
	   (and top-reduction-only r))
       #+debug(progn
		(debug-cgb "~&~3T~d reduction~:p" division-count)
		(when (endp r)
		  (debug-cgb " ---> 0")))
       (values (append (nreverse r) p)
	       division-count c))
    (declare (fixnum division-count) (integer c))
    (let ((i (position (lm p) fl :test #'monom-divisible-by-p :key #'lm)))
      (cond 
       (i				;division possible
	(incf division-count)
	(let* ((g (nth i fl))
	       (hash-tab (nth i cache))
	       (gcd (funcall (ring-gcd ring) (lc g) (lc p)))
	       (c1 (funcall (ring-/ ring) (lc g) gcd))
	       (c2 (funcall (ring-/ ring) (lc p) gcd))
	       (m (monom/ (lm p) (lm g))))
	  ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
	  (setf r (scalar-times-poly c1 r ring)
		c (funcall (ring-* ring) c c1)
		p (cached-grobner-op c2 c1 m p g pred ring hash-tab))))
       (t				;no division possible
	(setf r (cons (lt p) r)		;move lt(p) to remainder
	      p (rest p)))))))

(defun cached-grobner-op (c1 c2 m A B pred ring cache)
  "Returns the same result as (GROBNER-OP C1 C2 M A B PRED RING) but
it used variable CACHE which is described in CACHED-NORMAL-FORM for
efficient calculation of products of monomials with polynomials."
  (poly- (scalar-times-poly c2 A)
	 (scalar-times-poly c1 (cached-monom-times-poly m B cache) ring)
	 pred ring))

(defun cached-buchberger (F pred start top-reduction-only ring 
			    &aux (s (1- (length F))) B M CACHE)
  "It returns the same result as BUCHBERGER but it uses caching
of products of monomials with polynomials."
  (declare (fixnum s))
  (unless (plusp s) (return-from cached-buchberger F)) ;cut startup costs
  #+debug (debug-cgb "~&GROBNER BASIS - CACHED BUCHBERGER ALGORITHM")
  #+debug (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
  #+grobner-check  (when (plusp start)
		     (grobner-test (subseq F 0 start) (subseq F 0 start) pred ring))
  (setf B (nconc (makelist (list i j) (i 0 (1- start)) (j start s))
		 (makelist (list i j) (i start (1- s)) (j (1+ i) s)))
	M (make-hash-table :test #'equal)
	B (buchberger-sort-pairs B F pred ring))
  ;;Initialize cache 
  (dotimes (k (1+ s)) (push (make-hash-table :test #'equal) CACHE))
  ;;Initialize treated pairs
  (dotimes (i (1- start))
    (do ((j (1+ i) (1+ j))) ((>= j start))
      (setf (gethash (list i j) M) t)))
  (do ((G F))
      ((endp B)
       #+grobner-check(grobner-test G F pred ring)
       #+debug(debug-cgb "~&GROBNER END")
       G)
    (let ((pair (pop B)))
      (unless (or (Criterion-1 pair G)
		  (Criterion-2 pair G M))
	(let ((SP (cached-normal-form (cached-spoly (elt G (first pair))
						    (elt G (second pair))
						    pred
						    ring
						    (elt CACHE (first pair))
						    (elt CACHE (second pair)))
				      G pred top-reduction-only ring CACHE)))
	  (unless (poly-zerop SP)	;S-polynomial ---> 0;
	    (setf s (1+ s)
		  SP (grobner-primitive-part SP ring)
		  (cdr (last G)) (list SP)
		  (cdr (last Cache)) (list (make-hash-table :test #'equal))
		  B (funcall *buchberger-merge-pairs* B (makelist (list i s)
								  (i 0 (1- s)))
			     G pred ring))
	    #+debug (debug-cgb "~&Polynomials: ~d; Pairs left: ~d; Treated pairs: ~d"
			       (length G) (length B)
			       (hash-table-count M)))))
      (setf (gethash pair M) t))))

(defun cached-spoly (f g pred ring cache-f cache-g)
  "Uses CACHE to fetch products of monomials and polynomials."
  (declare (optimize (speed 3) (safety 0)))
  (let* ((lcm (monom-lcm (lm f) (lm g)))
	 (m1 (monom/ lcm (lm f)))
	 (m2 (monom/ lcm (lm g))))
    (let* ((c (funcall (ring-lcm ring) (lc f) (lc g)))
	   (cf (funcall (ring-/ ring) c (lc f)))
	   (cg (funcall (ring-/ ring) c (lc g))))
      (poly-
       (scalar-times-poly cf (cached-monom-times-poly m1 f cache-f) ring)
       (scalar-times-poly cg (cached-monom-times-poly m2 g cache-g) ring)
       pred
       ring))))


(defun cached-monom-times-poly (m f cache)
  (let ((prod (gethash m cache)))
    (cond
     (prod #+debug(debug-cgb "c") prod)
     (t
      #+debug(debug-cgb ".")
      (setf (gethash m cache) (monom-times-poly m f))))))
    
;;------------------------------------------------------------------------------------------------
;;
;; An implementation of the sugar strategy with product caching
;;
;;------------------------------------------------------------------------------------------------
(defun cached-normal-form-with-sugar (f fl pred top-reduction-only ring cache)
  "Normal form of the polynomial with sugar F with respect to
a list of polynomials with sugar FL. Assumes that FL is not empty.
Parameters PRED, TOP-REDUCTION-ONLY and RING have the same meaning as in NORMAL-FORM.
The parameter SUGAR-LIMIT should be set to a positive integer. If the sugar limit of
the partial remainder exceeds SUGAR-LIMIT then the calculation is stopped and the partial
remainder is returned, although it is not fully reduced with respect to FL."
  (declare (optimize (speed 3) (safety 0)))
  ;; Loop invariant: c*f=sum ai*fi+r+p
  (do ((r (cons nil 0))
       (c (ring-unit ring))
       (division-count 0)
       (p f))
      ((or (poly-with-sugar-zerop p)
	   (and top-reduction-only (not (poly-with-sugar-zerop r))))
       #+debug
       (progn
	 (debug-cgb "~&~3T~d reduction~:p" division-count)
	 (cond
	  ((poly-with-sugar-zerop r)
	   (debug-cgb " ---> 0"))))
       (values (poly-with-sugar-append (poly-with-sugar-nreverse r) p)
	       division-count c))
    (declare (fixnum division-count) (integer c))
    (let ((i (position (poly-with-sugar-lm p) fl
		   :test #'monom-divisible-by-p
		   :key #'poly-with-sugar-lm)))
      (cond 
       (i				;division possible
	(incf division-count)
	(let* ((g (nth i fl))
	       (hash-tab (nth i cache))
	       (gcd (funcall (ring-gcd ring) (poly-with-sugar-lc g) (poly-with-sugar-lc p)))
	       (c1 (funcall (ring-/ ring) (poly-with-sugar-lc g) gcd))
	       (c2 (funcall (ring-/ ring) (poly-with-sugar-lc p) gcd))
	       (m (monom/ (poly-with-sugar-lm p) (poly-with-sugar-lm g))))
	  ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
	  (setf r (scalar-times-poly-with-sugar c1 r ring)
		c (funcall (ring-* ring) c c1)
		p (scalar-times-poly-with-sugar c1 p ring)
		p (cached-poly-with-sugar-op p c2 m g pred ring hash-tab))))
       (t				;no division possible
	;;move lt(p) to remainder
	(setf r (cons (cons (poly-with-sugar-lt p) (poly-with-sugar-poly r))
		      (poly-with-sugar-sugar r))
	      p (poly-with-sugar-tail p)))))))

;; Parameter START is used when Grobner basis is calculated incrementally
;; and it signals that the elements of F from 0 to START-1 form a Grobner
;; basis already, so their pairs are not formed.
(defun cached-buchberger-with-sugar (F-no-sugar pred start top-reduction-only ring
						&aux (s (1- (length F-no-sugar))) B M F Cache)
  "Returns a Grobner basis of an ideal generated by a list of ordinary
polynomials F-NO-SUGAR. Adds initial sugar to the polynomials and
performs the Buchberger algorithm with ``sugar strategy''. It returns
an ordinary list of polynomials with no sugar. One of the most
effective algorithms for Grobner basis calculation."
  (declare (fixnum s))
  (unless (plusp s) (return-from cached-buchberger-with-sugar F-no-sugar)) ;cut startup costs
  #+debug (debug-cgb "~&GROBNER BASIS - BUCHBERGER WITH SUGER STRATEGY ALGORITHM")
  #+debug (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
  #+grobner-check  (when (plusp start)
		     (grobner-test (subseq F-no-sugar 0 start)
				   (subseq F-no-sugar 0 start) pred ring))
  (setf F (mapcar #'(lambda (x) (poly-add-sugar x ring)) F-no-sugar)
	B (nconc (makelist (list i j) (i 0 (1- start)) (j start s))
		 (makelist (list i j) (i start (1- s)) (j (1+ i) s)))
	M (make-hash-table :test #'equal))
  ;;Initialize cache 
  (dotimes (k (1+ s)) (push (make-hash-table :test #'equal) Cache))
  ;;Initialize treated pairs
  (dotimes (i (1- start))
    (do ((j (1+ i) (1+ j))) ((>= j start))
      (setf (gethash (list i j) M) t)))
  (setf B (buchberger-with-sugar-sort-pairs B F M pred ring))
  (do ((G F))
      ((endp B)
       (setf G (mapcar #'poly-with-sugar-poly G))
       ;;#+grobner-check(grobner-test G F-no-sugar pred ring)
       #+debug(debug-cgb "~&GROBNER END")
       G)
    (let ((pair (pop B)))
      (cond
       ((Criterion-1-with-sugar pair G))
       ((Criterion-2-with-sugar pair G M))
       (t
	(let ((SP (cached-normal-form-with-sugar
		   (cons (cached-spoly (poly-with-sugar-poly (elt G (first pair)))
				       (poly-with-sugar-poly (elt G (second pair)))
				       pred ring
				       (elt Cache (first pair))
				       (elt Cache (second pair)))
			 (third pair))
		   G pred top-reduction-only ring Cache)))
	  (cond
	   ((poly-with-sugar-zerop SP))
	   (t
	    (setf s (1+ s)
		  (poly-with-sugar-poly SP) (grobner-primitive-part
					     (poly-with-sugar-poly SP) ring)
		  (cdr (last G)) (list SP)
		  (cdr (last Cache)) (list (make-hash-table :test #'equal))
		  B (buchberger-with-sugar-merge-pairs B (makelist (list i s)
								   (i 0 (1- s)))
						       G M pred ring))
	    #+debug (debug-cgb "~&[~d] Polynomials: ~d; Pairs left: ~d; Treated pairs: ~d"
			       (third pair)
			       (length G) (length B)
			       (hash-table-count M)))))))
      (setf (cddr pair) nil
	    (gethash pair M) t))))

;; Implement f-term*g
(defun cached-poly-with-sugar-op (f c m g pred ring Cache)
  (poly-with-sugar- f (scalar-times-poly-with-sugar c (cached-monom-times-poly-with-sugar m g Cache) ring) pred ring))

(defun cached-monom-times-poly-with-sugar (m f Cache)
  (cons
   (cached-monom-times-poly m (poly-with-sugar-poly f) Cache)
   (+ (poly-with-sugar-sugar f) (monom-sugar m))))
@


1.7
log
@*** empty log message ***
@
text
@a86 3
(eval-when (compile)
  (proclaim '(inline grobner)))

a134 3
(eval-when (compile)
  (proclaim '(inline grobner-primitive-part grobner-content)))

a317 2
(eval-when (compile)
  (proclaim '(inline buchberger-sort-pairs)))
a346 6
(eval-when (compile)
  (proclaim '(inline
	      buchberger-merge-pairs-use-mock-spoly
	      buchberger-merge-pairs-smallest-lcm
	      buchberger-merge-pairs-use-smallest-degree)))

a441 3
(eval-when (compile)
  (proclaim '(inline Criterion-1 Criterion-2)))

a524 3
(eval-when (compile)
  (proclaim '(inline monom-depends-p term-depends-p poly-depends-p)))

a1300 3
(eval-when (compile)
  (proclaim '(inline Criterion-1-with-sugar Criterion-2-with-sugar)))

@


1.6
log
@*** empty log message ***
@
text
@a974 1
  (declare (ignore sorted))
@


1.5
log
@*** empty log message ***
@
text
@a1357 1
  (declare (ignore sorted))
@


1.4
log
@*** empty log message ***
@
text
@d57 2
a58 2
;;(proclaim '(optimize (speed 0) (debug 3)))
(proclaim '(optimize (speed 3) (debug 0)))
@


1.3
log
@*** empty log message ***
@
text
@d57 2
a58 1
(proclaim '(optimize (speed 0) (debug 3)))
@


1.2
log
@*** empty log message ***
@
text
@d57 1
a57 2
(eval-when (compile)
  (proclaim '(optimize (safety 0) (speed 3))))
@


1.1
log
@Initial revision
@
text
@d3 1
a3 1
	$Id: grobner.lisp,v 1.288 1998/01/18 09:06:58 marek Exp $
d15 1
a15 1
	"POLY-WITH-SUGAR" "LISP")
@