A HOLONOMIC SYSTEMS APPROACH TO SPECIAL FUNCTIONS IDENTITIES

Doron Zeilberger

Appeared in J. of Computational and Applied Math. 32 (1990), 321-368.

Certain functions appear so often that it is convenient to give them names. These functions are collectively called special functions. There are many examples and no single way of looking at them can illuminate all examples or even all the important properties of a single example of a special function.– (Richard Askey[Ask1])

ABSTRACT

We observe that many special functions are solutions of so-called holonomic systems. I. N. Bernstein’s deep theory of holonomic systems is then invoked to show that any identity involving sums and integrals of products of these special functions can be verified in a finite number of steps. This is partially substantiated by an algorithm that proves terminating hypergeometric series identities, and that is given both in ENGLISH and in MAPLE.

MR Classifications: Primary: 05A10, 05A19, 33A30, 33A99; Secondary: 05A30 , 33A65, 33-04, 35N, 39A10.

Supported in part by NSF grant DMS 8800663

1. Introduction

We are nowadays witnessing a spectacular comeback of the theory of special functions. This is expressed both by the various approaches ([Foa]),[Mill],[Ask1],[Ehr2]) that are successfully used to explain and give insight to large families of previously unrelated results, and by the dramatic applications of special function identities to the solution of longstanding open problems in pure mathematics ([deB], [Ask-Gas]; [vdP]).

This paper initiates yet another approach to special functions that is based on I.N.Bernstein’s theory of holonomic systems. Unlike the other approaches, the present approach does not give insight on any one particular identity. Instead it gives a kind of universal insight on many identities at once. It implies that a large class of special function identities, that includes all terminating hypergeometric (alias binomial coefficients) identities, is verifiable in a finite number of steps.

One of the most salient features of the classical sequences of orthogonal polynomials p_{n}(x) is that they satisfy both a second order linear differential equation and a second order linear recurrence equation, both with coefficients that are polynomials in n and x. For example, the Legendre polynomials P_{n}(x) satisfy ([Rai], ch.10):

(1 − x^{2})P′′_{n}(x) − 2xP′_{n}(x) + n(n + 1)P_{n}(x) = 0,  (1.1a)

(n + 2)P_{n+2}(x) − (2n + 3)xP_{n+1}(x) + (n + 1)P_{n}(x) = 0.  (1.1b)

Viewing the sequence P_{n}(x) as one function on N × R: F (n, x) := P_{n}(x) and introducing the differentiation and shift operators:

Df(n, x) := \frac{d}{dx}f(n, x) , Ef(n, x) := f(n + 1, x),

(1.1) can be rewritten as

[(1 − x^{2})D^{2} − 2xD + n(n + 1)]F (n, x) = 0,  (1.2a)

[(n + 2)E^{2} − (2n + 3)xE + (n + 1)]F (n, x) = 0,  (1.2b)

A system like (1.2) is called ” a maximally overdetermined system of linear differential-recurrence equations with polynomial coefficients on N × R”, or a holonomic system on N × R for short. Since (1.2a) and (1.2b) are independent,(in a certain technical sense to be made precise later), they uniquelydetermine F (n, x) given a finite number of initial conditions. (In this example these are: F (0, 0) = 1, F′(0, 0) = 0, F (1, 0) = 0, F′(1, 0) = 1.) A function like F (n, x) that is a solution of a holonomic system is called a holonomic function on N × R.

If we allow linear equations of any order, and generalize from one discrete and one continuous variable to several variables of each kind, we are lead to consider the class of holonomic functions f(x_{1}, ..., x_{n_{1}}, m_{1}, ..., m_{n_{2}}) on R^{n_{1}} × N^{n_{2}}. These are functions that satisfy ”as many (homogeneous ) linear (differential-recurrence) equations as possible (with polynomial coefficients)”. They are uniquely determined as solutions of a system of such equations subject to a finite number of initial conditions. The notion of ”satisfying as many equations as possible” is made precise in I.N.Bernstein’s theory of holonomic systems described in section 2.

The following two crucial facts make holonomic functions an ideal framework for special function identities:

(i) Every holonomic function can be described by a finite amount of information, and has a ”canonical holonomic representation”.

(ii) The product, addition, sum, and integral of holonomic functions is again holonomic. Furthermore, if one has canonical holonomic representations of f and g, then it is possible to find canonical holonomic representations for f + g, ,f − g, ,fg, , ∫ fdx_{1}, and ∑_{m_{1}} f, where x_{1} and m_{1} are any continuous and discrete variable respectively.

It follows that an expression like

P_{n_{1}}(x)P_{n_{2}}(x)P_{n_{3}}(x)

is holonomic in (n_{1}, n_{2}, n_{3}, x), and that its (definite) integral is holonomic in (n_{1}, n_{2}, n_{3}) . Since all the special functions that appear in the famous ”tableau d’Askey” ([Ask-Wil], [Lab]) are given by hypergeometric summation, whose summands are obviously holonomic in all their variables and parameters , it follows that they are all holonomic, both in their variable and their parameters, as are all expressions obtained from them by adding, multiplying, summing, and integrating. For example, the Jacobi polynomials P_{n}^{(α,β)} are not only holonomic in (n, x) for fixed (α, β), but are in fact holonomic in (n, x, α, β). It follows that an expression like

∫ e^{−α} ∑_{n} L_{n}^{(α)}(x)P_{n}^{(α,β)}(x)dα,

where L_{n}^{(α)}(x) are the generalized Laguerre polynomials, is holonomic in x and β.

The above two properties of holonomic functions enable us, at least in principle, to prove, or refute, any identity involving sums and integrals of products of holonomic functions. All one has to do is to bring everything to the left side and leave only 0 on the right side. Then a canonical holonomic representation of the expression on the left is evaluated step by step , and then it is determined whether such a representation is equivalent to 0. ( see section 4.) Note that canonical holonomic representations are not unique, but it is always possible to know when such a representation is equivalent to 0, and thus it is possible to know when two different canonical representations represent the same function. This is analogous to the fact that among the rational numbers a/b, with a and b integers, 0 has an infinite number of representations: 0/b, and we know that whenever the numerator is 0, the fraction represents 0. Because of this we can tell when two fractions a/b and c/d represent the same rational number: ad−bc = 0, since ad−bc is the numerator of a/b − c/d. Since the sum and product of rational numbers are again rational, and we know how to find ”canonical representations” for sums and products, it follows that any identity involving sums and products of (specific) rational numbers is provable. The same holds for polynomials and algebraic numbers.

Let us see how to prove an identity like

∑_{n=0}^{∞} P_{n}(x)t^{n} − (1 − 2xt + t^{2})^{−1/2} = 0.  (∗)

P_{n}(x) is holonomic in n and x, t^{n} is holonomic in t and n, thus the summand P_{n}(x)t^{n} is holonomic in (x, t, n). Summation with respect to n gets rid of the dependence on n, and the sum is holonomic in the surviving variables (x, t). Now (1−2xt+t^{2})^{−1/2} is clearly holonomic in x, t, as you can easily see by logarithmic differentiation with respect to x and t. Thus the left side of (*) is holonomic in (x, t), and it is possible to find a canonical holonomic representation for it, and determine whether it is indeed equivalent to 0.

An important special case of special functions identities is that of terminating hypergeometric summation. Thanks to the preaching of Richard Askey and his disciples ([Roy]), it is nowadays well known that the theory of combinatorial sums, traditionally pursued by combinatorialists, is just a special case of hypergeometric summation, and that many apparently distinct combinatorial identities are really equivalent qua hypergeometric sums.

I hope that more professional programmers and algorithmic designers will soon expand the rudimentary ideas in this paper and develop a symbolic software package to prove general special function identities. As a modest first step, I present, in sections 5 and 6, an explicit algorithm to verify any identity of the form

∑_{k} F (n, k) = a(n),

where F (n, k) is a quotient of products of factorial expressions of the form (an+bk+c)!, where a and b are integers, and c is a complex number or parameter. Since both (an + bk + c)! and its reciprocal are obviously holonomic in n and k, it follows that F (n, k) is holonomic, and thus that its sum is holonomic in n. But being holonomic in n means that it satisfies a certain (homogeneous)linear recurrence equation with polynomial coefficients. The algorithm of section 6 finds such a recurrence and then checks whether the right side a(n) also satisfies the same recurrence, and whether the initial conditions match.

A MAPLE program implementing this algorithm is given in the appendix. This program can be used by all those who have MAPLE, and who are willing to type it in, or who have e-mail, so that I can send it to them. In order to use it all you have to do is read the simple instructions in the appendix. Readers who are only interested in binomial coefficients identities should go directly to section 5.

I should warn the reader that my program requires a lot of memory, and it still remains to be seen whether it can prove all known identities with the memory available on a big computer. However, even on my AT&T 3B1 PC, the program was able to prove many non-trivial identities, some of which are given in section A.2 of the appendix.

The most time and space consuming part of the algorithms for proving special function identities is elimination. This elimination takes place in the non-commutative Weyl algebra of linear differentialrecurrence operators with polynomial coefficients. The algorithm that I describe is an adaptation of Sylvester’s classical dialytic elimination([VDW]) to the non-commutative context of the Weyl algebra. Sylvester’s method, while elegant and simple theoretically, is known to be computationally inefficient, especially when more that one variable is being eliminated. We are nowadays witnessing an ongoing revolution in computational commutative algebra that is centered around Buchberger’s ([Buc]) concept of Grobner bases. It turns out that, using Grobner bases, it is possible to perform elimination much faster.

Galligo[Gal] showed how to adapt the method of Grobner bases to the non-commutative Weyl algebra. This should result in a considerable improvement in the verification algorithms.

The connection between special functions and holonomic systems has been recently used by Takayama ([Tak1],[Tak2]), who also realized the importance of Grobner bases. Gelfand and his Moscow school ([Gel], [Gel-Gel],[Gel-Gra],[Gel-Zel]) are currently developing a general theory of hypergeometric functions in several variables, and they showed that they are holonomic functions.

Apology: In [Zei1] I attempted to develop a theory of special function identities by generalizing Sister Celine Fasenmayer’s ([Fas]) method. While her method is absolutely correct, my attempt to base a general theory on it was wrong. In addition, the definitions of ”multi-P-finite”, ”multi-P-recursive”, and ”special” given there are inaccurate. I wish to thank L.Lipshitz, R.Stanley, and D.E.Knuth for pointing out the shortcomings of [Zei1]. I am hereby retracting [Zei1] from my list of publications and beg that it be only considered as a historical document.

Limitations Of The Present Approach:

Besides the practical limitations of the verification algorithms, there is also a theoretical barrier. The present approach only shows that an identity with a fixed number of variables is verifiable in a finite number of steps. For example, the binomial theorem

(x + y)^{n} = ∑_{k=0}^{n}({n \atop k})x^{k}y^{n−k}

is trivially finitely verifiable. On the other hand, the multinomial theorem

(x_{1} + ... + x_{m})^{n} = ∑_{k_{1}+...+k_{m}=n} \frac{n!}{k_{1}!...k_{m}!}x_{1}^{k_{1}}...x_{m}^{k_{m}} ,

that has an indefinite number of variables, is beyond the scope of the present theory, and is not finitely verifiable, even not in principle. Of course for every specific m, say m = 1000, the identity is finitely verifiable, but not for an arbitrary number of variables m.

Other examples are Macdonald’s root-system conjectures ([Mac]). For every fixed root-system, the conjecture is verifiable in a finite number of steps, for example for all the exceptional root systems. But for the infinite families, not even zillion years will suffice.

Milne’s ([Miln]) deep theory of hypergeometric SU(N) identities is also not covered by the verification algorithms, since it involves an indefinite number of variables.

The theory and algorithms of this paper are easily extendible to the q-analogs of hypergeometric series and special functions. On the other hand, it is not extendible to bi-basic and multi-basic identities that are presently vigorously studied by Gasper and Rahman (see [Gas],[Gas-Rah1],[Gas-Rah2]).

Yet another example of the limitation of the present approach is furnished by the multi-variate generalizations of the finite form of the Rogers-Ramanujan identities found by Andrews and Bressoud (see [And1], p.30). The finite form of the Rogers-Ramanujan identity is provable by my algorithm, but its multi-variate generalization is not. Going in the other direction, toward the less general, the present method is incapable of proving the Rogers-Ramanujan identities directly, and we need a Schur or an Andrews to come up with a conjectured finite form that tends to the identities as n → ∞. It would be interesting if one could develop heuristics, in the spirit of chapter 10 of [And1], of conjecturing a finite form of a given q-series identity, that could then be proved by the present method.

Acknowledgement: I have already mentioned my indebtedness to Leonard Lipshitz, Richard Stanley, and Donald Knuth for pointing out the inadequacy of [Zei1]. I also wish to thank Richard Askey for his penetrating critique [Ask2] of the philosophical side of [Zei1] and for many helpful discussions.

It was Bernd Sturmfels who told me about Grobner bases and thus gave me hope that my clumsy algorithms may be improved considerably.

Special thanks are due to Joseph Bernstein (I.N.Bernstein) who spent two hours of his precious time explaining his beautiful theory to me, in a language that I understood and that enabled me to transmit the holonomic fire to my fellow mortals in the provinces of special functions and combinatorics.

I would also like to thank the four developers of MAPLE: Bruce Char, Keith Geddes, Gaston Gonnet, and Stephen Watt, for developing such a beautiful system, and Kurt Horton, the local UNIX guru, for informing me that MAPLE is available on the AT&T 3B1.

Many thanks are due to Gert Almkvist for his interest, encouragement, and many valuable suggestions. Readers who can read Swedish will enjoy Almkvist’s[Alm] delightful account of the present ideas.

Finally, I also benefited from numerous conversations with Leon Ehrenpreis, Dominique Foata, Ira Gessel, Yakar Kannai, Steve Milne, Ron Perline, Marci Perlstadt, Joel Spencer, Dennis Stanton, Jet Wimp, and Herbert Wilf.

2. A Short Course On Holonomic Systems

2.1 C-finite Functions and Sequences

2.1.1 C-finite functions

We all know that solutions of (homogeneous) linear ordinary differential equations with constant coefficients:

P (D)f ≡ 0,  (2.1)

where P is a polynomial in D, can be expressed as a finite linear combination of exponential polynomial solutions. Indeed, any solution of (2.1) can be expressed as

f(x) = ∑_{r=1}^{R} p_{r}(x)e^{λ_{r}x},  (2.2)

where λ_{r} are the roots of the characteristic equation P (z) = 0, and the degree of p_{r}(x), r = 1, …, R, is one less than the multiplicity of the root λ_{r}.

Let us call C-finite functions such solutions of constant coefficients linear ordinary differential equations. It is easily seen that conversely, any function that can be written in the form (2.2) is a C-finite function: P (D)f = 0, where

P (D) = ∏_{r=1}^{R} (D − λ_{r})^{degp_{r}+1}  (2.3)

In order to specify a C-finite function, one only needs a finite number of parameters. One way of coding such functions is by (2.2): specifying the complex numbers λ_{r} and the polynomials p_{r}(x). Another way is directly in terms of (2.1); specify the operator P(D), of order N, say, and the N initial conditions f(0), f′(0), ..., f^{(N−1)}(0).

In the foregoing we took P(D) to be with arbitrary complex coefficients, and we made the tacit assumption that it is always possible to ”write down” an arbitrary complex number. This is of course wrong, since real and complex numbers are merely fictitious entities, and infinite storage space is needed to specify an arbitrary real or complex number. In order to be rigorous we must restrict the coefficients of P(D) to be integer, rational, or at least to belong to some finite algebraic extension field of the rationals. Then the λ_{r} are algebraic numbers, and algebraic numbers are specifiable by a finite number of bits.

We will no longer dwell on this subtlety. All the results will be true for an arbitrary field of characteristic zero. However, whenever we talk about something being ”finite” we will tacitly assume that the members of the field are finitely specifiable.

The class of C-finite functions is an algebra : the sum and product of two expressions of the form (2.2) is again of the same form. This can be also proved directly from (2.1). It is easily seen that f is C-finite if and only if the vector space spanned by D^{i}f; i ≥ 0 is a finite dimensional vector space. The vector spaces spanned by D^{i}(f + g); i ≥ 0 and D^{i}(fg); i ≥ 0 are subspaces of the direct sum and (the homomorphic image of the) tensor product respectively of the finite dimensional vector spaces D^{i}f; i ≥ 0 and D^{i}g; i ≥ 0. Thus if f and g are C-finite so are f + g and fg.

Let us summarize the three definitions of C-finite functions:

Definition 0: f is C-finite if there exists a polynomial P such that P (D)f ≡ 0.

Definition 1 : f is C-finite if the vector space C[D]f = span{D^{i}f; i ≥ 0} is finite dimensional.

Definition 2 : f is C-finite if it can be expressed as a linear combination of exponential polynomial functions, as in (2.2).

