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