A Generalization of Gosper’s Algorithm to Bibasic Hypergeometric Summation

Axel Riese

Research Institute for Symbolic Computation Johannes Kepler University Linz A–4040 Linz, Austria

                  Axel.Riese@risc.uni-linz.ac.at

Submitted: May 9, 1996; Accepted: June 24, 1996.

Abstract An algebraically motivated generalization of Gosper’s algorithm to indefinite bibasic hypergeometric summation is presented. In particular, it is shown how Paule’s concept of greatest factorial factorization of polynomials can be extended to the bibasic case. It turns out that most of the bibasic hypergeometric summation identities from literature can be proved and even found this way. A Mathematica implementation of the algorithm is available from the author.

AMS Subject Classification. Primary 33D65, 68Q40; Secondary 33D20.

1 Introduction

Recently, Paule and Strehl [10] observed that the algorithm presented by Gosper [7] for indefinite hypergeometric summation extends quite naturally to the q-hypergeometric case by introducing a q-analogue of the canonical Gosper-Petkovšek (GP) representation for rational functions. Based on the new algebraic concept of greatest factorial factorization (GFF), Paule [8] developed an alternative but equivalent approach to hypergeometric telescoping. It was also shown by Paule (cf. Paule and Riese [9]) that the problem of q-hypergeometric telescoping can be treated along the same lines as the q = 1 case by making use of a q-version of GFF. Built on these concepts, a Mathematica implementation of q-analogues of Gosper’s as well as of Zeilberger’s [14] fast algorithm for definite q-hypergeometric summation has been carried out by the author (cf. Paule and Riese [9], and Riese [12]). The original approach to definite q-hypergeometric summation is due to Wilf and Zeilberger [13].

The object of this paper is to describe how the algorithm qTelescope presented in [9], a q-analogue of Gosper’s algorithm, generalizes to the bibasic hypergeometric case. In Section 2, the underlying theoretical background based on a bibasic extension of GFF is discussed, which leads to the bibasic counterpart of the algorithm qTelescope. In Section 3, the degree setting for solving the bibasic key equation is established. Applications are given in Section 4 to illustrate the usage of the newly developed Mathematica implementation which is available by

email request to Axel.Riese@risc.uni-linz.ac.at.

2 Theoretical Background

In this section, q-greatest factorial factorization (qGFF) of polynomials, which has been introduced by Paule (cf. Paule and Riese [9]) providing an algebraic explanation of q-hypergeometric telescoping, is extended to the bibasic hypergeometric case. It turns out that to this end the q-case argumentation can be carried over almost word by word.

2.1 Bibasic Greatest Factorial Factorization

Let ℤ denote the set of all integers, and ℕ the set of all non-negative integers. Let p, q, x, and y be fixed indeterminates. Assume K = L(κ_{1}, …, κ_{n}) to be the field of rational functions in a fixed number of indeterminates κ_{1}, …, κ_{n}, n ∈ ℕ, where p ≠ κ_{i} ≠ y and q ≠ κ_{i} ≠ x, 1 ≤ i ≤ n, over some computable field L of characteristic 0 and not containing p, q, x, and y. (For the sake of simplicity with regard to the implementation we will restrict ourselves to the case where L is the rational number field ℚ.) The transcendental extension of K by the indeterminates p and q is denoted by F, i.e., F = K(p, q).

For P ∈ F [x, y], let the bibasic shift operator ϵ be given by (ϵP)(x, y) = P (qx, py). The extension of this shift operator to the rational function field F (x, y), the quotient field of the polynomial ring F [x, y], will be also denoted by ϵ.

Definition 1. A polynomial P ∈ F [x] (resp. P ∈ F [y]) is called q-monic (resp. p-monic) if P (0) = 1. A polynomial P ∈ F [x, y] is called bibasic monic if P (x, 0) ≠ 0 ≠ P (0, y) and either P (0, 0) = 1, or P (0, 0) = 0 and the coefficients of P are relatively prime polynomials in F.^{†}

Example. (i) The following polynomials are bibasic monic:

P_{1}(x, y) = 1,  P_{2}(x, y) = 1 − apqx^{2}y^{3},  P_{3}(x, y) = (1 − q)^{2}x^{2} + py.

(ii) The following polynomials are not bibasic monic:

P_{4}(x, y) = q,  P_{5}(x, y) = xy − apqx^{2}y^{3},  P_{6}(x, y) = (1 − q)^{−1}px^{2} + py.

The properties of being q-monic, p-monic, and bibasic monic are clearly invariant with respect to the bibasic shift operator ϵ, i.e., if P is q-monic, p-monic, or bibasic monic, then the same holds true for ϵP. Furthermore, the product of two bibasic monic polynomials is again bibasic monic. Also note that a bibasic monic polynomial P satisfies gcd(x, P) = 1 = gcd(y, P).

Evidently, any non-zero polynomial P ∈ F [x, y] has a unique factorization, the bibasic monic decomposition, in the form

P = z · x^{α} · y^{β} · P^{∗},

where z ∈ F, α, β ∈ ℕ, and P^{∗} ∈ F [x, y] is bibasic monic.

The bibasic monic decomposition of a polynomial P ≠ 0 can be computed easily as follows. Define α := max{i ∈ ℕ : x^{i}|P }, β := max{j ∈ ℕ : y^{j}|P }, and put P̄ := x^{−α} · y^{−β} · P. If P̄(0, 0) ≠ 0 define z := P̄(0, 0), otherwise let l denote the least common multiple of all coefficient-denominators of P̄, let g denote the greatest common divisor of all coefficients of l · P̄, and define z := g/l. Then, for P^{∗} := z^{−1} · P̄, the bibasic monic decomposition of P is given by P = z · x^{α} · y^{β} · P^{∗}.

†In other words, P is assumed to be primitive over L[κ_{1}, …, κ_{n}, p, q] in this case, which will guarantee the uniqueness of the so-called bibasic monic decomposition of a polynomial as shown below.

Example. The bibasic monic decompositions of the polynomials P_{4}, P_{5}, and P_{6} from the example above are given by

P_{4} = q · x^{0} · y^{0} · 1, P_{5} = 1 · x · y · (1 − apqxy^{2}), P_{6} = \frac{p}{1 − q} · x^{0} · y^{0} · (x^{2} + (1 − q)y).

Moreover, we assume the result of any gcd computation over F [x, y] as being normalized in the following sense. If P_{1} = z_{1} · x^{α_{1}} · y^{β_{1}} · P_{1}^{∗} and P_{2} = z_{2} · x^{α_{2}} · y^{β_{2}} · P_{2}^{∗} are the bibasic monic decompositions of P_{1}, P_{2} ∈ F [x, y], we define

gcd_{p,q}(P_{1}, P_{2}) := gcd(x^{α_{1}}, x^{α_{2}}) · gcd(y^{β_{1}}, y^{β_{2}}) · gcd_{p,q}(P_{1}^{∗}, P_{2}^{∗}),

where the gcd_{p,q} of two bibasic monic polynomials is understood to be bibasic monic.

