| blob | a405937cf1ac3b3e3e125d30a484562b3af4ba4b |
1 ;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancements. ;;;;;
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 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
15 ;;; this is the definite integration package.
16 ;; defint does definite integration by trying to find an
17 ;;appropriate method for the integral in question. the first thing that
18 ;;is looked at is the endpoints of the problem.
19 ;;
20 ;; i(grand,var,a,b) will be used for integrate(grand,var,a,b)
22 ;; References are to "Evaluation of Definite Integrals by Symbolic
23 ;; Manipulation", by Paul S. Wang,
24 ;; (http://www.lcs.mit.edu/publications/pubs/pdf/MIT-LCS-TR-092.pdf;
25 ;; a better copy might be: https://maxima.sourceforge.io/misc/Paul_Wang_dissertation.pdf)
26 ;;
27 ;; nointegrate is a macsyma level flag which inhibits indefinite
28 ;;integration.
29 ;; abconv is a macsyma level flag which inhibits the absolute
30 ;;convergence test.
31 ;;
32 ;; $defint is the top level function that takes the user input
33 ;;and does minor changes to make the integrand ready for the package.
34 ;;
35 ;; next comes defint, which is the function that does the
36 ;;integration. it is often called recursively from the bowels of the
37 ;;package. defint does some of the easy cases and dispatches to:
38 ;;
39 ;; dintegrate. this program first sees if the limits of
40 ;;integration are 0,inf or minf,inf. if so it sends the problem to
41 ;;ztoinf or mtoinf, respectively.
42 ;; else, dintegrate tries:
43 ;;
44 ;; intsc1 - does integrals of sin's or cos's or exp(%i var)'s
45 ;; when the interval is 0,2 %pi or 0,%pi.
46 ;; method is conversion to rational function and find
47 ;; residues in the unit circle. [wang, pp 107-109]
48 ;;
49 ;; ratfnt - does rational functions over finite interval by
50 ;; doing polynomial part directly, and converting
51 ;; the rational part to an integral on 0,inf and finding
52 ;; the answer by residues.
53 ;;
54 ;; zto1 - i(x^(k-1)*(1-x)^(l-1),x,0,1) = beta(k,l) or
55 ;; i(log(x)*x^(x-1)*(1-x)^(l-1),x,0,1) = psi...
56 ;; [wang, pp 116,117]
57 ;;
58 ;; dintrad- i(x^m/(a*x^2+b*x+c)^(n+3/2),x,0,inf) [wang, p 74]
59 ;;
60 ;; dintlog- i(log(g(x))*f(x),x,0,inf) = 0 (by symmetry) or
61 ;; tries an integration by parts. (only routine to
62 ;; try integration by parts) [wang, pp 93-95]
63 ;;
64 ;; dintexp- i(f(exp(k*x)),x,a,inf) = i(f(x+1)/(x+1),x,0,inf)
65 ;; or i(f(x)/x,x,0,inf)/k. First case hold for a=0;
66 ;; the second for a=minf. [wang 96-97]
67 ;;
68 ;;dintegrate also tries indefinite integration based on certain
69 ;;predicates (such as abconv) and tries breaking up the integrand
70 ;;over a sum or tries a change of variable.
71 ;;
72 ;; ztoinf is the routine for doing integrals over the range 0,inf.
73 ;; it goes over a series of routines and sees if any will work:
74 ;;
75 ;; scaxn - sc(b*x^n) (sc stands for sin or cos) [wang, pp 81-83]
76 ;;
77 ;; ssp - a*sc^n(r*x)/x^m [wang, pp 86,87]
78 ;;
79 ;; zmtorat- rational function. done by multiplication by plog(-x)
80 ;; and finding the residues over the keyhole contour
81 ;; [wang, pp 59-61]
82 ;;
83 ;; log*rat- r(x)*log^n(x) [wang, pp 89-92]
84 ;;
85 ;; logquad0 log(x)/(a*x^2+b*x+c) uses formula
86 ;; i(log(x)/(x^2+2*x*a*cos(t)+a^2),x,0,inf) =
87 ;; t*log(a)/sin(t). a better formula might be
88 ;; i(log(x)/(x+b)/(x+c),x,0,inf) =
89 ;; (log^2(b)-log^2(c))/(2*(b-c))
90 ;;
91 ;; batapp - x^(p-1)/(b*x^n+a)^m uses formula related to the beta
92 ;; function [wang, p 71]
93 ;; there is also a special case when m=1 and a*b<0
94 ;; see [wang, p 65]
95 ;;
96 ;; sinnu - x^-a*n(x)/d(x) [wang, pp 69-70]
97 ;;
98 ;; ggr - x^r*exp(a*x^n+b)
99 ;;
100 ;; dintexp- see dintegrate
101 ;;
102 ;; ztoinf also tries 1/2*mtoinf if the integrand is an even function
103 ;;
104 ;; mtoinf is the routine for doing integrals on minf,inf.
105 ;; it too tries a series of routines and sees if any succeed.
106 ;;
107 ;; scaxn - when the integrand is an even function, see ztoinf
108 ;;
109 ;; mtosc - exp(%i*m*x)*r(x) by residues on either the upper half
110 ;; plane or the lower half plane, depending on whether
111 ;; m is positive or negative.
112 ;;
113 ;; zmtorat- does rational function by finding residues in upper
114 ;; half plane
115 ;;
116 ;; dintexp- see dintegrate
117 ;;
118 ;; rectzto%pi2 - poly(x)*rat(exp(x)) by finding residues in
119 ;; rectangle [wang, pp98-100]
120 ;;
121 ;; ggrm - x^r*exp((x+a)^n+b)
122 ;;
123 ;; mtoinf also tries 2*ztoinf if the integrand is an even function.
128 exp
132 ;;;rsn* is in comdenom. does a ratsimp of numerator.
133 ;expvar
135 ;impvar
137 context
138 ;;LIMITP T Causes $ASKSIGN to do special things
139 ;;For DEFINT like eliminate epsilon look for prin-inf
140 ;;take realpart and imagpart.
141 ))
146 "When @code{true}, definite integration tries to find poles in the integrand
149 ;; Currently, if true, $solvetrigwarn is set to true. No additional
150 ;; debugging information is displayed.
155 nil
156 "When NIL, print a message that the principal value of the integral has
160 nil
161 "When non-NIL, a divergent integral will throw to `divergent.
169 ;; Set to true when OSCIP-VAR returns true in DINTEGRATE.
187 ;; Distribute $defint over equations, lists, and matrices.
219 (merror (intl:gettext "defint: variable of integration must be a simple or subscripted variable.~%defint: found ~M") ivar))
248 ;; If antideriv can't do it, returns nil
249 ;; use limit to evaluate every answer returned by antideriv.
253 ;;;Hack the expression up for exponentials.
256 ;; Is this expr a candidate for SININT ?
265 t)))
266 ;;ERF type things go here.
271 t))))))
291 ;; We have some subscripted variable that's not our variable
292 ;; (I think), so it's deg-lessp.
293 ;;
294 ;; FIXME: Is this the best way to handle this? Are there
295 ;; other cases we're mising here?
296 t)))
306 ;; This routine tries to take a limit a couple of ways.
312 ans
315 ivar
316 0
325 ;; Test whether fun2 is inverse of fun1 at val.
331 ;; integration change of variable
339 ())
341 ;; This is a hack! If nv is of the form b*x^n+a, we can
342 ;; solve the equation manually instead of using solve.
343 ;; Why? Because solve asks us for the sign of yx and
344 ;; that's bogus.
346 ;; Solve yx = b*x^n+a, for x. Any root will do. So we
347 ;; have x = ((yx-a)/b)^(1/n).
349 d
355 ))))
372 ;; d: original variable (ivar) as a function of 'yx
373 ;; ind: boolean flag
374 ;; nv: new variable ('yx) as a function of original variable (ivar)
383 ;; converts limits of integration to values for new variable 'yx
386 ;; rewrites exp, the integrand in terms of ivar, the
387 ;; integrand in terms of 'yx, and returns the new
388 ;; integrand.
393 ll1 ul1)
402 ;; wrapper around limit, returns nil if
403 ;; limit not found (nounform returned), or undefined ($und or $ind)
410 ans))))
428 ;; D (not used? -- cwh)
441 ;; Look for integrals which involve log and exp functions.
442 ;; Maxima has a special algorithm to get general results.
453 %tan %sinh %cosh %tanh
454 %log %asin %acos %atan
455 %cot %acot %sec
456 %asec %csc %acsc
458 ;; Call ANTIDERIV with logabs disabled,
459 ;; because the Risch algorithm assumes
460 ;; the integral of 1/x is log(x), not log(abs(x)).
461 ;; Why not just assume logabs = false within RISCHINT itself?
462 ;; Well, there's at least one existing result which requires
463 ;; logabs = true in RISCHINT, so try to make a minimal change here instead.
504 ans)
507 ans)
522 ;;Should be a PROG inside of unwind-protect, but Multics has a compiler
523 ;;bug wrt. and I want to test this code now.
531 ;;;This seems((and (and (eq ul '$inf)
532 ;;;fairly losing (setq exp (subin (m+ ll ivar) exp))
533 ;;; (setq ll 0.))
534 ;;; (ztoinf exp ivar)))
539 ;; ((and (and (equal ul 1.)
540 ;; (equal ll 0.)) (zto1 exp)))
548 ;;;NOT COMPLETE for sin's and cos's.
570 ;; Call ANTIDERIV with logabs disabled,
571 ;; because the Risch algorithm assumes
572 ;; the integral of 1/x is log(x), not log(abs(x)).
573 ;; Why not just assume logabs = false within RISCHINT itself?
574 ;; Well, there's at least one existing result which requires
575 ;; logabs = true in RISCHINT, so try to make a minimal change here instead.
581 ;;;Recursion stopper
606 ;; adds up integrals of ranges between each pair of poles.
607 ;; checks if whole thing is divergent as limits of integration approach poles.
609 ;;; calling $logcontract causes antiderivative of 1/(1-x^5) to blow up
610 ;; (setq anti-deriv (cond ((involve anti-deriv '(%log))
611 ;; ($logcontract anti-deriv))
612 ;; (t anti-deriv)))
622 ivar))))
625 ;;Return section.
635 ;; I think this takes the pole-list and replaces $MINF with -PRIN-INF
636 ;; and $INF with PRIN-INF. The lower and upper integration limits
637 ;; (ll, ul) can also be modified to be -PRIN-INF and PRIN-INF. These
638 ;; special values are used in TAKE-PRINCIPAL.
656 ;; Assumes EXP is a rational expression with no polynomial part and
657 ;; converts the finite integration to integration over a half-infinite
658 ;; interval. The substitution y = (x-a)/(b-x) is used. Equivalently,
659 ;; x = (b*y+a)/(y+1).
660 ;;
661 ;; (I'm guessing CV means Change Variable here.)
664 ;; FIXME! This is a hack. We apply the transformation with
665 ;; symbolic limits and then substitute the actual limits later.
666 ;; That way method-by-limits (usually?) sees a simpler
667 ;; integrand.
668 ;;
669 ;; See Bugs 938235 and 941457. These fail because $FACTOR is
670 ;; unable to factor the transformed result. This needs more
671 ;; work (in other places).
675 ;; If the limit is a number, use $substitute so we simplify
676 ;; the result. Do we really want to do this?
684 ()))
686 ;; Integrate rational functions over a finite interval by doing the
687 ;; polynomial part directly, and converting the rational part to an
688 ;; integral from 0 to inf. This is evaluated via residues.
691 ;; PQR divides the rational expression and returns the quotient
692 ;; and remainder
699 ;; No polynomial part
702 ans
705 ;; Only polynomial part
708 ;; A non-zero quotient and remainder. Combine the results
709 ;; together.
713 ans))
718 ;; I think this takes a rational expression E, and finds the
719 ;; polynomial part. A cons is returned. The car is the quotient and
720 ;; the cdr is the remainder.
744 ;;;If leadop isn't plus don't do anything.
753 t)
759 t)
767 ;; May reorder the limits LL and UL so that LL <= UL. (See
768 ;; ORDER-LIMITS.) Hence, this function also returns the possibly
769 ;; updated values of LL and UL as additional values.
775 ())
808 ll ul))
822 ;; *Z* is a "zero parameter" for this package.
823 ;; Its also used to transform.
824 ;; limit(exp,var,val,dir) -- limit(exp,tvar,0,dir)
831 ;; EPSILON is used in principal value code to denote the familiar
832 ;; mathematical entity.
837 ;;; PRIN-INF Is a special symbol in the principal value code used to
838 ;;; denote an end-point which is proceeding to infinity.
849 ;; Order the limits LL and UL so that LL <= UL, as expected. Of
850 ;; course, this changes the sign of the integrand (in EXP), so that's
851 ;; also updated as well. Since the order can be changed, the possibly
852 ;; updated values of LL and UL are returned as additional values of
853 ;; this function.
867 ;; We have minf <= ll = ul <= inf
868 )
870 ;; We have minf <= ll < ul = inf
871 )
873 ;; We have minf = ll < ul < inf
874 ;;
875 ;; Now substitute
876 ;;
877 ;; ivar -> -ivar
878 ;; ll -> -ul
879 ;; ul -> inf
880 ;;
881 ;; so that minf < ll < ul = inf
887 ;; We have minf <= ul < ll
888 ;;
889 ;; Now substitute
890 ;;
891 ;; exp -> -exp
892 ;; ll <-> ul
893 ;;
894 ;; so that minf <= ll < ul
897 t))
898 ll ul))
912 ;; Substitute a and b into integral e
913 ;;
914 ;; Looks for discontinuties in integral, and works around them.
915 ;; For example, in
916 ;;
917 ;; integrate(x^(2*n)*exp(-(x)^2),x) ==>
918 ;; -gamma_incomplete((2*n+1)/2,x^2)*x^(2*n+1)*abs(x)^(-2*n-1)/2
919 ;;
920 ;; the integral has a discontinuity at x=0.
921 ;;
928 ivar a b)))))
942 total))))))
944 ;; look for terms with a negative exponent
945 ;;
946 ;; recursively traverses exp in order to find discontinuities such as
947 ;; erfc(1/x-x) at x=0
959 ;; returns list of places where exp might be discontinuous in ivar.
960 ;; list begins with ll and ends with ul, and include any values between
961 ;; ll and ul.
962 ;; return '$no or '$unknown if no discontinuities found.
979 ;; (multiplicity (cdar dummy)) ;; not used
984 ;; Carefully substitute the integration limits A and B into the
985 ;; expression E.
990 ;; Note: MAKE-DEFINT-ASSUMPTIONS may reorder the limits A
991 ;; and B, but I (rtoy) don't think that's should ever
992 ;; happen because the limits should already be in the
993 ;; correct order when this function is called. We don't
994 ;; check for that, though.
1003 ;; Try easy substitutions. Return NIL if we can't.
1006 ;; Infinite limits aren't easy
1007 nil)
1011 %gamma_incomplete %expintegral_ei)))
1013 ;; It's easy if we have a polynomial. I (rtoy) think
1014 ;; it's also easy if the denominator is free of the
1015 ;; integration variable and also if the expression
1016 ;; doesn't involve inverse functions.
1017 ;;
1018 ;; gamma_incomplete and expintegral_ie
1019 ;; included because of discontinuity in
1020 ;; gamma_incomplete(0, exp(%i*x)) and
1021 ;; expintegral_ei(exp(%i*x))
1022 ;;
1023 ;; XXX: Are there other cases we've forgotten about?
1024 ;;
1025 ;; So just try to substitute the limits into the
1026 ;; expression. If no errors are produced, we're done.
1031 ;; no-err-sub has returned T. An error was catched.
1032 nil)
1046 ;; check for divergent integral
1058 ;;;This function works only on things with ATAN's in them now.
1060 ;; POLES-IN-INTERVAL doesn't know about the poles of tan(x). Call
1061 ;; trigsimp to convert tan into sin/cos, which POLES-IN-INTERVAL
1062 ;; knows how to handle.
1063 ;;
1064 ;; XXX Should we fix POLES-IN-INTERVAL instead?
1065 ;;
1066 ;; XXX Is calling trigsimp too much? Should we just only try to
1067 ;; substitute sin/cos for tan?
1068 ;;
1069 ;; XXX Should the result try to convert sin/cos back into tan? (A
1070 ;; call to trigreduce would do it, among other things.)
1073 ;;POLES -> ((mlist) ((mequal) ((%atan) foo) replacement) ......)
1074 ;;We can then use $SUBSTITUTE
1079 ul-ans)
1081 nil)))
1095 ivar ll ul)))
1112 llist)))))))))))))))
1127 ;; There are no calls where fn1 is ever equal to
1128 ;; 'mtorat. Hence this case is never true. What is
1129 ;; this testing for?
1136 nil)
1155 ;; e is of form poly(x)*exp(m*%i*x)
1156 ;; s is degree of denominator
1157 ;; adds e to *bptu* or *bptd* according to sign of m
1174 ;; Check term is of form poly(x)*exp(m*%i*x)
1175 ;; n is degree of denominator.
1185 ;; Set updn to T if the coefficient of IVAR in the
1186 ;; polynomial is known to be positive. Otherwise set to NIL.
1187 ;; (What does updn really mean?)
1198 term)
1233 ;; exponentialize sin and cos
1239 '$%i
1261 exp)
1277 ;; computes integral from minf to inf for expressions of the form
1278 ;; exp(%i*m*x)*r(x) by residues on either the upper half
1279 ;; plane or the lower half plane, depending on whether
1280 ;; m is positive or negative. [wang p. 77]
1281 ;;
1282 ;; exponentializes sin and cos before applying residue method.
1283 ;; can handle some expressions with poles on real line, such as
1284 ;; sin(x)*cos(x)/x.
1294 ivar
1295 nil)
1298 ;;;Above tests for applicability of this method.
1307 ;; call to $expand can create rational terms
1308 ;; with no exp(m*%i*x)
1312 ;;;Function RIB sets up the values of BPTU and BPTD
1319 ;;;If there is BPTU then CSEMIUP must succeed.
1320 ;;;Likewise for BPTD.
1324 ;; The original integrand was already
1325 ;; determined to have no poles by initial-analysis.
1326 ;; If individual terms of the expansion have poles, the poles
1327 ;; must cancel each other out, so we can ignore them.
1330 ;; if integral of ratterms is divergent, ratans is nil,
1331 ;; and mtosc returns nil
1348 t)
1350 t)
1358 t)
1360 t)
1365 s nc dc
1367 nn-var dn-var)
1388 sd new-sd))
1389 result))
1392 sn nil)
1394 sd nil))
1397 ;;;
1398 ;;;New section.
1416 s
1419 ivar)))
1429 findout
1433 on
1456 s)
1460 ;; It's very important here to bind VAR
1461 ;; because the PLOG simplifier checks
1462 ;; for VAR. Without this, the
1463 ;; simplifier converts plog(-x) to
1464 ;; log(x)+%i*%pi, which messes up the
1465 ;; keyhole routine.
1467 d
1468 ivar)))
1471 ;; Let's not signal an error here. Return nil so that we
1472 ;; eventually return a noun form if no other algorithm gives
1473 ;; a result.
1476 nil)))
1478 ;;(setq *dflag* nil)
1507 nn-var dn-var)
1534 sd new-sd))
1535 result))
1571 ;; Test to see if exp is of the form p(x)*f(exp(x)). If so, set pp to
1572 ;; be p(x) and set pe to f(exp(x)).
1585 s
1588 ivar)))
1595 ;; We have an integral of the form p(x)*F(exp(x)), where
1596 ;; p(x) is a polynomial.
1598 ;; No polynomial
1602 ;; These limits must exist for the integral to converge.
1605 ;; This only handles the case when the F(z) is a
1606 ;; rational function.
1609 ;; If we get here, F(z) is not a rational function.
1610 ;; We transform it using the substitution x=log(y)
1611 ;; which gives us an integral of the form
1612 ;; p(log(y))*F(y)/y, which maxima should be able to
1613 ;; handle.
1616 ;; Give up. We don't know how to handle this.
1618 en
1632 exp)))
1634 ;;; given (b*x+a)^n+c returns (a b n c)
1642 nil)
1644 ;;;get high degree of poly
1646 ;;;diff down to linear.
1648 ;;;all the way to constant.
1651 ;;;get rid of factorial from differentiation.
1654 ;;;Sees if can be expressed as (a*x+b)^n + part freeof x.
1664 s)
1689 n
1690 d)))))
1693 on
1748 e nil)
1762 e))))))
1768 ;; Compute the integral(n/d,x,0,inf) by computing the negative of the
1769 ;; sum of residues of log(-x)*n/d over the poles of n/d inside the
1770 ;; keyhole contour. This contour is basically an disk with a slit
1771 ;; along the positive real axis. n/d must be a rational function.
1776 ;; Ok if not on the positive real axis.
1781 t)
1792 ;; Look at an expression e of the form sin(r*x)^k, where k is an
1793 ;; integer. Return the list (1 r k). (Not sure if the first element
1794 ;; of the list is always 1 because I'm not sure what partition is
1795 ;; trying to do here.)
1809 ;; Look at an expression e of the form sin(r*x) and return r.
1820 ;; integrate(a*sc(r*x)^k/x^n,x,0,inf).
1823 ;; Get the argument of the involved trig function.
1826 ;; I don't think this needs to be special.
1827 #+nil
1829 ;; Replace (1-cos(arg)^2) with sin(arg)^2.
1831 ;(m+t 1. (m- (m^t `((%cos) ,ivar) 2.)))
1832 ;; The code from above generates expressions with
1833 ;; a missing simp flag. Furthermore, the
1834 ;; substitution has to be done for the complete
1835 ;; argument of the trig function. (DK 02/2010)
1838 exp))
1843 ;; n is the power of the denominator.
1845 ;; The simple case.
1851 ;; Do this for a sum of such terms.
1853 c))))))))))
1855 ;; We have an integral of the form sin(r*x)^k/x^n. C is the list (1 r k).
1856 ;;
1857 ;; The substitution y=r*x converts this integral to
1858 ;;
1859 ;; r^(n-1)*integral(sin(y)^k/y^n,y,0,inf)
1860 ;;
1861 ;; (If r is negative, we need to negate the result.)
1862 ;;
1863 ;; The recursion Wang gives on p. 87 has a typo. The second term
1864 ;; should be subtracted from the first. This formula is given in G&R,
1865 ;; 3.82, formula 12.
1866 ;;
1867 ;; integrate(sin(x)^r/x^s,x) =
1868 ;; r*(r-1)/(s-1)/(s-2)*integrate(sin(x)^(r-2)/x^(s-2),x)
1869 ;; - r^2/(s-1)/(s-2)*integrate(sin(x)^r/x^(s-2),x)
1870 ;;
1871 ;; (Limits are assumed to be 0 to inf.)
1872 ;;
1873 ;; This recursion ends up with integrals with s = 1 or 2 and
1874 ;;
1875 ;; integrate(sin(x)^p/x,x,0,inf) = integrate(sin(x)^(p-1),x,0,%pi/2)
1876 ;;
1877 ;; with p > 0 and odd. This latter integral is known to maxima, and
1878 ;; it's value is beta(p/2,1/2)/2.
1879 ;;
1880 ;; integrate(sin(x)^2/x^2,x,0,inf) = %pi/2*binomial(q-3/2,q-1)
1881 ;;
1882 ;; where q >= 2.
1883 ;;
1885 ;; Compute sign(r)*r^(n-1)*integrate(sin(y)^k/y^n,y,0,inf)
1887 c
1893 recursion)
1894 ;; Recursion failed. Return the integrand
1895 ;; The following code generates expressions with a missing simp flag
1896 ;; for the sin function. Use better simplifying code. (DK 02/2010)
1897 ; (let ((integrand (div (pow `((%sin) ,(m* r ivar))
1898 ; k)
1899 ; (pow ivar n))))
1901 k)
1906 ;; integrate(sin(x)^n/x^2,x,0,inf) = pi/2*binomial(n-3/2,n-1).
1907 ;; Express in terms of Gamma functions, though.
1913 ;; integrate(sin(x)^n/x,x,0,inf) = beta((n+1)/2,1/2)/2, for n odd and
1914 ;; n > 0.
1917 half%pi)
1920 ;; This implements the recursion for computing
1921 ;; integrate(sin(y)^l/y^k,y,0,inf). (Note the change in notation from
1922 ;; the above!)
1928 ;; Integral diverges if l-(k-1) < 0.
1931 ;; If l + k is not even, return NIL. (Is this the right
1932 ;; thing to do?)
1933 nil)
1935 ;; We have integrate(sin(y)^l/y^2). Use sevn to evaluate.
1938 ;; We have integrate(sin(y)^l/y,y)
1944 ;; j = k-2, i = (k-1)*(k-2)
1945 ;;
1946 ;;
1947 ;; The main recursion:
1948 ;;
1949 ;; i(sin(y)^l/y^k)
1950 ;; = l*(l-1)/(k-1)/(k-2)*i(sin(y)^(l-2)/y^k)
1951 ;; - l^2/(k-1)/(k-1)*i(sin(y)^l/y^(k-2))
1959 ;; Returns the fractional part of a?
1963 ;; Why do we return 0 if a is a number? Perhaps we really
1964 ;; mean integer?
1967 ;; If we're here, this basically assumes a is a rational.
1968 ;; Compute the remainder and return the result.
1976 ;; Doesn't appear to be used anywhere in Maxima. Not sure what this
1977 ;; was intended to do.
1978 #+nil
1989 ;;; THE FOLLOWING FUNCTION IS FOR TRIG FUNCTIONS OF THE FOLLOWING TYPE:
1990 ;;; LOWER LIMIT=0 B A MULTIPLE OF %PI SCA FUNCTION OF SIN (X) COS (X)
1991 ;;; B<=%PI2
1996 ;; means there was no error
1999 ; returns cons of (integer_part . fractional_part) of a
2001 ;; I think we really want to compute how many full periods are in a
2002 ;; and the remainder.
2007 ; returns cons of (integer_part . fractional_part) of a
2009 ;; I think we really want to compute how many full periods are in a
2010 ;; and the remainder.
2017 ;; Return the integer part of r.
2019 ;; r - fpart(r)
2023 ;;;Try making exp(%i*ivar) --> yy, if result is rational then do integral
2024 ;;;around unit circle. Make corrections for limits of integration if possible.
2044 rat-form
2045 ivar)))
2050 ;;; Do integrals of sin and cos. this routine makes sure lower limit
2051 ;;; is zero.
2053 ;; integrate(e,var,a,b)
2069 ivar e)
2070 cc))
2080 ;; Exit if b or a is not a constant or if the integrand
2081 ;; doesn't appear to have a period of 2 pi.
2084 ;; Multiples of 2*%pi in limits.
2091 ;; The integral is not over a full period (2*%pi) or multiple
2092 ;; of a full period.
2094 ;; Wang p. 111: The integral integrate(f(x),x,a,b) can be
2095 ;; written as:
2096 ;;
2097 ;; n * integrate(f,x,0,2*%pi) + integrate(f,x,0,c)
2098 ;; - integrate(f,x,0,d)
2099 ;;
2100 ;; for some integer n and d >= 0, c < 2*%pi because there exist
2101 ;; integers p and q such that a = 2 * p *%pi + d and b = 2 * q
2102 ;; * %pi + c. Then n = q - p.
2104 ;; Compute q and c for the upper limit b.
2111 ;; Compute p and d for the lower limit a.
2113 ;; avoid an extra trip around the circle - helps skip principal values
2118 ;; Compute -integrate(f,x,0,d)
2123 ;; Compute n = q - p (stored in nzp2)
2125 out
2126 ;; Compute n*integrate(f,x,0,2*%pi)
2128 ;; n = 0
2131 ;; n is not zero, so compute
2132 ;; integrate(f,x,0,2*%pi)
2134 ;; Unable to compute integrate(f,x,0,2*%pi)
2136 ;; Compute integrate(f,x,0,c)
2140 ;; All three pieces succeeded.
2143 ;; Less than 1 full period, so intsc can integrate it.
2144 ;; Apply the substitution to make the lower limit 0.
2145 ;; This is last resort because substitution often causes intsc to fail.
2147 limit-diff ivar))
2148 ;; nothing worked
2151 ;; integrate(sc, var, 0, b), where sc is f(sin(x), cos(x)).
2152 ;; calls intsc with a wrapper to just return nil if integral is divergent,
2153 ;; rather than generating an error.
2158 nil
2159 ans)))
2161 ;; integrate(sc, ivar, 0, b), where sc is f(sin(x), cos(x)). I (rtoy)
2162 ;; think this expects b to be less than 2*%pi.
2165 0
2172 ;; Partition the integrand SC into the factors that do not
2173 ;; contain VAR (the car part) and the parts that do (the cdr
2174 ;; part).
2179 ;; integrate(sc, ivar, 0, b), where sc is f(sin(x), cos(x)).
2181 ;; Determine if sc is a product of sin's and cos's.
2185 ;; We have a product of sin's and cos's. We handle some
2186 ;; special cases here.
2188 ;; Wang p. 110, formula (1):
2189 ;; integrate(sin(x)^m*cos(x)^n, x, 0, %pi/2) =
2190 ;; gamma((m+1)/2)*gamma((n+1)/2)/2/gamma((n+m+2)/2)
2193 ;; Wang p. 110, near the bottom, says
2194 ;;
2195 ;; int(f(sin(x),cos(x)), x, 0, %pi) =
2196 ;; int(f(sin(x),cos(x)) + f(sin(x),-cos(x)),x,0,%pi/2)
2217 ;; Wang, p. 111 says
2218 ;;
2219 ;; int(f(sin(x),cos(x)),x,0,3*%pi/2) =
2220 ;; int(f(sin(x),cos(x)),x,0,%pi)
2221 ;; + int(f(-sin(x),-cos(x)),x,0,%pi/2)
2225 ;; We don't have a product of sin's and cos's.
2230 dn*)
2235 ;;;Is careful about substitution of limits where the denominator may be zero
2236 ;;;because of various assumptions made.
2257 nil
2258 ans)))
2264 nil
2265 ans)))
2267 ;; Determine whether E is of the form sin(x)^m*cos(x)^n and return the
2268 ;; list (m n).
2272 m n)
2291 ;; Says whether (m^t value exponent) has at least one real branch.
2292 ;; Only works for values of 1 and -1 now. Returns $yes $no
2293 ;; $unknown.
2304 ;; Compute beta((m+1)/2,(n+1)/2)/2.
2310 ;;Seems like Guys who call this don't agree on what it should return.
2323 ;; Check e for an expression of the form x^kk*(b*x^n+a)^l. If it
2324 ;; matches, Return the two values kk and (list l a n b).
2329 ;; We have f(x)^y. Look to see if f(x) has the desired
2330 ;; form. Then f(x)^y has the desired form too, with
2331 ;; suitably modified values.
2332 ;;
2333 ;; XXX: Should we ask for the sign of f(x) if y is not an
2334 ;; integer? This transformation we're going to do requires
2335 ;; that f(x)^y be real.
2337 e
2342 ;; Got a match. Adjust kk and cc appropriately.
2344 cc
2359 ;;(DEFUN BATAP (E)
2360 ;; (PROG (K C L)
2361 ;; (COND ((NOT (BATA0 E)) (RETURN NIL))
2362 ;; ((AND (EQUAL -1. (CADDDR C))
2363 ;; (EQ ($askSIGN (SETQ K (m+ 1. K)))
2364 ;; '$pos)
2365 ;; (EQ ($askSIGN (SETQ L (m+ 1. (CAR C))))
2366 ;; '$pos)
2367 ;; (ALIKE1 (CADR C)
2368 ;; (m^ UL (CADDR C)))
2369 ;; (SETQ E (CADR C))
2370 ;; (EQ ($askSIGN (SETQ C (CADDR C))) '$pos))
2371 ;; (RETURN (M// (m* (m^ UL (m+t K (m* C (m+t -1. L))))
2372 ;; `(($BETA) ,(SETQ NN* (M// K C))
2373 ;; ,(SETQ DN* L)))
2374 ;; C))))))
2377 ;; Integrals of the form i(log(x)^m*x^k*(a+b*x^n)^l,x,0,ul). There
2378 ;; are two cases to consider: One case has ul>0, b<0, a=-b*ul^n, k>-1,
2379 ;; l>-1, n>0, m a nonnegative integer. The second case has ul=inf, l < 0.
2380 ;;
2381 ;; These integrals are essentially partial derivatives of the Beta
2382 ;; function (i.e. the Eulerian integral of the first kind). Note
2383 ;; that, currently, with the default setting intanalysis:true, this
2384 ;; function might not even be called for some of these integrals.
2385 ;; However, this can be palliated by setting intanalysis:false.
2404 ;; We have i(x^kk/(d+cc*x^r)^s,x,0,inf) =
2405 ;; beta(aa,bb)/(cc^aa*d^bb*r). Compute this, and then
2406 ;; differentiate it m times to get the log term
2407 ;; incorporated.
2412 0
2413 m)))
2417 r))))
2426 ;; The result looks like B(k/n,l) ( ... ).
2427 ;; Perhaps, we should just $factor, instead of
2428 ;; pulling out beta like this.
2430 beta
2431 ($fullratsimp
2439 ivar m)
2441 beta))))))))))))
2444 ;;; If e is of the form given below, make the obvious change
2445 ;;; of variables (substituting ul*x^(1/n) for x) in order to reduce
2446 ;;; e to the usual form of the integrand in the Eulerian
2447 ;;; integral of the first kind.
2448 ;;; N. B: The old version of ZTO1 completely ignored this
2449 ;;; substitution; the log(x)s were just thrown in, which,
2450 ;;; of course would give wrong results.
2453 ;; Parse e
2457 ;; e=x^k*(a+b*x^n)^l
2459 c
2471 ;; Wang p. 71 gives the following formula for a beta function:
2472 ;;
2473 ;; integrate(x^(k-1)/(c*x^r+d)^s,x,0,inf)
2474 ;; = beta(a,b)/(c^a*d^b*r)
2475 ;;
2476 ;; where a = k/r > 0, b = s - a > 0, s > k > 0, r > 0, c*d > 0.
2477 ;;
2478 ;; This function matches this and returns k-1, d, r, c, a, b. And
2479 ;; also checks that all the conditions hold. If not, NIL is returned.
2480 ;;
2486 c
2503 ;; Handles beta integrals.
2511 nil)
2514 c
2515 ;; e = x^k*(d+c*x^al)^l.
2520 new-k)
2528 ;; Compute exp(d)*gamma((c+1)/b)/b/a^((c+1)/b). In essence, this is
2529 ;; the value of integrate(x^c*exp(d-a*x^b),x,0,inf).
2540 ;; Evaluates the contour integral of GRAND around the unit circle
2541 ;; using residues.
2548 ;; Is pt stricly inside the unit circle?
2554 ;; Is pt on the unit circle? (Use
2555 ;; the cached value computed
2556 ;; above.)
2577 ;; Wang 81-83. Unfortunately, the pdf version has page 82 as all
2578 ;; black, so here is, as best as I can tell, what Wang is doing.
2579 ;; Fortunately, p. 81 has the necessary hints.
2580 ;;
2581 ;; First consider integrate(exp(%i*k*x^n),x) around the closed contour
2582 ;; consisting of the real axis from 0 to R, the arc from the angle 0
2583 ;; to %pi/(2*n) and the ray from the arc back to the origin.
2584 ;;
2585 ;; There are no poles in this region, so the integral must be zero.
2586 ;; But consider the integral on the three parts. The real axis is the
2587 ;; integral we want. The return ray is
2588 ;;
2589 ;; exp(%i*%pi/2/n) * integrate(exp(%i*k*(t*exp(%i*%pi/2/n))^n),t,R,0)
2590 ;; = exp(%i*%pi/2/n) * integrate(exp(%i*k*t^n*exp(%i*%pi/2)),t,R,0)
2591 ;; = -exp(%i*%pi/2/n) * integrate(exp(-k*t^n),t,0,R)
2592 ;;
2593 ;; As R -> infinity, this last integral is gamma(1/n)/k^(1/n)/n.
2594 ;;
2595 ;; We assume the integral on the circular arc approaches 0 as R ->
2596 ;; infinity. (Need to prove this.)
2597 ;;
2598 ;; Thus, we have
2599 ;;
2600 ;; integrate(exp(%i*k*t^n),t,0,inf)
2601 ;; = exp(%i*%pi/2/n) * gamma(1/n)/k^(1/n)/n.
2602 ;;
2603 ;; Equating real and imaginary parts gives us the desired results:
2604 ;;
2605 ;; integrate(cos(k*t^n),t,0,inf) = G * cos(%pi/2/n)
2606 ;; integrate(sin(k*t^n),t,0,inf) = G * sin(%pi/2/n)
2607 ;;
2608 ;; where G = gamma(1/n)/k^(1/n)/n.
2609 ;;
2617 ;; Ok, we have cos(b*x^n) or sin(b*x^n), and we set e = (n
2618 ;; b)
2620 ;; n = 1. Give up. (Why not divergent?)
2621 nil)
2626 ;; s is the sign of b. Give up if it's zero.
2627 nil)
2629 ;; Give up if n-1 <= 0. (Why give up? Isn't the
2630 ;; integral divergent?)
2631 nil)
2633 ;; We can apply our formula now. g = gamma(1/n)/n/b^(1/n)
2640 e)))))))
2643 ;; this is the second part of the definite integral package
2648 ())
2650 ivar t)
2672 nil)
2674 t)))))
2684 n
2687 plm))))
2699 plm
2700 factors
2701 pl
2702 rl
2703 pl1
2704 rl1))))
2712 ans)
2714 ;; AFAICT, this call to PLOG doesn't need
2715 ;; to bind VAR.
2717 ans))))
2724 ;; AFAICT, this call to PLOG doesn't need
2725 ;; to bind VAR. An example where this is
2726 ;; used is
2727 ;; integrate(log(x)^2/(1+x^2),x,0,1) =
2728 ;; %pi^3/16.
2731 n)
2732 d
2733 plm)))))
2737 0
2740 ;; Handle integral(n(x)/d(x)*log(x)^m,x,0,inf). n and d are
2741 ;; polynomials.
2759 ans))))
2783 n
2787 ;;Makes a new poly such that np(x)-np(x+2*%i*%pi)=p(x).
2788 ;;Constructs general POLY of degree one higher than P with
2789 ;;arbitrary coeff. and then solves for coeffs by equating like powers
2790 ;;of the varibale of integration.
2791 ;;Can probably be made simpler now.
2799 ;;;Now, poly(x)-poly(x+2*%i*%pi)=p(x), P is the original poly.
2808 -1
2817 ;; Integrate(p(x)*f(exp(x))/g(exp(x)),x,minf,inf) by applying the
2818 ;; transformation y = exp(x) to get
2819 ;; integrate(p(log(y))*f(y)/g(y)/y,y,0,inf). This should be handled
2820 ;; by dintlog.
2827 ;; This implements Wang's algorithm in Chapter 5.2, pp. 98-100.
2828 ;;
2829 ;; This is a very brief description of the algorithm. Basically, we
2830 ;; have integrate(R(exp(x))*p(x),x,minf,inf), where R(x) is a rational
2831 ;; function and p(x) is a polynomial.
2832 ;;
2833 ;; We find a polynomial q(x) such that q(x) - q(x+2*%i*%pi) = p(x).
2834 ;; Then consider a contour integral of R(exp(z))*q(z) over a
2835 ;; rectangular contour. Opposite corners of the rectangle are (-R,
2836 ;; 2*%i*%pi) and (R, 0).
2837 ;;
2838 ;; Wang shows that this contour integral, in the limit, is the
2839 ;; integral of R(exp(x))*q(x)-R(exp(x))*q(x+2*%i*%pi), which is
2840 ;; exactly the integral we're looking for.
2841 ;;
2842 ;; Thus, to find the value of the contour integral, we just need the
2843 ;; residues of R(exp(z))*q(z). The only tricky part is that we want
2844 ;; the log function to have an imaginary part between 0 and 2*%pi
2845 ;; instead of -%pi to %pi.
2847 ;; We have R(exp(x))*p(x) represented as p(x)*pe(exp(x))/d(exp(x)).
2853 ;; At this point denom-exponential has converted d(exp(x)) to the
2854 ;; polynomial d(z), where z = exp(x).
2857 pe))
2859 ;; Find the poles of the denominator. denom-exponential is the
2860 ;; denominator of R(x).
2861 ;;
2862 ;; It seems as if polelist returns a list of several items.
2863 ;; The first element is a list consisting of the pole and (z -
2864 ;; pole). We don't care about this, so we take the rest of the
2865 ;; result.
2868 ;; The imaginary part is nonzero,
2869 ;; or the realpart is negative.
2873 ;; The realpart is not zero.
2875 ;; Not sure what this does.
2877 ;; No roots at all, so return
2881 ;; We have simple roots or roots in REGION1
2884 ;; The cadr of pl is the list of the simple poles of
2885 ;; denom-exponential. Take the log of them to find the
2886 ;; poles of the original expression. Then compute the
2887 ;; residues at each of these poles and sum them up and put
2888 ;; the result in B. (If no simple poles set B to 0.)
2894 ;; I think this handles the case of poles outside the
2895 ;; regions. The sum of these residues are placed in C.
2899 temp)))
2904 ;; We have the repeated poles of deonom-exponential, so we
2905 ;; need to convert them to the actual pole values for
2906 ;; R(exp(x)), by taking the log of the value of poles.
2911 ;; Compute the residues at all of these poles and sum
2912 ;; them up.
2915 poles))
2930 ;; Check to see if each term in exp that is of the form exp(k*x) has
2931 ;; an integer value for k.
2947 ;; Test to see if exp is of the form f(exp(x)), and if so, replace
2948 ;; exp(x) with 'z*. That is, basically return f(z*).
2951 exp)
2955 ivar))
2958 ;; Does e really need to be special here? Seems to be ok without
2959 ;; it; testsuite works.
2960 #+nil
2963 e)
2994 ;; I think this is looking at f(exp(x)) and tries to find some
2995 ;; rational function R and some number k such that f(exp(x)) =
2996 ;; R(exp(k*x)).
3013 ;; Integration by parts.
3014 ;;
3015 ;; integrate(u(x)*diff(v(x),x),x,a,b)
3016 ;; |b
3017 ;; = u(x)*v(x)| - integrate(v(x)*diff(u(x),x))
3018 ;; |a
3019 ;;
3021 ;;;SINCE ONLY CALLED FROM DINTLOG TO get RID OF LOGS - IF LOG REMAINS, QUIT
3025 nil)
3032 nil)
3039 ;; integrate(f(exp(k*x)),x,a,b), where f(z) is rational.
3040 ;;
3041 ;; See Wang p. 96-97.
3042 ;;
3043 ;; If the limits are minf to inf, we use the substitution y=exp(k*x)
3044 ;; to get integrate(f(y)/y,y,0,inf)/k. If the limits are 0 to inf,
3045 ;; use the substitution s+1=exp(k*x) to get
3046 ;; integrate(f(s+1)/(s+1),s,0,inf).
3052 ;; If we can integrate it directly, do so and take the
3053 ;; appropriate limits.
3054 )
3056 ;; ans is the list (f(x) exp(k*x)).
3059 ;; Use the substitution s + 1 = exp(k*x). The
3060 ;; integral becomes integrate(f(s+1)/(s+1),s,0,inf)
3063 ;; Use the substitution y=exp(k*x) because the
3064 ;; limits are minf to inf.
3066 ;; Apply the substitution and integrate it.
3069 ;; integrate(log(g(x))*f(x),x,0,inf)
3077 exp
3078 ivar))
3080 ;; Make the substitution y=1/x. If the integrand has
3081 ;; exactly the same form, the answer has to be 0.
3087 ;; It's easy if we have the antiderivative.
3088 ;; but intsubs sometimes gives results containing %limit
3090 ;; Ok, the easy cases didn't work. We now try integration by
3091 ;; parts. Set ANS to f(x).
3094 ;; Bad. f(x) contains a log term, so we give up.
3099 ivar ll ul)))
3100 ;; The arg of the log function is the same as the
3101 ;; integration variable. We can do something a little
3102 ;; simpler than integration by parts. We have something
3103 ;; like f(x)*log(x). Consider f(x)*x^z. If we
3104 ;; differentiate this wrt to z, the integrand becomes
3105 ;; f(x)*log(x)*x^z. When we evaluate this at z = 0, we
3106 ;; get the desired integrand.
3107 ;;
3108 ;; So we need f(0) to be 0 at 0. If we can integrate
3109 ;; f(x)*x^z, then we differentiate the result and
3110 ;; evaluate it at z = 0.
3113 ;; Try integration by parts.
3116 ;; Compute diff(e,ivar,n) at the point pt.
3120 ;;; GGR and friends
3122 ;; MAYBPC returns (COEF EXPO CONST)
3123 ;;
3124 ;; This basically picks off b*x^n+a and returns the list
3125 ;; (b n a).
3132 ;; At this point, e is of the form (a n b)
3137 ;; If we're here, b is complex, and n > 0. zn = imagpart(b).
3138 ;;
3139 ;; Set ivar to the same sign as zn.
3144 ;; zd = exp(ivar*%i*%pi*(1+nd)/(2*n). (ZD is special!)
3147 ;; Return zn, n, a, zd.
3153 ;; We're here if realpart(b) <= 0, and n >= 0. Then return -b, n, a.
3156 ;; Integrate x^m*exp(b*x^n+a), with realpart(m) > -1.
3157 ;;
3158 ;; See Wang, pp. 84-85.
3159 ;;
3160 ;; I believe the formula Wang gives is incorrect. The derivation is
3161 ;; correct except for the last step.
3162 ;;
3163 ;; Let J = integrate(x^m*exp(%i*k*x^n),x,0,inf), with real k.
3164 ;;
3165 ;; Consider the case for k < 0. Take a sector of a circle bounded by
3166 ;; the real line and the angle -%pi/(2*n), and by the radii, r and R.
3167 ;; Since there are no poles inside this contour, the integral
3168 ;;
3169 ;; integrate(z^m*exp(%i*k*z^n),z) = 0
3170 ;;
3171 ;; Then J = exp(-%pi*%i*(m+1)/(2*n))*integrate(R^m*exp(k*R^n),R,0,inf)
3172 ;;
3173 ;; because the integral along the boundary is zero except for the part
3174 ;; on the real axis. (Proof?)
3175 ;;
3176 ;; Wang seems to say this last integral is gamma(s/n/(-k)^s) where s =
3177 ;; (m+1)/n. But that seems wrong. If we use the substitution R =
3178 ;; (y/(-k))^(1/n), we end up with the result:
3179 ;;
3180 ;; integrate(y^((m+1)/n-1)*exp(-y),y,0,inf)/(n*k^((m+1)/n).
3181 ;;
3182 ;; or gamma((m+1)/n)/k^((m+1)/n)/n.
3183 ;;
3184 ;; Note that this also handles the case of
3185 ;;
3186 ;; integrate(x^m*exp(-k*x^n),x,0,inf);
3187 ;;
3188 ;; where k is positive real number. A simple change of variables,
3189 ;; y=k*x^n, gives
3190 ;;
3191 ;; integrate(y^((m+1)/n-1)*exp(-y),y,0,inf)/(n*k^((m+1)/n))
3192 ;;
3193 ;; which is the same form above.
3212 ;; e = (m b n a). That is, the integral is of the form
3213 ;; x^m*exp(b*x^n+a). I think we want to compute
3214 ;; gamma((m+1)/n)/b^((m+1)/n)/n.
3215 ;;
3216 ;; FIXME: If n > m + 1, the integral converges. We need
3217 ;; to check for this.
3219 e
3222 ;; The derivation only holds if b is purely real or
3223 ;; purely imaginary. Give up if it's not.
3225 ;; Check for convergence. If b is complex, we need n -
3226 ;; m > 1. If b is real, we need b < 0.
3235 ;; If we're here, b must be negative.
3238 ;; Complex b. Take the imaginary part
3240 n a))
3242 ;; FIXME: Why do we set %emode here? Shouldn't we just
3243 ;; bind it? And why do we want it bound to T anyway?
3244 ;; Shouldn't the user control that? The same goes for
3245 ;; dosimp.
3246 ;;(setq $%emode t)
3252 ;; Match x^m*exp(b*x^n+a). If it does, return (list m b n a).
3258 ;; We're looking at something like exp(f(ivar)). See if it's
3259 ;; of the form b*x^n+a, and return (list 0 b n a). (The 0 is
3260 ;; so we can graft something onto it if needed.)
3265 ;; E should be the product of exactly 2 terms
3267 ;; Check to see if one of the terms is of the form
3268 ;; ivar^p. If so, make sure the realpart of p > -1. If
3269 ;; so, check the other term has the right form via
3270 ;; another call to ggr1.
3281 ;; Both terms have the right form and nn* contains the ivar of
3282 ;; the exponential term. Put dn* as the car of nn*. The
3283 ;; result is something like (m b n a) when we have the
3284 ;; expression x^m*exp(b*x^n+a).
3288 ;; Match b*x^n+a. If a match is found, return the list (a n b).
3289 ;; Otherwise, return NIL
3302 ;; Match b*x^n. Return the list (n b) if found or NIL if not.
3312 ;; nn* should be the value of var. This is only called by bx**n with
3313 ;; the second arg of var.
3324 ;;; given (b*x^n+a)^m returns (m a n b)
3331 nil)
3338 ;;;Catches Unfactored forms.
3341 ivar)
3354 ;;;Is E = VAR raised to some power? If so return power or 0.
3360 ;;;Is E = VAR1 raised to some power? If so return power.
3381 ans)))))
3409 ;; Note: This doesn't seem be called from anywhere.
3434 ivar)
3462 ;;;Looks at the IMAGINARY part of a log and puts it in the interval 0 2*%pi.
3465 ;; We take the $rectform above to make sure that the log is
3466 ;; expanded out for the situations where simplifying plog itself
3467 ;; doesn't do it. This should probably be considered a bug in the
3468 ;; plog simplifier and should be fixed there.
3481 ;;; Temporary fix for a lacking in taylor, which loses with %i in denom.
3482 ;;; Besides doesn't seem like a bad thing to do in general.
3486 ;; Multiply the denominator by it's conjugate to get rid of
3487 ;; %i.
3492 ;; If the new denominator still contains %i, just give up.
3494 exp
3495 new-exp)))
3498 ;;; LL and UL must be real otherwise this routine return $UNKNOWN.
3499 ;;; Returns $no $unknown or a list of poles in the interval (ll ul)
3500 ;;; for exp w.r.t. ivar.
3501 ;;; Form of list ((pole . multiplicity) (pole1 . multiplicity) ....)
3518 ;; this clause handles cases where we can't find the exact roots,
3519 ;; but we know that they occur outside the interval of integration.
3520 ;; example: integrate ((1+exp(t))/sqrt(t+exp(t)), t, 0, 1);
3536 ;; (multiplicity (cdar dummy)) (not used? -- cwh)
3546 pole-list))))))))))))
3549 ;;;Returns $YES if there is no pole and $NO if there is one.
3562 ;;;Takes care of forms that the ratio test fails on.
3578 (%einvolve-var
3581 ivar))
3615 ;; real values for ll and ul; place can be imaginary.
3622 '$no
3627 ;; returns true or nil
3629 ;; real values for ll and ul; place can be imaginary.
3645 ;; If *failures is set, we may have missed some roots.
3646 ;; We still return the roots that we have found.
3651 rootlist)
3658 rootlist)))))))))
3661 ;;; Is x > y. X or Y can be $MINF or $INF, zeroA or zeroB.
3701 nil)
3704 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3705 ;;;
3706 ;;; Integrate Definite Integrals involving log and exp functions. The algorithm
3707 ;;; are taken from the paper "Evaluation of CLasses of Definite Integrals ..."
3708 ;;; by K.O.Geddes et. al.
3709 ;;;
3710 ;;; 1. CASE: Integrals generated by the Gamma function.
3711 ;;;
3712 ;;; inf
3713 ;;; /
3714 ;;; [ w m s - m - 1
3715 ;;; I t log (t) expt(- t ) dt = s signum(s)
3716 ;;; ]
3717 ;;; /
3718 ;;; 0
3719 ;;; !
3720 ;;; m !
3721 ;;; d !
3722 ;;; (--- (gamma(z))! )
3723 ;;; m !
3724 ;;; dz ! w + 1
3725 ;;; !z = -----
3726 ;;; s
3727 ;;;
3728 ;;; The integral converges for:
3729 ;;; s # 0, m = 0, 1, 2, ... and realpart((w+1)/s) > 0.
3730 ;;;
3731 ;;; 2. CASE: Integrals generated by the Incomplete Gamma function.
3732 ;;;
3733 ;;; inf !
3734 ;;; / m !
3735 ;;; [ w m s d s !
3736 ;;; I t log (t) exp(- t ) dt = (--- (gamma_incomplete(a, x ))! )
3737 ;;; ] m !
3738 ;;; / da ! w + 1
3739 ;;; x !z = -----
3740 ;;; s
3741 ;;; - m - 1
3742 ;;; s signum(s)
3743 ;;;
3744 ;;; The integral converges for:
3745 ;;; s # 0, m = 0, 1, 2, ... and realpart((w+1)/s) > 0.
3746 ;;; The shown solution is valid for s>0. For s<0 gamma_incomplete has to be
3747 ;;; replaced by gamma(a) - gamma_incomplete(a,x^s).
3748 ;;;
3749 ;;; 3. CASE: Integrals generated by the beta function.
3750 ;;;
3751 ;;; 1
3752 ;;; /
3753 ;;; [ m s r n
3754 ;;; I log (1 - t) (1 - t) t log (t) dt =
3755 ;;; ]
3756 ;;; /
3757 ;;; 0
3758 ;;; !
3759 ;;; ! !
3760 ;;; n m ! !
3761 ;;; d d ! !
3762 ;;; --- (--- (beta(z, w))! )!
3763 ;;; n m ! !
3764 ;;; dz dw ! !
3765 ;;; !w = s + 1 !
3766 ;;; !z = r + 1
3767 ;;;
3768 ;;; The integral converges for:
3769 ;;; n, m = 0, 1, 2, ..., s > -1 and r > -1.
3770 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3774 ;;; Recognize c*z^w*log(z)^m*exp(-t^s)
3786 ;;; Recognize c*z^r*log(z)^n*(1-z)^s*log(1-z)^m
3799 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3806 ;; var1 is used as a parameter for differentiation. Add var1>0 to the
3807 ;; database, to get the desired simplification of the differentiation of
3808 ;; the gamma_incomplete function.
3816 ;; The integrand matches the cases 1 and 2.
3835 ;; Case 1: Generated by the Gamma function.
3842 var1
3847 ; derivate the involved hypergeometric
3848 ; functions.
3853 ;; Case 2: Generated by the Incomplete Gamma function.
3869 ;; Case 3: Generated by the Beta function.
3895 var2 m)
3897 var1 n)
3903 ;; Simplify result and set $gamma_expand to global value
3906 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;