source: CGBLisp/src/colored-poly.lisp@ 1

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

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