source: branches/f4grobner/buchberger.lisp@ 3560

;;; -*-  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 "BUCHBERGER"

  (:use :cl :grobner-debug 

        :ring :polynomial :order :ring-and-order :division 

        :criterion :pair-queue :priority-queue

        )

  (:export "BUCHBERGER" "PARALLEL-BUCHBERGER"))



(in-package :buchberger)



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

;;

;; Buchberger Algorithm Implementation

;;

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



(defun buchberger (ring-and-order f 

                   &optional 

                     (start 0)

                     (top-reduction-only $poly_top_reduction_only)

                   &aux

                     (ring (ro-ring ring-and-order)))

  "An implementation of the Buchberger algorithm. Return Grobner basis

of the ideal generated by the polynomial list F.  Polynomials 0 to

START-1 are assumed to be a Grobner basis already, so that certain

critical pairs will not be examined. If TOP-REDUCTION-ONLY set, top

reduction will be preformed. This function assumes that all polynomials

in F are non-zero."

  (declare (type ring-and-order ring-and-order) (type fixnum start))

  (when (endp f) (return-from buchberger f)) ;cut startup costs

  (debug-cgb "~&GROBNER BASIS - BUCHBERGER ALGORITHM")

  (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))

  #+grobner-check  (when (plusp start)

                     (grobner-test ring-and-order (subseq f 0 start) (subseq f 0 start)))

  ;;Initialize critical pairs

  (let ((b (pair-queue-initialize (make-pair-queue)

                                  f start))

        (b-done (make-hash-table :test #'equal)))

    (declare (type priority-queue b) (type hash-table b-done))

    (dotimes (i (1- start))

      (do ((j (1+ i) (1+ j))) ((>= j start))

        (setf (gethash (list (elt f i) (elt f j)) b-done) t)))

    (do ()

        ((pair-queue-empty-p b)

         #+grobner-check(grobner-test ring-and-order f f)

         (debug-cgb "~&GROBNER END")

         f)

      (let ((pair (pair-queue-remove b)))

        (declare (type pair pair))

        (cond

          ((criterion-1 pair) nil)

          ((criterion-2 pair b-done f) nil)

          (t 

           (let ((sp (normal-form ring-and-order 

                                  (spoly ring-and-order

                                         (pair-first pair) 

                                         (pair-second pair))

                                  f top-reduction-only)))

             (declare (type poly sp))

             (cond

               ((poly-zerop sp)

                nil)

               (t

                (setf sp (poly-primitive-part ring sp)

                      f (nconc f (list sp)))

                ;; Add new critical pairs

                (dolist (h f)

                  (pair-queue-insert b (make-pair h sp)))

                (debug-cgb "~&Sugar: ~d Polynomials: ~d; Pairs left: ~d; Pairs done: ~d;"

                           (pair-sugar pair) (length f) (pair-queue-size b)

                           (hash-table-count b-done)))))))

        (setf (gethash (list (pair-first pair) (pair-second pair)) b-done)

              t)))))



(defun parallel-buchberger (ring-and-order f 

                            &optional 

                              (start 0)

                              (top-reduction-only $poly_top_reduction_only)

                            &aux

                              (ring (ro-ring ring-and-order)))

  "An implementation of the Buchberger algorithm. Return Grobner basis

of the ideal generated by the polynomial list F.  Polynomials 0 to

START-1 are assumed to be a Grobner basis already, so that certain

critical pairs will not be examined. If TOP-REDUCTION-ONLY set, top

reduction will be preformed."

  (declare (type ring-and-order ring-and-order)

           (ignore top-reduction-only)

           (type fixnum start))

  (when (endp f) (return-from parallel-buchberger f)) ;cut startup costs

  (debug-cgb "~&GROBNER BASIS - PARALLEL-BUCHBERGER ALGORITHM")

  (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))

  #+grobner-check  (when (plusp start)

                     (grobner-test ring-and-order (subseq f 0 start) (subseq f 0 start)))

  ;;Initialize critical pairs

  (let ((b (pair-queue-initialize (make-pair-queue) f start))

        (b-done (make-hash-table :test #'equal)))

    (declare (type priority-queue b)

             (type hash-table b-done))

    (dotimes (i (1- start))

      (do ((j (1+ i) (1+ j))) ((>= j start))

        (declare (type fixnum j))

        (setf (gethash (list (elt f i) (elt f j)) b-done) t)))

    (do ()

        ((pair-queue-empty-p b)

         #+grobner-check(grobner-test ring-and-order f f)

         (debug-cgb "~&GROBNER END")

         f)

      (let ((pair (pair-queue-remove b)))

        (when (null (pair-division-data pair))

          (setf (pair-division-data pair) (list (spoly ring-and-order

                                                       (pair-first pair)

                                                       (pair-second pair))

                                                (make-poly-zero)

                                                (funcall (ring-unit ring))

                                                0)))

        (cond

          ((criterion-1 pair) nil)

          ((criterion-2 pair b-done f) nil)

          (t

           (let* ((dd (pair-division-data pair))

                  (p (first dd))

                  (sp (second dd))

                  (c (third dd))

                  (division-count (fourth dd)))

             (cond

               ((poly-zerop p)          ;normal form completed

                (debug-cgb "~&~3T~d reduction~:p" division-count)

                (cond 

                  ((poly-zerop sp)

                   (debug-cgb " ---> 0")

                   nil)

                  (t

                   (setf sp (poly-nreverse sp)

                         sp (poly-primitive-part ring sp)

                         f (nconc f (list sp)))

                   ;; Add new critical pairs

                   (dolist (h f)

                     (pair-queue-insert b (make-pair h sp)))

                   (debug-cgb "~&Sugar: ~d Polynomials: ~d; Pairs left: ~d; Pairs done: ~d;"

                              (pair-sugar pair) (length f) (pair-queue-size b)

                              (hash-table-count b-done))))

                (setf (gethash (list (pair-first pair) (pair-second pair))

                               b-done) t))

               (t                               ;normal form not complete

                (do ()

                    ((cond

                       ((> (poly-sugar sp) (pair-sugar pair))

                        (debug-cgb "(~a)?" (poly-sugar sp))

                        t)

                       ((poly-zerop p)

                        (debug-cgb ".")

                        t)

                       (t nil))

                     (setf (first dd) p

                           (second dd) sp

                           (third dd) c

                           (fourth dd) division-count

                           (pair-sugar pair) (poly-sugar sp))

                     (pair-queue-insert b pair))

                  (multiple-value-setq (p sp c division-count)

                    (normal-form-step ring-and-order f p sp c division-count))))))))))))