Multibasic and Mixed Hypergeometric Gosper-Type Algorithms

Andrej Bauer

Department of Computer Science Carnegie Mellon University Pittsburgh PA 15213, U.S.A.

         Andrej.Bauer@cs.cmu.edu

Marko Petkovšek^{∗}

Department of Mathematics University of Ljubljana Jadranska 19, 1111 Ljubljana, Slovenia

                                           Marko.Petkovsek@fmf.uni-lj.si

Abstract

Gosper’s summation algorithm finds a hypergeometric closed form of an indefinite sum of hypergeometric terms, if such a closed form exists. We extend his algorithm to the case when the terms are simultaneously hypergeometric and multibasic hypergeometric. We also provide algorithms for finding polynomial as well as hypergeometric solutions of recurrences in the mixed case. We do not require the bases to be transcendental, but only that q_{1}^{k_{1}} ⋯ q_{m}^{k_{m}} ≠ 1 unless k_{1} = ⋯ = k_{m} = 0. Finally, we generalize the concept of greatest factorial factorization to the mixed hypergeometric case.

1 Introduction and notation

Let 𝔽 be a field of characteristic zero and ⟨t_{n}⟩_{n=0}^{∞} a sequence of elements from 𝔽 which is eventually non-zero. Call t_{n}:

• hypergeometric, if there are polynomials p_{1}, p_{2} ∈ 𝔽[x] such that p_{1}(n)t_{n+1} = p_{2}(n)t_{n} for all n;

• q-hypergeometric or basic hypergeometric, if there are polynomials p_{1}, p_{2} ∈ 𝔽[x] such that p_{1}(q^{n})t_{n+1} = p_{2}(q^{n})t_{n} for all n, where q ∈ 𝔽 \ {0} is the base;

• multibasic hypergeometric, if there are polynomials p_{1}, p_{2} ∈ 𝔽[y_{1}, …, y_{m}] such that p_{1}(q_{1}^{n}, …, q_{m}^{n})t_{n+1} = p_{2}(q_{1}^{n}, …, q_{m}^{n})t_{n} for all n, where q_{1}, …, q_{m} ∈ 𝔽 \ {0} are the bases;

• mixed hypergeometric, if there are polynomials p_{1}, p_{2} ∈ 𝔽[x, y_{1}, …, y_{m}] such that p_{1}(n, q_{1}^{n}, …, q_{m}^{n})t_{n+1} = p_{2}(n, q_{1}^{n}, …, q_{m}^{n})t_{n} for all n.

The well-known Gosper’s algorithm [8, 9] finds hypergeometric solutions f_{n} of the nonhomogeneous first-order recurrence

f_{n+1} − f_{n} = t_{n}

where t_{n} is a given hypergeometric sequence. Besides its obvious use for indefinite hypergeometric summation, it also plays a crucial role in the algorithms for definite hypergeometric summation, construction of annihilating recurrences, and automated verification of identities [25, 26, 23]. Therefore it is not surprising that analogous algorithms have been designed for many other settings, e.g., integration of hyperexponential functions [4], basic [24, 13, 17] and bibasic [20] hypergeometric summation. We generalize Gosper’s algorithm, as well as some related ones, to the mixed hypergeometric case.

The algebraic setting of the paper (with the exception of Section 8) is the rational-function field 𝔽(x, y) where 𝔽 is an arbitrary field of characteristic zero, together with an 𝔽-automorphism E which acts by Ex = x + 1 and Ey_{i} = q_{i}y_{i}. This is discussed in detail in Section 2. Some auxiliary algorithms used later as subroutines are sketched in Section 3, while in Sections 4 and 5 the necessary ingredients for Gosper’s algorithm are developed. Although there only first-order recurrences are checked for polynomial solutions, we provide in Section 4 algorithm MixedPoly^{1} which finds all polynomial solutions of a parametric nonhomogeneous polynomial-coefficient recurrence of any order. A mixed hypergeometric canonical form of rational functions is described in Section 5. After these preparations, we present in Section 6 an analogue of Gosper’s algorithm for the mixed hypergeometric case. Our algorithm MixedGosper^{2} is a common generalization of the algorithms presented in [9, 24, 20]. When specialized to the bibasic case, it essentially agrees with the algorithm given in [20]. However, looking at the case analysis in the computation of the multiplicities γ and δ [20, pp. 7–8], it is not immediately clear how to extend that to the multibasic case. In Section 7 we provide algorithm MixedHyper which finds all mixed hypergeometric solutions of a homogeneous polynomial-coefficient recurrence of any order. This is a common generalization of the algorithms presented in [19] and [3]. In Section 8 we extend the concept of greatest factorial factorization [16] to an arbitrary automorphism σ of the multivariate polynomial ring.

∗corresponding author 1available at http://www.cis.upenn.edu/˜ wilf/AeqB.html in the Mathematica package gosper.m as MixedPoly

Notation. The set of integers is denoted by ℤ, the set of nonnegative integers by ℕ_{0}, and the field of rational numbers by ℚ.

If n, m ∈ ℕ_{0} and a = (a_{1}, a_{2}, …, a_{m}), b = (b_{1}, b_{2}, …, b_{m}) are m-tuples of elements of a ring, we write ab for the componentwise product (a_{1}b_{1}, a_{2}b_{2}, …, a_{m}b_{m}), and a^{n} for the componentwise power (a_{1}^{n}, a_{2}^{n}, …, a_{m}^{n}). If α = (α_{1}, α_{2}, …, α_{m}) ∈ ℕ_{0}^{m} then we write a^{α} for the power product a_{1}^{α_{1}}a_{2}^{α_{2}} ⋯ a_{m}^{α_{m}}.

We say that two multivariate polynomials over a field are coprime if they do not have a non-constant common factor. When a and b are coprime, we write a ⊥ b. When S is a set of polynomials and a ⊥ b for all b ∈ S, we write a ⊥ S.

2 Algebraic preliminaries

Let 𝔽 be a field of characteristic zero. Let q_{1}, …, q_{m} ∈ 𝔽 \ {0}, and suppose that for any integers k_{1}, …, k_{m} ∈ ℤ,

q_{1}^{k_{1}}q_{2}^{k_{2}} ⋯ q_{m}^{k_{m}} = 1 ⟹ k_{1} = k_{2} = ⋯ = k_{m} = 0.  (1)

This generalizes the condition that q is not a root of unity in the q-hypergeometric case [3]. For example, if 𝔽 = ℝ, q_{1} = \sqrt[3]{2} and q_{2} = \sqrt[5]{2}, then q_{1}^{3}q_{2}^{−5} = 1 and we should have chosen q = \sqrt[15]{2} in the first place. On the other hand, q_{1} = 2 and q_{2} = 3 would be a legitimate choice in this case. We call q_{i}’s the bases, and write q = (q_{1}, q_{2}, …, q_{m}).

Let y = (y_{1}, y_{2}, …, y_{m}) be an m-tuple of variables, 𝔽[x, y] the ring of polynomials over 𝔽 in x and y, and 𝔽(x, y) the corresponding rational function field. We define an 𝔽-automorphism E of 𝔽(x, y) (i.e., E is a field automorphism of 𝔽(x, y) which fixes each element of 𝔽 ⊆ 𝔽(x, y)) by stipulating further that Ex = x + 1 and Ey_{k} = q_{k}y_{k} for k = 1, …, m. Then 𝔽(x, y) together with E is a difference field and 𝔽[x, y] is a difference subring of 𝔽(x, y) (see [6] for the relevant definitions).

Let M be the set of power products in y_{1}, y_{2}, …, y_{m}:

M = {y_{1}^{k_{1}}y_{2}^{k_{2}} ⋯ y_{m}^{k_{m}} | k_{i} ∈ ℕ_{0} for i = 1, …, m}.

If u = y_{1}^{k_{1}}y_{2}^{k_{2}} ⋯ y_{m}^{k_{m}} ∈ M, we write u(q) for the corresponding power product of the bases q_{1}^{k_{1}}q_{2}^{k_{2}} ⋯ q_{m}^{k_{m}}. Note that Eu = u(q)u for all u ∈ M.

As a multiplicative monoid, M is obviously isomorphic to ℕ_{0}^{m}, the direct product of m copies of the additive monoid ℕ_{0}. We denote by ⪯ an admissible term order in ℕ_{0}^{m}, which is a total order satisfying

1. 0 ⪯ α,

2. α ⪯ β ⇒ α + γ ⪯ β + γ

for all α, β, γ ∈ ℕ_{0}^{m}. An example of an admissible term order is the lexicographic order ⪯_{lex}, with α ≺_{lex} β when α ≠ β and α_{k} < β_{k} where k = min{i; α_{i} ≠ β_{i}}.

Definition 2.1 Let α, β ∈ ℕ_{0}^{m}. Then we write

α ⊆ β

whenever α_{i} ≤ β_{i} for all i between 1 and m.

Clearly, (ℕ_{0}^{m}, ⊆) is a partial order isomorphic to (M, |) where | denotes divisibility of power products, and is contained in any admissible term order:

α ⊆ β ⇒ α ⪯ β, for all α, β ∈ ℕ_{0}^{m}.

2available at http://www.cis.upenn.edu/˜ wilf/AeqB.html in the Mathematica package gosper.m as GosperSum

We adjoin to ℕ_{0}^{m} an absorbing bottom element^{3}, ⊥, such that for all α ∈ ℕ_{0}^{m}

⊥ ≺ α,
⊥ + α = α + ⊥ = ⊥.

Definition 2.2 Let p ∈ 𝔽[x, y]. Write

p(x, y) = ∑_{α∈ℕ_{0}^{m}} p_{α}(x)y^{α} = ∑_{i∈ℕ_{0}} c_{i}(y)x^{i}  (2)

where only finitely many p_{α} ∈ 𝔽[x] and c_{i} ∈ 𝔽[y] are non-zero.

1. We define the multidegree of p in y as

