[1969] | 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 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
| 23 | ;;
|
---|
[1997] | 24 | ;; Polynomials implemented in CLOS
|
---|
[1998] | 25 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
[1969] | 26 | ;;
|
---|
[1998] | 27 | ;; A polynomial is an collection of terms. A
|
---|
| 28 | ;; term has a monomial and a coefficient.
|
---|
| 29 | ;;
|
---|
| 30 | ;; A polynomial can be represented by an s-expp
|
---|
| 31 | ;; (EXPR . VARS) where EXPR is an arithmetical formula
|
---|
| 32 | ;; recursively built of the arithmetical operations,
|
---|
| 33 | ;; and VARS are the variables of the polynomial.
|
---|
| 34 | ;; If a subtree of this s-exp is not an arithmetical
|
---|
| 35 | ;; operator +, -, *, expt, and is not a member
|
---|
| 36 | ;; of VARS then it represents a scalar expression
|
---|
| 37 | ;; which the Lisp reader must know how to convert
|
---|
| 38 | ;; into an object for which can be multiplied by a variable,
|
---|
| 39 | ;; subject to commutativity and associativity rules.
|
---|
| 40 | ;;
|
---|
[1969] | 41 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
---|
| 42 |
|
---|
[1970] | 43 | (defpackage "POL"
|
---|
[1969] | 44 | (:use :cl :ring :ring-and-order :monom :order :term :termlist :infix)
|
---|
| 45 | (:export "POLY"
|
---|
| 46 | "POLY-TERMLIST"
|
---|
| 47 | "POLY-SUGAR"
|
---|
| 48 | "POLY-RESET-SUGAR"
|
---|
| 49 | "POLY-LT"
|
---|
| 50 | "MAKE-POLY-FROM-TERMLIST"
|
---|
| 51 | "MAKE-POLY-ZERO"
|
---|
| 52 | "MAKE-POLY-VARIABLE"
|
---|
| 53 | "POLY-UNIT"
|
---|
| 54 | "POLY-LM"
|
---|
| 55 | "POLY-SECOND-LM"
|
---|
| 56 | "POLY-SECOND-LT"
|
---|
| 57 | "POLY-LC"
|
---|
| 58 | "POLY-SECOND-LC"
|
---|
| 59 | "POLY-ZEROP"
|
---|
| 60 | "POLY-LENGTH"
|
---|
| 61 | "SCALAR-TIMES-POLY"
|
---|
| 62 | "SCALAR-TIMES-POLY-1"
|
---|
| 63 | "MONOM-TIMES-POLY"
|
---|
| 64 | "TERM-TIMES-POLY"
|
---|
| 65 | "POLY-ADD"
|
---|
| 66 | "POLY-SUB"
|
---|
| 67 | "POLY-UMINUS"
|
---|
| 68 | "POLY-MUL"
|
---|
| 69 | "POLY-EXPT"
|
---|
| 70 | "POLY-APPEND"
|
---|
| 71 | "POLY-NREVERSE"
|
---|
| 72 | "POLY-REVERSE"
|
---|
| 73 | "POLY-CONTRACT"
|
---|
| 74 | "POLY-EXTEND"
|
---|
| 75 | "POLY-ADD-VARIABLES"
|
---|
| 76 | "POLY-LIST-ADD-VARIABLES"
|
---|
| 77 | "POLY-STANDARD-EXTENSION"
|
---|
| 78 | "SATURATION-EXTENSION"
|
---|
| 79 | "POLYSATURATION-EXTENSION"
|
---|
| 80 | "SATURATION-EXTENSION-1"
|
---|
| 81 | "COERCE-COEFF"
|
---|
| 82 | "POLY-EVAL"
|
---|
| 83 | "POLY-EVAL-SCALAR"
|
---|
| 84 | "SPOLY"
|
---|
| 85 | "POLY-PRIMITIVE-PART"
|
---|
| 86 | "POLY-CONTENT"
|
---|
| 87 | "READ-INFIX-FORM"
|
---|
| 88 | "READ-POLY"
|
---|
| 89 | "STRING->POLY"
|
---|
| 90 | "POLY->ALIST"
|
---|
| 91 | "STRING->ALIST"
|
---|
| 92 | "POLY-EQUAL-NO-SUGAR-P"
|
---|
| 93 | "POLY-SET-EQUAL-NO-SUGAR-P"
|
---|
| 94 | "POLY-LIST-EQUAL-NO-SUGAR-P"
|
---|
| 95 | ))
|
---|
| 96 |
|
---|
[1972] | 97 | (in-package :pol)
|
---|
[1969] | 98 |
|
---|
[1973] | 99 | (proclaim '(optimize (speed 0) (space 0) (safety 3) (debug 3)))
|
---|
[1969] | 100 |
|
---|
[1979] | 101 | (defclass poly ()
|
---|
[2000] | 102 | ((termlist :initarg :termlist :accessor termlist))
|
---|
| 103 | (:default-initargs :termlist nil))
|
---|
[1969] | 104 |
|
---|
[2000] | 105 | (defun make-poly-from-termlist (termlist)
|
---|
[1982] | 106 | (make-instance 'poly :termlist termlist :sugar sugar))
|
---|
| 107 |
|
---|
| 108 | (defun make-poly-zero (&aux (termlist nil) (sugar -1))
|
---|
[2000] | 109 | (make-instance 'poly :termlist termlist))
|
---|
[1982] | 110 |
|
---|
[1974] | 111 | (defun make-poly-variable (ring nvars pos &optional (power 1)
|
---|
| 112 | &aux
|
---|
| 113 | (termlist (list
|
---|
[2000] | 114 | (make-term-variable ring nvars pos power))))
|
---|
| 115 | (make-instance 'poly :termlist termlist))
|
---|
[1974] | 116 |
|
---|
| 117 | (defun poly-unit (ring dimension
|
---|
| 118 | &aux
|
---|
[2000] | 119 | (termlist (termlist-unit ring dimension)))
|
---|
| 120 | (make-instance 'poly :termlist termlist))
|
---|
[1974] | 121 |
|
---|
| 122 |
|
---|
[2000] | 123 | (defmethod print-object ((self poly) stream)
|
---|
| 124 | (princ (slot-value self 'termlist)))
|
---|
[1986] | 125 |
|
---|
[2000] | 126 | (defmethod poly-termlist ((self poly))
|
---|
| 127 | (slot-value self 'termlist))
|
---|
[1987] | 128 |
|
---|
[1988] | 129 | (defmethod (setf poly-termlist) (new-value (poly poly))
|
---|
| 130 | (setf (slot-value poly 'termlist) new-value))
|
---|
[1987] | 131 |
|
---|
[2000] | 132 | (defmethod poly-add ((p poly) (q poly)))
|
---|
[1969] | 133 |
|
---|
[2000] | 134 | (defmethod poly-sub ((p poly) (q poly)))
|
---|
[1969] | 135 |
|
---|
[2000] | 136 | (defmethod poly-uminus ((self poly)))
|
---|
[1969] | 137 |
|
---|
[2000] | 138 | (defmethod poly-mul ((p poly) (poly q)))
|
---|
[1969] | 139 |
|
---|
[2000] | 140 | (defmethod poly-expt ((self poly) n))
|
---|
[1969] | 141 |
|
---|
[2000] | 142 | (defmethod initialize-instance :after ((self poly) expr vars)
|
---|
[1969] | 143 | "Evaluate Lisp form EXPR to a polynomial or a list of polynomials in
|
---|
| 144 | variables VARS. Return the resulting polynomial or list of
|
---|
| 145 | polynomials. Standard arithmetical operators in form EXPR are
|
---|
| 146 | replaced with their analogues in the ring of polynomials, and the
|
---|
| 147 | resulting expression is evaluated, resulting in a polynomial or a list
|
---|
| 148 | of polynomials in internal form. A similar operation in another computer
|
---|
| 149 | algebra system could be called 'expand' or so."
|
---|
| 150 | (cond
|
---|
| 151 | ((null expr) (error "Empty expression"))
|
---|
| 152 | ((eql expr 0) (make-poly-zero))
|
---|
| 153 | ((member expr vars :test #'equalp)
|
---|
| 154 | (let ((pos (position expr vars :test #'equalp)))
|
---|
| 155 | (make-poly-variable ring (length vars) pos)))
|
---|
| 156 | ((atom expr)
|
---|
[1996] | 157 | (scalar->poly ring expr vars))
|
---|
[1969] | 158 | ((eq (car expr) list-marker)
|
---|
| 159 | (cons list-marker (p-eval-list (cdr expr))))
|
---|
| 160 | (t
|
---|
| 161 | (case (car expr)
|
---|
| 162 | (+ (reduce #'p-add (p-eval-list (cdr expr))))
|
---|
| 163 | (- (case (length expr)
|
---|
| 164 | (1 (make-poly-zero))
|
---|
| 165 | (2 (poly-uminus ring (p-eval (cadr expr))))
|
---|
| 166 | (3 (poly-sub ring-and-order (p-eval (cadr expr)) (p-eval (caddr expr))))
|
---|
| 167 | (otherwise (poly-sub ring-and-order (p-eval (cadr expr))
|
---|
| 168 | (reduce #'p-add (p-eval-list (cddr expr)))))))
|
---|
| 169 | (*
|
---|
| 170 | (if (endp (cddr expr)) ;unary
|
---|
| 171 | (p-eval (cdr expr))
|
---|
| 172 | (reduce #'(lambda (p q) (poly-mul ring-and-order p q)) (p-eval-list (cdr expr)))))
|
---|
| 173 | (/
|
---|
| 174 | ;; A polynomial can be divided by a scalar
|
---|
| 175 | (cond
|
---|
| 176 | ((endp (cddr expr))
|
---|
| 177 | ;; A special case (/ ?), the inverse
|
---|
[1996] | 178 | (scalar->poly ring (apply (ring-div ring) (cdr expr)) vars))
|
---|
[1969] | 179 | (t
|
---|
| 180 | (let ((num (p-eval (cadr expr)))
|
---|
| 181 | (denom-inverse (apply (ring-div ring)
|
---|
| 182 | (cons (funcall (ring-unit ring))
|
---|
| 183 | (mapcar #'p-eval-scalar (cddr expr))))))
|
---|
| 184 | (scalar-times-poly ring denom-inverse num)))))
|
---|
| 185 | (expt
|
---|
| 186 | (cond
|
---|
| 187 | ((member (cadr expr) vars :test #'equalp)
|
---|
| 188 | ;;Special handling of (expt var pow)
|
---|
| 189 | (let ((pos (position (cadr expr) vars :test #'equalp)))
|
---|
| 190 | (make-poly-variable ring (length vars) pos (caddr expr))))
|
---|
| 191 | ((not (and (integerp (caddr expr)) (plusp (caddr expr))))
|
---|
| 192 | ;; Negative power means division in coefficient ring
|
---|
| 193 | ;; Non-integer power means non-polynomial coefficient
|
---|
[1996] | 194 | (scalar->poly ring expr vars))
|
---|
[1969] | 195 | (t (poly-expt ring-and-order (p-eval (cadr expr)) (caddr expr)))))
|
---|
| 196 | (otherwise
|
---|
[1996] | 197 | (scalar->poly ring expr vars)))))))
|
---|
[1969] | 198 |
|
---|
| 199 | (defun poly-eval-scalar (expr
|
---|
| 200 | &optional
|
---|
| 201 | (ring +ring-of-integers+)
|
---|
| 202 | &aux
|
---|
| 203 | (order #'lex>))
|
---|
| 204 | "Evaluate a scalar expression EXPR in ring RING."
|
---|
| 205 | (declare (type ring ring))
|
---|
| 206 | (poly-lc (poly-eval expr nil ring order)))
|
---|
| 207 |
|
---|
| 208 | (defun spoly (ring-and-order f g
|
---|
| 209 | &aux
|
---|
| 210 | (ring (ro-ring ring-and-order)))
|
---|
| 211 | "It yields the S-polynomial of polynomials F and G."
|
---|
| 212 | (declare (type ring-and-order ring-and-order) (type poly f g))
|
---|
| 213 | (let* ((lcm (monom-lcm (poly-lm f) (poly-lm g)))
|
---|
| 214 | (mf (monom-div lcm (poly-lm f)))
|
---|
| 215 | (mg (monom-div lcm (poly-lm g))))
|
---|
| 216 | (declare (type monom mf mg))
|
---|
| 217 | (multiple-value-bind (c cf cg)
|
---|
| 218 | (funcall (ring-ezgcd ring) (poly-lc f) (poly-lc g))
|
---|
| 219 | (declare (ignore c))
|
---|
| 220 | (poly-sub
|
---|
| 221 | ring-and-order
|
---|
| 222 | (scalar-times-poly ring cg (monom-times-poly mf f))
|
---|
| 223 | (scalar-times-poly ring cf (monom-times-poly mg g))))))
|
---|
| 224 |
|
---|
| 225 |
|
---|
| 226 | (defun poly-primitive-part (ring p)
|
---|
| 227 | "Divide polynomial P with integer coefficients by gcd of its
|
---|
| 228 | coefficients and return the result."
|
---|
| 229 | (declare (type ring ring) (type poly p))
|
---|
| 230 | (if (poly-zerop p)
|
---|
| 231 | (values p 1)
|
---|
| 232 | (let ((c (poly-content ring p)))
|
---|
| 233 | (values (make-poly-from-termlist
|
---|
| 234 | (mapcar
|
---|
| 235 | #'(lambda (x)
|
---|
| 236 | (make-term :monom (term-monom x)
|
---|
| 237 | :coeff (funcall (ring-div ring) (term-coeff x) c)))
|
---|
| 238 | (poly-termlist p))
|
---|
| 239 | (poly-sugar p))
|
---|
| 240 | c))))
|
---|
| 241 |
|
---|
| 242 | (defun poly-content (ring p)
|
---|
| 243 | "Greatest common divisor of the coefficients of the polynomial P. Use the RING structure
|
---|
| 244 | to compute the greatest common divisor."
|
---|
| 245 | (declare (type ring ring) (type poly p))
|
---|
| 246 | (reduce (ring-gcd ring) (mapcar #'term-coeff (rest (poly-termlist p))) :initial-value (poly-lc p)))
|
---|
| 247 |
|
---|
| 248 | (defun read-infix-form (&key (stream t))
|
---|
| 249 | "Parser of infix expressions with integer/rational coefficients
|
---|
| 250 | The parser will recognize two kinds of polynomial expressions:
|
---|
| 251 |
|
---|
| 252 | - polynomials in fully expanded forms with coefficients
|
---|
| 253 | written in front of symbolic expressions; constants can be optionally
|
---|
| 254 | enclosed in (); for example, the infix form
|
---|
| 255 | X^2-Y^2+(-4/3)*U^2*W^3-5
|
---|
| 256 | parses to
|
---|
| 257 | (+ (- (EXPT X 2) (EXPT Y 2)) (* (- (/ 4 3)) (EXPT U 2) (EXPT W 3)) (- 5))
|
---|
| 258 |
|
---|
| 259 | - lists of polynomials; for example
|
---|
| 260 | [X-Y, X^2+3*Z]
|
---|
| 261 | parses to
|
---|
| 262 | (:[ (- X Y) (+ (EXPT X 2) (* 3 Z)))
|
---|
| 263 | where the first symbol [ marks a list of polynomials.
|
---|
| 264 |
|
---|
| 265 | -other infix expressions, for example
|
---|
| 266 | [(X-Y)*(X+Y)/Z,(X+1)^2]
|
---|
| 267 | parses to:
|
---|
| 268 | (:[ (/ (* (- X Y) (+ X Y)) Z) (EXPT (+ X 1) 2))
|
---|
| 269 | Currently this function is implemented using M. Kantrowitz's INFIX package."
|
---|
| 270 | (read-from-string
|
---|
| 271 | (concatenate 'string
|
---|
| 272 | "#I("
|
---|
| 273 | (with-output-to-string (s)
|
---|
| 274 | (loop
|
---|
| 275 | (multiple-value-bind (line eof)
|
---|
| 276 | (read-line stream t)
|
---|
| 277 | (format s "~A" line)
|
---|
| 278 | (when eof (return)))))
|
---|
| 279 | ")")))
|
---|
| 280 |
|
---|
| 281 | (defun read-poly (vars &key
|
---|
| 282 | (stream t)
|
---|
| 283 | (ring +ring-of-integers+)
|
---|
| 284 | (order #'lex>))
|
---|
| 285 | "Reads an expression in prefix form from a stream STREAM.
|
---|
| 286 | The expression read from the strem should represent a polynomial or a
|
---|
| 287 | list of polynomials in variables VARS, over the ring RING. The
|
---|
| 288 | polynomial or list of polynomials is returned, with terms in each
|
---|
| 289 | polynomial ordered according to monomial order ORDER."
|
---|
| 290 | (poly-eval (read-infix-form :stream stream) vars ring order))
|
---|
| 291 |
|
---|
| 292 | (defun string->poly (str vars
|
---|
| 293 | &optional
|
---|
| 294 | (ring +ring-of-integers+)
|
---|
| 295 | (order #'lex>))
|
---|
| 296 | "Converts a string STR to a polynomial in variables VARS."
|
---|
| 297 | (with-input-from-string (s str)
|
---|
| 298 | (read-poly vars :stream s :ring ring :order order)))
|
---|
| 299 |
|
---|
| 300 | (defun poly->alist (p)
|
---|
| 301 | "Convert a polynomial P to an association list. Thus, the format of the
|
---|
| 302 | returned value is ((MONOM[0] . COEFF[0]) (MONOM[1] . COEFF[1]) ...), where
|
---|
| 303 | MONOM[I] is a list of exponents in the monomial and COEFF[I] is the
|
---|
| 304 | corresponding coefficient in the ring."
|
---|
| 305 | (cond
|
---|
| 306 | ((poly-p p)
|
---|
| 307 | (mapcar #'term->cons (poly-termlist p)))
|
---|
| 308 | ((and (consp p) (eq (car p) :[))
|
---|
| 309 | (cons :[ (mapcar #'poly->alist (cdr p))))))
|
---|
| 310 |
|
---|
| 311 | (defun string->alist (str vars
|
---|
| 312 | &optional
|
---|
| 313 | (ring +ring-of-integers+)
|
---|
| 314 | (order #'lex>))
|
---|
| 315 | "Convert a string STR representing a polynomial or polynomial list to
|
---|
| 316 | an association list (... (MONOM . COEFF) ...)."
|
---|
| 317 | (poly->alist (string->poly str vars ring order)))
|
---|
| 318 |
|
---|
| 319 | (defun poly-equal-no-sugar-p (p q)
|
---|
| 320 | "Compare polynomials for equality, ignoring sugar."
|
---|
| 321 | (declare (type poly p q))
|
---|
| 322 | (equalp (poly-termlist p) (poly-termlist q)))
|
---|
| 323 |
|
---|
| 324 | (defun poly-set-equal-no-sugar-p (p q)
|
---|
| 325 | "Compare polynomial sets P and Q for equality, ignoring sugar."
|
---|
| 326 | (null (set-exclusive-or p q :test #'poly-equal-no-sugar-p )))
|
---|
| 327 |
|
---|
| 328 | (defun poly-list-equal-no-sugar-p (p q)
|
---|
| 329 | "Compare polynomial lists P and Q for equality, ignoring sugar."
|
---|
| 330 | (every #'poly-equal-no-sugar-p p q))
|
---|
[1993] | 331 |
|
---|
| 332 |
|
---|