[1201] | 1 | ;;; -*- Mode: Lisp -*-
|
---|
[77] | 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 |
|
---|
[431] | 22 | (defpackage "POLYNOMIAL"
|
---|
[3055] | 23 | (:use :cl :utils :ring :monom :order :term #| :infix |# )
|
---|
[2596] | 24 | (:export "POLY"
|
---|
| 25 | "POLY-TERMLIST"
|
---|
[3016] | 26 | "POLY-TERM-ORDER"
|
---|
[3071] | 27 | "CHANGE-TERM-ORDER"
|
---|
[3099] | 28 | "STANDARD-EXTENSION"
|
---|
[3101] | 29 | "STANDARD-EXTENSION-1"
|
---|
[3094] | 30 | "SATURATION-EXTENSION"
|
---|
| 31 | "ALIST->POLY")
|
---|
[2522] | 32 | (:documentation "Implements polynomials"))
|
---|
[143] | 33 |
|
---|
[431] | 34 | (in-package :polynomial)
|
---|
| 35 |
|
---|
[1927] | 36 | (proclaim '(optimize (speed 3) (space 0) (safety 0) (debug 0)))
|
---|
[52] | 37 |
|
---|
[2442] | 38 | (defclass poly ()
|
---|
[2697] | 39 | ((termlist :initarg :termlist :accessor poly-termlist
|
---|
| 40 | :documentation "List of terms.")
|
---|
| 41 | (order :initarg :order :accessor poly-term-order
|
---|
| 42 | :documentation "Monomial/term order."))
|
---|
[2695] | 43 | (:default-initargs :termlist nil :order #'lex>)
|
---|
| 44 | (:documentation "A polynomial with a list of terms TERMLIST, ordered
|
---|
[2696] | 45 | according to term order ORDER, which defaults to LEX>."))
|
---|
[2442] | 46 |
|
---|
[2471] | 47 | (defmethod print-object ((self poly) stream)
|
---|
[2600] | 48 | (format stream "#<POLY TERMLIST=~A ORDER=~A>"
|
---|
[2595] | 49 | (poly-termlist self)
|
---|
| 50 | (poly-term-order self)))
|
---|
[2469] | 51 |
|
---|
[3015] | 52 | (defgeneric change-term-order (self other)
|
---|
[3012] | 53 | (:documentation "Change term order of SELF to the term order of OTHER.")
|
---|
[3010] | 54 | (:method ((self poly) (other poly))
|
---|
| 55 | (unless (eq (poly-term-order self) (poly-term-order other))
|
---|
| 56 | (setf (poly-termlist self) (sort (poly-termlist self) (poly-term-order other))
|
---|
| 57 | (poly-term-order self) (poly-term-order other)))
|
---|
[3012] | 58 | self))
|
---|
[3010] | 59 |
|
---|
[3095] | 60 | (defun alist->poly (alist &aux (poly (make-instance 'poly)))
|
---|
[3093] | 61 | "It reads polynomial from an alist formatted as ( ... (exponents . coeff) ...)."
|
---|
[3099] | 62 | (dolist (x alist poly)
|
---|
[3095] | 63 | (insert-item poly (make-instance 'term :exponents (car x) :coeff (cdr x)))))
|
---|
[3092] | 64 |
|
---|
| 65 |
|
---|
[2650] | 66 | (defmethod r-equalp ((self poly) (other poly))
|
---|
[2680] | 67 | "POLY instances are R-EQUALP if they have the same
|
---|
| 68 | order and if all terms are R-EQUALP."
|
---|
[2651] | 69 | (and (every #'r-equalp (poly-termlist self) (poly-termlist other))
|
---|
| 70 | (eq (poly-term-order self) (poly-term-order other))))
|
---|
[2650] | 71 |
|
---|
[2513] | 72 | (defmethod insert-item ((self poly) (item term))
|
---|
| 73 | (push item (poly-termlist self))
|
---|
[2514] | 74 | self)
|
---|
[2464] | 75 |
|
---|
[2513] | 76 | (defmethod append-item ((self poly) (item term))
|
---|
| 77 | (setf (cdr (last (poly-termlist self))) (list item))
|
---|
| 78 | self)
|
---|
[2466] | 79 |
|
---|
[52] | 80 | ;; Leading term
|
---|
[2442] | 81 | (defgeneric leading-term (object)
|
---|
| 82 | (:method ((self poly))
|
---|
[2525] | 83 | (car (poly-termlist self)))
|
---|
| 84 | (:documentation "The leading term of a polynomial, or NIL for zero polynomial."))
|
---|
[52] | 85 |
|
---|
| 86 | ;; Second term
|
---|
[2442] | 87 | (defgeneric second-leading-term (object)
|
---|
| 88 | (:method ((self poly))
|
---|
[2525] | 89 | (cadar (poly-termlist self)))
|
---|
| 90 | (:documentation "The second leading term of a polynomial, or NIL for a polynomial with at most one term."))
|
---|
[52] | 91 |
|
---|
| 92 | ;; Leading coefficient
|
---|
[2442] | 93 | (defgeneric leading-coefficient (object)
|
---|
| 94 | (:method ((self poly))
|
---|
[2526] | 95 | (r-coeff (leading-term self)))
|
---|
[2545] | 96 | (:documentation "The leading coefficient of a polynomial. It signals error for a zero polynomial."))
|
---|
[52] | 97 |
|
---|
| 98 | ;; Second coefficient
|
---|
[2442] | 99 | (defgeneric second-leading-coefficient (object)
|
---|
| 100 | (:method ((self poly))
|
---|
[2526] | 101 | (r-coeff (second-leading-term self)))
|
---|
[2906] | 102 | (:documentation "The second leading coefficient of a polynomial. It
|
---|
| 103 | signals error for a polynomial with at most one term."))
|
---|
[52] | 104 |
|
---|
| 105 | ;; Testing for a zero polynomial
|
---|
[2445] | 106 | (defmethod r-zerop ((self poly))
|
---|
| 107 | (null (poly-termlist self)))
|
---|
[52] | 108 |
|
---|
| 109 | ;; The number of terms
|
---|
[2445] | 110 | (defmethod r-length ((self poly))
|
---|
| 111 | (length (poly-termlist self)))
|
---|
[52] | 112 |
|
---|
[2483] | 113 | (defmethod multiply-by ((self poly) (other monom))
|
---|
[2501] | 114 | (mapc #'(lambda (term) (multiply-by term other))
|
---|
| 115 | (poly-termlist self))
|
---|
[2483] | 116 | self)
|
---|
[2469] | 117 |
|
---|
[2501] | 118 | (defmethod multiply-by ((self poly) (other scalar))
|
---|
[2502] | 119 | (mapc #'(lambda (term) (multiply-by term other))
|
---|
[2501] | 120 | (poly-termlist self))
|
---|
[2487] | 121 | self)
|
---|
| 122 |
|
---|
[2607] | 123 |
|
---|
[2761] | 124 | (defmacro fast-add/subtract (p q order-fn add/subtract-fn uminus-fn)
|
---|
[2755] | 125 | "Return an expression which will efficiently adds/subtracts two
|
---|
| 126 | polynomials, P and Q. The addition/subtraction of coefficients is
|
---|
| 127 | performed by calling ADD/SUBTRACT-METHOD-NAME. If UMINUS-METHOD-NAME
|
---|
| 128 | is supplied, it is used to negate the coefficients of Q which do not
|
---|
[2756] | 129 | have a corresponding coefficient in P. The code implements an
|
---|
| 130 | efficient algorithm to add two polynomials represented as sorted lists
|
---|
| 131 | of terms. The code destroys both arguments, reusing the terms to build
|
---|
| 132 | the result."
|
---|
[2742] | 133 | `(macrolet ((lc (x) `(r-coeff (car ,x))))
|
---|
| 134 | (do ((p ,p)
|
---|
| 135 | (q ,q)
|
---|
| 136 | r)
|
---|
| 137 | ((or (endp p) (endp q))
|
---|
| 138 | ;; NOTE: R contains the result in reverse order. Can it
|
---|
| 139 | ;; be more efficient to produce the terms in correct order?
|
---|
[2774] | 140 | (unless (endp q)
|
---|
[2776] | 141 | ;; Upon subtraction, we must change the sign of
|
---|
| 142 | ;; all coefficients in q
|
---|
[2774] | 143 | ,@(when uminus-fn
|
---|
[2775] | 144 | `((mapc #'(lambda (x) (setf x (funcall ,uminus-fn x))) q)))
|
---|
[2774] | 145 | (setf r (nreconc r q)))
|
---|
[2742] | 146 | r)
|
---|
| 147 | (multiple-value-bind
|
---|
| 148 | (greater-p equal-p)
|
---|
[2766] | 149 | (funcall ,order-fn (car p) (car q))
|
---|
[2742] | 150 | (cond
|
---|
| 151 | (greater-p
|
---|
| 152 | (rotatef (cdr p) r p)
|
---|
| 153 | )
|
---|
| 154 | (equal-p
|
---|
[2766] | 155 | (let ((s (funcall ,add/subtract-fn (lc p) (lc q))))
|
---|
[2742] | 156 | (cond
|
---|
| 157 | ((r-zerop s)
|
---|
| 158 | (setf p (cdr p))
|
---|
| 159 | )
|
---|
| 160 | (t
|
---|
| 161 | (setf (lc p) s)
|
---|
| 162 | (rotatef (cdr p) r p))))
|
---|
| 163 | (setf q (cdr q))
|
---|
| 164 | )
|
---|
| 165 | (t
|
---|
[2743] | 166 | ;;Negate the term of Q if UMINUS provided, signallig
|
---|
| 167 | ;;that we are doing subtraction
|
---|
[2908] | 168 | ,(when uminus-fn
|
---|
| 169 | `(setf (lc q) (funcall ,uminus-fn (lc q))))
|
---|
[2743] | 170 | (rotatef (cdr q) r q)))))))
|
---|
[2585] | 171 |
|
---|
[2655] | 172 |
|
---|
[2763] | 173 | (defmacro def-add/subtract-method (add/subtract-method-name
|
---|
[2752] | 174 | uminus-method-name
|
---|
| 175 | &optional
|
---|
[2913] | 176 | (doc-string nil doc-string-supplied-p))
|
---|
[2615] | 177 | "This macro avoids code duplication for two similar operations: ADD-TO and SUBTRACT-FROM."
|
---|
[2749] | 178 | `(defmethod ,add/subtract-method-name ((self poly) (other poly))
|
---|
[2615] | 179 | ,@(when doc-string-supplied-p `(,doc-string))
|
---|
[2769] | 180 | ;; Ensure orders are compatible
|
---|
[3015] | 181 | (change-term-order other self)
|
---|
[2772] | 182 | (setf (poly-termlist self) (fast-add/subtract
|
---|
| 183 | (poly-termlist self) (poly-termlist other)
|
---|
| 184 | (poly-term-order self)
|
---|
| 185 | #',add/subtract-method-name
|
---|
| 186 | ,(when uminus-method-name `(function ,uminus-method-name))))
|
---|
[2609] | 187 | self))
|
---|
[2487] | 188 |
|
---|
[2916] | 189 | (eval-when (:compile-toplevel :load-toplevel :execute)
|
---|
[2777] | 190 |
|
---|
| 191 | (def-add/subtract-method add-to nil
|
---|
| 192 | "Adds to polynomial SELF another polynomial OTHER.
|
---|
[2610] | 193 | This operation destructively modifies both polynomials.
|
---|
| 194 | The result is stored in SELF. This implementation does
|
---|
[2752] | 195 | no consing, entirely reusing the sells of SELF and OTHER.")
|
---|
[2609] | 196 |
|
---|
[2777] | 197 | (def-add/subtract-method subtract-from unary-minus
|
---|
[2753] | 198 | "Subtracts from polynomial SELF another polynomial OTHER.
|
---|
[2610] | 199 | This operation destructively modifies both polynomials.
|
---|
| 200 | The result is stored in SELF. This implementation does
|
---|
[2752] | 201 | no consing, entirely reusing the sells of SELF and OTHER.")
|
---|
[2610] | 202 |
|
---|
[2916] | 203 | )
|
---|
[2777] | 204 |
|
---|
[2916] | 205 |
|
---|
| 206 |
|
---|
[2691] | 207 | (defmethod unary-minus ((self poly))
|
---|
[2694] | 208 | "Destructively modifies the coefficients of the polynomial SELF,
|
---|
| 209 | by changing their sign."
|
---|
[2692] | 210 | (mapc #'unary-minus (poly-termlist self))
|
---|
[2683] | 211 | self)
|
---|
[52] | 212 |
|
---|
[2795] | 213 | (defun add-termlists (p q order-fn)
|
---|
[2794] | 214 | "Destructively adds two termlists P and Q ordered according to ORDER-FN."
|
---|
[2917] | 215 | (fast-add/subtract p q order-fn #'add-to nil))
|
---|
[2794] | 216 |
|
---|
[2800] | 217 | (defmacro multiply-term-by-termlist-dropping-zeros (term termlist
|
---|
[2927] | 218 | &optional (reverse-arg-order-P nil))
|
---|
[2799] | 219 | "Multiplies term TERM by a list of term, TERMLIST.
|
---|
[2792] | 220 | Takes into accound divisors of zero in the ring, by
|
---|
[2927] | 221 | deleting zero terms. Optionally, if REVERSE-ARG-ORDER-P
|
---|
[2928] | 222 | is T, change the order of arguments; this may be important
|
---|
[2927] | 223 | if we extend the package to non-commutative rings."
|
---|
[2800] | 224 | `(mapcan #'(lambda (other-term)
|
---|
[2907] | 225 | (let ((prod (r*
|
---|
[2923] | 226 | ,@(cond
|
---|
[2930] | 227 | (reverse-arg-order-p
|
---|
[2925] | 228 | `(other-term ,term))
|
---|
| 229 | (t
|
---|
| 230 | `(,term other-term))))))
|
---|
[2800] | 231 | (cond
|
---|
| 232 | ((r-zerop prod) nil)
|
---|
| 233 | (t (list prod)))))
|
---|
| 234 | ,termlist))
|
---|
[2790] | 235 |
|
---|
[2796] | 236 | (defun multiply-termlists (p q order-fn)
|
---|
[2787] | 237 | (cond
|
---|
[2917] | 238 | ((or (endp p) (endp q))
|
---|
| 239 | ;;p or q is 0 (represented by NIL)
|
---|
| 240 | nil)
|
---|
[2789] | 241 | ;; If p= p0+p1 and q=q0+q1 then p*q=p0*q0+p0*q1+p1*q
|
---|
[2787] | 242 | ((endp (cdr p))
|
---|
[2918] | 243 | (multiply-term-by-termlist-dropping-zeros (car p) q))
|
---|
| 244 | ((endp (cdr q))
|
---|
[2919] | 245 | (multiply-term-by-termlist-dropping-zeros (car q) p t))
|
---|
| 246 | (t
|
---|
[2948] | 247 | (cons (r* (car p) (car q))
|
---|
[2949] | 248 | (add-termlists
|
---|
| 249 | (multiply-term-by-termlist-dropping-zeros (car p) (cdr q))
|
---|
| 250 | (multiply-termlists (cdr p) q order-fn)
|
---|
| 251 | order-fn)))))
|
---|
[2793] | 252 |
|
---|
[2803] | 253 | (defmethod multiply-by ((self poly) (other poly))
|
---|
[3014] | 254 | (change-term-order other self)
|
---|
[2803] | 255 | (setf (poly-termlist self) (multiply-termlists (poly-termlist self)
|
---|
| 256 | (poly-termlist other)
|
---|
| 257 | (poly-term-order self)))
|
---|
| 258 | self)
|
---|
| 259 |
|
---|
[2939] | 260 | (defmethod r* ((poly1 poly) (poly2 poly))
|
---|
| 261 | "Non-destructively multiply POLY1 by POLY2."
|
---|
| 262 | (multiply-by (copy-instance POLY1) (copy-instance POLY2)))
|
---|
[2916] | 263 |
|
---|
[3044] | 264 | (defmethod left-tensor-product-by ((self poly) (other term))
|
---|
| 265 | (setf (poly-termlist self)
|
---|
| 266 | (mapcan #'(lambda (term)
|
---|
[3047] | 267 | (let ((prod (left-tensor-product-by term other)))
|
---|
[3044] | 268 | (cond
|
---|
| 269 | ((r-zerop prod) nil)
|
---|
| 270 | (t (list prod)))))
|
---|
[3048] | 271 | (poly-termlist self)))
|
---|
[3044] | 272 | self)
|
---|
| 273 |
|
---|
| 274 | (defmethod right-tensor-product-by ((self poly) (other term))
|
---|
[3045] | 275 | (setf (poly-termlist self)
|
---|
| 276 | (mapcan #'(lambda (term)
|
---|
[3046] | 277 | (let ((prod (right-tensor-product-by term other)))
|
---|
[3045] | 278 | (cond
|
---|
| 279 | ((r-zerop prod) nil)
|
---|
| 280 | (t (list prod)))))
|
---|
[3048] | 281 | (poly-termlist self)))
|
---|
[3045] | 282 | self)
|
---|
[3044] | 283 |
|
---|
[3062] | 284 | (defmethod left-tensor-product-by ((self poly) (other monom))
|
---|
| 285 | (setf (poly-termlist self)
|
---|
| 286 | (mapcan #'(lambda (term)
|
---|
| 287 | (let ((prod (left-tensor-product-by term other)))
|
---|
| 288 | (cond
|
---|
| 289 | ((r-zerop prod) nil)
|
---|
| 290 | (t (list prod)))))
|
---|
| 291 | (poly-termlist self)))
|
---|
| 292 | self)
|
---|
[3044] | 293 |
|
---|
[3062] | 294 | (defmethod right-tensor-product-by ((self poly) (other monom))
|
---|
| 295 | (setf (poly-termlist self)
|
---|
| 296 | (mapcan #'(lambda (term)
|
---|
| 297 | (let ((prod (right-tensor-product-by term other)))
|
---|
| 298 | (cond
|
---|
| 299 | ((r-zerop prod) nil)
|
---|
| 300 | (t (list prod)))))
|
---|
| 301 | (poly-termlist self)))
|
---|
| 302 | self)
|
---|
| 303 |
|
---|
| 304 |
|
---|
[3084] | 305 | (defun standard-extension (plist &aux (k (length plist)) (i 0))
|
---|
[2716] | 306 | "Calculate [U1*P1,U2*P2,...,UK*PK], where PLIST=[P1,P2,...,PK]
|
---|
[3060] | 307 | is a list of polynomials. Destructively modifies PLIST elements."
|
---|
[3061] | 308 | (mapc #'(lambda (poly)
|
---|
[3085] | 309 | (left-tensor-product-by
|
---|
| 310 | poly
|
---|
| 311 | (prog1
|
---|
| 312 | (make-monom-variable k i)
|
---|
| 313 | (incf i))))
|
---|
[3061] | 314 | plist))
|
---|
[52] | 315 |
|
---|
[3091] | 316 | (defmethod poly-dimension ((poly poly))
|
---|
| 317 | (cond ((r-zerop poly) -1)
|
---|
| 318 | (t (monom-dimension (leading-term poly)))))
|
---|
| 319 |
|
---|
[3087] | 320 | (defun standard-extension-1 (plist
|
---|
| 321 | &aux
|
---|
[3096] | 322 | (plist (standard-extension plist))
|
---|
[3087] | 323 | (nvars (poly-dimension (car plist))))
|
---|
[3081] | 324 | "Calculate [U1*P1-1,U2*P2-1,...,UK*PK-1], where PLIST=[P1,P2,...,PK].
|
---|
[3087] | 325 | Firstly, new K variables U1, U2, ..., UK, are inserted into each
|
---|
| 326 | polynomial. Subsequently, P1, P2, ..., PK are destructively modified
|
---|
[3105] | 327 | tantamount to replacing PI with UI*PI-1. It assumes that all
|
---|
| 328 | polynomials have the same dimension."
|
---|
[3089] | 329 | ;; Implementation note: we use STANDARD-EXTENSION and then subtract
|
---|
| 330 | ;; 1 from each polynomial; since UI*PI has no constant term,
|
---|
| 331 | ;; we just need to append the constant term at the end
|
---|
| 332 | ;; of each termlist.
|
---|
[3064] | 333 | (flet ((subtract-1 (p)
|
---|
[3104] | 334 | (append-item p (make-instance 'term :coeff -1 :dimension nvars))))
|
---|
[3083] | 335 | (setf plist (mapc #'subtract-1 plist)))
|
---|
[3077] | 336 | plist)
|
---|
[52] | 337 |
|
---|
[3096] | 338 | #|
|
---|
[52] | 339 |
|
---|
[3087] | 340 | (defun standard-sum (F plist
|
---|
[1475] | 341 | &aux
|
---|
| 342 | (k (length plist))
|
---|
[3079] | 343 | (d (+ k (monom-dimension (poly-lt (car plist)))))
|
---|
[1494] | 344 | ;; Add k variables to f
|
---|
[1493] | 345 | (f (poly-list-add-variables f k))
|
---|
[1495] | 346 | ;; Set PLIST to [U1*P1,U2*P2,...,UK*PK]
|
---|
[3077] | 347 | (plist (apply #'nconc (poly-standard-extension plist))))
|
---|
[3087] | 348 | "Calculate the polynomial U1*P1+U2*P2+...+UK*PK-1, where PLIST=[P1,P2,...,PK].
|
---|
| 349 | Firstly, new K variables, U1, U2, ..., UK, are inserted into each
|
---|
| 350 | polynomial. Subsequently, P1, P2, ..., PK are destructively modified
|
---|
| 351 | tantamount to replacing PI with UI*PI, and the resulting polynomials
|
---|
[3088] | 352 | are added. It should be noted that the term order is not modified,
|
---|
| 353 | which is equivalent to using a lexicographic order on the first K
|
---|
| 354 | variables."
|
---|
[1493] | 355 | (setf (cdr (last (poly-termlist plist)))
|
---|
[3087] | 356 | ;; Add -1 as the last term
|
---|
[1845] | 357 | (list (make-term :monom (make-monom :dimension d)
|
---|
| 358 | :coeff (funcall (ring-uminus ring) (funcall (ring-unit ring))))))
|
---|
[1493] | 359 | (append f (list plist)))
|
---|
[52] | 360 |
|
---|
[3076] | 361 |
|
---|
| 362 |
|
---|
[3096] | 363 |
|
---|
[1477] | 364 | (defun saturation-extension-1 (ring f p)
|
---|
[1497] | 365 | "Calculate [F, U*P-1]. It destructively modifies F."
|
---|
[1908] | 366 | (declare (type ring ring))
|
---|
[1477] | 367 | (polysaturation-extension ring f (list p)))
|
---|
[53] | 368 |
|
---|
| 369 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
| 370 | ;;
|
---|
| 371 | ;; Evaluation of polynomial (prefix) expressions
|
---|
| 372 | ;;
|
---|
| 373 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
| 374 |
|
---|
| 375 | (defun coerce-coeff (ring expr vars)
|
---|
| 376 | "Coerce an element of the coefficient ring to a constant polynomial."
|
---|
| 377 | ;; Modular arithmetic handler by rat
|
---|
[1908] | 378 | (declare (type ring ring))
|
---|
[1846] | 379 | (make-poly-from-termlist (list (make-term :monom (make-monom :dimension (length vars))
|
---|
| 380 | :coeff (funcall (ring-parse ring) expr)))
|
---|
[53] | 381 | 0))
|
---|
| 382 |
|
---|
[1046] | 383 | (defun poly-eval (expr vars
|
---|
| 384 | &optional
|
---|
[1668] | 385 | (ring +ring-of-integers+)
|
---|
[1048] | 386 | (order #'lex>)
|
---|
[1170] | 387 | (list-marker :[)
|
---|
[1047] | 388 | &aux
|
---|
| 389 | (ring-and-order (make-ring-and-order :ring ring :order order)))
|
---|
[1168] | 390 | "Evaluate Lisp form EXPR to a polynomial or a list of polynomials in
|
---|
[1208] | 391 | variables VARS. Return the resulting polynomial or list of
|
---|
| 392 | polynomials. Standard arithmetical operators in form EXPR are
|
---|
| 393 | replaced with their analogues in the ring of polynomials, and the
|
---|
| 394 | resulting expression is evaluated, resulting in a polynomial or a list
|
---|
[1209] | 395 | of polynomials in internal form. A similar operation in another computer
|
---|
| 396 | algebra system could be called 'expand' or so."
|
---|
[1909] | 397 | (declare (type ring ring))
|
---|
[1050] | 398 | (labels ((p-eval (arg) (poly-eval arg vars ring order))
|
---|
[1140] | 399 | (p-eval-scalar (arg) (poly-eval-scalar arg))
|
---|
[53] | 400 | (p-eval-list (args) (mapcar #'p-eval args))
|
---|
[989] | 401 | (p-add (x y) (poly-add ring-and-order x y)))
|
---|
[53] | 402 | (cond
|
---|
[1128] | 403 | ((null expr) (error "Empty expression"))
|
---|
[53] | 404 | ((eql expr 0) (make-poly-zero))
|
---|
| 405 | ((member expr vars :test #'equalp)
|
---|
| 406 | (let ((pos (position expr vars :test #'equalp)))
|
---|
[1657] | 407 | (make-poly-variable ring (length vars) pos)))
|
---|
[53] | 408 | ((atom expr)
|
---|
| 409 | (coerce-coeff ring expr vars))
|
---|
| 410 | ((eq (car expr) list-marker)
|
---|
| 411 | (cons list-marker (p-eval-list (cdr expr))))
|
---|
| 412 | (t
|
---|
| 413 | (case (car expr)
|
---|
| 414 | (+ (reduce #'p-add (p-eval-list (cdr expr))))
|
---|
| 415 | (- (case (length expr)
|
---|
| 416 | (1 (make-poly-zero))
|
---|
| 417 | (2 (poly-uminus ring (p-eval (cadr expr))))
|
---|
[989] | 418 | (3 (poly-sub ring-and-order (p-eval (cadr expr)) (p-eval (caddr expr))))
|
---|
| 419 | (otherwise (poly-sub ring-and-order (p-eval (cadr expr))
|
---|
[53] | 420 | (reduce #'p-add (p-eval-list (cddr expr)))))))
|
---|
| 421 | (*
|
---|
| 422 | (if (endp (cddr expr)) ;unary
|
---|
| 423 | (p-eval (cdr expr))
|
---|
[989] | 424 | (reduce #'(lambda (p q) (poly-mul ring-and-order p q)) (p-eval-list (cdr expr)))))
|
---|
[1106] | 425 | (/
|
---|
| 426 | ;; A polynomial can be divided by a scalar
|
---|
[1115] | 427 | (cond
|
---|
| 428 | ((endp (cddr expr))
|
---|
[1117] | 429 | ;; A special case (/ ?), the inverse
|
---|
[1119] | 430 | (coerce-coeff ring (apply (ring-div ring) (cdr expr)) vars))
|
---|
[1128] | 431 | (t
|
---|
[1115] | 432 | (let ((num (p-eval (cadr expr)))
|
---|
[1142] | 433 | (denom-inverse (apply (ring-div ring)
|
---|
| 434 | (cons (funcall (ring-unit ring))
|
---|
| 435 | (mapcar #'p-eval-scalar (cddr expr))))))
|
---|
[1118] | 436 | (scalar-times-poly ring denom-inverse num)))))
|
---|
[53] | 437 | (expt
|
---|
| 438 | (cond
|
---|
| 439 | ((member (cadr expr) vars :test #'equalp)
|
---|
| 440 | ;;Special handling of (expt var pow)
|
---|
| 441 | (let ((pos (position (cadr expr) vars :test #'equalp)))
|
---|
[1657] | 442 | (make-poly-variable ring (length vars) pos (caddr expr))))
|
---|
[53] | 443 | ((not (and (integerp (caddr expr)) (plusp (caddr expr))))
|
---|
| 444 | ;; Negative power means division in coefficient ring
|
---|
| 445 | ;; Non-integer power means non-polynomial coefficient
|
---|
| 446 | (coerce-coeff ring expr vars))
|
---|
[989] | 447 | (t (poly-expt ring-and-order (p-eval (cadr expr)) (caddr expr)))))
|
---|
[53] | 448 | (otherwise
|
---|
| 449 | (coerce-coeff ring expr vars)))))))
|
---|
| 450 |
|
---|
[1133] | 451 | (defun poly-eval-scalar (expr
|
---|
| 452 | &optional
|
---|
[1668] | 453 | (ring +ring-of-integers+)
|
---|
[1133] | 454 | &aux
|
---|
| 455 | (order #'lex>))
|
---|
| 456 | "Evaluate a scalar expression EXPR in ring RING."
|
---|
[1910] | 457 | (declare (type ring ring))
|
---|
[1133] | 458 | (poly-lc (poly-eval expr nil ring order)))
|
---|
| 459 |
|
---|
[1189] | 460 | (defun spoly (ring-and-order f g
|
---|
| 461 | &aux
|
---|
| 462 | (ring (ro-ring ring-and-order)))
|
---|
[55] | 463 | "It yields the S-polynomial of polynomials F and G."
|
---|
[1911] | 464 | (declare (type ring-and-order ring-and-order) (type poly f g))
|
---|
[55] | 465 | (let* ((lcm (monom-lcm (poly-lm f) (poly-lm g)))
|
---|
[2913] | 466 | (mf (monom-div lcm (poly-lm f)))
|
---|
| 467 | (mg (monom-div lcm (poly-lm g))))
|
---|
[55] | 468 | (declare (type monom mf mg))
|
---|
| 469 | (multiple-value-bind (c cf cg)
|
---|
| 470 | (funcall (ring-ezgcd ring) (poly-lc f) (poly-lc g))
|
---|
| 471 | (declare (ignore c))
|
---|
| 472 | (poly-sub
|
---|
[1189] | 473 | ring-and-order
|
---|
[55] | 474 | (scalar-times-poly ring cg (monom-times-poly mf f))
|
---|
| 475 | (scalar-times-poly ring cf (monom-times-poly mg g))))))
|
---|
[53] | 476 |
|
---|
| 477 |
|
---|
[55] | 478 | (defun poly-primitive-part (ring p)
|
---|
| 479 | "Divide polynomial P with integer coefficients by gcd of its
|
---|
| 480 | coefficients and return the result."
|
---|
[1912] | 481 | (declare (type ring ring) (type poly p))
|
---|
[55] | 482 | (if (poly-zerop p)
|
---|
| 483 | (values p 1)
|
---|
[2913] | 484 | (let ((c (poly-content ring p)))
|
---|
| 485 | (values (make-poly-from-termlist
|
---|
| 486 | (mapcar
|
---|
| 487 | #'(lambda (x)
|
---|
| 488 | (make-term :monom (term-monom x)
|
---|
| 489 | :coeff (funcall (ring-div ring) (term-coeff x) c)))
|
---|
| 490 | (poly-termlist p))
|
---|
| 491 | (poly-sugar p))
|
---|
| 492 | c))))
|
---|
[55] | 493 |
|
---|
| 494 | (defun poly-content (ring p)
|
---|
| 495 | "Greatest common divisor of the coefficients of the polynomial P. Use the RING structure
|
---|
| 496 | to compute the greatest common divisor."
|
---|
[1913] | 497 | (declare (type ring ring) (type poly p))
|
---|
[55] | 498 | (reduce (ring-gcd ring) (mapcar #'term-coeff (rest (poly-termlist p))) :initial-value (poly-lc p)))
|
---|
[1066] | 499 |
|
---|
[1091] | 500 | (defun read-infix-form (&key (stream t))
|
---|
[1066] | 501 | "Parser of infix expressions with integer/rational coefficients
|
---|
| 502 | The parser will recognize two kinds of polynomial expressions:
|
---|
| 503 |
|
---|
| 504 | - polynomials in fully expanded forms with coefficients
|
---|
| 505 | written in front of symbolic expressions; constants can be optionally
|
---|
| 506 | enclosed in (); for example, the infix form
|
---|
| 507 | X^2-Y^2+(-4/3)*U^2*W^3-5
|
---|
| 508 | parses to
|
---|
| 509 | (+ (- (EXPT X 2) (EXPT Y 2)) (* (- (/ 4 3)) (EXPT U 2) (EXPT W 3)) (- 5))
|
---|
| 510 |
|
---|
| 511 | - lists of polynomials; for example
|
---|
| 512 | [X-Y, X^2+3*Z]
|
---|
| 513 | parses to
|
---|
| 514 | (:[ (- X Y) (+ (EXPT X 2) (* 3 Z)))
|
---|
| 515 | where the first symbol [ marks a list of polynomials.
|
---|
| 516 |
|
---|
| 517 | -other infix expressions, for example
|
---|
| 518 | [(X-Y)*(X+Y)/Z,(X+1)^2]
|
---|
| 519 | parses to:
|
---|
| 520 | (:[ (/ (* (- X Y) (+ X Y)) Z) (EXPT (+ X 1) 2))
|
---|
| 521 | Currently this function is implemented using M. Kantrowitz's INFIX package."
|
---|
| 522 | (read-from-string
|
---|
| 523 | (concatenate 'string
|
---|
[2913] | 524 | "#I("
|
---|
| 525 | (with-output-to-string (s)
|
---|
| 526 | (loop
|
---|
| 527 | (multiple-value-bind (line eof)
|
---|
| 528 | (read-line stream t)
|
---|
| 529 | (format s "~A" line)
|
---|
| 530 | (when eof (return)))))
|
---|
| 531 | ")")))
|
---|
| 532 |
|
---|
[1145] | 533 | (defun read-poly (vars &key
|
---|
| 534 | (stream t)
|
---|
[1668] | 535 | (ring +ring-of-integers+)
|
---|
[1145] | 536 | (order #'lex>))
|
---|
[1067] | 537 | "Reads an expression in prefix form from a stream STREAM.
|
---|
[1144] | 538 | The expression read from the strem should represent a polynomial or a
|
---|
| 539 | list of polynomials in variables VARS, over the ring RING. The
|
---|
| 540 | polynomial or list of polynomials is returned, with terms in each
|
---|
| 541 | polynomial ordered according to monomial order ORDER."
|
---|
[1146] | 542 | (poly-eval (read-infix-form :stream stream) vars ring order))
|
---|
[1092] | 543 |
|
---|
[1146] | 544 | (defun string->poly (str vars
|
---|
[1164] | 545 | &optional
|
---|
[1668] | 546 | (ring +ring-of-integers+)
|
---|
[1146] | 547 | (order #'lex>))
|
---|
| 548 | "Converts a string STR to a polynomial in variables VARS."
|
---|
[1097] | 549 | (with-input-from-string (s str)
|
---|
[1165] | 550 | (read-poly vars :stream s :ring ring :order order)))
|
---|
[1095] | 551 |
|
---|
[1143] | 552 | (defun poly->alist (p)
|
---|
| 553 | "Convert a polynomial P to an association list. Thus, the format of the
|
---|
| 554 | returned value is ((MONOM[0] . COEFF[0]) (MONOM[1] . COEFF[1]) ...), where
|
---|
| 555 | MONOM[I] is a list of exponents in the monomial and COEFF[I] is the
|
---|
| 556 | corresponding coefficient in the ring."
|
---|
[1171] | 557 | (cond
|
---|
| 558 | ((poly-p p)
|
---|
| 559 | (mapcar #'term->cons (poly-termlist p)))
|
---|
| 560 | ((and (consp p) (eq (car p) :[))
|
---|
[1172] | 561 | (cons :[ (mapcar #'poly->alist (cdr p))))))
|
---|
[1143] | 562 |
|
---|
[1164] | 563 | (defun string->alist (str vars
|
---|
[2913] | 564 | &optional
|
---|
| 565 | (ring +ring-of-integers+)
|
---|
| 566 | (order #'lex>))
|
---|
[1143] | 567 | "Convert a string STR representing a polynomial or polynomial list to
|
---|
[1158] | 568 | an association list (... (MONOM . COEFF) ...)."
|
---|
[1166] | 569 | (poly->alist (string->poly str vars ring order)))
|
---|
[1440] | 570 |
|
---|
| 571 | (defun poly-equal-no-sugar-p (p q)
|
---|
| 572 | "Compare polynomials for equality, ignoring sugar."
|
---|
[1914] | 573 | (declare (type poly p q))
|
---|
[1440] | 574 | (equalp (poly-termlist p) (poly-termlist q)))
|
---|
[1559] | 575 |
|
---|
| 576 | (defun poly-set-equal-no-sugar-p (p q)
|
---|
| 577 | "Compare polynomial sets P and Q for equality, ignoring sugar."
|
---|
| 578 | (null (set-exclusive-or p q :test #'poly-equal-no-sugar-p )))
|
---|
[1560] | 579 |
|
---|
| 580 | (defun poly-list-equal-no-sugar-p (p q)
|
---|
| 581 | "Compare polynomial lists P and Q for equality, ignoring sugar."
|
---|
| 582 | (every #'poly-equal-no-sugar-p p q))
|
---|
[2456] | 583 | |#
|
---|