Index: branches/f4grobner/division.lisp
===================================================================
--- branches/f4grobner/division.lisp	(revision 4434)
+++ branches/f4grobner/division.lisp	(revision 4463)
@@ -43,21 +43,7 @@
 
 (defun grobner-op (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."
-  (declare (type monom m)
-	   (type poly f g))
-  (assert (universal-zerop
-	   (subtract
-	    (multiply c2 (leading-coefficient f))
-	    (multiply c1 (leading-coefficient g)))))
-  (assert (universal-equalp (leading-monomial f) (multiply m (leading-monomial g))))
-  ;; Note that below we can drop the leading terms of f ang g for the
-  ;; purpose of polynomial arithmetic.  
-  ;;
-  ;; TODO: Make sure that the sugar calculation is correct if leading
-  ;; terms are dropped.
-  (subtract
-   (multiply f c2)
-   (multiply g m c1)))
+  "Returns C2*F-C1*M*G, where F and G are polynomials M is a monomial."
+  (declare (type monom m) (type poly f g))
+  (subtract (multiply f c2) (multiply g m c1)))
 
 (defun check-loop-invariant (c f a fl r p &aux (p-zero (make-zero-for f)))
@@ -74,5 +60,5 @@
   |#
   (let* ((prod (inner-product a fl add multiply p-zero))
-	 (succeeded-p (universal-zerop (subtract (multiply f c) (add prod r p)))))
+	 (succeeded-p (universal-zerop (subtract (multiply f c) (add prod (make-instance 'poly :termlist (reverse r)) p)))))
     (unless succeeded-p 
       (error "#### Polynomial division Loop invariant failed ####:~%C=~A~%F=~A~%A=~A~%FL=~A~%R=~A~%P=~A~%" 
@@ -92,5 +78,5 @@
   (declare (type poly f) (list fl))
   ;; Loop invariant: c*f=sum ai*fi+r+p, where p must eventually become 0
-  (do ((r (make-zero-for f))
+  (do ((r nil)
        (c (make-unit-for (leading-coefficient f)))
        (a (make-list (length fl) :initial-element (make-zero-for f)))
@@ -100,6 +86,6 @@
        #+grobner-check(check-loop-invariant c f a fl r p)
        (debug-cgb "~&~3T~d reduction~:p" division-count)
-       (when (universal-zerop r) (debug-cgb " ---> 0"))
-       (values a r c division-count))
+       (when (null r) (debug-cgb " ---> 0"))
+       (values a (make-instance 'poly :termlist r) c division-count))
     (declare (fixnum division-count))
     ;; Check the loop invariant here
@@ -108,7 +94,6 @@
 	 (b a (rest b)))
 	((cond
-	   ((endp fl)				;no division occurred
-	    (setf r (add-to r (leading-term p)) ;move lt(p) to remainder
-		  p (subtract-from p (leading-term p)))	;remove lt(p) from p
+	   ((endp fl)			  ;no division occurred
+	    (push (poly-remove-term p) r) ;move lt(p) to remainder
 	    t)
 	   ((divides-p (leading-monomial (car fl)) (leading-monomial p)) ;division occurred
@@ -122,5 +107,5 @@
 			  (setf (car x) (multiply-by (car x) c1)))
 		      a)
-		(setf r (multiply-by r c1)
+		(setf r (mapc #'multiply-by r c1)
 		      c (multiply-by c c1)
 		      p (grobner-op c2 c1 m p (car fl)))
@@ -154,4 +139,5 @@
 				    :test #'divisible-by-p
 				    :key #'leading-monomial)))
+  ;; NOTE: Currently R is a list of terms of the remainder
   (cond
    (g					;division possible
@@ -162,12 +148,11 @@
       (let ((m (divide (leading-monomial p) (leading-monomial g))))
 	;; Multiply the equation c*f=sum ai*fi+r+p by cg.
-	(setf r (multiply-by r cg)
+	(setf r (mapc #'(lambda (trm) (multiply-by trm cg)) r)
 	      c (multiply-by c cg)
               ;; p := cg*p-cp*m*g
 	      p (grobner-op cp cg m p g))))
     (debug-cgb "/"))
-   (t							;no division possible
-    (setf r (add-to r (leading-term p))) ;move lt(p) to remainder
-    (setf p (subtract-from p (leading-term p))) ;move lt(p) to remainder
+   (t					   ;no division possible
+    (setf r (push (poly-remove-term p) r))	   ;move lt(p) to remainder
     (debug-cgb "+")))
   (values p r c division-count))
@@ -187,18 +172,17 @@
     ;; unit in the coefficient field.
     (return-from normal-form (values f nil 0)))
-  (do ((r (make-zero-for f))
+  (do ((r nil)
        (c (make-unit-for (leading-coefficient f)))
        (division-count 0))
       ((or (universal-zerop f)
 	   ;;(endp fl)
-	   (and top-reduction-only (not (universal-zerop r))))
+	   (and top-reduction-only (not (null r))))
        (progn
 	 (debug-cgb "~&~3T~D reduction~:P" division-count)
-	 (when (universal-zerop r)
+	 (when (null r)
 	   (debug-cgb " ---> 0")))
-       (setf f (add-to f r))
+       (setf (poly-termlist f) (nreconc r (poly-termlist f)))
        (values f c division-count))
-    (declare (fixnum division-count)
-	     (type poly r))
+    (declare (fixnum division-count))
     (multiple-value-setq (f r c division-count)
       (normal-form-step fl f r c division-count))))
