Gosper’s Algorithm, Accurate Summation, and the Discrete Newton-Leibniz Formula^{∗}

S. A. Abramov^{†}

Russian Academy of Sciences Dorodnicyn Computing Centre

Vavilova 40, 119991, Moscow GSP-1, Russia

sabramov@ccas.ru

M. Petkovšek^{‡}

Department of Mathematics University of Ljubljana

Jadranska 19, SI-1000 Ljubljana, Slovenia

     marko.petkovsek@fmf.uni-lj.si

ABSTRACT

Sufficient conditions are given for validity of the discrete Newton-Leibniz formula when the indefinite sum is obtained either by Gosper’s algorithm or by Accurate Summation algorithm. It is shown that sometimes a polynomial can be factored from the summand in such a way that the safe summation range is increased.

1. INTRODUCTION

Let K be a field of characteristic zero. A function t : I ↦ K defined on an interval of integers I ⊆ ℤ is a

• hypergeometric term if there are nonzero polynomials a_{0}(n), a_{1}(n) ∈ K[n] such that a_{1}(n)t(n + 1) + a_{0}(n)t(n) = 0 for all n such that n, n + 1 ∈ I;

• P -recursive sequence if there are polynomials a_{0}(n), a_{1}(n), …, a_{ρ}(n) ∈ K[n] such that a_{0}(n), a_{ρ}(n) are nonzero and a_{ρ}(n)t(n + ρ) + ⋯ + a_{0}(n)t(n) = 0 for all n such that n, …, n + ρ ∈ I.

Each hypergeometric term is, of course, a P -recursive sequence.

If t(n) is a hypergeometric term, one can use the wellknown Gosper’s algorithm [6] to find (if it exists) another hypergeometric term u(n) which satisfies the key equation

u(n + 1) − u(n) = t(n)  (1)

for all n ∈ I \ S where S is a finite set. Summing this equation on n from v to w we get the discrete analog of the Newton-Leibniz formula

∑_{n=v}^{w} t(n) = u(w + 1) − u(v)  (2)

∗The work is partially supported by the ECO-NET program of the French Foreign Affairs Ministry. †Partially supported by RFBR under grant 04-01-00757. ‡Partially supported by ARRS under grant P1-0294.

Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for pro�t or commercial advantage and that copies bear this notice and the full citation on the �rst page. To copy otherwise, to republish, to post on servers or to redistribute to lists, requires prior speci�c permission and/or a fee. ISSAC’05, July 24�27, 2005, Beijing, China. Copyright 2005 ACM 1-59593-095-705/0007 ...$5.00.

provided that [v, w] ∩ ℤ ⊆ I \ S.

In many existing implementations of Gosper’s algorithm, however, indiscriminate use of (2) sometimes results in wrong answers. Here is a case in point.

Example 1. Consider the sequence

t(n) = \frac{\binom{2n−3}{n}}{4^{n}} ,  (3)

which is defined for all n ∈ ℤ. This is a hypergeometric term which satisfies

2(n + 1)(n − 2)t(n + 1) = (2n − 1)(n − 1)t(n) (4)

for all n ∈ ℤ. Gosper’s algorithm succeeds with input t(n) and returns

u(n) = \frac{2n(n + 1)\binom{2n−3}{n}}{(n − 2)4^{n}} .

Summing equation (1) on n from 0 to m the left-hand side telescopes, and we obtain

∑_{n=0}^{m} t(n) =(?) u(m + 1) − u(0) = \frac{(m + 1)(m + 2)\binom{2m−1}{m+1}}{2(m − 1)4^{m}} .

(5)

But the expression on the right gives the true value of the sum only at m = 0. At m = 1 it is undefined, while at each m ≥ 2 its value is 3/8 less than the actual value of the sum. The problem here is that u(n) is undefined at n = 2, hence equation (1) does not hold for n ∈ {1, 2}, and summing it over a range including 1 or 2 may give a wrong answer.

This is not an isolated example: a similar phenomenon seems to occur with the sum

∑_{n=0}^{m} \frac{\binom{2n−p}{n}}{4^{n}}

for each positive integer p.

If t is a P -recursive sequence, then one can use Accurate Summation algorithm from [3], or its generalization in [5], to solve equation (1) (we discuss this algorithm in Section 5). Problems similar to those arising in Example 1 are possible when one uses the resulting Newton-Leibniz formula. Notice that one can apply Accurate Summation algorithm in the case ρ = 1 as an alternative to Gosper’s algorithm; then the incorrect formula (5) will appear again.

This common error is the discrete analogon of a wellknown error in definite integration committed by some of the early symbolic integrators: when attempting to evaluate I = ∫_{a}^{b} f(x)dx by computing first an antiderivative F (x) such that F′(x) = f(x), and then using the Newton-Leibniz formula I = F (b)−F (a), we may obtain an incorrect answer unless F (x) is continuous on [a, b]. For example, the actual value of

∫_{−1}^{1} \frac{x^{2} + 1}{x^{4} + x^{2} − 1} dx

is π, but using the antiderivative arctan(x − 1/x) in the Newton-Leibniz formula gives 0.

The obvious solution is to split the summation interval into several subintervals that do not contain the exceptional points from S. In this paper we analyze the exceptional set S that appears in Gosper’s algorithm when summing hypergeometric terms, and more generally, in the Accurate Summation algorithm [3] when summing P -recursive sequences.

Section 3 provides sufficient conditions for the Newton-Lebniz formula (2) to hold when the indefinite sum u(n) is obtained by Gosper’s algorithm, and Section 5 does the same for Accurate Summation. These conditions provide a bounding interval for the exceptional set S, and are of two kinds: a priori, which are weaker but readily available even before running the algorithms, as they are based on the singularities of the operator annihilating the summand; and a posteriori, which are stronger but available only after running the algorithms, as they are based on their output. On the other hand, in Section 4 we prove that for proper hypergeometric terms the discrete Newton-Leibniz formula is valid without restrictions. For general P -recursive sequences Section 6 shows that sometimes a polynomial can be factored from the summand in such a way that the size of the bounding interval in the a priori condition is decreased.

