Recurrences and Legendre Transform

Volker Strehl

Institut für Mathematische Maschinen und Datenverarbeitung (Informatik I) Universität Erlangen-Nürnberg, Germany

Abstract

A binomial identity ((1) below), which relates the famous Apéry numbers and the sums of cubes of binomial coefficients (for which Franel has established a recurrence relation almost 100 years ago), can be seen as a particular instance of a Legendre transform between sequences. A proof of this identity can be based on the more general fact that the Apéry and Franel recurrence relations themselves are conjugate via Legendre transform. This motivates a closer look at conjugacy of sequences satisfying linear recurrence relations with polynomial coefficients. The rôle of computer-aided proof and verification in the study of binomial identities and recurrence relations is illustrated, and potential applications of conjugacy in diophantine approximation are mentioned. This article is an expanded version of a talk given at the 29. meeting of the Séminaire Lothringien de Combinatoire, Thurnau, september 1992.

1 Introduction

In this article I will discuss some general aspects related to the following beautiful binomial identity:

∑_{k=0}^{n} \binom{n}{k}^{2}\binom{n + k}{k}^{2} = ∑_{k=0}^{n} \binom{n}{k}\binom{n + k}{k} ∑_{j=0}^{k} \binom{k}{j}^{3} (n ≥ 0)  (1)

An interesting aspect of this identity is the fact that it relates the famous Apéry numbers a_{n} := ∑_{k} \binom{n}{k}^{2}\binom{n+k}{k}^{2}, which played an important rôle in R. Apéry’s original proof of the irrationality of ζ(3) (see the entertaining and instructive article [20] by A. v.d. Poorten for an account of this), with the sums of cubes of binomial coefficients f_{n} := ∑_{k} \binom{n}{k}^{3}, for which J. Franel had found recurrences long ago ( [6], see [4], [8], [11], [5] ), and which for that reason will be called Franel numbers in this article.

This identity was brought to my attention early in 1992 as a conjecture by my colleague B. Voigt from Bielefeld. The origin of it lies in the following question put forward by A. L. Schmidt from Copenhagen:

Let numbers c_{k} (k ≥ 0), independent of n, be defined by

∑_{k=0}^{n} \binom{n}{k}^{2}\binom{n + k}{k}^{2} = ∑_{k=0}^{n} \binom{n}{k}\binom{n + k}{k} c_{k} (n ≥ 0)

[It is easy to see (by induction) that these rational numbers are uniquely determined.] Is it then true that these numbers are always integers? [This is not obvious from the definition, but an evaluation of the first few hundred values clearly suggests that.]

It was soon observed by W. Deuber, W. Thumser and B. Voigt that the beginning of the sequence

(c_{k})_{k≥0} = ( 1 , 2 , 10 , 56 , 346 , 2252 , 15184 , 104960 , 739162, …)

is the same as the beginning of the sequence of the Franel numbers f_{k} = ∑_{j=0}^{k} \binom{k}{j}^{3}. From this observation the above conjecture (1) was formulated. Since I could not spot the conjectured identity in the relevant literature (e.g. [12], [8]), I tried to find a proof for myself. In fact, there are now several different proofs available, none of them really ”trivial”, and I will briefly mention below these diverse approaches that I worked out. The ordering corresponds to the chronological ordering in which these proofs were discovered. A detailed discussion of this affair will be given in [18].

1. Identity (1) is hidden in a classical result, due to W.N. Bailey, from the theory of hypergeometric functions. In fact, it can be pulled out of Bailey’s bilinear generating function for the Jacobi polynomials ([2], [17], [19]) via a kind of diagonalization. The special case where the Jacobi polynomials degenerate into Legendre polynomials is sufficient for our identity, but the same trick can be played in the general situation, so that (1) turns out to be a degenerate case belonging to a family of binomial identities with two free parameters.

2. As has been shown by myself in [19], Bailey’s bilinear generating function can be proved from a suitable combinatorial interpretation. This is somewhat involved, and on the way of proving the full bilinear generating function I established a certain family of binomial identities. It turned out that these already include (1) (up to routine transformations) and its generalization mentioned before. This answers positively the question about a possible combinatorial proof of (1), but such a proof is more complicated than the rather simple appearance of this identity might suggest.

3. By using the method of Legendre inverse pairs, as exposed in Riordan’s book [12], the problem can be transformed into an equivalent one to which D. Zeilberger’s method of constructing recurrence operators for definite hypergeometric sums (in the case of single sums) can be applied. This gives a ”computer-supported” proof, which can be turned into a ”conventional” one - once the operators are known.

4. The desired result can be obtained by a series of applications of hypergeometric transforms (D. Stanton, Ch. Krattenthaler and G. Andrews have provided comments on that approach). Once found, such a proof is routine to check - finding the right track is the problem!

5. The Wilf-Zeilberger-method (see remarks below) for multiple hypergeometric sums can be applied directly to the identity. Doron Zeilberger provided such a completely ”machine-made” proof in june 1992.

6. Using the known recurrences for the Apéry numbers a_{n} = ∑_{k} \binom{n}{k}^{2}\binom{n+k}{k}^{2} and the Franel numbers f_{n} = ∑_{k} \binom{n}{k}^{3} one can give a short (partially computer supported) proof of (1), using no sophisticated algorithms, but only standard simplification capabilities of any reasonable computer algebra program. Such a proof is outlined in the next section, and indeed, this kind of proof was the motivation for writing the present note. Let me remark that a proof similar to this one was found independently by A. L. Schmidt.

The following sections are organized as follows.

After illustrating the use (neither the principles, nor the theory behind) of available implementations of algorithms for a computer-aided treatment of binomial sums and sequences satisfying linear recurrence relations with polynomial coefficients, I will present a short proof of (1) which is based on the (known) recurrence relations of the Apéry and Franel numbers. The crucial fact is this: identity (1) claims that Franel’s and Apéry’s sequences are related via Legendre transform. This concept has been used, with number-theoretic applications in mind, by A. L. Schmidt, see [14], [15], [16], and, in particular, his most recent article is closely related to the present one, and some of his results will be mentioned below. This turns out to be an immediate consequence of the stronger fact that the linear recurrences (of second order, with polynomial coefficients) generating those sequences are themselves conjugates via Legendre transform in a sense to be made precise below. Once guessed correctly, a proof of this conjugacy relation is remarkably simple (relative to a routine use of a machine-based simplification procedure).

In section 3 some notation will be developed in order to deal with this approach from a more general perspective.

The situation of (reduced) Legendre transform is then studied in some detail in section 4, and the question under which condition (and if so: how) recurrence relations are transformed into recurrence relations via conjugation is answered completely in a particular case. It should become clear, however, that similar results could be obtained by the same approach in related and more general situations. Thus these results should be taken as exemplary - no complete treatment was intended here, just an outline of how one may proceed.

