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

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

* empty log message *

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