close Warning: Can't synchronize with repository "(default)" (The repository directory has changed, you should resynchronize the repository with: trac-admin $ENV repository resync '(default)'). Look in the Trac log for more information.

source: branches/f4grobner/grobner.lisp@ 60

Last change on this file since 60 was 56, checked in by Marek Rychlik, 9 years ago

* empty log message *

File size: 57.0 KB
RevLine 
[1]1;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*-
2;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3;;;
[4]4;;; Copyright (C) 1999, 2002, 2009 Marek Rychlik <rychlik@u.arizona.edu>
[1]5;;;
6;;; This program is free software; you can redistribute it and/or modify
7;;; it under the terms of the GNU General Public License as published by
8;;; the Free Software Foundation; either version 2 of the License, or
9;;; (at your option) any later version.
10;;;
11;;; This program is distributed in the hope that it will be useful,
12;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
13;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14;;; GNU General Public License for more details.
15;;;
16;;; You should have received a copy of the GNU General Public License
17;;; along with this program; if not, write to the Free Software
18;;; Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
19;;;
20;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
21
22(in-package :maxima)
23
24(macsyma-module cgb-maxima)
25
26(eval-when
27 #+gcl (load eval)
28 #-gcl (:load-toplevel :execute)
29 (format t "~&Loading maxima-grobner ~a ~a~%"
[47]30 "$Revision: 2.0 $" "$Date: 2015/06/02 0:34:17 $"))
[1]31
32;;FUNCTS is loaded because it contains the definition of LCM
33($load "functs")
34
35
36
37
38
39
40;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
41;;
42;; Global switches
43;; (Can be used in Maxima just fine)
44;;
45;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
46
47(defmvar $poly_monomial_order '$lex
48 "This switch controls which monomial order is used in polynomial
49and Grobner basis calculations. If not set, LEX will be used")
50
51(defmvar $poly_coefficient_ring '$expression_ring
52 "This switch indicates the coefficient ring of the polynomials
53that will be used in grobner calculations. If not set, Maxima's
54general expression ring will be used. This variable may be set
55to RING_OF_INTEGERS if desired.")
56
57(defmvar $poly_primary_elimination_order nil
58 "Name of the default order for eliminated variables in elimination-based functions.
59If not set, LEX will be used.")
60
61(defmvar $poly_secondary_elimination_order nil
62 "Name of the default order for kept variables in elimination-based functions.
63If not set, LEX will be used.")
64
65(defmvar $poly_elimination_order nil
66 "Name of the default elimination order used in elimination calculations.
67If set, it overrides the settings in variables POLY_PRIMARY_ELIMINATION_ORDER
68and SECONDARY_ELIMINATION_ORDER. The user must ensure that this is a true
69elimination order valid for the number of eliminated variables.")
70
71(defmvar $poly_return_term_list nil
72 "If set to T, all functions in this package will return each polynomial as a
73list of terms in the current monomial order rather than a Maxima general expression.")
74
75(defmvar $poly_grobner_debug nil
76 "If set to TRUE, produce debugging and tracing output.")
77
78(defmvar $poly_grobner_algorithm '$buchberger
79 "The name of the algorithm used to find grobner bases.")
80
81(defmvar $poly_top_reduction_only nil
82 "If not FALSE, use top reduction only whenever possible.
83Top reduction means that division algorithm stops after the first reduction.")
84
85
86
87;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
88;;
89;; Coefficient ring operations
90;;
91;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
92;;
93;; These are ALL operations that are performed on the coefficients by
94;; the package, and thus the coefficient ring can be changed by merely
95;; redefining these operations.
96;;
97;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
98
99(defstruct (ring)
100 (parse #'identity :type function)
101 (unit #'identity :type function)
102 (zerop #'identity :type function)
103 (add #'identity :type function)
104 (sub #'identity :type function)
105 (uminus #'identity :type function)
106 (mul #'identity :type function)
107 (div #'identity :type function)
108 (lcm #'identity :type function)
109 (ezgcd #'identity :type function)
110 (gcd #'identity :type function))
111
112(defparameter *ring-of-integers*
113 (make-ring
114 :parse #'identity
115 :unit #'(lambda () 1)
116 :zerop #'zerop
117 :add #'+
118 :sub #'-
119 :uminus #'-
120 :mul #'*
121 :div #'/
122 :lcm #'lcm
123 :ezgcd #'(lambda (x y &aux (c (gcd x y))) (values c (/ x c) (/ y c)))
124 :gcd #'gcd)
125 "The ring of integers.")
126
127
128
129;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
130;;
131;; This is how we perform operations on coefficients
132;; using Maxima functions.
133;;
134;; Functions and macros dealing with internal representation structure
135;;
136;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
137
138
139
140
141;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
142;;
143;; Low-level polynomial arithmetic done on
144;; lists of terms
145;;
146;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
147
148(defmacro termlist-lt (p) `(car ,p))
149(defun termlist-lm (p) (term-monom (termlist-lt p)))
150(defun termlist-lc (p) (term-coeff (termlist-lt p)))
151
152(define-modify-macro scalar-mul (c) coeff-mul)
153
154(defun scalar-times-termlist (ring c p)
155 "Multiply scalar C by a polynomial P. This function works
156even if there are divisors of 0."
157 (mapcan
158 #'(lambda (term)
159 (let ((c1 (funcall (ring-mul ring) c (term-coeff term))))
160 (unless (funcall (ring-zerop ring) c1)
161 (list (make-term (term-monom term) c1)))))
162 p))
163
164
165(defun term-mul (ring term1 term2)
166 "Returns (LIST TERM) wheter TERM is the product of the terms TERM1 TERM2,
167or NIL when the product is 0. This definition takes care of divisors of 0
168in the coefficient ring."
169 (let ((c (funcall (ring-mul ring) (term-coeff term1) (term-coeff term2))))
170 (unless (funcall (ring-zerop ring) c)
171 (list (make-term (monom-mul (term-monom term1) (term-monom term2)) c)))))
172
173(defun term-times-termlist (ring term f)
174 (declare (type ring ring))
175 (mapcan #'(lambda (term-f) (term-mul ring term term-f)) f))
176
177(defun termlist-times-term (ring f term)
178 (mapcan #'(lambda (term-f) (term-mul ring term-f term)) f))
179
180(defun monom-times-term (m term)
181 (make-term (monom-mul m (term-monom term)) (term-coeff term)))
182
183(defun monom-times-termlist (m f)
184 (cond
185 ((null f) nil)
186 (t
187 (mapcar #'(lambda (x) (monom-times-term m x)) f))))
188
189(defun termlist-uminus (ring f)
190 (mapcar #'(lambda (x)
191 (make-term (term-monom x) (funcall (ring-uminus ring) (term-coeff x))))
192 f))
193
194(defun termlist-add (ring p q)
195 (declare (type list p q))
196 (do (r)
197 ((cond
198 ((endp p)
199 (setf r (revappend r q)) t)
200 ((endp q)
201 (setf r (revappend r p)) t)
202 (t
203 (multiple-value-bind
204 (lm-greater lm-equal)
205 (monomial-order (termlist-lm p) (termlist-lm q))
206 (cond
207 (lm-equal
208 (let ((s (funcall (ring-add ring) (termlist-lc p) (termlist-lc q))))
209 (unless (funcall (ring-zerop ring) s) ;check for cancellation
210 (setf r (cons (make-term (termlist-lm p) s) r)))
211 (setf p (cdr p) q (cdr q))))
212 (lm-greater
213 (setf r (cons (car p) r)
214 p (cdr p)))
215 (t (setf r (cons (car q) r)
216 q (cdr q)))))
217 nil))
218 r)))
219
220(defun termlist-sub (ring p q)
221 (declare (type list p q))
222 (do (r)
223 ((cond
224 ((endp p)
225 (setf r (revappend r (termlist-uminus ring q)))
226 t)
227 ((endp q)
228 (setf r (revappend r p))
229 t)
230 (t
231 (multiple-value-bind
232 (mgreater mequal)
233 (monomial-order (termlist-lm p) (termlist-lm q))
234 (cond
235 (mequal
236 (let ((s (funcall (ring-sub ring) (termlist-lc p) (termlist-lc q))))
237 (unless (funcall (ring-zerop ring) s) ;check for cancellation
238 (setf r (cons (make-term (termlist-lm p) s) r)))
239 (setf p (cdr p) q (cdr q))))
240 (mgreater
241 (setf r (cons (car p) r)
242 p (cdr p)))
243 (t (setf r (cons (make-term (termlist-lm q) (funcall (ring-uminus ring) (termlist-lc q))) r)
244 q (cdr q)))))
245 nil))
246 r)))
247
248;; Multiplication of polynomials
249;; Non-destructive version
250(defun termlist-mul (ring p q)
251 (cond ((or (endp p) (endp q)) nil) ;p or q is 0 (represented by NIL)
252 ;; If p=p0+p1 and q=q0+q1 then pq=p0q0+p0q1+p1q
253 ((endp (cdr p))
254 (term-times-termlist ring (car p) q))
255 ((endp (cdr q))
256 (termlist-times-term ring p (car q)))
257 (t
258 (let ((head (term-mul ring (termlist-lt p) (termlist-lt q)))
259 (tail (termlist-add ring (term-times-termlist ring (car p) (cdr q))
260 (termlist-mul ring (cdr p) q))))
261 (cond ((null head) tail)
262 ((null tail) head)
263 (t (nconc head tail)))))))
264
265(defun termlist-unit (ring dimension)
266 (declare (fixnum dimension))
267 (list (make-term (make-monom dimension :initial-element 0)
268 (funcall (ring-unit ring)))))
269
270(defun termlist-expt (ring poly n &aux (dim (monom-dimension (termlist-lm poly))))
271 (declare (type fixnum n dim))
272 (cond
273 ((minusp n) (error "termlist-expt: Negative exponent."))
274 ((endp poly) (if (zerop n) (termlist-unit ring dim) nil))
275 (t
276 (do ((k 1 (ash k 1))
277 (q poly (termlist-mul ring q q)) ;keep squaring
278 (p (termlist-unit ring dim) (if (not (zerop (logand k n))) (termlist-mul ring p q) p)))
279 ((> k n) p)
280 (declare (fixnum k))))))
281
282
283
284
285
286
287;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
288;;
289;; Debugging/tracing
290;;
291;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
292(defmacro debug-cgb (&rest args)
293 `(when $poly_grobner_debug (format *terminal-io* ,@args)))
294
295;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
296;;
297;; An implementation of Grobner basis
298;;
299;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
300
301
302
303
304;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
305;;
306;; An implementation of the division algorithm
307;;
308;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
309
310(defun grobner-op (ring c1 c2 m f g)
311 "Returns C2*F-C1*M*G, where F and G are polynomials M is a monomial.
[42]312Assume that the leading terms will cancel."
[1]313 #+grobner-check(funcall (ring-zerop ring)
[45]314 (funcall (ring-sub ring)
315 (funcall (ring-mul ring) c2 (poly-lc f))
[1]316 (funcall (ring-mul ring) c1 (poly-lc g))))
317 #+grobner-check(monom-equal-p (poly-lm f) (monom-mul m (poly-lm g)))
318 ;; Note that we can drop the leading terms of f ang g
319 (poly-sub ring
320 (scalar-times-poly-1 ring c2 f)
321 (scalar-times-poly-1 ring c1 (monom-times-poly m g))))
322
323(defun poly-pseudo-divide (ring f fl)
324 "Pseudo-divide a polynomial F by the list of polynomials FL. Return
325multiple values. The first value is a list of quotients A. The second
326value is the remainder R. The third argument is a scalar coefficient
327C, such that C*F can be divided by FL within the ring of coefficients,
328which is not necessarily a field. Finally, the fourth value is an
329integer count of the number of reductions performed. The resulting
330objects satisfy the equation: C*F= sum A[i]*FL[i] + R."
331 (declare (type poly f) (list fl))
332 (do ((r (make-poly-zero))
333 (c (funcall (ring-unit ring)))
334 (a (make-list (length fl) :initial-element (make-poly-zero)))
335 (division-count 0)
336 (p f))
337 ((poly-zerop p)
338 (debug-cgb "~&~3T~d reduction~:p" division-count)
339 (when (poly-zerop r) (debug-cgb " ---> 0"))
340 (values (mapcar #'poly-nreverse a) (poly-nreverse r) c division-count))
341 (declare (fixnum division-count))
342 (do ((fl fl (rest fl)) ;scan list of divisors
343 (b a (rest b)))
344 ((cond
345 ((endp fl) ;no division occurred
346 (push (poly-lt p) (poly-termlist r)) ;move lt(p) to remainder
347 (setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
348 (pop (poly-termlist p)) ;remove lt(p) from p
349 t)
350 ((monom-divides-p (poly-lm (car fl)) (poly-lm p)) ;division occurred
351 (incf division-count)
352 (multiple-value-bind (gcd c1 c2)
353 (funcall (ring-ezgcd ring) (poly-lc (car fl)) (poly-lc p))
354 (declare (ignore gcd))
355 (let ((m (monom-div (poly-lm p) (poly-lm (car fl)))))
356 ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
357 (mapl #'(lambda (x)
358 (setf (car x) (scalar-times-poly ring c1 (car x))))
359 a)
360 (setf r (scalar-times-poly ring c1 r)
361 c (funcall (ring-mul ring) c c1)
362 p (grobner-op ring c2 c1 m p (car fl)))
363 (push (make-term m c2) (poly-termlist (car b))))
364 t)))))))
365
366(defun poly-exact-divide (ring f g)
367 "Divide a polynomial F by another polynomial G. Assume that exact division
368with no remainder is possible. Returns the quotient."
369 (declare (type poly f g))
370 (multiple-value-bind (quot rem coeff division-count)
371 (poly-pseudo-divide ring f (list g))
372 (declare (ignore division-count coeff)
373 (list quot)
374 (type poly rem)
375 (type fixnum division-count))
376 (unless (poly-zerop rem) (error "Exact division failed."))
377 (car quot)))
378
379
380
381;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
382;;
383;; An implementation of the normal form
384;;
385;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
386
387(defun normal-form-step (ring fl p r c division-count
388 &aux (g (find (poly-lm p) fl
389 :test #'monom-divisible-by-p
390 :key #'poly-lm)))
391 (cond
392 (g ;division possible
393 (incf division-count)
[25]394 (multiple-value-bind (gcd cg cp)
[1]395 (funcall (ring-ezgcd ring) (poly-lc g) (poly-lc p))
396 (declare (ignore gcd))
397 (let ((m (monom-div (poly-lm p) (poly-lm g))))
398 ;; Multiply the equation c*f=sum ai*fi+r+p by cg.
399 (setf r (scalar-times-poly ring cg r)
400 c (funcall (ring-mul ring) c cg)
401 ;; p := cg*p-cp*m*g
402 p (grobner-op ring cp cg m p g))))
403 (debug-cgb "/"))
404 (t ;no division possible
405 (push (poly-lt p) (poly-termlist r)) ;move lt(p) to remainder
406 (setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
407 (pop (poly-termlist p)) ;remove lt(p) from p
408 (debug-cgb "+")))
409 (values p r c division-count))
410
411;; Merge it sometime with poly-pseudo-divide
412(defun normal-form (ring f fl &optional (top-reduction-only $poly_top_reduction_only))
413 ;; Loop invariant: c*f0=sum ai*fi+r+f, where f0 is the initial value of f
414 #+grobner-check(when (null fl) (warn "normal-form: empty divisor list."))
415 (do ((r (make-poly-zero))
416 (c (funcall (ring-unit ring)))
417 (division-count 0))
418 ((or (poly-zerop f)
419 ;;(endp fl)
420 (and top-reduction-only (not (poly-zerop r))))
421 (progn
422 (debug-cgb "~&~3T~d reduction~:p" division-count)
423 (when (poly-zerop r)
424 (debug-cgb " ---> 0")))
425 (setf (poly-termlist f) (nreconc (poly-termlist r) (poly-termlist f)))
426 (values f c division-count))
427 (declare (fixnum division-count)
428 (type poly r))
429 (multiple-value-setq (f r c division-count)
430 (normal-form-step ring fl f r c division-count))))
431
432
433
434;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
435;;
436;; These are provided mostly for debugging purposes To enable
437;; verification of grobner bases with BUCHBERGER-CRITERION, do
438;; (pushnew :grobner-check *features*) and compile/load this file.
439;; With this feature, the calculations will slow down CONSIDERABLY.
440;;
441;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
442
443(defun buchberger-criterion (ring g)
444 "Returns T if G is a Grobner basis, by using the Buchberger
445criterion: for every two polynomials h1 and h2 in G the S-polynomial
446S(h1,h2) reduces to 0 modulo G."
447 (every
448 #'poly-zerop
449 (makelist (normal-form ring (spoly ring (elt g i) (elt g j)) g nil)
450 (i 0 (- (length g) 2))
451 (j (1+ i) (1- (length g))))))
452
453(defun grobner-test (ring g f)
454 "Test whether G is a Grobner basis and F is contained in G. Return T
455upon success and NIL otherwise."
456 (debug-cgb "~&GROBNER CHECK: ")
457 (let (($poly_grobner_debug nil)
458 (stat1 (buchberger-criterion ring g))
459 (stat2
460 (every #'poly-zerop
461 (makelist (normal-form ring (copy-tree (elt f i)) g nil)
462 (i 0 (1- (length f)))))))
463 (unless stat1 (error "~&Buchberger criterion failed."))
464 (unless stat2
465 (error "~&Original polys not in ideal spanned by Grobner.")))
466 (debug-cgb "~&GROBNER CHECK END")
467 t)
468
469
470
471;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
472;;
473;; Pair queue implementation
474;;
475;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
476
477(defun sugar-pair-key (p q &aux (lcm (monom-lcm (poly-lm p) (poly-lm q)))
478 (d (monom-sugar lcm)))
479 "Returns list (S LCM-TOTAL-DEGREE) where S is the sugar of the S-polynomial of
480polynomials P and Q, and LCM-TOTAL-DEGREE is the degree of is LCM(LM(P),LM(Q))."
481 (declare (type poly p q) (type monom lcm) (type fixnum d))
482 (cons (max
483 (+ (- d (monom-sugar (poly-lm p))) (poly-sugar p))
484 (+ (- d (monom-sugar (poly-lm q))) (poly-sugar q)))
485 lcm))
486
487(defstruct (pair
488 (:constructor make-pair (first second
489 &aux
490 (sugar (car (sugar-pair-key first second)))
491 (division-data nil))))
492 (first nil :type poly)
493 (second nil :type poly)
494 (sugar 0 :type fixnum)
495 (division-data nil :type list))
496
497;;(defun pair-sugar (pair &aux (p (pair-first pair)) (q (pair-second pair)))
498;; (car (sugar-pair-key p q)))
499
500(defun sugar-order (x y)
501 "Pair order based on sugar, ties broken by normal strategy."
502 (declare (type cons x y))
503 (or (< (car x) (car y))
504 (and (= (car x) (car y))
505 (< (monom-total-degree (cdr x))
506 (monom-total-degree (cdr y))))))
507
508(defvar *pair-key-function* #'sugar-pair-key
509 "Function that, given two polynomials as argument, computed the key
510in the pair queue.")
511
512(defvar *pair-order* #'sugar-order
513 "Function that orders the keys of pairs.")
514
515(defun make-pair-queue ()
516 "Constructs a priority queue for critical pairs."
517 (make-priority-queue
518 :element-type 'pair
519 :element-key #'(lambda (pair) (funcall *pair-key-function* (pair-first pair) (pair-second pair)))
520 :test *pair-order*))
521
522(defun pair-queue-initialize (pq f start
523 &aux
524 (s (1- (length f)))
525 (b (nconc (makelist (make-pair (elt f i) (elt f j))
526 (i 0 (1- start)) (j start s))
527 (makelist (make-pair (elt f i) (elt f j))
528 (i start (1- s)) (j (1+ i) s)))))
529 "Initializes the priority for critical pairs. F is the initial list of polynomials.
530START is the first position beyond the elements which form a partial
531grobner basis, i.e. satisfy the Buchberger criterion."
532 (declare (type priority-queue pq) (type fixnum start))
533 (dolist (pair b pq)
534 (priority-queue-insert pq pair)))
535
536(defun pair-queue-insert (b pair)
537 (priority-queue-insert b pair))
538
539(defun pair-queue-remove (b)
540 (priority-queue-remove b))
541
542(defun pair-queue-size (b)
543 (priority-queue-size b))
544
545(defun pair-queue-empty-p (b)
546 (priority-queue-empty-p b))
547
548;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
549;;
550;; Buchberger Algorithm Implementation
551;;
552;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
553
554(defun buchberger (ring f start &optional (top-reduction-only $poly_top_reduction_only))
555 "An implementation of the Buchberger algorithm. Return Grobner basis
556of the ideal generated by the polynomial list F. Polynomials 0 to
557START-1 are assumed to be a Grobner basis already, so that certain
558critical pairs will not be examined. If TOP-REDUCTION-ONLY set, top
559reduction will be preformed. This function assumes that all polynomials
560in F are non-zero."
561 (declare (type fixnum start))
562 (when (endp f) (return-from buchberger f)) ;cut startup costs
563 (debug-cgb "~&GROBNER BASIS - BUCHBERGER ALGORITHM")
564 (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
565 #+grobner-check (when (plusp start)
566 (grobner-test ring (subseq f 0 start) (subseq f 0 start)))
567 ;;Initialize critical pairs
568 (let ((b (pair-queue-initialize (make-pair-queue)
569 f start))
570 (b-done (make-hash-table :test #'equal)))
571 (declare (type priority-queue b) (type hash-table b-done))
572 (dotimes (i (1- start))
573 (do ((j (1+ i) (1+ j))) ((>= j start))
574 (setf (gethash (list (elt f i) (elt f j)) b-done) t)))
575 (do ()
576 ((pair-queue-empty-p b)
577 #+grobner-check(grobner-test ring f f)
578 (debug-cgb "~&GROBNER END")
579 f)
580 (let ((pair (pair-queue-remove b)))
581 (declare (type pair pair))
582 (cond
583 ((criterion-1 pair) nil)
584 ((criterion-2 pair b-done f) nil)
585 (t
586 (let ((sp (normal-form ring (spoly ring (pair-first pair)
587 (pair-second pair))
588 f top-reduction-only)))
589 (declare (type poly sp))
590 (cond
591 ((poly-zerop sp)
592 nil)
593 (t
594 (setf sp (poly-primitive-part ring sp)
595 f (nconc f (list sp)))
596 ;; Add new critical pairs
597 (dolist (h f)
598 (pair-queue-insert b (make-pair h sp)))
599 (debug-cgb "~&Sugar: ~d Polynomials: ~d; Pairs left: ~d; Pairs done: ~d;"
600 (pair-sugar pair) (length f) (pair-queue-size b)
601 (hash-table-count b-done)))))))
602 (setf (gethash (list (pair-first pair) (pair-second pair)) b-done)
603 t)))))
604
605(defun parallel-buchberger (ring f start &optional (top-reduction-only $poly_top_reduction_only))
606 "An implementation of the Buchberger algorithm. Return Grobner basis
607of the ideal generated by the polynomial list F. Polynomials 0 to
608START-1 are assumed to be a Grobner basis already, so that certain
609critical pairs will not be examined. If TOP-REDUCTION-ONLY set, top
610reduction will be preformed."
611 (declare (ignore top-reduction-only)
612 (type fixnum start))
613 (when (endp f) (return-from parallel-buchberger f)) ;cut startup costs
614 (debug-cgb "~&GROBNER BASIS - PARALLEL-BUCHBERGER ALGORITHM")
615 (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
616 #+grobner-check (when (plusp start)
617 (grobner-test ring (subseq f 0 start) (subseq f 0 start)))
618 ;;Initialize critical pairs
619 (let ((b (pair-queue-initialize (make-pair-queue) f start))
620 (b-done (make-hash-table :test #'equal)))
621 (declare (type priority-queue b)
622 (type hash-table b-done))
623 (dotimes (i (1- start))
624 (do ((j (1+ i) (1+ j))) ((>= j start))
625 (declare (type fixnum j))
626 (setf (gethash (list (elt f i) (elt f j)) b-done) t)))
627 (do ()
628 ((pair-queue-empty-p b)
629 #+grobner-check(grobner-test ring f f)
630 (debug-cgb "~&GROBNER END")
631 f)
632 (let ((pair (pair-queue-remove b)))
633 (when (null (pair-division-data pair))
634 (setf (pair-division-data pair) (list (spoly ring
635 (pair-first pair)
636 (pair-second pair))
637 (make-poly-zero)
638 (funcall (ring-unit ring))
639 0)))
640 (cond
641 ((criterion-1 pair) nil)
642 ((criterion-2 pair b-done f) nil)
643 (t
644 (let* ((dd (pair-division-data pair))
645 (p (first dd))
646 (sp (second dd))
647 (c (third dd))
648 (division-count (fourth dd)))
649 (cond
650 ((poly-zerop p) ;normal form completed
651 (debug-cgb "~&~3T~d reduction~:p" division-count)
652 (cond
653 ((poly-zerop sp)
654 (debug-cgb " ---> 0")
655 nil)
656 (t
657 (setf sp (poly-nreverse sp)
658 sp (poly-primitive-part ring sp)
659 f (nconc f (list sp)))
660 ;; Add new critical pairs
661 (dolist (h f)
662 (pair-queue-insert b (make-pair h sp)))
663 (debug-cgb "~&Sugar: ~d Polynomials: ~d; Pairs left: ~d; Pairs done: ~d;"
664 (pair-sugar pair) (length f) (pair-queue-size b)
665 (hash-table-count b-done))))
666 (setf (gethash (list (pair-first pair) (pair-second pair))
667 b-done) t))
668 (t ;normal form not complete
669 (do ()
670 ((cond
671 ((> (poly-sugar sp) (pair-sugar pair))
672 (debug-cgb "(~a)?" (poly-sugar sp))
673 t)
674 ((poly-zerop p)
675 (debug-cgb ".")
676 t)
677 (t nil))
678 (setf (first dd) p
679 (second dd) sp
680 (third dd) c
681 (fourth dd) division-count
682 (pair-sugar pair) (poly-sugar sp))
683 (pair-queue-insert b pair))
684 (multiple-value-setq (p sp c division-count)
685 (normal-form-step ring f p sp c division-count))))))))))))
686
687;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
688;;
689;; Grobner Criteria
690;;
691;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
692
693(defun criterion-1 (pair)
694 "Returns T if the leading monomials of the two polynomials
695in G pointed to by the integers in PAIR have disjoint (relatively prime)
696monomials. This test is known as the first Buchberger criterion."
697 (declare (type pair pair))
698 (let ((f (pair-first pair))
699 (g (pair-second pair)))
700 (when (monom-rel-prime-p (poly-lm f) (poly-lm g))
701 (debug-cgb ":1")
702 (return-from criterion-1 t))))
703
704(defun criterion-2 (pair b-done partial-basis
705 &aux (f (pair-first pair)) (g (pair-second pair))
706 (place :before))
707 "Returns T if the leading monomial of some element P of
708PARTIAL-BASIS divides the LCM of the leading monomials of the two
709polynomials in the polynomial list PARTIAL-BASIS, and P paired with
710each of the polynomials pointed to by the the PAIR has already been
711treated, as indicated by the absence in the hash table B-done."
712 (declare (type pair pair) (type hash-table b-done)
713 (type poly f g))
714 ;; In the code below we assume that pairs are ordered as follows:
715 ;; if PAIR is (I J) then I appears before J in the PARTIAL-BASIS.
716 ;; We traverse the list PARTIAL-BASIS and keep track of where we
717 ;; are, so that we can produce the pairs in the correct order
718 ;; when we check whether they have been processed, i.e they
719 ;; appear in the hash table B-done
720 (dolist (h partial-basis nil)
721 (cond
722 ((eq h f)
723 #+grobner-check(assert (eq place :before))
724 (setf place :in-the-middle))
725 ((eq h g)
726 #+grobner-check(assert (eq place :in-the-middle))
727 (setf place :after))
728 ((and (monom-divides-monom-lcm-p (poly-lm h) (poly-lm f) (poly-lm g))
729 (gethash (case place
730 (:before (list h f))
731 ((:in-the-middle :after) (list f h)))
732 b-done)
733 (gethash (case place
734 ((:before :in-the-middle) (list h g))
735 (:after (list g h)))
736 b-done))
737 (debug-cgb ":2")
738 (return-from criterion-2 t)))))
739
740
741
742;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
743;;
744;; An implementation of the algorithm of Gebauer and Moeller, as
745;; described in the book of Becker-Weispfenning, p. 232
746;;
747;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
748
749(defun gebauer-moeller (ring f start &optional (top-reduction-only $poly_top_reduction_only))
750 "Compute Grobner basis by using the algorithm of Gebauer and
751Moeller. This algorithm is described as BUCHBERGERNEW2 in the book by
752Becker-Weispfenning entitled ``Grobner Bases''. This function assumes
753that all polynomials in F are non-zero."
754 (declare (ignore top-reduction-only)
755 (type fixnum start))
756 (cond
757 ((endp f) (return-from gebauer-moeller nil))
758 ((endp (cdr f))
759 (return-from gebauer-moeller (list (poly-primitive-part ring (car f))))))
760 (debug-cgb "~&GROBNER BASIS - GEBAUER MOELLER ALGORITHM")
761 (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
762 #+grobner-check (when (plusp start)
763 (grobner-test ring (subseq f 0 start) (subseq f 0 start)))
764 (let ((b (make-pair-queue))
765 (g (subseq f 0 start))
766 (f1 (subseq f start)))
767 (do () ((endp f1))
768 (multiple-value-setq (g b)
769 (gebauer-moeller-update g b (poly-primitive-part ring (pop f1)))))
770 (do () ((pair-queue-empty-p b))
771 (let* ((pair (pair-queue-remove b))
772 (g1 (pair-first pair))
773 (g2 (pair-second pair))
774 (h (normal-form ring (spoly ring g1 g2)
775 g
776 nil #| Always fully reduce! |#
777 )))
778 (unless (poly-zerop h)
779 (setf h (poly-primitive-part ring h))
780 (multiple-value-setq (g b)
781 (gebauer-moeller-update g b h))
782 (debug-cgb "~&Sugar: ~d Polynomials: ~d; Pairs left: ~d~%"
783 (pair-sugar pair) (length g) (pair-queue-size b))
784 )))
785 #+grobner-check(grobner-test ring g f)
786 (debug-cgb "~&GROBNER END")
787 g))
788
789(defun gebauer-moeller-update (g b h
790 &aux
791 c d e
792 (b-new (make-pair-queue))
793 g-new)
794 "An implementation of the auxillary UPDATE algorithm used by the
795Gebauer-Moeller algorithm. G is a list of polynomials, B is a list of
796critical pairs and H is a new polynomial which possibly will be added
797to G. The naming conventions used are very close to the one used in
798the book of Becker-Weispfenning."
799 (declare
800 #+allegro (dynamic-extent b)
801 (type poly h)
802 (type priority-queue b))
803 (setf c g d nil)
804 (do () ((endp c))
805 (let ((g1 (pop c)))
806 (declare (type poly g1))
807 (when (or (monom-rel-prime-p (poly-lm h) (poly-lm g1))
808 (and
809 (notany #'(lambda (g2) (monom-lcm-divides-monom-lcm-p
810 (poly-lm h) (poly-lm g2)
811 (poly-lm h) (poly-lm g1)))
812 c)
813 (notany #'(lambda (g2) (monom-lcm-divides-monom-lcm-p
814 (poly-lm h) (poly-lm g2)
815 (poly-lm h) (poly-lm g1)))
816 d)))
817 (push g1 d))))
818 (setf e nil)
819 (do () ((endp d))
820 (let ((g1 (pop d)))
821 (declare (type poly g1))
822 (unless (monom-rel-prime-p (poly-lm h) (poly-lm g1))
823 (push g1 e))))
824 (do () ((pair-queue-empty-p b))
825 (let* ((pair (pair-queue-remove b))
826 (g1 (pair-first pair))
827 (g2 (pair-second pair)))
828 (declare (type pair pair)
829 (type poly g1 g2))
830 (when (or (not (monom-divides-monom-lcm-p
831 (poly-lm h)
832 (poly-lm g1) (poly-lm g2)))
833 (monom-lcm-equal-monom-lcm-p
834 (poly-lm g1) (poly-lm h)
835 (poly-lm g1) (poly-lm g2))
836 (monom-lcm-equal-monom-lcm-p
837 (poly-lm h) (poly-lm g2)
838 (poly-lm g1) (poly-lm g2)))
839 (pair-queue-insert b-new (make-pair g1 g2)))))
840 (dolist (g3 e)
841 (pair-queue-insert b-new (make-pair h g3)))
842 (setf g-new nil)
843 (do () ((endp g))
844 (let ((g1 (pop g)))
845 (declare (type poly g1))
846 (unless (monom-divides-p (poly-lm h) (poly-lm g1))
847 (push g1 g-new))))
848 (push h g-new)
849 (values g-new b-new))
850
851
852
853;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
854;;
855;; Standard postprocessing of Grobner bases
856;;
857;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
858
859(defun reduction (ring plist)
860 "Reduce a list of polynomials PLIST, so that non of the terms in any of
861the polynomials is divisible by a leading monomial of another
862polynomial. Return the reduced list."
863 (do ((q plist)
864 (found t))
865 ((not found)
866 (mapcar #'(lambda (x) (poly-primitive-part ring x)) q))
867 ;;Find p in Q such that p is reducible mod Q\{p}
868 (setf found nil)
869 (dolist (x q)
870 (let ((q1 (remove x q)))
871 (multiple-value-bind (h c div-count)
872 (normal-form ring x q1 nil #| not a top reduction! |# )
873 (declare (ignore c))
874 (unless (zerop div-count)
875 (setf found t q q1)
876 (unless (poly-zerop h)
877 (setf q (nconc q1 (list h))))
878 (return)))))))
879
880(defun minimization (p)
881 "Returns a sublist of the polynomial list P spanning the same
882monomial ideal as P but minimal, i.e. no leading monomial
883of a polynomial in the sublist divides the leading monomial
884of another polynomial."
885 (do ((q p)
886 (found t))
887 ((not found) q)
888 ;;Find p in Q such that lm(p) is in LM(Q\{p})
889 (setf found nil
890 q (dolist (x q q)
891 (let ((q1 (remove x q)))
892 (when (member-if #'(lambda (p) (monom-divides-p (poly-lm x) (poly-lm p))) q1)
893 (setf found t)
894 (return q1)))))))
895
896(defun poly-normalize (ring p &aux (c (poly-lc p)))
897 "Divide a polynomial by its leading coefficient. It assumes
898that the division is possible, which may not always be the
899case in rings which are not fields. The exact division operator
900is assumed to be provided by the RING structure of the
901COEFFICIENT-RING package."
902 (mapc #'(lambda (term)
903 (setf (term-coeff term) (funcall (ring-div ring) (term-coeff term) c)))
904 (poly-termlist p))
905 p)
906
907(defun poly-normalize-list (ring plist)
908 "Divide every polynomial in a list PLIST by its leading coefficient. "
909 (mapcar #'(lambda (x) (poly-normalize ring x)) plist))
910
911
912
913;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
914;;
915;; Algorithm and Pair heuristic selection
916;;
917;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
918
919(defun find-grobner-function (algorithm)
920 "Return a function which calculates Grobner basis, based on its
921names. Names currently used are either Lisp symbols, Maxima symbols or
922keywords."
923 (ecase algorithm
924 ((buchberger :buchberger $buchberger) #'buchberger)
925 ((parallel-buchberger :parallel-buchberger $parallel_buchberger) #'parallel-buchberger)
926 ((gebauer-moeller :gebauer_moeller $gebauer_moeller) #'gebauer-moeller)))
927
928(defun grobner (ring f &optional (start 0) (top-reduction-only nil))
929 ;;(setf F (sort F #'< :key #'sugar))
930 (funcall
931 (find-grobner-function $poly_grobner_algorithm)
932 ring f start top-reduction-only))
933
934(defun reduced-grobner (ring f &optional (start 0) (top-reduction-only $poly_top_reduction_only))
935 (reduction ring (grobner ring f start top-reduction-only)))
936
937(defun set-pair-heuristic (method)
938 "Sets up variables *PAIR-KEY-FUNCTION* and *PAIR-ORDER* used
939to determine the priority of critical pairs in the priority queue."
940 (ecase method
941 ((sugar :sugar $sugar)
942 (setf *pair-key-function* #'sugar-pair-key
943 *pair-order* #'sugar-order))
944; ((minimal-mock-spoly :minimal-mock-spoly $minimal_mock_spoly)
945; (setf *pair-key-function* #'mock-spoly
946; *pair-order* #'mock-spoly-order))
947 ((minimal-lcm :minimal-lcm $minimal_lcm)
948 (setf *pair-key-function* #'(lambda (p q)
949 (monom-lcm (poly-lm p) (poly-lm q)))
950 *pair-order* #'reverse-monomial-order))
951 ((minimal-total-degree :minimal-total-degree $minimal_total_degree)
952 (setf *pair-key-function* #'(lambda (p q)
953 (monom-total-degree
954 (monom-lcm (poly-lm p) (poly-lm q))))
955 *pair-order* #'<))
956 ((minimal-length :minimal-length $minimal_length)
957 (setf *pair-key-function* #'(lambda (p q)
958 (+ (poly-length p) (poly-length q)))
959 *pair-order* #'<))))
960
961
962
963;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
964;;
965;; Operations in ideal theory
966;;
967;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
968
969;; Does the term depend on variable K?
970(defun term-depends-p (term k)
971 "Return T if the term TERM depends on variable number K."
972 (monom-depends-p (term-monom term) k))
973
974;; Does the polynomial P depend on variable K?
975(defun poly-depends-p (p k)
976 "Return T if the term polynomial P depends on variable number K."
977 (some #'(lambda (term) (term-depends-p term k)) (poly-termlist p)))
978
979(defun ring-intersection (plist k)
980 "This function assumes that polynomial list PLIST is a Grobner basis
981and it calculates the intersection with the ring R[x[k+1],...,x[n]], i.e.
982it discards polynomials which depend on variables x[0], x[1], ..., x[k]."
983 (dotimes (i k plist)
984 (setf plist
985 (remove-if #'(lambda (p)
986 (poly-depends-p p i))
987 plist))))
988
[21]989(defun elimination-ideal (ring flist k
990 &optional (top-reduction-only $poly_top_reduction_only) (start 0)
991 &aux (*monomial-order*
992 (or *elimination-order*
993 (elimination-order k))))
994 (ring-intersection (reduced-grobner ring flist start top-reduction-only) k))
995
996(defun colon-ideal (ring f g &optional (top-reduction-only $poly_top_reduction_only))
997 "Returns the reduced Grobner basis of the colon ideal Id(F):Id(G),
998where F and G are two lists of polynomials. The colon ideal I:J is
999defined as the set of polynomials H such that for all polynomials W in
1000J the polynomial W*H belongs to I."
1001 (cond
1002 ((endp g)
1003 ;;Id(G) consists of 0 only so W*0=0 belongs to Id(F)
1004 (if (every #'poly-zerop f)
[1]1005 (error "First ideal must be non-zero.")
1006 (list (make-poly
1007 (list (make-term
1008 (make-monom (monom-dimension (poly-lm (find-if-not #'poly-zerop f)))
1009 :initial-element 0)
1010 (funcall (ring-unit ring))))))))
1011 ((endp (cdr g))
1012 (colon-ideal-1 ring f (car g) top-reduction-only))
1013 (t
1014 (ideal-intersection ring
1015 (colon-ideal-1 ring f (car g) top-reduction-only)
1016 (colon-ideal ring f (rest g) top-reduction-only)
1017 top-reduction-only))))
1018
1019(defun colon-ideal-1 (ring f g &optional (top-reduction-only $poly_top_reduction_only))
1020 "Returns the reduced Grobner basis of the colon ideal Id(F):Id({G}), where
1021F is a list of polynomials and G is a polynomial."
1022 (mapcar #'(lambda (x) (poly-exact-divide ring x g)) (ideal-intersection ring f (list g) top-reduction-only)))
1023
1024
1025(defun ideal-intersection (ring f g &optional (top-reduction-only $poly_top_reduction_only)
1026 &aux (*monomial-order* (or *elimination-order*
1027 #'elimination-order-1)))
1028 (mapcar #'poly-contract
1029 (ring-intersection
1030 (reduced-grobner
1031 ring
1032 (append (mapcar #'(lambda (p) (poly-extend p (make-monom 1 :initial-element 1))) f)
1033 (mapcar #'(lambda (p)
1034 (poly-append (poly-extend (poly-uminus ring p)
1035 (make-monom 1 :initial-element 1))
1036 (poly-extend p)))
1037 g))
1038 0
1039 top-reduction-only)
1040 1)))
1041
1042(defun poly-lcm (ring f g)
1043 "Return LCM (least common multiple) of two polynomials F and G.
1044The polynomials must be ordered according to monomial order PRED
1045and their coefficients must be compatible with the RING structure
1046defined in the COEFFICIENT-RING package."
1047 (cond
1048 ((poly-zerop f) f)
1049 ((poly-zerop g) g)
1050 ((and (endp (cdr (poly-termlist f))) (endp (cdr (poly-termlist g))))
1051 (let ((m (monom-lcm (poly-lm f) (poly-lm g))))
1052 (make-poly-from-termlist (list (make-term m (funcall (ring-lcm ring) (poly-lc f) (poly-lc g)))))))
1053 (t
1054 (multiple-value-bind (f f-cont)
1055 (poly-primitive-part ring f)
1056 (multiple-value-bind (g g-cont)
1057 (poly-primitive-part ring g)
1058 (scalar-times-poly
1059 ring
1060 (funcall (ring-lcm ring) f-cont g-cont)
1061 (poly-primitive-part ring (car (ideal-intersection ring (list f) (list g) nil)))))))))
1062
1063;; Do two Grobner bases yield the same ideal?
1064(defun grobner-equal (ring g1 g2)
1065 "Returns T if two lists of polynomials G1 and G2, assumed to be Grobner bases,
1066generate the same ideal, and NIL otherwise."
1067 (and (grobner-subsetp ring g1 g2) (grobner-subsetp ring g2 g1)))
1068
1069(defun grobner-subsetp (ring g1 g2)
1070 "Returns T if a list of polynomials G1 generates
1071an ideal contained in the ideal generated by a polynomial list G2,
1072both G1 and G2 assumed to be Grobner bases. Returns NIL otherwise."
1073 (every #'(lambda (p) (grobner-member ring p g2)) g1))
1074
1075(defun grobner-member (ring p g)
1076 "Returns T if a polynomial P belongs to the ideal generated by the
1077polynomial list G, which is assumed to be a Grobner basis. Returns NIL otherwise."
1078 (poly-zerop (normal-form ring p g nil)))
1079
1080;; Calculate F : p^inf
1081(defun ideal-saturation-1 (ring f p start &optional (top-reduction-only $poly_top_reduction_only)
1082 &aux (*monomial-order* (or *elimination-order*
1083 #'elimination-order-1)))
1084 "Returns the reduced Grobner basis of the saturation of the ideal
1085generated by a polynomial list F in the ideal generated by a single
1086polynomial P. The saturation ideal is defined as the set of
1087polynomials H such for some natural number n (* (EXPT P N) H) is in the ideal
1088F. Geometrically, over an algebraically closed field, this is the set
1089of polynomials in the ideal generated by F which do not identically
1090vanish on the variety of P."
1091 (mapcar
1092 #'poly-contract
1093 (ring-intersection
1094 (reduced-grobner
1095 ring
1096 (saturation-extension-1 ring f p)
1097 start top-reduction-only)
1098 1)))
1099
1100
1101
1102;; Calculate F : p1^inf : p2^inf : ... : ps^inf
1103(defun ideal-polysaturation-1 (ring f plist start &optional (top-reduction-only $poly_top_reduction_only))
1104 "Returns the reduced Grobner basis of the ideal obtained by a
1105sequence of successive saturations in the polynomials
1106of the polynomial list PLIST of the ideal generated by the
1107polynomial list F."
1108 (cond
1109 ((endp plist) (reduced-grobner ring f start top-reduction-only))
1110 (t (let ((g (ideal-saturation-1 ring f (car plist) start top-reduction-only)))
1111 (ideal-polysaturation-1 ring g (rest plist) (length g) top-reduction-only)))))
1112
1113(defun ideal-saturation (ring f g start &optional (top-reduction-only $poly_top_reduction_only)
1114 &aux
1115 (k (length g))
1116 (*monomial-order* (or *elimination-order*
1117 (elimination-order k))))
1118 "Returns the reduced Grobner basis of the saturation of the ideal
1119generated by a polynomial list F in the ideal generated a polynomial
1120list G. The saturation ideal is defined as the set of polynomials H
1121such for some natural number n and some P in the ideal generated by G
1122the polynomial P**N * H is in the ideal spanned by F. Geometrically,
1123over an algebraically closed field, this is the set of polynomials in
1124the ideal generated by F which do not identically vanish on the
1125variety of G."
1126 (mapcar
1127 #'(lambda (q) (poly-contract q k))
1128 (ring-intersection
1129 (reduced-grobner ring
1130 (polysaturation-extension ring f g)
1131 start
1132 top-reduction-only)
1133 k)))
1134
1135(defun ideal-polysaturation (ring f ideal-list start &optional (top-reduction-only $poly_top_reduction_only))
1136 "Returns the reduced Grobner basis of the ideal obtained by a
1137successive applications of IDEAL-SATURATION to F and lists of
1138polynomials in the list IDEAL-LIST."
1139 (cond
1140 ((endp ideal-list) f)
1141 (t (let ((h (ideal-saturation ring f (car ideal-list) start top-reduction-only)))
1142 (ideal-polysaturation ring h (rest ideal-list) (length h) top-reduction-only)))))
1143
1144
1145
1146;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1147;;
1148;; Set up the coefficients to be polynomials
1149;;
1150;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1151
1152;; (defun poly-ring (ring vars)
1153;; (make-ring
1154;; :parse #'(lambda (expr) (poly-eval ring expr vars))
1155;; :unit #'(lambda () (poly-unit ring (length vars)))
1156;; :zerop #'poly-zerop
1157;; :add #'(lambda (x y) (poly-add ring x y))
1158;; :sub #'(lambda (x y) (poly-sub ring x y))
1159;; :uminus #'(lambda (x) (poly-uminus ring x))
1160;; :mul #'(lambda (x y) (poly-mul ring x y))
1161;; :div #'(lambda (x y) (poly-exact-divide ring x y))
1162;; :lcm #'(lambda (x y) (poly-lcm ring x y))
1163;; :ezgcd #'(lambda (x y &aux (gcd (poly-gcd ring x y)))
1164;; (values gcd
1165;; (poly-exact-divide ring x gcd)
1166;; (poly-exact-divide ring y gcd)))
1167;; :gcd #'(lambda (x y) (poly-gcd x y))))
1168
1169
1170
1171;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1172;;
1173;; Conversion from internal to infix form
1174;;
1175;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1176
1177(defun coerce-to-infix (poly-type object vars)
1178 (case poly-type
1179 (:termlist
1180 `(+ ,@(mapcar #'(lambda (term) (coerce-to-infix :term term vars)) object)))
1181 (:polynomial
1182 (coerce-to-infix :termlist (poly-termlist object) vars))
1183 (:poly-list
1184 `([ ,@(mapcar #'(lambda (p) (coerce-to-infix :polynomial p vars)) object)))
1185 (:term
1186 `(* ,(term-coeff object)
1187 ,@(mapcar #'(lambda (var power) `(expt ,var ,power))
1188 vars (monom-exponents (term-monom object)))))
1189 (otherwise
1190 object)))
1191
1192
1193
1194;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1195;;
1196;; Maxima expression ring
1197;;
1198;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1199
1200(defparameter *expression-ring*
1201 (make-ring
1202 ;;(defun coeff-zerop (expr) (meval1 `(($is) (($equal) ,expr 0))))
1203 :parse #'(lambda (expr)
1204 (when modulus (setf expr ($rat expr)))
1205 expr)
1206 :unit #'(lambda () (if modulus ($rat 1) 1))
[12]1207 :zerop #'(lambda (expr)
1208 ;;When is exactly a maxima expression equal to 0?
1209 (cond ((numberp expr)
[1]1210 (= expr 0))
1211 ((atom expr) nil)
1212 (t
1213 (case (caar expr)
1214 (mrat (eql ($ratdisrep expr) 0))
1215 (otherwise (eql ($totaldisrep expr) 0))))))
1216 :add #'(lambda (x y) (m+ x y))
1217 :sub #'(lambda (x y) (m- x y))
1218 :uminus #'(lambda (x) (m- x))
1219 :mul #'(lambda (x y) (m* x y))
1220 ;;(defun coeff-div (x y) (cadr ($divide x y)))
1221 :div #'(lambda (x y) (m// x y))
1222 :lcm #'(lambda (x y) (meval1 `((|$LCM|) ,x ,y)))
1223 :ezgcd #'(lambda (x y) (apply #'values (cdr ($ezgcd ($totaldisrep x) ($totaldisrep y)))))
1224 ;; :gcd #'(lambda (x y) (second ($ezgcd x y)))))
1225 :gcd #'(lambda (x y) ($gcd x y))))
1226
1227(defvar *maxima-ring* *expression-ring*
1228 "The ring of coefficients, over which all polynomials
1229are assumed to be defined.")
1230
1231
1232
1233;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1234;;
1235;; Maxima expression parsing
1236;;
1237;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1238
1239(defun equal-test-p (expr1 expr2)
1240 (alike1 expr1 expr2))
1241
1242(defun coerce-maxima-list (expr)
1243 "convert a maxima list to lisp list."
1244 (cond
1245 ((and (consp (car expr)) (eql (caar expr) 'mlist)) (cdr expr))
1246 (t expr)))
1247
1248(defun free-of-vars (expr vars) (apply #'$freeof `(,@vars ,expr)))
1249
1250(defun parse-poly (expr vars &aux (vars (coerce-maxima-list vars)))
1251 "Convert a maxima polynomial expression EXPR in variables VARS to internal form."
1252 (labels ((parse (arg) (parse-poly arg vars))
1253 (parse-list (args) (mapcar #'parse args)))
1254 (cond
1255 ((eql expr 0) (make-poly-zero))
1256 ((member expr vars :test #'equal-test-p)
1257 (let ((pos (position expr vars :test #'equal-test-p)))
1258 (make-variable *maxima-ring* (length vars) pos)))
1259 ((free-of-vars expr vars)
1260 ;;This means that variable-free CRE and Poisson forms will be converted
1261 ;;to coefficients intact
1262 (coerce-coeff *maxima-ring* expr vars))
1263 (t
1264 (case (caar expr)
1265 (mplus (reduce #'(lambda (x y) (poly-add *maxima-ring* x y)) (parse-list (cdr expr))))
1266 (mminus (poly-uminus *maxima-ring* (parse (cadr expr))))
1267 (mtimes
1268 (if (endp (cddr expr)) ;unary
1269 (parse (cdr expr))
1270 (reduce #'(lambda (p q) (poly-mul *maxima-ring* p q)) (parse-list (cdr expr)))))
1271 (mexpt
1272 (cond
1273 ((member (cadr expr) vars :test #'equal-test-p)
1274 ;;Special handling of (expt var pow)
1275 (let ((pos (position (cadr expr) vars :test #'equal-test-p)))
1276 (make-variable *maxima-ring* (length vars) pos (caddr expr))))
1277 ((not (and (integerp (caddr expr)) (plusp (caddr expr))))
1278 ;; Negative power means division in coefficient ring
1279 ;; Non-integer power means non-polynomial coefficient
1280 (mtell "~%Warning: Expression ~%~M~%contains power which is not a positive integer. Parsing as coefficient.~%"
1281 expr)
1282 (coerce-coeff *maxima-ring* expr vars))
1283 (t (poly-expt *maxima-ring* (parse (cadr expr)) (caddr expr)))))
1284 (mrat (parse ($ratdisrep expr)))
1285 (mpois (parse ($outofpois expr)))
1286 (otherwise
1287 (coerce-coeff *maxima-ring* expr vars)))))))
1288
1289(defun parse-poly-list (expr vars)
1290 (case (caar expr)
1291 (mlist (mapcar #'(lambda (p) (parse-poly p vars)) (cdr expr)))
1292 (t (merror "Expression ~M is not a list of polynomials in variables ~M."
1293 expr vars))))
1294(defun parse-poly-list-list (poly-list-list vars)
1295 (mapcar #'(lambda (g) (parse-poly-list g vars)) (coerce-maxima-list poly-list-list)))
1296
1297
1298
1299;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1300;;
1301;; Order utilities
1302;;
1303;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1304(defun find-order (order)
1305 "This function returns the order function bases on its name."
1306 (cond
1307 ((null order) nil)
1308 ((symbolp order)
1309 (case order
1310 ((lex :lex $lex) #'lex>)
1311 ((grlex :grlex $grlex) #'grlex>)
1312 ((grevlex :grevlex $grevlex) #'grevlex>)
1313 ((invlex :invlex $invlex) #'invlex>)
1314 ((elimination-order-1 :elimination-order-1 elimination_order_1) #'elimination-order-1)
1315 (otherwise
1316 (mtell "~%Warning: Order ~M not found. Using default.~%" order))))
1317 (t
1318 (mtell "~%Order specification ~M is not recognized. Using default.~%" order)
1319 nil)))
1320
1321(defun find-ring (ring)
1322 "This function returns the ring structure bases on input symbol."
1323 (cond
1324 ((null ring) nil)
1325 ((symbolp ring)
1326 (case ring
1327 ((expression-ring :expression-ring $expression_ring) *expression-ring*)
1328 ((ring-of-integers :ring-of-integers $ring_of_integers) *ring-of-integers*)
1329 (otherwise
1330 (mtell "~%Warning: Ring ~M not found. Using default.~%" ring))))
1331 (t
1332 (mtell "~%Ring specification ~M is not recognized. Using default.~%" ring)
1333 nil)))
1334
1335(defmacro with-monomial-order ((order) &body body)
1336 "Evaluate BODY with monomial order set to ORDER."
1337 `(let ((*monomial-order* (or (find-order ,order) *monomial-order*)))
1338 . ,body))
1339
1340(defmacro with-coefficient-ring ((ring) &body body)
1341 "Evaluate BODY with coefficient ring set to RING."
1342 `(let ((*maxima-ring* (or (find-ring ,ring) *maxima-ring*)))
1343 . ,body))
1344
1345(defmacro with-elimination-orders ((primary secondary elimination-order)
1346 &body body)
1347 "Evaluate BODY with primary and secondary elimination orders set to PRIMARY and SECONDARY."
[17]1348 `(let ((*primary-elimination-order* (or (find-order ,primary) *primary-elimination-order*))
1349 (*secondary-elimination-order* (or (find-order ,secondary) *secondary-elimination-order*))
1350 (*elimination-order* (or (find-order ,elimination-order) *elimination-order*)))
1351 . ,body))
[19]1352
[17]1353
[20]1354
[17]1355;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1356;;
[26]1357;; Conversion from internal form to Maxima general form
1358;;
[17]1359;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1360
[1]1361(defun maxima-head ()
1362 (if $poly_return_term_list
1363 '(mlist)
1364 '(mplus)))
1365
1366(defun coerce-to-maxima (poly-type object vars)
1367 (case poly-type
1368 (:polynomial
1369 `(,(maxima-head) ,@(mapcar #'(lambda (term) (coerce-to-maxima :term term vars)) (poly-termlist object))))
1370 (:poly-list
1371 `((mlist) ,@(mapcar #'(lambda (p) ($ratdisrep (coerce-to-maxima :polynomial p vars))) object)))
1372 (:term
1373 `((mtimes) ,($ratdisrep (term-coeff object))
1374 ,@(mapcar #'(lambda (var power) `((mexpt) ,var ,power))
1375 vars (monom-exponents (term-monom object)))))
1376 ;; Assumes that Lisp and Maxima logicals coincide
1377 (:logical object)
1378 (otherwise
1379 object)))
1380
1381
1382
1383;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1384;;
1385;; Macro facility for writing Maxima-level wrappers for
1386;; functions operating on internal representation
1387;;
1388;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1389
1390(defmacro with-parsed-polynomials (((maxima-vars &optional (maxima-new-vars nil new-vars-supplied-p))
1391 &key (polynomials nil)
1392 (poly-lists nil)
1393 (poly-list-lists nil)
1394 (value-type nil))
1395 &body body
1396 &aux (vars (gensym))
1397 (new-vars (gensym)))
1398 `(let ((,vars (coerce-maxima-list ,maxima-vars))
1399 ,@(when new-vars-supplied-p
1400 (list `(,new-vars (coerce-maxima-list ,maxima-new-vars)))))
1401 (coerce-to-maxima
1402 ,value-type
1403 (with-coefficient-ring ($poly_coefficient_ring)
1404 (with-monomial-order ($poly_monomial_order)
1405 (with-elimination-orders ($poly_primary_elimination_order
1406 $poly_secondary_elimination_order
1407 $poly_elimination_order)
[18]1408 (let ,(let ((args nil))
[1]1409 (dolist (p polynomials args)
1410 (setf args (cons `(,p (parse-poly ,p ,vars)) args)))
1411 (dolist (p poly-lists args)
1412 (setf args (cons `(,p (parse-poly-list ,p ,vars)) args)))
1413 (dolist (p poly-list-lists args)
1414 (setf args (cons `(,p (parse-poly-list-list ,p ,vars)) args))))
1415 . ,body))))
1416 ,(if new-vars-supplied-p
1417 `(append ,vars ,new-vars)
1418 vars))))
[18]1419
[1]1420(defmacro define-unop (maxima-name fun-name
1421 &optional (documentation nil documentation-supplied-p))
1422 "Define a MAXIMA-level unary operator MAXIMA-NAME corresponding to unary function FUN-NAME."
1423 `(defun ,maxima-name (p vars
1424 &aux
1425 (vars (coerce-maxima-list vars))
1426 (p (parse-poly p vars)))
1427 ,@(when documentation-supplied-p (list documentation))
1428 (coerce-to-maxima :polynomial (,fun-name *maxima-ring* p) vars)))
1429
1430(defmacro define-binop (maxima-name fun-name
1431 &optional (documentation nil documentation-supplied-p))
1432 "Define a MAXIMA-level binary operator MAXIMA-NAME corresponding to binary function FUN-NAME."
1433 `(defmfun ,maxima-name (p q vars
1434 &aux
1435 (vars (coerce-maxima-list vars))
1436 (p (parse-poly p vars))
1437 (q (parse-poly q vars)))
1438 ,@(when documentation-supplied-p (list documentation))
1439 (coerce-to-maxima :polynomial (,fun-name *maxima-ring* p q) vars)))
1440
1441
1442
1443;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1444;;
1445;; Maxima-level interface functions
1446;;
1447;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1448
1449;; Auxillary function for removing zero polynomial
1450(defun remzero (plist) (remove #'poly-zerop plist))
1451
1452;;Simple operators
1453
1454(define-binop $poly_add poly-add
1455 "Adds two polynomials P and Q")
1456
1457(define-binop $poly_subtract poly-sub
1458 "Subtracts a polynomial Q from P.")
1459
1460(define-binop $poly_multiply poly-mul
1461 "Returns the product of polynomials P and Q.")
1462
1463(define-binop $poly_s_polynomial spoly
1464 "Returns the syzygy polynomial (S-polynomial) of two polynomials P and Q.")
1465
1466(define-unop $poly_primitive_part poly-primitive-part
1467 "Returns the polynomial P divided by GCD of its coefficients.")
1468
1469(define-unop $poly_normalize poly-normalize
1470 "Returns the polynomial P divided by the leading coefficient.")
1471
1472;;Functions
1473
1474(defmfun $poly_expand (p vars)
1475 "This function is equivalent to EXPAND(P) if P parses correctly to a polynomial.
1476If the representation is not compatible with a polynomial in variables VARS,
1477the result is an error."
1478 (with-parsed-polynomials ((vars) :polynomials (p)
1479 :value-type :polynomial)
1480 p))
1481
1482(defmfun $poly_expt (p n vars)
1483 (with-parsed-polynomials ((vars) :polynomials (p) :value-type :polynomial)
1484 (poly-expt *maxima-ring* p n)))
1485
1486(defmfun $poly_content (p vars)
1487 (with-parsed-polynomials ((vars) :polynomials (p))
1488 (poly-content *maxima-ring* p)))
1489
1490(defmfun $poly_pseudo_divide (f fl vars
1491 &aux (vars (coerce-maxima-list vars))
[29]1492 (f (parse-poly f vars))
[1]1493 (fl (parse-poly-list fl vars)))
1494 (multiple-value-bind (quot rem c division-count)
1495 (poly-pseudo-divide *maxima-ring* f fl)
1496 `((mlist)
1497 ,(coerce-to-maxima :poly-list quot vars)
1498 ,(coerce-to-maxima :polynomial rem vars)
1499 ,c
1500 ,division-count)))
1501
1502(defmfun $poly_exact_divide (f g vars)
1503 (with-parsed-polynomials ((vars) :polynomials (f g) :value-type :polynomial)
1504 (poly-exact-divide *maxima-ring* f g)))
1505
1506(defmfun $poly_normal_form (f fl vars)
1507 (with-parsed-polynomials ((vars) :polynomials (f)
1508 :poly-lists (fl)
1509 :value-type :polynomial)
1510 (normal-form *maxima-ring* f (remzero fl) nil)))
1511
1512(defmfun $poly_buchberger_criterion (g vars)
1513 (with-parsed-polynomials ((vars) :poly-lists (g) :value-type :logical)
1514 (buchberger-criterion *maxima-ring* g)))
1515
1516(defmfun $poly_buchberger (fl vars)
1517 (with-parsed-polynomials ((vars) :poly-lists (fl) :value-type :poly-list)
1518 (buchberger *maxima-ring* (remzero fl) 0 nil)))
1519
1520(defmfun $poly_reduction (plist vars)
1521 (with-parsed-polynomials ((vars) :poly-lists (plist)
1522 :value-type :poly-list)
1523 (reduction *maxima-ring* plist)))
1524
1525(defmfun $poly_minimization (plist vars)
1526 (with-parsed-polynomials ((vars) :poly-lists (plist)
1527 :value-type :poly-list)
1528 (minimization plist)))
1529
1530(defmfun $poly_normalize_list (plist vars)
1531 (with-parsed-polynomials ((vars) :poly-lists (plist)
1532 :value-type :poly-list)
1533 (poly-normalize-list *maxima-ring* plist)))
1534
1535(defmfun $poly_grobner (f vars)
1536 (with-parsed-polynomials ((vars) :poly-lists (f)
1537 :value-type :poly-list)
1538 (grobner *maxima-ring* (remzero f))))
1539
1540(defmfun $poly_reduced_grobner (f vars)
1541 (with-parsed-polynomials ((vars) :poly-lists (f)
1542 :value-type :poly-list)
1543 (reduced-grobner *maxima-ring* (remzero f))))
1544
1545(defmfun $poly_depends_p (p var mvars
1546 &aux (vars (coerce-maxima-list mvars))
1547 (pos (position var vars)))
1548 (if (null pos)
1549 (merror "~%Variable ~M not in the list of variables ~M." var mvars)
1550 (poly-depends-p (parse-poly p vars) pos)))
1551
1552(defmfun $poly_elimination_ideal (flist k vars)
1553 (with-parsed-polynomials ((vars) :poly-lists (flist)
1554 :value-type :poly-list)
1555 (elimination-ideal *maxima-ring* flist k nil 0)))
1556
1557(defmfun $poly_colon_ideal (f g vars)
1558 (with-parsed-polynomials ((vars) :poly-lists (f g) :value-type :poly-list)
1559 (colon-ideal *maxima-ring* f g nil)))
1560
1561(defmfun $poly_ideal_intersection (f g vars)
1562 (with-parsed-polynomials ((vars) :poly-lists (f g) :value-type :poly-list)
1563 (ideal-intersection *maxima-ring* f g nil)))
1564
1565(defmfun $poly_lcm (f g vars)
1566 (with-parsed-polynomials ((vars) :polynomials (f g) :value-type :polynomial)
1567 (poly-lcm *maxima-ring* f g)))
1568
1569(defmfun $poly_gcd (f g vars)
1570 ($first ($divide (m* f g) ($poly_lcm f g vars))))
1571
1572(defmfun $poly_grobner_equal (g1 g2 vars)
1573 (with-parsed-polynomials ((vars) :poly-lists (g1 g2))
1574 (grobner-equal *maxima-ring* g1 g2)))
1575
1576(defmfun $poly_grobner_subsetp (g1 g2 vars)
1577 (with-parsed-polynomials ((vars) :poly-lists (g1 g2))
1578 (grobner-subsetp *maxima-ring* g1 g2)))
1579
1580(defmfun $poly_grobner_member (p g vars)
1581 (with-parsed-polynomials ((vars) :polynomials (p) :poly-lists (g))
1582 (grobner-member *maxima-ring* p g)))
1583
1584(defmfun $poly_ideal_saturation1 (f p vars)
1585 (with-parsed-polynomials ((vars) :poly-lists (f) :polynomials (p)
1586 :value-type :poly-list)
1587 (ideal-saturation-1 *maxima-ring* f p 0)))
1588
1589(defmfun $poly_saturation_extension (f plist vars new-vars)
1590 (with-parsed-polynomials ((vars new-vars)
1591 :poly-lists (f plist)
1592 :value-type :poly-list)
1593 (saturation-extension *maxima-ring* f plist)))
1594
1595(defmfun $poly_polysaturation_extension (f plist vars new-vars)
[26]1596 (with-parsed-polynomials ((vars new-vars)
1597 :poly-lists (f plist)
1598 :value-type :poly-list)
1599 (polysaturation-extension *maxima-ring* f plist)))
1600
1601(defmfun $poly_ideal_polysaturation1 (f plist vars)
1602 (with-parsed-polynomials ((vars) :poly-lists (f plist)
1603 :value-type :poly-list)
1604 (ideal-polysaturation-1 *maxima-ring* f plist 0 nil)))
1605
1606(defmfun $poly_ideal_saturation (f g vars)
1607 (with-parsed-polynomials ((vars) :poly-lists (f g)
1608 :value-type :poly-list)
1609 (ideal-saturation *maxima-ring* f g 0 nil)))
1610
1611(defmfun $poly_ideal_polysaturation (f ideal-list vars)
1612 (with-parsed-polynomials ((vars) :poly-lists (f)
1613 :poly-list-lists (ideal-list)
1614 :value-type :poly-list)
1615 (ideal-polysaturation *maxima-ring* f ideal-list 0 nil)))
1616
1617(defmfun $poly_lt (f vars)
1618 (with-parsed-polynomials ((vars) :polynomials (f) :value-type :polynomial)
1619 (make-poly-from-termlist (list (poly-lt f)))))
1620
1621(defmfun $poly_lm (f vars)
1622 (with-parsed-polynomials ((vars) :polynomials (f) :value-type :polynomial)
1623 (make-poly-from-termlist (list (make-term (poly-lm f) (funcall (ring-unit *maxima-ring*)))))))
1624
Note: See TracBrowser for help on using the repository browser.