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

source: branches/f4grobner/monom.lisp@ 3571

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

* empty log message *

File size: 20.3 KB
Line 
1;;; -*- Mode: Lisp -*-
2;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3;;;
4;;; Copyright (C) 1999, 2002, 2009, 2015 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(defpackage "MONOM"
23 (:use :cl :copy)
24 (:export "MONOM"
25 "EXPONENT"
26 "MONOM-DIMENSION"
27 "MONOM-EXPONENTS"
28 "MONOM-EQUALP"
29 "MONOM-ELT"
30 "MONOM-TOTAL-DEGREE"
31 "MONOM-SUGAR"
32 "MONOM-MULTIPLY-BY"
33 "MONOM-DIVIDE-BY"
34 "MONOM-COPY-INSTANCE"
35 "MONOM-MULTIPLY-2"
36 "MONOM-MULTIPLY"
37 "MONOM-DIVIDES-P"
38 "MONOM-DIVIDES-LCM-P"
39 "MONOM-LCM-DIVIDES-LCM-P"
40 "MONOM-LCM-EQUAL-LCM-P"
41 "MONOM-DIVISIBLE-BY-P"
42 "MONOM-REL-PRIME-P"
43 "MONOM-LCM"
44 "MONOM-GCD"
45 "MONOM-DEPENDS-P"
46 "MONOM-LEFT-TENSOR-PRODUCT-BY"
47 "MONOM-RIGHT-TENSOR-PRODUCT-BY"
48 "MONOM-LEFT-CONTRACT"
49 "MAKE-MONOM-VARIABLE"
50 "MONOM->LIST"
51 "LEX>"
52 "GRLEX>"
53 "REVLEX>"
54 "GREVLEX>"
55 "INVLEX>"
56 "REVERSE-MONOMIAL-ORDER"
57 "MAKE-ELIMINATION-ORDER-FACTORY")
58 (:documentation
59 "This package implements basic operations on monomials, including
60various monomial orders.
61
62DATA STRUCTURES: Conceptually, monomials can be represented as lists:
63
64 monom: (n1 n2 ... nk) where ni are non-negative integers
65
66However, lists may be implemented as other sequence types, so the
67flexibility to change the representation should be maintained in the
68code to use general operations on sequences whenever possible. The
69optimization for the actual representation should be left to
70declarations and the compiler.
71
72EXAMPLES: Suppose that variables are x and y. Then
73
74 Monom x*y^2 ---> (1 2) "))
75
76(in-package :monom)
77
78(proclaim '(optimize (speed 3) (space 0) (safety 0) (debug 0)))
79
80(deftype exponent ()
81 "Type of exponent in a monomial."
82 'fixnum)
83
84(defclass monom ()
85 ((exponents :initarg :exponents :accessor monom-exponents
86 :documentation "The powers of the variables."))
87 ;; default-initargs are not needed, they are handled by SHARED-INITIALIZE
88 ;;(:default-initargs :dimension 'foo :exponents 'bar :exponent 'baz)
89 (:documentation
90 "Implements a monomial, i.e. a product of powers
91of variables, like X*Y^2."))
92
93(defmethod print-object ((self monom) stream)
94 (print-unreadable-object (self stream :type t :identity t)
95 (with-accessors ((exponents monom-exponents))
96 self
97 (format stream "EXPONENTS=~A"
98 exponents))))
99
100(defmethod initialize-instance :after ((self monom)
101 &key
102 (dimension 0 dimension-supplied-p)
103 (exponents nil exponents-supplied-p)
104 (exponent 0)
105 &allow-other-keys
106 )
107 "The following INITIALIZE-INSTANCE method allows instance initialization
108of a MONOM in a style similar to MAKE-ARRAY, e.g.:
109
110 (MAKE-INSTANCE :EXPONENTS '(1 2 3)) --> #<MONOM EXPONENTS=#(1 2 3)>
111 (MAKE-INSTANCE :DIMENSION 3) --> #<MONOM EXPONENTS=#(0 0 0)>
112 (MAKE-INSTANCE :DIMENSION 3 :EXPONENT 7) --> #<MONOM EXPONENTS=#(7 7 7)>
113
114If both DIMENSION and EXPONENTS are supplied, they must be compatible,
115i.e. the length of EXPONENTS must be equal DIMENSION. If EXPONENTS
116is not supplied, a monom with repeated value EXPONENT is created.
117By default EXPONENT is 0, which results in a constant monomial.
118"
119 (cond
120 (exponents-supplied-p
121 (when (and dimension-supplied-p
122 (/= dimension (length exponents)))
123 (error "EXPONENTS (~A) must have supplied length DIMENSION (~A)"
124 exponents dimension))
125 (let ((dim (length exponents)))
126 (setf (slot-value self 'exponents) (make-array dim :initial-contents exponents))))
127 (dimension-supplied-p
128 ;; when all exponents are to be identical
129 (setf (slot-value self 'exponents) (make-array (list dimension)
130 :initial-element exponent
131 :element-type 'exponent)))
132 (t
133 (error "Initarg DIMENSION or EXPONENTS must be supplied."))))
134
135(defgeneric monom-dimension (m)
136 (:method ((m monom))
137 (length (monom-exponents m))))
138
139(defgeneric universal-equalp (object1 object2)
140 (:documentation "Returns T iff OBJECT1 and OBJECT2 are equal.")
141 (:method ((m1 monom) (m2 monom))
142 "Returns T iff monomials M1 and M2 have identical EXPONENTS."
143 (equalp (monom-exponents m1) (monom-exponents m2))))
144
145(defgeneric monom-elt (m index)
146 (:documentation "Return the power in the monomial M of variable number INDEX."
147 (:method ((m monom) index)
148 "Return the power in the monomial M of variable number INDEX."
149 (with-slots (exponents)
150 m
151 (elt exponents index))))
152
153(defgeneric (setf monom-elt) (new-value m index)
154 (:documentation "Set the power in the monomial M of variable number INDEX.")
155 (:method (new-value (m monom) index)
156 (with-slots (exponents)
157 m
158 (setf (elt exponents index) new-value))))
159
160(defgeneric total-degree (m &optional start end)
161 (:documentation "Return the total degree of a monomoal M. Optinally, a range
162of variables may be specified with arguments START and END.")
163 (:method ((m monom) &optional (start 0) (end (monom-dimension m)))
164 (declare (type fixnum start end))
165 (with-slots (exponents)
166 m
167 (reduce #'+ exponents :start start :end end))))
168
169(defgeneric sugar (m &optional start end)
170 (:documentation "Return the sugar of a monomial M. Optinally, a range
171of variables may be specified with arguments START and END.")
172 (:method ((m monom) &optional (start 0) (end (monom-dimension m)))
173 (declare (type fixnum start end))
174 (total-degree m start end)))
175
176(defgeneric multiply-by (self other)
177 (:documentation "Multiply SELF by OTHER, return SELF.")
178 (:method ((self monom) (other monom))
179 (with-slots ((exponents1 exponents))
180 self
181 (with-slots ((exponents2 exponents))
182 other
183 (unless (= (length exponents1) (length exponents2))
184 (error "Incompatible dimensions"))
185 (map-into exponents1 #'+ exponents1 exponents2)))
186 self))
187
188(defgeneric divide-by (self other)
189 (:documentation "Divide SELF by OTHER, return SELF.")
190 (:method ((self monom) (other monom))
191 (with-slots ((exponents1 exponents))
192 self
193 (with-slots ((exponents2 exponents))
194 other
195 (unless (= (length exponents1) (length exponents2))
196 (error "divide-by: Incompatible dimensions."))
197 (unless (every #'>= exponents1 exponents2)
198 (error "divide-by: Negative power would result."))
199 (map-into exponents1 #'- exponents1 exponents2)))
200 self))
201
202(defmethod copy-instance :around ((object monom) &rest initargs &key &allow-other-keys)
203 "An :AROUND method of COPY-INSTANCE. It replaces
204exponents with a fresh copy of the sequence."
205 (declare (ignore object initargs))
206 (let ((copy (call-next-method)))
207 (setf (monom-exponents copy) (copy-seq (monom-exponents copy)))
208 copy))
209
210(defun multiply-2 (object1 object2)
211 "Multiply OBJECT1 by OBJECT2"
212 (multiply-by (copy-instance object1) (copy-instance object2)))
213
214(defun multiply (&rest factors)
215 "Non-destructively multiply list FACTORS."
216 (reduce #'multiply-2 factors))
217
218(defun divide (numerator &rest denominators)
219 "Non-destructively divide object NUMERATOR by product of DENOMINATORS."
220 (divide-by (copy-instance numerator) (multiply denominators)))
221
222(defmethod monom-divides-p ((m1 monom) (m2 monom))
223 "Returns T if monomial M1 divides monomial M2, NIL otherwise."
224 (with-slots ((exponents1 exponents))
225 m1
226 (with-slots ((exponents2 exponents))
227 m2
228 (every #'<= exponents1 exponents2))))
229
230(defmethod monom-divides-lcm-p ((m1 monom) (m2 monom) (m3 monom))
231 "Returns T if monomial M1 divides LCM(M2,M3), NIL otherwise."
232 (every #'(lambda (x y z) (<= x (max y z)))
233 m1 m2 m3))
234
235(defmethod monom-lcm-divides-lcm-p ((m1 monom) (m2 monom) (m3 monom) (m4 monom))
236 "Returns T if monomial MONOM-LCM(M1,M2) divides MONOM-LCM(M3,M4), NIL otherwise."
237 (declare (type monom m1 m2 m3 m4))
238 (every #'(lambda (x y z w) (<= (max x y) (max z w)))
239 m1 m2 m3 m4))
240
241(defmethod monom-lcm-equal-lcm-p ((m1 monom) (m2 monom) (m3 monom) (m4 monom))
242 "Returns T if monomial LCM(M1,M2) equals LCM(M3,M4), NIL otherwise."
243 (with-slots ((exponents1 exponents))
244 m1
245 (with-slots ((exponents2 exponents))
246 m2
247 (with-slots ((exponents3 exponents))
248 m3
249 (with-slots ((exponents4 exponents))
250 m4
251 (every
252 #'(lambda (x y z w) (= (max x y) (max z w)))
253 exponents1 exponents2 exponents3 exponents4))))))
254
255(defgeneric divisible-by-p (object1 object2)
256 (:documentation "Return T if OBJECT1 is divisible by OBJECT2.")
257 (:method ((m1 monom) (m2 monom))
258 "Returns T if monomial M1 is divisible by monomial M2, NIL otherwise."
259 (with-slots ((exponents1 exponents))
260 m1
261 (with-slots ((exponents2 exponents))
262 m2
263 (every #'>= exponents1 exponents2)))))
264
265(defgeneric rel-prime-p (object1 object2)
266 "Returns T if objects OBJECT1 and OBJECT2 are relatively prime."
267 (:method ((m1 monom) (m2 monom))
268 "Returns T if two monomials M1 and M2 are relatively prime (disjoint)."
269 (with-slots ((exponents1 exponents))
270 m1
271 (with-slots ((exponents2 exponents))
272 m2
273 (every #'(lambda (x y) (zerop (min x y))) exponents1 exponents2)))))
274
275(defgeneric monom-lcm (object1 object2)
276 (:documentation "Returns the multiple of objects OBJECT1 and OBJECT2.")
277 (:method ((m1 monom) (m2 monom))
278 "Returns least common multiple of monomials M1 and M2."
279 (with-slots ((exponents1 exponents))
280 m1
281 (with-slots ((exponents2 exponents))
282 m2
283 (let* ((exponents (copy-seq exponents1)))
284 (map-into exponents #'max exponents1 exponents2)
285 (make-instance 'monom :exponents exponents))))))
286
287
288(defgeneric universal-gcd (object1 object2)
289 (:documentation "Returns GCD of objects OBJECT1 and OBJECT2")
290 (:method ((m1 monom) (m2 monom))
291 "Returns greatest common divisor of monomials M1 and M2."
292 (with-slots ((exponents1 exponents))
293 m1
294 (with-slots ((exponents2 exponents))
295 m2
296 (let* ((exponents (copy-seq exponents1)))
297 (map-into exponents #'min exponents1 exponents2)
298 (make-instance 'monom :exponents exponents))))))
299
300(defgeneric depends-p (object k)
301 (:documentation "Returns T iff object OBJECT depends on variable K.")
302 (:method ((m monom) k)
303 "Return T if the monomial M depends on variable number K."
304 (declare (type fixnum k))
305 (with-slots (exponents)
306 m
307 (plusp (elt exponents k)))))
308
309(defgeneric left-tensor-product-by (self other)
310 (:documentation "Returns a tensor product SELF by OTHER, stored into
311 SELF. Return SELF.")
312 (:method ((self monom) (other monom))
313 (with-slots ((exponents1 exponents))
314 self
315 (with-slots ((exponents2 exponents))
316 other
317 (setf exponents1 (concatenate 'vector exponents2 exponents1))))
318 self))
319
320(defgeneric right-tensor-product-by (self other)
321 (:documentation "Returns a tensor product of OTHER by SELF, stored
322 into SELF. Returns SELF.")
323 (:method ((self monom) (other monom))
324 (with-slots ((exponents1 exponents))
325 self
326 (with-slots ((exponents2 exponents))
327 other
328 (setf exponents1 (concatenate 'vector exponents1 exponents2))))
329 self))
330
331(defgeneric left-contract (self k)
332 (:documentation "Drop the first K variables in object SELF.")
333 (:method ((self monom) k)
334 "Drop the first K variables in monomial M."
335 (declare (fixnum k))
336 (with-slots (exponents)
337 self
338 (setf exponents (subseq exponents k)))
339 self))
340
341(defun make-monom-variable (nvars pos &optional (power 1)
342 &aux (m (make-instance 'monom :dimension nvars)))
343 "Construct a monomial in the polynomial ring
344RING[X[0],X[1],X[2],...X[NVARS-1]] over the (unspecified) ring RING
345which represents a single variable. It assumes number of variables
346NVARS and the variable is at position POS. Optionally, the variable
347may appear raised to power POWER. "
348 (declare (type fixnum nvars pos power) (type monom m))
349 (with-slots (exponents)
350 m
351 (setf (elt exponents pos) power)
352 m))
353
354(defmethod monom->list ((m monom))
355 "A human-readable representation of a monomial M as a list of exponents."
356 (coerce (monom-exponents m) 'list))
357
358
359;; pure lexicographic
360(defgeneric lex> (p q &optional start end)
361 (:documentation "Return T if P>Q with respect to lexicographic
362order, otherwise NIL. The second returned value is T if P=Q,
363otherwise it is NIL.")
364 (:method ((p monom) (q monom) &optional (start 0) (end (monom-dimension p)))
365 (declare (type fixnum start end))
366 (do ((i start (1+ i)))
367 ((>= i end) (values nil t))
368 (cond
369 ((> (monom-elt p i) (monom-elt q i))
370 (return-from lex> (values t nil)))
371 ((< (monom-elt p i) (monom-elt q i))
372 (return-from lex> (values nil nil)))))))
373
374;; total degree order, ties broken by lexicographic
375(defgeneric grlex> (p q &optional start end)
376 (:documentation "Return T if P>Q with respect to graded
377lexicographic order, otherwise NIL. The second returned value is T if
378P=Q, otherwise it is NIL.")
379 (:method ((p monom) (q monom) &optional (start 0) (end (monom-dimension p)))
380 (declare (type monom p q) (type fixnum start end))
381 (let ((d1 (monom-total-degree p start end))
382 (d2 (monom-total-degree q start end)))
383 (declare (type fixnum d1 d2))
384 (cond
385 ((> d1 d2) (values t nil))
386 ((< d1 d2) (values nil nil))
387 (t
388 (lex> p q start end))))))
389
390;; reverse lexicographic
391(defgeneric revlex> (p q &optional start end)
392 (:documentation "Return T if P>Q with respect to reverse
393lexicographic order, NIL otherwise. The second returned value is T if
394P=Q, otherwise it is NIL. This is not and admissible monomial order
395because some sets do not have a minimal element. This order is useful
396in constructing other orders.")
397 (:method ((p monom) (q monom) &optional (start 0) (end (monom-dimension p)))
398 (declare (type fixnum start end))
399 (do ((i (1- end) (1- i)))
400 ((< i start) (values nil t))
401 (declare (type fixnum i))
402 (cond
403 ((< (monom-elt p i) (monom-elt q i))
404 (return-from revlex> (values t nil)))
405 ((> (monom-elt p i) (monom-elt q i))
406 (return-from revlex> (values nil nil)))))))
407
408
409;; total degree, ties broken by reverse lexicographic
410(defgeneric grevlex> (p q &optional start end)
411 (:documentation "Return T if P>Q with respect to graded reverse
412lexicographic order, NIL otherwise. The second returned value is T if
413P=Q, otherwise it is NIL.")
414 (:method ((p monom) (q monom) &optional (start 0) (end (monom-dimension p)))
415 (declare (type fixnum start end))
416 (let ((d1 (monom-total-degree p start end))
417 (d2 (monom-total-degree q start end)))
418 (declare (type fixnum d1 d2))
419 (cond
420 ((> d1 d2) (values t nil))
421 ((< d1 d2) (values nil nil))
422 (t
423 (revlex> p q start end))))))
424
425(defgeneric invlex> (p q &optional start end)
426 (:documentation "Return T if P>Q with respect to inverse
427lexicographic order, NIL otherwise The second returned value is T if
428P=Q, otherwise it is NIL.")
429 (:method ((p monom) (q monom) &optional (start 0) (end (monom-dimension p)))
430 (declare (type fixnum start end))
431 (do ((i (1- end) (1- i)))
432 ((< i start) (values nil t))
433 (declare (type fixnum i))
434 (cond
435 ((> (monom-elt p i) (monom-elt q i))
436 (return-from invlex> (values t nil)))
437 ((< (monom-elt p i) (monom-elt q i))
438 (return-from invlex> (values nil nil)))))))
439
440(defun reverse-monomial-order (order)
441 "Create the inverse monomial order to the given monomial order ORDER."
442 #'(lambda (p q &optional (start 0) (end (monom-dimension q)))
443 (declare (type monom p q) (type fixnum start end))
444 (funcall order q p start end)))
445
446;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
447;;
448;; Order making functions
449;;
450;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
451
452;; This returns a closure with the same signature
453;; as all orders such as #'LEX>.
454(defun make-elimination-order-factory-1 (&optional (secondary-elimination-order #'lex>))
455 "It constructs an elimination order used for the 1-st elimination ideal,
456i.e. for eliminating the first variable. Thus, the order compares the degrees of the
457first variable in P and Q first, with ties broken by SECONDARY-ELIMINATION-ORDER."
458 #'(lambda (p q &optional (start 0) (end (monom-dimension p)))
459 (declare (type monom p q) (type fixnum start end))
460 (cond
461 ((> (monom-elt p start) (monom-elt q start))
462 (values t nil))
463 ((< (monom-elt p start) (monom-elt q start))
464 (values nil nil))
465 (t
466 (funcall secondary-elimination-order p q (1+ start) end)))))
467
468;; This returns a closure which is called with an integer argument.
469;; The result is *another closure* with the same signature as all
470;; orders such as #'LEX>.
471(defun make-elimination-order-factory (&optional
472 (primary-elimination-order #'lex>)
473 (secondary-elimination-order #'lex>))
474 "Return a function with a single integer argument K. This should be
475the number of initial K variables X[0],X[1],...,X[K-1], which precede
476remaining variables. The call to the closure creates a predicate
477which compares monomials according to the K-th elimination order. The
478monomial orders PRIMARY-ELIMINATION-ORDER and
479SECONDARY-ELIMINATION-ORDER are used to compare the first K and the
480remaining variables, respectively, with ties broken by lexicographical
481order. That is, if PRIMARY-ELIMINATION-ORDER yields (VALUES NIL T),
482which indicates that the first K variables appear with identical
483powers, then the result is that of a call to
484SECONDARY-ELIMINATION-ORDER applied to the remaining variables
485X[K],X[K+1],..."
486 #'(lambda (k)
487 (cond
488 ((<= k 0)
489 (error "K must be at least 1"))
490 ((= k 1)
491 (make-elimination-order-factory-1 secondary-elimination-order))
492 (t
493 #'(lambda (p q &optional (start 0) (end (monom-dimension p)))
494 (declare (type monom p q) (type fixnum start end))
495 (multiple-value-bind (primary equal)
496 (funcall primary-elimination-order p q start k)
497 (if equal
498 (funcall secondary-elimination-order p q k end)
499 (values primary nil))))))))
500
501(defclass term (monom)
502 ((coeff :initarg :coeff :accessor term-coeff))
503 (:default-initargs :coeff nil)
504 (:documentation "Implements a term, i.e. a product of a scalar
505and powers of some variables, such as 5*X^2*Y^3."))
506
507(defmethod print-object ((self term) stream)
508 (print-unreadable-object (self stream :type t :identity t)
509 (with-accessors ((exponents monom-exponents)
510 (coeff term-coeff))
511 self
512 (format stream "EXPONENTS=~A COEFF=~A"
513 exponents coeff))))
514
515(defmethod universal-equalp ((term1 term) (term2 term))
516 "Returns T if TERM1 and TERM2 are equal as MONOM, and coefficients
517are UNIVERSAL-EQUALP."
518 (and (call-next-method)
519 (universal-equalp (term-coeff term1) (term-coeff term2))))
520
521(defmethod update-instance-for-different-class :after ((old monom) (new term) &key)
522 (setf (term-coeff new) 1))
523
524(defmethod multiply-by :before ((self term) (other term))
525 "Destructively multiply terms SELF and OTHER and store the result into SELF.
526It returns SELF."
527 (setf (term-coeff self) (multiply-by (term-coeff self) (scalar-coeff other))))
528
529(defmethod term-left-tensor-product-by :before ((self term) (other term))
530 (setf (term-coeff self) (universal-multiply-by (term-coeff self) (term-coeff other))))
531
532(defmethod term-right-tensor-product-by :before ((self term) (other term))
533 (setf (term-coeff self) (multiply-by (term-coeff self) (term-coeff other))))
534
535(defmethod divide-by :before ((self term) (other term))
536 (setf (term-coeff self) (divide-by (term-coeff self) (term-coeff other))))
537
538(defmethod monom-unary-minus ((self term))
539 (setf (term-coeff self) (monom-unary-minus (term-coeff self)))
540 self)
541
542(defmethod monom-multiply ((term1 term) (term2 term))
543 "Non-destructively multiply TERM1 by TERM2."
544 (monom-multiply-by (copy-instance term1) (copy-instance term2)))
545
546(defmethod monom-multiply ((term1 number) (term2 monom))
547 "Non-destructively multiply TERM1 by TERM2."
548 (monom-multiply term1 (change-class (copy-instance term2) 'term)))
549
550(defmethod monom-zerop ((self term))
551 (zerop (term-coeff self)))
Note: See TracBrowser for help on using the repository browser.