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