source: branches/f4grobner/gebauer-moeller.lisp@ 1181

;;; -*-  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.  
;;;                                                                              
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; An implementation of the algorithm of Gebauer and Moeller, as
;; described in the book of Becker-Weispfenning, p. 232
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defpackage "GEBAUER-MOELLER"
  (:use :cl :grobner-debug 
        :division :ring :monomial :polynomial :order 
        :pair-queue :priority-queue
        )
  (:export "GEBAUER-MOELLER"))

(in-package :gebauer-moeller)

(defun gebauer-moeller (ring f start &optional (top-reduction-only $poly_top_reduction_only))
  "Compute Grobner basis by using the algorithm of Gebauer and
Moeller.  This algorithm is described as BUCHBERGERNEW2 in the book by
Becker-Weispfenning entitled ``Grobner Bases''. This function assumes
that all polynomials in F are non-zero."
  (declare (ignore top-reduction-only)
           (type fixnum start))
  (cond
   ((endp f) (return-from gebauer-moeller nil))
   ((endp (cdr f))
    (return-from gebauer-moeller (list (poly-primitive-part ring (car f))))))
   (debug-cgb "~&GROBNER BASIS - GEBAUER MOELLER ALGORITHM")
   (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
  #+grobner-check  (when (plusp start)
                     (grobner-test ring (subseq f 0 start) (subseq f 0 start)))
  (let ((b (make-pair-queue))
        (g (subseq f 0 start))
        (f1 (subseq f start)))
    (do () ((endp f1))
      (multiple-value-setq (g b)
        (gebauer-moeller-update g b (poly-primitive-part ring (pop f1)))))
    (do () ((pair-queue-empty-p b))
      (let* ((pair (pair-queue-remove b))
             (g1 (pair-first pair))
             (g2 (pair-second pair))
             (h (normal-form ring (spoly ring g1 g2)
                             g
                             nil #| Always fully reduce! |#
                             )))
        (unless (poly-zerop h)
          (setf h (poly-primitive-part ring h))
          (multiple-value-setq (g b)
            (gebauer-moeller-update g b h))
          (debug-cgb "~&Sugar: ~d Polynomials: ~d; Pairs left: ~d~%"
                     (pair-sugar pair) (length g) (pair-queue-size b))
          )))
    #+grobner-check(grobner-test ring g f)
    (debug-cgb "~&GROBNER END")
    g))

(defun gebauer-moeller-update (g b h
                 &aux
                 c d e
                 (b-new (make-pair-queue))
                 g-new)
  "An implementation of the auxillary UPDATE algorithm used by the
Gebauer-Moeller algorithm. G is a list of polynomials, B is a list of
critical pairs and H is a new polynomial which possibly will be added
to G. The naming conventions used are very close to the one used in
the book of Becker-Weispfenning."
  (declare
   #+allegro (dynamic-extent b)
   (type poly h)
   (type priority-queue b))
  (setf c g d nil) 
  (do () ((endp c))
    (let ((g1 (pop c)))
      (declare (type poly g1))
      (when (or (monom-rel-prime-p (poly-lm h) (poly-lm g1))
                (and
                 (notany #'(lambda (g2) (monom-lcm-divides-monom-lcm-p
                                         (poly-lm h) (poly-lm g2)
                                         (poly-lm h) (poly-lm g1)))
                         c)
                 (notany #'(lambda (g2) (monom-lcm-divides-monom-lcm-p
                                         (poly-lm h) (poly-lm g2)
                                         (poly-lm h) (poly-lm g1)))
                         d)))
        (push g1 d))))
  (setf e nil)
  (do () ((endp d))
    (let ((g1 (pop d)))
      (declare (type poly g1))
      (unless (monom-rel-prime-p (poly-lm h) (poly-lm g1))
        (push g1 e))))
  (do () ((pair-queue-empty-p b))
    (let* ((pair (pair-queue-remove b))
           (g1 (pair-first pair))
           (g2 (pair-second pair)))
      (declare (type pair pair)
               (type poly g1 g2))
      (when (or (not (monom-divides-monom-lcm-p
                      (poly-lm h)
                      (poly-lm g1) (poly-lm g2)))
                (monom-lcm-equal-monom-lcm-p
                 (poly-lm g1) (poly-lm h)
                 (poly-lm g1) (poly-lm g2))
                (monom-lcm-equal-monom-lcm-p
                 (poly-lm h) (poly-lm g2)
                 (poly-lm g1) (poly-lm g2)))
        (pair-queue-insert b-new (make-pair g1 g2)))))
  (dolist (g3 e)
    (pair-queue-insert b-new (make-pair h g3)))
  (setf g-new nil)
  (do () ((endp g))
    (let ((g1 (pop g)))
      (declare (type poly g1))
      (unless (monom-divides-p (poly-lm h) (poly-lm g1))
        (push g1 g-new))))
  (push h g-new)
  (values g-new b-new))