The polynomial degree in x and y of any P ∈ F [x, y] is denoted by deg_{x}(P) and deg_{y}(P), respectively.

Definition 2. For any bibasic monic polynomial P ∈ F [x, y] and k ∈ ℕ, the k-th falling bibasic factorial [P]\underline{k} ; p,q of P is defined as

[P]\underline{k} ; p,q := ∏_{i=0}^{k−1} ϵ^{−i}P.

Note that by the null convention ∏_{i∈∅} P_{i} := 1 we have [P]\underline{0} ; p,q = 1. In general, polynomials arising in bibasic hypergeometric summation have several different representations in terms of falling bibasic factorials. From all possibilities, we shall consider only the one taking care of maximal chains, which informally can be obtained as follows. One selects irreducible factors of P in such a way that their product, say

P_{k,1}(x, y) · P_{k,1}(q^{−1}x, p^{−1}y) ⋯ P_{k,1}(q^{−k+1}x, p^{−k+1}y),

forms a falling bibasic factorial [P_{k,1}]\underline{k} ; p,q of maximal length k. For the remaining irreducible factors of P this procedure is applied again in order to find all k-th falling factorial divisors [P_{k,1}]\underline{k} ; p,q, …, [P_{k,l}]\underline{k} ; p,q of that type. Then [P_{k}]\underline{k} ; p,q := [P_{k,1} ⋯ P_{k,l}]\underline{k} ; p,q forms the bibasic factorial factor of P of maximal length k. Iterating this procedure one gets a factorization of P in terms of “greatest” factorial factors.

Definition 3. We say that ⟨P_{1}, …, P_{k}⟩, P_{i} ∈ F [x, y], is a bibasic GFF-form of a bibasic monic polynomial P ∈ F [x, y], written as GFF_{p,q}(P) = ⟨P_{1}, …, P_{k}⟩, if the following conditions hold:

(GFF_{p,q} 1) P = [P_{1}]\underline{1} ; p,q ⋯ [P_{k}]\underline{k} ; p,q,

(GFF_{p,q} 2) each P_{i} is bibasic monic, and k > 0 implies P_{k} ≠ 1,

(GFF_{p,q} 3) for i ≤ j we have gcd_{p,q}([P_{i}]\underline{i} ; p,q, ϵP_{j}) = 1 = gcd_{p,q}([P_{i}]\underline{i} ; p,q, ϵ^{−j}P_{j}).

Note that GFF_{p,q}(1) = ⟨⟩. Condition (GFF_{p,q} 3) intuitively can be understood as prohibiting “overlaps” of bibasic factorials that violate length maximality. The following theorem states that, as in the q-hypergeometric case, the bibasic GFF-form is unique and thus provides a canonical form.

Theorem 1. If ⟨P_{1}, …, P_{k}⟩ and ⟨P′_{1}, …, P′_{l}⟩ are bibasic GFF-forms of a bibasic monic polynomial P ∈ F [x, y], then k = l and P_{i} = P′_{i} for all 1 ≤ i ≤ k.

Proof. The corresponding result for the ordinary hypergeometric case (p = q = 1) has been proved by Paule [8, Thm. 2.1]. The arguments used there extend immediately to the bibasic hypergeometric case proceeding by induction on d := deg_{x}(P) + deg_{y}(P).

From algorithmic point of view it is important to note that the bibasic GFF-form can be computed in an iterative manner essentially involving only gcd computations.

In q-hypergeometric summation, the normalized gcd of a polynomial P and its q-shift ϵP plays a fundamental role, as the gcd of P and its shift EP does in ordinary hypergeometric summation, where (EP)(x) = P (x + 1). The same is true for bibasic hypergeometric summation with respect to the bibasic shift operator ϵ. The mathematical and algorithmic essence lies in the following lemma.

Lemma 1 (Fundamental GFF_{p,q} Lemma). Let P ∈ F [x, y] be a bibasic monic polynomial with GFF_{p,q}(P) = ⟨P_{1}, …, P_{k}⟩. Then

gcd_{p,q}(P, ϵP) = [P_{1}]\underline{0} ; p,q ⋯ [P_{k}]\underline{k−1} ; p,q.

Proof. Due to the choice of the bibasic shift operator ϵ, the proof of the so-called Fundamental qGFF Lemma (cf. Paule and Riese [9, Lemma 1]) can be carried over to the bibasic hypergeometric case completely unchanged.

Thus, if GFF_{p,q}(P) = ⟨P_{1}, …, P_{k}⟩, then GFF_{p,q}(gcd_{p,q}(P, ϵP)) = ⟨P_{2}, …, P_{k}⟩. Consequently, dividing P with GFF_{p,q}(P) = ⟨P_{1}, …, P_{k}⟩ by ϵ^{−1} gcd_{p,q}(P, ϵP) or gcd_{p,q}(P, ϵP) results in separating the product of the first, respectively last, falling bibasic factorial entries, or in other words

\frac{P}{ϵ^{−1} gcd_{p,q}(P, ϵP)} = P_{1} · P_{2} ⋯ P_{k} and \frac{P}{gcd_{p,q}(P, ϵP)} = P_{1} · (ϵ^{−1}P_{2}) ⋯ (ϵ^{−k+1}P_{k}).

2.2 Bibasic Hypergeometric Telescoping

A sequence (f_{k})_{k∈ℤ} is said to be bibasic hypergeometric (see, e.g. Petkovšek, Wilf, and Zeilberger [11]) in p and q over F, if there exists a rational function ρ ∈ F (x, y) such that f_{k+1}/f_{k} = ρ(q^{k}, p^{k}) for all k where the quotient is well-defined.

Assume we are given a bibasic hypergeometric sequence (f_{k})_{k∈ℤ}. Then the problem of bibasic hypergeometric telescoping is to decide whether there exists a bibasic hypergeometric sequence (g_{k})_{k∈ℤ} such that

g_{k+1} − g_{k} = f_{k},  (1)

and if so, to determine (g_{k})_{k∈ℤ} with the motive that for a, b ∈ ℤ, a ≤ b,

∑_{k=a}^{b} f_{k} = g_{b+1} − g_{a},

which solves the indefinite summation problem.

For the rational function ρ, related to f_{k+1}/f_{k} as above, there exists a representation ρ(x, y) = z · x^{α} · y^{β} · A^{∗}(x, y)/B^{∗}(x, y) with bibasic monic A^{∗}, B^{∗} ∈ F [x, y], z ∈ F, and α, β ∈ ℤ, which we call a rational representation of the bibasic hypergeometric sequence (f_{k})_{k∈ℤ}. If additionally A^{∗} and B^{∗} are relatively prime, then ρ(x, y) is called the reduced rational representation of (f_{k})_{k∈ℤ}. For α ∈ ℤ, let α_{+} := max(α, 0) and α_{−} := max(−α, 0).

It will be shown below that bibasic hypergeometric telescoping can be decided constructively as follows.

