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: trunk/grobner.lisp@ 777

Last change on this file since 777 was 637, checked in by Marek Rychlik, 9 years ago
File size: 79.1 KB
Line 
1;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*-
2;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3;;;
4;;; Copyright (C) 1999, 2002, 2009 Marek Rychlik <rychlik@u.arizona.edu>
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~%"
30 "$Revision: 1.1 $" "$Date: 2008/09/08 21:40:10 $"))
31
32;;FUNCTS is loaded because it contains the definition of LCM
33($load "functs")
34
35;; Macros for making lists with iterators - an exammple of GENSYM
36;; MAKELIST-1 makes a list with one iterator, while MAKELIST accepts an
37;; arbitrary number of iterators
38
39;; Sample usage:
40;; Without a step:
41;; >(makelist-1 (* 2 i) i 0 10)
42;; (0 2 4 6 8 10 12 14 16 18 20)
43;; With a step of 3:
44;; >(makelist-1 (* 2 i) i 0 10 3)
45;; (0 6 12 18)
46
47;; Generate sums of squares of numbers between 1 and 4:
48;; >(makelist (+ (* i i) (* j j)) (i 1 4) (j 1 i))
49;; (2 5 8 10 13 18 17 20 25 32)
50;; >(makelist (list i j '---> (+ (* i i) (* j j))) (i 1 4) (j 1 i))
51;; ((1 1 ---> 2) (2 1 ---> 5) (2 2 ---> 8) (3 1 ---> 10) (3 2 ---> 13)
52;; (3 3 ---> 18) (4 1 ---> 17) (4 2 ---> 20) (4 3 ---> 25) (4 4 ---> 32))
53
54;; Evaluate expression expr with variable set to lo, lo+1,... ,hi
55;; and put the results in a list.
56(defmacro makelist-1 (expr var lo hi &optional (step 1))
57 (let ((l (gensym)))
58 `(do ((,var ,lo (+ ,var ,step))
59 (,l nil (cons ,expr ,l)))
60 ((> ,var ,hi) (reverse ,l))
61 (declare (fixnum ,var)))))
62
63(defmacro makelist (expr (var lo hi &optional (step 1)) &rest more)
64 (if (endp more)
65 `(makelist-1 ,expr ,var ,lo ,hi ,step)
66 (let* ((l (gensym)))
67 `(do ((,var ,lo (+ ,var ,step))
68 (,l nil (nconc ,l `,(makelist ,expr ,@more))))
69 ((> ,var ,hi) ,l)
70 (declare (fixnum ,var))))))
71
72;;----------------------------------------------------------------
73;; This package implements BASIC OPERATIONS ON MONOMIALS
74;;----------------------------------------------------------------
75;; DATA STRUCTURES: Monomials are represented as lists:
76;;
77;; monom: (n1 n2 ... nk) where ni are non-negative integers
78;;
79;; However, lists may be implemented as other sequence types,
80;; so the flexibility to change the representation should be
81;; maintained in the code to use general operations on sequences
82;; whenever possible. The optimization for the actual representation
83;; should be left to declarations and the compiler.
84;;----------------------------------------------------------------
85;; EXAMPLES: Suppose that variables are x and y. Then
86;;
87;; Monom x*y^2 ---> (1 2)
88;;
89;;----------------------------------------------------------------
90
91(deftype exponent ()
92 "Type of exponent in a monomial."
93 'fixnum)
94
95(deftype monom (&optional dim)
96 "Type of monomial."
97 `(simple-array exponent (,dim)))
98
99;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
100;;
101;; Construction of monomials
102;;
103;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
104
105(defmacro make-monom (dim &key (initial-contents nil initial-contents-supplied-p)
106 (initial-element 0 initial-element-supplied-p))
107 "Make a monomial with DIM variables. Additional argument
108INITIAL-CONTENTS specifies the list of powers of the consecutive
109variables. The alternative additional argument INITIAL-ELEMENT
110specifies the common power for all variables."
111 ;;(declare (fixnum dim))
112 `(make-array ,dim
113 :element-type 'exponent
114 ,@(when initial-contents-supplied-p `(:initial-contents ,initial-contents))
115 ,@(when initial-element-supplied-p `(:initial-element ,initial-element))))
116
117
118
119;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
120;;
121;; Operations on monomials
122;;
123;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
124
125(defmacro monom-elt (m index)
126 "Return the power in the monomial M of variable number INDEX."
127 `(elt ,m ,index))
128
129(defun monom-dimension (m)
130 "Return the number of variables in the monomial M."
131 (length m))
132
133(defun monom-total-degree (m &optional (start 0) (end (length m)))
134 "Return the todal degree of a monomoal M. Optinally, a range
135of variables may be specified with arguments START and END."
136 (declare (type monom m) (fixnum start end))
137 (reduce #'+ m :start start :end end))
138
139(defun monom-sugar (m &aux (start 0) (end (length m)))
140 "Return the sugar of a monomial M. Optinally, a range
141of variables may be specified with arguments START and END."
142 (declare (type monom m) (fixnum start end))
143 (monom-total-degree m start end))
144
145(defun monom-div (m1 m2 &aux (result (copy-seq m1)))
146 "Divide monomial M1 by monomial M2."
147 (declare (type monom m1 m2 result))
148 (map-into result #'- m1 m2))
149
150(defun monom-mul (m1 m2 &aux (result (copy-seq m1)))
151 "Multiply monomial M1 by monomial M2."
152 (declare (type monom m1 m2 result))
153 (map-into result #'+ m1 m2))
154
155(defun monom-divides-p (m1 m2)
156 "Returns T if monomial M1 divides monomial M2, NIL otherwise."
157 (declare (type monom m1 m2))
158 (every #'<= m1 m2))
159
160(defun monom-divides-monom-lcm-p (m1 m2 m3)
161 "Returns T if monomial M1 divides MONOM-LCM(M2,M3), NIL otherwise."
162 (declare (type monom m1 m2 m3))
163 (every #'(lambda (x y z) (declare (type exponent x y z)) (<= x (max y z))) m1 m2 m3))
164
165(defun monom-lcm-divides-monom-lcm-p (m1 m2 m3 m4)
166 "Returns T if monomial MONOM-LCM(M1,M2) divides MONOM-LCM(M3,M4), NIL otherwise."
167 (declare (type monom m1 m2 m3 m4))
168 (every #'(lambda (x y z w) (declare (type exponent x y z w)) (<= (max x y) (max z w))) m1 m2 m3 m4))
169
170(defun monom-lcm-equal-monom-lcm-p (m1 m2 m3 m4)
171 "Returns T if monomial MONOM-LCM(M1,M2) equals MONOM-LCM(M3,M4), NIL otherwise."
172 (declare (type monom m1 m2 m3 m4))
173 (every #'(lambda (x y z w) (declare (type exponent x y z w)) (= (max x y) (max z w))) m1 m2 m3 m4))
174
175(defun monom-divisible-by-p (m1 m2)
176 "Returns T if monomial M1 is divisible by monomial M2, NIL otherwise."
177 (declare (type monom m1 m2))
178 (every #'>= m1 m2))
179
180(defun monom-rel-prime-p (m1 m2)
181 "Returns T if two monomials M1 and M2 are relatively prime (disjoint)."
182 (declare (type monom m1 m2))
183 (every #'(lambda (x y) (declare (type exponent x y)) (zerop (min x y))) m1 m2))
184
185(defun monom-equal-p (m1 m2)
186 "Returns T if two monomials M1 and M2 are equal."
187 (declare (type monom m1 m2))
188 (every #'= m1 m2))
189
190(defun monom-lcm (m1 m2 &aux (result (copy-seq m1)))
191 "Returns least common multiple of monomials M1 and M2."
192 (declare (type monom m1 m2))
193 (map-into result #'max m1 m2))
194
195(defun monom-gcd (m1 m2 &aux (result (copy-seq m1)))
196 "Returns greatest common divisor of monomials M1 and M2."
197 (declare (type monom m1 m2))
198 (map-into result #'min m1 m2))
199
200(defun monom-depends-p (m k)
201 "Return T if the monomial M depends on variable number K."
202 (declare (type monom m) (fixnum k))
203 (plusp (elt m k)))
204
205(defmacro monom-map (fun m &rest ml &aux (result `(copy-seq ,m)))
206 `(map-into ,result ,fun ,m ,@ml))
207
208(defmacro monom-append (m1 m2)
209 `(concatenate 'monom ,m1 ,m2))
210
211(defmacro monom-contract (k m)
212 `(subseq ,m ,k))
213
214(defun monom-exponents (m)
215 (declare (type monom m))
216 (coerce m 'list))
217
218
219
220;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
221;;
222;; Implementations of various admissible monomial orders
223;;
224;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
225
226;; pure lexicographic
227(defun lex> (p q &optional (start 0) (end (monom-dimension p)))
228 "Return T if P>Q with respect to lexicographic order, otherwise NIL.
229The second returned value is T if P=Q, otherwise it is NIL."
230 (declare (type monom p q) (type fixnum start end))
231 (do ((i start (1+ i)))
232 ((>= i end) (values nil t))
233 (declare (type fixnum i))
234 (cond
235 ((> (monom-elt p i) (monom-elt q i))
236 (return-from lex> (values t nil)))
237 ((< (monom-elt p i) (monom-elt q i))
238 (return-from lex> (values nil nil))))))
239
240;; total degree order , ties broken by lexicographic
241(defun grlex> (p q &optional (start 0) (end (monom-dimension p)))
242 "Return T if P>Q with respect to graded lexicographic order, otherwise NIL.
243The second returned value is T if P=Q, otherwise it is NIL."
244 (declare (type monom p q) (type fixnum start end))
245 (let ((d1 (monom-total-degree p start end))
246 (d2 (monom-total-degree q start end)))
247 (cond
248 ((> d1 d2) (values t nil))
249 ((< d1 d2) (values nil nil))
250 (t
251 (lex> p q start end)))))
252
253
254;; total degree, ties broken by reverse lexicographic
255(defun grevlex> (p q &optional (start 0) (end (monom-dimension p)))
256 "Return T if P>Q with respect to graded reverse lexicographic order,
257NIL otherwise. The second returned value is T if P=Q, otherwise it is NIL."
258 (declare (type monom p q) (type fixnum start end))
259 (let ((d1 (monom-total-degree p start end))
260 (d2 (monom-total-degree q start end)))
261 (cond
262 ((> d1 d2) (values t nil))
263 ((< d1 d2) (values nil nil))
264 (t
265 (revlex> p q start end)))))
266
267
268;; reverse lexicographic
269(defun revlex> (p q &optional (start 0) (end (monom-dimension p)))
270 "Return T if P>Q with respect to reverse lexicographic order, NIL
271otherwise. The second returned value is T if P=Q, otherwise it is
272NIL. This is not and admissible monomial order because some sets do
273not have a minimal element. This order is useful in constructing other
274orders."
275 (declare (type monom p q) (type fixnum start end))
276 (do ((i (1- end) (1- i)))
277 ((< i start) (values nil t))
278 (declare (type fixnum i))
279 (cond
280 ((< (monom-elt p i) (monom-elt q i))
281 (return-from revlex> (values t nil)))
282 ((> (monom-elt p i) (monom-elt q i))
283 (return-from revlex> (values nil nil))))))
284
285
286(defun invlex> (p q &optional (start 0) (end (monom-dimension p)))
287 "Return T if P>Q with respect to inverse lexicographic order, NIL otherwise
288The second returned value is T if P=Q, otherwise it is NIL."
289 (declare (type monom p q) (type fixnum start end))
290 (do ((i (1- end) (1- i)))
291 ((< i start) (values nil t))
292 (declare (type fixnum i))
293 (cond
294 ((> (monom-elt p i) (monom-elt q i))
295 (return-from invlex> (values t nil)))
296 ((< (monom-elt p i) (monom-elt q i))
297 (return-from invlex> (values nil nil))))))
298
299
300
301;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
302;;
303;; Order making functions
304;;
305;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
306
307(defvar *monomial-order* #'lex>
308 "Default order for monomial comparisons")
309
310(defmacro monomial-order (x y)
311 `(funcall *monomial-order* ,x ,y))
312
313(defun reverse-monomial-order (x y)
314 (monomial-order y x))
315
316(defvar *primary-elimination-order* #'lex>)
317
318(defvar *secondary-elimination-order* #'lex>)
319
320(defvar *elimination-order* nil
321 "Default elimination order used in elimination-based functions.
322If not NIL, it is assumed to be a proper elimination order. If NIL,
323we will construct an elimination order using the values of
324*PRIMARY-ELIMINATION-ORDER* and *SECONDARY-ELIMINATION-ORDER*.")
325
326(defun elimination-order (k)
327 "Return a predicate which compares monomials according to the
328K-th elimination order. Two variables *PRIMARY-ELIMINATION-ORDER*
329and *SECONDARY-ELIMINATION-ORDER* control the behavior on the first K
330and the remaining variables, respectively."
331 (declare (type fixnum k))
332 #'(lambda (p q &optional (start 0) (end (monom-dimension p)))
333 (declare (type monom p q) (type fixnum start end))
334 (multiple-value-bind (primary equal)
335 (funcall *primary-elimination-order* p q start k)
336 (if equal
337 (funcall *secondary-elimination-order* p q k end)
338 (values primary nil)))))
339
340(defun elimination-order-1 (p q &optional (start 0) (end (monom-dimension p)))
341 "Equivalent to the function returned by the call to (ELIMINATION-ORDER 1)."
342 (declare (type monom p q) (type fixnum start end))
343 (cond
344 ((> (monom-elt p start) (monom-elt q start)) (values t nil))
345 ((< (monom-elt p start) (monom-elt q start)) (values nil nil))
346 (t (funcall *secondary-elimination-order* p q (1+ start) end))))
347
348
349
350;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
351;;
352;; Priority queue stuff
353;;
354;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
355
356(defparameter *priority-queue-allocation-size* 16)
357
358(defun priority-queue-make-heap (&key (element-type 'fixnum))
359 (make-array *priority-queue-allocation-size* :element-type element-type :fill-pointer 1
360 :adjustable t))
361
362(defstruct (priority-queue (:constructor priority-queue-construct))
363 (heap (priority-queue-make-heap))
364 test)
365
366(defun make-priority-queue (&key (element-type 'fixnum)
367 (test #'<=)
368 (element-key #'identity))
369 (priority-queue-construct
370 :heap (priority-queue-make-heap :element-type element-type)
371 :test #'(lambda (x y) (funcall test (funcall element-key y) (funcall element-key x)))))
372
373(defun priority-queue-insert (pq item)
374 (priority-queue-heap-insert (priority-queue-heap pq) item (priority-queue-test pq)))
375
376(defun priority-queue-remove (pq)
377 (priority-queue-heap-remove (priority-queue-heap pq) (priority-queue-test pq)))
378
379(defun priority-queue-empty-p (pq)
380 (priority-queue-heap-empty-p (priority-queue-heap pq)))
381
382(defun priority-queue-size (pq)
383 (fill-pointer (priority-queue-heap pq)))
384
385(defun priority-queue-upheap (a k
386 &optional
387 (test #'<=)
388 &aux (v (aref a k)))
389 (declare (fixnum k))
390 (assert (< 0 k (fill-pointer a)))
391 (loop
392 (let ((parent (ash k -1)))
393 (when (zerop parent) (return))
394 (unless (funcall test (aref a parent) v)
395 (return))
396 (setf (aref a k) (aref a parent)
397 k parent)))
398 (setf (aref a k) v)
399 a)
400
401
402(defun priority-queue-heap-insert (a item &optional (test #'<=))
403 (vector-push-extend item a)
404 (priority-queue-upheap a (1- (fill-pointer a)) test))
405
406(defun priority-queue-downheap (a k
407 &optional
408 (test #'<=)
409 &aux (v (aref a k)) (j 0) (n (fill-pointer a)))
410 (declare (fixnum k n j))
411 (loop
412 (unless (<= k (ash n -1))
413 (return))
414 (setf j (ash k 1))
415 (if (and (< j n) (not (funcall test (aref a (1+ j)) (aref a j))))
416 (incf j))
417 (when (funcall test (aref a j) v)
418 (return))
419 (setf (aref a k) (aref a j)
420 k j))
421 (setf (aref a k) v)
422 a)
423
424(defun priority-queue-heap-remove (a &optional (test #'<=) &aux (v (aref a 1)))
425 (when (<= (fill-pointer a) 1) (error "Empty queue."))
426 (setf (aref a 1) (vector-pop a))
427 (priority-queue-downheap a 1 test)
428 (values v a))
429
430(defun priority-queue-heap-empty-p (a)
431 (<= (fill-pointer a) 1))
432
433
434
435;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
436;;
437;; Global switches
438;; (Can be used in Maxima just fine)
439;;
440;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
441
442(defmvar $poly_monomial_order '$lex
443 "This switch controls which monomial order is used in polynomial
444and Grobner basis calculations. If not set, LEX will be used")
445
446(defmvar $poly_coefficient_ring '$expression_ring
447 "This switch indicates the coefficient ring of the polynomials
448that will be used in grobner calculations. If not set, Maxima's
449general expression ring will be used. This variable may be set
450to RING_OF_INTEGERS if desired.")
451
452(defmvar $poly_primary_elimination_order nil
453 "Name of the default order for eliminated variables in elimination-based functions.
454If not set, LEX will be used.")
455
456(defmvar $poly_secondary_elimination_order nil
457 "Name of the default order for kept variables in elimination-based functions.
458If not set, LEX will be used.")
459
460(defmvar $poly_elimination_order nil
461 "Name of the default elimination order used in elimination calculations.
462If set, it overrides the settings in variables POLY_PRIMARY_ELIMINATION_ORDER
463and SECONDARY_ELIMINATION_ORDER. The user must ensure that this is a true
464elimination order valid for the number of eliminated variables.")
465
466(defmvar $poly_return_term_list nil
467 "If set to T, all functions in this package will return each polynomial as a
468list of terms in the current monomial order rather than a Maxima general expression.")
469
470(defmvar $poly_grobner_debug nil
471 "If set to TRUE, produce debugging and tracing output.")
472
473(defmvar $poly_grobner_algorithm '$buchberger
474 "The name of the algorithm used to find grobner bases.")
475
476(defmvar $poly_top_reduction_only nil
477 "If not FALSE, use top reduction only whenever possible.
478Top reduction means that division algorithm stops after the first reduction.")
479
480
481
482;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
483;;
484;; Coefficient ring operations
485;;
486;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
487;;
488;; These are ALL operations that are performed on the coefficients by
489;; the package, and thus the coefficient ring can be changed by merely
490;; redefining these operations.
491;;
492;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
493
494(defstruct (ring)
495 (parse #'identity :type function)
496 (unit #'identity :type function)
497 (zerop #'identity :type function)
498 (add #'identity :type function)
499 (sub #'identity :type function)
500 (uminus #'identity :type function)
501 (mul #'identity :type function)
502 (div #'identity :type function)
503 (lcm #'identity :type function)
504 (ezgcd #'identity :type function)
505 (gcd #'identity :type function))
506
507(defparameter *ring-of-integers*
508 (make-ring
509 :parse #'identity
510 :unit #'(lambda () 1)
511 :zerop #'zerop
512 :add #'+
513 :sub #'-
514 :uminus #'-
515 :mul #'*
516 :div #'/
517 :lcm #'lcm
518 :ezgcd #'(lambda (x y &aux (c (gcd x y))) (values c (/ x c) (/ y c)))
519 :gcd #'gcd)
520 "The ring of integers.")
521
522
523
524;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
525;;
526;; This is how we perform operations on coefficients
527;; using Maxima functions.
528;;
529;; Functions and macros dealing with internal representation structure
530;;
531;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
532
533(defun make-term-variable (ring nvars pos
534 &optional
535 (power 1)
536 (coeff (funcall (ring-unit ring)))
537 &aux
538 (monom (make-monom nvars :initial-element 0)))
539 (declare (fixnum nvars pos power))
540 (incf (monom-elt monom pos) power)
541 (make-term monom coeff))
542
543(defstruct (term
544 (:constructor make-term (monom coeff))
545 ;;(:constructor make-term-variable)
546 ;;(:type list)
547 )
548 (monom (make-monom 0) :type monom)
549 (coeff nil))
550
551(defun term-sugar (term)
552 (monom-sugar (term-monom term)))
553
554(defun termlist-sugar (p &aux (sugar -1))
555 (declare (fixnum sugar))
556 (dolist (term p sugar)
557 (setf sugar (max sugar (term-sugar term)))))
558
559
560
561
562;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
563;;
564;; Low-level polynomial arithmetic done on
565;; lists of terms
566;;
567;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
568
569(defmacro termlist-lt (p) `(car ,p))
570(defun termlist-lm (p) (term-monom (termlist-lt p)))
571(defun termlist-lc (p) (term-coeff (termlist-lt p)))
572
573(define-modify-macro scalar-mul (c) coeff-mul)
574
575(defun scalar-times-termlist (ring c p)
576 "Multiply scalar C by a polynomial P. This function works
577even if there are divisors of 0."
578 (mapcan
579 #'(lambda (term)
580 (let ((c1 (funcall (ring-mul ring) c (term-coeff term))))
581 (unless (funcall (ring-zerop ring) c1)
582 (list (make-term (term-monom term) c1)))))
583 p))
584
585
586(defun term-mul (ring term1 term2)
587 "Returns (LIST TERM) wheter TERM is the product of the terms TERM1 TERM2,
588or NIL when the product is 0. This definition takes care of divisors of 0
589in the coefficient ring."
590 (let ((c (funcall (ring-mul ring) (term-coeff term1) (term-coeff term2))))
591 (unless (funcall (ring-zerop ring) c)
592 (list (make-term (monom-mul (term-monom term1) (term-monom term2)) c)))))
593
594(defun term-times-termlist (ring term f)
595 (declare (type ring ring))
596 (mapcan #'(lambda (term-f) (term-mul ring term term-f)) f))
597
598(defun termlist-times-term (ring f term)
599 (mapcan #'(lambda (term-f) (term-mul ring term-f term)) f))
600
601(defun monom-times-term (m term)
602 (make-term (monom-mul m (term-monom term)) (term-coeff term)))
603
604(defun monom-times-termlist (m f)
605 (cond
606 ((null f) nil)
607 (t
608 (mapcar #'(lambda (x) (monom-times-term m x)) f))))
609
610(defun termlist-uminus (ring f)
611 (mapcar #'(lambda (x)
612 (make-term (term-monom x) (funcall (ring-uminus ring) (term-coeff x))))
613 f))
614
615(defun termlist-add (ring p q)
616 (declare (type list p q))
617 (do (r)
618 ((cond
619 ((endp p)
620 (setf r (revappend r q)) t)
621 ((endp q)
622 (setf r (revappend r p)) t)
623 (t
624 (multiple-value-bind
625 (lm-greater lm-equal)
626 (monomial-order (termlist-lm p) (termlist-lm q))
627 (cond
628 (lm-equal
629 (let ((s (funcall (ring-add ring) (termlist-lc p) (termlist-lc q))))
630 (unless (funcall (ring-zerop ring) s) ;check for cancellation
631 (setf r (cons (make-term (termlist-lm p) s) r)))
632 (setf p (cdr p) q (cdr q))))
633 (lm-greater
634 (setf r (cons (car p) r)
635 p (cdr p)))
636 (t (setf r (cons (car q) r)
637 q (cdr q)))))
638 nil))
639 r)))
640
641(defun termlist-sub (ring p q)
642 (declare (type list p q))
643 (do (r)
644 ((cond
645 ((endp p)
646 (setf r (revappend r (termlist-uminus ring q)))
647 t)
648 ((endp q)
649 (setf r (revappend r p))
650 t)
651 (t
652 (multiple-value-bind
653 (mgreater mequal)
654 (monomial-order (termlist-lm p) (termlist-lm q))
655 (cond
656 (mequal
657 (let ((s (funcall (ring-sub ring) (termlist-lc p) (termlist-lc q))))
658 (unless (funcall (ring-zerop ring) s) ;check for cancellation
659 (setf r (cons (make-term (termlist-lm p) s) r)))
660 (setf p (cdr p) q (cdr q))))
661 (mgreater
662 (setf r (cons (car p) r)
663 p (cdr p)))
664 (t (setf r (cons (make-term (termlist-lm q) (funcall (ring-uminus ring) (termlist-lc q))) r)
665 q (cdr q)))))
666 nil))
667 r)))
668
669;; Multiplication of polynomials
670;; Non-destructive version
671(defun termlist-mul (ring p q)
672 (cond ((or (endp p) (endp q)) nil) ;p or q is 0 (represented by NIL)
673 ;; If p=p0+p1 and q=q0+q1 then pq=p0q0+p0q1+p1q
674 ((endp (cdr p))
675 (term-times-termlist ring (car p) q))
676 ((endp (cdr q))
677 (termlist-times-term ring p (car q)))
678 (t
679 (let ((head (term-mul ring (termlist-lt p) (termlist-lt q)))
680 (tail (termlist-add ring (term-times-termlist ring (car p) (cdr q))
681 (termlist-mul ring (cdr p) q))))
682 (cond ((null head) tail)
683 ((null tail) head)
684 (t (nconc head tail)))))))
685
686(defun termlist-unit (ring dimension)
687 (declare (fixnum dimension))
688 (list (make-term (make-monom dimension :initial-element 0)
689 (funcall (ring-unit ring)))))
690
691(defun termlist-expt (ring poly n &aux (dim (monom-dimension (termlist-lm poly))))
692 (declare (type fixnum n dim))
693 (cond
694 ((minusp n) (error "termlist-expt: Negative exponent."))
695 ((endp poly) (if (zerop n) (termlist-unit ring dim) nil))
696 (t
697 (do ((k 1 (ash k 1))
698 (q poly (termlist-mul ring q q)) ;keep squaring
699 (p (termlist-unit ring dim) (if (not (zerop (logand k n))) (termlist-mul ring p q) p)))
700 ((> k n) p)
701 (declare (fixnum k))))))
702
703
704
705;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
706;;
707;; Additional structure operations on a list of terms
708;;
709;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
710
711(defun termlist-contract (p &optional (k 1))
712 "Eliminate first K variables from a polynomial P."
713 (mapcar #'(lambda (term) (make-term (monom-contract k (term-monom term))
714 (term-coeff term)))
715 p))
716
717(defun termlist-extend (p &optional (m (make-monom 1 :initial-element 0)))
718 "Extend every monomial in a polynomial P by inserting at the
719beginning of every monomial the list of powers M."
720 (mapcar #'(lambda (term) (make-term (monom-append m (term-monom term))
721 (term-coeff term)))
722 p))
723
724(defun termlist-add-variables (p n)
725 "Add N variables to a polynomial P by inserting zero powers
726at the beginning of each monomial."
727 (declare (fixnum n))
728 (mapcar #'(lambda (term)
729 (make-term (monom-append (make-monom n :initial-element 0)
730 (term-monom term))
731 (term-coeff term)))
732 p))
733
734
735
736;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
737;;
738;; Arithmetic on polynomials
739;;
740;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
741
742(defstruct (poly
743 ;;BOA constructor, by default constructs zero polynomial
744 (:constructor make-poly-from-termlist (termlist &optional (sugar (termlist-sugar termlist))))
745 (:constructor make-poly-zero (&aux (termlist nil) (sugar -1)))
746 ;;Constructor of polynomials representing a variable
747 (:constructor make-variable (ring nvars pos &optional (power 1)
748 &aux
749 (termlist (list
750 (make-term-variable ring nvars pos power)))
751 (sugar power)))
752 (:constructor poly-unit (ring dimension
753 &aux
754 (termlist (termlist-unit ring dimension))
755 (sugar 0))))
756 (termlist nil :type list)
757 (sugar -1 :type fixnum))
758
759;; Leading term
760(defmacro poly-lt (p) `(car (poly-termlist ,p)))
761
762;; Second term
763(defmacro poly-second-lt (p) `(cadar (poly-termlist ,p)))
764
765;; Leading monomial
766(defun poly-lm (p) (term-monom (poly-lt p)))
767
768;; Second monomial
769(defun poly-second-lm (p) (term-monom (poly-second-lt p)))
770
771;; Leading coefficient
772(defun poly-lc (p) (term-coeff (poly-lt p)))
773
774;; Second coefficient
775(defun poly-second-lc (p) (term-coeff (poly-second-lt p)))
776
777;; Testing for a zero polynomial
778(defun poly-zerop (p) (null (poly-termlist p)))
779
780;; The number of terms
781(defun poly-length (p) (length (poly-termlist p)))
782
783(defun scalar-times-poly (ring c p)
784 (make-poly-from-termlist (scalar-times-termlist ring c (poly-termlist p)) (poly-sugar p)))
785
786(defun monom-times-poly (m p)
787 (make-poly-from-termlist (monom-times-termlist m (poly-termlist p)) (+ (poly-sugar p) (monom-sugar m))))
788
789(defun term-times-poly (ring term p)
790 (make-poly-from-termlist (term-times-termlist ring term (poly-termlist p)) (+ (poly-sugar p) (term-sugar term))))
791
792(defun poly-add (ring p q)
793 (make-poly-from-termlist (termlist-add ring (poly-termlist p) (poly-termlist q)) (max (poly-sugar p) (poly-sugar q))))
794
795(defun poly-sub (ring p q)
796 (make-poly-from-termlist (termlist-sub ring (poly-termlist p) (poly-termlist q)) (max (poly-sugar p) (poly-sugar q))))
797
798(defun poly-uminus (ring p)
799 (make-poly-from-termlist (termlist-uminus ring (poly-termlist p)) (poly-sugar p)))
800
801(defun poly-mul (ring p q)
802 (make-poly-from-termlist (termlist-mul ring (poly-termlist p) (poly-termlist q)) (+ (poly-sugar p) (poly-sugar q))))
803
804(defun poly-expt (ring p n)
805 (make-poly-from-termlist (termlist-expt ring (poly-termlist p) n) (* n (poly-sugar p))))
806
807(defun poly-append (&rest plist)
808 (make-poly-from-termlist (apply #'append (mapcar #'poly-termlist plist))
809 (apply #'max (mapcar #'poly-sugar plist))))
810
811(defun poly-nreverse (p)
812 (setf (poly-termlist p) (nreverse (poly-termlist p)))
813 p)
814
815(defun poly-contract (p &optional (k 1))
816 (make-poly-from-termlist (termlist-contract (poly-termlist p) k)
817 (poly-sugar p)))
818
819(defun poly-extend (p &optional (m (make-monom 1 :initial-element 0)))
820 (make-poly-from-termlist
821 (termlist-extend (poly-termlist p) m)
822 (+ (poly-sugar p) (monom-sugar m))))
823
824(defun poly-add-variables (p k)
825 (setf (poly-termlist p) (termlist-add-variables (poly-termlist p) k))
826 p)
827
828(defun poly-list-add-variables (plist k)
829 (mapcar #'(lambda (p) (poly-add-variables p k)) plist))
830
831(defun poly-standard-extension (plist &aux (k (length plist)))
832 "Calculate [U1*P1,U2*P2,...,UK*PK], where PLIST=[P1,P2,...,PK]."
833 (declare (list plist) (fixnum k))
834 (labels ((incf-power (g i)
835 (dolist (x (poly-termlist g))
836 (incf (monom-elt (term-monom x) i)))
837 (incf (poly-sugar g))))
838 (setf plist (poly-list-add-variables plist k))
839 (dotimes (i k plist)
840 (incf-power (nth i plist) i))))
841
842(defun saturation-extension (ring f plist &aux (k (length plist)) (d (monom-dimension (poly-lm (car plist)))))
843 "Calculate [F, U1*P1-1,U2*P2-1,...,UK*PK-1], where PLIST=[P1,P2,...,PK]."
844 (setf f (poly-list-add-variables f k)
845 plist (mapcar #'(lambda (x)
846 (setf (poly-termlist x) (nconc (poly-termlist x)
847 (list (make-term (make-monom d :initial-element 0)
848 (funcall (ring-uminus ring) (funcall (ring-unit ring)))))))
849 x)
850 (poly-standard-extension plist)))
851 (append f plist))
852
853
854(defun polysaturation-extension (ring f plist &aux (k (length plist))
855 (d (+ k (length (poly-lm (car plist))))))
856 "Calculate [F, U1*P1+U2*P2+...+UK*PK-1], where PLIST=[P1,P2,...,PK]."
857 (setf f (poly-list-add-variables f k)
858 plist (apply #'poly-append (poly-standard-extension plist))
859 (cdr (last (poly-termlist plist))) (list (make-term (make-monom d :initial-element 0)
860 (funcall (ring-uminus ring) (funcall (ring-unit ring))))))
861 (append f (list plist)))
862
863(defun saturation-extension-1 (ring f p) (polysaturation-extension ring f (list p)))
864
865
866
867
868;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
869;;
870;; Evaluation of polynomial (prefix) expressions
871;;
872;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
873
874(defun coerce-coeff (ring expr vars)
875 "Coerce an element of the coefficient ring to a constant polynomial."
876 ;; Modular arithmetic handler by rat
877 (make-poly-from-termlist (list (make-term (make-monom (length vars) :initial-element 0)
878 (funcall (ring-parse ring) expr)))
879 0))
880
881(defun poly-eval (ring expr vars &optional (list-marker '[))
882 (labels ((p-eval (arg) (poly-eval ring arg vars))
883 (p-eval-list (args) (mapcar #'p-eval args))
884 (p-add (x y) (poly-add ring x y)))
885 (cond
886 ((eql expr 0) (make-poly-zero))
887 ((member expr vars :test #'equalp)
888 (let ((pos (position expr vars :test #'equalp)))
889 (make-variable ring (length vars) pos)))
890 ((atom expr)
891 (coerce-coeff ring expr vars))
892 ((eq (car expr) list-marker)
893 (cons list-marker (p-eval-list (cdr expr))))
894 (t
895 (case (car expr)
896 (+ (reduce #'p-add (p-eval-list (cdr expr))))
897 (- (case (length expr)
898 (1 (make-poly-zero))
899 (2 (poly-uminus ring (p-eval (cadr expr))))
900 (3 (poly-sub ring (p-eval (cadr expr)) (p-eval (caddr expr))))
901 (otherwise (poly-sub ring (p-eval (cadr expr))
902 (reduce #'p-add (p-eval-list (cddr expr)))))))
903 (*
904 (if (endp (cddr expr)) ;unary
905 (p-eval (cdr expr))
906 (reduce #'(lambda (p q) (poly-mul ring p q)) (p-eval-list (cdr expr)))))
907 (expt
908 (cond
909 ((member (cadr expr) vars :test #'equalp)
910 ;;Special handling of (expt var pow)
911 (let ((pos (position (cadr expr) vars :test #'equalp)))
912 (make-variable ring (length vars) pos (caddr expr))))
913 ((not (and (integerp (caddr expr)) (plusp (caddr expr))))
914 ;; Negative power means division in coefficient ring
915 ;; Non-integer power means non-polynomial coefficient
916 (coerce-coeff ring expr vars))
917 (t (poly-expt ring (p-eval (cadr expr)) (caddr expr)))))
918 (otherwise
919 (coerce-coeff ring expr vars)))))))
920
921
922
923
924;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
925;;
926;; Debugging/tracing
927;;
928;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
929
930
931
932(defmacro debug-cgb (&rest args)
933 `(when $poly_grobner_debug (format *terminal-io* ,@args)))
934
935;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
936;;
937;; An implementation of Grobner basis
938;;
939;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
940
941(defun spoly (ring f g)
942 "It yields the S-polynomial of polynomials F and G."
943 (declare (type poly f g))
944 (let* ((lcm (monom-lcm (poly-lm f) (poly-lm g)))
945 (mf (monom-div lcm (poly-lm f)))
946 (mg (monom-div lcm (poly-lm g))))
947 (declare (type monom mf mg))
948 (multiple-value-bind (c cf cg)
949 (funcall (ring-ezgcd ring) (poly-lc f) (poly-lc g))
950 (declare (ignore c))
951 (poly-sub
952 ring
953 (scalar-times-poly ring cg (monom-times-poly mf f))
954 (scalar-times-poly ring cf (monom-times-poly mg g))))))
955
956
957(defun poly-primitive-part (ring p)
958 "Divide polynomial P with integer coefficients by gcd of its
959coefficients and return the result."
960 (declare (type poly p))
961 (if (poly-zerop p)
962 (values p 1)
963 (let ((c (poly-content ring p)))
964 (values (make-poly-from-termlist (mapcar
965 #'(lambda (x)
966 (make-term (term-monom x)
967 (funcall (ring-div ring) (term-coeff x) c)))
968 (poly-termlist p))
969 (poly-sugar p))
970 c))))
971
972(defun poly-content (ring p)
973 "Greatest common divisor of the coefficients of the polynomial P. Use the RING structure
974to compute the greatest common divisor."
975 (declare (type poly p))
976 (reduce (ring-gcd ring) (mapcar #'term-coeff (rest (poly-termlist p))) :initial-value (poly-lc p)))
977
978
979
980;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
981;;
982;; An implementation of the division algorithm
983;;
984;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
985
986(defun grobner-op (ring c1 c2 m f g)
987 "Returns C2*F-C1*M*G, where F and G are polynomials M is a monomial.
988Assume that the leading terms will cancel."
989 #+grobner-check(funcall (ring-zerop ring)
990 (funcall (ring-sub ring)
991 (funcall (ring-mul ring) c2 (poly-lc f))
992 (funcall (ring-mul ring) c1 (poly-lc g))))
993 #+grobner-check(monom-equal-p (poly-lm f) (monom-mul m (poly-lm g)))
994 (poly-sub ring
995 (scalar-times-poly ring c2 f)
996 (scalar-times-poly ring c1 (monom-times-poly m g))))
997
998(defun poly-pseudo-divide (ring f fl)
999 "Pseudo-divide a polynomial F by the list of polynomials FL. Return
1000multiple values. The first value is a list of quotients A. The second
1001value is the remainder R. The third argument is a scalar coefficient
1002C, such that C*F can be divided by FL within the ring of coefficients,
1003which is not necessarily a field. Finally, the fourth value is an
1004integer count of the number of reductions performed. The resulting
1005objects satisfy the equation: C*F= sum A[i]*FL[i] + R."
1006 (declare (type poly f) (list fl))
1007 (do ((r (make-poly-zero))
1008 (c (funcall (ring-unit ring)))
1009 (a (make-list (length fl) :initial-element (make-poly-zero)))
1010 (division-count 0)
1011 (p f))
1012 ((poly-zerop p)
1013 (debug-cgb "~&~3T~d reduction~:p" division-count)
1014 (when (poly-zerop r) (debug-cgb " ---> 0"))
1015 (values (mapcar #'poly-nreverse a) (poly-nreverse r) c division-count))
1016 (declare (fixnum division-count))
1017 (do ((fl fl (rest fl)) ;scan list of divisors
1018 (b a (rest b)))
1019 ((cond
1020 ((endp fl) ;no division occurred
1021 (push (poly-lt p) (poly-termlist r)) ;move lt(p) to remainder
1022 (setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
1023 (pop (poly-termlist p)) ;remove lt(p) from p
1024 t)
1025 ((monom-divides-p (poly-lm (car fl)) (poly-lm p)) ;division occurred
1026 (incf division-count)
1027 (multiple-value-bind (gcd c1 c2)
1028 (funcall (ring-ezgcd ring) (poly-lc (car fl)) (poly-lc p))
1029 (declare (ignore gcd))
1030 (let ((m (monom-div (poly-lm p) (poly-lm (car fl)))))
1031 ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
1032 (mapl #'(lambda (x)
1033 (setf (car x) (scalar-times-poly ring c1 (car x))))
1034 a)
1035 (setf r (scalar-times-poly ring c1 r)
1036 c (funcall (ring-mul ring) c c1)
1037 p (grobner-op ring c2 c1 m p (car fl)))
1038 (push (make-term m c2) (poly-termlist (car b))))
1039 t)))))))
1040
1041(defun poly-exact-divide (ring f g)
1042 "Divide a polynomial F by another polynomial G. Assume that exact division
1043with no remainder is possible. Returns the quotient."
1044 (declare (type poly f g))
1045 (multiple-value-bind (quot rem coeff division-count)
1046 (poly-pseudo-divide ring f (list g))
1047 (declare (ignore division-count coeff)
1048 (list quot)
1049 (type poly rem)
1050 (type fixnum division-count))
1051 (unless (poly-zerop rem) (error "Exact division failed."))
1052 (car quot)))
1053
1054
1055
1056;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1057;;
1058;; An implementation of the normal form
1059;;
1060;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1061
1062(defun normal-form-step (ring fl p r c division-count
1063 &aux (g (find (poly-lm p) fl
1064 :test #'monom-divisible-by-p
1065 :key #'poly-lm)))
1066 (cond
1067 (g ;division possible
1068 (incf division-count)
1069 (multiple-value-bind (gcd cg cp)
1070 (funcall (ring-ezgcd ring) (poly-lc g) (poly-lc p))
1071 (declare (ignore gcd))
1072 (let ((m (monom-div (poly-lm p) (poly-lm g))))
1073 ;; Multiply the equation c*f=sum ai*fi+r+p by cg.
1074 (setf r (scalar-times-poly ring cg r)
1075 c (funcall (ring-mul ring) c cg)
1076 ;; p := cg*p-cp*m*g
1077 p (grobner-op ring cp cg m p g))))
1078 (debug-cgb "/"))
1079 (t ;no division possible
1080 (push (poly-lt p) (poly-termlist r)) ;move lt(p) to remainder
1081 (setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
1082 (pop (poly-termlist p)) ;remove lt(p) from p
1083 (debug-cgb "+")))
1084 (values p r c division-count))
1085
1086;; Merge it sometime with poly-pseudo-divide
1087(defun normal-form (ring f fl &optional (top-reduction-only $poly_top_reduction_only))
1088 ;; Loop invariant: c*f0=sum ai*fi+r+f, where f0 is the initial value of f
1089 #+grobner-check(when (null fl) (warn "normal-form: empty divisor list."))
1090 (do ((r (make-poly-zero))
1091 (c (funcall (ring-unit ring)))
1092 (division-count 0))
1093 ((or (poly-zerop f)
1094 ;;(endp fl)
1095 (and top-reduction-only (not (poly-zerop r))))
1096 (progn
1097 (debug-cgb "~&~3T~d reduction~:p" division-count)
1098 (when (poly-zerop r)
1099 (debug-cgb " ---> 0")))
1100 (setf (poly-termlist f) (nreconc (poly-termlist r) (poly-termlist f)))
1101 (values f c division-count))
1102 (declare (fixnum division-count)
1103 (type poly r))
1104 (multiple-value-setq (f r c division-count)
1105 (normal-form-step ring fl f r c division-count))))
1106
1107
1108
1109;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1110;;
1111;; These are provided mostly for debugging purposes To enable
1112;; verification of grobner bases with BUCHBERGER-CRITERION, do
1113;; (pushnew :grobner-check *features*) and compile/load this file.
1114;; With this feature, the calculations will slow down CONSIDERABLY.
1115;;
1116;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1117
1118(defun buchberger-criterion (ring g)
1119 "Returns T if G is a Grobner basis, by using the Buchberger
1120criterion: for every two polynomials h1 and h2 in G the S-polynomial
1121S(h1,h2) reduces to 0 modulo G."
1122 (every
1123 #'poly-zerop
1124 (makelist (normal-form ring (spoly ring (elt g i) (elt g j)) g nil)
1125 (i 0 (- (length g) 2))
1126 (j (1+ i) (1- (length g))))))
1127
1128(defun grobner-test (ring g f)
1129 "Test whether G is a Grobner basis and F is contained in G. Return T
1130upon success and NIL otherwise."
1131 (debug-cgb "~&GROBNER CHECK: ")
1132 (let (($poly_grobner_debug nil)
1133 (stat1 (buchberger-criterion ring g))
1134 (stat2
1135 (every #'poly-zerop
1136 (makelist (normal-form ring (copy-tree (elt f i)) g nil)
1137 (i 0 (1- (length f)))))))
1138 (unless stat1 (error "~&Buchberger criterion failed."))
1139 (unless stat2
1140 (error "~&Original polys not in ideal spanned by Grobner.")))
1141 (debug-cgb "~&GROBNER CHECK END")
1142 t)
1143
1144
1145
1146;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1147;;
1148;; Pair queue implementation
1149;;
1150;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1151
1152(defun sugar-pair-key (p q &aux (lcm (monom-lcm (poly-lm p) (poly-lm q)))
1153 (d (monom-sugar lcm)))
1154 "Returns list (S LCM-TOTAL-DEGREE) where S is the sugar of the S-polynomial of
1155polynomials P and Q, and LCM-TOTAL-DEGREE is the degree of is LCM(LM(P),LM(Q))."
1156 (declare (type poly p q) (type monom lcm) (type fixnum d))
1157 (cons (max
1158 (+ (- d (monom-sugar (poly-lm p))) (poly-sugar p))
1159 (+ (- d (monom-sugar (poly-lm q))) (poly-sugar q)))
1160 lcm))
1161
1162(defstruct (pair
1163 (:constructor make-pair (first second
1164 &aux
1165 (sugar (car (sugar-pair-key first second)))
1166 (division-data nil))))
1167 (first nil :type poly)
1168 (second nil :type poly)
1169 (sugar 0 :type fixnum)
1170 (division-data nil :type list))
1171
1172;;(defun pair-sugar (pair &aux (p (pair-first pair)) (q (pair-second pair)))
1173;; (car (sugar-pair-key p q)))
1174
1175(defun sugar-order (x y)
1176 "Pair order based on sugar, ties broken by normal strategy."
1177 (declare (type cons x y))
1178 (or (< (car x) (car y))
1179 (and (= (car x) (car y))
1180 (< (monom-total-degree (cdr x))
1181 (monom-total-degree (cdr y))))))
1182
1183(defvar *pair-key-function* #'sugar-pair-key
1184 "Function that, given two polynomials as argument, computed the key
1185in the pair queue.")
1186
1187(defvar *pair-order* #'sugar-order
1188 "Function that orders the keys of pairs.")
1189
1190(defun make-pair-queue ()
1191 "Constructs a priority queue for critical pairs."
1192 (make-priority-queue
1193 :element-type 'pair
1194 :element-key #'(lambda (pair) (funcall *pair-key-function* (pair-first pair) (pair-second pair)))
1195 :test *pair-order*))
1196
1197(defun pair-queue-initialize (pq f start
1198 &aux
1199 (s (1- (length f)))
1200 (b (nconc (makelist (make-pair (elt f i) (elt f j))
1201 (i 0 (1- start)) (j start s))
1202 (makelist (make-pair (elt f i) (elt f j))
1203 (i start (1- s)) (j (1+ i) s)))))
1204 "Initializes the priority for critical pairs. F is the initial list of polynomials.
1205START is the first position beyond the elements which form a partial
1206grobner basis, i.e. satisfy the Buchberger criterion."
1207 (declare (type priority-queue pq) (type fixnum start))
1208 (dolist (pair b pq)
1209 (priority-queue-insert pq pair)))
1210
1211(defun pair-queue-insert (b pair)
1212 (priority-queue-insert b pair))
1213
1214(defun pair-queue-remove (b)
1215 (priority-queue-remove b))
1216
1217(defun pair-queue-size (b)
1218 (priority-queue-size b))
1219
1220(defun pair-queue-empty-p (b)
1221 (priority-queue-empty-p b))
1222
1223;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1224;;
1225;; Buchberger Algorithm Implementation
1226;;
1227;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1228
1229(defun buchberger (ring f start &optional (top-reduction-only $poly_top_reduction_only))
1230 "An implementation of the Buchberger algorithm. Return Grobner basis
1231of the ideal generated by the polynomial list F. Polynomials 0 to
1232START-1 are assumed to be a Grobner basis already, so that certain
1233critical pairs will not be examined. If TOP-REDUCTION-ONLY set, top
1234reduction will be preformed. This function assumes that all polynomials
1235in F are non-zero."
1236 (declare (type fixnum start))
1237 (when (endp f) (return-from buchberger f)) ;cut startup costs
1238 (debug-cgb "~&GROBNER BASIS - BUCHBERGER ALGORITHM")
1239 (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
1240 #+grobner-check (when (plusp start)
1241 (grobner-test ring (subseq f 0 start) (subseq f 0 start)))
1242 ;;Initialize critical pairs
1243 (let ((b (pair-queue-initialize (make-pair-queue)
1244 f start))
1245 (b-done (make-hash-table :test #'equal)))
1246 (declare (type priority-queue b) (type hash-table b-done))
1247 (dotimes (i (1- start))
1248 (do ((j (1+ i) (1+ j))) ((>= j start))
1249 (setf (gethash (list (elt f i) (elt f j)) b-done) t)))
1250 (do ()
1251 ((pair-queue-empty-p b)
1252 #+grobner-check(grobner-test ring f f)
1253 (debug-cgb "~&GROBNER END")
1254 f)
1255 (let ((pair (pair-queue-remove b)))
1256 (declare (type pair pair))
1257 (cond
1258 ((criterion-1 pair) nil)
1259 ((criterion-2 pair b-done f) nil)
1260 (t
1261 (let ((sp (normal-form ring (spoly ring (pair-first pair)
1262 (pair-second pair))
1263 f top-reduction-only)))
1264 (declare (type poly sp))
1265 (cond
1266 ((poly-zerop sp)
1267 nil)
1268 (t
1269 (setf sp (poly-primitive-part ring sp)
1270 f (nconc f (list sp)))
1271 ;; Add new critical pairs
1272 (dolist (h f)
1273 (pair-queue-insert b (make-pair h sp)))
1274 (debug-cgb "~&Sugar: ~d Polynomials: ~d; Pairs left: ~d; Pairs done: ~d;"
1275 (pair-sugar pair) (length f) (pair-queue-size b)
1276 (hash-table-count b-done)))))))
1277 (setf (gethash (list (pair-first pair) (pair-second pair)) b-done)
1278 t)))))
1279
1280(defun parallel-buchberger (ring f start &optional (top-reduction-only $poly_top_reduction_only))
1281 "An implementation of the Buchberger algorithm. Return Grobner basis
1282of the ideal generated by the polynomial list F. Polynomials 0 to
1283START-1 are assumed to be a Grobner basis already, so that certain
1284critical pairs will not be examined. If TOP-REDUCTION-ONLY set, top
1285reduction will be preformed."
1286 (declare (ignore top-reduction-only)
1287 (type fixnum start))
1288 (when (endp f) (return-from parallel-buchberger f)) ;cut startup costs
1289 (debug-cgb "~&GROBNER BASIS - PARALLEL-BUCHBERGER ALGORITHM")
1290 (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
1291 #+grobner-check (when (plusp start)
1292 (grobner-test ring (subseq f 0 start) (subseq f 0 start)))
1293 ;;Initialize critical pairs
1294 (let ((b (pair-queue-initialize (make-pair-queue) f start))
1295 (b-done (make-hash-table :test #'equal)))
1296 (declare (type priority-queue b)
1297 (type hash-table b-done))
1298 (dotimes (i (1- start))
1299 (do ((j (1+ i) (1+ j))) ((>= j start))
1300 (declare (type fixnum j))
1301 (setf (gethash (list (elt f i) (elt f j)) b-done) t)))
1302 (do ()
1303 ((pair-queue-empty-p b)
1304 #+grobner-check(grobner-test ring f f)
1305 (debug-cgb "~&GROBNER END")
1306 f)
1307 (let ((pair (pair-queue-remove b)))
1308 (when (null (pair-division-data pair))
1309 (setf (pair-division-data pair) (list (spoly ring
1310 (pair-first pair)
1311 (pair-second pair))
1312 (make-poly-zero)
1313 (funcall (ring-unit ring))
1314 0)))
1315 (cond
1316 ((criterion-1 pair) nil)
1317 ((criterion-2 pair b-done f) nil)
1318 (t
1319 (let* ((dd (pair-division-data pair))
1320 (p (first dd))
1321 (sp (second dd))
1322 (c (third dd))
1323 (division-count (fourth dd)))
1324 (cond
1325 ((poly-zerop p) ;normal form completed
1326 (debug-cgb "~&~3T~d reduction~:p" division-count)
1327 (cond
1328 ((poly-zerop sp)
1329 (debug-cgb " ---> 0")
1330 nil)
1331 (t
1332 (setf sp (poly-nreverse sp)
1333 sp (poly-primitive-part ring sp)
1334 f (nconc f (list sp)))
1335 ;; Add new critical pairs
1336 (dolist (h f)
1337 (pair-queue-insert b (make-pair h sp)))
1338 (debug-cgb "~&Sugar: ~d Polynomials: ~d; Pairs left: ~d; Pairs done: ~d;"
1339 (pair-sugar pair) (length f) (pair-queue-size b)
1340 (hash-table-count b-done))))
1341 (setf (gethash (list (pair-first pair) (pair-second pair))
1342 b-done) t))
1343 (t ;normal form not complete
1344 (do ()
1345 ((cond
1346 ((> (poly-sugar sp) (pair-sugar pair))
1347 (debug-cgb "(~a)?" (poly-sugar sp))
1348 t)
1349 ((poly-zerop p)
1350 (debug-cgb ".")
1351 t)
1352 (t nil))
1353 (setf (first dd) p
1354 (second dd) sp
1355 (third dd) c
1356 (fourth dd) division-count
1357 (pair-sugar pair) (poly-sugar sp))
1358 (pair-queue-insert b pair))
1359 (multiple-value-setq (p sp c division-count)
1360 (normal-form-step ring f p sp c division-count))))))))))))
1361
1362;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1363;;
1364;; Grobner Criteria
1365;;
1366;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1367
1368(defun criterion-1 (pair)
1369 "Returns T if the leading monomials of the two polynomials
1370in G pointed to by the integers in PAIR have disjoint (relatively prime)
1371monomials. This test is known as the first Buchberger criterion."
1372 (declare (type pair pair))
1373 (let ((f (pair-first pair))
1374 (g (pair-second pair)))
1375 (when (monom-rel-prime-p (poly-lm f) (poly-lm g))
1376 (debug-cgb ":1")
1377 (return-from criterion-1 t))))
1378
1379(defun criterion-2 (pair b-done partial-basis
1380 &aux (f (pair-first pair)) (g (pair-second pair))
1381 (place :before))
1382 "Returns T if the leading monomial of some element P of
1383PARTIAL-BASIS divides the LCM of the leading monomials of the two
1384polynomials in the polynomial list PARTIAL-BASIS, and P paired with
1385each of the polynomials pointed to by the the PAIR has already been
1386treated, as indicated by the absence in the hash table B-done."
1387 (declare (type pair pair) (type hash-table b-done)
1388 (type poly f g))
1389 ;; In the code below we assume that pairs are ordered as follows:
1390 ;; if PAIR is (I J) then I appears before J in the PARTIAL-BASIS.
1391 ;; We traverse the list PARTIAL-BASIS and keep track of where we
1392 ;; are, so that we can produce the pairs in the correct order
1393 ;; when we check whether they have been processed, i.e they
1394 ;; appear in the hash table B-done
1395 (dolist (h partial-basis nil)
1396 (cond
1397 ((eq h f)
1398 #+grobner-check(assert (eq place :before))
1399 (setf place :in-the-middle))
1400 ((eq h g)
1401 #+grobner-check(assert (eq place :in-the-middle))
1402 (setf place :after))
1403 ((and (monom-divides-monom-lcm-p (poly-lm h) (poly-lm f) (poly-lm g))
1404 (gethash (case place
1405 (:before (list h f))
1406 ((:in-the-middle :after) (list f h)))
1407 b-done)
1408 (gethash (case place
1409 ((:before :in-the-middle) (list h g))
1410 (:after (list g h)))
1411 b-done))
1412 (debug-cgb ":2")
1413 (return-from criterion-2 t)))))
1414
1415
1416
1417;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1418;;
1419;; An implementation of the algorithm of Gebauer and Moeller, as
1420;; described in the book of Becker-Weispfenning, p. 232
1421;;
1422;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1423
1424(defun gebauer-moeller (ring f start &optional (top-reduction-only $poly_top_reduction_only))
1425 "Compute Grobner basis by using the algorithm of Gebauer and
1426Moeller. This algorithm is described as BUCHBERGERNEW2 in the book by
1427Becker-Weispfenning entitled ``Grobner Bases''. This function assumes
1428that all polynomials in F are non-zero."
1429 (declare (ignore top-reduction-only)
1430 (type fixnum start))
1431 (cond
1432 ((endp f) (return-from gebauer-moeller nil))
1433 ((endp (cdr f))
1434 (return-from gebauer-moeller (list (poly-primitive-part ring (car f))))))
1435 (debug-cgb "~&GROBNER BASIS - GEBAUER MOELLER ALGORITHM")
1436 (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
1437 #+grobner-check (when (plusp start)
1438 (grobner-test ring (subseq f 0 start) (subseq f 0 start)))
1439 (let ((b (make-pair-queue))
1440 (g (subseq f 0 start))
1441 (f1 (subseq f start)))
1442 (do () ((endp f1))
1443 (multiple-value-setq (g b)
1444 (gebauer-moeller-update g b (poly-primitive-part ring (pop f1)))))
1445 (do () ((pair-queue-empty-p b))
1446 (let* ((pair (pair-queue-remove b))
1447 (g1 (pair-first pair))
1448 (g2 (pair-second pair))
1449 (h (normal-form ring (spoly ring g1 g2)
1450 g
1451 nil #| Always fully reduce! |#
1452 )))
1453 (unless (poly-zerop h)
1454 (setf h (poly-primitive-part ring h))
1455 (multiple-value-setq (g b)
1456 (gebauer-moeller-update g b h))
1457 (debug-cgb "~&Sugar: ~d Polynomials: ~d; Pairs left: ~d~%"
1458 (pair-sugar pair) (length g) (pair-queue-size b))
1459 )))
1460 #+grobner-check(grobner-test ring g f)
1461 (debug-cgb "~&GROBNER END")
1462 g))
1463
1464(defun gebauer-moeller-update (g b h
1465 &aux
1466 c d e
1467 (b-new (make-pair-queue))
1468 g-new)
1469 "An implementation of the auxillary UPDATE algorithm used by the
1470Gebauer-Moeller algorithm. G is a list of polynomials, B is a list of
1471critical pairs and H is a new polynomial which possibly will be added
1472to G. The naming conventions used are very close to the one used in
1473the book of Becker-Weispfenning."
1474 (declare
1475 #+allegro (dynamic-extent b)
1476 (type poly h)
1477 (type priority-queue b))
1478 (setf c g d nil)
1479 (do () ((endp c))
1480 (let ((g1 (pop c)))
1481 (declare (type poly g1))
1482 (when (or (monom-rel-prime-p (poly-lm h) (poly-lm g1))
1483 (and
1484 (notany #'(lambda (g2) (monom-lcm-divides-monom-lcm-p
1485 (poly-lm h) (poly-lm g2)
1486 (poly-lm h) (poly-lm g1)))
1487 c)
1488 (notany #'(lambda (g2) (monom-lcm-divides-monom-lcm-p
1489 (poly-lm h) (poly-lm g2)
1490 (poly-lm h) (poly-lm g1)))
1491 d)))
1492 (push g1 d))))
1493 (setf e nil)
1494 (do () ((endp d))
1495 (let ((g1 (pop d)))
1496 (declare (type poly g1))
1497 (unless (monom-rel-prime-p (poly-lm h) (poly-lm g1))
1498 (push g1 e))))
1499 (do () ((pair-queue-empty-p b))
1500 (let* ((pair (pair-queue-remove b))
1501 (g1 (pair-first pair))
1502 (g2 (pair-second pair)))
1503 (declare (type pair pair)
1504 (type poly g1 g2))
1505 (when (or (not (monom-divides-monom-lcm-p
1506 (poly-lm h)
1507 (poly-lm g1) (poly-lm g2)))
1508 (monom-lcm-equal-monom-lcm-p
1509 (poly-lm g1) (poly-lm h)
1510 (poly-lm g1) (poly-lm g2))
1511 (monom-lcm-equal-monom-lcm-p
1512 (poly-lm h) (poly-lm g2)
1513 (poly-lm g1) (poly-lm g2)))
1514 (pair-queue-insert b-new (make-pair g1 g2)))))
1515 (dolist (g3 e)
1516 (pair-queue-insert b-new (make-pair h g3)))
1517 (setf g-new nil)
1518 (do () ((endp g))
1519 (let ((g1 (pop g)))
1520 (declare (type poly g1))
1521 (unless (monom-divides-p (poly-lm h) (poly-lm g1))
1522 (push g1 g-new))))
1523 (push h g-new)
1524 (values g-new b-new))
1525
1526
1527
1528;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1529;;
1530;; Standard postprocessing of Grobner bases
1531;;
1532;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1533
1534(defun reduction (ring plist)
1535 "Reduce a list of polynomials PLIST, so that non of the terms in any of
1536the polynomials is divisible by a leading monomial of another
1537polynomial. Return the reduced list."
1538 (do ((q plist)
1539 (found t))
1540 ((not found)
1541 (mapcar #'(lambda (x) (poly-primitive-part ring x)) q))
1542 ;;Find p in Q such that p is reducible mod Q\{p}
1543 (setf found nil)
1544 (dolist (x q)
1545 (let ((q1 (remove x q)))
1546 (multiple-value-bind (h c div-count)
1547 (normal-form ring x q1 nil #| not a top reduction! |# )
1548 (declare (ignore c))
1549 (unless (zerop div-count)
1550 (setf found t q q1)
1551 (unless (poly-zerop h)
1552 (setf q (nconc q1 (list h))))
1553 (return)))))))
1554
1555(defun minimization (p)
1556 "Returns a sublist of the polynomial list P spanning the same
1557monomial ideal as P but minimal, i.e. no leading monomial
1558of a polynomial in the sublist divides the leading monomial
1559of another polynomial."
1560 (do ((q p)
1561 (found t))
1562 ((not found) q)
1563 ;;Find p in Q such that lm(p) is in LM(Q\{p})
1564 (setf found nil
1565 q (dolist (x q q)
1566 (let ((q1 (remove x q)))
1567 (when (member-if #'(lambda (p) (monom-divides-p (poly-lm x) (poly-lm p))) q1)
1568 (setf found t)
1569 (return q1)))))))
1570
1571(defun poly-normalize (ring p &aux (c (poly-lc p)))
1572 "Divide a polynomial by its leading coefficient. It assumes
1573that the division is possible, which may not always be the
1574case in rings which are not fields. The exact division operator
1575is assumed to be provided by the RING structure of the
1576COEFFICIENT-RING package."
1577 (mapc #'(lambda (term)
1578 (setf (term-coeff term) (funcall (ring-div ring) (term-coeff term) c)))
1579 (poly-termlist p))
1580 p)
1581
1582(defun poly-normalize-list (ring plist)
1583 "Divide every polynomial in a list PLIST by its leading coefficient. "
1584 (mapcar #'(lambda (x) (poly-normalize ring x)) plist))
1585
1586
1587
1588;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1589;;
1590;; Algorithm and Pair heuristic selection
1591;;
1592;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1593
1594(defun find-grobner-function (algorithm)
1595 "Return a function which calculates Grobner basis, based on its
1596names. Names currently used are either Lisp symbols, Maxima symbols or
1597keywords."
1598 (ecase algorithm
1599 ((buchberger :buchberger $buchberger) #'buchberger)
1600 ((parallel-buchberger :parallel-buchberger $parallel_buchberger) #'parallel-buchberger)
1601 ((gebauer-moeller :gebauer_moeller $gebauer_moeller) #'gebauer-moeller)))
1602
1603(defun grobner (ring f &optional (start 0) (top-reduction-only nil))
1604 ;;(setf F (sort F #'< :key #'sugar))
1605 (funcall
1606 (find-grobner-function $poly_grobner_algorithm)
1607 ring f start top-reduction-only))
1608
1609(defun reduced-grobner (ring f &optional (start 0) (top-reduction-only $poly_top_reduction_only))
1610 (reduction ring (grobner ring f start top-reduction-only)))
1611
1612
1613
1614;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1615;;
1616;; Operations in ideal theory
1617;;
1618;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1619
1620;; Does the term depend on variable K?
1621(defun term-depends-p (term k)
1622 "Return T if the term TERM depends on variable number K."
1623 (monom-depends-p (term-monom term) k))
1624
1625;; Does the polynomial P depend on variable K?
1626(defun poly-depends-p (p k)
1627 "Return T if the term polynomial P depends on variable number K."
1628 (some #'(lambda (term) (term-depends-p term k)) (poly-termlist p)))
1629
1630(defun ring-intersection (plist k)
1631 "This function assumes that polynomial list PLIST is a Grobner basis
1632and it calculates the intersection with the ring R[x[k+1],...,x[n]], i.e.
1633it discards polynomials which depend on variables x[0], x[1], ..., x[k]."
1634 (dotimes (i k plist)
1635 (setf plist
1636 (remove-if #'(lambda (p)
1637 (poly-depends-p p i))
1638 plist))))
1639
1640(defun elimination-ideal (ring flist k
1641 &optional (top-reduction-only $poly_top_reduction_only) (start 0)
1642 &aux (*monomial-order*
1643 (or *elimination-order*
1644 (elimination-order k))))
1645 (ring-intersection (reduced-grobner ring flist start top-reduction-only) k))
1646
1647(defun colon-ideal (ring f g &optional (top-reduction-only $poly_top_reduction_only))
1648 "Returns the reduced Grobner basis of the colon ideal Id(F):Id(G),
1649where F and G are two lists of polynomials. The colon ideal I:J is
1650defined as the set of polynomials H such that for all polynomials W in
1651J the polynomial W*H belongs to I."
1652 (cond
1653 ((endp g)
1654 ;;Id(G) consists of 0 only so W*0=0 belongs to Id(F)
1655 (if (every #'poly-zerop f)
1656 (error "First ideal must be non-zero.")
1657 (list (make-poly
1658 (list (make-term
1659 (make-monom (monom-dimension (poly-lm (find-if-not #'poly-zerop f)))
1660 :initial-element 0)
1661 (funcall (ring-unit ring))))))))
1662 ((endp (cdr g))
1663 (colon-ideal-1 ring f (car g) top-reduction-only))
1664 (t
1665 (ideal-intersection ring
1666 (colon-ideal-1 ring f (car g) top-reduction-only)
1667 (colon-ideal ring f (rest g) top-reduction-only)
1668 top-reduction-only))))
1669
1670(defun colon-ideal-1 (ring f g &optional (top-reduction-only $poly_top_reduction_only))
1671 "Returns the reduced Grobner basis of the colon ideal Id(F):Id({G}), where
1672F is a list of polynomials and G is a polynomial."
1673 (mapcar #'(lambda (x) (poly-exact-divide ring x g)) (ideal-intersection ring f (list g) top-reduction-only)))
1674
1675
1676(defun ideal-intersection (ring f g &optional (top-reduction-only $poly_top_reduction_only)
1677 &aux (*monomial-order* (or *elimination-order*
1678 #'elimination-order-1)))
1679 (mapcar #'poly-contract
1680 (ring-intersection
1681 (reduced-grobner
1682 ring
1683 (append (mapcar #'(lambda (p) (poly-extend p (make-monom 1 :initial-element 1))) f)
1684 (mapcar #'(lambda (p)
1685 (poly-append (poly-extend (poly-uminus ring p)
1686 (make-monom 1 :initial-element 1))
1687 (poly-extend p)))
1688 g))
1689 0
1690 top-reduction-only)
1691 1)))
1692
1693(defun poly-lcm (ring f g)
1694 "Return LCM (least common multiple) of two polynomials F and G.
1695The polynomials must be ordered according to monomial order PRED
1696and their coefficients must be compatible with the RING structure
1697defined in the COEFFICIENT-RING package."
1698 (cond
1699 ((poly-zerop f) f)
1700 ((poly-zerop g) g)
1701 ((and (endp (cdr (poly-termlist f))) (endp (cdr (poly-termlist g))))
1702 (let ((m (monom-lcm (poly-lm f) (poly-lm g))))
1703 (make-poly-from-termlist (list (make-term m (funcall (ring-lcm ring) (poly-lc f) (poly-lc g)))))))
1704 (t
1705 (multiple-value-bind (f f-cont)
1706 (poly-primitive-part ring f)
1707 (multiple-value-bind (g g-cont)
1708 (poly-primitive-part ring g)
1709 (scalar-times-poly
1710 ring
1711 (funcall (ring-lcm ring) f-cont g-cont)
1712 (poly-primitive-part ring (car (ideal-intersection ring (list f) (list g) nil)))))))))
1713
1714;; Do two Grobner bases yield the same ideal?
1715(defun grobner-equal (ring g1 g2)
1716 "Returns T if two lists of polynomials G1 and G2, assumed to be Grobner bases,
1717generate the same ideal, and NIL otherwise."
1718 (and (grobner-subsetp ring g1 g2) (grobner-subsetp ring g2 g1)))
1719
1720(defun grobner-subsetp (ring g1 g2)
1721 "Returns T if a list of polynomials G1 generates
1722an ideal contained in the ideal generated by a polynomial list G2,
1723both G1 and G2 assumed to be Grobner bases. Returns NIL otherwise."
1724 (every #'(lambda (p) (grobner-member ring p g2)) g1))
1725
1726(defun grobner-member (ring p g)
1727 "Returns T if a polynomial P belongs to the ideal generated by the
1728polynomial list G, which is assumed to be a Grobner basis. Returns NIL otherwise."
1729 (poly-zerop (normal-form ring p g nil)))
1730
1731;; Calculate F : p^inf
1732(defun ideal-saturation-1 (ring f p start &optional (top-reduction-only $poly_top_reduction_only)
1733 &aux (*monomial-order* (or *elimination-order*
1734 #'elimination-order-1)))
1735 "Returns the reduced Grobner basis of the saturation of the ideal
1736generated by a polynomial list F in the ideal generated by a single
1737polynomial P. The saturation ideal is defined as the set of
1738polynomials H such for some natural number n (* (EXPT P N) H) is in the ideal
1739F. Geometrically, over an algebraically closed field, this is the set
1740of polynomials in the ideal generated by F which do not identically
1741vanish on the variety of P."
1742 (mapcar
1743 #'poly-contract
1744 (ring-intersection
1745 (reduced-grobner
1746 ring
1747 (saturation-extension-1 ring f p)
1748 start top-reduction-only)
1749 1)))
1750
1751
1752
1753;; Calculate F : p1^inf : p2^inf : ... : ps^inf
1754(defun ideal-polysaturation-1 (ring f plist start &optional (top-reduction-only $poly_top_reduction_only))
1755 "Returns the reduced Grobner basis of the ideal obtained by a
1756sequence of successive saturations in the polynomials
1757of the polynomial list PLIST of the ideal generated by the
1758polynomial list F."
1759 (cond
1760 ((endp plist) (reduced-grobner ring f start top-reduction-only))
1761 (t (let ((g (ideal-saturation-1 ring f (car plist) start top-reduction-only)))
1762 (ideal-polysaturation-1 ring g (rest plist) (length g) top-reduction-only)))))
1763
1764(defun ideal-saturation (ring f g start &optional (top-reduction-only $poly_top_reduction_only)
1765 &aux
1766 (k (length g))
1767 (*monomial-order* (or *elimination-order*
1768 (elimination-order k))))
1769 "Returns the reduced Grobner basis of the saturation of the ideal
1770generated by a polynomial list F in the ideal generated a polynomial
1771list G. The saturation ideal is defined as the set of polynomials H
1772such for some natural number n and some P in the ideal generated by G
1773the polynomial P**N * H is in the ideal spanned by F. Geometrically,
1774over an algebraically closed field, this is the set of polynomials in
1775the ideal generated by F which do not identically vanish on the
1776variety of G."
1777 (mapcar
1778 #'(lambda (q) (poly-contract q k))
1779 (ring-intersection
1780 (reduced-grobner ring
1781 (polysaturation-extension ring f g)
1782 start
1783 top-reduction-only)
1784 k)))
1785
1786(defun ideal-polysaturation (ring f ideal-list start &optional (top-reduction-only $poly_top_reduction_only))
1787 "Returns the reduced Grobner basis of the ideal obtained by a
1788successive applications of IDEAL-SATURATION to F and lists of
1789polynomials in the list IDEAL-LIST."
1790 (cond
1791 ((endp ideal-list) f)
1792 (t (let ((h (ideal-saturation ring f (car ideal-list) start top-reduction-only)))
1793 (ideal-polysaturation ring h (rest ideal-list) (length h) top-reduction-only)))))
1794
1795
1796
1797;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1798;;
1799;; Set up the coefficients to be polynomials
1800;;
1801;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1802
1803;; (defun poly-ring (ring vars)
1804;; (make-ring
1805;; :parse #'(lambda (expr) (poly-eval ring expr vars))
1806;; :unit #'(lambda () (poly-unit ring (length vars)))
1807;; :zerop #'poly-zerop
1808;; :add #'(lambda (x y) (poly-add ring x y))
1809;; :sub #'(lambda (x y) (poly-sub ring x y))
1810;; :uminus #'(lambda (x) (poly-uminus ring x))
1811;; :mul #'(lambda (x y) (poly-mul ring x y))
1812;; :div #'(lambda (x y) (poly-exact-divide ring x y))
1813;; :lcm #'(lambda (x y) (poly-lcm ring x y))
1814;; :ezgcd #'(lambda (x y &aux (gcd (poly-gcd ring x y)))
1815;; (values gcd
1816;; (poly-exact-divide ring x gcd)
1817;; (poly-exact-divide ring y gcd)))
1818;; :gcd #'(lambda (x y) (poly-gcd x y))))
1819
1820
1821
1822;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1823;;
1824;; Conversion from internal to infix form
1825;;
1826;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1827
1828(defun coerce-to-infix (poly-type object vars)
1829 (case poly-type
1830 (:termlist
1831 `(+ ,@(mapcar #'(lambda (term) (coerce-to-infix :term term vars)) object)))
1832 (:polynomial
1833 (coerce-to-infix :termlist (poly-termlist object) vars))
1834 (:poly-list
1835 `([ ,@(mapcar #'(lambda (p) (coerce-to-infix :polynomial p vars)) object)))
1836 (:term
1837 `(* ,(term-coeff object)
1838 ,@(mapcar #'(lambda (var power) `(expt ,var ,power))
1839 vars (monom-exponents (term-monom object)))))
1840 (otherwise
1841 object)))
1842
1843
1844
1845;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1846;;
1847;; Maxima expression ring
1848;;
1849;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1850
1851(defparameter *expression-ring*
1852 (make-ring
1853 ;;(defun coeff-zerop (expr) (meval1 `(($is) (($equal) ,expr 0))))
1854 :parse #'(lambda (expr)
1855 (when modulus (setf expr ($rat expr)))
1856 expr)
1857 :unit #'(lambda () (if modulus ($rat 1) 1))
1858 :zerop #'(lambda (expr)
1859 ;;When is exactly a maxima expression equal to 0?
1860 (cond ((numberp expr)
1861 (= expr 0))
1862 ((atom expr) nil)
1863 (t
1864 (case (caar expr)
1865 (mrat (eql ($ratdisrep expr) 0))
1866 (otherwise (eql ($totaldisrep expr) 0))))))
1867 :add #'(lambda (x y) (m+ x y))
1868 :sub #'(lambda (x y) (m- x y))
1869 :uminus #'(lambda (x) (m- x))
1870 :mul #'(lambda (x y) (m* x y))
1871 ;;(defun coeff-div (x y) (cadr ($divide x y)))
1872 :div #'(lambda (x y) (m// x y))
1873 :lcm #'(lambda (x y) (meval1 `((|$LCM|) ,x ,y)))
1874 :ezgcd #'(lambda (x y) (apply #'values (cdr ($ezgcd ($totaldisrep x) ($totaldisrep y)))))
1875 ;; :gcd #'(lambda (x y) (second ($ezgcd x y)))))
1876 :gcd #'(lambda (x y) ($gcd x y))))
1877
1878(defvar *maxima-ring* *expression-ring*
1879 "The ring of coefficients, over which all polynomials
1880are assumed to be defined.")
1881
1882
1883
1884;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1885;;
1886;; Maxima expression parsing
1887;;
1888;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1889
1890(defun equal-test-p (expr1 expr2)
1891 (alike1 expr1 expr2))
1892
1893(defun coerce-maxima-list (expr)
1894 "convert a maxima list to lisp list."
1895 (cond
1896 ((and (consp (car expr)) (eql (caar expr) 'mlist)) (cdr expr))
1897 (t expr)))
1898
1899(defun free-of-vars (expr vars) (apply #'$freeof `(,@vars ,expr)))
1900
1901(defun parse-poly (expr vars &aux (vars (coerce-maxima-list vars)))
1902 "Convert a maxima polynomial expression EXPR in variables VARS to internal form."
1903 (labels ((parse (arg) (parse-poly arg vars))
1904 (parse-list (args) (mapcar #'parse args)))
1905 (cond
1906 ((eql expr 0) (make-poly-zero))
1907 ((member expr vars :test #'equal-test-p)
1908 (let ((pos (position expr vars :test #'equal-test-p)))
1909 (make-variable *maxima-ring* (length vars) pos)))
1910 ((free-of-vars expr vars)
1911 ;;This means that variable-free CRE and Poisson forms will be converted
1912 ;;to coefficients intact
1913 (coerce-coeff *maxima-ring* expr vars))
1914 (t
1915 (case (caar expr)
1916 (mplus (reduce #'(lambda (x y) (poly-add *maxima-ring* x y)) (parse-list (cdr expr))))
1917 (mminus (poly-uminus *maxima-ring* (parse (cadr expr))))
1918 (mtimes
1919 (if (endp (cddr expr)) ;unary
1920 (parse (cdr expr))
1921 (reduce #'(lambda (p q) (poly-mul *maxima-ring* p q)) (parse-list (cdr expr)))))
1922 (mexpt
1923 (cond
1924 ((member (cadr expr) vars :test #'equal-test-p)
1925 ;;Special handling of (expt var pow)
1926 (let ((pos (position (cadr expr) vars :test #'equal-test-p)))
1927 (make-variable *maxima-ring* (length vars) pos (caddr expr))))
1928 ((not (and (integerp (caddr expr)) (plusp (caddr expr))))
1929 ;; Negative power means division in coefficient ring
1930 ;; Non-integer power means non-polynomial coefficient
1931 (mtell "~%Warning: Expression ~%~M~%contains power which is not a positive integer. Parsing as coefficient.~%"
1932 expr)
1933 (coerce-coeff *maxima-ring* expr vars))
1934 (t (poly-expt *maxima-ring* (parse (cadr expr)) (caddr expr)))))
1935 (mrat (parse ($ratdisrep expr)))
1936 (mpois (parse ($outofpois expr)))
1937 (otherwise
1938 (coerce-coeff *maxima-ring* expr vars)))))))
1939
1940(defun parse-poly-list (expr vars)
1941 (case (caar expr)
1942 (mlist (mapcar #'(lambda (p) (parse-poly p vars)) (cdr expr)))
1943 (t (merror "Expression ~M is not a list of polynomials in variables ~M."
1944 expr vars))))
1945(defun parse-poly-list-list (poly-list-list vars)
1946 (mapcar #'(lambda (g) (parse-poly-list g vars)) (coerce-maxima-list poly-list-list)))
1947
1948
1949
1950;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1951;;
1952;; Order utilities
1953;;
1954;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1955(defun find-order (order)
1956 "This function returns the order function bases on its name."
1957 (cond
1958 ((null order) nil)
1959 ((symbolp order)
1960 (case order
1961 ((lex :lex $lex) #'lex>)
1962 ((grlex :grlex $grlex) #'grlex>)
1963 ((grevlex :grevlex $grevlex) #'grevlex>)
1964 ((invlex :invlex $invlex) #'invlex>)
1965 ((elimination-order-1 :elimination-order-1 elimination_order_1) #'elimination-order-1)
1966 (otherwise
1967 (mtell "~%Warning: Order ~M not found. Using default.~%" order))))
1968 (t
1969 (mtell "~%Order specification ~M is not recognized. Using default.~%" order)
1970 nil)))
1971
1972(defun find-ring (ring)
1973 "This function returns the ring structure bases on input symbol."
1974 (cond
1975 ((null ring) nil)
1976 ((symbolp ring)
1977 (case ring
1978 ((expression-ring :expression-ring $expression_ring) *expression-ring*)
1979 ((ring-of-integers :ring-of-integers $ring_of_integers) *ring-of-integers*)
1980 (otherwise
1981 (mtell "~%Warning: Ring ~M not found. Using default.~%" ring))))
1982 (t
1983 (mtell "~%Ring specification ~M is not recognized. Using default.~%" ring)
1984 nil)))
1985
1986(defmacro with-monomial-order ((order) &body body)
1987 "Evaluate BODY with monomial order set to ORDER."
1988 `(let ((*monomial-order* (or (find-order ,order) *monomial-order*)))
1989 . ,body))
1990
1991(defmacro with-coefficient-ring ((ring) &body body)
1992 "Evaluate BODY with coefficient ring set to RING."
1993 `(let ((*maxima-ring* (or (find-ring ,ring) *maxima-ring*)))
1994 . ,body))
1995
1996(defmacro with-elimination-orders ((primary secondary elimination-order)
1997 &body body)
1998 "Evaluate BODY with primary and secondary elimination orders set to PRIMARY and SECONDARY."
1999 `(let ((*primary-elimination-order* (or (find-order ,primary) *primary-elimination-order*))
2000 (*secondary-elimination-order* (or (find-order ,secondary) *secondary-elimination-order*))
2001 (*elimination-order* (or (find-order ,elimination-order) *elimination-order*)))
2002 . ,body))
2003
2004
2005
2006;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2007;;
2008;; Conversion from internal form to Maxima general form
2009;;
2010;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2011
2012(defun maxima-head ()
2013 (if $poly_return_term_list
2014 '(mlist)
2015 '(mplus)))
2016
2017(defun coerce-to-maxima (poly-type object vars)
2018 (case poly-type
2019 (:polynomial
2020 `(,(maxima-head) ,@(mapcar #'(lambda (term) (coerce-to-maxima :term term vars)) (poly-termlist object))))
2021 (:poly-list
2022 `((mlist) ,@(mapcar #'(lambda (p) ($ratdisrep (coerce-to-maxima :polynomial p vars))) object)))
2023 (:term
2024 `((mtimes) ,($ratdisrep (term-coeff object))
2025 ,@(mapcar #'(lambda (var power) `((mexpt) ,var ,power))
2026 vars (monom-exponents (term-monom object)))))
2027 ;; Assumes that Lisp and Maxima logicals coincide
2028 (:logical object)
2029 (otherwise
2030 object)))
2031
2032
2033
2034;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2035;;
2036;; Macro facility for writing Maxima-level wrappers for
2037;; functions operating on internal representation
2038;;
2039;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2040
2041(defmacro with-parsed-polynomials (((maxima-vars &optional (maxima-new-vars nil new-vars-supplied-p))
2042 &key (polynomials nil)
2043 (poly-lists nil)
2044 (poly-list-lists nil)
2045 (value-type nil))
2046 &body body
2047 &aux (vars (gensym))
2048 (new-vars (gensym)))
2049 `(let ((,vars (coerce-maxima-list ,maxima-vars))
2050 ,@(when new-vars-supplied-p
2051 (list `(,new-vars (coerce-maxima-list ,maxima-new-vars)))))
2052 (coerce-to-maxima
2053 ,value-type
2054 (with-coefficient-ring ($poly_coefficient_ring)
2055 (with-monomial-order ($poly_monomial_order)
2056 (with-elimination-orders ($poly_primary_elimination_order
2057 $poly_secondary_elimination_order
2058 $poly_elimination_order)
2059 (let ,(let ((args nil))
2060 (dolist (p polynomials args)
2061 (setf args (cons `(,p (parse-poly ,p ,vars)) args)))
2062 (dolist (p poly-lists args)
2063 (setf args (cons `(,p (parse-poly-list ,p ,vars)) args)))
2064 (dolist (p poly-list-lists args)
2065 (setf args (cons `(,p (parse-poly-list-list ,p ,vars)) args))))
2066 . ,body))))
2067 ,(if new-vars-supplied-p
2068 `(append ,vars ,new-vars)
2069 vars))))
2070
2071(defmacro define-unop (maxima-name fun-name
2072 &optional (documentation nil documentation-supplied-p))
2073 "Define a MAXIMA-level unary operator MAXIMA-NAME corresponding to unary function FUN-NAME."
2074 `(defun ,maxima-name (p vars
2075 &aux
2076 (vars (coerce-maxima-list vars))
2077 (p (parse-poly p vars)))
2078 ,@(when documentation-supplied-p (list documentation))
2079 (coerce-to-maxima :polynomial (,fun-name *maxima-ring* p) vars)))
2080
2081(defmacro define-binop (maxima-name fun-name
2082 &optional (documentation nil documentation-supplied-p))
2083 "Define a MAXIMA-level binary operator MAXIMA-NAME corresponding to binary function FUN-NAME."
2084 `(defmfun ,maxima-name (p q vars
2085 &aux
2086 (vars (coerce-maxima-list vars))
2087 (p (parse-poly p vars))
2088 (q (parse-poly q vars)))
2089 ,@(when documentation-supplied-p (list documentation))
2090 (coerce-to-maxima :polynomial (,fun-name *maxima-ring* p q) vars)))
2091
2092
2093
2094;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2095;;
2096;; Maxima-level interface functions
2097;;
2098;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2099
2100;; Auxillary function for removing zero polynomial
2101(defun remzero (plist) (remove #'poly-zerop plist))
2102
2103;;Simple operators
2104
2105(define-binop $poly_add poly-add
2106 "Adds two polynomials P and Q")
2107
2108(define-binop $poly_subtract poly-sub
2109 "Subtracts a polynomial Q from P.")
2110
2111(define-binop $poly_multiply poly-mul
2112 "Returns the product of polynomials P and Q.")
2113
2114(define-binop $poly_s_polynomial spoly
2115 "Returns the syzygy polynomial (S-polynomial) of two polynomials P and Q.")
2116
2117(define-unop $poly_primitive_part poly-primitive-part
2118 "Returns the polynomial P divided by GCD of its coefficients.")
2119
2120(define-unop $poly_normalize poly-normalize
2121 "Returns the polynomial P divided by the leading coefficient.")
2122
2123;;Functions
2124
2125(defmfun $poly_expand (p vars)
2126 "This function is equivalent to EXPAND(P) if P parses correctly to a polynomial.
2127If the representation is not compatible with a polynomial in variables VARS,
2128the result is an error."
2129 (with-parsed-polynomials ((vars) :polynomials (p)
2130 :value-type :polynomial)
2131 p))
2132
2133(defmfun $poly_expt (p n vars)
2134 (with-parsed-polynomials ((vars) :polynomials (p) :value-type :polynomial)
2135 (poly-expt *maxima-ring* p n)))
2136
2137(defmfun $poly_content (p vars)
2138 (with-parsed-polynomials ((vars) :polynomials (p))
2139 (poly-content *maxima-ring* p)))
2140
2141(defmfun $poly_pseudo_divide (f fl vars
2142 &aux (vars (coerce-maxima-list vars))
2143 (f (parse-poly f vars))
2144 (fl (parse-poly-list fl vars)))
2145 (multiple-value-bind (quot rem c division-count)
2146 (poly-pseudo-divide *maxima-ring* f fl)
2147 `((mlist)
2148 ,(coerce-to-maxima :poly-list quot vars)
2149 ,(coerce-to-maxima :polynomial rem vars)
2150 ,c
2151 ,division-count)))
2152
2153(defmfun $poly_exact_divide (f g vars)
2154 (with-parsed-polynomials ((vars) :polynomials (f g) :value-type :polynomial)
2155 (poly-exact-divide *maxima-ring* f g)))
2156
2157(defmfun $poly_normal_form (f fl vars)
2158 (with-parsed-polynomials ((vars) :polynomials (f)
2159 :poly-lists (fl)
2160 :value-type :polynomial)
2161 (normal-form *maxima-ring* f (remzero fl) nil)))
2162
2163(defmfun $poly_buchberger_criterion (g vars)
2164 (with-parsed-polynomials ((vars) :poly-lists (g) :value-type :logical)
2165 (buchberger-criterion *maxima-ring* g)))
2166
2167(defmfun $poly_buchberger (fl vars)
2168 (with-parsed-polynomials ((vars) :poly-lists (fl) :value-type :poly-list)
2169 (buchberger *maxima-ring* (remzero fl) 0 nil)))
2170
2171(defmfun $poly_reduction (plist vars)
2172 (with-parsed-polynomials ((vars) :poly-lists (plist)
2173 :value-type :poly-list)
2174 (reduction *maxima-ring* plist)))
2175
2176(defmfun $poly_minimization (plist vars)
2177 (with-parsed-polynomials ((vars) :poly-lists (plist)
2178 :value-type :poly-list)
2179 (minimization plist)))
2180
2181(defmfun $poly_normalize_list (plist vars)
2182 (with-parsed-polynomials ((vars) :poly-lists (plist)
2183 :value-type :poly-list)
2184 (poly-normalize-list *maxima-ring* plist)))
2185
2186(defmfun $poly_grobner (f vars)
2187 (with-parsed-polynomials ((vars) :poly-lists (f)
2188 :value-type :poly-list)
2189 (grobner *maxima-ring* (remzero f))))
2190
2191(defmfun $poly_reduced_grobner (f vars)
2192 (with-parsed-polynomials ((vars) :poly-lists (f)
2193 :value-type :poly-list)
2194 (reduced-grobner *maxima-ring* (remzero f))))
2195
2196(defmfun $poly_depends_p (p var mvars
2197 &aux (vars (coerce-maxima-list mvars))
2198 (pos (position var vars)))
2199 (if (null pos)
2200 (merror "~%Variable ~M not in the list of variables ~M." var mvars)
2201 (poly-depends-p (parse-poly p vars) pos)))
2202
2203(defmfun $poly_elimination_ideal (flist k vars)
2204 (with-parsed-polynomials ((vars) :poly-lists (flist)
2205 :value-type :poly-list)
2206 (elimination-ideal *maxima-ring* flist k nil 0)))
2207
2208(defmfun $poly_colon_ideal (f g vars)
2209 (with-parsed-polynomials ((vars) :poly-lists (f g) :value-type :poly-list)
2210 (colon-ideal *maxima-ring* f g nil)))
2211
2212(defmfun $poly_ideal_intersection (f g vars)
2213 (with-parsed-polynomials ((vars) :poly-lists (f g) :value-type :poly-list)
2214 (ideal-intersection *maxima-ring* f g nil)))
2215
2216(defmfun $poly_lcm (f g vars)
2217 (with-parsed-polynomials ((vars) :polynomials (f g) :value-type :polynomial)
2218 (poly-lcm *maxima-ring* f g)))
2219
2220(defmfun $poly_gcd (f g vars)
2221 ($first ($divide (m* f g) ($poly_lcm f g vars))))
2222
2223(defmfun $poly_grobner_equal (g1 g2 vars)
2224 (with-parsed-polynomials ((vars) :poly-lists (g1 g2))
2225 (grobner-equal *maxima-ring* g1 g2)))
2226
2227(defmfun $poly_grobner_subsetp (g1 g2 vars)
2228 (with-parsed-polynomials ((vars) :poly-lists (g1 g2))
2229 (grobner-subsetp *maxima-ring* g1 g2)))
2230
2231(defmfun $poly_grobner_member (p g vars)
2232 (with-parsed-polynomials ((vars) :polynomials (p) :poly-lists (g))
2233 (grobner-member *maxima-ring* p g)))
2234
2235(defmfun $poly_ideal_saturation1 (f p vars)
2236 (with-parsed-polynomials ((vars) :poly-lists (f) :polynomials (p)
2237 :value-type :poly-list)
2238 (ideal-saturation-1 *maxima-ring* f p 0)))
2239
2240(defmfun $poly_saturation_extension (f plist vars new-vars)
2241 (with-parsed-polynomials ((vars new-vars)
2242 :poly-lists (f plist)
2243 :value-type :poly-list)
2244 (saturation-extension *maxima-ring* f plist)))
2245
2246(defmfun $poly_polysaturation_extension (f plist vars new-vars)
2247 (with-parsed-polynomials ((vars new-vars)
2248 :poly-lists (f plist)
2249 :value-type :poly-list)
2250 (polysaturation-extension *maxima-ring* f plist)))
2251
2252(defmfun $poly_ideal_polysaturation1 (f plist vars)
2253 (with-parsed-polynomials ((vars) :poly-lists (f plist)
2254 :value-type :poly-list)
2255 (ideal-polysaturation-1 *maxima-ring* f plist 0 nil)))
2256
2257(defmfun $poly_ideal_saturation (f g vars)
2258 (with-parsed-polynomials ((vars) :poly-lists (f g)
2259 :value-type :poly-list)
2260 (ideal-saturation *maxima-ring* f g 0 nil)))
2261
2262(defmfun $poly_ideal_polysaturation (f ideal-list vars)
2263 (with-parsed-polynomials ((vars) :poly-lists (f)
2264 :poly-list-lists (ideal-list)
2265 :value-type :poly-list)
2266 (ideal-polysaturation *maxima-ring* f ideal-list 0 nil)))
2267
2268(defmfun $poly_lt (f vars)
2269 (with-parsed-polynomials ((vars) :polynomials (f) :value-type :polynomial)
2270 (make-poly-from-termlist (list (poly-lt f)))))
2271
2272(defmfun $poly_lm (f vars)
2273 (with-parsed-polynomials ((vars) :polynomials (f) :value-type :polynomial)
2274 (make-poly-from-termlist (list (make-term (poly-lm f) (funcall (ring-unit *maxima-ring*)))))))
2275
Note: See TracBrowser for help on using the repository browser.