source: branches/f4grobner/buchberger.lisp@ 77

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


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Buchberger Algorithm Implementation
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defun buchberger (ring f start &optional (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 ring (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 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 (spoly ring (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 f start &optional (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 ring (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 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
                                                       (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 f p sp c division-count))))))))))))