Algorithm Telescope_{p,q}. Input: a bibasic hypergeometric sequence (f_{k})_{k∈ℤ} specified by its reduced rational representation ρ = z · x^{α} · y^{β} · A^{∗}/B^{∗}; Output: a bibasic hypergeometric solution (g_{k})_{k∈ℤ} of (1); in case such a solution does not exist, the algorithm stops.

(i) Compute the bibasic GP form of (f_{k})_{k∈ℤ}, i.e.,

(a) determine unique bibasic monic polynomials P^{∗}, Q^{∗}, R^{∗} ∈ F [x, y] such that

\frac{A^{∗}}{B^{∗}} = \frac{ϵP^{∗}}{P^{∗}} · \frac{Q^{∗}}{ϵR^{∗}},  (2)

where gcd_{p,q}(P^{∗}, Q^{∗}) = 1 = gcd_{p,q}(P^{∗}, R^{∗}) and gcd_{p,q}(Q^{∗}, ϵ^{j}R^{∗}) = 1 for all j ≥ 1, and

(b) let a_{x}, b_{x}, a_{y}, and b_{y} denote the coefficients of the lowest occurring powers of x and y in A^{∗}(x, 0), B^{∗}(x, 0), A^{∗}(0, y), and B^{∗}(0, y), respectively. Define

(γ, δ) := \begin{cases}(φ, ψ) & if α = 0 = β and q^{φ} · p^{μ} · b_{y}/a_{y} = z = p^{ψ} · q^{ν} · b_{x}/a_{x} for φ, ψ ∈ ℕ and μ, ν ∈ ℤ, \\ (φ, 0) & if α = 0 ≠ β and z = q^{φ} · p^{μ} · b_{y}/a_{y} for φ ∈ ℕ, μ ∈ ℤ, \\ (0, ψ) & if α ≠ 0 = β and z = p^{ψ} · q^{ν} · b_{x}/a_{x} for ψ ∈ ℕ, ν ∈ ℤ, \\ (0, 0) & otherwise, \end{cases}

and put

P := x^{γ} · y^{δ} · P^{∗},
Q := z · q^{−γ} · p^{−δ} · x^{α_{+}} · y^{β_{+}} · Q^{∗},  (3)
ϵR := x^{α_{−}} · y^{β_{−}} · ϵR^{∗},

with the motive that then

ρ = \frac{ϵP}{P} · \frac{Q}{ϵR}.

(ii) Try to solve the bibasic key equation

P = Q · ϵY − R · Y  (4)

for a polynomial Y ∈ F [x, y].

(iii) If such a polynomial solution Y exists, then

g_{k} = \frac{R(q^{k}, p^{k}) · Y (q^{k}, p^{k})}{P (q^{k}, p^{k})} · f_{k}  (5)

is a bibasic hypergeometric solution of (1), otherwise no bibasic hypergeometric solution (g_{k})_{k∈ℤ} exists.

The steps of Algorithm Telescope_{p,q} are derived as follows. First, assume that a bibasic hypergeometric solution (g_{k})_{k∈ℤ} with rational representation g_{k+1}/g_{k} = σ(q^{k}, p^{k}) of (1) exists. Then evidently we have

g_{k} = τ(q^{k}, p^{k}) · f_{k},  (6)

where τ(x, y) = 1/(σ(x, y) − 1) ∈ F (x, y).

By relation (6), equation (1) is equivalent to

z · x^{α_{+}} · y^{β_{+}} · A^{∗} · ϵτ − x^{α_{−}} · y^{β_{−}} · B^{∗} · τ = x^{α_{−}} · y^{β_{−}} · B^{∗},  (7)

where the reduced rational representation of (f_{k})_{k∈ℤ} is given by ρ = z · x^{α} · y^{β} · A^{∗}/B^{∗}.

Vice versa, any rational solution τ ∈ F (x, y) of (7) gives rise to a bibasic hypergeometric solution g_{k} := τ(q^{k}, p^{k})·f_{k} of (1). This means, bibasic hypergeometric telescoping is equivalent to finding a rational solution τ of (7).

Any τ ∈ F (x, y) can be represented as the quotient of relatively prime polynomials in the form τ = U/V where U, V ∈ F [x, y] with V = x^{φ} · y^{ψ} · V^{∗} the bibasic monic decomposition of V. In case such a solution τ of (7) exists, assume we know V or a multiple V ∈ F [x, y] of V. Then by clearing denominators in

z · x^{α_{+}} · y^{β_{+}} · A^{∗} · \frac{ϵU}{ϵV} − x^{α_{−}} · y^{β_{−}} · B^{∗} · \frac{U}{V} = x^{α_{−}} · y^{β_{−}} · B^{∗},

the problem reduces further to finding a polynomial solution U ∈ F [x, y] of the resulting difference equation with polynomial coefficients,

z · x^{α_{+}} · y^{β_{+}} · A^{∗} · V · ϵU − x^{α_{−}} · y^{β_{−}} · B^{∗} · (ϵV) · U = x^{α_{−}} · y^{β_{−}} · B^{∗} · V · ϵV.  (8)

Note that at least one polynomial solution, namely U = U·V/V, exists. Furthermore, equations of that type simplify by canceling gcd_{p,q}’s. For instance, in order to get more information about the denominator V, let V_{i} := ϵ^{i}V/ gcd_{p,q}(V, ϵV), i ∈ {0, 1}. Then (7) is equivalent to

z · x^{α_{+}} · y^{β_{+}} · A^{∗} · V_{0} · ϵU − x^{α_{−}} · y^{β_{−}} · B^{∗} · V_{1} · U = x^{α_{−}} · y^{β_{−}} · B^{∗} · V_{0} · V_{1} · gcd_{p,q}(V, ϵV). (9)

Now, if ⟨P_{1}, …, P_{m}⟩, m ∈ ℕ, is the bibasic GFF-form of V^{∗}, it follows from gcd_{p,q}(U, V) = 1 = gcd_{p,q}(V_{0}, V_{1}) and the Fundamental GFF_{p,q} Lemma that

V_{0} = (ϵ^{0}P_{1}) ⋯ (ϵ^{−m+1}P_{m}) | B^{∗} and V_{1} = q^{φ} · p^{ψ} · (ϵP_{1}) ⋯ (ϵP_{m}) | A^{∗}.

This observation gives rise to a simple and straightforward algorithm for computing a multiple V^{∗}:= [P_{1}]\underline{1} ; p,q ⋯ [P_{n}]\underline{n} ; p,q of V^{∗}. For instance, if P_{1} := gcd_{p,q}(ϵ^{−1}A^{∗}, B^{∗}) then obviously P_{1}|P_{1}. Actually, one can iteratively extract bibasic monic P_{i}-multiples P_{i} such that ϵP_{i}|A^{∗} and ϵ^{−i+1}P_{i}|B^{∗} by the following algorithm.

