source: branches/f4grobner/division.lisp@ 4100

;;; -*-  Mode: Lisp -*- 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;                                                                              
;;;  Copyright (C) 1999, 2002, 2009, 2015 Marek Rychlik <rychlik@u.arizona.edu>          
;;;                                                                              
;;;  This program is free software; you can redistribute it and/or modify        
;;;  it under the terms of the GNU General Public License as published by        
;;;  the Free Software Foundation; either version 2 of the License, or           
;;;  (at your option) any later version.                                         
;;;                                                                              
;;;  This program is distributed in the hope that it will be useful,             
;;;  but WITHOUT ANY WARRANTY; without even the implied warranty of              
;;;  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               
;;;  GNU General Public License for more details.                                
;;;                                                                              
;;;  You should have received a copy of the GNU General Public License           
;;;  along with this program; if not, write to the Free Software                 
;;;  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  
;;;                                                                              
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defpackage "DIVISION"
  (:use :cl :copy :utils :monom :polynomial :grobner-debug :symbolic-polynomial)
  (:export "$POLY_TOP_REDUCTION_ONLY"
           "POLY-PSEUDO-DIVIDE"
           "POLY-EXACT-DIVIDE"
           "NORMAL-FORM-STEP"
           "NORMAL-FORM"
           "POLY-NORMALIZE"
           "POLY-NORMALIZE-LIST"
           "BUCHBERGER-CRITERION"
           "GROBNER-TEST"
           )
  (:documentation
   "An implementation of the division algorithm in the polynomial ring."))

(in-package :division)

