Séminaire Lotharingien de Combinatoire 56 (2007), Article B56b

SYMBOLIC SUMMATION ASSISTS COMBINATORICS

CARSTEN SCHNEIDER

Abstract. We present symbolic summation tools in the context of difference fields that help scientists in practical problem solving. Throughout this article we present multi-sum examples which are related to combinatorial problems.

1. Introduction

At the 56^{th} Séminaire Lotharingien de Combinatoire I presented in a series of talks the summation package Sigma [Sch04b]. Starting with Karr’s summation algorithm [Kar81, Kar85, Sch00], the discrete analogue of Risch’s integration algorithm [Ris69, Ris70], I developed various extensions and generalizations to tackle non-trivial multi-sum problems. In this survey article, which is based on this presentation, I illustrate all these new features. The article consists of three parts. Each of them can be considered independently.

In the first part we introduce the summation principles of Sigma: “telescoping”, “creative telescoping” [Zei91], and “solving recurrences”. Multi-sum examples from [AU85, GKP94, Zha99, FK00, DPSW06a, DPSW06b, PS07] illustrate all these techniques.

In the second part we explain how the summation principles are formulated in difference fields and we demonstrate how the underlying algorithms work. We present our algorithmic extensions [Sch01, Sch04c, Sch04a, Sch05e, Sch05a, Sch05d, Sch05b, KS06b, KS06a, Sch07].

In the third part we show how Sigma can be applied for even a wider class of multi-sums which covers big parts of ∂-finite and holonomic sequences [Zei90, CS98, Chy00]. Examples from [AS65, PWZ96, APS05] illustrate how our methods [Sch05c] find recurrences for such sums. We demonstrate that one can derive also differential equations with these techniques.

I am very grateful to Volker Strehl and Christian Krattenthaler who gave me the opportunity to present my work at the Séminaire Lotharingien de Combinatoire. I would like to thank the referee for his very careful reading and his valuable comments and suggestions.

Part 1. Symbolic summation methods and applications

Inspired by [PWZ96] the summation package Sigma is based on the paradigms of telescoping, creative telescoping, and solving recurrences. We show in the frame of Sigma how one can apply these summation principles to attack multi-sum problems.

2. The summation principles

Given an indefinite sum S(n) = ∑_{k=0}^{n} f(k) we are interested in solving the following problem.

Telescoping

Then summing (1.1) over k (under the assumption that (1.1) holds for all 0 ≤ k ≤ n) gives

g(n + 1) − g(0) = ∑_{k=0}^{n} f(k).

Given f(k); find g(k) such that

g(k + 1) − g(k) = f(k).  (1.1)

Supported by the SFB-grant F1305 and the grant P16613-N12 of the Austrian FWF.

There are various algorithms on the market for such input sequences f(k), like Gosper’s algorithm [Gos78] for hypergeometric terms^{1}; for theoretical insight see [Pau95] and for implementations see [PS95a, Koe95, ACGL04]. Similarly, there are variations for q-hy-pergeometric terms [Koo93, PR97, BK99], or generalizations for the mixed case [BP99].

More generally, the summation package Sigma can deal with rational expressions in terms of indefinite nested sums and products. The underlying telescoping algorithms [Sch04c, Sch05b, Sch07] extend Karr’s summation algorithm [Kar81].

Sigma session. Consider the indefinite sum S(n) = ∑_{k=m}^{n} f(k) with the summand f(k) = \left[{k \atop m}\right]_{q}q^{k} ∑_{i=1}^{k} \frac{q^{i}}{1−q^{i}}; \left[{k \atop m}\right]_{q} denotes the q-binomial. Then we can simplify S(n) as follows. After loading Sigma into the computer algebra system Mathematica

In[1]:= << Sigma.m

Sigma - A summation package by Carsten Schneider \begin{matrix}c \\ ◯ \end{matrix} RISC-Linz

we insert S(n)(=mySum):

In[2]:= mySum = SigmaSum[SigmaqBinomial[q, k, m]SigmaPower[q, k] SigmaSum[SigmaPower[q, i]/(1 − SigmaPower[q, i]), {i, 1, k}], {k, m, n}]

Out[2]= ∑_{k=m}^{n}\left[{k \atop m}\right]_{q}q^{k} ∑_{i=1}^{k} \frac{q^{i}}{1 − q^{i}}

Sigma manual. The basic functions SigmaSum and SigmaProduct are used to describe all indefinite nested sums and products that can be expressed in our setting. Additional functions, like SigmaPower, SigmaBinomial, or SigmaqBinomial are just shortcuts. E.g., our q-binomial can be described by \left[{k \atop m}\right]_{q} = ∏_{i=m+1}^{k} \frac{1−∏_{j=1}^{i} q}{1−∏_{j=1}^{i} q^{i}/q^{m}} for k ≥ m.

Then we find the following closed form by the function call

In[3]:= SigmaReduce[mySum]

Out[3]= \frac{q^{m}}{(qq^{m} − 1)^{2}} (q(q^{n} − q^{m}) + (qq^{m} − 1)(qq^{n} − 1) ∑_{i=1}^{n} \frac{q^{i}}{1 − q^{i}} )\left[{n \atop m}\right]_{q}

Internally, Sigma computes

g(k) = \frac{q^{m} − q^{k}}{(qq^{m} − 1)^{2}} ( − qq^{m} + (1 − qq^{m}) ∑_{i=1}^{k} \frac{q^{i}}{1 − q^{i}} )\left[{k \atop m}\right]_{q}

which satisfies (1.1) for all m ≤ k ≤ n. Note that the correctness can be verified independently of the computation steps. Namely, by using the relations ∑_{i=1}^{k+1} \frac{q^{i}}{1−q^{i}} =

∑_{i=1}^{k} \frac{q^{i}}{1−q^{i}} + \frac{qq^{k}}{1−qq^{k}} and \left[{k+1 \atop m}\right]_{q} = \frac{1−qq^{k}}{1−qq^{k}/q^{m}} \left[{k \atop m}\right]_{q} we can verify with simple polynomial arithmetic that (1.1) holds for all 0 ≤ m ≤ k ≤ n. Hence summing (1.1) over k produces

∑_{k=m}^{n} \left[{k \atop m}\right]_{q}q^{k} ∑_{i=1}^{k} \frac{q^{i}}{1 − q^{i}} = q^{m}\left[{n + 1 \atop m + 1}\right]_{q}( ∑_{i=1}^{n+1} \frac{q^{i}}{1 − q^{i}} − \frac{q^{m+1}}{1 − q^{m+1}} )

which is equivalent to Out[3]; compare identity [AU85, (2.5)]. Multiplying by 1 − q and sending q to 1 yields the identity [GKP94, (6.70)]

∑_{k=m}^{n} \binom{k}{m}H_{k} = \binom{n + 1}{m + 1}(H_{n+1} − \frac{1}{m + 1})

where H_{k} = ∑_{i=1}^{k} \frac{1}{i} denotes the harmonic numbers.

1f(k) is hypergeometric (resp. q-hypergeometric), if \frac{f(k+1)}{f(k)} is a rational function in k (resp. in q^{k}).

In most cases one fails to find such a g(k). If f(k) depends on an additional parameter n, we can apply a more flexible paradigm.

Zeilberger’s creative telescoping

Given f(n, k) and δ ∈ ℕ; find c_{0}(n), …, c_{δ}(n), free of k, and g(n, k) such that

g(n, k + 1) − g(n, k) = c_{0}(n)f(n, k) + ⋯ + c_{δ}(n)f(n + δ, k).  (1.2)

With creative telescoping we can attack definite sums as follows: Given a sum S(n) = ∑_{k=0}^{m} f(n, k) where m might depend linearly on n, find a solution c_{0}(n), …, c_{δ}(n) and g(n, k) of (1.2); here one usually starts with δ = 0 (which is nothing else than telescoping), and increases δ step by step. If one finds such a solution for (1.2) which holds for all 0 ≤ k ≤ n, then summing (1.2) over k gives a recurrence relation of the form

q(n) = c_{0}(n)S(n) + ⋯ + c_{δ}(n)S(n + δ).

As it turns out, all the telescoping algorithms mentioned above can be extended to creative telescoping. This was observed the first time by Zeilberger [Zei91] for Gosper’s algorithm; for theoretical insight see [PS95b] and for implementations and additional details see [PS95a, ACGL04, Ger04]. For the q-hypergeometric, the mixed and the holonomic case see [Koo93, PR97, BK99, BP99, Chy00]; different generalizations can be found in [M06].

Moreover, as recognized in [Sch00], Karr’s algorithm and all our extensions can be used for creative telescoping; see Sections 6.2 and 7.1.3.

Sigma session. Consider the definite sum S(n) given by

In[4]:= hSum = ∑_{k=0}^{n} \binom{n}{k}H_{k};

where the binomial \binom{n}{k} is interpreted as ∏_{i=1}^{k} \frac{n+1−i}{i}. By typing in

In[5]:= hRec = GenerateRecurrence[hSum]

Out[5]= {4(1 + n)SUM[n] − 2(3 + 2n)SUM[n + 1] + (2 + n)SUM[n + 2] == 1}

we compute the recurrence hRec for S(n) =SUM[n]. Note that by the proof certificate

In[6]:= CreativeTelescoping[hSum][[2]]

Out[6]= \{4 (1 + n), −2 (3 + 2n), 2 + n, \frac{(1+n) (−2+k−n+(2k−2k^{2}+kn)H_{k}) \binom{n}{k})}{(1−k+n) (2−k+n)}\}

we can conclude that the recurrence is correct: Given f(n, k) := \binom{n}{k}H_{k} and

\begin{aligned}c_{0}(n) & := 4 (1 + n), c_{1}(n) := −2 (3 + 2 n), c_{2}(n) := 2 + n, \\ g(n, k) & := \frac{(1+n) (−2+k−n+(2 k−2 k^{2}+k n)H_{k}) \binom{n}{k})}{(1−k+n) (2−k+n)}, \end{aligned}  (1.3)

we show that (1.2) with δ = 2 holds for all 0 ≤ k ≤ n as follows. Represent in (1.2) the expressions g(n, k + 1) and f(n + i, k) for i = 1, 2 in terms of H_{k} and \binom{n}{k}. This is possible by H_{k+1} = H_{k} + 1/(k + 1) and the fact that \binom{n}{k} is hypergeometric in n and k. E.g., we can write g(n, k + 1) = \frac{(n+1)(nH_{k}−2kH_{k}−1)}{n−k+1}. Given this representation, verify (1.2) by polynomial arithmetic. Summing (1.2) over k gives Out[5].

Similarly, we can compute for the following q-version of S(n) a recurrence:

In[7]:= qhSum = ∑_{k=0}^{n} q^{k(k−1)/2}\left[{n \atop k}\right]_{q} ∑_{i=1}^{k} \frac{q^{i}}{1 − q^{i}} ;

In[8]:= qhRec = GenerateRecurrence[qhSum]

Out[8]= {(q^{n+1} + q) (q^{n+1} − 1) (q^{n+1} + 1) SUM[n] −

             (q^{n+1} + 1) (2q^{n+2} − q − 1) SUM[n + 1] + (q^{n+2} − 1) SUM[n + 2] == −q^{2n+3}}

Remark. For the general input class of Sigma no criteria are known which guarantee the existence of a creative telescoping solution (1.2) for sufficiently large δ. So far there is only one way: try it out. For (q-)hypergeometric terms necessary and sufficient conditions are available; see [Abr03, AL05, CHM05]. Restricting to proper hypergeometric terms upper bounds for δ are known which guarantee solutions; see [PWZ96, Thm. 4.4.1] and [MZ05].

In many applications a recurrence is just the desired result, see e.g. Out[13]. If one is interested in a more explicit representation, the following attempt might help.

Solving recurrences

Given a recurrence

a_{0}(n)S(n) + ⋯ + a_{r}(n)S(n + r) = f(n);  (1.4)

find a particular solution p(n) and solutions of the homogeneous version, say h_{1}(n), …, h_{d}(n).

If we compute sufficiently^{2} many solutions, we can find constants k_{0}, …, k_{d} such that S(i) = k_{0}h_{0}(i) + ⋯ + k_{d}h_{d}(i) + p(i) holds for i ∈ {0, …, d − 1}. This implies that

S(n) = k_{0}h_{0}(n) + ⋯ + k_{d}h_{d}(n) + p(n), n ≥ 0.  (1.5)

Sigma session. We solve the recurrence hRec from Out[5] by typing in

In[9]:= recSol = SolveRecurrence[hRec[[1]], SUM[n]]

Out[9]= \{\{0, ∏_{i=1}^{n}2\}, \{0, ( ∏_{i=1}^{n}2) ∑_{i=1}^{n} \frac{1}{i}\}, \{1, −( ∏_{i=1}^{n}2) ∑_{i=1}^{n} \frac{1}{i ∏_{j=1}^{i}2}\}\}

Sigma manual. The 0s in the first two entries tell us that h_{0}(n) = 2^{n} and h_{1}(n) = 2^{n} ∑_{i=1}^{n} H_{n} are solutions of the homogeneous version of hRec, and the 1 in the third entry indicates that p(n) = −2^{n} ∑_{i=1}^{n} \frac{1}{2^{i}i} is a particular solution of hRec itself.

Finally, using the first two initial values we combine the solutions in the form (1.5):

In[10]:= FindLinearCombination[recSol, hSum, 2]

Out[10]= ( ∏_{i=1}^{n}2) ∑_{i=1}^{n} \frac{1}{i} − ( ∏_{i=1}^{n}2) ∑_{i=1}^{n} \frac{1}{i ∏_{j=1}^{i}2}

This leads to [PS03, Equ. (39)]

∑_{k=0}^{n} \binom{n}{k}H_{k} = 2^{n}H_{n} − 2^{n} ∑_{i=1}^{n} \frac{1}{2^{i}i}.  (1.6)

Completely analogously we solve the recurrence qhRec from Out[8] by the function call

In[11]:= recSol = SolveRecurrence[qhRec[[1]], SUM[n]]

Out[11]= \{\{0, ∏_{i=1}^{n} \frac{q + q^{i}}{q} \}, \{0, ( ∏_{i=1}^{n} \frac{q + q^{i}}{q} ) ∑_{i=1}^{n} \frac{q^{i}}{1 − q^{i}} \}, \{1, −( ∏_{i=1}^{n} \frac{q + q^{i}}{q} ) ∑_{i=1}^{n} \frac{q^{i}}{(1 − q^{i}) ∏_{j=1}^{i} \frac{q+q^{j}}{q}} \}\}

and find the following combination for our sum In[7]:

In[12]:= FindLinearCombination[recSol, qhSum, 2]

Out[12]= ( ∏_{i=1}^{n} \frac{q + q^{i}}{q} ) ∑_{i=1}^{n} \frac{q^{i}}{1 − q^{i}} − ( ∏_{i=1}^{n} \frac{q + q^{i}}{q} ) ∑_{i=1}^{n} \frac{q^{i}}{(1 − q^{i}) ∏_{j=1}^{i} \frac{q+q^{j}}{q}}

With the q-shifted factorial (−1, q)_{n} = ∏_{i=1}^{n} \frac{q+q^{i}}{q} we arrive at the identity

∑_{k=0}^{n} q^{\binom{k}{2}}\left[{n \atop k}\right]_{q} ∑_{i=1}^{k} \frac{q^{i}}{1 − q^{i}} = (−1, q)_{n} (∑_{i=1}^{n} \frac{q^{i}}{1 − q^{i}} − ∑_{i=1}^{n} \frac{q^{i}}{(−1, q)_{i}(1 − q^{i})});  (1.7)

note that we rediscover identity (1.6) by multiplying it with 1 − q and sending q to 1.

2If a_{0}(n), a_{r}(n) have only finitely many zeros, there are exactly r linearly independent solutions of (1.4)

with f(n) = 0, see [PWZ96, Thm. 8.2.1]; here we identify sequences, if they agree from some point on.

In Sigma the coefficients a_{i}(n) and the inhomogeneous part f(n) of the recurrence (1.4) can be rational expressions in terms of nested sums and products. Given such a recurrence, we can find the class of d’Alembertian solutions [AP94], a subclass of Liouvillian solutions [HS99]; typical examples are given in Out[9] and Out[11]. Here the crucial point is that those solutions are indefinite nested sums and products which fit in the input class of our telescoping algorithms. Hence, given such a representation (1.5), Sigma can help to simplify its right-hand side. If we are lucky, we can end up with a “closed form” for S(n).

The interaction of the different summation principles for a definite sum can be summarized by the Sigma-summation spiral [Sch04b]:

[figure: definite sum; creative telescoping; simplified solutions; combination of solutions; recurrence; solving; d’Alembertian solutions; telescoping]

3. Applications of the Sigma-summation spiral

Subsequently, we will illustrate all the aspects from Section 2 by concrete examples.

3.1. Quadratic Padé approximation. In [Wei05] A. Weideman looks for polynomials r_{n}(x), s_{n}(x) and t_{n}(x) with degree at most n such that

r_{n}(x) (log x)^{2} + s_{n}(x) log(x) + t_{n}(x) = O((x − 1)^{3n+2}).  (1.8)

He discovers that those polynomials can be written as a linear combination of the polynomials

A_{n}(x) = ∑_{k=0}^{n} \binom{n}{k}^{3}(−x)^{k}, B_{n}(x) = ∑_{k=0}^{n} [\frac{d}{dk} \binom{n}{k}^{3}](−x)^{k}, C_{n}(x) = ∑_{k=0}^{n} [\frac{d^{2}}{dk^{2}} \binom{n}{k}^{3}](−x)^{k}.

Evaluations at x = 1 and n = 0, 1, ... show empirically how the polynomials A_{n}(x), B_{n}(x) and C_{n}(x) must be combined to get r_{n}(x), s_{n}(x) and t_{n}(x). To obtain a rigorous proof for the guessed representation, it turns out that one has to show the key identity

∑_{k=0}^{n}(−1)^{k}\binom{n}{k}^{3}(3(H_{n−k} − H_{k})^{2} + H_{n−k}^{(2)} + H_{k}^{(2)}) = 0;  (1.9)

H_{k}^{(2)} = ∑_{i=1}^{k} \frac{1}{i^{2}} are the harmonic numbers of second order. In [DPSW06a] Sigma played the main role in proving this identity: Namely, we can compute the following recurrence