Algorithm V^{∗}MULT. Input: relatively prime and bibasic monic polynomials A^{∗}, B^{∗} ∈ F [x, y] that constitute the bibasic monic quotient of ρ = z · x^{α} · y^{β} · A^{∗}/B^{∗} ∈ F (x, y); Output: bibasic monic polynomials P_{1}, …, P_{n} such that V^{∗}:= [P_{1}]\underline{1} ; p,q ⋯ [P_{n}]\underline{n} ; p,q is a multiple of V^{∗}, the bibasic monic part of the denominator V = x^{φ} · y^{ψ} · V^{∗} of τ ∈ F (x, y).

(i) Compute n = min{j ∈ ℕ | gcd_{p,q}(ϵ^{−1}A^{∗}, ϵ^{k−1}B^{∗}) = 1 for all integers k > j}. (ii) Set A_{0} = A^{∗}, B_{0} = B^{∗}, and compute for i from 1 to n:

P_{i} = gcd_{p,q}(ϵ^{−1}A_{i−1}, ϵ^{i−1}B_{i−1}),
A_{i} = A_{i−1}/ϵP_{i},
B_{i} = B_{i−1}/ϵ^{−i+1}P_{i}.

A proof for the fact that the P_{i} are indeed multiples of the P_{i} has been worked out for the ordinary hypergeometric case by Paule [8, Lemma 5.1]. It can be carried over to the bibasic hypergeometric world almost word by word. Hence we leave the steps of the verification to the reader.

Note that in general step (i) of Algorithm V^{∗}MULT would be a rather time-consuming task involving resultant computations which could be solved by generalizing the univariate case (cf. Abramov, Paule, and Petkovšek [1]) in a straightforward way, for instance, as follows. Define R_{1}(v, w) := Res_{x}(A^{∗}(x, y), B^{∗}(vx, wy)) and R_{2}(v, w) := Res_{y}(A^{∗}(x, y), B^{∗}(vx, wy)), viewed as polynomials of v and w over F [y], respectively F [x]. Then n is the maximal positive integer such that R_{1}(q^{n}, p^{n}) · R_{2}(q^{n}, p^{n}) = 0 if such an integer exists, and n = 0 otherwise. However, in our implementation we make use of the fact that A^{∗} and B^{∗} already come in nicely factored form so that the computation of n boils down to a comparison of those factors.

Moreover, Algorithm V^{∗}MULT also delivers the constituents of the bibasic monic part of the GP representation (2) as stated in the following lemma.

Lemma 2. Let n, A_{n}, B_{n}, and the tuple ⟨P_{1}, …, P_{n}⟩ be computed as in Algorithm V^{∗}MULT. Then for P^{∗} = V^{∗}, Q^{∗} = A_{n}, and R^{∗} = ϵ^{−1}B_{n} we have

\frac{A^{∗}}{B^{∗}} = \frac{ϵP^{∗}}{P^{∗}} · \frac{Q^{∗}}{ϵR^{∗}},

where gcd_{p,q}(P^{∗}, Q^{∗}) = 1 = gcd_{p,q}(P^{∗}, R^{∗}) and gcd_{p,q}(Q^{∗}, ϵ^{j}R^{∗}) = 1 for all j ≥ 1.

For more details on GP representations in the q-hypergeometric case, see Abramov, Paule, and Petkovšek [1], or Paule and Strehl [10]. The results obtained there also apply in the bibasic hypergeometric case.

With the multiple V^{∗} of V^{∗} in hands, all what is left for solving (7), and thus the bibasic hypergeometric telescoping problem (1), is to determine appropriate multiplicities γ and δ such that

V = x^{γ} · y^{δ} · V^{∗} is a multiple of V = x^{φ} · y^{ψ} · V^{∗}.

For that we consider equation (9) again in the equivalent version

z · x^{α_{+}} · y^{β_{+}} · A^{∗} · V^{∗} · ϵU − x^{α_{−}} · y^{β_{−}} · B^{∗} · q^{φ} · p^{ψ} · (ϵV^{∗}) · U = x^{α_{−}} · y^{β_{−}} · B^{∗} · V^{∗} · ϵV, (10)

and distinguish the following cases corresponding to step (ib) of Algorithm Telescope_{p,q}.

(i) Assume that either α_{−} ≠ 0 or α_{+} ≠ 0. In the first case we have α_{+} = 0 and x^{α_{−}} | U, hence φ must be 0 because of gcd_{p,q}(U, V) = 1. This means, we can choose γ := 0. In the second case we have α_{−} = 0 and x^{min(α_{+},φ)} | U, because of ϵV = x^{φ} · y^{ψ} · q^{φ} · p^{ψ} · ϵV^{∗}. Again φ must be 0, and again we can choose γ := 0. Analogously, if β ≠ 0 we can choose δ := 0.

(ii) Assume that α = 0 and β ≠ 0, hence ψ = 0 by (i). For φ > 0, evaluating equation (10) at x = 0 results in

z · y^{β_{+}} · A^{∗}(0, y) · V^{∗}(0, y) · U(0, py) − y^{β_{−}} · B^{∗}(0, y) · q^{φ} · V^{∗}(0, py) · U(0, y) = 0. (11)

In order to evaluate (11) at y = 0, note that P ∈ F [x, y] being bibasic monic does not necessarily imply that P (0, y) ∈ F [y] is p-monic. To overcome this problem, let us consider the p-monic decompositions of U(0, y) and V^{∗}(0, y), say U(0, y) = u · y^{β_{u}} · Ū(y) and V^{∗}(0, y) = v · y^{β_{v}} · V̄(y), respectively. Now, dividing equation (11) by U(0, y) · V^{∗}(0, y) ≠ 0 leads to

z · y^{β_{+}} · A^{∗}(0, y) · p^{β_{u}} · \frac{Ū(py)}{Ū(y)} − y^{β_{−}} · B^{∗}(0, y) · q^{φ} · p^{β_{v}} · \frac{V̄(py)}{V̄(y)} = 0.  (12)

Additionally, let the p-monic decompositions of A^{∗}(0, y) and B^{∗}(0, y) be given by A^{∗}(0, y) = a_{y} · y^{β_{a}} · Ā(y) and B^{∗}(0, y) = b_{y} · y^{β_{b}} · B̄(y), respectively. Then the powers y^{β_{a}+β_{+}} and y^{β_{b}+β_{−}} must be equal, and after cancellation equation (12) at y = 0 turns into

z · a_{y} · p^{β_{u}} − b_{y} · q^{φ} · p^{β_{v}} = 0.

This means, we obtain as a condition for φ > 0 that z = q^{φ} · p^{μ} · b_{y}/a_{y} with μ ∈ ℤ. Hence, in this case we choose γ := φ, i.e., we set γ to this q-power if z has this particular form, and γ := 0 otherwise. Analogously, if α ≠ 0 and β = 0 we define δ := ψ > 0, if z = p^{ψ} · q^{ν} · b_{x}/a_{x} with ν ∈ ℤ, and δ := 0 otherwise.

(iii) Finally, for the case α = 0 = β similar reasoning as in case (ii) leads to the conditions

q^{φ} · p^{μ} · b_{y}/a_{y} = z = p^{ψ} · q^{ν} · b_{x}/a_{x},  (13)