In order to motivate holonomic systems, we will introduce some notions that will seem like ”abstract nonsense” in the present context.

C[D] is the ring of polynomials in D, (that is the ring of constant coefficients ordinary linear differential operators). Consider the following left ideal in C[D]:

I_{f} = {P ϵC[D] ; P (D)f = 0} .  (2.4)

Of course since C[D] is a principal ideal ring, we can write I_{f} = P_{0}(D)C[D], where P_{0}(D) is the minimal order operator that satisfies P_{0}(D)f = 0. Both C[D]f and C[D]/I_{f,} have obvious structure as ”C[D] modules”: Q(D)(P (D)f) = (Q(D)P (D))f and Q(D)P̄(D) = Q̄P (D), where P̄ is the image of P under the natural mapping , C[D] → C[D]/I_{f}. It is readily seen that these two C[D] modules are in fact isomorphic, since P_{1}(D)f = P_{2}(D)f if and only if P_{1} − P_{2} ϵI_{f}.

We thus have

Definition 1’: f is C-finite if C[D]/I_{f} is finite dimensional.

For every function f, let us introduce the ”variety associated with the ideal I_{f}”:

V_{f} := {zϵC; P (z) = 0 for every P ϵI_{f}}  (2.5)

Let us consider an arbitrary function f on the real line. Most likely it is not a solution of any constant-coefficients linear differential equation, in which case the ideal I_{f} is the trivial ideal consisting of 0 alone, and C[D]/I_{f} is equal to C[D], and V_{f} is the whole of C, and therefore it is a ”one dimensional variety”. The other extreme is f ≡ 0, in which case I_{f} = C[D] and C[D]/I_{f} is the zero module, and V_{f} is empty, i.e. it is a ” −∞ dimensional variety”.

Between these two extremes are the C-finite functions for which V_{f} is a finite set of points, i.e. a ”zero dimensional variety” . Thus we have

Definition 3: f is C-finite if V_{f} is a zero dimensional variety (i.e. it is a finite set of points).

Since the class of C-finite functions is an algebra, and every C-finite function is finitely-specifiable, it follows that every identity involving sums and products of C-finite functions is routinely verifiable.

Trivial Example: Prove that sin(x + a) = sin(x)cos(a) + cos(x)sin(a). Set f(x) := sin(x + a) − sin(x)cos(a) − cos(x)sin(a), we have to show that f(x) ≡ 0. First Solution:

f(x) = \frac{e^{i(x+a)} − e^{−i(x+a)}}{2i} − \frac{e^{ix} − e^{−ix}}{2i} \frac{e^{ia} + e^{−ia}}{2} − \frac{e^{ix} + e^{−ix}}{2} \frac{e^{ia} − e^{−ia}}{2i} = 0.

Second solution: f(x) satisfies (D^{2} + 1)f = 0, and f(0) = 0 , and f′(0) = 0, so f(x) ≡ 0. 2.1.2 C-Finite Sequences Analogous considerations apply to the class of sequences that are solutions of some (homogeneous) ordinary linear recurrence equation with constant coefficients:

P (E)a ≡ 0,  (2.6)

where P is a polynomial, and E is the shift operator: Ea_{n} := a_{n+1}. For example the Fibonacci numbers F_{n} constitute such a sequence since (E^{2} − E − 1)F ≡ 0. Let us call such sequences C-finite sequences. It is well known that a sequence is C-finite if and only it can be written in the form

a_{n} = ∑_{r=1}^{l} p_{r}(n)z_{r}^{n}  (2.7)

where the z_{r} are the roots of P (z) = 0 and degp_{r}= mult(z_{r}) − 1. Everything we said before has its discrete analog. In particular, we can define the ideal in C[E]:

I_{a} := {P ϵC[E]; P (E)a = 0}.

and the notion of C-finite sequences can be defined by either one of the following four equivalent definitions: Definition 0 : A sequence a is C-finite if there exists a polynomial P such that P (E)a = 0 . Definition 1 : a sequence a is C-finite if the vector space C[E]a = span{E^{i}a; i ≥ 0}is finite dimensional. Definition 1’ : a sequence a is C-finite if the C[E] module C[E]/I_{a} is a finite dimensional vector space. Definition 2 : a sequence a is C-finite if it is of the form (2.7). Definition 3 : a sequence a is C-finite if

V_{a} := {zϵC; P (z) = 0 for every P in I_{a}}

is a finite set of points (i.e. a ”zero dimensional variety”).

It follows similarly that the class of C-finite sequences is an algebra and that every identity involving sums and products of C-finite sequences is routinely verifiable.

Trivial Example: Prove Cassini’s identity

F_{n+1}F_{n−1} − F_{n}^{2} = (−1)^{n}

where F_{n} are the Fibonacci numbers defined by F_{0} = 0 , F_{1} = 1, and F_{n+2} = F_{n+1} + F_{n}.

Set a(n) := F_{n+1}F_{n−1} − F_{n}^{2} − (−1)^{n}, we have to prove that a(n) ≡ 0.

First solution: Express F_{n} in the form (2.7), in terms of the ”golden ratio”, and do the trivial, but tedious algebra.

Second Solution: It can be easily seen that both F_{n}^{2} and F_{n+1}F_{n−1} satisfy the same third order recurrence equation and therefore a(n) satisfies a certain fourth order linear recurrence equation with constant coefficients, whose exact form is irrelevant. It follows that it is enough to check the four initial conditions a(0)=0,a(1)=0,a(2)=0,a(3)=0, which is easily done by inspection.

2.2 Multi C-finite Functions and Sequences

2.2.1 Multi C-finite Functions

What is the several variable analog of C-finite functions? The straightforward analog of Definition 0 is ”a solution of a linear partial differential equation with constant coefficients”. However, it is soon realized that all the other definitions do not go through. Furthermore, the ”finiteness property” is not preserved. In order to specify, say, a solution of (D_{1}^{2} + D_{2}^{2})f = 0, one needs to state the ”Dirichlet boundary conditions” on the boundary of a closed region, i.e. one has to furnish an infinite amount of information to specify such a function.

The fact that every solution of an ordinary linear differential equation with constant coefficient is a linear combination of exponential polynomial solutions has a celebrated analog in several variables. It is called the Ehrenpreis-Palamadov ([Ehr1], [Pal]) theorem. This theorem implies that with an appropriate definition of convergence, every solution of a constant coefficients linear partial differential equation can be expressed as an (infinite) ”linear combination” of exponential polynomial solutions. This theorem generalizes to so-called over-determined systems

P_{i}(D_{1}, ..., D_{n})f ≡ 0, i = 1, …, L,  (2.8)

where P_{i}(D_{1}, ..., D_{n}) are polynomials in D_{1}, ..., D_{n}. The Ehrenpreis-Palamadov theorem implies that every solution of (2.8) can be expressed as a (usually infinite) ”linear combination” of exponentialpolynomial solutions.

An exponential polynomial solution f(x) = p_{λ}(x)e^{λx} = p_{λ}(x) e^{λ_{1}x_{1}+...λ_{n}x_{n}} is a solution of (2.8) only if

P_{i}(λ_{1}, ..., λ_{n}) = 0, i = 1, ..., L.

Thus the Ehrenpreis-Palamadov theorem asserts that every solution of the system (2.8) is a (usually infinite) ”linear combination” of exponential polynomial functions p_{λ}(x)e^{λx} where λ ranges over the algebraic variety

V := {λ ∈ C^{n}; P_{i}(λ) = 0, i = 1, ..., L}.

Guided by Definition 2 of C-finite functions, we want the ”usually infinite linear combination” to be a good old fashioned finite linear combination, and thus make our entities ”specifiable by a finite amount of information”. This prompts us to define

Definition 2: f is multi C-finite if it can be written as

f(x) = ∑_{r=1}^{R} p_{r}(x)e^{λ^{(r)}x}.  (2.9)

In this case the variety V consist of a finite set of points {λ^{(1)}, …, λ^{(R)}}, and so is a ”zerodimensional variety”. Thus we have a natural analog of definition 3: f is C-finite if it is a solution of a system (2.8) whose associated variety V_{f} is zero dimensional. More precisely let us define, as before,

I_{f} := {P ∈ C[D_{1}, ..., D_{n}]; P f = 0},  (2.10)

and

V_{f} := {λϵC^{n}; P (λ) = 0 for every P ∈ I_{f}},  (2.11)

then we have:

Definition 3 : f : R^{n} → C is multi C-finite if the algebraic variety V_{f} is zero dimensional , i.e. is a finite set of points.

Note that for an arbitrary f, I_{f} is usually the 0 ideal and then V_{f} is the whole of C^{n}, i.e. it is an ”n-dimensional variety”. In general the dimension of V_{f} is between 0 and n, and being multi C-finite means that the function satisfies as many independent constant coefficients linear partial differential equations as possible. Thus being multi C-finite is the same as being a solution of a ”maximally over-determined system”. Of course we have again excluded the trivial case of f ≡ 0 in which case I_{f} = C[D] and V_{f} is the empty set.

What about the analog of Definition 0? It is easily seen that now there are n differential operators P_{1}, ..., P_{n}, each of which is ”ordinary” in D_{1}, ...D_{n} respectively, such that:

P_{i}(D_{i})f = 0, i = 1, …, n .  (2.12)

They are given by:

P_{i}(D_{i}) := ∏_{r=1}^{R} (D_{i} − λ_{i}^{(r)})^{deg(p_{r})+1}  (2.13)

where

λ^{(r)} = (λ_{1}^{(r)}, …, λ_{n}^{(r)}).

Let P_{i}(D_{i}) be of order α_{i}. It follows from (2.12) that D_{i}^{α_{i}}f can be expressed as a linear combination of f, D_{i}f, ..., D_{i}^{α_{i}−1}f, and by iterating , D_{i}^{j}f can be thus expressed for any j ≥ α_{i}. Similarly, it is readily seen that every D_{1}^{i_{1}} ⋯ D_{n}^{i_{n}}f can be expressed as a linear combination of such entities for which i_{1} < α_{1}, ..., i_{n} < α_{n}. We thus have

Definition 1 : f is multi C-finite if C[D_{1}, ..., D_{n}]f is finite dimensional.

Definition 1’: f is multi C-finite if the C[D_{1}, ..., D_{n}]-module, C[D_{1}, ..., D_{n}]/I_{f}, is a finite dimensional vector space.

When we talked before about the ”dimension of the variety”, we meant the usual, analytical dimension, as a complex manifold. However there are several other, algebraic, notions of dimension, for example the Krull dimension . The notion of dimension that was used by I.N.Bernstein, and which is the one that we will need, is that of the Hilbert dimension of a graded C[z_{1}, ..., z_{n}] module, which we shall now introduce.

The ring of polynomials C[z_{1}, ..., z_{n}] with coefficients in C has a natural ”filtration”:

C[z_{1}, ..., z_{n}] = ∪_{v}F_{v}  (2.14)

where, F_{v} is the vector space of all polynomials of (total) degree ≤ v. Since a basis for this vector space is z_{1}^{α_{1}}...z_{n}^{α_{n}},; , α_{i} ≥ 0, α_{1} + ...α_{n} ≤ v, its dimension is the number of solutions in non-negative integers of α_{1} + ...α_{n} ≤ v which is well known, and easily seen, to be (n + v)!/n!v!. This is a polynomial in v of degree n with leading coefficient 1/n!.

A graded C[z_{1}, ..., z_{n}] module is a C[z_{1}, ..., z_{n}] module S which has the decomposition

S = ⨁ S(v),

(direct sum)

where S(v) are C-subspaces of S and x_{j}S(v) ⊂ S(v + 1) for 1 ≤ j ≤ n and all v. A famous theorem of Hilbert (e.g. [Bjo1], p.7) states that for v >> 0 ∑_{j≤v} dim_{C}(S(j)) is always a polynomial in v with rational coefficients (called the Hilbert polynomial). The degree of this polynomial is called the Hilbert dimension of S and is a measure of ”how big” S is. Of course C[z_{1}, ...z_{n}] itself is a graded C[z_{1}, ..., z_{n}] module with

S(j) = span{z_{1}^{α_{1}} ⋯ z_{n}^{α_{n}}; α_{1} + … + α_{n} = j}

We saw that ∑_{j≤v} dim(S(j)) = dim(F_{v}) = \binom{n+v}{n}, a polynomial of degree n in v. So the Hilbert dimension of C[z_{1}, ..., z_{n}] coincides with its natural dimension, n.

More generally, for any ideal I in C[z_{1}, ..., z_{n}], C[z_{1}, ..., z_{n}]/I has a natural structure as a graded C[z_{1}, ..., z_{n}] module, inherited from (2.14) by modding out by I. If the ideal I is the 0 ideal, then the Hilbert dimension is n. Excluding this trivial case, the dimension is between 0 and n-1. The Hilbert dimension of C[z_{1}, ..., z_{n}]/I is zero if and only if it is a finite dimensional C-vector space, because then ∑_{j≤v} dim_{C}(S(j)) is eventually constant. We thus have

Definition 1”: f is multi C-finite if the graded C[z_{1}, ..., z_{n}] module C[z_{1}, ..., z_{n}]/I_{f} has Hilbert dimension 0, which is the smallest possible.

Of course we have excluded f ≡ 0 for which the Hilbert polynomial is identically 0, which is a ” polynomial of degree −∞”.

It can be seen in a variety of ways that the class of multi C-finite functions is an algebra. The easiest way is by using definition 2. However, it is also possible to use Def 1, like we did with ordinary C-finite functions. C[D_{1}, ..., D_{n}](f + g) is a subspace of C[D_{1}, ..., D_{n}]f ⨁ C[D_{1}, ..., D_{n}]g and C[D_{1}, ..., D_{n}](fg) is a subspace of a homomorphic image of C[D_{1}, ..., D_{n}]f ⨂ C[D_{1}, ..., D_{n}]g and so are finite dimensional vector spaces.

Another obvious property, that we will generalize later on, is that the diagonal f̄(x̄) of a multi Cfinite f: f̄(x̄) := f(x̄, ..., x̄) is an ordinary C-finite function in the single variable x̄. This is immediate from Definition 2 but is also not hard to deduce using any of the other definitions. For example, from definition 1 it follows that the vector space span {D_{1}^{r}...D_{n}^{r}f; r ≥ 0} is finite dimensional and hence there is an operator P̄(D̄) in D̄ := D_{1}...D_{n} such that annihilates f(x_{1}, ..., x_{n}), and thus P (D̄_{x})f̄(x̄) = 0.

2.2.2 Multi C-finite Sequences

Everything in the previous section has an obvious discrete analog. Now we have multi sequences that are nothing but functions from Z^{n} to C, or N^{n} to C, or for that matter one can consider functions from any subset of Z^{n} to C. Given a function a : Z^{n} → C, we define the fundamental shift operators

E_{i}a(m_{1}, ..., m_{i}, ..., m_{n}) := a(m_{1}, ..., m_{i} + 1, ..., m_{n})  (2.15)

A linear partial recurrence operator with constant coefficients is a polynomial P (E_{1}, ..., E_{n}) in the fundamental shift operators. For example the operator a(m_{1}, m_{2}) → a(m_{1}, m_{2})+ 3a(m_{1}+1, m_{2}+2) is written I + 3E_{1}E_{2}^{2}. The ring of linear partial recurrence operators with constant coefficients is in fact C[E_{1}, ..., E_{n}], the commutative polynomial algebra in the n indeterminates E_{1}, ..., E_{n} over the field C. We leave it to the readers to do the obvious discrete analog of the previous subsection, and give the definitions of the notion ”multi C-finite” as applied to discrete functions. For example:

Definition 1: a : Z^{n} → C is multi C-finite if C[E_{1}, ..., E_{n}]a is a finite dimensional vector space.

2.2.3 Mixed Continuous and Discrete Functions

Consider functions f : R^{n_{1}} × Z^{n_{2}} → C. Then C[D_{1}, ..., D_{n_{1}}, E_{1}, ..., E_{n_{2}}] acts in a natural way on these functions, and everything goes through (except the definition of the diagonal).

2.3 P-finite Functions and Sequences

If, in Definition 0 of C-finite functions, we replace ”a solution of a constant coefficients linear ordinary differential equation” by ”a solution of a polynomial coefficients linear ordinary differential equation” we get the class of P-finite functions. Similarly, a P-finite sequence is a sequence that is a solution of a (homogeneous) linear ordinary recurrence equation with polynomial coefficients. For example, the Apéry numbers

a_{n} := ∑_{k=0}^{n}({n \atop k})^{2}({n + k \atop k})^{2}

of ζ(3) fame ([vdP]) is a P-finite sequence since it satisfies

n^{3}a_{n} − (34n^{3} − 51n^{2} + 27n − 5)a_{n−1} + (n − 1)^{3}a_{n−2} = 0.

