source: branches/f4grobner/division.lisp@ 63

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; 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))))

(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))))))