.EQ delim $$ .EN .nr Pi 10 .nr Pt 1 .LP \fBTHE METHOD OF DIFFERENTIATING UNDER THE INTEGRAL SIGN\fR .SP2 .LP \fIGert Almkvist* .FS* Department of Mathematics, University of Lund, Box 118, 22100 Lund, Sweden. .FE and Doron Zeilberger**\fR .FS** Department of Mathematics and Computer Science, Drexel University, Philadelphia, PA 19104. Supported in part by NSF grant DMS8800663. .FE .SP2 .DS (Submitted to J. Symbolic Computation, April 1989; Revised: Nov 1989) .DE .SP10 .DS \fI"I could never resist a definite integral"\fR (G. H. Hardy [Co]) .DE .SP2 .I .P "One thing I never did learn was contour integration. I had learned to do integrals by various methods shown in a book that my high-school physics teacher Mr. Bader had given me. ... .P The book also showed how to differentiate parameters under the integral sign-It's a certain operation. It turns out that's not taught very much in the universities; they don't emphasize it. But I caught on how to use that method, and I used that one damn tool again and again. So because I was self-taught using that book, I had peculiar methods of doing integrals. .P The result was that, when guys at MIT or Princeton had trouble doing a certain integral, it was because they couldn't do it with the standard methods they had learned in school. If it was contour integration, they would have found it; if it was a simple series expansion, they would have found it. Then I come along and try differentiating under the integral sign, and often it worked. So I got a great reputation for doing integrals, only because my box of tools was different from everybody else's, and they had tried all their tools on it before giving the problem to me." .SP1 .DS (Richard P. Feynman[F]) .DE .SK .H 1 "Introduction" .P The \fImethod of differentiating under the integral sign\fR can be described as follows. Given a function $F(x,y)$ of $x$ and $y$, one is interested in evaluating .SP1 .EQ(1) R(x):= int from {- inf} to inf F(x,y) dy~~~. .EN .SP1 .LP By differentiating $R(x)$ once or more with respect to x, and using integration by parts or change of variable, one gets a differential equation for $R(x)$. Then one evaluates $R(x)$ and its derivatives in an "easy" point, and then solves the differential equation. .P Although termed a "method", differentiating under the integral sign could hardly have been considered more than a trick, and the examples given in the few textbooks that treat it ([Ed]) are rather simple. One notable exception is Melzak's[Me] delightful book, that devotes a whole section to this method. .P In this paper we are going to present a unified theory of differentiation under the integral sign for the important and wide class of so-called \fIholonomic functions\fR. We will recall the "slow" elimination algorithm of [Z1], and present it in a form that will hopefully lead to a fast version using Buchberger's method of Grobner bases adapted to the Weyl algebra. For the subclass of "hyperexponential" functions we will give a fast algorithm, that given $F(x,y)$, finds a linear differential equation satisfied by $R(x)$. This "fast" algorithm is a continuous analog of Zeilberger's [Z2] [Z3] (see also [W-Z]) algorithm for hypergeometric summation, that uses Gosper's[Go] algorithm as a subroutine. In Appendix 1 (which is published in [A-Z]) we give a listing of a MAPLE program that implements this algorithm, while Appendix 2 contains sample inputs and outputs. .P While the problem of \fIindefinite integration\fR in closed form is completely solved and implemented ([Ri], [No], [Tr], [DST] ), there is currently no well-developed theory of definite integration ( see [No], 3.2). The present paper can be viewed as a first step toward developing a theory of definite integration. .H 1 "The Theoretical Foundation" .P In [Z1] Bernstein's class of holonomic functions is considered, and it is shown how the methodology of holonomic functions and systems forms a natural framework for special function identities. Roughly a function $F( x sub 1 , ... , x sub n )$ is holonomic if it "satisfies as many linear differential equations with polynomial coefficients as possible" (see appendix 4, footnote 1 for the precise definition). .P The algebra of linear differential operators with polynomial coefficients in $( x sub 1 , ... , x sub n ) $ is called the \fIWeyl algebra\fR. It is generated by the "indeterminates" $x sub 1 , ... , x sub n$, $D sub 1 , ... , D sub n$, that obey the commutation relations .SP1 .LP .EQ x sub i x sub j = x sub j x sub i~~~, D sub i D sub j = D sub j D sub i~~~, x sub i D sub j = D sub j x sub i ~~~(i != j) .EN .SP1 .LP and the "uncertainty principle" .SP1 .EQ D sub i x sub i - x sub i D sub i ^=^1~~~, .EN .SP1 .LP where $1$ is the identity operator. .P Excellent expositions on the Weyl algebra and Bernstein's theory can be found in [Bj] and [Eh]. It was proved in [Z1](section 6.4), using ideas of Bernstein[Be], that if $F(x,y)$ is holonomic and $P(x,y, D sub x , D sub y )$ and $Q(x,y, D sub x , D sub y )$ are operators that generate a "holonomic ideal" then one can "eliminate" $y$. In other words, it is possible to find operators $A(x,y, D sub x , D sub y )$ and $B(x,y, D sub x , D sub y )$ such that .SP1 .EQ S(x, D sub x , D sub y ) = A(x,y, D sub x , D sub y )P(x,y, D sub x , D sub y )^+^ B(x,y, D sub x , D sub y )Q(x,y, D sub x , D sub y ) .EN .SP1 .LP is independent of $y$. It follows that $F(x,y)$ is also annihilated by $S$: $S(x, D sub x , D sub y ) F(x,y) == 0$. One can write .SP1 .EQ(2) S(x, D sub x , D sub y ) = S bar (x , D sub x ) + D sub y S' (x, D sub x , D sub y )~~~. .EN .SP1 .LP It is possible to show that it is always possible to have $S bar != 0$. It follows that .SP1 .EQ(3) S bar (x, D sub x ) F(x,y)~ =~ - D sub y S' ( x, D sub x , D sub y ) F(x,y) ~~. .EN .SP1 .LP In the applications we have in mind, we will have that $lim sub {|y| -> inf } x sup n y sup m D sub x sup s D sub y sup t F(x,y)$ $=0$ for all integers $n,m,s$ and $t$. Therefore integrating with respect to $y$ from $- inf$ to $inf$, we get .SP1 .EQ S bar ( x , D sub x ) R(x) = 0~~~. .EN .P It was observed in [Z1] that the slow elimination algorithm for the discrete case, that was an adaptation of Sylvester's dialytic elimination method from commutative algebra, goes verbatim to the continuous case. The hope was raised that Buchberger's fruitful concept of Grobner bases could be used to get a fast elimination algorithm. We now reiterate this hope and make one more observation that should make the elimination even faster. .SP1 .LP \fBObservation\fR: In order for things to work we don't need that $S'$ in (2) be independent of $y$! Our previous insistence probably made $S bar$ of higher order than necessary. .SP1 .P This leads to the following .SP1 .LP \fBResearch Problem\fR: Develop a fast algorithm, using Grobner bases in the Weyl algebra ([Ga]), that inputs two operators $P(x,y, D sub x , D sub y )$ and $Q(x,y, D sub x , D sub y )$ and outputs operators $A,B,C, S bar $ such that .SP1 .EQ(4) S bar (x, D sub x ) := AP+BQ+ D sub y C .EN .SP1 .LP only depends on $x$ and $D sub x$. .P Basically we have to find an element $S bar (x, D sub x )$ that does not involve $y$ and $D sub y$, that belongs to the "ideal" generated by $P$,$Q$ and $D sub y$. But note that this ideal is neither "left" nor "right", and therefore not an ideal at all, in the usual sense of the word. It is "left" with respect to $P$ and $Q$ and "right" with respect to $D sub y$, and so one has to develop an ambidextrous Grobner bases algorithm in the Weyl algebra. .P A major step towards this goal was recently achieved by Takayama[Ta]. .P A related notion to holonomicity is that of D-finiteness, which was considered by Lipshitz[Li]. A formal power series $f( x sub 1 , ... , x sub n )$ is \fID-finite\fR if the vector space generated by all its partial derivatives is finite dimensional over the field of rational functions in $( x sub 1 , ... , x sub n )$. Every holonomic function is \fIipso facto\fR D-finite, but the converse is not known (see [Ta]). For our present purposes, the elegant elementary approach of Lipshitz[Li] suffices and there is no need to invoke the theory of holonomic functions. .H 1 "A Detailed Example" .P In order to make the ideas clear, we will use ad hoc elimination to treat the following integral. .SP1 .EQ R(x) := int from { - inf} to inf exp( - {x sup 2} over {y sup 2} - y sup 2 ) dy~~~. .EN .SP1 .LP With $F(x,y)=exp( - x sup 2 / y sup 2 - y sup 2 ) $ we have .SP1 .LP .EQ left "{" lpile{ { D sub x F = (-2x/ y sup 2 ) F} above { D sub y F = ( 2 x sup 2 / y sup 3 - 2y )F}} .EN .SP1 .LP or .SP1 .EQ left "{" lpile { { ( y sup 2 D sub x + 2x ) F^=^0} above { ( y sup 3 D sub y + 2 y sup 4 -2 x sup 2 )F^=^0} }~~~. .EN .SP1 .LP We want to apply $int from {- inf} to inf dy $ on these equations. Therefore we want to get rid of all terms containing $y$, except those of the type $D sub y G$. This depends on the fact that .SP1 .EQ int from { - inf } to inf D sub y G dy = [G] sub { - inf} sup inf = 0~~~, .EN .SP1 .LP since we assume that everything in sight tends to zero as $y -> +- inf$. This means that after integration we can ignore all terms starting with $D sub y$ ( i.e. put $D sub y =0$). Now .SP1 .EQ y sup 3 D sub y = D sub y y sup 3 - 3 y sup 2~~~, .EN .SP1 .LP hence the second equation becomes .SP1 .EQ QF=( D sub y y sup 3 - 3 y sup 2 + 2 y sup 4 - 2 x sup 2 )F =0~~. .EN .SP1 .LP Apply $D sub x$ on the left and use the first equation .SP1 .EQ y sup 2 D sub x F = - 2xF~~~ .EN .SP1 .LP several times. We obtain .SP1 .EQ D sub x Q^F^=^D sub y D sub x y sup 3 F - 3 D sub x y sup 2 F + 2 y sup 2 D sub x y sup 2 F - 2 D sub x x sup 2 F .EN .SP1 .EQ = ^D sub y D sub x y sup 3 F - 3 (-2x)F+ 2 y sup 2 (-2x) F + - 2 D sub x x sup 2 F~~~. .EN .SP2 .EQ D sub x sup 2 Q F= D sub y D sub x sup 2 y sup 3 F + 6 D sub x x F -4 y sup 2 ( D sub x x ) F - 2 D sub x sup 2 x sup 2F .EN .SP1 .EQ = D sub y D sub x sup 2 y sup 3 F + 6 D sub x x F -4 y sup 2 (x D sub x +1 ) F - 2 D sub x sup 2 x sup 2 F .EN .SP1 .EQ = D sub y D sub x sup 2 y sup 3 F + 6 D sub x x F + 8 x sup 2 F - 4 y sup 2 F - 2 D sub x sup 2 x sup 2 F ~~. .EN .SP2 .EQ D sub x sup 3 Q F = D sub y D sub x sup 3 y sup 3 F + 6 D sub x sup 2 x F + 8 D sub x x sup 2 F - 4 D sub x y sup 2 F - 2 D sub x sup 3 x sup 2 F = .EN .SP1 .EQ = D sub y D sub x sup 3 y sup 3 F + 6 D sub x sup 2 x F + 8 D sub x x sup 2 F + 8x F - 2 D sub x sup 3 x sup 2 F ~~~. .EN .SP1 .LP Differentiating .SP1 .EQ S(x):= int from {- inf} to inf F(x,y) dy .EN .SP1 .LP with respect to $x$, we obtain .SP1 .EQ D sub x sup r S(x) = int from {- inf} to inf D sub x sup r F(x,y) dy~~. .EN .SP1 .LP Integrating $D sub x sup 3 Q F$, the first term vanishes, and we get .SP1 .EQ (6 D sub x sup 2 x + 8 D sub x x sup 2 + 8x - 2 D sub x sup 3 x sup 2 ) S(x) == 0~~. .EN .SP1 .LP Using the rules of "multiplication" in the Weyl algebra, this factors to: .SP1 .EQ ( x D sub x + 3)( D sub x -2)( D sub x +2)S(x) = 0~~. .EN .SP1 .LP This differential equation has the general solution .SP1 .EQ S(x)= C sub 1 e sup 2x + C sub 2 e sup { -2x} + C sub 3 y sub 3~~, .EN .SP1 .LP where $y sub 3 $ satisfies .SP1 .EQ y'' sub 3 - 4 y sub 3 = A/ x sup 3 ~~ for~~ A != 0~~. .EN .SP1 .LP One can show that $y sub 3$ cannot be bounded at $0$. But .SP1 .EQ S(0)= int from {- inf} to inf e sup { - y sup 2 } dy = sqrt pi ~~, .EN .SP1 .LP hence $C sub 3 =0$. Assume that $x > 0$. Then $S(x) -> 0$ as $x -> inf $ implies that also $C sub 1 =0$. It follows that .SP1 .EQ S(x)= sqrt pi e sup -2x~~. .EN .SP1 .P The above example has only pedagogical value, since it is done much easier by performing the substitution $t= y-x/y$ on the "obvious" integral $int from {- inf} to inf exp( - t sup 2 ) = pi$ (see Appendix 4, footnote 2) or by an argument that combines differentiation under the integral sign and substitution,that is given in p. 220 of Edwards book [Ed] (reproduced in Appendix 4, footnote 3) . The advantage of the present method is that it is purely mechanical and has been programmed. We really hope that the experts on Grobner bases will develop a fast elimination algorithm. .H 1 "The Class of Hyperexponential Functions" .P The discrete analog of the exponential function $ e sup cx$ is the geometric sequence $ a sub n = c sup n $, that has the property that $a sub n+1 / a sub n = $ is constant. An important generalization of the class of geometric sequences is the class of \fIhypergeometric\fR sequences, that consist of sequences {$a sub n $} for which $a sub n+1 / a sub n $ is a \fIrational\fR function of $n$. Hypergeometric sequences are important because they are the terms of \fIhypergeometric series\fR, by which most of the classical special functions can be expressed. .P Since the exponential function can be defined as a function $f(x)$ for which $f'(x)/f(x)$ is a constant, it is natural to define, as the continuous analog of hypergeometric sequences, an \fIhyperexponential\fR function to be a function $f(x)$ for which the logarithmic derivative $f'(x)/f(x)$ is a rational function of $x$. .P A function of several variables $F( x sub 1 , ... , x sub n )$ is hyperexponential if the logarithmic derivatives $D sub x sub i F / F$, with respect to every variable $x sub i$, are rational functions of $( x sub 1 , ... , x sub n )$. .P Consider equation (3) where $F(x,y)$ is hyperexponential. Since a hyperexponential function is obviously D-finite (see [Li]), it follows from a result of Lipshitz ([Li], Lemma 3) that there exist operators $S bar$ and $S'$ such that (3) holds. Now set .SP1 .EQ(5) G(x,y):=~ -^ S' ( x, D sub x , D sub y ) F(x,y)~~. .EN .SP1 .LP It is easily proved by induction that for hyperexponential $F$, $D sub x sup i D sub j sup j F /F$ is a rational function for every nonnegative $i$ and $j$. It follows that $G(x,y)$ is a multiple of $F(x,y)$ by a certain rational function, and thus hyperexponential on its own right. Now equation (3) can be rewritten as .SP1 .EQ(6) S bar (x, D sub x ) F(x,y)~ =~ D sub y G(x,y)~~. .EN .SP1 Thus we can forget about $S'$ and concentrate our attention on finding the winning pair $S bar$ and $G$, now that we know that they exist. Suppose you already knew, or guessed, $S bar$. Then (6) asserts that $S bar ( x, D sub x )F(x,y)$ can be integrated in closed form with respect to $y$, and Risch's algorithm [Ri], that is implemented on many symbolic packages, can be used to find $G(x,y)$. However, since $S bar$ is not given, but has to be found, we chose to develop our own algorithm for indefinite integration, that is, in some sense, a special case of Risch's algorithm, but that we will be able to extend to find both $S bar$ and $G$ of (6). This algorithm, that is a continuous analog of Gosper's ([Go],see also [La], [GKP]) decision procedure for indefinite hypergeometric summation, will now be presented. .H 1 "A Decision Procedure For Indefinite Hyperexponential Integration or: Gosper's Algorithm Translated To The Continuous" .P This section follows closely Gosper's[Go] classical paper, and the numbering of the equations in this section correspond to that of Gosper's paper. To distinguish them from the other equations in this paper, they are prefixed with a "G". .SP1 .LP \fBInput:\fR A hyperexponential function $a(y)$, i.e. a rational function $P(y)/Q(y)$ such that $a'(y)/a(y)=P(y)/Q(y)$. (Throughout this section $a'(y)$ means $D sub y a(y) = d over dy a(y)$. ) .SP1 .LP \fBOutput:\fR: A hyperexponential function $S(y)$ such that $S'(y)=a(y)$ if such a hyperexponential function $S(y)$ exists, or the statement "$a(y)$ does not have a hyperexponential indefinite integral". .SP1 .P We first find polynomials $p(y), q(y), r(y)$ such that .SP1 .EQ(G5) (a(y)/p(y))' over (a(y)/p(y)) = q(y) over r(y)~~, .EN .SP1 .LP where $q(y)$ and $r(y)$ satisfy the condition .SP1 .EQ(G6) gcd(r(y), q(y)- jr'(y))=1 ~~for~every ~ non-negative~integer~j. .EN .SP1 .LP It is always possible to do this for if $gcd(r(y),q(y)-jr'(y))=g(y)$, then this common factor can be eliminated by .SP1 .EQ(G6`) ~ .EN .EQ p(y) ~ size +2 { <- }~ p(y) g(y) sup j~, ~ r(y) size +2 { <-} ~ r(y)/g(y)~, ~q(y) size +2 { <-} ~ j(r(y)/g(y))' ^+^(q(y)^-^j^r'(y)^)/g(y)~. .EN .SP1 .LP The values of $j$ for which such $g(y)'s$ exist can be readily detected as the non-negative integer roots of the resultant of $r(y)$ and $q(y)-jr'(y)$ with respect to $y$. .P We now write .SP1 .EQ(G7) S(y)=r(y)f(y)a(y)/p(y)~~~, .EN .SP1 where $f(y)$ is to be determined. By using $S'(y)=a(y)$, .SP1 .EQ f(y)=(S(y)/S'(y))p(y)/r(y)~~~, .EN .SP1 .LP so $f(y)$ is a rational function whenever $S(y)/S'(y)$ is, (i.e. whenever $S(y)$ is hyperexponential.) By substituting Eq. (G7) into $S'(y)=a(y)$, we get .SP1 .EQ a(y)=^r(y)f(y)(a(y)/p(y))'~+~r(y)f'(y)(a(y)/p(y))~+~r'(y)f(y)(a(y)/p(y)) .EN .SP1 .LP Multiplying this through by $p(y)/a(y)$, and using Eq. (G5), we have .SP1 .EQ(G8) p(y)=(q(y)+r'(y))f(y)+r(y)f'(y)~~~, .EN .SP1 the differential equation for $f(y)$. .SP1 .LP \fBTheorem:\fR If $S'(y)/S(y)$ is a rational function of $y$, then $f(y)$ is a polynomial. .SP1 .LP \fBProof\fR: We already know that $f(y)$ is a rational function when $S'(y)/S(y)$ is, so suppose .SP1 .EQ(G9) f(y)= c(y) over d(y)~~, .EN .SP1 .LP where $d(y)$ is a polynomial of positive degree, and .SP1 .EQ(G9a) gcd(c(y),d(y))=1~~~. .EN .SP1 .LP Substituting (G9) in (G8) we get .SP1 .LP .EQ(G10) qcd+r'cd+r(c'd-cd')= d sup 2 p~~. .EN .SP1 .LP It is always possible to write .SP1 .LP .EQ(G11) d(y) = A(y) sup k d bar (y)~~~, where~ gcd(A(y), d bar (y)) =1,~ gcd( A(y), A'(y))=1~~~. .EN .SP1 .LP Substituting (G11) in (G10) it is readily seen that $A(y)$ divides $r(y)$, i.e.. for some polynomial $r bar (y)$ .SP1 .LP .EQ(G12) r(y)= A(y) r bar (y)~~~~~. .EN .SP1 .LP Substituting (G12) in (G10) easily yields that $A(y)$ is also a factor of $q(y)~ -~ (k-1) r'(y)$, but this contradicts the assumption (G6), and thus $d(y)$ must be constant, QED. .P All that remains is to look for a polynomial $f(y)$ satisfying equation (G8), given $p(y)$, $q(y)$, and $r(y)$. To do this, we must find an upper bound $k$, for the degree of the polynomial $f(y)$. .SP1 .LP \fBCase 1:\fR .SP1 .EQ deg(r(y)) <= deg( q(y)+ r'(y))= l~~~. .EN .SP1 .LP Setting $f(y)= c sub k y sup k + O( y sup k-1 ) $, we get $p(y)= L c sub k y sup k+l + O( y sup k+l-1 ) $, where $L$ is a nonzero constant. Since the two sides of the equation must be of equal degree, .SP1 .EQ k= deg(f(y))= deg(p(y))-l~~~. .EN .SP1 .LP \fBCase 2:\fR .SP1 .LP .EQ deg(q(y)+r'(y)) < deg(r(y))=l~~~. .EN .SP1 .LP Estimating .SP1 .EQ f(y)= c sub k y sup k + c sub k-1 y sup k-1 + O( y sup k-2 ) .EN .SP1 .LP and .SP1 .EQ f'(y)= k c sub k y sup k-1 + (k-1) c sub k-1 y sup k-2 + O( y sup k-3) .EN .SP1 .LP in Eq. (G8) leads to .SP1 .EQ p(y)= L(k) c sub k y sup k+l-1 + O( y sup k+l-2 )~~, .EN .SP1 .LP where $L(k)$ is linear in $k$ and free of $c sub k$. Let $k sub 0 $ be the root of $L(k)$. Then, the largest $k$ for which $c sub k$ can be nonzero is .SP1 .EQ k=max "{" ~ k sub 0 ~,~ deg(p(y))-l+1 ~"}", ~if~k sub 0~ is~an~integer~,or .EN .SP1 .EQ k= deg(p(y))-l+1~,~if ~not~. .EN .SP1 .LP In either case, if $k < 0$, the indefinite integral $S(y)$ is not hyperexponential; otherwise we can substitute a $kth$ degree polynomial, with $k+1$ undetermined coefficients, for $f(y)$ in Eq. (G8). Equating like powers of $y$, we get a system of linear equations whose consistency is equivalent to the existence of a hyperexponential $S(y)$ such that $S'(y)=a(y)$ and whose solution provides $f(y)$, whereupon the indefinite integral, $S(y)$, is given by Eq. (G7). .LP .H 1 "A Fast Algorithm For Definite Hyperexponential Integration" .SP1 .P Consider (1) where $F(x,y)$ is hyperexponential. By the theory of holonomic systems, we know \fIfor sure\fR that there exist $S bar (x, D sub x )$ and a hyperexponential $G(x,y)$ such that .SP1 .EQ(6) S bar (x, D sub x ) F(x,y) = D sub y G(x,y)~~. .EN .SP1 .P We now consider $x$ as an auxiliary parameter and set .SP1 .LP .EQ S bar (x, D sub x ) = sum from 0 to m a sub i (x) D sub x sup i .EN .SP1 .LP where $m$ is the guessed order of $S bar$. Letting the $ a sub i$ be unknown, we evaluate the left of (6), \fIwith the undetermined\fR $a sub i (x)$, and apply the algorithm of section 5 to find the $G (x,y)$. In addition to the unknown coefficients of $f(y)$ (that now depend on $x$) in this algorithm, we also have the unknowns $a sub i (x)$. Since they appear linearly, we can modify the above algorithm and require that the resulting system of linear equations be consistent, and when it is, we get as the solution, not only the information about $G(x,y)$, but also the successful $S bar (x, D sub x )$. We first try $m=0$ and increase $m$ until we get a solution. Note that there is no "trial and error" involved here, and we are \fIguaranteed\fR to eventually get a successful $m$. Furthermore, it is not hard to derive an a priori upper bound for $m$. .P When we use the subroutine of section 5, the left side of (6) already has the form of a polynomial times a hyperexponential function, so the starting value of $p(y)$ is not $1$, and we already have a head start. .P A listing of a MAPLE program that implements both the algorithms of section 6 and 7 is given in appendix 1, and appendix 2 contains examples of input and output. .H 1 "The method of differencing under the integral sign" .P The same idea can be used to find a linear recurrence relation satisfied by .SP1 .LP .EQ a sub n := int from {- inf} to inf F sub n (y) dy~~, .EN .LP where both $F' sub n (y) / F sub n (y)$ and $F sub n+1 (y) / F sub n (y) $ are rational functions of $n$ and $y$, that possibly depend on other variables and parameters. We try a generic recurrence operator applied to $F sub n (y)$, and enquire, by the above algorithm, when it has a hyperexponential primitive function with respect to $y$. .SP1 .LP \fBSimple example (Courtesy of L. Euler)\fR: If $F sub n (y) = y sup n e sup -y $, we have .SP1 .EQ (n+1-N) F sub n (y) = D sub y (y F sub n (y) )~~~, .EN .SP1 .LP where $N a sub n := a sub n+1$. Integrating w.r.t to $y$ from $- inf$ to $inf$ yields $(n+1-N) a sub n =0 $, i.e. $a sub n+1 = (n+1) a sub n $, and so $ a sub n = n!$. .P Appendix 3 gives a computer-generated proof of the recurrence that leads to the evaluation of Euler's beta integral as well as four other, less trivial, recurrences. .SP1 .LP .H 1 "Applications To The Theory Of Special Functions" .SP1 .P It is perhaps surprising that while many deep hypergeometric sum identities are known([Ba]), all of which can be handled by the algorithm of [Z2] and [Z3], relatively few hyperexponential integral identities of the form (1) are known, where both $F(x,y)$ and $R(x)$ are hyperexponential. In fact, the only non-trivial cases known to us are those that are obtained by applying the change of variable .SP1 .EQ t= x - a over x .EN .SP1 .LP to the integral $int from {- inf} to inf exp( - t sup 2 ) dt = sqrt pi$, and iterating (see appendix 4, footnote 2). Furthermore, the differential equation that the algorithm outputs is never, in our experience, first order, so one has to resort to comparing initial values, that is much harder to do then in the discrete case([Z3]). In the discrete case the sums are usually terminating, and checking initial conditions is trivial. .P Many definite integral formulas that one can find in integral tables can be handled by the holonomic "slow" algorithm, but few fall into the hyperexponential class. On the other hand, the present algorithm is very useful for proving theorems in the theory of special functions. All the classical families of orthogonal polynomials may be defined in terms of a generating function .SP1 .EQ(7) sum from n=0 to inf P sub n (x) y sup n~=~ A(x,y)~~~, .EN .SP1 .LP where $A(x,y)$ turns out to be hyperexponential ( except for the Jacobi polynomials). Then (7) can be rewritten as .SP1 .EQ P sub n (x) = 1 over {2 pi i } int from {|y| = rho} A(x,y) over {y sup n+1 } dy~~~. .EN .SP1 .LP Since the integrand is hyperexponential, one can apply our algorithm to find, and prove, the differential equation satisfied by the sequence {$P sub n (x)$} defined by (7). Similarly, one can use our algorithm for differencing under the integral sign, outlined in the previous section, to find and prove the linear recurrence equation satisfied by the family. Examples are given in appendices 2 and 3. .P Another application is to proving Rodrigues-type formulas( e.g. [Ho], p.9). Here we use Cauchy's formula .SP1 .EQ d sup n over {dx sup n} f(x)= n! over {(2 pi i) sup n } int from { |y-x| = rho } f(y) over { (y-x) sup n+1 } dy~~~. .EN .SP1 .LP By using the algorithms of sections 6 and 7 we can get a differential equation and a recurrence equation respectively. .P Another application of the algorithm of section 7 is for finding a pure recurrence relation in $n$ for sequences {$a sub n $} given by .SP1 .LP .EQ a sub n =~ coeff.~ of~ x sup b ~in~ P(x) sup n~~~, .EN .SP1 .LP where $P(x)$ is a Laurent polynomial. This is useful in one-dimensional random walk problems. .SP1 .LP \fBAcknowledgement\fR: We would like to thank Viggo Kann for helpful suggestions. We also thank Shalosh B. Ekhad and Sol Tre for executing our programs. Finally we would like to thank the two referees for their numerous helpful comments. .SP1 .LP \fBAppendix 1: A Maple Program Implementing The Algorithm Of Section 6\fR .SP1 .LP \fBA.1.1 The program itself\fR: .P The listing for the program is given in [A-Z]. It is available from Zeilberger via e-mail (Please send your e-mail address via regular mail). Place that program in a file, and name it WHATYOUWILL. .SP1 .LP \fBA.1.2: A Typical Input File\fR: .P We will describe how to use the program in the UNIX operating system. Readers who use other operating systems will easily make the appropriate modifications. .P Place the following input file in a file called INPUT. You should replace WHATYOUWILL by the full path name of the program file. .SP1 .DS #Beginning Of Typical Input File read `WHATYOUWILL`: expli:=1: A:=exp(-x^2/y^2-y^2): B:=exp(-2*x): ORDER:=1: duis(expli,conj,A,B,ORDER,disorcon): quit: #End of Typical Input File .DE .SP1 .LP In order to use the program, you type on your terminal .SP1 .DS maple < INPUT > OUTPUT .DE .SP1 .LP After the execution is completed you can read the output file by entering "cat OUTPUT". .P The program can find both a differential equation and a recurrence equation. The default is the former case, where \fIdisorcon\fR is not set to anything. Then A(x,y) is the integrand, B(x) is the conjectured right side (if known), in which case expli is set to 1, and the program itself finds "conj", which is the first order differential operator annihilating $B(x)$. If $B(x)$ is not known explicitly, then we don't set anything for $B$, but set \fIexpli:=0\fR. If $B(x)$ is conjectured to satisfy a certain differential equation $conj(x,D)B(x)^=^0$, then we set $conj:=conj(x,D)$. If we don't know what differential equation $B(x)$ satisfies, then we set \fIconj:=1\fR. In every case the output is the differential operator $gorem(x,D)$ and the hyperexponential function $G(x,y)$ such that .SP1 .EQ gorem(x,D)conj(x,D)A(x,y)^=^D sub y ^G(x,y)~~~. .EN .SP1 .LP In other words $S bar (x,D)^=^gorem(x,D)conj(x,D)$. ORDER is the guessed order of $gorem(x,D)$ Note that $A(x,y)$ must always be a function of $x$ and $y$, and it is always assumed that $y$ is the integration variable (so $y$ is a global variable). .P If a recurrence equation is desired, then the integrand, $A(n,y)$, must depend on $n$ and $y$, and the program finds a recurrence operator, expressed in terms of the shift operator $N$ , $Nf(n):=f(n+1)$. The above input file is the same, only now one sets \fIdisorcon:=discrete\fR. In the discrete case set expli:=0 always. .SP1 .LP \fBAppendix 2: Some Sample Inputs and Outputs For Finding Differential Equations\fR .SP1 .LP 1. With the following input the program proves .SP1 .EQ int from {- inf} to inf exp( - x sup 2 / y sup 2 - y sup 2 ) dy~=~ sqrt pi e sup -2x~~~, .EN .SP1 .LP after the appropriate initial conditions are checked. .SP1 .DS #Begin input read `WHATYOUWILL`: expli:=1: conj:=1: A:=exp(-x^2/y^2-y^2): B:=exp(-2*x): ORDER:=1: duis(expli,conj,A,B,ORDER,disorcon): quit: #End input .DE .SP1 .LP The output was: .SP1 .DS the conjectured rhs B is, exp(- 2 x) conj, the first order operator annihilating B(x) is, D + 2 Let R(x) be the integral of the following w.r.t y .DE .SP1 .EQ ~~~~~~~~~~~~~~~~exp left ( { - { x sup 2 } over { y sup 2 } - y sup 2 } right ) .EN .SP1 .DS The operator annihilating R(x) is the following times conj - 2 + D G(x,y) is the integrand times 2 1/y .DE .SP1 .LP \fB2.\fR With the following input, the program finds, and then proves, that .SP1 .EQ R(x):= int from { - inf} to inf exp( - x sup 2 / y sup 2 + a y sup 2 ) dy~~ .EN .SP1 .LP satisfies the differential equation $(4a + D sup 2 )R(x) = 0$. .SP1 .DS #Begin input read `WHATYOUWILL`: expli:=0: conj:=1: A:=exp(-x^2/y^2+a*y^2): ORDER:=2: duis(expli,conj,A,B,ORDER,disorcon): quit: #End input .DE .SP1 .LP The output was: .DS conj is 1 Let R(x) be the integral of the following w.r.t y .DE .SP1 .EQ ~~~~~~~~~~~~~~~~exp left ( { - { x sup 2 } over { y sup 2 } + a y sup 2 } right ) .EN .SP1 .DS The operator annihilating R(x) is the following times conj .DE .SP1 .EQ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 4 a + D sup 2 .EN .SP1 .DS G(x,y) is the integrand times 2 1/y .DE .SP1 .LP \fB3.\fR .SP1 .LP .DS #Begin of input read `WHATYOUWILL`: expli:=0: conj:=1: .DE .DS A:=exp(-x^6/y^4-y^2): ORDER:=3: duis(expli,conj,A,B,ORDER,disorcon): quit: #End of input .DE .SP1 .LP The output was .DS conj is 1 Let R(x) be the integral of the following w.r.t y .DE .SP1 .EQ ~~~~~~~~~~~~~~~~exp left ( { - { x sup 6 } over { y sup 4 } - y sup 2 } right ) .EN .SP1 .DS The operator annihilating R(x) is the following times conj .DE .SP1 .EQ ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 216 x sup 5 + 7 D - 12 x D sup 2 + 4 x sup 2 D .EN .DS G(x,y) is the integrand times .DE .SP1 .EQ ~~~~~~~~~~~~~~~~~~~~~~~~~~~ -54 {( 4 x sup 6 - y sup 4 + 2 y sup 6 ) x sup 5 } over { y sup 7 } .EN .SP1 .LP \fB4.\fR With the following input, the program finds, and then proves, the differential equation satisfied by the Legendre polynomials, when the later are defined by their generating function: .SP1 .EQ sum from n=0 to inf P sub n (x) t sup n ~=~ (1-2xt+ t sup 2 ) sup -1/2 ~~~~. .EN .SP1 .LP .DS #Begin of input #The diff. eq. for Legendre out of the gen. fun. read `WHATYOUWILL`: expli:=0: A:=1/(y^(n+1)*(1-2*x*y+y^2)^(1/2)): ORDER:=2: duis(expli,conj,A,B,ORDER,disorcon): quit: #End of input .DE .SP1 .LP The output was: .DS #The diff. eq. for Legendre out of the gen. fun. conj is 1 Let R(x) be the integral of the following w.r.t y .DE .SP1 .EQ ~~~~~~~~~~~~~~~~~~1 over { y sup (n+1) (1~-~2~x~y~+y sup 2 ) sup 1/2 } .EN .SP1 .DS The operator annihilating R(x) is the following times conj .DE .EQ ~~~~~~~~~~~~~~~~~~~ n ~(n~ +~ 1)~ -~ 2~ x~ D ~-~ (x~ -~ 1) (x~ + ~1)~ D sup 2 .EN .DS G(x,y) is the integrand times .DE .SP1 .EQ ~~~~~~~~~~~~~~~~~~~~~ -~ { (-~ n ~- ~1 +~ 2~ n~ x ~y + ~x~ y ~-~ n ~y sup 2 ) ~y } over { (-~ 1 ~+~ 2~ x~ y~ -~ y sup 2 )~ n } .EN .SP1 .DS words used=156843, alloc=88064, time=29.37 .DE .SP1 .LP \fB5.\fR .SP1 .LP With the following input, the program proves the Mehler formula for Hermite polynomials (e.g. [Ho], p. 24, ex. 4). .P The Mehler formula came to the limelights when Foata[Fo1] gave a beautiful combinatorial proof that inaugurated the currently very active area of combinatorial special function theory. Unlike the computer-generated proof below, that is a \fIverification\fR, Foata's proof is a \fIderivation\fR, where every part of the generating function arises naturally. Furthermore, Foata's proof lead naturally to a proof of a multi-variate generalization of the Mehler formula that was accomplished by Foata and Garsia[FG](see also [Fo2]). The theory of holonomic systems is incapable, \fIeven in principle\fR, of proving multi-variate identities with an \fIarbitrary\fR number of variables. .SP1 .DS #Begin input #The Mehler formula for Hermite polynomials read `WHATYOUWILL`: expli:=0: A:=exp((2*x*y*z-y^2*(x^2+z^2))/(2*(1-y^2)))/((1-y^2)^(1/2)*y^(n+1)): ORDER:=2: duis(expli,conj,A,B,ORDER,disorcon): quit: #End of input .DE .SP1 .LP The output was: .DS #The Mehler formula for Hermite polynomials conj is 1 Let R(x) be the integral of the following w.r.t y .DE .EQ ~~~~~~~~~~~~~~~~~~~~~~~~ {exp( {2~x~y~z~-y sup 2 ~ ( x sup 2 ~ + z sup 2 )} over {2~-~2~y sup 2 } ) } over { (1~-~ y sup 2 ) sup 1/2 ~y sup {(n~+1) } } .EN .DS The operator annihilating R(x) is the following times conj .DE .SP1 .EQ ~~~~~~~~~~~~~~~~~~~~~~~~ - ~n ~+~ x~ D~ -~ D sup 2 .EN .SP1 .DS G(x,y) is the integrand times y words used=138796, alloc=79872, time=26.87 .DE .SP1 .LP \fB6.\fR With the following input the program proves the Doetsch identity (see [Fo2]). .SP1 .DS #To Prove the Doetsch formula for Hermite polynomials(see Foata,Adv. #Appl. Math. 2(1981),250-259.) read `WHATYOUWILL`: expli:=0: .DE .DS A:=y^(-n-1)*(1+2*y^2)^(-3/2)*(1+y*x+2*y^2)*exp(y^2*x^2/(1+2*y^2)): ORDER:=2: duis(expli,conj,A,B,ORDER,disorcon): quit: .DE .SP1 .LP The output was: .SP1 .DS conj is , 1 Let R(x) be the integral of the following w.r.t y .DE .SP1 .LP .EQ ~~~~~~~~~~~~~~~~ { y sup {(-~n~-1)} (1~+~y~x~+2~y sup 2 ) exp( { y sup 2 x sup 2 } over { 1~+~ 2 y ~ sup 2 } ) } over {(1~+~2~y sup 2 ) sup 3/2 } .EN .SP1 .DS The operator annihilating R(x) is the following times conj .DE .SP1 .EQ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~-n~+x~D~-~D sup 2 .EN .SP1 .DS G(x,y) is the integrand times y words used=194380, alloc=112640, time=40.10 .DE .SP1 .LP \fB7.\fR .br .LP With the following input, the program proves the Rodrigues formula for Jacobi polynomials ([Ra], p. 257) .SP1 .DS #Begin input #Proves Rodrigues' Formula for Jacobi polynomials by finding the #diff. eq read `WHATYOUWILL`: expli:=0: A:=(y^2-1)^n*(1-y)^a*(1+y)^b*(1-x)^(-a)*(1+x)^(-b)/(y-x)^(n+1): ORDER:=2: duis(expli,conj,A,B,ORDER,disorcon): quit: #End input .DE .SP1 .LP The output was: .SP1 .DS #Proves Rodrigues' Formula for Jacobi polynomials by finding the #diff. eq conj is the identity 1 Let R(x) be the integral of the following w.r.t y .DE .SP1 .EQ ~~~~~~~~ { ( y sup 2 ~-~1) sup n ~(1~-~y) sup a ~(1+~y) sup b ~(1~-x) sup {(- ~a)}~(1~+~x) sup {(-~b)} } over { (y-~x) sup {(n~+~1)} } .EN .SP1 .DS The operator annihilating R(x) is the following times conj .DE .SP1 .EQ - (n + a + 1 + b) n + (a + a x - b + x b + 2 x) D + (- 1 + x) (1 + x) D sup 2 .EN .SP1 .DS G(x,y) is the integrand times .DE .SP1 .EQ ~~~~~~~~~~~~~~~~~~~~~ { (n~ + ~1) (-~ 1~ +~ y) ~(1~ +~ y) } over {(- ~y~ +~ x)~ } .EN .SP1 .DS words used=397108, alloc=155648, time=96.80 .DE .SP1 .LP \fB8.\fR .SP1 .LP With the following input, the program proves Euler's integral representation for the hypergeometric function $"" sub 2 F sub 1 (a,b;c;x)$ ([Ra], Th. 16, p. 47), if the latter is defined by its differential equation, or alternatively, finds the differential equation if it is defined in terms of the integral representation. .SP1 .DS #Begin of input #Proves Euler's integral representation for 2F1 read `WHATYOUWILL`: expli:=0: A:=(1-x*y)^(-a)*y^(b-1)*(1-y)^(c-b-1): ORDER:=2: duis(expli,conj,A,B,ORDER,disorcon): quit: #End of input .DE .SP1 .LP The output was: .DS #Proves Euler's representation for 2F1 conj is , 1 Let R(x) be the integral of the following w.r.t y .DE .SP1 .EQ ~~~~~~~~~~~~~~~~~~~ (1~-~x~y) sup {(-~a)}~y sup {(b~-~1)}~(1-~y~) sup {(c~-~b~-~1)} .EN .SP1 .DS The operator annihilating R(x) is the following times conj .DE .SP1 .EQ ~~~~~~~~~~~~~~~~ b~ a~ +~ (x~ +~ a~ x~ + ~x ~b~ -~ c)~ D + ~x~ (x~ -~ 1) ~D sup 2 .EN .SP1 .DS G(x,y) is the integrand times .DE .SP1 .EQ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ { a y~ (- ~1~ + ~y) } over { { - ~1~ +~ x~ y} } .EN .SP1 .DS words used=132760, alloc=77824, time=25.07 .DE .SP1 .LP \fB9.\fR .SP1 .LP Rational functions can always be integrated either by contour-integration or by partial fractions, at least in principle. Here is an example how our program finds a differential equation satisfied by such an integral. .SP1 .LP .DS #Begin input read `WHATYOUWILL`: expli:=0: A:=x^2/((x^3+y^3)*(1+y^3)): ORDER:=2: duis(expli,conj,A,B,ORDER,disorcon): quit: #End input .DE .SP1 .LP The program outputted the following differential operator. We omit $G(x,y)$. .SP1 .EQ 3~x sup 2 ~+~D~+5~x sup 3 ~D~+~x sup 4 ~D sup 2~-x~D sup 2 .EN .SP1 .LP \fBAppendix 3: Sample Inputs And Outputs For Finding Recurrence Equations\fR .SP1 .P In all the following, the same program was used, and the input files have the same form, only now the parameter "disorcon" is set to \fIdiscrete\fR , and $A$ should involve $n$ and $y$. $N$ is the shift operator in the $n$ direction: $Nf(n):=f(n+1)$. .SP1 .LP \fB1.\fR .SP1 .LP With the following input the program finds, and proves, the recurrence equation satisfied by the Legendre polynomials, if the latter are defined by the Rodrigues formula. .SP1 .DS #Beginning of input read `WHATYOUWILL`: expli:=0: A:=(y^2-1)^n/(y-x)^(n+1): ORDER:=2: duis(expli,conj,A,B,ORDER,discrete): quit: #End of input .DE .SP1 .LP The output was: .DS conj is , 1 Let R(n) be the integral of the following w.r.t y .DE .SP1 .EQ ~~~~~~~~~~~~~~~~~~~~~~ {( y sup 2 ~-~1) sup n} over {(y~-~x) sup {(n~+~1)}} .EN .SP1 .DS The operator annihilating R(n) is the following times conj .DE .SP1 .EQ - ~4~ n~ - ~4~ +~ 2~ x ~(3~ +~ 2~ n)~ N ~ +~ (-~ n~ -~ 2)~ N sup 2 .EN .SP1 .DS G(n,y) is the integrand times .DE .SP1 .EQ ~~~~~~~~~~~~~~~~~~~~ - {(y~-~1)~(y~+~1)~(-~ y sup 2 ~+~2~y~x~-~1)} over {-~y~+~x} .EN .SP1 .DS words used=161099, alloc=94208, time=31.43 .DE .SP1 .LP \fB2.\fR .SP1 .LP With the following input, the program finds the recurrence equation satisfied by the Jacobi polynomials, when defined in terms of their Rodrigues's formula. .SP1 .DS #Begin of input #Finds, and proves, the recurrence equation for the Jacobi polynomials #when they are defined in terms of the Rodrigues formula. read `WHATYOUWILL`: expli:=0: A:=(y^2-1)^n*(1-y)^a*(1+y)^b*(1-x)^(-a)*(1+x)^(-b)/(y-x)^(n+1): ORDER:=2: duis(expli,conj,A,B,ORDER,discrete): quit: #End of input .DE .SP1 .LP The output is omitted( do it on your computer!). .SP1 .LP \fB3.\fR .SP1 .LP With the following input the program finds, and proves the recurrence satisfied by the sequence $R(n):=$constant term of $(y+1+1/y) sup n$. .SP1 .LP .DS #Begin of input #Finds the recurrence for the constant term of (1/y+1+y)^n: read `WHATYOUWILL`: expli:=0: A:=(1/y+1+y)^n/y: ORDER:=2: duis(expli,conj,A,B,ORDER,discrete): quit: #End of input .DE .SP1 .LP The output was: .DS #Finds the recurrence for the constant term of (1/y+1+y)^n: conj is , 1 Let R(n) be the integral of the following w.r.t y .DE .SP1 .EQ ~~~~~~~~~~~~~~~~~~~~~~~{(1/y~+~1~+~y) sup n} over y .EN .SP1 .DS The operator annihilating R(n) is the following times conj .DE .SP1 .EQ ~~~~~~~~~~~~~~~~~~~~ -~ 3~ n~ - ~3~ +~ (-~ 2 ~n ~- ~3)~ N ~+~ (n~ +~ 2)~ N sup 2 .EN .SP1 .DS G(n,y) is the integrand times .DE .SP1 .EQ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ {(y - 1) (y + 1) (1 + y + y sup 2 ) } over { y } .EN .SP1 .DS words used=130521, alloc=73728, time=23.02 .DE .SP1 .LP \fB4.\fR .SP1 .LP With the following input the program finds, and proves, the recurrence equation satisfied by the Legendre polynomials, when they are given in terms of their generating function([Ra], p. 157, (1)). .SP1 .DS #The recurrence eq. for Legendre out of the gen. fun. read `WHATYOUWILL`: expli:=0: A:=1/(y^(n+1)*(1-2*x*y+y^2)^(1/2)): ORDER:=2: duis(expli,conj,A,B,ORDER,discrete): quit: .DE .SP1 .LP The output was: .SP1 .DS #The recurrence eq. for Legendre out of the gen. fun. conj is , 1 Let R(n) be the integral of the following w.r.t y .DE .SP1 .EQ ~~~~~~~~~~~~~~~~~~~~~~~~~~ 1 over { y sup {(n~+~1)}~(1~-~2~x~y~+~y sup 2 ) sup 1/2 } .EN .SP1 .DS The operator annihilating R(n) is the following times conj .DE .SP1 .EQ ~~~~~~~~~~~~~~~~~~~~~~ n~ +~ 1~ -~ x~ (3~ +~ 2~ n)~ N~ +~ (n~ + ~2)~ N sup 2 .EN .SP1 .DS G(n,y) is the integrand times .DE .SP1 .EQ ~~~~~~~~~~~~~~~~~~~~~~~ { - ~1~ +~ 2~ x~ y~ -~ y sup 2 } over y .EN .SP1 .DS words used=117996, alloc=65536, time=21.25 .DE .SP1 .SP1 .LP \fB5. Euler's Beta integral\fR .P One of the referees suggested that we present the evaluation of the recurrence equation that leads to Euler's derivation of the beta function. Note that the recurrence, and the proof, hold when $m$ and $n$ are thought of as real arguments, and then by standard uniqueness theorems, factorials can be replaced by Gamma functions. .SP1 .LP The input file was: .SP1 .DS read `WHATYOUWILL` expli:=0: conj:=1: A:=y^n*(1-y)^m: ORDER:=1: duis(expli,conj,A,B,ORDER,discrete): quit: .DE .SP1 .LP The output was: .DS conj is , 1 Let R(n) be the integral of the following w.r.t y n m y (1 - y) The operator annihilating R(n) is the following times conj - n - 1 + (n + m + 2) N G(n,y) is the integrand times y (- 1 + y) words used=74417, alloc=43008, time=12.17 .DE .SP1 .LP \fBAppendix 4: Footnotes\fR .SP1 .LP \fB1.The definition of "holonomic function"\fR .SP1 .LP The precise definition of holonomicity is as follows([Bj],p. viii). Recall that the principal symbol of a linear partial differential operator .SP1 .EQ P(x,D) = sum from { | alpha | <= m } a sub alpha (x) D sup alpha .EN .SP1 .LP is the "highest homogeneous part in the D's": .SP1 .EQ sigma (P)(x, zeta ) = sum from { | alpha | = m } a sub alpha (x) zeta sup alpha .EN .SP1 .LP A function or distribution $F$ is holonomic if the ideal of linear partial differential operators with polynomial coefficients annihilating it .SP1 .EQ I sub F := "{" P(x,D)~;~ P(x,D) F =0 ~"}" .EN .SP1 has the property that the "variety of common zeroes of the principal symbols of the elements of $I sub F$" .SP1 .EQ "{" (x, zeta ) epsilon C sup 2n ~;~ sigma (P) (x, zeta ) =0~ ~,~P(x,D) epsilon I sub F ~"}" .EN has the minimal dimension it could have, namely $n$. (That such a variety cannot have dimension less than $n$ follows from the celebrated Bernstein inequality ([Bj]). .P As an application we will show that a wide class of hyperexponential functions is holonomic. .SP1 .LP \fBLemma A.4.1:\fR Let $S(x,y)$ be a rational function: $S(x,y)=P(x,y)/Q(x,y)$ (P, Q polynomials) such that the varieties {$(x,y) epsilon bold C sup 2 ~;~ P(x,y)=0, Q(x,y)=0$} and {$(x,y) epsilon bold C sup 2 ~;Q sub y (x,y)/G(x,y)=0$ $, Q sub x (x,y) /G(x,y)=0$} ($G:= gcd( Q sub x , Q sub y )$) and {$(x,y) epsilon bold C sup 2 ~;~ G(x,y)=0, Q(x,y)=0$} are zero-dimensional (i.e. consists of a finite set of points) then $F(x,y):=exp(S)$ is holonomic. .SP1 .LP \fBProof:\fR $F(x,y)=exp(S(x,y))$ is easily seen to satisfy the following differential equations .SP1 .EQ(i) ( Q sup 2 D sub x - ( P sub x Q - P Q sub x ) ) F(x,y) == 0 .EN .SP1 .EQ(ii) ( Q sup 2 D sub y - ( P sub y Q - P Q sub y ) ) F(x,y) == 0~~~. .EN .SP1 .LP $( P sub y Q - P Q sub y ) (i)-( P sub x Q - P Q sub x )(ii)$ yields that $F(x,y)$ also satisfies ( since $F(x,y)$ is $C sup inf$ we can "divide" by $Q$): .SP1 .EQ(iii) ( ( P sub y Q - P Q sub y ) D sub x -( P sub x Q - P Q sub x ) D sub y ) F(x,y) == 0~~~. .EN .SP1 .LP Thus the variety of common zeroes of the principal symbols of $I sub F$ is a subset of the variety .SP1 .EQ "{"( x , y, zeta sub 1 , zeta sub 2 ) epsilon bold C sup 4 ~;~ Q sup 2 zeta sub 1 =0 , Q sup 2 zeta 2 =0 , ( ( P sub y Q - P Q sub y ) zeta sub 1 -( P sub x Q - P Q sub x ) zeta sub 2 =0 "}"~~, .EN .SP1 .LP which is easily seen to be two dimensional. It follows that the variety of zeroes has dimension $<= 2$, and by Bernstein's inequality it is actually $2$ and $F$ is holonomic. .P The same method of proof yields that under the condition of the lemma $(P(x,y)/Q(x,y)) sup alpha$ is always holonomic, and, since holonomicity is closed under multiplication([Z1]) also any product of such functions. .P Let us remark that since holonomic functions are \fBD\fR-finite ([Z1], 4.1), the existence of (3) could have also been established by Lipshitz's[Li] method. .SP1 .LP \fB2. The Secret Behind Many Definite Integrals\fR .SP1 .P The secret behind many definite integrals is a lemma that can be found in [He]. It is probably not the earliest reference. It is likely that the trick was known to Euler or at least to Cauchy (Herr quotes Cauchy, among others , in the preface ). .SP1 .LP \fBLemma\fR: .SP1 .LP .EQ int from {- inf } to inf f(x) dx = int from {- inf } to inf f(x- a over x ) dx ~~~~, a>0 .EN .SP1 .LP \fBProof\fR: $t= x - a over x $ if and only if .SP1 .EQ x~= left "{" lpile { { t/2 + sqrt { t sup 2 / 4 + a } ~~~if x > 0 } above { t/2 - sqrt { t sup 2 / 4 + a } ~~~if x < 0 } } .EN .SP1 .LP so, .SP1 .EQ dx~= left "{" lpile { { 1/2 + t/~ sqrt { t sup 2 / 4 + a } ~~~if x > 0 } above { 1/2 - t/ ~sqrt { t sup 2 / 4 + a } ~~~if x < 0 } } .EN .SP1 .LP We thus have .SP1 .EQ int from 0 to inf f(x- a over x ) dx = int from { - inf } to inf f(t) (1/2 + t/~ sqrt ) dt~~~, .EN .SP1 .EQ int from {- inf } to 0 f(x- a over x ) dx = int from { - inf } to inf f(t) (1/2 - t/ ~ sqrt ) dt~~~. .EN .SP1 The proof of the lemma is completed by adding the above two equations. .SP1 .LP \fB3. How Feynman would have evaluated the integral of section 3\fR .SP1 .LP Since the derivation of the integral .SP1 .EQ R(x) := int from { - inf} to inf exp( - {x sup 2} over {y sup 2} - y sup 2 ) dy~~~, .EN .SP1 .LP given in ([Ed], p. 220) is "not as well known as it should be", and is "the argument that Feynman would have used directly", one of the referees suggested that we reproduce this derivation. It runs as follows. .SP1 .EQ R(x)= 2 int from 0 to inf exp( - {x sup 2} over {y sup 2} - y sup 2 ) dy~~~, .EN .SP1 .LP so .SP1 .EQ R'(x)= -4x int from 0 to inf exp( - {x sup 2} over {y sup 2} - y sup 2 ) dy over { y sup 2 } ~~~, .EN .SP1 .EQ =~ -4x int from inf to 0 exp( - t sup 2 - {x sup 2} over {t sup 2} ) ( -xdt over {t sup 2} ) {t sup 2} over { x sup 2 } ~~~, .EN .SP1 .EQ = -2 R(x) ~~~~. .EN .SP1 .LP This integrates to $R(x)=A e sup -2x $ and $x=0$ gives $R(0)= sqrt pi $. .SP1 .LP \fBReferences\fR .SP1 .LP [A-Z] Almkvist, G., and Zeilberger, D., \fIA MAPLE program that finds, and proves, recurrences and differential equations satisfied by hyperexponential definite integrals\fR, SIGSAM Bulletin, to appear. .SP1 .LP [Ba] Bailey, W. N. , "Generalized Hypergeometric Series", Cambridge Math. Tracts \fB32\fR, Cambridge University Press, London, 1935. (Reprinted: Hafner, New York, 1964 .) .SP1 .LP [Be] Bernstein, I. N. , \fIModules over a ring of differential operators, study of the fundamental solutions of equations with constant coefficients\fR, Functional Analysis and Its Applications, Vol \fB5\fR, No. 2(1971), Russian original :1-16, English translation:89-101. .SP1 .LP [Bj] Bjo\*:rk, J. -E., "\fIRings Of Differential Operators\fR", North-Holland, Amsterdam, 1979. .SP1 .LP [Bu] Buchberger ,B. \fIGrobner bases - an algorithmic method in polynomial ideal theory\fR, Chapter 6 in N.K.Bose (ed.): "Multidimensional System Theory", D.Reidel, Publ. Comp., 1985. .SP1 .LP [Co] Coxeter, H. S. M., "Twelve Geometric Essays", Southern Illinois University Press, Carbondale, 1968, preface. .SP1 .LP [DST] Davenport, J. H., Siret, Y., Tournier, E., "\fIComputer Algebra\fR", Academic Press, London, 1988. .SP1 .LP [Ed] Edwards, J., "\fIA treatise on the Integral Calculus\fR", vol. 2, Macmillan, New York, 1921. .SP1 .LP [Eh] Ehlers, F. \fIThe Weyl algebra\fR, in: "Algebraic D-modules", by A. Borel et. al., Academic Press, Boston, 1987. .SP1 .LP [F] Feynman, R. P., "\fISurely You are Joking, Mr. Feynman!\fR", pp. 71-72, paperback, Bantam books 1986. originally published by W. W. Norton, 1985 .SP1 .LP [Fo1] Foata, D., \fIA combinatorial proof of the Mehler formula\fR, J. Comb. Theory, ser. A, \fB24\fR(1978), 250-259. .SP1 .LP [Fo2] Foata, D., \fISome Hermite polynomial identities and their combinatorics\fR, Advances in App. Math. \fB2\fR(1981), 250-258. .SP1 .LP [FG] Foata, D., and Garsia, A. M., \fIA combinatorial approach to the Mehler formulas for Hermite polynomials\fR, in: "\fIRelations between combinatorics and other parts of mathematics\fR" (D. K. Ray-Chaudhury, Ed.), Proceedings of Symposia in Pure Math. No. 34, 163-179, Amer. Math. Soc., Providence, 1979. .SP1 .LP [Gal] Galligo A., \fISome algorithmic questions on ideals of differential operators\fR, in: Lecture Notes in Computer Science \fB204\fR, ( edited by Bob F.Caviness, EUROCAL '85, vol.2), Springer, Berlin. .SP1 .LP [Go] Gosper, R. W., Jr., \fIDecision Procedure of Indefinite Summation\fR, Proc. Natl. Acad. Sci. USA\fR \fB75\fR, 40-42 , 1978. .SP1 .LP [GKP] Graham, R., Patashnik, O., and Knuth, D. E., "\fIConcrete Mathematics\fR", Addison Wesley, Reading, 1989. .SP1 .LP [He] Herr, J. P., "Lehrbuch der hoeren Mathematik", Bd 2, p.356, Wien, 1857. .SP1 .LP [Ho]Hochstadt, H., "\fISpecial Functions of Mathematical Physics\fR", Holt, Rinehart and Winston, New York, 1961. .SP1 .LP [La] Lafon, J. C., \fISummation in finite terms\fR, in: "\fIComputer Algebra Symbolic and Algebraic Computation\fR, edited by B. Buchberger, G. E. Collins, and R. Loos , in cooperation with R. Albrecht, Springer-Verlag, Vienna, second edition, 1983. .SP1 .LP [Li] Lipshitz, L., \fIThe diagonal of a D-finite power series is D-finite\fR, J. Algebra \fB113\fR(1988), 373-378. .SP1 .LP [Me] Melzak, Z. A.,"\fICompanion to Concrete Mathematics\fR", John Wiley, New York, 1973. .SP1 .LP [No] Norman, A. C., \fIIntegration in finite terms\fR, in: "\fIComputer Algebra Symbolic and Algebraic Computation\fR, edited by B. Buchberger, G. E. Collins, and R. Loos , in cooperation with R. Albrecht, Springer-Verlag, Vienna, second edition, 1983. .SP1 .LP [Ra] Rainville, E. D., "\fISpecial Functions\fR", Chelsea, Bronx, 1971. Originally published by Macmillan, 1960. .SP1 .LP [Ri] Risch, R. H. , \fIThe solution of the problem of integrating in finite terms\fR, Bull. AMS \fB76\fR (1970). .SP1 .LP [Ta] Takayama, N., \fIAn approach to the zero recognition problem by the Buchberger algorithm\fR, preprint, Kobe University. .SP1 .LP [Tr] Trager, B. M.,"\fIOn the Integration of Algebraic Functions\fR", Ph. D. thesis, MIT, 1985. .SP1 .LP [W-Z] Wilf, H. S., and Zeilberger, D., \fITowards computerized proofs of identities\fR, Bulletin of the Amer. Math. Soc., to appear. .SP1 .LP [Z1] Zeilberger, D., \fIA holonomic systems approach to special function identities\fR, preprint, Drexel University. .SP1 .LP [Z2] Zeilberger, D., \fIA fast algorithm for proving terminating hypergeometric identities\fR, Discrete Math., to appear. .SP1 .LP [Z3] Zeilberger, D., \fIThe method of creative telescoping\fR, preprint, Drexel University.