close Warning: Can't synchronize with repository "(default)" (The repository directory has changed, you should resynchronize the repository with: trac-admin $ENV repository resync '(default)'). Look in the Trac log for more information.

source: branches/f4grobner/pol.lisp@ 2005

Last change on this file since 2005 was 2005, checked in by Marek Rychlik, 9 years ago

* empty log message *

File size: 6.2 KB
Line 
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;;
24;; Polynomials implemented in CLOS
25;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
26;;
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;;
41;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
42
43(defpackage "POL"
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
97(in-package :pol)
98
99(proclaim '(optimize (speed 0) (space 0) (safety 3) (debug 3)))
100
101(defclass poly ()
102 ((expr :initarg :expr :accessor expr))
103 ((vars :initarg :vars :accessor vars))
104 (:default-initargs :expr 0 :vars nil))
105
106(defmethod print-object ((self poly) stream)
107 (princ (slot-value self 'expr)))
108
109(defmethod poly-add ((p poly) (q poly)))
110
111(defmethod poly-sub ((p poly) (q poly)))
112
113(defmethod poly-uminus ((self poly)))
114
115(defmethod poly-mul ((p poly) (q poly)))
116
117(defmethod poly-expt ((self poly) n))
118
119(defun poly-eval (expr vars))
120 "Evaluate Lisp form EXPR to a polynomial or a list of polynomials in
121variables VARS. Return the resulting polynomial or list of
122polynomials. Standard arithmetical operators in form EXPR are
123replaced with their analogues in the ring of polynomials, and the
124resulting expression is evaluated, resulting in a polynomial or a list
125of polynomials in internal form. A similar operation in another computer
126algebra system could be called 'expand' or so."
127 (cond
128 ((null expr)
129 (error "Empty expression"))
130 ((eql expr 0)
131 ())
132 ((member expr vars :test #'equalp)
133 (let ((pos (position expr vars :test #'equalp)))
134 (make-poly-variable ring (length vars) pos)))
135 ((atom expr)
136 (scalar->poly ring expr vars))
137 ((eq (car expr) list-marker)
138 (cons list-marker (p-eval-list (cdr expr))))
139 (t
140 (case (car expr)
141 (+ (reduce #'p-add (p-eval-list (cdr expr))))
142 (- (case (length expr)
143 (1 (make-poly-zero))
144 (2 (poly-uminus ring (p-eval (cadr expr))))
145 (3 (poly-sub ring-and-order (p-eval (cadr expr)) (p-eval (caddr expr))))
146 (otherwise (poly-sub ring-and-order (p-eval (cadr expr))
147 (reduce #'p-add (p-eval-list (cddr expr)))))))
148 (*
149 (if (endp (cddr expr)) ;unary
150 (p-eval (cdr expr))
151 (reduce #'(lambda (p q) (poly-mul ring-and-order p q)) (p-eval-list (cdr expr)))))
152 (/
153 ;; A polynomial can be divided by a scalar
154 (cond
155 ((endp (cddr expr))
156 ;; A special case (/ ?), the inverse
157 (scalar->poly ring (apply (ring-div ring) (cdr expr)) vars))
158 (t
159 (let ((num (p-eval (cadr expr)))
160 (denom-inverse (apply (ring-div ring)
161 (cons (funcall (ring-unit ring))
162 (mapcar #'p-eval-scalar (cddr expr))))))
163 (scalar-times-poly ring denom-inverse num)))))
164 (expt
165 (cond
166 ((member (cadr expr) vars :test #'equalp)
167 ;;Special handling of (expt var pow)
168 (let ((pos (position (cadr expr) vars :test #'equalp)))
169 (make-poly-variable ring (length vars) pos (caddr expr))))
170 ((not (and (integerp (caddr expr)) (plusp (caddr expr))))
171 ;; Negative power means division in coefficient ring
172 ;; Non-integer power means non-polynomial coefficient
173 (scalar->poly ring expr vars))
174 (t (poly-expt ring-and-order (p-eval (cadr expr)) (caddr expr)))))
175 (otherwise
176 (scalar->poly ring expr vars))))))
177
178(defun poly-eval-scalar (expr
179 &optional
180 (ring +ring-of-integers+)
181 &aux
182 (order #'lex>))
183 "Evaluate a scalar expression EXPR in ring RING."
184 (declare (type ring ring))
185 (poly-lc (poly-eval expr nil ring order)))
186
187
188
189
Note: See TracBrowser for help on using the repository browser.