source: branches/f4grobner/buchberger.lisp@ 4197

;;; -*-  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 
        :polynomial :division 
        :criterion :pair-queue :priority-queue
        )
  (:export "BUCHBERGER" "PARALLEL-BUCHBERGER")
  (:documentation "Buchberger Algorithm Implementation."))

(in-package :buchberger)


(defun buchberger (f 
                   &optional 
                     (start 0)
                     (top-reduction-only $poly_top_reduction_only))
  "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 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 (subseq f 0 start) (subseq f 0 start)))
  ;;Initialize critical pairs
  (let ((b (make-critical-pair-queue *normal-strategy* f start))
        (b-done (make-hash-table :test #'equal)))
    (declare (type critical-pair-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 ()
        ((queue-empty-p b)
         #+grobner-check(grobner-test f f)
         (debug-cgb "~&GROBNER END")
         f)
      (let ((pair (dequeue b)))
        (declare (type critical-pair pair))
        (cond
          ((criterion-1 pair) nil)
          ((criterion-2 pair b-done f) nil)
          (t 
           (let ((sp (normal-form (s-polynomial
                                   (critical-pair-first pair) 
                                   (critical-pair-second pair))
                                  f top-reduction-only)))
             (declare (type poly sp))
             (cond
               ((universal-zerop sp)
                nil)
               (t
                (setf sp (poly-primitive-part sp)
                      f (nconc f (list sp)))
                ;; Add new critical pairs
                (dolist (h f)
                  (enqueue b (make-instance 'critical-pair :first h :second sp)))
                (debug-cgb "~&Polynomials: ~d; Pairs left: ~d; Pairs done: ~d;"
                           (length f) (queue-size b)
                           (hash-table-count b-done)))))))
        (setf (gethash (list (critical-pair-first pair) (critical-pair-second pair)) b-done)
              t)))))

(defun parallel-buchberger (f 
                            &optional 
                              (start 0)
                              (top-reduction-only $poly_top_reduction_only))
  "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 (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 (subseq f 0 start) (subseq f 0 start)))
  ;;Initialize critical pairs
  (let ((b (make-critical-pair-queue *normal-strategy* f start))
        (b-done (make-hash-table :test #'equal))
        (coeff-zero (make-zero-for (leading-coefficient (car f))))
        (coeff-unit (make-unit-for (leading-coefficient (car f)))))
    (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 ()
        ((queue-empty-p b)
         #+grobner-check(grobner-test f f)
         (debug-cgb "~&GROBNER END")
         f)
      (let ((pair (dequeue b)))
        (when (null (critical-pair-data pair))
          (setf (critical-pair-data pair) (list (s-polynomial
                                                          (critical-pair-first pair)
                                                          (critical-pair-second pair))
                                                         coeff-zero
                                                         coeff-unit
                                                         0)))
        (cond
          ((criterion-1 pair) nil)
          ((criterion-2 pair b-done f) nil)
          (t
           (let* ((dd (critical-pair-data pair))
                  (p (first dd))
                  (sp (second dd))
                  (c (third dd))
                  (division-count (fourth dd)))
             (cond
               ((universal-zerop p)     ;normal form completed
                (debug-cgb "~&~3T~d reduction~:p" division-count)
                (cond 
                  ((universal-zerop sp)
                   (debug-cgb " ---> 0")
                   nil)
                  (t
                   (setf sp (poly-primitive-part sp)
                         f (nconc f (list sp)))
                   ;; Add new critical pairs
                   (dolist (h f)
                     (enqueue b (make-instance 'critical-pair :first h :second sp)))
                   (debug-cgb "~&Polynomials: ~d; Pairs left: ~d; Pairs done: ~d;"
                              (length f) (queue-size b)
                              (hash-table-count b-done))))
                (setf (gethash (list (critical-pair-first pair) (critical-pair-second pair))
                               b-done) t))
               (t                               ;normal form not complete
                (do ()
                    ((cond
                       ((universal-zerop p)
                        (debug-cgb ".")
                        t)
                       (t nil))
                     (setf (first dd) p
                           (second dd) sp
                           (third dd) c
                           (fourth dd) division-count)
                     (enqueue b pair))
                  (multiple-value-setq (p sp c division-count)
                    (normal-form-step f p sp c division-count))))))))))))