In[13]:= GenerateRecurrence[∑_{k=0}^{n}(−1)^{k}\binom{n}{k}^{3}(3(H_{n−k} − H_{k})^{2} + H_{n−k}^{(2)} + H_{k}^{(2)})]

Out[13]= {(n + 2)^{2}SUM[n + 2] + 3(3n + 2)(3n + 4)SUM[n] == 0}

for the left-hand side of (1.9). Since this sum is zero for n = 0, 1, it must evaluate to zero for all n ≥ 0. We remark that together with the creative telescoping solution given in [DPSW06a] the correctness of the recurrence can be checked as for (1.3). Later C. Krattenthaler [Kra03] derived a non-algorithmic proof based on differentiation and hypergeometric transformations. W. Chu [Chu05] presents different techniques to show (1.9).

As illustrated in [DPSW06b] we cannot only prove identity (1.9), but we can discover it. To be more precise, we find with our machinery the identities

∑_{k=0}^{2n}(−1)^{k}\binom{2n}{k}^{3}H_{k}H_{2n−k} =\frac{(3n)!(−1)^{n}}{n!n!n!} \frac{1}{12}(3H_{n}^{2} − 6H_{n}H_{3n} + 3H_{3n}^{2} + H_{n}^{(2)}
+ 12H_{2n}(H_{2n} + H_{n} − H_{3n}) + 4H_{2n}^{(2)} − 3H_{3n}^{(2)}),  (1.10)
∑_{k=0}^{2n}(−1)^{k}\binom{2n}{k}^{3}H_{k}^{2} =\frac{(3n)!(−1)^{n}}{n!n!n!} \frac{1}{12}(3H_{n}^{2} − 6H_{n}H_{3n} + 3H_{3n}^{2} − H_{n}^{(2)}
+ 12H_{2n}(H_{2n} + H_{n} − H_{3n}) + 2H_{2n}^{(2)} − 3H_{3n}^{(2)}),  (1.11)
∑_{k=0}^{2n}(−1)^{k}\binom{2n}{k}^{3}H_{k}^{(2)} =\frac{1}{2} \frac{(3n)!(−1)^{n}}{n!n!n!} (H_{n}^{(2)} + H_{2n}^{(2)}).  (1.12)

Then it is easy to see that the right combination of these sums, see (1.9), evaluates to zero. E.g., we discover identity (1.10) as follows. Given the sum

In[14]:= mySum = ∑_{k=0}^{2n}(−1)^{k}\binom{2n}{k}^{3}H_{k}H_{2n−k};

we compute the recurrence relation

In[15]:= rec = GenerateRecurrence[mySum, SimplifyByExt → DepthNumber][[1]]

Out[15]= −18(n + 1)^{2}(n + 2)(2n + 1)^{2}(2n + 3)(3n + 1)^{2}(3n + 2)^{2}(108n^{3} + 495n^{2} + 752n + 378)SUM[n] +
      6(n + 1)^{3}(n + 2)(2n + 1)(2n + 3)(3888n^{7} + 29484n^{6} + 92250n^{5} + 153369n^{4} + 145192n^{3} + 77561n^{2} +
      21420n + 2316)SUM[n + 1] + 2(n + 1)^{3}(n + 2)^{4}(2n + 1)(2n + 3)^{2}(108n^{3} + 171n^{2} + 86n + 13)SUM[n + 2]
   == (2519424n^{11} + 26873856n^{10} + 126618552n^{9} + 347114322n^{8} + 613953513n^{7} + 734258088n^{6} +
        604816090n^{5}+342574260n^{4}+130558875n^{3}+31842320n^{2}+4469856n+273984) ∑_{k=0}^{2n}(−1)^{k}\binom{2n}{k}^{3}

Sigma manual. With the option SimplifyByExt→DepthNumber we search for a recurrence with sum extensions which have less objects in the summand than the input summand. In our case we find a recurrence of order two by using the sum e(n) = ∑_{k=0}^{2n}(−1)^{k}\binom{2n}{k}^{3}. Without this option, i.e., without using e(n), we compute a rather big recurrence relation of order 3.

Since the automatically found sum e(n) is definite, it does not fit into the input class of recurrences that Sigma can handle. But with our machinery, as described above, we find e(n) = \frac{(−1)^{n}(3n)!}{n!^{3}} which is a particular instance of Dixon’s identity [Dix03]; see also [PWZ96, Example 6.4.4]. Hence, we can simplify the recurrence by replacing e(n) with \frac{(−1)^{n}(3n)!}{n!^{3}}:

In[16]:= rec = rec/. ∑_{k=0}^{2n}(−1)^{k}\binom{2n}{k}^{3} → \frac{(−1)^{n}(3n)!}{n!^{3}};

This recurrence can be handled with Sigma. Namely, we compute the solutions

In[17]:= recSol = SolveRecurrence[rec, SUM[n]]

Out[17]= \{\{0, \frac{(−1)^{n}(3n)!}{n!^{3}}\}, \{0, \frac{(−1)^{n}(3n)!}{n!^{3}} ∑_{k=1}^{n} \frac{108k^{3} − 153k^{2} + 68k − 10}{k(2k − 1)(3k − 2)(3k − 1)}\}, \{1, \frac{1}{18}\frac{(−1)^{n}(3n)!}{n!^{3}}

∑_{k=1}^{n}^{\frac{−4860k^{6}+13770k^{5}−15849k^{4}+9504k^{3}−3148k^{2}+550k−40+(k(2k−1)(3k−2)(3k−1)(108k^{3}−153k^{2}+68k−10))∑_{i=1}^{k} \frac{108i^{3}−153i^{2}+68i−10}{i(2i−1)(3i−2)(3i−1)}}{k^{2}(2k−1)^{2}(3k−2)^{2}(3k−1)^{2}}} \}\}

Here it is important to mention that this type of solutions (d’Alembertian solutions) are just in the input class of Sigma’s telescoping algorithms; we get the following simplification:

In[18]:= recSol = SigmaReduce[recSol, n, SimplifyByExt → Depth]

Out[18]= \{\{0, \frac{(−1)^{n}(3n)!}{n!^{3}}\}, \{0, \frac{(−1)^{n}(3n)!}{n!^{3}} ∑_{k=1}^{n} \frac{108k^{3} − 153k^{2} + 68k − 10}{k(2k − 1)(3k − 2)(3k − 1)}\},

      \{1, \frac{1}{36}\frac{(−1)^{n}(3n)!}{n!^{3}}(( ∑_{i=1}^{n} \frac{108i^{3}−153i^{2}+68i−10}{i(2i−1)(3i−2)(3i−1)})^{2} + ∑_{k=1}^{n} \frac{1944k^{6}−5508k^{5}+6399k^{4}−3960k^{3}+1388k^{2}−260k+20}{k^{2}(2k−1)^{2}(3k−2)^{2}(3k−1)^{2}})\}\}

Séminaire Lotharingien de Combinatoire 56 (2007), Article B56b 7

Sigma manual. With the option SimplifyByExt→Depth we look for appropriate sum extensions in order to reduce the nested depth; see Section 7.1. E.g., the second sum in Out[18] has been computed in order to represent the expression in Out[17] by single nested sums. Using the first two initial values we combine the solutions to represent In[14] as follows.

In[19]:= closedForm = FindLinearCombination[recSol, mySum, 2]

Out[19]= \frac{1}{36}\frac{(−1)^{n}(3n)!}{n!^{3}}(( ∑_{i=1}^{n} \frac{108i^{3}−153i^{2}+68i−10}{i(2i−1)(3i−2)(3i−1)})^{2} + ∑_{k=1}^{n} \frac{1944k^{6}−5508k^{5}+6399k^{4}−3960k^{3}+1388k^{2}−260k+20}{k^{2}(2k−1)^{2}(3k−2)^{2}(3k−1)^{2}})

Next, we split our sums so that their denominators consist of irreducible polynomials:

In[20]:= closedForm = SigmaReduce[closedForm, n, SimpleSumRepresentation → True]

Out[20]= \frac{(−1)^{n}(3n)!}{n!^{3}}( − \frac{1}{4} ∑_{k=1}^{n} \frac{1}{(3k − 1)^{2}} − \frac{1}{4} ∑_{k=1}^{n} \frac{1}{(3k − 2)^{2}} + \frac{1}{3} ∑_{k=1}^{n} \frac{1}{(2k − 1)^{2}} + \frac{5}{36} ∑_{k=1}^{n} \frac{1}{k^{2}} +

            (\frac{5}{6} ∑_{k=1}^{n} \frac{1}{k} + ∑_{k=1}^{n} \frac{1}{2k − 1} − \frac{1}{2} ∑_{k=1}^{n} \frac{1}{2k − 2} − \frac{1}{12} ∑_{k=1}^{n} \frac{1}{2k − 3})^{2})

By some cosmetic rewriting we end up with the right-hand side of (1.10):

In[21]:= closedForm = SigmaReduce[closedForm, n, Tower → {H_{n}, H_{2n}, H_{n}^{(2)}, H_{2n}^{(2)}, H_{3n}, H_{3n}^{(2)}}]

Out[21]= \frac{(−1)^{n}(3n)!}{n!^{3}} \frac{1}{12}(3H_{n}^{2} + 12H_{2n}H_{n} + 12H_{2n}^{2} + 3H_{3n}^{2} + (−6H_{n} − 12H_{2n}) H_{3n} + H_{n}^{(2)} + 4H_{2n}^{(2)} − 3H_{3n}^{(2)})

3.2. A variation of Calkin’s sum. Alternating versions of Calkin’s identity [Cal94]

∑_{k=0}^{n} ( ∑_{j=0}^{k} \binom{n}{j})^{3} = \frac{n}{2} 8^{n} + 8^{n} − \frac{3 n}{4} 2^{n} \binom{2 n}{n}  (1.13)

have been considered in [Zha99]. We supplement this collection with the following sum.

In[22]:= mySum = ∑_{k=0}^{2n}(−1)^{k}( ∑_{j=0}^{k} \binom{2n}{j})^{3};

Following our Sigma-spiral we compute a recurrence

In[23]:= rec = GenerateRecurrence[mySum][[1]]

Out[23]= = (2n + 3)(11n + 3)(n + 1)^{2}SUM[n + 2] + (946n^{4} + 2799n^{3} + 2907n^{2} + 1252n + 180) SUM[n + 1] +
 96n(3n + 1)(3n + 2)(11n + 14)SUM[n] == 16 (5005n^{4} + 16897n^{3} + 20210n^{2} + 9884n + 1512) ( ∑_{k=0}^{2n}\binom{2n}{k})^{3}

In[24]:= rec = rec/.( ∑_{k=0}^{2n} \binom{2n}{k})^{3} → 64^{n};

solve the recurrence in terms of d’Alembertian solutions

In[25]:= recSol = SolveRecurrence[rec, SUM[n]]

Out[25]= \{\{0, ∏_{i=2}^{n} −\frac{32(i − 1)}{2i − 1}\}, \{0, ( ∏_{i=2}^{n} −\frac{32(i − 1)}{2i − 1}) ∑_{i=1}^{n} \frac{i^{3}(11i − 8) ∏_{j=1}^{i} \frac{3(2j−1)(3j−2)(3j−1)}{32j^{3}}}{(2i − 1)(3i − 2)(3i − 1)} \}, \{1, \frac{64^{n}}{2} \}\}

and combine the solutions

In[26]:= FindLinearCombination[recSol, mySum, n, 2, MinInitialValue → 1]

Out[26]= \frac{64^{n}}{2} + \frac{64}{3}( ∏_{i=2}^{n} −\frac{32(i − 1)}{2i − 1}) ∑_{i=1}^{n} \frac{i^{3}(11i − 8) ∏_{j=1}^{i} \frac{3(2j−1)(3j−2)(3j−1)}{32j^{3}}}{(2i − 1)(3i − 2)(3i − 1)}

Finally, after some rewriting we arrive, for n ≥ 1, at the identity

∑_{k=0}^{2 n}(−1)^{k}( ∑_{j=0}^{k} \binom{2 n}{j})^{3} = \frac{64^{n}}{2} − \frac{(−1)^{n}}{16 n} \frac{64^{n}}{\binom{2 n}{n}} ∑_{i=0}^{n−1}(3 + 11 i)\binom{2 i}{i}^{2}\binom{3 i}{i}64^{−i}. (1.14)

3.3. A problem from the SLC’06. At the SLC’06 Wenchang Chu showed in his talk various harmonic number identities of the type (1.10), (1.11), (1.12) which he produced by differentiating the Dougall-Dixon formula. Since he could not find an evaluation of the sum S(n) := ∑_{k=0}^{2n+1}(−1)^{k}\binom{2n+1}{k}^{3}H_{k}^{(2)} in terms of harmonic numbers, he posed this problem to Sigma. After some seconds of computations Sigma could give the following answer:

S(n) = \frac{(−1)^{n}(6n + 1)(n!)^{3}(6n)!}{2(2n + 1)^{2}((2n)!)^{3}(3n)!}( − 2 + ∑_{i=1}^{n} \frac{(72i^{2} + 36i + 5)((2i)!)^{3}((3i)!)^{2}}{2i(2i + 1)(6i + 1)(i!)^{6}(6i)!} ); (1.15)

observe that the sum on the right-hand side cannot be expressed by harmonic numbers. The derivation of the identity is based on the Sigma-spiral: First we compute a recurrence

− 12(2n + 3)(6n + 5)(6n + 7)(n + 1)^{3}S(n) − 4(2n + 3)^{3}(n + 1)^{3}S(n + 1)

= 3(3n + 1)(3n + 2) (72n^{2} + 180n + 113) \frac{(−1)^{n}(3n)!}{(n!)^{3}},

then we solve this recurrence and find its right-hand side.

3.4. A problem from rhombus tilings. Define

S_{n} := ∑_{k=0}^{n−1} \frac{(−1)^{n+k}(n + k + 4)!H_{k+1}}{(k + 2)!(k + 3)!(n − k − 1)!} and T_{n} := ∑_{k=0}^{n−1} \frac{(−1)^{k}(n + k + 4)!}{(k + 1)(k + 2)!^{2}(n − k − 1)!}.

Then in about four pages of highly non-trivial transformations the following evaluation has been found in [FK00, Lemma 26]:

S_{n} = \frac{−5 − 3n}{(1 + n)(2 + n)} − 2H_{n} + (−1)^{n}(\frac{5 + 2n − 2n^{2} − n^{3}}{(1 + n)(2 + n)} + 2(2 + n)H_{n}),
T_{n} = 1 − 9n − 9n^{2} − 2n^{3} + 2(1 + n)(2 + n)(3 + n)H_{n} − (−1)^{n}.

This finally shows that S_{n} + \frac{(1−(−1)^{n}(n+2))n!}{(n+3)!}T_{n} = (−1)^{n}(n + 2) − 2.

As illustrated in [Sch04b] Sigma finds these results in a straightforward manner by following the Sigma-spiral. Here it is worthwhile to mention that (−1)^{n} pops up with the algebraic relation ((−1)^{n})^{2} = 1; for more details see Sections 5.2 and 7.3.

3.5. Evaluation of a quadruple sum. In 2004 Doron Zeilberger sent an email to Robin Pemantle and Herbert Wilf with Cc to me:

Robin and Herb, I am willing to bet that Carsten Schneider’s SIGMA package for handling sums with harmonic numbers (among others) can do it in a jiffy. I am Cc-ing this to Carsten. Carsten: please do it, and Cc- the answer to me. -Doron

Of course, I and Sigma were eager to win the bet for Doron. So, “we” looked at Robin’s problem attached in this email which reads as follows:

I have a sum that, when I evaluate numerically, looks suspiciously like it comes out to exactly 1. Is there a way I can automatically decide this? The sum may be written in many ways, but one is:

S := ∑_{j,k=1}^{∞} \frac{H_{j}(H_{k+1} − 1)}{jk(k + 1)(j + k)}.

One week later I could reply [Sch06]: the sum is not exactly 1, but

S = −4ζ(2) − 2ζ(3) + 4ζ(2)ζ(3) + 2ζ(5) = 0.999222...  (1.16)

Whereas the full details are given in [PS07] we present here only the Sigma-part. Take the truncated version

S(a, b) = ∑_{k=1}^{b} \frac{H_{k+1} − 1}{k(k + 1)} ∑_{j=1}^{a} \frac{H_{j}}{j(j + k)},

Séminaire Lotharingien de Combinatoire 56 (2007), Article B56b 9

i.e., S = lim_{a,b→∞} S(a, b). Then Sigma can find a more appropriate sum representation of the inner sum following the Sigma-spiral. Namely, we compute a recurrence

In[27]:= rec = GenerateRecurrence[∑_{j=1}^{a} \frac{H_{j}}{j(j + k)}, k][[1]]

Out[27]= (k + 1)(a + k + 1)(a + k + 2)k^{2}SUM[k] − (k + 1)^{2}(a + k + 1)(a + k + 2)(2k + 1)SUM[k + 1] +

(k + 1)^{2}(k + 2)(a + k + 1)(a + k + 2)SUM[k + 2] == a(a + k + 2) + (−a − 1)(k + 1)H_{a}

solve the recurrence

In[28]:= recSol = SolveRecurrence[rec, SUM[k]]