P-finite functions and sequences were introduced and studied by Stanley[Stan]. Please note that our P-finite functions are called by him ”D-finite functions”, and our P-finite sequences are called by him ”P-recursive sequences”. In the interest of thawing the cold war between the discrete and the continuous I have decided to combine these two names into one.

So far we have introduced P-finite functions through the analog of Definition 0. Do the other definitions go through? Unfortunately, Definition 2 does not go through, but all the other definitions do. For example,

Definition 1: f is P-finite if C[D]f is a finite dimensional vector space over the field of rational functions C(x).

The analogs of definitions 1’,1”, and 3 will be discussed in the next subsection, in the more general context of holonomic functions.

Both P-finite functions and P-finite sequences are closed under addition and multiplication ([Stan]). This is easily seen from Definition 1, since the C(x) vector space C[D](f+g) is a subspace of C[D]f ⨁_{C(x)} C[D]g and the C(x) vector space C[D](fg) is a subspace of a homomorphic image of C[D]f ⨂_{C(x)} C[D]g , and thus are finite dimensional over C(x).

Once again, the discrete analog of this is immediate.

2.4 Holonomic Systems and Holonomic functions

After this somewhat lengthy introduction we are finally ready to introduce holonomic systems, that form the foundation of our approach to special function identities.

Holonomic systems were introduced and studied by I.N.Bernstein ([Ber1],[Ber2]) and were used by him to give an elementary proof of a famous conjecture of Gelfand concerning the existence of a meromorphic extension of the distribution-valued complex function λ → P^{λ}, where P is a polynomial of several variables in R^{n}. This theory is very deep and stands at the forefront of current research in analysis. Fortunately we will only need the most basic notions from this imposing theory, and everything we will need will be stated and explained. An excellent introduction to the theory of holonomic systems can be found in the introduction and chapters 1 and 7 of Bjork’s monograph [Bjo1]. Much of the following follows closely Bjork’s readable account on pages ix and x of [Bjo1].

Holonomic systems are obtained when the linear partial differential operators with constant coefficients that feature in the definition of multi C-finite functions are replaced by linear partial differential operator with polynomial coefficients. Thus instead of the commutative algebra C[D_{1}, ..., D_{n}] we will have to consider the so-called Weyl algebra A_{n}(C) = C < D_{1}, ...D_{n}, x_{1}, ..., x_{n} > which is the (non-commutative) algebra generated by the indeterminates D_{1}, ..., D_{n}, x_{1}, ..., x_{n}. Here x_{i} denotes the operator of multiplication by x_{i} : f → x_{i}f. These generators satisfy the following commutation relations: [x_{i}, x_{j}] = 0, [D_{i}, D_{j}] = 0, and from the product rule of differentiation [D_{i}, x_{j}] = δ_{i,j}, where δ_{i,j} is Kronecker’s delta function that equals 1 when i = j and equals 0 otherwise. (recall that ”1” means the identity operator : f → f).

The non-commutative structure is already so strong that the ring A_{n}(C) has no two-sided ideals, except the zero-ideal and the whole ring A_{n}(C). But of course there are one-sided ideals, left or right, and more generally left and right A_{n}(C) -modules.

A remarkable result, due to J.T.Stafford ([Staf], [Bjo1](1.7)), asserts that every left ideal in A_{n}(C) can be generated by 2 elements.

In particular A_{n}(C) is a left and right notherian ring. The main idea is to use filtrations on the ring A_{n}(C) and its modules.

The ring A_{n}(C) is equipped with the F -filtration= ∪F_{v} , where F_{v}= the finite dimensional vector space generated by the (x, D)-monomials x^{α}D^{β} whose total weight |α| + |β| ≤ v.

Thus F_{0} ⊂ F_{1} ⊂ F_{2} ⊂ ... is an increasing sequence of subspaces of A_{n}(C). Of course ∪F_{v} = A_{n}(C) and the inclusion F_{v}F_{k} ⊂ F_{v+k} for all pairs of non-negative integers v and k, follows from the definition of the ring product of A_{n}(C).

Let us now consider a left ideal L in the ring A_{n}(C). If v ≥ 0 we can consider the dimensions of the complex vector spaces F_{v}/(F_{v} ∩ L), which form the induced filtration on the left A_{n}(C) module A_{n}(C)/L. Let H(v) be the dimension of F_{v}/F_{v} ∩ L . Then H(v) = dim_{C}(F_{v}) − dim_{C}(F_{v} ∩ L). It turns out ([Bjo1], 1.3) that H(v) is a polynomial function a_{d}v^{d} + ... + a_{0} when v >> 0. Here a_{0}, ...ma_{d} are rational numbers and (d!)a_{d} is a positive integer that is denoted by e.

The degree d of the Hilbert polynomial a_{d}v^{d} + ... + a_{0} which computes H(v) when v >> 0, is called the F -dimension of the left A_{n}(C) module A_{n}(C)/L and is denoted by d_{F} (A_{n}/L).

The celebrated Bernstein inequality asserts that for every left ideal L , d_{F} (A_{n}/L) ≥ n, unless , of course L is the whole ring A_{n}.

The above notion of dimension can be extended to general finitely generated left A_{n}(C) modules of which modules of the form A_{n}(C)/L are special cases. One considers a certain class of ”good filtrations”, and shows that the notion of dimension makes sense for each such good filtration. Then it is shown that the dimension is independent of the particular filtration, and thus one is able to talk about the dimension of the A_{n}(C) module M, d(M), without reference to any particular filtration. Bernstein([Ber1]; [Bjo1], 1.4) proved the more general statement that d(M) ≥ n for every non-zero and finitely generated left A_{n}(C) module. An extremely short and elegant proof of this important inequality was given by Antony Joseph, and can be found in [Bjo2], and in [Gal](p. 417).

Let f be either a C^{∞} function, or a distribution, or a formal power series, or what have you(anything on which it is possible to differentiate and multiply by x_{i}), in n variables. Then one can consider all the elements of A_{n}(C) that annihilate f:

I_{f} := {P ∈ A_{n}(C); P f = 0},  (2.16)

Obviously I_{f} is a left ideal. By Bernstein’s inequality the dimension of the A_{n}(C) module is ≥ n, unless f is identically zero. We will call functions f for which the dimension of A_{n}(C)/I_{f} is the smallest possible, namely n, holonomic functions (or distributions, formal power series etc.) It can be easily seen that multi C-finite functions are ipso facto holonomic functions . This makes sense since in the passage from C[D_{1}, ..., D_{n}]/I_{f} to A_{n}(C)/I_{f} we have acquired n ”extra dimensions”, corresponding to the generators x_{1}, ..., x_{n}, and the ”gap” between numerator and denominator is still n.

Because of the analogy to multi C-finite functions and P-finite functions it would make sense to call holonomic functions ”multi P-finite functions”. Bjork [Bjo1] calls holonomic functions ”members of the Bernstein class”.

The Weyl algebra A_{n}(C) acts just as well on discrete functions f : Z^{n} → C, via x_{i} → E_{i}^{−1} and D_{i} → (n_{i} + 1)E_{i}, as can be easily seen by examining the action of x_{i} and D_{i} on a generic formal power series, and see how the coefficients get transformed. Of course this extends to mixed continuous-discrete functions f(x_{1}, ..., x_{n}, m_{1}, ..., m_{n′}) on which the Weyl algebra A_{n+n′}(C) acts naturally.

So let us state the formal definition of holonomic functions(or distribution, formal power series, etc.)

Definition 1” : Let f be a non-zero member of a family on which the Weyl algebra acts naturally. Let the ideal I_{f} be defined by (2.16) above. Then f is holonomic if the left A_{n}(C) module A_{n}(C)/I_{f} has its smallest possible dimension, namely n.

It is readily seen that this definition is the natural extension of definition 1” of multi C-finite functions, in section 2.2.1. Since A_{n}(C)/I_{f} is naturally isomorphic, as A_{n}(C) modules, to the A_{n}(C) module A_{n}(C)f. We can also state

Definition 1’ : f is holonomic if the A_{n}(C) module A_{n}(C)f has the smallest possible dimension, namely n.

Definition 3 can also be extended to the present context. But we first need some notation. A typical element of the Weyl algebra A_{n}(C) may be written

P = ∑_{|α|≤N} a_{α}(x)D^{α},

where α = (α_{1}, ..., α_{n}), D^{α} := D_{1}^{α_{1}}...D_{n}^{α_{n}}, and a_{α}(x) = a_{α}(x_{1}, ..., x_{n)} are polynomials. Its symbol is obtained by replacing D_{i} by the complex variables ζ_{i}, and replacing x_{i} by the complex variables z_{i}, obtaining a polynomial in the 2n complex variables z_{1}, ..., z_{n}, ζ_{1}, ..., ζ_{n}. The principal symbol, σ(P) is obtained by only retaining the leading terms for which |α| = N:

σ(P) = ∑_{|α|=N} a_{α}(z)D^{α} ,

For any left ideal L in A_{n}(C), lLet us consider the following variety in C^{n} × C^{n} = C^{2n}:

V_{L} = {(z, ζ)ϵC^{n} × C^{n}; σ(Q)(z, ζ) = 0, for every Q in L} = ∩{σ(Q)^{−1}(0) ; Q ∈ L}

It is a deep result (see [Bjo1],x, and chapter 2) that the (complex analytic) dimension of the variety V_{L} coincides with the previously defined dimension of the A_{n}(C) module A_{n}(C)/L, for any left ideal L in A_{n}(C). In particular we can define

Definition 3 : Let f be a non-zero member of a family on which the Weyl algebra acts naturally. Then f is holonomic if the variety

V_{f} := {(z, ζ) ∈ C^{n}×C^{n}; σ(Q)(z, ζ) = 0, for every Q such that Q(f) = 0} = ∩{σ(Q)^{−1}(0); Q(f) = 0}

has dimension n.

2.5 Examples of Holonomic Functions

In the next section we will show that the class of holonomic functions is an algebra: the sum and product of holonomic functions are also holonomic. We will also prove that it is closed under integration (or summation) with respect to a variable: If f(x_{1}, ..., x_{n}) is a holonomic function of n variables, then ∫_{−∞}^{∞} fdx_{n} is holonomic in the surviving n − 1 variables x_{1}, ..., x_{n−1}, if it is defined. Similarly if f = f(m_{1}, ..., m_{n}) : Z^{n} → C is holonomic then ∑_{m_{n}} f is holonomic in the surviving variables m_{1}, ..., m_{n−1}.

All polynomials and exponential-polynomial functions, being multi C-finite, are automatically holonomic. The simplest discrete function of a single variable that is P-finite but not C-finite is f(n) := n!, that satisfies (E − n − 1)f = 0. Similarly f(n) := 1/n! is P-finite since it satisfies ((n + 1)E − 1)f = 0 .

The function f := m_{1}! is a holonomic function of m_{1}, ..., m_{n}. One way of seeing this is noting that I_{f} contains the operators I −E_{i} , 2 ≤ i ≤ n, as well as E_{1} −m_{1} −1 . When translating this becomes 1 − z_{i}, 2 ≤ i ≤ m , and z_{1} − z_{1}D_{1} − 1, whose principal symbols are −z_{i} and −z_{1}ζ_{1}. The variety V_{f} is a subset of the set of common zeroes of −z_{i}, 2 ≤ i ≤ n and z_{1}ζ_{1}, which is a union of the two n dimensional varieties in C^{2n}: z_{1=}... = z_{n} = 0 and ζ_{1} = 0, z_{2} = ... = z_{n} = 0, and thus its dimension must be ≤ n, and by Bernstein’s inequality it should be equal to n, and thus is holonomic.

In a similar way we can prove that 1/m_{1}! is holonomic and of course this is true for all m_{i}! and 1/m_{i}! . Similarly (m_{1} +...+m_{n})! is holonomic and thus the multinomial coefficient (m_{1} +...+m_{n})!/ m_{1}!...m_{n!}. This would also follow from the following general result, upon taking P = 1−x_{1}−...−x_{n}, and looking at the coefficients of 1/P, viewed as a formal power series.

Proposition 2.1: Let P = P (x_{1}, ..., x_{n}) be a polynomial in n variables then 1/P is holonomic, whenever it is defined.

First Proof: Let us try and find as many linear differential operators with polynomial coefficients that annihilate f:=1/P as possible. Among the many such operators are P D_{1} + P^{(1)} and P^{(i)}D_{1} − P^{(1)}D_{i}, for 2 ≤ i ≤ n. (P^{(i)} is the partial derivative of P with respect to x_{i}: ∂P/∂x_{i}.) The principal symbols of these are P ξ_{1} and P^{(j)}ξ_{1} − P^{(1)}ξ_{j}. The variety V_{f} is a subset of the set of common zeroes of these n polynomials in x_{1}, ..., x_{n}, ξ_{1}, ..., ξ_{n}, which is easily seen to comprise an n dimensional variety. Thus the dimension of the variety V_{f} is ≤ n, and by Bernstein’s inequality it must be equal to n. It follows that 1/P is holonomic.

Second proof: A_{n}(C)(1/P) is a submodule of C[x_{1}, ..., x_{n}, P^{−1}], and this last A_{n}(C) module is shown to be a holonomic A_{n}(C)module in Bjork[Bjo1], Theorem I.5.5 (p.13). .

Since every polynomial Q is obviously holonomic it follows (from the closure under product to be proved in section 3) that all rational functions Q/P are holonomic. If P has a non-zero constant term, then Q/P makes sense as a formal power series, and its multi sequence of coefficients f defined by

Q/P = ∑ f(m_{1}, ..., m_{n})x_{1}^{m_{1}}...x_{n}^{m_{n}}

is a holonomic discrete function from N^{n} to C.

3 Operations That Preserve Holonomicity

3.1 Continuous Functions and Distributions

We will now show that the sum and product of holonomic functions are themselves holonomic , and that if f(x_{1}, ..., x_{n}) is holonomic on R^{n}, then f(x_{1}, ..., x_{n−1}, c) ( c a constant) and ∫ f(x_{1}, ..., x_{n})dx_{n} are holonomic functions of (x_{1}, ..., x_{n−1}).

Proposition 3.1 : Let f : R^{n} → C, g : R^{n} → C be holonomic, then f + g is also holonomic.

Proof: Let A_{n}(C) = ∪_{v=0}^{∞}F_{v} be the natural filtration of A_{n}(C), described above. Then the natural filtrations for A_{n}(C)f, A_{n}(C)g, and A_{n}(C)(f +g) are ∪_{v=0}^{∞}F_{v}f, ∪_{v=0}^{∞}F_{v}g, and ∪_{v=0}^{∞}F_{v}(f +g). Now for every v we have F_{v}(f + g) ⊂ F_{v}f ⨁ F_{v}g. Thus dimF_{v}(f + g) ≤ dim(F_{v}f) + dim(F_{v}g). Since f and g are holonomic, we know that for v >> 0 both dim(F_{v}f) and dim(F_{v}g) are polynomials in v of degree n, and the same is true for their sum. We do know a priori that for v >> 0 dimF_{v}(f +g) is a polynomial, and so it follows that its degree must be ≤ n, and therefore, by Bernstein’s inequality it is equal to n, and f + g is holonomic.

Proposition 3.2: If f and g are holonomic so is their product (if it is defined).

Proof: Let M_{1} and M_{2} be A_{n}(C) modules. Bernstein ([Ber2],p.277) defines the A_{n}(C) module M_{1} ⨂ M_{2} as follows: as an C[x_{1}, ..., x_{n}] module M_{1} × M_{2}= M_{1} ×_{C[x]} M_{2} and the operators ∂/∂x_{i} acts as follows

\frac{∂}{∂x_{i}} (f_{1} × f_{2}) = (\frac{∂}{∂x_{i}} f_{1}) × f_{2} + f_{1} × (\frac{∂}{∂x_{i}} f_{1}) , (f_{1}ϵm_{1}, f_{2}ϵM_{2}) .

Bernstein proved ([Ber2], p.278, theorem 3.2, part 3 ) that if M_{1} and M_{2} are holonomic A_{n}(C) modules, so is M_{1} × M_{2}.

Now, by Leibnitz’s rule

A_{n}(C)(fg) ⊂ (A_{n}(C)f) × (A_{n}(C)g).

It follows that A_{n}(C)(fg) is a holonomic A_{n}(C) module, since it is a non-zero submodule of a holonomic A_{n}(C) module.

Proposition 3.2 has some far-reaching consequences. It is obvious that the Dirac delta measure δ_{a}, concentrated at x = a (i.e. δ_{a}(f) = f(a), for every test function) is a holonomic distribution in a single variable, since it satisfies the equations (of order zero) (x − a)δ_{a} = 0. When δ_{a} is embedded in R^{n} as the Lebesgue measure that is supported at the hyperplane x_{n} = a, it is still holonomic, but now in R^{n}. Multiplication by δ_{a} from the right is nothing but the evaluation x_{n} = a:

