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

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

;;;                                                                              

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;



(in-package :ngrobner)



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;

;; An implementation of the algorithm of Gebauer and Moeller, as

;; described in the book of Becker-Weispfenning, p. 232

;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;



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