Out[28]= \{\{0, \frac{1}{k}\}, \{0, \frac{∑{k \atop i=2} \frac{1}{i−1}}{k} \}, \{1, \frac{1}{k} ∑_{j=2}^{k} \frac{∑{j \atop i=2} \frac{−(a+1)H_{a}(i−1)+a(a+i)}{(−1+i)(−1+a+i)(a+i)}}{j − 1} \}

simplify the solutions

In[29]:= SigmaReduce[recSol, k, SimpleSumRepresentation → True, SimplifyByExt → Depth]

Out[29]= \{\{0, \frac{1}{k}\}, \{0, \frac{∑{k \atop i=2} \frac{1}{i−1}}{k} \}, \{1, \frac{H_{a}}{(a + 1)k} − \frac{(kH_{a} − 1)}{k^{2}} ∑_{i=1}^{k} \frac{1}{a + i} − \frac{1}{k} ∑_{i=1}^{k} \frac{1}{i} ∑_{j=1}^{i} \frac{1}{a + j} +

           \frac{∑_{i=2}^{k} \frac{1}{(i−1)^{2}} + ( ∑_{i=2}^{k} \frac{1}{i−1})^{2}}{2k} \}

and combine the solutions to get the identity

∑_{j=1}^{a} \frac{H_{j}}{j(j + k)} = \frac{kH_{k}^{2} − 2H_{k} + kH_{k}^{(2)} + 2kH_{a}^{(2)}}{2k^{2}} − \frac{(kH_{a} − 1)}{k^{2}} ∑_{i=1}^{k} \frac{1}{a + i} − \frac{1}{k} ∑_{i=1}^{k} \frac{1}{i} ∑_{j=1}^{i} \frac{1}{a + j} .

By simple limit considerations it follows lim_{a,b→∞} S′(a, b) = S for

S′(a, b) := ∑_{k=1}^{b} \frac{H_{k+1} − 1}{k(k + 1)} \frac{kH_{k}^{2} − 2H_{k} + kH_{k}^{(2)} + 2kH_{a}^{(2)}}{2k^{2}} .

Sigma can simplify S′(a, b) further to S′(a, b) = A(a, b) + B(a, b) + C(a, b) where

A(a, b) :=\frac{1}{2(b + 1)^{2}} (6H_{b} + 4bH_{b} + 4H_{b}^{2} + 3bH_{b}^{2} + H_{b}^{3} + bH_{b}^{3} − 6bH_{a}^{(2)}
+ 2H_{b}H_{a}^{(2)} + 2bH_{b}H_{a}^{(2)} − 2H_{b}^{(2)} − 7bH_{b}^{(2)} + H_{b}H_{b}^{(2)} + bH_{b}H_{b}^{(2)}),
B(a, b) := − \frac{2b^{2}}{(b + 1)^{2}} (H_{a}^{(2)} + H_{b}^{(2)}),
C(a, b) :=(H_{a}^{(2)} − 1) ∑_{i=1}^{b} \frac{H_{i}}{i^{2}} − ∑_{i=1}^{b} \frac{H_{i}^{2}}{i^{3}} + \frac{1}{2} ∑_{i=1}^{b} \frac{H_{i}^{3}}{i^{2}} + \frac{1}{2} ∑_{i=1}^{b} \frac{H_{i}H_{i}^{(2)}}{i^{2}} .

Since lim_{a,b→∞} A(a, b) = 0 and lim_{a,b→∞} B(a, b) = −4ζ(2), we get

S = lim_{a,b→∞} S′(a, b) = −4ζ(2) + lim_{a,b→∞} C(a, b).

Now we can use, e.g. [BG96, FS98], and find ∑_{i=1}^{∞} \frac{H_{i}}{i^{2}} = 2ζ(3), ∑_{i=1}^{∞} \frac{H_{i}^{2}}{i^{3}} = −ζ(2)ζ(3) + \frac{7}{2}ζ(5), ∑_{i=1}^{∞} \frac{H_{i}^{3}}{i^{2}} = ζ(2)ζ(3)+10ζ(5), and ∑_{i=1}^{∞} \frac{H_{i}H_{i}^{(2)}}{i^{2}} = ζ(2)ζ(3)+ζ(5). This proves (1.16).

In [PP05] a computer-free proof is given. Conversely, in [PS07] we show that the evaluation of such sums can be handled by a more systematic approach to computer-proofs.

Part 2. Summation in difference fields

This part is directed to readers, who are curious how Sigma works. In Sections 4–6 we describe our (creative) telescoping method, which is the basis of all our extensions and variations. By concrete examples we focus on two main aspects: First, we explain how nested multi-sums can be formulated in the so-called ΠΣ^{∗}-fields. Given these notions we can specify precisely, what type of nested multi-sums Sigma can handle and where Sigma might run into problems. Second, we demonstrate how the telescoping algorithm works.

Finally, in Section 7 we give an overview of all summation problems that Sigma can handle.

4. Concrete examples for telescoping and creative telescoping

We illustrate the basic strategy of our (creative) telescoping method: (1) Rephrase the problem in terms of difference fields, (2) apply the corresponding algorithm in the given difference field, and (3) interpret this result to get a solution for the original problem. All other algorithms, see Section 7, follow the same strategy.

4.1. A telescoping example. Given f(k) = H_{k}, find a solution for (1.1). In order to accomplish this task, Sigma constructs a difference field in which the summation objects and the action of the shift operator S w.r.t. k can be described properly. In our case, take the rational function field 𝔽 := ℚ(k)(h) over the rational numbers ℚ where h represents H_{k}. Moreover, take the automorphism σ : 𝔽 → 𝔽 defined by σ(c) = c for all c ∈ ℚ and

\begin{aligned}σ(k) & = k + 1 & ⟷ & S k & = k + 1, \\ σ(h) & = h + \frac{1}{k + 1} & ⟷ & S H_{k} & = H_{k} + \frac{1}{k + 1}. \end{aligned}  (2.1)

By construction σ reflects the action of the shift operator on k and H_{k} in the field 𝔽. In a nutshell, our summation objects are represented in difference fields.

Definition. A difference field (𝔽, σ) is a field 𝔽 plus a field automorphism σ : 𝔽 → 𝔽; in this article we restrict to fields of characteristic zero.

Next, we solve (1.1) in our difference field (𝔽, σ), i.e., we look for a g ∈ 𝔽 such that

σ(g) − g = h.  (2.2)

In Section 6.1 we show how Sigma computes the solution

g = (h − 1)k ∈ 𝔽.  (2.3)

Rephrasing the result in terms of our summation objects, we obtain the solution g(k) = (H_{k} − 1)k for (1.1). By telescoping we get ∑_{k=1}^{n} H_{k} = (H_{n+1} − 1)(n + 1).

Summarizing, telescoping (1.1) can be formulated in a difference field (𝔽, σ) as follows:

Telescoping in difference fields

Given f ∈ 𝔽 where f represents f(k); find g ∈ 𝔽 with

σ(g) − g = f.  (2.4)

4.2. A creative telescoping example. In order to get the solution (1.3), we proceed as follows. First, we construct a difference field (𝔽, σ) in which we can describe f(n, k) = \binom{n}{k}H_{k} with the shift S in k. Namely, take the rational function field 𝔽 := ℚ(n)(k)(h)(b) with the automorphism σ : 𝔽 → 𝔽 given by σ(c) = c for all c ∈ ℚ(n) and

\begin{aligned}σ(k) & = k + 1 & ⟷ & S k & = k + 1, \\ σ(h) & = h + \frac{1}{k + 1} & ⟷ & S H_{k} & = H_{k} + \frac{1}{k + 1}, \\ σ(b) & = \frac{n − k}{k + 1} b & ⟷ & S \binom{n}{k} & = \frac{n − k}{k + 1} \binom{n}{k}. \end{aligned}  (2.5)

Now observe that one can also represent f(n + i, k) for i ≥ 0 in 𝔽, e.g.,

\begin{aligned}f(n, k) = H_{k} \binom{n}{k} & ⟷ h b =: f_{0} \\ f(n + 1, k) = \frac{(n + 1) H_{k} \binom{n}{k}}{n + 1 − k} & ⟷ \frac{(n + 1) h b}{n + 1 − k} =: f_{1} \\ f(n + 2, k) = \frac{(n + 1) (n + 2) H_{k} \binom{n}{k}}{(n + 1 − k) (n + 2 − k)} & ⟷ \frac{(n + 1) (n + 2) h b}{(n + 1 − k) (n + 2 − k)} =: f_{2}. \end{aligned} (2.6)

Then the creative telescoping problem (1.2) with δ = 2 can be formulated in (𝔽, σ) as follows: FIND c_{0}, c_{1}, c_{2} ∈ ℚ(n) and g ∈ 𝔽 such that

σ(g) − g = c_{0}f_{0} + c_{1}f_{1} + c_{2}f_{2}.  (2.7)

Sigma computes the solution

c_{0} := 4 (1 + n), c_{1} := −2 (3 + 2 n), c_{2} := 2 + n, g := \frac{(1+n) (−2+k−n+(2 k−2 k^{2}+k n)h) b)}{(1−k+n) (2−k+n)};

representing this result in terms of \binom{n}{k} and H_{k} gives (1.3).

Observe that the shift operator S and automorphism σ keep n fixed, i.e., n is a constant. More generally, let (𝔽, σ) be a difference field. Then the set of constants is defined by

const_{σ}𝔽 := {c ∈ 𝔽|σ(c) = c}.

Since const_{σ}𝔽 is a sub-field of 𝔽, we call const_{σ}𝔽 also the constant field of (𝔽, σ); it is easy to see that the rational numbers ℚ must be contained in const_{σ}𝔽. E.g., for the difference field (ℚ(k)(h), σ) the constant field is ℚ and for (ℚ(n)(k)(h), σ) the constant field is ℚ(n); see Section 5.

For a given difference field (𝔽, σ) with constant field 𝕂 creative telescoping reads as follows:

Creative telescoping in difference fields

Given f_{0}, …, f_{δ} ∈ 𝔽 where f_{i} corresponds to f(n + i, k); find c_{i} ∈ 𝕂, not all zero, and g ∈ 𝔽 with

σ(g) − g = c_{0} f_{0} + ⋯ + c_{δ} f_{δ}.  (2.8)

Subsequently we explain in details how the ideas from above can be handled by the computer: we show how one can represent f(k) (resp. f(n + i, k) for i = 0 … δ) in a difference field and we present an algorithm for (creative) telescoping. As it turns out, the construction of an appropriate difference field and telescoping are closely related.

5. A constructive theory of ΠΣ^{∗}-extensions

In Sigma the automatic construction of difference fields relies on the fact that indefinite nested sums and products can be adjoined step by step in form of difference field extensions.

Definition. Let (𝔼, σ) and (𝔽, σ′) be difference fields. We call (𝔼, σ) a difference field extension of (𝔽, σ′) if 𝔼 is a field extension of 𝔽 and σ′(g) = σ(g) for all g ∈ 𝔽. Since σ′ and σ agree on 𝔽, we do not distinguish between the two automorphisms σ and σ′ anymore.

We can construct a difference field for f(k) (resp. for f(n, k)) where we follow the rule that inner objects of a sum or product come first. E.g., for f(n, k) = H_{k}\binom{n}{k} we start with

(ℚ(n), σ) where const_{σ}ℚ(n) = ℚ(n) and we adjoin the objects in the order →^{(1)} k →^{(2)} H_{k} →^{(3)} \binom{n}{k} as follows.

(1) k with Sk = k + 1: Take the rational function field ℚ(n)(k) and extend σ to the automorphism σ : ℚ(n)(k) → ℚ(n)(k) with σ(k) = k + 1. We obtain the difference field (ℚ(n)(k), σ).

(2) H_{k} with SH_{k} = H_{k} + \frac{1}{k+1}: Take the rational function field ℚ(n)(k)(h) and extend the automorphism from σ : ℚ(n)(k) → ℚ(n)(k) to σ : ℚ(n)(k)(h) → ℚ(n)(k)(h) with σ(h) = h + \frac{1}{k+1}. We get the difference field (ℚ(n)(k)(h), σ).

(3) \binom{n}{k} with S\binom{n}{k} = \frac{n−k}{k+1} \binom{n}{k}: Take the rational function field ℚ(n)(k)(h)(b) and construct the difference field extension (ℚ(n)(k)(h)(b), σ) of (ℚ(n)(k)(h), σ) with σ(b) = \frac{n−k}{k+1} b. We arrive at the difference field (ℚ(n)(k)(h)(b), σ) with (2.5) and represent f(n, k) with h b. Following this construction, one obtains difference fields in which one can formulate

rational expressions in terms of nested sums and products. There is only one drawback:

Life gets difficult when one wishes to tackle telescoping and creative telescoping; see Technical remark 2. Thus, we refine our construction with Karr’s ΠΣ^{∗}-theory [Kar81, Kar85]: during the construction one is not allowed to extend the constant field.

Definition. A difference field extension (𝔽(t), σ) of (𝔽, σ) is called a Σ^{∗}-extension (resp. a Π -extension) if t is transcendental over 𝔽, σ(t) = t + a (resp. σ(t) = a t) for some a ∈ 𝔽^{∗}, and const_{σ}𝔽(t) = const_{σ}𝔽. A Π Σ^{∗}-extension is a Π- or a Σ^{∗}-extension. A Π Σ^{∗}-field (𝕂(t_{1}) … (t_{e}), σ) over 𝕂 is a tower of ΠΣ^{∗}-extensions starting with the constant field 𝕂.

As it turns out, the extensions in (ℚ(n)(k)(h)(b), σ) are all ΠΣ^{∗}-extensions (see below). This means that const_{σ}ℚ(n)(k)(h)(b) = ℚ(n). In particular, (ℚ(n)(k)(h)(b), σ) is a ΠΣ^{∗}- field over ℚ(n). Similarly, (ℚ(k)(h)(b), σ) with (2.5) is a ΠΣ^{∗}-field over ℚ. Moreover, all the examples in Part 1, except in Sections 3.2 and 3.4, can be formalized in ΠΣ^{∗}-fields.

Technical remark 1. Karr defines Σ-extensions [Kar81, Kar85] which are of the form σ(t) = α t + β and which must satisfy rather technical side conditions. Restricting to the sum case (α = 1), we obtain exactly the class of Σ^{∗}-extensions. □

Subsequently, we present methods that construct a ΠΣ^{∗}-field for a given expression in terms of nested sums and products.

5.1. Σ^{∗}-extensions. The following beautiful result is a direct consequence of Karr’s theory [Kar81]. For an explicit proof see [Sch01, Cor. 2.2.4].

Σ-Theorem. Let (𝔽(t), σ) be a difference field extension of (𝔽, σ) with σ(t) = t + f for some f ∈ 𝔽; note that t might be algebraic or transcendental over 𝔽. Then this extension is a Σ^{∗}-extension if and only if there is no g ∈ 𝔽 with σ(g) = g + f.

With this result we can easily show that the constructed difference field (ℚ(k)(h), σ) with (2.1) is a ΠΣ^{∗}-field. Consider the difference field extension (ℚ(k), σ) of (ℚ, σ). Since ℚ is the constant field, there is no g ∈ ℚ with σ(g) = g + 1. Thus, (ℚ(k), σ) is a Σ^{∗}- extension of (ℚ, σ). Similarly, by using e.g. Gosper’s algorithm or the Sigma-package, one can check that there is no g ∈ ℚ(k) with σ(g) = g + \frac{1}{k+1}. Hence, (ℚ(k)(h), σ) is a Σ^{∗}- extension of (ℚ(k)(h), σ). Analogously, one can check that (ℚ(n)(k)(h), σ) with (2.1) is a tower of Σ^{∗}-extensions.

Here is an example where the construction from above fails: We cannot construct a Σ^{∗}- extension (ℚ(k)(h)(t), σ) of (ℚ(k)(h), σ) with the shift-behavior σ(t) = t + h. This follows by the Σ-Theorem and the fact that there is the solution (2.3) for σ(g) − g = h. Luckily, there is no need to adjoin a new variable t: if we need t with σ(t) = t + h, then take t := g.

In general, suppose we are given a sum T (k) = ∑_{i=0}^{k} F (i) and a difference field (𝔽, σ) in which we can represent F (k) with f ∈ 𝔽. Due to the telescoping algorithm for ΠΣ^{∗}-fields (see Section 6), only the following two situations can occur during the construction of a ΠΣ^{∗}-field in which one can represent the sum T (k) with

T (k + 1) = T (k) + F (k + 1) :

(1) We find a t ∈ 𝔽 with

σ(t) = t + σ(f).  (2.9)

Then we can represent T (k) by t in 𝔽 with the desired shift-behavior (2.9).

(2) We show that there is no t ∈ 𝔽 with (2.9). Then by the Σ-Theorem we can adjoin the sum T (k) in form of the Σ^{∗}-extension (𝔽(t), σ) of (𝔽, σ) with (2.9).

SUMMARY: We can always represent a rational expression of indefinite nested sums in form of a tower of Σ^{∗}-extensions. Together with the telescoping algorithm for ΠΣ^{∗}-fields, see Section 6, the whole construction mechanism turns out to be algorithmic.

5.2. Π-extensions. For Π-extensions Karr provides the following result [Kar81, Thm. 2].

Π-Theorem. Let (𝔽(t), σ) be a difference field extension of (𝔽, σ) with t ∉ 𝔽 and σ(t) = a t for some a ∈ 𝔽^{∗}; note that t might be algebraic or transcendental over 𝔽. Then this extension is a Π-extension if and only if there are no r > 0 and g ∈ 𝔽^{∗} with

σ(g) = a^{r}g.  (2.10)

Motivated by this result, Karr developed an algorithm [Kar81] which solves the following problem: Given a ΠΣ^{∗}-field (𝔽, σ), decide if there are r > 0 and g ∈ 𝔽^{∗} with (2.10). Summarizing, one can check algorithmically if an extension over a given ΠΣ^{∗}-field is a Π-extension.

E.g., consider the difference field (ℚ(n)(k)(h)(b), σ) from the beginning of Section 5. Since (ℚ(n)(k)(h), σ) is a ΠΣ^{∗}-field, we can apply Karr’s algorithm and show that there is no g ∈ ℚ(n)(k)(h) and r > 0 with

σ(g) = (\frac{n − k}{k + 1} )^{r} g.

Thus, by the Π-Theorem it follows that the difference field extension (ℚ(n)(k)(h)(b), σ) of (ℚ(n)(k)(h), σ) is a Π-extension. Summarizing, we can represent our summand f(n, k) = \binom{n}{k}H_{k} by b h in the ΠΣ^{∗}-field (ℚ(n)(k)(h)(b), σ) in a completely automatic fashion.

Finally, we compute the shift relations in (2.6) as follows. From (2.5) we get

S\binom{n + 1}{k} = \frac{n + 1 − k}{k + 1} \binom{n + 1}{k}.

With the algorithm from Section 6 we compute for

σ(g) = \frac{n + 1 − k}{k + 1} g.  (2.11)

the general solution g = c \frac{(n+1)}{n+1−k}b with c ∈ ℚ(n). Hence, we get \binom{n+1}{k} = c \frac{(n+1)}{n+1−k} \binom{n}{k}; checking the initial value k = 0, shows that c = 1. Therefore we can represent \binom{n+1}{k} by g, and we find \frac{(n+1)}{n+1−k} \binom{n}{k}H_{k} for f(n + 1, k). Similarly, we can proceed with f(n + 2, k). Summarizing, we derive the ΠΣ^{∗}-field (ℚ(n)(k)(h)(b), σ) with the representations (2.6) completely algorithmically.

SUMMARY: Express a product in a given ΠΣ^{∗}-field, or, if this is not possible, try to adjoin it in the form of a Π-extension. But, be careful: This construction might fail!

E.g. we cannot express (−1)^{k} with a Π-extension: for g = 1 and r = 2 we have σ(g) = (−1)^{r}g. Hence, by the Π-Theorem there is no Π-extension (ℚ(P), σ) of (ℚ, σ) with σ(P) = −P. Note that such elements can be only expressed in rings, since we have zero-divisors, like

(1 − (−1)^{k})(1 + (−1)^{k}) = (1 − (−1)^{2k}) = 0.

Luckily, a big class of products can be expressed by Π-extensions. E.g., in [Sch05d] we show that any hypergeometric term can be represented by a Π-extension over the constant field; there is only one exceptional case: a hypergeometric term which can be written in the form r(n)α^{n} where r(n) is a rational function and α is a root of unity; see, e.g., identity (2.22).

6. The basic algorithm

We present a simplified version [Sch05e] of Karr’s algorithm [Kar81].

6.1. Telescoping. Our algorithm finds the solution (2.3) for (2.2) in three reduction steps.

1. Denominator bounding: COMPUTE a polynomial d ∈ ℚ(k)[h]^{∗} such that

∀g ∈ ℚ(k)(h) : σ(g) − g = h ⇒ g d ∈ ℚ(k)[h].

Given such a denominator bound d ∈ ℚ(k)[h]^{∗}, one only has to find the “numerator”, i.e., to find g′ ∈ ℚ(k)[h] with σ(\frac{g′}{d}) − \frac{g′}{d} = h. Observe that this is equivalent to finding all solutions g′ of the first order linear difference equation

a_{1}σ(g′) + a_{0}g′ = h  (2.12)

for the given a_{1} = \frac{1}{σ(d)} and a_{0} = −\frac{1}{d}. In our concrete case we compute the denominator bound d = 1; see [Bro00, Sch04a]. Thus (2.12) is nothing else than our original problem (2.2) with g replaced by g′. We proceed with the second step.

2. Degree bounding: COMPUTE b ≥ 0 with the following property:

∀g ∈ ℚ(k)[h] : σ(g) − g = h ⇒ deg(g) ≤ b.

In our particular case we compute the degree bound b = 2; for more details see [Sch05a].

Technical remark 2. At this point we heavily depend on the fact that (ℚ(k)(h), σ) is a Σ^{∗}-extension of (ℚ(k), σ). We motivate this fact by the following considerations. Suppose we are given a difference field extension (𝔽(t), σ) of (𝔽, σ) where t is transcendental over 𝔽 and σ(t) = α t + β for some α ∈ 𝔽^{∗} and β ∈ 𝔽; let f ∈ 𝔽(t). Then the following holds: If our extension is not a ΠΣ^{∗}-extension and if there is a solution g_{0} ∈ 𝔽(t) for (2.4), then there is no denominator bound and degree bound for (2.4). This can be seen as follows. Since our extension is not a ΠΣ^{∗}-extension, we can take a v ∈ const_{σ}𝔽(t) \ const_{σ}𝔽, i.e., v ∈ 𝔽(t) \ 𝔽 with σ(v) = v; observe that σ(v^{r}) = σ(v)^{r} = v^{r} for any r ≥ 0. Hence,

σ(g_{0} + v^{r}) − (g_{0} + v^{r}) = σ(g_{0}) − g_{0} = f.

Therefore we can increase in g_{0} + v^{r} either the degree of the numerator or the degree of the denominator (or both) by increasing r. Thus, either the denominator bound or the degree bound cannot exist. Similarly, by taking 1/v we can guarantee that neither bound exists. Conversely, if (𝔽(t), σ) is a ΠΣ^{∗}-extension of (𝔽, σ), such bounds exist. Furthermore, they can be computed, if (𝔽, σ) is a ΠΣ^{∗}-field; see [Kar81, Bro00, Sch04a, Sch05a]. □

3. Polynomial solutions: FIND all the coefficients in

g = g_{2} h^{2} + g_{1} h + g_{0} ∈ ℚ(k)[h]  (2.13)

such that (2.2) holds. Graphically the computations can be illustrated as follows:

[figure: [σ(g_{2})(h + \frac{1}{k+1})^{2} + σ(g_{1}h + g_{0})] − [g_{2} h^{2} + g_{1}h + g_{0}] = h; coeff. comp.; σ(g_{2}) − g_{2} = 0; g_{2} = c ∈ ℚ; σ(g_{1} h + g_{0}) − (g_{1} h + g_{0}) = h − c[\frac{2h(k+1)+1}{(k+1)^{2}}]; coeff. comp.; σ(g_{1}) − g_{1} = 1 − c \frac{2}{k+1}; c = 0, \begin{matrix}g_{1} = k + d \\ d ∈ ℚ \end{matrix}; d = 0, \begin{matrix}g_{0} = −k + e \\ e ∈ ℚ \end{matrix} ← σ(g_{0}) − g_{0} = −1 − d\frac{1}{k+1}]

First, we plug the possible solution (2.13) into (2.2) and obtain the equation

[σ(g_{2})(h + \frac{1}{k + 1})^{2} + σ(g_{1}h + g_{0})] − [g_{2} h^{2} + g_{1}h + g_{0}] = h.  (2.14)

By coefficient comparison of the leading terms we get the condition

σ(g_{2}) − g_{2} = 0  (2.15)

for g_{2}. This means that g_{2} = c ∈ ℚ where the constant c is not determined yet. We substitute this partial result into (2.14) and get

σ(g_{1} h + g_{0}) − (g_{1} h + g_{0}) = h − c[\frac{2h(k + 1) + 1}{(k + 1)^{2}}].  (2.16)

Thus we have to find g_{0}, g_{1} ∈ ℚ(k) and c ∈ ℚ with (2.16). We repeat this strategy: By coefficient comparison of the leading terms in (2.16) we get the constraint

σ(g_{1}) − g_{1} = 1 − c \frac{2}{k+1}  (2.17)

for g_{1} ∈ ℚ(k). Solving this problem we find the generic solution c = 0 and g_{1} = k + d for some constant d ∈ ℚ. If we repeat these ideas, see the diagram, we obtain the constraint

σ(g_{0}) − g_{0} = −1 − d\frac{1}{k + 1}  (2.18)

for g_{0} ∈ ℚ(k). This gives d = 0 and g_{0} = −k + e for some constant e ∈ ℚ. Hence, we obtain the general solution g = kh − k + e with e ∈ ℚ; by setting e = 0, we arrive at (2.3).

Since we want to turn this reduction strategy into an algorithm, we need algorithms for the “coefficient problems” (2.15),(2.17), and (2.18) in the difference field (ℚ(k), σ).

More generally, suppose we are given a ΠΣ^{∗}-extension (𝔽(t), σ) of (𝔽, σ). Then solving the telescoping problem (2.4) with our three reduction steps leads to coefficient problems of the following form.

Given f_{1}, …, f_{δ} ∈ 𝔽, a_{0}, a_{1} ∈ 𝔽; find all (c_{1}, …, c_{δ}) ∈ 𝕂^{δ} and h ∈ 𝔽 with

a_{1} σ(h) − a_{0} h = c_{1} f_{1} + ⋯ + c_{δ} f_{δ}.  (2.19)

Technical remark 3. (1) Since 𝕍 := {(c_{1}, …, c_{δ}, g) ∈ 𝕂^{δ} × 𝔽|(2.19) holds} is a vector space over 𝕂 with dimension ≤ δ + 1, problem (2.19) can be solved by finding a basis of 𝕍.

(2) Let (𝔽, σ) be itself a ΠΣ^{∗}-field over 𝕂, i.e., 𝔽 = 𝕂(t_{1}) … (t_{e}) is a tower of ΠΣ^{∗}- extensions. Then in [Sch05e] we work out that problem (2.19) can be solved as above: Compute a denominator bound d ∈ 𝕂(t_{1}) … (t_{e−1})[t_{e}], see [Bro00, Sch04a], then bound the degree of the possible numerators, see [Sch05a], and afterwards extract the coefficients of the numerator by solving problems of the type (2.19) in 𝕂(t_{1}) … (t_{e−1}). Hence we can reduce problem (2.19) in 𝕂(t_{1}) … (t_{e}) to several problems of (2.19) in the smaller field 𝕂(t_{1}) … (t_{e−1}). By recursion, we end up at the base case (2.19) with 𝔽 = 𝕂 which can be solved with linear algebra. □

SUMMARY: There is an algorithm that solves problem (2.19) for a given ΠΣ^{∗}-field (𝔽, σ).

6.2. Creative telescoping. The key observation is that problem (2.19) covers creative telescoping (2.8). Hence, applying our algorithm from Section 6.1 we get a creative telescoping algorithm for ΠΣ^{∗}-fields.

6.3. Solving first order linear difference equations. Solving first order linear difference equations is contained in problem (2.19). E.g., we can compute with our algorithm the solution g = \frac{(n+1)}{n+1−k}b for the homogeneous equation (2.11). More generally, our algorithm can solve recurrences of the form (1.4) with order r = 1 where a_{0}(n), a_{1}(n) and f(n) can be indefinite nested sum-product expressions.

7. Generalizations

By variations and generalizations of our telescoping algorithm we can solve the following summation problems in Sigma.

7.1. Refined (creative) telescoping.

Refined Telescoping

Given f(k); find g(k) and f′(k) such that

f(k) = g(k + 1) − g(k) + f′(k)  (2.20)

where f′(k) is simpler than f(k).

Given such a solution, one finds (under the assumption that (2.20) holds for all 0 ≤ k ≤ n) the identity

∑_{k=0}^{n} f(k) = g(n + 1) − g(0) + ∑_{k=0}^{n} f′(k).

Subsequently, we suppose that we are given a ΠΣ^{∗}-field (𝔽(t), σ) where we can represent f(k) by f ∈ 𝔽(t); see Section 5. Then we can handle the following variations of “simpler”.

7.1.1. Degree optimal w.r.t the top extension.

Given f ∈ 𝔽(t); find (f′, g) ∈ 𝔽(t)^{2} such that

σ(g) − g + f′ = f  (2.21)

where in f′ the degrees of the numerator and denominator polynomials are minimal.

Reinterpreting f′ and g as sequences f′(k) and g(k) solves (2.20) where the sum or product t occurs with optimal degree in the numerator and denominator of f′(k).

The rational case (𝔽(t) = 𝕂(t) with σ(t) = t + 1) has been considered in [Abr75, Pau95, PS95c]. In Sigma the general case of ΠΣ^{∗}-fields can be handled; see [Sch07]. E.g., with the option SimplifyByExt→DepthNumberDegree our solver finds the following simplifications:

Séminaire Lotharingien de Combinatoire 56 (2007), Article B56b 17

In[30]:= SigmaReduce[∑_{k=2}^{n} \frac{2 − kH_{k} + H_{k}^{4} − kH_{k}^{5}}{H_{k} − kH_{k}^{2}} , SimplifyByExt → DepthNumberDegree]

Out[30]= ∑_{k=2}^{n} \frac{k^{2} + 2H_{k}}{k^{2}H_{k}} + (n + 1)H_{n}^{3} − (2n + 1)(\frac{3}{2}H_{n}^{2} − 3H_{n}) − \frac{3}{2}(4n + 1) + \frac{1}{H_{n}}

In[31]:= SigmaReduce[∑_{k=1}^{n} H_{k}^{4}, SimplifyByExt → DepthNumberDegree]

Out[31]= ∑_{k=1}^{n} \frac{2H_{k}k − 2k − 1}{k^{3}} + (n + 1)H_{n}^{4} − 2(2n + 1)H_{n}^{3} + 6(2n + 1)H_{n}^{2} − 12(2n + 1)H_{n} + 24n

Note that the found sums can be simplified further to

∑_{k=2}^{n} \frac{k^{2} + 2H_{k}}{k^{2}H_{k}} = ∑_{k=1}^{n} \frac{1}{H_{k}} + 2H_{n}^{(2)} − 3, ∑_{k=1}^{n} \frac{2H_{k}k − 2k − 1}{k^{3}} = −H_{n}^{(3)} − 2H_{n}^{(2)} + 2 ∑_{k=1}^{n} \frac{H_{k}}{k^{2}} .

Remark: The analogous problem for Π-extensions has been considered in [Sch05d] which generalizes hypergeometric results from [AP02, ALP03]:

Given f ∈ 𝔽(t); find (f′, g) ∈ 𝔽(t)^{2} with f = \frac{σ(g)}{g}f′ where in f′ the degrees of the numerator and denominator polynomials are minimal.

Reinterpreting f′ and g as sequences f′(k) and g(k), we get, with some mild extra conditions, the product representation ∏_{k=1}^{n} f(k) = \frac{g(n+1)}{g(1)} ∏_{k=1}^{n} f′(k). Examples are

∏_{k=1}^{n} \frac{(−k−1)(k+7)}{(k+4)^{2}} = \frac{4}{35} \frac{(n+5)(n+6)(n+7)}{(n+2)(n+3)(n+4)} (−1)^{n},  (2.22)

∏_{k=1}^{n} \frac{(k+3)(H_{k}(k+1)+1)^{2}(H_{k}(k+2)(k+1)+2k+3)}{(k+1)^{2}H_{k}(H_{k}(k+3)(k+2)(k+1)+3(k+4)k+11)} = \frac{11}{6} \frac{(n+3)(n+2)(H_{n}(n+1)+1)^{2}}{(n+1)(H_{n}(n+3)(n+2)(n+1)+3(n+4)n+11)} ∏_{k=1}^{n} H_{k}.

7.1.2. Simpler w.r.t. the depth.

Given f ∈ 𝔽; find (f′, g) ∈ 𝔽^{2} with (2.21) where the nested depth of the sums and products in f′ is minimal.

Reinterpreting f′ and g as sequences f′(k) and g(k) solves (2.20) where only those sums and products of f(k) occur in f′(k) which make the depth of f′(k) optimal. This mechanism [Sch04c] is activated by setting the option SimplifyByExt→Depth:

In[32]:= SigmaReduce[∑_{k=1}^{n} H_{k}^{2}H_{k}^{(2)}, SimplifyByExt → Depth]

Out[32]= \frac{1}{3} ∑_{i=1}^{n} \frac{1}{i^{3}} − \frac{1}{3}H_{n}^{3} + ((n + 1)H_{n}^{(2)} + 1) H_{n}^{2} + (2n + 1) (1 − H_{n}) H_{n}^{(2)} − 2H_{n}

Note that in Out[32] we find the sum extension H_{n}^{(3)} = ∑_{i=1}^{n} \frac{1}{i^{3}} in order to represent ∑_{k=1}^{n} H_{k}^{2}H_{k}^{(2)} in terms of sums with depth 1. For further examples see In[18], In[29], or

∑_{k=1}^{n} H_{k}^{3} = (n + 1)H_{n}^{3} − \frac{3}{2}(2n + 1)H_{n}^{2} + 3(2n + 1)H_{n} − 6n + \frac{1}{2}H_{n}^{(2)},

∑_{k=0}^{n} (∑_{i=0}^{k} \binom{m}{i})^{2} = (m − n)\binom{m}{n} ∑_{i=0}^{n} \binom{m}{i} − \frac{m − 2n − 2}{2} (∑_{i=0}^{n} \binom{m}{i})^{2} − \frac{m}{2}∑_{i=0}^{n} \binom{m}{i}^{2}.

Note that the last expression can be simplified further to \frac{n4^{n}}{2} + 4^{n} − \frac{n}{2} \binom{2n}{n} if m = n; see [Hir96, AP99].

More generally, in [Sch05b] we solve the following problem:

Given f ∈ 𝔽; find a tower of ΠΣ^{∗}-extensions^{a} (𝔽(t_{1}) … (t_{e}), σ) over (𝔽, σ) and (f′, g) ∈ 𝔽(t_{1}) … (t_{e})^{2} with (2.21) where f′ has minimal depth.

aThe nested depths of the extensions t_{i} are smaller than the depth of any element in f.

Reinterpreting f′ and g as sequences f′(k) and g(k) solves (2.20) where the nested depth of the sums and products occurring in f′(k) are optimal. With Sigma we find

∑_{k=1}^{n} ( ∑_{j=1}^{k} ^{\frac{H_{j}^{(2)}}{j^{3}}} )^{2} = −(H_{n}^{(2)}^{2} + H_{n}^{(4)}) ∑_{j=1}^{n} ^{\frac{H_{j}^{(2)}}{j^{3}}} + (n + 1)( ∑_{j=1}^{n} ^{\frac{H_{j}^{(2)}}{j^{3}}} )^{2} + ∑_{j=1}^{n} ^{\frac{H_{j}^{(2)}((jH_{j}^{(2)})^{2}−H_{j}^{(2)}+j^{2}H_{j}^{(4)})}{j^{5}}},

∑_{k=1}^{n} \frac{1}{k^{3}} ∑_{j=1}^{k} \frac{H_{j}}{j^{2}} = H_{n}^{(3)} ∑_{j=1}^{n} \frac{H_{j}}{j^{2}} − ∑_{j=1}^{n} \frac{H_{j}(j^{3}H_{j}^{(3)} − 1)}{j^{5}}.  (2.23)

E.g., for (2.23) we compute the extensions H_{n}^{(3)} and ∑_{j=1}^{n} H_{j}(j^{3}H_{j}^{(3)} − 1)/j^{5} in order to reduce the 3-nested sum on the left-hand side to an expression with at most 2-nested sums.

7.1.3. Creative telescoping. The refined telescoping algorithms from [Sch04c, Sch05b] can be carried over to creative telescoping:

Given δ ≥ 0 and f_{i} ∈ 𝔽; find (f′, g) ∈ 𝔽^{2} and c_{i} ∈ const_{σ}𝔽 such that

σ(g) − g + f′ = c_{0}f_{0} + ⋯ + c_{δ}f_{δ}

where f′ is nicer than the f_{i}.

E.g., by using the option SimplifyByExt->DepthNumber in GenerateRecurrence one

looks for an f′ where the number of used objects are smaller than the objects occurring in the f_{i}. Typical examples can be found in In[15] or in [Sch02, PS03, DPSW06b].

7.2. Solving linear difference equations of higher order. The recurrence solver in Sigma works as follows. Given a recurrence (1.4), Sigma represents the coefficients a_{i}(n) and the inhomogeneous part f(n) in a ΠΣ^{∗}-field (𝔽, σ) with a_{i}, f ∈ 𝔽; see Section 5. We call (𝔽, σ) also the coefficient field of the given recurrence. Then there are various options how to continue.

7.2.1. Solutions in the coefficient field.

Given f, a_{0}, …, a_{r} ∈ 𝔽; find all g ∈ 𝔽 with

a_{r} σ^{r}(g) + ⋯ + a_{0} g = f.  (2.24)

Then rephrasing the elements g ∈ 𝔽 to sequences g(k) produces solutions for the recurrence (1.4). There are such solvers for the rational case and the q-rational case, see [Abr89a, Abr89b, Pet92, ABP95, Abr95, PWZ96, Hoe98, APP98]. Also in Sigma such an efficient solver is available: a typical example is the first homogeneous solution in the result Out[28].

More generally, Sigma contains methods for ΠΣ^{∗}-fields [Sch05a]. With our solver we compute, e.g., the particular solution in the result Out[25]. Similarly, given

In[33]:= rec = −n(n + 1)H_{n}((n + 1)H_{n} + 1)F[n] + n(n + 1)(2H_{n} + 1)(n + 1)H_{n} + 1)F[n + 1] −