f(x_{1}, ..., x_{n−1}, x_{n})δ_{a} = f(x_{1}, ..., x_{n−1}, a) , when x_{n} = a, and 0 otherwise.

It follows from proposition 3.2 that if f is holonomic and fδ_{a} is defined (which is the case, for example, when f is continuous) then it is holonomic. In other words, if f(x_{1}, ..., x_{n−1}, x_{n}) is holonomic in R^{n} then f(x_{1}, ..., x_{n−1}, a) is holonomic in R^{n−1}.

It turns out ([Bjo1],7.4.1) that it is always possible to define the evaluation x_{n} = a on any holonomic distribution. This is certainly not possible for an arbitrary distribution, so the class of holonomic distributions is indeed privileged.

The differentiation operators \frac{1}{i} \frac{∂}{∂x_{j}} and the ”multiplication by x_{j}” operators x_{j} are dual to each other under the Fourier transform, so it is easily seen ([Bjo1], p.293) that the Fourier transform of a holonomic entity on R^{n}, if it exists (and it does if it is a tempered distribution ), is also holonomic, and vice versa. The Fourier transform sends the act of multiplication into the act of convolution: ˆ(fg) = f̂∗ ĝ, so we have as an immediate consequence of proposition 3.2 :

Proposition 3.2: If f and g are holonomic so is their convolution (if it is defined).

This proposition can be also easily proved directly, by mimicking the proof of proposition 3.2.

Let δ_{n}^{(n−1)} be the Lebesgue measure of the ”line” x_{1} = ... = x_{n−1} = 0 in R^{n} , i.e. for any test function phi:

δ_{n}^{(n−1)}(ϕ(x_{1}, ..., x_{n})) = ∫ ϕ(0, …, 0, x_{n})dx_{n} .

It is holonomic because it satisfies the system:

x_{i}δ_{n}^{(n−1)} = 0, i = 1, …, n − 1

D_{n}δ_{n}^{(n−1)} = 0.

Now

f̄(x_{1}, ..., x_{n−1}) := ∫ f(x_{1}, ..., x_{n−1}, x_{n})dx_{n}  (3.1)

is nothing but the convolution f ∗δ_{n}^{(n−1)}, and so it follows from proposition 3.2 that if f is holonomic in R^{n}, then f̄, as given in (3.1) is holonomic in R^{n−1}.

Summarizing we have:

Proposition 3.3: Let f = f(x_{1}, ..., x_{n}) be holonomic in R^{n}, and let a be any constant, then f(x_{1}, ..., x_{n−1}, a) is holonomic in R^{n−1}.

Proposition 3.4: Let f(x_{1}, ..., x_{n}) be holonomic in R^{n}, then ∫ f(x_{1}, ..., x_{n})dx_{n} is holonomic in R^{n−1}, if it is defined.

It is obvious that the class of holonomic functions (or distributions) is closed under differentiation and indefinite integration. Furthermore, the closure under indefinite integration and evaluation implies closure under definite integration, so we have

Proposition 3.5: If f(x_{1}, ..., x_{n}) is holonomic, then, for any constants a and b,

∫_{a}^{b} f(x_{1}, ..., x_{n−1}, x_{n})dx_{n}

is also holonomic, if it is defined. 3.2 Discrete Functions Everything that we said about the continuous realm goes over smoothly to the discrete realm. Given a discrete function f : Z^{n} → C, we can talk about its Fourier transform

f̂(t_{1}, ..., t_{n}) = ∑_{mϵZ^{n}} f(m)e^{im_{1}t_{1}+...+im_{n}t_{n}}

that lives on the torus T^{n}=[−π, π)^{n}. Alternatively, by taking e^{it_{j}} → z_{j}, (j = 1, ..., n), we get the ”z-transform”

f̂(z_{1}, ..., z_{n}) = ∑_{mϵZ^{n}} f(m)z_{1}^{m_{1}}...z_{n}^{m_{n}}  (3.2)

(3.2) always makes sense as a ”formal Laurent series”, which may be viewed as ”distributions” with respect to the Laurent polynomials as ”test functions”. If f(m) is of polynomial growth then f̂ is a distribution on T^{n}, if f(m) is of exponential growth, then f̂ is analytic in some multi-annulus that contains T^{n}, and if f(m) is supported in N^{n}, then fˆ is a nice and honest formal power series. Let us apply the operators z \frac{d}{dz} and z on a single variable formal Laurent series

f̂(z) = ∑_{−∞}^{∞} f(m)z^{m}.

We get

z \frac{d}{dz} f̂(z) = ∑_{−∞}^{∞} mf(m)z^{m}.

and

zf̂(z) = ∑_{−∞}^{∞} f(m)z^{m+1} = ∑_{−∞}^{∞} f(m − 1)z^{m}

so the z-transform sends the operator z \frac{d}{dz} to the operator ”multiplication by m” , and it sends the operator ”z” into the ”backward shift” E^{−1} ( E^{−1}f(m) = f(m − 1) ). In several variables we have

z_{j} ↔ E_{j}^{−1} , z_{j}D_{j} ↔ m_{j}

where, as always m_{j} is shorthand for the operator ”multiplication by m_{j}”. A linear partial recurrence operator with polynomial coefficients is written

P (m_{1}, ..., m_{n}, E_{1}^{−1}, ..., E_{n}^{−1}),

and this corresponds to the differential operator

P (z_{1}D_{1}, ..., z_{n}D_{n}, z_{1}, ..., z_{n}).

Conversely any linear differential operator with polynomial coefficients can be written as

z_{1}^{−a_{1}}...z_{n}^{−a_{n}}P (z_{1}D_{1}, ..., z_{n}D_{n}, z_{1}, ..., z_{n})

for some non-negative integers a_{1}, ..., a_{n}. So it is obvious that f(m) is holonomic if and only if fˆ (z) is holonomic. In fact this last sentence is a tautology since we have basically defined f(m) to be holonomic when fˆ (z) is holonomic. The same reasoning as in the continuous case yields: Proposition 3.1’: the sum of two holonomic discrete functions from (a subset of) Z^{n} to C is also holonomic. Proposition 3.2’: the product of two holonomic discrete functions from (a subset of) Z^{n} to C is also holonomic. Proposition 3.2′ : If f and g are holonomic discrete functions so is their convolution (if it is defined). The characteristic function of the discrete hyperplane m_{n} = a is clearly holonomic, so we have: Proposition 3.3’: Let f = f(m_{1}, ..., m_{n}) be holonomic in Z^{n}, and let a be any constant, then f(m_{1}, ..., m_{n−1}, a) is holonomic in Z^{n−1}. The characteristic function of the discrete line m_{1} = ... = m_{n−1} = 0 is obviously holonomic so we have : Proposition 3.4’: Let f(m_{1}, ..., m_{n}) be holonomic in Z^{n}, then

∑_{m_{n}} f(m_{1}, ..., m_{n−1}, m_{n})

is holonomic in Z^{n−1}. (If it is defined.) Since the property of being holonomic is preserved under addition we have Proposition 3.5’: If f(m_{1}, ..., m_{n}) is holonomic in Z^{n}, then, for any constants a and b,

∑_{m_{n}=a}^{b} f(m_{1}, ..., m_{n−1}, m_{n})

is holonomic in Z^{n−1} (in the surviving variables m_{1}, ..., m_{n−1}).

3.2.1 The Diagonal of a Holonomic Formal Power Series is D-finite

If

f̂(z_{1}, ..., z_{n}) = ∑_{mϵN^{n}} f(m_{1}, ..., m_{n})z_{1}^{m_{1}}...z_{n}^{m_{n}}

is a formal power series, we define the diagonal as

f̂_{D}(w) = ∑_{m} f(m, ..., m)w^{m}.

In other words, we restrict f(m_{1}, ..., m_{n}) to the diagonal and replace z_{1}...z_{n} by a new variable w.

The z-transform of the characteristic set of the diagonal (m, ..., m), m ≥ 0 is the rational formal power series 1/(1 − (z_{1}...z_{n})) , which, by proposition 2.1 (section 2.5) is holonomic (this can be also checked directly). If f̂ is a holonomic formal power series then f itself must be a holonomic discrete function and it follows by proposition 3.2’ that the product of f with the characteristic function of the diagonal (m, ..., m), m ≥ 0 is holonomic. But this product is the restriction f(m, ..., m) of f to the diagonal, and so f(m, ..., m) is a single variable holonomic discrete function ( what Stanley[Stan] calls P-recursive), and its z-transform f_{ˆ D}(w) must be a holonomic formal power series in one variable, what Stanley[Stan] calls D-finite.

By proposition 2.1, a rational formal power series Q(z_{1}, ..., z_{n})/P (z_{1}, ..., z_{n}) ( P (0, ..., 0) ≠ 0) is holonomic, so it follows that the diagonal of a rational formal power series is D-finite. This fact was conjectured by Stanley[Stan], and incompletely proved in [Zei1] and [Ges1]. An elementary proof was given by Lipshitz[Lip1].

3.2 Mixed Continuous and Discrete Functions

Our discussion extends immediately to mixed continuous-discrete functions f : R^{n_{1}} × Z^{n_{2}} → C (in fact the continuous component may be a distribution, i.e. f(−, m) may be a distribution on R^{n_{1}} for every m in Z^{n_{2}}.)

Holonomic entities on R^{n_{1}} × Z^{n_{2}} are closed under addition, and if well defined, also under multiplication and convolution.

The same argument as above shows that if f : R^{n_{1}}×Z^{n_{2}} → C is holonomic then f(x_{1}, ..., x_{n_{1}−1}, b, m_{1}, ..., m_{n_{2}}) and ∫ fdx_{n_{1}} are holonomic in R^{n_{1}−1} × Z^{n_{2}} and f(x_{1}, ..., x_{n_{1}}, m_{1}, ..., m_{n_{2}−1}, a) and ∑_{m_{n_{2}}} f are holonomic as functions on R^{n_{1}} × Z^{n_{2}−1}.

3.3 Applications To Combinatorics

The general problem in enumerative combinatorics is that of counting the number of elements a(n) in a family A(n) of combinatorial objects, parametrized by the (single or multi) discrete variable n. For example the number of subsets of an n-element set is 2^{n}. Of course it is nice to get a simple explicit expression for a(n), but this is very rarely possible. But what is an explicit answer? ([Wil].) One narrow definition is that a(n + 1)/a(n) be a rational function in n. In other words, we know the answer explicitly if a(n) satisfies a first order recurrence equation with polynomial coefficients. One way of broadening our definition of ”explicit” is not to insist on the recurrence being first order. In other words, given a combinatorial sequence a(n), to which there is no simple explicit form, it is of interest to know whether it is P-finite, i.e. is a solution of a linear recurrence equation with polynomial coefficients (see Stanley’s paper [Stan]).

We saw in section 2.5 that the factorial and multinomial coefficients are holonomic in their variables. Many times in combinatorics , there is an explicit expression for a(n) in terms of a huge sum of summands that are obviously holonomic. The general scenario is that the the summand depends on, say, k parameters, and the sum contains k − 1 sigmas. It follows then from an iteration of propositions 3.4’ and/or 3.5’ that a(n) is holonomic in the surviving variable n. But being holonomic in a single variable is nothing but being P-finite! We have thus a powerful way of proving that sequences of combinatorial interest are P-finite.

One example is t_{k}(n), the number of standard Young tableaux with n cells and height ≤ k. For every fixed k , this is a sum of f_{λ} over the set

{(m_{1}, ...m_{k}); m_{1} + ... + m_{k} = n, m_{1} ≥ m_{2} ≥ ... ≥ m_{k} ≥ 0}.

f_{λ} is the number of standard Young tableaux of shape λ = (m_{1}, ..., m_{k}), and is given explicitly by the famous Young-Frobenius formula:

f_{λ} = [ ∏_{1≤i<j≤k}(m_{i} − m_{j} + j − i)]\frac{(m_{1} + ... + m_{k})!}{(m_{1} + k − 1)!...(m_{k})!}.

Since f_{λ} is obviously holonomic in its variables, it follows that for every fixed k, t_{k}(n) is P-finite. (Explicit expressions for t_{k}(n) for small values of k were given by A. Regev[Reg](k = 3), and D. Gouyou-Beauchamps[Gouy](k = 4, 5), and I. Goulden[Goul] developed an intersting symmetric functions method for handling these and related sums.)

Many other combinatorial applications were found by Gessel [Ges2],[Ges3]. L. Lipshitz[Lip1][Lip2] gave elegant elementary proofs of many of the results in this section for D-finite formal power series.

4 How To ”Write Down” Holonomic Functions And How To Verify Identities

4.1 Continuous Holonomic Functions

In order to actually work with holonomic functions we should have some way to ”write them down”. A holonomic function is completely determined by a left ideal J of A_{n}(C) that annihilates it, and by ”initial conditions”. We will now show that every holonomic function in n variables has n ”ordinary” operators P_{i}(D_{i}, x_{1}, ..., x_{n}) , i = 1, ..., n that annihilate it. This would follow from the following lemma, whose proof was kindly shown to me by Joseph Bernstein.

Lemma: Let L be a left ideal in A_{n}(C) such that A_{n}(C)/L is a holonomic A_{n}(C) module. For every n + 1 generators out of the 2n generators < x_{1}, ..., x_{n}, D_{1}, ..., D_{n} > of A_{n}(C) there is a non-zero member of L that only depends on these n+1 generators.

Proof: For the sake of definiteness let us take the n+1 generators to be x_{1}, ..., x_{n}, D_{1}. The proof in general is similar. Consider the mapping

ϕ : C < x_{1}, ..., x_{n}, D_{1} >→ A_{n}(C)/J,

given by ϕ(f) = f mod J. Let ∪F_{v}^{(0)} and ∪F_{v}/J be the natural filtrations of C < x_{1}, ..., x_{n}, D_{1} > and A_{n}(C)/J respectively. (Recall that F_{v} is the vector space of all members of A_{n}(C) of degree ≤ v .) Since A_{n}(C)/J is holonomic, we have that dim(F_{v}/J) is a polynomial in v of degree n for v >> 0. Obviously dim(F_{v}^{(0))} is a polynomial of degree n + 1 , so it follows that there is a v such that dim(F_{v}^{(0)}) > dim(F_{v}/J). The restriction of the linear map phi to the finite-dimensional vector space F_{v}^{(0)} is therefore a linear transformation from a higher dimensional vector space to a lower dimensional vector space, and its kernel must therefore be non-zero. But this kernel is precisely J ∩ C < x_{1}, ..., x_{n}, D_{1} >.

The above lemma shows that it is always possible to eliminate any n − 1 of the generators. The above proof is an existence proof. In section 5, I will present an algorithm that actually constructs the eliminated operators.

In order to write down a holonomic function f in ”holonomic notation”, we can give any set of generators of its corresponding ideal, with the appropriate initial conditions. However, in general it is not clear how many initial conditions are required to uniquely specify the function f. It is therefore necessary to introduce what I call ”canonical holonomic representation” as follows. we find the n ”ordinary” operators guaranteed by the above lemma

P_{i}(D_{i}, x_{1}, …, x_{n}), i = 1, …, n,  (4.1a)

of order α_{i} (in D_{i}), say, that annihilate f, and give the ”initial conditions”

D_{1}^{i_{1}}...D_{n}^{i_{n}}f(x_{0}), 0 ≤ i_{1} < α_{1}, ..., 0 ≤ i_{n} < α_{n}  (4.1b)

where x_{0} is any point that is not on the ”characteristic set” of the system (4.1a) (the characteristic set of a system is the set of common zeroes of the leading coefficients of the operators P_{i}.) The discussion below assumes, for the sake of simplicity, that the characteristic set of (4.1a) is finite. If it is not, we have to take some other operators out of the annihilation ideal I_{f}, until the characteristic set is finite, and then give the appropriate initial conditions for the enlarged set of equations. A good ”canonical holonomic representation”, that always works, is in terms of Grobner bases.

In the last section we proved that if f and g are holonomic then so are their sum and product. In order to manipulate concrete holonomic functions we should know how to find the ” canonical holonomic representation” (4.1) of f + g and fg given the canonical representations of f and g.