deg_{y} p(x, y) = \{ \begin{matrix}max_{⪯}\{α ∈ ℕ_{0}^{m}; p_{α} ≠ 0\}, & p ≠ 0, \\ ⊥, & p = 0. \end{matrix}

2. Similarly,

mindeg_{y}p(x, y) = \{ \begin{matrix}min_{⪯}\{α ∈ ℕ_{0}^{m}; p_{α} ≠ 0\}, & p ≠ 0, \\ ⊥, & p = 0. \end{matrix}

3. We write [y^{α}] p(x, y) for p_{α}(x) and [x^{i}] p(x, y) for c_{i}(y) in (2).

4. When deg_{y} p ≺ α we write p = o(y^{α}).

5. Let δ = deg_{y}p. We call p mixed monic when [y^{δ}] p(x, y) is monic as a univariate polynomial in x.

Note that the concepts of multidegree and mixed monicity are relative to the chosen term order ⪯. By convention, gcd(a, b) always denotes a mixed monic greatest common divisor of a, b ∈ 𝔽[x, y].

We need the following well-known result from the theory of linear recurrent sequences.

Lemma 2.1 Let 𝔽 be a field of characteristic zero and r_{1}, …, r_{k} ∈ 𝔽 \ {0}, with r_{i} ≠ r_{j} for i ≠ j. Let d_{1}, d_{2}, …, d_{k} ∈ ℕ_{0} and d = d_{1} + d_{2} + ⋯ + d_{k}. Then the d functions g_{ij} : ℕ_{0} → 𝔽, defined by g_{ij}(n) = n^{j} r_{i}^{n}, for i = 1, …, k, j = 0, …, d_{i} − 1, are linearly independent in the vector space ℕ_{0} → 𝔽 over 𝔽.

For a proof, see, e.g., [22, Thm. 4.1.1].

The main object of our interest is the ring of sequences ℕ_{0} → 𝔽. To simplify notation, we denote the sequence ⟨0, 1, 2, …⟩ by n, and ⟨1, q_{i}, q_{i}^{2}, …⟩ by q_{i}^{n}, for i = 1, 2, …, m. We write 𝔽[n, q^{n}] for the subring of ℕ_{0} → 𝔽 generated by n, q_{1}^{n}, …, q_{m}^{n} and the constant sequences, and call its elements the polynomial sequences. This is justified by the following theorem.

Theorem 2.2 Let Φ : 𝔽[x, y] → 𝔽[n, q^{n}] be the ring homomorphism mapping x ↦ n and y_{i} ↦ q_{i}^{n}. Then Φ is an isomorphism between the ring of polynomials 𝔽[x, y] and the ring of polynomial sequences 𝔽[n, q^{n}].

Proof: It is obvious that Φ is an epimorphism. We show that it is a monomorphism. Let f ∈ 𝔽[x, y]. Write f as

f = ∑_{i=1}^{k} p_{i}u_{i},

where p_{1}, …, p_{k} ∈ 𝔽[x], u_{1}, …, u_{k} ∈ M, and u_{i} ≠ u_{j} for i ≠ j. Suppose Φf = 0:

0 = Φf = ∑_{i=1}^{k} p_{i}(n)u_{i}(q)^{n}.

Because q_{1}, …, q_{m} satisfy condition (1), u_{i}(q) ≠ u_{j}(q) for i ≠ j. The result now follows from Lemma 2.1.

3not to be confused with the notation for coprime polynomials introduced at the end of Section 1

As a consequence, 𝔽[n, q^{n}] is an integral domain, and its field of fractions 𝔽(n, q^{n}) whose elements we call the rational sequences is isomorphic to the rational function field 𝔽(x, y). The map Φ : 𝔽[x, y] → 𝔽[n, q^{n}] defined in Theorem 2.2 can be naturally extended to a map from 𝔽(x, y) to 𝔽(n, q^{n}).

We define a homomorphism S on ℕ_{0} → 𝔽 by setting (Sa)(n) = a(n + 1) for all a: ℕ_{0} → 𝔽. This makes 𝔽(n, q^{n}) into a difference field and 𝔽[n, q^{n}] into a difference subring of 𝔽(n, q^{n}). As Φ ∘ E = S ∘ Φ, we see that Φ extended to 𝔽(x, y) is a difference isomorphism of the two fields 𝔽(x, y) and 𝔽(n, q^{n}), as well as of the two rings 𝔽[x, y] and 𝔽[n, q^{n}]. This allows us to work in 𝔽(x, y) resp. 𝔽[x, y] instead of in 𝔽(n, q^{n}) resp. 𝔽[n, q^{n}] whenever suitable.

We conclude this section by two simple technical lemmas.

Lemma 2.3 Let p ∈ 𝔽[x, y] \ {0}, k ∈ ℤ \ {0}, and a ∈ 𝔽. Then E^{k}p = ap if and only if p = ru for some r ∈ 𝔽, u ∈ M, and a = u(q)^{k}.

Proof: Sufficiency is obvious. Suppose E^{k}p = ap. Write p as

p(x, y) = ∑_{i=1}^{n} p_{i}(x) u_{i}(y),

where u_{1}, …, u_{n} ∈ M are pairwise different and p_{1}, …, p_{n} ∈ 𝔽[x] \ {0}. It follows that

∑_{i=1}^{n} a p_{i}(x) u_{i}(y) = ap = E^{k}p = ∑_{i=1}^{n} p_{i}(x + k) u_{i}(q)^{k} u_{i}(y).

Hence, for i = 1, …, n

a p_{i}(x) = u_{i}(q)^{k} p_{i}(x + k).

By comparing the leading coefficients in the above equation, we conclude that a = u_{i}(q)^{k} for i = 1, …, n. However, if it were the case that u_{i}(q)^{k} = u_{j}(q)^{k} for some i ≠ j, condition (1) would be violated. It follows that n = 1, and p(x, y) = r(x) u(y) for some r ∈ 𝔽[x] \ {0} and u ∈ M. From E^{k}p = ap we get r(x + k) = r(x), which is only possible if r is a constant.

Definition 2.3 For 1 ≤ i ≤ m, we denote by π_{i} the endomorphism of 𝔽[x, y] which substitutes 0 for y_{i}.

Lemma 2.4 The endomorphisms π_{i}, 1 ≤ i ≤ m, commute with E and E^{−1}.

Proof: Let Y_{i} = {y_{1}, y_{2}, …, y_{m}} \ {y_{i}} and p ∈ 𝔽[x, y]. Consider p to be a polynomial in 𝔽[x, Y_{i}][y_{i}]. It is easy to check that Eπ_{i}p = π_{i}Ep and E^{−1}π_{i}p = π_{i}E^{−1}p.

3 Algorithmic preliminaries

The main algorithmic subproblems that we encounter are the following:

1. (disp) For polynomials a, b ∈ 𝔽[x, y] \ {0} such that a, b ⊥ M, compute the dispersion set

D(a, b) = {n ∈ ℕ_{0}; a ⊥̸ E^{n}b}

containing all nonnegative integers n such that a and E^{n}b have a non-constant common divisor.

2. (introot) Find the set of all nonnegative integer roots n of P (n, q^{n}) = 0 where P ∈ 𝔽[x, y] \ {0}.

3. (qmon) Given a ∈ 𝔽 \ {0}, find integers k_{1}, …, k_{m} (if any) such that a = q_{1}^{k_{1}} ⋯ q_{m}^{k_{m}}. Note that by (1) such integers are unique.

In Section 3.1 we reduce disp to introot. In Sections 3.2 and 3.3 we show how to solve introot in two important special cases when 𝔽 = ℚ(q_{1}, …, q_{m}) is a purely transcendental extension of ℚ, and when 𝔽 = ℚ (and hence q_{1}, …, q_{m} ∈ ℚ), respectively. We do not elaborate on qmon, because in the two special cases of transcendental resp. rational bases it is rather obvious how to solve it.

3.1 Computing the dispersion set

Define polynomials R_{1}, R_{2}, …, R_{m}, R ∈ 𝔽[x, y][ξ, η] as polynomial resultants

R_{i}(ξ, η) = Res_{y_{i}}(a(x, y), b(x + ξ, ηy)) (1 ≤ i ≤ m),
R(ξ, η) = Res_{x}(a(x, y), b(x + ξ, ηy)).

Here ξ is a variable and η = (η_{1}, η_{2}, …, η_{m}) is an m-tuple of variables. Let

P (ξ, η) = R(ξ, η) ∏_{i=1}^{m} R_{i}(ξ, η).  (3)

The following lemma leads to an algorithm for computing D(a, b): Lemma 3.1 D(a, b) = {n ∈ ℕ_{0}; P (n, q^{n}) = 0}. Proof: For n ∈ ℕ_{0}, let ϕ_{n} : 𝔽[x, y, ξ, η] → 𝔽[x, y] be the evaluation homomorphism which substitutes n for ξ and q^{n} for η. It is easy to see that for any non-zero polynomial p ∈ 𝔽[x, y], the homomorphic image ϕ_{n}(p(x + ξ, ηy)) = p(x + n, q^{n}y) = E^{n}p(x, y) is non-zero. Therefore, by the Homomorphism Lemma for resultants (see, e.g., [14, Lemma 7.3.1]),

R_{i}(n, q^{n}) = ϕ_{n}(R_{i}(ξ, η)) = Res_{y_{i}}(ϕ_{n}(a(x, y)), ϕ_{n}(b(x + ξ, ηy))) = Res_{y_{i}}(a, E^{n}b)
(1 ≤ i ≤ m),
R(n, q^{n}) = ϕ_{n}(R(ξ, η)) = Res_{x}(ϕ_{n}(a(x, y)), ϕ_{n}(b(x + ξ, ηy))) = Res_{x}(a, E^{n}b).

Thus we have the following chain of equivalences:

n ∈ D(a, b) ⟺ one of deg_{x} gcd(a, E^{n}b), deg_{y_{i}} gcd(a, E^{n}b) is positive
⟺ one of Res_{x}(a, E^{n}b), Res_{y_{i}}(a, E^{n}b) vanishes
⟺ one of R(n, q^{n}), R_{i}(n, q^{n}) vanishes
⟺ R(n, q^{n}) ∏_{i=1}^{m} R_{i}(n, q^{n}) = 0
⟺ P (n, q^{n}) = 0.  (4)

The second equivalence above follows from the well-known properties of polynomial resultants.

Next we show how to find integral solutions n of equation (4) in two special cases.

3.2 Transcendental bases

Let 𝔽 = ℚ(q_{1}, …, q_{m}) where q_{1}, …, q_{m} are algebraically independent over ℚ. Let p ∈ 𝔽[x, y] \ {0}. We look for n ∈ ℕ_{0} such that

p(n, q_{1}^{n}, …, q_{m}^{n}) = 0.  (5)

We present a recursive algorithm for finding an upper bound for n. Once the bound is known, all integers between zero and the bound can be checked.

In equation (5), the coefficients are elements of 𝔽, which are rational functions of q_{1}, …, q_{m}. We can clear the denominators and obtain an equation in which q_{i} occur polynomially:

r(n, q_{1}, …, q_{m}, q_{1}^{n}, …, q_{m}^{n}) = 0,  (6)

where r ∈ ℚ[x, z_{1}, …, z_{m}, y_{1}, …, y_{m}] \ {0}. We show how to reduce recursively the problem of finding an upper bound for solutions of (6). Consider all terms of r with lowest degree of y_{m}, and let that degree be j. Among these terms, consider the one with the lowest degree of z_{m}, and let d be that degree. The term has the form sz_{m}^{d}y_{m}^{j} for some s ∈ ℚ[x, z_{1}, …, z_{m−1}, y_{1}, …, y_{m−1}] \ {0}. Let M be an upper bound on natural solutions of equation

s(n, q_{1}, …, q_{m−1}, q_{1}^{n}, …, q_{m−1}^{n}) = 0,  (7)

which we can get recursively. Then max(M, d) is an upper bound for solutions of (6). Suppose n > max(M, d). Then n is not a solution of (7), and the lowest power of q_{m} that occurs in (6) is d + nj. Since this power occurs only in the term s(n)q_{m}^{d}q_{m}^{nj}, the term does not cancel, and n is not a solution of (6).

The base case of the recursion is an equation r(n) = 0, where r ∈ ℚ[x] \ {0}. This can be handled easily, since any natural solution of this equation must divide the constant term (after we have cleared the denominators).

3.3 Rational bases

Suppose q_{1}, …, q_{m} ∈ ℚ. Let p ∈ ℚ[x, y] \ {0}. Again we consider the problem of finding n ∈ ℕ_{0} such that

p(n, q_{1}^{n}, …, q_{m}^{n}) = 0.  (8)

Write p as

p = ∑_{i=1}^{k} p_{i}u_{i},

where p_{1}, …, p_{k} ∈ ℚ[x] \ {0}, u_{1}, …, u_{k} ∈ M, and u_{i} ≠ u_{j} for i ≠ j. Equation (8) can be written as

∑_{i=1}^{k} p_{i}(n)u_{i}(q)^{n} = 0.  (9)

Because bases q_{1}, …, q_{m} satisfy condition (1), |u_{i}(q)| ≠ |u_{j}(q)| for i ≠ j. Let s_{i} = u_{i}(q) for i = 1, …, k. We may assume that |s_{1}| < |s_{2}| < ⋯ < |s_{k}|. Suppose p_{k}(x) = a_{d}x^{d} + a_{d−1}x^{d−1} + ⋯ + a_{0}. Equation (9) is equivalent to

a_{d}n^{d}s_{k}^{n} + ∑_{i=0}^{d−1} a_{i}n^{i}s_{k}^{n} + ∑_{i=1}^{k−1} p_{i}(n) s_{i}^{n} = 0.  (10)

The first term in (10) dominates the sum of the others. We only need a lower bound on n, such that the absolute value of the first term is greater than the absolute value of the sum of the other terms. Then we can check all integers between zero and the lower bound.

Let dom(a, b, k) be a function which gives an integer lower bound, such that for all n ≥ dom(a, b, k) it is true that a^{n} > bn^{k}. Here a > 1, b > 0 and k ∈ ℤ.

Let δ = 1/(d + k). For i = 0, …, d − 1, define

M_{i} = ⌈| \frac{a_{i}}{a_{d}δ} |^{\frac{1}{d−i}} ⌉ .

Let K_{i} be the maximum absolute value of the coefficients of p_{i}. For i = 1, …, k − 1, define

N_{i} = dom (| \frac{s_{k}}{s_{i}} | , | \frac{2K_{i}}{δa_{d}} | , deg p_{i} − d) .

Let N = max(2, M_{0}, …, M_{d−1}, N_{1}, …, N_{k−1}). The choice of M_{i} ensures that

|δa_{d}n^{d}s_{k}^{n}| > |a_{i}n^{i}s_{k}^{n}|

for all n ≥ N. The choice of K_{i} ensures that |p_{i}(n)| < 2K_{i}n^{deg p_{i}} for all n ≥ 2. Therefore,

|δa_{d}n^{d}s_{k}^{n}| > |p_{i}(n)s_{i}^{n}|

for all n > N. This means that equation (8) does not have any solutions larger than N. We can find all solutions of (8) by checking all integers between 0 and N.

4 Polynomial solutions

In this section we present an algorithm for finding all polynomial solutions f ∈ 𝔽[x, y] of parametric nonhomogeneous equations of the form

Lf = g + ∑_{j=1}^{s} λ_{j}h_{j}  (11)

where

L = ∑_{i=0}^{ρ} p_{i}E^{i}  (12)

is a linear recurrence operator with polynomial coefficients p_{i} ∈ 𝔽[x, y], λ_{j} are free parameters (ranging over 𝔽) to be determined, and g, h_{j} ∈ 𝔽[x, y] are given polynomials. More precisely, the problem is to compute a basis of the affine space L_{p}^{−1}(g) where L_{p} : 𝔽[x, y] ⊕ 𝔽^{s} → 𝔽[x, y] and L_{p} : (f, λ) ↦ Lf − ∑_{j=1}^{s} λ_{j}h_{j} for f ∈ 𝔽[x, y] and λ ∈ 𝔽^{s}. Thus by a solution of (11) we mean a pair (f, λ) with f ∈ 𝔽[x, y] and λ ∈ 𝔽^{s} such that (11) is satisfied.

As a special case, (11) includes nonhomogeneous equations without parameters (when all h_{j} = 0) as well as homogeneous equations (when also g = 0). The ability to solve parametric nonhomogeneous equations is crucial if one wants to apply Zeilberger’s Creative Telescoping algorithm [26] in the mixed hypergeometric case. Another reason for allowing linear parameters in the equation is the nature of our algorithm which finds the terms of the solution one by one, introducing new free parameters into the right-hand side at each step.

Let f(x, y) be a polynomial solution of (11). Write

α = max_{0≤i≤ρ} deg_{y} p_{i},  (13)
p_{i,α}(x) = [y^{α}] p_{i},  (14)
d = max_{0≤i≤ρ} deg_{x} p_{i,α}(x),  (15)
p_{i,α,d} = [x^{d}] p_{i,α}(x),  (16)
rhs(λ) = g + ∑_{j=1}^{s} λ_{j}h_{j},  (17)
φ = deg_{y} f(x, y),  (18)
t(x) = [y^{φ}] f(x, y),  (19)
β = deg_{y} rhs(λ),  (20)
r_{β} = [y^{β}] rhs(λ),  (21)

where t ∈ 𝔽[x] \ {0}, p_{i,α} ∈ 𝔽[x] and p_{i,α,d} ∈ 𝔽. In (20) we regard λ_{j}’s as variables over 𝔽(x), and rhs(λ) as belonging to 𝔽(x, λ)[y]. This means that after the parameters λ_{j} are given specific values λ′_{j} ∈ 𝔽, the multidegree of rhs(λ′) = g + ∑_{j=1}^{s} λ′_{j}h_{j} in y can be lower than β.

Lemma 4.1 Let L, α, p_{i,α,d} and φ be as given in (12)–(18). If deg_{y} Lf ≺ α + φ, then φ = deg_{y} f satisfies P (q^{φ}) = 0 where

P (x) = ∑_{i=0}^{ρ} p_{i,α,d}x^{i}  (22)

is the characteristic polynomial of L.

Proof: From (19), E^{i}f = t(x + i)q^{iφ}y^{φ} + o(y^{φ}), so Lf = T (x)y^{α+φ} + o(y^{α+φ}) where

T (x) = ∑_{i=0}^{ρ} p_{i,α}(x)q^{iφ}t(x + i).  (23)

If deg_{y} Lf ≺ α + φ then T = 0. This is an ordinary recurrence relation with non-zero polynomial solution t(x). As the coefficient of x^{d+deg_{x} t} in T (x) must vanish, ∑_{i=0}^{ρ} p_{i,α,d}q^{iφ} = 0 as claimed.

Let R denote the set of exponents of those roots of the characteristic polynomial (22) (if any) which are power products of the bases:

R = {σ ∈ ℕ_{0}^{m}; P (q^{σ}) = 0}.  (24)

When R is empty we take max R = ⊥. The following lemma gives rise to an algorithm for finding all polynomial solutions of equation (11).

Lemma 4.2 Let (f, λ′) be a solution of (11) with f ∈ 𝔽[x, y] and λ′ ∈ 𝔽^{s}.

1. If α + max R ≻ β then deg_{y} f ⪯ max R.

2. Let α + max R ⪯ β.

(a) If α ⊆ β then deg_{y} f ⪯ β − α.

(b) If α ⊈ β then deg_{y} rhs(λ′) ≺ β.

Proof: Let φ = deg_{y} f. Let T be as in (23).

1. α + max R ≻ β If T = 0 then deg_{y} Lf ≺ α + φ and φ ∈ R, by Lemma 4.1. If T ≠ 0 then deg_{y} Lf = α + φ. As deg_{y} rhs(λ′) ⪯ β, it follows that α + φ ⪯ β ≺ α + max R, so φ ≺ max R. In either case, φ ⪯ max R as claimed.

2. α + max R ⪯ β

(a) α ⊆ β If T = 0 then deg_{y} Lf ≺ α + φ and φ ∈ R by Lemma 4.1, so φ ⪯ max R and therefore α + φ ⪯ β. If T ≠ 0 then deg_{y} Lf = α + φ. As deg_{y} rhs(λ′) ⪯ β, it follows that α + φ ⪯ β. In either case, φ ⪯ β − α as claimed.

(b) α ⊈ β Assume that deg_{y} rhs(λ′) = β. If T = 0 then deg_{y} Lf ≺ α + φ and φ ∈ R by Lemma 4.1, so deg_{y} Lf ≺ α + φ ⪯ α + max R ⪯ β = deg_{y} rhs(λ′), a contradiction. If T ≠ 0 then α + φ = β which implies that α ⊆ β, contrary to the assumption. Both cases lead to contradiction, so deg_{y} rhs(λ′) ≺ β as claimed.

Based on Lemma 4.2, we can find the general solution (f, λ) of equation (11) as follows: First compute the set R as given in (24). Then distinguish three cases:

1. α + max R ≻ β Set φ = max R and look for f in the form

f = t(x)y^{φ} + f_{1}  (25)

where f_{1} = o(y^{φ}). To find t(x), apply the algorithm of [2] to T = 0 (an ordinary homogeneous recurrence relation). Then remove max R from R and find f_{1} recursively by solving

Lf_{1} = rhs(λ) − L (t(x)y^{φ}) .  (26)

2. α + max R ⪯ β

(a) α ⊆ β Set φ = β − α and look for f in the form (25). To find t(x), apply the algorithm of [2] to T = r_{β} (an ordinary parametric nonhomogeneous recurrence relation). Then remove max R from R (only in case that α + max R = β), and find f_{1} recursively by solving (26).

(b) α ⊈ β Let λ = λ′ be the solution of the system of linear algebraic equations for the free parameters λ obtained by equating the coefficients of powers of x in r_{β} to zero. Then find f recursively by solving Lf = rhs(λ′).

Remarks:

1. Note that in steps 1 and 2(a), t(x) can contain new free parameters which are then joined with the existing ones. This explains the need for allowing parameters in the right-hand side of the equation.

2. In step 2(b), the number of free parameters will drop by the rank of the linear system to be solved.

3. If the ordinary recurrence in steps 1 or 2(a) has no polynomial solution, or the linear system in step 2(b) is unsolvable, then the original parametric recurrence has no polynomial solution, and the algorithm terminates unsuccessfully.

4. At each step, either the cardinality of the set R drops, or else it stays the same but the multidegree β = deg_{y} rhs(λ) decreases in the admissible term order ≺. It follows that the pair (cardR, β) decreases in the lexicographic ordering of ℕ_{0} × ℕ_{0}^{m} which uses < in the first component and ≺ in the second component. As every admissible term order is a well-order, this assures termination of the algorithm.

5. Unless the algorithm terminates unsuccessfully, eventually R becomes empty and rhs(λ) becomes 0. Then the only polynomial solution of (11) is f = 0.

An iterative version of this tail-recursive algorithm called MixedPoly is given in Appendix A.

5 A canonical form for rational functions

Let r ∈ 𝔽(x, y) \ {0}. Write r as

r = \frac{u}{v} · \frac{a_{0}}{b_{0}} ,

where u, v ∈ M, a_{0}, b_{0} ∈ 𝔽[x, y], a_{0}b_{0} ⊥ M, ua_{0} ⊥ vb_{0}, and b_{0} is mixed monic (Def. 2.2).

There are finitely many values h ∈ ℕ_{0} such that a_{0} ⊥̸ E^{h}b_{0}. These values are the elements of the dispersion set D(a_{0}, b_{0}) which can be found as described in Section 3. So let D(a_{0}, b_{0}) = {h_{1}, h_{2}, …, h_{N}} where 0 ≤ h_{1} < h_{2} < ⋯ < h_{N}.

Lemma 5.1 Consider the algorithm CanonicalForm in Appendix B. Define h_{N+1} = ∞, and let 0 ≤ k ≤ i, j ≤ N, h ∈ ℕ_{0} and h < h_{k+1}. Then a_{i} ⊥ E^{h}b_{j}.

Proof: Let S = {h_{1}, …, h_{N}}. Suppose h ∉ S. Since a_{i} | a_{0} and b_{j} | b_{0} and a_{0} ⊥ E^{h}b_{0}, it follows that a_{i} ⊥ E^{h}b_{j}.

To prove the lemma for h ∈ S, we use induction on k. When k = 0, there is nothing to prove because there is no h ∈ S such that h < h_{1}. Assume that the lemma holds for all h ∈ S, h < h_{k}. We show that it holds for h = h_{k}. Since a_{i} | a_{k} and b_{j} | b_{k}, it follows that gcd(a_{i}, E^{h_{k}}b_{j}) divides gcd(a_{k}, E^{h_{k}}b_{k}). Furthermore,

gcd(a_{k}, E^{h_{k}}b_{k}) = gcd( \frac{a_{k−1}}{s_{k}} , \frac{E^{h_{k}}b_{k−1}}{s_{k}}) = 1

by the definition of a_{k}, b_{k} and s_{k} in algorithm CanonicalForm. This completes the proof.

Theorem 5.2 Let r ∈ 𝔽(x, y) \ {0}. There exist polynomials a, b, c ∈ 𝔽[x, y] \ {0} such that

1. b, c are mixed monic,

2. c ⊥ M,

3. a ⊥ E^{k}b for all k ∈ ℕ_{0},

4. a ⊥ c,

5. b ⊥ Ec, and

r = \frac{a}{b} · \frac{Ec}{c} .  (27)

Proof: Let a, b, c be constructed by the algorithm CanonicalForm from Appendix B. Conditions 1 and 2 are satisfied by construction, and condition 3 follows from Lemma 5.1 by taking i = j = k = N. Identity (27) is verified directly,

\frac{a}{b} · \frac{Ec}{c} = \frac{u · a_{N}}{v · b_{N}} · ∏_{i=1}^{N} ∏_{j=1}^{h_{i}} \frac{E^{−j+1}s_{i}}{E^{−j}s_{i}} =

= \frac{u · a_{0}}{∏_{i=1}^{N} s_{i}} · \frac{∏_{i=1}^{N} E^{−h_{i}}s_{i}}{v · b_{0}} · ∏_{i=1}^{N} \frac{s_{i}}{E^{−h_{i}}s_{i}} = \frac{u · a_{0}}{v · b_{0}} = r.

Proof of 4: Suppose a ⊥̸ c. Then also a_{N} ⊥̸ E^{−j}s_{i} for some i and j such that 1 ≤ i ≤ N and 1 ≤ j ≤ h_{i}. By definition E^{h_{i}−j}b_{i−1} = E^{h_{i}−j}b_{i} · E^{−j}s_{i}, so it follows that a_{N} ⊥̸ E^{h_{i}−j}b_{i−1}. Since h_{i} − j < h_{i}, this contradicts Lemma 5.1.

Proof of 5: Suppose b ⊥̸ Ec. Then also b_{N} ⊥̸ E^{−j}s_{i} for some i and j such that 1 ≤ i ≤ N and 0 ≤ j ≤ h_{i} − 1. By definition E^{−j}a_{i−1} = E^{−j}a_{i} · E^{−j}s_{i}, so it follows that a_{i−1} ⊥̸ E^{j}b_{N}. Since j < h_{i}, this contradicts Lemma 5.1.

Lemma 5.3 Let a, b, c, A, B, C ∈ 𝔽[x, y] \ {0} be polynomials such that a ⊥ c, b ⊥ Ec, c ⊥ M, and A ⊥ E^{k}B for all k ∈ ℕ_{0}. If

\frac{a}{b} · \frac{Ec}{c} = \frac{A}{B} · \frac{EC}{C} ,  (28)

then c divides C.

Proof: Let g = gcd(c, C), d = c/g, and D = C/g. Then d ⊥ D, a ⊥ d, and b ⊥ Ed. Clear denominators in (28) and cancel g · Eg on both sides. The result

A · b · d · ED = a · B · D · Ed

shows that d | B · Ed and Ed | A · d. Using these two relations repeatedly, we find that

d | B · EB ⋯ E^{k−1}B · E^{k}d
d | E^{−1}A · E^{−2}A ⋯ E^{−k}A · E^{−k}d

for all k ∈ ℕ_{0}. Because d ⊥ M, and 𝔽 has characteristic zero, d ⊥ E^{k}d and d ⊥ E^{−k}d for large enough k. It follows that d divides both B · EB ⋯ E^{k−1}B and E^{−1}A · E^{−2}A ⋯ E^{−k}A for large enough k. But these two polynomials are coprime by assumption, so d must be a constant. Hence, c divides C.

Theorem 5.4 Let r ∈ 𝔽[x, y] \ {0}. Then the factorization of r described in Theorem 5.2 is unique.

Proof: Suppose that a, b, c and A, B, C are two factorizations of r, as described in Theorem 5.2. Then

r = \frac{a}{b} · \frac{Ec}{c} = \frac{A}{B} · \frac{EC}{C} .

By Lemma 5.3, c divides C, and vice versa. As c and C are mixed monic they must be equal, hence a/b = A/B. As a ⊥ b, A ⊥ B, and b, B are mixed monic, it follows that b = B and a = A as well.

The factorization of non-zero rational functions described in Theorem 5.2 is thus a canonical form. We introduce special notation for it.

Definition 5.1 Let r ∈ 𝔽(x, y) \ {0} be a non-zero rational function. We write

(a, b, c) = C.f.(r)

to denote the unique polynomials a, b, c ∈ 𝔽[x, y] which satisfy the conditions of Theorem 5.2.

Theorem 5.5 Let a, b ∈ 𝔽[x, y]\{0}, and (A, B, C) = C.f.(b/a). The homogeneous first-order recurrence

a · Ef − bf = 0  (29)

has a non-zero polynomial solution f ∈ 𝔽[x, y] if and only if A/B = u(q) for some u ∈ M. In that case, f = λ · u · C for some λ ∈ 𝔽 \ {0}.

Proof: Suppose (29) has a non-zero solution f ∈ 𝔽[x, y]. Write f = λ · u · g where λ ∈ 𝔽 \ {0}, u ∈ M and g ⊥ M is mixed monic. Then C.f.(Ef/f) = (u(q), 1, g). Since from (29)

\frac{Ef}{f} = \frac{b}{a} = \frac{A}{B} · \frac{EC}{C} ,

(A, B, C) = C.f.(Ef/f) as well. By Theorem 5.4 it follows that A = u(q), B = 1 and C = g, so A/B = u(q) and f = λ · u · C.

Conversely, if A/B = u(q) for some u ∈ M, then f = u · C is a non-zero solution of (29).

We remark that our canonical form differs from the Paule/Riese/Strehl form (PRS, for short) described in [18, 17] for the basic and in [20] for the bibasic case, in the following three respects:

1. In the PRS form the monomial factors of the numerator and denominator of r are listed separately while in our form they are included with A resp. B.

2. In the PRS form all polynomials either have unit constant terms or else are primitive and the overall constant factor is listed separately, while in our form B and C have unit leading coefficients in the chosen term order and the overall constant factor is included with A.

3. In the PRS form the polynomial corresponding to our B is given by a constant multiple of E^{−1}B.

6 Mixed Gosper’s algorithm

Let S_{n} = ∑_{k=0}^{n−1} t_{k}. Clearly substituting S_{n} for s_{n} satisfies the first-order recurrence

s_{n+1} − s_{n} = t_{n}.  (30)

Conversely, any solution s_{n} of (30) differs from S_{n} only by an additive constant – more precisely, S_{n} = s_{n} − s_{0}. Therefore we consider the following problem:

Given a sequence t_{n}, decide if equation (30) has a mixed hypergeometric solution s_{n}, and if so, find it.

Let s_{n} and t_{n} satisfy (30), with

\frac{s_{n+1}}{s_{n}} =: T_{n} ∈ 𝔽(n, q^{n}).

Then the two quotients

r_{n} := \frac{t_{n+1}}{t_{n}} = \frac{s_{n+2} − s_{n+1}}{s_{n+1} − s_{n}} = \frac{T_{n+1} − 1}{1 − 1/T_{n}}

and

R_{n} := \frac{s_{n}}{t_{n}} = \frac{s_{n}}{s_{n+1} − s_{n}} = \frac{1}{T_{n} − 1}

both belong to 𝔽(n, q^{n}). So t_{n} must be mixed hypergeometric itself, and s_{n} is a rational multiple of t_{n}: s_{n} = R_{n}t_{n}. Using this, (30) yields a recurrence for the unknown rational sequence R_{n},

r_{n}R_{n+1} − R_{n} = 1.  (31)

By Theorem 2.2, equation (31) is equivalent to

r · ER − R = 1,  (32)

where r, R ∈ 𝔽(x, y) are Φ^{−1}-isomorphic images of r_{n} and R_{n}, respectively.

Next we show how to find rational solutions R ∈ 𝔽(x, y) of equation (32). The following theorem provides a multiple of the denominator and a divisor of the numerator of R. The missing factor in the numerator can then be found using algorithm MixedPoly of Section 4.

Definition 6.1 Let r ∈ 𝔽(x, y) \ {0} be a non-zero rational function, and (a, b, c) = C.f.(r). For 1 ≤ i ≤ m define exponents e_{i}(r) as follows: If π_{i}(a)π_{i}(b) ≠ 0, let (a_{i}, b_{i}, c_{i}) = C.f.(π_{i}(b)/π_{i}(a)). If there are v, w ∈ M such that v ⊥ w and a_{i}/b_{i} = v(q)/w(q), then e_{i}(r) = deg_{y_{i}} w. If not, or if π_{i}(a)π_{i}(b) = 0, then e_{i}(r) = 0.

Theorem 6.1 Let R = f/(ug) be a rational solution of (32) with f, g ∈ 𝔽[x, y], u ∈ M, g ⊥ M, and f ⊥ ug. Then

1. g | c where (a, b, c) = C.f.(r),

2. deg_{y_{i}} u ≤ e_{i}(r),

3. E^{−1}b | f.

Proof:

1. From (32),

r = \frac{R + 1}{ER} = \frac{(f + ug)u(q)}{Ef} · \frac{Eg}{g} .  (33)

On the other hand, (a, b, c) = C.f.(r), so

r = \frac{a}{b} · \frac{Ec}{c} .  (34)

As Ef ⊥ Eg, g ⊥ (f + ug)u(q), g ⊥ M, and a ⊥ E^{k}b for all k ∈ ℕ_{0}, it follows by Lemma 5.3 that g divides c.

2. Write F = fc/g ∈ 𝔽[x, y]. Then R = f/(ug) = F/(uc). Combining this with (33) and (34), we find

that

(F + uc) · u(q) · b = a · EF.  (35)

Now assume that y_{i} | u. Then applying π_{i} (see Def. 2.3) to Eqn. (35) and rearranging yields

π_{i}(a) · Eπ_{i}(F) − u(q) · π_{i}(b) · π_{i}(F) = 0.  (36)

Because F | fc, f ⊥ u and c ⊥ M, it follows that y_{i} |̸ F and π_{i}(F) ≠ 0. Assume that π_{i}(a) = 0. Then y_{i} | a and from (35), y_{i} | b · F. But y_{i} |̸ b because a ⊥ b, so y_{i} | F. This contradiction shows that π_{i}(a) ≠ 0. In an analogous way we conclude that π_{i}(b) ≠ 0. Let (a_{i}, b_{i}, c_{i}) = C.f.(π_{i}(b)/π_{i}(a)). Then ((u(q)·a_{i}, b_{i}, c_{i}) = C.f.(u(q)·π_{i}(b)/π_{i}(a)). Since equation (36) has a non-zero polynomial solution π_{i}(F), it follows by Theorem 5.5 that there is u_{1} ∈ M such that u(q) · a_{i}/b_{i} = u_{1}(q), and that π_{i}(F) = λu_{1}c_{i} for some λ ∈ 𝔽. Then a_{i}/b_{i} = u_{1}(q)/u(q) is a quotient of two monomials. Write u_{1} = v · t and u = w · t where t, v, w ∈ M and v ⊥ w. By Definiton 6.1, e_{i}(r) = deg_{y_{i}} w. As t | u_{1} | π_{i}(F) it follows that t ⊥ y_{i}, so

deg_{y_{i}} u = deg_{y_{i}} w = e_{i}(r).

We have shown that deg_{y_{i}} u is either 0 or e_{i}(r), so in either case deg_{y_{i}} u ≤ e_{i}(r).

3. From (35) it follows that b | a·EF. As a ⊥ b, we have that b | EF | EfEc. But b ⊥ Ec, so E^{−1}b | f.

From Theorem 6.1 it follows that we can look for R in the form

R = \frac{E^{−1}b · p}{u · c}  (37)

where (a, b, c) = C.f.(r) and u = ∏_{i=1}^{m} y_{i}^{e_{i}(r)} are known while p ∈ 𝔽[x, y] is an unknown polynomial. Inserting (37) and (34) into (32) yields

a · Ep − u(q) · E^{−1}b · p = u(q)u · c,  (38)

an nonhomogeneous first-order linear recurrence relation with polynomial coefficients satisfied by p. Algorithm MixedPoly of Section 4 can now be applied to find a polynomial solution p of Eqn. (38). The full algorithm is given in Appendix C.

We conclude this section by giving some examples of sums which can be evaluated automatically by MixedGosper. We write (a; q)_{n} = ∏_{i=0}^{n−1}(1 − aq^{i}).

Many bibasic examples can be found in [7] and [20]. An indefinite multibasic summation formula (too big to reproduce it here) is proved in [21]. The formula contains an arbitrary number, k, of bases. Such formulæ cannot be proved by our algorithm. However, any specialization of this formula in which k is replaced by a specific natural number can be, at least in principle, not only proved, but also derived by MixedGosper. In [21], it is shown that several well-known basic and bibasic summation formulæ can be obtained as specializations of this k-basic master formula.

The following two examples are due to Gosper [10].

Example 6.1 In this tribasic example 𝔽 = ℚ(a, b, c, p, q, r) where a, b, c are parameters and p, q, r are the bases.

∑_{k=0}^{n} \frac{(−a)^{k}p^{\binom{k}{2}}(1 − abp^{k}q^{k})(1 − acp^{k}r^{k})(b; q)_{k}(c; r)_{k}}{(ap; p)_{k}(abcpqr; pqr)_{k}}

= (a − 1)(abc − 1) + \frac{(−a)^{n}ap^{\binom{n+1}{2}}(b; q)_{n+1}(c; r)_{n+1}}{(ap; p)_{n}(abcpqr; pqr)_{n}} .

Example 6.2 In this quadbasic example 𝔽 = ℚ(a, b, c, d, p, q, r, s) where a, b, c, d are parameters and p, q, r, s are the bases.

∑_{k=0}^{n} \frac{b^{k}q^{\binom{k+1}{2}}(1 − \frac{ap^{k}}{bq^{k}})(1 − abp^{k}q^{k})(1 − \frac{ds^{k}}{cr^{k}})(1 − cdr^{k}s^{k})( \frac{a}{c}; \frac{p}{r})_{k}(ac; pr)_{k}( \frac{d}{b}; \frac{s}{q})_{k}(bd; qs)_{k}}{a^{k}p^{\binom{k+1}{2}}( \frac{bq}{cr}; \frac{q}{r})_{k}(bcqr; qr)_{k}( \frac{ds}{ap}; \frac{s}{p})_{k}(adps; ps)_{k}}

= (b − c)(1 − \frac{1}{bc})(a − d)(ad − 1) +

\frac{b^{n}q^{\binom{n}{2}}(cr^{n} − ap^{n})(acp^{n}r^{n} − 1)(bq^{n} − ds^{n})(bdq^{n}s^{n} − 1)( \frac{a}{c}; \frac{p}{r})_{n}(ac; pr)_{n}( \frac{d}{b}; \frac{s}{q})_{n}(bd; qs)_{n}}{a^{n}bcp^{\binom{n+1}{2}}r^{n}( \frac{bq}{cr}; \frac{q}{r})_{n}(bcqr; qr)_{n}( \frac{ds}{ap}; \frac{s}{p})_{n}(adps; ps)_{n}} .

Example 6.3 For a simple mixed hypegeometric example, define

φ_{q,n}(a, b, c) = ∏_{i=0}^{n−1}(a + bi + cq^{i}).

Then, using MixedGosper, one obtains

∑_{k=0}^{n}(bk + cq^{k})φ_{q,k}(1, b, c) = φ_{q,n+1}(1, b, c) − 1.  (39)

As φ_{q,n}(1, 1, 0) = n! and φ_{q,n}(1, 0, −a) = (a; q)_{n}, both well-known summation formulæ

∑_{k=0}^{n} k k! = (n + 1)! − 1,  ∑_{k=0}^{n}(a; q)_{k}q^{k} = \frac{1 − (a; q)_{n+1}}{a} ,

turn out to be special cases of (39). – It is not hard to imagine how more complex mixed hypergeometric formulæ could be built using φ or similar functions.

7 Mixed hypergeometric solutions

In this section we derive algorithm MixedHyper for finding all mixed hypergeometric solutions f of Lf = 0 where L is as in (12). Let Ef = rf where r ∈ 𝔽(x, y), then E^{i}f = ∏_{j=0}^{i−1} (E^{j}r) f. The crucial idea is to look for r in the canonical form described in Theorem 5.2. More precisely, we use a slightly modified canonical form

r = z \frac{a}{b} \frac{Ec}{c}  (40)

where z ∈ 𝔽 \ {0}, a ∈ 𝔽[x, y] is mixed monic, and a, b, c satisfy conditions 1–5 of Theorem 5.2. After inserting (40) into Lf = 0, clearing denominators and cancelling f we obtain

∑_{i=0}^{ρ} z^{i}P_{i} · E^{i}c = 0  (41)

where

P_{i} = p_{i} ∏_{j=0}^{i−1} E^{j}a ∏_{j=i}^{ρ−1} E^{j}b.

Since all terms in (41) except for i = 0 contain a as an explicit factor, it follows that a divides p_{0}c ∏_{j=0}^{ρ−1} E^{j}b. Because of properties 3 and 4 of the canonical form, a divides p_{0}. Similarly, all terms in (41) except for i = ρ contain E^{ρ−1}b as an explicit factor, therefore E^{ρ−1}b divides z^{ρ}p_{ρ}E^{ρ}c ∏_{j=0}^{ρ−1} E^{j}a. Because of properties 3 and 5 of the canonical form, E^{ρ−1}b divides p_{ρ}. Thus we have a finite number of choices for a and b: they are mixed monic factors of p_{0} and E^{1−ρ}p_{ρ}, respectively.

For each choice of a and b, equation (41) is a linear recurrence with polynomial coefficients satisfied by the unknown polynomial c(x, y). However, z ∈ 𝔽 \ {0} is also unknown. To find z(Ec/c), write α = min_{0≤i≤ρ} mindeg_{y}P_{i}, p_{i,α}(x) = [y^{α}] P_{i}, d = max_{0≤i≤ρ} deg_{x} p_{i,α}(x), p_{i,α,d} = [x^{d}] p_{i,α}(x), and φ = mindeg_{y}c(x, y). Looking at the coefficient of y^{α+φ} in (41), we find that P (zq^{φ}) = 0 where

P (x) = ∑_{i=0}^{ρ} p_{i,α,d}x^{i}.  (42)

Write zq^{φ} = τ where τ is a root of P. Then z = τq^{−φ}, hence (41) can be rewritten as

∑_{i=0}^{ρ} τ^{i}q^{−iφ}P_{i} · E^{i}c = 0.

Dividing by y^{φ} we obtain finally

∑_{i=0}^{ρ} τ^{i}P_{i} · E^{i}(c/y^{φ}) = 0.

Suppose that we know c̄ = c/y^{φ}. Then

τ \frac{Ec̄}{c̄} = τq^{−φ} \frac{Ec}{c} = z \frac{Ec}{c}

which is just what we are looking for. It remains to see how to find all Laurent polynomial solutions c̄ ∈ 𝔽[x, y, y^{−1}] of

∑_{i=0}^{ρ} τ^{i}P_{i} · E^{i}(c̄) = 0.  (43)

Obviously, if we know lower bounds b_{j} for the degrees of y_{j} in c̄, then we can substitute c̄ = d̄y_{1}^{b_{1}} ⋯ y_{m}^{b_{m}} in (43) and use MixedPoly to find polynomial solutions d̄ of the resulting equation. Observe that any lexicographic order ⪯ on ℤ^{m} is total (though not well-founded) and satisfies α ⪯ β ⇒ α + γ ⪯ β + γ for all α, β, γ ∈ ℤ^{m}. To obtain b_{j}, we order the terms lexicographically with y_{j} as the first variable, and write once again α_{j} = min_{0≤i≤ρ} mindeg_{y}P_{i}, p_{i,α_{j}} (x) = [y^{α^{j}}] τ^{i}P_{i}, d_{j} = max_{0≤i≤ρ} deg_{x} p_{i,α_{j}} (x), p_{i,α_{j} ,d_{j}} = [x^{d_{j}}] p_{i,α_{j}} (x), and φ_{j} = mindeg_{y}c̄(x, y). Then P^{(j)}(q^{φ^{j}}) = 0 where

P^{(j)}(x) = ∑_{i=0}^{ρ} p_{i,α_{j} ,d_{j}} x^{i}.

The bound b_{j} can now be read off as the j-th component of φ_{j}.

In summary, we find the factors of r = τ(a/b)(Ec̄/c̄) as follows:

1. a is a mixed monic factor of p_{0},

2. b is a mixed monic factor of E^{1−ρ}p_{ρ},

3. τ is a root of polynomial P (x) defined in (42),

4. c̄ is a non-zero Laurent polynomial solution of (43).

Checking each admissible triple a, b, τ for Laurent polynomial solutions c̄ of recurrence (43) constitutes algorithm MixedHyper which is given in appendix D.

8 Greatest factorial factorization

The concept of greatest factorial factorization of polynomials (GFF, for short) which is an analogue of the well-known square-free factorization (SFF), plays a fundamental role in symbolic summation. It has been introduced by Paule in [16] for the hypergeometric case, and subsequently extended to the basic [18, 17] as well as bibasic cases [20]. Here we sketch an extension of the GFF concept to an arbitrary polynomial ring with an automorphism σ, including the hypergeometric, basic, bibasic, multibasic, and mixed hypergeometric GFF, as well as SFF, as special cases.

Let 𝔽 be a field of characteristic zero and 𝔽[x] = 𝔽[x_{1}, x_{2}, …, x_{n}] the ring of n-variate polynomials over 𝔽. For p, q ∈ 𝔽[x], we write p ∼ q if there is an a ∈ 𝔽 \ {0} such that p = aq. Such p and q are called associated.

Let σ be an 𝔽-automorphism of 𝔽[x] (i.e., a ring automorphism of 𝔽[x] which fixes each element of 𝔽 ⊆ 𝔽[x]). To specify σ it suffices to give the n polynomials σx_{1}, …, σx_{n}. Note that σ preserves irreducibility of polynomials, and so for any irreducible p ∈ 𝔽[x] \ 𝔽, either σp ∼ p or σp ⊥ p.

In analogy to Moenck [15] we write

[p]\underline{k} ; σ = ∏_{i=0}^{k−1} σ^{−i}p

for the k-th falling σ-factorial of p.

Definition 8.1 Let p ∈ 𝔽[x] \ {0}. Then

σ-span(p) = max\{k ∈ ℕ_{0}; [q]\underline{k} ; σ divides p, for some q ∈ 𝔽[x] \ 𝔽\}.

Note that 0 ≤ σ-span(p) ≤ max_{1≤i≤n} deg_{x_{i}}p.

Definition 8.2 Let p ∈ 𝔽[x] \ {0}. A list of polynomials ⟨p_{1}, …, p_{k}⟩ from 𝔽[x] where k ≥ 0 is a σ-GFF of p if the following conditions are satisfied:

(σ-GFF1) p ∼ [p_{1}]\underline{1} ; σ ⋯ [p_{k}]\underline{k} ; σ,

(σ-GFF2) if k > 0 then p_{k} ∈ 𝔽[x] \ 𝔽,

(σ-GFF3) for all 1 ≤ i ≤ k, either p_{i} = 1, or p_{i} ∈ 𝔽[x] \ 𝔽 and σ-span ([p_{1}]\underline{1} ; σ ⋯ [p_{i}]\underline{i} ; σ) = i.

It is clear that every polynomial p ∈ 𝔽[x] \ {0} has a σ-GFF. To find it, factor p into non-constant irreducible factors, combine those which compose falling σ-factorials of length k = σ-span(p) into [p_{k}]\underline{k} ; σ, then repeat this procedure with p/[p_{k}]\underline{k} ; σ. Note that the σ-GFF of a constant polynomial p ∈ 𝔽 \ \{0\} is ⟨ ⟩, the empty list.

Example 8.1 1. When σ is the identity automorphism, σ-GFF agrees with SFF.

2. When the polynomial ring is 𝔽[x] and σx = x + 1, σ-GFF agrees, up to a constant normalization factor, with GFF as defined in [16].

3. Let 𝔽 = 𝕂(q) where 𝕂 is a field of characteristic zero and q is transcendental over 𝕂. When the polynomial ring is 𝔽[x] and σx = qx, σ-GFF agrees (on polynomials which are not divisible by x), up to a constant normalization factor, with qGFF as defined in [17].

4. Let 𝔽 = 𝕂(p, q) where 𝕂 is a field of characteristic zero and p, q are algebraically independent over 𝕂. When the polynomial ring is 𝔽[x, y] and σx = qx, σy = py, σ-GFF agrees (on polynomials divisible by neither x nor y), up to a constant normalization factor, with GFF_{p,q} as defined in [20].

5. When the polynomial ring is 𝔽[x, y], σ = E as defined in Section 2, and q_{1}, …, q_{m} satisfy (1), σ-GFF provides GFF for the mixed hypergeometric case.

Note that cases 2, 3 and 4 are all contained in case 5.

As shown by the following example, σ-GFFs are in general not unique.

Example 8.2 Let σ be the 𝔽-automorphism of 𝔽[x] defined by σx = −x. Obviously both ⟨1, 1 + x⟩ and ⟨1, 1 − x⟩ are σ-GFFs of 1 − x^{2} in this case.

In order to have unique σ-GFF’s (up to associated factors), we restrict our attention to a subclass of automorphisms σ.

Definition 8.3 An 𝔽-automorphism σ of 𝔽[x] is aperiodic if for every irreducible p ∈ 𝔽[x] \ 𝔽 either a) σp ∼ p, or

b) σ^{m}p ⊥ p, for all m ∈ ℤ \ {0}.

Example 8.3 1. The identity automorphism is aperiodic because in this case σp ∼ p for all p ∈ 𝔽[x].

2. Let q ∈ 𝔽 be a primitive m-th root of unity. Then the 𝔽-automorphism σ of 𝔽[x] defined by σx = qx is not aperiodic because, e.g., σ(x + 1) = qx + 1 ≁ x + 1 while σ^{m}(x + 1) = q^{m}x + 1 ∼ x + 1.

3. If q_{1}, …, q_{m} ∈ 𝔽 satisfy (1) it follows from Lemma 2.3 that the 𝔽-automorphism σ = E of 𝔽[x] as defined in Section 2 is aperiodic. This includes cases 2 – 5 of Example 8.1.

Lemma 8.1 Assume that σ is an aperiodic 𝔽-automorphism of 𝔽[x]. Let q ∈ 𝔽[x], and let p_{1}, …, p_{k−1} ∈ 𝔽[x] \ {0}, p_{k} ∈ 𝔽[x] \ 𝔽 satisfy

(D1) σp_{j}, σ^{−j}p_{j} ⊥ [p_{i}]\underline{i} ; σ, for any 1 ≤ i < j ≤ k,

(D2) rσ^{j}r does not divide p_{i}, for any r ∈ 𝔽[x] \ 𝔽 and 1 ≤ j ≤ i ≤ k. If [q]\underline{k} ; σ divides [p_{1}]\underline{1} ; σ ⋯ [p_{k}]\underline{k} ; σ, then q divides p_{k}.

Proof: Assume first that q is irreducible. As σ is aperiodic, we distinguish two cases.

a) σq ∼ q In this case [q]\underline{k} ; σ ∼ q^{k} and so q^{k} | [p_{1}]\underline{1} ; σ ⋯ [p_{k}]\underline{k} ; σ. Let i ≤ k be minimal such that q | p_{i}. Assume that i < k. If there is j > i such that q | p_{j} then σp_{j} ⊥̸ [p_{i}]\underline{i} ; σ, contrary to (D1). Otherwise q^{k} | [p_{i}]\underline{i} ; σ and hence q^{2} | p_{i}, contrary to (D2). So in case a) we must have i = k. It follows that q | p_{k}.