In the final section I present some consequences in the field of diophantine approximation which are directly related to identity (1) and the conjugacy relation behind it. These results are taken from A. L. Schmidt’s recent article [16] mentioned above - here they should be taken as a motivation showing that the study of binomial (or hypergeometric) identities and recurrences related to them might lead to interesting applications in number theory.

2 A short proof of identity (1)

Two basic facts are used for a proof of (1) as presented here:

• [Franel, 1895] The Franel numbers f_{n} = ∑_{k} \binom{n}{k}^{3} satify the following secondorder recurrence with polynomial coefficients ^{1}

(n + 1)^{2}f_{n+1} − (7n^{2} + 7n + 2)f_{n} − 8n^{2}f_{n−1} = 0 (n ≥ 0)  (2)

• [Apéry, 1978] The Apéry numbers f_{n} = ∑_{k} \binom{n}{k}^{2}\binom{n+k}{k}^{2} satify the following second-order recurrence with polynomial coefficients^{2}

(n + 1)^{3}a_{n+1} − ((n + 1)^{3} + n^{3} + 4(2n + 1)^{3}) a_{n} + n^{3}a_{n−1} = 0 (n ≥ 0) (3)

Note that even if one knows the recurrences beforehand, it can be rather tedious to verify them. To illustrate this situation, let me cite (with a slight change in notation)^{3} from A. v.d. Poorten’s article [20]:

... To convince ourselves of the validity of Apéry’s proof we need only complete the following exercise: Let a_{n} = ∑_{k} \binom{n}{k}^{2}\binom{n+k}{k}^{2}, then a_{0} = 1, a_{1} = 5 and the sequence \{ a_{n} \} satisfies the recurrence (3). ... Neither Cohen nor I had been able to prove this in the intervening two months ...

Then v.d. Poorten mentions that D. Zagier was able to provide a solution to this problem with irritating speed by virtue of a technique called creative telescoping. Actually, this technique, combined with the powerful algorithm for indefinite hypergeometric summation due to W. Gosper [7], is at the heart of D. Zeilberger’s method. Today, finding (and proving!) recurrences like (2) and (3) is a routine application of Zeilberger’s algorithm for simple terminating hypergeometric sums, in fact a matter of very few seconds (in these cases! In a case like the one mentioned above in 5., i.e. directly verifying that the right hand side of (1) satisfies the same recurrence as the Apéry numbers do, application of Zeilberger’s algorithm for hypergeometric multisums leads to nonnegligible computing time, even on fast machines^{4}. To illustrate this, I include input and output from an implementation of Zeilberger’s algorithm which was programmed by J. Hornegger [9] in the context of the AXIOM^{5} computer algebra system.

1Franel has also found a similar recurrence for ∑_{k=0}^{n} \binom{n}{k}^{4}. It took quite a while until corresponding recurrence relations for higher exponents were found, see the recent articles by Perlstadt [11] and Cusick [5].

2Apéry has also found a similar recurrence for ∑_{k=0}^{n} \binom{n}{k}^{2}\binom{n+k}{k}. See the article by van der Poorten [20] and an article [1] by R. Askey and J. Wilson for illuminating remarks from the hypergeometric viewpoint and a generalization.

3Actually, v.d. Poorten writes about two solutions of the recurrence, the other one, call it { b_{n} } specified by initial values b_{0} = 0 and b_{1} = 6. This second solution will only appear in the last section of this article. Note that the rôle of a and b has been interchanged w.r.t. [20].

4D. Zeilberger’s first direct verification of (1) took about 1500 sec.

Recurrences for ∑_{k=0}^{n} \binom{n}{k}^{e} (e = 1, 2, 3)

k1 : Symbol := k
p1 : SMP(FRAC INT, Symbol)      := n
p2 : SMP(FRAC INT, Symbol)      := k
arg: FRAC SMP(FRAC INT, Symbol) := 1

Exponent e = 1:

(69)  r1:= findrec(bin(p1, p2, n1, k1, 1), 1, arg, 1, n1, k1)
  (69)  - B + 2, 1
Type: Union(List FRAC SMP(FRAC INT,Symbol),...)
             Time: 0.08 (IN) + 0.13 (EV) + 0.04 (OT) = 0.25 sec

The result consists of the recurrence operator (written as a polynomial in the shift operator B and n) and the certificate provided by the method (plus informations about typing and timing).

Exponent e = 2:

(71)  r2:= findrec(bin(p1, p2, n1, k1, 2), 1, arg, 1, n1, k1)
  (71)  (B - 4)n + B - 2, - 3n + 2k - 1
Type:  ...
            Time: 0.08 (IN) + 0.31 (EV) + 0.05 (OT) = 0.44 sec

This result expresses the fact that ∑_{k=0}^{n} \binom{n}{k}^{2} = \binom{2n}{n}, because both sides are annihilated by the recurrence operator.

Exponent e = 3:

(72)  r3:= findrec(bin(p1, p2, n1, k1, 3), 1, arg, 1, n1, k1)
  (72)
   2            2       2                    2
 (B  - 7 B - 8)n  + (4 B  - 21 B - 16)n + 4 B  - 16 B - 8,
              5              4         2               3
         - 14n  + (27k - 75)n  + (- 18k  + 111k - 161)n
       +
            3      2               2      3      2
         (4k  - 54k  + 171k - 173)n  + (8k  - 54k  + 117k - 93)n
       +
           3      2
         4k  - 18k  + 30k - 20

5AXIOM is a trademark of NAG, Numerical Algoriths Group Ltd.,Oxford, England.

    /
        3              2      2               3     2
       n  + (- 3k + 3)n  + (3k  - 6k + 3)n - k  + 3k  - 3k + 1
Type:  ...
         Time: 0.08 (IN) + 13.30 (EV) + 0.38 (OT) +4.15 (GC) = 17.91 sec
(73)  factorCoefficients r3.1
                2 2      2        16             2
   (73)  (n + 2) B  - 7(n  + 3n + --)B - 8(n + 1)
                                  7
Type: SUP Factored SMP(FRAC INT,Symbol)
         Time: 0.01 (IN) + 0.11 (EV) + 0.09 (OT) = 0.21 sec

This is Franel’s recurence. The certificate in factorized form:

(74)  factorCertifyingFunction r3.2
   (74)
                 2
      - 14(n + 1)
   *
      3      27     47  2    9  2   57     53     2  3   9  2  15     10
    (n  + (- -- k + --)n  + (- k  - -- k + --)n - - k  + - k - -- k + --)
             14     14       7      14     14     7      7      7      7
 /
               3
    (n - k + 1)
Type: FRAC Factored SMP(FRAC INT,Symbol)
        Time: 0.04 (IN) + 0.78 (EV) + 0.48 (OT) + 1.50 (GC) = 2.80 sec

Now for Apéry’s recurrence (omitting the certificate):

a1: SMP(FRAC INT, Symbol):= n
a2: SMP(FRAC INT, Symbol):= k
a3: SMP(FRAC INT, Symbol):= a1 + a2
r:= recurrence  findrec(bin(a1, a2, n1, k1, 2),
     bin(a3, a2, n1, k1, 2), 1, 2, n1, k1)
  3     2            2         3       2
(n  + 6n  + 12n + 8)B  + (- 34n  - 153n  - 231n - 117)B +
       3     2
      n  + 3n  + 3n + 1
Type: SUP SMP(FRAC INT, Symbol)
        Time: 0.20 (IN) + 13.71 (EV) + 0.23 (OT) = 14.14 sec

Here I will not give any further explanations w.r.t. Zeilberger’s method. The interested reader may look into the numerous articles available (e.g. [25], [26], [21], [22], [24], [23], [3]), [10], [9]) describing the method, the algorithms, and their applications.