Suppose P (D_{1}, x_{1}, ..., x_{n})f = 0, Q(D_{1}, x_{1}, ..., x_{n})g = 0, where P has order α and Q has order β (in D_{1}). We can find an operator R = R(D_{1}, x_{1}, ..., x_{n}) of order α + β that annihilates f + g as follows. The equation P f = 0 can be used to express D_{1}^{α}f as a linear combination, with coefficients that are rational functions in x_{1}, ..., x_{n}, of f, , D_{1}f, , D_{1}^{2}f, , , ..., , , D_{1}^{α−1}f. By successively applying D_{1} , using the the product rule of differentiation, and replacing D_{1}^{α}f by the above mentioned linear combination, we can express D_{1}^{i}f, for 0 ≤ i ≤ α + β in terms of D_{1}^{i}f for 0 ≤ i < α. Similarly we can express D_{1}^{i}g for 0 ≤ i ≤ α + β in terms of D_{1}^{i}g for 0 ≤ i < β. It follows that we can express each of the α + β + 1 quantities D_{1}^{i}(f + g), 0 ≤ i ≤ α + β , as linear combinations of the α + β quantities D_{1}^{i}f, 0 ≤ i < α, and D_{1}^{j}g, 0 ≤ j < β. It follows that D^{i}(f + g), 0 ≤ i ≤ α + β are linearly dependent (over the rational functions in x_{1}, ..., x_{n}). The linear relation between them, after clearing denominators, is exactly a differential equation in D_{1}, satisfied by f + g. Now we repeat the process for D_{i} and get a system of the form (4.1) satisfied by f + g. To get the initial conditions we just use the information about f and g. To get the higher order initial conditions for f and g, out of those that are given, we iteratively use the differential equations and plug in the point where the initial conditions are taken.

Trivial Example: The holonomic functions (in the single variable x) f = e^{−x} and g = e^{−x^{2}} have the following canonical holonomic representations:

f = [D + 1; f(0) = 1] , g = [D + 2x; g(0) = 1] .

Now

(f + g) = f + g , D(f + g) = −f − 2xg , D^{2}(f + g) = f + (4x^{2} − 2)g .

Eliminating f and g yields the following differential equations for (f + g):

(−2x + 1)(f + g)′′ + (−4x^{2} + 3)(f + g)′ + (−4x^{2} + 2x + 2)(f + g) = 0.

Since (f +g) satisfies a second order differential equation, we need to know (f +g)(0) and (f +g)′(0). Of course (f +g)(0) = f(0)+g(0) = 2, but we don’t know (f +g)′(0) right away, since f′(0) and g′(0) are not part of the canonical representation of f and g respectively. However, we can find f′(0) and g′(0) from their respective differential equations:f′(0) = −f(0) = −1 and g′(0) = −2(0)g(0) = 0, so (f +g)′(0) = −1+0 = −1. It follows that the canonical holonomic representation of h := (f +g) is

h = [(−2x + 1)D^{2} + (−4x^{2} + 3)D + (−4x^{2} + 2x + 2); h(0) = 2, h′(0) = −1].

A similar process applies to fg. If P has order α and Q has order β, then it is possible to find a differential equation of order αβ satisfied by fg. Now we use Leibnitz’s rule to express the αβ + 1 quantities D_{1}^{i}(fg), 0 ≤ i ≤ αβ in terms of the αβ quantities (D_{1}^{i}f)(D_{1}^{j}g), for 0 ≤ i < α and 0 ≤ j < β.

The canonical representation (4.1) , with its accompanying initial conditions is not unique. Given such a representation, we can left multiply by any operator in D_{i}, get higher order equations, and add the appropriate number of initial conditions. It is therefore necessary to know when two holonomic functions, given in terms of their canonical representation (4.1), are in fact the same. If f and g are holonomic functions suspected to be the same, which are given in terms of canonical representations, we can use the above method to find a canonical representation of f − g. Then we find the appropriate initial values of f − g and check that they are all 0. Unfortunately, we actually have to find a concrete canonical representation, because otherwise we would not know the orders of the P_{i}, and the characteristic set of the system, both of which are needed to find how many initial conditions are required.

4.2 Discrete Holonomic Systems

Analogous considerations apply in the discrete realm. Every holonomic function f of the n discrete variables m_{1}, , , ..., , m_{n} can be specified by giving n ”ordinary” recurrence equations:

P_{i}(E_{i}, m_{1}, ..., m_{n})f = 0, i = 1, …, n  (4.2)

with appropriate initial conditions. If we write

P_{i} = ∑_{j=0}^{α_{i}} p_{j}^{(i)}E_{i}^{j}

where α_{i} is the order of P_{i} as a recurrence operator in E_{i}, then the initial conditions needed are

f(a_{1} + b_{1}, ..., a_{n} + b_{n}), 0 ≤ b_{i} < α_{i}, i = 1, ..., n.

where (a_{1}, ..., a_{n}) is any point in Z^{n}. In addition we need the values of f in the characteristic set of all common zeroes, on Z^{n}, of the leading coefficients p_{α_{i}}(m_{1}, ..., m_{n}). If this set is not finite, then we must add to the ”ordinary” equations in (4.2) some more ”partial” equations out of its annihilation ideal I_{f}, with the appropriate initial conditions.

4.3 Continuous-Discrete Functions

Everything applies equally well to mixed continuous-discrete functions on R^{N} × Z^{M}, with the obvious analogy.

4.4 How To Compute Indefinite Integrals And Sums

We know from Propositions 3.5 and 3.5’ that if f is holonomic in R^{N} × Z^{M}, then

F_{1} := ∫ fdx_{N} and F_{2} := ∑_{m_{M}} f,

are holonomic in R^{N−1} × Z^{M} and R^{N} × Z^{M−1} respectively. However, in order to compute with holonomic functions we would need a way to find canonical holonomic representations of F_{1} and F_{2} out of a canonical holonomic representation of f.

It can be seen that the ideal I_{F_{1}} in A_{N−1,M}(C) that annihilates F_{1} is nothing but the ideal obtained by setting D_{N} = 0 in the ”elimination ideal” I_{f} ∩ C < x_{1}, ..., x_{N−1}, D_{1}, ... >. Similarly I_{F_{2}} is obtained by setting E_{M} = 1 in the elimination ideal I_{f} ∩ C < ..., ..., ..., m_{1}, ..., m_{M−1} >.

In the next section we will see how to adapt Sylvester’s classical dialytic elimination method ([VDW],II)to the Weyl algebra. Galligo[Gal] showed how to adapt Buchberger’s powerful method of Grobner bases to the Weyl algebra which implies, in particular, efficient algorithms for elimination.

4.5 How To Verify Holonomic Functions Identities

We are now capable, at least in principle, of verifying any identity that involves a fixed number of sums, differences, products, integrals, and sigmas of holonomic functions in a fixed set of variables. Simply bring it to the form in which the right side is zero, and compute a canonical holonomic representation of the right. Finally compute the necessary initial conditions and verify that they are all zero.

4.6 Special Functions Identities

Many special functions, in particular all those in the famous ”tableau d’Askey”([Ask-Wil],[Lab]), are holonomic in all their variables and parameters. For example the Legendre polynomials P_{n}(x), when viewed as a function of n and x, is holonomic, and the Jacobi polynomials P_{n}^{(α,β)}(x) when viewed as a function of n, x, α, and β is a holonomic fucnction on N × R^{3}. In general, any function that is given by hypergeometric summation, whose entries and arguments are linear combinations of the variables and parameters, is holonomic. This follows from the fact the summand (w.r.t. to k, say) is clearly holonomic in all its variables and parameters, and it follows from proposition 3.5’ that the sum itself is a holonomic function.

It follows from the above comments that, at least in principle, we can verify any identity involving a finite number of sums, products, sigmas and integrals of these special functions. However, at present I am only able to implement successfully the verification of single sum terminating hypergeometric identities, which is carried out in sections 5 and 6 and in the appendix . I am afraid that it is a long way before the Askey-Gasper identity ([Ask-Gas], (2.30)) can be machine-proved. On the other hand Apery’s recurrence can be verified using the program. Of course, we still need the human to conjecture the identities, but once conjectured, the machine can save us the chore of proving them.

5.Sylvester’s Dialytic Elimination in the Weyl Algebra

¯5.0

In the previous section we saw that if P (n, k, E_{n}, E_{k}) and Q(n, k, E_{n}, E_{k}) are independent in the sense that they generate a holonomic system, then we can find an operator R(n, E_{n}, E_{k}), independent of k, in the ideal I generated by P and Q. The proof given there however was an existence proof (using the pigeonhole principle). In this section we will present an algorithm that inputs P (n, k, E_{n}, E_{k}) and Q(n, k, E_{n}, E_{k}) and outputs an operator R(n, E_{n}, E_{k}), with k missing , and two operators A and B such that R = AP + BQ. We will then show how this implies an effective algorithm for verifying any binomial coefficients identity (≡ terminating hypergeometric sum).

Given F (n, k) that is holonomic , proposition 3.4’ guarantees that

a(n) := ∑_{k} F (n, k)

is holonomic in the single variable n, wherever a(n) is defined (i.e. the sum converges). But being holonomic in the single variable n means that a(n) satisfies a linear recurrence equation with polynomial coefficients. Once we found an operator R(n, E_{n}, E_{k}) that annihilates F (n, k), i.e. RF ≡ 0, we can find an equation satisfied by a(n) as follows.

Write

R(n, E_{n}, E_{k}) = (1 − E_{k})^{α}R̃(n, E_{n}, E_{k}),

where α is maximal, i.e. R̃(n, E_{n}, I) ≠ 0.

We get

(1 − E_{k})^{α}{R̃(n, E_{n}, E_{k})F (n, k)} = 0.

But this means that

G(n, k) := R̃(n, E_{n}, E_{k})F (n, k)

is a polynomial in k for every fixed n. If G(n, k) were a non-zero polynomial in k, then ∑_{k} G(n, k) would have diverged, a contradiction, since the convergence of ∑_{k} F (n, k) implies the convergence of ∑_{k} G(n, k). So G(n, k) ≡ 0 and we have

R̃(n, E_{n}, E_{k})F (n, k) = 0.

Now write

R̃(n, E_{n}, E_{k}) = S(n, E_{n}) + (1 − E_{k})R̃′(n, E_{n}, E_{k})

where S(n, E_{n}) := R̃(n, E_{n}, I) (≠ 0).

We have

0 = ∑_{k} R̃(n, E_{n}, E_{k})F (n, k) = ∑_{k} S(n, E_{n})F (n, k) + ∑_{k}(I − E_{k})[R̃′(n, E_{n}, E_{k})F (n, k)].

The second term vanishes by telescoping, so this is equal to

S(n, E_{n}) ∑_{k} F (n, k) = S(n, E_{n})a(n).

We have just proved

Theorem:Let F : Z^{2} → C be such that

∑_{k} F (n, k)

converges for every n (this happens, in particular, if F (n) has finite support, like in all terminating hypergeometric series). If F (n, k) is a solution of a linear partial recurrence equation with polynomial coefficients of the form

R(n, E_{n}, E_{k})F (n, k) ≡ 0

with k missing from R, then

a(n) := ∑_{k} F (n, k)

satisfies the following ordinary linear recurrence equation with polynomial coefficients

S(n, E_{n})a(n) ≡ 0, , ,

where S(n, E_{n}) = R̃(n, E_{n}, I) and R̃ is obtained from R by dividing by the highest possible power of (I − E_{k}).

5.1 Sylvester’s Dialytic Elimination in Commutative Algebra

We will now review Sylvester’s classic dialytic elimination method [Syl] (see also [vdW]). In the next section we will describe how to modify it to the context of the (non-commutative) algebra of linear partial recurrence operators with polynomial coefficients. Sylvester’s method basically consists of reducing algebraic elimination to linear elimination, so we will start with linear elimination.

5.1.1 Linear (Gaussian) Elimination: Let R be any commutative ring, and let x be an indeterminate that does not necessarily commute with the elements of R. Consider the two linear affine forms P = ax + b and Q = cx + d. In order to eliminate x, we multiply P by c, Q by a, and substruct: cP − aQ = cb − ad. So the ring element cb − ad is in the ”ideal generated by P and Q”. All we needed was that a and c commute with each other. The procedure is still valid if b and d do not commute between themselves or with a and c.

More generally, consider the n affine-linear forms in the indeterminates x_{1}, ..., x_{n−1}:

P_{j} = ∑_{i=1}^{n−1} a_{i,j}x_{j} + b_{j}, j = 1, …, n ,

and consider the matrix of coefficients

( \begin{matrix}a_{1,1} & … & a_{n−1,1} & b_{1} \\ … & … & … & … \\ a_{1,n} & … & a_{n−1,b} & b_{n} \end{matrix} ) .

Let Δ_{i} be the cofactor of b_{i} in the above matrix A. Then obviously

Δ_{1}P_{1} + ...Δ_{n}P_{n} = det(A) .

and det(A), independent of the x_{i}, is in the ”ideal generated by P_{1}, ..., P_{n}.” It is important for the sequel to note that all we need is that the a_{i,j} be pairwise commutative, but the b_{i} neither have to commute with the a_{i,j}, nor among themselves. Of course , in that case we must be careful and evaluate the determinant with respect to the last column.

5.1.2. Sylvester’s Dialytic Elimination: Following Sylvester[Syl], we will first do an example. Let

P = ax^{2} + bx + c; Q = a′x^{2} + b′x + c.

Form xP and xQ, and consider P, Q, xP, xQ as affine-linear forms in x, x^{2}, x^{3}:

( \begin{matrix}c & b & a & 0 \\ 0 & c & b & a \\ c′ & b‘ & a′ & 0 \\ 0 & c′ & b′ & a′ \end{matrix} ) ( \begin{matrix}1 \\ x \\ x^{2} \\ x^{3} \end{matrix} ) = ( \begin{matrix}0 \\ 0 \\ 0 \\ 0 \end{matrix} )  (5.1)

Let A be the matrix on the right side of (5.1) and let Δ_{1}, Δ_{2},Δ_{3},and Δ_{4}, be the cofactors of the entries of the last column in the above matrix, then

(Δ_{1} + Δ_{2}x)P + (Δ_{3} + Δ_{4}x)Q = det(A).

It follows that det(A) belongs to the ideal in R[x] generated by the two polynomials P (x) and Q(x), and is independent of x, i.e. an element of the ring of coefficients R. This is called the resultant of P and Q. If det(A) is equal to zero then it means that P and Q have a common factor ([vdW] I,p.83). So if P and Q are ”independent” in the sense that they don’t have a common factor, (or, equivalently, a common zero), then the resultant must be non-zero.

More generally, to eliminate x out of

P = a_{m}x^{m} + a_{m−1}x^{m−1} + ... + a_{0}
Q = b_{n}x^{n} + b_{n−1}x^{n−1} + ... + b_{0},

we form P, xP, ..., x^{n−1}P, and Q, xQ, ..., x^{m−1}Q and get m + n linear forms in the indeterminates 1, x, …, x^{m+n−1}:

(P, …, x^{n−1}P, Q, …, x^{m−1}Q)^{T} = A(x^{m+n−1}, …, 1)^{T}  (5.2)

where A is the appropriate matrix. Let Δ_{i} (i = 1, ..., n) and Δ′_{j}(j = 1, ..., m) be the m+n cofactors of the entries of the last column of A, then we have

(Δ_{1} + Δ_{2}x + … + Δ_{n}x^{n−1})P + (Δ′_{1} + Δ′_{2}x + ...Δ′_{m}x^{m−1})Q = det(A).

det(A) is called the resultant since it is an expression in the coefficients of P and Q that determines whether they have a common factor: det(A) is zero if and only if P and Q have a common factor ([vdW],I, section 27, p.83). det(A) belongs to ”the ideal generated by P and Q”, and is independent of x.

5.2 Dialytic Elimination In The Weyl Algebra: Two Simple Examples

Given two partial linear recurrence operators with polynomial coefficients P (N, K, n, k) and Q(N, K, n, k) we would like to ”eliminate k”, i.e. find two operators, A and B, such that R := AP + BQ is independent of k, and, of course, non-zero. Trying to emulate Sylvester’s dialytic method, the first thought that comes to mind is to multiply P and Q by powers of k. However, we cannot take the determinant, since the ”ring of coefficients” is C < N, K, n >, which is not commutative, so we don’t get rid of k. Our twist on Sylvester’s method is to multiply P and Q by various monomials n^{i}k^{j}, get many more ”affine-linear forms” in the ”indeterminates” n^{α}k^{β}, (β > 0) with coefficients in the commutative ring C[N, K], with the constant part, i.e. the coefficient of ”1”, belonging to C < N, K, n >. Only now can we take the determinant of the resulting system, and get the eliminated operator R. I will first present two examples.

Let us define the shift operators N and K by Nf(n, k) := f(n + 1, k), Kf(n, k) := f(n, k + 1). A partial recurrence operator with polynomial coefficientsis any polynomial in the four indeterminates n,k,N,K. They satisfy the commutation relations: Nn = (n+1)N, Kk = (k +1)K, Nk = kN, nK = Kn, nk = kn.

By using the commutation relations it is easy to see that every operator can be written in the form

∑_{r,s} p_{r,s}(N, K)n^{r}k^{s}.  (1)

Trivial Example : F (n, k) := n!/k!(n − k)!, find R(n, N, K) such that RF = 0.