b) σ^{m}q ⊥ q, for all m ∈ ℤ \ {0} Let j ≤ k be maximal such that [q]\underline{k} ; σ ⊥̸ [p_{j}]\underline{j} ; σ. Then there are u and v, 0 ≤ u < k and 0 ≤ v < j, such that σ^{−u}q | σ^{−v}p_{j}. Hence σ^{v−u}q | p_{j}. We are going to prove that u = v and j = k, by distinguishing two subcases.

b1) u > v

In this case σ^{v−u+1} | [q]\underline{k} ; σ | [p_{1}]\underline{1} ; σ ⋯ [p_{j}]\underline{j} ; σ. If σ^{v−u+1} | [p_{l}]\underline{l} ; σ for some l < j then σp_{j} ⊥̸ [p_{l}]\underline{l} ; σ, contrary

to (D1). Otherwise σ^{v−u+1} | [p_{j}]\underline{j} ; σ and so σ^{v−u+1} | σ^{−w}p_{j} for some w such that 0 ≤ w < j. Hence r, σ^{w+1}r | p_{j} where r = σ^{v−u}q. As 1 ≤ w + 1 ≤ j, it follows from b) that r ⊥ σ^{w+1}r whence rσ^{w+1}r | p_{j}, contrary to (D2). So this case is not possible.

