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