Solution: NF/F = (n+1)/(n−k+1), KF/F = (n−k)/(k+1), so P = (n−k+1)N −(n+1), Q = (k + 1)K − (n − k). Converting to the format of (1) we have:

{P = −N(k) \atop Q = (K + 1)(k)} {+(Nn − n − 1)(1) \atop −n(1)}.

In this case the only indeterminates are ”k” and ”1”, and we already have as many equations as unknowns, so we don’t have to multiply by any monomials n^{r}k^{s}. The matrix of coefficients is

({−N \atop K + 1} {Nn − n − 1 \atop −n})

whose determinant is R(n, N, K) := (−N)(−n) − (K + 1)(Nn − n − 1)= (n + 1)(−NK + K + 1). Thus if n ≠ −1, F (n, k) = \binom{n}{k} satisfies

(−NK + K + 1)\binom{n}{k} = 0,

and we have rediscovered Pascal’s triangle. In particular R(n, N, I) = (n + 1)(−N + 2I), and it follows that if

a(n) := ∑_{k} F (n, k),

then (n + 1)(−N + 2I)a ≡ 0 which in everyday language means (n + 1)(−a(n + 1) + 2a(n)) = 0, whose solution is C2^{n}, for some constant C, which is easily seen to be 1 by plugging in n=0. We have just proved the deep binomial coefficients identity

∑_{k} \binom{n}{k} = 2^{n}.

A Simple But Not Quiet So Trivial Example Let

F (n, k) = ({n \atop k})({n + k \atop k}) = \frac{(n + k)!}{k!^{2}(n − k)!}

Here NF/F = (n + k + 1)/(n − k + 1), NKF/F = (n + k + 1)(n + k + 2)/(k + 1)^{2}. The relevant operators, in standard form are

P = (NK − 1)(k^{2}) − 2(nk) − 3(k) − (n + 1)(n + 2),
Q = −(N + 1)(k) + (N − 1)n − 1.

Now we have four indeterminates : k^{2}, k, nk, and 1, where the indeterminate ”1” is allowed to have coefficients in C < n, N, K >. However, we only have two equations. Is it possible to get more equations without introducing more indeterminates? Usually not, because creating more equations usually brings in more indeterminates and all you can do is reduce the excess ”#indeterminates − #equations, until you get 0. However, whenever one of the operators is linear in k and n, like Q is in this example, then we can get away without increasing the number of indeterminates.

Now multiply Q by n and k respectively, convert to the format (1), and you get

kQ = −(N + 1)k^{2} + (N − 1)nk − k
nQ = −(N + 1)nk + Nk + N(n^{2} − n) − (n^{2} + n)

So

( \begin{matrix}P \\ Q \\ kQ \\ nQ \end{matrix} ) = ( \begin{matrix}NK − 1 & −2 & −3 & −(n + 1)(n + 2) \\ 0 & 0 & −(N + 1) & (N − 1)n − 1 \\ −(N + 1) & N − 1 & −1 & 0 \\ 0 & −(N + 1) & N & N(n^{2} − n) − (n^{2} + n) \end{matrix} ) ( \begin{matrix}k^{2} \\ nk \\ k \\ 1 \end{matrix} )

Now take the determinant of this, and get R(n, N, K). Plugging in K = I gives a linear recurrence operator that annihilates

a(n) := ∑_{k} F (n, k) = ∑_{k}({n \atop k})({n + k \atop k}).

It turns out to be

n^{2}N^{4} + (−5n^{2} − 10n − 3)N^{3} + (2 − 5n^{2})N^{2} + (n^{2} + 2n + 1)N.

This operator is not the minimal operator annihilating a(n), so this method only produces an operator annihilating ∑_{k} F (n, k). However , a minimal operator can be found empirically and then one can use the Euclidean algorithm (adapted to linear ordinary recurrence operators with polynomial coefficients) to show that the conjectured operator is a right factor of the operator that the method outputed.

5.3 Sylvester Dialytic Elimination In The Weyl Algebra: The General Method

Consider two operators P (N, K, n, k) and Q(N, K, n, k) of degrees b_{1} and b_{2} respectively in k, and maximal degree c in n, and let’s write them in the standard form, and with the k-free part set apart at the end:

P = ∑_{i=1}^{b_{1}} ∑_{j=0}^{c} p_{i,j}(N, K)k^{i}n^{j} + R_{P} (n, N, K),

and

Q = ∑_{i=1}^{b_{2}} ∑_{j=0}^{c} q_{i,j}(N, K)k^{i}n^{j} + R_{Q}(n, N, K).

We consider all operators as polynomials in the variables n and k with coefficients that belong to the commutative ring C[N, K]. The coefficients mutually commute and the only non-commutativity that arises is from the interaction between the coefficients and the variables. As in the commutative case, let us form

P, kP, …, k^{b_{2}−1}P,

and

Q, kQ, …, k^{b_{1}−1}Q.

We must bring them to the standard form as above, using k^{i}N^{α}K^{β} = N^{α}K^{β}(k − β)^{i}.

We now consider these as linear forms in the ”indeterminates”

k^{α}n^{β}1 ≤ α ≤ b_{1} + b_{2} − 1, 0 ≤ β ≤ c.

Alas, now we have far too few linear forms, or equivalently, far too many indeterminates. Namely, we have b_{1} + b_{2} linear forms and (b_{1} + b_{2} − 1)(c + 1) indeterminates, and the parameter ”#linear forms − #indeterminates” is 1 − c(b_{1} + b_{2} − 1). We need this parameter to be 1.

To set the balance right, we will multiply by successive powers of n. At the r^{th} stage, when we multiply the b_{1} + b_{2} linear forms

{P, kP, ..., k^{b_{2}−1}P, Q, ..., k^{b_{1}−1}Q, }  (5.3)

by the power n^{r} we gain b_{1} + b_{2} new linear forms, but are only burdened with b_{1} + b_{2} − 1 new indeterminates: k^{α}n^{r+c}, for α = 1, ..., b_{1} + b_{2} − 1. It follows that the budgetary deficit ”#indeterminates − #linear forms” gets reduced by one for every power of n by which we multiply. So if we multiply (5.3) by n^{r} for r = 1, ..., c(b_{1} + b_{2} − 1), we will finally have one more linear form than we have indeterminates.

Let A_{f}(i, j) be the ”k-free” part of k^{i}n^{j}P (i = 0, ..., b_{2} − 1, j = 0, ..., c(b_{1} + b_{2} − 1)) , and B_{f}(i, j) be the ”k-free” part of k^{i}n^{j}Q (i = 0, ..., b_{1}−1, j = 0, ...c(b_{1}+b_{2}−1)) . Also let A(i, j) = k^{i}n^{j}P −A_{f}(i, j) and B(i, j) = k^{i}n^{j}Q − B_{f}(i, j). By linear algebra, A(i, j) ∪ B(i, j) are linearly dependent, so there exist a_{i,j}(N, K), and b_{i,j}(N, K), not all zero such that

∑ a_{i,j}(N, K)A(i, j) + ∑ b_{i,j}(N, K)B(i, j) = 0.  (5.4)

It follows that

R := (∑ a_{i,j}(N, K)k^{i}n^{j})P +(∑ b_{i,j}(N, K)k^{i}n^{j})Q = ∑ a_{i,j}(N, K)A_{f}(i, j)+∑ b_{i,j}(N, K)B_{f}(i, j)

is independent of k, i.e. R = R(n, N, K). Obviously R is an operator that belongs to the left ideal generated by P and Q, and so if F (n, k) is a solution of the two linear partial recurrence equations P (n, k, N, K)F ≡ 0 and Q(n, k, N, K)F ≡ 0, then F automatically satisfies R(n, N, K)F ≡ 0. In other words we have succeeded in eliminating k from P and Q.

We have to make sure that R is not the zero operator. Of course if P and Q are dependent, for example if they are both left multiples of the same operator: P = P′A , Q = Q′A, then R will be zero. The technical condition of ”independence” is exactly that P, Q generate a holonomic system, i.e. that the ideal I generated by P and Q is such that A_{2}(C)/I is a holonomic A_{2}(C) module.

Suppose that R was the zero operator. This means that we can find two operators A(k, n, K, N) of degree < b_{2} in k and B(k, n, K, N) of degree < b_{1} in k such that AP = BQ. Now taking the ”principal symbols” as operators in A_{2}(C) will show that the images of P and Q in the commutative graded ring are dependent, which is impossible if P and Q generate a holonomic system.

To find the a_{i,j} and b_{i,j} of (5.4) we can take the cofactors of the last column in the determinant formed by the matrix of coefficients of the k^{i}n^{j}P, k^{i}n^{j}Q expressed as linear combinations of the monomials n^{α}k^{β}, for β > 0 and the ”monomial” k^{0}. If all this cofactors are zero, than we would have to take the cofactors of an appropriate non-singular subdeterminant.

5.4 Implementation

In the preceding discussion we worked in terms of the k-degrees in k, b_{1}, and b_{2}, of P and Q respectively, and the highest n-degree , c. Thus it applies to operators whose ”supports”:

{(i, j); p_{i,j} ≠ 0}, {(i, j); q_{i,j} ≠ 0},

are included in the rectangles

{(i, j); 0 ≤ i ≤ b_{1}, 0 ≤ j ≤ c},

and

{(i, j); , 0 ≤ i ≤ b_{2}0 ≤ j ≤ c}, , .

It turns out that most operators that occur ”in nature” have their ”support” in the shape of a triangle:

{(i, j); 0 ≤ i, 0 ≤ j, i + j ≤ d},

so it is a waste to consider the smallest rectangle containing their support, since half of it will be zeroes. In short, it is more natural to consider the total degree in k and n. This reduces the size of the determinant considerably. In section 6 I will give an explicit algorithm that uses the total degree.

5.5 Elimination In General

A similar method works for the continuous case. Just bring all operators to the form

P = ∑_{i=1}^{b_{1}} ∑_{j=0}^{c} p_{i,j}(D_{x}, D_{y})y^{i}x^{j} + R_{P} (n, N, K),

multiply from the left by monomials y^{α}x^{β}, and convert back to that form using

xD_{x}^{i}D_{y}^{j} = D_{x}^{i}D_{y}^{j}x − iD_{x}^{i−1}D_{y}^{j}

and

yD_{x}^{i}D_{y}^{j} = D_{x}^{i}D_{y}^{j}y − jD_{x}^{i}D_{y}^{j−1}.

The q-analog is just as straightforward. For the discrete q-analog you replace n by q^{n}, k by q^{k}, and the commutation rules nN = N(n − 1), kK = K(k − 1), by q^{n}N = Nq^{n}/q, q^{k}K = Kq^{k}/q. For the continuous q-analog, featuring the q-dilation operators Q_{x}f(x, y) = f(qx, y) and Q_{y}f(x, y) = f(x, qy), we have the commutation rules xQ_{x} = Q_{x}x/q and yQ_{y} = Q_{y}y/q.

Everything extends to the elimination of several variables from a general annihilation ideal, as in the classical case ([vdW]), but in this cases Buchberger’s method of Grobner bases, as adapted by Galligo[Gal], is far superior.

6. An Algorithm For Proving Binomial Coefficients Identities

When someone finds a new identity, there are not many people who get excited about it any more, except the discoverer. (Donald E.Knuth [Knu],p.53)

6.1 The Algorithm

We saw that, in principle, it is possible to verify any (holonomic) special function identity involving sums, integrals, and products. However, as Richard Askey said so aptly in ([Ask2]): ”Many things that can be done in theory cannot be done in practice”. In this section I will show that at least some things can be done in practice, by giving an explicit algorithm for proving, or refuting, binomial coefficients identities (≡ terminating hypergeometric identities) of the form

∑_{k} F (n, k) = rhs(n).  (6.1)

A MAPLE implementation of this algorithm is given in the appendix.

In (6.1) it is assumed that F (n, k) has the form

F (n, k) = \frac{∏_{i=1}^{A}(a_{i}n + a′_{i}k + a′′_{i})!}{∏_{i=1}^{B}(b_{i}n + b′_{i}k + b′′_{i})!}z^{k}  (6.2a)

where the a_{i}, a′_{i}, b_{i}, b′_{i} have to be constant, specific, (positive or negative) integers, but z, a′′_{i}, and b′′_{i} can be any complex numbers or parameters, and x! means Γ(x + 1).

The right side of (6.1), rhs(n), may be given explicitly, in the form

rhs(n) = C \frac{∏_{i=1}^{A}(ā_{i}n + ā′_{i})!}{∏_{i=1}^{B}(b̄_{i}n + b̄′_{i})!}x̄^{n}  (6.2b)

where ā_{i} and b̄_{i} are specific (positive or negative) integers, and ā′_{i}, b̄′_{i} ,C, and x̄, are complex numbers or parameters. Another possibility is that the right side is given implicitly in terms of a minimal (ordinary) linear recurrence operator with polynomial coefficients, conj(n, N) that annihilates rhs(n), together with the appropriate initial conditions.

The algorithm consists of 10 steps.

Step I: (a) Find operators p(k, n, N, K) and q(k, n, N, K) that annihilate F (n, k). [ In a more general situation p and q may be given from the outset, in which case you go to step 2 directly. ]

You do this by computing the rational functions F (n + 1, k)/F (n, k) and F (n, k + 1)/F (n, k):

\frac{F (n + 1, k)}{F (n, k)} = \frac{P (n, k)}{Q(n, k)} , \frac{F (n, k + 1)}{F (n, k)} = \frac{P′(n, k)}{Q′(n, k)}.  (6.4)

where P, Q, P′, Q′ are polynomials in n and k and, of course, gcd(P, Q) = 1 and gcd(P′, Q′) = 1. The operators p and q are given by

p(k, n, N, K) := NQ(n − 1, k) − P (n, k), , , , q(k, n, N, K) := KQ′(n, k − 1) − P′(n, k), , , .

[p and q both annihilate F (n, k) ].

(b) If the right side of (6.1) is given explicitly, find the (minimal) linear recurrence operator with polynomial coefficients that annihilates rhs(n). You do this by computing rhs(n + 1)/rhs(n) = P̄(n)/Q̄(n), say, and set

conj(n, N) := NQ̄(n − 1) − P̄(n), , , , .

Step II : Find the degree, in k, of p(k, n, N, K) and q(k, n, N, K), say α and β respectively. [It is readily seen that because the form (6.2a) of F (n, k) the (generic) degree of p in n equals to the degree in k, which equals to the total degree in (n, k), and similarly for q. We will view all operators as polynomials in n and k with coefficients that belong to C[N, K] . ]

Step III : Let l1 and l2 be the following sets of pairs of integers (i, j):

l1 := {(i, j); 0 ≤ i ≤ β − 1j ≥ 0, i + j ≤ α(β − 1)}

l2 := {(i, j); 0 ≤ i ≤ α − 1, j ≥ 0, i + j ≤ (α − 1)β}.

For each pair (i, j) in l1 form the operator k^{i}n^{j}p, convert it to standard form, and expand in terms of the monomials k^{i′}n^{j′}. Similarly, for q and l2. So let

k^{i}n^{j}p = ∑ A_{i′j′}^{(i,j)}(K, N)k^{i′}n^{j′}(i, j) ∈ l1  (6.4a)

k^{i}n^{j}q = ∑ B_{i′j′}^{(i,j)}(K, N)k^{i′}n^{j′}(i, j) ∈ l2.  (6.4b)

[ A simple subroutine left-multiplies any operator by a monomial and converts it to standard form, using k^{i}n^{j}K^{a}N^{b} = K^{a}N^{b}(k − a)^{i}(n − b)^{j} .]

Step IV : [ We now consider (6.4) as affine-linear forms in the ”indeterminates” k^{i′}n^{j′}. It is readily seen that the monomials that feature in (6.4) are those belonging to the following set

l3 := {(i′, j′); 1 ≤ i′ ≤ α + β − 1, j′ ≥ 0, i′ + j′ ≤ αβ}.

with coefficients that are in C[N, K] . The ”affine part” (i.e. ”the coefficient of k^{0}”) is in C < N, K, n >. Note that |l1| + |l2| = |l3| + 1 = (2αβ − α − β + 2)(α + β − 1)/2 + 1 . ]

Let

SIZE := (2αβ − α − β + 2)(α + β − 1)/2 + 1.

Form the SIZE × SIZE matrix M as follows. First convert the sets l1, l2, l3, into lists:

l1 := l1[r], , r = 1, ..., |l1|, l2 := l2[r], r = 1, ..., |l2|, l3 := l3[r], r = 1, ..., |l3|.

Now, for 1 ≤ r ≤ |l1|, and 1 ≤ s ≤ SIZE − 1,set

M[r, s] := A_{i′j′}^{(i,j)}(N, K)(of (6.4a)),

where (i, j) := l1[r] and (i′, j′) := l3[s].

