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

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

* empty log message *

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