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

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

* empty log message *

File size: 10.2 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 :ring)
24 (:export "MONOM"
25 "EXPONENT"
26 "MONOM-DIMENSION"
27 "MONOM-EXPONENTS"
28 "MAKE-MONOM-VARIABLE")
29 (:documentation
30 "This package implements basic operations on monomials.
31DATA STRUCTURES: Conceptually, monomials can be represented as lists:
32
33 monom: (n1 n2 ... nk) where ni are non-negative integers
34
35However, lists may be implemented as other sequence types, so the
36flexibility to change the representation should be maintained in the
37code to use general operations on sequences whenever possible. The
38optimization for the actual representation should be left to
39declarations and the compiler.
40
41EXAMPLES: Suppose that variables are x and y. Then
42
43 Monom x*y^2 ---> (1 2) "))
44
45(in-package :monom)
46
47(proclaim '(optimize (speed 3) (space 0) (safety 0) (debug 0)))
48
49(deftype exponent ()
50 "Type of exponent in a monomial."
51 'fixnum)
52
53(defclass monom ()
54 ((exponents :initarg :exponents :accessor monom-exponents
55 :documentation "The powers of the variables."))
56 ;; default-initargs are not needed, they are handled by SHARED-INITIALIZE
57 ;;(:default-initargs :dimension 'foo :exponents 'bar :exponent 'baz)
58 (:documentation
59 "Implements a monomial, i.e. a product of powers
60of variables, like X*Y^2."))
61
62(defmethod print-object ((self monom) stream)
63 (print-unreadable-object (self stream :type t :identity t)
64 (with-accessors ((exponents monom-exponents))
65 self
66 (format stream "EXPONENTS=~A"
67 exponents))))
68
69;; The following INITIALIZE-INSTANCE method allows instance
70;; initialization in a style similar to MAKE-ARRAY, e.g.
71;;
72;; (MAKE-INSTANCE :EXPONENTS '(1 2 3)) --> #<MONOM EXPONENTS=#(1 2 3)>
73;; (MAKE-INSTANCE :DIMENSION 3) --> #<MONOM EXPONENTS=#(0 0 0)>
74;; (MAKE-INSTANCE :DIMENSION 3 :EXPONENT 7) --> #<MONOM EXPONENTS=#(7 7 7)>
75;;
76(defmethod initialize-instance :after ((self monom)
77 &key
78 (dimension 0 dimension-supplied-p)
79 (exponents nil exponents-supplied-p)
80 (exponent 0)
81 &allow-other-keys
82 )
83 (cond
84 (exponents-supplied-p
85 (when dimension-supplied-p
86 (warn "Ignoring initarg DIMENSION."))
87 (let ((dim (length exponents)))
88 (setf (slot-value self 'exponents) (make-array dim :initial-contents exponents))))
89 (dimension-supplied-p
90 ;; when all exponents are to be identical
91 (setf (slot-value self 'exponents) (make-array (list dimension)
92 :initial-element exponent
93 :element-type 'exponent)))
94 (t
95 (error "Initarg DIMENSION or EXPONENTS must be supplied."))))
96
97(defmethod monom-dimension ((m monom))
98 (length (monom-exponents m)))
99
100(defmethod r-equalp ((m1 monom) (m2 monom))
101 "Returns T iff monomials M1 and M2 have identical
102EXPONENTS."
103 (equalp (monom-exponents m1) (monom-exponents m2)))
104
105(defmethod r-coeff ((m monom))
106 "A MONOM can be treated as a special case of TERM,
107where the coefficient is 1."
108 1)
109
110(defmethod r-elt ((m monom) index)
111 "Return the power in the monomial M of variable number INDEX."
112 (with-slots (exponents)
113 m
114 (elt exponents index)))
115
116(defmethod (setf r-elt) (new-value (m monom) index)
117 "Return the power in the monomial M of variable number INDEX."
118 (with-slots (exponents)
119 m
120 (setf (elt exponents index) new-value)))
121
122(defmethod r-total-degree ((m monom) &optional (start 0) (end (monom-dimension m)))
123 "Return the todal degree of a monomoal M. Optinally, a range
124of variables may be specified with arguments START and END."
125 (declare (type fixnum start end))
126 (with-slots (exponents)
127 m
128 (reduce #'+ exponents :start start :end end)))
129
130
131(defmethod r-sugar ((m monom) &aux (start 0) (end (monom-dimension m)))
132 "Return the sugar of a monomial M. Optinally, a range
133of variables may be specified with arguments START and END."
134 (declare (type fixnum start end))
135 (r-total-degree m start end))
136
137(defmethod multiply-by ((self monom) (other monom))
138 (with-slots ((exponents1 exponents))
139 self
140 (with-slots ((exponents2 exponents))
141 other
142 (unless (= (length exponents1) (length exponents2))
143 (error "Incompatible dimensions"))
144 (map-into exponents1 #'+ exponents1 exponents2)))
145 self)
146
147(defmethod divide-by ((self monom) (other monom))
148 (with-slots ((exponents1 exponents))
149 self
150 (with-slots ((exponents2 exponents))
151 other
152 (unless (= (length exponents1) (length exponents2))
153 (error "Incompatible dimensions"))
154 (map-into exponents1 #'- exponents1 exponents2)))
155 self)
156
157(defmethod copy-instance :around ((object monom) &rest initargs &key &allow-other-keys)
158 "An :AROUNT method for COPY-INSTANCE. The primary method is a shallow copy,
159 while for monomials we typically need a fresh copy of the
160 exponents."
161 (declare (ignore object initargs))
162 (let ((copy (call-next-method)))
163 (setf (monom-exponents copy) (copy-seq (monom-exponents copy)))
164 copy))
165
166(defmethod r* ((m1 monom) (m2 monom))
167 "Non-destructively multiply monomial M1 by M2."
168 (multiply-by (copy-instance m1) (copy-instance m2)))
169
170(defmethod r/ ((m1 monom) (m2 monom))
171 "Non-destructively divide monomial M1 by monomial M2."
172 (divide-by (copy-instance m1) (copy-instance m2)))
173
174(defmethod r-divides-p ((m1 monom) (m2 monom))
175 "Returns T if monomial M1 divides monomial M2, NIL otherwise."
176 (with-slots ((exponents1 exponents))
177 m1
178 (with-slots ((exponents2 exponents))
179 m2
180 (every #'<= exponents1 exponents2))))
181
182
183(defmethod r-divides-lcm-p ((m1 monom) (m2 monom) (m3 monom))
184 "Returns T if monomial M1 divides LCM(M2,M3), NIL otherwise."
185 (every #'(lambda (x y z) (<= x (max y z)))
186 m1 m2 m3))
187
188
189(defmethod r-lcm-divides-lcm-p ((m1 monom) (m2 monom) (m3 monom) (m4 monom))
190 "Returns T if monomial MONOM-LCM(M1,M2) divides MONOM-LCM(M3,M4), NIL otherwise."
191 (declare (type monom m1 m2 m3 m4))
192 (every #'(lambda (x y z w) (<= (max x y) (max z w)))
193 m1 m2 m3 m4))
194
195(defmethod r-lcm-equal-lcm-p (m1 m2 m3 m4)
196 "Returns T if monomial LCM(M1,M2) equals LCM(M3,M4), NIL otherwise."
197 (with-slots ((exponents1 exponents))
198 m1
199 (with-slots ((exponents2 exponents))
200 m2
201 (with-slots ((exponents3 exponents))
202 m3
203 (with-slots ((exponents4 exponents))
204 m4
205 (every
206 #'(lambda (x y z w) (= (max x y) (max z w)))
207 exponents1 exponents2 exponents3 exponents4))))))
208
209(defmethod r-divisible-by-p ((m1 monom) (m2 monom))
210 "Returns T if monomial M1 is divisible by monomial M2, NIL otherwise."
211 (with-slots ((exponents1 exponents))
212 m1
213 (with-slots ((exponents2 exponents))
214 m2
215 (every #'>= exponents1 exponents2))))
216
217(defmethod r-rel-prime-p ((m1 monom) (m2 monom))
218 "Returns T if two monomials M1 and M2 are relatively prime (disjoint)."
219 (with-slots ((exponents1 exponents))
220 m1
221 (with-slots ((exponents2 exponents))
222 m2
223 (every #'(lambda (x y) (zerop (min x y))) exponents1 exponents2))))
224
225
226(defmethod r-lcm ((m1 monom) (m2 monom))
227 "Returns least common multiple of monomials M1 and M2."
228 (with-slots ((exponents1 exponents))
229 m1
230 (with-slots ((exponents2 exponents))
231 m2
232 (let* ((exponents (copy-seq exponents1)))
233 (map-into exponents #'max exponents1 exponents2)
234 (make-instance 'monom :exponents exponents)))))
235
236
237(defmethod r-gcd ((m1 monom) (m2 monom))
238 "Returns greatest common divisor of monomials M1 and M2."
239 (with-slots ((exponents1 exponents))
240 m1
241 (with-slots ((exponents2 exponents))
242 m2
243 (let* ((exponents (copy-seq exponents1)))
244 (map-into exponents #'min exponents1 exponents2)
245 (make-instance 'monom :exponents exponents)))))
246
247(defmethod r-depends-p ((m monom) k)
248 "Return T if the monomial M depends on variable number K."
249 (declare (type fixnum k))
250 (with-slots (exponents)
251 m
252 (plusp (elt exponents k))))
253
254(defmethod left-tensor-product-by ((self monom) (other monom))
255 (with-slots ((exponents1 exponents))
256 self
257 (with-slots ((exponents2 exponents))
258 other
259 (setf exponents1 (concatenate 'vector exponents2 exponents1))))
260 self)
261
262(defmethod right-tensor-product-by ((self monom) (other monom))
263 (with-slots ((exponents1 exponents))
264 self
265 (with-slots ((exponents2 exponents))
266 other
267 (setf exponents1 (concatenate 'vector exponents1 exponents2))))
268 self)
269
270(defmethod left-contract ((self monom) k)
271 "Drop the first K variables in monomial M."
272 (declare (fixnum k))
273 (with-slots (exponents)
274 self
275 (setf exponents (subseq exponents k)))
276 self)
277
278(defun make-monom-variable (nvars pos &optional (power 1)
279 &aux (m (make-instance 'monom :dimension nvars)))
280 "Construct a monomial in the polynomial ring
281RING[X[0],X[1],X[2],...X[NVARS-1]] over the (unspecified) ring RING
282which represents a single variable. It assumes number of variables
283NVARS and the variable is at position POS. Optionally, the variable
284may appear raised to power POWER. "
285 (declare (type fixnum nvars pos power) (type monom m))
286 (with-slots (exponents)
287 m
288 (setf (elt exponents pos) power)
289 m))
290
291(defmethod r->list ((m monom))
292 "A human-readable representation of a monomial M as a list of exponents."
293 (coerce (monom-exponents m) 'list))
294
295(defmethod r-dimension ((self monom))
296 (monom-dimension self))
297
298(defmethod r-exponents ((self monom))
299 (monom-exponents self))
Note: See TracBrowser for help on using the repository browser.