For |l1| + 1 ≤ r ≤ SIZE, and 1 ≤ s ≤ SIZE − 1,set

M[r, s] := B_{i′j′}^{(i,j)}(N, K), (of (6.4b)),

where (i, j) := l2[r − |l1|] and (i′, j′) := l3[s].

Finally, for 1 ≤ r ≤ |l1|,

M[r, SIZE] := ∑_{j′} A_{0j′}^{(i,j)}(N, K)(of (6.4a)),

where (i, j) := l1[r]. In other words M[r, SIZE] is the coefficient of k^{0} in k^{i}n^{j}p. Similarly, for |l1| ≤ r ≤ SIZE,

M[r, SIZE] := ∑_{j′} B_{0j′}^{(i,j)}(N, K)(of (6.4b)))

where (i, j) := l2[r − |l1|]. In other words M[r, SIZE] is the coefficient of k^{0} in k^{i}n^{j}q.

Step V: Let

R(n, N, K) := determinant(M).

If R is identically zero (a very rare event), then take an appropriate subdeterminant, making sure that the last column gets chosen. [ More precisely, find the rank of the matrix A without the last column, pick that many linearly independent rows and and one less column and adjoin the (relevant part of the) last column .] [R(n, N, K) annihilates F (n, k)].

Step VI: Find R̄(n, N, K) such that

R(n, N, K) = (1 − K)^{g}R̄(n, N, K),

where g is as big as possible. [In section 5 we showed that R̄(n, N, K) also annihilates F (n, k)].

Step VII: Let S(n, N) be the ordinary linear recurrence operator with polynomial coefficients (in n) obtained by substituting K = 1 in R̄(n, N, K):

S(n, N) := R̄(n, N, 1).

[S(n, N) annihilates a(n) := ∑_{k} F (n, k) . ]

Step VIII: Use the Euclidean algorithm (adapted to the algebra of ordinary linear recurrence operators) to find operators T (n, N) and rem(n, N) such that

s(n, N) = T (n, N)conj(n, N) + rem(n, N) ,

where the degree of rem(n, N) in N (i.e. its order as a recurrence operator) is smaller than that of conj(n, N).

Step IX : If rem(n, N) ≠ 0 then the identity is false, i.e. the left side is not annihilated by conj(n, N), while the left side is (by definition). This follows from the fact that conj(n, N) was taken to be the lowest order (non-zero) operator that annihilates the right side. If the left side would have been equal to the right side, then rem(n, N) would also annihilate the right side, contradicting the minimality of conj(n, N).

If rem(n, N) ≡ 0, then the identity is true provided it is true at some initial points, which are determined as follows. We have

T (n, N)[conj(n, N)a(n)] ≡ 0.

Let the order (alias degree in N) of T (n, N) be t, and let the leading coefficient, of N^{t}, be mekadem(n). Then it is obvious that conj(n, N)a(n) ≡ 0 provided it is 0 at n = 0, ..., t − 1, and at the (the usually empty) set of positive integer zeroes of mekadem(n).

Step X : Now we know that

a(n) := ∑_{k} F (n, k)

and rhs(n) are both solutions of the linear recurrence equation

conj(n, N)s(n) ≡ 0.

In order to infer that a(n) and rhs(n) are identically equal we have to check that they match at the first r values n = 0, ..., r − 1, where r is the order of rhs(n, N), and at the ”characteristic set” of positive integer zeroes of the leading coefficient of rhs(n, N).

6.2 Discussion

The above algorithm is capable of proving any identity of the form (6.1), with the summand of the form (6.2a). This includes all terminating hypergeometric identities in Bailey[Bai], as well as those of Gosper’s[Gos] list, and almost all those of Gessel and Stanton[Ges-Sta]. A MAPLE program implementing the algorithm is given in the appendix. This program requires a lot of memory. For example, on my AT&T 3B1 PC, I was only able to do those identities for which SIZE, the size of the matrix whose determinant we take, is smaller than 13. In other words, my computer can handle the following values for (α, β) : {(α, 1);α ≤ 5}, {(1, β); β ≤ 5} and (2, 2). A sampling of the identities proved by my computer, using the program, is given in the appendix. For example, my computer can prove the Saalschutz identity, for which α = 2, β = 2, but it ran out of memory for Dixon’s identity, for which α = 2, β = 3, and hence SIZE = 19. Among the identities in [Ges-Sta], my program proved (1.1) (originally due to Andrews[And2]), (1.5), (5.24), and (5.25). I am sure that with a bigger computer, and even more importantly, with a better programmer, the algorithm would be able to handle much deeper identities. Incidentally, one may use the parameter SIZE as a measure of the depth of an identity.

In any case, the main virtue of the above algorithm is theoretical, because it tells you that you are guaranteed to find an operator R(n, N, K) annihilating F (n, k), and from its expression as a determinant, it is very easy to derive a priory bounds for the degrees in K,N, and n. Now that we know that there is such an operator, we can easily try out a generic form:

R(n, N, K) := ∑ a_{i,j,l}n^{i}N^{j}K^{l}

require that R(n, N, K)F (n, k) ≡ 0, plug in specific values of (n, k) and solve for the coefficients a_{i,j,l}. Having conjectured such an operator we can easily verify that indeed R(n, N, K) annihilates F (n, k) and obtain S(n, N) as before. Alternatively, we can use the method of creative telescoping.

6.3 The Method Of Creative Telescoping

The term creative telescoping was coined by van der Poorten in his charming account [vdP] of Apery’s proof of the irrationality of ζ(3). The method consists in ”creating” a certain multiple of F (n, k) by a rational function P (n, k)/Q(n, k):

G(n, k) := \frac{P (n, k)}{Q(n, k)}F (n, k)  (6.5)

with the property that

G(n, k + 1) − G(n, k)(= (K − 1)G(n, k)) = conj(n, N)F (n, k).  (6.6)

which immediately implies that a(n) := ∑_{k} F (n, k) satisfies conj(n, N)a(n) ≡ 0.

Perlstadt[Per] used creative telescoping and MACSYMA to find recurrences for sums of powers (up to the 6^{th}) of binomial coefficients.

In order to implement creative telescoping we write P (n, k) and Q(n, k) in generic form (guessing their degrees) , plug into (6.6), divide by F (n, k), clear denominators, and compare coefficients, getting linear equations that can be solved.

Can you always use creative telescoping? The algorithm of section 6.1 guarantees that we have an operator R̄(n, N, K) that annihilates F (n, k), and S(n, N) := R̄(n, N, 1) is non-zero. Then we have

R̄(n, N, K) − S(n, N) = (1 − K)L(n, N, K).

It follows that

S(n, N)F (n, k) = (1 − K)[L(n, N, K)F (n, k)].

We saw above that if the identity is true then S(n, N) must be a left multiple of conj(n, N). Taking

G(n, k) := L(n, N, K)F (n, k),

we see that (6.6) is always true, but sometimes with conj(n, N) replaced by some left multiple of it. So we now have an explanation why creative telescoping works. First try to find G(n, k) in (6.6) with the conjectured operator conj(n, N), if you fail, try to do it over with some left multiple of it.

Note added at the revised revision: The ideas in this section were extended to a much faster algorithm([Zei2],[Zei3]). This lead to the notions of ”WZ pair” and ”Rational function certification” ([W-Z1], [W-Z2], see also [Cip]).

6.4 Continuous and Discrete-Continuous Analogues

The algorithm of (6.1) can be easily adapted to the continuous case of eliminating y out of p(x, y, D_{x}, D_{y}) and q(x, y, D_{x}, D_{y}) , and proving identities of the form

∫_{−∞}^{∞} F (x, y)dy = a(x),

where F (x, y) has the form

F (x, y) = \frac{A(x, y)}{B(x, y)}exp(C(x, y)/D(x, y)),

for some polynomials A, B, C, D in x and y, and a(x) is either given explicitly, or in terms of a linear differential equation with polynomial coefficients that it satisfies.

We can also easily adapt the algorithm to handle discrete-continuous identities like

∑_{n=0}^{∞} P_{n}(x)t^{n} = 1/\sqrt{(1 − 2xt + t2)}.

Note added at the revised version: the ideas in this section were extended to a much faster algorithm in [A-Z].

6.5 The q-Analogue

The algorithm of section 6.1 can be easily adapted to the proof of any q-binomial identity of the form

∑_{k} F (n, k) = a(n)

where F (n, k) has the form

F (n, k) = \frac{∏_{i=1}^{A}(q^{a′′_{i}}; q)_{a_{i}n+a′_{i}k}}{∏_{i=1}^{B}(q^{b′′_{i}}; q)_{b_{i}n+b′_{i}k}} z^{k}

where the a_{i}, a′_{i}, b_{i}, b′_{i} have to be constant, specific, (positive or negative) integers, but z, a′′_{i}, and b′′_{i} can be any complex numbers or parameters. (Note added at the revised version: a much faster algorithm is given in [Zei4]).

More generally, the whole theory of holonomic systems can be q-analogized, where instead of C < E_{1}, ..., E_{n}, m_{1}, ..., m_{n} > we have C < E_{1}, ..., E_{n}, q^{m_{1}}, ..., q^{m_{n}} >, and the commutation rule is E_{i}(q^{m_{i}}) = q(q^{m_{i}}E_{i}). Also, instead of C < D_{1}, ..., D_{n}, x_{1}, ..., x_{n} >, we have C < Q_{1}, ..., Q_{n}, x_{1}, ..., x_{n} >, where Q_{i} are the q-dilation operators: Q_{i}f(x_{1}, ..., x_{i}, ..., x_{n}) = f(x_{1}, ..., qx_{i}, ..., x_{n}).

APPENDIX: A MAPLE IMPLEMENTATION OF THE ALGORITHM IN SECTION 6.1 THAT PROVES BINOMIAL COEFFICIENTS IDENTITIES

In a world in which the price of calculation continues to decrease rapidly, but the price of theorem proving continues to hold steady or increase, elementary economics indicates that we ought to spend a larger and larger fraction of our time on calculation. (John W.Tukey, The American Statistician, 40(1986)p.74)

A.1 The Program

The following is a MAPLE program that implements the algorithm in section 6.1. In order to use it, either type it into a file or get it from me via electronic mail. Call the program file verifier (or any name of your choice). Then in the UNIX operating system, that has MAPLE, you type the following on your terminal:

maple < verifier

The input given in the program below is for proving the following deep binomial coefficients identity

∑_{k} 1/(k!(n − k)!) = 2^{n}/n!, , .

In order to use this program for other identities

∑_{k} F (n, k) = a(n), , ,

you must first bring it to the form (6.2), i.e. make sure that the left has no terms depending on n only, and express everything in terms of factorials. Then you replace the 1/(k!(n − k)!) in the INPUT below by F (n, k) and the 2^{n}/n! by a(n). A Sampling of inputs for which this program ran successfully (i.e. did not run out of memory) on my AT&T 3B1 PC (that has 2 Meg of RAM, 67 Meg of hard disk, and 60free disk space), is given in section A.2.

If the right side a(n) is given implicitly, in terms of a recurrence, then you should change ”expli:=1” to ”expli:=0”, erase the line with rhs, and replace it by ” conj:= ...”, where the right side is the conjectured operator conj(n, N). See examples 1,3,6-10, in section A.2.

#THE IDENTITY VERIFIER #This program proves binomial coefficients identities(alias #terminating hypergeometric summation) of the #form described in section 6.1. The input consists of

#F(n,k), the summand of the l.h.s, and the conjectured r.h.s, #both of which should be quotients of products of factorials #of the form (an+bk+c)! for some integers a and b,where b is not 0. # The output consists of a human readable proof or refutation of the #identity

#THIS IS THE INPUT expli:=1: summand:=1/(k!*(n-k)!): rhs:=2n/n!: #END OF INPUT

if expli=1 then rat:=subs(n=n+1,rhs)/rhs: rat:=expand(rat): conj:=denom(rat)*N-numer(rat): fi: rat:=subs(n=n+1,summand)/summand: rat:=expand(rat): p1:=subs(n=n-1,denom(rat))*N-numer(rat): rat:=subs(k=k+1,summand)/summand: rat:=expand(rat): q1:=subs(k=k-1,denom(rat))*K-numer(rat): dp1:=degree(p1,k): dq1:=degree(q1,k): p1:=expand(p1): q1:=expand(q1):

#to generate the lists of exponents mapa, maq #by which p1 and q1 should be multiplied respectively, and khezka #of exponents that feature, given the degree of p1, #dp1, and the degree of q1, dq1 kamap:=0: for i from 0 to dq1-1 do for j from 0 to (dq1-1)*dp1-i do kamap:=kamap+1: mapa[kamap,1]:=i; mapa[kamap,2]:=j: od: od:

kamaq:=0: for i from 0 to dp1-1 do for j from 0 to (dp1-1)*dq1-i do kamaq:=kamaq+1: maq[kamaq,1]:=i: maq[kamaq,2]:=j: od: od: kamakh:=0: for i from 1 to dp1+dq1-1 do for j from 0 to dp1*dq1-i do kamakh:=kamakh+1: khezka[kamakh,1]:=i; khezka[kamakh,2]:=j; od: od;

with(linalg): #The following function converts an operator p(n,x) from #Nj*ni form to ni*Nj form” convert1:=proc(p) local i,p1,pa,r,term: p1:=expand(p): pa:=0: r:=degree(p,N): for i from 0 to r do term:=coeff(p1,N,i): term:=expand(term): term:=subs(n=n+i,term): term:=expand(term): pa:=pa+term*Ni: od: RETURN(expand(pa)): end:

timesn:= proc(p): RETURN(expand(n*p-N*diff(p,N))): end: timesk:= proc(p): RETURN(expand(k*p-K*diff(p,K))): end:

timeskn:=proc(p,exk,exn) local pa,i : pa:=p: for i from 1 to exk do pa:=timesk(pa): pa:=expand(pa): od: for i from 1 to exn do pa:=timesn(pa): pa:=expand(pa): od: RETURN(pa): end:

p1:=expand(p1): q1:=expand(q1): kama:=kamap+kamaq: ma:=array(1..kama,1..kama): for i from 1 to kamap do p11:=timeskn(p1,mapa[i,1],mapa[i,2]): p11:=expand(p11): for j from 1 to kama-1 do tem:=coeff(p11,k,khezka[j,1]): tem:=expand(tem): tem:=coeff(tem,n,khezka[j,2]): ma[i,j]:=expand(tem): od: p11:=expand(p11): tem:=coeff(p11,k,0): ma[i,kama]:=expand(tem): od:

for i from kamap+1 to kamap+kamaq do p11:=timeskn(q1,maq[i-kamap,1],maq[i-kamap,2]): p11:=expand(p11): for j from 1 to kama-1 do tem:=coeff(p11,k,khezka[j,1]): tem:=expand(tem): tem:=coeff(tem,n,khezka[j,2]): ma[i,j]:=expand(tem): od: tem:=coeff(p11,k,0): ma[i,kama]:=expand(tem): od:

equ:=det(ma): if equ=0 then print(‘the determinant vanishes, take an appropriate subdeterminant‘): fi: equ:=subs(K=K+1,equ): equ:=expand(equ): equ1:=coeff(equ,K,0): while equ1=0 do equ:=expand(equ/K): equ1:=coeff(equ,K,0): od:

equ:=convert1(equ): equ:=subs(K=K-1,equ): print(‘F(n,k) is annihilated by the following, kfree, operator(check!)‘): print(‘(Warning: it should be in the standard form in which n is in)‘): print(‘front of N, if it is not , it is still meant that way‘): print(equ): equ1:=convert1(equ1): equ1:=expand(equ1): print(‘Therefore, putting K=1, we get an operator A(N,n) annihilating the rhs‘): print(equ1): timesx:= proc(p,nup) local gu: gu:=expand(subs(n=n+nup,p))*Nnup: RE-TURN(gu): end:

aa:=1: sa:=0: ra:=equ1: qaq:=conj: degq:=degree(qaq,N): while degree(ra,N)¿=degq do degr:=degree(ra,N): nu:=degr-degq: qaqa:=expand(timesx(qaq,nu)): leadr:=coeff(expand(ra),N,degr): leadqa:=coeff(qaqa,N,degr): gcd(leadr,leadqa,’leadr’,’leadqa’): aa:=aa*leadqa: sa:=sa*leadqa: sa:=sa+leadr*Nnu: sa:=expand(sa): ra:=leadqa*ra-leadr*qaqa: ra:=expand(ra): od:

