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