Context Navigation

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@ 55

Last change on this file since 55 was 54, checked in by Marek Rychlik, 10 years ago
* empty log message *
File size: 58.3 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: