%begin macros
\parskip=8pt
\magnification=\magstephalf
\parindent=0pt
\overfullrule=0in
%end macros
\rm

{\bf THE METHOD OF DIFFERENTIATING UNDER THE INTEGRAL SIGN}

\it 
Gert Almkvist
\footnote{*}{\rm
Department of Mathematics, University of Lund, Box 118,
22100 Lund, Sweden.
}
and Doron Zeilberger\rm 
\footnote{**}{\rm
Department of Mathematics and Computer Science, Drexel University,
Philadelphia, PA 19104. Supported in part by NSF grant DMS8800663.
Address in 1999: Dept. of Mathematics, Temple Univ.,
Philadelphia, PA 19122, zeilberg@math.temple.edu
}

{\bf Appeared in J. Symbolic Computation 10 (1990), 571-591}


{\it "I could never resist a definite integral"}
                                 (G. H. Hardy [Co])

\it

"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. ...

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.

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."
\rm
                              (Richard P. Feynman[F])

\medskip

{\bf Introduction}

\rm

The {\it method of differentiating under the integral sign} can
be described as follows. Given a function $F(x,y)$
of $x$ and $y$, one is interested in evaluating

$$
 R(x):= \int_{- \infty}^{\infty} F(x,y) dy \, \, \,. 
\eqno(1)
$$          

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.

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.

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 {\it holonomic functions}. 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.

While the problem of {\it indefinite integration} 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.

{\bf The Theoretical Foundation}

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_1 , ... , x_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).

The algebra of linear differential operators with polynomial coefficients
in $( x_1 , ... , x_n ) $ is called the {\it Weyl algebra}. It
is generated by the "indeterminates" 
$x_1 , ... , x_n$, $D_1 , ... , D_n$, 
that obey the commutation relations

$$
  x_i x_j = x_j x_{i \, \, \,,} D_i D_j = D_j D_{i \, \, \,,} x_i D_j = D_j x_i  \, \, \,(i \neq  j) 
$$

and the "uncertainty principle"

$$
  D_i x_i - x_i D_i =1 \, \, \,, 
$$

where $1$ is the identity operator.

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_x , D_y )$ and $Q(x,y, D_x , D_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_x , D_y )$ and
$B(x,y, D_x , D_y )$ such that

$$
  S(x, D_x , D_y ) = A(x,y, D_x , D_y )P(x,y, D_x , D_y )+ B(x,y, D_x , D_y )Q(x,y, D_x , D_y ) 
$$

is independent of $y$. It follows that $F(x,y)$ is also
annihilated by $S$: $S(x, D_x , D_y ) F(x,y) \equiv 0$.
One can write

$$
 S(x, D_x , D_y ) = \bar{S} (x , D_x ) + D_y S' (x, D_x , D_y ) \, \, \,. 
\eqno(2)
$$         

It is possible to show that it is always possible to have 
$\bar{S} \neq  0$. It follows that

$$
 \bar{S} (x, D_x ) F(x,y) \, = \, - D_y S' ( x, D_x , D_y ) F(x,y)  \, \,. 
\eqno(3)
$$         

In the applications we have in mind, we will have that
$lim_{|y| \rightarrow \infty } x^n y^m D_x^sD_y^tF(x,y)$
$=0$ for all integers $n,m,s$ and $t$.  Therefore integrating with
respect to $y$ from $- \infty$ to $\infty$, we get

$$
  \bar{S} ( x , D_x ) R(x) = 0 \, \, \,. 
$$

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.

{\bf Observation}: In order for things to work we don't need
that $S'$ in (2) be independent of $y$! Our previous insistence
probably made $\bar{S}$ of higher order than necessary.

This leads to the following

{\bf Research Problem}: Develop a fast algorithm, using Grobner
bases in the Weyl algebra ([Ga]), that inputs two operators
$P(x,y, D_x , D_y )$ and $Q(x,y, D_x , D_y )$ and
outputs operators $A,B,C, \bar{S} $ such that

$$
 \bar{S} (x, D_x ) := AP+BQ+ D_y C 
\eqno(4)
$$          

only depends on $x$ and $D_x$.

Basically we have to find an element $\bar{S} (x, D_x )$ that
does not involve $y$ and $D_y$, that belongs to the
"ideal" generated by $P$,$Q$ and $D_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_y$, and so one
has to develop an ambidextrous Grobner bases algorithm in the
Weyl algebra. 

A major step towards this goal was recently achieved by
Takayama[Ta].

A related notion to holonomicity is that of D-finiteness, which
was considered by Lipshitz[Li]. A formal power series
$f( x_1 , ... , x_n )$ is {\it D-finite} if the vector
space generated by all its partial derivatives is finite dimensional
over the field of rational functions in 
$( x_1 , ... , x_n )$.
Every holonomic function is
{\it ipso facto} 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.

{\bf A Detailed Example}

In order to make the ideas clear, we will use ad hoc elimination
to treat the following integral.

$$
  R(x) := \int_{ - \infty}^{\infty} exp( - {{x^2}\over{y^2}} - y^2 ) dy \, \, \,. 
$$

With $F(x,y)=exp( - x^2 / y^2 - y^2 ) $ we have

