Changeset 4049 for branches/f4grobner

Timestamp:

2016-05-31T17:29:55-07:00 (9 years ago)

Author:

Marek Rychlik

Message:

* empty log message *

File:

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

branches/f4grobner/division.lisp

@@ -21 +21 @@
 
 (defpackage "DIVISION"
-  (:use :cl :utils :ring :monom :polynomial :grobner-debug :term :ring-and-order)
+  (:use :cl :utils :monom :polynomial :grobner-debug)
   (:export "$POLY_TOP_REDUCTION_ONLY"
            "POLY-PSEUDO-DIVIDE"
@@ -55 +55 @@
                    (multiply c2 (leading-coefficient f))
                    (multiply c1 (leading-coefficient g))))
-  #+grobner-check(monom-equal-p (poly-lm f) (monom-mul m (poly-lm g)))
+  #+grobner-check(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.  
@@ -61 +61 @@
   ;; TODO: Make sure that the sugar calculation is correct if leading
   ;; terms are dropped.
-  (poly-sub ring-and-order
-            (scalar-times-poly-1 ring c2 f)
-            (scalar-times-poly-1 ring c1 (monom-times-poly m g))))
-
-(defun check-loop-invariant (ring-and-order c f a fl r p
+  (subtract
+   (multiply c2 f)
+   (multiply c1 (multiply m g))))
+
+(defun check-loop-invariant (c f a fl r p
                              &aux 
-                               (ring (ro-ring ring-and-order))
                                (p-zero (make-poly-zero))
                                (a (mapcar #'poly-reverse a))
@@ -84 +83 @@
           c f a fl r p)
   |#
-  (flet ((p-add (x y) (poly-add ring-and-order x y))
-         (p-sub (x y) (poly-sub ring-and-order x y))
-         (p-mul (x y) (poly-mul ring-and-order x y)))
-    (let* ((prod (inner-product a fl p-add p-mul p-zero))
-           (succeeded-p 
-            (poly-zerop
-             (p-sub
-              (scalar-times-poly ring c f)
-              (reduce #'p-add (list prod r p))))))
-      (unless succeeded-p 
-          (error "#### Polynomial division Loop invariant failed ####:~%C=~A~%F=~A~%A=~A~%FL=~A~%R=~A~%P=~A~%" 
-                 c f a fl r p))
-      succeeded-p)))
+  (let* ((prod (inner-product a fl #'add #'multiply 0))
+         (succeeded-p 
+          (universal-zerop
+           (subtract
+            (multiply c f)
+            (reduce #'add (list prod r p))))))
+    (unless succeeded-p 
+      (error "#### Polynomial division Loop invariant failed ####:~%C=~A~%F=~A~%A=~A~%FL=~A~%R=~A~%P=~A~%" 
+             c f a fl r p))
+    succeeded-p))
   
 
-(defun poly-pseudo-divide (ring-and-order f fl
-                           &aux
-                             (ring (ro-ring ring-and-order)))
+(defun poly-pseudo-divide (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
@@ -112 +106 @@
   (declare (type poly f) (list fl))
   ;; Loop invariant: c*f=sum ai*fi+r+p, where p must eventually become 0
-  (do ((r (make-poly-zero))
-       (c (funcall (ring-unit ring)))
-       (a (make-list (length fl) :initial-element (make-poly-zero)))
+  (do ((r 0)
+       (c 1)
+       (a (make-list (length fl) :initial-element 0))
        (division-count 0)
        (p f))
-      ((poly-zerop p)
-       #+grobner-check(check-loop-invariant ring-and-order c f a fl r p)
+      ((universal-zerop p)
+       #+grobner-check(check-loop-invariant c f a fl r p)
        (debug-cgb "~&~3T~d reduction~:p" division-count)
-       (when (poly-zerop r) (debug-cgb " ---> 0"))
+       (when (universal-zerop r) (debug-cgb " ---> 0"))
        ;; We obtained the terms in reverse order, so must fix that
        (setf a (mapcar #'poly-nreverse a)
              r (poly-nreverse r))
        ;; Initialize the sugar of the quotients
-       (mapc #'poly-reset-sugar a)
+       ;; (mapc #'poly-reset-sugar a) ;; TODO: Sugar is currently unimplemented
        (values a r c division-count))
     (declare (fixnum division-count))
     ;; Check the loop invariant here
-    #+grobner-check(check-loop-invariant ring-and-order c f a fl r p)
+    #+grobner-check(check-loop-invariant c f a fl r p)
     (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))))
+            (push (leading-term 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
+           ((monom-divides-p (leading-monomial (car fl)) (leading-monomial p)) ;division occurred
             (incf division-count)
             (multiple-value-bind (gcd c1 c2)
-                (funcall (ring-ezgcd ring) (poly-lc (car fl)) (poly-lc p))
+                (universal-ezgcd (leading-coefficient (car fl)) (leading-coefficient p))
               (declare (ignore gcd))
-              (let ((m (monom-div (poly-lm p) (poly-lm (car fl)))))
+              (let ((m (divide (leading-monomial p) (leading-monomial (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))))
+                          (setf (car x) (multiply c1 (car x))))
                       a)
-                (setf r (scalar-times-poly ring c1 r)
-                      c (funcall (ring-mul ring) c c1)
-                      p (grobner-op ring-and-order c2 c1 m p (car fl)))
+                (setf r (multiply c1 r)
+                      c (multiply c c1)
+                      p (grobner-op c2 c1 m p (car fl)))
                 (push (make-term :monom m :coeff c2) (poly-termlist (car b))))
               t))))
       )))
 
-(defun poly-exact-divide (ring-and-order f g)
+(defun poly-exact-divide (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) (type ring-and-order ring-and-order))
+  (declare (type poly f g))
   (multiple-value-bind (quot rem coeff division-count)
-      (poly-pseudo-divide ring-and-order f (list g))
+      (poly-pseudo-divide 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."))
+    (unless (universal-zerop rem) (error "Exact division failed."))
     (car quot)))
 
@@ -174 +168 @@
 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
-(defun normal-form-step (ring-and-order fl p r c division-count
+(defun normal-form-step (fl p r c division-count
                          &aux 
-                           (ring (ro-ring ring-and-order))
-                           (g (find (poly-lm p) fl
+                           (g (find (leading-monomial p) fl
                                     :test #'monom-divisible-by-p
-                                    :key #'poly-lm)))
+                                    :key #'leading-monomial)))
   (cond
    (g                                   ;division possible
     (incf division-count)
     (multiple-value-bind (gcd cg cp)
-        (funcall (ring-ezgcd ring) (poly-lc g) (poly-lc p))
+        (universal-ezgcd (leading-coefficient g) (leading-coefficient p))
       (declare (ignore gcd))
-      (let ((m (monom-div (poly-lm p) (poly-lm g))))
+      (let ((m (divide (leading-monomial p) (leading-monomial 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)
+        (setf r (multiply cg r)
+              c (multiply c cg)
               ;; p := cg*p-cp*m*g
-              p (grobner-op ring-and-order cp cg m p g))))
+              p (grobner-op 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))))
+    (push (leading-term 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 "+")))
@@ -207 +200 @@
 ;; method would be to test NORMAL-FORM using POLY-PSEUDO-DIVIDE.
 ;;
-(defun normal-form (ring-and-order f fl
-                    &optional 
-                      (top-reduction-only $poly_top_reduction_only)
-                      (ring (ro-ring ring-and-order)))
+(defun normal-form (f fl
+                    &optional
+                      (top-reduction-only $poly_top_reduction_only))
   #+grobner-check(when (null fl) (warn "normal-form: empty divisor list."))
-  (do ((r (make-poly-zero))
-       (c (funcall (ring-unit ring)))
+  (do ((r 0)
+       (c 1)
        (division-count 0))
-      ((or (poly-zerop f)
+      ((or (universal-zerop f)
            ;;(endp fl)
-           (and top-reduction-only (not (poly-zerop r))))
+           (and top-reduction-only (not (universal-zerop r))))
        (progn
          (debug-cgb "~&~3T~D reduction~:P" division-count)
-         (when (poly-zerop r)
+         (when (universal-zerop r)
            (debug-cgb " ---> 0")))
        (setf (poly-termlist f) (nreconc (poly-termlist r) (poly-termlist f)))
@@ -227 +219 @@
              (type poly r))
     (multiple-value-setq (f r c division-count)
-      (normal-form-step ring-and-order fl f r c division-count))))
+      (normal-form-step fl f r c division-count))))
 
 (defun buchberger-criterion (ring-and-order g)