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

Last change on this file since 294 was 29, checked in by Marek Rychlik, 12 years ago

* empty log message *

File size: 80.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(defun set-pair-heuristic (method)
1613 "Sets up variables *PAIR-KEY-FUNCTION* and *PAIR-ORDER* used
1614to determine the priority of critical pairs in the priority queue."
1615 (ecase method
1616 ((sugar :sugar $sugar)
1617 (setf *pair-key-function* #'sugar-pair-key
1618 *pair-order* #'sugar-order))
1619; ((minimal-mock-spoly :minimal-mock-spoly $minimal_mock_spoly)
1620; (setf *pair-key-function* #'mock-spoly
1621; *pair-order* #'mock-spoly-order))
1622 ((minimal-lcm :minimal-lcm $minimal_lcm)
1623 (setf *pair-key-function* #'(lambda (p q)
1624 (monom-lcm (poly-lm p) (poly-lm q)))
1625 *pair-order* #'reverse-monomial-order))
1626 ((minimal-total-degree :minimal-total-degree $minimal_total_degree)
1627 (setf *pair-key-function* #'(lambda (p q)
1628 (monom-total-degree
1629 (monom-lcm (poly-lm p) (poly-lm q))))
1630 *pair-order* #'<))
1631 ((minimal-length :minimal-length $minimal_length)
1632 (setf *pair-key-function* #'(lambda (p q)
1633 (+ (poly-length p) (poly-length q)))
1634 *pair-order* #'<))))
1635
1636
1637
1638;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1639;;
1640;; Operations in ideal theory
1641;;
1642;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1643
1644;; Does the term depend on variable K?
1645(defun term-depends-p (term k)
1646 "Return T if the term TERM depends on variable number K."
1647 (monom-depends-p (term-monom term) k))
1648
1649;; Does the polynomial P depend on variable K?
1650(defun poly-depends-p (p k)
1651 "Return T if the term polynomial P depends on variable number K."
1652 (some #'(lambda (term) (term-depends-p term k)) (poly-termlist p)))
1653
1654(defun ring-intersection (plist k)
1655 "This function assumes that polynomial list PLIST is a Grobner basis
1656and it calculates the intersection with the ring R[x[k+1],...,x[n]], i.e.
1657it discards polynomials which depend on variables x[0], x[1], ..., x[k]."
1658 (dotimes (i k plist)
1659 (setf plist
1660 (remove-if #'(lambda (p)
1661 (poly-depends-p p i))
1662 plist))))
1663
1664(defun elimination-ideal (ring flist k
1665 &optional (top-reduction-only $poly_top_reduction_only) (start 0)
1666 &aux (*monomial-order*
1667 (or *elimination-order*
1668 (elimination-order k))))
1669 (ring-intersection (reduced-grobner ring flist start top-reduction-only) k))
1670
1671(defun colon-ideal (ring f g &optional (top-reduction-only $poly_top_reduction_only))
1672 "Returns the reduced Grobner basis of the colon ideal Id(F):Id(G),
1673where F and G are two lists of polynomials. The colon ideal I:J is
1674defined as the set of polynomials H such that for all polynomials W in
1675J the polynomial W*H belongs to I."
1676 (cond
1677 ((endp g)
1678 ;;Id(G) consists of 0 only so W*0=0 belongs to Id(F)
1679 (if (every #'poly-zerop f)
1680 (error "First ideal must be non-zero.")
1681 (list (make-poly
1682 (list (make-term
1683 (make-monom (monom-dimension (poly-lm (find-if-not #'poly-zerop f)))
1684 :initial-element 0)
1685 (funcall (ring-unit ring))))))))
1686 ((endp (cdr g))
1687 (colon-ideal-1 ring f (car g) top-reduction-only))
1688 (t
1689 (ideal-intersection ring
1690 (colon-ideal-1 ring f (car g) top-reduction-only)
1691 (colon-ideal ring f (rest g) top-reduction-only)
1692 top-reduction-only))))
1693
1694(defun colon-ideal-1 (ring f g &optional (top-reduction-only $poly_top_reduction_only))
1695 "Returns the reduced Grobner basis of the colon ideal Id(F):Id({G}), where
1696F is a list of polynomials and G is a polynomial."
1697 (mapcar #'(lambda (x) (poly-exact-divide ring x g)) (ideal-intersection ring f (list g) top-reduction-only)))
1698
1699
1700(defun ideal-intersection (ring f g &optional (top-reduction-only $poly_top_reduction_only)
1701 &aux (*monomial-order* (or *elimination-order*
1702 #'elimination-order-1)))
1703 (mapcar #'poly-contract
1704 (ring-intersection
1705 (reduced-grobner
1706 ring
1707 (append (mapcar #'(lambda (p) (poly-extend p (make-monom 1 :initial-element 1))) f)
1708 (mapcar #'(lambda (p)
1709 (poly-append (poly-extend (poly-uminus ring p)
1710 (make-monom 1 :initial-element 1))
1711 (poly-extend p)))
1712 g))
1713 0
1714 top-reduction-only)
1715 1)))
1716
1717(defun poly-lcm (ring f g)
1718 "Return LCM (least common multiple) of two polynomials F and G.
1719The polynomials must be ordered according to monomial order PRED
1720and their coefficients must be compatible with the RING structure
1721defined in the COEFFICIENT-RING package."
1722 (cond
1723 ((poly-zerop f) f)
1724 ((poly-zerop g) g)
1725 ((and (endp (cdr (poly-termlist f))) (endp (cdr (poly-termlist g))))
1726 (let ((m (monom-lcm (poly-lm f) (poly-lm g))))
1727 (make-poly-from-termlist (list (make-term m (funcall (ring-lcm ring) (poly-lc f) (poly-lc g)))))))
1728 (t
1729 (multiple-value-bind (f f-cont)
1730 (poly-primitive-part ring f)
1731 (multiple-value-bind (g g-cont)
1732 (poly-primitive-part ring g)
1733 (scalar-times-poly
1734 ring
1735 (funcall (ring-lcm ring) f-cont g-cont)
1736 (poly-primitive-part ring (car (ideal-intersection ring (list f) (list g) nil)))))))))
1737
1738;; Do two Grobner bases yield the same ideal?
1739(defun grobner-equal (ring g1 g2)
1740 "Returns T if two lists of polynomials G1 and G2, assumed to be Grobner bases,
1741generate the same ideal, and NIL otherwise."
1742 (and (grobner-subsetp ring g1 g2) (grobner-subsetp ring g2 g1)))
1743
1744(defun grobner-subsetp (ring g1 g2)
1745 "Returns T if a list of polynomials G1 generates
1746an ideal contained in the ideal generated by a polynomial list G2,
1747both G1 and G2 assumed to be Grobner bases. Returns NIL otherwise."
1748 (every #'(lambda (p) (grobner-member ring p g2)) g1))
1749
1750(defun grobner-member (ring p g)
1751 "Returns T if a polynomial P belongs to the ideal generated by the
1752polynomial list G, which is assumed to be a Grobner basis. Returns NIL otherwise."
1753 (poly-zerop (normal-form ring p g nil)))
1754
1755;; Calculate F : p^inf
1756(defun ideal-saturation-1 (ring f p start &optional (top-reduction-only $poly_top_reduction_only)
1757 &aux (*monomial-order* (or *elimination-order*
1758 #'elimination-order-1)))
1759 "Returns the reduced Grobner basis of the saturation of the ideal
1760generated by a polynomial list F in the ideal generated by a single
1761polynomial P. The saturation ideal is defined as the set of
1762polynomials H such for some natural number n (* (EXPT P N) H) is in the ideal
1763F. Geometrically, over an algebraically closed field, this is the set
1764of polynomials in the ideal generated by F which do not identically
1765vanish on the variety of P."
1766 (mapcar
1767 #'poly-contract
1768 (ring-intersection
1769 (reduced-grobner
1770 ring
1771 (saturation-extension-1 ring f p)
1772 start top-reduction-only)
1773 1)))
1774
1775
1776
1777;; Calculate F : p1^inf : p2^inf : ... : ps^inf
1778(defun ideal-polysaturation-1 (ring f plist start &optional (top-reduction-only $poly_top_reduction_only))
1779 "Returns the reduced Grobner basis of the ideal obtained by a
1780sequence of successive saturations in the polynomials
1781of the polynomial list PLIST of the ideal generated by the
1782polynomial list F."
1783 (cond
1784 ((endp plist) (reduced-grobner ring f start top-reduction-only))
1785 (t (let ((g (ideal-saturation-1 ring f (car plist) start top-reduction-only)))
1786 (ideal-polysaturation-1 ring g (rest plist) (length g) top-reduction-only)))))
1787
1788(defun ideal-saturation (ring f g start &optional (top-reduction-only $poly_top_reduction_only)
1789 &aux
1790 (k (length g))
1791 (*monomial-order* (or *elimination-order*
1792 (elimination-order k))))
1793 "Returns the reduced Grobner basis of the saturation of the ideal
1794generated by a polynomial list F in the ideal generated a polynomial
1795list G. The saturation ideal is defined as the set of polynomials H
1796such for some natural number n and some P in the ideal generated by G
1797the polynomial P**N * H is in the ideal spanned by F. Geometrically,
1798over an algebraically closed field, this is the set of polynomials in
1799the ideal generated by F which do not identically vanish on the
1800variety of G."
1801 (mapcar
1802 #'(lambda (q) (poly-contract q k))
1803 (ring-intersection
1804 (reduced-grobner ring
1805 (polysaturation-extension ring f g)
1806 start
1807 top-reduction-only)
1808 k)))
1809
1810(defun ideal-polysaturation (ring f ideal-list start &optional (top-reduction-only $poly_top_reduction_only))
1811 "Returns the reduced Grobner basis of the ideal obtained by a
1812successive applications of IDEAL-SATURATION to F and lists of
1813polynomials in the list IDEAL-LIST."
1814 (cond
1815 ((endp ideal-list) f)
1816 (t (let ((h (ideal-saturation ring f (car ideal-list) start top-reduction-only)))
1817 (ideal-polysaturation ring h (rest ideal-list) (length h) top-reduction-only)))))
1818
1819
1820
1821;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1822;;
1823;; Set up the coefficients to be polynomials
1824;;
1825;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1826
1827;; (defun poly-ring (ring vars)
1828;; (make-ring
1829;; :parse #'(lambda (expr) (poly-eval ring expr vars))
1830;; :unit #'(lambda () (poly-unit ring (length vars)))
1831;; :zerop #'poly-zerop
1832;; :add #'(lambda (x y) (poly-add ring x y))
1833;; :sub #'(lambda (x y) (poly-sub ring x y))
1834;; :uminus #'(lambda (x) (poly-uminus ring x))
1835;; :mul #'(lambda (x y) (poly-mul ring x y))
1836;; :div #'(lambda (x y) (poly-exact-divide ring x y))
1837;; :lcm #'(lambda (x y) (poly-lcm ring x y))
1838;; :ezgcd #'(lambda (x y &aux (gcd (poly-gcd ring x y)))
1839;; (values gcd
1840;; (poly-exact-divide ring x gcd)
1841;; (poly-exact-divide ring y gcd)))
1842;; :gcd #'(lambda (x y) (poly-gcd x y))))
1843
1844
1845
1846;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1847;;
1848;; Conversion from internal to infix form
1849;;
1850;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1851
1852(defun coerce-to-infix (poly-type object vars)
1853 (case poly-type
1854 (:termlist
1855 `(+ ,@(mapcar #'(lambda (term) (coerce-to-infix :term term vars)) object)))
1856 (:polynomial
1857 (coerce-to-infix :termlist (poly-termlist object) vars))
1858 (:poly-list
1859 `([ ,@(mapcar #'(lambda (p) (coerce-to-infix :polynomial p vars)) object)))
1860 (:term
1861 `(* ,(term-coeff object)
1862 ,@(mapcar #'(lambda (var power) `(expt ,var ,power))
1863 vars (monom-exponents (term-monom object)))))
1864 (otherwise
1865 object)))
1866
1867
1868
1869;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1870;;
1871;; Maxima expression ring
1872;;
1873;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1874
1875(defparameter *expression-ring*
1876 (make-ring
1877 ;;(defun coeff-zerop (expr) (meval1 `(($is) (($equal) ,expr 0))))
1878 :parse #'(lambda (expr)
1879 (when modulus (setf expr ($rat expr)))
1880 expr)
1881 :unit #'(lambda () (if modulus ($rat 1) 1))
1882 :zerop #'(lambda (expr)
1883 ;;When is exactly a maxima expression equal to 0?
1884 (cond ((numberp expr)
1885 (= expr 0))
1886 ((atom expr) nil)
1887 (t
1888 (case (caar expr)
1889 (mrat (eql ($ratdisrep expr) 0))
1890 (otherwise (eql ($totaldisrep expr) 0))))))
1891 :add #'(lambda (x y) (m+ x y))
1892 :sub #'(lambda (x y) (m- x y))
1893 :uminus #'(lambda (x) (m- x))
1894 :mul #'(lambda (x y) (m* x y))
1895 ;;(defun coeff-div (x y) (cadr ($divide x y)))
1896 :div #'(lambda (x y) (m// x y))
1897 :lcm #'(lambda (x y) (meval1 `((|$LCM|) ,x ,y)))
1898 :ezgcd #'(lambda (x y) (apply #'values (cdr ($ezgcd ($totaldisrep x) ($totaldisrep y)))))
1899 ;; :gcd #'(lambda (x y) (second ($ezgcd x y)))))
1900 :gcd #'(lambda (x y) ($gcd x y))))
1901
1902(defvar *maxima-ring* *expression-ring*
1903 "The ring of coefficients, over which all polynomials
1904are assumed to be defined.")
1905
1906
1907
1908;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1909;;
1910;; Maxima expression parsing
1911;;
1912;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1913
1914(defun equal-test-p (expr1 expr2)
1915 (alike1 expr1 expr2))
1916
1917(defun coerce-maxima-list (expr)
1918 "convert a maxima list to lisp list."
1919 (cond
1920 ((and (consp (car expr)) (eql (caar expr) 'mlist)) (cdr expr))
1921 (t expr)))
1922
1923(defun free-of-vars (expr vars) (apply #'$freeof `(,@vars ,expr)))
1924
1925(defun parse-poly (expr vars &aux (vars (coerce-maxima-list vars)))
1926 "Convert a maxima polynomial expression EXPR in variables VARS to internal form."
1927 (labels ((parse (arg) (parse-poly arg vars))
1928 (parse-list (args) (mapcar #'parse args)))
1929 (cond
1930 ((eql expr 0) (make-poly-zero))
1931 ((member expr vars :test #'equal-test-p)
1932 (let ((pos (position expr vars :test #'equal-test-p)))
1933 (make-variable *maxima-ring* (length vars) pos)))
1934 ((free-of-vars expr vars)
1935 ;;This means that variable-free CRE and Poisson forms will be converted
1936 ;;to coefficients intact
1937 (coerce-coeff *maxima-ring* expr vars))
1938 (t
1939 (case (caar expr)
1940 (mplus (reduce #'(lambda (x y) (poly-add *maxima-ring* x y)) (parse-list (cdr expr))))
1941 (mminus (poly-uminus *maxima-ring* (parse (cadr expr))))
1942 (mtimes
1943 (if (endp (cddr expr)) ;unary
1944 (parse (cdr expr))
1945 (reduce #'(lambda (p q) (poly-mul *maxima-ring* p q)) (parse-list (cdr expr)))))
1946 (mexpt
1947 (cond
1948 ((member (cadr expr) vars :test #'equal-test-p)
1949 ;;Special handling of (expt var pow)
1950 (let ((pos (position (cadr expr) vars :test #'equal-test-p)))
1951 (make-variable *maxima-ring* (length vars) pos (caddr expr))))
1952 ((not (and (integerp (caddr expr)) (plusp (caddr expr))))
1953 ;; Negative power means division in coefficient ring
1954 ;; Non-integer power means non-polynomial coefficient
1955 (mtell "~%Warning: Expression ~%~M~%contains power which is not a positive integer. Parsing as coefficient.~%"
1956 expr)
1957 (coerce-coeff *maxima-ring* expr vars))
1958 (t (poly-expt *maxima-ring* (parse (cadr expr)) (caddr expr)))))
1959 (mrat (parse ($ratdisrep expr)))
1960 (mpois (parse ($outofpois expr)))
1961 (otherwise
1962 (coerce-coeff *maxima-ring* expr vars)))))))
1963
1964(defun parse-poly-list (expr vars)
1965 (case (caar expr)
1966 (mlist (mapcar #'(lambda (p) (parse-poly p vars)) (cdr expr)))
1967 (t (merror "Expression ~M is not a list of polynomials in variables ~M."
1968 expr vars))))
1969(defun parse-poly-list-list (poly-list-list vars)
1970 (mapcar #'(lambda (g) (parse-poly-list g vars)) (coerce-maxima-list poly-list-list)))
1971
1972
1973
1974;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1975;;
1976;; Order utilities
1977;;
1978;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1979(defun find-order (order)
1980 "This function returns the order function bases on its name."
1981 (cond
1982 ((null order) nil)
1983 ((symbolp order)
1984 (case order
1985 ((lex :lex $lex) #'lex>)
1986 ((grlex :grlex $grlex) #'grlex>)
1987 ((grevlex :grevlex $grevlex) #'grevlex>)
1988 ((invlex :invlex $invlex) #'invlex>)
1989 ((elimination-order-1 :elimination-order-1 elimination_order_1) #'elimination-order-1)
1990 (otherwise
1991 (mtell "~%Warning: Order ~M not found. Using default.~%" order))))
1992 (t
1993 (mtell "~%Order specification ~M is not recognized. Using default.~%" order)
1994 nil)))
1995
1996(defun find-ring (ring)
1997 "This function returns the ring structure bases on input symbol."
1998 (cond
1999 ((null ring) nil)
2000 ((symbolp ring)
2001 (case ring
2002 ((expression-ring :expression-ring $expression_ring) *expression-ring*)
2003 ((ring-of-integers :ring-of-integers $ring_of_integers) *ring-of-integers*)
2004 (otherwise
2005 (mtell "~%Warning: Ring ~M not found. Using default.~%" ring))))
2006 (t
2007 (mtell "~%Ring specification ~M is not recognized. Using default.~%" ring)
2008 nil)))
2009
2010(defmacro with-monomial-order ((order) &body body)
2011 "Evaluate BODY with monomial order set to ORDER."
2012 `(let ((*monomial-order* (or (find-order ,order) *monomial-order*)))
2013 . ,body))
2014
2015(defmacro with-coefficient-ring ((ring) &body body)
2016 "Evaluate BODY with coefficient ring set to RING."
2017 `(let ((*maxima-ring* (or (find-ring ,ring) *maxima-ring*)))
2018 . ,body))
2019
2020(defmacro with-elimination-orders ((primary secondary elimination-order)
2021 &body body)
2022 "Evaluate BODY with primary and secondary elimination orders set to PRIMARY and SECONDARY."
2023 `(let ((*primary-elimination-order* (or (find-order ,primary) *primary-elimination-order*))
2024 (*secondary-elimination-order* (or (find-order ,secondary) *secondary-elimination-order*))
2025 (*elimination-order* (or (find-order ,elimination-order) *elimination-order*)))
2026 . ,body))
2027
2028
2029
2030;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2031;;
2032;; Conversion from internal form to Maxima general form
2033;;
2034;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2035
2036(defun maxima-head ()
2037 (if $poly_return_term_list
2038 '(mlist)
2039 '(mplus)))
2040
2041(defun coerce-to-maxima (poly-type object vars)
2042 (case poly-type
2043 (:polynomial
2044 `(,(maxima-head) ,@(mapcar #'(lambda (term) (coerce-to-maxima :term term vars)) (poly-termlist object))))
2045 (:poly-list
2046 `((mlist) ,@(mapcar #'(lambda (p) ($ratdisrep (coerce-to-maxima :polynomial p vars))) object)))
2047 (:term
2048 `((mtimes) ,($ratdisrep (term-coeff object))
2049 ,@(mapcar #'(lambda (var power) `((mexpt) ,var ,power))
2050 vars (monom-exponents (term-monom object)))))
2051 ;; Assumes that Lisp and Maxima logicals coincide
2052 (:logical object)
2053 (otherwise
2054 object)))
2055
2056
2057
2058;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2059;;
2060;; Macro facility for writing Maxima-level wrappers for
2061;; functions operating on internal representation
2062;;
2063;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2064
2065(defmacro with-parsed-polynomials (((maxima-vars &optional (maxima-new-vars nil new-vars-supplied-p))
2066 &key (polynomials nil)
2067 (poly-lists nil)
2068 (poly-list-lists nil)
2069 (value-type nil))
2070 &body body
2071 &aux (vars (gensym))
2072 (new-vars (gensym)))
2073 `(let ((,vars (coerce-maxima-list ,maxima-vars))
2074 ,@(when new-vars-supplied-p
2075 (list `(,new-vars (coerce-maxima-list ,maxima-new-vars)))))
2076 (coerce-to-maxima
2077 ,value-type
2078 (with-coefficient-ring ($poly_coefficient_ring)
2079 (with-monomial-order ($poly_monomial_order)
2080 (with-elimination-orders ($poly_primary_elimination_order
2081 $poly_secondary_elimination_order
2082 $poly_elimination_order)
2083 (let ,(let ((args nil))
2084 (dolist (p polynomials args)
2085 (setf args (cons `(,p (parse-poly ,p ,vars)) args)))
2086 (dolist (p poly-lists args)
2087 (setf args (cons `(,p (parse-poly-list ,p ,vars)) args)))
2088 (dolist (p poly-list-lists args)
2089 (setf args (cons `(,p (parse-poly-list-list ,p ,vars)) args))))
2090 . ,body))))
2091 ,(if new-vars-supplied-p
2092 `(append ,vars ,new-vars)
2093 vars))))
2094
2095(defmacro define-unop (maxima-name fun-name
2096 &optional (documentation nil documentation-supplied-p))
2097 "Define a MAXIMA-level unary operator MAXIMA-NAME corresponding to unary function FUN-NAME."
2098 `(defun ,maxima-name (p vars
2099 &aux
2100 (vars (coerce-maxima-list vars))
2101 (p (parse-poly p vars)))
2102 ,@(when documentation-supplied-p (list documentation))
2103 (coerce-to-maxima :polynomial (,fun-name *maxima-ring* p) vars)))
2104
2105(defmacro define-binop (maxima-name fun-name
2106 &optional (documentation nil documentation-supplied-p))
2107 "Define a MAXIMA-level binary operator MAXIMA-NAME corresponding to binary function FUN-NAME."
2108 `(defmfun ,maxima-name (p q vars
2109 &aux
2110 (vars (coerce-maxima-list vars))
2111 (p (parse-poly p vars))
2112 (q (parse-poly q vars)))
2113 ,@(when documentation-supplied-p (list documentation))
2114 (coerce-to-maxima :polynomial (,fun-name *maxima-ring* p q) vars)))
2115
2116
2117
2118;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2119;;
2120;; Maxima-level interface functions
2121;;
2122;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2123
2124;; Auxillary function for removing zero polynomial
2125(defun remzero (plist) (remove #'poly-zerop plist))
2126
2127;;Simple operators
2128
2129(define-binop $poly_add poly-add
2130 "Adds two polynomials P and Q")
2131
2132(define-binop $poly_subtract poly-sub
2133 "Subtracts a polynomial Q from P.")
2134
2135(define-binop $poly_multiply poly-mul
2136 "Returns the product of polynomials P and Q.")
2137
2138(define-binop $poly_s_polynomial spoly
2139 "Returns the syzygy polynomial (S-polynomial) of two polynomials P and Q.")
2140
2141(define-unop $poly_primitive_part poly-primitive-part
2142 "Returns the polynomial P divided by GCD of its coefficients.")
2143
2144(define-unop $poly_normalize poly-normalize
2145 "Returns the polynomial P divided by the leading coefficient.")
2146
2147;;Functions
2148
2149(defmfun $poly_expand (p vars)
2150 "This function is equivalent to EXPAND(P) if P parses correctly to a polynomial.
2151If the representation is not compatible with a polynomial in variables VARS,
2152the result is an error."
2153 (with-parsed-polynomials ((vars) :polynomials (p)
2154 :value-type :polynomial)
2155 p))
2156
2157(defmfun $poly_expt (p n vars)
2158 (with-parsed-polynomials ((vars) :polynomials (p) :value-type :polynomial)
2159 (poly-expt *maxima-ring* p n)))
2160
2161(defmfun $poly_content (p vars)
2162 (with-parsed-polynomials ((vars) :polynomials (p))
2163 (poly-content *maxima-ring* p)))
2164
2165(defmfun $poly_pseudo_divide (f fl vars
2166 &aux (vars (coerce-maxima-list vars))
2167 (f (parse-poly f vars))
2168 (fl (parse-poly-list fl vars)))
2169 (multiple-value-bind (quot rem c division-count)
2170 (poly-pseudo-divide *maxima-ring* f fl)
2171 `((mlist)
2172 ,(coerce-to-maxima :poly-list quot vars)
2173 ,(coerce-to-maxima :polynomial rem vars)
2174 ,c
2175 ,division-count)))
2176
2177(defmfun $poly_exact_divide (f g vars)
2178 (with-parsed-polynomials ((vars) :polynomials (f g) :value-type :polynomial)
2179 (poly-exact-divide *maxima-ring* f g)))
2180
2181(defmfun $poly_normal_form (f fl vars)
2182 (with-parsed-polynomials ((vars) :polynomials (f)
2183 :poly-lists (fl)
2184 :value-type :polynomial)
2185 (normal-form *maxima-ring* f (remzero fl) nil)))
2186
2187(defmfun $poly_buchberger_criterion (g vars)
2188 (with-parsed-polynomials ((vars) :poly-lists (g) :value-type :logical)
2189 (buchberger-criterion *maxima-ring* g)))
2190
2191(defmfun $poly_buchberger (fl vars)
2192 (with-parsed-polynomials ((vars) :poly-lists (fl) :value-type :poly-list)
2193 (buchberger *maxima-ring* (remzero fl) 0 nil)))
2194
2195(defmfun $poly_reduction (plist vars)
2196 (with-parsed-polynomials ((vars) :poly-lists (plist)
2197 :value-type :poly-list)
2198 (reduction *maxima-ring* plist)))
2199
2200(defmfun $poly_minimization (plist vars)
2201 (with-parsed-polynomials ((vars) :poly-lists (plist)
2202 :value-type :poly-list)
2203 (minimization plist)))
2204
2205(defmfun $poly_normalize_list (plist vars)
2206 (with-parsed-polynomials ((vars) :poly-lists (plist)
2207 :value-type :poly-list)
2208 (poly-normalize-list *maxima-ring* plist)))
2209
2210(defmfun $poly_grobner (f vars)
2211 (with-parsed-polynomials ((vars) :poly-lists (f)
2212 :value-type :poly-list)
2213 (grobner *maxima-ring* (remzero f))))
2214
2215(defmfun $poly_reduced_grobner (f vars)
2216 (with-parsed-polynomials ((vars) :poly-lists (f)
2217 :value-type :poly-list)
2218 (reduced-grobner *maxima-ring* (remzero f))))
2219
2220(defmfun $poly_depends_p (p var mvars
2221 &aux (vars (coerce-maxima-list mvars))
2222 (pos (position var vars)))
2223 (if (null pos)
2224 (merror "~%Variable ~M not in the list of variables ~M." var mvars)
2225 (poly-depends-p (parse-poly p vars) pos)))
2226
2227(defmfun $poly_elimination_ideal (flist k vars)
2228 (with-parsed-polynomials ((vars) :poly-lists (flist)
2229 :value-type :poly-list)
2230 (elimination-ideal *maxima-ring* flist k nil 0)))
2231
2232(defmfun $poly_colon_ideal (f g vars)
2233 (with-parsed-polynomials ((vars) :poly-lists (f g) :value-type :poly-list)
2234 (colon-ideal *maxima-ring* f g nil)))
2235
2236(defmfun $poly_ideal_intersection (f g vars)
2237 (with-parsed-polynomials ((vars) :poly-lists (f g) :value-type :poly-list)
2238 (ideal-intersection *maxima-ring* f g nil)))
2239
2240(defmfun $poly_lcm (f g vars)
2241 (with-parsed-polynomials ((vars) :polynomials (f g) :value-type :polynomial)
2242 (poly-lcm *maxima-ring* f g)))
2243
2244(defmfun $poly_gcd (f g vars)
2245 ($first ($divide (m* f g) ($poly_lcm f g vars))))
2246
2247(defmfun $poly_grobner_equal (g1 g2 vars)
2248 (with-parsed-polynomials ((vars) :poly-lists (g1 g2))
2249 (grobner-equal *maxima-ring* g1 g2)))
2250
2251(defmfun $poly_grobner_subsetp (g1 g2 vars)
2252 (with-parsed-polynomials ((vars) :poly-lists (g1 g2))
2253 (grobner-subsetp *maxima-ring* g1 g2)))
2254
2255(defmfun $poly_grobner_member (p g vars)
2256 (with-parsed-polynomials ((vars) :polynomials (p) :poly-lists (g))
2257 (grobner-member *maxima-ring* p g)))
2258
2259(defmfun $poly_ideal_saturation1 (f p vars)
2260 (with-parsed-polynomials ((vars) :poly-lists (f) :polynomials (p)
2261 :value-type :poly-list)
2262 (ideal-saturation-1 *maxima-ring* f p 0)))
2263
2264(defmfun $poly_saturation_extension (f plist vars new-vars)
2265 (with-parsed-polynomials ((vars new-vars)
2266 :poly-lists (f plist)
2267 :value-type :poly-list)
2268 (saturation-extension *maxima-ring* f plist)))
2269
2270(defmfun $poly_polysaturation_extension (f plist vars new-vars)
2271 (with-parsed-polynomials ((vars new-vars)
2272 :poly-lists (f plist)
2273 :value-type :poly-list)
2274 (polysaturation-extension *maxima-ring* f plist)))
2275
2276(defmfun $poly_ideal_polysaturation1 (f plist vars)
2277 (with-parsed-polynomials ((vars) :poly-lists (f plist)
2278 :value-type :poly-list)
2279 (ideal-polysaturation-1 *maxima-ring* f plist 0 nil)))
2280
2281(defmfun $poly_ideal_saturation (f g vars)
2282 (with-parsed-polynomials ((vars) :poly-lists (f g)
2283 :value-type :poly-list)
2284 (ideal-saturation *maxima-ring* f g 0 nil)))
2285
2286(defmfun $poly_ideal_polysaturation (f ideal-list vars)
2287 (with-parsed-polynomials ((vars) :poly-lists (f)
2288 :poly-list-lists (ideal-list)
2289 :value-type :poly-list)
2290 (ideal-polysaturation *maxima-ring* f ideal-list 0 nil)))
2291
2292(defmfun $poly_lt (f vars)
2293 (with-parsed-polynomials ((vars) :polynomials (f) :value-type :polynomial)
2294 (make-poly-from-termlist (list (poly-lt f)))))
2295
2296(defmfun $poly_lm (f vars)
2297 (with-parsed-polynomials ((vars) :polynomials (f) :value-type :polynomial)
2298 (make-poly-from-termlist (list (make-term (poly-lm f) (funcall (ring-unit *maxima-ring*)))))))
2299
Note: See TracBrowser for help on using the repository browser.