As a side remark: I would like to take the occasion to point out another approach, experimental in character, which can be helpful in situations similar to the one discussed here. There is now a nice set of tools available with which given the first few coefficients of a generatin function, one may cleverly guess the generating function (if it has a ”nice” closed form), or a linear differential equation satisfied by it, or a linear recurrence relation satisfied by its coefficients. The mathematical model behind it is that of holonomic functions of one variable (or equivalently: D-finite power series, or P-recursive sequences), with important special cases such as rational and algebraic generating functions. Even if these tools only provide guesses and not proofs in the strict sense, they often can be used to obtain enough information about a problem so that a strict solution becomes feasible. These tools, written by B. Salvy and P. Zimmerman (INRIA, Paris) are available under the Maple^{6} computer algebra system as a package named gfun in the Maple share library ^{7}. A rather comprehensive description [13] will soon be available from the authors.

As an illustration, let us guess the recurrences (2) and (3) using gfun:

    |\^/|      MAPLE V
._|\|   |/|_.  Copyright (c) 1981-1990 by the University of Waterloo.
 \  MAPLE  /   All rights reserved.  MAPLE is a registered trademark of
 <________>   Waterloo Maple Software.
      |        Type ? for help.
> readlib(gfun):with(gfun):

The Apéry numbers are defined:

a := ⟨n ↦ ∑_{k=0}^{n} \binom{n}{k}^{2}\binom{n + k}{k}^{2}⟩

A list of the first few values:

> [a(i)$i=0..10];
       [1, 5, 73, 1445, 33001, 819005, 21460825, 584307365,
       16367912425, 468690849005, 13657436403073]

gfun does not find a ”nice” generating function:

> guessgf([a(i)$i=0..15],t);
                                      FAIL

Eleven initial values are not sufficient for the recurrence:

6Maple is a trademark of Waterloo Maple Software, Waterloo, Canada.

7The Maple share library be accessed using anonymous ftp from daisy.uwaterloo.ca in Waterloo (Canada) or neptune.inf.ethz.ch in Zurich (Switzerland).

> listtorec([a(i)$i=0..10],A(n));
                                      FAIL

but sixteen initial values are:

> listtorec([a(i)$i=0..15],A(n));

[ {n^{3}A(n) + (−5 − 27 n − 51 n^{2} − 34 n^{3}) A(n + 1)
+ (1 + 3 n + 3 n^{2} + n^{3}) A(n + 2), A(1) = 1, A(2) = 5} , ogf ]

We can get a differential equation from the recurrence:

> rectodiffeq("[1],A(n),aa(z));

[ D(aa)(0) = 1, D^{(2)}(aa)(0) = 10,

(3 z^{4} − 51 z^{3}) \frac{d^{2}}{dz^{2}} aa(z) + (z + z^{3} − 10 z^{2}) \frac{d}{dz} aa(z) + (5 z − 1) aa(z)

−5 A(0)z + A(0) + z^{5} \frac{d^{3}}{dz^{3}} aa(z) + (\frac{d^{3}}{dz^{3}} aa(z)) z^{3} − 34 (\frac{d^{3}}{dz^{3}} aa(z)) z^{4} , ogf ]

Now we consider sums of cubes of binomial coefficients:

f := ⟨n ↦ ∑_{k=0}^{n} \binom{n}{k}^{3}⟩

> [f(i)$i=0..10];
[1, 2, 10, 56, 346, 2252, 15184, 104960, 739162, 5280932, 38165260]

Too few initial values will not lead to a differential equation:

> listtodiffeq([f(i)$i=0..10],F(z));
                                      FAIL
> listtodiffeq([f(i)$i=0..15],F(z));

\{ (−2 − 8 z) F (z) + (1 − 14 z − 24 z^{2}) \frac{d}{dz} F (z) + (z − 7 z^{2} − 8 z^{3}) \frac{d^{2}}{dz^{2}} F (z),

D(F)(0) = 2, F (0) = 1 }

From the differential equation we get the Franel recurrence:

> diffeqtorec(",F(z),ff(n));

{ (−7n − 2 − 7n^{2}) ff (n) + (n^{2} + 2n + 1) ff (n + 1) − 8ff (n − 1)n^{2},
ff (0) = 1, ff (1) = 2 }

Let us now return to the proposed proof of (1). It will be presented in matrix form. To begin with, let us define a doubly infinite tridiagonal matrix F = (f_{i,j})_{i,j≥0} representing the difference operator of the Franel recurrence (2) :

f_{i,j} = f_{i−j}(i) where f_{k}(z) := \{ \begin{matrix}(z + 1)^{2} & if k = −1 \\ −(7z^{2} + 7z + 2) & if k = 0 \\ −8z^{2} & if k = 1 \\ 0 & if |k| > 1 \end{matrix}

so that

F = ( \begin{matrix}f_{0}(0) & f_{−1}(0) & f_{−2}(0) & … \\ f_{1}(1) & f_{0}(1) & f_{−1}(1) & … \\ f_{2}(2) & f_{1}(2) & f_{0}(2) & … \\ ⋮ & ⋮ & ⋮ & ... \end{matrix} ) = ( \begin{matrix}−2 & 1 & 0 & 0 & … \\ −8 & −16 & 4 & 0 & … \\ 0 & −32 & −44 & 9 & … \\ 0 & 0 & −72 & −86 & … \\ ⋮ & ⋮ & ⋮ & ⋮ & ... \end{matrix} )

If we denote by f = (f_{n})_{n≥0} the infinite column vector with the Franel numbers f_{n} = ∑_{k} \binom{n}{k}^{3} as entries, then Franel’s result (2) is equivalent to

F · f = 0

where 0 denotes the zero column vector.

Similarly, the difference operator occuring in the Apéry recurrence (3) has matrix form

A = ( \begin{matrix}a_{0}(0) & a_{−1}(0) & a_{−2}(0) & … \\ a_{1}(1) & a_{0}(1) & a_{−1}(1) & … \\ a_{2}(2) & a_{1}(2) & a_{0}(2) & … \\ ⋮ & ⋮ & ⋮ & ... \end{matrix} ) = ( \begin{matrix}−5 & 1 & 0 & 0 & … \\ 1 & −117 & 8 & 0 & … \\ 0 & 8 & −535 & 27 & … \\ 0 & 0 & 27 & −1463 & … \\ ⋮ & ⋮ & ⋮ & ⋮ & ... \end{matrix} )

where

a_{i,j} = a_{i−j}(i) where a_{k}(z) := \{ \begin{matrix}(z + 1)^{3} & if k = −1 \\ −(z + 1)^{3} − z^{3} − 4(2z + 1)^{3} & if k = 0 \\ z^{3} & if k = 1 \\ 0 & if |k| > 1 \end{matrix}

Then Apéry’s recurrence can be written as

A · a = 0

where a = (a_{n})_{n≥0} is the vector of the Apéry numbers, a_{n} = ∑_{k} \binom{n}{k}^{2}\binom{n+k}{k}^{2}.

Finally we introduce the doubly infinite, lower triangular matrix P = (p_{i,j})_{i,j≥0} of the Legendre transform:

P = (\binom{i}{j}\binom{i + j}{j})_{i,j≥0} = ( \begin{matrix}1 & 0 & 0 & 0 & … \\ 1 & 2 & 0 & 0 & … \\ 1 & 6 & 6 & 0 & … \\ 1 & 12 & 30 & 20 & … \\ ⋮ & ⋮ & ⋮ & ⋮ & ... \end{matrix} )

Note that the n-th row of this matrix contains the coefficients of the Legendre polynomials P_{n}, if written as

P_{n}(z) = ∑_{k=0}^{n} \binom{n}{k}\binom{n + k}{k} (\frac{z − 1}{2})^{k}

Using this notation, it becomes clear that A.L. Schmidt’s question mentioned in the beginning asks for the inverse image of the sequence of Apéry’s numbers under the Legendre transform, and that the conjectured identity (1) claims that Apéry’s sequence a is the image of Franel’s sequence f under the Legendre transform, i.e. a = P · f. Put into matrix terms, what we would like to show is that

A · P · f = 0

Obviously, we would be done if we could exhibit a matrix X such that

A · P = X · F

Computing initial segments of this unknown (infinite) matrix X leads to the following (surprising?) guess:

X = D · P and hence A = D · P · F · P^{−1}  (4)

where D = (d_{i,j})_{i,j≥0} is a diagonal matrix given by d_{i,i} = 4i + 2 (i ≥ 0).

It remains to prove this claim which not only says that the Franel and the Apéry sequences are related via the Legendre transform, but also that the associated difference operators are related in the sense of conjugation via Legendre transform (up to multiplication with a diagonal matrix). Interestingly, even though we have to deal with infinite matrices containing binomial coefficients, i.e. nonrational terms, the proof can be established by rational arithmetic alone. For this to see, let us write the assertion (4) to be proved in the form:

A · P = D · P · F

Consider the (i, j)-entry on both sides of this equation. Since both A and F are tridiagonal matrices, every such term involves three summands only. Write this as

lhs = a_{1}(i) p_{i−1,j} + a_{0}(i) p_{i,j} + a_{−1}(i) p_{i+1,j}
rhs = (4i + 2) [p_{i,j−1} f_{−1}(j − 1) + p_{i,j} f_{0}(j) + p_{i,j+1} f_{1}(j + 1)]

where now i and j are treated as variables. Now ask for simplification of lhs − rhs. The Maple command

     expand(lhs - rhs)

gives back an expression of considerable size. Of course, simplification could be done by hand, but this tedious task is better accomplished by your computer algebra system^{8}:

8The subtle point here is this: the simplification algorithm must recognize the fact that even though there are binomial coefficients with symbolic entries around, the problem is essentially one of rational arithmetic because of the hypergeometric nature of binomial coefficients: walking along rows and columns in the Pascal triangle can locally be performed via multiplication with suitable rational functions of the locations, hence rational normalization is possible as soon as only bounded neighbourhoods have to be taken into account. Here the tridiagonal structure of the matrices A and F is important.

     simplify(expand(lhs - rhs))
                                 0

To conclude this section, let us summarize the above discussion

Theorem 1 The Franel recurrence operator F and the Apéry recurrence operator A are Legendre conjugates in the following sense:

A = D · P · F · P^{−1}

with A, F, D, P as above.

The binomial identity (1) we started with is just one of the consequences of this general fact - perhaps the most interesting one. In section 5 it will be shown how this conjugacy relation can be further exploited.

Note that in the proof of the theorem we essentially used the fact that the matrices A and F are tridiagonal matrices with ”polynomial values along the diagonals”, and that the matrix of the Legendre transform P has a similar structure^{9}. It is therefore natural to ask under which conditions a conjugation operation (with a nonsingular, lower triangular matrix such as P) will transform matrices representing linear recurrence operators with polynomial coefficients into matrices of the same kind. I will not try to give a definitive answer to this question in generality. Instead, I will look into the case of the (reduced) Legendre transform and show how, at least for recurrence operators of low order and with low degree polynomial coefficients, such a question can be attacked.

3 Notation and some generalities

Let Φ = (ϕ_{i}(z))_{i∈Z} be a family of functions defined on the integers. With this family we will associate a doubly-infinite matrix D_{Φ} = (d_{i,j})_{i,j≥0} = (ϕ_{i−j}(i))_{i,j≥0}, i.e.

D_{Φ} = ( \begin{matrix}ϕ_{0}(0) & ϕ_{−1}(0) & ϕ_{−2}(0) & … \\ ϕ_{1}(1) & ϕ_{0}(1) & ϕ_{−1}(1) & … \\ ϕ_{2}(2) & ϕ_{1}(2) & ϕ_{0}(2) & … \\ ⋮ & ⋮ & ⋮ & ... \end{matrix} )

In the previous section, we have seen the matrices A and F belonging to the Apéry and the Franel recurrence, respectively, when written in this form, are defined by polynomial families (a_{i}(z))_{i∈Z} and (f_{i}(z))_{i∈Z} (with three nonzero members in each case). The matrix P of the Legendre transform does not belong to a polynomial family, simply because i ↦ \binom{2i}{i} does not grow polynomially. (One might use the Γfunction in order to represent the matrix by a family of analytic functions. Instead, we will pass to a related polynomial family below).

9No quite, actually: the values ”along the diagonals” are no longer polynomial, but still hypergeometric. That is why in section 4 we will pass to a related transform, called reduced Legendre transform, which is polynomial and thus behaves a bit nicer.

In particular, if Φ = (t_{i}(z))_{i∈Z} is such that t_{i} = 0 for i < 0 and t_{0}(j) ≠ 0 for all j ∈ N, i.e. D_{Φ} is an invertible lower triangular triangular matrix, then its inverse matrix (D_{Φ})^{−1} can be written as D_{Ψ} where the family Ψ = (t′_{i}(z))_{i∈Z} is given by t′_{i}(z) = 0 for i < 0, t′_{0}(j) = 1/t_{0}(j) for all j ∈ N, and for i > 0, (−1)^{i}t′_{i}(z) is given by:

\frac{det ( \begin{matrix}t_{1}(z − i + 1) & t_{0}(z − i + 1) & 0 & 0 & … & 0 \\ t_{2}(z − i + 2) & t_{1}(z − i + 2) & t_{0}(z − i + 2) & 0 & … & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ... & ⋮ \\ t_{i−1}(z − 1) & t_{i−2}(z − 1) & … & … & … & t_{0}(z − 1) \\ t_{i}(z) & t_{i−1}(z) & t_{i−2}(z) & t_{i−3}(z) & … & t_{1}(z) \end{matrix} )}{t_{0}(z − i) · . . . · t_{0}(z)}

This shows in particular: if Φ is a polynomial family in the sense that all t_{i}(z) are polynomials, and if in addition t_{0}(z) is a non-vanishing constant, then the inverse family Ψ is of the same kind. Moreover, if deg t_{i}(z) ≤ i for all i ≥ 0, then expansion of the determinant shows that t′_{i}(z) is a sum of terms of degree ≤ i, hence deg t′_{i}(z) ≤ i for i ≥ 0.

Now take a function g(z) and δ ∈ Z. The singleton family Γ_{g,δ} = (g_{i}(z))_{i∈Z} is given by g_{δ}(z) = g(z) and g_{i}(z) = 0 for i ≠ δ. Then consider the family Λ = (h_{d}(z))_{d∈Z} given by conjugation

D_{Φ} · D_{Γ} · (D_{Φ})^{−1} = D_{Λ}

with a ‘triangular’ family Φ = (t_{i}(z))_{i∈Z} as before. It is easy to see that h_{d}(z) = 0 for d < δ, and that for d ≥ δ one gets

h_{d}(z) = ∑_{0≤k≤d−δ} t_{k}(z) · g(z − k) · t′_{d−δ−k}(z − k − δ)  (5)

which can be written in determinantal form (and thus avoiding explicit reference to the inverse family t′_{i}(z)) as

h_{d}(z) = \frac{det M(Φ, g, δ, d)}{t_{0}(z − d) · . . . · t_{0}(z − δ)}

where the square matrix M(Φ, g, δ, d) of size d − δ + 1 is given by

( \begin{matrix}t_{d−δ}(z) · g(z − d + δ) & t_{d−δ−1}(z) · g(z − d + δ + 1) & … & t_{0}(z) · g(z) \\ t_{1}(z − d + 1) & t_{0}(z − d + 1) & … & 0 \\ t_{2}(z − d + 2) & t_{1}(z − d + 2) & … & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ t_{d−δ}(z − δ) & t_{d−δ−1}(z − δ) & … & t_{0}(z − δ) \end{matrix} )

Either of these presentations shows that for a polynomial g(z) and δ ∈ Z the family Λ = (h_{i}(z))_{i∈Z} obtained via conjugation with a polynomial family Φ = (t_{i}(z))_{i∈Z} is again polynomial (and vanishes for indices d < δ).

By superposition and linearity, if Γ = (g_{i}(z))_{i∈Z} is a family of polynomials such that g_{i}(z) = 0 for all i less than some integer K, and if Φ is as before, than conjugation

D_{Φ} · D_{Γ} · (D_{Φ})^{−1}

is well-defined and gives a family of polynomials Λ = (h_{d}(z))_{d∈Z} which vanishes for d < K, and for d ≥ K its members are given by

h_{d}(z) = ∑_{K≤δ≤d} ∑_{0≤k≤d−δ} t_{k}(z) · g_{δ}(z − k) · t′_{d−δ−k}(z − k − δ)

In general, if the family Γ = (g_{i}(z))_{i∈Z} has a finite number of nonzero terms - and hence corresponds to a linear difference operator with polynomial coefficients - the same will not be true for the members of the conjugate family Λ = (h_{i}(z))_{i∈Z}, as the case of a singleton family Γ_{g,δ} clearly shows. We also expect the degree of the conjugate polynomials h_{i}(z) to be higher than the degrees of the g_{j}(z). Thus it will come rather as an exception if the conjugate under Φ-transform of a finite polynomial family Γ (in the above sense) is again of the same type.

In the next section we will study this situation with respect to a variant of the Legendre transform.

4 (Reduced) Legendre transform preserves degrees

We consider now the matrix P = (p_{i,j})_{i,j≥0} belonging to the Legendre transform, i.e.

p_{i,j} = \binom{i}{j}\binom{i + j}{j} = \binom{2j}{j}\binom{i + j}{i − j}

which corresponds to a family Π = (p_{k}(z))_{k∈Z} whose members do not behave polynomially. Since it is slightly more convenient to work with polynomials only, we will write

P = Q · Δ

where

Q = (q_{i,j})_{i,j≥0} with q_{i,j} = \binom{i + j}{i − j} and Δ = diag (\binom{2i}{i})_{i≥0}  (6)

The matrix Q belongs to the polynomial family Φ = (q_{k}(z))_{k∈Z} where q_{k}(z) = \binom{2z−k}{k} (k ≥ 0). Note that the family Φ′ = (q′_{k}(z))_{k∈Z} belonging to the inverse matrix Q^{−1} is given by

q′_{k}(z) = (−1)^{k} [\binom{2z}{k} − \binom{2z}{k − 1}] = (−1)^{k}\binom{2z}{k}\frac{2z − 2k + 1}{2z − k + 1}  (7)

so that q_{0}(z) = q′_{0}(z) = 1 and deg q_{k}(z) = deg q′_{k}(z) = k (k ≥ 0). The transformation provided by the matrix Q and its inverse will be called reduced Legendre transform. Since it differs from the Legendre transform just by multiplication with the diagonal matrix Δ, it will be easy to carry over results from one situation to the other.

Let us now take singleton families Γ_{g,δ} as above, with polynomial g(z) and δ ∈ Z, and study their reduced Legendre transform. Let Λ = (h_{d}(z))_{d∈Z} be the conjugate family obtained from (reduced) Legendre transform, i.e.

D_{Λ} = D_{Φ} · D_{Γ_{g,δ}} · (D_{Φ})^{−1}

Below I will sketch a proof of

Proposition 2 With the notions introduced, one has deg h_{d}(z) ≤ deg g(z) for all d ≥ δ, and hence the same result holds by linearity for any (summable) combination of such singleton families

The proof will only given as a sketch, technical details will be left out. Note that the assertion is robust under linear transformations of the z-variable, so free use can be made of this fact.

From (5),(6), and (7) we know that we have to evaluate

h_{d}(z) = ∑_{k=0}^{n} g(z − k)\binom{2z − k}{k}(−1)^{n−k}\binom{2(z − k − δ)}{n − k}\frac{2(z − δ − n) + 1}{2(z − δ) − n − k + 1}

where n = d − δ ≥ 0. By expanding and linear transformation it turns out that it is sufficient to show that for t ≥ 0 a sum like

∑_{k=0}^{n}(−1)^{n−k} k^{t} \binom{z + d − k}{k}\binom{z − 2k}{n − k}\frac{z − 2n + 1}{z − n − k + 1}

is a polynomial of degree ≤ min(n, t) in z. Now, this follows from

Lemma 3 For each integer t ≥ 0

r_{n,t}(z, d) = ∑_{k=0}^{n}(−1)^{n−k}k(k − 1) … (k − t + 1)\binom{z + d − k}{k}\binom{z − 2k}{n − k}\frac{z − 2n + 1}{z − n − k + 1}

is a polynomial in z and d with:

1. r_{n,t}(z, d) = 0 for 0 ≤ n < t (trivially);

2. for n ≥ t, r_{n,t}(z, d) is a polynomial of total degree n, and its terms of highest degree are given by

\frac{1}{(n − t)!}(z + d)^{t} · d^{n−t}

3. the degree in z of r_{n,t}(z, d) is always ≤ t, and in particular

r_{n,0}(z, d) = \binom{n + d − 1}{n}

For our purposes, we are only interested in the assertion about the z-degree. The proof of this lemma can be reduced to the proof of a similar assertion for sums like

∑_{k=0}^{n}(−1)^{k} k^{t} \binom{z + d − k}{k}\binom{z − 2k}{n − k}

or

∑_{k=0}^{n}(−1)^{k}\binom{k}{t}\binom{z + d − k}{k}\binom{z − 2k}{n − k}

for nonnegative integer values of t.

Let us now consider a bit more generally sums like

∑_{k=0}^{n}(−1)^{k}\binom{k}{t}\binom{x + λ k}{k}\binom{y + μ (n − k)}{n − k}  (8)

Note that the above sums show up when specializing:

x = z + d , y = z − 2n , λ = −1 , μ = 2

For t > n it is clear, that the sum vanishes, since every summand equals zero. For 0 < t < n, the first t terms of (8) vanish, and this together with rewriting of binomial coefficients and linear transformation of the variables makes it possible to replace the sum by a similar one in which the range 0 … n has been replaced by 0 … n − t. (Some care has to be taken to see that the inductive argument goes through in this situation). The critical case which remains to be treated and which can be eventually be made responsible for the cancellation phenomenon is then the case t = 0. Here we can use:

Lemma 4

∑_{k=0}^{n}(−1)^{k}\binom{x + λk}{k}\binom{y + μ(n − k)}{n − k}

can be written as a polynomial in (x − y), provided λ + μ = 1.

In particular, if we replace x by z + a and similarly y by z + b, for some a, b independent of z and k, the the sum is a constant w.r.t. z.

Let us now take a tridiagonal matrix G associated with polynomials of second degree, i.e. G = D_{Γ} where Γ = (g_{i}(z))_{i∈Z} s. th.

\begin{aligned}g_{−1}(z) & = & a_{0} + a_{1} z + a_{2} z^{2} \\ g_{0}(z) & = & b_{0} + b_{1} z + b_{2} z^{2} \\ g_{1}(z) & = & c_{0} + c_{1} z + c_{2} z^{2} \\ g_{k}(z) & = & 0 if |k| > 1 \end{aligned}  (9)

Let Λ = (h_{d}(z))_{d∈Z} be the conjugate family obtained from Γ via reduced Legendre transform. We know from the above considerations, that for d ≥ −1 the polynomials h_{d}(z) are polynomials of degree 2 if we take the coefficients a_{0}, …, c_{2} as indeterminates. What we are interested in is the question: under which conditions will the matrix D_{Λ} again be a tridiagonal matrix, (i.e. h_{d}(z) = 0 for |d| > 1) ? It should be clear that the occurence of this phenomenon is expressible as a linear condition in the indeterminates. So let us use the above formulas and calculate the first few values of h_{d}(z) = λ_{d,0} + λ_{d,1} z + λ_{d,2} z^{2}

d  λ_{d,0}  λ_{d,1}  λ_{d,2}
−1  a_{0}  a_{1}  a_{2}
0  −2a_{0} + a_{1} − a_{2} −4a_{1} + 4a_{2} + b_{1}  −6a_{2} + b_{2}
1 a_{0} − 5a_{1} + 11a_{2}  7a_{1} − 26a_{2}  17a_{2} − 4b_{2} + c_{2}
+b_{1} − b_{2} + c_{0} −2b_{1} + 4b_{2} + c_{1}
2 12a_{1} − 54a_{2} −8a_{1} + 48a_{2} −32a_{2} + 8b_{2} − 2c_{2}
−3b_{1} + 9b_{2}  +2b_{1} − 18b_{2}
2c_{0} + c_{1} − c_{2}  +4c_{2}
3 −20a_{1} + 170a_{2} 8a_{1} − 188a_{2} 48a_{2} − 12b_{2} + 3c_{2}
+5b_{1} − 35b_{2}  −2b_{1} + 44b_{2}
+3c_{0} − c_{1} + 7c_{2}  c_{1} − 10c_{2}
4  28a_{1} − 406a_{2} −8a_{1} + 340a_{2} −64a_{2} + 16b_{2} − 4c_{2}
−7b_{1} + 91b_{2}  +2b_{1} − 82b_{2}
+4c_{0} + 2c_{1} − 20c_{2}  +20c_{2}
5 −36a_{1} + 810a_{2} 8a_{1} − 540a_{2} 80a_{2} − 20b_{2} + 5c_{2}
+9b_{1} − 189b_{2}  −2b_{1} + 132b_{2}
+5c_{0} − 2c_{1}  c_{1} − 32c_{2}

Either via manipulation of formulas (tedious) or interpolation from the first few values we can find the following explicit general form for the polynomials h_{d}(z) as functions of the indeterminates a_{0}, …, c_{2}:

h_{d}(z) = (−1)^{d} \{[−16a_{2} + 4b_{2} − c_{2}] d z^{2} + + [24a_{2} − 6b_{2} + \frac{3}{2}c_{2}] d^{2} z

+ [−16a_{2} + 4b_{2} − c_{2}] d z + [−8a_{1} + 20a_{2} + 2b_{1} − 2b_{2}] z

+ [−8a_{2} + 2b_{2} − \frac{1}{2}c_{2}] d^{3} + [12a_{2} − 3b_{2} + \frac{4}{2}c_{2}] d^{2}

+ [8a_{1} − 24a_{2} − 2b_{1} + 3b_{2} + (−1)^{d}c_{0} + \frac{1}{2}c_{1}] d

+ [−4a_{1} + 10a_{2} + b_{1} − b_{2}]\} + \frac{1 − (−1)^{d}}{2} [c_{1} + \frac{1}{2}c_{2}] (z − \frac{1}{2})

From this we conclude that for h(d) to vanish identically for d ≥ 2 it is necessary and sufficient to have

\begin{aligned}c_{2} & = & −16a_{2} + 4b_{2} \\ c_{1} & = & −\frac{1}{2}c_{2} = 8a_{2} − 2b_{2} \\ c_{0} & = & 0 \\ b_{1} & = & 4a_{1} + 10a_{2} + b_{2} \end{aligned}  (10)

and the parameters a_{0}, a_{1}, a_{2}, b_{0}, b_{2} may be chosen independently. (Actually, rank consideration show that the vanishing of h_{2}(z) and h_{3}(z) is sufficient to impose these conditions, whereas vanishing of h_{2}(z) or h_{3}(z) alone is not).

Let us note that under the mentioned conditions the polynomials h_{−1}(z), h_{0}(z) and h_{1}(z) turn out to be:

\begin{aligned}h_{−1}(z) & = & g_{−1}(z) = a_{1} + a_{1} z + a_{2} z^{2} \\ h_{0}(z) & = & (−6a_{2} + b_{2})z^{2} + (−6a_{2} + b_{2})z − 2a_{0} + a_{1} + a_{2} + b_{0} \\ & = & (−6a_{2})z^{2} + (4a_{2} − 4a_{1})z − 2a_{0} + a_{1} − a_{2} + g_{0}(z) \\ h_{1}(z) & = & a_{2} z^{2} + (2a_{2} − a_{1})z + a_{0} − a_{1} + a_{2} \\ & = & a_{2}(z + 1)^{2} − a_{1}(z + 1) + a_{0} = g_{1}(−z − 1) \end{aligned}  (11)

Thus we have shown:

Theorem 5 For any tridiagonal matrix G = D_{Γ} as specified in (9), the matrix H = D_{Λ} obtained via reduced Legendre transform, i.e.

H = Q · G · Q^{−1}

is again tridiagonal if and only if the conditions in (10) are satisfied. In this case, the family Λ is given by (11).

Note that a verification of the if-part of this result is just as easy as the verificational proof of identity (4). Again, we can verify

H · Q = Q · G

by comparing the (i, j)-entries of the matrices on both sides with i and j as indeterminates. The same remarks as made at the end of section 3 apply.

The result of theorem 5 can be rephrased in terms of the original Legendre transform. Let

G̃ = Δ^{−1} · G · Δ

with Δ = diag (\binom{2k}{k})_{k≥0}. Then G̃ = DΓ̃ where Γ̃ = (g̃_{k}(z))_{k∈Z} is given by

g̃_{k}(z) = \{ \begin{matrix}\frac{4z+2}{z+1} g_{−1}(z) & if k = −1 \\ g_{0}(z) & if k = 0 \\ \frac{z}{4z−2} g_{1}(z) & if k = 1 \\ 0 & else \end{matrix}

If we impose conditions (10) and furthermore require that z + 1 | g_{−1}(z) (which introduces a further linear constraint: a_{0} = a_{1} − a_{2}), then we find:

Corollary 6 Legendre transform

H = P · G̃ · P^{−1}

maps the polynomials

g̃_{−1}(z) = A z^{2} + B z + C
g̃_{0}(z) = (A + D) z^{2} + (−A + B + D) z + E
g̃_{1}(z) = D z^{2}

onto

h_{−1}(z) = \frac{z + 1}{4z + 2} (A z^{2} + B z + c)
h_{0}(z) = \frac{1}{2} [(2D − A) z(z + 1) + 2E − C]
h_{1}(z) = \frac{z}{4z + 2} [A z^{2} + (2A − B) z + (A − B + C)]

Note that the translation is given by A = 4a_{2}, B = −2a_{2} + 4a_{1}, C = 2(a_{1} − a_{2}), D = c_{2}/4 and E = b_{0}.

This is essentially the contents^{10} of theorem 1 in [16], where a rather different proof using creative telescoping is given. Note that the situation of the Franel-Apéry conjugacy is given by the particular values A = C = 1, B = 2, D = −8, E = 2. Schmidt’s article contains a number of further examples illustrating this transformation.

5 A brief look at diophantine approximation

In this final section I would like to draw attention to an application of the foregoing to diophantine approximation. The results mentioned here are adapted from [16].

As mentioned in the beginning, Apéry’s surprising proof of the irrationality of ζ(3) is based on the recurrence relation for the Apéry numbers. In a few words, his idea of proof goes as follows (see [20]):

10The equivalence is not complete, because the additional linear relation introduced here reduces the number of free parameters by one. But this is a technical matter that can be dealt with.

Consider the two sequences

(a_{n})_{n≥0} = ( 1 , 5 , 73 , 1445 , 33001 , …)
(b_{n})_{n≥0} = ( 0 , 6 , \frac{351}{4} , \frac{62531}{36} , \frac{11424695}{288} , …)

satisfying both Apéry’s linear recurrence relation

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

The surprising fact about (a_{n})_{n≥0} is that all its terms are integers, which is contained in (3), but by no means obvious from (12) alone. The terms of (b_{n})_{n≥0} are rational numbers, where the denominator of b_{n} divides 2 lcm(1, 2, …, n)^{3}. This follows from

b_{n} = ∑_{k=0}^{n} \binom{n}{k}^{2}\binom{n + k}{k}^{2} c_{n,k}

where

c_{n,k} = ∑_{m=1}^{n} \frac{1}{m^{3}} + ∑_{m=1}^{k} \frac{(−1)^{m−1}}{2 m^{3}\binom{n}{m}\binom{n+m}{m}}

and this also shows that

lim_{n→∞} \frac{b_{n}}{a_{n}} = ζ(3)

The joint recurrence (12) for numerator and denominator of this recurrence allows for an estimate of the rate of convergence:

ζ(3) − \frac{b_{n}}{a_{n}} = ∑_{k=n+1}^{∞} \frac{a_{k} b_{k+1} − a_{k+1} b_{k}}{a_{k} a_{k+1}} = O (\frac{1}{a_{n}^{2}} )

since

a_{n−1} b_{n} − a_{n} b_{n−1} = \frac{a_{0} b_{1} − a_{1} b_{0}}{n^{3}} = \frac{6}{n^{3}}

The sequences (p_{n})_{n≥0} and (q_{n})_{n≥0} with

p_{n} = 2 lcm(1, 2, …, n)^{3} b_{n} and q_{n} = 2 lcm(1, 2, …, n)^{3} a_{n} (n ≥ 0)

respectively, are integer sequences with

lim \frac{p_{n}}{q_{n}} = lim \frac{b_{n}}{a_{n}} = ζ(3)

and from an estimate of the growth rate of (a_{n})_{n≥0} (which is obtained from co-

efficients of the recurrence (12)) one can show that the convergence of (\frac{p_{n}}{q_{n}} )_{n≥0}

towards ζ(3) is rapidly enough (in terms of the growth of the q_{n}) so as to establish non-rationality of the limit.

In terms of our matrix notation, writing A for the matrix representing the Apéry recurrence, a for the (column)-vector of the a_{n} (n ≥ 0) and similarly for b, we have

A · a = 0 and A · b = \left[{6 \atop 0}\right]

where on the right we have a vector with a 0 in all positions, except a 6 in the top position.

As to the Franel recurrence, it has been remarked by Cusick (see a footnote in [20]) that both sequences

(f_{n})_{n≥0} = (1, 2, 10, 56, 346, …)
(g_{n})_{n≥0} = (0, 3, 12, \frac{208}{3}, \frac{1280}{3}, …)

satisfying the Franel recurrence, which we may write as

F · f = 0 and F · g = \left[{3 \atop 0}\right]

lead to a diophantine approximation:

lim \frac{g_{n}}{f_{n}} = \frac{π^{2}}{8}  (13)

We have seen above that the equation

A · a = 0 and F · f = 0

are related via Legendre conjugacy. Now conjugacy provides us with a solution of an inhomogeneous Franel recurrence associated to the second solution b of the (homogeneous) Apéry recurrence:

A · b = \left[{6 \atop 0}\right] vs. F · h = k

where h and k are vectors belonging to the sequences

(h_{n})_{n≥0} = ( 0 , 3 , \frac{93}{8} , \frac{1217}{18} , \frac{239429}{576} , …)
(k_{n})_{n≥0} = ( 3 , −\frac{3}{2} , 1 , −\frac{4}{3} , …, (−1)^{n} \frac{3}{n + 1} , …)

Similarly, the second solution g of the (homogeneous) Franel recurrence provides us with a solution of an inhomogeneous Apéry recurrence:

F · g = \left[{3 \atop 0}\right] vs. A · c = d

where c and d are vectors belonging to the sequences

(c_{n})_{n≥0} = ( 0 , 6 , 90 , \frac{5348}{3} , \frac{122140}{3} …)
(d_{n})_{n≥0} = ( 6 , 18 , 30 , …, 3 (4n + 2) , …)

As a result, we get the approximations

lim_{n→∞} \frac{h_{n}}{f_{n}} = ζ(3) and lim_{n→∞} \frac{c_{n}}{a_{n}} = \frac{π^{2}}{8}

i.e. one approximation of ζ(3) in terms of Franel-recursive sequences, and one of π^{2}/8 in terms of Apéry-recursive sequences. From this, one may hope for a proof of independence of { 1 , π^{2}/8 , ζ(3) } over the rationals.

References

[1] R. Askey and J. A. Wilson. A recurrence relation generalizing those of Apéry. J. Austral. Math. Soc. (Series A), 36:267–278, 1984.

[2] W. N. Bailey. The generating function for Jacobi polynomials. J. London Math. Soc., 13:243–246, 1938.

[3] P. Cartier. Démonstration ”automatique” d’identités et fonctions hypergéometriques [d’après D. Zeilberger]. In Séminaire Bourbaki, volume 746, 1991.

[4] L. Comtet. Advanced Combinatorics. Reidel, Dordrecht, 1974.

[5] T. Cusick. Recurrences for sums of powers of binomial coefficients. J. Comb. Theory (A), 52:77–83, 1989.

[6] J. Franel. notes published in: L’Intermédiaire des Mathématiciens, vol. 1 (1894), 45-47, and vol. 2 (1895), 33-35.

[7] R. W. Gosper. Decision procedure for indefinite hypergeometric summation. Proc. Nat. Acad. Sci. USA, 75:40–42, 1978.

[8] H. W. Gould. Combinatorial Identities: A Standardized Set of Tables Listing 500 Binomial Coefficient Identities (rev. ed.). Morgantown, W. VA., 1972.

[9] J. Hornegger. Hypergeometrische Summation und polynomiale Rekursion. Diplomarbeit, IMMD-I, Erlangen, 1992.

[10] T. H. Koornwinder. On Zeilberger’s algorithm and its q-analogue: a rigorous description. Technical Report AM-R9207, Centre for Mathematics and Computer Science, Amsterdam, 1992.

[11] M. A. Perlstadt. Some recurrences for sums of powers of binomial coefficients. J. Number Theory, 27:304–309, 1987.

[12] J. Riordan. Combinatorial Identities. John Wiley & Sons, Inc., New York-London-Sidney, 1968.

[13] B. Salvy and P. Zimmerman. Gfun: a Maple package for the manipulation of generating and holonomic functions in one variable. preprint, October 1992.

[14] A. L. Schmidt. Generalized Legendre polynomials. J. Reine Angew. Math., 404:192–202, 1990.

[15] A. L. Schmidt. Generalized Legendre and q-Legendre polynomials. In VII Simposium sobre polinomios ortogonales y applicaciones, Granada 1991, 1991. (to appear in: J. Comp. Appl. Math.).

[16] A. L. Schmidt. Legendre transforms and Apéry’s sequences. preprint series 14, Matematisk Institut, University of Copenhagen, August 1992.

[17] D. Stanton. A short proof of a generating function for the Jacobi polynomials. Proc. Amer. Math. Soc., 80:398–400, 1980.

[18] V. Strehl. Binomial identities: combinatorial and algorithmic aspects. (in preparation).

[19] V. Strehl. Zykel-Enumeration bei lokal-strukturierten Funktionen. Habilitationsschrift, Erlangen, 1990.

[20] A. J. van der Poorten. A proof that Euler missed ... Apéry’s proof of the irrationality of ζ(3). Math. Intelligencer, 1:195–203, 1978.

[21] H. S. Wilf and D. Zeilberger. Rational functions certify combinatorial identities. J. Amer. Math. Soc., 3:147–158, 1990.

[22] H. S. Wilf and D. Zeilberger. Towards computerized proofs of combinatorial identities. Bull. Amer. Math. Soc. (New Series), pages 77–83, 1990.

[23] H. S. Wilf and D. Zeilberger. An algorithmic proof theory for hypergeometric (ordinary and ”q”) multisum/integral identities. Invent. Math., 108:575–633, 1992.

[24] H. S. Wilf and D. Zeilberger. Rational function certification of multisum/integral/”q” identites. Bull. Amer. Math. Soc. (New Series), 27:148–153, 1992.

[25] D. Zeilberger. A fast algorithm for proving terminating hypergeometric identities. Discrete Mathematics, 80:207–211, 1990.

[26] D. Zeilberger. The method of creative telescoping. J. Symb. Comp., 11:195– 205, 1991.