n(n + 1)H_{n} (n + (n + 1)H_{n} + 2) F[n + 2] == H_{n}((n + 1)H_{n} + 1);

we can compute a particular solution:

In[34]:= SolveRecurrence[rec, F[n], Extension → None]

Out[34]= \{\{1, \frac{nH_{n} − 1}{n}\}\}

Technical remark 4. In [Sch05a] we generalize the telescoping algorithm from Section 6 to solve problem (2.24) for a ΠΣ^{∗}-field (𝔽(t), σ). Here the following remarks are in place:

Reduction 1: Denominator bounds can be computed for Σ^{∗}-extensions. For Π-extensions

it can be determined up to a factor of the form t^{l} with l ∈ ℕ_{0}; see [Bro00, Sch04a]. Reduction 2: Degree bounds can be computed for several special cases [Sch05a]. In [Sch01]

a method has been developed that can compute degree bounds for Σ^{∗}-extensions. So far, I did not find a proof for termination.

Reduction 3: The coefficient problems are of the following type.

Given a_{0}, …, a_{r} ∈ 𝔽 and f_{1}, …, f_{δ} ∈ 𝔽; find all g ∈ 𝔽 and (c_{1}, …, c_{δ}) ∈ 𝕂^{δ} such that

a_{r} σ^{r}(g) + ⋯ + a_{0} g = c_{1} f_{1} + ⋯ + c_{δ} f_{δ}.  (2.25)