for φ > 0 or ψ > 0, and μ, ν ∈ ℤ. Thus, if both conditions (13) are satisfied we choose γ := φ and δ := ψ, and otherwise γ = δ := 0.

The remaining steps of Algorithm Telescope_{p,q} now are explained as follows. Once again, employing the GP representation for the bibasic monic quotient of ρ,

\frac{A^{∗}}{B^{∗}} = \frac{ϵP^{∗}}{P^{∗}} · \frac{Q^{∗}}{ϵR^{∗}},

it is easily seen that equation (8) can be written as

z · q^{−γ} · p^{−δ} · x^{α_{+}} · y^{β_{+}} · \frac{Q^{∗}}{ϵR^{∗}} · ϵU − x^{α_{−}} · y^{β_{−}} · U = x^{γ+α_{−}} · y^{δ+β_{−}} · P^{∗}.  (14)

Because of relative primeness of certain polynomials, we observe that x^{α_{−}}| U, y^{β_{−}}| U, and ϵR^{∗} | ϵU. Hence by defining Y by the relation

U = x^{α_{−}} · y^{β_{−}} · q^{−α_{−}} · p^{−β_{−}} · R^{∗} · Y,

the task to solve equation (8) for U reduces to solve

z · q^{−γ} · p^{−δ} · x^{α_{+}} · y^{β_{+}} · Q^{∗} · ϵY − x^{α_{−}} · y^{β_{−}} · q^{−α_{−}} · p^{−β_{−}} · R^{∗} · Y = x^{γ} · y^{δ} · P^{∗}  (15)

for Y ∈ F [x, y]. By definition (3) of P, Q, and R, equation (15) immediately turns into the bibasic key equation (4),

Q · ϵY − R · Y = P.

Finally, from U/V = R · Y/P, again by definition (3), it follows directly that

g_{k} = \frac{R(q^{k}, p^{k}) · Y (q^{k}, p^{k})}{P (q^{k}, p^{k})} · f_{k}

as in (5) actually is a solution of the bibasic hypergeometric telescoping problem (1). This completes the proof of the correctness of Algorithm Telescope_{p,q}.

3 Degree Setting for Solving the Bibasic Key Equation

To solve the bibasic key equation

P = Q · ϵY − R · Y  (16)

we first have to determine degree bounds d_{1} and d_{2}, say, for the solution polynomial Y ∈ F [x, y] with respect to x and y, respectively, as shown in Theorem 2 below. Then we put

Y (x, y) := ∑_{i=0}^{d_{1}} ∑_{j=0}^{d_{2}} y_{i,j} · x^{i} · y^{j}

with undetermined y_{i,j} and solve (16) for the y_{i,j} by equating to zero all coefficients of x^{i}y^{j} in the equation

P − Q · ϵY + R · Y = 0,

which corresponds to solving a system of linear equations.

Theorem 2. Let l_{Q}^{x}(y), l_{Q}^{y}(x), l_{R}^{x}(y), and l_{R}^{y}(x) denote the leading coefficient polynomials of Q and R with respect to x and y, respectively. Let QR^{+} := Q + R and QR^{−} := Q − R. Then bounds for deg_{x}(Y) and deg_{y}(Y) are given by:

(A_{x}) If deg_{x}(QR^{+}) ≠ deg_{x}(QR^{−}), then

deg_{x}(Y) ≤ max{deg_{x}(P) − max{deg_{x}(QR^{+}), deg_{x}(QR^{−})}, 0}.

(A_{y}) If deg_{y}(QR^{+}) ≠ deg_{y}(QR^{−}), then

deg_{y}(Y) ≤ max{deg_{y}(P) − max{deg_{y}(QR^{+}), deg_{y}(QR^{−})}, 0}.

(B_{x}) If deg_{x}(QR^{+}) = deg_{x}(QR^{−}), then

(B1_{x}) if deg_{x}(Q) ≠ deg_{x}(R), then
deg_{x}(Y) = deg_{x}(P) − deg_{x}(QR^{+}),
(B2_{x}) if deg_{x}(Q) = deg_{x}(R), then

(B2a_{x}) if l_{R}^{x}(y)/l_{Q}^{x}(y) is of the form p^{μ} · q^{ν} · r(y) with μ, ν ∈ ℕ, and r(y) a rational function with r(0) = 1, then

deg_{x}(Y) ≤ max{deg_{x}(P) − deg_{x}(QR^{+}), ν},

(B2b_{x}) otherwise

deg_{x}(Y) = deg_{x}(P) − deg_{x}(QR^{+}).

(B_{y}) If deg_{y}(QR^{+}) = deg_{y}(QR^{−}), then

(B1_{y}) if deg_{y}(Q) ≠ deg_{y}(R), then
deg_{y}(Y) = deg_{y}(P) − deg_{y}(QR^{+}),

(B2_{y}) if deg_{y}(Q) = deg_{y}(R), then

(B2a_{y}) if l_{R}^{y}(x)/l_{Q}^{y}(x) is of the form p^{μ} · q^{ν} · r(x) with μ, ν ∈ ℕ, and r(x) a rational function with r(0) = 1, then

deg_{y}(Y) ≤ max{deg_{y}(P) − deg_{y}(QR^{+}), μ},

(B2b_{y}) otherwise

deg_{y}(Y) = deg_{y}(P) − deg_{y}(QR^{+}).

Proof. We rewrite the key equation to obtain

2 P = QR^{+} · (ϵY − Y) + QR^{−} · (ϵY + Y).  (17)

Cases (A_{x}) and (A_{y}) follow immediately. Note that it might happen that

deg_{x}(QR^{+}) > deg_{x}(P) and deg_{x}(QR^{−}) = deg_{x}(P),

and simultaneously

deg_{y}(QR^{+}) > deg_{y}(P) and deg_{y}(QR^{−}) = deg_{y}(P).

In this case, setting deg_{x}(Y) = deg_{y}(Y) = 0 could yield a solution, since ϵY − Y = 0 then.

For Case (B1_{x}) let a := deg_{x}(Q), c := deg_{x}(Y), and let l_{Y}^{x} (y) denote the leading coefficient polynomial of Y with respect to x. Assume that deg_{x}(Q) > deg_{x}(R). Then (17) gives

2 P (x, y) = (l_{Q}^{x}(y) x^{a} + …) · [(l_{Y}^{x} (py) q^{c} − l_{Y}^{x} (y)) x^{c} + …]
+ (l_{Q}^{x}(y) x^{a} + …) · [(l_{Y}^{x} (py) q^{c} + l_{Y}^{x} (y)) x^{c} + …]
= 2 l_{Q}^{x}(y) l_{Y}^{x} (py) q^{c} x^{a+c} + ….  (18)

Clearly, the coefficient of x^{a+c} in (18) will never vanish. Therefore we have

deg_{x}(Y) = deg_{x}(P) − deg_{x}(Q).

Including the case deg_{x}(Q) < deg_{x}(R), we obtain

deg_{x}(Y) = deg_{x}(P) − max{deg_{x}(Q), deg_{x}(R)} = deg_{x}(P) − deg_{x}(QR^{+}).