$$
   \{ { { D_x F = (-2x/ y^2 ) F}      \atop        { D_y F = ( 2 x^2 / y^3 - 2y )F}} 
$$

or

$$
   \{  { { ( y^2 D_x + 2x ) F=0}      \atop        { ( y^3 D_y + 2 y^4 -2 x^2 )F=0} } \, \, \,. 
$$

We want to apply $\int_{- \infty}^{\infty}  dy $ on these equations.
Therefore we want to get rid of all terms containing $y$, except those
of the type $D_y G$. This depends on the fact that

$$
  \int_{ - \infty }^{\infty} D_y G dy = [G]_{ - \infty}^{\infty} = 0 \, \, \,, 
$$

since we assume that everything in sight tends to zero as 
$y \rightarrow \pm \infty$. This means that after integration we can ignore
all terms starting with $D_y$ ( i.e. put $D_y =0$).
Now 

$$
  y^3 D_y = D_y y^3 - 3 y^{2 \, \, \,,} 
$$

hence the second equation becomes

$$
  QF=( D_y y^3 - 3 y^2 + 2 y^4 - 2 x^2 )F =0 \, \,. 
$$

Apply $D_x$ on the left and use the first equation

$$
  y^2 D_x F = - 2xF \, \, \, 
$$

several times. We obtain

$$
  D_x QF=D_y D_x y^3 F  - 3 D_x y^2 F + 2 y^2 D_x y^2 F - 2 D_x x^2 F 
$$

$$
  = D_y D_x y^3 F  - 3 (-2x)F+ 2 y^2 (-2x) F +  - 2 D_x x^2  F \, \, \,. 
$$

$$
  D_x^2Q F= D_y D_x^2y^3 F + 6 D_x x F -4 y^2 ( D_x x ) F - 2 D_x^2x^{2F}  
$$

$$
  = D_y D_x^2y^3 F + 6 D_x x F -4 y^2 (x D_x +1 ) F - 2 D_x^2x^2 F  
$$

$$
  = D_y D_x^2y^3 F + 6 D_x x F + 8 x^2 F - 4 y^2 F -  2 D_x^2x^2 F  \, \,. 
$$

$$
  D_x^3Q F = D_y D_x^3y^3 F + 6 D_x^2x F + 8 D_x x^2 F - 4 D_x y^2 F - 2 D_x^3 x^2 F = 
$$

$$
   = D_y D_x^3y^3 F + 6 D_x^2x F + 8 D_x x^2 F + 8x F - 2 D_x^3 x^2 F  \, \, \,. 
$$

Differentiating

$$
  S(x):= \int_{- \infty}^{\infty} F(x,y) dy 
$$

with respect to $x$, we obtain

$$
  D_x^rS(x) = \int_{- \infty}^{\infty} D_x^rF(x,y) dy \, \,. 
$$

Integrating $D_x^3Q F$, the first term vanishes, and we
get

$$
  (6 D_x^2x + 8 D_x x^2 + 8x - 2 D_x^3x^2 ) S(x) \equiv 0 \, \,. 
$$

Using the rules of "multiplication" in the Weyl algebra, this factors
to:

$$
  ( x D_x + 3)( D_x -2)( D_x +2)S(x) = 0 \, \,. 
$$

This differential equation has the general solution

$$
  S(x)= C_1 e^{2x} + C_2 e^{ -2x} + C_3 y_{3 \, \,,} 
$$

where $y_3 $ satisfies

$$
  y''_3 - 4 y_3 = A/ x^3  \, \, for \, \, A \neq  0 \, \,. 
$$

One can show that $y_3$ cannot be bounded at $0$.
But

$$
  S(0)= \int_{- \infty}^{\infty} e^{ - y^2 } dy = \sqrt \pi  \, \,, 
$$

hence $C_3 =0$. Assume that $x > 0$. Then
$S(x) \rightarrow 0$ as $x \rightarrow \infty $ implies that also $C_1 =0$.
It follows that

$$
  S(x)= \sqrt \pi e^{-2x \, \,}. 
$$

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_{- \infty}^{\infty} exp( - t^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.

{\bf The Class of Hyperexponential Functions}

The discrete analog of the exponential function $ e^{cx}$ is
the geometric sequence $ a_n = c^n $, that has the property
that $a_{n+1} / a_n = $ is constant. An important generalization
of the class of geometric sequences is the 
class of {\it hypergeometric} sequences, that consist of sequences 
{$a_n $} for which $a_{n+1} / a_n $ is a {\it rational}
function of $n$. Hypergeometric sequences are important because
they are the terms of  {\it hypergeometric series}, by which most of
the classical special functions can be expressed.

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 {\it hyperexponential} function to be a function
$f(x)$ for which the logarithmic derivative $f'(x)/f(x)$ is a
rational function of $x$. 

A function of several variables $F( x_1 , ... , x_n )$ is 
hyperexponential if the logarithmic derivatives
$D_{x_i} F / F$, with respect to every variable $x_i$, are
rational functions of $( x_1 , ... , x_n )$.

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 $\bar{S}$ and $S'$ such that (3) holds.
Now set

$$
 G(x,y):= \, - S' ( x, D_x , D_y ) F(x,y) \, \,. 
\eqno(5)
$$          

It is easily proved by induction that for hyperexponential $F$,
$D_x^iD_j^jF /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

$$
 \bar{S} (x, D_x ) F(x,y) \, = \,  D_y G(x,y) \, \,. 
\eqno(6)
$$          

Thus we can forget about $S'$ and concentrate our attention on
finding the winning pair $\bar{S}$ and $G$, now that we know that
they exist.  Suppose you already
knew, or guessed, $\bar{S}$. Then (6) asserts that 
$\bar{S} ( x, D_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 $\bar{S}$ 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 $\bar{S}$ 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.

or: Gosper's Algorithm Translated To The Continuous"

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".

{\bf Input:} 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_y a(y) = {{{d}\over{dy}}} a(y). )$

{\bf Output:}: 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".

We first find polynomials $p(y), q(y), r(y)$ such that

$$
 {{(a(y)/p(y))'}\over {(a(y)/p(y))}} =  {{q(y)}\over {r(y)}} \, \,, 
\eqno(G5)
$$          

where $q(y)$ and $r(y)$ satisfy the condition

$$
 gcd(r(y), q(y)- jr'(y))=1  \, \,for \,every  \, non-negative \,integer \,j. 
\eqno(G6)
$$          

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

$$
  \, 
$$
$$
  p(y)  \,  { \leftarrow } \, p(y) g(y)^{j \,,}  \, r(y)  { \leftarrow}  \, r(y)/g(y) \,,   \,q(y)  { \leftarrow}  \, j(r(y)/g(y))' +(q(y)-jr'(y))/g(y) \,. 
\eqno(G6`)
$$       

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$.

We now write

$$
 S(y)=r(y)f(y)a(y)/p(y) \, \, \,, 
\eqno(G7)
$$          

where $f(y)$ is to be determined. By using $S'(y)=a(y)$,

$$
  f(y)=(S(y)/S'(y))p(y)/r(y) \, \, \,, 
$$

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

$$
  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)) 
$$

Multiplying this through by $p(y)/a(y)$, and using Eq. (G5), we have

$$
 p(y)=(q(y)+r'(y))f(y)+r(y)f'(y) \, \, \,, 
\eqno(G8)
$$          

the differential equation for $f(y)$.

{\bf Theorem:} If $S'(y)/S(y)$ is a rational function of $y$, then
$f(y)$ is a polynomial.

{\bf Proof}: We already know that $f(y)$ is a rational function
when $S'(y)/S(y)$ is, so suppose

$$
 f(y)= {{c(y)}\over {d(y)}} \, \,, 
\eqno(G9)
$$          

where $d(y)$ is a polynomial of positive degree, and

$$
 gcd(c(y),d(y))=1 \, \, \,. 
\eqno(G9a)
$$          

Substituting (G9) in (G8) we get

$$
 qcd+r'cd+r(c'd-cd')= d^2 p \, \,. 
\eqno(G10)
$$          

It is always possible to write

$$
 d(y) = A(y)^k \bar{d} (y) \, \, \,, where \, gcd(A(y), \bar{d} (y)) =1, \, gcd( A(y), A'(y))=1 \, \, \,. 
\eqno(G11)
$$         

Substituting (G11) in (G10) it is readily seen that $A(y)$ divides
$r(y)$, i.e.. for some polynomial $\bar{r} (y)$

$$
 r(y)= A(y) \bar{r} (y) \, \, \, \, \,. 
\eqno(G12)
$$          

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.

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)$.

{\bf Case 1:} 

$$
  deg(r(y)) \leq deg( q(y)+ r'(y))= l \, \, \,. 
$$

Setting $f(y)= c_k y^k + O( y^{k-1} ) $, we get
$p(y)= L c_k y^{k+l} + O( y^{k+l-1} ) $, where $L$ is a
nonzero constant. Since the two sides of the equation must be of
equal degree,

$$
  k= deg(f(y))= deg(p(y))-l \, \, \,. 
$$

{\bf Case 2:}

$$
  deg(q(y)+r'(y)) < deg(r(y))=l \, \, \,. 
$$

Estimating

$$
  f(y)= c_k y^k + c_{k-1} y^{k-1} + O( y^{k-2} ) 
$$

and

$$
  f'(y)= k c_k y^{k-1} + (k-1) c_{k-1} y^{k-2} + O( y^{k-3)} 
$$

in Eq. (G8) leads to 

$$
  p(y)= L(k) c_k y^{k+l-1} + O( y^{k+l-2} ) \, \,, 
$$

where $L(k)$ is linear in $k$ and free of $c_k$. Let
$k_0 $ be the root of $L(k)$. Then, the largest $k$ for
which $c_k$ can be nonzero is

$$
  k=max \{  \, k_0  \,, \, deg(p(y))-l+1  \,\},  \,if \,k_{0 \,} is \,an \,integer \,,or 
$$

$$
  k= deg(p(y))-l+1 \,, \,if  \,not \,. 
$$

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).

{\bf A Fast Algorithm For Definite Hyperexponential Integration}

Consider (1) where $F(x,y)$ is hyperexponential. By 
the theory of holonomic systems, we know {\it for sure} that there
exist $\bar{S} (x, D_x )$ and a hyperexponential $G(x,y)$ such that

$$
 \bar{S} (x, D_x ) F(x,y) =  D_y G(x,y) \, \,. 
\eqno(6)
$$          

We now consider $x$ as an auxiliary parameter and set

$$
  \bar{S} (x, D_x ) = \sum_{0}^{m} a_i (x) D_x^i
$$

where $m$ is the guessed order of $\bar{S}$. Letting the $ a_i$ be
unknown, we evaluate the left of (6), {\it with the undetermined}
$a_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_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 
$\bar{S} (x, D_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 {\it guaranteed} to eventually get a successful $m$.
Furthermore, it is not hard to derive an a priori upper bound for
$m$.

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.

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.

{\bf The method of differencing under the integral sign}

The same idea can be used to find a linear recurrence relation
satisfied by

$$
  a_n := \int_{- \infty}^{\infty} F_n (y) dy \, \,, 
$$

where both $F'_n (y) / F_n (y)$ and 
$F_{n+1} (y) / F_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_n (y)$, and enquire,
by the above algorithm, when it has a hyperexponential primitive
function with respect to $y$.

{\bf Simple example (Courtesy of L. Euler)}: 
If $F_n (y) = y^n e^{-y} $,
we have

$$
  (n+1-N) F_n (y) = D_y (y F_n (y) ) \, \, \,, 
$$

where $N a_n := a_{n+1}$. Integrating w.r.t to $y$ from
$- \infty$ to $\infty$ yields $(n+1-N) a_n =0 $, i.e.
$a_{n+1} = (n+1) a_n $, and so $ a_n = n!$.

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.

{\bf Applications To The Theory Of Special Functions}

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

$$
  t= x - {{a}\over {x}} 
$$

to the integral 
$\int_{- \infty}^{\infty} exp( - t^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.

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

$$
 \sum_{n=0}^{\infty} P_n (x) y^{n} \,=  A(x,y) \, \, \,, 
\eqno(7)
$$          

where $A(x,y)$ turns out to be hyperexponential ( except for the
Jacobi polynomials).  Then (7) can be rewritten as

$$
  P_n (x) = {{1}\over{2 \pi i }}  \int_{|y| = \rho}  
{{A(x,y)} \over {y^{n+1} } } dy \, \, \,. 
$$

Since the integrand is hyperexponential, one can apply our algorithm
to find, and prove, the differential equation satisfied by the sequence
{$P_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.

Another application is to proving Rodrigues-type formulas( e.g. [Ho],
p.9).  Here we use Cauchy's formula

$$
  {{d^n}\over{dx^n}} f(x)= {{n!}\over{(2 \pi i)^n }} \int_{ |y-x| = \rho } 
{{f(y)} \over { (y-x)^{n+1} } }dy \, \, \,. 
$$

By using  the algorithms of sections 6 and 7 we can
get a differential equation and a recurrence equation respectively.

Another application of the algorithm of section 7 is for finding
a pure recurrence relation in $n$ for sequences {$a_n $} given by

$$
  a_n = \, coeff. \, of \, x^b  \,in \, P(x)^{n \, \, \,,} 
$$

where $P(x)$ is a Laurent polynomial. This is useful in 
one-dimensional random walk problems. 

{\bf Acknowledgement}: 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.

\bf 
Appendix 1: A Maple Program Implementing The Algorithm Of
Section 6\rm 

{\bf A.1.1 The program itself}:

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.

{\bf A.1.2: A Typical Input File}:

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.

Place the following input file in a file called INPUT.
You should replace WHATYOUWILL by the full path name of the program
file. 

\#Beginning Of Typical Input File
read `WHATYOUWILL`:
expli:=1:
A:=exp(-x2/y2-y2):
B:=exp(-2*x):
ORDER:=1:
duis(expli,conj,A,B,ORDER,disorcon):
quit:
\#End of Typical Input File

In order to use the program, you type on your terminal

maple $<$ INPUT $>$ OUTPUT

After the execution is completed you can read the output file
by entering  "cat OUTPUT".

The program can find both a differential equation and a recurrence
equation. The default is the former case, where {\it disorcon} 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 {\it expli:=0}. 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 {\it conj:=1}. In every case the
output is the differential operator $gorem(x,D)$ and the hyperexponential
function $G(x,y)$ such that

$$
  gorem(x,D)conj(x,D)A(x,y)=D_y G(x,y) \, \, \,. 
$$

In other words $\bar{S} (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).

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 {\it disorcon:=discrete}.
In the discrete case set expli:=0  always.

\bf 
Appendix 2: Some Sample Inputs and Outputs For Finding Differential
Equations\rm 

1. With the following input the program proves 

$$
  \int_{- \infty}^{\infty} exp( - x^2 / y^2 - y^2 ) dy \,= \, \sqrt \pi e^{-2x \, \, \,,} 
$$

after the appropriate initial conditions are checked.

\#Begin input
read `WHATYOUWILL`:
expli:=1:
conj:=1:
A:=exp(-x2/y2-y2):
B:=exp(-2*x):
ORDER:=1:
duis(expli,conj,A,B,ORDER,disorcon):
quit:
\#End input

The output was:

                      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

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \,exp  ( { - {{ x^2 }\over{ y^2 }}                               - y^2 }  ) 
$$

           The operator annihilating R(x) is the following times conj
                                    - 2 + D
                         G(x,y) is the integrand times
                                     2 1/y

{\bf 2.} With the following input, 
the program finds, and then proves, that

$$
  R(x):= \int_{ - \infty}^{\infty} exp( - x^2 / y^2 + a y^2 ) dy \, \, 
$$

satisfies the differential equation $(4a + D^2 )R(x) = 0$.

\#Begin input
read `WHATYOUWILL`:
expli:=0:
conj:=1:
A:=exp(-x2/y2+a*y2):
ORDER:=2:

duis(expli,conj,A,B,ORDER,disorcon):
quit:
\#End input

The output was:

                             conj is  1
               Let R(x) be the integral of the following w.r.t y

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \,exp  ( { - {{ x^2 }\over{ y^2 }}                               + a y^2 }  ) 
$$

           The operator annihilating R(x) is the following times conj

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \,                                               4 a + D^2 
$$

                         G(x,y) is the integrand times
                                     2 1/y

{\bf 3.}

\#Begin of input
read `WHATYOUWILL`:
expli:=0:
conj:=1:

A:=exp(-x6/y4-y2):
ORDER:=3:
duis(expli,conj,A,B,ORDER,disorcon):
quit:
\#End of input

The output was

                             conj is 1
               Let R(x) be the integral of the following w.r.t y

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \,exp  ( { - {{ x^6 }\over{ y^4 }}                               -  y^2 }  ) 
$$

           The operator annihilating R(x) is the following times conj

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \,                         216 x^5 + 7 D - 12 x D^2  + 4 x^2  D  
$$

                      G(x,y) is the integrand times

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, -54 {{( 4 x^6 - y^4 + 2 y^6 ) x^5 }\over{ y^7 }} 
$$

{\bf 4.} 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:

$$
  \sum_{n=0}^{\infty} P_n (x) t^n  \,= \, (1-2xt+ t^2 )^{-1/2}  \, \, \, \,. 
$$

\#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+y2)(1/2)):
ORDER:=2:
duis(expli,conj,A,B,ORDER,disorcon):
quit:
\#End of input

The output was:

\#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

$$
{ {(1 - 2 x y +y^2 )^{1/2} } \over
 {y^{n+1}} }
$$

           The operator annihilating R(x) is the following times conj

$$
n  \,(n \, + \, 1) \, - \, 2 \, x \, D  \,- \,  
(x \, - \, 1) (x \, +  \,1) \, D^2 
$$

                         G(x,y) is the integrand times

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \,                                            - \, {{  (- \, n  \,-  \,1 + \, 2 \, n \, x  \,y +   \,x \, y  \,- \, n  \,y^2 )  \,y }\over{            (- \, 1  \,+ \, 2 \, x \, y \, - \, y^2 ) \, n        }}  
$$

words used=156843, alloc=88064, time=29.37

{\bf 5.}

With the following input, the program proves the Mehler
formula for Hermite polynomials (e.g. [Ho], p. 24, ex. 4).

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 {\it verification}, Foata's proof is a 
{\it derivation}, 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, {\it even in principle}, of proving multi-variate
identities with an {\it arbitrary} number of variables.

\#Begin input
\#The Mehler formula for Hermite polynomials
read `WHATYOUWILL`:
expli:=0:
A:=exp((2*x*y*z-y2*(x2+z2))/(2*(1-y2)))/((1-y2)(1/2)*y(n+1)):
ORDER:=2:

duis(expli,conj,A,B,ORDER,disorcon):
quit:
\#End of input

The output was:

\#The Mehler formula for Hermite polynomials
                             conj is 1
               Let R(x) be the integral of the following w.r.t y

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, {exp( {{2 \,x \,y \,z \,-y^2  \, ( x^2  \, + z^2 )}\over{2 \,- \,2 \,y^2 }} ) } 
\over { (1 \,- \, y^2 )^{1/2}  \,y^{(n \,+1) } } 
$$

           The operator annihilating R(x) is the following times conj

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \,                                                                                -  \,n  \,+ \, x \, D \, - \, D^2 
$$

                         G(x,y) is the integrand times
                                       y
words used=138796, alloc=79872, time=26.87

{\bf 6.} With the following input the program proves the Doetsch identity
(see [Fo2]).

\#To Prove the Doetsch formula for Hermite polynomials(see Foata,Adv.
\#Appl. Math. 2(1981),250-259.)
read `WHATYOUWILL`:
expli:=0:

A:=y(-n-1)*(1+2*y2)(-3/2)*(1+y*x+2*y2)*exp(y2*x2/(1+2*y2)):
ORDER:=2:
duis(expli,conj,A,B,ORDER,disorcon):
quit:

The output was:

                                  conj is , 1
               Let R(x) be the integral of the following w.r.t y

$$
{
{ y^{(- \,n \,-1)} (1 \,+ \,y \,x \,+2 \,y^2 ) 
exp( {{ y^2 x^2 }\over{ 1 \,+ \, 2 y  \,^2 }} )  
} 
\over 
{(1 \,+ \,2 \,y^2 )^{3/2} } 
}
$$

           The operator annihilating R(x) is the following times conj

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \,-n \,+x \,D \,- \,D^2 
$$

                         G(x,y) is the integrand times
                                      y
words used=194380, alloc=112640, time=40.10

{\bf 7.}

With the following input, the program proves the Rodrigues formula
for Jacobi polynomials ([Ra], p. 257)

\#Begin input
\#Proves Rodrigues' Formula for Jacobi polynomials by finding the
\#diff. eq
read `WHATYOUWILL`:
expli:=0:
A:=(y2-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

The output was:

\#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

$$
   \, \, \, \, \, \, \, \, {{ ( y^2  \,- \,1)^n  \,(1 \,- \,y)^a  \,(1+ \,y)^b  \,(1 \,-x)^{(-  \,a)} \,(1 \,+ \,x)^{(- \,b)} } \over { (y- \,x)^{(n \,+ \,1)} }} 
$$

           The operator annihilating R(x) is the following times conj

$$
      - (n + a + 1 + b) n + (a + a x - b + x b + 2 x) D +  (- 1 + x) (1 + x) D^2 
$$

                         G(x,y) is the integrand times

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \,                      {{      (n \, +  \,1) (- \, 1 \, + \, y)  \,(1 \, + \, y)  }\over{(-  \,y \, + \, x) \,  }} 
$$

words used=397108, alloc=155648, time=96.80

{\bf 8.}

With the following input, the program proves Euler's
integral representation for the hypergeometric function
${}_2 F_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.

\#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

The output was:

\#Proves Euler's representation for 2F1
                                  conj is , 1
               Let R(x) be the integral of the following w.r.t y

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, (1 \,- \,x \,y)^{(- \,a)} \,y^{(b \,- \,1)} \,(1- \,y \,)^{(c \,- \,b \,- \,1)} 
$$

           The operator annihilating R(x) is the following times conj

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \,                    b \, a \, + \, (x \, + \, a \, x \, +  \,x  \,b \, - \, c) \, D +   \,x \, (x \, - \, 1)  \,D^2 
$$

                         G(x,y) is the integrand times

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \,                                 {{ a y \, (-  \,1 \, +  \,y) }\over{                                 { -  \,1 \, + \, x \, y} }} 
$$

words used=132760, alloc=77824, time=25.07

{\bf 9.}

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.

\#Begin input
read `WHATYOUWILL`:
expli:=0:
A:=x2/((x3+y3)*(1+y3)):
ORDER:=2:

duis(expli,conj,A,B,ORDER,disorcon):
quit:
\#End input

The program outputted the following differential operator. We omit
$G(x,y)$.

$$
  3 \,x^2  \,+ \,D \,+5 \,x^3  \,D \,+ \,x^4  \,D^{2} \,-x \,{D}^2 
$$

\bf 
Appendix 3: Sample Inputs And Outputs For Finding Recurrence 
Equations\rm 

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
{\it discrete} , and $A$ should involve $n$ and $y$. $N$ is the shift
operator in the $n$ direction: $Nf(n):=f(n+1)$.

{\bf 1.}

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.

\#Beginning of input
read `WHATYOUWILL`:
expli:=0:
A:=(y2-1)n/(y-x)(n+1):
ORDER:=2:
duis(expli,conj,A,B,ORDER,discrete):
quit:
\#End of input

The output was:

                                  conj is , 1
               Let R(n) be the integral of the following w.r.t y

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, {{( y^2  \,- \,1)^n}\over{(y \,- \,x)^{(n \,+ \,1)}}} 
$$

           The operator annihilating R(n) is the following times conj

$$
                                              -  \,4 \, n \, -  \,4 \, + \, 2 \, x  \,(3 \, + \, 2 \, n) \, N  \, + \, (- \, n \, - \, 2) \, N^2 
$$

                         G(n,y) is the integrand times

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, - {{(y \,- \,1) \,(y \,+ \,1) \,(- \, y^2  \,+ \,2 \,y \,x \,- \,1)}\over{- \,y \,+ \,x}} 
$$

words used=161099, alloc=94208, time=31.43

{\bf 2.}

With the following input, the program finds the recurrence equation
satisfied by the Jacobi polynomials, when defined in terms of their
Rodrigues's formula. 

\#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:=(y2-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 

The output is omitted( do it on your computer!).

{\bf 3.}

With the following input the program finds, and proves the
recurrence satisfied by the sequence
$R(n):=$constant term of $(y+1+1/y)^n$.

\#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

The output was:

\#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

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \,{{(1/y \,+ \,1 \,+ \,y)^n}\over{y}} 
$$

           The operator annihilating R(n) is the following times conj

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \,                       - \, 3 \, n \, -  \,3 \, + \, (- \, 2  \,n  \,-  \,3) \, N   \,+ \, (n \, + \, 2) \, N^2 
$$

                         G(n,y) is the integrand times

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \,                                                                                {{(y - 1) (y + 1) (1 + y + y^2 ) }\over{  y }} 
$$

words used=130521, alloc=73728, time=23.02

{\bf 4.}

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)).

\#The recurrence eq. for Legendre out of the gen. fun.
read `WHATYOUWILL`:
expli:=0:
A:=1/(y(n+1)*(1-2*x*y+y2)(1/2)):
ORDER:=2:
duis(expli,conj,A,B,ORDER,discrete):
quit:

The output was:

\#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

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, {{1} \over { y^{(n \,+ \,1)} \,(1 \,- \,2 \,x \,y \,+ \,y^2 )^{1/2} } }
$$

           The operator annihilating R(n) is the following times conj

$$
    \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \,                                                  n \, + \, 1 \, - \, x \, (3 \, + \, 2 \, n) \, N \, + \, (n \, +  \,2) \, N^2 
$$

                         G(n,y) is the integrand times

$$
   \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \,    {{    -  \,1 \, + \, 2 \, x \, y \, - \, y^2 }\over{y}} 
$$

words used=117996, alloc=65536, time=21.25

{\bf 5. Euler's Beta integral}

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.

The input file was:

read `WHATYOUWILL`
expli:=0:
conj:=1:
A:=yn*(1-y)m:
ORDER:=1:

duis(expli,conj,A,B,ORDER,discrete):
quit:

The output was:

                                  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

{\bf Appendix 4: Footnotes}

{\bf 1.The definition of "holonomic function"}

The precise definition of holonomicity
is as follows([Bj],p. viii).  Recall that the principal
symbol of a linear partial differential operator

$$
  P(x,D) = \sum_{ | \alpha | \leq m } a_{\alpha} (x) D^{\alpha} 
$$

is the "highest homogeneous part in the D's":

$$
  \sigma (P)(x, \zeta ) = \sum_{ | \alpha | = m }  a_{\alpha} (x) 
\zeta^{\alpha} 
$$

A function or
distribution $F$ is holonomic if the ideal of linear partial
differential operators with polynomial coefficients annihilating
it

$$
  I_F := \{ P(x,D) \,; \, P(x,D) F =0  \,\} 
$$

has the property that the "variety of common zeroes of the principal
symbols of the elements of $I_F$"

$$
  \{ (x, \zeta )  \epsilon  C^{2n}  \,; \, \sigma (P) (x, \zeta ) =0 \,  \,, \,P(x,D)  \epsilon  I_F  \,\} 
$$
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]).

As an application we will show that a wide class of hyperexponential
functions is holonomic.  

{\bf Lemma A.4.1:} 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   C^2  \,; \, P(x,y)=0, Q(x,y)=0$} and
{$(x,y)  \epsilon   C^2  \,;Q_y (x,y)/G(x,y)=0$
$, Q_x (x,y) /G(x,y)=0$} ($G:= gcd( Q_x , Q_y )$) 
and
{$(x,y)  \epsilon   C^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.

{\bf Proof:} $F(x,y)=exp(S(x,y))$ is easily seen to satisfy the following
differential equations

$$
 (  Q^2 D_x - ( P_x Q - P Q_x ) ) F(x,y) \equiv 0 
\eqno(i)
$$          

$$
 (  Q^2 D_y - ( P_y Q - P Q_y ) ) F(x,y) \equiv 0 \, \, \,. 
\eqno(ii)
$$          

$( P_y Q - P Q_y ) (i)-( P_x Q - P Q_x )(ii)$ yields
that $F(x,y)$ also satisfies ( since $F(x,y)$ is $C^{\infty}$ we can
"divide" by $Q$):

$$
 ( ( P_y Q - P Q_y ) D_x -( P_x Q - P Q_x ) D_y ) F(x,y) \equiv 0 \, \, \,. 
\eqno(iii)
$$         

Thus the variety of common zeroes of the principal symbols of
$I_F$ is a subset of the variety

$$
  \{( x , y, \zeta_1 , \zeta_2 )  \epsilon  {\bf C}^4  \,; \, Q^2 
\zeta_1 =0 , Q^2 \zeta 2 =0 ,  ( ( P_y Q - P Q_y ) 
\zeta_1 -( P_x Q - P Q_x )  \zeta_2 =0 \} \, \,, 
$$

which is easily seen to be two dimensional.  It follows that
the variety of zeroes has dimension $\leq 2$, and by Bernstein's
inequality it is actually $2$ and $F$ is holonomic.

The same method of proof yields that under the condition of the lemma
$(P(x,y)/Q(x,y))^{\alpha}$ is always holonomic, and, since
holonomicity is closed under multiplication([Z1]) also any
product of such functions.

Let us remark that since
holonomic functions are 
{\bf D}-finite ([Z1], 4.1), the
existence of (3) could have also been established by 
Lipshitz's[Li] method.

{\bf 2. The Secret Behind Many Definite Integrals}

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 ).

{\bf Lemma}:

$$
  \int_{- \infty }^{\infty} f(x) dx = \int_{- \infty }^{\infty} f(x- {{a}\over {x}} ) dx  \, \, \, \,, a>0 
$$

{\bf Proof}: $t= x - {{{a}\over{x}}}$  if and only if

$$
  x \,=  \{  { { t/2 + \sqrt { t^2 / 4 + a }  \, \, \,if x > 0 } \atop { t/2 - \sqrt { t^2 / 4 + a }  \, \, \,if x < 0 } } 
$$

so,

$$
   dx \,=  \{  { { 1/2 + t/ \, \sqrt { t^2 / 4 + a }  \, \, \,if x > 0 } \atop { 1/2 - t/  \,\sqrt { t^2 / 4 + a }  \, \, \,if x < 0 } } 
$$

We thus have

$$
  \int_{0}^{\infty} f(x- {{a}\over {x}} ) dx = \int_{ - \infty }^{\infty} f(t) (1/2 + t/ \, \sqrt ) dt \, \, \,, 
$$

$$
  \int_{- \infty }^{0} f(x- {{a}\over {x}} ) dx = \int_{ - \infty }^{\infty} f(t) (1/2 - t/  \, \sqrt ) dt \, \, \,. 
$$

The proof of the lemma is completed by adding the above two
equations.

{\bf 3. How Feynman would have evaluated the integral of section 3}

Since the derivation of the integral

$$
  R(x) := \int_{ - \infty}^{\infty} exp( - {{x^2}\over{y^2}} - y^2 ) dy \, \, \,, 
$$

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.

$$
  R(x)= 2 \int_{0}^{\infty} exp( - {{x^2}\over{y^2}} - y^2 ) dy \, \, \,, 
$$

so

$$
  R'(x)= -4x \int_{0}^{\infty} exp( - {{x^2}\over{y^2}} - y^2 ) {{dy}\over{ y^2 }}  \, \, \,, 
$$

$$
  = \, -4x \int_{\infty}^{0} exp(  - t^2 - {{x^2}\over{t^2}}  ) ( {{-xdt}\over{t^2}} ) {{t^2}\over{ x^2 }}  \, \, \,, 
$$

$$
  = -2 R(x)  \, \, \, \,. 
$$

This integrates to $R(x)=A e^{-2x} $ and $x=0$ gives 
$R(0)= \sqrt \pi $.

{\bf References}

[A-Z] Almkvist, G., and Zeilberger, D., 
\it 
A MAPLE program that finds, and proves, recurrences and differential
equations satisfied by hyperexponential  definite 
integrals\rm 
, SIGSAM Bulletin, to appear.

[Ba] Bailey, W. N. , "Generalized Hypergeometric Series",
Cambridge Math. Tracts {\bf 32}, Cambridge University Press, London,
1935. (Reprinted: Hafner, New York, 1964 .)

[Be] Bernstein, I. N. , \it 
Modules over a ring of
differential operators, study of the fundamental solutions of
equations with constant coefficients\rm 
, Functional Analysis and Its
Applications, Vol {\bf 5}, No. 2(1971), Russian original :1-16,
English translation:89-101.

[Bj] Bjo\*:rk, J. -E., "{\it Rings Of Differential Operators}",
North-Holland, Amsterdam, 1979.

[Bu] Buchberger ,B. \it 
Grobner bases - an algorithmic method
in polynomial ideal theory\rm 
, Chapter 6 in N.K.Bose (ed.):
"Multidimensional System Theory", D.Reidel, Publ. Comp., 1985.

[Co] Coxeter, H. S. M., "Twelve Geometric Essays",
Southern Illinois University Press, Carbondale, 1968, preface.

[DST] Davenport, J. H., Siret, Y., Tournier, E., 
"{\it Computer Algebra}", Academic Press, London, 1988.

[Ed] Edwards, J., "{\it A treatise on the Integral Calculus}",
vol. 2, Macmillan, New York, 1921.

[Eh] Ehlers, F. {\it The Weyl algebra}, in: "Algebraic D-modules",
by A. Borel et. al., Academic Press, Boston, 1987.

[F] Feynman, R. P., "{\it Surely You are Joking, Mr. Feynman!}",
pp. 71-72, paperback, Bantam books 1986. originally published by
W. W. Norton, 1985

[Fo1] Foata, D., {\it A combinatorial proof of the Mehler formula},
J. Comb. Theory, ser. A, {\bf 24}(1978), 250-259.

[Fo2] Foata, D., \it 
Some Hermite polynomial identities and their
combinatorics\rm 
, Advances in App. Math. {\bf 2}(1981), 250-258.

[FG] Foata, D., and Garsia, A. M., \it 
A combinatorial approach to
the Mehler formulas for Hermite polynomials\rm 
, in: "\it 
Relations
between combinatorics and other parts of mathematics\rm 
" (D. K.
Ray-Chaudhury, Ed.), Proceedings of Symposia in Pure Math. No. 34,
163-179, Amer. Math. Soc., Providence, 1979.

[Gal] Galligo A., \it 
Some algorithmic questions on ideals of
differential operators\rm 
, in: Lecture Notes in Computer Science
{\bf 204}, ( edited by Bob F.Caviness, EUROCAL '85, vol.2),
Springer, Berlin.

[Go] Gosper, R. W., Jr., {\it Decision Procedure of Indefinite Summation},
Proc. Natl. Acad. Sci. USA\rm 
 {\bf 75}, 40-42 , 1978.

[GKP] Graham, R., Patashnik, O., and Knuth, D. E., "\it 
Concrete
Mathematics\rm 
", Addison Wesley, Reading, 1989.

[He] Herr, J. P., "Lehrbuch der hoeren Mathematik", Bd 2, p.356,
Wien, 1857.

[Ho]Hochstadt, H., "{\it Special Functions of Mathematical Physics}",
Holt, Rinehart and Winston, New York, 1961.

[La] Lafon, J. C., {\it Summation in finite terms}, in:
"{\it Computer Algebra Symbolic and Algebraic Computation}, edited by
B. Buchberger, G. E. Collins, and R. Loos , in cooperation with
R. Albrecht, Springer-Verlag, Vienna, second edition, 1983.

[Li] Lipshitz, L., \it 
The diagonal of a D-finite power series is
D-finite\rm 
, J. Algebra {\bf 113}(1988), 373-378.

[Me] Melzak, Z. A.,"{\it Companion to Concrete Mathematics}", John
Wiley, New York, 1973.

[No] Norman, A. C., {\it Integration in finite terms}, in:
"{\it Computer Algebra Symbolic and Algebraic Computation}, edited by
B. Buchberger, G. E. Collins, and R. Loos , in cooperation with
R. Albrecht, Springer-Verlag, Vienna, second edition, 1983.

[Ra] Rainville, E. D., "{\it Special Functions}", Chelsea, Bronx,
1971. Originally published by Macmillan, 1960.

[Ri] Risch, R. H. , \it 
The solution of the problem of integrating in
finite terms\rm 
, Bull. AMS {\bf 76} (1970).

[Ta] Takayama, N., \it 
An approach to the zero recognition problem
by the Buchberger algorithm\rm 
, preprint, Kobe University.

[Tr] Trager, B. M.,"{\it On the Integration of Algebraic Functions}",
Ph. D. thesis, MIT, 1985.

[W-Z] Wilf, H. S., and Zeilberger, D., \it 
Towards computerized
proofs of identities\rm 
, Bulletin of the Amer. Math. Soc., to
appear.

[Z1] Zeilberger, D., \it 
A holonomic systems approach to special
function identities\rm 
, preprint, Drexel University.

[Z2] Zeilberger, D., \it 
A fast algorithm for proving terminating
hypergeometric identities\rm 
, Discrete Math., to appear.

[Z3] Zeilberger, D., {\it The method of creative telescoping}, 
preprint, Drexel University.
\end