Note that 𝕍 := {(c_{1}, …, c_{δ}, g) ∈ 𝕂^{δ} × 𝔽|(2.25) holds} is a vector space over 𝕂 with dimension ≤ δ + r. Notice that problem (2.19) is a special case of (2.25) with r = 1. E.g., in Out[9],Out[11],Out[17], Out[28] we output bases of such vector spaces with δ = 1, r = 2.

In order to solve problem (2.25), we apply the three reduction steps recursively as in the telescoping algorithm. Here we emphasize the following result: Although there are open subproblems in the reduction steps 1 and 2, it has been shown in [Sch05e] that there is a recursive enumeration procedure that eventually outputs all solutions for a given equation (2.25). Further investigations are going on [Bro05, AP06] to overcome these open problems. □

SUMMARY: The methods for the “master problem” (2.25) are the algorithmic heart of Sigma. All the other problems treated here, such as (creative) telescoping and solving recurrences, can be reduced to it.

As one can see in Out[34], we missed the homogeneous solutions of the recurrence. The problem is that we searched for solutions only in the coefficient field “ℚ(n, H_{n})” given by In[33]. In order to extend the search space, the following possibilities are available.

7.2.2. Manual extensions. The coefficient field can be extended manually by using the option Tower → {ext_{1}, …, ext_{e}}. This feature might be useful, if one has additional insight, i.e., one expects that certain sums or products should occur in the solution.

7.2.3. Automatic extensions. Sigma finds certain type of solutions for (1.4), or it outputs that solutions of such type do not exist. We focus on the following problem.

Find solutions by extensions

Given (1.4); find all solutions of the form

h(n) ∑_{k_{1}=0}^{n} b_{1}(k_{1}) ∑_{k_{2}=0}^{k_{2}} b_{2}(k_{2}) ⋯ ∑_{k_{s}=0}^{k_{s−1}} b_{s}(k_{s}).  (2.26)

• Sum solutions. Sum solutions are of the form (2.26) where the b_{i}(k_{i}) and h(n) can be represented in the given coefficient field (𝔽, σ). In other words, the b_{i}(k_{i}) and h(n) are expressions in terms of the objects given in a_{i}(n) and f(n); for examples see Out[17] and Out[28].

Technical remark 5. (1) Any solution, that can be represented in a tower of Σ^{∗}- extensions, can be represented by a sum solution of the form (2.26); see [Sch01, Thm. 4.5.4].

(2) If there exists a sum solution (2.26), then the expression h(n) must be a solution of the homogeneous version of (1.4); see [Sch01, Thm. 4.5.1].

(3) Sum-solutions are obtained by factorizing the linear difference equation as much as possible into first order linear right factors over the given difference field/ring. Then each factor corresponds basically to one indefinite summation quantifier; see [AP94, Sch01]. If one fails to split off such a first order factor, there is still hope to continue. Namely, if one finds a product solution of the remaining difference equation, then this corresponds exactly to one additional factor. Hence, product extensions can lead to a refined factorization, and therefore can produce additional solutions of a given linear difference equation. □

SUMMARY: We find all solutions in terms of indefinite sum expressions by looking for all sum solutions. If there is no homogeneous solution in the coefficient field (𝔽, σ), like in Out[34], there is no sum solution at all. Luckily, product extensions can contribute to finding additional sum solutions, see Out[9], Out[11], or Out[25].

• Product solutions. Given a recurrence with rational coefficients in 𝕂(n), there are algorithms, like [Pet92] or [Hoe99], that can compute the so-called hypergeometric solutions; for q-hypergeometric solutions see [APP98]. Typical examples are the first entry in Out[9] for the hypergeometric case and in Out[11] for the q-hypergeometric case.

A generalized version of the algorithms [Pet92, APP98] is implemented in Sigma for ΠΣ^{∗}-fields. E.g. for our recurrence from In[33] we compute the following product solution:

In[35]:= SolveRecurrence[rec, F[n], Extension → PRODS]

Out[35]= \{\{1, \frac{nH_{n} − 1}{n}\}, \{0, ∏_{i=2}^{n} \frac{−1 + iH_{i}}{−1 + i + iH_{i}} \}\}

Using this product extension we can completely solve the recurrence with sum solutions:

In[36]:= recSol = SolveRecurrence[rec, F[n], Tower → \{∏_{i=2}^{n} \frac{−1+iH_{i}}{−1+i+iH_{i}} \}, Extension → SUMS]

Out[36]= \{\{1, \frac{nH_{n} − 1}{n}\}, \{0, ∏_{i=2}^{n} \frac{−1 + iH_{i}}{−1 + i + iH_{i}} \}, \{0, ∏_{i=2}^{n} \frac{−1 + iH_{i}}{−1 + i + iH_{i}} ∑_{i=1}^{n} \frac{−1 + H_{i}i}{−1 + i + H_{i}i} ∏_{j=2}^{i} \frac{−1 + j + jH_{j}}{−1 + jH_{j}} \}\}

• d’Alembertian solutions. In Out[36] we have computed solutions of the type (2.26) where the b_{i}(k_{i}) and h(n) can be either elements from the coefficient field, or products over such elements from 𝔽. Such type of solutions are also called d’Alembertian solutions [AP94], a subclass of Liouvillian solutions [HS99].

With the option Extension→dAlembert we can compute all such solutions in one stroke. E.g., we solve our recurrence In[33] at once with

In[37]:= SolveRecurrence[rec, F[n], Extension → dAlembert]

Out[37]= \{\{1, \frac{nH_{n} − 1}{n}\}, \{0, ∏_{i=2}^{n} \frac{−1 + iH_{i}}{−1 + i + iH_{i}} \}, \{0, ∏_{i=2}^{n} \frac{−1 + iH_{i}}{−1 + i + iH_{i}} ∑_{i=1}^{n} \frac{−1 + H_{i}i}{−1 + i + H_{i}i} ∏_{j=2}^{i} \frac{−1 + j + jH_{j}}{−1 + jH_{j}} \}\}

Note that Extension→dAlembert is the default option in SolveRecurrence, i.e., in the computation steps In[9], In[11],In[17],In[25],In[28] the d’Alembertian-machinery was activated.

Technical remark 6. Since problem (2.25) occurs as a subproblem in the algorithms for sum-solutions and product-solutions, the open problems listed in Technical remark 4 are relevant here; further investigations are going on [Bro05, AP06]. Note that we can find all such solutions by recursive enumeration. □

SUMMARY: Sigma can find all d’Alembertian solutions. If Sigma fails to find any product solution (including solutions in the coefficient field 𝔽) for a given recurrence, then there does not exist a d’Alembertian solution at all. In this case Sigma’s weapons are exhausted.

7.3. Algebraic extensions. In various summation problems, like in Sections 3.2 and 3.4, the algebraic term (−1)^{n} with ((−1)^{n})^{2} = 1 pops up. As shown in Section 5.2, such an object can be formulated only in rings with zero-divisors. In Sigma our methods for problem (2.25) have been extended for such algebraic extensions; see [Sch01, Section 3.6].

7.4. Unspecified sequences. In joint work with Manuel Kauers [KS06b, KS06a] Sigma has been extended as follows. Our summation objects can be represented in a tower of ΠΣ^{∗}-extensions over a free difference field [Coh65]. More precisely, take the field 𝔽 := 𝕂(…, x_{−1}, x_{0}, x_{1}, …) with infinitely many variables x_{i} and define the field automorphism σ : 𝔽 → 𝔽 by σ(c) = c for all c ∈ 𝕂 and σ(x_{i}) = x_{i+1} for i ∈ ℤ. Then given this difference field (𝔽, σ) we extend it with a tower of ΠΣ^{∗}-extensions, say (𝔽(t_{1}) … (t_{e}), σ).

For such a difference field we managed to carry over all the summation algorithms mentioned earlier. In order to apply this machinery, we have to load in

In[38]:= << Free.m

Free.m - Solver for PLDEs over the free difference field by Manuel Kauers \begin{matrix}c \\ ◯ \end{matrix} RISC-Linz

As carried out in detail in [KS06b] we can find the identity [KS06b, Equ. (5)]

∑_{k=1}^{n} k^{2} ∑_{i=1}^{k} X_{i} = \frac{1}{6}(n(n + 1)(2n + 1) ∑_{k=1}^{n} X_{k} − ∑_{k=1}^{n} kX_{k} + 3 ∑_{k=1}^{n} k^{2}X_{k} − 2 ∑_{k=1}^{n} k^{3}X_{k})

where X_{i} stand for a generic/unspecified sequence. More precisely, we compute the righthand side by simply executing our telescoping-solver:

In[39]:= SigmaReduce[∑_{k=1}^{n} k^{2} ∑_{i=1}^{k} X[i], SimpleSumRepresentation → True]

Out[39]= \frac{1}{6}(n(n + 1)(2n + 1) ∑_{k=1}^{n}X[k] − ∑_{k=1}^{n}kX[k] + 3 ∑_{k=1}^{n}k^{2}X[k] − 2 ∑_{k=1}^{n}k^{3}X[k])

As observed in [KS06b] we can now specialize this identity and get, for instance with X_{k} = \frac{1}{n+k} and H_{n+k} = H_{n} + ∑_{i=1}^{k} \frac{1}{n+i} the identity [GKP94, Bonus problem 6.69]:

∑_{k=1}^{n} k^{2}H_{n+k} = \frac{1}{3}n(n + \frac{1}{2})(n + 1)(2H_{2n} − H_{n}) − \frac{1}{36}(10n^{2} + 9n − 1).

Similarly, we find

∑_{k=1}^{n} a^{k} ∑_{j=1}^{k} X_{j} = \frac{1}{a − 1}(a^{n+1} ∑_{k=1}^{n} X_{k} − ∑_{k=1}^{n} a^{k}X_{k}), a ≠ 1

which generalizes the identity [KS06b, Equ. (7)]. With X_{j} := \frac{1}{j} we get

∑_{k=1}^{n} a^{k}H_{k} = \frac{1}{a − 1}[a^{n+1}H_{n} − ∑_{k=1}^{n} \frac{a^{k}}{k} ]

and with X_{j} = \binom{m}{j−1}, a = −1, and n := m + 1 we rediscover [Zha99]:

∑_{k=0}^{m}(−1)^{k+1} ∑_{j=0}^{k} \binom{m}{j} = \frac{1}{2}(−1)^{m+1}2^{m}.

Observe that for a = 1 we derive a different version, namely [KS06b, Equ. (5)]:

∑_{k=1}^{n} ∑_{j=1}^{k} X_{j} = (n + 1) ∑_{k=1}^{n} X_{k} − ∑_{k=1}^{n} kX_{k}.

For X_{j} := \frac{1}{j^{2}} we find ∑_{k=1}^{n} H_{k}^{(2)} = (n + 1)H_{n}^{(2)} − H_{k} and for X_{j} = \binom{m}{j} we get

∑_{k=0}^{n} ∑_{j=0}^{k} \binom{m}{j} = (n + 1) ∑_{k=0}^{n} \binom{m}{k} − ∑_{k=0}^{n} k\binom{m}{k} = \frac{1}{2}(m − n)\binom{m}{n} + (2n − m + 2) ∑_{i=0}^{n} \binom{m}{i}.

If m = n, we get \frac{n+2}{2}2^{n}; see [Hir96, AP99]. Further identities are given in [KS06b, KS06a].

Part 3. Multi-summation and applications

There are various approaches for multi-summation available. As illustrated in the first two parts, the difference field approach enables one to handle a rather general class of nested multi-sums. Zeilberger’s holonomic systems approach [Zei90] was an important breakthrough for another class of multi-sums. This work forms a common framework for the Sister Celine/WZ-method and the holonomic/∂-finite function approach.

The Sister Celine/WZ-Method: Following Sister Celine Fasenmyer’s PhD–thesis [Fas45] and D. Zeilberger/H. Wilf [WZ92] one computes suitable recurrences for hypergeometric summands by setting up a system of linear equations; the summand-recurrence can then be transformed to a recurrence for the corresponding hypergeometric multi-sum. An efficient machinery has been developed by Wegschaider [Weg97] where ideas of Sister Celine/Wilf/Zeilberger are combined in a non-trivial manner with results of Verbaeten [Ver74] and its simplification presented in [Hor92]. For a q-version see [Rie03]. Related approaches are [CHM06, AZ06].

The holonomic/∂-finite function approach: Pioneering work of the holonomic/∂-finite approach has been done in [CS98]. In particular, in [Chy00] Zeilberger’s algorithm [Zei91] has been generalized to general holonomic and ∂-finite functions. This method treats also multiple sum (and multiple integration) problems.

A new Sigma approach: In [Sch05c] it has been shown that Chyzak’s approach [Chy00] can be substantially simplified, if one attacks multi-sum problems in a slightly restricted way. Moreover, it turns out that one can bring Karr’s ΠΣ^{∗}-world and Chyzak’s approach under one umbrella. This leads to a rather general and surprisingly efficient machinery which has been implemented in Sigma. Subsequently, we shall illustrate our “Sigma approach” by various concrete examples. In addition to the results [Sch05c], we show in Section 9.3 that our method can be also used to compute differential equations for such general multi-sums.

8. The basic idea for telescoping

We consider the following problem from [BPP^{+}06] which arose in joint cooperation with the JKU-Finite Element group. Find a closed form for the hypergeometric multi-sum

S(n) = ∑_{k=1}^{n} \frac{2k + 1}{k + 1} ∑_{j=0}^{k} \frac{(−k)_{j}(k + 1)_{j}(2)_{k}}{j!k!(2)_{j}} (\frac{1 − x}{2})^{j};

here we use the standard Pochhammer symbol (a)_{j} = ∏_{i=1}^{j}(a + i − 1). Subsequently, we denote the inner sum in S(n) with P (k); note that P (k) are the Jacobi-polynomials P_{k}^{(α,β)}(x) for the specific choice (α, β) = (1, −1).

First notice that we cannot apply our (creative) telescoping algorithm presented earlier to the summand \frac{2k+1}{k+1}P (k). To see this, we compute a recurrence for P (k):

In[40]:= recP = GenerateRecurrence[∑_{j=0}^{k} \frac{(−k)_{j}(k + 1)_{j}(2)_{k}}{j!k!(2)_{j}} (\frac{1 − x}{2})^{j}][[1]]/.SUM → P

Out[40]= P[k + 2] == \frac{(2k + 3)x}{k + 2}P[k + 1] − \frac{k}{k + 1}P[k]

Then we show with Sigma that there are no d’Alembertian solutions, i.e., we cannot represent P (k) in terms of a ΠΣ^{∗}-field. In particular, we cannot simplify Out[40] to a first order recurrence. Summarizing, we cannot handle the sum S(n) with the tools presented so far.

Motivated by such examples, we extended Sigma in order to handle also multi-sums where the inner sum is described by recurrences of higher order. E.g., using Out[40] we can simplify

In[41]:= mySum = ∑_{k=1}^{n} \frac{2k + 1}{k + 1}P[k];

with our generalized telescoping solver as follows.

In[42]:= SigmaReduce[mySum, {recP, P[k]}]