b2) u ≤ v

Unless u = v and j = k, we have σ^{v−u−j} | [q]\underline{k} ; σ | [p_{1}]\underline{1} ; σ ⋯ [p_{j}]\underline{j} ; σ. If σ^{v−u−j} | [p_{l}]\underline{l} ; σ for some l < j then

σ^{−j}p_{j} ⊥̸ [p_{l}]\underline{l} ; σ, contrary to (D1). Otherwise σ^{v−u−j} | [p_{j}]\underline{j} ; σ and so σ^{v−u−j} | σ^{−w}p_{j} for some w such that 0 ≤ w < j. Hence r, σ^{j−w}r | p_{j} where r = σ^{w+v−u−j}q. As 1 ≤ j − w ≤ j, it follows from b) that r ⊥ σ^{j−w}r whence rσ^{j−w}r | p_{j}, contrary to (D2). So in case b) we must have u = v and j = k. It follows that q | p_{k}.

Finally, if q is reducible write q = q_{1} ⋯ q_{m} where m > 1 and the q_{i}’s are irreducible. As [q_{m}]\underline{k} ; σ | [q]\underline{k} ; σ | [p_{1}]\underline{1} ; σ ⋯ [p_{k}]\underline{k} ; σ and q_{m} is irreducible, we already know that q_{m} | p_{k}. Therefore