Analogous reasoning leads to Case (B1_{y}).

For Case (B2_{x}) we similarly observe that

2 P (x, y) = [(l_{Q}^{x}(y) + l_{R}^{x}(y)) x^{a} + …] · [(l_{Y}^{x} (py) q^{c} − l_{Y}^{x} (y)) x^{c} + …]
+ [(l_{Q}^{x}(y) − l_{R}^{x}(y)) x^{a} + …] · [(l_{Y}^{x} (py) q^{c} + l_{Y}^{x} (y)) x^{c} + …]
= 2 [l_{Q}^{x}(y) l_{Y}^{x} (py) q^{c} − l_{R}^{x}(y) l_{Y}^{x} (y)] x^{a+c} + ….  (19)

Now we no longer have the guarantee that the coefficient of x^{a+c} in (19) does not vanish, but it is easily seen that this happens only for

q^{c} = \frac{l_{R}^{x}(y)}{l_{Q}^{x}(y)} · \frac{l_{Y}^{x} (y)}{l_{Y}^{x} (py)}.  (20)

Note that l_{Y}^{x} (y) is actually not known. However, for any non-zero polynomial h(y) = h_{0} + h_{1} y + ⋯ + h_{d} y^{d}, the quotient h(y)/h(py) is of the form p^{−m} · s(y), where s(y) is a rational function with s(0) = 1 and m is the zero-root multiplicity of h(y). Hence, the rightmost fraction in (20) may eliminate only positive integer powers of p and a rational function of y but never introduce a power of q. This proves Case (B2a_{x}), and after interchanging x and p with y and q, respectively, also Case (B2a_{y}).

On the other hand, if the coefficient of x^{a+c} in (19) does not vanish, we obtain Case (B2b_{x}) and analogously Case (B2b_{y}).

4 Applications

In this section we shall illustrate the method of bibasic hypergeometric telescoping using the author’s Mathematica implementation qTelescope, which is a bibasic extension of a q-ana-logue of Gosper’s algorithm originally described in Paule and Riese [9].

Let the q-shifted factorial of a ∈ F be defined as usual (see, e.g. Gasper and Rahman [5]) by

(a; q)_{k} := \begin{cases}(1 − a)(1 − aq) ⋯ (1 − aq^{k−1}), & if k > 0, \\ 1, & if k = 0, \\ [(1 − aq^{−1})(1 − aq^{−2}) ⋯ (1 − aq^{k})]^{−1}, & if k < 0, \end{cases}

and

(a; q)_{∞} := ∏_{k=0}^{∞}(1 − aq^{k}),

where products of q-shifted factorials will be abbreviated by

(a_{1}, a_{2}, …, a_{n}; q)_{k} := (a_{1}; q)_{k} (a_{2}; q)_{k} ⋯ (a_{n}; q)_{k}.

In the present implementation we allow as summand any bibasic hypergeometric sequence (f_{k})_{k∈ℤ} of the form

f_{k} = \frac{∏_{r}(C_{r} q^{(c_{r}i_{r})k+d_{r}}; q^{i_{r}})_{a_{r}k+b_{r}}}{∏_{s}(D_{s} q^{(v_{s}j_{s})k+w_{s}}; q^{j_{s}})_{t_{s}k+u_{s}}} · \frac{∏_{r}(C′_{r} p^{(c′_{r}i′_{r})k+d′_{r}}; p^{i′_{r}})_{a′_{r}k+b′_{r}}}{∏_{s}(D′_{s} p^{(v′_{s}j′_{s})k+w′_{s}}; p^{j′_{s}})_{t′_{s}k+u′_{s}}}
× R(q^{k}, p^{k}) · q^{α\binom{k}{2}} · p^{β\binom{k}{2}} · z^{k},

with C_{r}, D_{s} power products in K(p), C′_{r}, D′_{s} power products in K(q), a_{r}, t_{s}, a′_{r}, t′_{s} specific integers (i.e., integers free of any parameters), b_{r}, u_{s}, b′_{r}, u′_{s} integer parameters free of k, or ±∞ if a_{r} (resp. t_{s}, a′_{r}, t′_{s}) = 0, c_{r}, v_{s}, c′_{r}, v′_{s} specific integers, d_{r}, w_{s}, d′_{r}, w′_{s} integer parameters free of k, i_{r}, j_{s}, i′_{r}, j′_{s} specific non-zero integers, R a rational function in F (q^{k}, p^{k}) such that the denominator factors completely

into a product of terms of the form (1 − D q^{vk+w}) and (1 − D′p^{v′k+w′}),

α, β specific integers, and z a rational function in F.

For the actual computation of the GP representation let ρ(x, y) denote the possibly nonreduced rational representation of the summand f_{k}. It is obvious from the input specification that ρ can always be converted into the form

ρ(x, y) = \frac{(ϵP̄)(x, y)}{P̄(x, y)} · \frac{∏_{i}(1 − Γ_{i} x^{γ_{i}})}{∏_{j}(1 − Δ_{j} x^{δ_{j}})} · \frac{∏_{i}(1 − Γ′_{i} y^{γ′_{i}})}{∏_{j}(1 − Δ′_{j} y^{δ′_{j}})} · xᾱ · yβ̄ · z̄

= \frac{(ϵP̄)(x, y)}{P̄(x, y)} · \frac{Ā(x, y)}{B̄(x, y)} · xᾱ · yβ̄ · z̄,

where P̄ ∈ F [x, y] is bibasic monic and satisfies gcd_{p,q}(P̄, Ā) = 1 = gcd_{p,q}(ϵP̄, B̄); the Γ_{i}, Δ_{j}, Γ′_{i}, Δ′_{j} are power products in F, the γ_{i}, δ_{j}, γ′_{i}, δ′_{j} are positive integers, ᾱ, β̄ ∈ ℤ, and z̄ ∈ F.

Concerning Algorithm V^{∗}MULT, it is clear from above that any P̄ ≠ 1 will actually contribute to [P_{1}]\underline{1} ; p,q and thus can be treated separately. Due to our input restrictions — this is the reason for admitting only power products instead of arbitrary rational functions — it is possible to find n in step (i) of Algorithm V^{∗}MULT simply by comparing all factors in Ā and B̄ as already mentioned.

Furthermore, since Ā and B̄ are both products of a q-monic and a p-monic polynomial, they will never contribute to b_{x}/a_{x} and b_{y}/a_{y}. Thus, b_{x}/a_{x} and b_{y}/a_{y} are in any case integer powers of q and p, respectively, coming from ϵP̄/P̄. Therefore, they do not take influence on the computation of γ and δ at all.

4.1 Bibasic Summation Formulas

In 1989, Gasper [3] derived the indefinite bibasic summation formula

∑_{k=0}^{n} f_{k} = ∑_{k=0}^{n} \frac{(1 − ap^{k}q^{k})(1 − bp^{k}q^{−k})}{(1 − a)(1 − b)} \frac{(a, b; p)_{k} (c, a/bc; q)_{k}}{(q, aq/b; q)_{k} (ap/c, bcp; p)_{k}} q^{k}