Out[42]= \frac{(3x − 2)P[1]}{2(x − 1)} − \frac{P[2]}{x − 1} − \frac{nP[n]}{(n + 1)(x − 1)} + \frac{P[n + 1]}{x − 1}

With P (1) = x + 1, P (2) = \frac{3}{2}x(x + 1) we arrive at the identity [BPP^{+}06, Equ. (16)]

S(n) = −\frac{x + 1}{x − 1} − \frac{nP (n)}{(n + 1)(x − 1)} + \frac{P (n + 1)}{x − 1},  n ≥ 1.  (3.1)

8.1. Verification. The result Out[42] can be produced by summing the telescoping equation

g(k + 1) − g(k) = \frac{2k + 1}{k + 1}P (k)  (3.2)

with the computed solution

g(k) = \frac{1 + k − x − 2kx}{(x − 1)(k + 1)}P (k) + \frac{1}{x − 1}P (k + 1).  (3.3)

Note that the correctness of (3.2) can be verified independently of the computational steps: Represent P (k + 2) in g(k + 1) as a linear combination of P (k) and P (k + 1) by using the relation Out[40]. Then verify (3.2) by polynomial arithmetic.

8.2. The method. We consider the following telescoping problem: Given Out[40], find g(k) = g_{0}(k)P (k) + g_{1}(k)P (k + 1) with unknown coefficients g_{0}(k) and g_{1}(k) such that (3.2) holds. This problem is equivalent to finding g_{0}(k) and g_{1}(k) with

[g_{0}(k+1)P (k+1)+g_{1}(k+1)P (k+2)]−[g_{0}(k)P (k)+g_{1}(k)P (k+1)] = \frac{2k + 1}{k + 1}P (k). (3.4)

Applying the relation Out[40] and collecting terms w.r.t. P (k) and P (k + 1) we get equivalently

P (k)[ − \frac{k}{k + 1}g_{1}(k + 1) − g_{0}(k) − \frac{2k + 1}{k + 1}]

+ P (k + 1)[g_{0}(k + 1) + \frac{(2k + 3)x}{k + 2}g_{1}(k + 1) − g_{1}(k)] = 0. (3.5)

Hence, if g_{0}(k) and g_{1}(k) satisfy

g_{0}(k) = −\frac{k}{k + 1}g_{1}(k + 1) − \frac{2k + 1}{k + 1},  (3.6)
g_{0}(k + 1) + \frac{(2k + 3)x}{k + 2}g_{1}(k + 1) − g_{1}(k) = 0,  (3.7)

then (3.4) and (3.5) hold. Note that the other direction might not hold in general, but in our concrete case it does; see Technical remark 7.1. Finally, by taking the shifted version of (3.6) we can rewrite (3.7) in the form of the linear difference equation

−\frac{k + 1}{k + 2}g_{1}(k + 2) + \frac{(2k + 3)x}{k + 2}g_{1}(k + 1) − g_{1}(k) = \frac{2k + 3}{k + 2}.  (3.8)

Summarizing, if g_{0}(k) and g_{1}(k) satisfy (3.6) and (3.8), g(k) = g_{0}(k)P (k)+g_{1}(k)P (k +1) is a solution of (3.2) and (3.4). Now observe that (3.8) is a linear recurrence in g_{1}(k). Hence we can run Sigma’s recurrence solver and compute the solution g_{1}(k) = \frac{1}{x−1}. Finally, g_{0}(k) = \frac{1+k−x−2kx}{(x−1)(k+1)} is determined by (3.6). We end up with our solution (3.3).

In general, we obtain a method for the following telescoping problem:

Telescoping with a recurrence

Given f(k) := h_{0}(k)P (k) + ⋯ + h_{s}(k)P (k + s) and

P (k + s + 1) = a_{0}(k)P (k) + ⋯ + a_{s}(k)P (k + s);  (3.9)

find a solution for (1.1) which is of the form

g(k) = g_{0}(k)P (k) + ⋯ + g_{s}(k)P (k + s).  (3.10)

Namely, by inserting (3.10) with the unknown coefficients g_{r}(k) in (1.1) and doing coefficient comparison we find, as above, the following coupled system; see [Sch05c, Lemma 1]:

g_{0}(k) = a_{0}(k) g_{s}(k + 1) − h_{0}(k),  (3.11)
g_{r}(k) = g_{r−1}(k + 1) + a_{r}(k) g_{s}(k + 1) − h_{r}(k),  1 ≤ r ≤ s.  (3.12)

This means that any solution g_{0}(k), …, g_{s}(k) of this system produces a solution (3.10) for (1.1). Now the crucial step is that this system can be uncoupled, once and for all, and we obtain, in addition, the following linear difference equation for g_{s}(k) [Sch05c, Lemma 2]:

∑_{j=0}^{s} a_{s−j}(k + j)g_{s}(k + j + 1) − g_{s}(k) = ∑_{j=0}^{s} h_{s−j}(k + j)  (3.13)

Summarizing, we arrive at the following method:

1. FIND a solution g_{s}(k) for (3.13); this is a particular instance of problem (2.24).

2. COMPUTE g_{0}(k) by (3.11); then compute the remaining g_{1}(k), …, g_{s−1}(k) by (3.12).

Technical remark 7. (1) If g(k) = g_{0}P (k) + ⋯ + g_{s}(k)P (k + s) is a solution of (1.1), we cannot guarantee that the g_{r}(k) satisfy the system (3.11), (3.12), (3.13). Hence, our method might fail, although there exists a solution of the form (3.10). This can be only guaranteed, if the recurrence order of (3.9) is minimal. E.g., the recurrence Out[40] has minimal order. This implies that P (k) and P (k + 1) are linearly independent. Thus (3.5) implies (3.6) and (3.7).

(2) We can apply this method for all input terms h_{i}(k) and a_{i}(k) for which one has solvers for (3.13) or equivalently for (2.24). In the Sigma-package the h_{i}(k) and a_{i}(k) can be any expression represented in a ΠΣ^{∗}-field. Furthermore, in [Sch05c] we allow that the recurrence (3.9) might have an inhomogeneous part. □

9. The basic idea for creative telescoping

In [PWZ96] the following identity pops up:

∑_{k=0}^{n} ∑_{j=0}^{n}(−1)^{n+k+j}\binom{n}{k}\binom{n}{j}\binom{n + k}{k}\binom{n + j}{j}\binom{2n − j − k}{n} = ∑_{k=0}^{n} \binom{n}{k}^{4}.  (3.14)

As observed in [Sch05c] Sigma can prove this identity by a slight generalization of the techniques presented in Section 8. First, we compute a recurrence in k for the inner sum P (n, k):

In[43]:= innerSum = ∑_{j=0}^{n}(−1)^{n+k+j}\binom{n}{k}\binom{n}{j}\binom{n + k}{k}\binom{n + j}{j}\binom{2n − j − k}{n};

In[44]:= recK = GenerateRecurrence[innerSum, k][[1]]/.SUM → P

Out[44]= P[k + 2] == \frac{(n−k)^{3}(1+k+n)(2+k+n)}{(1+k)^{2}(2+k)^{2}(k−3n)}P[k] + \frac{(1+k)^{2}(2+k+n)(k+2k^{2}−3n−6kn+3n^{2})}{(1+k)^{2}(2+k)^{2}(k−3n)}P[k + 1]

Besides this we compute a recurrence with one shift in n and the remaining shifts in k:

In[45]:= recKN = GenerateRecurrence[innerSum, k, OneShiftIn → n][[1]]/.SUM → P

Out[45]= P[n + 1, k] == \frac{(1 + k)^{2}(−1 + k − 3n)(6 − 8k + 3k^{2} + 12n − 8kn + 6n^{2})}{(1 − k + n)^{3}(1 + n)^{2}}P[k + 1] +

      \frac{−(1+k+n)(−5k+12k^{2}−10k^{3}+3k^{4}+3n−32kn+42k^{2}n−16k^{3}n+15n^{2}−57kn^{2}+33k^{2}n^{2}+21n^{3}−30kn^{3}+9n^{4})}{(1−k+n)^{3}(1+n)^{2}}P[k]

Remark. Given the summand f(n, k, j), set up the creative telescoping equation

g(n, k, j + 1) − g(n, k, j) = c_{0}(n, k)f(n, k, j) + c_{1}(n, k)f(n, k + 1, j) + c_{2}(n, k)f(n + 1, k, j) (3.15)

and solve the underlying problem (2.8) where f_{0}, f_{1}, f_{2} correspond to f(n, k, j), f(n, k + 1, j), f(n+1, k, j), respectively. Summing the resulting equation (3.15) over k gives Out[45].

Finally, we compute a recurrence for the sum on the left-hand side of (3.14) as follows:

In[46]:= GenerateRecurrence[∑_{k=0}^{n} P[k], n, recK, P[k], recKN]

Out[46]= {−4(1 + n)(3 + 4n)(5 + 4n)SUM[n] − 2(3 + 2n)(7 + 9n + 3n^{2})SUM[1 + n] + (2 + n)^{3}SUM[2 + n] == 0}

We remark that Chyzak’s general holonomic approach (which takes 2300s) and Wegschaider’s implementation (which takes 510s) are much slower on that; we need only 12s on the same machine.

To this end, we can compute with Sigma the same recurrence for the right-hand side of (3.14). Since both sides of (3.14) agree at n = 0, 1, equality follows for all n ≥ 0.

9.1. Verification. The correctness of Out[46] follows by

g(n, k + 1) − g(n, k) = c_{0}(n)P (n, k) + c_{1}(n)P (n + 1, k) + c_{2}(n)P (n + 2, k)  (3.16)

and the proof certificate

\begin{aligned}c_{0}(n) & = 4(n + 1)^{3}(4n + 3)(4n + 5), & c_{2}(n) & = −(n + 1)^{2}(n + 2)^{3}, \\ c_{1}(n) & = 2(n + 1)^{2}(2n + 3) (3n^{2} + 9n + 7) , \end{aligned}  (3.17)

and

g(n, k) = g_{0}(n, k)P (n, k) + g_{1}(n, k)P (n, k + 1)  (3.18)

where the rational functions g_{i}(n, k) in n and k are derived in Example 9.4; the explicit expression can be found in [Sch05c, p. 763]. The correctness of (3.16) can be verified for all 0 ≤ k ≤ n as follows. Rewrite the expression P (n, k + 2) in g(n, k + 1) in terms of P (n, k) and P (n, k + 1) by using the relation Out[44]. Similarly, rewrite the expression P (n + 1, k) in (3.16) in the form of a linear combination of P (n, k) and P (n, k + 1) by using the relation Out[45]. Moreover, express P (n + 2, k) in (3.16) by a linear combination in P (n + 1, k) and P (n + 1, k + 1) which itself can be expressed by a linear combination in P (n, k) and P (n, k + 1) by using the “rewrite rules” Out[44] and Out[45]. Then the correctness of (3.16) follows by polynomial arithmetic. Summing (3.16) over k produces Out[46].

9.2. A method for recurrences. We present the following strategy from [Sch05c] to compute a recurrence for a hypergeometric double sum

S(n) = ∑_{k} ∑_{j} h(n, k, j)

where h(n, k, j) is hypergeometric in n, k and j. We start with the inner sum P (n, k) = ∑_{j} h(n, k, j). If one is lucky, one can compute with Sigma not only a recurrence in k, say

P (n, k + s + 1) = a_{0}(n, k)P (n, k) + ⋯ + a_{s}(n, k)P (n, k + s),  (3.19)

but also a recurrence with one shift in n and the remaining shifts in k, like

P (n + 1, k) = b_{0}(n, k)P (n, k) + ⋯ + b_{s}(n, k)P (n, k + s);  (3.20)

examples are In[44] and In[45]. In [Pau04] an existence theory is presented which closely relates to the situation of Zeilberger’s algorithm. This question is analysed further in [PS04].

Given such a recurrence system (3.19),(3.20), our creative telescoping problem reads as follows.

Creative telescoping with a recurrence system

Given δ ∈ ℕ and f(n, k) := h_{0}(n, k)P (n, k) + ⋯ + h_{s}(n, k)P (n, k + s) plus (3.19) and (3.20); find a solution c_{0}(n), …, c_{δ}(n), not all zero, and g(n, k) for (1.2) where g(n, k) is of the form

g(n, k) = g_{0}(n, k)P (n, k) + ⋯ + g_{s}(n, k)P (n, k + s).  (3.21)

Summing (3.21) over k gives, with some extra conditions, a recurrence of the form (1.4).

Example 9.1. Take f(n, k) = P (n, k) with (3.19) and (3.20) where the a_{i}(k) and b_{i}(k) are given by Out[44] and Out[45], respectively. We look for c_{0}(n), c_{1}(n), c_{2}(n) and (3.18) with (3.16). □

Given the two relations (3.19) and (3.20), the terms P (n + i, k), …, P (n + i, k + s) in f(n + i, k) can be expressed by a linear combination in P (n, k), …, P (n, k + s). E.g., if s = 1 and i = 2,

P (n + 2, k) =^{(3.20)} b_{0}(n + 1, k)P (n + 1, k) + b_{1}(n + 1, k)P (n + 1, k + 1)
=^{(3.20)}b_{0}(n + 1, k)[b_{0}(n, k)P (n, k) + b_{1}(n, k)P (n, k + 1)]
+ b_{1}(n + 1, k)[b_{0}(n, k + 1)P (n, k + 1) + b_{1}(n, k + 1)P (n, k + 2)]
=^{(3.19)}P (n, k)[b_{0}(n + 1, k)b_{0}(n, k) + b_{1}(n + 1, k)b_{1}(n, k + 1)a_{0}(n, k)]

+ P (n, k + 1)[b_{0}(n + 1, k)b_{1}(n, k) + b_{1}(n + 1, k)(b_{0}(n, k + 1) + b_{1}(n, k + 1)a_{1}(n, k))]. (3.22)

Subsequently, denote f(n + i, k) by f_{i}(k); from now on we suppress the parameter n. Then by the above considerations the expressions f_{0}(k) := f(n, k), …, f_{δ}(k) := f(n + δ, k) can be set up in the form

f_{0}(k) := h_{0}^{(0)}(k)P (k) + ⋯ + h_{s}^{(0)}(k)P (k + s) ; ⋮ ; f_{δ}(k) := h_{0}^{(δ)}(k)P (k) + ⋯ + h_{s}^{(δ)}(k)P (k + s).  (3.23)

Example 9.2 (Cont.). Write f_{0}(k) := f(n, k), f_{1}(k) := f(n+1, k) and f_{2}(k) := f(n+2, k) in the form (3.23), i.e., set h_{0}^{(0)}(k) := 1 and h_{1}^{(0)}(k) := 0, take h_{i}^{(1)}(k) := b_{i}(n, k) for i = 0, 1, and define h_{0}^{(2)}(k) and h_{1}^{(2)}(k) by the coefficients of P (n, k) and P (n, k + 1) in (3.22). □

Consequently, our creative telescoping problem can be stated in the following form. Given (3.23) and (3.9); find g(k) = g_{0}(k)P (k) + ⋯ + g_{s}(k)P (k + s) and c_{0}, …, c_{δ} with

g(k + 1) − g(k) = c_{0}f_{0}(k) + ⋯ + c_{δ}f_{δ}(k).  (3.24)

Example 9.3 (Cont.). Finding solutions c_{0}, c_{1}, c_{2} and g(k) = g_{0}(k)P (k) + g_{1}(k)P (k + 1) for (3.24) is equivalent to looking for solutions c_{0}(n), c_{1}(n), c_{2}(n) and (3.18) for (3.16). □

KEY OBSERVATION: By replacing f(k) with c_{0}f_{0}(k) + ⋯ + c_{δ}f(k) in our telescoping method from Section 8.2 we get the following method to find a solution for (3.24).

1. FIND a solution c_{0}, …, c_{δ} and g_{s}(k) for

∑_{j=0}^{s} a_{s−j}(k + j)g_{s}(k + j + 1) − g_{s}(k) = ∑_{i=0}^{δ} c_{i} ∑_{j=0}^{s} h_{s−j}^{(i)}(k + j);

this is a particular instance of problem (2.25).

2. COMPUTE the remaining g_{0}(k), …, g_{s−1}(k) by

g_{0}(k) = a_{0}(k)g_{s}(k + 1) − ∑_{i=0}^{δ} c_{i}h_{0}^{(i)}(k),

g_{r}(k) = g_{r−1}(k + 1) + a_{r}(k)g_{s}(k + 1) − ∑_{i=0}^{δ} c_{i}h_{r}^{(i)}(k), 0 < r < s.

Example 9.4 (Cont.). We set up our parameterized linear difference equation

− \frac{(k−n+1)^{3}(k+n+2)(k+n+3)}{(k+2)^{2}(k+3)^{2}(k−3n+1)}g_{1}(k + 2) + \frac{(k+n+2)(2k^{2}−6nk+k+3n^{2}−3n)}{(k+2)^{2}(k−3n)}g_{1}(k + 1) − g_{1}(k)

= c_{0}ϕ_{0}(k) + c_{1}ϕ_{1}(k) + c_{2}ϕ_{2}(k)

where ϕ_{i}(n, k) := ∑_{j=0}^{s} h_{s−j}^{(i)}(k + j) are rational functions in n, k. Sigma computes (3.17) and

g_{1}(k) = (k^{2}(k + 1)^{2}(−k + 3n + 1)(17n^{6} − 121kn^{5} + 161n^{5} + 225k^{2}n^{4} − 944kn^{4}

+ 625n^{4} − 177k^{3}n^{3} + 1389k^{2}n^{3} − 2901kn^{3} + 1271n^{3} + 61k^{4}n^{2} − 808k^{3}n^{2} + 3174k^{2}n^{2}
− 4386kn^{2} + 1426n^{2} − 7k^{5}n + 182k^{4}n − 1220k^{3}n + 3183k^{2}n − 3260kn + 836n − 10k^{5}

+ 136k^{4} − 610k^{3} + 1182k^{2} − 952k + 200))/((−k + n + 1)^{3}(−k + n + 2)^{3});

we set g_{0}(k) := a_{0}(k)g_{1}(k + 1) − ∑_{i=0}^{2} c_{i}h_{0}^{(i)}(k) and get the solution (3.18) for (3.16). □

