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