source: branches/f4grobner/division.lisp@ 461

;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;                                                                              
;;;  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 :ring :monomial :polynomial)
  (:export "POLY-PSEUDO-DIVIDE"
           "POLY-EXACT-DIVIDE"
           "NORMAL-FORM"
           "POLY-NORMALIZE"
           ))

(in-package :division)

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


(defun poly-normalize (ring p &aux (c (poly-lc 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 of the
COEFFICIENT-RING package."
  (mapc #'(lambda (term)
            (setf (term-coeff term) (funcall (ring-div ring) (term-coeff term) c)))
        (poly-termlist p))
  p)

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