close Warning: Can't synchronize with repository "(default)" (The repository directory has changed, you should resynchronize the repository with: trac-admin $ENV repository resync '(default)'). Look in the Trac log for more information.

source: branches/f4grobner/division.lisp@ 4528

Last change on this file since 4528 was 4497, checked in by Marek Rychlik, 8 years ago

* empty log message *

File size: 9.5 KB
Line 
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(defpackage "DIVISION"
23 (:use :cl :copy :utils :monom :polynomial :grobner-debug :symbolic-polynomial)
24 (:export "$POLY_TOP_REDUCTION_ONLY"
25 "POLY-PSEUDO-DIVIDE"
26 "POLY-EXACT-DIVIDE"
27 "NORMAL-FORM-STEP"
28 "NORMAL-FORM"
29 "POLY-NORMALIZE"
30 "POLY-NORMALIZE-LIST"
31 "BUCHBERGER-CRITERION"
32 "GROBNER-TEST"
33 )
34 (:documentation
35 "An implementation of the division algorithm in the polynomial ring."))
36
37(in-package :division)
38
39(proclaim '(optimize (speed 3) (space 0) (safety 0) (debug 0)))
40
41(defvar $poly_top_reduction_only nil
42 "If not FALSE, use top reduction only whenever possible.
43Top reduction means that division algorithm stops after the first reduction.")
44
45(defmacro grobner-op (c1 c2 m f g)
46 "Returns C2*F-C1*M*G, where F and G are polynomials M is a monomial."
47 `(subtract (multiply ,f ,c2) (multiply ,g ,m ,c1)))
48
49(defun check-loop-invariant (c f a fl r p &aux (p-zero (make-zero-for f)))
50 "Check loop invariant of division algorithms, when we divide a
51polynomial F by the list of polynomials FL. The invariant is the
52identity C*F=SUM AI*FI+R+P, where F0 is the initial value of F, A is
53the list of partial quotients, R is the intermediate value of the
54remainder, and P is the intermediate value which eventually becomes
550."
56 #|
57 (format t "~&----------------------------------------------------------------~%")
58 (format t "#### Loop invariant check ####:~%C=~A~%F=~A~%A=~A~%FL=~A~%R=~A~%P=~A~%"
59 c f a fl r p)
60 |#
61 (let* ((prod (inner-product a fl add multiply p-zero))
62 (succeeded-p (universal-zerop (subtract (multiply f c) (add prod (make-instance 'poly :termlist (reverse r)) p)))))
63 (unless succeeded-p
64 (error "#### Polynomial division Loop invariant failed ####:~%C=~A~%F=~A~%A=~A~%FL=~A~%R=~A~%P=~A~%"
65 c f a fl r p))
66 succeeded-p))
67
68
69(defun poly-pseudo-divide (f fl)
70 "Pseudo-divide a polynomial F by the list of polynomials FL. Return
71multiple values. The first value is a list of quotients A. The second
72value is the remainder R. The third argument is a scalar coefficient
73C, such that C*F can be divided by FL within the ring of coefficients,
74which is not necessarily a field. Finally, the fourth value is an
75integer count of the number of reductions performed. The resulting
76objects satisfy the equation: C*F= sum A[i]*FL[i] + R. The sugar of
77the quotients is initialized to default."
78 (declare (type poly f) (list fl))
79 ;; Loop invariant: c*f=sum ai*fi+r+p, where p must eventually become 0
80 (do ((r nil)
81 (c (make-unit-for (leading-coefficient f)))
82 (a (make-list (length fl) :initial-element (make-zero-for f)))
83 (division-count 0)
84 (p f))
85 ((universal-zerop p)
86 #+grobner-check(check-loop-invariant c f a fl r p)
87 (debug-cgb "~&~3T~d reduction~:p" division-count)
88 (when (null r) (debug-cgb " ---> 0"))
89 (values a (make-instance 'poly :termlist (nreverse r)) c division-count))
90 (declare (fixnum division-count))
91 ;; Check the loop invariant here
92 #+grobner-check(check-loop-invariant c f a fl r p)
93 (do ((fl fl (rest fl)) ;scan list of divisors
94 (b a (rest b)))
95 ((cond
96 ((endp fl) ;no division occurred
97 (push (poly-remove-term p) r) ;move lt(p) to remainder
98 t)
99 ((divides-p (leading-monomial (car fl)) (leading-monomial p)) ;division occurred
100 (incf division-count)
101 (multiple-value-bind (gcd c1 c2)
102 (universal-ezgcd (leading-coefficient (car fl)) (leading-coefficient p))
103 (declare (ignore gcd))
104 (let ((m (divide (leading-monomial p) (leading-monomial (car fl)))))
105 ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
106 (mapl #'(lambda (x)
107 (setf (car x) (multiply-by (car x) c1)))
108 a)
109 (setf r (mapc #'multiply-by r c1)
110 c (multiply-by c c1)
111 p (grobner-op c2 c1 m p (car fl)))
112 (setf (car b) (add (car b)
113 (change-class m 'term :coeff c2))))
114 t))))
115 )))
116
117(defun poly-exact-divide (f g)
118 "Divide a polynomial F by another polynomial G. Assume that exact division
119with no remainder is possible. Returns the quotient."
120 (declare (type poly f g))
121 (multiple-value-bind (quot rem coeff division-count)
122 (poly-pseudo-divide f (list g))
123 (declare (ignore division-count coeff)
124 (list quot)
125 (type poly rem)
126 (type fixnum division-count))
127 (unless (universal-zerop rem) (error "Exact division failed."))
128 (car quot)))
129
130;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
131;;
132;; An implementation of the normal form
133;;
134;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
135
136(defun normal-form-step (fl p r c division-count
137 &aux
138 (g (find (leading-monomial p) fl
139 :test #'divisible-by-p
140 :key #'leading-monomial)))
141 ;; NOTE: Currently R is a list of terms of the remainder
142 (cond
143 (g ;division possible
144 (incf division-count)
145 (multiple-value-bind (gcd cg cp)
146 (universal-ezgcd (leading-coefficient g) (leading-coefficient p))
147 (declare (ignore gcd))
148 (let ((m (divide (leading-monomial p) (leading-monomial g))))
149 ;; Multiply the equation c*f=sum ai*fi+r+p by cg.
150 (setf r (mapc #'(lambda (trm) (multiply-by trm cg)) r)
151 c (multiply-by c cg)
152 ;; p := cg*p-cp*m*g
153 p (grobner-op cp cg m p g))))
154 (debug-cgb "/"))
155 (t ;no division possible
156 (setf r (push (poly-remove-term p) r)) ;move lt(p) to remainder
157 (debug-cgb "+")))
158 (values p r c division-count))
159
160;;
161;; Merge NORMAL-FORM someday with POLY-PSEUDO-DIVIDE.
162;;
163;; TODO: It is hard to test normal form as there is no loop invariant,
164;; like for POLY-PSEUDO-DIVIDE. Is there a testing strategy? One
165;; method would be to test NORMAL-FORM using POLY-PSEUDO-DIVIDE.
166;;
167(defun normal-form (f fl &optional (top-reduction-only $poly_top_reduction_only))
168 #+grobner-check(when (null fl) (warn "normal-form: empty divisor list."))
169 (when (universal-zerop f)
170 #+grobner-check(when (null fl) (warn "normal-form: Dividend is zero."))
171 ;; NOTE: When the polynomial F is zero, we cannot constuct the
172 ;; unit in the coefficient field.
173 (return-from normal-form (values f nil 0)))
174 (do ((r nil)
175 (c (make-unit-for (leading-coefficient f)))
176 (division-count 0))
177 ((or (universal-zerop f)
178 ;;(endp fl)
179 (and top-reduction-only (not (null r))))
180 (progn
181 (debug-cgb "~&~3T~D reduction~:P" division-count)
182 (when (null r)
183 (debug-cgb " ---> 0")))
184 (setf (poly-termlist f) (nreconc r (poly-termlist f)))
185 (values f c division-count))
186 (declare (fixnum division-count))
187 (multiple-value-setq (f r c division-count)
188 (normal-form-step fl f r c division-count))))
189
190(defun buchberger-criterion (g)
191 "Returns T if G is a Grobner basis, by using the Buchberger
192criterion: for every two polynomials h1 and h2 in G the S-polynomial
193S(h1,h2) reduces to 0 modulo G."
194 (every #'universal-zerop
195 (makelist (normal-form (s-polynomial (elt g i) (elt g j)) g nil)
196 (i 0 (- (length g) 2))
197 (j (1+ i) (1- (length g))))))
198
199
200(defun poly-normalize (p &aux (c (leading-coefficient p)))
201 "Divide a polynomial by its leading coefficient. It assumes
202that the division is possible, which may not always be the
203case in rings which are not fields. The exact division operator
204is assumed to be provided by the RING structure."
205 (mapc #'(lambda (term)
206 (setf (term-coeff term) (divide (term-coeff term) c)))
207 (poly-termlist p))
208 p)
209
210(defun poly-normalize-list (plist)
211 "Divide every polynomial in a list PLIST by its leading coefficient. "
212 (mapcar #'(lambda (x) (poly-normalize x)) plist))
213
214(defun grobner-test (g f)
215 "Test whether G is a Grobner basis and F is contained in G. Return T
216upon success and NIL otherwise. The function GROBNER-TEST is provided
217primarily for debugging purposes. To enable verification of grobner
218bases with BUCHBERGER-CRITERION, do
219(pushnew :grobner-check *features*) and compile/load this file. With
220this feature, the calculations will slow down CONSIDERABLY."
221 (debug-cgb "~&GROBNER CHECK: ")
222 (let (($poly_grobner_debug nil)
223 (stat1 (buchberger-criterion g))
224 (stat2
225 (every #'universal-zerop
226 (makelist (normal-form (copy-instance (elt f i)) (mapcar #'copy-instance g) nil)
227 (i 0 (1- (length f)))))))
228 (unless stat1
229 (error "~&Buchberger criterion failed, not a grobner basis: ~A" g))
230 (unless stat2
231 (error "~&Original polynomials not in ideal spanned by Grobner basis: ~A" f)))
232 (debug-cgb "~&GROBNER CHECK END")
233 t)
Note: See TracBrowser for help on using the repository browser.