= \frac{(ap, bp; p)_{n} (cq, aq/bc; q)_{n}}{(q, aq/b; q)_{n} (ap/c, bcp; p)_{n}} = g_{n}  (21)

by showing that g_{k} is a bibasic hypergeometric solution of the equation f_{k} = g_{k} − g_{k−1}, however, without revealing how to come up with g_{k}. With our implementation the job of finding g_{k} is left to the computer.

In[1]:= <<qTelescope.m
Out[1]= Axel Riese's qTelescope implementation version 2.0 loaded
In[2]:= qTelescope[(1-a p^k q^k) (1-b p^k/q^k) qfac[a,p,k] qfac[b,p,k] qfac[c,q,k] *
                   qfac[a/b/c,q,k] q^k / ((1-a) (1-b) qfac[q,q,k] qfac[a q/b,q,k] *
                   qfac[a p/c,p,k] qfac[b c p,p,k]), {k, 0, n}]
                                             a q
        qfac[a p, p, n] qfac[b p, p, n] qfac[---, q, n] qfac[c q, q, n]
                                             b c
Out[2]= ---------------------------------------------------------------
             a p                                             a q
        qfac[---, p, n] qfac[b c p, p, n] qfac[q, q, n] qfac[---, q, n]
              c                                               b

Applying the same argumentation, Gasper and Rahman [4] generalized (21) to

∑_{k=−m}^{n}(1 − adp^{k}q^{k})(1 − bp^{k}/dq^{k}) \frac{(a, b; p)_{k} (c, ad^{2}/bc; q)_{k}}{(dq, adq/b; q)_{k} (adp/c, bcp/d; p)_{k}} q^{k}

= \frac{(1 − a)(1 − b)(1 − c)(1 − ad^{2}/bc)}{d(1 − c/d)(1 − ad/bc)}

× \{\frac{(ap, bp; p)_{n} (cq, ad^{2}q/bc; q)_{n}}{(dq, adq/b; q)_{n} (adp/c, bcp/d; p)_{n}} − \frac{(c/ad, d/bc; p)_{m+1} (1/d, b/ad; q)_{m+1}}{(1/c, bc/ad^{2}; q)_{m+1} (1/a, 1/b; p)_{m+1}} \}. (22)

Obviously, (21) is the case d = 1, m = 0 of (22). Since the output of qTelescope for identity (22) is quite lengthy, here we shall consider only the case m = −1 after dividing the summand by the constant fraction on the right hand side. Of course, the algorithm works for symbolic m as well.

In[3]:= qTelescope[(1-a d p^k q^k) (1-b/d p^k/q^k) qfac[a,p,k] qfac[b,p,k] *
                   qfac[c,q,k] qfac[a d^2/b/c,q,k] q^k d (1-c/d) (1-a d/b/c) /
                  (qfac[d q,q,k] qfac[a d q/b,q,k] qfac[a d p/c,p,k] *
                   qfac[b c p/d,p,k] (1-a) (1-b) (1-c) (1-a d^2/b/c)), {k, 1, n}]
                                                                      2
                                                                   a d  q
Out[3]= -1 + (qfac[a p, p, n] qfac[b p, p, n] qfac[c q, q, n] qfac[------, q, n]) /
                                                                    b c
         b c p             a d p                             a d q
   (qfac[-----, p, n] qfac[-----, p, n] qfac[d q, q, n] qfac[-----, q, n])
           d                 c                                 b

4.2 Bibasic Matrix Inverses

Al-Salam and Verma [2] showed that the triangular matrices H = (h_{n,k}) and G = (g_{k,n}), where

h_{n,k} = \frac{(−1)^{n+k} (hqp^{n}; q)_{n−1} (1 − hq^{k}p^{k})}{(p; p)_{n−k} (hqp^{n}; q)_{k}}

and

g_{k,n} = \frac{(hp^{n}q^{n}; q)_{k−n}}{(p; p)_{k−n}} p^{\binom{k−n}{2}}

are inverse to each other. This result is equivalent to the fact that

∑_{k=m}^{n} h_{n,k} · g_{k,m} = δ_{n,m},  (23)

where δ_{n,m} denotes the Kronecker symbol. Running the algorithm we obtain:

In[4]:= qTelescope[(-1)^(n+k) qfac[h q p^n,q,n-1] (1-h q^k p^k) *
                   qfac[h p^m q^m,q,k-m] p^Binomial[k-m,2] /
                  (qfac[p,p,n-k] qfac[h q p^n,q,k] qfac[p,p,k-m]), {k, m, n}]
Out[4]= {0, {-m + n != 0}}

This means, we algorithmically proved identity (23) for m ≠ n, but evaluation failed for m = n. However, it is easily seen that h_{n,n} · g_{n,n} = 1, which completes the proof.

These matrices were used in a slightly modified form also by Gessel and Stanton [6] in the derivation of a family of q-Lagrange inversion formulas.

Al-Salam and Verma [2] employed the fact that the n-th q-difference of a polynomial of degree less than n is equal to zero, to show that

(1 − \frac{a}{q}) ∑_{k=0}^{n} \frac{(−1)^{k} (ap^{k}; q)_{n−1}}{(p; p)_{k} (p; p)_{n−k}} p^{\binom{k}{2}} = δ_{n,0}.  (24)

Unfortunately, for d_{k} := (ap^{k}; q)_{n−1}, we find that

\frac{d_{k+1}}{d_{k}} = \frac{(1 − ap^{k+1})(1 − ap^{k+1}q) ⋯ (1 − ap^{k+1}q^{n−2})}{(1 − ap^{k})(1 − ap^{k}q) ⋯ (1 − ap^{k}q^{n−2})}

is a rational function of q^{k} and p^{k} only for fixed n. Therefore d_{k} is not a valid input for the algorithm. To overcome the problem, we replace k, n, and a in (24) by k − m, n − m, and a^{−1}p^{m}q^{1−n}, respectively, such that (24) turns into the orthogonality relation

c_{n,m} ∑_{k=m}^{n} a_{n,k} · b_{k,m} = δ_{n,m}  (25)

with

c_{n,m} = (1 − a^{−1}p^{m}q^{−n}) a^{1+m−n} q^{\binom{m+1}{2}−\binom{n}{2}},
a_{n,k} = \frac{(ap^{−k}; q)_{n}}{(p; p)_{n−k}} (−1)^{1+k+n} p^{\binom{n−k}{2}},
b_{k,m} = \frac{p^{−k(m+1)}}{(p; p)_{k−m} (ap^{−k}; q)_{m+1}} .

Note that a_{n,k} and b_{k,m} still do not fit into the input specification of the algorithm. For A = (a_{n,k}), B = (b_{k,m}), and C = (c_{n,m}), relation (25) could be rewritten as A·B = diag(C)^{−1}, showing that the matrix diag(C) · A = (c_{n,n} · a_{n,k}) is inverse to the matrix B. Since inverse matrices commute, (25) is equivalent to

