| 1 | ;;; -*- Mode: Lisp; Syntax: Common-Lisp; Package: Grobner; Base: 10 -*-
|
|---|
| 2 | #|
|
|---|
| 3 | $Id: grobner.lisp,v 1.8 2009/01/23 10:40:26 marek Exp $
|
|---|
| 4 | *--------------------------------------------------------------------------*
|
|---|
| 5 | | Copyright (C) 1994, Marek Rychlik (e-mail: rychlik@math.arizona.edu) |
|
|---|
| 6 | | Department of Mathematics, University of Arizona, Tucson, AZ 85721 |
|
|---|
| 7 | | |
|
|---|
| 8 | | Everyone is permitted to copy, distribute and modify the code in this |
|
|---|
| 9 | | directory, as long as this copyright note is preserved verbatim. |
|
|---|
| 10 | *--------------------------------------------------------------------------*
|
|---|
| 11 | |#
|
|---|
| 12 | (defpackage "GROBNER"
|
|---|
| 13 | (:use "ORDER" "PARSE" "PRINTER" "DIVISION" "MONOM"
|
|---|
| 14 | "MAKELIST" "COEFFICIENT-RING" "TERM" "POLY"
|
|---|
| 15 | "POLY-WITH-SUGAR" "COMMON-LISP")
|
|---|
| 16 | (:export
|
|---|
| 17 | spoly
|
|---|
| 18 | grobner
|
|---|
| 19 | reduction
|
|---|
| 20 | normal-form
|
|---|
| 21 | reduced-grobner
|
|---|
| 22 | normalize-poly
|
|---|
| 23 | normalize-basis
|
|---|
| 24 | elimination-ideal
|
|---|
| 25 | ring-intersection
|
|---|
| 26 | poly-contract
|
|---|
| 27 | ideal-intersection
|
|---|
| 28 | poly-lcm
|
|---|
| 29 | grobner-gcd
|
|---|
| 30 | pseudo-divide
|
|---|
| 31 | buchberger-criterion
|
|---|
| 32 | grobner-test
|
|---|
| 33 | grobner-equal
|
|---|
| 34 | grobner-subsetp
|
|---|
| 35 | grobner-member
|
|---|
| 36 | ideal-equal
|
|---|
| 37 | ideal-subsetp
|
|---|
| 38 | ideal-member
|
|---|
| 39 | ideal-saturation
|
|---|
| 40 | ideal-polysaturation
|
|---|
| 41 | ideal-saturation-1
|
|---|
| 42 | ideal-polysaturation-1
|
|---|
| 43 | grobner-primitive-part
|
|---|
| 44 | grobner-content
|
|---|
| 45 | colon-ideal
|
|---|
| 46 | *grobner-debug*
|
|---|
| 47 | buchberger-set-pair-heuristic
|
|---|
| 48 | gebauer-moeller-set-pair-heuristic
|
|---|
| 49 | select-grobner-algorithm
|
|---|
| 50 | spoly-with-sugar
|
|---|
| 51 | normal-form-with-sugar
|
|---|
| 52 | grobner-with-sugar
|
|---|
| 53 | ))
|
|---|
| 54 |
|
|---|
| 55 | (in-package "GROBNER")
|
|---|
| 56 |
|
|---|
| 57 | #+debug(proclaim '(optimize (speed 0) (debug 3)))
|
|---|
| 58 | #-debug(proclaim '(optimize (speed 3) (debug 0)))
|
|---|
| 59 |
|
|---|
| 60 | #+debug
|
|---|
| 61 | (defvar *grobner-debug* nil
|
|---|
| 62 | "If T, produce debugging and tracing output.")
|
|---|
| 63 |
|
|---|
| 64 | (defvar *buchberger-merge-pairs* 'buchberger-merge-pairs-smallest-lcm
|
|---|
| 65 | "A variable holding the predicate used to sort critical pairs
|
|---|
| 66 | in the order of decreasing priority for the Buchberger algorithm.")
|
|---|
| 67 |
|
|---|
| 68 | (defvar *gebauer-moeller-merge-pairs* 'gebauer-moeller-merge-pairs-use-mock-spoly
|
|---|
| 69 | "A variable holding the predicate used to sort critical pairs
|
|---|
| 70 | in the order of decreasing priority for the Gebauer-Moeller algorithm")
|
|---|
| 71 |
|
|---|
| 72 | (defvar *grobner-function* 'buchberger
|
|---|
| 73 | "The name of the function used to calculate Grobner basis")
|
|---|
| 74 |
|
|---|
| 75 | (defun select-grobner-algorithm (algorithm)
|
|---|
| 76 | "Select one of the several implementations of the Grobner basis
|
|---|
| 77 | algorithm."
|
|---|
| 78 | (setf *grobner-function*
|
|---|
| 79 | (ecase algorithm
|
|---|
| 80 | (:buchberger 'buchberger)
|
|---|
| 81 | (:cached-buchberger 'cached-buchberger)
|
|---|
| 82 | (:gebauer-moeller 'gebauer-moeller)
|
|---|
| 83 | (:sugar 'buchberger-with-sugar)
|
|---|
| 84 | (:cached-sugar 'cached-buchberger-with-sugar)
|
|---|
| 85 | (:gebauer-moeller-with-sugar 'gebauer-moeller-with-sugar))))
|
|---|
| 86 |
|
|---|
| 87 | (defun grobner (F &optional (pred #'lex>)
|
|---|
| 88 | (start 0)
|
|---|
| 89 | (top-reduction-only nil)
|
|---|
| 90 | (ring *coefficient-ring*))
|
|---|
| 91 | "Return Grobner basis of an ideal generated by polynomial list F.
|
|---|
| 92 | Assume that the monomials of each polynomial are ordered according to
|
|---|
| 93 | the admissible monomial order PRED. The result will be a list of
|
|---|
| 94 | polynomials ordered as well according to PRED. The polynomials from 0
|
|---|
| 95 | to START-1 in the list F are assumed to be a Grobner basis, and thus
|
|---|
| 96 | certain critical pairs do not have to be calculated. If
|
|---|
| 97 | TOP-REDUCTION-ONLY is not NIL then only top reductions will be
|
|---|
| 98 | performed, instead of full division, and thus the returned Grobner
|
|---|
| 99 | basis will have its polynomials only top-reduced with respect to the
|
|---|
| 100 | preceding polynomials. RING should be set to the RING structure (see
|
|---|
| 101 | COEFFICIENT-RING package) corresponding to the coefficients of the
|
|---|
| 102 | polynomials in the list F. This function is currently just an
|
|---|
| 103 | interface to several algorithms which fine Grobner bases. Use
|
|---|
| 104 | SELECT-GROBNER-ALGORITHM in order to change the algorithm."
|
|---|
| 105 | (funcall *grobner-function* F pred start top-reduction-only ring))
|
|---|
| 106 |
|
|---|
| 107 |
|
|---|
| 108 | #+debug
|
|---|
| 109 | (defmacro debug-cgb (&rest args)
|
|---|
| 110 | "This macro is used in printing debugging and tracing messages
|
|---|
| 111 | and its arguments are analogous to the FORMAT function, except
|
|---|
| 112 | that the stream argument is omitted. By default, the output is
|
|---|
| 113 | sent to *trace-output*, which can be set to produce output
|
|---|
| 114 | in a file."
|
|---|
| 115 | `(when *grobner-debug* (format *trace-output* ,@args)))
|
|---|
| 116 |
|
|---|
| 117 | (defun spoly (f g pred ring)
|
|---|
| 118 | "Return the S-polynomial of F and G, which should
|
|---|
| 119 | be two polynomials with terms ordered by PRED and
|
|---|
| 120 | coefficients in a ring whose operations are the slots
|
|---|
| 121 | of the RING structure."
|
|---|
| 122 | (declare (optimize (speed 3) (safety 0)))
|
|---|
| 123 | (let* ((lcm (monom-lcm (lm f) (lm g)))
|
|---|
| 124 | (m1 (monom/ lcm (lm f)))
|
|---|
| 125 | (m2 (monom/ lcm (lm g))))
|
|---|
| 126 | (let* ((c (funcall (ring-lcm ring) (lc f) (lc g)))
|
|---|
| 127 | (cf (funcall (ring-/ ring) c (lc f)))
|
|---|
| 128 | (cg (funcall (ring-/ ring) c (lc g))))
|
|---|
| 129 | (poly-
|
|---|
| 130 | (term-times-poly (cons m1 cf) (rest f) ring)
|
|---|
| 131 | (term-times-poly (cons m2 cg) (rest g) ring)
|
|---|
| 132 | pred
|
|---|
| 133 | ring))))
|
|---|
| 134 |
|
|---|
| 135 |
|
|---|
| 136 | (defun grobner-primitive-part (p ring)
|
|---|
| 137 | "Divide polynomial P with integer coefficients by gcd of its
|
|---|
| 138 | coefficients and return the result. Use the RING structure from the
|
|---|
| 139 | COEFFICIENT-RING package to perform the operations on the
|
|---|
| 140 | coefficients."
|
|---|
| 141 | (if (poly-zerop p)
|
|---|
| 142 | (values p 1)
|
|---|
| 143 | (let ((c (grobner-content p ring)))
|
|---|
| 144 | (values (mapcar
|
|---|
| 145 | #'(lambda (x)
|
|---|
| 146 | (cons (car x)
|
|---|
| 147 | (funcall
|
|---|
| 148 | (ring-/ ring)
|
|---|
| 149 | (cdr x) c)))
|
|---|
| 150 | p)
|
|---|
| 151 | c))))
|
|---|
| 152 |
|
|---|
| 153 |
|
|---|
| 154 | (defun grobner-content (p ring)
|
|---|
| 155 | "Greatest common divisor of the coefficients of the polynomial P. Use the RING structure
|
|---|
| 156 | to compute the greatest common divisor."
|
|---|
| 157 | (let ((c (apply (ring-gcd ring) (mapcar #'cdr p))))
|
|---|
| 158 | (when (/= (funcall (ring-signum ring) (cdar p))
|
|---|
| 159 | (funcall (ring-signum ring) c))
|
|---|
| 160 | (setf c (funcall (ring-- ring) c)))
|
|---|
| 161 | c))
|
|---|
| 162 |
|
|---|
| 163 | (defun normal-form (f fl pred top-reduction-only ring)
|
|---|
| 164 | "Calculate the normal form of the polynomial F with respect to
|
|---|
| 165 | the polynomial list FL. Destructive to F; thus, copy-tree should be used
|
|---|
| 166 | where F needs to be preserved. It returns multiple values.
|
|---|
| 167 | The first value is the polynomial R which is reduced (or top-reduced
|
|---|
| 168 | if TOP-REDUCTION-ONLY is not NIL). The second value is an integer
|
|---|
| 169 | which is the number of reductions actually performed. The third value
|
|---|
| 170 | is the coefficient by which F needs to be multiplied in order to be
|
|---|
| 171 | able to perform the division without actually having to divide in the
|
|---|
| 172 | coefficient ring (a form of pseudo-division is used). All operations
|
|---|
| 173 | on the coefficients are performed using the RING structure from the
|
|---|
| 174 | COEFFICIENT-RING package."
|
|---|
| 175 | (declare (optimize (speed 3) (safety 0)))
|
|---|
| 176 | ;; Loop invariant: c*f=sum ai*fi+r+p
|
|---|
| 177 | (do (r (c (ring-unit ring))
|
|---|
| 178 | (division-count 0)
|
|---|
| 179 | (p f))
|
|---|
| 180 | ((or (endp p)
|
|---|
| 181 | (endp fl)
|
|---|
| 182 | (and top-reduction-only r))
|
|---|
| 183 | #+debug(progn
|
|---|
| 184 | (debug-cgb "~&~3T~d reduction~:p" division-count)
|
|---|
| 185 | (when (endp r)
|
|---|
| 186 | (debug-cgb " ---> 0")))
|
|---|
| 187 | (values (append (nreverse r) p)
|
|---|
| 188 | division-count c))
|
|---|
| 189 | (declare (fixnum division-count) (integer c))
|
|---|
| 190 | (let ((g (find (lm p) fl :test #'monom-divisible-by-p :key #'lm)))
|
|---|
| 191 | (cond
|
|---|
| 192 | (g ;division possible
|
|---|
| 193 | (incf division-count)
|
|---|
| 194 | (let* ((gcd (funcall (ring-gcd ring) (lc g) (lc p)))
|
|---|
| 195 | (c1 (funcall (ring-/ ring) (lc g) gcd))
|
|---|
| 196 | (c2 (funcall (ring-/ ring) (lc p) gcd))
|
|---|
| 197 | (m (monom/ (lm p) (lm g))))
|
|---|
| 198 | ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
|
|---|
| 199 | (setf r (scalar-times-poly c1 r ring)
|
|---|
| 200 | c (funcall (ring-* ring) c c1)
|
|---|
| 201 | p (grobner-op c2 c1 m (rest p) (rest g) pred ring))))
|
|---|
| 202 | (t ;no division possible
|
|---|
| 203 | (setf r (cons (lt p) r) ;move lt(p) to remainder
|
|---|
| 204 | p (rest p))))))) ;remove lt(p) from p
|
|---|
| 205 |
|
|---|
| 206 | (defun buchberger (F pred start top-reduction-only ring
|
|---|
| 207 | &aux (s (1- (length F))) B M)
|
|---|
| 208 | "An implementation of the Buchberger algorithm. Return Grobner
|
|---|
| 209 | basis of the ideal generated by the polynomial list F.
|
|---|
| 210 | The terms of each polynomial are ordered by the admissible
|
|---|
| 211 | monomial order PRED. Polynomials 0 to START-1 are assumed to
|
|---|
| 212 | be a Grobner basis already, so that certain critical pairs
|
|---|
| 213 | will not be examined. If TOP-REDUCTION-ONLY set, top reduction
|
|---|
| 214 | will be preformed. The RING structure is used to perform operations
|
|---|
| 215 | on coefficents (see COEFFICIENT-RING package)."
|
|---|
| 216 | (declare (fixnum s))
|
|---|
| 217 | (unless (plusp s) (return-from buchberger F)) ;cut startup costs
|
|---|
| 218 | #+debug (debug-cgb "~&GROBNER BASIS - BUCHBERGER ALGORITHM")
|
|---|
| 219 | #+debug (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
|
|---|
| 220 | #+grobner-check (when (plusp start)
|
|---|
| 221 | (grobner-test (subseq F 0 start) (subseq F 0 start) pred ring))
|
|---|
| 222 | (setf B (nconc (makelist (list i j) (i 0 (1- start)) (j start s))
|
|---|
| 223 | (makelist (list i j) (i start (1- s)) (j (1+ i) s)))
|
|---|
| 224 | M (make-hash-table :test #'equal)
|
|---|
| 225 | B (buchberger-sort-pairs B F pred ring))
|
|---|
| 226 | ;;Initialize treated pairs
|
|---|
| 227 | (dotimes (i (1- start))
|
|---|
| 228 | (do ((j (1+ i) (1+ j))) ((>= j start))
|
|---|
| 229 | (setf (gethash (list i j) M) t)))
|
|---|
| 230 | (do ((G F))
|
|---|
| 231 | ((endp B)
|
|---|
| 232 | #+grobner-check(grobner-test G F pred ring)
|
|---|
| 233 | #+debug(debug-cgb "~&GROBNER END")
|
|---|
| 234 | G)
|
|---|
| 235 | (let ((pair (pop B)))
|
|---|
| 236 | (unless (or (Criterion-1 pair G)
|
|---|
| 237 | (Criterion-2 pair G M))
|
|---|
| 238 | (let ((SP (normal-form (spoly (elt G (first pair))
|
|---|
| 239 | (elt G (second pair))
|
|---|
| 240 | pred
|
|---|
| 241 | ring)
|
|---|
| 242 | G pred top-reduction-only ring)))
|
|---|
| 243 | (unless (poly-zerop SP) ;S-polynomial ---> 0;
|
|---|
| 244 | (setf s (1+ s)
|
|---|
| 245 | SP (grobner-primitive-part SP ring)
|
|---|
| 246 | G (nconc G (list SP))
|
|---|
| 247 | B (funcall *buchberger-merge-pairs* B (makelist (list i s)
|
|---|
| 248 | (i 0 (1- s)))
|
|---|
| 249 | G pred ring))
|
|---|
| 250 | #+debug (debug-cgb "~&Polynomials: ~d; Pairs left: ~d; Treated pairs: ~d"
|
|---|
| 251 | (length G) (length B)
|
|---|
| 252 | (hash-table-count M)))))
|
|---|
| 253 | (setf (gethash pair M) t))))
|
|---|
| 254 |
|
|---|
| 255 | (defun grobner-op (c1 c2 m A B pred ring &aux n A1 A2)
|
|---|
| 256 | "Perform a calculation equivalent to
|
|---|
| 257 | (POLY- (SCALAR-TIMES-POLY C2 A)
|
|---|
| 258 | (SCALAR-TIMES-POLY C1 (MONOM-TIMES-POLY M B))
|
|---|
| 259 | PRED)
|
|---|
| 260 | but more efficiently; it destructively modifies A and stores the
|
|---|
| 261 | result in it; A is returned."
|
|---|
| 262 | (dolist (x A) (setf (cdr x) (funcall (ring-* ring) c2 (cdr x))))
|
|---|
| 263 | (unless B
|
|---|
| 264 | (return-from grobner-op A))
|
|---|
| 265 | (unless A
|
|---|
| 266 | (return-from grobner-op
|
|---|
| 267 | (mapcar
|
|---|
| 268 | #'(lambda (x) (cons (monom* m (car x))
|
|---|
| 269 | (funcall (ring-- ring)
|
|---|
| 270 | (funcall (ring-* ring)
|
|---|
| 271 | c1 (cdr x)))))
|
|---|
| 272 | B)))
|
|---|
| 273 | (setf A1 (cons nil A) A2 A1)
|
|---|
| 274 | (tagbody
|
|---|
| 275 | begin-outer-loop
|
|---|
| 276 | (unless B
|
|---|
| 277 | (return-from grobner-op (cdr A2)))
|
|---|
| 278 | (setf n (monom* m (caar B)))
|
|---|
| 279 | (tagbody
|
|---|
| 280 | begin-inner-loop
|
|---|
| 281 | (unless (cdr A1)
|
|---|
| 282 | (setf (cdr A1)
|
|---|
| 283 | (mapcar #'(lambda (x) (cons (monom* m (car x))
|
|---|
| 284 | (funcall (ring-- ring)
|
|---|
| 285 | (funcall (ring-* ring)
|
|---|
| 286 | c1 (cdr x))))) B))
|
|---|
| 287 | (return-from grobner-op (cdr A2)))
|
|---|
| 288 | (multiple-value-bind (lm-greater lm-equal)
|
|---|
| 289 | (funcall pred (caadr A1) n)
|
|---|
| 290 | (cond
|
|---|
| 291 | (lm-equal
|
|---|
| 292 | (let ((s (funcall (ring-- ring)
|
|---|
| 293 | (cdadr A1)
|
|---|
| 294 | (funcall (ring-* ring) c1 (cdar B)))))
|
|---|
| 295 | (if (funcall (ring-zerop ring) s) ;check for cancellation
|
|---|
| 296 | (setf (cdr A1) (cddr A1))
|
|---|
| 297 | (setf (cdadr A1) s))
|
|---|
| 298 | (setf B (cdr B))
|
|---|
| 299 | (go begin-outer-loop)))
|
|---|
| 300 | (lm-greater
|
|---|
| 301 | (setf A1 (cdr A1))
|
|---|
| 302 | (go begin-inner-loop))
|
|---|
| 303 | (t ;(lm B) is greater than (lm A)
|
|---|
| 304 | (setf (cdr A1)
|
|---|
| 305 | (cons
|
|---|
| 306 | (cons
|
|---|
| 307 | n
|
|---|
| 308 | (funcall (ring-- ring)
|
|---|
| 309 | (funcall (ring-* ring) c1 (cdar B))))
|
|---|
| 310 | (cdr A1))
|
|---|
| 311 | B (cdr B)
|
|---|
| 312 | A1 (cdr A1))
|
|---|
| 313 | (go begin-outer-loop)))))))
|
|---|
| 314 |
|
|---|
| 315 | ;;----------------------------------------------------------------
|
|---|
| 316 | ;; Pairs
|
|---|
| 317 | ;;----------------------------------------------------------------
|
|---|
| 318 |
|
|---|
| 319 | (defun buchberger-sort-pairs (B G pred ring)
|
|---|
| 320 | "Sort critical pairs B by the function which is
|
|---|
| 321 | the value of the vairable *buchberger-merge-pairs*.
|
|---|
| 322 | G is the list of polynomials which the elemnts
|
|---|
| 323 | of B point into. PRED is the admissible monomial order.
|
|---|
| 324 | The RING structure holds operations on the coefficients
|
|---|
| 325 | as described in the COEFFICIENT-RING package."
|
|---|
| 326 | (funcall *buchberger-merge-pairs* nil B G pred ring))
|
|---|
| 327 |
|
|---|
| 328 |
|
|---|
| 329 | (defun mock-spoly (f g pred)
|
|---|
| 330 | "Returns an upper-bound on the leading
|
|---|
| 331 | monomial of the S-polynomial of F and G, which are two polynomials
|
|---|
| 332 | sorted by the admissible monomial order PRED."
|
|---|
| 333 | (let* ((lcm (monom-lcm (lm f) (lm g)))
|
|---|
| 334 | (m1 (monom/ lcm (lm f)))
|
|---|
| 335 | (m2 (monom/ lcm (lm g))))
|
|---|
| 336 | (cond
|
|---|
| 337 | ((endp (second f))
|
|---|
| 338 | (monom* m2 (car (second g))))
|
|---|
| 339 | ((endp (second g))
|
|---|
| 340 | (monom* m1 (car (second f))))
|
|---|
| 341 | (t (let ((n1 (monom* m2 (car (second g))))
|
|---|
| 342 | (n2 (monom* m1 (car (second f)))))
|
|---|
| 343 | (if (funcall pred n1 n2)
|
|---|
| 344 | n1
|
|---|
| 345 | n2))))))
|
|---|
| 346 |
|
|---|
| 347 | ;;----------------------------------------------------------------
|
|---|
| 348 | ;; A stragegy based on a predicted approximate S-polynomials
|
|---|
| 349 | ;;----------------------------------------------------------------
|
|---|
| 350 | (defun buchberger-merge-pairs-use-mock-spoly (B C G pred ring)
|
|---|
| 351 | "Merges lists of critical pairs B and C into one list.
|
|---|
| 352 | The pairs are assumed to be ordered according to the
|
|---|
| 353 | increasing value of MOCK-SPOLY, as determined by the admissible
|
|---|
| 354 | monomial order PRED. G is the list of polynomials which the pairs
|
|---|
| 355 | of integers in B and C point into, and RING is the ring structure
|
|---|
| 356 | defined in the COEFFICIENT-RING package."
|
|---|
| 357 | (declare (ignore ring))
|
|---|
| 358 | ;; Delete pairs consisting of two single monomials
|
|---|
| 359 | (setf C (delete-if
|
|---|
| 360 | #'(lambda (p) (and (endp (cdr (elt G (first p))))
|
|---|
| 361 | (endp (cdr (elt G (second p))))))
|
|---|
| 362 | C))
|
|---|
| 363 | (labels ((pair-key (p) (mock-spoly (elt G (first p))
|
|---|
| 364 | (elt G (second p)) pred))
|
|---|
| 365 | (pair-order (x y) (funcall pred y x)))
|
|---|
| 366 | (merge 'list B (sort C #'pair-order :key #'pair-key)
|
|---|
| 367 | #'pair-order :key #'pair-key)))
|
|---|
| 368 |
|
|---|
| 369 |
|
|---|
| 370 | ;;----------------------------------------------------------------
|
|---|
| 371 | ;; A strategy based on the smallest lcm of the leading
|
|---|
| 372 | ;; monomials of the two polynomials - so called normal strategy
|
|---|
| 373 | ;;----------------------------------------------------------------
|
|---|
| 374 | (defun buchberger-merge-pairs-smallest-lcm (B C G pred ring)
|
|---|
| 375 | "Merges lists of critical pairs B and C into a single ordered
|
|---|
| 376 | list. Implements a strategy of ordering pairs according to the
|
|---|
| 377 | smallest lcm of the leading monomials of the two polynomials - so
|
|---|
| 378 | called normal strategy."
|
|---|
| 379 | (declare (ignore ring))
|
|---|
| 380 | (labels ((pair-key (p) (monom-lcm (lm (elt G (first p)))
|
|---|
| 381 | (lm (elt G (second p)))))
|
|---|
| 382 | (pair-order (x y) (funcall pred y x)))
|
|---|
| 383 | (merge 'list B (sort C #'pair-order :key #'pair-key)
|
|---|
| 384 | #'pair-order :key #'pair-key)))
|
|---|
| 385 |
|
|---|
| 386 | ;;----------------------------------------------------------------
|
|---|
| 387 | ;; A strategy based on the smallest degree of the lcm of the leading
|
|---|
| 388 | ;; monomials of the two polynomials
|
|---|
| 389 | ;;----------------------------------------------------------------
|
|---|
| 390 | (defun buchberger-merge-pairs-use-smallest-degree (B C G pred ring)
|
|---|
| 391 | "Mergest lists B and C of critical pairs. This strategy is based on
|
|---|
| 392 | ordering the pairs according to the smallest degree of the lcm of the
|
|---|
| 393 | leading monomials of the two polynomials."
|
|---|
| 394 | (declare (ignore pred ring))
|
|---|
| 395 | (labels ((pair-key (p) (total-degree
|
|---|
| 396 | (monom-lcm (lm (elt G (first p)))
|
|---|
| 397 | (lm (elt G (second p)))))))
|
|---|
| 398 | (merge 'list B (sort C #'< :key #'pair-key) #'< :key #'pair-key)))
|
|---|
| 399 |
|
|---|
| 400 | (defun buchberger-merge-pairs-use-smallest-length (B C G pred ring)
|
|---|
| 401 | "Mergest lists B and C of critical pairs. This strategy is based on
|
|---|
| 402 | ordering the pairs according to the smallest total length of the
|
|---|
| 403 | two polynomials, where length is the number of terms."
|
|---|
| 404 | (declare (ignore pred ring))
|
|---|
| 405 | (labels ((pair-key (p) (+
|
|---|
| 406 | (length (elt G (first p)))
|
|---|
| 407 | (length (elt G (second p))))))
|
|---|
| 408 | (merge 'list B (sort C #'< :key #'pair-key) #'< :key #'pair-key)))
|
|---|
| 409 |
|
|---|
| 410 | (defun buchberger-merge-pairs-use-smallest-coefficient-length (B C G pred ring)
|
|---|
| 411 | "Mergest lists B and C of critical pairs. This strategy is based on
|
|---|
| 412 | ordering the pairs according to the smallest combined length of the
|
|---|
| 413 | coefficients of the two polynomials."
|
|---|
| 414 | (declare (ignore pred))
|
|---|
| 415 | (labels ((coefficient-length (poly)
|
|---|
| 416 | (apply #'+ (mapcar (ring-length ring)
|
|---|
| 417 | (mapcar #'cdr poly))))
|
|---|
| 418 | (pair-key (p) (+
|
|---|
| 419 | (coefficient-length (elt G (first p)))
|
|---|
| 420 | (coefficient-length (elt G (second p))))))
|
|---|
| 421 | (merge 'list B (sort C #'< :key #'pair-key) #'< :key #'pair-key)))
|
|---|
| 422 |
|
|---|
| 423 | (defun buchberger-set-pair-heuristic
|
|---|
| 424 | (method
|
|---|
| 425 | &aux
|
|---|
| 426 | (strategy-fn
|
|---|
| 427 | (ecase method
|
|---|
| 428 | (:minimal-mock-spoly #'buchberger-merge-pairs-use-mock-spoly)
|
|---|
| 429 | (:minimal-lcm #'buchberger-merge-pairs-smallest-lcm)
|
|---|
| 430 | (:minimal-total-degree #'buchberger-merge-pairs-use-smallest-degree)
|
|---|
| 431 | (:minimal-length #'buchberger-merge-pairs-use-smallest-length)
|
|---|
| 432 | (:minimal-coefficient-length #'buchberger-merge-pairs-use-smallest-coefficient-length))))
|
|---|
| 433 | "Simply sets the variable *buchberger-merge-pairs* to the heuristic
|
|---|
| 434 | function STRATEGY-FN implementing one of several strategies introduces
|
|---|
| 435 | before of selecting the most promising. METHOD should be one of the
|
|---|
| 436 | listed keywords."
|
|---|
| 437 | (setf *buchberger-merge-pairs* strategy-fn))
|
|---|
| 438 |
|
|---|
| 439 | ;;----------------------------------------------------------------
|
|---|
| 440 | ;; Grobner Criteria
|
|---|
| 441 | ;;----------------------------------------------------------------
|
|---|
| 442 | (defun Criterion-1 (pair G &aux (i (first pair)) (j (second pair)))
|
|---|
| 443 | "Returns T if the leading monomials of the two polynomials
|
|---|
| 444 | in G pointed to by the integers in PAIR have disjoint (relatively prime)
|
|---|
| 445 | monomials. This test is known as the first Buchberger criterion."
|
|---|
| 446 | (monom-rel-prime (lm (elt G i))
|
|---|
| 447 | (lm (elt G j))))
|
|---|
| 448 |
|
|---|
| 449 |
|
|---|
| 450 | (defun Criterion-2 (pair G M &aux (i (first pair)) (j (second pair)))
|
|---|
| 451 | "Returns T if the leading monomial of some element P of G divides
|
|---|
| 452 | the LCM of the leading monomials of the two polynomials in the
|
|---|
| 453 | polynomial list G, pointed to by the two integers in PAIR, and P
|
|---|
| 454 | paired with each of the polynomials pointed to by the the PAIR has
|
|---|
| 455 | already been treated, as indicated by the hash table M, which stores
|
|---|
| 456 | value T for every treated pair so far."
|
|---|
| 457 | (labels ((make-pair (u v) (if (< u v) (list u v) (list v u))))
|
|---|
| 458 | (dotimes (k (length G) nil)
|
|---|
| 459 | (when (and (/= k i)
|
|---|
| 460 | (/= k j)
|
|---|
| 461 | (monom-divides-p
|
|---|
| 462 | (lm (elt G k))
|
|---|
| 463 | (monom-lcm (lm (elt G i)) (lm (elt G j))))
|
|---|
| 464 | (gethash (make-pair i k) M)
|
|---|
| 465 | (gethash (make-pair k j) M))
|
|---|
| 466 | #+debug(when *grobner-debug* (format t ":B2"))
|
|---|
| 467 | (return-from Criterion-2 t)))))
|
|---|
| 468 |
|
|---|
| 469 | ;;----------------------------------------------------------------
|
|---|
| 470 | ;; End of GROBNER implementation
|
|---|
| 471 | ;;----------------------------------------------------------------
|
|---|
| 472 |
|
|---|
| 473 | (defun normalize-poly (p ring)
|
|---|
| 474 | "Divide a polynomial by its leading coefficient. It assumes
|
|---|
| 475 | that the division is possible, which may not always be the
|
|---|
| 476 | case in rings which are not fields. The exact division operator
|
|---|
| 477 | is assumed to be provided by the RING structure of the
|
|---|
| 478 | COEFFICIENT-RING package."
|
|---|
| 479 | (mapcar #'(lambda (term)
|
|---|
| 480 | (cons
|
|---|
| 481 | (car term)
|
|---|
| 482 | (funcall (ring-/ ring) (cdr term) (lc p)))) p))
|
|---|
| 483 |
|
|---|
| 484 | (defun normalize-basis (plist ring)
|
|---|
| 485 | "Divide every polynomial in a list PLIST by its leading coefficient.
|
|---|
| 486 | Use RING structure to perform the division of the coefficients."
|
|---|
| 487 | (mapcar #'(lambda (x) (normalize-poly x ring)) plist))
|
|---|
| 488 |
|
|---|
| 489 | (defun reduction (P pred ring)
|
|---|
| 490 | "Reduce a list of polynomials P, so that non of the terms in any of the
|
|---|
| 491 | polynomials is divisible by a leading monomial of another polynomial.
|
|---|
| 492 | Return the reduced list."
|
|---|
| 493 | (do ((Q P)
|
|---|
| 494 | (found t))
|
|---|
| 495 | ((not found)
|
|---|
| 496 | (mapcar #'(lambda (x) (grobner-primitive-part x ring)) Q))
|
|---|
| 497 | ;;Find p in Q such that p is reducible mod Q\{p}
|
|---|
| 498 | (setf found nil)
|
|---|
| 499 | (dolist (x Q)
|
|---|
| 500 | (let ((Q1 (remove x Q)))
|
|---|
| 501 | (multiple-value-bind (h div-count)
|
|---|
| 502 | (normal-form x Q1
|
|---|
| 503 | pred
|
|---|
| 504 | nil ;not a top reduction
|
|---|
| 505 | ring)
|
|---|
| 506 | (unless (zerop div-count)
|
|---|
| 507 | (setf found t Q Q1)
|
|---|
| 508 | (when h
|
|---|
| 509 | (setf Q (nconc Q1 (list h))))
|
|---|
| 510 | (return)))))))
|
|---|
| 511 |
|
|---|
| 512 |
|
|---|
| 513 | (defun reduced-grobner (F &optional (pred #'lex>)
|
|---|
| 514 | (start 0)
|
|---|
| 515 | (top-reduction-only nil)
|
|---|
| 516 | (ring *coefficient-ring*))
|
|---|
| 517 | "Return the reduced Grobner basis of the ideal generated by a polynomial
|
|---|
| 518 | list F. This combines calls to two functions: GROBNER and REDUCTION.
|
|---|
| 519 | The parameters have the same meaning as in GROBNER or BUCHBERGER."
|
|---|
| 520 | (reduction (grobner F pred start top-reduction-only ring) pred ring))
|
|---|
| 521 |
|
|---|
| 522 | #+kcl
|
|---|
| 523 | (defvar *vars* nil "Keep track of variables for tracing SPOLY.")
|
|---|
| 524 |
|
|---|
| 525 | ;; Does the monomial P depend on variable K
|
|---|
| 526 | (defun monom-depends-p (m k)
|
|---|
| 527 | "Return T if the monomial M depends on variable number K."
|
|---|
| 528 | (plusp (elt m k)))
|
|---|
| 529 |
|
|---|
| 530 | ;; Does the term depend on variable K
|
|---|
| 531 | (defun term-depends-p (term k)
|
|---|
| 532 | "Return T if the term TERM depends on variable number K."
|
|---|
| 533 | (monom-depends-p (car term) k))
|
|---|
| 534 |
|
|---|
| 535 | ;; Does the polynomial P depend on variable K
|
|---|
| 536 | (defun poly-depends-p (p k)
|
|---|
| 537 | "Return T if the term polynomial P depends on variable number K."
|
|---|
| 538 | (some #'(lambda (term) (term-depends-p term k)) p))
|
|---|
| 539 |
|
|---|
| 540 | (defun ring-intersection (plist k &key (key #'identity))
|
|---|
| 541 | "This function assumes that polynomial list PLIST is a Grobner basis
|
|---|
| 542 | and it calculates the intersection with the ring R[x[k+1],...,x[n]], i.e.
|
|---|
| 543 | it discards polynomials which depend on variables x[0], x[1], ..., x[k]."
|
|---|
| 544 | (dotimes (i k plist)
|
|---|
| 545 | (setf plist
|
|---|
| 546 | (remove-if #'(lambda (p)
|
|---|
| 547 | (poly-depends-p p (funcall key i)))
|
|---|
| 548 | plist))))
|
|---|
| 549 |
|
|---|
| 550 | (defun elimination-ideal (flist k
|
|---|
| 551 | &key
|
|---|
| 552 | (primary-order #'lex>)
|
|---|
| 553 | (secondary-order #'lex>)
|
|---|
| 554 | (start 0)
|
|---|
| 555 | (order (elimination-order
|
|---|
| 556 | k
|
|---|
| 557 | :primary-order primary-order
|
|---|
| 558 | :secondary-order secondary-order
|
|---|
| 559 | ))
|
|---|
| 560 | (top-reduction-only nil)
|
|---|
| 561 | (ring *coefficient-ring*))
|
|---|
| 562 | "Returns the K-th elimination ideal of the ideal generated by polynomial list FLIST.
|
|---|
| 563 | Thus, a Grobner basis of the ideal generated by FLIST is calculated and
|
|---|
| 564 | polynomials depending on variables 0-K are discarded. The monomial order in the
|
|---|
| 565 | calculation is an elimination order formed by combining two orders: PRIMARY-ORDER
|
|---|
| 566 | used on variables 0-K and SECONDARY-ORDER used on variables starting from variable K+1.
|
|---|
| 567 | Thus, if both orders are set to the lexicographic order #'LEX> then the resulting
|
|---|
| 568 | order is simply equivalent to #'LEX>. But one could use arbitrary two orders in this
|
|---|
| 569 | function, for instance, two copies of #'GREVLEX> (see ORDER package). When doing so,
|
|---|
| 570 | the caller must ensure that the same order has been used to sort the terms of
|
|---|
| 571 | the polynomials FLIST."
|
|---|
| 572 | (ring-intersection
|
|---|
| 573 | (reduced-grobner flist order start top-reduction-only ring)
|
|---|
| 574 | k))
|
|---|
| 575 |
|
|---|
| 576 |
|
|---|
| 577 | (defun ideal-intersection (F G pred top-reduction-only ring)
|
|---|
| 578 | "Return the Grobner basis of the intersection of the ideals
|
|---|
| 579 | generated by the polynomial lists F and G. PRED is an admissible
|
|---|
| 580 | monomial order with respect to which the terms of F and G have been
|
|---|
| 581 | sorted. This order is going to be used in the Grobner basis
|
|---|
| 582 | calculation. If TOP-REDUCTION-ONLY is not NIL, internally the Grobner
|
|---|
| 583 | basis algorithm will perform top reduction only. The RING parameter,
|
|---|
| 584 | as usual, should be set to a structure defined in the package
|
|---|
| 585 | COEFFICIENT-RING, and it will be used to perform all operations on
|
|---|
| 586 | coefficients."
|
|---|
| 587 | (mapcar
|
|---|
| 588 | #'poly-contract
|
|---|
| 589 | (ring-intersection
|
|---|
| 590 | (reduced-grobner
|
|---|
| 591 | (append
|
|---|
| 592 | (mapcar #'(lambda (p) (poly-extend p (list 1))) F)
|
|---|
| 593 | (mapcar #'(lambda (p)
|
|---|
| 594 | (append (poly-extend (minus-poly p ring) (list 1))
|
|---|
| 595 | (poly-extend p)))
|
|---|
| 596 | G))
|
|---|
| 597 | (elimination-order-1 pred)
|
|---|
| 598 | 0
|
|---|
| 599 | top-reduction-only
|
|---|
| 600 | ring)
|
|---|
| 601 | 1)))
|
|---|
| 602 |
|
|---|
| 603 |
|
|---|
| 604 | (defun poly-contract (f &optional (k 1))
|
|---|
| 605 | "Return a polynomial obtained from polynomial F by dropping the
|
|---|
| 606 | first K variables, if F does not depend on them. Note that this
|
|---|
| 607 | operation makes the polynomial incompatible for certain arithmetical
|
|---|
| 608 | operations with the original polynomial F. Calling this function on a polynomial F
|
|---|
| 609 | which does depend on variables 0-K will result in error."
|
|---|
| 610 | (when (some #'(lambda (term) (find-if-not #'zerop (car term) :end k)) f)
|
|---|
| 611 | (error "Contracting a polynomial which depends on contracted variables"))
|
|---|
| 612 | (mapcar #'(lambda (term) (cons (nthcdr k (car term)) (cdr term))) f))
|
|---|
| 613 |
|
|---|
| 614 |
|
|---|
| 615 | (defun poly-lcm (f g &optional (pred #'lex>) (ring *coefficient-ring*))
|
|---|
| 616 | "Return LCM (least common multiple) of two polynomials F and G.
|
|---|
| 617 | The polynomials must be ordered according to monomial order PRED
|
|---|
| 618 | and their coefficients must be compatible with the RING structure
|
|---|
| 619 | defined in the COEFFICIENT-RING package."
|
|---|
| 620 | (cond
|
|---|
| 621 | ((endp f) nil)
|
|---|
| 622 | ((endp g) nil)
|
|---|
| 623 | ((and (endp (cdr f)) (endp (cdr g)))
|
|---|
| 624 | (list (cons (monom-lcm (caar f) (caar g))
|
|---|
| 625 | (lcm (cdar f) (cdar g)))))
|
|---|
| 626 | (t
|
|---|
| 627 | (multiple-value-bind (f f-cont)
|
|---|
| 628 | (grobner-primitive-part f ring)
|
|---|
| 629 | (multiple-value-bind (g g-cont)
|
|---|
| 630 | (grobner-primitive-part g ring)
|
|---|
| 631 | (scalar-times-poly
|
|---|
| 632 | (funcall (ring-lcm ring) f-cont g-cont)
|
|---|
| 633 | (grobner-primitive-part (car (ideal-intersection (list f) (list g) pred nil ring))
|
|---|
| 634 | ring)
|
|---|
| 635 | ring))))))
|
|---|
| 636 |
|
|---|
| 637 | ;; GCD of two polynomials, using Grobner bases
|
|---|
| 638 | (defun grobner-gcd (f g &optional (pred #'lex>) (ring *coefficient-ring*))
|
|---|
| 639 | "Return GCD (greatest common divisor) of two polynomials F and G.
|
|---|
| 640 | The polynomials must be ordered according to monomial order PRED
|
|---|
| 641 | and their coefficients must be compatible with the RING structure
|
|---|
| 642 | defined in the COEFFICIENT-RING package."
|
|---|
| 643 | (poly-exact-divide (poly* f g pred ring)
|
|---|
| 644 | (poly-lcm f g pred ring)
|
|---|
| 645 | pred ring))
|
|---|
| 646 |
|
|---|
| 647 | ;; Do two Grobner bases yield the same ideal?
|
|---|
| 648 | (defun grobner-equal (g1 g2 &optional (pred #'lex>) (ring *coefficient-ring*))
|
|---|
| 649 | "Returns T if two lists of polynomials G1 and G2, assumed to be Grobner bases,
|
|---|
| 650 | generate the same ideal, and NIL otherwise."
|
|---|
| 651 | (and (grobner-subsetp g1 g2 pred ring)
|
|---|
| 652 | (grobner-subsetp g2 g1 pred ring)))
|
|---|
| 653 |
|
|---|
| 654 | (defun grobner-subsetp (g1 g2 &optional (pred #'lex>) (ring *coefficient-ring*))
|
|---|
| 655 | "Returns T if a list of polynomials G1 generates
|
|---|
| 656 | an ideal contained in the ideal generated by a polynomial list G2,
|
|---|
| 657 | both G1 and G2 assumed to be Grobner bases. Returns NIL otherwise."
|
|---|
| 658 | (every #'(lambda (p) (grobner-member p g2 pred ring)) g1))
|
|---|
| 659 |
|
|---|
| 660 | (defun grobner-member (p g &optional (pred #'lex>) (ring *coefficient-ring*))
|
|---|
| 661 | "Returns T if a polynomial P belongs to the ideal generated by the
|
|---|
| 662 | polynomial list G, which is assumed to be a Grobner basis. Returns NIL otherwise."
|
|---|
| 663 | (endp (normal-form (copy-tree p) g pred nil ring)))
|
|---|
| 664 |
|
|---|
| 665 | (defun ideal-equal (f1 f2 pred top-reduction-only ring)
|
|---|
| 666 | "Returns T if two ideals generated by polynomial lists F1 and F2 are identical.
|
|---|
| 667 | Returns NIL otherwise."
|
|---|
| 668 | (grobner-equal (grobner f1 pred 0 top-reduction-only ring)
|
|---|
| 669 | (grobner f2 pred 0 top-reduction-only ring)
|
|---|
| 670 | pred ring))
|
|---|
| 671 |
|
|---|
| 672 | (defun ideal-subsetp (f1 f2 &optional (pred #'lex>) (ring *coefficient-ring*))
|
|---|
| 673 | "Returns T if the ideal spanned by the polynomial list F1 is contained
|
|---|
| 674 | in the ideal spanned by the polynomial list F2. Returns NIL otherwise."
|
|---|
| 675 | (grobner-subsetp
|
|---|
| 676 | (grobner f1 pred 0 nil ring)
|
|---|
| 677 | (grobner f2 pred 0 nil ring)
|
|---|
| 678 | pred ring))
|
|---|
| 679 |
|
|---|
| 680 | (defun ideal-member (p plist pred top-reduction-only ring)
|
|---|
| 681 | "Returns T if the polynomial P belongs to the ideal spanned by the polynomial
|
|---|
| 682 | list PLIST, and NIL otherwise."
|
|---|
| 683 | (endp (normal-form (copy-tree p) (grobner plist pred 0 top-reduction-only ring)
|
|---|
| 684 | pred nil ring)))
|
|---|
| 685 |
|
|---|
| 686 | ;; Calculate F : p^inf
|
|---|
| 687 | (defun ideal-saturation-1 (F p pred start top-reduction-only ring
|
|---|
| 688 | &aux (pred (elimination-order-1 pred)))
|
|---|
| 689 | "Returns the reduced Grobner basis of the saturation of the ideal
|
|---|
| 690 | generated by a polynomial list F in the ideal generated by a single
|
|---|
| 691 | polynomial P. The saturation ideal is defined as the set of
|
|---|
| 692 | polynomials H such for some natural number n (* (EXPT P N) H) is in the ideal
|
|---|
| 693 | F. Geometrically, over an algebraically closed field, this is the set
|
|---|
| 694 | of polynomials in the ideal generated by F which do not identically
|
|---|
| 695 | vanish on the variety of P."
|
|---|
| 696 | (mapcar
|
|---|
| 697 | #'poly-contract
|
|---|
| 698 | (ring-intersection
|
|---|
| 699 | (reduced-grobner
|
|---|
| 700 | (saturation-extension-1 F p ring)
|
|---|
| 701 | pred
|
|---|
| 702 | start top-reduction-only ring)
|
|---|
| 703 | 1)))
|
|---|
| 704 |
|
|---|
| 705 |
|
|---|
| 706 | (defun add-variables (plist n)
|
|---|
| 707 | "Add N new variables, adn the first N variables, to every polynomial in polynomial list PLIST.
|
|---|
| 708 | The resulting polynomial will not depend on these N variables."
|
|---|
| 709 | (mapcar
|
|---|
| 710 | #'(lambda (g)
|
|---|
| 711 | (mapcar #'(lambda (q) (cons (append (make-list n :initial-element 0)
|
|---|
| 712 | (car q))
|
|---|
| 713 | (cdr q))) g))
|
|---|
| 714 | plist))
|
|---|
| 715 |
|
|---|
| 716 | ;; Calculate <u1*p1,u2*p2,...,uk*pk>
|
|---|
| 717 | (defun extend-polynomials (plist &aux (k (length plist)))
|
|---|
| 718 | "Returns polynomial list {U1*P1', U2*P2', ... , UK*PK'}
|
|---|
| 719 | where Ui are new variables and PLIST={P1, P2, ... , PK} is a polynomial
|
|---|
| 720 | list. PI' is obtained from PI by adding new variables U1, U2, ..., UN
|
|---|
| 721 | UK as the first K variables. Thus, the resulting list is consists of polynomials whose
|
|---|
| 722 | terms are ordered by the lexicographic order on variables UI, with
|
|---|
| 723 | ties resolved by the monomial order which was used to order the terms
|
|---|
| 724 | of the polynomials in PLIST. We note that the monomial order does not
|
|---|
| 725 | explicitly participate in this calculation, and neither do the
|
|---|
| 726 | variable names."
|
|---|
| 727 | (labels ((incf-power (g i)
|
|---|
| 728 | (dolist (x g g)
|
|---|
| 729 | (incf (nth i (car x))))))
|
|---|
| 730 | (setf plist (add-variables plist k))
|
|---|
| 731 | (dotimes (i k plist)
|
|---|
| 732 | (setf (nth i plist) (incf-power (nth i plist) i)))))
|
|---|
| 733 |
|
|---|
| 734 | ;; Calculate <F, u1*p1-1,u2*p2-1,...,uk*pk-1>
|
|---|
| 735 | (defun saturation-extension (F plist ring &aux (k (length plist))
|
|---|
| 736 | (d (+ k (length (caaar plist)))))
|
|---|
| 737 | "Returns F' union {U1*P1-1,U2*P2-1,...,UK*PK-1} where Ui are new variables
|
|---|
| 738 | and F' is a polynomial list obtained from F by adding variables U1, U2, ..., Uk
|
|---|
| 739 | as the first K variables to each polynomial in F. Thus, the resulting list is
|
|---|
| 740 | consists of polynomials whose terms are ordered by the lexicographic order on
|
|---|
| 741 | variables Ui, with ties resolved by the monomial order which was used to order the
|
|---|
| 742 | terms of the polynomials in F. We note that the monomial order does not explicitly
|
|---|
| 743 | participate in this calculation, and neither do the variable names."
|
|---|
| 744 | (setf F (add-variables F k)
|
|---|
| 745 | plist (mapcar #'(lambda (x)
|
|---|
| 746 | (nconc x (list (cons (make-list d :initial-element 0)
|
|---|
| 747 | (funcall (ring-- ring)
|
|---|
| 748 | (ring-unit ring))))))
|
|---|
| 749 | (extend-polynomials plist)))
|
|---|
| 750 | (append F plist))
|
|---|
| 751 |
|
|---|
| 752 | ;; Calculate <F, u1*p1+u2*p2+...+uk*pk-1>
|
|---|
| 753 | (defun polysaturation-extension (F plist ring &aux (k (length plist))
|
|---|
| 754 | (d (+ k (length (caaar plist)))))
|
|---|
| 755 | "Returns F' union {U1*P1+U2*P2+UK*PK-1} where Ui are new variables
|
|---|
| 756 | and F' is a polynomial list obtained from F by adding variables U1, U2, ..., Uk
|
|---|
| 757 | as the first K variables to each polynomial in F. Thus, the resulting list is
|
|---|
| 758 | consists of polynomials whose terms are ordered by the lexicographic order on
|
|---|
| 759 | variables Ui, with ties resolved by the monomial order which was used to order the
|
|---|
| 760 | terms of the polynomials in F. We note that the monomial order does not explicitly
|
|---|
| 761 | participate in this calculation, and neither do the variable names."
|
|---|
| 762 | (setf F (add-variables F k)
|
|---|
| 763 | plist (apply #'append (extend-polynomials plist))
|
|---|
| 764 | (cdr (last plist)) (list (cons (make-list d :initial-element 0)
|
|---|
| 765 | (funcall (ring-- ring)
|
|---|
| 766 | (ring-unit ring)))))
|
|---|
| 767 | (append F (list plist)))
|
|---|
| 768 |
|
|---|
| 769 | (defun saturation-extension-1 (F p ring)
|
|---|
| 770 | "Return the list F' union {U*P-1} where U is a new variable,
|
|---|
| 771 | where F' is obtained from the polynomial list F by adding U as the first variable.
|
|---|
| 772 | The variable U will be added as the first variable, so that it can be easily
|
|---|
| 773 | eliminated."
|
|---|
| 774 | (saturation-extension F (list p) ring))
|
|---|
| 775 |
|
|---|
| 776 | ;; Calculate F : p1^inf : p2^inf : ... : ps^inf
|
|---|
| 777 | (defun ideal-polysaturation-1 (F plist pred start top-reduction-only ring)
|
|---|
| 778 | "Returns the reduced Grobner basis of the ideal obtained by a
|
|---|
| 779 | sequence of successive saturations in the polynomials
|
|---|
| 780 | of the polynomial list PLIST of the ideal generated by the
|
|---|
| 781 | polynomial list F."
|
|---|
| 782 | (cond
|
|---|
| 783 | ((endp plist) (reduced-grobner F pred start top-reduction-only ring))
|
|---|
| 784 | (t (let ((G (ideal-saturation-1 F (car plist) pred start top-reduction-only ring)))
|
|---|
| 785 | (ideal-polysaturation-1 G (rest plist) pred (length G) top-reduction-only ring)))))
|
|---|
| 786 |
|
|---|
| 787 | (defun ideal-saturation (F G pred start top-reduction-only ring
|
|---|
| 788 | &aux (k (length G))
|
|---|
| 789 | (pred (elimination-order k
|
|---|
| 790 | :primary-order #'lex>
|
|---|
| 791 | :secondary-order pred)))
|
|---|
| 792 | "Returns the reduced Grobner basis of the saturation of the ideal
|
|---|
| 793 | generated by a polynomial list F in the ideal generated a polynomial
|
|---|
| 794 | list G The saturation ideal is defined as the set of polynomials H
|
|---|
| 795 | such for some natural number n and some P in the ideal generated by G
|
|---|
| 796 | the polynomial P**N * H is in the ideal spanned by F. Geometrically,
|
|---|
| 797 | over an algebraically closed field, this is the set of polynomials in
|
|---|
| 798 | the ideal generated by F which do not identically vanish on the
|
|---|
| 799 | variety of G."
|
|---|
| 800 | (mapcar
|
|---|
| 801 | #'(lambda (q) (poly-contract q k))
|
|---|
| 802 | (ring-intersection
|
|---|
| 803 | (reduced-grobner (polysaturation-extension F G ring)
|
|---|
| 804 | pred start
|
|---|
| 805 | top-reduction-only ring)
|
|---|
| 806 | k)))
|
|---|
| 807 |
|
|---|
| 808 | (defun ideal-polysaturation (F ideal-list pred start top-reduction-only ring)
|
|---|
| 809 | "Returns the reduced Grobner basis of the ideal obtained by a
|
|---|
| 810 | successive applications of IDEAL-SATURATIONS to F and
|
|---|
| 811 | lists of polynomials in the list IDEAL-LIST."
|
|---|
| 812 | (cond
|
|---|
| 813 | ((endp ideal-list) F)
|
|---|
| 814 | (t (let ((H (ideal-saturation F (car ideal-list) pred start top-reduction-only ring)))
|
|---|
| 815 | (ideal-polysaturation H (rest ideal-list) pred (length H) top-reduction-only ring)))))
|
|---|
| 816 |
|
|---|
| 817 | (defun buchberger-criterion (G pred ring)
|
|---|
| 818 | "Returns T if G is a Grobner basis, by using the Buchberger
|
|---|
| 819 | criterion: for every two polynomials h1 and h2 in G the S-polynomial
|
|---|
| 820 | S(h1,h2) reduces to 0 modulo G."
|
|---|
| 821 | (every
|
|---|
| 822 | #'endp
|
|---|
| 823 | (makelist (normal-form (spoly (elt G i) (elt G j) pred ring) G pred nil ring)
|
|---|
| 824 | (i 0 (- (length G) 2))
|
|---|
| 825 | (j (1+ i) (1- (length G))))))
|
|---|
| 826 |
|
|---|
| 827 |
|
|---|
| 828 | (defun grobner-test (G F pred ring)
|
|---|
| 829 | "Test whether G is a Grobner basis and F is contained in G. Return T
|
|---|
| 830 | upon success and NIL otherwise."
|
|---|
| 831 | #+debug(debug-cgb "~&GROBNER CHECK: ")
|
|---|
| 832 | (let ((*grobner-debug* nil)
|
|---|
| 833 | (stat1 (buchberger-criterion G pred ring))
|
|---|
| 834 | (stat2
|
|---|
| 835 | (every #'endp
|
|---|
| 836 | (makelist (normal-form (copy-tree (elt F i)) G pred nil ring)
|
|---|
| 837 | (i 0 (1- (length F)))))))
|
|---|
| 838 | (unless stat1 (error "~&Buchberger criterion failed."))
|
|---|
| 839 | (unless stat2
|
|---|
| 840 | (error "~&Original polys not in ideal spanned by Grobner.")))
|
|---|
| 841 | #+debug(debug-cgb "~&GROBNER CHECK END")
|
|---|
| 842 | T)
|
|---|
| 843 |
|
|---|
| 844 |
|
|---|
| 845 | (defun minimization (P pred)
|
|---|
| 846 | "Returns a sublist of the polynomial list P spanning the same
|
|---|
| 847 | monomial ideal as P but minimal, i.e. no leading monomial
|
|---|
| 848 | of a polynomial in the sublist divides the leading monomial
|
|---|
| 849 | of another polynomial."
|
|---|
| 850 | (declare (ignore pred))
|
|---|
| 851 | (do ((Q P)
|
|---|
| 852 | (found t))
|
|---|
| 853 | ((not found) Q)
|
|---|
| 854 | ;;Find p in Q such that lm(p) is in LM(Q\{p})
|
|---|
| 855 | (setf found nil
|
|---|
| 856 | Q (dolist (x Q Q)
|
|---|
| 857 | (let ((Q1 (remove x Q)))
|
|---|
| 858 | (when (member-if #'(lambda (p) (monom-divides-p (lm x) (lm p))) Q1)
|
|---|
| 859 | (setf found t)
|
|---|
| 860 | (return Q1)))))))
|
|---|
| 861 |
|
|---|
| 862 |
|
|---|
| 863 | (defun add-minimized (f Gred pred)
|
|---|
| 864 | "Adds a polynomial f to GRED, reduced Grobner basis, preserving the
|
|---|
| 865 | property described in the documentation for MINIMIZATION."
|
|---|
| 866 | (if (member-if #'(lambda (p) (monom-divides-p (lm p) (lm f))) Gred)
|
|---|
| 867 | Gred
|
|---|
| 868 | (merge 'list (list f)
|
|---|
| 869 | (delete-if #'(lambda (p) (monom-divides-p (lm f) (lm p))) Gred)
|
|---|
| 870 | pred :key #'lm)))
|
|---|
| 871 |
|
|---|
| 872 | (defun colon-ideal (F G pred top-reduction-only ring)
|
|---|
| 873 | "Returns the reduced Grobner basis of the colon ideal Id(F):Id(G), where
|
|---|
| 874 | F and G are two lists of polynomials. The colon ideal I:J is defined as the
|
|---|
| 875 | set of polynomials H such that there is a polynomial W in J for which
|
|---|
| 876 | W*H belongs to I."
|
|---|
| 877 | (cond
|
|---|
| 878 | ((endp G)
|
|---|
| 879 | (if (every #'endp F) nil
|
|---|
| 880 | (list (list
|
|---|
| 881 | (cons
|
|---|
| 882 | (make-list (length (caar (find-if-not #'endp F)))
|
|---|
| 883 | :initial-element 0)
|
|---|
| 884 | 1)))))
|
|---|
| 885 | ((endp (cdr G))
|
|---|
| 886 | (colon-ideal-1 F (car G) pred top-reduction-only ring))
|
|---|
| 887 | (t
|
|---|
| 888 | (ideal-intersection (colon-ideal-1 F (car G) pred top-reduction-only ring)
|
|---|
| 889 | (colon-ideal F (rest G) pred top-reduction-only ring)
|
|---|
| 890 | pred
|
|---|
| 891 | top-reduction-only
|
|---|
| 892 | ring))))
|
|---|
| 893 |
|
|---|
| 894 | (defun colon-ideal-1 (F g pred top-reduction-only ring)
|
|---|
| 895 | "Returns the reduced Grobner basis of the colon ideal Id(F):Id({G}), where
|
|---|
| 896 | F is a list of polynomials and G is a polynomial."
|
|---|
| 897 | (mapcar #'(lambda (x) (poly-exact-divide x g pred ring))
|
|---|
| 898 | (ideal-intersection F (list g) pred top-reduction-only ring)))
|
|---|
| 899 |
|
|---|
| 900 |
|
|---|
| 901 | (defun pseudo-divide (f fl pred ring)
|
|---|
| 902 | "Pseudo-divide a polynomial F by the list of polynomials FL. Return
|
|---|
| 903 | multiple values. The first value is a list of quotients A.
|
|---|
| 904 | The second value is the remainder R. The third value is an integer
|
|---|
| 905 | count of the number of reductions performed. Finally, the fourth
|
|---|
| 906 | argument is a scalar coefficient C, such that C*F can be divided by FL
|
|---|
| 907 | within the ring of coefficients, which is not necessarily a field.
|
|---|
| 908 | The resulting objects satisfy the equation:
|
|---|
| 909 | C*F= sum A[i]*FL[i] + R
|
|---|
| 910 | "
|
|---|
| 911 | (declare (optimize (speed 3) (safety 0)))
|
|---|
| 912 | (do (r
|
|---|
| 913 | (c (ring-unit ring))
|
|---|
| 914 | (a (make-list (length fl)))
|
|---|
| 915 | (division-count 0)
|
|---|
| 916 | (p f))
|
|---|
| 917 | ((endp p)
|
|---|
| 918 | #+debug(debug-cgb "~&~3T~d reduction~:p" division-count)
|
|---|
| 919 | #+debug(when (endp r) (debug-cgb " ---> 0"))
|
|---|
| 920 | (values a (nreverse r) division-count c))
|
|---|
| 921 | (declare (fixnum division-count) (integer c))
|
|---|
| 922 | (do ((fl fl (rest fl)) ;scan list of divisors
|
|---|
| 923 | (b a (rest b)))
|
|---|
| 924 | ((cond
|
|---|
| 925 | ((endp fl) ;no division occurred
|
|---|
| 926 | (setf r (cons (lt p) r) ;move lt(p) to remainder
|
|---|
| 927 | p (rest p)) ;remove lt(p) from p
|
|---|
| 928 | t)
|
|---|
| 929 | ((monom-divides-p (lm (car fl)) (lm p)) ;division occurred
|
|---|
| 930 | (incf division-count)
|
|---|
| 931 | (let* ((gcd (funcall (ring-gcd ring) (lc (car fl)) (lc p)))
|
|---|
| 932 | (c1 (funcall (ring-/ ring) (lc (car fl)) gcd))
|
|---|
| 933 | (c2 (funcall (ring-/ ring) (lc p) gcd))
|
|---|
| 934 | (quot (cons (monom/ (lm p) (lm (car fl))) c2)))
|
|---|
| 935 | ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
|
|---|
| 936 | (mapl #'(lambda (x)
|
|---|
| 937 | (setf (car x) (scalar-times-poly c1 (car x) ring)))
|
|---|
| 938 | a)
|
|---|
| 939 | (setf p (scalar-times-poly c1 (rest p) ring) ;drop first term
|
|---|
| 940 | r (scalar-times-poly c1 r ring)
|
|---|
| 941 | c (funcall (ring-* ring) c c1)
|
|---|
| 942 | p (poly-op p quot (rest (car fl)) pred ring)
|
|---|
| 943 | (car b) (nconc (car b) (list quot))))
|
|---|
| 944 | t))))))
|
|---|
| 945 |
|
|---|
| 946 | ;;----------------------------------------------------------------
|
|---|
| 947 | ;; An implementation of the algorithm of Gebauer and Moeller,
|
|---|
| 948 | ;; as described in the book of Becker-Weispfenning, p. 232
|
|---|
| 949 | ;;----------------------------------------------------------------
|
|---|
| 950 | (defun gebauer-moeller (F pred start top-reduction-only ring
|
|---|
| 951 | &aux B G F1)
|
|---|
| 952 | "Compute Grobner basis by using the algorithm of Gebauer and Moeller.
|
|---|
| 953 | This algorithm is described as BUCHBERGERNEW2 in the book by
|
|---|
| 954 | Becker-Weispfenning entitled ``Grobner Bases''"
|
|---|
| 955 | ;;cut startup costs if (length F)<=1
|
|---|
| 956 | (cond
|
|---|
| 957 | ((endp F) (return-from gebauer-moeller nil))
|
|---|
| 958 | ((endp (cdr F))
|
|---|
| 959 | (return-from gebauer-moeller (list (grobner-primitive-part (car F) ring)))))
|
|---|
| 960 | #+debug (debug-cgb "~&GROBNER BASIS - GEBAUER MOELLER ALGORITHM")
|
|---|
| 961 | #+debug (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
|
|---|
| 962 | #+grobner-check (when (plusp start)
|
|---|
| 963 | (grobner-test (subseq F 0 start) (subseq F 0 start)
|
|---|
| 964 | pred ring))
|
|---|
| 965 | (setf B nil
|
|---|
| 966 | G (subseq F 0 start)
|
|---|
| 967 | F1 (subseq F start))
|
|---|
| 968 | (do () ((endp F1))
|
|---|
| 969 | (multiple-value-setq (G B)
|
|---|
| 970 | (update G B (grobner-primitive-part (pop F1) ring) pred ring)))
|
|---|
| 971 | (do () ((endp B))
|
|---|
| 972 | (let* ((pair (pop B))
|
|---|
| 973 | (g1 (first pair))
|
|---|
| 974 | (g2 (second pair))
|
|---|
| 975 | (h (normal-form (spoly g1 g2 pred ring)
|
|---|
| 976 | G pred top-reduction-only ring)))
|
|---|
| 977 | (unless (poly-zerop h)
|
|---|
| 978 | (setf h (grobner-primitive-part h ring))
|
|---|
| 979 | (multiple-value-setq (G B)
|
|---|
| 980 | (update G B h pred ring))
|
|---|
| 981 | #+debug (debug-cgb "~&Polynomials: ~d; Pairs left: ~d~%"
|
|---|
| 982 | (length G) (length B))
|
|---|
| 983 | )))
|
|---|
| 984 | #+grobner-check(grobner-test G F pred ring)
|
|---|
| 985 | #+debug(debug-cgb "~&GROBNER END")
|
|---|
| 986 | G)
|
|---|
| 987 |
|
|---|
| 988 | (defun update (G B h pred ring &aux C D E B-new G-new pair)
|
|---|
| 989 | "An implementation of the auxillary UPDATE algorithm used by the
|
|---|
| 990 | Gebauer-Moeller algorithm. G is a list of polynomials, B is a list of
|
|---|
| 991 | critical pairs and H is a new polynomial which possibly will be added
|
|---|
| 992 | to G. The naming conventions used are very close to the one used in
|
|---|
| 993 | the book of Becker-Weispfenning."
|
|---|
| 994 | (setf C G D nil)
|
|---|
| 995 | (do (g1) ((endp C))
|
|---|
| 996 | (setf g1 (pop C))
|
|---|
| 997 | (when (or (monom-rel-prime (lm h) (lm g1))
|
|---|
| 998 | (and
|
|---|
| 999 | (not (some
|
|---|
| 1000 | #'(lambda (g2) (monom-divides-p (monom-lcm (lm h) (lm g2))
|
|---|
| 1001 | (monom-lcm (lm h) (lm g1))))
|
|---|
| 1002 | C))
|
|---|
| 1003 | (not (some
|
|---|
| 1004 | #'(lambda (g2) (monom-divides-p (monom-lcm (lm h) (lm g2))
|
|---|
| 1005 | (monom-lcm (lm h) (lm g1))))
|
|---|
| 1006 | D))))
|
|---|
| 1007 | (push g1 D)))
|
|---|
| 1008 | (setf E nil)
|
|---|
| 1009 | (do (g1) ((endp D))
|
|---|
| 1010 | (setf g1 (pop D))
|
|---|
| 1011 | (unless (monom-rel-prime (lm h) (lm g1))
|
|---|
| 1012 | (push g1 E)))
|
|---|
| 1013 | (setf B-new nil)
|
|---|
| 1014 | (do (g1 g2) ((endp B))
|
|---|
| 1015 | (setf pair (pop B)
|
|---|
| 1016 | g1 (first pair)
|
|---|
| 1017 | g2 (second pair))
|
|---|
| 1018 | (when (or (not (monom-divides-p (lm h) (monom-lcm (lm g1) (lm g2))))
|
|---|
| 1019 | (equal (monom-lcm (lm g1) (lm h))
|
|---|
| 1020 | (monom-lcm (lm g1) (lm g2)))
|
|---|
| 1021 | (equal (monom-lcm (lm h) (lm g2))
|
|---|
| 1022 | (monom-lcm (lm g1) (lm g2))))
|
|---|
| 1023 | (setf B-new (funcall *gebauer-moeller-merge-pairs* B-new (list (list g1 g2)) pred ring))))
|
|---|
| 1024 | (dolist (g3 E)
|
|---|
| 1025 | (setf B-new (funcall *gebauer-moeller-merge-pairs* B-new (list (list h g3)) pred ring)))
|
|---|
| 1026 | (setf G-new nil)
|
|---|
| 1027 | (do (g1) ((endp G))
|
|---|
| 1028 | (setf g1 (pop G))
|
|---|
| 1029 | (unless (monom-divides-p (lm h) (lm g1))
|
|---|
| 1030 | (push g1 G-new)))
|
|---|
| 1031 | (push h G-new)
|
|---|
| 1032 | (values G-new B-new))
|
|---|
| 1033 |
|
|---|
| 1034 | (defun gebauer-moeller-merge-pairs-use-mock-spoly (B C pred ring)
|
|---|
| 1035 | "The merging strategy used by the Gebauer-Moeller algorithm, based
|
|---|
| 1036 | on MOCK-SPOLY; see the documentation of
|
|---|
| 1037 | BUCHBERGER-MERGE-PAIRS-USE-MOCK-SPOLY."
|
|---|
| 1038 | (declare (ignore ring))
|
|---|
| 1039 | ;; Delete pairs consisting of two single monomials
|
|---|
| 1040 | (setf C (delete-if
|
|---|
| 1041 | #'(lambda (p) (and (endp (cdr (first p)))
|
|---|
| 1042 | (endp (cdr (second p)))))
|
|---|
| 1043 | C))
|
|---|
| 1044 | (labels ((pair-key (p) (mock-spoly (first p) (second p) pred))
|
|---|
| 1045 | (pair-order (x y) (funcall pred y x)))
|
|---|
| 1046 | (merge 'list B C #'pair-order :key #'pair-key)))
|
|---|
| 1047 |
|
|---|
| 1048 | (defun gebauer-moeller-merge-pairs-smallest-lcm (B C pred ring)
|
|---|
| 1049 | "The merging strategy based on the smallest-lcm (normal strategy); see
|
|---|
| 1050 | the documentation of BUCHBERGER-MERGE-PAIRS-SMALLEST-LCM."
|
|---|
| 1051 | (declare (ignore ring))
|
|---|
| 1052 | ;; Delete pairs consisting of two single monomials
|
|---|
| 1053 | (labels ((pair-key (p) (monom-lcm (lm (first p))
|
|---|
| 1054 | (lm (second p))))
|
|---|
| 1055 | (pair-order (x y) (funcall pred y x)))
|
|---|
| 1056 | (merge 'list B C #'pair-order :key #'pair-key)))
|
|---|
| 1057 |
|
|---|
| 1058 |
|
|---|
| 1059 | (defun gebauer-moeller-merge-pairs-use-smallest-degree (B C pred ring)
|
|---|
| 1060 | "The merging strategy based on the smallest-lcm (normal strategy); see
|
|---|
| 1061 | the documentation of BUCHBERGER-MERGE-PAIRS-USE-SMALLEST-DEGREE."
|
|---|
| 1062 | (declare (ignore ring))
|
|---|
| 1063 | (declare (ignore pred))
|
|---|
| 1064 | (labels ((pair-key (p) (total-degree
|
|---|
| 1065 | (monom-lcm (lm (first p))
|
|---|
| 1066 | (lm (second p)))))
|
|---|
| 1067 | (pair-order (x y) (funcall #'< x y)))
|
|---|
| 1068 | (merge 'list B C #'pair-order :key #'pair-key)))
|
|---|
| 1069 |
|
|---|
| 1070 |
|
|---|
| 1071 | (defun gebauer-moeller-merge-pairs-use-smallest-length (B C pred ring)
|
|---|
| 1072 | (declare (ignore ring))
|
|---|
| 1073 | "The merging strategy based on the smallest-lcm (normal strategy); see
|
|---|
| 1074 | the documentation of BUCHBERGER-MERGE-PAIRS-USE-SMALLEST-LENGTH."
|
|---|
| 1075 | (declare (ignore pred))
|
|---|
| 1076 | (labels ((pair-key (p) (+
|
|---|
| 1077 | (length (first p))
|
|---|
| 1078 | (length (second p))))
|
|---|
| 1079 | (pair-order (x y) (funcall #'< x y)))
|
|---|
| 1080 | (if B
|
|---|
| 1081 | (merge 'list B C #'pair-order :key #'pair-key)
|
|---|
| 1082 | (sort C #'pair-order :key #'pair-key))))
|
|---|
| 1083 |
|
|---|
| 1084 |
|
|---|
| 1085 | (defun gebauer-moeller-merge-pairs-use-smallest-coefficient-length (B C pred ring)
|
|---|
| 1086 | "The merging strategy based on the smallest-lcm (normal strategy); see
|
|---|
| 1087 | the documentation of BUCHBERGER-MERGE-PAIRS-USE-SMALLEST-COEFFICIENT-LENGTH."
|
|---|
| 1088 | (declare (ignore pred))
|
|---|
| 1089 | (labels ((coefficient-length (poly)
|
|---|
| 1090 | (apply #'+ (mapcar (ring-length ring) (mapcar #'cdr poly))))
|
|---|
| 1091 | (pair-key (p) (+
|
|---|
| 1092 | (coefficient-length (first p))
|
|---|
| 1093 | (coefficient-length (second p))))
|
|---|
| 1094 | (pair-order (x y) (funcall #'< x y)))
|
|---|
| 1095 | (merge 'list B C #'pair-order :key #'pair-key)))
|
|---|
| 1096 |
|
|---|
| 1097 | (defun gebauer-moeller-set-pair-heuristic
|
|---|
| 1098 | (method
|
|---|
| 1099 | &aux (strategy-fn
|
|---|
| 1100 | (ecase method
|
|---|
| 1101 | (:minimal-mock-spoly #'gebauer-moeller-merge-pairs-use-mock-spoly)
|
|---|
| 1102 | (:minimal-lcm #'gebauer-moeller-merge-pairs-smallest-lcm)
|
|---|
| 1103 | (:minimal-total-degree #'gebauer-moeller-merge-pairs-use-smallest-degree)
|
|---|
| 1104 | (:minimal-length #'gebauer-moeller-merge-pairs-use-smallest-length)
|
|---|
| 1105 | (:minimal-coefficient-length #'gebauer-moeller-merge-pairs-use-smallest-coefficient-length))))
|
|---|
| 1106 | (setf *gebauer-moeller-merge-pairs* strategy-fn))
|
|---|
| 1107 |
|
|---|
| 1108 | ;;------------------------------------------------------------------------------------------------
|
|---|
| 1109 | ;;
|
|---|
| 1110 | ;; An implementation of the sugar strategy
|
|---|
| 1111 | ;;
|
|---|
| 1112 | ;;------------------------------------------------------------------------------------------------
|
|---|
| 1113 |
|
|---|
| 1114 | ;; Calculate sugar of the S-polynomial of f and g
|
|---|
| 1115 | (defun spoly-sugar (f-with-sugar g-with-sugar ring)
|
|---|
| 1116 | "Calculate the ``sugar'' of the S-polynomial of two polynomials with
|
|---|
| 1117 | sugar. A polynomial with sugar is simply a pair (P . SUGAR) where
|
|---|
| 1118 | SUGAR is an integer constant defined according to the following
|
|---|
| 1119 | algorighm: the sugar of sum or difference of two polynomials with
|
|---|
| 1120 | sugar is the MAX of the sugars of those two polynomials. The sugar of
|
|---|
| 1121 | a product of a term and a polynomial is the sum of the degree of the
|
|---|
| 1122 | term and the sugar of the polynomial. The idea is to accumulate sugar
|
|---|
| 1123 | as we perform the arithmetic operations, and that polynomials with
|
|---|
| 1124 | small (little) sugar should be given priority in the calculations.
|
|---|
| 1125 | Thus, the ``sugar strategy'' of the critical pair selection is to
|
|---|
| 1126 | select a pair with the smallest value of sugar of the resulting
|
|---|
| 1127 | S-polynomial."
|
|---|
| 1128 | (let* ((f (poly-with-sugar-poly f-with-sugar))
|
|---|
| 1129 | (g (poly-with-sugar-poly g-with-sugar))
|
|---|
| 1130 | (lcm (monom-lcm (lm f) (lm g)))
|
|---|
| 1131 | (m1 (monom/ lcm (lm f)))
|
|---|
| 1132 | (m2 (monom/ lcm (lm g))))
|
|---|
| 1133 | (let* ((c (funcall (ring-lcm ring) (lc f) (lc g)))
|
|---|
| 1134 | (cf (funcall (ring-/ ring) c (lc f)))
|
|---|
| 1135 | (cg (funcall (ring-/ ring) c (lc g))))
|
|---|
| 1136 | (max
|
|---|
| 1137 | (+ (term-sugar (cons m1 cf) ring) (poly-with-sugar-sugar f-with-sugar))
|
|---|
| 1138 | (+ (term-sugar (cons m2 cg) ring) (poly-with-sugar-sugar g-with-sugar))))))
|
|---|
| 1139 |
|
|---|
| 1140 | ;; Calculate the S-polynomial of f and g with sugar
|
|---|
| 1141 | (defun spoly-with-sugar (f-with-sugar g-with-sugar pred ring)
|
|---|
| 1142 | "The S-polynomials of two polynomials with SUGAR strategy."
|
|---|
| 1143 | (let* ((f (poly-with-sugar-poly f-with-sugar))
|
|---|
| 1144 | (g (poly-with-sugar-poly g-with-sugar))
|
|---|
| 1145 | (lcm (monom-lcm (lm f) (lm g)))
|
|---|
| 1146 | (m1 (monom/ lcm (lm f)))
|
|---|
| 1147 | (m2 (monom/ lcm (lm g))))
|
|---|
| 1148 | (let* ((c (funcall (ring-lcm ring) (lc f) (lc g)))
|
|---|
| 1149 | (cf (funcall (ring-/ ring) c (lc f)))
|
|---|
| 1150 | (cg (funcall (ring-/ ring) c (lc g))))
|
|---|
| 1151 | (poly-with-sugar-
|
|---|
| 1152 | (term-times-poly-with-sugar (cons m1 cf) (poly-with-sugar-tail f-with-sugar) ring)
|
|---|
| 1153 | (term-times-poly-with-sugar (cons m2 cg) (poly-with-sugar-tail g-with-sugar) ring)
|
|---|
| 1154 | pred
|
|---|
| 1155 | ring))))
|
|---|
| 1156 |
|
|---|
| 1157 |
|
|---|
| 1158 | (defun normal-form-with-sugar (f fl pred top-reduction-only ring)
|
|---|
| 1159 | "Normal form of the polynomial with sugar F with respect to
|
|---|
| 1160 | a list of polynomials with sugar FL. Assumes that FL is not empty.
|
|---|
| 1161 | Parameters PRED, TOP-REDUCTION-ONLY and RING have the same meaning as in NORMAL-FORM.
|
|---|
| 1162 | The parameter SUGAR-LIMIT should be set to a positive integer. If the sugar limit of
|
|---|
| 1163 | the partial remainder exceeds SUGAR-LIMIT then the calculation is stopped and the partial
|
|---|
| 1164 | remainder is returned, although it is not fully reduced with respect to FL."
|
|---|
| 1165 | (declare (optimize (speed 3) (safety 0)))
|
|---|
| 1166 | ;; Loop invariant: c*f=sum ai*fi+r+p
|
|---|
| 1167 | (do ((r (cons nil 0))
|
|---|
| 1168 | (c (ring-unit ring))
|
|---|
| 1169 | (division-count 0)
|
|---|
| 1170 | (p f))
|
|---|
| 1171 | ((or (poly-with-sugar-zerop p)
|
|---|
| 1172 | (and top-reduction-only (not (poly-with-sugar-zerop r))))
|
|---|
| 1173 | #+debug
|
|---|
| 1174 | (progn
|
|---|
| 1175 | (debug-cgb "~&~3T~d reduction~:p" division-count)
|
|---|
| 1176 | (cond
|
|---|
| 1177 | ((poly-with-sugar-zerop r)
|
|---|
| 1178 | (debug-cgb " ---> 0"))))
|
|---|
| 1179 | (values (poly-with-sugar-append (poly-with-sugar-nreverse r) p)
|
|---|
| 1180 | division-count c))
|
|---|
| 1181 | (declare (fixnum division-count) (integer c))
|
|---|
| 1182 | (let ((g (find (poly-with-sugar-lm p) fl
|
|---|
| 1183 | :test #'monom-divisible-by-p
|
|---|
| 1184 | :key #'poly-with-sugar-lm)))
|
|---|
| 1185 | (cond
|
|---|
| 1186 | (g ;division possible
|
|---|
| 1187 | (incf division-count)
|
|---|
| 1188 | (let* ((gcd (funcall (ring-gcd ring) (poly-with-sugar-lc g) (poly-with-sugar-lc p)))
|
|---|
| 1189 | (c1 (funcall (ring-/ ring) (poly-with-sugar-lc g) gcd))
|
|---|
| 1190 | (c2 (funcall (ring-/ ring) (poly-with-sugar-lc p) gcd))
|
|---|
| 1191 | (m (monom/ (poly-with-sugar-lm p) (poly-with-sugar-lm g))))
|
|---|
| 1192 | ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
|
|---|
| 1193 | (setf r (scalar-times-poly-with-sugar c1 r ring)
|
|---|
| 1194 | c (funcall (ring-* ring) c c1)
|
|---|
| 1195 | p (scalar-times-poly-with-sugar c1 p ring)
|
|---|
| 1196 | p (poly-with-sugar-op p (cons m c2) g pred ring))))
|
|---|
| 1197 | (t ;no division possible
|
|---|
| 1198 | ;;move lt(p) to remainder
|
|---|
| 1199 | (setf r (cons (cons (poly-with-sugar-lt p) (poly-with-sugar-poly r))
|
|---|
| 1200 | (poly-with-sugar-sugar r))
|
|---|
| 1201 | p (poly-with-sugar-tail p)))))))
|
|---|
| 1202 |
|
|---|
| 1203 | ;; Parameter START is used when Grobner basis is calculated incrementally
|
|---|
| 1204 | ;; and it signals that the elements of F from 0 to START-1 form a Grobner
|
|---|
| 1205 | ;; basis already, so their pairs are not formed.
|
|---|
| 1206 | (defun buchberger-with-sugar (F-no-sugar pred start top-reduction-only ring
|
|---|
| 1207 | &aux (s (1- (length F-no-sugar))) B M F)
|
|---|
| 1208 | "Returns a Grobner basis of an ideal generated by a list of ordinary
|
|---|
| 1209 | polynomials F-NO-SUGAR. Adds initial sugar to the polynomials and
|
|---|
| 1210 | performs the Buchberger algorithm with ``sugar strategy''. It returns
|
|---|
| 1211 | an ordinary list of polynomials with no sugar. One of the most
|
|---|
| 1212 | effective algorithms for Grobner basis calculation."
|
|---|
| 1213 | (declare (fixnum s))
|
|---|
| 1214 | (unless (plusp s) (return-from buchberger-with-sugar F-no-sugar)) ;cut startup costs
|
|---|
| 1215 | #+debug (debug-cgb "~&GROBNER BASIS - BUCHBERGER WITH SUGER STRATEGY ALGORITHM")
|
|---|
| 1216 | #+debug (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
|
|---|
| 1217 | #+grobner-check (when (plusp start)
|
|---|
| 1218 | (grobner-test (subseq F-no-sugar 0 start)
|
|---|
| 1219 | (subseq F-no-sugar 0 start) pred ring))
|
|---|
| 1220 | (setf F (mapcar #'(lambda (x) (poly-add-sugar x ring)) F-no-sugar)
|
|---|
| 1221 | B (nconc (makelist (list i j) (i 0 (1- start)) (j start s))
|
|---|
| 1222 | (makelist (list i j) (i start (1- s)) (j (1+ i) s)))
|
|---|
| 1223 | M (make-hash-table :test #'equal))
|
|---|
| 1224 | ;;Initialize treated pairs
|
|---|
| 1225 | (dotimes (i (1- start))
|
|---|
| 1226 | (do ((j (1+ i) (1+ j))) ((>= j start))
|
|---|
| 1227 | (setf (gethash (list i j) M) t)))
|
|---|
| 1228 | (setf B (buchberger-with-sugar-sort-pairs B F M pred ring))
|
|---|
| 1229 | (do ((G F))
|
|---|
| 1230 | ((endp B)
|
|---|
| 1231 | (setf G (mapcar #'poly-with-sugar-poly G))
|
|---|
| 1232 | ;;#+grobner-check(grobner-test G F-no-sugar pred ring)
|
|---|
| 1233 | #+debug(debug-cgb "~&GROBNER END")
|
|---|
| 1234 | G)
|
|---|
| 1235 | (let ((pair (pop B)))
|
|---|
| 1236 | (cond
|
|---|
| 1237 | ((Criterion-1-with-sugar pair G))
|
|---|
| 1238 | ((Criterion-2-with-sugar pair G M))
|
|---|
| 1239 | (t
|
|---|
| 1240 | (let ((SP (normal-form-with-sugar
|
|---|
| 1241 | (cons (spoly (poly-with-sugar-poly (elt G (first pair)))
|
|---|
| 1242 | (poly-with-sugar-poly (elt G (second pair))) pred ring)
|
|---|
| 1243 | (third pair))
|
|---|
| 1244 | G pred top-reduction-only ring)))
|
|---|
| 1245 | (cond
|
|---|
| 1246 | ((poly-with-sugar-zerop SP))
|
|---|
| 1247 | (t
|
|---|
| 1248 | (setf s (1+ s)
|
|---|
| 1249 | (poly-with-sugar-poly SP) (grobner-primitive-part
|
|---|
| 1250 | (poly-with-sugar-poly SP) ring)
|
|---|
| 1251 | (cdr (last G)) (list SP)
|
|---|
| 1252 | B (buchberger-with-sugar-merge-pairs B (makelist (list i s)
|
|---|
| 1253 | (i 0 (1- s)))
|
|---|
| 1254 | G M pred ring))
|
|---|
| 1255 | #+debug (debug-cgb "~&[~d] Polynomials: ~d; Pairs left: ~d; Treated pairs: ~d"
|
|---|
| 1256 | (third pair)
|
|---|
| 1257 | (length G) (length B)
|
|---|
| 1258 | (hash-table-count M)))))))
|
|---|
| 1259 | (setf (cddr pair) nil
|
|---|
| 1260 | (gethash pair M) t))))
|
|---|
| 1261 |
|
|---|
| 1262 |
|
|---|
| 1263 | (defun buchberger-with-sugar-merge-pairs (B C G M pred ring)
|
|---|
| 1264 | "Merges lists of critical pairs. It orders pairs according to
|
|---|
| 1265 | increasing sugar, with ties broken by smaller MOCK-SPOLY value of the
|
|---|
| 1266 | two polynomials. In this function B must already be sorted."
|
|---|
| 1267 | ;;Remove pairs consisting of two monomials
|
|---|
| 1268 | ;;These pairs produce zero S-polynomial
|
|---|
| 1269 | (setf C (delete-if
|
|---|
| 1270 | #'(lambda (p)
|
|---|
| 1271 | (when (and (endp (cdr (poly-with-sugar-poly (elt G (first p)))))
|
|---|
| 1272 | (endp (cdr (poly-with-sugar-poly (elt G (second p))))))
|
|---|
| 1273 | ;;Record the pair as treated
|
|---|
| 1274 | (setf (gethash p M) t)
|
|---|
| 1275 | t))
|
|---|
| 1276 | C))
|
|---|
| 1277 | (labels ((pair-sugar (p)
|
|---|
| 1278 | (spoly-sugar (elt G (first p)) (elt G (second p)) ring))
|
|---|
| 1279 | (pair-mock-spoly (p)
|
|---|
| 1280 | (mock-spoly (poly-with-sugar-poly (elt G (first p)))
|
|---|
| 1281 | (poly-with-sugar-poly (elt G (second p))) pred))
|
|---|
| 1282 | (pair-order (p q)
|
|---|
| 1283 | (let ((sugar-p (third p))
|
|---|
| 1284 | (sugar-q (third q)))
|
|---|
| 1285 | (cond
|
|---|
| 1286 | ((< sugar-p sugar-q) t)
|
|---|
| 1287 | ((> sugar-p sugar-q) nil)
|
|---|
| 1288 | (t
|
|---|
| 1289 | (funcall pred (fourth q) (fourth p)))))))
|
|---|
| 1290 | ;;Add sugar and mock S-polynomial as the third and fourth element
|
|---|
| 1291 | ;;of new pair in C
|
|---|
| 1292 | (dolist (p C)
|
|---|
| 1293 | (setf (cdr (last p))
|
|---|
| 1294 | (list (pair-sugar p) (pair-mock-spoly p))))
|
|---|
| 1295 | (merge 'list B (sort C #'pair-order) #'pair-order)))
|
|---|
| 1296 |
|
|---|
| 1297 | (defun buchberger-with-sugar-sort-pairs (C G M pred ring)
|
|---|
| 1298 | "Sorts critical pairs C according to sugar strategy."
|
|---|
| 1299 | (buchberger-with-sugar-merge-pairs nil C G M pred ring))
|
|---|
| 1300 |
|
|---|
| 1301 | (defun Criterion-1-with-sugar (pair G &aux (i (first pair)) (j (second pair)))
|
|---|
| 1302 | "An implementation of Buchberger's first criterion for polynomials
|
|---|
| 1303 | with sugar. See the documentation of BUCHBERGER-CRITERION-1."
|
|---|
| 1304 | (monom-rel-prime (poly-with-sugar-lm (elt G i))
|
|---|
| 1305 | (poly-with-sugar-lm (elt G j))))
|
|---|
| 1306 |
|
|---|
| 1307 | (defun Criterion-2-with-sugar (pair G M &aux (i (first pair)) (j (second pair)))
|
|---|
| 1308 | "An implementation of Buchberger's first criterion for polynomials
|
|---|
| 1309 | with sugar. See the documentation of BUCHBERGER-CRITERION-2."
|
|---|
| 1310 | (labels ((make-pair (u v) (if (< u v) (list u v) (list v u))))
|
|---|
| 1311 | (dotimes (k (length G) nil)
|
|---|
| 1312 | (when (and (/= k i)
|
|---|
| 1313 | (/= k j)
|
|---|
| 1314 | (monom-divides-p
|
|---|
| 1315 | (poly-with-sugar-lm (elt G k))
|
|---|
| 1316 | (monom-lcm (poly-with-sugar-lm (elt G i))
|
|---|
| 1317 | (poly-with-sugar-lm (elt G j))))
|
|---|
| 1318 | (gethash (make-pair i k) M)
|
|---|
| 1319 | (gethash (make-pair k j) M))
|
|---|
| 1320 | #+debug(when *grobner-debug* (format t ":B2"))
|
|---|
| 1321 | (return-from Criterion-2-with-sugar t)))))
|
|---|
| 1322 |
|
|---|
| 1323 |
|
|---|
| 1324 | ;;----------------------------------------------------------------
|
|---|
| 1325 | ;; An implementation of the algorithm of Gebauer and Moeller,
|
|---|
| 1326 | ;; as described in the book of Becker-Weispfenning, p. 232
|
|---|
| 1327 | ;; with sugar strategy
|
|---|
| 1328 | ;;----------------------------------------------------------------
|
|---|
| 1329 | (defun gebauer-moeller-with-sugar (F pred start top-reduction-only ring
|
|---|
| 1330 | &aux B G F1)
|
|---|
| 1331 | "Compute Grobner basis by using the algorithm of Gebauer and Moeller.
|
|---|
| 1332 | This algorithm is described as BUCHBERGERNEW2 in the book by
|
|---|
| 1333 | Becker-Weispfenning entitled ``Grobner Bases''"
|
|---|
| 1334 | ;;cut startup costs if (length F)<=1
|
|---|
| 1335 | (cond
|
|---|
| 1336 | ((endp F) (return-from gebauer-moeller-with-sugar nil))
|
|---|
| 1337 | ((endp (cdr F))
|
|---|
| 1338 | (return-from gebauer-moeller-with-sugar (list (grobner-primitive-part (car F) ring)))))
|
|---|
| 1339 | #+debug (debug-cgb "~&GROBNER BASIS - GEBAUER MOELLER ALGORITHM WITH SUGAR STRATEGY")
|
|---|
| 1340 | #+debug (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
|
|---|
| 1341 | #+grobner-check (when (plusp start)
|
|---|
| 1342 | (grobner-test (subseq F 0 start) (subseq F 0 start)
|
|---|
| 1343 | pred ring))
|
|---|
| 1344 | (setf B nil
|
|---|
| 1345 | G (subseq F 0 start)
|
|---|
| 1346 | F1 (subseq F start))
|
|---|
| 1347 | ;;Initialize sugar
|
|---|
| 1348 | (setf G (mapcar #'(lambda (x) (poly-add-sugar x ring)) G)
|
|---|
| 1349 | F1 (mapcar #'(lambda (x) (poly-add-sugar x ring)) F1))
|
|---|
| 1350 | (do () ((endp F1))
|
|---|
| 1351 | (multiple-value-setq (G B)
|
|---|
| 1352 | (update-with-sugar G B (grobner-primitive-part-with-sugar (pop F1) ring) pred ring)))
|
|---|
| 1353 | (do ()
|
|---|
| 1354 | ((endp B))
|
|---|
| 1355 | (let* ((pair (pop B))
|
|---|
| 1356 | (g1 (first pair))
|
|---|
| 1357 | (g2 (second pair))
|
|---|
| 1358 | (h (normal-form-with-sugar (spoly-with-sugar g1 g2 pred ring)
|
|---|
| 1359 | G pred top-reduction-only ring)))
|
|---|
| 1360 | (cond
|
|---|
| 1361 | ((poly-with-sugar-zerop h))
|
|---|
| 1362 | (t
|
|---|
| 1363 | (setf h (grobner-primitive-part-with-sugar h ring))
|
|---|
| 1364 | (multiple-value-setq (G B)
|
|---|
| 1365 | (update-with-sugar G B h pred ring))
|
|---|
| 1366 | #+debug (debug-cgb "~&[[~d]] Polynomials: ~d; Pairs left: ~d~%"
|
|---|
| 1367 | (poly-with-sugar-sugar h)
|
|---|
| 1368 | (length G) (length B))
|
|---|
| 1369 | ))))
|
|---|
| 1370 | (setf G (mapcar #'poly-with-sugar-poly G))
|
|---|
| 1371 | ;;#+grobner-check(grobner-test G F pred ring)
|
|---|
| 1372 | #+debug(debug-cgb "~&GROBNER END")
|
|---|
| 1373 | G)
|
|---|
| 1374 |
|
|---|
| 1375 | (defun update-with-sugar (G B h pred ring &aux C D E B-new G-new pair)
|
|---|
| 1376 | "An implementation of the auxillary UPDATE algorithm used by the
|
|---|
| 1377 | Gebauer-Moeller algorithm. G is a list of polynomials, B is a list of
|
|---|
| 1378 | critical pairs and H is a new polynomial which possibly will be added
|
|---|
| 1379 | to G. The naming conventions used are very close to the one used in
|
|---|
| 1380 | the book of Becker-Weispfenning. Operates on polynomials with sugar."
|
|---|
| 1381 | (setf C G D nil)
|
|---|
| 1382 | (do (g1) ((endp C))
|
|---|
| 1383 | (setf g1 (pop C))
|
|---|
| 1384 | (when (or (monom-rel-prime (poly-with-sugar-lm h) (poly-with-sugar-lm g1))
|
|---|
| 1385 | (and
|
|---|
| 1386 | (not (some
|
|---|
| 1387 | #'(lambda (g2) (monom-divides-p (monom-lcm (poly-with-sugar-lm h)
|
|---|
| 1388 | (poly-with-sugar-lm g2))
|
|---|
| 1389 | (monom-lcm (poly-with-sugar-lm h)
|
|---|
| 1390 | (poly-with-sugar-lm g1))))
|
|---|
| 1391 | C))
|
|---|
| 1392 | (not (some
|
|---|
| 1393 | #'(lambda (g2) (monom-divides-p (monom-lcm (poly-with-sugar-lm h)
|
|---|
| 1394 | (poly-with-sugar-lm g2))
|
|---|
| 1395 | (monom-lcm (poly-with-sugar-lm h)
|
|---|
| 1396 | (poly-with-sugar-lm g1))))
|
|---|
| 1397 | D))))
|
|---|
| 1398 | (push g1 D)))
|
|---|
| 1399 | (setf E nil)
|
|---|
| 1400 | (do (g1) ((endp D))
|
|---|
| 1401 | (setf g1 (pop D))
|
|---|
| 1402 | (unless (monom-rel-prime (poly-with-sugar-lm h) (poly-with-sugar-lm g1))
|
|---|
| 1403 | (push g1 E)))
|
|---|
| 1404 | (setf B-new nil)
|
|---|
| 1405 | (do (g1 g2) ((endp B))
|
|---|
| 1406 | (setf pair (pop B)
|
|---|
| 1407 | g1 (first pair)
|
|---|
| 1408 | g2 (second pair))
|
|---|
| 1409 | (when (or (not (monom-divides-p (poly-with-sugar-lm h) (monom-lcm (poly-with-sugar-lm g1)
|
|---|
| 1410 | (poly-with-sugar-lm g2))))
|
|---|
| 1411 | (equal (monom-lcm (poly-with-sugar-lm g1) (poly-with-sugar-lm h))
|
|---|
| 1412 | (monom-lcm (poly-with-sugar-lm g1) (poly-with-sugar-lm g2)))
|
|---|
| 1413 | (equal (monom-lcm (poly-with-sugar-lm h) (poly-with-sugar-lm g2))
|
|---|
| 1414 | (monom-lcm (poly-with-sugar-lm g1) (poly-with-sugar-lm g2))))
|
|---|
| 1415 | (setf B-new (gebauer-moeller-with-sugar-merge-pairs B-new (list (list g1 g2)) pred ring))))
|
|---|
| 1416 | (dolist (g3 E)
|
|---|
| 1417 | (setf B-new (gebauer-moeller-with-sugar-merge-pairs B-new (list (list h g3)) pred ring)))
|
|---|
| 1418 | (setf G-new nil)
|
|---|
| 1419 | (do (g1) ((endp G))
|
|---|
| 1420 | (setf g1 (pop G))
|
|---|
| 1421 | (unless (monom-divides-p (poly-with-sugar-lm h) (poly-with-sugar-lm g1))
|
|---|
| 1422 | (push g1 G-new)))
|
|---|
| 1423 | (push h G-new)
|
|---|
| 1424 | (values G-new B-new))
|
|---|
| 1425 |
|
|---|
| 1426 | (defun gebauer-moeller-with-sugar-merge-pairs (B C pred ring)
|
|---|
| 1427 | "Merges lists of critical pairs. It orders pairs according to
|
|---|
| 1428 | increasing sugar, with ties broken by smaller lcm of head monomials In
|
|---|
| 1429 | this function B must already be sorted. Operates on polynomials with sugar."
|
|---|
| 1430 | (setf C (delete-if
|
|---|
| 1431 | #'(lambda (p) (and (endp (cdr (first p)))
|
|---|
| 1432 | (endp (cdr (second p)))))
|
|---|
| 1433 | C))
|
|---|
| 1434 | (labels ((pair-sugar (p)
|
|---|
| 1435 | (spoly-sugar (first p) (second p) ring))
|
|---|
| 1436 | (pair-mock-spoly (p)
|
|---|
| 1437 | (mock-spoly (poly-with-sugar-poly (first p))
|
|---|
| 1438 | (poly-with-sugar-poly (second p)) pred))
|
|---|
| 1439 | (pair-order (p q)
|
|---|
| 1440 | (let ((sugar-p (pair-sugar p))
|
|---|
| 1441 | (sugar-q (pair-sugar q)))
|
|---|
| 1442 | (cond
|
|---|
| 1443 | ((< sugar-p sugar-q) t)
|
|---|
| 1444 | ((> sugar-p sugar-q) nil)
|
|---|
| 1445 | (t
|
|---|
| 1446 | (funcall pred (pair-mock-spoly q) (pair-mock-spoly p)))))))
|
|---|
| 1447 | (merge 'list B (sort C #'pair-order) #'pair-order)))
|
|---|
| 1448 |
|
|---|
| 1449 | (defun grobner-primitive-part-with-sugar (h ring)
|
|---|
| 1450 | (cons (grobner-primitive-part (car h) ring) (cdr h)))
|
|---|
| 1451 |
|
|---|
| 1452 | ;;----------------------------------------------------------------
|
|---|
| 1453 | ;;
|
|---|
| 1454 | ;; An implementation of the algorithms using product caching
|
|---|
| 1455 | ;;
|
|---|
| 1456 | ;;----------------------------------------------------------------
|
|---|
| 1457 | (defun cached-normal-form (f fl pred top-reduction-only ring cache)
|
|---|
| 1458 | "This version of NORMAL-FORM uses a cache which stores products of monomials
|
|---|
| 1459 | and polynomials found in FL. The cache is a hash table storing pairs (G . G-HASH-TABLE)
|
|---|
| 1460 | where G is a polynomial in FL serving as a key. The hash table G-HASH-TABLE is a list
|
|---|
| 1461 | of pairs (MONOMIAL . PRODUCT-POLY) where PRODUCT-POLY is (POLY* MONOMIAL G)."
|
|---|
| 1462 |
|
|---|
| 1463 | (declare (optimize (speed 3) (safety 0)))
|
|---|
| 1464 | ;; Loop invariant: c*f=sum ai*fi+r+p
|
|---|
| 1465 | (do (r (c (ring-unit ring))
|
|---|
| 1466 | (division-count 0)
|
|---|
| 1467 | (p f))
|
|---|
| 1468 | ((or (endp p)
|
|---|
| 1469 | (endp fl)
|
|---|
| 1470 | (and top-reduction-only r))
|
|---|
| 1471 | #+debug(progn
|
|---|
| 1472 | (debug-cgb "~&~3T~d reduction~:p" division-count)
|
|---|
| 1473 | (when (endp r)
|
|---|
| 1474 | (debug-cgb " ---> 0")))
|
|---|
| 1475 | (values (append (nreverse r) p)
|
|---|
| 1476 | division-count c))
|
|---|
| 1477 | (declare (fixnum division-count) (integer c))
|
|---|
| 1478 | (let ((i (position (lm p) fl :test #'monom-divisible-by-p :key #'lm)))
|
|---|
| 1479 | (cond
|
|---|
| 1480 | (i ;division possible
|
|---|
| 1481 | (incf division-count)
|
|---|
| 1482 | (let* ((g (nth i fl))
|
|---|
| 1483 | (hash-tab (nth i cache))
|
|---|
| 1484 | (gcd (funcall (ring-gcd ring) (lc g) (lc p)))
|
|---|
| 1485 | (c1 (funcall (ring-/ ring) (lc g) gcd))
|
|---|
| 1486 | (c2 (funcall (ring-/ ring) (lc p) gcd))
|
|---|
| 1487 | (m (monom/ (lm p) (lm g))))
|
|---|
| 1488 | ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
|
|---|
| 1489 | (setf r (scalar-times-poly c1 r ring)
|
|---|
| 1490 | c (funcall (ring-* ring) c c1)
|
|---|
| 1491 | p (cached-grobner-op c2 c1 m p g pred ring hash-tab))))
|
|---|
| 1492 | (t ;no division possible
|
|---|
| 1493 | (setf r (cons (lt p) r) ;move lt(p) to remainder
|
|---|
| 1494 | p (rest p)))))))
|
|---|
| 1495 |
|
|---|
| 1496 | (defun cached-grobner-op (c1 c2 m A B pred ring cache)
|
|---|
| 1497 | "Returns the same result as (GROBNER-OP C1 C2 M A B PRED RING) but
|
|---|
| 1498 | it used variable CACHE which is described in CACHED-NORMAL-FORM for
|
|---|
| 1499 | efficient calculation of products of monomials with polynomials."
|
|---|
| 1500 | (poly- (scalar-times-poly c2 A)
|
|---|
| 1501 | (scalar-times-poly c1 (cached-monom-times-poly m B cache) ring)
|
|---|
| 1502 | pred ring))
|
|---|
| 1503 |
|
|---|
| 1504 | (defun cached-buchberger (F pred start top-reduction-only ring
|
|---|
| 1505 | &aux (s (1- (length F))) B M CACHE)
|
|---|
| 1506 | "It returns the same result as BUCHBERGER but it uses caching
|
|---|
| 1507 | of products of monomials with polynomials."
|
|---|
| 1508 | (declare (fixnum s))
|
|---|
| 1509 | (unless (plusp s) (return-from cached-buchberger F)) ;cut startup costs
|
|---|
| 1510 | #+debug (debug-cgb "~&GROBNER BASIS - CACHED BUCHBERGER ALGORITHM")
|
|---|
| 1511 | #+debug (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
|
|---|
| 1512 | #+grobner-check (when (plusp start)
|
|---|
| 1513 | (grobner-test (subseq F 0 start) (subseq F 0 start) pred ring))
|
|---|
| 1514 | (setf B (nconc (makelist (list i j) (i 0 (1- start)) (j start s))
|
|---|
| 1515 | (makelist (list i j) (i start (1- s)) (j (1+ i) s)))
|
|---|
| 1516 | M (make-hash-table :test #'equal)
|
|---|
| 1517 | B (buchberger-sort-pairs B F pred ring))
|
|---|
| 1518 | ;;Initialize cache
|
|---|
| 1519 | (dotimes (k (1+ s)) (push (make-hash-table :test #'equal) CACHE))
|
|---|
| 1520 | ;;Initialize treated pairs
|
|---|
| 1521 | (dotimes (i (1- start))
|
|---|
| 1522 | (do ((j (1+ i) (1+ j))) ((>= j start))
|
|---|
| 1523 | (setf (gethash (list i j) M) t)))
|
|---|
| 1524 | (do ((G F))
|
|---|
| 1525 | ((endp B)
|
|---|
| 1526 | #+grobner-check(grobner-test G F pred ring)
|
|---|
| 1527 | #+debug(debug-cgb "~&GROBNER END")
|
|---|
| 1528 | G)
|
|---|
| 1529 | (let ((pair (pop B)))
|
|---|
| 1530 | (unless (or (Criterion-1 pair G)
|
|---|
| 1531 | (Criterion-2 pair G M))
|
|---|
| 1532 | (let ((SP (cached-normal-form (cached-spoly (elt G (first pair))
|
|---|
| 1533 | (elt G (second pair))
|
|---|
| 1534 | pred
|
|---|
| 1535 | ring
|
|---|
| 1536 | (elt CACHE (first pair))
|
|---|
| 1537 | (elt CACHE (second pair)))
|
|---|
| 1538 | G pred top-reduction-only ring CACHE)))
|
|---|
| 1539 | (unless (poly-zerop SP) ;S-polynomial ---> 0;
|
|---|
| 1540 | (setf s (1+ s)
|
|---|
| 1541 | SP (grobner-primitive-part SP ring)
|
|---|
| 1542 | (cdr (last G)) (list SP)
|
|---|
| 1543 | (cdr (last Cache)) (list (make-hash-table :test #'equal))
|
|---|
| 1544 | B (funcall *buchberger-merge-pairs* B (makelist (list i s)
|
|---|
| 1545 | (i 0 (1- s)))
|
|---|
| 1546 | G pred ring))
|
|---|
| 1547 | #+debug (debug-cgb "~&Polynomials: ~d; Pairs left: ~d; Treated pairs: ~d"
|
|---|
| 1548 | (length G) (length B)
|
|---|
| 1549 | (hash-table-count M)))))
|
|---|
| 1550 | (setf (gethash pair M) t))))
|
|---|
| 1551 |
|
|---|
| 1552 | (defun cached-spoly (f g pred ring cache-f cache-g)
|
|---|
| 1553 | "Uses CACHE to fetch products of monomials and polynomials."
|
|---|
| 1554 | (declare (optimize (speed 3) (safety 0)))
|
|---|
| 1555 | (let* ((lcm (monom-lcm (lm f) (lm g)))
|
|---|
| 1556 | (m1 (monom/ lcm (lm f)))
|
|---|
| 1557 | (m2 (monom/ lcm (lm g))))
|
|---|
| 1558 | (let* ((c (funcall (ring-lcm ring) (lc f) (lc g)))
|
|---|
| 1559 | (cf (funcall (ring-/ ring) c (lc f)))
|
|---|
| 1560 | (cg (funcall (ring-/ ring) c (lc g))))
|
|---|
| 1561 | (poly-
|
|---|
| 1562 | (scalar-times-poly cf (cached-monom-times-poly m1 f cache-f) ring)
|
|---|
| 1563 | (scalar-times-poly cg (cached-monom-times-poly m2 g cache-g) ring)
|
|---|
| 1564 | pred
|
|---|
| 1565 | ring))))
|
|---|
| 1566 |
|
|---|
| 1567 |
|
|---|
| 1568 | (defun cached-monom-times-poly (m f cache)
|
|---|
| 1569 | (let ((prod (gethash m cache)))
|
|---|
| 1570 | (cond
|
|---|
| 1571 | (prod #+debug(debug-cgb "c") prod)
|
|---|
| 1572 | (t
|
|---|
| 1573 | #+debug(debug-cgb ".")
|
|---|
| 1574 | (setf (gethash m cache) (monom-times-poly m f))))))
|
|---|
| 1575 |
|
|---|
| 1576 | ;;------------------------------------------------------------------------------------------------
|
|---|
| 1577 | ;;
|
|---|
| 1578 | ;; An implementation of the sugar strategy with product caching
|
|---|
| 1579 | ;;
|
|---|
| 1580 | ;;------------------------------------------------------------------------------------------------
|
|---|
| 1581 | (defun cached-normal-form-with-sugar (f fl pred top-reduction-only ring cache)
|
|---|
| 1582 | "Normal form of the polynomial with sugar F with respect to
|
|---|
| 1583 | a list of polynomials with sugar FL. Assumes that FL is not empty.
|
|---|
| 1584 | Parameters PRED, TOP-REDUCTION-ONLY and RING have the same meaning as in NORMAL-FORM.
|
|---|
| 1585 | The parameter SUGAR-LIMIT should be set to a positive integer. If the sugar limit of
|
|---|
| 1586 | the partial remainder exceeds SUGAR-LIMIT then the calculation is stopped and the partial
|
|---|
| 1587 | remainder is returned, although it is not fully reduced with respect to FL."
|
|---|
| 1588 | (declare (optimize (speed 3) (safety 0)))
|
|---|
| 1589 | ;; Loop invariant: c*f=sum ai*fi+r+p
|
|---|
| 1590 | (do ((r (cons nil 0))
|
|---|
| 1591 | (c (ring-unit ring))
|
|---|
| 1592 | (division-count 0)
|
|---|
| 1593 | (p f))
|
|---|
| 1594 | ((or (poly-with-sugar-zerop p)
|
|---|
| 1595 | (and top-reduction-only (not (poly-with-sugar-zerop r))))
|
|---|
| 1596 | #+debug
|
|---|
| 1597 | (progn
|
|---|
| 1598 | (debug-cgb "~&~3T~d reduction~:p" division-count)
|
|---|
| 1599 | (cond
|
|---|
| 1600 | ((poly-with-sugar-zerop r)
|
|---|
| 1601 | (debug-cgb " ---> 0"))))
|
|---|
| 1602 | (values (poly-with-sugar-append (poly-with-sugar-nreverse r) p)
|
|---|
| 1603 | division-count c))
|
|---|
| 1604 | (declare (fixnum division-count) (integer c))
|
|---|
| 1605 | (let ((i (position (poly-with-sugar-lm p) fl
|
|---|
| 1606 | :test #'monom-divisible-by-p
|
|---|
| 1607 | :key #'poly-with-sugar-lm)))
|
|---|
| 1608 | (cond
|
|---|
| 1609 | (i ;division possible
|
|---|
| 1610 | (incf division-count)
|
|---|
| 1611 | (let* ((g (nth i fl))
|
|---|
| 1612 | (hash-tab (nth i cache))
|
|---|
| 1613 | (gcd (funcall (ring-gcd ring) (poly-with-sugar-lc g) (poly-with-sugar-lc p)))
|
|---|
| 1614 | (c1 (funcall (ring-/ ring) (poly-with-sugar-lc g) gcd))
|
|---|
| 1615 | (c2 (funcall (ring-/ ring) (poly-with-sugar-lc p) gcd))
|
|---|
| 1616 | (m (monom/ (poly-with-sugar-lm p) (poly-with-sugar-lm g))))
|
|---|
| 1617 | ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
|
|---|
| 1618 | (setf r (scalar-times-poly-with-sugar c1 r ring)
|
|---|
| 1619 | c (funcall (ring-* ring) c c1)
|
|---|
| 1620 | p (scalar-times-poly-with-sugar c1 p ring)
|
|---|
| 1621 | p (cached-poly-with-sugar-op p c2 m g pred ring hash-tab))))
|
|---|
| 1622 | (t ;no division possible
|
|---|
| 1623 | ;;move lt(p) to remainder
|
|---|
| 1624 | (setf r (cons (cons (poly-with-sugar-lt p) (poly-with-sugar-poly r))
|
|---|
| 1625 | (poly-with-sugar-sugar r))
|
|---|
| 1626 | p (poly-with-sugar-tail p)))))))
|
|---|
| 1627 |
|
|---|
| 1628 | ;; Parameter START is used when Grobner basis is calculated incrementally
|
|---|
| 1629 | ;; and it signals that the elements of F from 0 to START-1 form a Grobner
|
|---|
| 1630 | ;; basis already, so their pairs are not formed.
|
|---|
| 1631 | (defun cached-buchberger-with-sugar (F-no-sugar pred start top-reduction-only ring
|
|---|
| 1632 | &aux (s (1- (length F-no-sugar))) B M F Cache)
|
|---|
| 1633 | "Returns a Grobner basis of an ideal generated by a list of ordinary
|
|---|
| 1634 | polynomials F-NO-SUGAR. Adds initial sugar to the polynomials and
|
|---|
| 1635 | performs the Buchberger algorithm with ``sugar strategy''. It returns
|
|---|
| 1636 | an ordinary list of polynomials with no sugar. One of the most
|
|---|
| 1637 | effective algorithms for Grobner basis calculation."
|
|---|
| 1638 | (declare (fixnum s))
|
|---|
| 1639 | (unless (plusp s) (return-from cached-buchberger-with-sugar F-no-sugar)) ;cut startup costs
|
|---|
| 1640 | #+debug (debug-cgb "~&GROBNER BASIS - BUCHBERGER WITH SUGER STRATEGY ALGORITHM")
|
|---|
| 1641 | #+debug (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
|
|---|
| 1642 | #+grobner-check (when (plusp start)
|
|---|
| 1643 | (grobner-test (subseq F-no-sugar 0 start)
|
|---|
| 1644 | (subseq F-no-sugar 0 start) pred ring))
|
|---|
| 1645 | (setf F (mapcar #'(lambda (x) (poly-add-sugar x ring)) F-no-sugar)
|
|---|
| 1646 | B (nconc (makelist (list i j) (i 0 (1- start)) (j start s))
|
|---|
| 1647 | (makelist (list i j) (i start (1- s)) (j (1+ i) s)))
|
|---|
| 1648 | M (make-hash-table :test #'equal))
|
|---|
| 1649 | ;;Initialize cache
|
|---|
| 1650 | (dotimes (k (1+ s)) (push (make-hash-table :test #'equal) Cache))
|
|---|
| 1651 | ;;Initialize treated pairs
|
|---|
| 1652 | (dotimes (i (1- start))
|
|---|
| 1653 | (do ((j (1+ i) (1+ j))) ((>= j start))
|
|---|
| 1654 | (setf (gethash (list i j) M) t)))
|
|---|
| 1655 | (setf B (buchberger-with-sugar-sort-pairs B F M pred ring))
|
|---|
| 1656 | (do ((G F))
|
|---|
| 1657 | ((endp B)
|
|---|
| 1658 | (setf G (mapcar #'poly-with-sugar-poly G))
|
|---|
| 1659 | ;;#+grobner-check(grobner-test G F-no-sugar pred ring)
|
|---|
| 1660 | #+debug(debug-cgb "~&GROBNER END")
|
|---|
| 1661 | G)
|
|---|
| 1662 | (let ((pair (pop B)))
|
|---|
| 1663 | (cond
|
|---|
| 1664 | ((Criterion-1-with-sugar pair G))
|
|---|
| 1665 | ((Criterion-2-with-sugar pair G M))
|
|---|
| 1666 | (t
|
|---|
| 1667 | (let ((SP (cached-normal-form-with-sugar
|
|---|
| 1668 | (cons (cached-spoly (poly-with-sugar-poly (elt G (first pair)))
|
|---|
| 1669 | (poly-with-sugar-poly (elt G (second pair)))
|
|---|
| 1670 | pred ring
|
|---|
| 1671 | (elt Cache (first pair))
|
|---|
| 1672 | (elt Cache (second pair)))
|
|---|
| 1673 | (third pair))
|
|---|
| 1674 | G pred top-reduction-only ring Cache)))
|
|---|
| 1675 | (cond
|
|---|
| 1676 | ((poly-with-sugar-zerop SP))
|
|---|
| 1677 | (t
|
|---|
| 1678 | (setf s (1+ s)
|
|---|
| 1679 | (poly-with-sugar-poly SP) (grobner-primitive-part
|
|---|
| 1680 | (poly-with-sugar-poly SP) ring)
|
|---|
| 1681 | (cdr (last G)) (list SP)
|
|---|
| 1682 | (cdr (last Cache)) (list (make-hash-table :test #'equal))
|
|---|
| 1683 | B (buchberger-with-sugar-merge-pairs B (makelist (list i s)
|
|---|
| 1684 | (i 0 (1- s)))
|
|---|
| 1685 | G M pred ring))
|
|---|
| 1686 | #+debug (debug-cgb "~&[~d] Polynomials: ~d; Pairs left: ~d; Treated pairs: ~d"
|
|---|
| 1687 | (third pair)
|
|---|
| 1688 | (length G) (length B)
|
|---|
| 1689 | (hash-table-count M)))))))
|
|---|
| 1690 | (setf (cddr pair) nil
|
|---|
| 1691 | (gethash pair M) t))))
|
|---|
| 1692 |
|
|---|
| 1693 | ;; Implement f-term*g
|
|---|
| 1694 | (defun cached-poly-with-sugar-op (f c m g pred ring Cache)
|
|---|
| 1695 | (poly-with-sugar- f (scalar-times-poly-with-sugar c (cached-monom-times-poly-with-sugar m g Cache) ring) pred ring))
|
|---|
| 1696 |
|
|---|
| 1697 | (defun cached-monom-times-poly-with-sugar (m f Cache)
|
|---|
| 1698 | (cons
|
|---|
| 1699 | (cached-monom-times-poly m (poly-with-sugar-poly f) Cache)
|
|---|
| 1700 | (+ (poly-with-sugar-sugar f) (monom-sugar m))))
|
|---|
![(please configure the [header_logo] section in trac.ini)](/trac/CGBLisp/chrome/site/your_project_logo.png){kind=logo}