[q_{1} ⋯ q_{m−1}]\underline{k} ; σ | [p_{1}]\underline{1} ; σ ⋯ [p_{k−1}]\underline{k−1} ; σ[ \frac{p_{k}}{q_{m}}]\underline{k} ; σ.

Clearly, replacing p_{k} by p_{k}/q_{m} invalidates neither (D1) nor (D2). Hence, by induction on m, q_{1} ⋯ q_{m−1} | p_{k}/q_{m}. It follows that q = q_{1} ⋯ q_{m} | p_{k}.

We can now give a characterization of σ-GFFs which is akin to the definition of GFF in [16]. Corollary 8.2 Assume that σ is an aperiodic 𝔽-automorphism of 𝔽[x]. Let p = [p_{1}]\underline{1} ; σ ⋯ [p_{k}]\underline{k} ; σ where p_{1}, …, p_{k−1} ∈ 𝔽[x] \ {0} and p_{k} ∈ 𝔽[x] \ 𝔽. Then ⟨p_{1}, …, p_{k}⟩ is a σ-GFF for p if and only if p_{1}, …, p_{k} satisfy conditions (D1) and (D2) of Lemma 8.1. Proof: (⇒) Let ⟨p_{1}, …, p_{k}⟩ be a σ-GFF for p.

