1 | ;;; -*- Mode: Lisp -*-
|
---|
2 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
3 | ;;;
|
---|
4 | ;;; Copyright (C) 1999, 2002, 2009, 2015 Marek Rychlik <rychlik@u.arizona.edu>
|
---|
5 | ;;;
|
---|
6 | ;;; This program is free software; you can redistribute it and/or modify
|
---|
7 | ;;; it under the terms of the GNU General Public License as published by
|
---|
8 | ;;; the Free Software Foundation; either version 2 of the License, or
|
---|
9 | ;;; (at your option) any later version.
|
---|
10 | ;;;
|
---|
11 | ;;; This program is distributed in the hope that it will be useful,
|
---|
12 | ;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
13 | ;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
14 | ;;; GNU General Public License for more details.
|
---|
15 | ;;;
|
---|
16 | ;;; You should have received a copy of the GNU General Public License
|
---|
17 | ;;; along with this program; if not, write to the Free Software
|
---|
18 | ;;; Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
|
---|
19 | ;;;
|
---|
20 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
21 |
|
---|
22 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
23 | ;;
|
---|
24 | ;; An implementation of the algorithm of Gebauer and Moeller, as
|
---|
25 | ;; described in the book of Becker-Weispfenning, p. 232
|
---|
26 | ;;
|
---|
27 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
28 |
|
---|
29 | (defpackage "GEBAUER-MOELLER"
|
---|
30 | (:use :cl :grobner-debug
|
---|
31 | :division :ring :monom :polynomial :order :ring-and-order
|
---|
32 | :pair-queue :priority-queue
|
---|
33 | )
|
---|
34 | (:export "GEBAUER-MOELLER"))
|
---|
35 |
|
---|
36 | (in-package :gebauer-moeller)
|
---|
37 |
|
---|
38 | (defun gebauer-moeller (ring-and-order f
|
---|
39 | &optional
|
---|
40 | (start 0)
|
---|
41 | (top-reduction-only $poly_top_reduction_only)
|
---|
42 | &aux
|
---|
43 | (ring (ro-ring ring-and-order)))
|
---|
44 | "Compute Grobner basis by using the algorithm of Gebauer and
|
---|
45 | Moeller. This algorithm is described as BUCHBERGERNEW2 in the book by
|
---|
46 | Becker-Weispfenning entitled ``Grobner Bases''. This function assumes
|
---|
47 | that all polynomials in F are non-zero."
|
---|
48 | (declare (ignore top-reduction-only)
|
---|
49 | (type fixnum start)
|
---|
50 | (type ring-and-order ring-and-order))
|
---|
51 | (cond
|
---|
52 | ((endp f) (return-from gebauer-moeller nil))
|
---|
53 | ((endp (cdr f))
|
---|
54 | (return-from gebauer-moeller (list (poly-primitive-part ring (car f))))))
|
---|
55 | (debug-cgb "~&GROBNER BASIS - GEBAUER MOELLER ALGORITHM")
|
---|
56 | (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
|
---|
57 | #+grobner-check (when (plusp start)
|
---|
58 | (grobner-test ring-and-order (subseq f 0 start) (subseq f 0 start)))
|
---|
59 | (let ((b (make-pair-queue))
|
---|
60 | (g (subseq f 0 start))
|
---|
61 | (f1 (subseq f start)))
|
---|
62 | (do () ((endp f1))
|
---|
63 | (multiple-value-setq (g b)
|
---|
64 | (gebauer-moeller-update g b (poly-primitive-part ring (pop f1)))))
|
---|
65 | (do () ((pair-queue-empty-p b))
|
---|
66 | (let* ((pair (pair-queue-remove b))
|
---|
67 | (g1 (pair-first pair))
|
---|
68 | (g2 (pair-second pair))
|
---|
69 | (h (normal-form ring-and-order (spoly ring-and-order g1 g2)
|
---|
70 | g
|
---|
71 | nil #| Always fully reduce! |#
|
---|
72 | )))
|
---|
73 | (unless (poly-zerop h)
|
---|
74 | (setf h (poly-primitive-part ring h))
|
---|
75 | (multiple-value-setq (g b)
|
---|
76 | (gebauer-moeller-update g b h))
|
---|
77 | (debug-cgb "~&Sugar: ~d Polynomials: ~d; Pairs left: ~d~%"
|
---|
78 | (pair-sugar pair) (length g) (pair-queue-size b))
|
---|
79 | )))
|
---|
80 | #+grobner-check(grobner-test ring-and-order g f)
|
---|
81 | (debug-cgb "~&GROBNER END")
|
---|
82 | g))
|
---|
83 |
|
---|
84 | (defun gebauer-moeller-update (g b h
|
---|
85 | &aux
|
---|
86 | c d e
|
---|
87 | (b-new (make-pair-queue))
|
---|
88 | g-new)
|
---|
89 | "An implementation of the auxillary UPDATE algorithm used by the
|
---|
90 | Gebauer-Moeller algorithm. G is a list of polynomials, B is a list of
|
---|
91 | critical pairs and H is a new polynomial which possibly will be added
|
---|
92 | to G. The naming conventions used are very close to the one used in
|
---|
93 | the book of Becker-Weispfenning."
|
---|
94 | (declare
|
---|
95 | #+allegro (dynamic-extent b)
|
---|
96 | (type poly h)
|
---|
97 | (type priority-queue b))
|
---|
98 | (setf c g d nil)
|
---|
99 | (do () ((endp c))
|
---|
100 | (let ((g1 (pop c)))
|
---|
101 | (declare (type poly g1))
|
---|
102 | (when (or (monom-rel-prime-p (poly-lm h) (poly-lm g1))
|
---|
103 | (and
|
---|
104 | (notany #'(lambda (g2) (monom-lcm-divides-monom-lcm-p
|
---|
105 | (poly-lm h) (poly-lm g2)
|
---|
106 | (poly-lm h) (poly-lm g1)))
|
---|
107 | c)
|
---|
108 | (notany #'(lambda (g2) (monom-lcm-divides-monom-lcm-p
|
---|
109 | (poly-lm h) (poly-lm g2)
|
---|
110 | (poly-lm h) (poly-lm g1)))
|
---|
111 | d)))
|
---|
112 | (push g1 d))))
|
---|
113 | (setf e nil)
|
---|
114 | (do () ((endp d))
|
---|
115 | (let ((g1 (pop d)))
|
---|
116 | (declare (type poly g1))
|
---|
117 | (unless (monom-rel-prime-p (poly-lm h) (poly-lm g1))
|
---|
118 | (push g1 e))))
|
---|
119 | (do () ((pair-queue-empty-p b))
|
---|
120 | (let* ((pair (pair-queue-remove b))
|
---|
121 | (g1 (pair-first pair))
|
---|
122 | (g2 (pair-second pair)))
|
---|
123 | (declare (type pair pair)
|
---|
124 | (type poly g1 g2))
|
---|
125 | (when (or (not (monom-divides-monom-lcm-p
|
---|
126 | (poly-lm h)
|
---|
127 | (poly-lm g1) (poly-lm g2)))
|
---|
128 | (monom-lcm-equal-monom-lcm-p
|
---|
129 | (poly-lm g1) (poly-lm h)
|
---|
130 | (poly-lm g1) (poly-lm g2))
|
---|
131 | (monom-lcm-equal-monom-lcm-p
|
---|
132 | (poly-lm h) (poly-lm g2)
|
---|
133 | (poly-lm g1) (poly-lm g2)))
|
---|
134 | (pair-queue-insert b-new (make-pair g1 g2)))))
|
---|
135 | (dolist (g3 e)
|
---|
136 | (pair-queue-insert b-new (make-pair h g3)))
|
---|
137 | (setf g-new nil)
|
---|
138 | (do () ((endp g))
|
---|
139 | (let ((g1 (pop g)))
|
---|
140 | (declare (type poly g1))
|
---|
141 | (unless (monom-divides-p (poly-lm h) (poly-lm g1))
|
---|
142 | (push g1 g-new))))
|
---|
143 | (push h g-new)
|
---|
144 | (values g-new b-new))
|
---|