[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.
|
---|
| 38 | Top 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.
|
---|
| 50 | Assume 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
|
---|
| 73 | polynomial F by the list of polynomials FL. The invariant is the
|
---|
[1242] | 74 | identity C*F=SUM AI*FI+R+P, where F0 is the initial value of F, A is
|
---|
[1238] | 75 | the list of partial quotients, R is the intermediate value of the
|
---|
[1242] | 76 | remainder, and P is the intermediate value which eventually becomes
|
---|
[1269] | 77 | 0. A thing to remember is that the terms of polynomials in A and
|
---|
| 78 | the polynomial R have their terms in reversed order. Hence, before
|
---|
| 79 | the 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
|
---|
| 102 | multiple values. The first value is a list of quotients A. The second
|
---|
| 103 | value is the remainder R. The third argument is a scalar coefficient
|
---|
| 104 | C, such that C*F can be divided by FL within the ring of coefficients,
|
---|
| 105 | which is not necessarily a field. Finally, the fourth value is an
|
---|
| 106 | integer count of the number of reductions performed. The resulting
|
---|
[1220] | 107 | objects satisfy the equation: C*F= sum A[i]*FL[i] + R. The sugar of
|
---|
[1221] | 108 | the 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
|
---|
| 155 | with 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
|
---|
| 223 | criterion: for every two polynomials h1 and h2 in G the S-polynomial
|
---|
| 224 | S(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
|
---|
| 233 | that the division is possible, which may not always be the
|
---|
| 234 | case in rings which are not fields. The exact division operator
|
---|
[1197] | 235 | is 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))
|
---|