(D1) Assume that q | [p_{i}]\underline{i} ; σ where q ∈ 𝔽[x] \ 𝔽 is irreducible and 1 ≤ i < j ≤ k. If q | σp_{j} then

[σ^{−1}q]\underline{j} ; σ | [p_{j}]\underline{j} ; σ, so [q]\underline{j+1} ; σ | [p_{i}]\underline{i} ; σ[p_{j}]\underline{j} ; σ. If q | σ^{−j}p_{j} then [σ^{j}q]\underline{j} ; σ | [p_{j}]\underline{j} ; σ, so [σ^{j}q]\underline{j+1} ; σ | [p_{i}]\underline{i} ; σ[p_{j}]\underline{j} ; σ. In either case

σ-span ([p_{1}]\underline{1} ; σ ⋯ [p_{j}]\underline{j} ; σ) ≥ j + 1, contrary to (σ-GFF3).

(D2) Assume that q σ^{j}q | p_{i} for some irreducible q ∈ 𝔽[x]\𝔽 and 1 ≤ j ≤ i ≤ k. Then [q]\underline{i} ; σ[σ^{j}q]\underline{i} ; σ | [p_{i}]\underline{i} ; σ. But σq | [σ^{j}q]\underline{i} ; σ, so [σq]\underline{i+1} ; σ | [p_{i}]\underline{i} ; σ. Hence σ-span ([p_{1}]\underline{1} ; σ ⋯ [p_{i}]\underline{i} ; σ) ≥ i + 1, contrary to (σ-GFF3).

