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