1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancements. ;;;;;
5 ;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
6 ;;; All rights reserved ;;;;;
7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
8 ;;; (c) Copyright 1980 Massachusetts Institute of Technology ;;;
9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
13 (macsyma-module ratmac macro
)
18 ;; A polynomial in a single variable is stored as a sparse list of coefficients,
19 ;; whose first element is the polynomial's variable. The rest of the elements
20 ;; form a plist with keys the powers (in descending order) and values the
23 ;; For example, 42*x^2 + 1 might be stored as ($x 2 42 0 1). If, say,
24 ;; x*sin(x)+x^2 is respresented as a polynomial in x, we might expect it to come
25 ;; out as something like
27 ;; ($x 2 1 1 ((%sin) $x)),
29 ;; but to make it easier to work with polynomials we don't allow arbitrary
30 ;; conses as coefficients. What actually happens is that the expression is
31 ;; thought of as a polynomial in two variables x and "sin(x)". More on that
34 ;; Multivariate polynomials are stored in basically the same way as single
35 ;; variable polynomials, using the observation that a polynomial in X and Y with
36 ;; coefficients in K is the same as a polynomial in X with coefficients in K[Y].
38 ;; Specifically, the coefficient terms can be polynomials themselves (in other
39 ;; variables). So x^2 + x*y could be rperesented as (($x 2 1 1 ($y 1 1))) or
40 ;; alternatively as (($y 1 ($x 1 1) 0 ($x 2 1))), depending on whether x or y
41 ;; was taken as the primary variable. If you add together two polynomials in
42 ;; different variables (say x+1 and y+2) in the rational function code, then it
43 ;; decides on the main variable using POINTERGP. This only works if the
44 ;; variables have already been given a numbering by the rest of the rational
47 ;; In the x*sin(x) + x^2 example above, the expression can be represented as
48 ;; something like ($x 2 1 1 (sinx 1 1)). When passed around as expressions
49 ;; outside of the core rational function code, polynomials come with some header
50 ;; information that explains what the variables are. In this case, it would be
51 ;; responsible for remembering that "sinx" means sin(x).
53 ;; As a slightly special case, a polynomial can also be an atom, in which case
54 ;; it is treated as a degree zero polynomial in no particular variable. Test for
55 ;; this using the pcoefp macro defined below.
57 ;; There are accessor macros for the parts of a polynomial defined below: p-var,
58 ;; p-terms, p-lc, p-le and p-red (which extract the primary variable, the list
59 ;; of powers and coefficients, the leading coefficient, the leading exponent and
60 ;; the list of powers and coefficients except the leading coefficient,
63 ;; Rational expressions
64 ;; ====================
66 ;; A rational expression is just a quotient of two polynomials p/q and, as such,
67 ;; is stored as a cons pair (p . q), where p and q are in the format described
68 ;; above. Since bare coefficients are allowed as polynomials, we can represent
69 ;; zero as (0 . 1), which literally means 0/1.
71 ;; These format are also documented in the "Introduction to Polynomials" page of
76 ;; Returns true if E (which is hopefully a polynomial expression) should be
77 ;; thought of as a bare coefficient.
78 (defmacro pcoefp
(e) `(atom ,e
))
82 ;; Return true iff the polynomial X is syntactically the zero polynomial. This
83 ;; only happens when the polynomial is a bare coefficient and that coefficient
85 (declaim (inline pzerop
))
88 ((fixnump x
) (zerop x
))
90 ((floatp x
) (zerop x
))))
94 ;; A simple macro that evaluates to the zero polynomial.
99 ;; TERMS should be a list of terms of a polynomial. Returns T if that list is
100 ;; empty, so the polynomial has no terms.
101 (defmacro ptzerop
(terms) `(null ,terms
))
105 ;; A simple macro that evaluates to an empty list of polynomial terms,
106 ;; representing the zero polynomial.
107 (defmacro ptzero
() '())
111 ;; Return the negation of a coefficient, which had better be numeric.
112 (defmacro cminus
(c) `(- ,c
))
116 ;; Return T if the coefficient C is negative. Only works if C is a real number.
117 (defmacro cminusp
(c) `(minusp ,c
))
122 ;; Retrieve a stored value from the given symbol, stored by VALPUT. This is used
123 ;; in the rational function code, which uses it to store information on gensyms
124 ;; that represent variables.
126 ;; Historical note from 2000 (presumably wfs):
128 ;; The PDP-10 implementation used to use the PRINTNAME of the gensym as a
129 ;; place to store a VALUE. Somebody changed this to value-cell instead, even
130 ;; though using the value-cell costs more. Anyway, in NIL I want it to use the
131 ;; property list, as thats a lot cheaper than creating a value cell. Actually,
132 ;; better to use the PACKAGE slot, a kludge is a kludge right?
133 (defmacro valget
(item)
134 `(symbol-value ,item
))
138 ;; Store a value on the given symbol, which can be later retrieved by
139 ;; valget. This is used by the rational function code.
140 (defmacro valput
(item val
)
141 `(setf (symbol-value ,item
) ,val
))
145 ;; Test whether one symbol should occur before another in a canonical ordering.
147 ;; A historical note from Richard Fateman, on the maxima list, 2006/03/17:
149 ;; "The name pointergp comes from the original hack when we wanted a bunch of
150 ;; atoms that could be ordered fast, we just generated, say, 10 gensyms. Then
151 ;; we sorted them by the addresses of the symbols in memory. Then we
152 ;; associated them with x,y,z,.... This meant that pointergp was one or two
153 ;; instructions on a PDP-10, in assembler."
155 ;; "That version of pointergp turned out to be more trouble than it was worth
156 ;; because we sometimes had to interpolate between two gensym "addresses" and
157 ;; to do that we had to kind of renumber too much of the universe. Or maybe
158 ;; we just weren't clever enough to do it without introducing bugs."
160 ;; Richard Fateman also says pointergp needs to be fast because it's called a
161 ;; lot. So if you get an error from pointergp, it's probably because someone
162 ;; forgot to initialize things correctly.
163 (declaim (inline pointergp
))
164 (defun pointergp (a b
)
165 (> (symbol-value a
) (symbol-value b
)))
169 ;; V should be a symbol. If V has an "algebraic value" (stored in the TELLRAT
170 ;; property) then return it, provided that the $ALGEBRAIC flag is
171 ;; true. Otherwise, return NIL.
173 `(and $algebraic
(get ,v
'tellrat
)))
177 ;; Expands to the zero rational expression (literally 0/1)
178 (defmacro rzero
() ''(0 .
1))
182 ;; Test whether a rational expression is zero. A quotient p/q is zero iff p is
184 (defmacro rzerop
(a) `(pzerop (car ,a
)))
188 ;; Calculate the primitive part of a polynomial. This is the polynomial divided
189 ;; through by the greatest common divisor of its coefficients.
190 (defmacro primpart
(p) `(second (oldcontent ,p
)))
194 ;; A convenience macro for constructing polynomials. VAR is the main variable
195 ;; and should be a symbol. With no other arguments, it constructs the polynomial
196 ;; representing the linear polynomial "VAR".
198 ;; A single optional argument is taken to be a coefficient list and, if the list
199 ;; is known at compile time, some tidying up is done with PSIMP.
201 ;; With two optional arguments, they are taken to be an exponent/coefficient
202 ;; pair. So (make-poly 'x 3 2) represents 2*x^3.
204 ;; With all three optional arguments, the first two are interpreted as above,
205 ;; but this coefficient is prepended to an existing list of terms passed in the
208 ;; Note: Polynomials are normally assumed to have terms listed in descending
209 ;; order of exponent. MAKE-POLY does not ensure this, so
210 ;; (make-poly 'x 1 2 '(2 1)) results in '(x 1 2 2 1), for example.
211 (defmacro make-poly
(var &optional
(terms-or-e nil options?
) (c nil e-c?
) (terms nil terms?
))
212 (cond ((null options?
) `(cons ,var
'(1 1)))
213 ((null e-c?
) `(psimp ,var
,terms-or-e
))
214 ((null terms?
) `(list ,var
,terms-or-e
,c
))
215 (t `(psimp ,var
(list* ,terms-or-e
,c
,terms
)))))
219 ;; Extract the main variable from the polynomial P. Note: this does not work for
220 ;; a bare coefficient.
221 (defmacro p-var
(p) `(car ,p
))
225 ;; Extract the list of terms from the polynomial P. Note: this does not work for
226 ;; a bare coefficient.
227 (defmacro p-terms
(p) `(cdr ,p
))
231 ;; Extract the leading coefficient of the polynomial P. Note: this does not work for
232 ;; a bare coefficient.
233 (defmacro p-lc
(p) `(caddr ,p
))
237 ;; Extract the leading exponent or degree of the polynomial P. Note: this does
238 ;; not work for a bare coefficient.
239 (defmacro p-le
(p) `(cadr ,p
))
243 ;; Extract the terms of the polynomial P, save the leading term.
244 (defmacro p-red
(p) `(cdddr ,p
))
248 ;; Extract the leading coefficient from TERMS, a list of polynomial terms.
249 (defmacro pt-lc
(terms) `(cadr ,terms
))
253 ;; Extract the leading exponent (or degree) from TERMS, a list of polynomial
255 (defmacro pt-le
(terms) `(car ,terms
))
259 ;; Return all but the leading term of TERMS, a list of polynomial terms.
260 (defmacro pt-red
(terms) `(cddr ,terms
))
264 ;; Sum one or more rational expressions with RATPL
267 (t `(ratpl ,r
(r+ ,@l
)))))
271 ;; Take the product of one or more rational expressions with RATTI.
274 (t `(ratti ,r
(r* ,@l
) t
))))
278 ;; Subtract the sum of the second and following rational expressions from the
279 ;; first, using RATDIF.
281 (cond ((null l
) `(ratminus (ratfix ,r
)))
282 (t `(ratdif (ratfix ,r
) (r+ ,@l
)))))