We emphasize that the sketched double-sum approach can be carried over to general multisums; for more details see [Sch05c, Section 4].

9.3. A method for differential equations. By a slight variation of our method [Sch05c], see Section 9.2, one can compute linear differential equations for a sum

S(z) = ∑_{k} ∑_{j} h(z, k, j)  (3.25)

where h(z, k, j) is hypergeometric in k and j and \frac{d}{dz}h(z, k, j)/h(z, k, j) is a rational function in z, k, j. E.g., consider the following identities which will pop up in Section 10.2:

\begin{matrix}∑_{k=0}^{∞}(4k + 1)\frac{(2k)!}{k!^{2}2^{2k}} P (z, k) = 1, ∑_{k=0}^{∞} −(4k + 1)k \frac{(2k + 1)!}{k!^{2}2^{2k−1}}P (z, k) = −z^{2}, \\ ∑_{k=0}^{∞}(4k + 1)k(2k − 1)\frac{(2k + 2)!}{k!^{2}2^{2k+1}}P (z, k) = \frac{1}{4}z^{4} + z^{2} \end{matrix}  (3.26)

where

P (z, k) = ∑_{i=0}^{∞} \frac{\sqrt{π}}{2} \frac{(−z^{2}/4)^{i}(z/2)^{2k}}{i!(2k + i + 1/2)!}.  (3.27)

Note that P (z, k) is equal to j_{2k}(z) where j_{k}(z) are the Spherical Bessel functions of the first kind; see [AS65]. In order to prove these identities, we compute differential equations for the left-hand sides of (3.26) as follows. First, we derive a recurrence of the form

P (z, k + s + 1) = a_{0}(z, k)P (z, k) + ⋯ + a_{s}(z, k)P (z, k + s).  (3.28)

28 CARSTEN SCHNEIDER

With Sigma we get for our sum P (z, k)

In[47]:= sumP = ∑_{i=0}^{∞} \frac{\sqrt{π}}{2} \frac{(−z^{2}/4)^{i}(z/2)^{2k}}{i!(2k + i + 1/2)!};

the recurrence relation

In[48]:= recP = GenerateRecurrence[sumP, k][[1]]/.SUM → P

Out[48]= P[k + 2] == −\frac{4k + 7}{4k + 3}P[k] + \frac{(4k + 5)(16k^{2} + 40k − 2z^{2} + 21)}{(4k + 3)z^{2}}P[k + 1]

Next, we look for a difference-differential equation of the form

\frac{d}{dz} P (z, k) = b_{0}(z, k)P (z, k) + ⋯ + b_{s}(z, k)P (z, k + s).  (3.29)

This can be accomplished by the function call.

In[49]:= recOneDiff = GenerateRecurrence[sumP, k, OneDiffIn → z][[1]]/.SUM → P Out[49]= P^{(0,1)}[z, k] == \frac{8k^{2} + 6k − z^{2}}{(4k + 3)z}P[k] − \frac{z}{4k + 3}P[k + 1]

Remark. Given the summand f(z, k, j), set up the creative telescoping equation

g(z, k, i+1)−g(z, k, i) = c_{0}(z, k)f(z, k, i)+c_{1}(z, k)f(z, k+1, i)+c_{2}(z, k)\frac{d}{dz} f(z, k, i) (3.30)

and solve the underlying problem (2.8) where f_{0}, f_{1}, and f_{2} correspond to f(z, k, i), f(z, k+ 1, i), and \frac{d}{dz}f(z, k, i), respectively. Summing the result (3.30) over k gives Out[49]. Using Out[48] we can compute a differential equation for the left-hand side of (3.26):

In[50]:= GenerateDE[∑_{k=0}^{∞} −(4k + 1)k\frac{(2k + 1)!}{k!^{2}2^{2k−1}}P[k], z, \{recP, P[k]\}, recOneDiff] Out[50]= {2SUM[z] − zSUM′[z] == 0}

Internally, Sigma solves the problem as follows: Given Out[48] and Out[49], we look for a solution g(z, k) = g_{0}(z, k)P (z, k) + g_{1}(z, k)P (z, k + 1) with

g(z, k + 1) − g(z, k) =c_{0}(z)P (z, k) + c_{1}(z)\frac{d}{dz} P (z, k)
=^{Out[49]}c_{0}(z)P (z, k) + c_{1}(z)[\frac{8k^{2} + 6k − z^{2}}{(4k + 3)z} P (z, k) − \frac{z}{4k + 3}P (z, k + 1)].

Together with Out[48] this is nothing else than a certain instance of problem (3.24). Therefore we continue as in Section 9.2 and can compute the solution c_{0} = 2, c_{1} = −z and

g(z, k) = \frac{(2k)!}{k!^{2}2^{2k}} \frac{2(k − 1)k(2k + 1)}{4k + 3}[(z^{2} − 16k^{2} − 16k − 3)P (z, k) + z^{2}P (z, k + 1)].

The correctness can be checked along the lines of Section 9.1. Hence, summing this equation over k gives Out[50]. Completely analogously we find differential equations for the remaining two sums in (3.26):

In[51]:= GenerateDE[∑_{k=0}^{∞}(4k + 1)\frac{(2k)!}{k!^{2}2^{2k}} P[k], z, \{recP, P[k]\}, recOneDiff] Out[51]= {SUM′[z] == 0}

In[52]:= GenerateDE[∑_{k=0}^{∞}(4k + 1)k(2k − 1)\frac{(2k + 2)!}{k!^{2}2^{2k+1}}P[k], z, \{recP, P[k]\}, recOneDiff] Out[52]= {−4(z^{2} + 2)SUM[z] + z(z^{2} + 4)SUM′[z] == 0}

By looking at the solutions of these differential equations, the identities in (3.26) follow.

Summarizing, we propose the following strategy to find a differential equation for (3.25). First, try to compute a difference-differential system (3.28), (3.29) for the inner sum P (z, k) = ∑_{j} h(z, k, j) of (3.25). Then we have methods in hand for the following problem.

Creative telescoping with a difference-differential system

Given δ ∈ ℕ and f(z, k) := h_{0}(z, k)P (z, k) + ⋯ + h_{s}(z, k)P (z, k + s) plus (3.19) and (3.29); find a solution c_{0}(z), …, c_{δ}(z), not all zero, and g(z, k) for

g(z, k + 1) − g(z, k) = c_{0}(z)f(z, k) + c_{1}(z)\frac{d}{dz} f(z, k) + ⋯ + c_{δ}(z)\frac{d^{δ}}{dz^{δ}} f(z, k)  (3.31)

where g(z, k) is of the form g(z, k) = g_{0}(z, k)P (z, k) + ⋯ + g_{s}(z, k)P (z, k + s).

More precisely, with (3.19) and (3.29) we can rewrite

f_{0}(k) := f(z, k), …, f_{δ}(k) := \frac{d^{δ}}{dz^{δ}} f(z, k)

in the form (3.23), respectively. Summarizing, we can reduce problem (3.31) to problem (3.24). A solution of (3.24) by our method from Section 9.2 will provide us with a solution for problem (3.31). Then summing such a solution (3.31) over k gives (under the assumption that (3.31) holds for all 0 ≤ k ≤ n) a differential equation of the form

q(z) = c_{0}(z)S(z) + c_{1}(z)S′(z) + ⋯ + c_{δ}(z)S^{(δ)}(z)

for some function q(z). We emphasize that the sketched double-sum approach for differential equations can be carried over to general multi-sums along the lines of [Sch05c, Section 4].

10. Applications

10.1. Stembridge’s TSPP Theorem. In [APS05] we derived a computer-assisted proof of Stembridge’s Totally Symmetric Plane Partition (TSPP) Theorem [Ste95]. Consider a plane partition with largest part ≤ n, i.e., a matrix

n ≥ a_{11} ≥ a_{12} ≥ a_{13} ≥ …  a_{1r}
∨|  ∨|  ∨|
a_{21} ≥ a_{22}  …  a_{2r}
∨|  ∨|
⋮  ⋮
∨|  ∨|
a_{s1} ≥ a_{s2} ≥ a_{s3} ≥ ⋯ ≥ a_{s,r} ≥ 0

where in each row and column the positive integers are weakly decreasing. Typically, one can represent a plane partition in 3-d as follows:

3 2 1 2 1 0 1 0 0 ⟷

Then a TSPP is such a 3-d cube which can be rotated and reflected without changing the graphical picture. In the beginning of the eighties G.E. Andrews, I.G. Macdonald and R.P. Stanley conjectured that for the number T_{n} of TSPPs with largest part ≤ n the following identity holds:

T_{n} = ∏_{1≤i≤j≤k≤n} \frac{i + j + k − 1}{i + j + k − 2},  n ≥ 1;  (3.32)

cf. [Sta86, case 4] and [Ste95]. Moreover, in [Oka89] a matrix M(n) = (μ(i, j))_{0≤i,j≤n−1} with explicit expressions μ(i, j) was given such that

T_{n}^{2} = det (M).

Since the expressions μ(i, j) are rather complicated, there was no chance to simplify the determinant evaluation in order to show (3.32). Finally, in the beginning of the nineties

G. Andrews could guess a highly non-trivial matrix W = (\begin{matrix}∗ & … & ∗ \\ & ... & ⋮ \\ 0 & & ∗ \end{matrix}) with det(W) = 1

and

MW = (\begin{matrix}∗ & & 0 \\ ⋮ & ... \\ ∗ & … & ∗ \end{matrix}) =: U.  (3.33)

Since det(M) = det(M) det(W) = det(M W) = det(U), one only has to take the product of the diagonal elements of U, which leads exactly to (3.32). The only remaining task is to show (3.33). Unfortunately, the expressions in W (triple sums are involved) look even more complicated than the expressions in M (only double sums are involved). Hence verifying (3.33) leads to identities with quadruple sums! Summarizing, at that time a proof of (3.33) was out of scope – even by using the computer. Finally, J. Stembridge found an elegant proof [Ste95] which combines the combinatorics of Pfaffians and reduction of such to known determinant representations.

Luckily, G. Andrews did not forget his fascinating attempt and sent his problem (3.33) to RISC in 2003. After several weeks of hard work, Peter Paule and I won the battle. In a first step P. Paule managed to simplify in a non-trivial way the underlying identities involving quadruple sums to identities with triple sums. E.g., define

A_{0}(i, m) := ∑_{k=0}^{2m} \binom{i + k − 3}{i − 2}h(k, m) and A_{2}(i, m) := ∑_{k=i}^{2m}(−1)^{k}h(k, m) where

h(k, m) := ∑_{s=0}^{⌊ \frac{2m−k}{2} ⌋−1}\frac{k}{m − s}\binom{m − s}{2m − 2s − k}\frac{(−1)^{s+k}}{2m 4^{s}} ∑_{r=0}^{s} \frac{(m − r)(m)_{r}(−3m − 1)_{r}}{r!(\frac{1}{2} − 2m)_{r}} .

Then one among the many identities looks like follows: for all m ≥ 1 and all 3 ≤ i ≤ 2m+1,

2h(i − 2, m) − 5h(i − 1, m) − A_{0}(i, m) + 6(−1)^{i}A_{2}(i, m) − 3(−1)^{i} ∏_{s=1}^{2m−1} \frac{2(m+s−1)}{2m+s−2} = 0. (3.34)

Here Sigma enters the game: we compute for each of the ingredients recurrences in i. E.g., for A_{0}(i, m) we get a recurrence as follows:

In[53]:= hSum = ∑_{s=0}^{⌊ \frac{2m−k}{2} ⌋−1}\frac{k}{m − s}\binom{m − s}{2m − 2s − k}\frac{(−1)^{s+k}}{2m 4^{s}} ∑_{r=0}^{s} \frac{(m − r)(m)_{r}(−3m − 1)_{r}}{r!(\frac{1}{2} − 2m)_{r}} ;

In[54]:= rec = GenerateRecurrence[hSum, k, FiniteSupport → True][[1]]/.SUM → h Out[54]= −(k − 2m + 2)(k + 2m + 3)h(k + 3)(k + 1)^{2} + 2(k + 2)^{2}(k − 2m)(k + 2m + 1)h[k] +

(−5k^{4} − 29k^{3} + 12m^{2}k^{2} + 6mk^{2} − 58k^{2} + 40m^{2}k + 20mk − 46k + 24m^{2} + 12m − 12)h[k + 1] +

(4k^{4} + 26k^{3} − 12m^{2}k^{2} − 6mk^{2} + 59k^{2} − 28m^{2}k − 14mk + 55k − 12m^{2} − 6m + 18)h[k + 2] == 0 In[55]:= recA0 = GenerateRecurrence[∑_{k=0}^{2m} \binom{i + k − 3}{i − 2}h[k], i, \{rec, h[k]\}, FiniteSupport → True] Out[55]= (i^{2} + i + 2)(i + 2m − 1)(i − 2(m + 1))SUM[i] − (i + 3)(i^{3} − i^{2} + 2i + 4m^{2} + 2m − 2)SUM[i + 1] −

(i − 3)(i^{3} + i^{2} + 2i − 4m^{2} − 2m + 2)SUM[i + 2] + (i^{2} − i + 2)(i − 2m + 1)(i + 2m + 2)SUM[i + 3] == 0

Given all the recurrences, see [APS05], we combine them to one recurrence by using the Mathematica package GeneratingFunctions [Mal96], which is based on the ideas of [SZ94]. By checking initial values we show that all the sums combined in (3.34) evaluate to zero. We conclude our result with Zeilberger’s Opinion 65 [Zei05]:

Seeing all the details, (that nowadays can (and should!) be easily relegated to the computer), even if they are extremely hairy, is a hang-up that traditional mathematicians should learn to wean themselves from. A case in point is the excellent but unnecessarily long-winded recent article [APS05]. It is a new, computer-assisted proof, of John Stembridge’s celebrated TSPP theorem. It is so long because they insisted on showing explicitly all the hairy details, and easilyreproducible-by-the-reader ”proof certificates”. It would have been much better if they would have first applied their method to a much simpler case, that the reader can easily follow, that would take one page, and then state that the same method was applied to the complicated case of Stembridge’s theorem and the result was TRUE. For those poor people who are unable or unwilling to run the program themselves, they could have posted the computer output on their websites, but please, have mercy on the rain forest! You don’t need 30 pages, and frankly all this EXPLICIT LANGUAGE of hairy computer output is almost pornographic.

Here I would like to mention that our TSPP-proof is indeed hairy and highly non-trivial. The challenge was to illustrate that such non-trivial problems can be proven completely rigorously with the computer. As a consequence, we derived proof certificates that do not fill 30 pages, as Doron mentioned, but 80 pages :-) Interestingly enough, there was quite some human interaction necessary, e.g., to avoid summation over poles. Exactly this kind of problems have been checked carefully in the extended version [APS04].

10.2. Lost proofs of the Handbook of Mathematical Functions. In spring 2005 Frank Olver asked Peter Paule if the algorithms of the RISC combinatorics group can provide proofs of about twelve identities in the Handbook of Mathematical Functions [AS65]. The real challenge was that the original proofs have been lost and no alternative proofs were known. After a long weekend the Comb-group could find computer proofs for each of the identities [GKO^{+}06]. One of the identities is [AS65, Equ. (10.1.48)]

J_{0}(z sin θ) = ∑_{k=0}^{∞}(4k + 1)\frac{(2k)!}{2^{2k}k!^{2}} j_{2k}(z)P_{2k}(cos θ)  (3.35)

where P_{k}(z) are the Legendre polynomials, J_{k}(z) are the Bessel functions of the first kind, and j_{k}(z) are Spherical Bessel functions of the first kind.

When I have seen this identity for the first time, I wondered myself, what Sigma could have in common with all these functions. After having a closer look, it was clear: The problem can be transformed exactly to the input class of Sigma. First we apply the substitution t := cos θ. Hence with sin(θ) = \sqrt{1 − t^{2}} our identity reads as

J_{0}(z\sqrt{1 − t^{2}}) = ∑_{k=0}^{∞}(4k + 1)\frac{(2k)!}{2^{2k}k!^{2}} j_{2k}(z)P_{2k}(t).

Moreover, by hypergeometric series representations from [AS65] we get

J_{0}(z\sqrt{1 − t^{2}}) = ∑_{n=0}^{∞} \frac{(−\frac{1}{2}z^{2})^{n}}{n!^{2}} (\frac{1 − t^{2}}{2} )^{n} , P_{2k}(t) = ∑_{i=0}^{∞} \frac{(−2k)_{i}(2k + 1)_{i}}{i!^{2}} (\frac{1 − t}{2})^{i} ,

and P (k) := P (z, k) = j_{2k}(z) from (3.27). Summarizing, we have to show

∑_{n=0}^{∞} a_{n} (\frac{1 − t^{2}}{2} )^{n} = ∑_{n=0}^{∞} b_{n} (\frac{1 − t}{2})^{n}  (3.36)

where

a_{n} = \frac{(−\frac{1}{2}z^{2})^{n}}{n!^{2}} ,  and  b_{n} = ∑_{k=0}^{∞}(4k + 1)\frac{(2k)!}{2^{2k}k!} \frac{(−2k)_{n}(2k + 1)_{n}}{n!^{2}} P (k).

Given this representation, our RISC-packages can finish the job. First we derive recurrences for a_{k} and b_{k}; the recurrence z^{2}a_{n} + 2(n + 1)^{3}a_{n+1} = 0 is immediate. We get a recurrence for b_{k} by the following function call; here we use the recurrence Out[48] for P (k).

In[56]:= recB = GenerateRecurrence[

∑_{k=0}^{∞}(4k + 1)\frac{(2k)!}{2^{2k}k!} \frac{(−2k)_{n}(2k + 1)_{n}}{n!^{2}} P[k], n, \{recP, P[k]\}][[1]]/.SUM → b

Out[56]= 4(n + 3)z^{2}b[n] − 2(2n + 5)z^{2}b[n + 1] − (n + 2)(n^{2} + 5n + 6 − z^{2})b[n + 2] + (n + 2)(n + 3)^{2}b[n + 3] == 0

Next, we use closure properties of holonomic functions [SZ94] in order to compute differential equations for both sides of (3.36). The package [Mal96] (which is inspired by [SZ94]) helps here.

