source: CGBLisp/trunk/src/colored-poly.lisp@ 20

Last change on this file since 20 was 20, checked in by Marek Rychlik, 15 years ago

* empty log message *
File size: 37.0 KB
Line 
1#|
2 $Id: colored-poly.lisp,v 1.12 2009/01/24 11:07:32 marek Exp $
3 *--------------------------------------------------------------------------*
4 | Copyright (C) 1994, Marek Rychlik (e-mail: rychlik@math.arizona.edu) |
5 | Department of Mathematics, University of Arizona, Tucson, AZ 85721 |
6 | |
7 | Everyone is permitted to copy, distribute and modify the code in this |
8 | directory, as long as this copyright note is preserved verbatim. |
9 *--------------------------------------------------------------------------*
10|#
11(defpackage "COLORED-POLY"
12 (:use "MONOM" "ORDER" "PARSE" "PRINTER" "MAKELIST"
13 "GROBNER" "DIVISION" "POLY" "COEFFICIENT-RING" "COMMON-LISP")
14 (:export make-colored-poly colored-poly-print colored-poly-print-list
15 colored-poly-to-poly
16 cond-print cond-system-print determine string-determine
17 color-poly color-poly-list
18 grobner-system totally-green-p cond-hm cond-normal-form
19 cond-spoly make-colored-poly-list
20 string-grobner-system string-cond string-cover
21 *colored-poly-debug*
22 parse-to-colored-poly-list
23 ))
24(in-package "COLORED-POLY")
25
26(proclaim '(optimize (speed 0) (debug 3)))
27
28#+debug(defvar *colored-poly-debug* nil "If true debugging output is on.")
29#+debug
30(defmacro debug-cgb (&rest args)
31 `(when *colored-poly-debug* (format *trace-output* ,@args)))
32
33(defun make-colored-poly (poly k &key (key #'identity)
34 (main-order #'lex>)
35 (parameter-order #'lex>)
36 &aux l)
37 "Colored poly is represented as a list
38 (TERM1 TERM2 ... TERMS)
39where each term is a triple
40 (MONOM . (POLY . COLOR))
41where monoms and polys have common number of variables while color is
42one of the three: :RED, :GREEN or :WHITE. This function translates an
43ordinary polynomial into a colored one by dividing variables into K
44and N-K, where N is the total number of variables in the polynomial
45poly; the function KEY can be called to select variables in arbitrary
46order."
47 (when (endp poly) (return-from make-colored-poly))
48 (setf l (length (caar poly)))
49 (labels
50 ((monom-split-variables (monom)
51 (values
52 (makelist (elt monom (funcall key i)) (i 0 (1- k)))
53 (makelist (elt monom (funcall key i)) (i k (1- l)))))
54 (term-split-variables (term)
55 (multiple-value-bind (main par)
56 (monom-split-variables (car term))
57 (cons main (cons par (cdr term)))))
58 (collect-terms (p)
59 (do ((p p (rest p))
60 (q (mapcar #'(lambda (x) (cons x nil))
61 (remove-duplicates (mapcar #'car p)
62 :test #'equal))))
63 ((endp p) q)
64 (push (cdar p) (cdr (assoc (caar p) q :test #'equal)))
65 )))
66 (mapcar #'(lambda (term)
67 (cons (car term)
68 (cons (sort (cdr term) parameter-order :key #'car)
69 :white)))
70 (collect-terms
71 (sort (mapcar #'term-split-variables poly)
72 main-order :key #'car)))))
73
74(defun make-colored-poly-list (plist &rest rest)
75 "Translate a list of polynomials PLIST into a list of colored polynomials
76by calling MAKE-COLORED-POLY. Returns the resulting list."
77 (mapcar #'(lambda (p) (apply #'make-colored-poly (cons p rest)))
78 plist))
79
80(defun color-poly-list (flist &optional (cond (list nil nil)))
81 "Add colors to an ordinary list of polynomials FLIST, according to a
82condition COND. A condition is a pair of polynomial lists. Each
83polynomial in COND is a polynomial in parameters only. The list
84(FIRST COND) is called the ``green list'' and it consists of
85polynomials which vanish for the parameters associated with the
86condition. The list (SECOND COND) is called the ``red list'' and it
87consists of the polynomials which do not vanish for the parameters
88associated with the condition."
89 (mapcar #'(lambda (f) (color-poly f cond)) flist))
90
91(defun color-poly (f &optional (cond (list nil nil)))
92 "Add color to a single polynomial F, according to condition COND.
93See the documentation of COLOR-POLY-LIST."
94 (mapcar
95 #'(lambda (term)
96 (cons (car term)
97 (cons (cadr term)
98 (cond
99 ((member (cadr term) (car cond)
100 :test #'tree-equal) :green)
101 ((member (cadr term) (cadr cond)
102 :test #'tree-equal) :red)
103 (t :white)))))
104 (cdr f)))
105
106;;Conversion to ordinary polynomial
107(defun colored-poly-to-poly (cpoly)
108 "For a given colored polynomial CPOLY, removes the colors and
109it returns the polynomial as an ordinary polynomial with
110coefficients which are polynomials in parameters."
111 (mapcan #'(lambda (term)
112 (mapcar #'(lambda (x) (cons (append (car term) (car x)) (cdr x)))
113 (cadr term)))
114 (cdr cpoly)))
115
116(defun colored-poly-print (poly vars
117 &key (stream t) (beg t)
118 (print-green-part nil)
119 (mark-coefficients nil))
120 "Print a colored polynomial POLY. Use variables VARS to represent
121the variables. Some of the variables are going to be used as
122parameters, according to the length of the monomials in the main
123monomial and coefficient part of each term in POLY. The key variable
124STREAM may be used to redirect the output. If parameter
125PRINT-GREEN-PART is set then the coefficients which have color :GREEN
126will be printed, otherwise they are discarded silently. If
127MARK-COEFFICIENTS is not NIL then every coefficient will be marked
128according to its color, for instance G(U-1) would mean that U-1 is
129in the green list. Returns P."
130 (when (null poly)
131 (when beg (format stream "0") )
132 (return-from colored-poly-print))
133 (let* ((term (car poly))
134 (color (cddr term))
135 (n (length (car term))))
136 (unless (and (not print-green-part) (eq color :green))
137 (unless beg (format stream " + "))
138 (if mark-coefficients
139 (format stream "~c("
140 (case color (:red #\R) (:green #\G) (:white #\W)))
141 (format stream "("))
142 (poly-print (cadr term) (subseq vars n))
143 (format stream ")")
144 (print-monom (car term) (subseq vars 0 n) stream)
145 (setf beg nil)))
146 (colored-poly-print (rest poly) vars :stream stream :beg beg
147 :print-green-part print-green-part
148 :mark-coefficients mark-coefficients))
149
150(defun colored-poly-print-list (poly-list vars
151 &key (stream t)
152 (beg t)
153 (print-green-part nil)
154 (mark-coefficients nil))
155 "Pring a list of colored polynomials via a call to
156COLORED-POLY-PRINT."
157 (when (endp poly-list)
158 (when beg (format stream "["))
159 (format stream " ]")
160 (return-from colored-poly-print-list))
161 (unless (or (endp (car poly-list))
162 (and (not print-green-part) (endp (cond-part (car poly-list)))))
163 (when beg (format stream "[ "))
164 (unless beg (format stream ", "))
165 (colored-poly-print (car poly-list) vars
166 :stream stream
167 :print-green-part print-green-part
168 :mark-coefficients mark-coefficients)
169 (setf beg nil))
170 (colored-poly-print-list (rest poly-list) vars
171 :stream stream :beg beg
172 :print-green-part print-green-part
173 :mark-coefficients mark-coefficients))
174
175(defun determine (F &optional (cond (list nil nil)) (order #'lex>) (ring *coefficient-ring*))
176 "This function takes a list of colored polynomials F and a condition
177COND, and it returns a list of pairs (COND' F') such that COND' cover
178COND and F' is a ``determined'' version of the colored polynomial list
179F, i.e. every polynomial has its leading coefficient determined. This
180means that some of the initial coefficients in each polynomial in F'
181are in the green list of COND, and the first non-green coefficient is
182in the red list of COND. We note that F' differs from F only by
183different colors: some of the terms marked :WHITE are now marked
184either :GREEN or :RED. Coloring is done either by explicitly checking
185membership in red or green list of COND, or implicitly by performing
186Grobner basis calculations in the polynomial ring over the
187parameters. The admissible monomial order ORDER is used only in the
188parameter space. Also, the ring structure RING is used only for
189calculations on polynomials of the parameters only."
190 (cond
191 ((endp F) (list (list cond nil)))
192 (t
193 (let ((gs (mapcan #'(lambda (p &aux (cond1 (car p)) (F1 (cadr p)))
194 (determine-1 cond1 (car F) F1 nil order ring))
195 (determine (rest F) cond order ring))))
196 #+debug(let ((l (length gs)))
197 (when (> l 1)
198 (debug-cgb "~&Determined ~d new case~:p." (1- l))))
199 gs))))
200
201(defun determine-1 (cond P end GP order ring)
202 "Determine a single colored polynomial P according to condition COND.
203Prepend green part GP to P. Cons the result with END, which should be
204a list of colored polynomials, and return the resulting list of polynomials.
205This is an auxillary function of DETERMINE."
206 (cond
207 ((endp P) (list (list cond (append (list GP) end))))
208 ((eq (cddar P) :green)
209 (determine-1 cond (rest P) end (append GP (list (car P))) order ring))
210 ((eq (cddar P) :red)
211 (list (list cond (append (list (append GP P)) end))))
212 ;; white cases
213 (t (determine-white-term cond (car P) (rest P) end GP order ring))))
214
215#-use-saturation
216(defun determine-white-term (cond term restP end GP order ring)
217 "This is an auxillary function of DETERMINE.
218In this function the parameter COND is a condition.
219The parameters TERM, RESTP and GP are three parts of a polynomial being
220processed, where TERM is colored :WHITE. By testing membership
221in the red and green list of COND we try to determine whether the term
222is :RED or :GREEN. If we are successful, we simply change the color
223of the term and return the list ((COND P)) where P is obtained by
224appending GP, (LIST TERM) and RESTP. If we cannot determine whether
225TERM is :RED or :GREEN, we return the list ((COND' P') (COND'' P ''))
226where COND' is obtained by adding the coefficient of TERM to
227the red list of COND and P' is obtained by appending
228GP, (LIST TERM) and RESTP. COND'' is obtained by putting
229the coefficient of TERM into the green list of COND and
230P'' is obtaind by a recursive call to DETERMINE-1 on
231RESTP, together with GP and TERM which was marked :GREEN."
232 (cond
233 ((member (cadr term) (car cond) :test #'tree-equal) ;green
234 (determine-1 cond restP end
235 (append GP (list (list* (car term) (cadr term) :green)))
236 order ring))
237 ((or (member (cadr term) (cadr cond) :test #'tree-equal) ;red
238 (poly-constant-p (cadr term)))
239 (list
240 (list cond
241 (append (list
242 (append GP ;green part
243 (list (list* (car term) (cadr term) :red));red term
244 restP)) ;other terms
245 end))))
246 (t ;white
247 (cons
248 (list (list (car cond) (cons (cadr term) (cadr cond)))
249 (append
250 (list
251 (append
252 GP
253 (list (cons (car term) (cons (cadr term) :red)))
254 restP))
255 end))
256 (determine-1
257 (list (cons (cadr term) (car cond)) (cadr cond))
258 restP
259 end
260 (append GP (list (cons (car term) (cons (cadr term) :green))))
261 order ring)))))
262
263;; In this version, push on red list literally and test membership
264;; But keep green list to be a Grobner basis
265#+use-saturation
266(defun determine-white-term (cond term restP end GP order ring)
267 "This is an auxillary function of DETERMINE. In this function the
268parameter COND is a condition. The parameters TERM, RESTP and GP are
269three parts of a polynomial being processed, where TERM is colored
270:WHITE. We test the membership in the red and green list of COND we
271try to determine whether the term is :RED or :GREEN. This is done by
272performing ideal membership tests in the polynomial ring. Let C be the
273coefficient of TERM. Thus, C is a polynomial in parameters. We find
274whether C is in the green list by performing a plain ideal membership
275test. However, to test properly whether C is in the red list, one
276needs a different strategy. In fact, we test whether adding C to the
277red list would produce a non-empty set of parameters in some algebraic
278extension. The test is whether 1 belongs to the saturation ideal of
279(FIRST COND) in (CONS C (SECOND COND)). Thus, we use POLY-SATURATION.
280If we are successful in determining the color of TERM, we simply
281change the color of the term and return the list ((COND P)) where P is
282obtained by appending GP, (LIST TERM) and RESTP. If we cannot
283determine whether TERM is :RED or :GREEN, we return the list ((COND'
284P') (COND'' P'')) where COND' is obtained by adding the coefficient of
285TERM to the red list of COND and P' is obtained by appending GP, (LIST
286TERM) and RESTP. COND'' is obtained by putting the coefficient of TERM
287into the green list of COND and P'' is obtaind by a recursive call to
288DETERMINE-1 on RESTP, together with GP and TERM which was marked
289:GREEN."
290 (if (member (cadr term) (cadr cond) :test #'tree-equal)
291 ;;Paint red and return
292 (list
293 (list
294 cond
295 (append
296 (list (append GP (list (list* (car term) (cadr term) :red)) restP))
297 end)))
298 (let ((green-sat (ideal-saturation-1
299 (car cond)
300 (cadr term) order (length (car cond))
301 nil ring)))
302 (if (some #'poly-constant-p green-sat) ; in the radical of green
303 ;; paint the term green and determine the rest
304 (determine-1
305 cond restP end
306 (append GP (list (list* (car term) (cadr term) :green))) ;green term
307 order ring)
308 ;; else it does not contradict green and thus may be added to red list
309 ;; so it should be added to either red or green list
310 (cons
311 ;; Add to red
312 (list
313 (list green-sat (append (cadr cond) (list (cadr term))))
314 (append
315 (list (append
316 GP
317 (list (list* (car term) (cadr term) :red))
318 restP))
319 end))
320 ;; Add to green
321 (let ((sat (ideal-polysaturation-1
322 (append (car cond) (list (cadr term)))
323 (cadr cond)
324 order (length (car cond))
325 nil ring)))
326 (unless (some #'poly-constant-p sat) ;contradiction after all
327 (determine-1
328 (list sat (cadr cond))
329 restP end
330 (append GP (list (list* (car term) (cadr term) :green)))
331 order ring))))))))
332
333
334
335;; Print a conditional system, i.e. a list of pairs (gamma colored-poly-list)
336(defun cond-system-print (system vars params
337 &key (suppress-value t)
338 (print-green-part nil)
339 (mark-coefficients nil)
340 &aux (label 0))
341 "A conditional system SYSTEM is a list of pairs (COND PLIST), where COND
342is a condition (a pair (GREEN-LIST RED-LIST)) and PLIST is a list
343of colored polynomials. This function pretty-prints this list of
344pairs. A conditional system is the data structure returned by
345GROBNER-SYSTEM. This function returns SYSTEM, if SUPPRESS-VALUE is non-NIL and
346no value otherwise. If MARK-COEFFICIENTS is non-NIL coefficients will be marked
347as in G(u-1)*x+R(2)*y, which means that u-1 is :GREEN and 2 is :RED."
348 (dolist (pair system (if suppress-value (values) system))
349 (let ((cond (car pair))
350 (basis (cadr pair)))
351 (format t "~&------------------- CASE ~d -------------------"
352 (incf label))
353 (cond-print cond params)
354 (format t "~&~1TBasis: ")
355 (colored-poly-print-list basis (append vars params)
356 :print-green-part print-green-part
357 :mark-coefficients mark-coefficients))))
358
359;; Print a condition
360(defun cond-print (cond params)
361 "Pretty-print a condition COND, using symbol list PARAMS as parameter names."
362 (format t "~&Condition:")
363 (format t "~&~2TGreen list: ")
364 (poly-print (cons '[ (first cond)) params)
365 (format t "~&~2TRed list: ")
366 (poly-print (cons '[ (second cond)) params))
367
368
369(defun add-pairs (gs pred)
370 "The parameter GS shoud be a Grobner system, i.e. a set of pairs
371(CONDITION POLY-LIST) This functions adds the third component: the
372list of initial critical pairs (I J), as in the ordinary Grobner basis
373algorithm. In addition, it adds the length of of the POLY-LIST, less
3741, as the fourth component. The resulting list of quadruples is
375returned."
376 #-reorder-pairs(declare (ignore pred))
377 (mapcar
378 #'(lambda (gb &aux (n (length (cadr gb))))
379 (let ((B (makelist (list i j) (i 0 (- n 2)) (j (1+ i) (1- n)))))
380 #+reorder-pairs
381 (setf B (reorder-pairs B nil (cadr gb) pred t))
382 (append gb (list B (1- n)))))
383 gs))
384
385(defun cond-part (p)
386 "Find the part of a colored polynomial P starting with the first
387 non-green term."
388 (member :green p :test-not #'eq :key #'cddr))
389
390(defun cond-hm (p)
391 "Return the conditional head monomial of a colored polynomial P."
392 (let ((cp (cond-part p)))
393 (cond
394 ((endp cp) (error "Zero conditional part."))
395 ((eq (cddar cp) :red) (car cp))
396 (t (error "Head not determined.")))))
397
398(defun delete-green-polys (gamma)
399 "Delete totally green polynomials from in a grobner system GAMMA."
400 (dolist (gb gamma gamma)
401 (setf (cadr gb) (delete-if-not #'cond-part (cadr gb)))))
402
403;; B is a cover (i.e. a list of conditions)
404;; flist is a list of colored polynomials
405(defun grobner-system (F &key
406 (cover (list '(nil nil)))
407 (main-order #'lex>)
408 (parameter-order #'lex>)
409 (reduce t)
410 (green-reduce t)
411 (top-reduction-only nil)
412 (ring *coefficient-ring*)
413 &aux
414 (cover #-use-saturation cover
415 #+use-saturation
416 (saturate-cover cover parameter-order ring))
417 (gamma (delete-green-polys
418 (mapcan #'(lambda (cond)
419 (determine F cond
420 parameter-order ring))
421 cover))))
422 "This function returns a grobner system, given a list of colored polynomials F,
423Other parameters are:
424A cover COVER, i.e. a list of conditions, i.e. pairs of the form (GREEN-LIST RED-LIST), where
425GREEN-LIST and RED-LIST are to lists of ordinary polynomials in parameters.
426A monomial order MAIN-ORDER used on main variables (not parameters).
427A monomial order PARAMETER-ORDER used in calculations with parameters only.
428REDUCE, a flag deciding whether COLORED-REDUCTION will be performed on the resulting
429grobner system.
430GREEN-REDUCE, a flag deciding whether the green list of each condition will be reduced in
431a form of a reduced Grobner basis.
432TOP-REDUCTION-ONLY, a flag deciding whether in the internal calculations in the space of parameters
433top reduction only will be used.
434RING, a structure as in the package COEFFICIENT-RING, used in operations on the coefficients
435of the polynomials in parameters."
436 #+debug(debug-cgb "~&Initially ~d open case~:p." (length gamma))
437 (do ((open (add-pairs gamma main-order))
438 closed)
439 ((endp open)
440 ;; Post-process Grobner system
441 (tidy-grobner-system
442 (mapcar #'(lambda (gp) (butlast gp 2)) closed)
443 main-order parameter-order reduce green-reduce ring))
444 #+debug(debug-cgb "~&Currently ~d open case~:p." (length open))
445 (let* ((gb (pop open)) (cond (car gb)) (G (cadr gb)) (B (caddr gb))
446 (s (cadddr gb)))
447 (declare (fixnum s))
448 (assert (= (length G) (1+ s)))
449 #+debug(debug-cgb "~&Colored case of ~d polynomials and ~d pairs."
450 (1+ s) (length B))
451 (cond
452 ((endp B) ;no more pairs in this tuple
453 (push gb closed))
454 ((let* ((pair (car B))
455 (i (car pair))
456 (j (cadr pair)))
457 (declare (fixnum i j))
458 (or
459 ;;Buchberger criterion 1 or 2 succeeds
460 (colored-Criterion-1 i j G)
461 (colored-Criterion-2 i j G (rest B) s)))
462 (push (list cond G (rest B) s) open))
463 (t ;Grobner step - S-polynomial
464 (do* ((pair (car B))
465 (i (car pair))
466 (j (cadr pair))
467 (h (cond-spoly (elt G i) (elt G j)
468 main-order parameter-order ring))
469 (SP (cond-normal-form h G main-order parameter-order top-reduction-only ring))
470 (delta (determine (list SP) cond parameter-order ring)
471 (rest delta)))
472 ((endp delta))
473 (declare (fixnum i j))
474 (let ((cond1 (caar delta))
475 (SP1 (caadar delta)))
476 (cond
477 ((cond-part SP1) ;SP1 is not green
478 (let* ((G1 (append G (list SP1)))
479 (s1 (1+ s))
480 (Bnew (makelist (list k s1) (k 0 (1- s1))))
481 B1)
482 (assert (= (length G1) (1+ s1)))
483 #+reorder-pairs
484 (setf B1 (reorder-pairs (rest B) Bnew G1 main-order nil))
485 #-reorder-pairs
486 (setf B1 (append (rest B) Bnew))
487 (push (list cond1 G1 B1 s1) open)))
488 (t ;SP1 is totally green
489 (assert (= (length G) (1+ s)))
490 (push (list cond1 G (rest B) s) open))))))))))
491
492;; This is destructive to B and Bnew
493#+reorder-pairs
494(defun reorder-pairs (B Bnew G pred &optional (sort-first nil))
495 "Reorder pairs according to some heuristic. The heuristic at this time is ad hoc,
496in the future it should be replaced with sugar strategy and a mechanism for implementing
497new heuristic strategies, as in the GROBNER package."
498 (let ((order
499 #'(lambda (p1 p2)
500 (let* ((m1 (monom-lcm (cond-lm (elt G (car p1)))
501 (cond-lm (elt G (cadr p1)))))
502 (m2 (monom-lcm (cond-lm (elt G (car p2)))
503 (cond-lm (elt G (cadr p2)))))
504 (d1 (total-degree m1))
505 (d2 (total-degree m2)))
506 (cond ((< d1 d2) t)
507 ((= d1 d2) (funcall pred m2 m1))
508 (t nil))))))
509 (when sort-first (setf B (sort (copy-list B) order)))
510 (if Bnew
511 (setf B (merge 'list (sort Bnew order) (copy-list B) order))
512 B)))
513
514(defun colored-Criterion-1 (i j F)
515 "Buchberger criterion 1 for colored polynomials."
516 (declare (fixnum i j))
517 (let ((v (monom-rel-prime (cond-lm (elt F i))
518 (cond-lm (elt F j)))))
519 #+debug(when v (debug-cgb "~&~2TColored Buchberger1 succeded."))
520 v))
521
522(defun colored-Criterion-2 (i j F B s)
523 "Buchberger criterion 2 for colored polynomials."
524 (declare (fixnum i j s))
525 (labels ((pair (i j)
526 (declare (fixnum i j))
527 (if (< i j) (list i j) (list j i))))
528 (do ((k 1 (1+ k)))
529 ((>= k s))
530 (when (and (/= k i)
531 (/= k j)
532 (not (member (pair i k) B :test #'equal))
533 (not (member (pair j k) B :test #'equal))
534 (monom-divides-p (cond-lm (elt F k))
535 (monom-lcm (cond-lm (elt F i))
536 (cond-lm (elt F j)))))
537 #+debug(debug-cgb "~&~2TColored Buchberger2 succeded.")
538 (return-from colored-Criterion-2 t)))))
539
540
541(defun cond-normal-form (f fl main-order parameter-order top-reduction-only ring)
542 "Returns the conditional normal form of a colored polynomial F with respect to
543the list of colored polynomials FL. The list FL is assumed to consist of determined
544polynomials, i.e. such that the first term which is not marked :GREEN is
545:RED."
546 ;; Remove all zero (i.e totally green) polys from plist
547 ;; (setf fl (remove nil fl :key #'cond-part))
548 (do (r (division-count 0)
549 (p f))
550 ((or (endp p) (and top-reduction-only r))
551 #+debug(debug-cgb "~&~3T~d conditional reductions" division-count)
552 #+debug(when (endp r) (debug-cgb " ---> 0"))
553 (values (reverse r) division-count))
554 (cond
555 ((eq (cddar p) :green)
556 (setf r (cons (car p) r)
557 p (rest p)))
558 (t
559 ;; Find a divisor
560 (do ((fl fl (rest fl))) ;scan list of divisors
561 ((cond
562 ((endp fl)
563 ;; no division occurred
564 (setf r (cons (car p) r) ;move term to remainder
565 p (rest p)) ;remove car from p
566 t)
567 ((monom-divides-p (cond-lm (car fl))
568 (caar p))
569 (incf division-count)
570 (let* (#-colored-poly-use-grobner
571 (c1 (cons (make-list (length (caaar fl))
572 :initial-element 0)
573 (cond-lc (car fl))))
574 #+colored-poly-use-grobner
575 (lcm (poly-lcm (car (cond-lc (car fl))) (cadar p)
576 parameter-order ring))
577 #+colored-poly-use-grobner
578 (c1 (cons (make-list (length (cond-lm (car fl)))
579 :initial-element 0)
580 (cons (poly-exact-divide lcm (cadar p) parameter-order ring) :red)))
581 #+colored-poly-use-grobner
582 (c2 (poly-exact-divide
583 lcm
584 (car (cond-lc (car fl)))
585 parameter-order ring))
586 #-colored-poly-use-grobner
587 (c2 (cadar p))
588 ;; This works for both
589 (quot (cons (monom/ (caar p) (cond-lm (car fl)))
590 (cons c2 (cddr p)))))
591 ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
592 (setf r (colored-term-times-poly c1 r parameter-order ring)
593 p (colored-poly-
594 (colored-term-times-poly c1 p parameter-order ring)
595 (colored-term-times-poly quot (car fl) parameter-order ring)
596 main-order parameter-order ring)))
597 t))))))))
598
599
600
601(defun cond-spoly (f g main-order parameter-order ring)
602 "Returns the conditional S-polynomial of two colored polynomials F and G.
603Both polynomials are assumed to be determined."
604 (let* ((lcm (monom-lcm (cond-lm f) (cond-lm g)))
605 (m1 (monom/ lcm (cond-lm f)))
606 (m2 (monom/ lcm (cond-lm g))))
607 #-colored-poly-use-grobner
608 (colored-poly-
609 (colored-term-times-poly (cons m1 (cond-lc g)) (rest f)
610 parameter-order ring)
611 (colored-term-times-poly (cons m2 (cond-lc f)) (rest g)
612 parameter-order ring)
613 main-order parameter-order ring)
614
615 #+colored-poly-use-grobner
616 (let* ((lcm-2 (poly-lcm (car (cond-lc f)) (car (cond-lc g))
617 parameter-order ring))
618 (lcf (cond-lc f))
619 (lcg (cond-lc g))
620 (cf (cons (poly-exact-divide lcm-2 (car lcf) parameter-order ring)
621 :red))
622 (cg (cons (poly-exact-divide lcm-2 (car lcg) parameter-order ring)
623 :red)))
624 (colored-poly-
625 (colored-term-times-poly (cons m1 cf) f parameter-order ring)
626 (colored-term-times-poly (cons m2 cg) g parameter-order ring)
627 main-order parameter-order ring))
628 ))
629
630(defun cond-lm (f)
631 "Returns the conditional leading monomial of a colored polynomial F,
632which is assumed to be determined."
633 (car (cond-hm f)))
634
635;; Conditional leading coefficient; (poly . color)
636(defun cond-lc (f)
637 "Returns the conditional leading coefficient of a colored polynomial F,
638which is assumed to be determined."
639 (cdr (cond-hm f)))
640
641
642(defun colored-term-times-poly (term f order ring)
643 "Returns the product of a colored term TERM and a colored polynomial F."
644 (mapcar #'(lambda (x) (colored-term* term x order ring)) f))
645
646
647(defun colored-scalar-times-poly (c f ring)
648 "Returns the product of an element of the coefficient ring C a colored polynomial F."
649 (mapcar #'(lambda (x)
650 (cons (car x)
651 (cons (scalar-times-poly c (cadr x) ring)
652 (cddr x))))
653 f))
654
655(defun colored-term* (term1 term2 order ring)
656 "Returns the product of two colored terms TERM1 and TERM2."
657 (cons
658 (monom* (car term1) (car term2))
659 (cons (poly* (cadr term1) (cadr term2) order ring)
660 (color* (cddr term1) (cddr term2)))))
661
662(defun color* (c1 c2)
663 "Returns a product of two colores. Rules:
664:red * :red yields :red
665any * :green yields :green
666otherwise the result is :white."
667 (cond
668 ((and (eq c1 :red) (eq c2 :red)) :red)
669 ((or (eq c1 :green) (eq c2 :green)) :green)
670 (t :white)))
671
672(defun color+ (c1 c2)
673 "Returns a sum of colors. Rules:
674:green + :green yields :green,
675:red + :green yields :red
676any other result is :white."
677 (cond
678 ((and (eq c1 :green) (eq c2 :green)) :green)
679 ((and (eq c1 :red) (eq c2 :green)) :red)
680 ((and (eq c2 :red) (eq c1 :green)) :red)
681 (t :white)))
682
683(defun color- (c1 c2)
684 "Identical to COLOR+."
685 (color+ c1 c2))
686
687(defun colored-poly+ (p q main-order parameter-order ring)
688 "Returns the sum of colored polynomials P and Q."
689 (cond
690 ((endp p) q)
691 ((endp q) p)
692 (t
693 (multiple-value-bind
694 (mgreater mequal)
695 (funcall main-order (caar p) (caar q))
696 (cond
697 (mequal
698 (let ((s (poly+ (cadar p) (cadar q) parameter-order ring)))
699 (if (endp s) ;check for cancellation
700 (colored-poly+ (cdr p) (cdr q) main-order parameter-order ring)
701 (cons (cons (caar p) (cons s (color+ (cddar p) (cddar q))))
702 (colored-poly+ (cdr p) (cdr q) main-order
703 parameter-order ring)))))
704 (mgreater
705 (cons (car p)
706 (colored-poly+ (cdr p) q main-order parameter-order ring)))
707 (t (cons (cons (caar q) (cons (cadar q) (cddar q)))
708 (colored-poly+ p (cdr q) main-order parameter-order ring))))))))
709
710(defun colored-poly- (p q main-order parameter-order ring)
711 "Returns the difference of colored polynomials P and Q."
712 (do (r)
713 (nil)
714 (cond
715 ((endp p) (return (nreconc r (colored-minus-poly q ring))))
716 ((endp q) (return (nreconc r p)))
717 (t
718 (multiple-value-bind
719 (mgreater mequal)
720 (funcall main-order (caar p) (caar q))
721 (cond
722 (mequal
723 (let ((s (poly- (cadar p) (cadar q) parameter-order ring)))
724 (unless (endp s) ;check for cancellation
725 (setf r (cons (cons (caar p) (cons s (color- (cddar p) (cddar q)))) r)))
726 (setf p (cdr p) q (cdr q))))
727 (mgreater
728 (setf r (cons (car p) r)
729 p (cdr p)))
730 (t (setf r (cons (cons (caar q) (cons (minus-poly (cadar q) ring) (cddar q))) r)
731 q (cdr q)))))))))
732
733(defun colored-term-uminus (term ring)
734 "Returns the negation of a colored term TERM."
735 (cons (car term) (cons (minus-poly (cadr term) ring) (cddr term))))
736
737(defun colored-minus-poly (p ring)
738 "Returns the negation of a colored polynomial P."
739 (mapcar #'(lambda (x) (colored-term-uminus x ring)) p))
740
741(defun string-grobner-system (F vars params
742 &key (cover (list (list "[]" "[]")))
743 (main-order #'lex>)
744 (parameter-order #'lex>)
745 (ring *coefficient-ring*)
746 (suppress-value t)
747 (suppress-printing nil)
748 (mark-coefficients nil)
749 (reduce t)
750 (green-reduce t)
751 &aux
752 (F (parse-to-colored-poly-list
753 F vars params main-order parameter-order))
754 (cover (string-cover cover params parameter-order)))
755 "An interface to GROBNER-SYSTEM in which polynomials can be specified in infix notations
756as strings. Lists of polynomials are comma-separated list marked by a matchfix operators []"
757 (let ((gs (grobner-system F :cover cover :main-order main-order
758 :parameter-order parameter-order
759 :ring ring
760 :green-reduce green-reduce
761 :reduce reduce)))
762 (unless suppress-printing
763 (cond-system-print gs vars params :mark-coefficients mark-coefficients))
764 (if suppress-value (values) gs)))
765
766
767(defun string-cond (cond params &optional (order #'lex>))
768 "Return the internal representation of a condition COND, specified
769as pairs of strings (GREEN-LIST RED-LIST). GREEN-LIST and RED-LIST in
770the input are assumed to be strings which parse to two lists of
771polynomials with respect to variables whose names are in the list of
772symbols PARAMS. ORDER is the predicate used to sort the terms of
773the polynomials."
774 (list
775 (rest (parse-string-to-sorted-alist (car cond) params order))
776 (rest (parse-string-to-sorted-alist (cadr cond) params order))))
777
778(defun string-cover (cover params &optional (order #'lex>))
779 "Returns the internal representation of COVER, given in the form of
780a list of conditions. See STRING-COND for description of a condition."
781 (cond
782 ((endp cover) nil)
783 (t (cons (string-cond (car cover) params order)
784 (string-cover (cdr cover) params order)))))
785
786(defun saturate-cover (cover order ring)
787 "Brings every condition of a list of conditions COVER to the form (G R)
788where G is saturated with respect to R and G is a Grobner basis
789We could reduce R so that the elements of R are relatively prime,
790but this is not currently done."
791 (remove nil (mapcar #'(lambda (cond) (saturate-cond cond order ring)) cover)))
792
793(defun saturate-cond (cond order ring)
794 "Saturate a single condition COND. An auxillary function of SATURATE-COVER."
795 (let* ((green-sat (ideal-polysaturation-1 (car cond) (cadr cond) order 0 nil ring)))
796 (if (some #'poly-constant-p green-sat)
797 nil
798 (list green-sat (cadr cond)))))
799
800(defun string-determine (F vars params
801 &key (cond '("[]" "[]"))
802 (main-order #'lex>)
803 (parameter-order #'lex>)
804 (suppress-value t)
805 (suppress-printing nil)
806 (mark-coefficients nil)
807 (ring *coefficient-ring*)
808 &aux
809 (F (parse-to-colored-poly-list
810 F vars params main-order parameter-order))
811 (cond (string-cond cond params parameter-order)))
812 "A string interface to DETERMINE. See the documentation of STRING-GROBNER-SYSTEM."
813 (let ((gs (determine F cond parameter-order ring)))
814 (unless suppress-printing
815 (cond-system-print gs vars params :mark-coefficients mark-coefficients))
816 (if suppress-value (values) gs)))
817
818(defun tidy-grobner-system (gs main-order parameter-order
819 reduce green-reduce ring)
820 "Apply TIDY-PAIR to every pair of a Grobner system."
821 (mapcan #'(lambda (p) (tidy-pair p main-order parameter-order reduce green-reduce
822 ring))
823 gs))
824
825
826(defun tidy-pair (pair main-order parameter-order
827 reduce green-reduce ring
828 &aux gs)
829 "Make the output of Grobner system more readable by performing
830certain simplifications on an element of a Grobner system.
831If REDUCE is non-NIL then COLORED-reduction will be performed.
832In addition TIDY-COND is called on the condition part of the pair PAIR."
833 (if reduce
834 (setf gs (colored-reduction (car pair) (cadr pair)
835 main-order parameter-order ring))
836 (setf gs (list pair)))
837 (setf gs
838 (mapcar #'(lambda (pair)
839 (list (tidy-cond (car pair) parameter-order
840 ring)
841 (cadr pair)))
842 gs))
843 (when green-reduce
844 (setf gs (cond-system-green-reduce gs parameter-order ring)))
845 gs)
846
847(defun tidy-cond (cond order ring)
848 "Currently saturates condition COND and does RED-REDUCTION on the red list."
849 (let ((cond1 (saturate-cond cond order ring)))
850 (list (reduction (car cond1) order ring)
851 (red-reduction (cadr cond1) order ring))))
852
853(defun colored-reduction (cond P main-order parameter-order ring
854 &aux (open (list (list cond nil P))) closed)
855 "Reduce a list of colored polynomials P. The difficulty as compared
856to the usual Buchberger algorithm is that the polys may have the same
857leading monomial which may result in cancellations and polynomials
858which may not be determined. Thus, when we find those, we will have to
859split the condition by calling determine. Returns a list of pairs
860(COND' P') where P' is a reduced grobner basis with respect to any
861parameter choice compatible with condition COND'. Moreover, COND' form
862a cover of COND."
863 ;;We form a list of tripples (CONDITION G U) where G is reduced and U
864 ;;is the unreduced part of the basis and operate on these
865 (do ()
866 ((endp open) closed)
867 (let* ((tuple (pop open))
868 (cond1 (car tuple))
869 (G (cadr tuple))
870 (U (caddr tuple)))
871 (cond
872 ((endp U) ;no-more undetermined
873 (push (list cond1 G)
874 closed))
875 (t
876 (let ((f (car U)))
877 (multiple-value-bind (k div-count)
878 (cond-normal-form f (append G (cdr U))
879 main-order parameter-order nil ring)
880 (cond
881 ((zerop div-count)
882 (push (list cond1 (append G (list f)) (cdr U)) open))
883 (t
884 (do ((delta (determine (list k) cond1 parameter-order ring)
885 (rest delta)))
886 ((endp delta))
887 (let* ((eps (caar delta))
888 (k1 (caadar delta)))
889 (cond
890 ((cond-part k1) ;does not reduce to 0
891 ;; Replace f with k1 and start all over
892 (push (list eps nil (append G (list k1) (cdr U)))
893 open))
894 (t
895 ;; f reduces to 0 so just drop f
896 (push (list eps G (cdr U))
897 open))))))))))))))
898
899
900(defun green-reduce-colored-poly (cond f parameter-order ring)
901 "It takes a colored polynomial F and it returns a modified
902polynomial obtained by reducing coefficient of F modulo green list of
903the condition COND."
904 (dotimes (i (length f) f)
905 (multiple-value-bind (nf division-count c)
906 (normal-form (cadr (nth i f)) (car cond) parameter-order nil ring)
907 (declare (ignore division-count))
908 (unless (endp nf)
909 (setf f (colored-scalar-times-poly c f ring)))
910 (setf (cadr (nth i f)) nf))))
911
912
913(defun green-reduce-colored-list (cond fl parameter-order ring)
914 "Apply GREEN-REDUCE-COLORED-POLY to a list of polynomials FL."
915 (remove-if
916 #'endp
917 (cond
918 ((endp fl) nil)
919 (t
920 (cons
921 (green-reduce-colored-poly cond (car fl) parameter-order ring)
922 (green-reduce-colored-list cond (rest fl) parameter-order ring))))))
923
924(defun cond-system-green-reduce (gs parameter-order ring)
925 "Apply GREEN-REDUCE-COLORED-LIST to every pair of
926a grobner system GS."
927 (cond
928 ((endp gs) nil)
929 (t (cons
930 (list (caar gs) (green-reduce-colored-list
931 (caar gs) (cadar gs) parameter-order ring))
932 (cond-system-green-reduce (rest gs) parameter-order ring)))))
933
934
935(defun parse-to-colored-poly-list (F vars params main-order parameter-order
936 &aux (k (length vars))
937 (vars-params (append vars params)))
938 "Parse a list of polynomials F, given as a string, with respect to
939a list of variables VARS, given as a list of symbols, to the internal
940representation of a colored polynomial. The polynomials will be properly
941sorted by MAIN-ORDER, with the coefficients, which are polynomials in
942parameters, sorted by PARAMETER-ORDER. Both orders must be admissible
943monomial orders. This form is suitable for parsing polynomials with integer
944coefficients."
945 (make-colored-poly-list
946 (rest (parse-string-to-alist F vars-params))
947 k :main-order main-order :parameter-order parameter-order))
948(defun red-reduction (P pred ring
949 &aux
950 (P (remove-if #'poly-constant-p P)))
951 "Takes a family of polynomials and produce a list whose prime factors
952are the same but they are relatively prime
953Repetitively used the following procedure: it finds two elements f, g of
954P which are not relatively prime; it replaces f and g with f/GCD(f,g),
955 g/ GCD(f,f) and GCD(f,g)."
956 (when (endp P) (return-from red-reduction))
957 (do ((found t))
958 ((or (endp (cdr P)) (not found))
959 (mapcar #'(lambda (x) (grobner-primitive-part x ring)) P))
960 (setf found nil)
961 (tagbody
962 (do ((Q1 P (rest Q1)))
963 ((endp Q1))
964 (do ((Q2 (rest Q1) (rest Q2)))
965 ((endp Q2))
966 (let* ((f (car Q1))
967 (g (car Q2))
968 (h (grobner-gcd f g pred ring)))
969 (unless (poly-constant-p h)
970 (setf found t
971 P (remove f P)
972 P (remove G P)
973 P (cons h P))
974 (let ((f1 (poly-exact-divide f h pred ring))
975 (g1 (poly-exact-divide g h pred ring)))
976 (unless (poly-constant-p f1) (push f1 P))
977 (unless (poly-constant-p g1) (push g1 P)))
978 (go found)))))
979 found)))
Note: See TracBrowser for help on using the repository browser.