HYPERG : a Maple package for manipulating hypergeometric series

Bruno Gauthier ^{∗}

Abstract This article describes a Maple package for the manipulation of generalized hyper-

geometric series. This package provides tools for automatic proofs of identities. The aim of this article is to present the package and to give some results obtained with it. It can be obtained by either e-mail (Bruno.Gauthier@univ-mlv.fr) or on the web

  site (http://www-igm.univ-mlv.fr/~gauthier/HYPERG.html).

Introduction

HYPERG is a Maple package for the manipulation and identification of hypergeometric series and identities. It provides several tools for:

• transforming sums into hypergeometric notation ;

• manipulating factorial expressions ;

• applying summation and transformation formulas ;

• applying contiguous relations ;

• computing formal limits of hypergeometric expressions ;

• manipulating linear recurrence equations ;

• finding polynomial, rational and hypergeometric solutions of recurrences ;

• using the Gosper and Zeilberger algorithms. In the following sections, we describe Maple basic objects, principal features and applica-

tions.

∗Bruno.Gauthier@univ-mlv.fr, Université Marne-la-Vallée, Institut Gaspard Monge, Cité Descartes, 5 bd Descartes, Champs sur marne, 77454 Marne-la-Vallée Cedex 2, France.

Basic objects

First, we recall that a sequence (U_{k}) is hypergeometric if the ratio of two consecutive terms is a rational function of the summation index k:

\frac{U_{k+1}}{U_{k}} = \frac{p(k)}{q(k)} ,  (1)

where p and q are polynomials in k. Then we define two new hypergeometric objects in Maple: the rising factorial (also known as Pochhammer symbol) RF[a,k]:

(a)_{k} := \{ \begin{matrix}a (a + 1) … (a + k − 1) & k > 0, \\ 1 & k = 0, \\ [(a − 1) (a − 2) … (a + k)]^{−1} & k < 0, \end{matrix}

and also the hypergeometric series

HYP[[a_1,...,a_r],[b_1,...,b_s],z]:

rF_{s} [{a_{1}, …, a_{r} \atop b_{1}, …, b_{s}}; z] = ∑_{k=0}^{∞} \frac{(a_{1})_{k} … (a_{r})_{k}}{(b_{1})_{k} … (b_{s})_{k}} \frac{z^{k}}{k!} .

The hypergeometric series were already defined in the Maple computer algebra system by

hypergeom([a_1,...,a_r],[b_1,...,b_s],z), but we prefered redefine them to avoid Maple

automatic simplifications.

It is important to recognize when a given series is hypergeometric, because the general theory of hypergeometric functions is very powerful. So, what we want to do is to show how one may identify a given hypergeometric series with a particular _{r}F_{s} :

> SumToHyp( Sum(x^k/k!,k) ); # exponential series
                                 HYP[[], [], x]
> SumToHyp( Sum(1/(2*k-1)/(2*k+1)!,k=0..infinity) );
                         - HYP[[-1/2], [3/2, 1/2], 1/4]
> SumToHyp( Sum(4*k^2+1,k) ); # hypergeometric series with complex parameters
              HYP[[1 - 1/2 I, 1 + 1/2 I, 1], [- 1/2 I, 1/2 I], 1]
> SumToHyp( Sum(k^2+k!,k) );
Error, (in SumToHyp) Not a hypergeometric series

Database of formulas

This package contains a lot of summation, transformation and contiguous formulas involving hypergeometric series. (The numbering of these formulas is similar as for the Krattenthaler’s Mathematica package HYP [6].) Here, we give examples showing how these classical hypergeometric formulas can be used.

From Watson’s theorem, we derive Gauss’s second theorem :

> eg := HypSumPrint( HYP[[a,b,c],[d,e],z], 33 );
# z must be equal to 1
# 1/2+1/2*a+1/2*b must be equal to d
# 2*c must be equal to e
eg := HYP[[a, b, c], [d, e], z] =
   1/2
 Pi    GAMMA(1/2 + c) GAMMA(1/2 + 1/2 a + 1/2 b) GAMMA(1/2 - 1/2 a - 1/2 b + c)
    /
   /  (GAMMA(1/2 + 1/2 a) GAMMA(1/2 + 1/2 b) GAMMA(1/2 - 1/2 a + c)
  /
 GAMMA(1/2 - 1/2 b + c))
> eg := subs( z=1, d=1/2+1/2*a+1/2*b, e=2*c, eg );
eg := HYP[[a, b, c], [1/2 + 1/2 a + 1/2 b, 2 c], 1] =
   1/2
 Pi    GAMMA(1/2 + c) GAMMA(1/2 + 1/2 a + 1/2 b) GAMMA(1/2 - 1/2 a - 1/2 b + c)
    /
   /  (GAMMA(1/2 + 1/2 a) GAMMA(1/2 + 1/2 b) GAMMA(1/2 - 1/2 a + c)
  /
 GAMMA(1/2 - 1/2 b + c))
> eg := Lim( eg, c=infinity );
eg :=
                                               1/2
                                             Pi    GAMMA(1/2 + 1/2 a + 1/2 b)
HYP[[a, b], [1/2 + 1/2 a + 1/2 b], 1/2] = -------------------------------------
                                          GAMMA(1/2 + 1/2 a) GAMMA(1/2 + 1/2 b)

Ramanujan hypergeometric identity

1 − (\frac{1}{2})^{3} + (\frac{1/2 3/2}{2})^{3} + ⋯

= _{3}F_{2} [{1/2, 1/2, 1/2 \atop 1, 1}; −1]

= \frac{Γ(\frac{9}{8})^{2}}{Γ(\frac{5}{4})^{2} Γ(\frac{7}{8})^{2}} .

> H := HYP[[1/2,1/2,1/2],[1,1],-1]:

> Ht := HypTransf( H, 34 );
                                                      2
                        Ht := HYP[[1/4, 1/4], [1], -1]
> HypSum( Ht, 4 ); # Kummer's theorem
                                     1/2            2
                               (2 + 2   ) GAMMA(3/4)
                           1/4 ----------------------
                                               4
                                     GAMMA(7/8)

Formal limits

This section deals with formal limits of hypergeometrc expression.

> #obviously
> expr := RF[-n,r]/RF[1+a+b-c-n,r];
                                       RF[- n, r]
                        expr := ------------------------
                                RF[1 + a + b - c - n, r]
> Lim( expr, n=infinity );
                                       1
> #confluent principle
> Lim( HYP[[a,2*b],[c],z/b], b=infinity );
                              HYP[[a], [c], 2 z]
> # From Wilson polynomials, we get Continuous dual Hahn polynomials
> W(n):= HYP[[-n,n+a+b+c+d-1,a+I*x,a-I*x],[a+b,a+c,a+d],1]
         * RF[a+b,n] * RF[a+c,n] * RF[a+d,n]:
> Lim( W(n)/RF[a+d,n], d=infinity );
GAMMA(a + b + n) GAMMA(a + c + n)
    HYP[[- n, a + I x, a - I x], [a + b, a + c], 1]/(GAMMA(a + b) GAMMA(a + c))

Manipulating recurrence equations

The package HYPERG is fully compatible with gfun [9], therefor we need to add gfun’s features. In fact we are interested in the tools handling recurrence equations and differential equations. We state the existence of the rectodiffeq and diffeqtorec functions which compute the recurrence associated with the differential equation and vice versa. For instance:

> with(gfun):
> rectodiffeq( {(5*n+10)*u(n)+k*u(n+1)-u(n+2),u(0)=0,u(1)=0}, u(n), f(t) );

                           2                      3
                    (- 10 t  - k t + 1) f(t) - 5 t  D(f)(t)

The others important tools are closure operations such as:

  borel            cauchyproduct    invborel         rec+rec          rec*rec
  diffeq+diffeq    diffeq*diffeq    hadamardproduct

In this example we biuld from holonomic recurrences of 2^{n} and n! the recurrence satisfied by 2^{n} + n! using two methodes : closure by addition and the computation of the following casoratian determinant

| \begin{matrix}y(n) & n! & 2^{n} \\ y(n + 1) & (n + 1)! & 2^{n+1} \\ y(n + 2) & (n + 2)! & 2^{n+2} \end{matrix} | .

> with(HYPERG):
> with(gfun):
> rec1 := HypergToRec( 2^n, U(n) );
                         rec1 := 2 U(n) - U(n + 1) = 0
> rec2 := HypergToRec( n!, U(n) );
                      rec2 := U(n) (n + 1) - U(n + 1) = 0
> `rec+rec`( rec1, rec2, U(n) ); # closure by addition
        2                             2
   {(2 n  + 2 n) U(n) + (- 3 n + 2 - n ) U(n + 1) + (n - 1) U(n + 2),
       U(0) = _C[0], U(1) = _C[1], U(2) = 2 _C[1], U(3) = 4 _C[0] + 2 _C[1]}
> GenRec( {2^n,n!}, U(n) );      # casoratian determinant
                                       2
      2 n (n + 1) U(n) + (- 3 n + 2 - n ) U(n + 1) + (n - 1) U(n + 2) = 0

Solving recurrence relation

It is possible to manipulate and find polynomial [8], rational [1] and hypergeometric [7] solutions of linear recurrence relation

> PolySolRec( 3*U(n+2)-n*U(n+1)+(n-1)*U(n)=0, U(n) );
                                                2
                              _x2 (27 - 11 n + n )
> PolySolRec( r(n+2)-2*r(n+1)+r(n)=2, r(n), {r(0)=-1,r(1)=m} );
                                              2
                                 - 1 + m n + n
> RatioSolRec( (n+4)*(2*n+1)*(n+2)*U(n+3)-(2*n+3)*(n+3)*(n+1)*U(n+2)

>           +n*(n+2)*(2*n-3)*U(n+1)-(n-1)*(2*n-1)*(n+1)*U(n)=0, U(n) );
                                   _x2 (2 n - 3)
                              1/2 ---------------
                                  (n + 1) (n - 1)

The following HypSolRec function doesn’t return the hypergeometric solution but the ratio of two consecutive terms:

> HypSolRec( (n+4)*a(n+2)+a(n+1)-(n+1)*a(n)=0,a(n) );
                             (n + 1) (2 n + 5)  n + 1
                          {- -----------------, -----}
                             (n + 3) (2 n + 3)  n + 3

Here, we verify that the recurrence of the previous section is satisfied by 2^{n} and n!:

> HypSolRec( 2*n*(n+1)*U(n)+(-n^2-3*n+2)*U(n+1)+(n-1)*U(n+2) = 0, U(n) );
                                   {2, n + 1}

Finding a recurrence

Gosper algorithm

The Gosper algorithm [4] deals with the question of finding an upward antidifference T (k) for given U(k), i.e. a sequence T (k) for which

T (k + 1) − T (k) = U(k) ,

in the particular case when U(k) is a hypergeometric term.

> T := Gosper( k^2, k );
                                               2        3
                       T := _x0 + 1/6 k - 1/2 k  + 1/3 k
> U := subs(k=k+1,T) - T;
                                  2              3        2        3
            U := 1/6 - 1/2 (k + 1)  + 1/3 (k + 1)  + 1/2 k  - 1/3 k
> simplify( U );
                                        2
                                       k
> Gosper( k*k!, k );
                                       k!
> Gosper( (k/2)!, k );
Error, (in Gosper) (1/2*k)!, is not hypergeometric in, k

Remark: Gosper’s algorithm is in the computational heart of Zeilberger’s fast algorithm described in the next section.

Zeilberger algorithm

The (fast) Zeilberger algorithm [11, 12] deals with the definite summation of hypergeometric terms. It is a direct application of the Gosper algorithm. It finds a linear homogeneous recurrence relation for the sum

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

where F (n, k) is a hypergeometric term in n and k.

The first step is finding a recurrence relation for the summand F (n, k). It is an extension of the Gosper algorithm which finds polynomials a_{j}(n) and G(n, k) such that:

∑_{j=0}^{J} a_{j}(n) F (n + j, k) = G(n, k + 1) − G(n, k) ,  (2)

ith G(n, k) = R(n, k) F (n, k) where R(n, k) is a rational function called “certificate”. Summing (2) over all integers k, we get by creative telescoping [13]

∑_{j=0}^{J} a_{j}(n) s(n + j) = G(n, ∞) − G(n, −∞) ,  (3)

and checking that G(n, ±∞) = 0, the recurrence over the initial sum s(n) is proved:

a_{0}(n) s(n) + a_{1}(n) s(n + 1) + ⋯ + a_{J}(n) s(n + J) = 0 .  (4)

> SummandToRec( binomial(n,k), k, n, F, G, 'cert' );
               - 2 F(n, k) + F(n + 1, k) = G(n, k + 1) - G(n, k)
> cert;
                                       k
                                  -----------
                                  - n - 1 + k

Applications

Sum of binomial coefficients

This is a first trivial example which evals the sum of binomial coefficients:

> rec := SumToRec( Sum(binomial(n,k),k), n, S );
                         rec := - 2 S(n) + S(n + 1) = 0
> HypSolRec( rec, S(n) );
                                      {2}

Terminating hypergeometric series

Following suggestions of Koornwinder [5], we give a new proof of Kummer’s theorem involving Zeilberger’s algorithm, the WZ method and asymptotic estimates.

In the following Kummer’s theorem

2F_{1} [{a, b \atop 1 + a − b}; −1] = \frac{Γ(1 + \frac{a}{2}) Γ(1 + a − b)}{Γ(1 + a) Γ(1 + \frac{a}{2} − b)} ,  (5)

let us replace a by a + 2 n where n is a new variable that stands for a nonnegative integer. The resulting formula

2F_{1} [{a + 2 n, b \atop 1 + a + 2 n − b}; −1] = \frac{Γ(1 + \frac{a}{2} + n) Γ(1 + a + 2 n − b)}{Γ(1 + a + 2 n) Γ(1 + \frac{a}{2} + n − b)}  (6)

is clearly equivalent to (5). The introduction of the free parameter n is the key to this new computer-assisted proof.

Formula (6) can be written as:

\frac{∑_{k} f(n, k)}{S(n)} = 1 ,  (7)

where

f(n, k) = \frac{(a + 2 n)_{k} (b)_{k}}{(1 + a + 2 n − b)_{k}} \frac{(−1)^{k}}{k!}

and

S(n) = \frac{Γ(1 + \frac{a}{2} + n) Γ(1 + a + 2 n − b)}{Γ(1 + a + 2 n) Γ(1 + \frac{a}{2} + n − b)} .

(The coefficient 2 in the substitution a → a + 2 n makes S(n) hypergeometric in n.) Formula (7) can be rewritten as:

∑_{k} F (n, k) = 1  where F (n, k) = \frac{f(n, k)}{S(n)}.  (8)

We shall actually prove this last formula. Using HYPERG, we apply Zeilberger’s algorithm to obtain:

F (n, k) − F (n + 1, k) = G(n, k + 1) − G(n, k) ,  (9)

with G(n, k) = F (n, k) C(n, k) and C(n, k) = −\frac{(b − 1) k}{(1 + a + 2 n − b + k) (a + 2 n)}.

This C(n, k) is the so-called certificate of the WZ-pair (F, G). After summing (9) for k from 0 to K, the right-hand side telescopes:

∑_{k=0}^{K} F (n, k) − ∑_{k=0}^{K} F (n + 1, k) = ∑_{k=0}^{K} G(n, k + 1) − ∑_{k=0}^{K} G(n, k)  (10)
= G(n, K + 1) − G(n, 0)  (11)
= F (n, K + 1) C(n, K + 1) .  (12)

We now let K tend to infinity. Clearly,

lim_{K→+∞} C(n, K + 1) = −\frac{b − 1}{a + 2 n} .

To determine lim_{K→+∞} F (n, K + 1), we use the well-known estimate (see e.g. [5]):

\frac{Γ(a + k)}{Γ(b + k)} ∼ k^{a−b} as k → ∞ .  (13)

First,

F (n, k) = \frac{(a + 2 n)_{k} (b)_{k}}{(1 + a + 2 n − b)_{k}} \frac{(−1)^{k}}{k!} \frac{1}{S(n)} ,
= \frac{Γ(a + 2 n + k) Γ(b + k) (−1)^{k}}{Γ(1 + a + 2 n − b + k) Γ(k + 1)} T (n) ,

where T (n) = \frac{(a + 2 n) Γ(1 + \frac{a}{2} + n − b)}{Γ(b) Γ(1 + \frac{a}{2} + n)} does not depend on k.

Then,

|F (n, k)| ∼ k^{a+2n−(1+a+2n−b)} k^{b−1} T (n) ,
∼ k^{2b−2} T (n) .

We know that Re(b) < 1 (necessary condition for the convergence of the series in (5)). Therefore

lim_{K→+∞} F (n, K + 1) = 0 .

By a dominated convergence argument [3], we can let K to infinity in (12) to obtain :

∑_{k=0}^{∞} F (n, k) − ∑_{k=0}^{∞} F (n + 1, k) = 0 .

Hence, the sum ∑_{k=0}^{∞} F (n, k) is independent of n. Then

∑_{k=0}^{∞} F (n, k) = lim_{n→+∞} ∑_{k=0}^{∞} F (n, k) ,
= ∑_{k=0}^{∞} lim_{n→+∞} F (n, k) ,

where the last equality is justified by dominated convergence. Moreover,

∑_{k=0}^{∞} F (n, k) = ∑_{k=0}^{∞} lim_{n→+∞}(F (n, k) ,
= ∑_{k=0}^{∞} 2^{b} \frac{(b)_{k} (−1)^{k}}{(k)!}
= 2^{b} ∑_{k=0}^{∞} \frac{(b)_{k} (−1)^{k}}{k!} ,
= 2^{b} (1 − (−1))^{b} = 1 ,  by binomial theorem.

This is precisely formula (8).

References

[1] S. A. Abramov, Rational solutions of linear differential and q-difference equations with polynomial coefficients, Proceedings ISSAC ’95, ACM Press, New York, (1995) 285-289.

[2] B. Gauthier, HYPERG, Maple package, User’s reference manual, Version 1.0, (182 pages), version électronique disponible à l’URL

http://www-igm.univ-mlv.fr/~gauthier, (1999).

[3] B. Gauthier, A proof of Kummer’s theorem, e-Print archive, math.CA/9904061, (1999) 10 pages.

[4] R. W. Gosper, Jr, Decision procedure for indefinite hypergeometric summation, Proceedings of the national academy of sciences of the United States of America, 75 (1978) 40–42.

[5] T. H. Koornwinder, Identities of nonterminating series by Zeilberger’s algorithm, Journal of computational and applied mathematics, 99 (1998) 449–461.

[6] C. Krattenthaler, HYP and HYPQ, Mathematica packages for the manipulation of binomial sums and hypergeometric series, respectively q-binomial sums and basic hypergeometric series, Journal of symbolic computation, 20 (1995) 737– 744.

[7] M. Petkovšek, Hypergeometric solutions of linear recurrences with polynomial coefficients, Journal of symbolic computation, 14 (1992) 243–264.

[8] M. Petkovšek, H. S. Wilf, D. Zeilberger, A = B, A.K. Peters, Wellesley, Massachusetts, (1996). (Les programmes sont disponibles sur le site WWW

http://www.cis.upenn.edu/~wilf/AeqB.html).

[9] B. Salvy, P. Zimmermann, Gfun : a maple package for the manipulation of generating and holonomic functions in one variable, ACM Transactions on Mathematical Software, 20, 2 (1994) 163–177.

[10] L. J. Slater, Generalized hypergeometric functions, Cambridge University Press, Cambridge, (1966).

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

[12] D. Zeilberger, Three recitations on holonomic systems and hypergeometric series, Séminaire Lotharingien de Combinatoire, B24a (1990).

[13] D. Zeilberger, The method of creative telescoping, Journal of symbolic computation, 11 (1991) 195–204.