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 1982 Massachusetts Institute of Technology ;;; 9;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 10 11(in-package :maxima) 12 13(macsyma-module series) 14 15(declare-top (special var *n *a *m *c *index $cauchysum *gcd* 16 nn* dn* $ratsimpexpons *infsumsimp *roots *failures 17 *ratexp *var usexp $verbose ans *trigred 18 *form indl *noexpand $ratexpand)) 19 20(load-macsyma-macros rzmac) 21 22;;****************************************************************************** 23;; driver stage 24;;****************************************************************************** 25 26(defmfun $powerseries (expr var pt) 27 (when (and (atom var) (not (symbolp var))) 28 (improper-arg-err var '$powerseries)) 29 (powerseries expr var pt)) 30 31;; An error that can be raised deep in the bowels of POWERSERIES that we'll 32;; catch and return a noun form. 33(define-condition powerseries-expansion-error (error) 34 ((message :initarg :message :initform nil))) 35 36(defun powerseries-expansion-error (&optional message &rest arguments) 37 (error 'powerseries-expansion-error 38 :message (when message (cons message arguments)))) 39 40;; The top-level routine for $powerseries, which calls this function after 41;; spotting invalid arguments. 42(defun powerseries (expr var pt) 43 (handler-case 44 (if (and (signp e pt) (symbolp var)) 45 (seriesexpand* expr) 46 (cond 47 ((signp e pt) 48 (sbstpt expr 'x var var)) 49 ((eq pt '$inf) 50 (sbstpt expr (m// 1 'x) var (div* 1 var))) 51 ((eq pt '$minf) 52 (sbstpt expr (m- (m// 1 'x)) var (m- (div* 1 var)))) 53 (t 54 (sbstpt expr (m+ 'x pt) var (simplifya (m- var pt) nil))))) 55 56 ;; If expansion fails, print a message in verbose mode but then return a 57 ;; noun form. 58 (powerseries-expansion-error (e) 59 (when $verbose 60 (mtell (intl:gettext "Failed to expand ~M with respect to ~M at ~M.~%") 61 expr var pt) 62 (with-slots (message) e 63 (when message 64 (terpri) 65 (finish-output) 66 (apply #'mtell message)))) 67 `((%powerseries) ,expr ,var ,pt)))) 68 69;; Return a list of the terms in the expression that are used as integration or 70;; differentiation variables. 71(defun intdiff-vars-in-expr (expr) 72 (cond 73 ((atom expr) nil) 74 75 ;; Arguably, if this is a definite integral, the variable isn't free so we 76 ;; shouldn't return it. But maxima-substitute doesn't know about bound 77 ;; variables and this is a helper function to avoid silly substitutions. 78 ((eq (caar expr) '%integrate) (list (caddr expr))) 79 80 ((eq (caar expr) '%derivative) 81 (loop for x in (cddr expr) by #'cddr collecting x)) 82 83 (t 84 (reduce #'union (cdr expr) :key #'intdiff-vars-in-expr)))) 85 86;; Series expand EXP after substituting in SEXP for VAR1, the main 87;; variable. SEXP should be an expression in 'X and 'X will be replaced by a 88;; gensym for the calculation. When the expansion is done, substitute USEXP for 89;; the gensym. 90(defun sbstpt (exp sexp var1 usexp) 91 (let* ((var (gensym)) 92 (sub-exp (maxima-substitute (subst var 'x sexp) var1 exp)) 93 (expanded (seriesexpand* sub-exp))) 94 ;; If we've ended up with diff(foo, var) then we give up (we can't 95 ;; substitute arbitrary expressions for the differentiation / integration 96 ;; variable). 97 (if (memq var (intdiff-vars-in-expr expanded)) 98 (powerseries-expansion-error 99 (intl:gettext "~ 100Couldn't make substitution to evaluate at the given point because the~%~ 101power series expansion contained the expansion variable as an~%~ 102integration / differentiation variable.")) 103 (maxima-substitute usexp var expanded)))) 104 105(defun seriesexpand* (x) 106 (let ((*index (gensumindex)) *n *a *m *c 107 ($cauchysum t) ($ratsimpexpons t) 108 $ratexpand *infsumsimp *ratexp *trigred *noexpand) 109 (meval `(($declare) ,*index $integer)) 110 111 (sp2expand (seriespass1 x)))) 112 113(defun out-of (e) 114 (let ((e (cond ((and (boundp '*var) *var) 115 (subst (list '(mexpt) *var *gcd*) var e)) 116 (t e))) 117 (var (cond ((and (boundp '*var) *var)) (t var)))) 118 (cond ((and (boundp 'usexp) usexp) 119 (subst usexp var e)) 120 (t e)))) 121 122(defun show-exp (e) 123 (mtell "~%~%~M~%~%" 124 (list '(mlabel) nil (out-of e)))) 125 126(defun seriespass1 (e) 127 (let ((w (sratsimp (sp1 e)))) 128 (when $verbose 129 (terpri) 130 (finish-output) 131 (mtell (intl:gettext "powerseries: first simplification returned ~%")) 132 (show-exp w)) 133 w)) 134 135;; 136;;***************************************************************************** 137;; pass two expansion phase 138;;***************************************************************************** 139;; 140 141(defun sp2expand (exp) 142 (cond ((or (free exp var) (atom exp)) exp) 143 ((mbagp exp) (cons (car exp) (mapcar #'sp2expand (cdr exp)))) 144 ((eq (caar exp) 'mplus) (m+l (mapcar #'sp2expand (cdr exp)))) ;was below 'mtimes test--fixes powerseries(1+x^n,x,0) 145 ((sratp exp var) (ratexp exp)) 146 ((eq (caar exp) 'mexpt) (sp2expt exp)) 147 ((zl-get (caar exp) 'sp2) (sp2sub (sp2trig exp) (cadr exp))) 148 ((eq (caar exp) 'mtimes) (m*l (mapcar #'sp2expand (cdr exp)))) 149 ((eq (caar exp) '%log) (sp2log (cadr exp))) 150 ((eq (caar exp) '%derivative) (sp2diff (cadr exp) (cddr exp))) 151 ((eq (caar exp) '%integrate) 152 (sp2integ (cadr exp) (caddr exp) (cdddr exp))) 153 ((member (caar exp) '(%product %sum) :test #'eq) 154 (list* (car exp) (sp2expand (cadr exp)) (cddr exp))) 155 (t (list '(%sum) 156 (m* (m^ var *index) 157 (m^ (list '(mfactorial) *index) -1) 158 (list '(%at) (list '(%derivative) exp var *index) 159 (list '(mequal) var 0))) 160 *index 0 '$inf)))) 161 162(defun sp2sub (s exp) 163 (unless (smono exp var) 164 (powerseries-expansion-error 165 (intl:gettext "Can't substitute the value~%~M~%~ 166 into another expansion because it isn't a monomial in the~ 167 expansion variable.") 168 (out-of exp))) 169 (maxima-substitute exp 'sp2var (simplify s))) 170 171(defun ratexp (exp) 172 (let (*gcd*) 173 (if $verbose 174 (mtell (intl:gettext "powerseries: attempt rational function expansion of~%~M") 175 (list '(mlabel) nil exp))) 176 (multiple-value-call #'sratexpnd (numden exp)))) 177 178(defun sratexpnd (n d) 179 (let ((ans (list nil)) 180 (*splist*) 181 ;; A pattern that matches cc*(c*x^m + a)^n 182 (linpat 183 '((mtimes) ((coefftt) (cc not-zero-free var)) 184 ((mexpt) ((mplus) ((coeffpt) 185 (w m1 ((mexpt) (x equal var) 186 (m not-zero-free var))) 187 (c freevar)) 188 ((coeffpp) (a freevar))) 189 (n not-zero-free var))))) 190 191 (declare (special *splist*)) 192 (cond ((and (not (equal n 1)) (smono n var)) 193 (m* n (sratexpnd 1 d))) 194 ((free d var) 195 (cond ((poly? n var) 196 (m// n d)) 197 ((m1 n linpat) 198 (m// (srbinexpnd (cdr ans)) d)) 199 (t 200 (powerseries-expansion-error)))) 201 ((smonop d var) 202 (cond ((mplusp n) 203 (m+l (mapcar #'(lambda (q) (div* q d)) (cdr n)))) 204 (t (m// n d)))) 205 ((not (equal 1 (setq *gcd* (ggcd (nconc (exlist n) (exlist d)))))) 206 (sratsubst *gcd* n d)) 207 208 ;; Rational expansion theorem for distinct roots. The main point of 209 ;; the POLY? test here is that it guarantees $HIPOW will give the 210 ;; correct degrees (in particular, it only allows known non-negative 211 ;; integer exponents) 212 ((and (poly? n var) 213 (poly? d var) 214 (> ($hipow d var) ($hipow n var)) 215 (has-distinct-nonzero-roots-p d var)) 216 (expand-distinct-roots n d)) 217 218 ;; The SRBINEXPND call above dealt with expressions of the form 219 ;; cc*(c*x^m+a)^n. Here, we deal with b/(cc*(c*x^m+a)^n). If you 220 ;; explicitly write a polynomial like that, the simplifier will 221 ;; rewrite it as (..)^(-n) before it gets here, but things like 222 ;; 1/sqrt(x+1) won't have been rewritten successfully, so we catch 223 ;; them here. 224 ((and (free n var) 225 (prog2 (setq d (let (($ratfac t)) 226 (ratdisrep ($rat (factor d) var)))) 227 (m1 d linpat))) 228 ;; We had num/den and LINPAT matched den. We need to replace cc 229 ;; with num/cc and n with -n. 230 (setf (cdr (assoc 'n ans)) (m- (cdr (assoc 'n ans))) 231 (cdr (assoc 'cc ans)) (m// n (cdr (assoc 'cc ans)))) 232 (srbinexpnd (cdr ans))) 233 234 (t 235 ;; *RATEXP is set by RATEXPAND1, which we call to do the general 236 ;; expansion below. This check makes sure we can't end up recursing 237 ;; infinitely. 238 (when *ratexp 239 (powerseries-expansion-error)) 240 241 ;; Look for a power of var (a zero root) as a term. We already know 242 ;; that d isn't a monomial, so we can only have a zero root if d is 243 ;; a product. We also know that d is not atomic because otherwise 244 ;; we'd have taken the (free d var) clause above. 245 (let ((zero-root 246 (and (eq (caar d) 'mtimes) 247 (find-if (lambda (factor) 248 (or (eq factor var) 249 (and (mexptp factor) 250 (eq (cadr factor) var)))) 251 (cdr d))))) 252 (if (not zero-root) 253 (ratexpand1 n d) 254 (let ((other-factors (remove zero-root (cdr d) 255 :test #'eq :count 1))) 256 (m* (sratexpnd n (m*l other-factors)) 257 `((mexpt) ,zero-root -1))))))))) 258 259; is a sum with index and bounds from psp2form 260(defun psp2formp (exp) 261 (and (listp exp) 262 (listp (car exp)) 263 (eq (caar exp) '%sum) 264 (eq (caddr exp) *index) 265 (eql (cadddr exp) 0) 266 (eq (cadr (cdddr exp)) '$inf))) 267 268;; A stripped down version of (sumcontract (intosum EXP)). Coalesce occurrences 269;; of ((%SUM) <FOO> *INDEX 0 '$INF), which are allowed to start with 270;; multiplicative constants. So something like 271;; 272;; sum(a(i), i, 0, inf) + 2*sum(b(i), i, 0, inf) + c 273;; 274;; turns into 275;; 276;; sum(a(i) + 2*b(i), i, 0, inf) + c 277;; 278;; This is good enough to collect up the multiple infinite sums that you get 279;; when expanding a power series by partial fractions. It's (obviously) not 280;; clever enough to rename index variables or adjust ranges. 281 282(defun psp2foldsum (exp) 283 (when $verbose 284 (mtell (intl:gettext "powerseries: preparing to fold sums~%")) 285 (show-exp exp)) 286 287 (multiple-value-bind (sums others) 288 (partition-by (lambda (e) 289 (or (psp2formp e) 290 (and (mtimesp e) 291 (psp2formp (caddr e))))) 292 (cdr exp)) 293 (if (null sums) 294 ;; If there were no sums, just return the original expression. 295 exp 296 ;; Otherwise contract the sums 297 (let ((contracted 298 (list '(%sum) 299 (m+l (mapcar (lambda (e) 300 (if (eq (caar e) 'mtimes) 301 (m* (cadr e) (cadr (caddr e))) 302 (cadr e))) 303 sums)) 304 *index 0 '$inf))) 305 (if (null others) 306 contracted 307 (m+l (cons contracted others))))))) 308 309; solve returns a list: (soln mult soln mult ...) 310; distinct-nonzero-roots-p returns true if every 311; soln is not nonzero and every mult is 1 312(defun distinct-nonzero-roots-p (roots) 313 (or (null roots) 314 (and 315 ;; Check all roots have multiplicity one and are nonzero. 316 (loop 317 for (root mult) on roots by #'cddr 318 always (and (eql 1 mult) 319 (eq ($sign `((mabs) ,(caddr root))) '$pos))) 320 321 ;; Now check that the roots are genuinely different from one another 322 ;; (eg. if I have two symbolic roots, A and B, then I don't need to know 323 ;; the values of A and B, but I do need to know that they aren't equal) 324 (loop 325 for rest on roots by #'cddr 326 always (loop 327 for b in (cddr rest) by #'cddr 328 with a = (car rest) 329 always (eq '$pos 330 ($sign `((mabs) ,(m- (caddr b) (caddr a)))))))))) 331 332; returns t if polynomial poly in variable var has all distinct roots 333(defun has-distinct-nonzero-roots-p (poly var) 334 (let ((*roots nil) 335 (*failures nil)) 336 (solve poly var 1) 337 (and (not *failures) 338 (distinct-nonzero-roots-p *roots)))) 339 340; Rational Expansion Theorem for Distinct Roots 341; Graham, Knuth, Patashnik, "Concrete Math" 2nd ed, p 340 342; 343; If R(z) = P(z)/Q(z), where Q(z) = q0*(1-r_1*z)*...*(1-r_l*z) and the 344; numbers (r_1 ... r_l) are distinct, and if P(z) is a polynomial of degree less 345; than l, then 346; [z^n]R(z) = a_1*r_1^n + ... + a_l*r_l^n, where a_k = -r_k*P(1/r_k)/Q'(1/r_k) 347(defun expand-distinct-roots (num den) 348 (let ((*roots nil) 349 (*failures nil)) 350 (solve den var 1) 351 (cond (*failures (error "EXPAND-DISTINCT-ROOTS: failed to solve for roots.")) 352 ((distinct-nonzero-roots-p *roots) 353 (psp2form (m+l (mapcar #'(lambda (r) 354 (and $verbose 355 (prog2 356 (mtell (intl:gettext "powerseries: for root:~%")) 357 (show-exp r) 358 (mtell (intl:gettext "powerseries: numerator at root =~%")) 359 (show-exp (maxima-substitute r var num)) 360 (mtell (intl:gettext "powerseries: first derivative of denominator at root =~%")) 361 (show-exp (maxima-substitute r var ($diff den var))))) 362 (m* -1 363 (m// 1 r) 364 (maxima-substitute r var num) 365 (m// 1 (maxima-substitute r var ($diff den var))) 366 (m^ (m// 1 r) *index))) 367 (mapcar #'caddr (remove-mult *roots)))) 368 *index 0)) 369 (t (error "EXPAND-DISTINCT-ROOTS: roots are not distinct.~%"))))) 370 371;; Try to expand N/D as a power series in VAR about 0. 372;; 373;; Can safely assume that D has no zero roots (we remove them in SRATEXPND 374;; before calling this function). 375(defun ratexpand1 (n d) 376 (when $verbose 377 (mtell (intl:gettext 378 "powerseries: attempt partial fraction expansion of ")) 379 (show-exp (list '(mquotient) n d)) 380 (terpri) 381 (finish-output)) 382 383 ;; *RATEXP serves as a recursion guard: if SRATEXPND is about to call us 384 ;; *recursively, the flag will be set and it gives up instead. 385 (let* ((*ratexp t) 386 (fractions ($partfrac (div* n d) var))) 387 (cond 388 ;; If $PARTFRAC succeeds, it will return a sum of terms. In that case, 389 ;; expand each one (which we assume is going to be easier than what we 390 ;; started with) and try to fold the result into a single power series sum 391 ;; again. 392 ((mplusp fractions) 393 (when $verbose 394 (mtell (intl:gettext "which is ~%")) 395 (show-exp fractions)) 396 (psp2foldsum (m+l (mapcar #'ratexp (cdr fractions))))) 397 398 ;; If that didn't work, maybe it's because the numerator messed things 399 ;; up. If we're really lucky, the numerator is actually a polynomial 400 ;; though. In that case, factor it out and do an expansion of 1/d on its 401 ;; own. 402 ((poly? n var) 403 (when $verbose 404 (mtell (intl:gettext 405 "powerseries: partial fraction expansion failed, ~ 406 expanding denominator only.~%"))) 407 (m* n (ratexp (m// 1 d)))) 408 409 ;; If n is complicated, there's not really much we can do to make further 410 ;; progress, so give up and return a noun form. 411 (t (powerseries-expansion-error 412 (intl:gettext "Partial fraction expansion failed")))))) 413 414(defun sratsubst (gcd num den) 415 (and $verbose 416 (prog2 (mtell (intl:gettext "powerseries: substituting for the occurrences of")) 417 (show-exp (list '(mexpt) var gcd)) 418 (mtell (intl:gettext "in")) 419 (show-exp (list '(mquotient) num den)) 420 (terpri) 421 (finish-output))) 422 (funcall #'(lambda (var* *var) 423 (setq num (maxima-substitute (m^ var* (m^ gcd -1)) *var num) 424 den (maxima-substitute (m^ var* (m^ gcd -1)) *var den)) 425 (maxima-substitute (m^ *var gcd) var* 426 (funcall #'(lambda (var) (sratexpnd num den)) var*))) 427 (gensym) var)) 428 429(defun ggcd (l) 430 (cond ((null l) 1) 431 ((null (cdr l)) (car l)) 432 ((equal 1 (car l)) 1) 433 (t ($gcd (ggcd (cdr l)) (car l))))) 434 435(defun exlist (exp) 436 (cond ((null exp) nil) 437 ((atom exp) 438 (and (eq exp var) (ncons 1))) 439 ((and (not (atom (car exp))) (eq (caar exp) 'mplus)) 440 (exlist (cdr exp))) 441 ((smono (car exp) var) 442 (cond ((equal *n 0) (exlist (cdr exp))) 443 (t (cons *n (exlist (cdr exp)))))) 444 (t (exlist (cdr exp))))) 445 446;; Perform a binomial expansion of (c x^m + a)^n. 447;; 448;; If n isn't known to be an integer, we'd like to write a sum of terms 449;; 450;; x^(mk) c^k a^(n-k) / ((n+1) beta(n-k+1, k+1)) 451;; 452;; But we have to be a bit careful: it only works if c, x and a are 453;; nonzero. (Otherwise the simplifier quite reasonably simplifies each of the 454;; terms to zero, because the only nonzero term has the undefined term 0^0 in 455;; it). 456;; 457;; If we're not sure whether cx is zero, we can just split out the first term 458;; from the sum. If a is definitely zero, we can just write down a single 459;; term. If we're not sure about a, it's a bit more difficult: if we know that n 460;; is a positive integer, the sum is finite (and has at least two terms) and we 461;; can just split off the last term. Otherwise, give up. 462(defun srbinexpnd (ans) 463 (alist-bind (n a m c cc x) ans 464 (m* cc 465 (if (and (integerp n) (minusp n)) 466 (srintegexpd (neg n) a m c) 467 (let ((sgn-a ($sign (list '(mabs) a))) 468 (sgn-cx ($sign `((mabs) ((mtimes) ,c ,x)))) 469 (sgn-n ($sign n)) 470 (general-term 471 (m// (m* (m^ var (m* m *index)) 472 (m^ c *index) 473 (m^ a (m- n *index))) 474 (m* (list '($beta) (m- n (m1- *index)) (m1+ *index)) 475 (m1+ n))))) 476 (cond 477 ((eq sgn-n '$zero) 1) 478 ((eq sgn-cx '$zero) (m^ a n)) 479 ((eq sgn-a '$zero) (m* (m^ c n) (m^ x (m* m n)))) 480 481 ((and (eq sgn-a '$pos) (eq sgn-cx '$pos)) 482 (if (and ($featurep n '$integer) 483 (memq sgn-n '($pos $pz))) 484 `((%sum) ,general-term ,*index 0 ,n) 485 `((%sum) ,general-term ,*index 0 $inf))) 486 487 ((eq sgn-a '$pos) 488 (m+ (m^ a n) `((%sum) ,general-term ,*index 1 $inf))) 489 490 ((and ($featurep n '$integer) (eq sgn-n '$pos)) 491 (m+ `((%sum) ,general-term ,*index 0 ,(m1- n)) 492 (m* (m^ c n) (m^ x (m* n m))))) 493 494 (t 495 (powerseries-expansion-error 496 (intl:gettext 497 "Couldn't expand binomial~%~M~%~ 498 as we didn't know which terms were nonzero."))))))))) 499 500(defun psp2form (coeff exp bas) 501 (list '(%sum) (m* coeff (m^ var exp)) *index bas '$inf)) 502 503(defun srintegexpd (n a m c) 504 (and $verbose 505 (prog2 (mtell (intl:gettext "powerseries: apply rule for expressions of ~%")) 506 (show-exp '((mexpt) ((mplus) $a ((mtimes) $c ((mexpt) $var $m))) 507 ((mminus) $n))) 508 (mtell (intl:gettext "powerseries: here we have")) 509 (show-exp (list '(mlist) (list '(mequal) '$n n) (list '(mequal) '$a a) 510 (list '(mequal) '$c c) (list '(mequal) '$m m))))) 511 (cond ((= n 1) 512 (psp2form 513 (m* (m^ a (m* -1 (m+ 1 *index))) 514 (m^ (m* -1 c) *index)) 515 (if (equal m 1) *index (m* *index m)) 516 0)) 517 ((= 2 n) 518 (psp2form (m* (m+ 1 *index) 519 (m^ a (m* -1 (m+ 2 *index))) 520 (m^ (m* -1 c) *index)) 521 (if (equal m 1) *index (m* *index m)) 522 0)) 523 (t (psp2form (m* (do ((nn (1- n) (1- nn)) 524 (l nil (cons (list '(mplus) *index nn) l))) 525 ((= nn 0) 526 (m*l (cons (m// 1 (factorial (1- n))) l)))) 527 (m^ (m* -1 c (m^ a -1)) *index) 528 (m^ a (- n))) 529 (if (equal m 1) *index (m* *index m)) 530 0)))) 531 532(defun sratp (a var) 533 (cond ((atom a) t) 534 ((member (caar a) '(mplus mtimes) :test #'eq) (sandmap (cdr a))) 535 ((eq (caar a) 'mexpt) 536 (cond ((free (cadr a) var) (free (caddr a) var)) 537 ((smono a var) t) 538 ((and (free (caddr a) var) (sratp (cadr a) var))))) 539 (t (and (free (mop a) var) (every #'(lambda (s) (free s var)) (margs a)))))) 540 541(defun sandmap (l) (or (null l) (and (sratp (car l) var) (sandmap (cdr l))))) 542 543(defun sp2trig (exp) (subst *index '*index (zl-get (caar exp) 'sp2))) 544 545;; Take an expression, EXPR, and try to write it as a + b*VAR^c. On success, 546;; returns (VALUES A B C). On failure, raise a powerseries-expansion-error. 547;; 548;; One way to do this would be to call $EXPAND and look at the results, but that 549;; might be rather slow - if the expression is something like (1+x)^10, expand 550;; will take ages and won't help very much. Another approach would be to put it 551;; into CRE form, but that would be a mistake if the input was something like 552;; (1+x)^100... 553;; 554;; Instead, we're a bit stupider: we walk through the expression tree. If we see 555;; anything other than +,* and ^, we give up. If we find we have more than one 556;; different exponent in our terms, we give up. 557;; 558;; Obviously, there are always examples where this won't work, but $EXPAND will 559;; (something like (x+1)^2 - x^2, for example), but I'm assuming that this won't 560;; be something we encounter in practice. 561(defun split-two-term-poly (expr) 562 (cond 563 ((atom expr) 564 (if (eq expr var) 565 (values 0 1 1) 566 (values expr 0 1))) 567 568 ((freeof var expr) 569 (values expr 0 1)) 570 571 ((eq (caar expr) 'mplus) 572 (let ((a 0) (b) (c)) 573 (dolist (arg (cdr expr)) 574 (multiple-value-bind (aa bb cc) (split-two-term-poly arg) 575 (setf a (m+ a aa)) 576 (unless (eql bb 0) 577 (cond 578 ;; Is this the first nonzero power we've seen? 579 ((not b) 580 (setf b bb c cc)) 581 ;; Is this another term with the same power as one we've seen 582 ;; already? 583 ((eql c cc) 584 (setf b (m+ b bb))) 585 ;; Otherwise, give up. 586 (t 587 (powerseries-expansion-error)))))) 588 (values a b c))) 589 590 ((eq (caar expr) 'mtimes) 591 (let ((prod 1) a b c) 592 (dolist (arg (cdr expr)) 593 (multiple-value-bind (aa bb cc) (split-two-term-poly arg) 594 (cond 595 ((eql bb 0) 596 (setf prod (m* prod aa))) 597 ((not a) 598 (setf a aa b bb c cc)) 599 ;; This is the second term (a + b*x^c), so we know we'll end up 600 ;; with mixed terms. (We don't try to spot e.g. (1-x)*(1+x)) 601 (t 602 (powerseries-expansion-error))))) 603 604 (if a 605 (values (m* a prod) (m* b prod) c) 606 (values prod 0 1)))) 607 608 ((eq (caar expr) 'mexpt) 609 ;; We know that EXPR isn't free of VAR. Check that VAR isn't in the 610 ;; exponent. 611 (unless (freeof var (caddr expr)) 612 (powerseries-expansion-error)) 613 (multiple-value-bind (a b c) (split-two-term-poly (cadr expr)) 614 (cond 615 ((eql a 0) 616 (values 0 (m^ b (caddr expr)) (m* c (caddr expr)))) 617 ((eql b 0) 618 (values (m^ a (caddr expr)) 0 1)) 619 (t 620 (powerseries-expansion-error))))) 621 622 (t 623 (powerseries-expansion-error)))) 624 625;; Try to expand log(e) using the series for log(1+x). 626;; 627;; The basic idea is that if a is nonzero then 628;; 629;; log(a + b*x^k) = log(a*(1 + b/a*x^k)) 630;; = log(a) + log(1 + b/a*x^k) 631;; 632;; and we know a series for that. 633(defun sp2log (e) 634 (cond 635 ((or (atom e) (free e var)) 636 (list '(%log) e)) 637 (t 638 (multiple-value-bind (a b k) 639 (split-two-term-poly (specdisrep e)) 640 ;; If a is zero, we can't expand 641 (unless (eq '$positive 642 (asksign (list '(mabs) a))) 643 (powerseries-expansion-error)) 644 645 (let* ((coeff (m* b (m^ a -1))) 646 (coeff-sign ($sign coeff)) 647 (negate-coeff-p)) 648 ;; If we know that the coefficient is not positive, switch the series 649 ;; around (which gets rid of some ugly minus signs). If we're not sure, 650 ;; but the cofficient "looks negative" (so is something like -7*k), 651 ;; switch it around too. 652 (cond 653 ((member coeff-sign '($neg $nz)) 654 (setf negate-coeff-p t 655 coeff (m- coeff))) 656 ((and (not (member coeff-sign '($pos $pz))) 657 (mtimesp coeff) 658 (numberp (cadr coeff)) 659 (minusp (cadr coeff))) 660 (setq negate-coeff-p t 661 coeff (simptimes 662 (list* (car coeff) (- (cadr coeff)) (cddr coeff)) 663 1 t)))) 664 (list '(%sum) 665 (m* -1 666 (m^ (if negate-coeff-p 667 (m* coeff (m^ var k)) 668 (m* -1 coeff (m^ var k))) 669 *index) 670 (m^ *index -1)) 671 *index 1 '$inf)))))) 672 673(defun sp2expt (exp) 674 (let ((base (cadr exp)) 675 (power (caddr exp))) 676 (cond 677 ;; Do (a^b)^c = a^(bc) when c is a number 678 ((and (numberp power) (mexptp base)) 679 (sp2expt (m^ (cadr base) 680 (m* power (caddr base))))) 681 682 ;; If a positive power which is free of the expansion variable and less 683 ;; than the maximum expansion power then expand the base and do the 684 ;; multiplication explicitly. 685 ((and (free power var) 686 (signp g power) 687 (< power $maxposex)) 688 (m*l (dup (sp2expand base) power))) 689 690 ;; If the base is free of VAR, rewrite as a^b = e^(b log(a)) if necessary 691 ;; and then use the tabulated expansion of exp(x). If POWER is a sum then 692 ;; do the rewrite a^(b+c) = a^b a^c 693 ((free base var) 694 (let ((expansion 695 (subst *index '*index 696 (if (eq base '$%e) 697 (get 'mexpt 'sp2) 698 (subst `((mtimes) ((mlog) ,base) sp2var) 'sp2var 699 (get 'mexpt 'sp2)))))) 700 (cond 701 ((not (mplusp power)) 702 (sp2sub expansion power)) 703 704 (t 705 (m*l (mapcar (lambda (summand) (sp2sub expansion summand)) 706 (cdr power))))))) 707 708 (t (powerseries-expansion-error))))) 709 710(defun dup (x %n) 711 (if (= %n 1) 712 (ncons x) 713 (cons x (dup x (1- %n))))) 714 715(defun sp2diff (exp l) 716 (cond 717 ((free exp var) 718 (sp3form (sp2expand exp) (cons '(%derivative) l))) 719 (t 720 ;; If the expression isn't free of VAR, we know how to expand the result if 721 ;; the orders of differentiation are all explicit non-negative 722 ;; integers. Otherwise, give up. 723 (let ((indl) (remaining-derivs)) 724 (loop 725 for (y order) on l by #'cddr 726 do 727 (cond 728 ((not (typep order '(integer 0))) 729 (powerseries-expansion-error)) 730 731 ((not (eq y var)) 732 (setf remaining-derivs (list* order y remaining-derivs))) 733 734 (t 735 (dotimes (i order) 736 (setf indl nil) 737 (setf exp (sp2diff1 (sp2expand exp) nil nil)))))) 738 739 (if remaining-derivs 740 (sp3form exp `((%derivative) 741 ,@(nreverse remaining-derivs))) 742 exp))))) 743 744(defun sp2diff1 (exp ind lol) ;ind is a list of the indices 745 ;lol is a list of the lower limits 746 (cond ((atom exp) (sdiff exp var)) 747 ((eq (caar exp) 'mplus) 748 (cons '(mplus) 749 (mapcar #'(lambda (q) (sp2diff1 q ind lol)) 750 (cdr exp)))) 751 ((eq (caar exp) '%sum) 752 (setq indl (cons (copy-list (cddr exp)) indl)) 753 (sp2diff1 (cadr exp) 754 (cons (caddr exp) ind) 755 (cons (cadddr exp) lol))) 756 (t (sp2diff2 exp ind lol)))) 757 758(defun sp2diff2 (exp ind lol) 759 (let (e fr) 760 (setq e (m2 exp '((mtimes) ((coefftt) (fr freevar)) 761 ((coefftt) (e true)))) 762 fr (cdr (assoc 'fr e :test #'eq)) 763 e (cdr (assoc 'e e :test #'eq))) 764 (sp3reconst 765 (cond ((and (mexptp e) (eq (cadr e) var)) 766 (cond ((equal 0 (mbinding (ind lol) 767 (meval (m* fr (caddr e))))) 768 (m* (sp3substp1 ind ind (m* fr (caddr e))) e)) 769 ((mgrp 1 (mbinding (ind lol) 770 (simplify (mevalatoms (caddr e))))) 771 (m* fr (caddr e) (m^ (cadr e) (m- (caddr e) 1)))) 772 (t (sdiff exp var)))) 773 (t (sdiff exp var)))))) 774 775(defun sp2integ (exp v l) 776 (if (null l) 777 (if (eq var v) 778 (sp2integ1 ($expand (sp2expand exp))) 779 (sp3form (sp2expand exp) (list '(%integrate) v))) 780 (sp2integ2 exp v (car l) (cadr l)))) 781 782(defun sp2integ1 (exp) 783 (let (pair) 784 (cond ((ratp exp var) (ratint exp var)) 785 ((eq (caar exp) 'mplus) 786 (cons '(mplus) (mapcar #'sp2integ1 (cdr exp)))) 787 ((and (eq (caar exp) 'mtimes) 788 (prog2 (setq pair (partition exp var 1)) 789 (not (equal (car pair) 1)))) 790 (mul2* (car pair) (sp2integ1 (cdr pair)))) 791 ((and (eq (caar exp) 'mtimes) 792 (prog2 (setq exp ($intosum exp)) nil))) 793 ((or (not (eq (caar exp) '%sum)) (not (isinop (cadr exp) '%sum))) 794 (sinint exp var)) 795 (t (let ((indl (ncons (cddr exp)))) 796 (sp2integ12 (cadr exp) (ncons (caddr exp)) (ncons (cadddr exp)))))))) 797 798(defun sp2integ12 (exp ind lol) 799 (cond ((atom exp) 800 (sp3reconst (ratint exp var))) 801 ((eq (caar exp) 'mplus) 802 (sp3reconst 803 (m+l (mapcar #'(lambda (q) (sp2integ13 q ind lol)) 804 (cdr exp))))) 805 ((eq (caar exp) '%sum) 806 (setq indl (cons (cddr exp) indl)) 807 (sp2integ12 (cadr exp) 808 (cons (caddr exp) ind) 809 (cons (cadddr exp) lol))) 810 (t (sp3reconst (sp2integ13 exp ind lol))))) 811 812(defun sp2integ13 (exp ind lol) 813 (let (e fr) 814 (setq e (m2 exp '((mtimes) ((coefftt) (fr freevar)) 815 ((coefftt) (e true)))) 816 fr (cdr (assoc 'fr e :test #'eq)) 817 e (cdr (assoc 'e e :test #'eq))) 818 (cond ((and (mexptp e) (eq (cadr e) var)) 819 (cond ((mgrp -1 (mbinding (ind lol) 820 (meval (caddr e)))) 821 (m* (sp3substpn ind ind (m* fr (caddr e)) -1) e)) 822 (t (sinint exp var)))) 823 (t (sinint exp var))))) 824 825;; Substitute NEW for OLD in the expression tree EXPRESSION. This is a bit like 826;; MAXIMA-SUBSTITUTE, except it assumes OLD is a symbol. But the important bit 827;; is that it is bright enough to avoid substituting for bound variables in 828;; %INTEGRATE and %AT forms: even though the symbol might have the same name, 829;; it's thought of as a logically different variable. 830;; 831;; The function returns two values: the new expression and a flag that says 832;; whether that expression has changed. If the flag is true on a recursive call, 833;; int-diff-substitute resimplifies the result. 834(defun int-diff-substitute (new old expression) 835 (cond 836 ((eq expression old) (values new t)) 837 ((atom expression) (values expression nil)) 838 (t 839 (let ((op (caar expression)) 840 (args (cdr expression))) 841 (if (or (and (eq op '%integrate) (eq old (second args))) 842 (and (eq op '%at) 843 (not (atom (second args))) 844 (if (eq (caar (second args)) 'mlist) 845 ;; (second args) looks like ((mlist) ((mequal) ...) ...) 846 (memq old (mapcar #'second (rest (second args)))) 847 ;; (second args) looks like ((mequal) ...) 848 (eq old (second (second args)))))) 849 (values expression nil) 850 (let* ((some-changed-p nil) 851 (new-args 852 (mapcar (lambda (x) 853 (multiple-value-bind (new-val changed-p) 854 (int-diff-substitute new old x) 855 (when changed-p 856 (setf some-changed-p t)) 857 new-val)) 858 (cdr expression)))) 859 (if (not some-changed-p) 860 (values expression nil) 861 (values 862 (simplifya (cons (list op) new-args) nil) t)))))))) 863 864;; Expand a definite integral, integrate(exp, v, lo, hi). 865(defun sp2integ2 (exp v lo hi) 866 ;; Deal with the case where the integration variable is also the expansion 867 ;; variable. Something's a bit odd if we're doing this, but we assume that the 868 ;; aliasing was accidental and that the (bound) integration variable should be 869 ;; called something else. 870 (when (eq v var) 871 (setq v (gensym) 872 exp (subst v var exp))) 873 874 (when (boundp '*sp2integ-recursion-guard*) 875 (powerseries-expansion-error 876 (intl:gettext "Recursion when trying to expand the definite integral: ~M") 877 (out-of (symbol-value '*sp2integ-recursion-guard*)))) 878 879 (cond 880 ((not (and (free lo var) (free hi var))) 881 ;; Suppose one or both of the endpoints involves VAR. We'll give up unless 882 ;; they are monomials in VAR (because substituting a non-monomial into a 883 ;; power series is more difficult). If they *are* monomials in VAR, then we 884 ;; try to expand the integral as a power series and find an antiderivative. 885 (let ((lo-exp (sp2expand lo)) 886 (hi-exp (sp2expand hi)) 887 (*sp2integ-recursion-guard* exp)) 888 (declare (special *sp2integ-recursion-guard*)) 889 890 (unless (and (smono lo-exp var) (smono hi-exp var)) 891 (powerseries-expansion-error 892 (intl:gettext "Endpoints of definite integral ~M are not monomials in ~ 893 the expansion variable.") (out-of exp))) 894 895 ;; Since this is a formal powerseries calculation, we needn't concern 896 ;; ourselves with radii of convergence here. Just expand the integrand 897 ;; about zero. (Is there something cleverer we could do?) 898 (let ((anti-derivative (sinint ($powerseries exp v 0) v))) 899 ;; If the expansion had failed, we'd have thrown an error to top-level, 900 ;; so we can assume that ANTI-DERIVATIVE is a sum of monomials in 901 ;; V. Substituting in LO-EXP and HI-EXP, we'll get two sums of 902 ;; monomials in VAR. The difference of the two expressions isn't quite 903 ;; of the right form, but hopefully it's close enough for any other 904 ;; code that uses it. 905 (m- (int-diff-substitute hi-exp v anti-derivative) 906 (int-diff-substitute lo-exp v anti-derivative))))) 907 908 ((free exp var) 909 (list '(%integrate) (subst var v exp) var lo hi)) 910 911 (t 912 (sp3form (sp2expand exp) 913 (list '(%integrate) v lo hi))))) 914 915;;************************************************************************************ 916;; phase three miscellaneous garbage and final simplification 917;;************************************************************************************ 918 919(defun sp3reconst (e) 920 (do ((l indl (cdr l)) (e e (list* '(%sum) e (car l)))) 921 ((null l) e))) 922 923(defun sp3substpn (vars vals exp n) 924 (sp3subst vars (mapcar #'(lambda (q) (add2* q n)) vals) exp)) 925 926(defun sp3substp1 (vars vals exp) (sp3substpn vars vals exp 1)) 927 928(defun sp3subst (vars vals exp) 929 (simplify (sublis (mapcar #'cons (cdr vars) (cdr vals)) 930 (subst (car vals) (car vars) exp)))) 931 932(defun sp3form (e *form) (sp3form1 e)) 933 934(defun sp3form1 (e) 935 (cond ((atom e) (list* (car *form) e (cdr *form))) 936 ((eq (caar e) 'mplus) 937 (cons '(mplus) (mapcar #'sp3form1 (cdr e)))) 938 ((eq (caar e) '%sum) 939 (list* '(%sum) (sp3form1 (cadr e)) (cddr e))) 940 (t (list* (car *form) e (cdr *form))))) 941 942;; These are the series expansions for circular functions 943 944(defprop %sin 945 ((%sum) ((mtimes) 946 ((mexpt) -1 *index) 947 ((mexpt) ((mfactorial) ((mplus) ((mtimes) 2 *index) 1)) -1) 948 ((mexpt) sp2var ((mplus) ((mtimes) 2 *index) 1))) 949 *index 0 $inf) 950 sp2) 951 952(defprop %cos 953 ((%sum) ((mtimes) ((mexpt) -1 *index) 954 ((mexpt) ((mfactorial) ((mtimes) 2 *index)) -1) 955 ((mexpt) sp2var ((mtimes) 2 *index))) 956 *index 0 $inf) 957 sp2) 958 959(defprop %tan 960 ((%sum) ((mtimes) ((mexpt) -1 ((mplus) *index -1)) 961 ((mexpt) 2 ((mtimes) 2 *index)) 962 ((mplus) ((mexpt) 2 ((mtimes) 2 *index)) -1) 963 ((mexpt) ((mfactorial) ((mtimes) 2 *index)) -1) 964 (($bern) ((mtimes) 2 *index)) 965 ((mexpt) sp2var ((mplus) ((mtimes) 2 *index) -1))) 966 *index 0 $inf) 967 sp2) 968 969(defprop %csc 970 ((%sum) ((mtimes) 2 971 ((mexpt) -1 ((mplus) *index -1)) 972 ((mplus) ((mexpt) 2 ((mplus) ((mtimes) 2 *index) -1)) -1) 973 ((mexpt) ((mfactorial) ((mtimes) 2 *index)) -1) 974 (($bern) ((mtimes) 2 *index)) 975 ((mexpt) sp2var ((mplus) ((mtimes) 2 *index) -1))) 976 *index 0 $inf) 977 sp2) 978 979(defprop %cot 980 ((%sum) ((mtimes) 981 ((mexpt) -1 *index) 982 ((mexpt) 2 ((mtimes) 2 *index)) 983 ((mexpt) ((mfactorial) ((mtimes) 2 *index)) -1) 984 (($bern) ((mtimes) 2 *index)) 985 ((mexpt) sp2var ((mplus) ((mtimes) 2 *index) -1))) 986 *index 0 $inf) 987 sp2) 988 989(defprop %sec 990 ((%sum) ((mtimes) ((mexpt) -1 *index) 991 ((mexpt) ((mfactorial) ((mtimes) 2 *index)) -1) 992 (($euler) ((mtimes) 2 *index)) 993 ((mexpt) sp2var ((mtimes) 2 *index))) 994 *index 0 $inf) 995 sp2) 996 997;; These are the series definitions of exponential functions. 998 999(defprop mexpt 1000 ((%sum) 1001 ((mtimes) ((mexpt) ((mfactorial) *index) -1) ((mexpt) sp2var *index)) 1002 *index 0 $inf) 1003 sp2) 1004 1005(defprop %sinh 1006 ((%sum) ((mtimes) 1007 ((mexpt) ((mfactorial) ((mplus) ((mtimes) 2 *index) 1)) -1) 1008 ((mexpt) sp2var ((mplus) ((mtimes) 2 *index) 1))) 1009 *index 0 $inf) 1010 sp2) 1011 1012(defprop %cosh 1013 ((%sum) ((mtimes) 1014 ((mexpt) ((mfactorial) ((mtimes) 2 *index)) -1) 1015 ((mexpt) sp2var ((mtimes) 2 *index))) 1016 *index 0 $inf) 1017 sp2) 1018 1019(defprop %tanh 1020 ((%sum) 1021 ((mtimes) ((mexpt) 4 *index) 1022 ((mplus) ((mexpt) 4 *index) -1) 1023 (($bern) ((mtimes) 2 *index)) 1024 ((mexpt) sp2var ((mplus) ((mtimes) 2 *index) -1)) 1025 ((mexpt) 1026 ((mfactorial) ((mtimes) 2 *index)) 1027 -1)) 1028 *index 0 $inf) 1029 sp2) 1030 1031(defprop %coth 1032 ((%sum) 1033 ((mtimes) ((mexpt) 4 *index) 1034 (($bern) ((mtimes) 2 *index)) 1035 ((mexpt) ((mfactorial) ((mtimes) 2 *index)) -1) 1036 ((mexpt) sp2var ((mplus) ((mtimes) 2 *index) -1))) 1037 *index 0 $inf) 1038 sp2) 1039 1040(defprop %sech 1041 ((%sum) 1042 ((mtimes) (($euler) ((mtimes) 2 *index)) 1043 ((mexpt) ((mfactorial) ((mtimes) 2 *index)) -1) 1044 ((mexpt) sp2var ((mtimes) 2 *index))) 1045 *index 0 $inf) 1046 sp2) 1047 1048(defprop %csch 1049 ((%sum) 1050 ((mtimes) -2 ((mplus) ((mexpt) 2 ((mplus) ((mtimes) 2 *index) -1)) -1) 1051 ((mexpt) ((mfactorial) ((mtimes) *index 2)) -1) 1052 (($bern) ((mtimes) 2 *index)) 1053 ((mexpt) sp2var ((mplus) ((mtimes) 2 *index) -1))) 1054 *index 0 $inf) 1055 sp2) 1056 1057;;arc trigonometric function expansions 1058 1059(defprop %asin 1060 ((%sum) 1061 ((mtimes) ((%genfact) ((mplus) ((mtimes) 2 *index) -1) *index 2) 1062 ((mexpt) ((%genfact) ((mtimes) 2 *index) *index 2) -1) 1063 ((mexpt) ((mplus) ((mtimes) 2 *index) 1) -1) 1064 ((mexpt) sp2var ((mplus) ((mtimes) 2 *index) 1))) 1065 *index 0 $inf) 1066 sp2) 1067 1068(defprop %atan 1069 ((%sum) 1070 ((mtimes) ((mexpt) -1 *index) 1071 ((mexpt) ((mplus) ((mtimes) 2 *index) 1) -1) 1072 ((mexpt) sp2var ((mplus) ((mtimes) 2 *index) 1))) 1073 *index 0 $inf) 1074 sp2) 1075