Changeset 4463 for branches/f4grobner/division.lisp

Timestamp:

2016-06-14T19:35:34-07:00 (8 years ago)

Author:

Marek Rychlik

Message:

 

File:

• branches/f4grobner/division.lisp (modified) (9 diffs)

branches/f4grobner/division.lisp

@@ -43 +43 @@
 
 (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 +60 @@
   |#
   (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 +78 @@
   (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 +86 @@
        #+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 +94 @@
          (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 +107 @@
                           (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 +139 @@
                                     :test #'divisible-by-p
                                     :key #'leading-monomial)))
+  ;; NOTE: Currently R is a list of terms of the remainder
   (cond
    (g                                   ;division possible
@@ -162 +148 @@
       (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 +172 @@
     ;; 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))))