(⇐) Let p_{1}, …, p_{k} satisfy conditions (D1) and (D2) of Lemma 8.1. To prove that they satisfy (σ-GFF3) as well, assume that p_{i} ∈ 𝔽[x] \ 𝔽 and σ-span ([p_{1}]\underline{1} ; σ ⋯ [p_{i}]\underline{i} ; σ) = s > i, for some i such that

1 ≤ i ≤ k. Then there is an irreducible q ∈ 𝔽[x]\𝔽 such that [q]\underline{s} ; σ | [p_{1}]\underline{1} ; σ ⋯ [p_{i}]\underline{i} ; σ. Hence [q]\underline{i} ; σ | [p_{1}]\underline{1} ; σ ⋯ [p_{i}]\underline{i} ; σ and [σ^{−1}q]\underline{i} ; σ | [p_{1}]\underline{1} ; σ ⋯ [p_{i}]\underline{i} ; σ. By Lemma 8.1, q | p_{i} and σ^{−1}q | p_{i}. We distinguish two cases.

a) σq ∼ q If there is j < i such that q | p_{j} then σp_{i} ⊥̸ [p_{j}]\underline{j} ; σ, contrary to (D1). Otherwise q^{s} | [p_{i}]\underline{i} ; σ. As s > i, it follows that q^{2} | p_{i}, contrary to (D2).

b) σ^{m}q ⊥ q, for all m ∈ ℤ \ {0} In this case q ⊥ σ^{−1}q and so q σ^{−1}q | p_{i}, contrary to (D2).

Example 8.4 Let q ∈ 𝔽 be a primitive third root of unity, and let σ be the 𝔽-automorphism of 𝔽[x] defined by σx = qx. The polynomials p_{1} = 1 + x, p_{2} = 1, p_{3} = 1 + qx satisfy conditions (D1) and (D2) of Lemma 8.1. However, ⟨p_{1}, p_{2}, p_{3}⟩ is not a σ-GFF of p = [p_{1}]\underline{1} ; σ[p_{2}]\underline{2} ; σ[p_{3}]\underline{3} ; σ because p = [1 + x]\underline{4} ; σ and σ-span(p) = 4. This example shows that Corollary 8.2 can fail when σ is not aperiodic.

The next corollary shows uniqueness of σ-GFFs (up to associated factors). Corollary 8.3 Assume that σ is an aperiodic 𝔽-automorphism of 𝔽[x]. If ⟨p_{1}, …, p_{k}⟩ and ⟨q_{1}, …, q_{l}⟩ are σ-GFFs for the same p ∈ 𝔽[x], then k = l and p_{i} ∼ q_{i} for 1 ≤ i ≤ k. Proof: By (σ-GFF3), k = σ-span(p) = l.

We prove the rest of the corollary by induction on k.

If k = 0 the assertion holds vacuously. Let k > 0. By Corollary 8.2, the p_{i}’s as well as the q_{i}’s satisfy conditions (D1) and (D2) of Lemma 8.1. Therefore p_{k} | q_{k} and q_{k} | p_{k}, hence p_{k} ∼ q_{k} and

[p_{1}]\underline{1} ; σ ⋯ [p_{k−1}]\underline{k−1} ; σ ∼ [q_{1}]\underline{1} ; σ ⋯ [q_{k−1}]\underline{k−1} ; σ. Let r = [p_{1}]\underline{1} ; σ ⋯ [p_{k−1}]\underline{k−1} ; σ and m = σ-span(r). Then ⟨p_{1}, …, p_{m}⟩

and ⟨q_{1}, …, q_{m}⟩ are σ-GFFs for r, and p_{i} = q_{i} = 1 for m < i < k. By the inductive hypothesis, p_{i} ∼ q_{i} for 1 ≤ i ≤ m.

Finally, the formulæ

gcd(p, σp) ∼ [p_{1}]\underline{0} ; σ ⋯ [p_{k}]\underline{k−1} ; σ gcd(p_{1}, σp_{1}) ⋯ gcd(σ^{−k+1}p_{k}, σp_{k}),

\frac{p}{gcd(p, σ^{−1}p)} ∼ \frac{p_{1}p_{2} ⋯ p_{k}}{gcd(p_{1}, σ^{−1}p_{1}) ⋯ gcd(p_{k}, σ^{−k}p_{k})} ,
\frac{p}{gcd(p, σp)} ∼ \frac{p_{1}σ^{−1}p_{2} ⋯ σ^{−k+1}p_{k}}{gcd(p_{1}, σp_{1}) ⋯ gcd(σ^{−k+1}p_{k}, σp_{k})} .

can be proved in much the same way as the corresponding formulæ in [16].

We remark that analogously to the hypergeometric, basic and bibasic cases, σ-GFF could be used to derive and explain Gosper’s algorithm in the mixed hypergeometric case.

9 Concluding remarks

We have shown how to compute the hypergeometric canonical form of rational functions, how to perform Gosper’s algorithm, and how to find polynomial as well as hypergeometric solutions of recurrences, all in the mixed hypergeometric case. We have also indicated how to extend the concept of GFF to this case.

It remains to provide mixed hypergeometric generalizations of algorithms for finding rational solutions of recurrences [1, 11] and of algorithms for factorization of the corresponding operators [5]. The more efficient algorithm of van Hoeij for finding hypergeometric solutions [12] should also admit of a generalization to the mixed hypergeometric case.

Acknowledgements