∑_{k=m}^{n} c_{k,k} · b_{n,k} · a_{k,m} = δ_{n,m},

or, in other words

∑_{k=m}^{n} \frac{(−1)^{k+m} (1 − ap^{−k}q^{k}) (ap^{−m}; q)_{k}}{(p; p)_{n−k} (p; p)_{k−m} (ap^{−n}; q)_{k+1}} p^{\binom{k−m}{2}−n(k+1)+k(m+1)} = δ_{n,m}.  (26)

In[5]:= qTelescope[(-1)^(k+m) (1-a q^k/p^k) qfac[a/p^m,q,k] *
                   p^(Binomial[k-m,2]-n(k+1)+k(m+1)) /
                   (qfac[p,p,n-k] qfac[p,p,k-m] qfac[a/p^n,q,k+1]), {k, m, n}]
Out[5]= {0, {-m + n != 0}}

For m = 0, (26) reduces to

∑_{k=0}^{n} \frac{(−1)^{k} (1 − ap^{−k}q^{k}) (a; q)_{k}}{(p; p)_{n−k} (p; p)_{k} (ap^{−n}; q)_{k+1}} p^{\binom{n−k}{2}} = δ_{n,0}.

Out[6]= {0, {n != 0}}

Proceeding in the same way, we can prove the bibasic identity (cf. Gasper [3])

(1 − \frac{a}{q})(1 − \frac{b}{q} ) ∑_{k=0}^{n} \frac{(−1)^{k} (ap^{k}, bp^{−k}; q)_{n−1} (1 − ap^{2k}/b)}{(p; p)_{k} (p; p)_{n−k} (ap^{k}/b; p)_{n+1}} p^{k(n−1)+\binom{n−k}{2}} = δ_{n,0},

In[6]:= qTelescope[(-1)^k (1-a q^k/p^k) qfac[a,q,k] p^Binomial[n-k,2] /

(qfac[p,p,n-k] qfac[p,p,k] qfac[a/p^n,q,k+1]), {k, 0, n}]

by transforming it into the equivalent version

(1 − \frac{b}{a}) ∑_{k=0}^{n}(1 − ap^{−k}q^{k})(1 − bp^{k}q^{k}) \frac{(−1)^{k} (a, b; q)_{k} (bp^{k+1}/a; p)_{n−1}}{(p; p)_{k} (p; p)_{n−k} (ap^{−n}, bp^{n}; q)_{k+1}} p^{\binom{n−k}{2}} = δ_{n,0}.

In[7]:= qTelescope[(1-b/a) (1-a q^k/p^k) (1-b p^k q^k) (-1)^k qfac[a,q,k] *
                   qfac[b,q,k] qfac[b/a p^(k+1),p,n-1] p^Binomial[n-k,2] /
                  (qfac[p,p,k] qfac[p,p,n-k] qfac[a/p^n,q,k+1] qfac[b p^n,q,k+1]),
                   {k, 0, n}]
Out[7]= {0, {n != 0}}

4.3 Open Problems

With the input specification described above we actually have not taken into account that a bibasic hypergeometric summand f_{k} could involve q-shifted factorials with mixed bases such as (a; p^{i}q^{j})_{k} for i, j ∈ ℤ as well. However, since to our knowledge applications of this type have not arisen in practice up to now, this feature has not been implemented yet.

For the sake of simplicity we restricted ourselves to discuss in detail the bibasic case. Nevertheless, the presented approach should easily extend to the multi-basic case, i.e., to sequences being hypergeometric in independent bases q_{1}, …, q_{m}.

So far we found only one single bibasic example in the literature which we could not handle with our machinery, namely Gasper’s [3] transformation formulas

∑_{k=0}^{∞} \frac{1 − ap^{k}q^{k}}{1 − a} \frac{(a; p)_{k} (c/b; q)_{k}}{(q; q)_{k} (abp; p)_{k}} b^{k}

= \frac{1 − c}{1 − b} ∑_{k=0}^{∞} \frac{(ap; p)_{k} (c/b; q)_{k}}{(q; q)_{k} (abp; p)_{k}} (bq)^{k}

= \frac{1 − c}{1 − abp} ∑_{k=0}^{∞} \frac{(ap; p)_{k} (cq/b; q)_{k}}{(q; q)_{k} (abp^{2}; p)_{k}} b^{k}

= \frac{(1 − c) (ap; p)_{∞}}{(1 − b) (abp; p)_{∞}} ∑_{k=0}^{∞} \frac{(b; p)_{k} (cqp^{k}; q)_{∞}}{(p; p)_{k} (bqp^{k}; p)_{∞}} (ap)^{k},

when max(|p|, |q|, |ap|, |b|) < 1.

Acknowledgment. I wish to thank Peter Paule for his cooperation and comments.

References

[1] S.A. Abramov, P. Paule, and M. Petkovšek, q-Hypergeometric solutions of q-difference equations, preprint, 1995.

[2] W. Al-Salam and A. Verma, On quadratic transformations of basic series, SIAM J. Math. Anal., 15 (1984), 414–421.

[3] G. Gasper, Summation, transformation, and expansion formulas for bibasic series, Trans. Amer. Math. Soc., 312 (1989), 257–277.

[4] G. Gasper and M. Rahman, An indefinite bibasic summation formula and some quadratic, cubic, and quartic summation and transformation formulas, Canad. J. Math., 42 (1990), 1–27.

[5] G. Gasper and M. Rahman, Basic Hypergeometric Series, Encyclopedia of Mathematics and its Applications, 35 (G.-C. Rota, ed.), Cambridge University Press, London and New York, 1990.

[6] I.M. Gessel and D. Stanton, Applications of q-Lagrange inversion to basic hypergeometric series, Trans. Amer. Math. Soc., 277 (1983), 173–201.

[7] R.W. Gosper, Decision procedures for indefinite hypergeometric summation, Proc. Natl. Acad. Sci. U.S.A., 75 (1978), 40–42.

[8] P. Paule, Greatest factorial factorization and symbolic summation, J. Symbolic Computation, 20 (1995), 235–268.

[9] P. Paule and A. Riese, A Mathematica q-analogue of Zeilberger’s algorithm based on an algebraically motivated approach to q-hypergeometric telescoping, preprint, to appear in Fields Proceedings of the Workshop on “Special Functions, q-Series and Related Topics”, organized by the Fields Institute for Research in Mathematical Sciences at University College, 12–23 June 1995, Toronto, Ontario.

[10] P. Paule and V. Strehl, Symbolic summation — some recent developments, Computeralgebra in Science and Engineering — Algorithms, Systems, Applications (J. Fleischer, J. Grabmeier, F. Hehl, and W. Küchlin, eds.), pp. 138–162, World Scientific, Singapore, 1995.

[11] M. Petkovšek, H.S. Wilf, and D. Zeilberger, A = B, A.K. Peters, 1996.

[12] A. Riese, A Mathematica q-analogue of Zeilberger’s algorithm for proving q-hypergeometric identities, Diploma thesis, J. Kepler University, Linz, 1995.

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

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