A thorough analysis of the relationship between hypergeometric terms as syntactic objects and their analytic meaning in the context of summation has been provided by M. Schorn in [8]. The solution proposed there for evaluation of sums such as the one in Example 1 is by means of suitably chosen limiting processes.

2. PRELIMINARIES

Definition 1. Following conventional notation, the rising factorial power (α)_{n} and its reciprocal 1/(β)_{n} are defined for α, β ∈ K and n ∈ ℤ by

(α)_{n} = \{ \begin{matrix}∏_{k=0}^{n−1}(α + k), & n ≥ 0, \\ ∏_{k=1}^{|n|} \frac{1}{α − k} , & n < 0, α ≠ 1, 2, …, |n|, \\ undefined, & otherwise; \end{matrix}

\frac{1}{(β)_{n}} = \{ \begin{matrix}∏_{k=0}^{n−1} \frac{1}{β + k} , & n ≥ 0, β ≠ 0, −1, …, 1 − n, \\ ∏_{k=1}^{|n|}(β − k), & n < 0, \\ undefined, & otherwise. \end{matrix}

Note that if (α)_{n} resp. 1/(β)_{n} is defined for some n ∈ ℤ, then (α)_{n+1} resp. 1/(β)_{n−1} is defined for that n as well.

More precisely, if α ∈ ℤ and α ≥ 1 then (α)_{n} is defined on [−α + 1, ∞) ∩ ℤ, otherwise it is defined on all ℤ. Similarly, if β ∈ ℤ and β ≤ 0 then 1/(β)_{n} is defined on (−∞, −β] ∩ ℤ, otherwise it is defined on all ℤ. Thus (α)_{n} and 1/(β)_{n} are hypergeometric terms which satisfy

(α)_{n+1} = (α + n)(α)_{n}, (β + n)/(β)_{n+1} = 1/(β)_{n} (6)

whenever (α)_{n} and 1/(β)_{n+1} are defined.

If I ⊆ ℤ is an infinite interval of integers we denote

I^{+} = \{ \begin{matrix}(−∞, a + 1] ∩ ℤ, & if I = (−∞, a] ∩ ℤ, \\ I, & otherwise; \end{matrix}

I^{−} = \{ \begin{matrix}(−∞, a − 1] ∩ ℤ, & if I = (−∞, a] ∩ ℤ, \\ I, & otherwise. \end{matrix}

We use E to denote the shift operator w.r.t. n, so that E t(n) = t(n + 1). Since juxtaposition can mean not only operator application but also composition of operators, we use ∘ to denote the latter in case of ambiguity, so that, e.g., E ∘ t(n) = t(n + 1) ∘ E = t(n + 1)E. Sometimes we use parentheses to denote operator application, writing, e.g., E(1) = E 1 = 1.

Definition 2. For a linear difference operator

L = a_{ρ}(n)E^{ρ} + a_{ρ−1}(n)E^{ρ−1} + ⋯ + a_{0}(n)  (7)

where ρ ≥ 1, a_{ρ}(n), …, a_{0}(n) ∈ K[n], a_{ρ}(n)a_{0}(n) ≢ 0 and gcd(a_{0}(n), …, a_{ρ}(n)) = 1, we define the sets S_{Ll} of leading and S_{Lt} of trailing integer singularities by

S_{Ll} = {x ∈ ℤ; a_{ρ}(x − ρ) = 0},
S_{Lt} = {x ∈ ℤ; a_{0}(x) = 0}.

We call

• m_{Ll} = min(S_{Ll} ∪ {+∞}) the minimal leading singularity of L,

• M_{Ll} = max(S_{Ll} ∪ {−∞}) the maximal leading singularity of L,

• m_{Lt} = min(S_{Lt} ∪ {+∞}) the minimal trailing singularity of L,

• M_{Lt} = max(S_{Lt} ∪ {−∞}) the maximal trailing singularity of L.

Proposition 1. Let L be as in (7) and b(n) ∈ K[n]. If a rational function y(n) ∈ K(n) satisfies

a_{ρ}(n)y(n + ρ) + ⋯ + a_{0}(n)y(n) = b(n),  (8)

then y(n) has no integer poles outside the interval (possibly empty) [m_{Ll}, M_{Lt}].

For a proof, see [1].

3. WHEN CAN GOSPER’S ALGORITHM BE USED TO SUM HYPERGEOMETRIC TERMS?

We denote Gosper’s algorithm hereafter by GA. Consider the case when (8) has the form

a_{1}(n)t(n + 1) + a_{0}(n)t(n) = 0  (9)

and set L = a_{1}(n)E + a_{0}(n). Let a hypergeometric term t(n) satisfy equation (9). Given a_{0}(n), a_{1}(n) as input, GA tries to construct r(n) ∈ K(n) such that

a_{0}(n)r(n + 1) + a_{1}(n)r(n) = −a_{1}(n)  (10)

(this can also be done by the algorithms from [1] or [2]). If such r(n) exists then u(n) = r(n)t(n) satisfies the key equation (1), possibly with finitely many exceptions. We now give two kinds of sufficient conditions for this u(n) to satisfy equation (1) and for the discrete Newton-Leibniz formula in the form

∑_{k=v}^{w} t(k) = u(w) − u(v) + t(w)  (11)

to be valid:

1. an a posteriori condition, depending on the poles of r(n) (Proposition 2),

2. an a priori condition, depending only on the integer singularities of L (Theorem 1).

In both, we make the following assumptions:

• L = a_{1}(n)E + a_{0}(n) is an operator of type (7) with ρ = 1,

• r(n) ∈ K(n) is a rational function which satisfies (10) as an equation in K(n),

• v, w are integers such that v ≤ w,

• I_{1} := [v, w − 1] ∩ ℤ,

• t(n) is a K-valued sequence which is defined for all n ∈ [v, w] ∩ ℤ and satisfies (9) for all n ∈ I_{1},

• u(n) is a K-valued sequence such that u(n) = r(n)t(n) whenever both r(n) and t(n) are defined.

Remark 1. Since u(n) = r(n)t(n) it is clear that, in general, formula (11) should be used instead of (2), because the latter formula needs values of the summand lying outside the summation interval which however may be undefined. A nice example is provided by the sum

∑_{k=0}^{2n}(−1)^{k}\binom{4n}{2k}/\binom{2n}{k}

whose evaluation was posed as Problem 10494 in the Amer. Math. Monthly in 1996. Here GA succeeds (see [7]), but the summand is undefined everywhere outside the summation interval.

Proposition 2. (a posteriori condition for GA) If r(n) has no integer poles in [v, w], then the key equation (1) holds for all n ∈ I_{1}, and the discrete Newton-Leibniz formula (11) is valid.

Proof: By assumption, t(n), t(n + 1), r(n), r(n + 1), u(n), u(n + 1) are defined for all n ∈ I_{1}, and (9), (10) are valid on I_{1}. Therefore, for all n ∈ I_{1},

a_{0}(n)u(n + 1) = a_{0}(n)r(n + 1)t(n + 1)
= −a_{1}(n)(1 + r(n))t(n + 1) (by (10))
= a_{0}(n)(1 + r(n))t(n) (by (9))
= a_{0}(n)u(n) + a_{0}(n)t(n),

or equivalently,

a_{0}(n)(u(n + 1) − u(n)) = a_{0}(n)t(n).  (12)

Pick an n ∈ I_{1}. If a_{0}(n) ≠ 0 then (12) implies (1). If a_{0}(n) = 0 then, since gcd(a_{0}, a_{1}) = 1, we have a_{1}(n) ≠ 0. Hence (9) implies t(n+1) = 0 and (10) implies r(n)+1 = 0. Therefore u(n + 1) = r(n + 1)t(n + 1) = 0 and u(n) + t(n) = (r(n) + 1)t(n) = 0, so (1) holds for all n ∈ I_{1}. Summing (1) over I_{1} yields (11).

Theorem 1. (a priori condition for GA) If [v, w]∩[m_{Lt} + 1, M_{Ll}−1] = ∅, then the key equation (1) holds for all n ∈ I_{1}, and the discrete Newton-Leibniz formula (11) is valid.

Proof: Since r(n) satisfies (10), Proposition 1 implies that r(n) has no integer poles outside the interval [α, β] where

α = min({x ∈ ℤ; a_{0}(x − 1) = 0} ∪ {+∞})
= m_{Lt} + 1,
β = max({x ∈ ℤ; a_{1}(x) = 0} ∪ {−∞})
= M_{Ll} − 1.

By assumption, the interval [v, w] is disjoint from [α, β]. Hence r(n) has no poles in [v, w], and the assertion follows from Proposition 2.

In practice, one would run GA and then check if the a posteriori condition of Proposition 2 is satisfied, i.e., if r(n) has any integer poles in the summation interval. If yes, this interval would be split into several subintervals in order to guarantee correct evaluation. But it may be useful to check the a priori condition of Theorem 1 first, because this will, in general, restrict the relevant domain to check for poles of r(n).

Example 2. For the hypergeometric term t(n) = \binom{2n−3}{n}/4^{n} of Example 1, we have a_{0}(n) = −(2n − 1)(n − 1), a_{1}(n) = 2(n + 1)(n − 2), r(n) = 2n(n + 1)/(n − 2), and u(n) = 2n(n + 1)\binom{2n−3}{n}/((n − 2)4^{n}). Thus S_{Lt} = \{1\}, S_{Ll} = {0, 3}, m_{Lt} = 1, M_{Ll} = 3, [m_{Lt} + 1, M_{Ll} − 1] = [2, 2], and the only integer pole of r(n) is n = 2. In this case both the a priori and the a posteriori conditions give the same set {2} to be avoided by the summation interval. As predicted by either condition, the key equation (1) fails at n = 1 and n = 2 because u(n) or u(n + 1) are undefined there. One is tempted to absorb the denominator factor n − 2 into the binomial coefficient, and replace u(n) by the sequence ū(n) = n(n + 1)\binom{2n−1}{n}/((2n − 1)4^{n}) which is defined everywhere, and agrees with u(n) for all n ≠ 1, 2. But then equation

ū(n + 1) − ū(n) = t(n)  (13)

fails at n = 0 and n = 1.

Example 3. Let

t(n) = \{ \begin{matrix}(n − 2)(n − 3)(n − 5)(n − 1)!, & n ≥ 2 \\ (n − 2)(n − 3)(n − 5) \frac{(−1)^{n}}{(−n)!} , & n ≤ 1 \end{matrix}

where we define as usual 1/k! = 0 when k is a negative integer. This is a hypergeometric term which satisfies

(n − 5)(n − 3) t(n + 1) = (n − 4)(n − 1)n t(n)

for all n ∈ ℤ. Here we have a_{0}(n) = −(n − 4)(n − 1)n, a_{1}(n) = (n − 5)(n − 3), r(n) = (n − 6)/((n − 2)(n − 3)), and u(n) = r(n)t(n). Thus S_{Lt} = {0, 1, 4}, S_{Ll} = {4, 6}, m_{Lt} = 0, M_{Ll} = 6, [m_{Lt} + 1, M_{Ll} − 1] = [1, 5], and the set of integer poles of r(n) is {2, 3}. The set to be avoided by the summation interval given by the a priori condition is {1, 2, 3, 4, 5}, while the analogous set given by the a posteriori condition is {2, 3}. The key equation (1) fails at n = 1, 2, 3 as predicted by the a posteriori condition, because u(n) or u(n + 1) are undefined there. One can try cancelling the factor (n − 2)(n − 3) and replace u(n) by the sequence

ū(n) = \{ \begin{matrix}(n − 5)(n − 6)(n − 1)!, & n ≥ 2 \\ (n − 5)(n − 6) \frac{(−1)^{n}}{(−n)!} , & n ≤ 1 \end{matrix}

which is defined everywhere, and agrees with u(n) for all n ≠ 2, 3. But equation (13) still fails at n = 1.

4. SUMMATION OF PROPER HYPERGE-OMETRIC TERMS

It is clear that the a priori condition given in Theorem 1 is, in general, too cautious: e.g., if the summand is a polynomial sequence then the integer singularities of the corresponding recurrence present no obstacles to validity of the discrete Newton-Leibniz formula (11). The following example shows that even the a posteriori condition given in Proposition 2 can sometimes be too pessimistic.

Example 4. Let t(n) = (n − 2)(−1/2)_{n}/(4n!). This hypergeometric term is defined for all n ∈ ℤ (note that t(n) = 0 for n < 0) and satisfies Lt(n) = 0 for all n ∈ ℤ where L is the same operator as in Example 2. Thus both Theorem 1 and Proposition 2 require the point n = 2 to be excluded from the summation interval. Equation (1) indeed fails at n = 1 and n = 2 because u(n) = r(n)t(n) is undefined at n = 2. But if we cancel the factor n − 2 in the product r(n)t(n) and replace u(n) by the sequence

ū(n) = n(n + 1) \frac{(−1/2)_{n}}{2n!} ,

then equation (13) holds for all n ∈ ℤ, and the discrete Newton-Leibniz formula

∑_{n=v}^{w} t(n) = ū(w + 1) − ū(v)  (14)

is valid for all v ≤ w.

This example also shows that, thanks to possible singularities, a hypergeometric term (or a P -recursive sequence) is, in general, not uniquely defined by its annihilating operator and the appropriate number of initial values. In fact, it is shown in [4] that every natural number is the dimension of the kernel of some operator of type (7) with ρ = 1 in the space of sequences t : ℤ → K.

The hypergeometric term t(n) from Example 4 is an instance of a proper term which we are going to define now. Then we show in Theorem 2 that there are no restrictions on the validity of the discrete Newton-Leibniz formula for proper terms.

Definition 3. A hypergeometric term t(n) defined on an infinite interval I of integers is proper if there are

• a constant z ∈ K,

• a polynomial p(n) ∈ K[n], • nonnegative integers q, r,

• constants α_{1}, …, α_{q}, β_{1}, …, β_{r} ∈ K

such that

t(n) = p(n)z^{n} \frac{∏_{i=1}^{q}(α_{i})_{n}}{∏_{j=1}^{r}(β_{j})_{n}}  (15)

for all n ∈ I.

Theorem 2. Let t(n) be a proper hypergeometric term defined on an infinite interval I of integers and given by (15). Denote a(n) = z ∏_{i=1}^{q}(n + α_{i}) and b(n) = ∏_{j=1}^{r}(n + β_{j}). If a polynomial y(n) ∈ K[n] satisfies

a(n)y(n + 1) − b(n − 1)y(n) = p(n)  (16)

and if

ū(n) = y(n)z^{n} \frac{∏_{i=1}^{q}(α_{i})_{n}}{∏_{j=1}^{r}(β_{j})_{n−1}}

for all n ∈ I^{+} (see Section 2 for notation), then equation (13) holds for all n ∈ I, and the discrete Newton-Leibniz formula (14) is valid whenever [v, w] ∩ ℤ ⊆ I.

Proof: By assumption, (16) holds for all n ∈ I. Multiplying it by z^{n} ∏_{i=1}^{q}(α_{i})_{n}/ ∏_{j=1}^{r}(β_{j})_{n} yields

z^{n} \frac{∏_{i=1}^{q}(α_{i})_{n}}{∏_{j=1}^{r}(β_{j})_{n}} a(n)y(n + 1)

− z^{n} \frac{∏_{i=1}^{q}(α_{i})_{n}}{∏_{j=1}^{r}(β_{j})_{n}} b(n − 1)y(n)
= t(n).  (17)

Since (α_{i})_{n} and 1/(β_{j})_{n} are defined for all n ∈ I, (α_{i})_{n+1} and 1/(β_{j})_{n−1} are defined there too. By (6), a(n) ∏_{i=1}^{q}(α_{i})_{n} = z ∏_{i=1}^{q}(α_{i})_{n+1} and b(n−1)/ ∏_{j=1}^{r}(β_{j})_{n} = 1/ ∏_{j=1}^{r}(β_{j})_{n−1}. Hence (17) is the same as (13), and (14) follows by summing (13) over [v, w] ∩ ℤ.

Example 5. Even though the hypergeometric term (3) from Example 1 defined on I = ℤ can be written in terms of rising factorials as

t(n) = \frac{(n − 2)_{n}}{4^{n}(1)_{n}} ,

one can show that it is not a proper term on ℤ. However, it can also be written as

t(n) = \{ \begin{matrix}2 t^{∗}(n), & n < 2, \\ t^{∗}(n), & n ≥ 2, \end{matrix}

where

t^{∗}(n) = (2 − n) \frac{(−1/2)_{n}}{4(1)_{n}}

is a proper term (namely the one discussed in Example 4). So to evaluate ∑_{n=v}^{w} t(n) one can first split the summation range at n = 2, then use Theorem 2 on both subranges.

5. WHEN CAN ACCURATE SUMMATION BE USED TO SUM P-RECURSIVE SE-QUENCES?

By Accurate Summation algorithm (hereafter denoted by AS) we mean a specific version of the general Accurate Integration algorithm given in [3] for integration/summation of solutions of Ore equations. This version, which is adapted for sequences that satisfy equations of the form (8) with b(n) = 0, solves the following problem: Let a minimal annihilator L of the form (7) be known for a K-valued sequence t(n). Determine if there exists a sequence u(n) which satisfies (1), and has a minimal annihilator L̃ of order ρ.

It is shown in [3] that if such a u exists then it can be expressed as R t where R is an operator of order ρ − 1 with rational-function coefficients. AS constructs R if it exists. (GA solves this problem when ρ = 1.) In order to analyze the validity of the discrete Newton-Leibniz formula (11) in this case, we need to express explicitly the quotient and the remainder of a linear difference operator when divided by the first-order operator E − 1 from the left. The notion of the adjoint difference operator is useful here.

Definition 4. Let

L = ∑_{k=0}^{ρ} b_{k}(n)E^{k}

be an operator in K(n)[E]. Its adjoint L^{∗} ∈ K(n)[E^{−1}] is defined as

L^{∗} = ∑_{k=0}^{ρ} E^{−k} ∘ b_{k}(n) = ∑_{k=0}^{ρ} b_{k}(n − k)E^{−k}.

It is straightforward to verify that (L_{1} ∘ L_{2})^{∗} = L_{2}^{∗} ∘ L_{1}^{∗}.

Lemma 1. Let P = ∑_{k=0}^{ρ} b_{k}(n)E^{k}, R = ∑_{k=0}^{ρ−1} c_{k}(n)E^{k} be operators from K(n)[E], and p(n) ∈ K(n) a rational function such that P = (E − 1) ∘ R + p. Then c_{k}(n) = ∑_{j=1}^{ρ−k} b_{k+j}(n − j) and p = P^{∗}(1).

Proof: Let

Q = ∑_{k=0}^{ρ−1} ∑_{j=1}^{ρ−k} b_{k+j}(n − j)E^{k}

and q = P^{∗}(1). Then

(E − 1) ∘ Q + q
= ∑_{k=0}^{ρ−1} ∑_{j=1}^{ρ−k} b_{k+j}(n − j + 1)E^{k+1} − Q + P^{∗}(1)

= ∑_{k=1}^{ρ} ∑_{j=1}^{ρ−k+1} b_{k+j−1}(n − j + 1)E^{k}

− ∑_{k=0}^{ρ−1} ∑_{j=2}^{ρ−k+1} b_{k+j−1}(n − j + 1)E^{k} + ∑_{j=0}^{ρ} b_{j}(n − j)

= ∑_{k=0}^{ρ} ∑_{j=1}^{ρ−k+1} b_{k+j−1}(n − j + 1)E^{k}

− ∑_{k=0}^{ρ−1} ∑_{j=2}^{ρ−k+1} b_{k+j−1}(n − j + 1)E^{k}

= b_{ρ}(n)E^{ρ} + ∑_{k=0}^{ρ−1} b_{k}(n)E^{k} = P.

As the quotient and remainder in operator division are unique, it follows that R = Q and p = q.

Note that just to find the remainder p = P^{∗}(1), it suffices to take adjoints on both sides of equation P = (E −1)∘R+p which results in P^{∗} = R^{∗} ∘(E −1)^{∗} +p = R^{∗} ∘(E^{−1} −1)+p, and apply this to 1.

Remark 2. Let r(n) ∈ K(n) be a rational function, and L a difference operator as in (7). By Lemma 1, the remainder of 1 − rL when divided by E − 1 from the left is equal to

(1 − rL)^{∗}(1) = (1 − L^{∗} ∘ r)(1) = 1 − (L^{∗} ∘ r)(1) = 1 − L^{∗}r.

Hence there exists an operator R such that 1 − rL = (E − 1) ∘ R, if and only if L^{∗}r = 1. This observation forms the basis of Accurate Summation.

Algorithm AS

Input: L = ∑_{k=0}^{ρ} a_{k}(n)E^{k} ∈ K[n, E]. Output: r(n) ∈ K(n) and R ∈ K(n)[E]

such that 1 − rL = (E − 1) ∘ R, if they exist.

if there exists r(n) ∈ K(n) such that L^{∗}r = 1 then

for k := 0, 1, …, ρ − 1 do

c_{k}(n) := − ∑_{j=1}^{ρ−k} r(n − j)a_{k+j}(n − j);

R := ∑_{k=0}^{ρ−1} c_{k}(n)E^{k}; return (r(n), R)

else such r(n) and R do not exist.

We can find a rational-function solution r(n) of L^{∗}r = 1 using, e.g., the algorithm from [1] or the algorithm from [2]. A generalization of [3] was given in [5]; however, the approach taken in [3] has the advantage of simplicity, as it only uses the adjoint operator and algorithms for finding rational solutions. This simplifies the investigation of solutions that are obtained by AS, and enables us to formulate a priori conditions for AS, similar to Theorem 1 (see Theorems 3 and 5 below).

Assume that AS succeeds with L, returning r and R. It is shown in [3] that

• if L is a minimal annihilator for t, then a minimal annihilator for u = R t is L̃ = 1 − R ∘ (E − 1) (note that L̃ has the same order ρ as L);

• the sequence u(n) = R t(n) satisfies (1), possibly with finitely many exceptions, for any sequence t such that Lt = 0 (not only for those t’s whose minimal annihilator is L).

We now give two sufficient conditions for this u(n) to satisfy equation (1) and for the discrete Newton-Leibniz formula (11) to be valid:

1. an a posteriori condition, depending on the poles of r(n) and of the coefficients of R (Proposition 3),

2. an a priori condition, depending only on the integer singularities of L (Theorem 3).

In either case, we make the following assumptions:

• L ∈ K[n, E] is an operator of type (7),

• r(n) ∈ K(n) is a rational function which satisfies L^{∗}r = 1 as an equation in K(n),

• R ∈ K(n)[E] is an operator of order ρ − 1 which satisfies 1 − rL = (E − 1) ∘ R in K(n)[E],

• v, w are integers such that v ≤ w − ρ,

• I_{ρ} := [v, w − ρ] ∩ ℤ,

• t(n) is a K-valued sequence which is defined for all n ∈ [v, w] ∩ ℤ and satisfies L t(n) = 0 for all n ∈ I_{ρ},

• u(n) is a K-valued sequence such that u(n) = R t(n) whenever R t(n) is defined.

Proposition 3. (a posteriori condition for AS) If r(n) has no poles in I_{ρ} and the coefficients of R have no integer poles in [v, w − ρ + 1], then equation (1) holds for all n ∈ I_{ρ}, and the discrete Newton-Leibniz formula

∑_{k=v}^{w} t(k) = u(w − ρ + 1) − u(v) + ∑_{k=1}^{ρ} t(w − ρ + k) (18)

is valid.

Proof: By assumptions on t and R, u(n) and u(n + 1) are defined for all n ∈ I_{ρ}. As r(n) has no poles in I_{ρ},

u(n + 1) − u(n) = (E ∘ R) t(n) − R t(n)
= ((E − 1) ∘ R) t(n) = (1 − rL) t(n)
= t(n) − r(n) L t(n) = t(n)

for every n ∈ I_{ρ}. Thus (1) holds for all n ∈ I_{ρ}, and summing it over I_{ρ} yields (18).

Lemma 2. Let α, β ∈ ℤ ∪ {−∞, ∞}. If r(n) has no integer poles in [α, β] then the coefficients of R have

(i) no integer poles in [α + ρ, β + 1], and also

(ii) no integer poles in [α, β − ρ + 1].

Proof: Write R = ∑_{k=0}^{ρ−1} c_{k}(n)E^{k}. By Lemma 1,

c_{k}(n) = − ∑_{j=1}^{ρ−k} r(n − j)a_{k+j}(n − j)

for 0 ≤ k ≤ ρ − 1.

(i) By assumption, r(n − j) has no integer poles in [α + j, β + j], hence c_{k}(n) has no integer poles in ⋂_{1≤j≤ρ−k}[α + j, β + j] = [max_{1≤j≤ρ−k}(α + j), min_{1≤j≤ρ−k}(β+j)] = [α+ρ−k, β+1], and the coefficients of R have no integer poles in ⋂_{0≤k≤ρ−1}[α+ρ− k, β+1] = [max_{0≤k≤ρ−1}(α+ρ−k), β+1] = [α+ρ, β+1].

(ii) To prove the second assertion, we need to express the coefficients c_{k}(n) in a different way. Since 1 − rL = (E − 1) ∘ R, it follows from Remark 2 that L^{∗}r = ∑_{j=0}^{ρ} r(n − j)a_{j}(n − j) = 1. Shifting this k times we find that ∑_{j=0}^{ρ} r(n+k−j)a_{j}(n+k−j) = ∑_{j=−k}^{ρ−k} r(n− j)a_{k+j}(n − j) = 1. Therefore

c_{k}(n) = − ∑_{j=1}^{ρ−k} r(n − j)a_{k+j}(n − j)

+ ∑_{j=−k}^{ρ−k} r(n − j)a_{k+j}(n − j) − 1

= ∑_{j=0}^{k} r(n + j)a_{k−j}(n + j) − 1

for 0 ≤ k ≤ ρ − 1. By assumption, r(n + j) has no integer poles in [α − j, β − j], hence c_{k}(n) has no integer poles in ⋂_{0≤j≤k}[α − j, β − j] = [max_{0≤j≤k}(α − j), min_{0≤j≤k}(β − j)] = [α, β − k], and the coefficients of R have no integer poles in ⋂_{0≤k≤ρ−1}[α, β − k] = [α, min_{0≤k≤ρ−1}(β − k)] = [α, β − ρ + 1].

Theorem 3. (a priori condition for AS) If [v, w − ρ] ∩ [m_{Lt}, M_{Ll} − ρ] = ∅, then equation (1) holds for all n ∈ I_{ρ}, and the discrete Newton-Leibniz formula (18) is valid.

Proof: Rewrite L^{∗}r = 1 in the equivalent form L′r = 1 where L′ = E^{ρ} ∘ L^{∗} = ∑_{k=0}^{ρ} a_{ρ−k}(n + k)E^{k} ∈ K[n, E]. By Lemma 1, r(n) has no integer poles outside [m_{L′l}, M_{L′t}]. But S_{L′l} = S_{Lt} and S_{L′t} = S_{Ll} − ρ, therefore m_{L′l} = m_{Lt} and M_{L′t} = M_{Ll} −ρ, hence r(n) has no integer poles outside [m_{Lt}, M_{Ll} − ρ].

If m_{Lt} ≤ M_{Ll} − ρ, then both intervals [v, w − ρ] and [m_{Lt}, M_{Ll} − ρ] are nonempty, hence either w − ρ < m_{Lt} or M_{Ll} − ρ < v. In the former case, r(n) has no integer poles in (−∞, w − ρ], so by Lemma 2(i), the coefficients of R have no integer poles in (−∞, w−ρ+1]. In the latter case, r(n) has no integer poles in [v, ∞), so by Lemma 2(ii), the coefficients of R have no integer poles in [v, ∞). In either case, the result follows from Proposition 3.

If m_{Lt} > M_{Ll} −ρ then r(n) has no integer poles at all. By Lemma 2, the coefficients of R also have no integer poles, and the result again follows from Proposition 3.

A similar remark to the one stated immediately after the proof of Theorem 1 about the use of the a priori and a posteriori conditions in practice applies here as well.

Example 6. Let L = (n−3)(n−2)(n+1)E^{2}−(n−3)(n^{2}− 2n−1)E−(n−2)^{2}. Define t(n) by the initial values t(2) = a, t(3) = 0, t(4) = b, t(5) = c where a, b, c are arbitrary fixed complex numbers, and by the recurrence L t(n) = 0 when n ≤ 1 or n ≥ 6. Then it can be checked that L t(n) = 0 for all n ∈ ℤ. Algorithm AS succeeds with input L and returns r(n) = −1/((n−2)(n−3)), R = nE +1/(n−3). In this case S_{Lt} = {2}, m_{Lt} = 2, S_{Ll} = {1, 4, 5}, M_{Ll} − ρ = 5 − 2 = 3. So both the a posteriori and the a priori conditions reduce to 3 ∉ [v, w − 1]. This is the best possible, as the sequence u(n) = R t(n) is undefined at n = 3, and equation (1) does not hold for n = 2, 3. It can be verified that except in the special case b+4c = 0, there is no way to define u(3) so that (1) would hold for all n ∈ ℤ.

Remark 3. When ρ = 1 and v ≤ w − 1, Theorem 3 implies Theorem 1 in the following way. If r(n) satisfies equation (10) then it is easy to verify that r̄(n) := −r(n + 1)/a_{1}(n) satisfies L^{∗}r̄ = 1, and 1 − r̄L = (E − 1) ∘ R where R = r(n) is an operator of order 0. Thus u(n) = r(n)t(n) of Theorem 1 agrees with u(n) = R t(n) of Theorem 3. By the assumption of Theorem 1, [v, w] ∩ [m_{Lt} + 1, M_{Ll} − 1] = ∅.

If m_{Lt} + 1 ≤ M_{Ll} − 1 then either w ≤ m_{Lt} or M_{Ll} − 1 ≤ v − 1, hence [v, w − 1] ∩ [m_{Lt}, M_{Ll} − 1] = ∅, and the conclusion follows by Theorem 3. If m_{Lt} + 1 > M_{Ll} − 1 then m_{Lt} ≥ M_{Ll} − 1. Again we distinguish two cases: If m_{Lt} > M_{Ll} − 1, the conclusion follows by Theorem 3. If m_{Lt} = M_{Ll} − 1, then a_{0}(n) and a_{1}(n) have a common zero since a_{0}(m_{Lt}) = a_{1}(M_{Ll} − 1) = 0. But this contradicts the assumption of relative primality of the coefficients of L.

Note that when ρ ≥ 2, polynomials a_{0}(n) and a_{ρ}(n) need not be relatively prime.

6. TAKING ADVANTAGE OF POLYNO-MIAL FACTORS

In this section we show how polynomial factors of hypergeometric terms (even non-proper ones, such as the one in Example 3) and of P -recursive sequences can be used to strengthen the statements of Theorems 1 and 3 (i.e., to weaken the a priori conditions for validity of the discrete Newton-Leibniz formula).

Theorem 4. Assume that

• L = a_{1}(n)E + a_{0}(n) and L̄ = ā_{1}(n)E + ā_{0}(n) are operators of type (7) with ρ = 1,

• t(n), t̄(n) are K-valued sequences with infinitely many nonzero values, defined on an infinite interval of integers I such that I ∩ [m_{L̄t} + 1, M_{L̄l} − 1] = ∅, and satisfying L t(n) = L̄ t̄(n) = 0 on I^{−},

• p(n) ∈ K[n] is a polynomial such that t(n) = p(n)t̄(n) on I,

• r(n) ∈ K(n) is a rational function which satisfies (10) as an equation in K(n),

• r̄(n) = p(n)r(n) in K(n),

• ū(n) is a K-sequence such that ū(n) = r̄(n)t̄(n) whenever both r̄(n) and t̄(n) are defined,

• v, w are integers such that v ≤ w and [v, w] ∩ ℤ ⊆ I.

Then equation (13) holds for all n ∈ I^{−}, and the discrete Newton-Leibniz formula

∑_{k=v}^{w} t(k) = ū(w) − ū(v) + t(w)  (19)

is valid.

Proof: By assumption, we have for all n ∈ I^{−},

ā_{1}(n)t̄(n + 1) + ā_{0}(n)t̄(n) = 0, (20)
a_{1}(n)p(n + 1)t̄(n + 1) + a_{0}(n)p(n)t̄(n) = 0. (21)

Multiplying (20) by a_{1}(n)p(n + 1), (21) by ā_{1}(n) and subtracting, we find ā_{0}(n)t̄(n)a_{1}(n)p(n + 1) = a_{0}(n)p(n)t̄(n)ā_{1}(n). Since t̄(n) has infinitely many nonzero values on I^{−}, this implies that

a_{0}(n)p(n)ā_{1}(n) = ā_{0}(n)a_{1}(n)p(n + 1)  (22)

holds infinitely often, hence also as an equation in K[n]. Multiplying (10) by ā_{1}(n)p(n), (22) by r(n+1), subtracting, and cancelling a_{1}(n) in K(n), we obtain

ā_{0}(n)r̄(n + 1) + ā_{1}(n)r̄(n) = −ā_{1}(n)p(n)  (23)

as an equation in K(n). It follows from Proposition 1 that r̄(n) ∈ K(n) has no integer poles outside the interval [m_{L̄t} + 1, M_{L̄l}−1]. Therefore r̄(n) and ū(n) are defined on I, r̄(n+1) and ū(n + 1) are defined on I^{−}, and (23) is valid for all n ∈ I^{−}.

Pick an n ∈ I^{−}. Multiplying (20) by r̄(n + 1), (23) by t̄(n) and subtracting, we obtain

ā_{1}(n)r̄(n + 1)t̄(n + 1) − ā_{1}(n)r̄(n)t̄(n) = ā_{1}(n)p(n)t̄(n).

If ā_{1}(n) ≠ 0 this reduces to (13). If ā_{1}(n) = 0 then, by assumption, ā_{0}(n) ≠ 0, so (20) implies t̄(n) = t(n) = ū(n) = 0 and (23) implies r̄(n + 1) = ū(n + 1) = 0, hence (13) holds in this case as well. So (13) holds for all n ∈ I^{−}, and the second assertion follows by summing (13) on n from v to w − 1.

Example 7. Let t(n) be the hypergeometric term from Example 3 which satisfies L t(n) = 0 for all n ∈ ℤ where L = (n − 5)(n − 3) E − (n − 4)(n − 1)n. Define p(n) = (n − 5)(n − 3)(n − 2) and

t̄(n) = \{ \begin{matrix}(n − 1)!, & n ≥ 2, \\ \frac{(−1)^{n}}{(−n)!} , & n ≤ 1. \end{matrix}

Then L̄ t̄(n) = 0 for all n ≠ 1 where L̄ = E − n. So, in the notation of Theorem 4, the maximal possible interval I is either (−∞, 1] ∩ ℤ or [2, ∞) ∩ ℤ. As L̄ has no leading singularities, M_{L̄l} = −∞ and [m_{L̄t} + 1, M_{L̄l} − 1] = ∅. With r(n) = (n − 6)/((n − 2)(n − 3)), r̄(n) = (n − 5)(n − 6) and ū(n) = r̄(n)t̄(n), all the assumptions of Theorem 4 are satisfied, and it follows that formula (19) is valid provided that w ≤ 1 or v ≥ 2.

Now we consider the general case with ρ ≥ 1.

Proposition 4. Let (r(n), R) be the result of applying AS to an operator L ∈ K[n, E] of type (7), and let r(n) = s(n)/q(n) where s(n), q(n) ∈ K[n]. Then there exist p(n) ∈ K[n] and L̄ ∈ K[n, E] such that

L ∘ p = qL̄,  (24)

(E − 1) ∘ R ∘ p = p − s L̄, and R ∘ p ∈ K[n, E].

Proof: Let d ∈ K[n] be a polynomial and B ∈ K[n, E] an operator such that

E^{ρ} ∘ L^{∗} ∘ \frac{1}{q} = \frac{1}{d} B.

Then

\frac{1}{q} L ∘ E^{−ρ} = B^{∗} ∘ \frac{1}{d} .

Therefore L ∘ E^{−ρ} ∘ d = qB^{∗}. Multiplying this by E^{ρ} on the right gives L ∘ E^{−ρ} ∘ d ∘ E^{ρ} = L ∘ d(n − ρ) = qB^{∗} ∘ E^{ρ}. Take p(n) = d(n − ρ) and L̄ = B^{∗} ∘ E^{ρ}. Then (24) is satisfied and p−s L̄ = p−rqL̄ = p−rL∘p = (1−rL)∘p = (E −1)∘R∘p. Hence the operator R ∘ p is the left quotient of p − rqL̄ by E −1 and, as a consequence, has polynomial coefficients.

Theorem 5. Let

• L, L̄, R, r, p, q be such as in Proposition 4,

• v, w ∈ ℤ be such that v ≤ w,

• I_{ρ} = [v, w − ρ] ∩ ℤ,

• t̄(n) be a K-valued sequence defined for all n ∈ [v, w] ∩ ℤ such that L̄t̄(n) = 0 for all n ∈ I_{ρ}.

Then the K-valued sequence t(n) = p(n)t̄(n) satisfies Lt(n) = 0 for all n ∈ I_{ρ}, and the discrete Newton-Leibniz formula (18) can be applied to t(n) with u(n) = (R ∘ p)t̄(n).

Proof: By (24),

Lt(n) = (L ∘ p)t̄(n) = qL̄t̄(n) = 0

for all n ∈ I_{ρ}. Also, u(n) = (R ∘ p)t̄(n) and u(n + 1) = (E ∘ R ∘ p)t̄(n) are defined for all n ∈ I_{ρ}. Therefore, by Proposition 4,

u(n + 1) − u(n) = ((E − 1) ∘ R ∘ p)t̄(n) = (p − sL̄)t̄(n)
= p(n)t̄(n) − s(n)L̄t̄(n) = t(n)

for all n ∈ I_{ρ}, and (18) follows by summing this over I_{ρ}.

Example 8. Take

L = 2(n + 1)(n − 2)E − (2n − 1)(n − 1).

Here m_{Lt} = 1 and M_{Ll} = 3, so by Theorem 3, we can apply formula (18) if [v, w − 1] ∩ [1, 2] = ∅. Using the algorithm from [1] or the algorithm from [2] we compute the solution r(n) = −(n + 2)/((n − 1)(n − 2)) of L^{∗}y = 1, and set q(n) = (n − 1)(n − 2). We have

E ∘ L^{∗} ∘ \frac{1}{(n − 2)(n − 1)} = \frac{1}{n − 1} (−(2n + 1)E + 2(n + 1)),

therefore

d(n) = n − 1,
B = −(2n + 1)E + 2(n + 1),
L̄ = B^{∗} ∘ E = 2(n + 1)E − (2n − 1), (25)
p(n) = n − 2.

Finally we get

L ∘ (n − 2) = (n − 1)(n − 2)L̄

with [m_{L̄t}, M_{L̄l} − ρ] = ∅.

We compute R = 2n(n + 1)/(n − 2) and u(n) = 2n(n + 1)t̄(n). If w ≥ v we have

∑_{n=v}^{w−1} t(n) = 2w(w + 1)t̄(w) − 2v(v + 1)t̄(v)  (26)

for any sequence t(n) annihilated by L which has the form (n − 2)t̄(n) with t̄(n) annihilated by (25) such that t̄(n) is defined for all n ∈ [v, w]. The equation L̄y = 0 has solutions of the form t̄(n) = c \frac{(−1/2)_{n}}{(1)_{n}} where c is an arbitrary constant. Each term of this form is defined for all n ∈ ℤ and we can use the formula (26) to sum, e.g., t(n) = (2 − n) \frac{(−1/2)_{n}}{2(1)_{n}},

with t̄(n) = \frac{−(−1/2)_{n}}{2(1)_{n}}.

Acknowledgements. The authors wish to express their thanks to Jürgen Gerhard for providing (a close variant of) the hypergeometric term in Example 1, and to anonymous referees for their helpful remarks.

7. REFERENCES

[1] S. A. Abramov, Rational solutions of linear difference and differential equations with polynomial coefficients, USSR Comput. Math. Phys. 29 (1989), 7–12. Transl. from Zh. vychisl. mat. mat. fyz. 29 (1989), 1611–1620.

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

[3] S. A. Abramov, M. van Hoeij, Integration of solutions of linear functional equations, Integral Transforms and Special Functions 8 (1999), 3–12.

[4] S. A. Abramov, M. Petkovšek, Solution spaces of H-systems and the Ore-Sato theorem, Proceedings FPSAC ’05, to appear.

[5] F. Chyzak, An extension of Zeilberger’s fast algorithm to general holonomic functions, Discrete Math. 217 (2000), 115–134.

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

[7] I. Nemes, M. Petkovšek, H. S. Wilf, D. Zeilberger, How to do Monthly problems with your computer, Amer. Math. Monthly 104 (1997), 505–519.

[8] M. Schorn, Contributions to Symbolic Summation, RISC, J. Kepler University, Diploma Thesis, December 1995.