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