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 irinte) 14 15(load-macsyma-macros rzmac) 16 17(declare-top (special checkcoefsignlist var zerosigntest productcase)) 18 19(defun hasvar (exp) (not (freevar exp))) 20 21(defun zerp (a) (equal a 0)) 22 23(defun integerpfr (a) (if (not (maxima-integerp a)) (integerp1 a))) 24 25(defun nonzerp (a) (not (equal a 0))) 26 27(defun freevnz (a) (and (freevar a) (not (equal a 0)))) 28 29(defun inte (funct x) 30 (let ((checkcoefsignlist nil) 31 (*globalcareflag* nil) 32 ($radexpand t)) 33 (declare (special checkcoefsignlist *globalcareflag* $radexpand)) 34 (intir-ref funct x))) 35 36(defun intir-ref (fun x) 37 (prog (a) 38 (when (setq a (intir1 fun x)) (return a)) 39 (when (setq a (intir2 fun x)) (return a)) 40 (return (intir3 fun x)))) 41 42(defun intir1 (fun x) 43 (prog (assoclist e0 r0 e1 e2 r1 r2 d p) 44 (setq assoclist (factpow (specrepcheck fun) x)) 45 (setq e1 (cdras 'e1 assoclist) 46 e2 (cdras 'e2 assoclist)) 47 (cond ((null assoclist)(return nil))) 48 (setq d (cdras 'd assoclist) 49 p (cdras 'p assoclist) 50 e0 (cdras 'e0 assoclist) 51 r0 (cdras 'r0 assoclist) 52 r1 (cdras 'r1 assoclist) 53 r2 (cdras 'r2 assoclist)) 54 (cond ((floatp e0) 55 (setq e0 (rdis (ration1 e0))))) 56 (cond ((floatp e1) 57 (setq e1 (rdis (ration1 e1))))) 58 (cond ((floatp e2) 59 (setq e2 (rdis (ration1 e2))))) 60 (return (intir1-ref d p r0 e0 r1 e1 r2 e2 x)))) 61 62(defun intir2 (funct x) 63 (let ((res (intir funct x))) 64 (if res 65 res 66 (intirfactoroot funct x)))) 67 68(defun intir3 (exp x) 69 (prog ((assoclist (elliptquad exp x)) e f g r0) 70 (cond (assoclist 71 (setq e (cdras 'e assoclist) f (cdras 'f assoclist) 72 g (cdras 'g assoclist) r0 (cdras 'r0 assoclist)) 73 (assume `(($notequal) ,e 0)) 74 (return (intir3-r0test assoclist x e f g r0)))) 75 (return nil))) 76 77(defun intir3-r0test (assoclist x e f g r0) 78 (if (root+anything r0 x) 79 nil 80 (intir3-ref assoclist x e f g r0))) 81 82;; Handle integrals of the form d*p(x)*r1(x)^e1*r2(x)^e2*r0(x)^e0, 83;; where p(x) is a polynomial, e1 and e2 are both half an odd integer, 84;; and e3 is an integer. 85(defun intir1-ref (d p r0 e0 r1 e1 r2 e2 x) 86 (let ((nume1 (cadr e1)) ;; nume1 = 2*e1 87 (nume2 (cadr e2))) ;; nume2 = 2*e2 88 ;; I think what this does is try to rationalize r1(x)^e1 or 89 ;; r2(x)^e2 so we end up with a new integrand of the form 90 ;; d*p'(x)*r0'(x)^e0*ra(x)^ea, where p'(x) is a new polynomial 91 ;; obtained from rationaling one term and r0'(x) is a more 92 ;; complicated term. 93 (cond ((and (plusp nume1) (plusp nume2)) 94 (pp-intir1 d p r0 e0 r1 e1 r2 e2 x)) 95 ((and (minusp nume1) (minusp nume2)) 96 (mm-intir1 d p r0 e0 r1 e1 r2 e2 x)) 97 ((plusp nume1) 98 (pm-intir1 d p r0 e0 r1 e1 r2 e2 x)) 99 (t 100 (pm-intir1 d p r0 e0 r2 e2 r1 e1 x))))) 101 102(defun pp-intir1 (d p r0 e0 r1 e1 r2 e2 x) 103 (if (> (cadr e1) (cadr e2)) 104 (pp-intir1-exec d p r0 e0 r1 e1 r2 e2 x) 105 (pp-intir1-exec d p r0 e0 r2 e2 r1 e1 x))) 106 107;; Handle integrals of the form d*p(x)*r0(x)^e0*r1(x)^e1*r2(x)^e2 108;; where p(x) is a polynomial, e1 < 0, and e2 < 0 and are both half an 109;; odd integer, and e3 is an integer. 110(defun mm-intir1 (d p r0 e0 r1 e1 r2 e2 x) 111 (if (> (cadr e1) (cadr e2)) 112 (mm-intir1-exec d p r0 e0 r1 e1 r2 e2 x) 113 (mm-intir1-exec d p r0 e0 r2 e2 r1 e1 x))) 114 115;; Handle integrals of the form d*p(x)*r0(x)^e0*r1(x)^e1*r2(x)^e2 116;; where p(x) is a polynomial, e1 > 0, and e2 < 0 and are both half an 117;; odd integer, and e3 is an integer. 118;; 119(defun pm-intir1 (d p r0 e0 rofpos epos rofneg eneg x) 120 (let ((numepos (cadr epos)) ;; numepos = 2*epos = 2*e1 121 (numul-1eneg (mul -1 (cadr eneg)))) ;; numul-1eneg = -2*eneg = -2*e2 122 (cond ((> numepos numul-1eneg) 123 (mm-intir1 d (mul p (power rofpos (sub epos eneg))) 124 r0 e0 rofpos eneg rofneg eneg x)) 125 ((or (equal e0 0) (plusp e0)) 126 (pp-intir1 d (mul p (power rofneg (sub eneg epos))) 127 r0 e0 rofpos epos rofneg epos x)) 128 (t 129 (mm-intir1 d (mul p (power rofpos (sub epos eneg))) 130 r0 e0 rofpos eneg rofneg eneg x))))) 131 132(defun pp-intir1-exec (d p r0 e0 rofmax emax rofmin emin x) 133 (intir (mul d p (if (equal e0 0) 1 (power r0 e0)) 134 (power rofmax (add emax (mul -1 emin))) 135 (power ($expand (mul rofmax rofmin)) emin)) x)) 136 137;; Handle integrals of the form d*p(x)*r0(x)^e0*r1(x)^e1*r2(x)^e2 138;; where p(x) is a polynomial, e1 < 0, and e2 < 0 and are both half an 139;; odd integer, and e3 is an integer. And e2 > e1. 140(defun mm-intir1-exec (d p r0 e0 rofmin emin rofmax emax x) 141 (intir (mul d p 142 (if (equal e0 0) 1 143 (power r0 e0)) 144 (power rofmax (add emax (mul -1 emin))) 145 (power ($expand (mul rofmax rofmin)) emin)) x)) 146 147;; Integrating the form (e*x^2+f*x+g)^m*r0(x)^e0. 148 149(defun intir3-ref (assoclist x e f g r0) 150 (let ((signdisc (signdiscr e f g)) 151 (d (cdras 'd assoclist)) 152 (p (cdras 'p assoclist)) 153 (e0 (cdras 'e0 assoclist))) 154 (cond ((eq signdisc '$positive) 155 (pns-intir3 x e f g d p r0 e0)) 156 ((eq signdisc '$negative) 157 (ns-intir3 x e f g d p r0 e0)) 158 (t 159 (zs-intir3 x e f d p r0 e0))))) 160 161(defun root+anything (exp var) 162 (m2 exp '((mplus) 163 ((coeffpt) (c nonzerp) ((mexpt) (u hasvar) (v integerpfr))) 164 ((coeffpp) (c true))))) 165 166;; Handle d*p(x)/(e*x^2+f*x+g)*r0(x)^e0. We know that e*x^2+f*x+g has 167;; no repeated roots. Let D be the discriminant of this quadratic, 168;; sqrt(f^2-4*e*g). (If we're here, we already know that f^2-4*e*g > 169;; 0). Thus, we can factor this quadratic as 170;; (2*e*x+f-D)*(2*e*x+f+D)/(4*e). Thus, the original integrand 171;; becomes 172;; 173;; 4*e*d/(2*e*x+f-D)/(2*e*x+f+D)*p(x)*r0(x)^e0. 174;; 175;; We can separate this as a partial fraction to get 176;; 177;; (2*d*e/D/(2*e*x+f-D) - 2*d*e/D/(2*e*x+f+D))*p(x)*r0(x)^e0. 178;; 179;; So we have separated this into two "simpler" integrals. 180(defun pns-intir3 (x e f g d p r0 e0) 181 (let* ((discr (power (sub (mul f f) (mul 4 e g)) 1//2)) ;; Compute discriminant of 182 (p*r0^e0 (mul p (power r0 e0))) ;; quadratic: sqrt(f^2-4*e*g) 183 (2*e*x+f (add (mul 2 e x) f)) 184 (2*e*d*invdisc (mul 2 e d (inv discr)))) 185 (mul 2*e*d*invdisc 186 (sub (intir2 (mul (inv (sub 2*e*x+f discr)) p*r0^e0) x) 187 (intir2 (mul (inv (add 2*e*x+f discr)) p*r0^e0) x))))) 188 189;; Handle d*p(x)/(e*x^2+f*x+g)*r0(x)^e0. We know that e*x^2+f*x+g has 190;; repeated roots. 191;; 192(defun zs-intir3 (x e f d p r0 e0) 193 ;; Since e*x^2+f*x+g has repeated roots, it can be written as e*(x+r)^2. 194 ;; We easily see that r = f/(2*e), so rewrite the integrand as 195 ;; 196 ;; d*p(x)/e/(x+r)^2*r0(x)^e0. 197 (intir2 (mul d p (inv e) 198 (power (add x (div f (add e e))) -2) 199 (power r0 e0)) 200 x)) 201 202;; Handle d*p(x)/(e*x^2+f*x+g)*r0(x)^e0. We know that e*x^2+f*x+g has 203;; no real roots. 204;; 205;; G&R 2.252 shows how we can handle these integrals, but I'm too lazy 206;; to implement them right now, so return NIL to indicate we don't 207;; know what to do. But whatever it is we do, it's definitely not 208;; calling intir or intir2 like zs-intir3 or pns-intir3 do because 209;; they eventually call inti which only handles linear forms (e = 0.) 210;; We'll need to write custom versions. 211(defun ns-intir3 (xx ee fff gg dd pp r0 e0) 212 (declare (ignore xx ee fff gg dd pp r0 e0)) 213 nil) 214 215(defun cdras (a b) 216 (cdr (assoc a b :test #'equal))) 217 218(defun intir (funct x) 219 (inti funct x (jmaug (specrepcheck funct) x))) 220 221;; Integrate d*p(x)*(f*x+e)^m*(a*x^2+b*x+c)^n. p(x) is a polynomial, 222;; m is an integer, n is a number(?). a, b, c, e, and f are 223;; expressions independent of x (var). 224(defun inti (funct x assoclist) 225 (prog (met n expr f e #+nil denom) 226 (setq n (cdras 'n assoclist)) 227 (when (or (null assoclist) (maxima-integerp n)) 228 (return nil)) 229 (setq f (cdras 'f assoclist) 230 e (cdras 'e assoclist)) 231 ;; If e is 0 (or not given, we don't have to do the 232 ;; transformation. Just integrate it and return. 233 (when (or (equal e 0) (null e)) 234 (return (intira funct x))) 235 236 ;; (unless (numberp f) (go jump)) 237 ;; (when (plusp f) (go jump)) 238 ;; I (rtoy) think this is the case where f is a negative number. 239 ;; I think this is trying to convert f*x+e to -f*x-e to make the 240 ;; coefficient of x positive. And if I'm right, the code below 241 ;; isn't doing it correctly, except when m = 1 or m = -1. 242 ;; (setq denom (add (mul f x) e) 243 ;; f (mul -1 f) 244 ;; e (mul -1 e) 245 ;; funct (mul -1 (div (meval (mul denom funct)) 246 ;; (add (mul f x) e)))) 247 248 ;; jump 249 250 ;; Apply the linear substitution y = f*x+e. That is x = (y-e)/f. 251 ;; Then use INTIRA to integrate this. The integrand becomes 252 ;; something like p(y)*y^m and other terms. 253 (setq expr (intira (distrexpandroot 254 (cdr ($substitute 255 (mul (inv f) 256 (add (setq met (make-symbol (symbol-name '#:yannis))) 257 (mul -1 e))) 258 x funct))) 259 met)) 260 (setq expr (and expr (mul (inv f) expr))) 261 (return ($expand ($substitute (add (mul f x) e) met expr))))) 262 263(defun distrexpandroot (expr) 264 (if (null expr) 265 1 266 (mul (expandroot (car expr)) 267 (distrexpandroot (cdr expr))))) 268 269(defun expandroot (expr) 270 (if (atom expr) 271 expr 272 (if (and (eq (caar expr) 'mexpt) 273 (integerpfr (caddr expr))) 274 ($expand expr) 275 expr))) 276 277(defun intirfactoroot (expr x) 278 (let ((factors (timestest expr))) 279 (flet ((try-inti () 280 (let* ((e (distrfactor factors x)) 281 (alist (jmaug e x))) 282 (when alist (inti e x alist))))) 283 (or (try-inti) 284 (let ((*globalcareflag* t)) 285 (declare (special *globalcareflag*)) 286 (try-inti)))))) 287 288;; DISTRFACTOR 289;; 290;; Apply FACTOROOT to each element in the list, FACTORS. Accumulate the results 291;; by multiplication, associating to the right. 292(defun distrfactor (factors x) 293 (if (null factors) 294 1 295 (mul (factoroot (first factors) x) 296 (distrfactor (rest factors) x)))) 297 298;; FACTOROOT 299;; 300;; If EXPR is of the form A^B and is not free of VAR, use CAREFULFACTOR to try 301;; to factor it. Otherwise just return EXPR. 302(defun factoroot (expr var) 303 (if (and (mexptp expr) 304 (hasvar expr) 305 (integerpfr (caddr expr))) 306 (carefulfactor expr var) 307 expr)) 308 309;; CAREFULFACTOR 310;; 311;; Try to factor an expression of the form A^B. If *GLOBALCAREFLAG* is NIL, this 312;; is exactly the same as $FACTOR. Otherwise, use $FACTOR on (A/x)^B and then 313;; restore the missing x^B afterwards using RESTOREX. 314(defun carefulfactor (expr x) 315 (declare (special *globalcareflag*)) 316 (if (null *globalcareflag*) 317 ($factor expr) 318 (restorex ($factor (power (div (cadr expr) x) (caddr expr))) 319 x (caddr expr)))) 320 321;; RESTOREX 322;; 323;; Multiply EXPR by VAR^EXPT, trying to insert the factor of VAR inside terms in 324;; a product if possible. 325(defun restorex (expr var expt) 326 (cond 327 ((and (mexptp expr) 328 (equal expt (caddr expr))) 329 (power (restorex (cadr expr) var 1) (caddr expr))) 330 331 ((mtimesp expr) 332 (distrestorex (cdr expr) var expt)) 333 334 (t 335 (mul (power var expt) expr)))) 336 337;; DISTRESTOREX 338;; 339;; Takes a list of factors. Returns an expression that equals the product of 340;; these factors, but after multiplying one of them by VAR to try to multiply 341;; the entire product by VAR^EXPT. If it's not possible to multiply factors to 342;; do so, add a factor of VAR^EXPT to the end. 343(defun distrestorex (factors var expt) 344 (if (null factors) 345 (power var expt) 346 (let ((start (first factors))) 347 (if (and (mexptp start) 348 (equal expt (caddr start))) 349 350 (lmul 351 (cons (power ($expand (mul var (cadr start))) 352 (caddr start)) 353 (rest factors))) 354 355 (mul start (distrestorex (rest factors) var expt)))))) 356 357(defun timestest (expr) 358 (if (mtimesp expr) 359 (cdr expr) 360 (list expr))) 361 362;; Integrate a function of the form d*p(y)*y^m*(a*y^2+b*x+c)^n. 363;; n is half of an integer. 364(defun intira (funct x) 365 (prog (a b c *ec-1* d m n (assoclist (jmaug (specrepcheck funct) x)) 366 pluspowfo1 pluspowfo2 minuspowfo 367 polfact signn poszpowlist negpowlist) 368 (declare (special *ec-1*)) 369 (setq n (cdras 'n assoclist)) 370 ;; r12 1//2) 371 ;; (format t "n = ~A~%" n) 372 (when (or (null assoclist) 373 (maxima-integerp n)) 374 (return nil)) 375 (when (floatp n) 376 (setq n (rdis (ration1 n)))) 377 (setq d (cdras 'd assoclist)) 378 (when (equal d 0) (return 0)) 379 (setq c (cdras 'a assoclist)) 380 (when (equal c 0) (return nil)) 381 (setq m (cdras 'm assoclist) 382 polfact (cdras 'p assoclist) 383 ;; We know that n is of the form s/2, so just make n = s, 384 ;; and remember that the actual exponent needs to be 385 ;; divided by 2. 386 n (cadr n) 387 signn (checksigntm n) 388 *ec-1* (power c -1) 389 b (cdras 'b assoclist) 390 a (cdras 'c assoclist) 391 ;; pluspowfo1 = 1/2*(n-1), That is, the original exponent - 1/2. 392 pluspowfo1 (mul 1//2 (+ n -1)) 393 ;; minupowfo = 1/2*(n+1), that is, the original exponent + 1/2. 394 minuspowfo (mul 1//2 (+ n 1)) 395 ;; pluspowfo2 = -1/2*(n+1), that is, the negative of 1/2 396 ;; plus the original exponent. 397 pluspowfo2 (* -1 minuspowfo)) 398 (destructuring-bind (pos &optional neg) 399 (powercoeflist polfact m x) 400 (setf poszpowlist pos) 401 (setf negpowlist neg)) 402 403 #+nil (progn 404 (format t "n = ~A~%" n) 405 (format t "pluspowfo1 = ~A~%" pluspowfo1) 406 (format t "minuspowfo = ~A~%" minuspowfo) 407 (format t "pluspowfo2 = ~A~%" pluspowfo2)) 408 409 ;; I (rtoy) think powercoeflist computes p(x)/x^m as a Laurent 410 ;; series. POSZPOWLIST is a list of coefficients of the positive 411 ;; powers and NEGPOWLIST is a list of the negative coefficients. 412 (when (and (null negpowlist) 413 (not (null poszpowlist))) 414 ;; Only polynomial parts. 415 (when (eq signn '$positive) 416 (return (augmult (mul d 417 (nummnumn poszpowlist 418 pluspowfo1 419 minuspowfo c b a x))))) 420 (return (augmult (mul d 421 (nummdenn poszpowlist 422 pluspowfo2 c b a x))))) 423 (when (and (null poszpowlist) 424 (not (null negpowlist))) 425 ;; No polynomial parts 426 (when (eq signn '$positive) 427 (return (augmult (mul d 428 (denmnumn negpowlist minuspowfo c b a x))))) 429 (return (augmult (mul d 430 (denmdenn negpowlist pluspowfo2 c b a x))))) 431 (when (and (not (null negpowlist)) 432 (not (null poszpowlist))) 433 ;; Positive and negative powers. 434 (when (eq signn '$positive) 435 (return (add (augmult (mul d 436 (nummnumn poszpowlist 437 pluspowfo1 438 minuspowfo c b a x))) 439 (augmult (mul d 440 (denmnumn negpowlist 441 minuspowfo c b a x)))))) 442 (return (add (augmult (mul d 443 (nummdenn poszpowlist 444 pluspowfo2 c b a x))) 445 (augmult (mul d 446 (denmdenn negpowlist 447 pluspowfo2 c b a x)))))))) 448 449;; Match d*p(x)*(f*x+e)^m*(a*x^2+b*x+c)^n. p(x) is a polynomial, m is 450;; an integer, n is half of an integer. a, b, c, e, and f are 451;; expressions independent of x (var). 452(defun jmaug (exp var) 453 (m2 exp '((mtimes) 454 ((coefftt) (d freevar)) 455 ((coefftt) (p polyp)) 456 ((mexpt) ((mplus) ((coeffpt) (f freevar) (x varp)) 457 ((coeffpp)(e freevar))) 458 (m maxima-integerp)) 459 ((mexpt) ((mplus) 460 ((coeffpt) (a freevar) ((mexpt) (x varp) 2)) 461 ((coeffpt) (b freevar) (x varp)) 462 ((coeffpp) (c freevar))) 463 (n integerp1))))) 464 465;; Match d*p(x)*r1(x)^e1*r2(x)^e2*r0(x)^e0, where p(x) is a 466;; polynomial, e1 and e2 are both half an odd integer, and e3 is an 467;; integer. 468(defun factpow (exp var) 469 (m2 exp '((mtimes) ((coefftt) (d freevar)) 470 ((coefftt) (p polyp)) 471 ((mexpt) (r1 hasvar) 472 (e1 integerpfr)) 473 ((mexpt) (r2 hasvar) 474 (e2 integerpfr)) 475 ((mexpt) (r0 hasvar) 476 (e0 maxima-integerp))))) 477 478;; Match EXP to the form 479;; d*p/(e*x^2+f*x+g)*r0(x)^e0. p is a polynomial in x. 480(defun elliptquad (exp var) 481 (m2 exp '((mtimes) 482 ((coefftt) (d freevar)) 483 ((coefftt) (p polyp)) 484 ((mexpt) ((mplus) ((coeffpt) (e freevnz) ((mexpt) (x varp) 2)) 485 ((coeffpt) (f freevar) (x varp)) 486 ((coeffpp) (g freevar))) 487 -1) 488 ((mexpt) (r0 hasvar) 489 (e0 integerpfr))))) 490 491;; From the set of coefficients, generate the polynomial c*x^2+b*x+a. 492(defun polfoo (c b a x) 493 (add (mul c x x) 494 (mul b x) 495 a)) 496 497;; I think this is computing the list of coefficients of fun(x)/x^m, 498;; where fun(x) is a polynomial and m is a non-negative integer. The 499;; result is a list of two lists. The first list contains the 500;; polynomial part of fun(x)/x^m. The second is the principal part 501;; containing negative powers. 502;; 503;; Each of the lists is itself a list of power and coefficient itself. 504;; 505;; Thus (x+3)^2/x^2 = 1 + 6/x + 9/x^2 returns 506;; 507;; '(((0 1)) ((1 6) (2 9))) 508(defun powercoeflist (fun m var) 509 (prog ((expanfun (unquote ($expand (mul (prevconstexpan fun var) (power var m))))) 510 maxpowfun powfun coef poszpowlist negpowlist) 511 (when (and (equal fun 1) (plusp m)) 512 (return (cons nil (list (list (cons m (list 1))))))) 513 (when (and (equal fun 1) (minusp m)) 514 (return (cons nil (list (list (cons (- m) (list 1))))))) 515 (when (equal expanfun 1) 516 (return (cons (list (cons 0 (list 1))) (list nil)))) 517 (setq maxpowfun ($hipow expanfun var) 518 powfun ($lopow expanfun var)) 519 loop (setq coef ($coeff expanfun (power var powfun))) 520 (when (numberp coef) (go testjump)) 521 (go nojump) 522 testjump (when (and (not (zerop powfun)) (zerop coef)) 523 (go jump)) 524 nojump (when (plusp powfun) 525 (setq poszpowlist (append poszpowlist 526 (list (cons powfun (list coef)))))) 527 (when (zerop powfun) 528 (setq poszpowlist 529 (append poszpowlist 530 (list (cons 0 (list (consterm (cdr expanfun) var))))))) 531 (when (minusp powfun) 532 (setq negpowlist (append negpowlist (list (cons (- powfun) (list coef)))))) 533 (when (= powfun maxpowfun) 534 (return (list poszpowlist (reverse negpowlist)))) 535 jump (incf powfun) 536 (go loop))) 537 538(defun consterm (fun var) 539 (cond ((null fun) 0) 540 ((freeof var (car fun)) 541 (add (car fun) (consterm (cdr fun) var))) 542 (t (consterm (cdr fun) var)))) 543 544(defun prevconstexpan (fun var) 545 (cond ((atom fun) fun) 546 ((eq (caar fun) 'mplus) 547 (cond ((and (freeof var fun) 548 (not (inside fun 'mexpt))) 549 (list '(mquote) fun)) 550 ((and (freeof var fun) (inside fun 'mexpt)) 551 (list '(mquote) 552 (distrinplusprev (cdr fun) var))) 553 ((inside fun 'mexpt) 554 (distrinplusprev (cdr fun) var)) 555 (t fun))) 556 ((eq (caar fun) 'mtimes) 557 (distrintimesprev (cdr fun) var)) 558 ((and (not (inside (cdr fun) var)) 559 (eq (caar fun) 'mexpt)) 560 (power (prevconstexpan (cadr fun) var) (caddr fun))) 561 (t fun))) 562 563(defun distrinplusprev (fun var) 564 (cond ((null fun) 0) 565 (t (add (prevconstexpan (car fun) var) 566 (distrinplusprev (cdr fun) var))))) 567 568(defun distrintimesprev (fun var) 569 (cond ((null fun) 1) 570 (t (mul (prevconstexpan (car fun) var) 571 (distrintimesprev (cdr fun) var))))) 572 573(defun inside (fun arg) 574 (cond ((atom fun)(equal fun arg)) 575 ((inside (car fun) arg) t) 576 (t (inside (cdr fun) arg)))) 577 578(defun unquote (fun) 579 (cond ((not (inside fun 'mquote)) fun) 580 (t (unquote (meval fun))))) 581 582(defun checksigntm (expr) 583 (prog (aslist quest zerosigntest productcase) 584 (setq aslist checkcoefsignlist) 585 (cond ((atom expr) (go loop))) 586 (cond ((eq (caar expr) 'mtimes)(setq productcase t))) 587 loop (cond ((null aslist) 588 (setq checkcoefsignlist 589 (append checkcoefsignlist 590 (list (cons expr 591 (list 592 (setq quest (checkflagandact expr))))))) 593 (return quest))) 594 (cond ((equal (caar aslist) expr) (return (cadar aslist)))) 595 (setq aslist (cdr aslist)) 596 (go loop))) 597 598(defun checkflagandact (expr) 599 (cond (productcase 600 (setq productcase nil) 601 (findsignoftheirproduct (findsignofactors (cdr expr)))) 602 (t (asksign ($realpart expr))))) 603 604(defun findsignofactors (listofactors) 605 (cond ((null listofactors) nil) 606 ((eq zerosigntest '$zero) '$zero) 607 (t (append (list (setq zerosigntest (checksigntm (car listofactors)))) 608 (findsignofactors (cdr listofactors)))))) 609 610(defun findsignoftheirproduct (llist) 611 (prog (sign) 612 (cond ((eq llist '$zero) (return '$zero))) 613 (setq sign '$positive) 614 loop (cond ((null llist) (return sign))) 615 (cond ((eq (car llist) '$positive) 616 (setq llist (cdr llist)) 617 (go loop))) 618 (cond ((eq (car llist) '$negative) 619 (setq sign (changesign sign) llist (cdr llist)) 620 (go loop))) 621 (return '$zero))) 622 623(defun changesign (sign) 624 (if (eq sign '$positive) 625 '$negative 626 '$positive)) 627 628;; Integrate 1/sqrt(c*x^2+b*x+a). 629;; 630;; G&R 2.261 gives the following, where R = c*x^2+b*x+a and D = 631;; 4*a*b-b^2: 632;; 633;; c > 0: 634;; 1/sqrt(c)*log(2*sqrt(c*R)+2*c*x+b) 635;; 636;; c > 0, D > 0: 637;; 1/sqrt(c)*asinh((2*c*x+b)/sqrt(D)) 638;; 639;; c < 0, D < 0: 640;; -1/sqrt(-c)*asin((2*c*x+b)/sqrt(-D)) 641;; 642;; c > 0, D = 0: 643;; 1/sqrt(c)*log(2*c*x+b) 644;; 645(defun den1 (c b a x) 646 (let* ((expr (add (mul 2 c x) b)) ;; expr = 2*c*x+b 647 (signc (checksigntm (power c -1))) 648 (signb (checksigntm (power b 2))) 649 (signdiscrim (signdis2 c b a signc signb))) 650 (when (and (eq signc '$positive) 651 (eq signdiscrim '$negative)) 652 ;; c > 0, D > 0 653 (return-from den1 (augmult (mul* (power c -1//2) 654 `((%asinh) 655 ,(mul expr 656 (power (add (mul 4 c a) 657 (mul -1 b b)) 658 -1//2))))))) 659 (when (and (eq signc '$positive) 660 (eq signdiscrim '$zero)) 661 ;; c > 0, D = 0 662 (return-from den1 (augmult (mul* (power -1 expr) 663 (power c -1//2) 664 `((%log) ,expr))))) 665 (when (eq signc '$positive) 666 ;; c > 0 667 (return-from den1 (augmult (mul* (power c -1//2) 668 `((%log) 669 ,(add (mul 2 670 (power c 1//2) 671 (power (polfoo c b a x) 1//2)) 672 expr)))))) 673 (when (and (eq signc '$negative) 674 (eq signdiscrim '$positive)) 675 ;; c < 0, D > 0 676 (return-from den1 (augmult (mul* -1 677 (power (mul -1 c) -1//2) 678 `((%asin) 679 ,(mul expr 680 (power (add (mul b b) 681 (mul -4 c a)) 682 -1//2))))))) 683 (when (eq signc '$negative) 684 ;; c < 0. We try again, but flip the sign of the 685 ;; polynomial, and multiply by -%i. 686 (return-from den1 (augmult (mul (power -1 -1//2) 687 (den1 (mul -1 c) 688 (mul -1 b) 689 (mul -1 a) 690 x))))))) 691 692;; Compute sign of discriminant of the quadratic c*x^2+b*x+a. That 693;; is, the sign of b^2-4*c*a. 694(defun signdiscr (c b a) 695 (checksigntm (simplifya (add (power b 2) (mul -4 c a)) nil))) 696 697(defun askinver (a) 698 (checksigntm (inv a))) 699 700(defun signdis1 (c b a) 701 (cond ((equal (mul b a) 0) 702 (if (and (equal b 0) (equal a 0)) 703 '$zero 704 '$nonzero)) 705 (t 706 ;; Why are we checking the sign of (b^2-4*a*c)^2? 707 (checksigntm (power (add (mul b b) (mul -4 c a)) 2))))) 708 709;; Check sign of discriminant of c*x^2+b*x+a, but also taking into 710;; account the sign of c and b. 711(defun signdis2 (c b a signc signb) 712 (cond ((equal signb '$zero) 713 (cond ((equal a 0) '$zero) 714 (t (let ((askinv (askinver a))) 715 (if (or (and (eq signc '$positive) 716 (eq askinv '$negative)) 717 (and (eq signc '$negative) 718 (eq askinv '$positive))) 719 '$positive 720 '$negative))))) 721 (t (if (equal a 0) 722 '$positive 723 (signdiscr c b a))))) 724 725(defun signdis3 (c b a signa) 726 (declare (special *ec-1*)) 727 (cond ((equal b 0) 728 (if (equal (checksigntm *ec-1*) signa) 729 '$negative 730 '$positive)) 731 (t (signdiscr c b a)))) 732 733;; Integrate things like x^m*R^(p-1/2), p > 0, m > 0. 734;; 735;; I think pluspowfo1 = p - 1. 736(defun nummnumn (poszpowlist pluspowfo1 p c b a x) 737 (declare (special *ec-1*)) 738 (let ((expr (power (polfoo c b a x) (add p 1//2))) ;; expr = R^(p+1/2) 739 (expo *ec-1*) ;; expo = 1/c 740 (ex (power c -2))) ;; ex = 1/c^2 741 (prog ((result 0) 742 (controlpow (caar poszpowlist)) 743 (coef (cadar poszpowlist)) 744 count res1 res2 m partres) 745 #+nil (format t "p = ~A~%pluspowfo1 = ~A~%" p pluspowfo1) 746 (when (zerop controlpow) 747 ;; Integrate R^(p-1/2) 748 (setq result (augmult (mul coef (numn pluspowfo1 c b a x))) 749 count 1) 750 (go loop)) 751 752 jump1 753 ;; Handle x*R^(p-1/2) 754 ;; 755 ;; G&R 2.260, m = 1 756 ;; 757 ;; integrate(x*R^(2*p-1),x) = 758 ;; R^(p+1/2)/(2*p+1)/c 759 ;; - b/2/c*integrate(sqrt(R^(2*p-1)),x) 760 (setq res1 (add (augmult (mul expr expo 761 (power (+ p p 1) -1))) 762 (augmult (mul -1 b 1//2 expo 763 (numn pluspowfo1 c b a x))))) 764 (when (equal controlpow 1) 765 (setq result (add result (augmult (mul coef res1))) 766 count 2) 767 (go loop)) 768 769 jump2 770 ;; Handle x^2*R^(p-1/2) 771 ;; 772 ;; G&R 2.260, m = 2 773 ;; 774 ;; integrate(x^2*R^(2*p-1),x) = 775 ;; x*R^(p+1/2)/(2*p+2)/c 776 ;; - (2*p+3)*b/2/(2*p+2)/c*integrate(x*sqrt(R^(2*p-1)),x) 777 ;; - a/(2*p+2)/c*integrate(sqrt(R^(2*p-1)),x) 778 (setq res2 (add (augmult (mul* x expr expo 779 (inv (+ p p 2)))) 780 (augmult (mul* b (+ p p 3) 781 '((rat) -1 4) 782 ex 783 (inv (+ p p p 1 784 (* p p) 785 (* p p))) 786 expr)) 787 (augmult (mul (inv (1+ p)) 788 ex 789 '((rat simp) 1 8) 790 (add (mul (power b 2) 791 (+ p p 3)) 792 (mul -4 a c)) 793 (numn pluspowfo1 c b a x))))) 794 (when (equal controlpow 2) 795 (setq result (add result (augmult (mul coef res2))) 796 count 3) 797 (go loop)) 798 799 jump3 800 (setq count 4 801 m 3) 802 jump 803 ;; The general case: x^m*R^(p-1/2) 804 (setq partres 805 (let ((pro (inv (+ m p p)))) 806 ;; pro = 1/(m+2*p) 807 ;; 808 ;; G&R 2.260, m = 2 809 ;; 810 ;; integrate(x^m*R^(2*p-1),x) = 811 ;; x^(m-1)*R^(p+1/2)/(m+2*p)/c 812 ;; - (2*m+2*p-1)*b/2/(m+2*p)/c*integrate(x^(m-1)*sqrt(R^(2*p-1)),x) 813 ;; - (m-1)*a/(m+2*p)/c*integrate(x^(m-2)*sqrt(R^(2*p-1)),x) 814 (add (augmult (mul (power x (1- m)) 815 expr expo pro)) 816 (augmult (mul -1 b (+ p p m m -1) 817 1//2 expo pro res2)) 818 (augmult (mul -1 a (1- m) 819 expo pro res1))))) 820 (incf m) 821 (when (> m controlpow) 822 (setq result (add result (augmult (mul coef partres)))) 823 (go loop)) 824 825 jump4 826 (setq res1 res2 827 res2 partres) 828 (go jump) 829 830 loop 831 (setq poszpowlist (cdr poszpowlist)) 832 (when (null poszpowlist) (return result)) 833 (setq coef (cadar poszpowlist)) 834 (setq controlpow (caar poszpowlist)) 835 (when (equal count 4) (go jump4)) 836 (when (equal count 1) (go jump1)) 837 (when (equal count 2) (go jump2)) 838 (go jump3)))) 839 840;; Integrate R^(p+1/2) 841(defun numn (p c b a x) 842 (declare (special *ec-1*)) 843 (let ((exp1 *ec-1*) ;; exp1 = 1/c 844 (exp2 (add b (mul 2 c x))) ;; exp2 = b+2*c*x 845 (exp4 (add (mul 4 a c) (mul -1 b b))) ;; exp4 = 4*a*c-b^2 846 (exp5 (div 1 (1+ p)))) ;; exp5 = 1/(p+1) 847 (if (zerop p) 848 ;; integrate(sqrt(R),x) 849 ;; 850 ;; G&R 2.262 says 851 ;; 852 ;; integrate(sqrt(R),x) = 853 ;; (2*c*x+b)*sqrt(R)/4/c + del/8/c*integrate(1/sqrt(R),x) 854 ;; 855 ;; del = 4*a*c-b^2 856 (add (augmult (mul '((rat simp) 1 4) 857 exp1 exp2 858 (power (polfoo c b a x) 1//2))) 859 (augmult (mul '((rat simp) 1 8) 860 exp1 exp4 861 (den1 c b a x)))) 862 863 ;; The general case 864 ;; 865 ;; G&R 2.260, eq. 2: 866 ;; 867 ;; integrate(sqrt(R^(2*p+1)),x) = 868 ;; (2*c*x+b)/4/(p+1)/c*R^(p+1/2) 869 ;; + (2*p+1)*del/8/(p+1)/c*integrate(sqrt(R^(2*p-1)),x) 870 (add (augmult (mul '((rat simp) 1 4) 871 exp1 exp5 exp2 872 (power (polfoo c b a x) (add p 1//2)))) 873 (augmult (mul '((rat simp) 1 8) 874 exp1 exp5 (+ p p 1) 875 exp4 876 (numn (1- p) c b a x))))))) 877 878(defun augmult (x) 879 ($multthru (simplifya x nil))) 880 881;; Integrate things like 1/x^m/R^(p+1/2), p > 0. 882(defun denmdenn (negpowlist p c b a x) 883 (let ((exp1 (power (polfoo c b a x) (add 1//2 (- p))))) ;; exp1 = 1/R^(p-1/2) 884 (prog (result controlpow coef count res1 res2 m partres 885 (signa (checksigntm (simplifya a nil))) 886 ea-1) 887 (when (eq signa '$zero) 888 (return (noconstquad negpowlist p c b x))) 889 (setq result 0 890 controlpow (caar negpowlist) 891 ea-1 (power a -1)) 892 (setq coef (cadar negpowlist)) 893 (when (zerop controlpow) 894 ;; I'm not sure we ever get here because m = 0 is 895 ;; usually handled elsewhere. 896 (setq result (augmult (mul coef (denn p c b a x))) 897 count 1) 898 (go loop)) 899 900 jump1 901 ;; Handle 1/x/R^(p+1/2) 902 (setq res1 (den1denn p c b a x)) 903 (when (equal controlpow 1) 904 (setq result (add result (augmult (mul coef res1))) 905 count 2) 906 (go loop)) 907 908 jump2 909 ;; Handle 1/x^2/R^(p+1/2) 910 ;; 911 ;; G&R 2.268, m = 2 912 ;; 913 ;; integrate(1/x^2/R^(p+1/2),x) = 914 ;; -1/a/x/sqrt(R^(2*p-1)) 915 ;; -(2*p+1)*b/2/a*integrate(1/x/sqrt(R^(2*p+1)),x) 916 ;; -2*p*c/a*integrate(1/sqrt(R^(2*p+1)),x) 917 (setq res2 (add (augmult (mul -1 ea-1 (power x -1) exp1)) 918 (augmult (mul -1 b (+ 1 p p) 1//2 919 ea-1 (den1denn p c b a x))) 920 (augmult (mul -2 p c ea-1 (denn p c b a x))))) 921 (when (equal controlpow 2) 922 (setq result (add result (augmult (mul coef res2))) 923 count 3) 924 (go loop)) 925 926 jump3 927 (setq count 4 928 m 3) 929 930 jump 931 ;; General case 1/x^m/R^(p+1/2) 932 ;; 933 ;; G&R 2.268 934 ;; 935 ;; integrate(1/x^2/R^(p+1/2),x) = 936 ;; -1/(m-1)/a/x^(m-1)/sqrt(R^(2*p-1)) 937 ;; -(2*p+2*m-3)*b/2/(m-1)/a*integrate(1/x^(m-1)/sqrt(R^(2*p+1)),x) 938 ;; -(2*n+m-2)*c/(m-1)/a*integrate(1/x^(m-2)/sqrt(R^(2*p+1)),x) 939 (setq partres 940 (let ((exp2 (div -1 (1- m)))) 941 ;; exp2 = -1/(m-1) 942 (add (augmult (mul exp2 ea-1 943 (power x (1+ (- m))) 944 exp1)) 945 (augmult (mul b (+ p p m m -3) 1//2 946 ea-1 exp2 res2)) 947 (augmult (mul c ea-1 exp2 948 (+ p p m -2) res1))))) 949 (incf m) 950 (when (> m controlpow) 951 (setq result (add result (augmult (mul coef partres)))) 952 (go loop)) 953 954 jump4 955 (setq res1 res2 res2 partres) 956 (go jump) 957 958 loop 959 (setq negpowlist (cdr negpowlist)) 960 (when (null negpowlist) (return result)) 961 (setq coef (cadar negpowlist) 962 controlpow (caar negpowlist)) 963 (when (equal count 4) (go jump4)) 964 (when (equal count 1) (go jump1)) 965 (when (equal count 2) (go jump2)) 966 (go jump3)))) 967 968;; Integral of 1/(c*x^2+b*x+a)^(n), n > 0. p = n + 1/2. 969;; 970;; See G&R 2.263 formula 3. 971;; 972;; Let R = c*x^2+b*x+a. 973(defun denn (p c b a x) 974 (let ((signdisc (signdis1 c b a)) 975 (exp1 (add b (mul 2 c x))) ;; exp1 = b + 2*c*x; 976 (exp2 (add (mul 4 a c) (mul b b -1))) ;; exp2 = (4*a*c-b^2) 977 (exp3 (inv (+ p p -1))) ;; exp3 = 1/(2*p-1) 978 (*ec-1* (inv c))) 979 (declare (special *ec-1*)) 980 #+nil (format t "signdisc = ~A~%p = ~A~%" signdisc p) 981 (cond ((and (eq signdisc '$zero) (zerop p)) 982 ;; 1/sqrt(R), and R has duplicate roots. That is, we have 983 ;; 1/sqrt(c*(x+b/(2c))^2) = 1/sqrt(c)/sqrt((x+b/2/c)^2). 984 ;; 985 ;; We return 1/sqrt(c)*log(x+b/2/c). Shouldn't we return 986 ;; 1/c*log(|x+b/2/c|)? 987 (augmult (mul* (power *ec-1* 1//2) 988 `((%log) ,(add x (mul b 1//2 *ec-1*)))))) 989 ((and (eq signdisc '$zero) (plusp p)) 990 ;; 1/sqrt(R^(2*p+1)), with duplicate roots. 991 ;; 992 ;; That is, 1/sqrt((c*(x+b/2/c)^(2)^(2*p+1))), or 993 ;; 1/c^(p+1/2)/(x+b/2/c)^(2*p+1). So the result is 994 ;; -1/2/p*c^(-p-1/2)/(x+b/2/c)^(2*p) 995 (augmult (mul (div -1 (+ p p)) 996 (power c (mul -1//2 (+ p p 1))) 997 (power (add x (mul b 1//2 *ec-1*)) (* -2 p))))) 998 ((zerop p) 999 ;; 1/sqrt(R) 1000 (den1 c b a x)) 1001 ((equal p 1) 1002 ;; 1/sqrt(R^3). 1003 ;; 1004 ;; G&R 2.264 eq. 5 says 1005 ;; 1006 ;; 2*(2*c*x+b)/del/sqrt(R). 1007 (augmult (mul 2 exp1 (inv exp2) 1008 (power (polfoo c b a x) -1//2)))) 1009 (t 1010 ;; The general case. G&R 2.263 eq. 3. 1011 ;; 1012 ;; integrate(1/sqrt(R^(2*p+1)),x) = 1013 ;; 2*(2*c*x+b)/(2*p-1)/c/sqrt(R^(2*p-1)) 1014 ;; + 8*(p-1)*c/(2*p-1)/del*integrate(1/sqrt(R^(2*p-1)),x) 1015 ;; 1016 ;; where del = 4*a*c-b^2. 1017 (add (augmult (mul 2 exp1 1018 exp3 (inv exp2) 1019 (power (polfoo c b a x) 1020 (add 1//2 (- p))))) 1021 (augmult (mul 8 c (1- p) 1022 exp3 (inv exp2) 1023 (denn (1- p) c b a x)))))))) 1024 1025;; Integral of 1/x/(c*x^2+b*x+a)^(p+1/2), p > 0. 1026(defun den1denn (p c b a x) 1027 (let ((signa (checksigntm (power a 2))) ;; signa = sign of a^2 1028 (ea-1 (inv a))) ;; ea-1 = 1/a 1029 (cond ((eq signa '$zero) 1030 ;; This is wrong because noconstquad expects a 1031 ;; powercoeflist for th first arg. But I don't think 1032 ;; there's any way to get here from anywhere. I'm pretty 1033 ;; sure den1denn is never called with a equal to 0. 1034 (noconstquad 1 p c b x)) 1035 ((zerop p) 1036 ;; 1/x/sqrt(c*x^2+b*x+a) 1037 (den1den1 c b a x)) 1038 (t 1039 ;; The general case. See G&R 2.268: 1040 ;; 1041 ;; R=(c*x^2+b*x+a) 1042 ;; 1043 ;; integrate(1/x/sqrt(R^(2*p+1)),x) = 1044 ;; 1045 ;; 1/(2*p-1)/a/sqrt(R^(2*p-1)) 1046 ;; - b/2/a*integrate(1/sqrt(R^(2*p+1)),x) 1047 ;; + 1/a*integrate(1/x/sqrt(R^(2*p-1)),x) 1048 (add (augmult (mul (inv (+ p p -1)) 1049 ea-1 1050 (power (polfoo c b a x) 1051 (add 1//2 (- p))))) 1052 (augmult (mul ea-1 (den1denn (1- p) c b a x))) 1053 (augmult (mul -1 1//2 ea-1 b (denn p c b a x)))))))) 1054 1055;; Integral of 1/x/sqrt(c*x^2+b*x+a). 1056;; 1057;; G&R 2.266 gives the following results, where del is the 1058;; discriminant 4*a*c-b^2. (I think this is the opposite of what we 1059;; compute below for the discriminant.) 1060;; 1061(defun den1den1 (c b a x) 1062 (let ((exp2 (add (mul b x) a a)) ; exp2 = b*x+2*a 1063 (exp3 (inv (simplify (list '(mabs) x))))) ; exp3 = 1/abs(x) 1064 (prog (signdiscrim 1065 (condition (add (mul b x) a a)) 1066 (signa (checksigntm (simplifya a nil))) 1067 exp1) 1068 (when (eq signa '$zero) 1069 (return (noconstquad '((1 1)) 0 c b x))) 1070 (setq signdiscrim (signdis3 c b a signa) 1071 exp1 (power a (inv -2))) 1072 #+nil (format t "signa = ~A~%signdiscrim = ~A~%" signa signdiscrim) 1073 1074 (when (and (eq signa '$positive) 1075 (eq signdiscrim '$negative)) 1076 ;; G&R case a > 0, del > 0 1077 ;; 1078 ;; -1/sqrt(a)*asinh((2*a+b*x)/x/sqrt(del) 1079 (return (mul* -1 exp1 1080 `((%asinh) 1081 ,(augmult (mul exp2 exp3 1082 (power (add (mul 4 a c) 1083 (mul -1 b b)) 1084 -1//2))))))) 1085 (when (and (eq signdiscrim '$zero) 1086 (eq signa '$positive)) 1087 ;; G&R case del = 0, a > 0 1088 ;; 1089 ;; 1/sqrt(a)*log(x/(2*a+b*x)) 1090 (return (mul* (power -1 condition) 1091 -1 exp1 1092 `((%log) ,(augmult (mul exp3 exp2)))))) 1093 (when (eq signa '$positive) 1094 ;; G&R case a > 0 1095 ;; 1096 ;; -1/sqrt(a)*log((2*a+b*x+2*sqrt(a*R))/x) 1097 ;; 1098 ;; R = c*x^2+b*x+a. 1099 (return (mul* -1 exp1 1100 `((%log) 1101 ,(add b (mul 2 a exp3) 1102 (mul 2 exp3 1103 (power a 1//2) 1104 (power (polfoo c b a x) 1//2))))))) 1105 (when (and (eq signa '$negative) 1106 (eq signdiscrim '$positive)) 1107 ;; G&R case a < 0, del < 0 1108 ;; 1109 ;; 1/sqrt(-a)*asin((2*a+b*x)/x/sqrt(b^2-4*a*c)) 1110 (return (mul* (power (mul -1 a) -1//2) 1111 `((%asin) 1112 ,(augmult (mul exp2 exp3 1113 (power (add (mul b b) (mul -4 a c)) -1//2))))))) 1114 ;; I think this is the case of a < 0. We flip the sign of 1115 ;; coefficients of the quadratic to make a positive, and 1116 ;; multiply the result by 1/%i. 1117 ;; 1118 ;; Why can't we use the case a < 0 in G&R 2.266: 1119 ;; 1120 ;; 1/sqrt(-a)*atan((2*a+b*x)/2/sqrt(-a)/sqrt(R) 1121 ;; 1122 ;; FIXME: Why the multiplication by -1? 1123 (return (mul #+nil -1 1124 (power -1 1//2) 1125 (den1den1 (mul -1 c) (mul -1 b) (mul -1 a) x)))))) 1126 1127;; Integral of d/x^m/(c*x^2+b*x)^(p+1/2), p > 0. The values of m and 1128;; d are in NEGPOWLIST. 1129(defun noconstquad (negpowlist p c b x) 1130 (let ((exp1 (inv (+ p p 1))) ;; exp1 = 1/(2*p+1) 1131 (exp2 (inv x)) ;; exp2 = 1/x 1132 (exp3 (- p))) ;; exp3 = -p 1133 (prog (result controlpow coef count res1 signb m partres eb-1) 1134 (setq signb (checksigntm (power b 2))) 1135 (when (eq signb '$zero) 1136 (return (trivial1 negpowlist p c x))) 1137 (setq result 0 1138 controlpow (caar negpowlist) 1139 coef (cadar negpowlist) 1140 eb-1 (inv b)) 1141 (when (zerop controlpow) 1142 ;; Not sure if we ever actually get here. The case of 1143 ;; m=0 is usually handled at a higher level. 1144 (setq result (augmult (mul coef (denn p c b 0 x))) 1145 count 1) 1146 (go loop)) 1147 jump1 1148 ;; Handle 1/x/R^(p+1/2) 1149 ;; 1150 ;; G&R 2.268, a = 0, m = 1 1151 ;; 1152 ;; integrate(1/x^m/sqrt(R^(2*p+1)),x) = 1153 ;; -2/(2*p+1)/b/x/sqrt(R^(2*p-1)) 1154 ;; -4*p*c/(2*p+1)/b*integrate(1/sqrt(R^(2*p+1)),x) 1155 (setq res1 (add (augmult (mul -2 exp1 eb-1 exp2 1156 (power (polfoo c b 0 x) 1157 (add 1//2 exp3)))) 1158 (augmult (mul -4 p c exp1 eb-1 (denn p c b 0 x))))) 1159 (when (equal controlpow 1) 1160 (setq result (add result (augmult (mul coef res1))) 1161 count 2) 1162 (go loop)) 1163 jump2 (setq count 3 m 2) 1164 jump 1165 ;; Handle general case 1/x^m/R^(p+1/2) 1166 ;; 1167 ;; G&R 2.268, a = 0 1168 ;; 1169 ;; integrate(1/x^m/sqrt(R^(2*p+1)),x) = 1170 ;; -2/(2*p+2*m-1)/b/x^m/sqrt(R^(2*p+1)) 1171 ;; -(4*p+2*m-2)*c/(2*p+2*m-1)/b*integrate(1/x^(m-1)/sqrt(R^(2*p+1)),x) 1172 (setq partres 1173 (add (augmult (mul -2 (inv (+ p p m m -1)) 1174 eb-1 1175 (power x (mul -1 m)) 1176 (power (polfoo c b 0 x) 1177 (add 1//2 exp3)))) 1178 (augmult (mul -2 c (+ p p m -1) 1179 eb-1 1180 (inv (+ p p m m -1)) 1181 res1)))) 1182 (incf m) 1183 (when (> m controlpow) 1184 (setq result (add result (augmult (mul coef partres)))) 1185 (go loop)) 1186 jump3 1187 (setq res1 partres) 1188 (go jump) 1189 loop 1190 (setq negpowlist (cdr negpowlist)) 1191 (when (null negpowlist) (return result)) 1192 (setq coef (cadar negpowlist) 1193 controlpow (caar negpowlist)) 1194 (when (= count 3) (go jump3)) 1195 (when (= count 1) (go jump1)) 1196 (go jump2)))) 1197 1198;; The trivial case of d/x^m/(c*x^2+b*x+a)^(p+1/2), p > 0, and a=b=0. 1199(defun trivial1 (negpowlist p c x) 1200 (cond ((null negpowlist) 0) 1201 (t 1202 ;; d/x^m/c^(p+1/2)/x^(2*p+1) = d/c^(p+1/2)/x^(m+2*p+1) 1203 ;; The integral is obviously 1204 ;; 1205 ;; -d/c^(p+1/2)/x^(m+2*p)/(m+2*p) 1206 (add (augmult (mul (power x 1207 (add (* -2 p) 1208 (mul -1 (caar negpowlist)))) 1209 (cadar negpowlist) 1210 (power c (add (- p) -1//2)) 1211 (inv (add (* -2 p) 1212 (mul -1 (caar negpowlist)))))) 1213 (trivial1 (cdr negpowlist) p c x))))) 1214 1215;; Integrate pl(x)/(c*x^2+b*x+a)^(p+1/2) where pl(x) is a polynomial 1216;; and p > 0. The polynomial is given in POSZPOWLIST. 1217(defun nummdenn (poszpowlist p c b a x) 1218 (declare (special *ec-1*)) 1219 (let ((exp1 (inv (+ p p -1))) ;; exp1 = 1/(2*p-1) 1220 (exp2 (power (polfoo c b a x) (add 1//2 (- p)))) ;; exp2 = (a*x^2+b*x+c)^(p-1/2) 1221 (exp3 (add (mul 4 a c) (mul -1 b b))) ;; exp3 = (4*a*c-b^2) (negative of the discriminant) 1222 (exp4 (add x (mul b 1//2 *ec-1*))) ;; exp4 = x+b/2/c 1223 (exp5 (power c -2)) ;; exp5 = 1/c^2 1224 (exp6 (+ 2 (* -2 p))) ;; exp6 = -2*p+2 1225 (exp7 (1+ (* -2 p)))) ;; exp7 = -2*p+1 1226 (prog (result controlpow coef count res1 res2 m partres signdiscrim) 1227 ;; Let S=R^(p+1/2). 1228 ;; 1229 ;; We are trying to integrate pl(x)/S 1230 ;; = (p0 + p1*x + p2*x^3 + ...)/S 1231 ;; 1232 ;; So the general term is pm*x^m/S. This integral is given by 1233 ;; G&R 2.263, eq.1: 1234 ;; 1235 ;; integrate(x^m/sqrt(R^(2*p+1)),x) = 1236 ;; 1237 ;; x^(m-1)/(m-2*n)/sqrt(R^(2*p-1)) 1238 ;; - (2*m-2*n-1)*b/2/(m-2*n)/c*integrate(x^(m-1)/sqrt(R^(2*p+1)),x) 1239 ;; - (m-1)*a/(m-2*n)/c*integrate(x^(m-2)/sqrt(R^(2*p+1)),x) 1240 ;; 1241 ;; Thus the integration of x^m/S involves x^(m-1)/S and x^(m-2)/S. 1242 ;; 1243 ;; I think what this loop does is integrate x^(m-1)/S and 1244 ;; x^(m-2)/S, and remember them so that we can integrate x^m/S 1245 ;; without having to do all the integrals again. Thus we 1246 ;; start with the constant term and work our way up to the 1247 ;; highest term. 1248 ;; 1249 ;; I think this would be much simpler if we used memoization 1250 ;; and start with the highest power. Then all the 1251 ;; intermediate forms will have been computed, and we can just 1252 ;; simply integrate all the lower terms by looking them up. 1253 (setq result 0 1254 controlpow (caar poszpowlist)) 1255 (setq coef (cadar poszpowlist) 1256 signdiscrim (signdis1 c b a)) 1257 ;; We're looking at coef*x^controlpow/R^(p+1/2) right now. 1258 (when (zerop controlpow) 1259 ;; Actually it's coef/R^(p+1/2), so handle that now, go 1260 ;; to the next coefficient. 1261 (setq result (augmult (mul coef (denn p c b a x))) 1262 count 1) 1263 (go loop)) 1264 1265 jump1 1266 ;; 1267 ;; This handles the case coef*x/R^(p+1/2) 1268 ;; 1269 ;; res1 = -1/c/(2*p-1)*R^(p-1/2) 1270 ;; -b/2/c*integrate(R^(p+1/2),x) 1271 ;; 1272 (setq res1 1273 (add (augmult (mul -1 *ec-1* exp1 exp2)) 1274 (augmult (mul b -1//2 *ec-1* (denn p c b a x))))) 1275 (when (= controlpow 1) 1276 ;; Integrate coef*x/R^(p+1/2). 1277 ;; 1278 ;; x/R^(p+1/2) is in res1. 1279 (setq result (add result (augmult (mul coef res1))) 1280 count 2) 1281 (go loop)) 1282 jump2 1283 ;; This handles the case coef*x^2/R^(p+1/2) 1284 (when (and (plusp p) (not (eq signdiscrim '$zero))) 1285 ;; p > 0, no repeated roots 1286 (setq res2 1287 (add (augmult (mul *ec-1* exp1 (inv exp3) exp2 1288 (add (mul 2 a b) 1289 (mul 2 b b x) 1290 (mul -4 a c x)))) 1291 (augmult (mul *ec-1* (inv exp3) exp1 1292 (add (mul 4 a c) 1293 (mul 2 b b p) 1294 (mul -3 b b)) 1295 (denn (+ p -1) 1296 c b a x)))))) 1297 (when (and (zerop p) (not (eq signdiscrim '$zero))) 1298 ;; x^2/sqrt(R), no repeated roots. 1299 ;; 1300 ;; G&R 2.264, eq. 3 1301 ;; 1302 ;; integrate(x^2/sqrt(R),x) = 1303 ;; (x/2/c-3*b/4/c^2)*sqrt(R) 1304 ;; + (3*b^2/8/c^2-a/2/c)*integrate(1/sqrt(R),x) 1305 ;; 1306 ;; = (2*c*x-3*b)/4/c^2*sqrt(R) 1307 ;; + (3*b^2-4*a*c)/8/c^2*integrate(1/sqrt(R),x) 1308 (setq res2 1309 (add (augmult (mul '((rat simp) 1 4) 1310 exp5 1311 (add (mul 2 c x) 1312 (mul -3 b)) 1313 (power (polfoo c b a x) 1314 1//2))) 1315 (augmult (mul '((rat simp) 1 8) 1316 exp5 1317 (add (mul 3 b b) 1318 (mul -4 a c)) 1319 (den1 c b a x)))))) 1320 (when (and (zerop p) (eq signdiscrim '$zero)) 1321 ;; x^2/sqrt(R), repeated roots 1322 ;; 1323 ;; With repeated roots, R is really of the form 1324 ;; c*x^2+b*x+b^2/4/c = c*(x+b/2/c)^2. So we have 1325 ;; 1326 ;; x^2/sqrt(c)/(x+b/2/c) 1327 ;; 1328 ;; And the integral of this is 1329 ;; 1330 ;; b^2*log(x+b/2/c)/4/c^(5/2) + x^2/2/sqrt(c) - b*x/2/c^(3/2). 1331 ;; 1332 (setq res2 1333 ;; (add (augmult (mul* b b (list '(rat) 1 4) 1334 ;; (power c -3) 1335 ;; (list '(%log) exp4))) 1336 ;; (augmult (mul *ec-1* 1//2 (power exp4 2))) 1337 ;; (augmult (mul -1 b x exp5))) 1338 (add (augmult (mul* b b '((rat) 1 4) 1339 (power c (div -5 2)) 1340 `((%log) ,exp4))) 1341 (augmult (mul (power c -1//2) 1//2 (power x 2))) 1342 (augmult (mul -1//2 b x (power c (div -3 2))))))) 1343 1344 (when (and (= p 1) (eq signdiscrim '$zero)) 1345 ;; x^2/sqrt(R^3), repeated roots 1346 ;; 1347 ;; As above, we have c*(x+b/2/c)^2, so 1348 ;; 1349 ;; x^2/sqrt(R^3) = x^2/sqrt(c^3)/(x+b/2/c)^3 1350 ;; 1351 ;; And the integral is 1352 ;; 1353 ;; log(x+b/2/c)/c^(3/2) + z*(3*z+4*x)/2/c^(3/2)/(z+x)^2 1354 ;; 1355 ;; where z = b/2/c. 1356 (setq res2 1357 ;; (add (augmult (mul* *ec-1* (list '(%log) exp4))) 1358 ;; (augmult (mul b exp5 (power exp4 -1))) 1359 ;; (augmult (mul (list '(rat) -1 8) 1360 ;; (power c -3) 1361 ;; b b (power exp4 -2)))) 1362 (add (augmult (mul* (power c (div -3 2)) `((%log) ,exp4))) 1363 (augmult (mul b x (power c (div -5 2)) (power exp4 -2))) 1364 (augmult (mul (div 3 8) 1365 (power c (div -7 2)) 1366 b b (power exp4 -2)))))) 1367 1368 (when (and (eq signdiscrim '$zero) (> p 1)) 1369 ;; x^2/R^(p+1/2), repeated roots, p > 1 1370 ;; 1371 ;; As above, we have R=c*(x+b/2/c)^2, so the integral is 1372 ;; 1373 ;; x^2/(x+b/2/c)^(2*p+1)/c^(p+1/2). 1374 ;; 1375 ;; Let d = b/2/c. Then write x^2 = 1376 ;; (x+d)^2-2*d*(x+d)+d^2. The integrand becomes 1377 ;; 1378 ;; 1/(x+d)^(2*p-1) - 2*d/(x+d)^(2*p) + d^2/(x+d)^(2*p+1) 1379 ;; 1380 ;; whose integral is 1381 ;; 1382 ;; 1/(2*p-2)/(x+d)^(2*p-2) - 2*d/(2*p-1)/(x+d)^(2*p-1) 1383 ;; + d^2/(2*p)/(x+d)^(2*p) 1384 ;; 1385 ;; And don't forget the factor c^(-p-1/2). Finally, 1386 ;; 1387 ;; 1/c^(p+1/2)/(2*p-2)/(x+d)^(2*p-2) 1388 ;; - b/c^(p+3/2)/(2*p-1)/(x+d)^(2*p-1) 1389 ;; + b^2/8/c^(p+5/2)/p/(x+d)^(2*p) 1390 (setq res2 1391 ;; (add (augmult (mul *ec-1* 1392 ;; (power exp4 exp6) 1393 ;; (inv exp6))) 1394 ;; (augmult (mul -1 b exp5 (inv exp7) 1395 ;; (power exp4 exp7))) 1396 ;; (augmult (mul b b (list '(rat) -1 8) 1397 ;; (power c -3) 1398 ;; (inv p) 1399 ;; (power exp4 1400 ;; (* -2 p))))) 1401 (add (augmult (mul (inv (power c (add p 1//2))) 1402 (power exp4 exp6) 1403 (inv exp6))) 1404 (augmult (mul -1 b 1405 (inv (power c (add p (div 3 2)))) 1406 (inv exp7) 1407 (power exp4 exp7))) 1408 (augmult (mul b b '((rat simp) -1 8) 1409 (inv (power c (add p (div 5 2)))) 1410 (inv p) 1411 (power exp4 1412 (* -2 p))))))) 1413 (when (= controlpow 2) 1414 ;; x^2/R^(p+1/2) 1415 ;; 1416 ;; We computed this result above, so just multiply by 1417 ;; the coefficient and add it to the result. 1418 (setq result (add result (augmult (mul coef res2))) 1419 count 3) 1420 (go loop)) 1421 jump3 1422 (setq count 4 1423 m 3) 1424 jump 1425 ;; coef*x^m/R^(p+1/2). m >= 3 1426 (setq partres 1427 (let ((denom (+ p p (- m)))) 1428 ;; denom = 2*p-m 1429 1430 ;; G&R 2.263, eq 1: 1431 ;; 1432 ;; integrate(x^m/sqrt(R^(2*p+1)),x) = 1433 ;; x^(m-1)/c/(m-2*p)/sqrt(R^(2*p-1)) 1434 ;; - b*(2*m-2*p-1)/2/(m-2*p)*integrate(x^(m-1)/sqrt(R^(2*p+1)),x) 1435 ;; + (m-1)*a/(m-2*p)/c*integrate(x^(m-2)/sqrt(R^(2*p+1)),x) 1436 ;; 1437 ;; The two integrals here were computed above in res2 1438 ;; and res1, respectively. 1439 (add (augmult (mul* (power x (1- m)) 1440 *ec-1* (div -1 denom) 1441 (power (polfoo c b a x) 1442 (add 1//2 (- p))))) 1443 (augmult (mul b (+ p p 1 (* -2 m)) 1444 -1//2 1445 *ec-1* (inv denom) res2)) 1446 (augmult (mul a (1- m) *ec-1* (inv denom) res1))))) 1447 on 1448 ;; Move on to next higher power 1449 (incf m) 1450 (when (> m controlpow) 1451 (setq result (add result (augmult (mul coef partres)))) 1452 (go loop)) 1453 jump4 1454 (setq res1 res2 1455 res2 partres) 1456 (when (= m (+ p p)) 1457 (setq partres 1458 (let ((expr (nummdenn (list (list (1- m) 1)) p c b a x))) 1459 (add (mul x expr) 1460 (mul -1 (distrint (cdr ($expand expr)) x))))) 1461 (go on)) 1462 (go jump) 1463 loop 1464 ;; Done with first term of polynomial. Exit if we're done. 1465 (setq poszpowlist (cdr poszpowlist)) 1466 (when (null poszpowlist) (return result)) 1467 (setq coef (cadar poszpowlist) controlpow (caar poszpowlist)) 1468 (when (= count 4) (go jump4)) 1469 (when (= count 1) (go jump1)) 1470 (when (= count 2) (go jump2)) 1471 (go jump3)))) 1472 1473;; Integrate functions of the form coef*R^(pow-1/2)/x^m. NEGPOWLIST 1474;; contains the list of coef's and m's. 1475(defun denmnumn (negpowlist pow c b a x) 1476 (let ((exp1 (inv x)) ;; exp1 = 1/x 1477 (exp2 (+ pow pow -1))) ;; exp2 = 2*pow-1 1478 (prog (result controlpow coef count res1 res2 m partres signa ea-1 1479 (p (+ pow pow -1))) ;; p = 2*pow-1. 1480 ;; NOTE: p is not the same here as in other routines! 1481 ;; Why is there this special case for negpowlist? 1482 ;; CASE1 calls this in this way. 1483 (when (eq (car negpowlist) 't) 1484 (setq negpowlist (cdr negpowlist)) 1485 (go there)) 1486 (setq signa (checksigntm (power a 2))) 1487 (when (eq signa '$zero) 1488 (return (nonconstquadenum negpowlist p c b x))) 1489 (setq ea-1 (inv a)) 1490 there 1491 (setq result 0 1492 controlpow (caar negpowlist) 1493 coef (cadar negpowlist)) 1494 (when (zerop controlpow) 1495 ;; integrate(sqrt(R)). 1496 ;; I don't think we can normally get here. 1497 (setq result (augmult (mul coef 1498 (numn (add (mul p 1//2) 1//2) 1499 c b a x))) 1500 count 1) 1501 (go loop)) 1502 jump1 1503 ;; Handle integrate(sqrt(R^(2*pow-1))/x),x 1504 (setq res1 (den1numn pow c b a x)) 1505 (when (equal controlpow 1) 1506 (setq result (add result (augmult (mul coef res1))) 1507 count 2) 1508 (go loop)) 1509 jump2 1510 ;; Handle integrate(sqrt(R^(2*pow-1))/x^2,x) 1511 (unless (= p 1) 1512 ;; integrate(sqrt(R^(2*pow-1))/x^2,x) 1513 ;; 1514 ;; We can use integration by parts to get 1515 ;; 1516 ;; integrate(sqrt(R^(2*pow-1))/x^2,x) = 1517 ;; -R^(pow-1/2)/x 1518 ;; + (2*pow-1)*b/2*integrate(sqrt(R^(2*pow-3))/x,x) 1519 ;; + (2*pow-1)*c*integrate(sqrt(R^(2*pow-3)),x) 1520 (setq res2 1521 (add (augmult (mul -1 exp1 1522 (power (polfoo c b a x) 1523 (add pow -1//2)))) 1524 (augmult (mul b (div exp2 2) 1525 (den1numn (1- pow) c b a x))) 1526 (augmult (mul c exp2 (numn (- pow 2) c b a x)))))) 1527 (when (= p 1) 1528 ;; integrate(sqrt(R)/x^2,x) 1529 ;; 1530 ;; G&R 2.267, eq. 2 1531 ;; 1532 ;; integrate(sqrt(R)/x^2,x) = 1533 ;; -sqrt(R)/x 1534 ;; + b/2*integrate(1/x/sqrt(R),x) 1535 ;; + c*integrate(1/sqrt(R),x) 1536 ;; 1537 (setq res2 (add (augmult (mul -1 (power (polfoo c b a x) 1//2) 1538 exp1)) 1539 (augmult (mul b 1//2 (den1den1 c b a x))) 1540 (augmult (mul c (den1 c b a x)))))) 1541 (when (equal controlpow 2) 1542 (setq result (add result (augmult (mul coef res2))) 1543 count 3) 1544 (go loop)) 1545 jump3 1546 (setq count 4 1547 m 3) 1548 jump 1549 ;; The general case sqrt(R^(2*p-1))/x^m 1550 ;; 1551 ;; G&R 2.265 1552 ;; 1553 ;; integrate(sqrt(R^(2*p-1))/x^m,x) = 1554 ;; -sqrt(R^(2*p+1))/(m-1)/a/x^(m-1) 1555 ;; + (2*p-2*m+3)*b/2/(m-1)/a*integrate(sqrt(R^(2*p-3))/x^(m-1),x) 1556 ;; + (2*p-m+2)*c/(m-1)/a*integrate(sqrt(R^(2*n-3))/x^(m-2),x) 1557 ;; 1558 ;; NOTE: The p here is 2*pow-1. And we're 1559 ;; integrating R^(pow-1/2). 1560 1561 (setq partres 1562 (add (augmult (mul (div -1 (1- m)) 1563 ea-1 1564 (power x (1+ (- m))) 1565 (power (polfoo c b a x) 1566 (add (div p 2) 1)))) 1567 (augmult (mul (inv (+ m m -2)) 1568 ea-1 b 1569 (+ p 4 (* -2 m)) 1570 res2)) 1571 (augmult (mul c ea-1 1572 (+ p 3 (- m)) 1573 (inv (1- m)) res1)))) 1574 (incf m) 1575 (when (> m controlpow) 1576 (setq result (add result (augmult (mul coef partres)))) 1577 (go loop)) 1578 jump4 1579 (setq res1 res2 1580 res2 partres) 1581 (go jump) 1582 loop 1583 (setq negpowlist (cdr negpowlist)) 1584 (when (null negpowlist) (return result)) 1585 (setq coef (cadar negpowlist) 1586 controlpow (caar negpowlist)) 1587 (when (= count 4) 1588 (go jump4)) 1589 (when (= count 1) 1590 (go jump1)) 1591 (when (= count 2) 1592 (go jump2)) 1593 (go jump3)))) 1594 1595;; Like denmnumn, but a = 0. 1596(defun nonconstquadenum (negpowlist p c b x) 1597 (prog (result coef m) 1598 (cond ((equal p 1) 1599 (return (case1 negpowlist c b x)))) 1600 (setq result 0) 1601 loop 1602 (setq m (caar negpowlist) 1603 coef (cadar negpowlist)) 1604 (setq result (add result (augmult (mul coef (casegen m p c b x)))) 1605 negpowlist (cdr negpowlist)) 1606 (cond ((null negpowlist) (return result))) 1607 (go loop))) 1608 1609;; Integrate (c*x^2+b*x)^(p-1/2)/x^m 1610(defun casegen (m p c b x) 1611 (let ((exp1 (power (polfoo c b 0 x) (div p 2))) ;; exp1 = R^(p/2) 1612 (exp3 (power x (1+ (- m))))) ;; exp3 = 1/x^(m-1) 1613 (cond ((= p 1) 1614 ;; sqrt(c*x^2+b*x)/x^m 1615 (case1 (list (list m 1)) c b x)) 1616 ((zerop m) 1617 ;; (c*x^2+b*x)^(p-1/2) 1618 (case0 p c b x)) 1619 ((= m (1+ p)) 1620 ;; (c*x^2+b*x)^(p-1/2)/x^(p+1) 1621 (add (augmult (mul -1 exp1 (inv (1- m)) exp3)) 1622 (augmult (mul b 1//2 (casegen (1- m) (- p 2) c b x))) 1623 (augmult (mul c (casegen (- m 2) (- p 2) c b x))))) 1624 ((= m 1) 1625 ;; (c*x^2+b*x)^(p-1/2)/x 1626 ;; 1627 (add (augmult (mul (inv p) exp1)) 1628 (augmult (mul b 1//2 (case0 (- p 2) c b x))))) 1629 (t 1630 ;; (c*x^2+b*x)^(p-1/2)/x^m 1631 (add (augmult (mul -1 exp1 (inv (- m (1+ p))) exp3)) 1632 (augmult (mul -1 p b 1//2 (inv (- m (1+ p))) 1633 (casegen (1- m) (- p 2) c b x)))))))) 1634 1635;; Integrate things like sqrt(c*x^2+b*x))/x^m. 1636(defun case1 (negpowlist c b x) 1637 (declare (special *ec-1*)) 1638 (let ((exp1 (power c -1//2)) ;; exp1 = 1/sqrt(c) 1639 (eb-1 (inv b))) ;; eb-1 = 1/b 1640 (prog ((result 0) (controlpow (caar negpowlist)) (coef (cadar negpowlist)) 1641 m1 count res1 res2 m signc signb partres res) 1642 (setq m1 (- controlpow 2)) 1643 (when (zerop controlpow) 1644 (setq result (augmult (mul coef (case0 1 c b x))) 1645 count 1) 1646 (go loop)) 1647 jump1 1648 ;; sqrt(R)/x 1649 (when (= controlpow 1) 1650 (setq result 1651 (add result 1652 (augmult (mul coef (den1numn 1 c b 0 x)))) 1653 count 2) 1654 (go loop)) 1655 jump2 1656 ;; sqrt(R)/x^2 1657 (when (= controlpow 2) 1658 (setq result 1659 (add result 1660 (augmult (mul coef 1661 (denmnumn '(t (2 1)) 1 c b 0 x)))) 1662 count 3) 1663 (go loop)) 1664 jump3 1665 (setq signb (checksigntm (power b 2))) 1666 (when (eq signb '$zero) 1667 (setq count 5) 1668 (go jump5)) 1669 (setq count 4 1670 m 0 1671 signc (checksigntm *ec-1*)) 1672 (when (eq signc '$positive) 1673 (setq res 1674 (augmult (mul* 2 exp1 1675 `((%log) 1676 ,(add (power (mul c x) 1//2) 1677 (power (add b (mul c x)) 1//2)))))) 1678 (go jump4)) 1679 (setq res 1680 (augmult (mul* 2 exp1 1681 `((%atan) 1682 ,(power (mul c x 1683 (inv (add b (mul -1 c x)))) 1684 1//2))))) 1685 jump4 1686 (incf m) 1687 (setq res (add (augmult (mul -2 (power (polfoo c b 0 x) 1//2) 1688 eb-1 (inv (+ m m -1)) 1689 (power x (- m)))) 1690 (augmult (mul (div -2 (+ m m -1)) 1691 c (1- m) eb-1 res)))) 1692 (when (= m m1) 1693 (setq res2 res) 1694 (go jump4)) 1695 (when (= (1- m) m1) 1696 (if (null res2) 1697 (return nil)) 1698 (setq res1 res 1699 partres (add (augmult (mul -1 1700 (power (polfoo c b 0 x) 1//2) 1701 (r1m m) 1702 (power x (- m)))) 1703 (augmult (mul b 1//2 (r1m m) res1)) 1704 (augmult (mul c (r1m m) res2)))) 1705 (go on)) 1706 (go jump4) 1707 jump5 1708 (setq m controlpow) 1709 (when (zerop m) 1710 (setq partres (mul* exp1 `((%log) ,x))) 1711 (go on)) 1712 (setq partres (mul -1 exp1 (power x (- m)) (r1m m))) 1713 on 1714 (setq result (add result (augmult (mul coef partres)))) 1715 loop 1716 (setq negpowlist (cdr negpowlist)) 1717 (when (null negpowlist) (return result)) 1718 (setq coef (cadar negpowlist) 1719 controlpow (caar negpowlist)) 1720 (when (= count 5) (go jump5)) 1721 (when (= count 1) (go jump1)) 1722 (when (= count 2) (go jump2)) 1723 (when (= count 3) (go jump3)) 1724 (setq m1 (- controlpow 2)) 1725 (when (= m1 m) 1726 (setq res2 res1)) 1727 (go jump4)))) 1728 1729(defun r1m (m) 1730 (div 1 m)) 1731 1732;; Integrate (c*x^2+b*x)^(p-1/2) 1733(defun case0 (power c b x) 1734 (let ((exp1 '((rat simp) 1 4)) 1735 (exp2 (add b (mul 2 c x))) 1736 (exp3 (power c '((rat simp) -3 2))) 1737 (exp4 (add (mul 2 c x) (mul -1 b)))) 1738 ;; exp1 = 1/4 1739 ;; exp2 = b+2*c*x 1740 ;; exp3 = 1/c^(3/2) 1741 ;; exp4 = 2*c*x-b 1742 (declare (special *ec-1*)) 1743 (prog (signc p result) 1744 (setq signc (checksigntm *ec-1*) 1745 p 1) 1746 ;; sqrt(c*x^2+b*x) 1747 ;; 1748 ;; This could be handled by numn. Why don't we? 1749 (when (eq signc '$positive) 1750 (setq result 1751 (add (augmult (mul exp1 *ec-1* exp2 1752 (power (polfoo c b 0 x) 1//2))) 1753 (augmult (mul* b b '((rat) -1 8) 1754 exp3 1755 `((%log) 1756 ,(add exp2 1757 (mul 2 1758 (power c 1//2) 1759 (power (polfoo c b 0 x) 1//2))))))))) 1760 (when (eq signc '$negative) 1761 (setq result 1762 (add (augmult (mul exp1 *ec-1* exp4 1763 (power (polfoo (mul -1 c) b 0 x) 1//2))) 1764 (augmult (mul* b b '((rat) 1 8) 1765 exp3 1766 `((%asin) ,(mul (inv b) exp4))))))) 1767 loop 1768 (when (equal power p) (return result)) 1769 (incf p 2) 1770 1771 ;; integrate(sqrt(R^(2*n+1)),x) = 1772 ;; (2*c*x+b)/4/(n+1)/c*sqrt(R^(2*n+1)) 1773 ;; + (2*n+1)*del/8/(n+1)/c*integrate(sqrt(R^(2*n-1)),x) 1774 1775 (setq result (add (augmult (mul 1//2 *ec-1* (inv (1+ p)) exp2 1776 (power (polfoo c b 0 x) 1777 (div p 2)))) 1778 (augmult (mul p b b '((rat simp) -1 4) 1779 *ec-1* (inv (1+ p)) result)))) 1780 (go loop)))) 1781 1782;; Integrate R^(p-1/2)/x, p >= 1. 1783(defun den1numn (p c b a x) 1784 (cond ((= p 1) 1785 ;; integrate(sqrt(R)/x,x) 1786 ;; 1787 ;; G&R 2.267 eq. 1 1788 ;; 1789 ;; integrate(sqrt(R)/x,x) = 1790 ;; sqrt(R) 1791 ;; + a*integrate(1/x/sqrt(R),x) 1792 ;; + b/2*integrate(1/sqrt(R),x) 1793 (add (power (polfoo c b a x) 1//2) 1794 (augmult (mul a (den1den1 c b a x))) 1795 (augmult (mul b 1//2 (den1 c b a x))))) 1796 (t 1797 ;; General case 1798 ;; 1799 ;; G&R 2.265 1800 ;; 1801 ;; integrate(sqrt(R^(2*p-1)/x,x) = 1802 ;; R^(p-1/2)/(2*p-1) 1803 ;; + b/2*integrate(sqrt(R^(2*p-3)),x) 1804 ;; + a*integrate(sqrt(2^(2*p-3))/x,x) 1805 (add (augmult (mul (power (polfoo c b a x) 1806 (add p -1//2)) 1807 (inv (+ p p -1)))) 1808 (augmult (mul a (den1numn (+ p -1) c b a x))) 1809 (augmult (mul b 1//2 (numn (+ p -2) c b a x))))))) 1810 1811;; L is a list of expressions that INTIRA should be applied to. 1812;; Sum up the results of applying INTIRA to each. 1813(defun distrint (l x) 1814 (addn (mapcar #'(lambda (e) 1815 (let ((ie (intira e x))) 1816 (if ie 1817 ie 1818 `((%integrate simp) ,e ,x)))) 1819 l) 1820 t)) 1821