(defvar $poly_top_reduction_only nil
  "If not FALSE, use top reduction only whenever possible.
Top reduction means that division algorithm stops after the first reduction.")


(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))
  #+grobner-check(universal-zerop
                  (subtract
                   (multiply c2 (leading-coefficient f))
                   (multiply c1 (leading-coefficient 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.  
  ;;
  ;; TODO: Make sure that the sugar calculation is correct if leading
  ;; terms are dropped.
  (subtract
   (multiply f c2)
   (multiply m g c1)))

(defun check-loop-invariant (c f a fl r p
                             &aux 
                               (p-zero (make-zero-for f))
                               (a (mapcar #'poly-reverse a))
                               (r (poly-reverse r)))
  "Check loop invariant of division algorithms, when we divide a
polynomial F by the list of polynomials FL. The invariant is the
identity C*F=SUM AI*FI+R+P, where F0 is the initial value of F, A is
the list of partial quotients, R is the intermediate value of the
remainder, and P is the intermediate value which eventually becomes
0. A thing to remember is that the terms of polynomials in A and
the polynomial R have their terms in reversed order. Hence, before
the arithmetic is performed, we need to fix the order of terms"
  #|
  (format t "~&----------------------------------------------------------------~%")
  (format t "#### Loop invariant check ####:~%C=~A~%F=~A~%A=~A~%FL=~A~%R=~A~%P=~A~%" 
          c f a fl r p)
  |#
  (let* ((prod (inner-product a fl add multiply p-zero))
         (succeeded-p (universal-zerop (subtract (multiply f c) (add 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 (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. The sugar of 
the quotients is initialized to default."
  (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))
       (c 1)
       (a (make-list (length fl) :initial-element (make-zero-for f)))
       (division-count 0)
       (p f))
      ((universal-zerop p)
       #+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"))
       ;; We obtained the terms in reverse order, so must fix that
       (setf a (mapcar #'poly-reverse a)
             r (poly-reverse r))
       (values a r c division-count))
    (declare (fixnum division-count))
    ;; Check the loop invariant here
    #+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
            t)
           ((divides-p (leading-monomial (car fl)) (leading-monomial p)) ;division occurred
            (incf division-count)
            (multiple-value-bind (gcd c1 c2)
                (universal-ezgcd (leading-coefficient (car fl)) (leading-coefficient p))
              (declare (ignore gcd))
              (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) (multiply (car x) c1)))
                      a)
                (setf r (multiply r c1)
                      c (multiply c c1)
                      p (grobner-op c2 c1 m p (car fl)))
                (setf (car b) (add (car b)
                                   (change-class m 'term :coeff c2))))
              t))))
      )))

(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))
  (multiple-value-bind (quot rem coeff division-count)
      (poly-pseudo-divide f (list g))
    (declare (ignore division-count coeff)
             (list quot)
             (type poly rem)
             (type fixnum division-count))
    (unless (universal-zerop rem) (error "Exact division failed."))
    (car quot)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; An implementation of the normal form
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defun normal-form-step (fl p r c division-count
                         &aux 
                           (g (find (leading-monomial p) fl
                                    :test #'divisible-by-p
                                    :key #'leading-monomial)))
  (cond
   (g                                   ;division possible
    (incf division-count)
    (multiple-value-bind (gcd cg cp)
        (universal-ezgcd (leading-coefficient g) (leading-coefficient p))
      (declare (ignore gcd))
      (let ((m (divide (leading-monomial p) (leading-monomial g))))
        ;; Multiply the equation c*f=sum ai*fi+r+p by cg.
        (setf r (multiply r cg)
              c (multiply c cg)
              ;; p := cg*p-cp*m*g
              p (grobner-op cp cg m p g))))
    (debug-cgb "/"))
   (t                                                   ;no division possible
    (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 "+")))
  (values p r c division-count))

;;
;; Merge NORMAL-FORM someday with POLY-PSEUDO-DIVIDE.
;;
;; TODO: It is hard to test normal form as there is no loop invariant,
;; like for POLY-PSEUDO-DIVIDE.  Is there a testing strategy? One
;; method would be to test NORMAL-FORM using POLY-PSEUDO-DIVIDE.
;;
(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-zero-for f))
       (c 1)
       (division-count 0))
      ((or (universal-zerop f)
           ;;(endp fl)
           (and top-reduction-only (not (universal-zerop r))))
       (progn
         (debug-cgb "~&~3T~D reduction~:P" division-count)
         (when (universal-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 fl f r c division-count))))

(defun spoly (f g)
  "It yields the S-polynomial of polynomials F and G."
  (declare (type poly f g))
  (let* ((lcm (universal-lcm (leading-monomial f) (leading-monomial g)))
          (mf (divide lcm (leading-monomial f)))
          (mg (divide lcm (leading-monomial g))))
    (declare (type monom mf mg))
    (multiple-value-bind (c cf cg)
        (universal-ezgcd (leading-coefficient f) (leading-coefficient g))
      (declare (ignore c))
      (subtract
       (multiply (multiply mf f) cg)
       (multiply (multiply mg g) cf)))))

(defun buchberger-criterion (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 #'universal-zerop
         (makelist (normal-form (spoly (elt g i) (elt g j)) g nil)
                   (i 0 (- (length g) 2))
                   (j (1+ i) (1- (length g))))))


(defun poly-normalize (p &aux (c (leading-coefficient p)))
  "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."
  (mapc #'(lambda (term)
            (setf (term-coeff term) (divide (term-coeff term) c)))
        (poly-termlist p))
  p)

(defun poly-normalize-list (plist)
  "Divide every polynomial in a list PLIST by its leading coefficient. "
  (mapcar #'(lambda (x) (poly-normalize x)) plist))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; The function GROBNER-TEST is provided primarily for debugging purposes. To
;; enable verification of grobner bases with BUCHBERGER-CRITERION, do
;; (pushnew :grobner-check *features*) and compile/load this file.
;; With this feature, the calculations will slow down CONSIDERABLY.
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defun grobner-test (g f)
  "Test whether G is a Grobner basis and F is contained in G. Return T
upon success and NIL otherwise."
  (debug-cgb "~&GROBNER CHECK: ")
  (let (($poly_grobner_debug nil)
        (stat1 (buchberger-criterion g))
        (stat2
          (every #'universal-zerop
                 (makelist (normal-form (copy-instance (elt f i)) g nil)
                           (i 0 (1- (length f)))))))
    (unless stat1 (error "~&Buchberger criterion failed, not a grobner basis: ~A" g))
    (unless stat2
      (error "~&Original polynomials not in ideal spanned by Grobner basis: ~A" f)))
  (debug-cgb "~&GROBNER CHECK END")
  t)