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 sin) 14 15;;; Reference: J. Moses, Symbolic Integration, MIT-LCS-TR-047, 12-1-1967. 16;;; http://www.lcs.mit.edu/publications/pubs/pdf/MIT-LCS-TR-047.pdf. 17;;;; 18;;;; Unfortunately, some important pages in the scan are all black. 19;;;; 20;;;; A version with the missing pages is available (2008-12-14) from 21;;;; http://www.softwarepreservation.org/projects/LISP/MIT 22 23(declare-top (special $radexpand $%e_to_numlog ans arcpart coef 24 aa powerlist *a* *b* k stack e w y expres arg var 25 *powerl* *c* *d* exp varlist genvar $liflag $opsubst)) 26 27(defvar *debug-integrate* nil 28 "Enable debugging for the integrator routines.") 29 30;; When T do not call the risch integrator. This flag can be set by the risch 31;; integrator to avoid endless loops when calling the integrator from risch. 32(defvar *in-risch-p* nil) 33 34(defmacro op (frob) 35 `(get ,frob 'operators)) 36 37;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 38 39;;; Predicate functions 40 41(declaim (inline varp)) 42(defun varp (x) 43 (alike1 x var)) 44 45(defun integerp1 (x) 46 "Returns 2*x if 2*x is an integer, else nil" 47 (integerp2 (mul2* 2 x))) 48 49(defun integerp2 (x) 50 "Returns x if x is an integer, else false" 51 (let (u) 52 (cond ((not (numberp x)) nil) 53 ((not (floatp x)) x) 54 ((prog2 (setq u (maxima-rationalize x)) 55 (equal (cdr u) 1)) (car u))))) 56 57;; This predicate is used with m2 pattern matcher. 58;; A rational expression in var. 59(defun rat8 (ex) 60 (cond ((or (varp ex) (freevar ex)) 61 t) 62 ((member (caar ex) '(mplus mtimes) :test #'eq) 63 (do ((u (cdr ex) (cdr u))) 64 ((null u) t) 65 (if (not (rat8 (car u))) 66 (return nil)))) 67 ((not (eq (caar ex) 'mexpt)) 68 nil) 69 ((integerp (caddr ex)) 70 (rat8 (cadr ex))))) 71 72(defun rat8prime (c) 73 (and (rat8 c) 74 (or (not (mnump c)) 75 (not (zerop1 c))))) 76 77(defun elem (a) 78 (cond ((freevar a) t) 79 ((atom a) nil) 80 ((m2 a expres) t) 81 (t (every #'elem (cdr a))))) 82 83(defun freevar (a) 84 (cond ((atom a) (not (eq a var))) 85 ((varp a) nil) 86 ((and (not (atom (car a))) 87 (member 'array (cdar a) :test #'eq)) 88 (cond ((freevar (cdr a)) t) 89 (t (merror "~&FREEVAR: variable of integration appeared in subscript.")))) 90 (t (and (freevar (car a)) (freevar (cdr a)))))) 91 92;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 93 94;; possibly a bug: For var = x and *d* =3, we have expand(?subst10(x^9 * (x+x^6))) --> x^5+x^4, but 95;; ?subst10(expand(x^9 * (x+x^6))) --> x^5+x^3. (Barton Willis) 96 97(defun subst10 (ex) 98 (cond ((atom ex) ex) 99 ((and (eq (caar ex) 'mexpt) (eq (cadr ex) var)) 100 (list '(mexpt) var (integerp2 (quotient (caddr ex) *d*)))) 101 (t (cons (remove 'simp (car ex)) 102 (mapcar #'(lambda (c) (subst10 c)) (cdr ex)))))) 103 104(defun rationalizer (x) 105 (let ((ex (simplify ($factor x)))) 106 (if (not (alike1 ex x)) ex))) 107 108;; Like FIND-IF, but calls FUNC on elements of SEQ in turn until one returns 109;; non-NIL. At that point, return the result (rather than the input, which is 110;; what you'd get from FIND-IF) 111 112(defun map-find (func seq) 113 (map nil (lambda (x) 114 (let ((result (funcall func x))) 115 (when result (return-from map-find result)))) 116 seq)) 117 118;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 119 120;;; Stage II of the Integrator 121;;; 122;;; Check if the problem can be transformed or solved by special methods. 123;;; 11 Methods are implemented by Moses, some more have been added. 124 125(defun intform (expres &aux w) 126 (cond ((freevar expres) nil) 127 ((atom expres) nil) 128 129 ;; Map the function intform over the arguments of a sum or a product 130 ((member (caar expres) '(mplus mtimes) :test #'eq) 131 (map-find #'intform (cdr expres))) 132 133 ((or (eq (caar expres) '%log) 134 (arcp (caar expres))) 135 (cond 136 ;; Method 9: Rational function times a log or arctric function 137 ((setq arg (m2 exp 138 `((mtimes) ((,(caar expres)) (b rat8)) 139 ((coefftt) (c rat8prime))))) 140 ;; Integrand is of the form R(x)*F(S(x)) where F is a log, or 141 ;; arctric function and R(x) and S(x) are rational functions. 142 (ratlog exp var (cons (cons 'a expres) arg))) 143 (t 144 (prog (y z) 145 (cond 146 ((setq y (intform (cadr expres))) (return y)) 147 148 ;; Method 10: Rational function times log(b*x+a) 149 ((and (eq (caar expres) '%log) 150 (setq z (m2-b*x+a (cadr expres))) 151 (setq y (m2 exp 152 '((mtimes) 153 ((coefftt) (c rat8)) 154 ((coefftt) (d elem)))))) 155 (return 156 (let ((a (cdr (assoc 'a z :test #'eq))) 157 (b (cdr (assoc 'b z :test #'eq))) 158 (c (cdr (assoc 'c y :test #'eq))) 159 (d (cdr (assoc 'd y :test #'eq))) 160 (newvar (gensym "intform"))) 161 ;; keep var from appearing in questions to user 162 (putprop newvar t 'internal) 163 ;; Substitute y = log(b*x+a) and integrate again 164 (substint 165 expres 166 newvar 167 (integrator 168 (muln 169 (list (maxima-substitute 170 `((mquotient) ((mplus) ((mexpt) $%e ,newvar) 171 ((mtimes) -1 ,a)) 172 ,b) 173 var 174 c) 175 `((mquotient) ((mexpt) $%e ,newvar) ,b) 176 (maxima-substitute newvar expres d)) 177 nil) 178 newvar))))) 179 (t (return nil))))))) 180 181 ;; We have a special function with an integral on the property list. 182 ;; After the integral property was defined for the trig functions, 183 ;; in rev 1.52, need to exclude trig functions here. 184 ((and (not (atom (car expres))) 185 (not (optrig (caar expres))) 186 (not (eq (caar expres) 'mexpt)) 187 (get (caar expres) 'integral)) 188 (when *debug-integrate* 189 (format t "~&INTFORM: found 'INTEGRAL on property list~%")) 190 (cond 191 ((setq arg 192 (m2 exp `((mtimes) ((,(caar expres)) (b rat8)) ((coefftt) (c rat8prime))))) 193 ;; A rational function times the special function. 194 ;; Integrate with the method integration-by-parts. 195 (partial-integration (cons (cons 'a expres) arg) var)) 196 ;; The method of integration-by-parts can not be applied. 197 ;; Maxima tries to get a clue for the argument of the function which 198 ;; allows a substitution for the argument. 199 ((intform (cadr expres))) 200 (t nil))) 201 202 ;; Method 6: Elementary function of trigonometric functions 203 ((optrig (caar expres)) 204 (cond ((not (setq w (m2-b*x+a (cadr expres)))) 205 (intform (cadr expres))) 206 (t 207 (prog2 208 (setq *powerl* t) 209 (monstertrig exp var (cadr expres)))))) 210 211 ((and (eq (caar expres) '%derivative) 212 (eq (caar exp) (caar expres)) 213 (checkderiv exp))) 214 215 ;; Stop intform if we have not a power function. 216 ((not (eq (caar expres) 'mexpt)) nil) 217 218 ;; Method 2: Substitution for an integral power 219 ((integerp (caddr expres)) (intform (cadr expres))) 220 221 ;; Method 1: Elementary function of exponentials 222 ((freevar (cadr expres)) 223 (cond ((setq w (m2-b*x+a (caddr expres))) 224 (superexpt exp var (cadr expres) w)) 225 ((intform (caddr expres))) 226 ((and (eq '$%e (cadr expres)) 227 (isinop (caddr expres) '%log)) 228 ;; Found something like exp(r*log(x)) 229 (let* (($%e_to_numlog t) 230 ($radexpand nil) ; do not simplify sqrt(x^2) -> abs(x) 231 (nexp (resimplify exp))) 232 (cond ((alike1 exp nexp) nil) 233 (t (integrator (setq exp nexp) var))))) 234 (t nil))) 235 236 ;; The base is not a rational function. Try to get a clue for the base. 237 ((not (rat8 (cadr expres))) 238 (intform (cadr expres))) 239 240 ;; Method 3: Substitution for a rational root 241 ((and (setq w (m2-ratrootform (cadr expres))) ; e*(a*x+b) / (c*x+d) 242 (denomfind (caddr expres))) ; expon is ratnum 243 (or (progn 244 (setq *powerl* t) 245 (ratroot exp var (cadr expres) w)) 246 (inte exp var))) 247 248 ;; Method 4: Binomial - Chebyschev 249 ((not (integerp1 (caddr expres))) ; 2*exponent not integer 250 (cond ((m2-chebyform exp) 251 (chebyf exp var)) 252 (t (intform (cadr expres))))) 253 254 ;; Method 5: Arctrigonometric substitution 255 ((setq w (m2-c*x^2+b*x+a (cadr expres))) ; sqrt(c*x^2+b*x+a) 256 #+nil 257 (format t "expres = sqrt(c*x^2+b*x+a)~%") 258 ;; I think this is method 5, arctrigonometric substitutions. 259 ;; (Moses, pg 80.) The integrand is of the form 260 ;; R(x,sqrt(c*x^2+b*x+a)). This method first eliminates the b 261 ;; term of the quadratic, and then uses an arctrig substitution. 262 (inte exp var)) 263 264 ;; Method 4: Binomial - Chebyschev 265 ((m2-chebyform exp ) 266 (chebyf exp var)) 267 268 ;; Expand expres. 269 ;; Substitute the expanded factor into the integrand and try again. 270 ((not (m2 (setq w ($expand (cadr expres))) 271 (cadr expres))) 272 (prog2 273 (setq exp (maxima-substitute w (cadr expres) exp)) 274 (intform (simplify (list '(mexpt) w (caddr expres)))))) 275 276 ;; Factor expres. 277 ;; Substitute the factored factor into the integrand and try again. 278 ((setq w (rationalizer (cadr expres))) 279 ;; The forms below used to have $radexpand set to $all. But I 280 ;; don't think we really want to do that here because that makes 281 ;; sqrt(x^2) become x, which might be totally wrong. This is one 282 ;; reason why we returned -4/3 for the 283 ;; integrate(sqrt(x+1/x-2),x,0,1). We were replacing 284 ;; sqrt((x-1)^2) with x - 1, which is totally wrong since 0 <= x 285 ;; <= 1. 286 (setq exp (let (($radexpand $radexpand)) 287 (maxima-substitute w (cadr expres) exp))) 288 (intform (let (($radexpand '$all)) 289 (simplify (list '(mexpt) w (caddr expres)))))))) 290 291;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 292 293(defun separc (ex) 294 (cond ((arcfuncp ex) (setq arcpart ex coef 1)) 295 ((and (consp ex) (eq (caar ex) 'mtimes)) 296 (arclist (cdr ex)) 297 (setq coef (cond ((null (cdr coef)) (car coef)) 298 (t (setq coef (cons (car ex) coef)))))))) 299 300(defun arclist (list) 301 (cond ((null list) nil) 302 ((and (arcfuncp (car list)) (null arcpart)) 303 (setq arcpart (car list)) (arclist (cdr list))) 304 (t (setq coef (cons (car list) coef)) 305 (arclist (cdr list))))) 306 307(defun arcfuncp (ex) 308 (and (not (atom ex)) 309 (or (arcp (caar ex)) 310 (eq (caar ex) '%log) ; Experimentally treat logs also. 311 (and (eq (caar ex) 'mexpt) 312 (integerp2 (caddr ex)) 313 (> (integerp2 (caddr ex)) 0) 314 (arcfuncp (cadr ex)))))) 315 316;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 317 318;;; Five pattern for the Integrator and other routines. 319 320;; This is matching the pattern e*(a*x+b)/(c*x+d), where 321;; a, b, c, d, and e are free of x, and x is the variable of integration. 322(defun m2-ratrootform (expr) 323 (m2 expr 324 `((mtimes) 325 ((coefftt) (e freevar)) 326 ((mplus) 327 ((coeffpt) (a freevar) (var varp)) 328 ((coeffpt) (b freevar))) 329 ((mexpt) 330 ((mplus) 331 ((coeffpt) (c freevar) (var varp)) 332 ((coeffpt) (d freevar))) 333 -1)))) 334 335;; This is for matching the pattern a*x^r1*(c1+c2*x^q)^r2. 336(defun m2-chebyform (expr) 337 (m2 expr 338 `((mtimes) 339 ((mexpt) (var varp) (r1 numberp)) 340 ((mexpt) 341 ((mplus) 342 ((mtimes) 343 ((coefftt) (c2 freevar)) 344 ((mexpt) (var varp) (q free1))) 345 ((coeffpp) (c1 freevar))) 346 (r2 numberp)) 347 ((coefftt) (a freevar))))) 348 349;; Pattern to match b*x + a 350(defun m2-b*x+a (expr) 351 (m2 expr 352 `((mplus) 353 ((coeffpt) (b freevar) (x varp)) 354 ((coeffpt) (a freevar))))) 355 356;; This is the pattern c*x^2 + b * x + a. 357(defun m2-c*x^2+b*x+a (expr) 358 (m2 expr 359 `((mplus) 360 ((coeffpt) (c freevar) ((mexpt) (x varp) 2)) 361 ((coeffpt) (b freevar) (x varp)) 362 ((coeffpt) (a freevar))))) 363 364;; This is the pattern (a*x+b)*(c*x+d) 365(defun m2-a*x+b/c*x+d (expr) 366 (m2 expr 367 `((mtimes) 368 ((mplus) 369 ((coeffpt) (a freevar) (var varp)) 370 ((coeffpt) (b freevar))) 371 ((mplus) 372 ((coeffpt) (c freevar) (var varp)) 373 ((coeffpt) (d freevar)))))) 374 375;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 376 377;;; This is the main integration routine. It is called from sinint. 378 379(defun integrator (exp var) 380 (prog (y arg *powerl* const *b* w arcpart coef integrand result) 381 (declare (special *integrator-level*)) 382 ;; Increment recursion counter 383 (incf *integrator-level*) 384 385 ;; Trivial case. exp is not a function of var. 386 (if (freevar exp) (return (mul2* exp var))) 387 388 ;; Remove constant factors 389 (setq w (partition exp var 1)) 390 (setq const (car w)) 391 (setq exp (cdr w)) 392 #+nil 393 (progn 394 (format t "w = ~A~%" w) 395 (format t "const = ~A~%" const) 396 (format t "exp = ~A~%" exp)) 397 398 (cond ;; First stage, Method I: Integrate a sum. 399 ((mplusp exp) 400 (return (mul2* const (integrate1 (cdr exp))))) 401 402 ;; Convert atan2(a,b) to atan(a/b) and try again. 403 ((setq w (isinop exp '$atan2)) 404 (setq exp 405 (maxima-substitute (take '(%atan) (div (cadr w) (caddr w))) 406 w 407 exp)) 408 (return (mul* const 409 (integrator exp var)))) 410 411 ;; First stage, Method II: Integrate sums. 412 ((and (not (atom exp)) 413 (eq (caar exp) '%sum)) 414 (return (mul2* const (intsum exp var)))) 415 416 ;; First stage, Method III: Try derivative-divides method. 417 ;; This is the workhorse that solves many integrals. 418 ((setq y (diffdiv exp var)) 419 (return (mul2* const y)))) 420 421 ;; At this point, we have EXP as a product of terms. Make Y a 422 ;; list of the terms of the product. 423 (setq y (cond ((mtimesp exp) 424 (cdr exp)) 425 (t 426 (list exp)))) 427 428 ;; Second stage: 429 ;; We're looking at each term of the product and check if we can 430 ;; apply one of the special methods. 431 loop 432 #+nil 433 (progn 434 (format t "car y =~%") 435 (maxima-display (car y))) 436 (cond ((rat8 (car y)) 437 #+nil 438 (format t "In loop, go skip~%") 439 (go skip)) 440 ((and (setq w (intform (car y))) 441 ;; Do not return a noun form as result at this point, because 442 ;; we would like to check for further special integrals. 443 ;; We store the result for later use. 444 (setq result w) 445 (not (isinop w '%integrate))) 446 #+nil 447 (format t "In loop, case intform~%") 448 (return (mul2* const w))) 449 (t 450 #+nil 451 (format t "In loop, go special~%") 452 ;; Store a possible partial result 453 (setq result w) 454 (go special))) 455 skip 456 (setq y (cdr y)) 457 (cond ((null y) 458 ;; Method 8: Rational functions 459 (return (mul2* const (cond ((setq y (powerlist exp var)) y) 460 (t (ratint exp var))))))) 461 (go loop) 462 463 special 464 ;; Third stage: Try more general methods 465 466 ;; SEPARC SETQS ARCPART AND COEF SUCH THAT 467 ;; COEF*ARCEXP=EXP WHERE ARCEXP IS OF THE FORM 468 ;; ARCFUNC^N AND COEF IS ITS ALGEBRAIC COEFFICIENT 469 (separc exp) 470 471 #+nil 472 (progn 473 (format t "arcpart = ~A~%" arcpart) 474 (format t "coef =~%") 475 (maxima-display coef)) 476 (cond ((and (not (null arcpart)) 477 (do ((stacklist stack (cdr stacklist))) 478 ((null stacklist) t) 479 (cond ((alike1 (car stacklist) coef) 480 (return nil)))) 481 (not (isinop (setq w (let ((stack (cons coef stack))) 482 (integrator coef var))) 483 '%integrate)) 484 (setq integrand (mul2 w (sdiff arcpart var))) 485 (do ((stacklist stack (cdr stacklist))) 486 ((null stacklist) t) 487 (cond ((alike1 (car stacklist) integrand) 488 (return nil)))) 489 (not (isinop 490 (setq y (let ((stack (cons integrand stack)) 491 (integ integrand)) 492 (integrator integ var))) 493 '%integrate))) 494 (return (add* (list '(mtimes) const w arcpart) 495 (list '(mtimes) -1 const y)))) 496 (t 497 (return 498 (mul* const 499 (cond ((setq y (scep exp var)) 500 (cond ((cddr y) 501 #+nil 502 (progn 503 (format t "cddr y =~%") 504 (maxima-display (cddr y))) 505 (integrator ($trigreduce exp) var)) 506 (t (sce-int (car y) (cadr y) var)))) 507 ;; I don't understand why we do this. This 508 ;; causes the stack overflow in Bug 509 ;; 1487703, because we keep expanding exp 510 ;; into a form that matches the original 511 ;; and therefore we loop forever. To 512 ;; break this we keep track how how many 513 ;; times we've tried this and give up 514 ;; after 4 (arbitrarily selected) times. 515 ((and (< *integrator-level* 4) 516 (not (alike1 exp (setq y ($expand exp))))) 517 #+nil 518 (progn 519 (format t "exp = ~A~%" exp) 520 (maxima-display exp) 521 (format t "y = ~A~%" y) 522 (maxima-display y) 523 (break)) 524 (integrator y var)) 525 ((and (not *powerl*) 526 (setq y (powerlist exp var))) 527 y) 528 ((and (not *in-risch-p*) ; Not called from rischint 529 (setq y (rischint exp var)) 530 ;; rischint has not found an integral but 531 ;; returns a noun form. Do not return that 532 ;; noun form as result at this point, but 533 ;; store it for later use. 534 (setq result y) 535 (not (isinop y '%integrate))) 536 y) 537 ((setq y (integrate-exp-special exp var)) 538 ;; Maxima found an integral for a power function 539 y) 540 (t 541 ;; Integrate-exp-special has not found an integral 542 ;; We look for a previous result obtained by 543 ;; intform or rischint. 544 (if result 545 result 546 (list '(%integrate) exp var)))))))))) 547 548(defun optrig (x) 549 (member x '(%sin %cos %sec %tan %csc %cot) :test #'eq)) 550 551;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 552 553;;; Stage I 554;;; Implementation of Method 1: Integrate a sum 555 556;;after finding a non-integrable summand usually better to pass rest to risch 557(defun integrate1 (exp) 558 (do ((terms exp (cdr terms)) (ans)) 559 ((null terms) (addn ans nil)) 560 (let ($liflag) ; don't gen li's for 561 (push (integrator (car terms) var) ans)) ; parts of integrand 562 (when (and (not *in-risch-p*) ; Not called from rischint 563 (not (free (car ans) '%integrate)) 564 (cdr terms)) 565 (return (addn (cons (rischint (cons '(mplus) terms) var) (cdr ans)) 566 nil))))) 567 568(defun scep (expr var &aux trigl exp) ; Product of SIN, COS, EXP 569 (and (mtimesp expr) ; of linear args. 570 (loop for fac in (cdr expr) do 571 (cond ((atom fac) (return nil)) 572 ((trig1 (car fac)) 573 (if (linearp (cadr fac) var) (push fac trigl) 574 (return nil))) 575 ((and (mexptp fac) 576 (eq (cadr fac) '$%e) 577 (linearp (caddr fac) var)) 578 ;; should be only one exponential factor 579 (setq exp fac)) 580 (t (return nil))) 581 finally (return (cons exp trigl))))) 582 583;; Integrates exponential * sin or cos, all with linear args. 584 585(defun sce-int (exp s-c var) ; EXP is non-trivial 586 (let* ((e-coef (car (islinear (caddr exp) var))) 587 (sc-coef (car (islinear (cadr s-c) var))) 588 (sc-arg (cadr s-c)) 589 (abs-val (add (power e-coef 2) (power sc-coef 2)))) 590 (if (zerop1 abs-val) 591 ;; The numerator is zero. Exponentialize the integrand and try again. 592 (integrator ($exponentialize (mul exp s-c)) var) 593 (mul (div exp abs-val) 594 (add (mul e-coef s-c) 595 (if (eq (caar s-c) '%sin) 596 (mul* (neg sc-coef) `((%cos) ,sc-arg)) 597 (mul* sc-coef `((%sin) ,sc-arg)))))))) 598 599(defun checkderiv (expr) 600 (checkderiv1 (cadr expr) (cddr expr) () )) 601 602;; CHECKDERIV1 gets called on the expression being differentiated, 603;; an alternating list of variables being differentiated with 604;; respect to and powers thereof, and a reversed list of the latter 605;; that have already been examined. It returns either the antiderivative 606;; or (), saying this derivative isn't wrt the variable of integration. 607 608(defun checkderiv1 (expr wrt old-wrt) 609 (cond ((varp (car wrt)) 610 (if (equal (cadr wrt) 1) ;Power = 1? 611 (if (null (cddr wrt)) ;single or partial 612 (if (null old-wrt) 613 expr ;single 614 `((%derivative), expr ;partial in old-wrt 615 ,.(nreverse old-wrt))) 616 `((%derivative) ,expr ;Partial, return rest 617 ,.(nreverse old-wrt) 618 ,@(cddr wrt))) 619 `((%derivative) ,expr ;Higher order, reduce order 620 ,.(nreverse old-wrt) 621 ,(car wrt) ,(add* (cadr wrt) -1) 622 ,@ (cddr wrt)))) 623 ((null (cddr wrt)) () ) ;Say it doesn't apply here 624 (t (checkderiv1 expr (cddr wrt) ;Else we check later terms 625 (list* (cadr wrt) (car wrt) old-wrt))))) 626 627(defun integrallookups (exp) 628 (let (form dummy-args real-args) 629 (cond 630 ((eq (caar exp) 'mqapply) 631 ;; Transform to functional form and try again. 632 ;; For example: 633 ;; ((MQAPPLY SIMP) (($PSI SIMP ARRAY) 1) $X) 634 ;; => (($PSI) 1 $X) 635 (integrallookups `((,(caaadr exp)) ,@(cdadr exp) ,@(cddr exp)))) 636 637 ;; Lookup algorithm for integral of a special function. 638 ;; The integral form is put on the property list, and can be a 639 ;; lisp function of the args. If the form is nil, or evaluates 640 ;; to nil, then return noun form unevaluated. 641 ((and (not (atom (car exp))) 642 (setq form (get (caar exp) 'integral)) 643 (setq dummy-args (car form)) 644 (setq real-args (cdr exp)) 645 ;; search through the args of exp and find the arg containing var 646 ;; look up the integral wrt this arg from form 647 (setq form 648 (do ((x real-args (cdr x)) 649 (y (cdr form) (cdr y))) 650 ((or (null x) (null y)) nil) 651 (if (not (freevar (car x))) (return (car y))))) 652 ;; If form is a function then evaluate it with actual args 653 (or (not (functionp form)) 654 (setq form (apply form real-args)))) 655 (when *debug-integrate* 656 (format t "~&INTEGRALLOOKUPS: Found integral ~A.~%" (caar exp))) 657 (substitutel real-args dummy-args form)) 658 659 ((eq (caar exp) 'mplus) 660 (muln (list '((rat simp) 1 2) exp exp) nil)) 661 662 (t nil)))) 663 664;; Define the antiderivatives of some elementary special functions. 665;; This may not be the best place for this definition, but it is close 666;; to the original code. 667;; Antiderivatives that depend on the logabs flag are defined further below. 668(defprop %log ((x) ((mplus) ((mtimes) x ((%log) x)) ((mtimes) -1 x))) integral) 669(defprop %sin ((x) ((mtimes) -1 ((%cos) x))) integral) 670(defprop %cos ((x) ((%sin) x)) integral) 671(defprop %sinh ((x) ((%cosh) x)) integral) 672(defprop %cosh ((x) ((%sinh) x)) integral) 673;; No need to take logabs into account for tanh(x), because cosh(x) is positive. 674(defprop %tanh ((x) ((%log) ((%cosh) x))) integral) 675(defprop %sech ((x) ((%atan) ((%sinh) x))) integral) 676(defprop %asin ((x) ((mplus) ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) x 2))) ((rat) 1 2)) ((mtimes) x ((%asin) x)))) integral) 677(defprop %acos ((x) ((mplus) ((mtimes) -1 ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) x 2))) ((rat) 1 2))) ((mtimes) x ((%acos) x)))) integral) 678(defprop %atan ((x) ((mplus) ((mtimes) x ((%atan) x)) ((mtimes) ((rat) -1 2) ((%log) ((mplus) 1 ((mexpt) x 2)))))) integral) 679(defprop %acsc ((x) ((mplus) ((mtimes) ((rat) -1 2) ((%log) ((mplus) 1 ((mtimes) -1 ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) x -2))) ((rat) 1 2)))))) ((mtimes) ((rat) 1 2) ((%log) ((mplus) 1 ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) x -2))) ((rat) 1 2))))) ((mtimes) x ((%acsc) x)))) integral) 680(defprop %asec ((x) ((mplus) ((mtimes) ((rat) 1 2) ((%log) ((mplus) 1 ((mtimes) -1 ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) x -2))) ((rat) 1 2)))))) ((mtimes) ((rat) -1 2) ((%log) ((mplus) 1 ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) x -2))) ((rat) 1 2))))) ((mtimes) x ((%asec) x)))) integral) 681(defprop %acot ((x) ((mplus) ((mtimes) x ((%acot) x)) ((mtimes) ((rat) 1 2) ((%log) ((mplus) 1 ((mexpt) x 2)))))) integral) 682(defprop %asinh ((x) ((mplus) ((mtimes) -1 ((mexpt) ((mplus) 1 ((mexpt) x 2)) ((rat) 1 2))) ((mtimes) x ((%asinh) x)))) integral) 683(defprop %acosh ((x) ((mplus) ((mtimes) -1 ((mexpt) ((mplus) -1 ((mexpt) x 2)) ((rat) 1 2))) ((mtimes) x ((%acosh) x)))) integral) 684(defprop %atanh ((x) ((mplus) ((mtimes) x ((%atanh) x)) ((mtimes) ((rat) 1 2) ((%log) ((mplus) 1 ((mtimes) -1 ((mexpt) x 2))))))) integral) 685(defprop %acsch ((x) ((mplus) ((mtimes) ((rat) -1 2) ((%log) ((mplus) -1 ((mexpt) ((mplus) 1 ((mexpt) x -2)) ((rat) 1 2))))) ((mtimes) ((rat) 1 2) ((%log) ((mplus) 1 ((mexpt) ((mplus) 1 ((mexpt) x -2)) ((rat) 1 2))))) ((mtimes) x ((%acsch) x)))) integral) 686(defprop %asech ((x) ((mplus) ((mtimes) -1 ((%atan) ((mexpt) ((mplus) -1 ((mexpt) x -2)) ((rat) 1 2)))) ((mtimes) x ((%asech) x)))) integral) 687(defprop %acoth ((x) ((mplus) ((mtimes) x ((%acoth) x)) ((mtimes) ((rat) 1 2) ((%log) ((mplus) 1 ((mtimes) -1 ((mexpt) x 2))))))) integral) 688 689;; Define a little helper function to be used in antiderivatives. 690;; Depending on the logabs flag, it either returns log(x) or log(abs(x)). 691(defun log-or-logabs (x) 692 (take '(%log) (if $logabs (take '(mabs) x) x))) 693 694;; Define the antiderivative of tan(x), taking logabs into account. 695(defun integrate-tan (x) 696 (log-or-logabs (take '(%sec) x))) 697(putprop '%tan `((x) ,#'integrate-tan) 'integral) 698 699;; ... the same for csc(x) ... 700(defun integrate-csc (x) 701 (mul -1 (log-or-logabs (add (take '(%csc) x) (take '(%cot) x))))) 702(putprop '%csc `((x) ,#'integrate-csc) 'integral) 703 704;; ... the same for sec(x) ... 705(defun integrate-sec (x) 706 (log-or-logabs (add (take '(%sec) x) (take '(%tan) x)))) 707(putprop '%sec `((x) ,#'integrate-sec) 'integral) 708 709;; ... the same for cot(x) ... 710(defun integrate-cot (x) 711 (log-or-logabs (take '(%sin) x))) 712(putprop '%cot `((x) ,#'integrate-cot) 'integral) 713 714;; ... the same for coth(x) ... 715(defun integrate-coth (x) 716 (log-or-logabs (take '(%sinh) x))) 717(putprop '%coth `((x) ,#'integrate-coth) 'integral) 718 719;; ... the same for csch(x) ... 720(defun integrate-csch (x) 721 (log-or-logabs (take '(%tanh) (mul '((rat simp) 1 2) x)))) 722(putprop '%csch `((x) ,#'integrate-csch) 'integral) 723 724;; integrate(x^n,x) = if n # -1 then x^(n+1)/(n+1) else log-or-logabs(x). 725(defun integrate-mexpt-1 (x n) 726 (let ((n-is-minus-one ($askequal n -1))) 727 (cond ((eq '$yes n-is-minus-one) 728 (log-or-logabs x)) 729 (t 730 (setq n (add n 1)) 731 (div (take '(mexpt) x n) n))))) 732 733;; integrate(a^x,x) = a^x/log(a). 734(defun integrate-mexpt-2 (a x) 735 (div (take '(mexpt) a x) (take '(%log) a))) 736 737(putprop 'mexpt `((x n) ,#'integrate-mexpt-1 ,#'integrate-mexpt-2) 'integral) 738 739(defun rat10 (ex) 740 (cond ((freevar ex) t) 741 ((varp ex) nil) 742 ((eq (caar ex) 'mexpt) 743 (if (varp (cadr ex)) 744 (if (integerp2 (caddr ex)) 745 (setq powerlist (cons (caddr ex) powerlist))) 746 (and (rat10 (cadr ex)) (rat10 (caddr ex))))) 747 ((member (caar ex) '(mplus mtimes) :test #'eq) 748 (do ((u (cdr ex) (cdr u))) ((null u) t) 749 (if (not (rat10 (car u))) (return nil)))))) 750 751(defun integrate5 (ex var) 752 (if (rat8 ex) 753 (ratint ex var) 754 (integrator ex var))) 755 756(defun denomfind (x) 757 (cond ((ratnump x) (caddr x)) 758 ((not (numberp x)) nil) 759 ((not (floatp x)) 1) 760 (t (cdr (maxima-rationalize x))))) 761 762;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 763;;; 764;;; Stage II 765;;; Implementation of Method 1: Elementary function of exponentials 766;;; 767;;; The following examples are integrated with this method: 768;;; 769;;; integrate(exp(x)/(2+3*exp(2*x)),x) 770;;; integrate(exp(x+1)/(1+exp(x)),x) 771;;; integrate(10^x*exp(x),x) 772 773(let ((bas nil) ; The common base. 774 (pow nil) ; The common power of the form b*x+a. The values are 775 ; stored in a list which is returned from m2. 776 (exptflag nil)) ; When T, the substitution is not possible. 777 778 (defun superexpt (exp var bas1 pow1) 779 (prog (y (new-var (gensym "NEW-VAR-"))) 780 (setq bas bas1 781 pow pow1 782 exptflag nil) 783 ;; Transform the integrand. At this point resimplify, because it is not 784 ;; guaranteed, that a correct simplified expression is returned. 785 ;; Use a new variable to prevent facts on the old variable to be wrongly used. 786 (setq y (resimplify (maxima-substitute new-var var (elemxpt exp)))) 787 (when exptflag (return nil)) 788 ;; Integrate the transformed integrand and substitute back. 789 (return 790 ($multthru 791 (substint (list '(mexpt) bas 792 (list '(mplus) (cdras 'a pow) 793 (list '(mtimes) (cdras 'b pow) var))) 794 new-var 795 (integrator (div y 796 (mul new-var 797 (cdras 'b pow) 798 (take '(%log) bas))) new-var)))))) 799 800 ;; Transform expressions like g^(b*x+a) to the common base bas and 801 ;; do the substitution y = bas^(b*x+a) in the expr. 802 (defun elemxpt (expr &aux w) 803 (cond ((freevar expr) expr) 804 ;; var is the base of a subexpression. The transformation fails. 805 ((atom expr) (setq exptflag t)) 806 ((not (eq (caar expr) 'mexpt)) 807 (cons (car expr) 808 (mapcar #'(lambda (c) (elemxpt c)) (cdr expr)))) 809 ((not (freevar (cadr expr))) 810 (list '(mexpt) 811 (elemxpt (cadr expr)) 812 (elemxpt (caddr expr)))) 813 ;; Transform the expression to the common base. 814 ((not (eq (cadr expr) bas)) 815 (elemxpt (list '(mexpt) 816 bas 817 (mul (power (take '(%log) bas) -1) 818 (take '(%log) (cadr expr)) 819 (caddr expr))))) 820 ;; The exponent must be linear in the variable of integration. 821 ((not (setq w (m2-b*x+a (caddr expr)))) 822 (list (car expr) bas (elemxpt (caddr expr)))) 823 ;; Do the substitution y = g^(b*x+a). 824 (t 825 (setq w (cons (cons 'bb (cdras 'b pow)) w)) 826 (setq w (cons (cons 'aa (cdras 'a pow)) w)) 827 (setq w (cons (cons 'bas bas) w)) 828 (subliss w '((mtimes) 829 ((mexpt) bas a) 830 ((mexpt) 831 bas 832 ((mquotient) 833 ((mtimes) -1 aa b) bb)) 834 ((mexpt) 835 x 836 ((mquotient) b bb))))))) 837) ; End of let 838 839;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 840;;; 841;;; Stage II 842;;; Implementation of Method 3: 843;;; Substitution for a rational root of a linear fraction of x 844;;; 845;;; This method is applicable when the integrand is of the form: 846;;; 847;;; / 848;;; [ a x + b n1/m1 a x + b n1/m2 849;;; I R(x, (-------) , (-------) , ...) dx 850;;; ] c x + d c x + d 851;;; / 852;;; 853;;; Substitute 854;;; 855;;; (1) t = ((a*x+b)/(c*x+d))^(1/k), or 856;;; 857;;; (2) x = (b-d*t^k)/(c*t^k-a) 858;;; 859;;; where k is the least common multiplier of m1, m2, ... and 860;;; 861;;; (3) dx = k*(a*d-b*c)*t^(k-1)/(a-c*t^k)^2 * dt 862;;; 863;;; First, the algorithm calls the routine RAT3 to collect the roots of the 864;;; form ((a*x+b)/(c*x+d))^(n/m) in the list *ROOTLIST*. 865;;; search for the least common multiplier of m1, m2, ... then the 866;;; substitutions (2) and (3) are done and the new problem is integrated. 867;;; As always, W is an alist which associates to the coefficients 868;;; a, b... (and to VAR) their values. 869 870(defvar *ratroot* nil) ; Expression of the form (a*x+b)/(c*x+d) 871(defvar *rootlist* nil) ; List of powers of the expression *ratroot*. 872 873(defun ratroot (exp var *ratroot* w) 874 (prog (*rootlist* k y w1) 875 ;; Check if the integrand has a chebyform, if so return the result. 876 (when (setq y (chebyf exp var)) (return y)) 877 ;; Check if the integrand has a suitably form and collect the roots 878 ;; in the global special variable *ROOTLIST*. 879 (unless (rat3 exp t) (return nil)) 880 ;; Get the least common multiplier of m1, m2, ... 881 (setq k (apply #'lcm *rootlist*)) 882 (setq w1 (cons (cons 'k k) w)) 883 ;; Substitute for the roots. 884 (setq y 885 (subst41 exp 886 (subliss w1 887 '((mquotient) 888 ((mplus) ((mtimes) b e) 889 ((mtimes) -1 d ((mexpt) var k))) 890 ((mplus) ((mtimes) c ((mexpt) var k)) 891 ((mtimes) -1 e a)))) 892 var)) 893 ;; Integrate the new problem. 894 (setq y 895 (integrator 896 (mul y 897 (subliss w1 898 '((mquotient) 899 ((mtimes) e 900 ((mplus) 901 ((mtimes) a d k 902 ((mexpt) var ((mplus) -1 k))) 903 ((mtimes) -1 904 ((mtimes) b c k 905 ((mexpt) var ((mplus) -1 k)))))) 906 ((mexpt) ((mplus) 907 ((mtimes) c ((mexpt) var k)) 908 ((mtimes) -1 a e)) 909 2)))) 910 var)) 911 ;; Substitute back and return the result. 912 (return (substint (power *ratroot* (power k -1)) var y)))) 913 914(defun rat3 (ex ind) 915 (cond ((freevar ex) t) 916 ((atom ex) ind) 917 ((member (caar ex) '(mtimes mplus) :test #'eq) 918 (do ((u (cdr ex) (cdr u))) 919 ((null u) t) 920 (if (not (rat3 (car u) ind)) 921 (return nil)))) 922 ((not (eq (caar ex) 'mexpt)) 923 (rat3 (car (margs ex)) t)) 924 ((freevar (cadr ex)) 925 (rat3 (caddr ex) t)) 926 ((integerp (caddr ex)) 927 (rat3 (cadr ex) ind)) 928 ((and (m2 (cadr ex) *ratroot*) 929 (denomfind (caddr ex))) 930 (setq *rootlist* (cons (denomfind (caddr ex)) *rootlist*))) 931 (t (rat3 (cadr ex) nil)))) 932 933(let ((rootform nil) ; Expression of the form x = (b*e-d*t^k)/(c*t^k-e*a). 934 (rootvar nil)) ; The variable we substitute for the root. 935 936 (defun subst4 (ex) 937 (cond ((freevar ex) ex) 938 ((atom ex) rootform) 939 ((not (eq (caar ex) 'mexpt)) 940 (mapcar #'(lambda (u) (subst4 u)) ex)) 941 ((m2 (cadr ex) *ratroot*) 942 (list (car ex) rootvar (integerp2 (timesk k (caddr ex))))) 943 (t (list (car ex) (subst4 (cadr ex)) (subst4 (caddr ex)))))) 944 945 (defun subst41 (exp a b) 946 (setq rootform a 947 rootvar b) 948 ;; At this point resimplify, because it is not guaranteed, that a correct 949 ;; simplified expression is returned. 950 (resimplify (subst4 exp))) 951) ; End of let 952 953;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 954 955;;; Stage II 956;;; Implementation of Method 4: Binomial Chebyschev 957 958;; exp = a*t^r1*(c1+c2*t^q)^r2, where var = t. 959;; 960;; G&S 2.202 has says this integral can be expressed by elementary 961;; functions ii: 962;; 963;; 1. q is an integer 964;; 2. (r1+1)/q is an integer 965;; 3. (r1+1)/q+r2 is an integer. 966;; 967;; I (rtoy) think that for this code to work, r1, r2, and q must be numbers. 968(defun chebyf (exp var) 969 (prog (r1 r2 d1 d2 n1 n2 w q) 970 ;; Return NIL if the expression doesn't match. 971 (when (not (setq w (m2-chebyform exp))) 972 (return nil)) 973 #+nil 974 (format t "w = ~A~%" w) 975 (when (zerop1 (cdr (assoc 'c1 w :test #'eq))) 976 ;; rtoy: Is it really possible to be in this routine with c1 = 977 ;; 0? 978 (return 979 (mul* 980 ;; This factor is locally constant as long as t and 981 ;; c2*t^q avoid log's branch cut. 982 (subliss w '((mtimes) a ((mexpt) var ((mtimes) -1 q r2)) 983 ((mexpt) ((mtimes) c2 ((mexpt) var q)) r2))) 984 (integrator 985 (subliss w '((mexpt) var ((mplus) r1 ((mtimes) q r2)))) var)))) 986 (setq q (cdr (assoc 'q w :test #'eq))) 987 ;; Reset parameters. a = a/q, r1 = (1 - q + r1)/q 988 (setq w 989 (list* (cons 'a (div* (cdr (assoc 'a w :test #'eq)) q)) 990 (cons 991 'r1 992 (div* (addn (list 1 (neg (simplify q)) (cdr (assoc 'r1 w :test #'eq))) nil) q)) 993 w)) 994 #+nil 995 (format t "new w = ~A~%" w) 996 (setq r1 (cdr (assoc 'r1 w :test #'eq)) 997 r2 (cdr (assoc 'r2 w :test #'eq))) 998 #+nil 999 (progn 1000 (format t "new r1 = ~A~%" r1) 1001 (format t "r2 = ~A~%" r2)) 1002 ;; Write r1 = d1/n1, r2 = d2/n2, if possible. Update w with 1003 ;; these values, if so. If we can't, give up. I (rtoy) think 1004 ;; this only happens if r1 or r2 can't be expressed as rational 1005 ;; numbers. Hence, r1 and r2 have to be numbers, not variables. 1006 (cond 1007 ((not (and (setq d1 (denomfind r1)) 1008 (setq d2 (denomfind r2)) 1009 (setq n1 (integerp2 (timesk r1 d1))) 1010 (setq n2 (integerp2 (timesk r2 d2))) 1011 (setq w (list* (cons 'd1 d1) (cons 'd2 d2) 1012 (cons 'n1 n1) (cons 'n2 n2) 1013 w)))) 1014 #+nil 1015 (progn 1016 (format t "cheby can't find one of d1,d2,n1,n2:~%") 1017 (format t " d1 = ~A~%" d1) 1018 (format t " d2 = ~A~%" d2) 1019 (format t " n1 = ~A~%" n1) 1020 (format t " n2 = ~A~%" n2)) 1021 (return nil)) 1022 ((and (integerp2 r1) (> r1 0)) 1023 #+nil (format t "integer r1 > 0~%") 1024 ;; (r1+q-1)/q is positive integer. 1025 ;; 1026 ;; I (rtoy) think we are using the substitution z=(c1+c2*t^q). 1027 ;; Maxima thinks the resulting integral should then be 1028 ;; 1029 ;; a/q*c2^(-r1/q-1/q)*integrate(z^r2*(z-c1)^(r1/q+1/q-1),z) 1030 ;; 1031 (return 1032 (substint 1033 (subliss w '((mplus) c1 ((mtimes) c2 ((mexpt) var q)))) 1034 var 1035 (integrator 1036 (expands (list (subliss w 1037 ;; a*t^r2*c2^(-r1-1) 1038 '((mtimes) 1039 a 1040 ((mexpt) var r2) 1041 ((mexpt) 1042 c2 1043 ((mtimes) 1044 -1 1045 ((mplus) r1 1)))))) 1046 (cdr 1047 ;; (t-c1)^r1 1048 (expandexpt (subliss w 1049 '((mplus) 1050 var 1051 ((mtimes) -1 c1))) 1052 r1))) 1053 var)))) 1054 ((integerp2 r2) 1055 #+nil (format t "integer r2~%") 1056 ;; I (rtoy) think this is using the substitution z = t^(q/d1). 1057 ;; 1058 ;; The integral (as maxima will tell us) becomes 1059 ;; 1060 ;; a*d1/q*integrate(z^(n1/q+d1/q-1)*(c1+c2*z^d1)^r2,z) 1061 ;; 1062 ;; But be careful because the variable A in the code is 1063 ;; actually a/q. 1064 (return 1065 (substint (subliss w '((mexpt) var ((mquotient) q d1))) 1066 var 1067 (ratint (simplify (subliss w 1068 '((mtimes) 1069 d1 a 1070 ((mexpt) 1071 var 1072 ((mplus) 1073 n1 d1 -1)) 1074 ((mexpt) 1075 ((mplus) 1076 ((mtimes) 1077 c2 1078 ((mexpt) 1079 var d1)) 1080 c1) 1081 r2)))) 1082 var)))) 1083 ((and (integerp2 r1) (< r1 0)) 1084 #+nil (format t "integer r1 < 0~%") 1085 ;; I (rtoy) think this is using the substitution 1086 ;; 1087 ;; z = (c1+c2*t^q)^(1/d2) 1088 ;; 1089 ;; With this substitution, maxima says the resulting integral 1090 ;; is 1091 ;; 1092 ;; a/q*c2^(-r1/q-1/q)*d2* 1093 ;; integrate(z^(n2+d2-1)*(z^d2-c1)^(r1/q+1/q-1),z) 1094 (return 1095 (substint (subliss w 1096 ;; (c1+c2*t^q)^(1/d2) 1097 '((mexpt) 1098 ((mplus) 1099 c1 1100 ((mtimes) c2 ((mexpt) var q))) 1101 ((mquotient) 1 d2))) 1102 var 1103 (ratint (simplify (subliss w 1104 ;; This is essentially 1105 ;; the integrand above, 1106 ;; except A and R1 here 1107 ;; are not the same as 1108 ;; derived above. 1109 '((mtimes) 1110 a d2 1111 ((mexpt) 1112 c2 1113 ((mtimes) 1114 -1 1115 ((mplus) 1116 r1 1))) 1117 ((mexpt) 1118 var 1119 ((mplus) 1120 n2 d2 -1)) 1121 ((mexpt) 1122 ((mplus) 1123 ((mexpt) 1124 var d2) 1125 ((mtimes) -1 c1)) 1126 r1)))) 1127 var)))) 1128 ((integerp2 (add* r1 r2)) 1129 #+nil (format t "integer r1+r2~%") 1130 ;; If we're here, (r1-q+1)/q+r2 is an integer. 1131 ;; 1132 ;; I (rtoy) think this is using the substitution 1133 ;; 1134 ;; z = ((c1+c2*t^q)/t^q)^(1/d1) 1135 ;; 1136 ;; With this substitution, maxima says the resulting integral 1137 ;; is 1138 ;; 1139 ;; a*d2/q*c1^(r2+r1/q+1/q)* 1140 ;; integrate(z^(d2*r2+d2-1)*(z^d2-c2)^(-r2-r1/q-1/q-1),z) 1141 (return 1142 (substint (let (($radexpand '$all)) 1143 ;; Setting $radexpand to $all here gets rid of 1144 ;; ABS in the subtitution. I think that's ok in 1145 ;; this case. See Bug 1654183. 1146 (subliss w 1147 '((mexpt) 1148 ((mquotient) 1149 ((mplus) 1150 c1 1151 ((mtimes) c2 ((mexpt) var q))) 1152 ((mexpt) var q)) 1153 ((mquotient) 1 d1)))) 1154 var 1155 (ratint (simplify (subliss w 1156 '((mtimes) 1157 -1 a d1 1158 ((mexpt) 1159 c1 1160 ((mplus) 1161 r1 r2 1)) 1162 ((mexpt) 1163 var 1164 ((mplus) 1165 n2 d1 -1)) 1166 ((mexpt) 1167 ((mplus) 1168 ((mexpt) 1169 var d1) 1170 ((mtimes) 1171 -1 c2)) 1172 ((mtimes) 1173 -1 1174 ((mplus) 1175 r1 r2 1176 2)))))) 1177 var)))) 1178 (t (return (list '(%integrate) exp var)))))) 1179 1180(defun greaterratp (x1 x2) 1181 (cond ((and (numberp x1) (numberp x2)) 1182 (> x1 x2)) 1183 ((ratnump x1) 1184 (greaterratp (quotient (float (cadr x1)) 1185 (caddr x1)) 1186 x2)) 1187 ((ratnump x2) 1188 (greaterratp x1 1189 (quotient (float (cadr x2)) 1190 (caddr x2)))))) 1191 1192(defun trig1 (x) 1193 (member (car x) '(%sin %cos) :test #'eq)) 1194 1195(defun supertrig (exp) 1196 (declare (special *notsame* *trigarg*)) 1197 (cond ((freevar exp) t) 1198 ((atom exp) nil) 1199 ((member (caar exp) '(mplus mtimes) :test #'eq) 1200 (and (supertrig (cadr exp)) 1201 (or (null (cddr exp)) 1202 (supertrig (cons (car exp) 1203 (cddr exp)))))) 1204 ((eq (caar exp) 'mexpt) 1205 (and (supertrig (cadr exp)) 1206 (supertrig (caddr exp)))) 1207 ((eq (caar exp) '%log) 1208 (supertrig (cadr exp))) 1209 ((member (caar exp) 1210 '(%sin %cos %tan %sec %cot %csc) :test #'eq) 1211 (cond ((m2 (cadr exp) *trigarg*) t) 1212 ((m2-b*x+a (cadr exp)) 1213 (and (setq *notsame* t) nil)) 1214 (t (supertrig (cadr exp))))) 1215 (t (supertrig (cadr exp))))) 1216 1217(defun subst2s (ex pat) 1218 (cond ((null ex) nil) 1219 ((m2 ex pat) var) 1220 ((atom ex) ex) 1221 (t (cons (subst2s (car ex) pat) 1222 (subst2s (cdr ex) pat))))) 1223 1224;; Match (c*x+b), where c and b are free of x 1225(defun simple-trig-arg (exp) 1226 (m2 exp '((mplus) ((mtimes) 1227 ((coefftt) (c freevar)) 1228 ((coefftt) (v varp))) 1229 ((coeffpp) (b freevar))))) 1230 1231;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1232 1233;;; Stage II 1234;;; Implementation of Method 6: Elementary function of trigonometric functions 1235 1236(defun monstertrig (exp var *trigarg*) 1237 (declare (special *trigarg*)) 1238 (when (and (not (atom *trigarg*)) 1239 ;; Do not exute the following code when called from rischint. 1240 (not *in-risch-p*)) 1241 (let ((arg (simple-trig-arg *trigarg*))) 1242 (cond (arg 1243 ;; We have trig(c*x+b). Use the substitution y=c*x+b to 1244 ;; try to compute the integral. Why? Because x*sin(n*x) 1245 ;; takes longer and longer as n gets larger and larger. 1246 ;; This is caused by the Risch integrator. This is a 1247 ;; work-around for this issue. 1248 (let* ((c (cdras 'c arg)) 1249 (b (cdras 'b arg)) 1250 (new-var (gensym "NEW-VAR-")) 1251 (new-exp (maxima-substitute (div (sub new-var b) c) 1252 var exp))) 1253 (if (every-trigarg-alike new-exp new-var) 1254 ;; avoid endless recursion when more than one 1255 ;; trigarg exists or c is a float 1256 (return-from monstertrig 1257 (maxima-substitute 1258 *trigarg* 1259 new-var 1260 (div (integrator new-exp new-var) c)))))) 1261 (t 1262 (return-from monstertrig (rischint exp var)))))) 1263 (prog (*notsame* w a b y d) 1264 (declare (special *notsame*)) 1265 (cond 1266 ((supertrig exp) (go a)) 1267 ((null *notsame*) (return nil)) 1268 ;; Check for an expression like a*trig1(m*x)*trig2(n*x), 1269 ;; where trig1 and trig2 are sin or cos. 1270 ((not (setq y (m2 exp 1271 '((mtimes) 1272 ((coefftt) (a freevar)) 1273 (((b trig1)) 1274 ((mtimes) 1275 (x varp) 1276 ((coefftt) (m freevar)))) 1277 (((d trig1)) 1278 ((mtimes) 1279 (x varp) 1280 ((coefftt) (n freevar)))))))) 1281 (go b)) 1282; This check has been done with the pattern match. 1283; ((not (and (member (car (setq b (cdr (assoc 'b y :test #'eq)))) '(%sin %cos) :test #'eq) 1284; (member (car (setq d (cdr (assoc 'd y :test #'eq)))) '(%sin %cos) :test #'eq))) 1285; (return nil)) 1286 ((progn 1287 ;; The tests after this depend on values of b and d being 1288 ;; set. Set them here unconditionally, before doing the 1289 ;; tests. 1290 (setq b (cdras 'b y)) 1291 (setq d (cdras 'd y)) 1292 (and (eq (car b) '%sin) 1293 (eq (car d) '%sin))) 1294 ;; We have a*sin(m*x)*sin(n*x). 1295 ;; The integral is: a*(sin((m-n)*x)/(2*(m-n))-sin((m+n)*x)/(2*(m+n)) 1296 (return (subliss y 1297 '((mtimes) a 1298 ((mplus) 1299 ((mquotient) 1300 ((%sin) ((mtimes) ((mplus) m ((mtimes) -1 n)) x)) 1301 ((mtimes) 2 ((mplus) m ((mtimes) -1 n)))) 1302 ((mtimes) -1 1303 ((mquotient) 1304 ((%sin) ((mtimes) ((mplus) m n) x)) 1305 ((mtimes) 2 ((mplus) m n))))))))) 1306 ((and (eq (car b) '%cos) (eq (car d) '%cos)) 1307 ;; We have a*cos(m*x)*cos(n*x). 1308 ;; The integral is: a*(sin((m-n)*x)/(2*(m-n))+sin((m+n)*x)/(2*(m+n)) 1309 (return (subliss y 1310 '((mtimes) a 1311 ((mplus) 1312 ((mquotient) 1313 ((%sin) ((mtimes) ((mplus) m ((mtimes) -1 n)) x)) 1314 ((mtimes) 2 1315 ((mplus) m ((mtimes) -1 n)))) 1316 ((mquotient) 1317 ((%sin) ((mtimes) ((mplus) m n) x)) 1318 ((mtimes) 2 ((mplus) m n)))))))) 1319 ((or (and (eq (car b) '%cos) 1320 ;; The following (destructively!) swaps the values of 1321 ;; m and n if first trig term is sin. I (rtoy) don't 1322 ;; understand why this is needed. The formula 1323 ;; doesn't depend on that. 1324 (setq w (cdras 'm y )) 1325 (rplacd (assoc 'm y) (cdras 'n y)) 1326 (rplacd (assoc 'n y) w)) 1327 t) 1328 ;; We have a*cos(n*x)*sin(m*x). 1329 ;; The integral is: -a*(cos((m-n)*x)/(2*(m-n))+cos((m+n)*x)/(2*(m+n)) 1330 (return (subliss y 1331 '((mtimes) -1 a 1332 ((mplus) 1333 ((mquotient) 1334 ((%cos) ((mtimes) ((mplus) m ((mtimes) -1 n)) x)) 1335 ((mtimes) 2 ((mplus) m ((mtimes) -1 n)))) 1336 ((mquotient) 1337 ((%cos) ((mtimes) ((mplus) m n) x)) 1338 ((mtimes) 2 ((mplus) m n))))))))) 1339 b ;; At this point we have trig functions with different arguments, 1340 ;; but not a product of sin and cos. 1341 (cond ((not (setq y (prog2 1342 (setq *trigarg* var) 1343 (m2 exp 1344 '((mtimes) 1345 ((coefftt) (a freevar)) 1346 (((b trig1)) 1347 ((mtimes) 1348 (x varp) 1349 ((coefftt) (n integerp2)))) 1350 ((coefftt) (c supertrig))))))) 1351 (return nil))) 1352 ;; We have a product of trig functions: trig1(n*x)*trig2(y). 1353 ;; trig1 is sin or cos, where n is a numerical integer. trig2 is not a sin 1354 ;; or cos. The cos or sin function is expanded. 1355 (return 1356 (integrator 1357 ($expand 1358 (list '(mtimes) 1359 (cdras 'a y) ; constant factor 1360 (cdras 'c y) ; trig functions 1361 (cond ((eq (car (cdras 'b y)) '%cos) ; expand cos(n*x) 1362 (maxima-substitute var 1363 'x 1364 (supercosnx (cdras 'n y)))) 1365 (t ; expand sin(x*x) 1366 (maxima-substitute var 1367 'x 1368 (supersinx (cdras 'n y))))))) 1369 var)) 1370 a ;; A product of trig functions and all trig functions have the same 1371 ;; argument *trigarg*. Maxima substitutes *trigarg* with the variable var 1372 ;; of integration and calls trigint to integrate the new problem. 1373 (setq w (subst2s exp *trigarg*)) 1374 (setq b (cdras 'b (m2-b*x+a *trigarg*))) 1375 (setq a (substint *trigarg* var (trigint (div* w b) var))) 1376 (return (if (isinop a '%integrate) 1377 (list '(%integrate) exp var) 1378 a)))) 1379 1380(defun trig2 (x) 1381 (member (car x) '(%sin %cos %tan %cot %sec %csc) :test #'eq)) 1382 1383(defun supersinx (n) 1384 (let ((i (if (< n 0) -1 1))) 1385 ($expand (list '(mtimes) i (sinnx (timesk i n)))))) 1386 1387(defun supercosnx (n) 1388 ($expand (cosnx (timesk (if (< n 0) -1 1) n)))) 1389 1390(defun sinnx (n) 1391 (if (equal n 1) 1392 '((%sin) x) 1393 (list '(mplus) 1394 (list '(mtimes) '((%sin) x) (cosnx (1- n))) 1395 (list '(mtimes) '((%cos) x) (sinnx (1- n)))))) 1396 1397(defun cosnx (n) 1398 (if (equal n 1) 1399 '((%cos) x) 1400 (list '(mplus) 1401 (list '(mtimes) '((%cos) x) (cosnx (1- n))) 1402 (list '(mtimes) -1 '((%sin) x) (sinnx (1- n)))))) 1403 1404(defun poseven (x) 1405 (and (even x) (> x -1))) 1406 1407(defun trigfree (x) 1408 (if (atom x) 1409 (not (member x '(sin* cos* sec* tan*) :test #'eq)) 1410 (and (trigfree (car x)) (trigfree (cdr x))))) 1411 1412(defun rat1 (exp) 1413 (prog (*b1* *notsame*) 1414 (declare (special *yy* *b1* *notsame*)) 1415 (when (and (numberp exp) (zerop exp)) 1416 (return nil)) 1417 (setq *b1* (subst *b* 'b '((mexpt) b (n even)))) 1418 (return (prog2 1419 (setq *yy* (rats exp)) 1420 (cond ((not *notsame*) *yy*)))))) 1421 1422(defun rats (exp) 1423 (prog (y) 1424 (declare (special *notsame* *b1*)) 1425 (return 1426 (cond ((eq exp *a*) 'x) 1427 ((atom exp) 1428 (cond ((member exp '(sin* cos* sec* tan*) :test #'eq) 1429 (setq *notsame* t)) 1430 (t exp))) 1431 ((setq y (m2 exp *b1*)) 1432 (f3 y)) 1433 (t (cons (car exp) (mapcar #'(lambda (g) (rats g)) (cdr exp)))))))) 1434 1435(defun f3 (y) 1436 (maxima-substitute *c* 1437 'c 1438 (maxima-substitute (quotient (cdr (assoc 'n y :test #'eq)) 2) 1439 'n 1440 '((mexpt) 1441 ((mplus) 1442 1 1443 ((mtimes) 1444 c 1445 ((mexpt) x 2))) 1446 n)))) 1447 1448(defun odd1 (n) 1449 (declare (special *yz*)) 1450 (cond ((not (numberp n)) nil) 1451 ((not (equal (rem n 2) 0)) 1452 (setq *yz* 1453 (maxima-substitute *c* 1454 'c 1455 (list '(mexpt) 1456 '((mplus) 1 ((mtimes) c ((mexpt) x 2))) 1457 (quotient (1- n) 2))))) 1458 (t nil))) 1459 1460(defun subvar (x) 1461 (maxima-substitute var 'x x)) 1462 1463(defun subvardlg (x) 1464 (mapcar #'(lambda (m) 1465 (cons (maxima-substitute var 'x (car m)) (cdr m))) 1466 x)) 1467 1468;; This appears to be the implementation of Method 6, pp.82 in Moses' thesis. 1469 1470(defun trigint (exp var) 1471 (prog (y repl y1 y2 *yy* z m n *c* *yz* *a* *b* ) 1472 (declare (special *yy* *yz*)) 1473 ;; Transform trig(x) into trig* (for simplicity?) Convert cot to 1474 ;; tan and csc to sin. 1475 (setq y2 1476 (subliss (subvardlg '((((%sin) x) . sin*) 1477 (((%cos) x) . cos*) 1478 (((%tan) x) . tan*) 1479 (((%cot) x) . ((mexpt) tan* -1)) 1480 (((%sec) x) . sec*) 1481 (((%csc) x) . ((mexpt) sin* -1)))) 1482 exp)) 1483 1484 (when *debug-integrate* 1485 (format t "~& in TRIGINT:~%") 1486 (format t "~& : y2 = ~A~%" y2)) 1487 1488 ;; Now transform tan to sin/cos and sec to 1/cos. 1489 (setq y1 (setq y (subliss '((tan* . ((mtimes) sin* 1490 ((mexpt) cos* -1))) 1491 (sec* . ((mexpt) cos* -1))) 1492 y2))) 1493 1494 (when *debug-integrate* (format t "~& : y = ~A~%" y)) 1495 1496 (when (null (setq z 1497 (m2 y 1498 '((mtimes) 1499 ((coefftt) (b trigfree)) 1500 ((mexpt) sin* (m poseven)) 1501 ((mexpt) cos* (n poseven)))))) 1502 ;; Go if y is not of the form sin^m*cos^n for positive even m and n. 1503 (go l1)) 1504 1505 ;; Case III: 1506 ;; Handle the case of sin^m*cos^n, m, n both non-negative and even. 1507 1508 (setq m (cdras 'm z)) 1509 (setq n (cdras 'n z)) 1510 (setq *a* (integerp2 (* 0.5 (if (< m n) 1 -1) (+ n (* -1 m))))) 1511 (setq z (cons (cons 'a *a*) z)) 1512 (setq z (cons (cons 'x var) z)) 1513 1514 (when *debug-integrate* 1515 (format t "~& CASE III:~%") 1516 (format t "~& : m, n = ~A ~A~%" m n) 1517 (format t "~& : a = ~A~%" *a*) 1518 (format t "~& : z = ~A~%" z)) 1519 1520 ;; integrate(sin(y)^m*cos(y)^n,y) is transformed to the following form: 1521 ;; 1522 ;; m < n: integrate((sin(2*y)/2)^n*(1/2+1/2*cos(2*y)^((n-m)/2),y) 1523 ;; m >= n: integrate((sin(2*y)/2)^n*(1/2-1/2*cos(2*y)^((m-n)/2),y) 1524 (return 1525 (mul (cdras 'b z) 1526 (div 1 2) 1527 (substint 1528 (mul 2 var) 1529 var 1530 (integrator 1531 (cond ((< m n) 1532 (subliss z 1533 '((mtimes) 1534 ((mexpt) 1535 ((mtimes) ((rat simp) 1 2) ((%sin) x)) 1536 m) 1537 ((mexpt) 1538 ((mplus) 1539 ((rat simp) 1 2) 1540 ((mtimes) 1541 ((rat simp) 1 2) ((%cos) x))) a)))) 1542 (t 1543 (subliss z 1544 '((mtimes) 1545 ((mexpt) 1546 ((mtimes) ((rat simp) 1 2) ((%sin) x)) 1547 n) 1548 ((mexpt) 1549 ((mplus) 1550 ((rat simp) 1 2) 1551 ((mtimes) 1552 ((rat simp) -1 2) 1553 ((%cos) x))) a))))) 1554 var)))) 1555 l1 1556 ;; Case IV: 1557 ;; I think this is case IV, working on the expression in terms of 1558 ;; sin and cos. 1559 ;; 1560 ;; Elem(x) means constants, x, trig functions of x, log and 1561 ;; inverse trig functions of x, and which are closed under 1562 ;; addition, multiplication, exponentiation, and substitution. 1563 ;; 1564 ;; Elem(f(x)) is the same as Elem(x), but f(x) replaces x in the 1565 ;; definition. 1566 1567 (when *debug-integrate* (format t "~& Case IV:~%")) 1568 1569 (setq *c* -1) 1570 (setq *a* 'sin*) 1571 (setq *b* 'cos*) 1572 (when (and (m2 y '((coeffpt) (c rat1) ((mexpt) cos* (n odd1)))) 1573 (setq repl (list '(%sin) var))) 1574 ;; The case cos^(2*n+1)*Elem(cos^2,sin). Use the substitution z = sin. 1575 (go getout)) 1576 (setq *a* *b*) 1577 (setq *b* 'sin*) 1578 (when (and (m2 y '((coeffpt) (c rat1) ((mexpt) sin* (n odd1)))) 1579 (setq repl (list '(%cos) var))) 1580 ;; The case sin^(2*n+1)*Elem(sin^2,cos). Use the substitution z = cos. 1581 (go get3)) 1582 1583 ;; Case V: 1584 ;; Transform sin and cos to tan and sec to see if the integral is 1585 ;; of the form Elem(tan, sec^2). If so, use the substitution z = tan. 1586 1587 (when *debug-integrate* (format t "~& Case V:~%")) 1588 1589 (setq y (subliss '((sin* (mtimes) tan* ((mexpt) sec* -1)) 1590 (cos* (mexpt) sec* -1)) 1591 y2)) 1592 (setq *c* 1) 1593 (setq *a* 'tan*) 1594 (setq *b* 'sec*) 1595 (when (and (rat1 y) (setq repl (list '(%tan) var))) 1596 (go get1)) 1597 (setq *a* *b*) 1598 (setq *b* 'tan*) 1599 (when (and (m2 y '((coeffpt) (c rat1) ((mexpt) tan* (n odd1)))) 1600 (setq repl (list '(%sec) var))) 1601 (go getout)) 1602 (when (not (alike1 (setq repl ($expand exp)) exp)) 1603 (return (integrator repl var))) 1604 (setq y (subliss '((sin* (mtimes) 2 x 1605 ((mexpt) ((mplus) 1 ((mexpt) x 2)) -1)) 1606 (cos* (mtimes) 1607 ((mplus) 1 ((mtimes) -1 ((mexpt) x 2))) 1608 ((mexpt) ((mplus) 1 ((mexpt) x 2)) -1))) 1609 y1)) 1610 (setq y (list '(mtimes) 1611 y 1612 '((mtimes) 2 ((mexpt) ((mplus) 1 ((mexpt) x 2)) -1)))) 1613 (setq repl (subvar '((mquotient) ((%sin) x) ((mplus) 1 ((%cos) x))))) 1614 (go get2) 1615 get3 1616 (setq y (list '(mtimes) -1 *yy* *yz*)) 1617 (go get2) 1618 get1 1619 (setq y (list '(mtimes) '((mexpt) ((mplus) 1 ((mexpt) x 2)) -1) *yy*)) 1620 (go get2) 1621 getout 1622 (setq y (list '(mtimes) *yy* *yz*)) 1623 get2 1624 (when *debug-integrate* 1625 (format t "~& Call the INTEGRATOR with:~%") 1626 (format t "~& : y = ~A~%" y) 1627 (format t "~& : repl = ~A~%" repl)) 1628 ;; See Bug 2880797. We want atan(tan(x)) to simplify to x, so 1629 ;; set $triginverses to '$all. 1630 (return 1631 ;; Do not integrate for the global variable VAR, but substitute it. 1632 ;; This way possible assumptions on VAR are no longer present. The 1633 ;; algorithm of DEFINT depends on this behavior. See Bug 3085498. 1634 (let (($triginverses '$all) (newvar (gensym))) 1635 (substint repl 1636 newvar 1637 (integrator (maxima-substitute newvar 'x y) newvar)))))) 1638 1639(defmvar $integration_constant_counter 0) 1640(defmvar $integration_constant '$%c) 1641 1642;; This is the top level of the integrator 1643(defun sinint (exp var) 1644 ;; *integrator-level* is a recursion counter for INTEGRATOR. See 1645 ;; INTEGRATOR for more details. Initialize it here. 1646 (let ((*integrator-level* 0)) 1647 (declare (special *integrator-level*)) 1648 1649 ;; Sanity checks for variables 1650 (when (mnump var) 1651 (merror (intl:gettext "integrate: variable must not be a number; found: ~:M") var)) 1652 (when ($ratp var) (setf var (ratdisrep var))) 1653 (when ($ratp exp) (setf exp (ratdisrep exp))) 1654 1655 (cond 1656 ;; Distribute over lists and matrices 1657 ((mxorlistp exp) 1658 (cons (car exp) 1659 (mapcar #'(lambda (y) (sinint y var)) (cdr exp)))) 1660 1661 ;; The symbolic integration code doesn't really deal very well with 1662 ;; subscripted variables, so if we have one then replace occurrences of var 1663 ;; with an atomic gensym and recurse. 1664 ((and (not (atom var)) 1665 (member 'array (cdar var))) 1666 (let ((dummy-var (gensym))) 1667 (maxima-substitute var dummy-var 1668 (sinint (maxima-substitute dummy-var var exp) dummy-var)))) 1669 1670 ;; If exp is an equality, integrate both sides and add an integration 1671 ;; constant 1672 ((mequalp exp) 1673 (list (car exp) (sinint (cadr exp) var) 1674 (add (sinint (caddr exp) var) 1675 ($concat $integration_constant (incf $integration_constant_counter))))) 1676 1677 ;; If var is an atom which occurs as an operator in exp, then return a noun form. 1678 ((and (atom var) 1679 (isinop exp var)) 1680 (list '(%integrate) exp var)) 1681 1682 ((zerop1 exp) ;; special case because 0 will not pass sum-of-intsp test 1683 0) 1684 1685 ((let ((ans (simplify 1686 (let ($opsubst varlist genvar stack) 1687 (integrator exp var))))) 1688 (if (sum-of-intsp ans) 1689 (list '(%integrate) exp var) 1690 ans)))))) 1691 1692;; SUM-OF-INTSP 1693;; 1694;; This is a heuristic that SININT uses to work out whether the result from 1695;; INTEGRATOR is worth returning or whether just to return a noun form. If this 1696;; function returns T, then SININT will return a noun form. 1697;; 1698;; The logic, as I understand it (Rupert 01/2014): 1699;; 1700;; (1) If I integrate f(x) wrt x and get an atom other than x or 0, either 1701;; something's gone horribly wrong, or this is part of a larger 1702;; expression. In the latter case, we can return T here because hopefully 1703;; something else interesting will make the top-level return NIL. 1704;; 1705;; (2) If I get a sum, return a noun form if every one of the summands is no 1706;; better than what I started with. (Presumably this is where the name 1707;; came from) 1708;; 1709;; (3) If this is a noun form, it doesn't convey much information on its own, 1710;; so return T. 1711;; 1712;; (4) If this is a product, something interesting has probably happened. But 1713;; I still want a noun form if the result is like 2*'integrate(f(x),x), so 1714;; I'm only interested in terms in the product that are not free of 1715;; VAR. If one of those terms is worthy of report, that's great: return 1716;; NIL. Otherwise, return T if we saw at least two things (eg splitting an 1717;; integral into a product of two integrals) 1718;; 1719;; (5) If the result is free of VAR, we're in a similar position to (1). 1720;; 1721;; (6) Otherwise something interesting (and hopefully useful) has 1722;; happened. Return NIL to tell SININT to report it. 1723(defun sum-of-intsp (ans) 1724 (cond ((atom ans) 1725 ;; Result of integration should never be a constant other than zero. 1726 ;; If the result of integration is zero, it is either because: 1727 ;; 1) a subroutine inside integration failed and returned nil, 1728 ;; and (mul 0 nil) yielded 0, meaning that the result is wrong, or 1729 ;; 2) the original integrand was actually zero - this is handled 1730 ;; with a separate special case in sinint 1731 (not (eq ans var))) 1732 ((mplusp ans) (every #'sum-of-intsp (cdr ans))) 1733 ((eq (caar ans) '%integrate) t) 1734 ((mtimesp ans) 1735 (let ((int-factors 0)) 1736 (not (or (dolist (factor (cdr ans)) 1737 (unless (freeof var factor) 1738 (if (sum-of-intsp factor) 1739 (incf int-factors) 1740 (return t)))) 1741 (<= 2 int-factors))))) 1742 ((freeof var ans) t) 1743 (t nil))) 1744 1745;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1746 1747;;; Stage I 1748;;; Implementation of Method 2: Integrate a summation 1749 1750(defun intsum (form var) 1751 (prog (exp idx ll ul pair val) 1752 (setq exp (cadr form) 1753 idx (caddr form) 1754 ll (cadddr form) 1755 ul (car (cddddr form))) 1756 (if (or (not (atom var)) 1757 (not (free idx var)) 1758 (not (free ll var)) 1759 (not (free ul var))) 1760 (return (list '(%integrate) form var))) 1761 (setq pair (partition exp var 1)) 1762 (when (and (mexptp (cdr pair)) 1763 (eq (caddr pair) var)) 1764 (setq val (maxima-substitute ll idx (cadddr pair))) 1765 (cond ((equal val -1) 1766 (return (add (integrator (maxima-substitute ll idx exp) var) 1767 (intsum1 exp idx (add 1 ll) ul var)))) 1768 ((mlsp val -1) 1769 (return (list '(%integrate) form var))))) 1770 (return (intsum1 exp idx ll ul var)))) 1771 1772(defun intsum1 (exp idx ll ul var) 1773 (assume (list '(mgeqp) idx ll)) 1774 (if (not (eq ul '$inf)) 1775 (assume (list '(mgeqp) ul idx))) 1776 (simplifya (list '(%sum) (integrator exp var) idx ll ul) t)) 1777 1778(defun finds (x) 1779 (if (atom x) 1780 (member x '(%log %integrate %atan) :test #'eq) 1781 (or (finds (car x)) (finds (cdr x))))) 1782 1783;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1784 1785;;; Stage II 1786;;; Implementation of Method 9: 1787;;; Rational function times a log or arctric function 1788 1789;;; ratlog is called for an expression containing a log or arctrig function 1790;;; The integrand is like log(x)*f'(x). To obtain the result the technique of 1791;;; partial integration is applied: log(x)*f(x)-integrate(1/x*f(x),x) 1792 1793(defun ratlog (exp var form) 1794 (prog (b c d y z w) 1795 (setq y form) 1796 (setq b (cdr (assoc 'b y :test #'eq))) 1797 (setq c (cdr (assoc 'c y :test #'eq))) 1798 (setq y (integrator c var)) 1799 (when (finds y) (return nil)) 1800 (setq d (sdiff (cdr (assoc 'a form :test #'eq)) var)) 1801 1802 (setq z (integrator (mul2* y d) var)) 1803 (setq d (cdr (assoc 'a form :test #'eq))) 1804 (return (simplify (list '(mplus) 1805 (list '(mtimes) y d) 1806 (list '(mtimes) -1 z)))))) 1807 1808;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1809 1810;;; partial-integration is an extension of the algorithm of ratlog to support 1811;;; the technique of partial integration for more cases. The integrand 1812;;; is like g(x)*f'(x) and the result is g(x)*f(x)-integrate(g'(x)*f(x),x). 1813;;; 1814;;; Adding integrals properties for elementary functions led to infinite recursion 1815;;; with integrate(z*expintegral_shi(z),z). This was resolved by limiting the 1816;;; recursion depth. *integrator-level* needs to be at least 3 to solve 1817;;; o integrate(expintegral_ei(1/sqrt(x)),x) 1818;;; o integrate(sqrt(z)*expintegral_li(z),z) 1819;;; while a value of 4 causes testsuite regressions with 1820;;; o integrate(z*expintegral_shi(z),z) 1821(defun partial-integration (form var) 1822 (declare (special *integrator-level*)) 1823 (let ((g (cdr (assoc 'a form))) ; part g(x) 1824 (df (cdr (assoc 'c form))) ; part f'(x) 1825 (f nil)) 1826 (setq f (integrator df var)) ; integrate f'(x) wrt var 1827 (cond 1828 ((or (isinop f '%integrate) ; no result or 1829 (isinop f (caar g)) ; g in result 1830 (> *integrator-level* 3)) 1831 nil) ; we return nil 1832 (t 1833 ;; Build the result: g(x)*f(x)-integrate(g'(x)*f(x)) 1834 (add (mul f g) 1835 (mul -1 (integrator (mul f (sdiff g var)) var))))))) 1836 1837;; returns t if argument of every trig operation in y matches arg 1838(defun every-trigarg-alike (y arg) 1839 (cond ((atom y) t) 1840 ((optrig (caar y)) (alike1 arg (cadr y))) 1841 (t (every (lambda (exp) 1842 (every-trigarg-alike exp arg)) 1843 (cdr y))))) 1844 1845;; return argument of first trig operation encountered in y 1846(defun find-first-trigarg (y) 1847 (cond ((atom y) nil) 1848 ((optrig (caar y)) (cadr y)) 1849 (t (some (lambda (exp) 1850 (find-first-trigarg exp)) 1851 (cdr y))))) 1852 1853;; return constant factor that makes elements of alist match elements of blist 1854;; or nil if no match found 1855;; (we could replace this using rat package to divide alist and blist) 1856(defun matchsum (alist blist) 1857 (prog (r s *c* *d*) 1858 (setq s (m2 (car alist) ;; find coeff for first term of alist 1859 '((mtimes) 1860 ((coefftt) (a freevar)) 1861 ((coefftt) (c true))))) 1862 (setq *c* (cdr (assoc 'c s :test #'eq))) 1863 (cond ((not (setq r ;; find coeff for first term of blist 1864 (m2 (car blist) 1865 (cons '(mtimes) 1866 (cons '((coefftt) (b free1)) 1867 (cond ((mtimesp *c*) 1868 (cdr *c*)) 1869 (t (list *c*)))))))) 1870 (return nil))) 1871 (setq *d* (simplify (list '(mtimes) 1872 (subliss s 'a) 1873 (list '(mexpt) 1874 (subliss r 'b) 1875 -1)))) 1876 (cond ((m2 (cons '(mplus) alist) ;; check that all terms match 1877 (timesloop *d* blist)) 1878 (return *d*)) 1879 (t (return nil))))) 1880 1881(defun timesloop (a b) 1882 (cons '(mplus) (mapcar #'(lambda (c) (mul2* a c)) b))) 1883 1884(defun expands (aa b) 1885 (addn (mapcar #'(lambda (c) (timesloop c aa)) b) nil)) 1886 1887(defun powerlist (exp var) 1888 (prog (y *c* *d* powerlist *b*) 1889 (setq y (m2 exp 1890 '((mtimes) 1891 ((mexpt) (var varp) (c integerp2)) 1892 ((coefftt) (a freevar)) 1893 ((coefftt) (b true))))) 1894 (setq *b* (cdr (assoc 'b y :test #'eq))) 1895 (setq *c* (cdr (assoc 'c y :test #'eq))) 1896 (unless (rat10 *b*) (return nil)) 1897 (setq *d* (apply #'gcd (cons (1+ *c*) powerlist))) 1898 (when (or (eql 1 *d*) (zerop *d*)) (return nil)) 1899 (return 1900 (substint 1901 (list '(mexpt) var *d*) 1902 var 1903 (integrate5 (simplify (list '(mtimes) 1904 (power* *d* -1) 1905 (cdr (assoc 'a y :test #'eq)) 1906 (list '(mexpt) var (1- (quotient (1+ *c*) *d*))) 1907 (subst10 *b*))) 1908 var))))) 1909 1910;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1911 1912;;; Stage I 1913;;; Implementation of Method 3: Derivative-divides algorithm 1914 1915;; This is the derivative-divides algorithm of Moses. 1916;; 1917;; / 1918;; [ 1919;; Look for form I c * op(u(x)) * u'(x) dx 1920;; ] 1921;; / 1922;; 1923;; where: c is a constant 1924;; u(x) is an elementary expression in x 1925;; u'(x) is its derivative 1926;; op is an elementary operator: 1927;; - the indentity, or 1928;; - any function that can be integrated by INTEGRALLOOKUPS 1929;; 1930;; The method of solution, once the problem has been determined to 1931;; posses the form above, is to look up OP in a table and substitute 1932;; u(x) for each occurrence of x in the expression given in the table. 1933;; In other words, the method performs an implicit substitution y = u(x), 1934;; and obtains the integral of op(y)dy by a table look up. 1935;; 1936(defun diffdiv (exp var) 1937 (prog (y *a* x v *d* z w r) 1938 (cond ((and (mexptp exp) 1939 (mplusp (cadr exp)) 1940 (integerp (caddr exp)) 1941 (< (caddr exp) 6) 1942 (> (caddr exp) 0)) 1943 (return (integrator (expandexpt (cadr exp) (caddr exp)) var)))) 1944 1945 ;; If not a product, transform to a product with one term 1946 (setq exp (cond ((mtimesp exp) exp) (t (list '(mtimes) exp)))) 1947 1948 ;; Loop over the terms in exp 1949 (setq z (cdr exp)) 1950 a (setq y (car z)) 1951 1952 ;; This m2 pattern matches const*(exp/y) 1953 (setq r (list '(mplus) 1954 (cons '(coeffpt) 1955 (cons '(c free1) 1956 (remove y (cdr exp) :count 1))))) 1957 (cond 1958 ;; Case u(var) is the identity function. y is a term in exp. 1959 ;; Match if diff(y,var) == c*(exp/y). 1960 ;; This even works when y is a function with multiple args. 1961 ((setq w (m2 (sdiff y var) r)) 1962 (return (muln (list y y (power* (mul2* 2 (cdr (assoc 'c w :test #'eq))) -1)) nil)))) 1963 1964 ;; w is the arg in y. 1965 (let ((arg-freevar)) 1966 (setq w 1967 (cond 1968 ((or (atom y) (member (caar y) '(mplus mtimes) :test #'eq)) y) 1969 ;; Take the argument of a function with one value. 1970 ((= (length (cdr y)) 1) (cadr y)) 1971 ;; A function has multiple args, and exactly one arg depends on var 1972 ((= (count-if #'null (setq arg-freevar (mapcar #'freevar (cdr y)))) 1) 1973 (do ((args (cdr y) (cdr args)) 1974 (argf arg-freevar (cdr argf))) 1975 ((if (not (car argf)) (return (car args)))))) 1976 (t 0)))) 1977 1978 (cond 1979 ((setq w (cond ((and (setq x (sdiff w var)) 1980 (mplusp x) 1981 (setq *d* (remove y (cdr exp) :count 1)) 1982 (setq v (car *d*)) 1983 (mplusp v) 1984 (not (cdr *d*))) 1985 (cond ((setq *d* (matchsum (cdr x) (cdr v))) 1986 (list (cons 'c *d*))) 1987 (t nil))) 1988 (t (m2 x r)))) 1989 (return (cond ((null (setq x (integrallookups y))) nil) 1990 ((eq w t) x) 1991 (t (mul2* x (power* (cdr (assoc 'c w :test #'eq)) -1))))))) 1992 (setq z (cdr z)) 1993 (when (null z) (return nil)) 1994 (go a))) 1995 1996(defun subliss (alist expr) 1997 "Alist is an alist consisting of a variable (symbol) and its value. expr is 1998 an expression. For each entry in alist, substitute the corresponding 1999 value into expr." 2000 (let ((x expr)) 2001 (dolist (a alist x) 2002 (setq x (maxima-substitute (cdr a) (car a) x))))) 2003 2004(defun substint (x y expres) 2005 (if (and (not (atom expres)) (eq (caar expres) '%integrate)) 2006 (list (car expres) exp var) 2007 (substint1 (maxima-substitute x y expres)))) 2008 2009(defun substint1 (exp) 2010 (cond ((atom exp) exp) 2011 ((and (eq (caar exp) '%integrate) 2012 (null (cdddr exp)) 2013 (not (symbolp (caddr exp))) 2014 (not (free (caddr exp) var))) 2015 (simplify (list '(%integrate) 2016 (mul2 (cadr exp) (sdiff (caddr exp) var)) 2017 var))) 2018 (t (recur-apply #'substint1 exp)))) 2019 2020;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2021;;; 2022;:; Extension of the integrator for more integrals with power functions 2023;;; 2024;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2025 2026;;; Recognize (a^(c*(z^r)^p+d)^v 2027 2028(defun m2-exp-type-1a (expr) 2029 (m2 expr 2030 '((mexpt) 2031 ((mexpt) 2032 (a freevar0) 2033 ((mplus) 2034 ;; The order of the pattern is critical. If we change it, 2035 ;; we do not get the expected match. 2036 ((coeffpp) (d freevar)) 2037 ((coefft) (c freevar0) 2038 ((mexpt) 2039 ((mexpt) (z varp) (r freevar0)) 2040 (p freevar))))) 2041 (v freevar)))) 2042 2043;;; Recognize z^v*a^(b*z^r+d) 2044 2045(defun m2-exp-type-2 (expr) 2046 (m2 expr 2047 '((mtimes) 2048 ((mexpt) (z varp) (v freevar0)) 2049 ((mexpt) 2050 (a freevar0) 2051 ((mplus) 2052 ((coeffpp) (d freevar)) 2053 ((coefft) (b freevar0) ((mexpt) (z varp) (r freevar0)))))))) 2054 2055;;; Recognize z^v*%e^(a*z^r+b)^u 2056 2057(defun m2-exp-type-2-1 (expr) 2058 (m2 expr 2059 '((mtimes) 2060 ((mexpt) (z varp) (v freevar0)) 2061 ((mexpt) 2062 ((mexpt) 2063 $%e 2064 ((mplus) 2065 ((coeffpp) (b freevar)) 2066 ((coefft) (a freevar0) ((mexpt) (z varp) (r freevar0))))) 2067 (u freevar))))) 2068 2069;;; Recognize (a*z+b)^p*%e^(c*z+d) 2070 2071(defun m2-exp-type-3 (expr) 2072 (m2 expr 2073 '((mtimes) 2074 ((mexpt) 2075 ((mplus) 2076 ((coefft) (a freevar0) (z varp)) 2077 ((coeffpp) (b freevar))) 2078 (p freevar0)) 2079 ((mexpt) 2080 $%e 2081 ((mplus) 2082 ((coefft) (c freevar0) (z varp)) 2083 ((coeffpp) (d freevar))))))) 2084 2085;;; Recognize d^(a*z^2+b/z^2+c) 2086 2087(defun m2-exp-type-4 (expr) 2088 (m2 expr 2089 '((mexpt) 2090 (d freevar0) 2091 ((mplus) 2092 ((coefft) (a freevar0) ((mexpt) (z varp) 2)) 2093 ((coefft) (b freevar0) ((mexpt) (z varp) -2)) 2094 ((coeffpp) (c freevar)))))) 2095 2096;;; Recognize z^(2*n)*d^(a*z^2+b/z^2+c) 2097 2098(defun m2-exp-type-4-1 (expr) 2099 (m2 expr 2100 '((mtimes) 2101 ((mexpt) (z varp) (n freevar0)) 2102 ((mexpt) 2103 (d freevar0) 2104 ((mplus) 2105 ((coefft) (a freevar0) ((mexpt) (z varp) 2)) 2106 ((coefft) (b freevar0) ((mexpt) (z varp) -2)) 2107 ((coeffpp) (c freevar))))))) 2108 2109;;; Recognize z^n*d^(a*z^2+b*z+c) 2110 2111(defun m2-exp-type-5 (expr) 2112 (m2 expr 2113 '((mtimes) 2114 ((mexpt) (z varp) (n freevar)) 2115 ((mexpt) 2116 (d freevar0) 2117 ((mplus) 2118 ((coeffpt) (a freevar) ((mexpt) (z varp) 2)) 2119 ((coeffpt) (b freevar) (z varp)) 2120 ((coeffpp) (c freevar))))))) 2121 2122;;; Recognize z^n*(%e^(a*z^2+b*z+c))^u 2123 2124(defun m2-exp-type-5-1 (expr) 2125 (m2 expr 2126 '((mtimes) 2127 ((mexpt) (z varp) (n freevar0)) 2128 ((mexpt) 2129 ((mexpt) 2130 $%e 2131 ((mplus) 2132 ((coeffpp) (c freevar)) 2133 ((coefft) (a freevar0) ((mexpt) (z varp) 2)) 2134 ((coefft) (b freevar0) (z varp)))) 2135 (u freevar))))) 2136 2137;;; Recognize z^n*d^(a*sqrt(z)+b*z+c) 2138 2139(defun m2-exp-type-6 (expr) 2140 (m2 expr 2141 '((mtimes) 2142 ((mexpt) (z varp) (n freevar0)) 2143 ((mexpt) 2144 (d freevar0) 2145 ((mplus) 2146 ((coefft) (a freevar0) ((mexpt) (z varp) ((rat) 1 2))) 2147 ((coefft) (b freevar0) (z varp)) 2148 ((coeffpp) (c freevar))))))) 2149 2150;;; Recognize z^n*(%e^(a*sqrt(z)+b*z+c))^u 2151 2152(defun m2-exp-type-6-1 (expr) 2153 (m2 expr 2154 '((mtimes) 2155 ((mexpt) (z varp) (n freevar0)) 2156 ((mexpt) 2157 ((mexpt) 2158 $%e 2159 ((mplus) 2160 ((coeffpp) (c freevar)) 2161 ((coefft) (a freevar0) ((mexpt) (z varp) ((rat) 1 2))) 2162 ((coefft) (b freevar0) (z varp)))) 2163 (u freevar))))) 2164 2165;;; Recognize z^n*a^(b*z^r+e)*h^(c*z^r+g) 2166 2167(defun m2-exp-type-7 (expr) 2168 (m2 expr 2169 '((mtimes) 2170 ((mexpt) (z varp) (n freevar)) 2171 ((mexpt) 2172 (a freevar0) 2173 ((mplus) 2174 ((coefft) 2175 (b freevar0) 2176 ((mexpt) (z varp) (r freevar0))) 2177 ((coeffpp) (e freevar)))) 2178 ((mexpt) 2179 (h freevar0) 2180 ((mplus) 2181 ((coefft) 2182 (c freevar0) 2183 ((mexpt) (z varp) (r1 freevar0))) 2184 ((coeffpp) (g freevar))))))) 2185 2186;;; Recognize z^v*(%e^(b*z^r+e))^q*(%e^(c*z^r+g))^u 2187 2188(defun m2-exp-type-7-1 (expr) 2189 (m2 expr 2190 '((mtimes) 2191 ((mexpt) (z varp) (v freevar)) 2192 ((mexpt) 2193 ((mexpt) 2194 $%e 2195 ((mplus) 2196 ((coeffpp) (e freevar)) 2197 ((coefft) (b freevar0) ((mexpt) (z varp) (r freevar0))))) 2198 (q freevar)) 2199 ((mexpt) 2200 ((mexpt) 2201 $%e 2202 ((mplus) 2203 ((coeffpp) (g freevar)) 2204 ((coefft) (c freevar0) ((mexpt) (z varp) (r1 freevar0))))) 2205 (u freevar))))) 2206 2207;;; Recognize a^(b*sqrt(z)+d*z+e)*h^(c*sqrt(z)+f*z+g) 2208 2209(defun m2-exp-type-8 (expr) 2210 (m2 expr 2211 '((mtimes) 2212 ((mexpt) 2213 (a freevar0) 2214 ((mplus) 2215 ((coeffpt) (b freevar) ((mexpt) (z varp) ((rat) 1 2))) 2216 ((coeffpt) (d freevar) (z varp)) 2217 ((coeffpp) (e freevar)))) 2218 ((mexpt) 2219 (h freevar0) 2220 ((mplus) 2221 ((coeffpt) (c freevar) ((mexpt) (z varp) ((rat) 1 2))) 2222 ((coeffpt) (f freevar) (z varp)) 2223 ((coeffpp) (g freevar))))))) 2224 2225;;; Recognize (%e^(b*sqrt(z)+d*z+e))^u*(%e^(c*sqrt(z)+f*z+g))^v 2226 2227(defun m2-exp-type-8-1 (expr) 2228 (m2 expr 2229 '((mtimes) 2230 ((mexpt) 2231 ((mexpt) 2232 $%e 2233 ((mplus) 2234 ((coeffpp) (e freevar)) 2235 ((coeffpt) (b freevar) ((mexpt) (z varp) ((rat) 1 2))) 2236 ((coeffpt) (d freevar) (z varp)))) 2237 (u freevar)) 2238 ((mexpt) 2239 ((mexpt) 2240 $%e 2241 ((mplus) 2242 ((coeffpp) (g freevar)) 2243 ((coeffpt) (c freevar) ((mexpt) (z varp) ((rat) 1 2))) 2244 ((coeffpt) (f freevar) (z varp)))) 2245 (v freevar))))) 2246 2247;;; Recognize (%e^(b*z^r+e))^u*(%e^(c*z^r+g))^v 2248 2249(defun m2-exp-type-8-2 (expr) 2250 (m2 expr 2251 '((mtimes) 2252 ((mexpt) 2253 ((mexpt) 2254 $%e 2255 ((mplus) 2256 ((coeffpp) (e freevar)) 2257 ((coefft) (b freevar) ((mexpt) (z varp) (r freevar0))))) 2258 (u freevar)) 2259 ((mexpt) 2260 ((mexpt) 2261 $%e 2262 ((mplus) 2263 ((coeffpp) (g freevar)) 2264 ((coefft) (c freevar) ((mexpt) (z varp) (r1 freevar0))))) 2265 (v freevar))))) 2266 2267;;; Recognize z^n*a^(b*z^2+d*z+e)*h^(c*z^2+f*z+g) 2268 2269(defun m2-exp-type-9 (expr) 2270 (m2 expr 2271 '((mtimes) 2272 ((mexpt) (z varp) (n freevar)) 2273 ((mexpt) 2274 (a freevar0) 2275 ((mplus) 2276 ((coeffpt) (b freevar) ((mexpt) (z varp) 2)) 2277 ((coeffpt) (d freevar) (z varp)) 2278 ((coeffpp) (e freevar)))) 2279 ((mexpt) 2280 (h freevar0) 2281 ((mplus) 2282 ((coeffpt) (c freevar) ((mexpt) (z varp) 2)) 2283 ((coeffpt) (f freevar) (z varp)) 2284 ((coeffpp) (g freevar))))))) 2285 2286;;; Recognize z^n*(%e^(b*z^2+d*z+e))^q*(%e^(c*z^2+f*z+g))^u 2287 2288(defun m2-exp-type-9-1 (expr) 2289 (m2 expr 2290 '((mtimes) 2291 ((mexpt) (z varp) (n freevar)) 2292 ((mexpt) 2293 ((mexpt) 2294 $%e 2295 ((mplus) 2296 ((coeffpp) (e freevar)) 2297 ((coeffpt) (b freevar) ((mexpt) (z varp) 2)) 2298 ((coeffpt) (d freevar) (z varp)))) 2299 (q freevar)) 2300 ((mexpt) 2301 ((mexpt) 2302 $%e 2303 ((mplus) 2304 ((coeffpp) (g freevar)) 2305 ((coeffpt) (c freevar) ((mexpt) (z varp) 2)) 2306 ((coeffpt) (f freevar) (z varp)))) 2307 (u freevar))))) 2308 2309;;; Recognize z^n*a^(b*sqrt(z)+d*z+e)*h^(c*sqrt(z+)f*z+g) 2310 2311(defun m2-exp-type-10 (expr) 2312 (m2 expr 2313 '((mtimes) 2314 ((mexpt) (z varp) (n freevar)) 2315 ((mexpt) 2316 (a freevar0) 2317 ((mplus) 2318 ((coeffpt) (b freevar) ((mexpt) (z varp) ((rat) 1 2))) 2319 ((coeffpt) (d freevar) (z varp)) 2320 ((coeffpp) (e freevar)))) 2321 ((mexpt) 2322 (h freevar0) 2323 ((mplus) 2324 ((coeffpt) (c freevar) ((mexpt) (z varp) ((rat) 1 2))) 2325 ((coeffpt) (f freevar) (z varp)) 2326 ((coeffpp) (g freevar))))))) 2327 2328;;; Recognize z^n*(%e^(b*sqrt(z)+d*z+e))^q*(%e^(c*sqrt(z)+f*z+g))^u 2329 2330(defun m2-exp-type-10-1 (expr) 2331 (m2 expr 2332 '((mtimes) 2333 ((mexpt) (z varp) (n freevar)) 2334 ((mexpt) 2335 ((mexpt) 2336 $%e 2337 ((mplus) 2338 ((coeffpp) (e freevar)) 2339 ((coeffpt) (b freevar) ((mexpt) (z varp) ((rat) 1 2))) 2340 ((coeffpt) (d freevar) (z varp)))) 2341 (q freevar)) 2342 ((mexpt) 2343 ((mexpt) 2344 $%e 2345 ((mplus) 2346 ((coeffpp) (g freevar)) 2347 ((coeffpt) (c freevar) ((mexpt) (z varp) ((rat) 1 2))) 2348 ((coeffpt) (f freevar) (z varp)))) 2349 (u freevar))))) 2350 2351;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2352 2353(defun integrate-exp-special (expr var &aux w const) 2354 2355 ;; First factor the expression. 2356 (setq expr ($factor expr)) 2357 2358 ;; Remove constant factors. 2359 (setq w (partition expr var 1)) 2360 (setq const (car w)) 2361 (setq expr (cdr w)) 2362 2363 (schatchen-cond w 2364 ((m2-exp-type-1a (facsum-exponent expr)) 2365 (a c d r p v) 2366 (when *debug-integrate* 2367 (format t "~&Type 1a: (a^(c*(z^r)^p+d)^v : w = ~A~%" w)) 2368 2369 (mul -1 2370 const 2371 ;; 1/(p*r*(a^(c*v*(var^r)^p))) 2372 (inv (mul p r (power a (mul c v (power (power var r) p))))) 2373 var 2374 ;; (a^(d+c*(var^r)^p))^v 2375 (power (power a (add d (mul c (power (power var r) p)))) v) 2376 ;; gamma_incomplete(1/(p*r), -c*v*(var^r)^p*log(a)) 2377 (take '(%gamma_incomplete) 2378 (inv (mul p r)) 2379 (mul -1 c v (power (power var r) p) (take '(%log) a))) 2380 ;; (-c*v*(var^r)^p*log(a))^(-1/(p*r)) 2381 (power (mul -1 c v (power (power var r) p) (take '(%log) a)) 2382 (div -1 (mul p r))))) 2383 2384 ((m2-exp-type-2 (facsum-exponent expr)) 2385 (a b d v r) 2386 2387 (when *debug-integrate* 2388 (format t "~&Type 2: z^v*a^(b*z^r+d) : w = ~A~%" w)) 2389 2390 (mul 2391 const 2392 (div -1 r) 2393 (power a d) 2394 (power var (add v 1)) 2395 ($gamma_incomplete 2396 (div (add v 1) r) 2397 (mul -1 b (power var r) ($log a))) 2398 (power 2399 (mul -1 b (power var r) ($log a)) 2400 (mul -1 (div (add v 1) r))))) 2401 2402 ((m2-exp-type-2-1 (facsum-exponent expr)) 2403 (a b v r u) 2404 (when *debug-integrate* 2405 (format t "~&Type 2-1: z^v*(%e^(a*z^r+b))^u : w = ~A~%" w)) 2406 2407 (mul const 2408 -1 2409 (inv r) 2410 (power '$%e (mul -1 a u (power var r))) 2411 (power (power '$%e (add (mul a (power var r)) b)) u) 2412 (power var (add v 1)) 2413 (power (mul -1 a u (power var r)) (div (mul -1 (add v 1)) r)) 2414 (take '(%gamma_incomplete) 2415 (div (add v 1) r) 2416 (mul -1 a u (power var r))))) 2417 2418 ((m2-exp-type-3 (facsum-exponent expr)) 2419 (a b c d p) 2420 (when *debug-integrate* 2421 (format t "~&Type 3: (a*z+b)^p*%e^(c*z+d) : w = ~A~%" w)) 2422 (mul 2423 const 2424 (div -1 a) 2425 (power '$%e (sub d (div (mul b c) a))) 2426 (power (add b (mul a var)) (add p 1)) 2427 ($expintegral_e (mul -1 p) (mul (div -1 a) c (add b (mul a var)))))) 2428 2429 ((m2-exp-type-4 expr) 2430 (a b c d) 2431 (let (($trigsign nil)) ; Do not simplify erfc(-x) ! 2432 (when *debug-integrate* 2433 (format t "~&Type 4: d^(a*z^2+b/z^2+c) : w = ~A~%" w)) 2434 2435 (mul 2436 const 2437 (div 1 (mul 4 (power (mul -1 a ($log d)) (div 1 2)))) 2438 (mul 2439 (power d c) 2440 (power '$%pi (div 1 2)) 2441 (power '$%e 2442 (mul -2 2443 (power (mul -1 a ($log d)) (div 1 2)) 2444 (power (mul -1 b ($log d)) (div 1 2)))) 2445 (add 2446 ($erfc 2447 (add 2448 (div (power (mul -1 b ($log d)) (div 1 2)) var) 2449 (mul -1 var (power (mul -1 a ($log d)) (div 1 2))))) 2450 (mul -1 2451 (power '$%e 2452 (mul 4 2453 (power (mul -1 a ($log d)) (div 1 2)) 2454 (power (mul -1 b ($log d)) (div 1 2)))) 2455 ($erfc 2456 (add 2457 (mul var (power (mul -1 a ($log d)) (div 1 2))) 2458 (div (power (mul -1 b ($log d)) (div 1 2)) var))))))))) 2459 2460 ((and (m2-exp-type-4-1 expr) 2461 (poseven (cdras 'n w)) ; only for n a positive, even integer 2462 (symbolp (cdras 'a w))) ; a has to be a symbol 2463 (a b c d n) 2464 (let (($trigsign nil)) ; Do not simplify erfc(-x) ! 2465 2466 (when *debug-integrate* 2467 (format t "~&Type 4-1: z^(2*n)*d^(a*z^2+b/z^2+c) : w = ~A~%" w)) 2468 2469 (setq n (div n 2)) 2470 2471 (mul const 2472 (div 1 4) 2473 (power d c) 2474 (power '$%pi (div 1 2)) 2475 (simplify (list '(%derivative) 2476 (div 2477 (sub 2478 (mul 2479 (power ($log d) (mul -1 n)) 2480 (add 2481 (mul 2482 (power 2483 '$%e 2484 (mul -2 2485 (power (mul -1 a ($log d)) (div 1 2)) 2486 (power (mul -1 b ($log d)) (div 1 2)))) 2487 ($erfc 2488 (sub 2489 (div 2490 (power (mul -1 b ($log d)) (div 1 2)) 2491 var) 2492 (mul var (power (mul -1 ($log d)) (div 1 2)))))))) 2493 (mul 2494 (power 2495 '$%e 2496 (mul 2 2497 (power (mul -1 a ($log d)) (div 1 2)) 2498 (power (mul -1 b ($log d)) (div 1 2)))) 2499 ($erfc 2500 (add 2501 (power (mul -1 a ($log d)) (div 1 2)) 2502 (div (power (mul -1 b ($log d)) (div 1 2)) var))))) 2503 (power (mul -1 a ($log d)) (div 1 2))) 2504 a n))))) 2505 2506 ((and (m2-exp-type-5 (facsum-exponent expr)) 2507 (maxima-integerp (cdras 'n w)) 2508 (eq ($sign (cdras 'n w)) '$pos)) 2509 (a b c d n) 2510 2511 (when *debug-integrate* 2512 (format t "~&Type 5: z^n*d^(a*z^2+b*z+c) : w = ~A~%" w)) 2513 2514 (mul 2515 const 2516 (div -1 (mul 2 (power (mul a ($log d)) (div 1 2)))) 2517 (mul 2518 (power d (sub c (div (mul b b) (mul 4 a)))) 2519 (let ((index (gensumindex)) 2520 ($simpsum t)) 2521 (mfuncall '$sum 2522 (mul 2523 (power 2 (sub index n)) 2524 ($binomial n index) 2525 ($gamma_incomplete 2526 (div (add index 1) 2) 2527 (mul 2528 (div -1 (mul 4 a)) 2529 (power (add b (mul 2 a var)) 2) 2530 ($log d))) 2531 (power (mul a ($log d)) (mul -1 (add n (div 1 2)))) 2532 (power (mul -1 b ($log d)) (sub n index)) 2533 (power (mul (add b (mul 2 a var)) ($log d)) (add index 1)) 2534 (power 2535 (mul (div -1 a) (power (add b (mul 2 a var)) 2) ($log d)) 2536 (mul (div -1 2) (add index 1)))) 2537 index 0 n))))) 2538 2539 ((and (m2-exp-type-5-1 (facsum-exponent expr)) 2540 (maxima-integerp (cdras 'n w)) 2541 (eq ($sign (cdras 'n w)) '$pos)) 2542 (a b c u n) 2543 (when *debug-integrate* 2544 (format t "~&Type 5-1: z^n*(%e^(a*z^2+b*z+c))^u : w = ~A~%" w)) 2545 2546 (mul const 2547 (div -1 2) 2548 (power '$%e 2549 (add (mul -1 (div (mul b b u) (mul 4 a))) 2550 (mul -1 u (add (mul a var var) (mul b var))))) 2551 (power a (mul -1 (add n 1))) 2552 (power (power '$%e 2553 (add (mul a var var) (mul b var) c)) 2554 u) 2555 (let ((index (gensumindex)) 2556 ($simpsum t)) 2557 (dosum 2558 (mul (power 2 (sub index n)) 2559 (power (mul -1 b) (sub n index)) 2560 (power (add b (mul 2 a var)) (add index 1)) 2561 (power (div (mul -1 u (power (add b (mul 2 a var)) 2)) a) 2562 (mul (div -1 2) (add index 1))) 2563 (take '(%binomial) n index) 2564 (take '(%gamma_incomplete) 2565 (div (add index 1) 2) 2566 (div (mul -1 u (power (add b (mul 2 a var)) 2)) 2567 (mul 4 a)))) 2568 index 0 n t)))) 2569 2570 ((and (m2-exp-type-6 (facsum-exponent expr)) 2571 (maxima-integerp (cdras 'n w)) 2572 (eq ($sign (cdras 'n w)) '$pos)) 2573 (a b c d n) 2574 (when *debug-integrate* 2575 (format t "~&Type 6: z^n*d^(a*sqrt(z)+b*z+c) : w = ~A~%" w)) 2576 2577 (mul 2578 const 2579 (power 2 (mul -1 (add n 1))) 2580 (power d (sub c (div (mul a a) (mul 4 b)))) 2581 (power (mul b ($log d)) (mul -2 (add n 1))) 2582 (let ((index1 (gensumindex)) 2583 (index2 (gensumindex)) 2584 ($simpsum t)) 2585 (mfuncall '$sum 2586 (mfuncall '$sum 2587 (mul 2588 (power -1 (sub index1 index2)) 2589 (power 4 index1) 2590 ($binomial index1 index2) 2591 ($binomial n index1) 2592 ($log d) 2593 (power (mul a ($log d)) (sub (mul 2 n) (add index1 index2))) 2594 (power 2595 (mul (add a (mul 2 b (power var (div 1 2)))) ($log d)) 2596 (add index1 index2)) 2597 (power 2598 (mul 2599 (div -1 b) 2600 (power (add a (mul 2 b (power var (div 1 2)))) 2) 2601 ($log d)) 2602 (mul (div -1 2) (add index1 index2 1))) 2603 (add 2604 (mul 2 b 2605 (power 2606 (mul 2607 (div -1 b) 2608 (power (add a (mul 2 b (power var (div 1 2)))) 2) 2609 ($log d)) 2610 (div 1 2)) 2611 ($gamma_incomplete 2612 (div (add index1 index2 2) 2) 2613 (mul 2614 (div -1 (mul 4 b)) 2615 (power (add a (mul 2 b (power var (div 1 2)))) 2) 2616 ($log d)))) 2617 (mul a 2618 (add a (mul 2 b (power var (div 1 2)))) 2619 ($log d) 2620 ($gamma_incomplete 2621 (div (add index1 index2 1) 2) 2622 (mul 2623 (div -1 (mul 4 b)) 2624 (power (add a (mul 2 b (power var (div 1 2)))) 2) 2625 ($log d)))))) 2626 index2 0 index1) 2627 index1 0 n)))) 2628 2629 ((and (m2-exp-type-6-1 (facsum-exponent expr)) 2630 (maxima-integerp (cdras 'n w)) 2631 (eq ($sign (cdras 'n w)) '$pos)) 2632 (a b c u n) 2633 (when *debug-integrate* 2634 (format t "~&Type 6-1: z^n*(%e^(a*sqrt(z)+b*z+c))^u : w = ~A~%" w)) 2635 2636 (mul const 2637 (power 2 (mul -1 (add (mul 2 n) 1))) 2638 (power '$%e 2639 (add (div (mul -1 u a a) (mul 4 b)) 2640 (mul u (add (mul a (power var (div 1 2))) 2641 (mul b var) 2642 c)))) 2643 (power b (mul -2 (add n 1))) 2644 (power (power '$%e 2645 (add (mul a (power var (div 1 2))) 2646 (mul b var))) 2647 u) 2648 (let ((index1 (gensumindex)) 2649 (index2 (gensumindex)) 2650 ($simpsum t)) 2651 (dosum 2652 (dosum 2653 (mul (power -1 (sub index1 index2)) 2654 (power 4 index1) 2655 (power a (add (neg index2) (neg index1) (mul 2 n))) 2656 (power (add a (mul 2 b (power var (div 1 2)))) 2657 (add index1 index2)) 2658 (power (div (mul -1 u 2659 (power (add a 2660 (mul 2 2661 b 2662 (power var (div 1 2)))) 2663 2)) 2664 b) 2665 (mul (div -1 2) (add index1 index2 1))) 2666 (take '(%binomial) index1 index2) 2667 (take '(%binomial) n index1) 2668 (add (mul a 2669 (add a (mul 2 b (power var (div 1 2)))) 2670 (take '(%gamma_incomplete) 2671 (div (add index1 index2 1) 2) 2672 (div (mul -1 u 2673 (power (add a 2674 (mul 2 b 2675 (power var 2676 (div 1 2)))) 2677 2)) 2678 (mul 4 b)))) 2679 (mul (inv u) 2680 (power (div (mul -1 u 2681 (power (add a 2682 (mul 2 b 2683 (power var 2684 (div 1 2)))) 2685 2)) 2686 b) 2687 (div 1 2)) 2688 (mul 2 b) 2689 (take '(%gamma_incomplete) 2690 (div (add index1 index2 2) 2) 2691 (div (mul -1 u 2692 (power (add a 2693 (mul 2 b 2694 (power var (div 1 2)))) 2695 2)) 2696 (mul 4 b)))))) 2697 index2 0 index1 t) 2698 index1 0 n t)))) 2699 2700 ((and (m2-exp-type-7 (facsum-exponent expr)) 2701 (eq ($sign (sub (cdras 'r w) (cdras 'r1 w))) '$zero)) 2702 (a b c e g h r n) 2703 (when *debug-integrate* 2704 (format t "~&Type 7: z^n*a^(b*z^r+e)*h^(c*z^r+g) : w = ~A~%" w)) 2705 2706 (setq n (add n 1)) 2707 2708 (mul 2709 const 2710 (power var n) 2711 (div -1 r) 2712 (power a e) 2713 (power h g) 2714 (power 2715 (mul -1 2716 (power var r) 2717 (add (mul b ($log a)) (mul c ($log h)))) 2718 (div (mul -1 n) r)) 2719 ($gamma_incomplete 2720 (div n r) 2721 (mul -1 (power var r) (add (mul b ($log a)) (mul c ($log h))))))) 2722 2723 ((and (m2-exp-type-7-1 (facsum-exponent expr)) 2724 (eq ($sign (sub (cdras 'r w) (cdras 'r1 w))) '$zero)) 2725 (b c e g r v q u) 2726 (when *debug-integrate* 2727 (format t "~&Type 7-1: z^v*(%e^(b*z^r+e))^q*(%e^(c*z^r+g))^u : w = ~A~%" w)) 2728 2729 (mul const 2730 (div -1 r) 2731 (power '$%e (mul -1 (power var r) (add (mul b q) (mul c u)))) 2732 (power (power '$%e (add e (mul b (power var r)))) q) 2733 (power (power '$%e (add g (mul c (power var r)))) u) 2734 (power var (add v 1)) 2735 (power (mul -1 (power var r) (add (mul b q) (mul c u))) 2736 (div (mul -1 (add v 1)) r)) 2737 (take '(%gamma_incomplete) 2738 (div (add v 1) r) 2739 (mul -1 (power var r) (add (mul b q) (mul c u)))))) 2740 2741 ((m2-exp-type-8 (facsum-exponent expr)) 2742 (a b c d e f g h) 2743 (when *debug-integrate* 2744 (format t "~&Type 8: a^(b*sqrt(z)+d*z+e)*h^(c*sqrt(z)+f*z+g)") 2745 (format t "~& : w = ~A~%" w)) 2746 2747 (mul 2748 const 2749 (div 1 2) 2750 (power a e) 2751 (power h g) 2752 (add 2753 (mul 2 2754 (power a (add (mul b (power var (div 1 2))) (mul d var))) 2755 (power h (add (mul c (power var (div 1 2))) (mul f var))) 2756 (div 1 (add (mul d ($log a)) (mul f ($log h))))) 2757 (mul -1 2758 (power '$%pi (div 1 2)) 2759 (power '$%e 2760 (mul -1 2761 (div 2762 (power (add (mul b ($log a)) (mul c ($log h))) 2) 2763 (mul 4 (add (mul d ($log a)) (mul f ($log h))))))) 2764 ($erfi 2765 (div 2766 (add 2767 (mul b ($log a)) 2768 (mul c ($log h)) 2769 (mul 2 2770 (power var (div 1 2)) 2771 (add (mul d ($log a)) (mul f ($log h))))) 2772 (mul 2 2773 (power (add (mul d ($log a)) (mul f ($log h))) (div 1 2))))) 2774 (add (mul b ($log a)) (mul c ($log h))) 2775 (power (add (mul d ($log a)) (mul f ($log h))) (div -3 2)))))) 2776 2777 ((m2-exp-type-8-1 (facsum-exponent expr)) 2778 (b c d e f g u v) 2779 (when *debug-integrate* 2780 (format t "~&Type 8-1: (%e^(b*sqrt(z)+d*z+e))^u*(%e^(c*sqrt(z)+f*z+g))^v") 2781 (format t "~& : w = ~A~%" w)) 2782 2783 (mul const 2784 (div 1 2) 2785 (power (add (mul d u) (mul f v)) (div -3 2)) 2786 (mul (power '$%e 2787 (mul -1 2788 (power (add (mul b u) 2789 (mul 2 d u (power var (div 1 2))) 2790 (mul v (add c (mul 2 f (power var (div 1 2)))))) 2791 2) 2792 (inv (mul 4 (add (mul d u) (mul f v)))))) 2793 (power (power '$%e 2794 (add (mul b (power var (div 1 2))) 2795 e 2796 (mul d var))) 2797 u) 2798 (power (power '$%e 2799 (add (mul c (power var (div 1 2))) 2800 g 2801 (mul f var))) 2802 v) 2803 (add (mul 2 2804 (power '$%e 2805 (mul (power (add (mul b u) 2806 (mul 2 d u (power var (div 1 2))) 2807 (mul v (add c (mul 2 f (power var (div 1 2)))))) 2808 2) 2809 (inv (mul 4 (add (mul d u) (mul f v)))))) 2810 (power (add (mul d u) (mul f v)) (div 1 2))) 2811 (mul -1 2812 (power '$%pi (div 1 2)) 2813 (add (mul b u) (mul c v)) 2814 (take '(%erfi) 2815 (div (add (mul b u) 2816 (mul 2 d u (power var (div 1 2))) 2817 (mul c v) 2818 (mul 2 f v (power var (div 1 2)))) 2819 (mul 2 2820 (power (add (mul d u) (mul f v)) 2821 (div 1 2)))))))))) 2822 2823 ((and (m2-exp-type-8-2 (facsum-exponent expr)) 2824 (eq ($sign (sub (cdras 'r w) (cdras 'r1 w))) '$zero)) 2825 (b c e g r u v) 2826 (when *debug-integrate* 2827 (format t "~&Type 8-2: (%e^(b*z^r+e))^u*(%e^(c*z^r+g))^v") 2828 (format t "~& : w = ~A~%" w)) 2829 2830 (mul const 2831 -1 2832 (inv r) 2833 (power '$%e 2834 (mul -1 2835 (power var r) 2836 (add (mul b u) (mul c v)))) 2837 (power (power '$%e 2838 (add (power var r) e)) 2839 u) 2840 (power (power '$%e 2841 (add (power var r) g)) 2842 v) 2843 var 2844 (power (mul -1 2845 (power var r) 2846 (add (mul b u) (mul c v))) 2847 (div -1 r)) 2848 (take '(%gamma_incomplete) 2849 (div 1 r) 2850 (mul -1 (power var r) (add (mul b u) (mul c v)))))) 2851 2852 ((and (m2-exp-type-9 (facsum-exponent expr)) 2853 (maxima-integerp (cdras 'n w)) 2854 (eq ($sign (cdras 'n w)) '$pos) 2855 (or (not (eq ($sign (cdras 'b w)) '$zero)) 2856 (not (eq ($sign (cdras 'c w)) '$zero)))) 2857 (a b c d e f g h n) 2858 (when *debug-integrate* 2859 (format t "~&Type 9: z^n*a^(b*z^2+d*z+e)*h^(c*z^2+f*z+g)") 2860 (format t "~& : w = ~A~%" w)) 2861 2862 (mul 2863 const 2864 (div -1 2) 2865 (power a e) 2866 (power h g) 2867 (power '$%e 2868 (div 2869 (power (add (mul d ($log a)) (mul f ($log h))) 2) 2870 (mul -4 (add (mul b ($log a)) (mul c ($log h)))))) 2871 (power (add (mul b ($log a)) (mul c ($log h))) (mul -1 (add n 1))) 2872 (let ((index (gensumindex)) 2873 ($simpsum t)) 2874 (mfuncall '$sum 2875 (mul 2876 (power 2 (sub index n)) 2877 ($binomial n index) 2878 (power 2879 (add (mul -1 d ($log a)) (mul -1 f ($log h))) 2880 (sub n index)) 2881 (power 2882 (add 2883 (mul (add d (mul 2 b var)) ($log a)) 2884 (mul (add f (mul 2 c var)) ($log h))) 2885 (add index 1)) 2886 (power 2887 (mul -1 2888 (div 2889 (power 2890 (add 2891 (mul (add d (mul 2 b var)) ($log a)) 2892 (mul (add f (mul 2 c var)) ($log h))) 2893 2) 2894 (add (mul b ($log a)) (mul c ($log h))))) 2895 (div (add index 1) -2)) 2896 ($gamma_incomplete 2897 (div (add index 1) 2) 2898 (mul -1 2899 (div 2900 (power 2901 (add 2902 (mul (add d (mul 2 b var)) ($log a)) 2903 (mul (add f (mul 2 c var)) ($log h))) 2904 2) 2905 (mul 4 (add (mul b ($log a)) (mul c ($log h)))))))) 2906 index 0 n)))) 2907 2908 ((and (m2-exp-type-9-1 (facsum-exponent expr)) 2909 (maxima-integerp (cdras 'n w)) 2910 (eq ($sign (cdras 'n w)) '$pos) 2911 (or (not (eq ($sign (cdras 'b w)) '$zero)) 2912 (not (eq ($sign (cdras 'c w)) '$zero)))) 2913 (b c d e f g q u n) 2914 (when *debug-integrate* 2915 (format t "~&Type 9-1: z^n*(%e^(b*z^2+d*z+e))^q*(%e^(c*z^2+f*z+g))^u") 2916 (format t "~& : w = ~A~%" w)) 2917 2918 (mul const 2919 (div -1 2) 2920 (power (add (mul b q) (mul c u)) (div -1 2)) 2921 (power '$%e 2922 (add (div (power (add (mul d q) (mul f u)) 2) 2923 (mul -4 (add (mul b q) (mul c u)))) 2924 (mul -1 var 2925 (add (mul d q) 2926 (mul b q var) 2927 (mul f u) 2928 (mul c u var))))) 2929 (power (power '$%e 2930 (add e 2931 (mul var (add d (mul b var))))) 2932 q) 2933 (power (power '$%e 2934 (add g 2935 (mul var (add f (mul c var))))) 2936 u) 2937 (let ((index (gensumindex)) 2938 ($simpsum t)) 2939 (dosum 2940 (mul (power 2 (sub index n)) 2941 (power (add (mul b q) (mul c u)) (neg (add n (div 1 2)))) 2942 (power (add (neg (mul d q)) (neg (mul f u))) 2943 (sub n index)) 2944 (power (add (mul d q) 2945 (mul f u) 2946 (mul 2 var (add (mul b q) (mul c u)))) 2947 (add index 1)) 2948 (power (div (power (add (mul d q) 2949 (mul f u) 2950 (mul 2 2951 (add (mul b q) 2952 (mul c u)) 2953 var)) 2954 2) 2955 (neg (add (mul b q) (mul c u)))) 2956 (mul (div -1 2) (add index 1))) 2957 (take '(%binomial) n index) 2958 (take '(%gamma_incomplete) 2959 (div (add index 1) 2) 2960 (div (power (add (mul d q) 2961 (mul f u) 2962 (mul 2 2963 (add (mul b q) 2964 (mul c u)) 2965 var)) 2966 2) 2967 (mul -4 (add (mul b q) (mul c u)))))) 2968 index 0 n t)))) 2969 2970 ((and (m2-exp-type-10 (facsum-exponent expr)) 2971 (maxima-integerp (cdras 'n w)) 2972 (eq ($sign (cdras 'n w)) '$pos) 2973 (or (not (eq ($sign (cdras 'b w)) '$zero)) 2974 (not (eq ($sign (cdras 'c w)) '$zero)))) 2975 (a b c d e f g h n) 2976 (when *debug-integrate* 2977 (format t "~&Type 10: z^n*a^(b*sqrt(z)+d*z+e)*h^(c*sqrt(z)+f*z+g)") 2978 (format t "~& : w = ~A~%" w)) 2979 2980 (mul const 2981 (power 2 (add (mul -2 n) -1)) 2982 (power a e) 2983 (power h g) 2984 (power '$%e 2985 (div (power (add (mul b ($log a)) (mul c ($log h))) 2) 2986 (mul -4 (add (mul d ($log a)) (mul f ($log h)))))) 2987 (power (add (mul d ($log a)) (mul f ($log h))) (mul -2 (add n 1))) 2988 (let ((index1 (gensumindex)) 2989 (index2 (gensumindex)) 2990 ($simpsum t)) 2991 (dosum 2992 (dosum 2993 (mul (power -1 (sub index1 index2)) 2994 (power 4 index1) 2995 ($binomial index1 index2) 2996 ($binomial n index1) 2997 (power (add (mul b ($log a)) (mul c ($log h))) 2998 (sub (mul 2 n) (add index1 index2))) 2999 (power (add (mul b ($log a)) 3000 (mul c ($log h)) 3001 (mul 2 3002 (power var (div 1 2)) 3003 (add (mul d ($log a)) (mul f ($log h))))) 3004 (add index1 index2)) 3005 (power (mul -1 3006 (div (power (add (mul b ($log a)) 3007 (mul c ($log h)) 3008 (mul 2 3009 (power var (div 1 2)) 3010 (add (mul d ($log a)) 3011 (mul f ($log h))))) 3012 2) 3013 (add (mul d ($log a)) (mul f ($log h))))) 3014 (mul (div -1 2) (add index1 index2 1))) 3015 (add (mul ($gamma_incomplete (mul (div 1 2) 3016 (add index1 index2 1)) 3017 (mul (div -1 4) 3018 (div (power (add (mul b ($log a)) 3019 (mul c ($log h)) 3020 (mul 2 3021 (power var (div 1 2)) 3022 (add (mul d ($log a)) (mul f ($log h))))) 3023 2) 3024 (add (mul d ($log a)) (mul f ($log h)))))) 3025 (add (mul b ($log a)) (mul c ($log h))) 3026 (add (mul b ($log a)) 3027 (mul c ($log h)) 3028 (mul 2 3029 (power var (div 1 2)) 3030 (add (mul d ($log a)) (mul f ($log h)))))) 3031 (mul 2 3032 ($gamma_incomplete (mul (div 1 2) 3033 (add index1 index2 2)) 3034 (mul (div -1 4) 3035 (div (power (add (mul b ($log a)) 3036 (mul c ($log h)) 3037 (mul 2 3038 (power var (div 1 2)) 3039 (add (mul d ($log a)) 3040 (mul f ($log h))))) 3041 2) 3042 (add (mul d ($log a)) 3043 (mul f ($log h)))))) 3044 (add (mul d ($log a)) (mul f ($log h))) 3045 (power (mul -1 3046 (div (power (add (mul b ($log a)) 3047 (mul c ($log h)) 3048 (mul 2 3049 (power var (div 1 2)) 3050 (add (mul d ($log a)) 3051 (mul f ($log h))))) 3052 2) 3053 (add (mul d ($log a)) 3054 (mul f ($log h))))) 3055 (div 1 2))))) 3056 index2 0 index1 t) 3057 index1 0 n t)))) 3058 3059 ((and (m2-exp-type-10-1 (facsum-exponent expr)) 3060 (maxima-integerp (cdras 'n w)) 3061 (eq ($sign (cdras 'n w)) '$pos) 3062 (or (not (eq ($sign (cdras 'b w)) '$zero)) 3063 (not (eq ($sign (cdras 'c w)) '$zero)))) 3064 (b c d e f g q u n) 3065 (let ((bq+cu (add (mul b q) (mul c u))) 3066 (dq+fu (add (mul d q) (mul f u)))) 3067 (when *debug-integrate* 3068 (format t "~&Type 10-1: z^n*(%e^(b*sqrt(z)+d*z+e))^q*(%e^(c*sqrt(z)+f*z+g))^u") 3069 (format t "~& : w = ~A~%" w)) 3070 3071 (mul const 3072 (power 2 (mul -1 (add (mul 2 n) 1))) 3073 (power '$%e 3074 (add (div (mul -1 (power bq+cu 2)) (mul 4 dq+fu)) 3075 (mul -1 d var q) 3076 (mul -1 b (power var (div 1 2)) q) 3077 (mul -1 f var u) 3078 (mul -1 c (power var (div 1 2)) u))) 3079 (power (power '$%e 3080 (add (mul b (power var (div 1 2))) 3081 (mul d var) 3082 e)) 3083 q) 3084 (power (power '$%e 3085 (add (mul c (power var (div 1 2))) 3086 (mul f var) 3087 g)) 3088 u) 3089 (power dq+fu (mul -2 (add n 1))) 3090 (let ((index1 (gensumindex)) 3091 (index2 (gensumindex)) 3092 ($simpsum t)) 3093 (dosum 3094 (dosum 3095 (mul (power -1 (sub index1 index2)) 3096 (power 4 index1) 3097 (power bq+cu 3098 (add (neg index1) (neg index2) (mul 2 n))) 3099 (power (add bq+cu 3100 (mul 2 (power var (div 1 2)) dq+fu)) 3101 (add index1 index2)) 3102 (power (div (power (add bq+cu 3103 (mul 2 3104 (power var (div 1 2)) 3105 dq+fu)) 3106 2) 3107 (mul -1 dq+fu)) 3108 (mul (div -1 2) 3109 (add index1 index2 1))) 3110 (take '(%binomial) index1 index2) 3111 (take '(%binomial) n index1) 3112 (add (mul bq+cu 3113 (add bq+cu 3114 (mul 2 3115 (power var (div 1 2)) 3116 dq+fu)) 3117 (take '(%gamma_incomplete) 3118 (mul (div 1 2) 3119 (add index1 index2 1)) 3120 (div (power (add (mul b q) 3121 (mul c u) 3122 (mul 2 3123 (power var (div 1 2)) 3124 dq+fu)) 3125 2) 3126 (mul -4 3127 dq+fu)))) 3128 (mul 2 3129 (power (div (power (add bq+cu 3130 (mul 2 3131 (power var (div 1 2)) 3132 dq+fu)) 3133 2) 3134 (mul 1 dq+fu)) 3135 (div 1 2)) 3136 dq+fu 3137 (take '(%gamma_incomplete) 3138 (mul (div 1 2) 3139 (add index1 index2 2)) 3140 (div (power (add bq+cu 3141 (mul 2 3142 (power var (div 1 2)) 3143 dq+fu)) 3144 2) 3145 (mul -4 3146 dq+fu)))))) 3147 index2 0 index1 t) 3148 index1 0 n t))))))) 3149 3150;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3151 3152;;; Do a facsum for the exponent of power functions. 3153;;; This is necessary to integrate all general forms. The pattern matcher is 3154;;; not powerful enough to do the job. 3155 3156(defun facsum-exponent (expr) 3157 ;; Make sure that expr has the form ((mtimes) factor1 factor2 ...) 3158 (when (not (mtimesp expr)) (setq expr (list '(mtimes) expr))) 3159 (do ((result nil) 3160 (l (cdr expr) (cdr l))) 3161 ((null l) (cons (list 'mtimes) result)) 3162 (cond 3163 ((mexptp (car l)) 3164 ;; Found an power function. Factor the exponent with facsum. 3165 (let* ((fac (mfuncall '$facsum (caddr (car l)) var)) 3166 (num ($num fac)) 3167 (den ($denom fac))) 3168 (setq result 3169 (cons (cons (list 'mexpt) 3170 (cons (cadr (car l)) 3171 (if (equal 1 den) 3172 (list num) 3173 (list ($multthru (inv den) num))))) 3174 result)))) 3175 (t 3176 ;; Nothing to do. 3177 (setq result (cons (car l) result)))))) 3178 3179;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3180 3181