The authors wish to thank Bill Gosper for providing the interesting tribasic and quadbasic examples just in the nick of time. They also thank an anonymous referee for very careful reading of a preliminary version as well as for providing additional references to the relevant literature.

A Algorithm MixedPoly

INPUT: p_{0}, …, p_{ρ}, g, h_{1}, …, h_{s} ∈ 𝔽[x, y], p_{0}, p_{ρ} ≠ 0 OUTPUT: general solution (f, λ) ∈ 𝔽[x, y] × 𝔽^{s} of Lf = g + ∑_{j=1}^{s} λ_{j}h_{j} where L = ∑_{i=0}^{ρ} p_{i} · E^{i}

CALLING SEQUENCE: MixedPoly(eqn, unknown, params)
EXTERNAL ALGORITHMS USED:

Poly(e, t, λ) returns general solution (t, λ) of the parametric nonhomogeneous ordinary recurrence e (see [2]) LinSolve(e, x, λ) returns general solution λ of the linear algebraic equations resulting from equating the coefficients of like powers of x on both sides of the polynomial equation e

α := max_{0≤i≤ρ} deg_{y} p_{i} p_{i,α}(x) := [y^{α}] p_{i} d := max_{0≤i≤ρ} deg_{x} p_{i,α}(x) p_{i,α,d} := [x^{d}] p_{i,α}(x) P (x) := ∑_{i=0}^{ρ} p_{i,α,d}x^{i} R := {σ ∈ ℕ_{0}^{m}; P (q^{σ}) = 0} rhs := g + ∑_{j=1}^{s} λ_{j}h_{j} f := 0

while R ≠ ∅ or rhs ≠ 0 do

if R ≠ ∅ then μ := max_{≺} R else μ := ⊥ if rhs ≠ 0 then β := deg_{y} rhs else β := ⊥

    if α + μ ≻ β then

φ := μ
(t′, λ′) := Poly(∑_{i=0}^{ρ} p_{i,α}(x)q^{iφ}t(x + i) = 0, t, λ)
f := f + t′y^{φ}

rhs := rhs|_{λ←λ′} − L (t′y^{φ})
R := R \ {φ}

    end
    else if α ⊆ β then

φ := β − α

(t′, λ′) := Poly(∑_{i=0}^{ρ} p_{i,α}(x)q^{iφ}t(x + i) = [y^{β}] rhs, t, λ)

f := f + t′y^{φ}
rhs := rhs|_{λ←λ′} − L (t′y^{φ})
if φ = μ then R := R \ {φ}

    end
    else

λ′ := LinSolve([y^{β}] rhs = 0, x, λ)
rhs := rhs|_{λ←λ′}

    end
end
return f.

NB: If either Poly or LinSolve fails then MixedPoly fails as well.

B Algorithm CanonicalForm

INPUT: r ∈ 𝔽(x, y) \ {0} OUTPUT: canonical form of r

EXTERNAL ALGORITHMS USED:

Resultant(a, b, x) returns the resultant of polynomials a, b w.r.t. x GCD(a, b) returns the mixed monic gcd of polynomials a, b

let r = (u/v) · (a_{0}/b_{0}) where a_{0}b_{0} ⊥ M, a_{0}u ⊥ b_{0}v, u, v ∈ M, and b_{0} is mixed monic P (ξ, η) := Resultant(a(x, y), b(x + ξ, ηy), x) for i = 1, …, m do

P (ξ, η) := P (ξ, η)·Resultant(a(x, y), b(x + ξ, ηy), y_{i})

end

let h_{1} < h_{2} < ⋯ < h_{N} be the roots h ∈ ℕ_{0} of p(h) = P (h, q^{h}) c_{0} := 1 for i = 1, …, N do

s_{i} := GCD(a_{i−1}, E^{h_{i}}b_{i−1})

a_{i} := a_{i−1}/s_{i}
b_{i} := b_{i−1}/E^{−h_{i}}s_{i}
c_{i} := c_{i−1} ∏_{j=1}^{h_{i}} E^{−j}s_{i}

end

a := u · a_{N} b := v · b_{N} c := c_{N} return (a, b, c).

C Algorithm MixedGosper

INPUT: mixed hypergeometric sequence t_{n}
OUTPUT: mixed hypergeometric sequence s_{n} such that s_{n+1} − s_{n} = t_{n}, if it exists

r := Φ(t_{n+1}/t_{n}) (see Theorem 2.2)

(a, b, c) := CanonicalForm(r)

u := 1

for i = 1 to m do

if π_{i}(ab) = 0

then

e_{i} := 0

else

(a_{i}, b_{i}, c_{i}) := CanonicalForm(π_{i}(b/a))

if ∃v, w ∈ M : (v ⊥ w and a_{i}/b_{i} = v(q)/w(q))

    then

e_{i} := deg_{y_{i}} w else e_{i} := 0

end

u := u · y_{i}^{e_{i}}

end

(p, ⟨ ⟩) := MixedPoly(a · Ep − u(q) · E^{−1}b · p = u(q)u · c, p, ⟨ ⟩) R := (p · E^{−1}b)/(u · c) return Φ^{−1}(R) · t_{n}.

NB: If MixedPoly fails then MixedGosper fails as well.

D Algorithm MixedHyper

INPUT: p_{0}, …, p_{ρ} ∈ 𝔽[x, y], p_{0}, p_{ρ} ≠ 0 OUTPUT: mixed hypergeometric solution f of ∑_{i=0}^{ρ} p_{i} · E^{i}f = 0, if it exists

for all mixed monic factors a of p_{0} and b of E^{1−ρ}p_{ρ} do

for i = 1 to ρ do P_{i} := p_{i} ∏_{j=0}^{i−1} E^{j}a ∏_{j=i}^{ρ−1} E^{j}b α := min_{0≤i≤ρ} mindeg_{y}P_{i} for i = 0 to ρ do p_{i,α} := [y^{α}] P_{i} d := max_{0≤i≤ρ} deg_{x} p_{i,α} for i = 0 to ρ do p_{i,α,d} := [x^{d}] p_{i,α}

for all τ such that ∑_{i=1}^{ρ} p_{i,α,d}τ^{i} = 0 do
    for j = 1 to m do
      order terms lexicographically with y_{j} first

α_{j} := min_{0≤i≤ρ} mindeg_{y}P_{i}

for i = 0 to ρ do p_{i,α_{j}} := [y^{α^{j}}] P_{i}

d_{j} := max_{0≤i≤ρ} deg_{x} p_{i,α_{j}}

for i = 0 to ρ do p_{i,α_{j} ,d_{j}} := [x^{d}] p_{i,α_{j}}

let σ_{1}, …, σ_{s} be the roots of ∑_{i=0}^{ρ} p_{i,α_{j} ,d_{j}} σ^{i} = 0

      let b_{j} be the minimum of the j-th components of σ_{1}, . . . , σ_{s}

u := y_{1}^{b_{1}} ⋯ y_{m}^{b_{m}}

(d̄, ⟨ ⟩) := MixedPoly(∑_{i=0}^{ρ} τ^{i}u(q)^{i}P_{i} · E^{i}d̄= 0, d̄, ⟨ ⟩)

    if d̄≠ 0 then

r := τu(q) \frac{a}{b} \frac{Ed̄}{d̄}

      return a non-zero solution f of Ef = rf and stop
    end
end
end
stop.

NB: If the original term order is lexicographic with y_{1} as the first variable, then in the innermost loop j runs from 2 to m and we set b_{1} = 0.

References

[1] S. A. Abramov, Rational solutions of linear difference and q-difference equations with polynomial coefficients, Programming and Comput. Software 21 (1995) 273–278.

[2] S. A. Abramov, M. Bronstein and M. Petkovšek, On polynomial solutions of linear operator equations, in: T. Levelt, ed., Proc. ISSAC ’95 (ACM Press, New York, 1995) 290–296.

[3] S. A. Abramov, P. Paule and M. Petkovšek, q-Hypergeometric solutions of q-difference equations, Discrete Math. 180 (1998) 3–22.

[4] G. Almkvist and D. Zeilberger, The method of differentiating under the integral sign, J. Symbolic Comp. 10 (1990) 571–591.

[5] M. Bronstein and M. Petkovšek, An introduction to pseudo-linear algebra, Theoret. Comput. Sci. 157 (1996) 3–33.

[6] R. M. Cohn, Difference Algebra, Interscience Publishers, New York, 1965.

[7] G. Gasper and M. Rahman, Basic Hypergeometric Series, vol. 35 of Encyclopedia of Mathematics and its Applications, Cambridge University Press, Cambridge, 1990.

[8] R. W. Gosper, Jr., Indefinite hypergeometric sums in MACSYMA, in: Proc. ’77 MACSYMA Users Conference, Berkeley CA, 1977, 237–252.

[9] R. W. Gosper, Jr., Decision procedure for indefinite hypergeometric summation, Proc. Natl. Acad. Sci. USA 75 (1978) 40–42.

[10] R. W. Gosper, Jr., private communication, 1998.

[11] M. van Hoeij, Rational solutions of linear difference equations, Proc. ISSAC ’98, to appear.

[12] M. van Hoeij, Finite singularities and hypergeometric solutions of linear recurrence equations, submitted to MEGA ’98.

[13] T. H. Koornwinder, On Zeilberger’s algorithm and its q-analogue, J. Comp. Appl. Math. 48 (1993) 91–111.

[14] B. Mishra, Algorithmic Algebra, Springer-Verlag, New York 1993.

[15] R. Moenck, On computing closed forms for summations, in: Proc. ’77 MACSYMA Users Conference, Berkeley CA, 1977, 225–236.

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

[17] P. Paule and A. Riese, A Mathematica q-analogue of Zeilberger’s algorithm based on an algebraically motivated approach to q-hypergeometric telescoping, Fields Inst. Comm. 14 (1997) 179–210.

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

[19] M. Petkovšek, Hypergeometric solutions of linear recurrences with polynomial coefficients, J. Symbolic Comp. 14 (1992) 243–264.

[20] A. Riese, A generalization of Gosper’s algorithm to bibasic hypergeometric summation, Electr. J. Comb. 3 (1996) #R19.

[21] U. B. Singh, Some summation formula for multibasic hypergeometric series, Glasnik Mat. 31(51) (1996) 263–268.

[22] R. P. Stanley, Enumerative Combinatorics I, Wadsworth & Brooks/Cole, Monterey CA, 1986.

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

[24] D. Zeilberger, q-analog of Gosper’s algorithm, unpublished manuscript, 1990.

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

[26] D. Zeilberger, The method of creative telescoping, J. Symbolic Comp. 11 (1991) 195–204.