In[57]:= << GeneratingFunctions.m GeneratingFunctions Package by Christian Mallinger \begin{matrix}c \\ ◯ \end{matrix} RISC-Linz

Namely, given the recurrence for b_{k}, we get a differential equation for ∑_{k=0}^{∞} b_{k}t^{k} by the function call

In[58]:= deB = RE2DE[recB, b[k], B[t]]

Out[58]= 12z^{2}B[t] + 10(2tz^{2} − z^{2})B′[t] + (4t^{2}z^{2} − 4tz^{2} + z^{2} − 6)B′′[t] − 3(2t − 1)B^{(3)}[t] − t(t − 1)B^{(4)}[t] == 0

Then by the substitution t → (1 − t)/2 we compute

In[59]:= deB = ACompose[deB, B[t] == (1 − t)/2, B[t]]

Out[59]= −3z^{2}B[t] − 5tz^{2}B′[t] − (t^{2}z^{2} − 6)B′′(t) + 6tB^{(3)}[t] + (t + 1)(t − 1)B^{(4)}[t] == 0

for B(t) = ∑_{k=0}^{∞} b_{k} (\frac{1−t}{2})^{k}. Similarly, we obtain

In[60]:= deA = −t^{3}z^{2}A[t] + (t^{2} + 1)A′[t] + t(t + 1)(t − 1)A′′[t] == 0;

for A(t) = ∑_{k=0}^{∞} a_{k} (\frac{1−t^{2}}{2} )^{k}. With

In[61]:= DEPlus[(deA/.A → C), (deB/.B → C), C[t]]

Out[61]= −3z^{2}C[t] − 5tz^{2}C′[t] − (t^{2}z^{2} − 6)C′′(t) + 6tC^{(3)}[t] + (t + 1)(t − 1)C^{(4)}[t] == 0

we find a differential equation for C(t) := A(t) − B(t). By inspection we get the initial conditions b_{0} = a_{0}, b_{1} = 2a_{1}, and b_{2} = −2a_{1} + 4a_{2}. Notice that these are the already proven identities (3.26). This completes the proof of (3.36) and therefore of (3.35).

11. Conclusion

I illustrated how Sigma can handle non-trivial summation problems, most of them related to combinatorial questions. As a conclusion I want to emphasize that the title can be reversed: “Combinatorics assists symbolic summation”. Namely, most of the examples in this survey article were important case studies to improve the summation package Sigma. Even more, challenging problems, like the TSPP-problem, were the source to extend Sigma.

I am looking forward to see how symbolic summation and combinatorics will inspire each other in the future.

The first public release of Sigma is planned for summer 2007.

References

[ABP95] S.A. Abramov, M. Bronstein, and M. Petkovšek. On polynomial solutions of linear operator equations. In T. Levelt, editor, Proc. ISSAC’95, pages 290–296. ACM Press, 1995.

[Abr75] S.A. Abramov. The rational component of the solution of a first-order linear recurrence relation with a rational right-hand side. U.S.S.R. Comput. Maths. Math. Phys., 15:216–221, 1975. Transl. from Zh. vychisl. mat. mat. fiz. 15, pp. 1035–1039, 1975.

[Abr89a] S.A. Abramov. Problems in computer algebra that are connected with a search for polynomial solutions of linear differential and difference equations. Moscow Univ. Comput. Math. Cybernet., 3:63–68, 1989.

[Abr89b] S.A. Abramov. Rational solutions of linear differential and difference equations with polynomial coefficients. U.S.S.R. Comput. Math. Math. Phys., 29(6):7–12, 1989.

[Abr95] S.A. Abramov. Rational solutions of linear difference and q-difference equations with polynomial coefficients. In T. Levelt, editor, Proc. ISSAC’95, pages 285–289. ACM Press, 1995.

[Abr03] S.A. Abramov. When does Zeilberger’s algorithm succeed? Adv. in Appl. Math., 30:424–441, 2003.

[ACGL04] S.A. Abramov, J.J. Carette, K.O. Geddes, and H.Q. Le. Telescoping in the context of symbolic summation in Maple. J. Symbolic Comput., 38:1303–1326, 2004.

[AL05] S.A. Abramov and H.Q. Le. On the order of the recurrence produced by the method of creative telescoping. Discrete Math., 298:2–17, 2005.

[ALP03] S.A. Abramov, H.Q. Le, and M. Petkovsek. Rational canonical forms and efficient representations of hypergeometric terms. In J.R. Sendra, editor, Proc. ISSAC’03, pages 7–14. ACM Press, 2003.

[AP94] S.A. Abramov and M. Petkovšek. D’Alembertian solutions of linear differential and difference equations. In J. von zur Gathen, editor, Proc. ISSAC’94, pages 169–174. ACM Press, 1994.

[AP99] G.E. Andrews and P. Paule. MacMahon’s partition analysis IV: Hypergeometric multisums. Sém. Lothar. Combin., B42i:1–24, 1999.

[AP02] S.A. Abramov and M. Petkovšek. Rational normal forms and minimal decompositions of hypergeometric terms. J. Symbolic Comput., 33(5):521–543, 2002.

[AP06] S.A. Abramov and M. Petkovšek, 2006. Private Communication.

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

[APS04] G.E. Andrews, P. Paule, and C. Schneider. Plane partitions VI: Stembridge’s TSPP Theorem — a detailed algorithmic proof. Technical Report 04-08, RISC-Linz, J. Kepler University, 2004.

[APS05] G.E. Andrews, P. Paule, and C. Schneider. Plane Partitions VI: Stembridge’s TSPP Theorem. Advances in Applied Math., 34(4):709–739, 2005.

[AS65] M. Abramowitz and I.A. Stegun. Handbook of Mathematical Functions. Dover Publications, Inc, New York, 1965.

[AU85] G.E. Andrews and K. Uchimura. Identities in combinatorics IV: Differentiation and harmonic numbers. Util. Math., 28:265–269, 1985.

[AZ06] M. Apagodu and D. Zeilberger. Multi-variable Zeilberger and Almkvist–Zeilberger algorithms and the sharpening of Wilf–Zeilberger theory. Advances in Applied Math., 37:139–152, 2006.

[BG96] J.M. Borwein and R. Girgensohn. Evaluation of triple Euler sums. Electron. J. Combin., 3:1–27, 1996.

[BK99] H. Böing and W. Koepf. Algorithms for q-hypergeometric summation in computer algebra. J. Symbolic Comput., 28:777–799, 1999.

[BP99] A. Bauer and M. Petkovšek. Multibasic and mixed hypergeometric Gosper-type algorithms. J. Symbolic Comput., 28(4–5):711–736, 1999.

[BPP^{+}06] A. Bećirović, P. Paule, V. Pillwein, A. Riese, C. Schneider, and J. Schöberl. Hypergeometric

summation algorithms for high order finite elements. Computing, 78(3):235–249, 2006.

[Bro00] M. Bronstein. On solutions of linear ordinary difference equations in their coefficient field. J. Symbolic Comput., 29(6):841–877, 2000.

[Bro05] M. Bronstein, 2005. Private Communication.

[Cal94] N.J. Calkin. A curious binomial identity. Discrete Math., 131(1-3):335–337, 1994.

[CHM05] W.Y.C. Chen, Q.H. Hou, and Y.P. Mu. Applicability of the q-analogue of Zeilberger’s algorithm. J. Symbolic Comput., 39:155–170, 2005.

[CHM06] W.Y.C. Chen, Q.H. Hou, and Y.P. Mu. A telescoping method for double summations. J. Comput. Appl. Math., 196(2):553–566, 2006.

[Chu05] W. Chu. Harmonic number identities and Hermite-Padé approximations to the logarithm function. Journal of Approximation Theory, 137(1):42 – 56, 2005.

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

[Coh65] R.M. Cohn. Difference Algebra. Interscience Publishers, John Wiley & Sons, 1965.

[CS98] F. Chyzak and B. Salvy. Non-commutative elimination in ore algebras proves multivariate identities. J. Symbolic Comput., 26(2):187–227, 1998.

[Dix03] A.C. Dixon. Summation of a certain series. Proc. London Math. Soc. (1), 35:285–289, 1903.

[DPSW06a] K. Driver, H. Prodinger, C. Schneider, and J.A.C. Weideman. Padé approximations to the logarithm II: Identities, recurrences, and symbolic computation. Ramanujan J., 11(2):139– 158, 2006.

[DPSW06b] K. Driver, H. Prodinger, C. Schneider, and J.A.C. Weideman. Padé approximations to the logarithm III: Alternative methods and additional results. Ramanujan J., 12(3):299–314, 2006.

[Fas45] Sister Mary Celine Fasenmyer. Some generalized hypergeometric polynomials. PhD thesis, University of Michigan, 1945.

[FK00] M. Fulmek and C. Krattenthaler. The number of rhombus tilings of a symmetric hexagon which contains a fixed rhombus on the symmetric axis, II. European J. Combin., 21(5):601– 640, 2000.

[FS98] P. Flajolet and B. Salvy. Euler sums and contour integral representations. Experimental Math., 7(1):15–35, 1998.

[Ger04] J. Gerhard. Modular algorithms in symbolic summation and symbolic integration. Lecture Notes in Computer Science. Springer, 2004.

[GKO^{+}06] S. Gerhold, M. Kauers, F.W.J. Olver, P. Paule, C. Schneider, and B. Zimmermann. Computer proofs of some identities involving derivatives of bessel functions with respect to order. In preparation, 2006.

[GKP94] R.L. Graham, D.E. Knuth, and O. Patashnik. Concrete Mathematics: a foundation for computer science. Addison-Wesley Publishing Company, Amsterdam, 2nd edition, 1994.

[Gos78] R.W. Gosper. Decision procedures for indefinite hypergeometric summation. Proc. Nat. Acad. Sci. U.S.A., 75:40–42, 1978.

[Hir96] M. Hirschhorn. Calkin’s binomial identity. Discrete Math., 273-278, 1996.

[Hoe98] M. van Hoeij. Rational solutions of linear difference equations. In O. Gloor, editor, Proc. ISSAC’98, pages 120–123. ACM Press, 1998.

[Hoe99] M. van Hoeij. Finite singularities and hypergeometric solutions of linear recurrence equations. J. Pure Appl. Algebra, 139(1-3):109–131, 1999.

[Hor92] J. Hornegger. Hypergeometrische Summation und polynomiale Rekursion. Master’s thesis, University Erlangen, 1992.

[HS99] P.A. Hendriks and M.F. Singer. Solving difference equations in finite terms. J. Symbolic Comput., 27(3):239–259, 1999.

[Kar81] M. Karr. Summation in finite terms. J. ACM, 28:305–350, 1981.

[Kar85] M. Karr. Theory of summation in finite terms. J. Symbolic Comput., 1:303–315, 1985.

[Koe95] W. Koepf. Algorithms for m-fold hypergeometric summation. J. Symbolic Comput., 20(399-417), 1995.

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

[Kra03] C. Krattenthaler, 2003. Private Communication.

[KS06a] M. Kauers and C. Schneider. Application of unspecified sequences in symbolic summation. In J.G. Dumas, editor, Proc. ISSAC’06., pages 177–183. ACM Press, 2006.

[KS06b] M. Kauers and C. Schneider. Indefinite summation with unspecified summands. Discrete Math., 306(17):2021–2140, 2006.

[M06] Kauers M. Sum Cracker – A Package for Manipulating Symbolic Sums and Related Objects. J. Symbolic Computat., 41(9):1039–1057, 2006.

[Mal96] C. Mallinger. Algorithmic manipulations and transformations of univariate holonomic functions and sequences. Master’s thesis, RISC, J. Kepler University, Linz, 1996.

[MZ05] M. Mohammed and D. Zeilberger. Sharp upper bounds for the orders of the recurrences outputted by the Zeilberger and q-Zeilberger algorithms. J. Symbolic Comput., 39:201–207, 2005.

[Oka89] S. Okada. On the generating functions for certain classes of plane partitions. J. Combin. Theory Ser. A, 51:1–23, 1989.

[Pau95] P. Paule. Greatest factorial factorization and symbolic summation. J. Symbolic Comput., 20(3):235–268, 1995.

[Pau04] P. Paule. Contiguous relations and creative telescoping. Preprint, 2004.

[Pet92] M. Petkovšek. Hypergeometric solutions of linear recurrences with polynomial coefficients. J. Symbolic Comput., 14(2-3):243–264, 1992.

[PP05] A. Panholzer and H. Prodinger. Computer-free evaluation of an infinite double sum. Sém. Lothar. Combin., 55:1–3, 2005.

[PR97] P. Paule and A. Riese. A Mathematica q-analogue of Zeilberger’s algorithm based on an algebraically motivated aproach to q-hypergeometric telescoping. In M. Ismail and M. Rahman, editors, Special Functions, q-Series and Related Topics, volume 14, pages 179–210. Fields Institute Toronto, AMS, 1997.

[PS95a] P. Paule and M. Schorn. A Mathematica version of Zeilberger’s algorithm for proving binomial coefficient identities. J. Symbolic Comput., 20(5-6):673–698, 1995.

[PS95b] P. Paule and V. Strehl. Symbolic summation - some recent developments. In J. Fleischer et al., editor, Computer Algebra in Science and Engineering - Algorithms, Systems, and Applications, pages 138–162. World Scientific, Singapore, 1995.

[PS95c] R. Pirastu and V. Strehl. Rational summation and Gosper-Petkovšek representation. J. Symbolic Comput., 20:617–635, 1995.

[PS03] P. Paule and C. Schneider. Computer proofs of a new family of harmonic number identities. Adv. in Appl. Math., 31(2):359–378, 2003.

[PS04] P. Paule and C. Schneider. Creative telescoping for hypergeometric double sums. Preprint, 2004.

[PS07] R. Pemantle and C. Schneider. When is 0.999... equal to 1? To appear in Amer. Math. Monthly, April 2007.

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

[Rie03] A. Riese. qMultisum - A package for proving q-hypergeometric multiple summation identities. J. Symbolic Comput., 35:349–377, 2003.

[Ris69] R. Risch. The problem of integration in finite terms. Trans. Amer. Math. Soc., 139:167–189, 1969.

[Ris70] R. Risch. The solution to the problem of integration in finite terms. Bull. Amer. Math. Soc., 76:605–608, 1970.

[Sch00] C. Schneider. An implementation of Karr’s summation algorithm in Mathematica. Sém. Lothar. Combin., S43b:1–10, 2000.

[Sch01] C. Schneider. Symbolic summation in difference fields. Technical Report 01-17, RISC-Linz, J. Kepler University, November 2001. PhD Thesis.

[Sch02] C. Schneider. How one can play with sums. In H. Kredel and W. Seiler, editors, Proceedings of the 8th Rhine Workshop on Computer Algebra, pages 73–83, Mannheim, Germany, 2002.

[Sch04a] C. Schneider. A collection of denominator bounds to solve parameterized linear difference equations in ΠΣ-extensions. In D. Petcu et al., editor, Proc. SYNASC04, pages 269–282. Mirton Publishing, 2004.

[Sch04b] C. Schneider. The summation package Sigma: Underlying principles and a rhombus tiling application. Discrete Math. Theor. Comput. Sci., 6(2):365–386, 2004.

[Sch04c] C. Schneider. Symbolic summation with single-nested sum extensions. In J. Gutierrez, editor, Proc. ISSAC’04, pages 282–289. ACM Press, 2004.

[Sch05a] C. Schneider. Degree bounds to find polynomial solutions of parameterized linear difference equations in ΠΣ-fields. Appl. Algebra Engrg. Comm. Comput., 16(1):1–32, 2005.

[Sch05b] C. Schneider. Finding telescopers with minimal depth for indefinite nested sum and product expressions. In M. Kauers, editor, Proc. ISSAC’05, pages 285–292. ACM, 2005.

[Sch05c] C. Schneider. A new Sigma approach to multi-summation. Advances in Applied Math., 34(4):740–767, 2005.

[Sch05d] C. Schneider. Product representations in ΠΣ-fields. Annals of Combinatorics, 9(1):75–99, 2005.

[Sch05e] C. Schneider. Solving parameterized linear difference equations in terms of indefinite nested sums and products. J. Differ. Equations Appl., 11(9):799–821, 2005.

[Sch06] C. Schneider. Some notes on “When is 0.999... equal to 1?”. In Thierry Coquand et al., editor, Mathematics, Algorithms, Proofs, number 05021 in Dagstuhl Seminar Proceedings. Internationales Begegnungs- und Forschungszentrum (IBFI), Schloss Dagstuhl, Germany, 2006.

[Sch07] C. Schneider. Simplifying sums in ΠΣ^{∗}-extensions. To appear in J. Algebra Appl., 2007.

[Sta86] R.P. Stanley. Symmetries of plane partitions. J. Combin. Theory Ser. A, 43:103–113, 1986.

[Ste95] J. Stembridge. The enumeration of totally symmetric plane partitions. Advances in Math., 111:227–243, 1995.

[SZ94] B. Salvy and P. Zimmermann. Gfun: A package for the manipulation of generating and holonomic functions in one variable. ACM Trans. Math. Software, 20:163–177, 1994.

[Ver74] P. Verbaeten. The automatic construction of pure recurrence equations. ACM-SIGSAM Bulletin, 8:96–98, 1974.

[Weg97] K. Wegschaider. Computer generated proofs of binomial multi-sum identities. Diploma thesis, RISC Linz, Johannes Kepler University, 1997.

[Wei05] J.A.C. Weideman. Padé approximations to the logarithm I: Derivation via differential equations. Quaestiones Mathematicae, 28:375–390, 2005.

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

[Zei90] D. Zeilberger. A holonomic systems approach to special functions identities. J. Comput. Appl. Math., 32:321–368, 1990.

[Zei91] D. Zeilberger. The method of creative telescoping. J. Symbolic Comput., 11:195–204, 1991.

[Zei05] D. Zeilberger. Opinion 65. http://www.math.rutgers.edu/ zeilberg/Opinion65.html, 2005.

[Zha99] Z. Zhang. A kind of binomial identity. Discrete Math., 196(1-3):291–298, 1999.

Research Institute for Symbolic Computation, J. Kepler University Linz, A-4040 Linz, Austria

  E-mail address: Carsten.Schneider@risc.uni-linz.ac.at