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

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
	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 C2F-C1M*G, where F and G are polynomials M is a monomial.
[42]	312	Assume 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
	325	multiple values. The first value is a list of quotients A. The second
	326	value is the remainder R. The third argument is a scalar coefficient
	327	C, such that C*F can be divided by FL within the ring of coefficients,
	328	which is not necessarily a field. Finally, the fourth value is an
	329	integer count of the number of reductions performed. The resulting
	330	objects satisfy the equation: CF= 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 cf=sum aifi+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
	368	with 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 cf=sum aifi+r+p by cg.
	399	(setf r (scalar-times-poly ring cg r)
	400	c (funcall (ring-mul ring) c cg)
	401	;; p := cgp-cpm*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: cf0=sum aifi+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
	445	criterion: for every two polynomials h1 and h2 in G the S-polynomial
	446	S(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
	455	upon 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
	480	polynomials 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
	510	in 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.
	530	START is the first position beyond the elements which form a partial
	531	grobner 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
	556	of the ideal generated by the polynomial list F. Polynomials 0 to
	557	START-1 are assumed to be a Grobner basis already, so that certain
	558	critical pairs will not be examined. If TOP-REDUCTION-ONLY set, top
	559	reduction will be preformed. This function assumes that all polynomials
	560	in 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
	607	of the ideal generated by the polynomial list F. Polynomials 0 to
	608	START-1 are assumed to be a Grobner basis already, so that certain
	609	critical pairs will not be examined. If TOP-REDUCTION-ONLY set, top
	610	reduction 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
	695	in G pointed to by the integers in PAIR have disjoint (relatively prime)
	696	monomials. 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
	708	PARTIAL-BASIS divides the LCM of the leading monomials of the two
	709	polynomials in the polynomial list PARTIAL-BASIS, and P paired with
	710	each of the polynomials pointed to by the the PAIR has already been
	711	treated, 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
	751	Moeller. This algorithm is described as BUCHBERGERNEW2 in the book by
	752	Becker-Weispfenning entitled ``Grobner Bases''. This function assumes
	753	that 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
	795	Gebauer-Moeller algorithm. G is a list of polynomials, B is a list of
	796	critical pairs and H is a new polynomial which possibly will be added
	797	to G. The naming conventions used are very close to the one used in
	798	the 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
	861	the polynomials is divisible by a leading monomial of another
	862	polynomial. 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
	882	monomial ideal as P but minimal, i.e. no leading monomial
	883	of a polynomial in the sublist divides the leading monomial
	884	of 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
	898	that the division is possible, which may not always be the
	899	case in rings which are not fields. The exact division operator
	900	is assumed to be provided by the RING structure of the
	901	COEFFICIENT-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
	921	names. Names currently used are either Lisp symbols, Maxima symbols or
	922	keywords."
	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
	939	to 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
	981	and it calculates the intersection with the ring R[x[k+1],...,x[n]], i.e.
	982	it 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),
	998	where F and G are two lists of polynomials. The colon ideal I:J is
	999	defined as the set of polynomials H such that for all polynomials W in
	1000	J 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
	1021	F 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.
	1044	The polynomials must be ordered according to monomial order PRED
	1045	and their coefficients must be compatible with the RING structure
	1046	defined 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,
	1066	generate 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
	1071	an ideal contained in the ideal generated by a polynomial list G2,
	1072	both 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
	1077	polynomial 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
	1085	generated by a polynomial list F in the ideal generated by a single
	1086	polynomial P. The saturation ideal is defined as the set of
	1087	polynomials H such for some natural number n (* (EXPT P N) H) is in the ideal
	1088	F. Geometrically, over an algebraically closed field, this is the set
	1089	of polynomials in the ideal generated by F which do not identically
	1090	vanish 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
	1105	sequence of successive saturations in the polynomials
	1106	of the polynomial list PLIST of the ideal generated by the
	1107	polynomial 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
	1119	generated by a polynomial list F in the ideal generated a polynomial
	1120	list G. The saturation ideal is defined as the set of polynomials H
	1121	such for some natural number n and some P in the ideal generated by G
	1122	the polynomial P*N H is in the ideal spanned by F. Geometrically,
	1123	over an algebraically closed field, this is the set of polynomials in
	1124	the ideal generated by F which do not identically vanish on the
	1125	variety 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
	1137	successive applications of IDEAL-SATURATION to F and lists of
	1138	polynomials 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
	1229	are 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.
	1476	If the representation is not compatible with a polynomial in variables VARS,
	1477	the 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.

Download in other formats: