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@ 2624

Last change on this file since 2624 was 1965, checked in by Marek Rychlik, 9 years ago

* empty log message *

File size: 11.3 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 :utils :ring :monom :polynomial :grobner-debug :term :ring-and-order)
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
35(in-package :division)
36
37(defvar $poly_top_reduction_only nil
38 "If not FALSE, use top reduction only whenever possible.
39Top reduction means that division algorithm stops after the first reduction.")
40
41;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
42;;
43;; An implementation of the division algorithm
44;;
45;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
46
47(defun grobner-op (ring-and-order c1 c2 m f g
48 &aux
49 (ring (ro-ring ring-and-order)))
50 "Returns C2*F-C1*M*G, where F and G are polynomials M is a monomial.
51Assume that the leading terms will cancel."
52 (declare (type ring-and-order ring-and-order)
53 (type monom m)
54 (type poly f g))
55 #+grobner-check(funcall (ring-zerop ring)
56 (funcall (ring-sub ring)
57 (funcall (ring-mul ring) c2 (poly-lc f))
58 (funcall (ring-mul ring) c1 (poly-lc g))))
59 #+grobner-check(monom-equal-p (poly-lm f) (monom-mul m (poly-lm g)))
60 ;; Note that below we can drop the leading terms of f ang g for the
61 ;; purpose of polynomial arithmetic.
62 ;;
63 ;; TODO: Make sure that the sugar calculation is correct if leading
64 ;; terms are dropped.
65 (poly-sub ring-and-order
66 (scalar-times-poly-1 ring c2 f)
67 (scalar-times-poly-1 ring c1 (monom-times-poly m g))))
68
69(defun check-loop-invariant (ring-and-order c f a fl r p
70 &aux
71 (ring (ro-ring ring-and-order))
72 (p-zero (make-poly-zero))
73 (a (mapcar #'poly-reverse a))
74 (r (poly-reverse r)))
75 "Check loop invariant of division algorithms, when we divide a
76polynomial F by the list of polynomials FL. The invariant is the
77identity C*F=SUM AI*FI+R+P, where F0 is the initial value of F, A is
78the list of partial quotients, R is the intermediate value of the
79remainder, and P is the intermediate value which eventually becomes
800. A thing to remember is that the terms of polynomials in A and
81the polynomial R have their terms in reversed order. Hence, before
82the arithmetic is performed, we need to fix the order of terms"
83 #|
84 (format t "~&----------------------------------------------------------------~%")
85 (format t "#### Loop invariant check ####:~%C=~A~%F=~A~%A=~A~%FL=~A~%R=~A~%P=~A~%"
86 c f a fl r p)
87 |#
88 (flet ((p-add (x y) (poly-add ring-and-order x y))
89 (p-sub (x y) (poly-sub ring-and-order x y))
90 (p-mul (x y) (poly-mul ring-and-order x y)))
91 (let* ((prod (inner-product a fl p-add p-mul p-zero))
92 (succeeded-p
93 (poly-zerop
94 (p-sub
95 (scalar-times-poly ring c f)
96 (reduce #'p-add (list prod r p))))))
97 (unless succeeded-p
98 (error "#### Polynomial division Loop invariant failed ####:~%C=~A~%F=~A~%A=~A~%FL=~A~%R=~A~%P=~A~%"
99 c f a fl r p))
100 succeeded-p)))
101
102
103(defun poly-pseudo-divide (ring-and-order f fl
104 &aux
105 (ring (ro-ring ring-and-order)))
106 "Pseudo-divide a polynomial F by the list of polynomials FL. Return
107multiple values. The first value is a list of quotients A. The second
108value is the remainder R. The third argument is a scalar coefficient
109C, such that C*F can be divided by FL within the ring of coefficients,
110which is not necessarily a field. Finally, the fourth value is an
111integer count of the number of reductions performed. The resulting
112objects satisfy the equation: C*F= sum A[i]*FL[i] + R. The sugar of
113the quotients is initialized to default."
114 (declare (type poly f) (list fl))
115 ;; Loop invariant: c*f=sum ai*fi+r+p, where p must eventually become 0
116 (do ((r (make-poly-zero))
117 (c (funcall (ring-unit ring)))
118 (a (make-list (length fl) :initial-element (make-poly-zero)))
119 (division-count 0)
120 (p f))
121 ((poly-zerop p)
122 #+grobner-check(check-loop-invariant ring-and-order c f a fl r p)
123 (debug-cgb "~&~3T~d reduction~:p" division-count)
124 (when (poly-zerop r) (debug-cgb " ---> 0"))
125 ;; We obtained the terms in reverse order, so must fix that
126 (setf a (mapcar #'poly-nreverse a)
127 r (poly-nreverse r))
128 ;; Initialize the sugar of the quotients
129 (mapc #'poly-reset-sugar a)
130 (values a r c division-count))
131 (declare (fixnum division-count))
132 ;; Check the loop invariant here
133 #+grobner-check(check-loop-invariant ring-and-order c f a fl r p)
134 (do ((fl fl (rest fl)) ;scan list of divisors
135 (b a (rest b)))
136 ((cond
137 ((endp fl) ;no division occurred
138 (push (poly-lt p) (poly-termlist r)) ;move lt(p) to remainder
139 (setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
140 (pop (poly-termlist p)) ;remove lt(p) from p
141 t)
142 ((monom-divides-p (poly-lm (car fl)) (poly-lm p)) ;division occurred
143 (incf division-count)
144 (multiple-value-bind (gcd c1 c2)
145 (funcall (ring-ezgcd ring) (poly-lc (car fl)) (poly-lc p))
146 (declare (ignore gcd))
147 (let ((m (monom-div (poly-lm p) (poly-lm (car fl)))))
148 ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
149 (mapl #'(lambda (x)
150 (setf (car x) (scalar-times-poly ring c1 (car x))))
151 a)
152 (setf r (scalar-times-poly ring c1 r)
153 c (funcall (ring-mul ring) c c1)
154 p (grobner-op ring-and-order c2 c1 m p (car fl)))
155 (push (make-term :monom m :coeff c2) (poly-termlist (car b))))
156 t))))
157 )))
158
159(defun poly-exact-divide (ring-and-order f g)
160 "Divide a polynomial F by another polynomial G. Assume that exact division
161with no remainder is possible. Returns the quotient."
162 (declare (type poly f g) (type ring-and-order ring-and-order))
163 (multiple-value-bind (quot rem coeff division-count)
164 (poly-pseudo-divide ring-and-order f (list g))
165 (declare (ignore division-count coeff)
166 (list quot)
167 (type poly rem)
168 (type fixnum division-count))
169 (unless (poly-zerop rem) (error "Exact division failed."))
170 (car quot)))
171
172;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
173;;
174;; An implementation of the normal form
175;;
176;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
177
178(defun normal-form-step (ring-and-order fl p r c division-count
179 &aux
180 (ring (ro-ring ring-and-order))
181 (g (find (poly-lm p) fl
182 :test #'monom-divisible-by-p
183 :key #'poly-lm)))
184 (cond
185 (g ;division possible
186 (incf division-count)
187 (multiple-value-bind (gcd cg cp)
188 (funcall (ring-ezgcd ring) (poly-lc g) (poly-lc p))
189 (declare (ignore gcd))
190 (let ((m (monom-div (poly-lm p) (poly-lm g))))
191 ;; Multiply the equation c*f=sum ai*fi+r+p by cg.
192 (setf r (scalar-times-poly ring cg r)
193 c (funcall (ring-mul ring) c cg)
194 ;; p := cg*p-cp*m*g
195 p (grobner-op ring-and-order cp cg m p g))))
196 (debug-cgb "/"))
197 (t ;no division possible
198 (push (poly-lt p) (poly-termlist r)) ;move lt(p) to remainder
199 (setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
200 (pop (poly-termlist p)) ;remove lt(p) from p
201 (debug-cgb "+")))
202 (values p r c division-count))
203
204;;
205;; Merge NORMAL-FORM someday with POLY-PSEUDO-DIVIDE.
206;;
207;; TODO: It is hard to test normal form as there is no loop invariant,
208;; like for POLY-PSEUDO-DIVIDE. Is there a testing strategy? One
209;; method would be to test NORMAL-FORM using POLY-PSEUDO-DIVIDE.
210;;
211(defun normal-form (ring-and-order f fl
212 &optional
213 (top-reduction-only $poly_top_reduction_only)
214 (ring (ro-ring ring-and-order)))
215 #+grobner-check(when (null fl) (warn "normal-form: empty divisor list."))
216 (do ((r (make-poly-zero))
217 (c (funcall (ring-unit ring)))
218 (division-count 0))
219 ((or (poly-zerop f)
220 ;;(endp fl)
221 (and top-reduction-only (not (poly-zerop r))))
222 (progn
223 (debug-cgb "~&~3T~D reduction~:P" division-count)
224 (when (poly-zerop r)
225 (debug-cgb " ---> 0")))
226 (setf (poly-termlist f) (nreconc (poly-termlist r) (poly-termlist f)))
227 (values f c division-count))
228 (declare (fixnum division-count)
229 (type poly r))
230 (multiple-value-setq (f r c division-count)
231 (normal-form-step ring-and-order fl f r c division-count))))
232
233(defun buchberger-criterion (ring-and-order g)
234 "Returns T if G is a Grobner basis, by using the Buchberger
235criterion: for every two polynomials h1 and h2 in G the S-polynomial
236S(h1,h2) reduces to 0 modulo G."
237 (every #'poly-zerop
238 (makelist (normal-form ring-and-order (spoly ring-and-order (elt g i) (elt g j)) g nil)
239 (i 0 (- (length g) 2))
240 (j (1+ i) (1- (length g))))))
241
242
243(defun poly-normalize (ring p &aux (c (poly-lc p)))
244 "Divide a polynomial by its leading coefficient. It assumes
245that the division is possible, which may not always be the
246case in rings which are not fields. The exact division operator
247is assumed to be provided by the RING structure."
248 (mapc #'(lambda (term)
249 (setf (term-coeff term) (funcall (ring-div ring) (term-coeff term) c)))
250 (poly-termlist p))
251 p)
252
253(defun poly-normalize-list (ring plist)
254 "Divide every polynomial in a list PLIST by its leading coefficient. "
255 (mapcar #'(lambda (x) (poly-normalize ring x)) plist))
256
257;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
258;;
259;; The function GROBNER-CHECK is provided primarily for debugging purposes. To
260;; enable verification of grobner bases with BUCHBERGER-CRITERION, do
261;; (pushnew :grobner-check *features*) and compile/load this file.
262;; With this feature, the calculations will slow down CONSIDERABLY.
263;;
264;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
265
266(defun grobner-test (ring-and-order g f)
267 "Test whether G is a Grobner basis and F is contained in G. Return T
268upon success and NIL otherwise."
269 (debug-cgb "~&GROBNER CHECK: ")
270 (let (($poly_grobner_debug nil)
271 (stat1 (buchberger-criterion ring-and-order g))
272 (stat2
273 (every #'poly-zerop
274 (makelist (normal-form ring-and-order (copy-tree (elt f i)) g nil)
275 (i 0 (1- (length f)))))))
276 (unless stat1 (error "~&Buchberger criterion failed, not a grobner basis: ~A" g))
277 (unless stat2
278 (error "~&Original polynomials not in ideal spanned by Grobner basis: ~A" f)))
279 (debug-cgb "~&GROBNER CHECK END")
280 t)
Note: See TracBrowser for help on using the repository browser.