1;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;; 2;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3;;; The data in this file contains enhancments. ;;;;; 4;;; ;;;;; 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;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 10 11(in-package :maxima) 12 13(macsyma-module ratmac macro) 14 15;; Polynomials 16;; =========== 17;; 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 21;; coefficients. 22;; 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 26;; 27;; ($x 2 1 1 ((%sin) $x)), 28;; 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 32;; below. 33;; 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]. 37;; 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 45;; function code. 46;; 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). 52;; 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. 56;; 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, 61;; respectively). 62;; 63;; Rational expressions 64;; ==================== 65;; 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. 70;; 71;; These format are also documented in the "Introduction to Polynomials" page of 72;; the manual. 73 74;; PCOEFP 75;; 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)) 79 80;; PZEROP 81;; 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 84;; is zero. 85(declaim (inline pzerop)) 86(defun pzerop (x) 87 (cond 88 ((fixnump x) (zerop x)) 89 ((consp x) nil) 90 ((floatp x) (zerop x)))) 91 92;; PZERO 93;; 94;; A simple macro that evaluates to the zero polynomial. 95(defmacro pzero () 0) 96 97;; PTZEROP 98;; 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)) 102 103;; PTZERO 104;; 105;; A simple macro that evaluates to an empty list of polynomial terms, 106;; representing the zero polynomial. 107(defmacro ptzero () '()) 108 109;; CMINUS 110;; 111;; Return the negation of a coefficient, which had better be numeric. 112(defmacro cminus (c) `(- ,c)) 113 114;; CMINUSP 115;; 116;; Return T if the coefficient C is negative. Only works if C is a real number. 117(defmacro cminusp (c) `(minusp ,c)) 118 119 120;; VALGET 121;; 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. 125;; 126;; Historical note from 2000 (presumably wfs): 127;; 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)) 135 136;; VALPUT 137;; 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)) 142 143;; POINTERGP 144;; 145;; Test whether one symbol should occur before another in a canonical ordering. 146;; 147;; A historical note from Richard Fateman, on the maxima list, 2006/03/17: 148;; 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." 154;; 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." 159;; 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))) 166 167;; ALGV 168;; 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. 172(defmacro algv (v) 173 `(and $algebraic (get ,v 'tellrat))) 174 175;; RZERO 176;; 177;; Expands to the zero rational expression (literally 0/1) 178(defmacro rzero () ''(0 . 1)) 179 180;; RZEROP 181;; 182;; Test whether a rational expression is zero. A quotient p/q is zero iff p is 183;; zero. 184(defmacro rzerop (a) `(pzerop (car ,a))) 185 186;; PRIMPART 187;; 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))) 191 192;; MAKE-POLY 193;; 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". 197;; 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. 200;; 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. 203;; 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 206;; third argument. 207;; 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))))) 216 217;; P-VAR 218;; 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)) 222 223;; P-TERMS 224;; 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)) 228 229;; P-LC 230;; 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)) 234 235;; P-LE 236;; 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)) 240 241;; P-RED 242;; 243;; Extract the terms of the polynomial P, save the leading term. 244(defmacro p-red (p) `(cdddr ,p)) 245 246;; PT-LC 247;; 248;; Extract the leading coefficient from TERMS, a list of polynomial terms. 249(defmacro pt-lc (terms) `(cadr ,terms)) 250 251;; PT-LE 252;; 253;; Extract the leading exponent (or degree) from TERMS, a list of polynomial 254;; terms. 255(defmacro pt-le (terms) `(car ,terms)) 256 257;; PT-RED 258;; 259;; Return all but the leading term of TERMS, a list of polynomial terms. 260(defmacro pt-red (terms) `(cddr ,terms)) 261 262;; R+ 263;; 264;; Sum one or more rational expressions with RATPL 265(defmacro r+ (r . l) 266 (cond ((null l) r) 267 (t `(ratpl ,r (r+ ,@l))))) 268 269;; R* 270;; 271;; Take the product of one or more rational expressions with RATTI. 272(defmacro r* (r . l) 273 (cond ((null l) r) 274 (t `(ratti ,r (r* ,@l) t)))) 275 276;; R- 277;; 278;; Subtract the sum of the second and following rational expressions from the 279;; first, using RATDIF. 280(defmacro r- (r . l) 281 (cond ((null l) `(ratminus (ratfix ,r))) 282 (t `(ratdif (ratfix ,r) (r+ ,@l))))) 283