print(‘the conjectured r.h.s is annihilated by the following operator B(N,n) ‘): print(conj): if( ra=0) then print(‘Now A(N,n)=S(N,n)*B(N,n) (check!), where S(N,n) is‘): print(sa): gu:=degree(sa,N): print(‘Since S has order‘,gu,‘in N, THE IDENTITY IS TRUE provided‘): print(‘it is true for n=0,1,...,‘,gu-1,‘ and at the positive integer‘): print(‘roots ot the leading term of S, which is‘): sac:=coeff(sa,N,gu): print(factor(sac)): else print(‘since A(N,n) is not a left multiple of B(N,n), THE IDENTITY‘): print(‘IS FALSE‘): fi: quit A.2 Ten Sample Inputs To The Program The above program ran successfully (on my PC) with the following ten inputs. To test them yourself just replace the part between ”#THIS IS THE INPUT” and ”#END OF INPUT” in the program by the ones given below.

1) An identity of George Andrews[And2], ((1.1) of [Ges-Sta]): #THIS IS THE INPUT expli:=0: summand:=(-1)**k*(n+3*a+k-1)!*(a+k-1)!*3k/ ((n-k)!*(3*a/2+k-1)!*(3*a/2+k-1/2)!*k!*4k): conj:=(n+3)*N**3-(3*a+n): #END OF INPUT

2) An identity of Gessel and Stanton((5.24) of [Ges-Sta]). #THIS IS THE INPUT expli:=1: summand:=(-1)**k*(n+k)!*8k/((2*n-k)!*(1/3+k)!*k!*9k): rhs:=(- 1)n*(n-1/2)!*(1/6)!*n!/((-1/2)!*(1/6+n)!*(2*n)!*(1/3)!*3n): #END OF INPUT

3) With this input the program will prove that ([vdP],p.202)

a(n) := ∑_{k}({n \atop k})({n + k \atop k})

satisfies the recurrence

(n + 2)a(n + 2) − (6n + 9)a(n + 1) + (n + 1)a(n) = 0

#THIS IS THE INPUT expli:=0: summand:=(n+k)!/(k!2*(n-k)!): conj:=(n+2)*N**2-(6*n+9)*N+(n+1): #END OF INPUT

4) Saalchtuz’s identity (e.g. Knuth[Knu],1.2.6., ex. 31) #THIS IS THE INPUT expli:=1: summand:=(b+k)!/(k!*(a-b+c-k)!*(n-k)!*(b-c+k)!*(b+k-a-n)!): rhs:=b!*c!*(a+n)!/(a!*(b-a)!*n!*(c-n)!*(a-b+c)!*(n+b-c)!): #END OF INPUT

5) The Vandermonde-Chu identity (≡, terminating case of Gauss’s ””_{2}F_{1}(1)). #THIS IS THE INPUT expli:=1: summand:=1/(k!*(n-k)!*(c-k)!*(b-c+k)!): rhs:=(n+b)!/(n!*b!*c!*(n+bc)!): #END OF INPUT 6)It follows from formula (21) in section 8.2 of Riordan’s book [Rio] (page 201) that the so-called ”reduced straight Menage numbers”, V_{n}, satisfy the recurrence

V_{n+3} − (n + 2)V_{n+2} − (n + 2)V_{n+1} − V_{n} = 0.  (∗)

The explicit expression for V_{n} is given by (take t = 0, k → n − k in formula (6) of section 8.2 (p.197) of [Rio])

V_{n} = ∑_{k}({n + k \atop 2k})k!(−1)^{n−k}  (∗∗)

The program, with the following input, proves that indeed (**) implies (*).

#THIS IS THE INPUT expli:=0: summand:=(-1)(n-k)*(n+k)!*k!/((2*k)!*(n-k)!): conj:=N**3- (n+2)*N**2-(n+2)*N-1: #END OF INPUT

7) With the following input the program proves a recurrence obtained by Sister Celine Fasenmayer ([Rai],pp.234-235,(2) and (11)).

#THIS IS THE INPUT expli:=0: summand:=(-x)**k*(n+k)!/(k!2*(k-1/2)!*(n-k)!): conj:= (n+3)*N**3- (3*n+7-4*x)*N**2+ (3*n+4*x+5)*N-(n+1): #END OF INPUT

8) With the following input the program proves another recurrence of Sister Celine Fasenmayer ([Rai],p.236, (14) and (22)). Note that ”al” stands for α and ”be” stands for β.

#THIS IS THE INPUT expli:=0: summand:=(-x)**k*(be+k)!/(k!2*(al+k)!*(n-k)!): conj:= (n+3)2*(al+n+3)*N**3- (3*(n+3)2-3*(n+3)+1+al*(2*n+5)- (be+n+3)*x)*N2+ (al+3*n+6-x)*N-1: #END OF INPUT

9) With the following input the program proves that the Jacobi polynomials satisfy their recurrence ([Rai], (4) on p.255 and (1) of section 137 on p. 263). Once again ”al” stands for α and ”be” stands for β.

#THIS IS THE INPUT #Recurrence or Jacobi expli:=0: summand:=(al+be+n+k)!*((x-1)/2)**k/(k!*(n-k)!*(al+k)!): conj:=(2*n+4)*(al+be+2*n+2)*(al+n+1)*(al+n+2)*N2- (al+be+2*n+3)*((al2-be2)+x*(al+be+2*n+4)*(al+be+2*n+2))*(al+n+1)*N +2*(al+n+1)*(be+n+1)*(al+be+2*n+4)*(al+be+n+1): #END OF INPUT

10) With the following input the program proves that the Hahn polynomials indeed satisfy their three term recurrence ([Lab]). Here ”a” stands for α, ”b” stands for β, and ”N1” stands for N. (We cannot use ”N” since it has a special meaning , namely the shift operator in the n direction).

#THIS IS THE INPUT #Recurrence for Hahn expli:=0: summand:=(-1)**k*(n+a+b+k)!*(N1-k)!/((n-k)!*(x-k)!*(a+k)!*k!): conj:=(n+a+2)*(N1-n-1)*(n+2)*(2*n+a+b+2)*N**2- ((n+a+b+2)*(n+a+2)*(N1n-1)*(2*n+a+b+2)+(n+1)* (n+b+1)* (n+a+b+N1+2)* (2*n+a+b+4)-(2*n+a+b+2)* (2*n+a+b+3)* (2*n+a+b+4)*x)* N + (n+b+1)* (n+a+b+N1+2)* (n+a+b+1)* (2*n+a+b+4): #END OF IN-PUT

A.3 A Sample Output

1)When the program was run with the original (trivial) input, the output was as follows: (I have deleted the comments and the MAPLE logo)

F(n,k) is annihilated by the following, k-free, operator(check!) (Warning: it should be in the standard form in which n is in) front of N, if it is not , it is still meant that way - K - 1 + n N (K - 1) + (K - 1) N + n N + N Therefore, putting K=1, we get an operator A(N,n) annihilating the rhs n N + N - 2 the conjectured r.h.s is annihilated by the following operator B(N,n) (n + 1) N - 2 Now A(N,n)=S(N,n)*B(N,n) (check!), where S(N,n) is 1 Since S has order, 0, in N, THE IDENTITY IS TRUE provided it is true for n=0,1,...,, -1, and at the positive integer roots ot the leading term of S, which is 1 words used=57169, alloc=36864, time=8.82

REFERENCES

[Alm] Almkvist (Gert), Zeilbergers bevismaskin, Elementa, to appear, (in Swedish).

[A-Z] Almkvist (Gert) and Zeilberger (Doron), The method of differentiation under the integral sign, J. Symbolic Computation, to appear.

[And1] Andrews (George E.), q-Series:Their Development And Applications in Analysis, Number Theory, Combinatorics, Physics, and Computer Algebra, CBMS series 66, Amer. Math. Soc., Providence, 1986.

[And2] Andrews (George E.), Connection Coefficient problems and partitions, in ”AMS Proc. Symposia in Pure Mathematics 34 ”, D.Ray-Chaudhuri, ed., American Mathematical Society, Providence RI, 1979, 1-24.

[Ask1] Askey (Richard A.), Preface, in ”Special Functions: Group Theoretical Aspects and Applications”, ed. by R.A.Askey, T.H.Koornwinder, and W.Schempp, D.Reidel Publ. Co., Dordrecht, 1984.

[Ask2] Askey (Richard A.), Review of [Zei1], Mathematical Reviews MR 83f:33001.

[Ask-Gas] Askey (Richard A.) and Gasper (George), Inequalities for polynomials, in ”The Bieberbach Conjecture”, edited by A.Baernstein et. al., American Mathematical Society, Providence, 1986.

[Ask-Wil] Askey (Richard A.), and Wilson (James A.), Some basic hypergeometric orthogonal polynomials that generalize Jacobi polynomials, Memoirs of the Amer. Math. Soc. 318, 1985.

[Bai] Bailey (Wifrid N.), ”Generalized Hypergeometric Series”, Cambridge Math. Tracts 32, Cambridge University Press, London, 1935. (Reprinted: Hafner, New York, 1964 .)

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

[Ber2] Bernstein (I.N. [Joseph N.]),The analytic continuation of generalized functions with respect to a parameter, Functional Analysis and its Applications, Vol 6 No.4 (1972), Russian original: 26-40, English translation:273-285.

[Bjo1] Bjork (Jan-Erik), ”Rings Of Differential Operators”, North-Holland, Amsterdam, 1979.

[Bjo2] Bjork (Jan-Erik), Differential systems on algebraic manifolds, preprint, Univ. of Stockholm, 1984.

[deB] de Brange (Louis), A proof of the Bieberbach conjecture, Acta Math. 154(1985), 137-152.

[Buc] Buchberger (Bruno), Grobner bases - an algorithmic method in polynomial ideal theory, Chapter 6 in N.K.Bose (ed.): ”Multidimensional System Theory”, D.Reidel, Publ. Comp., 1985.

[Cip] Cipra (Barry), How the grinch stole mathematics, Science 245, 595 (11 August 1989).

[Ehr1] Ehrenpreis (Leon), ”Fourier Analysis In Several Variables”, North Holland, Amsterdam, 1973.

[Ehr2] Ehrenpreis (Leon), Hypergeometric functions, preprint, Temple University.

[Fas] Fasenmayer (Sister M. Celine), Some generalized hypergeometric polynomials, Bull. Amer. Math. Soc. 53 (1947), 806-812.

[Foa] Foata (Dominique),Combinatoire des identites sur les polynomes orthogonaux,in ”Proc.Inter.Congress of Math. [Warsaw, Aug. 16-24, 1983]”, 1541-1553, Varsovie, 1983.

[Gal]Galligo (Andre), Some algorithmic questions on ideals of differential operators, in: Lecture Notes in Computer Science 204, (edited by Bob F.Cavines, EUROCAL ’85, vol.2), Springer, Berlin.

[Gas] Gasper (George), Summation, transformation, and expansion formulas for bibasic series, preprint, Northwestern University.

[Gas-Rah1] Gasper (George), and Rahman (Mizan),An indefinite bibasic summation formula and some quadratic, cubic, and quartic summation and transformation formulas, preprint.

[Gas-Rah2] Gasper (George), and Rahman (Mizan),”Basic Hypergeometric Series”, Cambridge University Press, to appear.

[Gel] Gel’fand (I.M.), General theory of hypergeometric functions, Dokl. Akad. Nauk. USSR, 286(1986), No.1 , 14-18,translated in: Soviet Math. Doklady 33 (1986), No.3, 573-577.

[Gel-Gel] Gel’fand (I.M.) and Gel’fand (S.I.), Generalized hypergeometric equations, Dokl. Akad. Nauk. USSR 288 (1986), No.2, 279-282, translated in: Soviet Math. Doklady 33 (1986), No.3, 643-646.

[Gel-Gra] Gel’fand (I.M.) and Graev (M.I.), A duality theorem for general hypergeometric functions, Dokl. Akad. Nauk. USSR, 289(1986), No.1, translated in: Soviet Math. Doklady 34 (1987),No.1, 9-13.

[Gel-Zel] Gel’fand (I.M.) and Zelevinskii (A.V.), Algebraic and combinatorial aspects of the general theory of hypergeometric functions, Functional Analysis and Its Applications, 20, no.3 (1986), 17-34, translation:20 (1986), 183-197.

[Ges1] Gessel (Ira), Two theorems on rational power series, Utilitas Math., 19, 1981, 247-254.

[Ges2] Gessel (Ira), Counting Latin rectangles, Bull.(New Series) of the Amer. Math. Soc. 16 (1987), 79-82.

[Ges3] Gessel (Ira), Symmetric functions and P-recursiveness, preprint, Brandeis University.

[Ges-Sta] Gessel (Ira), and Stanton (Dennis),Strange evaluations of hypergeometric series, SIAM J. Math. Anal. 13 (1982), 295-308.

[Gos] Gosper (Bill), Letters to R.Askey and D.Stanton, private communication.

[Goul] Goulden (Ian P.), Exact values for degree sums over strips of Young diagrams, preprint, University of Waterloo.

[Gouy] Gouyou-Beauchamps (Dominique), Standard Young tableaux of height 4 and 5, preprint.

[Knu] Knuth (Donald E.), ”Fundamental Algorithms”, The Art of Computer programming vol.1, second edition, Addison-Wesley, Reading, MA, 1973.

[Lab] Labelle (Jacques), Tableau d’Askey, in ”Polynomes orthogonaux et applications”[ Bar-le-Duc 1984 : C. Brezinski et. al., eds.] XXXVI-XXXVII, Lecture Notes in Math. 1171, Springer, Berlin, 1985.

[Lip1] Lipshitz (Leonard), The diagonal of a D-finite power series is D-finite, J.Algebra 113(1988), 373-378.

[Lip2] Lipshitz (Leonard), D-finite power series, J. Alg, to appear.

[Mac1] Macdonald (Ian G.), Some conjectures for root systems, SIAM J. Math. Anal. 13 (1982), 91-143.

[Mill] Miller (Willard), ”Symmetry and Separation of Variables”, ”Encyclopedia of Mathematics and its Applications”, vol.4, Addison-Wesley, London, 1977.

[Miln] Milne (Steven C.), Hypergeometric series well poised in SU(n) and a generalization of Biedenharn’s G function, Advances in Math. 36(1980), 169-211.

[Pal] Palamadov (V.P.), ”Linear Differential Operators With Constant Coefficients”, Grundlehren 168, Springer, Berlin.

[Per] Perlstadt (Marci), Recurrences for sums of powers of binomial coefficients, J. Number Theory 27(1987), 304-309.

[vdP] van der Poorten (Alfred), A proof that Euler missed... Apery’s proof of the irrationality of ζ(3), Math. Intell. 1 (1979), 195-203.

[Rai] Rainville (Earl D.), ”Special Functions”, Macmillan, New York, 1960.(Reprinted: Chelsea, Bronx, NY, 1971).

[Reg] Regev (Amitai), Asymptotic values for degrees associated with strips of Young diagrams, Adv. in Math. 41 (1981), 115-136.

[Rio] Riordan (John), ”An Introduction to Combinatorial Analysis”, Princeton University Press, Princeton,1980. (Originally published by John Wiley, 1958).

[Roy] Roy (R.), Binomial identities and hypergeometric series, Amer. Math. Monthly, 94(1987), 37-46.

[Stan] Stanley (Richard), Differentiably finite power series, European J. Combin. 1 (1980), 175-188.

[Staf] Stafford (J.T.), Module structure of Weyl algebras, J. London Math. Soc. (2), 18(1978), 429-442.

[Syl] Sylvester (James J.),A method of determining by mere inspection the derivatives from two equations of any degree, Philosophical Magazine 16 (1840), 132-135. [≡ Collected Works I , 54-57].

[Tak1] Takayama (Nobuki), Holonomic solutions of Weisner’s operator, preprint, Tokushima University, Japan.

[Tak2] Takayama (Nobuki),Grobner basis and the problem of contiguous relations,preprint, Tokushima University, Japan.

[Vdw] van der Waerden (B.L.), ”Modern Algebra”, Volumes I and II, Frederik Ungar Publ. Co., New York.

[Wil] Wilf (Herbert S.), What is an answer?, Amer. Math. Monthly, 89(1982), 289-292. [W-Z1] Wilf (Herbert S.) and Zeilberger (Doron), Towards computerized proofs of identities, Bulletin of the Amer. Math. Soc., to appear. [W-Z2] Wilf (Herbert S.) and Zeilberger (Doron), Rational functions certify combinatorial identities, J. Amer. Math. Soc., to appear.

[Zei1] Zeilberger (Doron), Sister Celine’s technique and its generalizations, J. Math. Anal. Anal. Appl. 85(1982), 114-145.

[Zei2] Zeilberger (Doron), A fast algorithm for proving terminating hypergeometric identities, preprint.

[Zei3] Zeilberger (Doron), The method of creative telescoping, preprint.

[Zei4] Zeilberger (Doron), The method of creative telescoping for q-series, in preparation. Department of Mathematics and Computer Science, Drexel University, Philadelphia, PA19104, USA. June 23, 1988. Revised: Oct 19, 1989. TeX Version: Dec. 15, 1999. Address in 1999: Department of Mathematics, Temple University, Philadelphia, PA19122,USA.