# Reduction-Based Creative Telescoping for Definite Summation of D-Finite Functions

Hadrien Brochet, Bruno Salvy

► To cite this version:

Hadrien Brochet, Bruno Salvy. Reduction-Based Creative Telescoping for Definite Summation of D-Finite Functions. Journal of Symbolic Computation, 2024, 125, pp.102329. 10.1016/j.jsc.2024.102329. hal-04295759

HAL Id: hal-04295759

https://hal.science/hal-04295759v1

Submitted on 20 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# REDUCTION-BASED CREATIVE TELESCOPING FOR DEFINITE SUMMATION OF D-FINITE FUNCTIONS

HADRIEN BROCHET AND BRUNO SALVY

Dedicated to the memory of Marko Petkovšek

ABSTRACT. Creative telescoping is an algorithmic method initiated by Zeilberger to compute definite sums by synthesizing summands that telescope, called certificates. We describe a creative telescoping algorithm that computes telescopes for definite sums of D-finite functions as well as the associated certificates in a compact form. The algorithm relies on a discrete analogue of the generalized Hermite reduction, or equivalently, a generalization of the Abramov-Petkovšek reduction. We provide a Maple implementation with good timings on a variety of examples.

## 1. INTRODUCTION

The algorithmic computation of definite sums originates in Zeilberger's algorithm in the 1990's [36, 38, 35]. Initially designed to deal with hypergeometric sums, his method of creative telescoping has been extended to differential settings [1, 37, 30, 29] and next generalized to the large class of D-finite functions by Chyzak [16]. In order to compute a definite sum of $F(t, x_1, \dots, x_m)$ with respect to $t$, where each $x_i$ is a variable with respect to which one can apply a linear operator $\partial_i$ (generally, differentiation or shift or q-shift operator), the creative telescoping algorithm constructs identities of the form

$$ (1) \qquad \sum_{\alpha} c_{\alpha}(x_1, \dots, x_m) \partial^{\alpha}(F) = G(t+1, x_1, \dots, x_m) - G(t, x_1, \dots, x_m). $$

Here, the sum is over a finite number of multi-indices $\alpha$ and we use the multi-exponent notation $\partial^{\alpha} = \partial_1^{\alpha_1} \cdots \partial_m^{\alpha_m}$. In the original version for hypergeometric summation, the monomials $\partial^{\alpha}(F)$ are simply successive shifts $F(t, n), F(t, n+1), F(t, n+2), \dots$ of a hypergeometric sequence $F(t, n)$. Identities obtained that way can often be summed over $t$. The right-hand side telescopes by design. Since the coefficients $c_\alpha$ do not depend on the variable $t$, the left-hand side results in an operator applied to the definite sum of $F$. From there, other algorithms can be applied to compute information on the sum. The operator in the left-hand side of Eq. (1) is called a telescoper of $F$ and the function $G$ in the right-hand side is the corresponding certificate. Chyzak's algorithm also deals with the differential analogue of Eq. (1) where the right-hand side is a derivative; it is used to compute information on definite integrals. Chyzak's algorithm, like Zeilberger's, looks for telescopes with an increasing number of monomials $\partial^{\alpha}$ with indeterminate coefficients $c_\alpha$ and determines $c_\alpha$ such that a certificate $G$ exists in the vector space generated by the $\partial^{\beta}(F)$ for $\beta \in \mathbb{N}^{m+1}$ over the field of rational functions. The conditions of being D-finite is that this vector space has finite dimension, which allows for the existence of algorithms based on linear algebra. If no certificate exists, then the support is increased and the process is iterated. This stops either when sufficiently many operators have been found or when a prescribed bound on the orders is reached. (In the original hypergeometric case, no bound on the order is fixed a priori and termination is guaranteed for the family of proper hypergeometric terms [35].) Efficiency issues with this approach have led to the development of heuristics and a very useful Mathematica implementation by Koutschan [26].

The most recent approach to deal with the efficiency issues with creative telescoping was initiated by Bostan, Chen, Chyzak and many co-authors who developed a class of reduction-based algorithms [4, 5, 12, 7, 13, 14]. These algorithms avoid the computation of potentially large certificates. In the differential case, where the right-hand side of Eq. (1) is replaced by a derivative $\partial_t(G)$, the principle is to use a variant of Hermite reduction to compute an additive decomposition of each monomial in the form

$$ (2) \qquad \partial^{\alpha}(F) = R_{\alpha}(t, x_1, \dots, x_m)F + \partial_t(G_{\alpha}), $$
where $R_\alpha$ is a rational function with a certain minimality property. A telescoper is found by looking for a linear dependency between these rational functions for a family of monomials $\partial^\alpha$. The computation of the rational function $R_\alpha$ by Hermite reduction works by getting rid of multiple poles and isolating a polynomial part. This was first done for the integration of bivariate rational functions [4], of hyperexponential functions [5] and of mixed hypergeometric-hyperexponential functions [7]. In these three cases, the vector space generated by the functions $\partial^\beta(F)$ for $\beta \in \mathbb{N}^2$ has dimension only 1 over the rational functions. For summation, the analogous problem for bivariate hypergeometric sequences was solved by replacing the Hermite reduction by a modified Abramov-Petkovšek reduction, thereby providing a faster variant of Zeilberger's algorithm [12]. For bivariate problems of dimension larger than 1, the method was extended to the integration of bivariate algebraic functions [13], of Fuchsian functions [14] and more recently of P-recursive sequences [10, 20] by means of suitable integral bases. An extension to the integration of purely differential bivariate D-finite functions in arbitrary dimension was first achieved by turning the differential equations satisfied by the function to be integrated into first-order differential systems; then, a variant of Hermite reduction can be designed at the level of vectors of rational functions [31]. This approach generalizes to purely differential D-finite functions in more variables [33].

Another method relies on cyclic vectors and allows the integration of arbitrary D-finite functions [6]. Without loss of generality, we assume that $F$ is a cyclic vector for $\partial_t$, which means that all monomials $\partial^\alpha(F)$ rewrite as $M_\alpha(F)$ with $M_\alpha$ a linear operator in $\partial_t$ only. (If $F$ is not a cyclic vector, one finds a cyclic vector $G$, $F = M_F(G)$ for some linear operator $M_F$ in $\partial_t$ only and the rest of the reasoning is unchanged.) Next, for any rational function $u$ and any linear operator $M$ in $\partial_t$, repeated integration by parts implies Lagrange's identity

$$uM(F) - M^*(u)F = \partial_t(P_M(F, u)),$$

where $M^*$ is the adjoint of $M$ and $P_M$ is linear in $F, \partial_t(F), \dots$ and $u, \partial_t(u), \dots$. Thus, a reduction-based algorithm is obtained by the additive decomposition of Eq. (2) with $R_\alpha$ a solution of the generalized Hermite reduction

$$R_{\alpha} = M_{\alpha}^{*}(1) \mod \mathrm{Im}(L^{*}).$$

The tools used in the reduction modulo the image of the linear differential operator $L^*$ are classical techniques used when looking for rational solutions of linear differential equations.

This method of integration using generalized Hermite reduction extends to other contexts. This has been done in high generality in a preprint by van der Hoeven [32]. Our approach here is different: we focus on the case of summation only and give a simple self-contained presentation of the corresponding algorithm; our algorithm returns operators of minimal order<fnref n="1" />; we make the choice to avoid algebraic extensions when possible; we present a Maple implementation that performs well in practice. Note that while in terms of complexity, minimal operators cannot be computed in polynomial time in general, in practice this does not seem to be an obstacle.

## 2. EXAMPLE

The multiplication theorem for Bessel functions of the first kind $J_\nu$ states that [27, 10.23.1]

$$(3) \qquad J_{\nu}(\lambda z) = \lambda^{\nu} \sum_{n=0}^{\infty} \frac{(-1)^n (\lambda^2 - 1)^n (z/2)^n}{n!} J_{\nu+n}(z).$$

This can be proved automatically by showing that the left-hand side and the right-hand side satisfy the same set of mixed differential-difference equations with sufficiently many identical initial conditions.

<footnote n="1">¹Remark 5.6 in [32] seems to allude to a way of doing this, but the relevant space $E$ may contain rational functions that are not in the image of $L$. For example take $L = 1/z + 1/(z-1)\sigma^{-1}$, $\alpha = 0$, and $A = \{\alpha\}$, then one can check that $1/z \notin \mathrm{Im}(L)$ but $1/z \in E$.</footnote>
We write $F$ for the summand in Eq. (3). It is a function of the four variables $\nu, n, z, \lambda$. Basic properties of the Bessel function give the following four equations:

$$ (4) \qquad (\lambda^2 - 1)zS_{\nu}(F) + 2(n+1)S_{n}(F) = 0, $$

$$ (5) \qquad (\lambda^2 - 1)\partial_{\lambda}(F) - 2n\lambda F = 0, $$

$$ (6) \qquad (1 - \lambda^2)z\partial_z(F) + 2(n+1)S_n(F) + (\lambda^2 - 1)(2n + \nu)F = 0, $$

$$ (7) \qquad 4(n+1)(n+2)S_n^2(F) + 4(\lambda^2 - 1)(n+1)(n+\nu+1)S_n(F) + z^2(\lambda^2 - 1)^2 F = 0, $$

where $S_\nu$ denotes the shift with respect to $\nu$: $S_\nu : G(\nu) \mapsto G(\nu + 1)$ and similarly for $S_n$, while $\partial_z$ and $\partial_\lambda$ denote partial derivatives. These equations show that any shift or derivative $S_\nu^a \partial_\lambda^b \partial_z^c S_n^d F$ of $F$ with nonnegative integers $a, b, c, d$ rewrites as a $\mathbb{Q}(\nu, \lambda, z, n)$-linear combination of $F$ and $S_n(F)$. In particular, this implies that $F$ is D-finite with respect to these variables. The aim of creative telescoping is to find a similar set of equations, in the variables $\nu, \lambda, z$ only, for the sum in Eq. (3).

Let $\Delta_n$ be the difference operator $\Delta_n = S_n - 1$. Any product $\phi(n)S_n F$ with $\phi \in \mathbb{Q}(\nu, \lambda, z, n)$ can be rewritten $\phi(n-1)F + \Delta_n(\phi(n-1)F)$, i.e., as the sum of a rational function times $F$ plus a difference, that would telescope under summation. Consequently, any $S_\nu^a \partial_\lambda^b \partial_z^c S_n^d F$ with nonnegative integers $a, b, c, d$ rewrites in the form

$$ (8) \qquad \partial^{\alpha}(F) = R_{\alpha}F + \Delta_{n}(G_{\alpha}), $$

with $R_\alpha$ a rational function in $\mathbb{Q}(\nu, \lambda, z, n)$. This reduction as a sum of a rational function plus a difference is a general phenomenon (see Section 3.2). Reduction-based creative telescoping works by reducing this rational function further by pulling out parts that can be incorporated into the difference $\Delta_n(G_\alpha)$. Denote by $L$ the recurrence operator such that Eq. (7) is $LF = 0$. The adjoint of $L$ (see Proposition 2) is

$$ L^* = 4(n-1)nS_n^{-2} + 4(\lambda^2 - 1)n(n + \nu)S_n^{-1} + z^2(\lambda^2 - 1)^2 $$

where $S_n^{-1} : g(n) \mapsto g(n-1)$. Proposition 3 shows that a rational function $R$ is of the form $\Delta_n(M(F))$ for a recurrence operator $M(S_n)$ if and only if $R$ is in the image $L^*(\mathbb{Q}(\nu, \lambda, z, n))$. This is the basis for the computation of relations of the form (8) where $R_\alpha$ is now a reduced rational function (in a sense made precise in Definition 1). The next step is to look for linear combinations of these rational functions that yield telescopes.

The starting point is the monomial 1, which decomposes as

$$ (9) \qquad 1 \cdot F = 1 \cdot F + \Delta_n(0). $$

Using Eq. (5), the monomial $\partial_\lambda$ rewrites

$$ (10) \qquad \partial_{\lambda}(F) = \frac{2n\lambda}{\lambda^2 - 1} F + \Delta_{n}(0) $$

and the rational function is reduced. Taking the derivative of this equation and using Eq. (5) again gives a similar equation for $\partial_\lambda^2(F)$:

$$ \begin{aligned} \partial_{\lambda}^{2}(F) &= -\frac{2n(\lambda^{2} + 1)}{(\lambda^{2} - 1)^{2}}F + \frac{2n\lambda}{\lambda^{2} - 1}\partial_{\lambda}(F) + \Delta_{n}(0), \\ &= -\frac{2n(\lambda^{2} + 1)}{(\lambda^{2} - 1)^{2}}F + \frac{(2n\lambda)^{2}}{(\lambda^{2} - 1)^{2}}F + \Delta_{n}(0). \end{aligned} $$

This time, a reduction is possible. Indeed, Proposition 3 implies that $L^*(1)F$ is a difference $\Delta(A_n)$ (where $A_n$ can be computed explicitly). Since

$$ L^*(1) = 4\lambda^2 n^2 + 4((\lambda^2 - 1)\nu - 1)n + z^2(\lambda^2 - 1)^2, $$

we can eliminate the term in $n^2$ in the expression of $\partial_\lambda^2(F)$ to get

$$ (11) \qquad \partial_{\lambda}^{2}(F) = -\frac{2n(2\nu + 1)}{(\lambda^{2} - 1)}F - z^{2}F + \Delta_{n}\left(\frac{A_{n}}{(\lambda^{2} - 1)^{2}}\right). $$
A simple linear combination of Eqs. (9) to (11) then eliminates the term in $n$, showing that $F$ satisfies the equation

$$ \lambda \partial_{\lambda}^{2} F + (2\nu + 1)\partial_{\lambda} F + \lambda z^{2} F = \Delta_{n} \left( \frac{-\lambda A_{n}}{(\lambda^{2} - 1)^{2}} \right). $$

The left-hand side is a telescoper. The right-hand side is a certificate. It can be written more explicitly as

$$ \frac{-\lambda A_n}{(\lambda^2 - 1)^2} = -\frac{4(n\lambda^2 + \lambda^2\nu - \nu - 1)n\lambda}{(\lambda^2 - 1)^2} F - \frac{4\lambda(n+1)n}{(\lambda^2 - 1)^2} S_n(F). $$

In general, summation and telescoping of the certificate requires verification. Here, we first observe that the certificate does not have integer poles and thus is well defined at all points over which it is summed. Next, the certificate evaluates to zero at $n = 0$. Finally, it tends to zero when $n$ tends to infinity, as $J_{\nu+n}(z)$ decreases fast as $n \to \infty$ [27, 10.19.1].

In summary, we have obtained that the sum $S$ in the right-hand side of Eq. (3) satisfies

$$ \lambda \partial_{\lambda}^{2}(S) + (2\nu + 1)\partial_{\lambda}(S) + \lambda z^{2} S = 0. $$

Proceeding similarly with Eqs. (4) and (6), one gets the equations

$$ z\lambda S_{\nu}(S) + \partial_{\lambda}(S) = 0, \quad z\partial_{z}(S) - \lambda\partial_{\lambda}(S) - \nu S = 0. $$

Injecting $T = J_{\nu}(\lambda z)/\lambda^{\nu}$ in these equations and using basic equations for $J_{\nu}$ shows that it is a solution of this system too. The proof of the multiplication theorem is concluded by checking the equality of the initial conditions for $T$ and for the sum on the right-hand side of Eq. (3). As $\nu$ is associated to the shift, we need to check initial conditions for any $\nu$ satisfying $0 \le \text{Re}(\nu) < 1$. Indeed, both term of the identity equal $J_{\nu}(1)$ at $z = 1, \lambda = 1$, and $\nu \in [0, 1)$ and both their derivatives with respect to $\lambda$ equal $-J_{\nu+1}(1)$, which proves the identity.

### 3. BACKGROUND

In this section, we recall the basic framework for reducing creative telescoping to the generalized Abramov-Petkovšek decomposition. Most of this section is identical to the differential case [6, Sec. 4], except for the existence and computation of the cyclic vector and the use of the recurrence variant of Lagrange's identity [3]. More gentle introductions to Ore algebras, creative telescoping and their applications can be found in the references [17, 19].

#### 3.1. Telescoping ideal.

**Ore algebras.** Let $\mathbf{k}$ be a field of characteristic 0, $x_0, \dots, x_m$ be variables used to form the fields of rational functions $\mathbb{K} = \mathbf{k}(x_1, \dots, x_m)$ and $\hat{\mathbb{K}} = \mathbb{K}(x_0)$. The *Ore algebra* $\mathbb{A} = \hat{\mathbb{K}}\langle\partial_0, \dots, \partial_m\rangle$ is a polynomial ring over $\hat{\mathbb{K}}$, with $\partial_i\partial_j = \partial_j\partial_i$, and a commutation between the $\partial_i$ and the elements of $\hat{\mathbb{K}}$ ruled by relations

$$ (12) \qquad \partial_i R = \sigma_i(R)\partial_i + \delta_i(R), \quad R \in \hat{\mathbb{K}}, $$

with $\sigma_i$ a ring morphism of $\hat{\mathbb{K}}$ and $\delta_i$ a $\sigma_i$-derivation, which means that $\delta_i(ab) = \sigma_i(a)\delta_i(b) + \delta_i(a)b$ for all $a, b$ in $\hat{\mathbb{K}}$ [9, 19]. The typical cases are when $\partial_i$ is the differentiation $d/dx_i$ (then $\sigma_i$ is the identity and $\delta_i = d/dx_i$) and the shift operator $x_i \mapsto x_i + 1$ (then $\sigma_i(a) = a|_{x_i \leftarrow x_i+1}$ and $\delta_i = 0$).

**Annihilating and D-finite ideals.** For a given function $f$ in a left $\mathbb{A}$-module, the annihilating ideal of $f$ is the left ideal $\text{ann}\,f \subseteq \mathbb{A}$ of elements of $\mathbb{A}$ that annihilate $f$. A left ideal $\mathcal{I}$ of $\mathbb{A}$ is D-finite when the quotient $\mathbb{A}/\mathcal{I}$ is a finite dimensional $\mathbb{K}$-vector space. A function is called D-finite when its annihilating ideal is D-finite.

**Telescoping ideal.** As we focus here on summation, from now on, when we use $n$ and $S_n$, they stand for $x_0$ and the corresponding shift operator $\partial_0 : x_0 \mapsto x_0 + 1$.

The telescoping ideal $\mathcal{T}_{\mathcal{I}}$ of the left ideal $\mathcal{I} \subset \mathbb{A}$ with respect to $n$ is

$$ \mathcal{T}_{\mathcal{I}} = (\mathcal{I} + \Delta_n(\mathbb{A})) \cap \mathbb{K}\langle\partial_1, \dots, \partial_m\rangle, \quad \text{where} \quad \Delta_n = S_n - 1. $$

In other words, if $\mathcal{I} = \text{ann}\,F$, the telescoping ideal $\mathcal{T}_{\mathcal{I}}$ is the set of operators $T \in \mathbb{K}\langle\partial_1, \dots, \partial_m\rangle$ such that there exists $G \in \mathbb{A}$ such that $T + \Delta_n G \in \mathcal{I}$, or equivalently, such that Eq. (1) holds (with $t = n$).

#### 3.2. Cyclic vector and Lagrange identity.
**Cyclic vector.** Let $\mathcal{I}$ be a D-finite ideal of $\mathbb{A}$ and let $r$ be the dimension of the $\hat{\mathbb{K}}$-vector space $\mathbb{B} := \mathbb{A}/\mathcal{I}$. An element $\gamma \in \mathbb{B}$ is called cyclic with respect to $\partial_0$ if $\{\gamma, \dots, \partial_0^{r-1}\gamma\}$ is a basis of $\mathbb{B}$. In the differential case $(\partial_0 = d/dx_0)$, such a vector always exists and can be computed efficiently when $\mathcal{I}$ is D-finite [15]. In the shift case $(\partial_0 : x_0 \mapsto x_0 + 1)$, even for a D-finite ideal $\mathcal{I}$, it is not the case that there always exists a cyclic vector: in general, $\mathbb{B}$ decomposes as the sum of a vector space where $\partial_0$ is nilpotent and a part where it is cyclic [24]. However, we have the following.

**Proposition 1.** [23, Thm. B2] With the notation above, in the case when $\partial_0$ is the shift operator $x_0 \mapsto x_0+1$, let $E = (e_1, \dots, e_r)^\top$ be a basis of the vector space $\mathbb{B} = \mathbb{A}/\mathcal{I}$ and $A_0 \in \hat{\mathbb{K}}^{r \times r}$ be defined by $\partial_0 E = A_0 E$. If $A_0$ is invertible, then there exists a cyclic vector with respect to $\partial_0$ of the form $v = a_1 e_1 + \dots + a_r e_r$ with polynomial coefficients $a_i \in \mathbb{Z}[x_0]$ of degree at most $r-1$, and coefficients all in $\{0, \dots, r\}$.

Sufficient conditions for the matrix $A_0$ to be invertible are that $\mathcal{I} = \text{ann} f$ with $f$ in a $\hat{\mathbb{K}}[\partial_0, \partial_0^{-1}]$-module [23] or that $\mathcal{I}$ be a reflexive ideal [32]. In practice, this condition on $A_0$ can be checked from the input and appears to be always satisfied in the examples we have tried. From this proposition, the computation of a cyclic vector follows the same lines as that of the differential case [15]. Most often, $e_1 = 1$ is a cyclic vector, which simplifies the rest of the computation.

**Lagrange's identity.** For our purpose, the shift version of Lagrange's identity can be viewed as giving an explicit form of the result of the left Euclidean division by the difference operator $\Delta_n$, when applied to a left multiplication by a rational function.

**Proposition 2.** [3] Let $u \in \hat{\mathbb{K}}$ and let $L = \sum_{i=0}^r a_i S_n^i$ be an operator of order $r$ with $a_i$ in $\hat{\mathbb{K}}$. The adjoint operator $L^*$ of $L$ is defined as $L^* = \sum_{i=0}^r a_i (n-i) S_n^{-i}$ and it satisfies

$$ (13) \qquad uL - L^*(u) = \Delta_n P_L(u) $$

where

$$ (14) \qquad P_L(u(n)) = \sum_{i=0}^{r-1} \left( \sum_{j=i+1}^{r} a_j (n+i-j)u(n+i-j) \right) S_n^i. $$

Note that the term $L^*(u)$ denotes the evaluation of the operator $L^*$ at the rational function $u$, rather than the product of $L^*$ by $u$.

Let $\gamma$ be a cyclic vector. Then any element of $\mathbb{B}$ is of the form $A\gamma$ with $A \in \hat{\mathbb{K}}\langle S_n \rangle$. Applying Lagrange's identity with $u = 1$, $L = A$ and multiplying on the right by $\gamma$ shows that this is a rational multiple of $\gamma$ up to a difference:

$$ (15) \qquad A\gamma = A^*(1)\gamma + \Delta_n P_A(1)\gamma. $$

As in the differential case, all computations in $\mathbb{B}$ then reduce to $\hat{\mathbb{K}}$-linear operations on single rational functions, rather than vectors of them, by the following analogue of [6, Prop. 4.2].

**Proposition 3.** With the notation above, let $\gamma$ be a cyclic vector of $\mathbb{B} = \mathbb{A}/\mathcal{I}$ and for all $i = 0, \dots, m$, let $B_i \in \hat{\mathbb{K}}\langle S_n \rangle$ be such that $\partial_i \gamma = B_i \gamma$. Then for all $R \in \hat{\mathbb{K}}$,

$$ (16) \qquad \partial_i R \gamma = \varphi_i(R) \gamma + \Delta_n Q_i(R) \gamma, \\ \text{with} \quad \begin{cases} \varphi_i(R) = B_i^*(R(x_i+1)), & Q_i(R) = P_{B_i}(R(x_i+1)) & \text{if } \partial_i : x_i \mapsto x_i+1; \\ \varphi_i(R) = B_i^*(R) + \frac{d}{dx_i}(R), & Q_i(R) = P_{B_i}(R(x_i)) & \text{if } \partial_i = d/dx_i. \end{cases} $$

*Proof.* Multiplying Eq. (12) by $\gamma$ on the right and using the definition of $B_i$ gives

$$ \partial_i R \gamma = \sigma_i(R) B_i \gamma + \delta_i(R) \gamma. $$

The conclusion follows from Lagrange's identity (13) applied with $L = B_i$ and $u = \sigma_i(R)$. $\square$
**Algorithm 1** Creative Telescoping Algorithm

**Input** Generators of a D-finite ideal $\mathcal{I}$

**Output** A basis of the telescoping ideal $\mathcal{T}_{\mathcal{I}}$ of $\mathcal{I}$

$\gamma \leftarrow$ a cyclic vector of $\mathcal{B} := \mathbb{A}/\mathcal{I}$ with respect to $S_n$ \quad \{See Section 3.2\}
$L \leftarrow$ a minimal order operator in $S_n$ annihilating $\gamma$ \quad \{See Section 3.3\}
$\varphi_1, \dots, \varphi_n$ the maps described in Proposition 3
Initialize CanonicalForm \quad \{See Section 4\}
$R_0 \leftarrow$ CanonicalForm$(A_1^*(1), L^*)$ \quad \{See Section 3.4\}
$\mathcal{L} \leftarrow$ [1] \quad \{List of monomials in $\partial_1, \dots, \partial_m$ to visit\}
$\mathcal{G} \leftarrow$ {} \quad \{Gröbner basis\}
$\mathcal{Q} \leftarrow$ {} \quad \{Basis of the quotient\}
$\mathcal{R} \leftarrow$ {} \quad \{Set of reducible monomials\}

while $\mathcal{L} \neq \emptyset$ do

Remove the first element $\partial^{\alpha}$ of $\mathcal{L}$

if $\partial^{\alpha}$ is not a multiple of an element of $\mathcal{R}$ then

if $\partial^{\alpha} \neq 1$ then

Take $i$ such that $\partial^{\alpha}/\partial_i \in \mathcal{Q}$; $R_{\alpha} \leftarrow$ CanonicalForm($\varphi_j(R_{\partial^{\alpha}/\partial_i}), L^*$)

if there exists a $\mathbb{K}$-linear relation between $R_{\alpha}$ and $\{R_{\beta} \mid R_{\beta} \in \mathcal{Q}\}$ then

($\lambda_{\beta})_{R_{\beta} \in \mathcal{Q}} \leftarrow$ coefficients of the relation $R_{\alpha} = \sum_{R_{\beta} \in \mathcal{Q}} \lambda_{\beta} R_{\beta}$

Add $\partial^{\alpha} - \sum_{R_{\beta} \in \mathcal{Q}} \lambda_{\beta} \partial^{\beta}$ to $\mathcal{G}$; Add $\partial^{\alpha}$ to $\mathcal{R}$

else

Add $\partial^{\alpha}$ to $\mathcal{Q}$

for $j = 1$ to $m$ do Append the monomial $\partial_j \partial^{\alpha}$ to $\mathcal{L}$

return $\mathcal{G}$

**3.3. Canonical Form.** Proposition 3 shows how, given a cyclic vector $\gamma$, all elements of $\mathbb{B}$ can be reduced to the product of $\gamma$ by a rational function, up to a difference in $\Delta_n \mathbb{B}$. The starting point of the reduction-based creative telescoping is that one can actually identify those multiples that belong to $\Delta_n \mathbb{B}$.

**Proposition 4.** *With the same hypotheses as in Proposition 3, let $L$ be a minimal-order operator in $\hat{\mathbb{K}}\langle S_n \rangle$ annihilating $\gamma$, ie, the product $L\gamma$ is 0 in $\mathbb{A}/\mathcal{I}$, $L$ has order $r$ and no operator of order $r-1$ has that property. Then for all $R \in \hat{\mathbb{K}}$, $R\gamma \in \Delta_n(\mathbb{B}) \iff R \in L^*(\hat{\mathbb{K}})$.*

*Proof.* First, if $R = L^*(R')$ with $R' \in \hat{\mathbb{K}}$, Lagrange's identity (13) with $u = R'$ and $L\gamma = 0$ implies that $R\gamma = L^*(R')\gamma = \Delta_n G\gamma$ for $G = -P_L(R')$. Conversely, if $R\gamma \in \Delta_n(\mathbb{B})$, there exists $M \in \hat{\mathbb{K}}\langle S_n \rangle$ such that $R\gamma = \Delta_n M\gamma$. The operator $\Delta_n M - R$ annihilates $\gamma$. By minimality of $L$, there exists $N \in \hat{\mathbb{K}}\langle S_n \rangle$ such that $\Delta_n M - R = NL$. Taking the adjoint and evaluating at 1 gives $R = M^*\Delta_n^* - L^*N^*$ and finally $R = -L^*(N^*(1))$. $\square$

This proposition motivates the following.

**Definition 1.** [6] A *canonical form* associated to $L^*$ is a $\mathbb{K}$-linear map $[\cdot] : \mathbb{K}(n) \to \mathbb{K}(n)$ such that for all $R \in \mathbb{K}(n)$, $[L^*(R)] = 0$ and $R - [R] \in L^*(\mathbb{K}(n))$. A rational function $R \in \mathbb{K}(n)$ is called *reduced* when $[R] = R$.

The computation of canonical forms is the object of Section 4.

**3.4. Creative Telescoping Algorithm via Canonical Forms.** With the notation above, the creative telescoping algorithm from [6] applies verbatim. It is given in Algorithm 1. Its principle is to iterate on every monomial of the form $\partial_1^{\alpha_1} \dots \partial_m^{\alpha_m}$ by increasing order for some monomial order, e.g., the grevlex order, and to compute the reduced rational functions $R_{\alpha}$ such that

$$ (17) \qquad \partial^{\alpha} = R_{\alpha}\gamma + \Delta_n G_{\alpha} \mod \mathcal{I}. $$
The rational function $R_\alpha$ is obtained by

$$ (18) \qquad \begin{cases} R_{0,...,0} = [A_1^*(1)], \\ R_\alpha = [\varphi_i(R_\beta)] & \text{if } \partial^\alpha = \partial_i \partial^\beta, \end{cases} $$

where $A_1^*$ is the adjoint of the operator $A_1$ verifying $1 = A_1(\gamma)$. When a monomial $\alpha$ is dealt with, two situations are possible. The corresponding $R_\alpha$ can be a linear combination of the previous $R_\beta$. In that case, that linear combination makes the corresponding linear combination of $\partial^\alpha$ and the $\partial^\beta$ a newly discovered element of the telescoping ideal $\mathcal{T}_I$ and then it is not necessary to visit the multiples of this monomial. Otherwise, $\partial^\alpha$ is free from the previous ones and thus a new generator of $\mathbb{B}$ has been found. The algorithm terminates when there are no more monomials to visit.

The only difference with the differential case lies in the definition of the canonical form $[\cdot]$ associated to the adjoint $L^*$ of the minimal-order operator $L \in \mathbb{K}(n)\langle S_n \rangle$ annihilating $\gamma$. By Propositions 3 and 4 and the definition of canonical form, Eq. (17) is satisfied and the following equivalence holds:

$$ (19) \qquad a_1 R_1 + \dots + a_s R_s = 0 \quad \text{iff} \quad a_1 D^1 + \dots + a_s D^s \in \mathcal{T}_I. $$

The following result follows, with the same proof as in [6].

**Theorem 1.** Given as input the generators of a D-finite ideal $\mathcal{I}$ and a cyclic vector $\gamma$ for $S_n$, Algorithm 1 terminates if and only if $\mathcal{T}_I$ is D-finite. Then, it outputs a Gröbner basis of $\mathcal{T}_I$ for the grevlex order.

As in the differential situation [6], one can modify the algorithm to compute all elements of $\mathcal{T}_I$ up to a given bound on their degree, or to return as soon as one telescoper is found, thus allowing to recover a generating family of a sub-ideal of $\mathcal{T}_F$.

### 4. GENERALIZED ABRAMOV-PETKOVŠEK DECOMPOSITION

The main contribution of this article is an algorithm for the computation of canonical forms as in Definition 1 for the operator

$$ L^* = \sum_{i=0}^{r} p_i(n) S_n^{-i} $$

with polynomial coefficients $p_i \in \mathbb{K}[n]$. The modified Abramov-Petkovšek decomposition [12] is a special case of this reduction when $L$ has order 1 and once the shell [12] has been removed [6, Sec. 3.5.3].

The starting point is a decomposition of any rational function $R \in \mathbb{K}(n)$ in the form

$$ (20) \qquad R(n) = P_\infty(n) + \sum_{i,h} \frac{c_{i,h}(n)}{Q_i(n+h)^{\ell_{i,h}}}, $$

with $\ell_{i,h} \in \mathbb{N}^*$, polynomials $P_\infty, Q_i$ and $c_{i,h}$ in $\mathbb{K}[n]$ such that $\deg c_{i,h} < \ell_{i,h} \deg Q_i$ and $\gcd(Q_i(n+k), Q_j(n)) = 1$ for all $k \in \mathbb{Z}$ when $i \neq j$. This is discussed in Section 4.1.

The vector spaces

$$ \mathbb{K}_{Q_i}(n) \stackrel{\text{def}}{=} \text{Vect}_{\mathbb{K}} \left( \frac{n^\ell}{Q_i(n+h)^j} \mid h \in \mathbb{Z}, j \in \mathbb{N}^*, \ell < j \deg(Q_i) \right) $$

are in direct sum for distinct $Q_i$ and are left invariant by $L^*$ modulo $\mathbb{K}[n]$. This allows the reduction algorithm to operate in each of the $\mathbb{K}_{Q_i}(n)$ independently. This is described in Sections 4.2 and 4.3, before the reduction of the remaining polynomial part in Sections 4.4 and 4.5.

**Notation 1.** For two integers $a, b$ with $a \le b$, we write $\llbracket a; b \rrbracket$ for the set $\{a, a+1, \dots, b\}$.

4.1. **Decomposition of rational functions.** Recall that a polynomial $Q$ is square-free when it does not have multiple nontrivial factors. It is shift-free when $\gcd(Q(n), Q(n+k)) = 1$ for all $k \in \mathbb{Z}^*$.

A shiftless decomposition of a polynomial $Q$ is a factorization of the form

$$ Q = \prod_{i=1}^{v} \prod_{j=1}^{n_i} Q_i(n + h_{i,j})^{e_{i,j}}, $$
**Algorithm 2** Weak reduction of poles [·]Q

**Input** R, Q

**Output** a reduced form of R

while there exists $j < 0$ (minimal) such that $Q(n-j) \mid \text{den}(R)$ **do**

$$R \leftarrow R - L^* \left( \frac{A(n)}{Q(n-j)^{\text{ord}_j(\text{den}(R)) + \text{ord}_j(p_0)}} \right) \text{ with } A \text{ as in Eq. (22)}$$

while there exists $j \ge r$ (maximal) such that $Q(n-j) \mid \text{den}(R)$ **do**

$$R \leftarrow R - L^* \left( \frac{A'(n+r)}{Q(n-j+r)^{\text{ord}_j(\text{den}(R)) + \text{ord}_j(p_r)}} \right) \text{ with } A' \text{ as in Eq. (23)}$$

**return** R

with $e_{i,j} \in \mathbb{N}^*$, $h_{i,j} \in \mathbb{Z}$, and $Q_i \in \mathbb{K}[n]$ are such that each $Q_i$ is square-free and $\text{gcd}(Q_i(n+k), Q_j(n)) = 1$ for all $i, j$ and all $k \in \mathbb{Z}$ unless $i = j$ and $k = 0$. Such a factorization can be computed using only gcds, resultants and integer root finding [21].

Note that shiftless decompositions are not unique in general. One can be refined when a $Q_i$ is not irreducible, by splitting this factor further. In particular, the linear factors of the $Q_i$ can be isolated and dealt with more easily.

A polynomial $Q$ is refined with respect to a polynomial $P$ when it is such that for each $h \in \mathbb{Z}$, there exists $\ell \in \mathbb{N}$ such that $\text{gcd}(P, Q(n+h)^{\ell+1}) = Q(n+h)^{\ell}$. A shiftless decomposition is called refined with respect to $P$ when each $Q_i$ is. This refinement can be computed using gcds only and will be used with $P = p_0$ and $P = p_r$, the extreme coefficients of $L^*$.

From a shiftless decomposition, the partial fraction decomposition of Eq. (20) is then obtained by standard algorithms [34, 5.11].

## 4.2. Weak reduction of the polar part.

**Lemma 1.** Let $Q \in \mathbb{K}[n]$ be square-free, shift-free and refined with respect to the coefficients $p_0$ and $p_r$ of $L^*$. Given a rational function $R \in \mathbb{K}_Q(n)$, Algorithm 2 computes a rational function $[R]_Q \in \mathbb{K}_Q(n)$ with all its poles at zeros of $Q(n-j)$ such that $j \in [0; r-1]$ and $R - [R]_Q = P + L^*(T)$ for some $P \in \mathbb{K}[n]$ and $T \in \mathbb{K}_Q(n)$. The algorithm is $\mathbb{K}$-linear.

*Proof.* Assume that $R$ decomposes as

$$(21) \qquad R = \sum_{j \in J} \frac{\lambda_j(n)}{Q(n-j)^{s_j}} \quad \text{with } \deg(\lambda_j(n)) < s_j \deg(Q).$$

Let $j_m = \min(J)$ and $\text{ord}_j(p_0)$ be the largest integer $\ell$ such that $Q(n-j)^\ell \mid p_0$. Then,

$$Q(n - j_m)^{s_{jm}} L^* \left( \frac{1}{Q(n - j_m)^{s_{jm} + \text{ord}_{jm}(p_0)}} \right) = \tilde{p}_0(n) \mod Q(n - j_m)^{s_{jm}},$$

where $\tilde{p}_0(n)$ is the remainder in the Euclidean division of $p_0/Q(n-j_m)^{\text{ord}_{jm}(p_0)}$ by $Q(n-j_m)^{s_{jm}}$. The poles of this rational function are at zeros of $Q(n-j)$ with $j \in \tilde{J} := J \setminus \{j_m\} \cup (j_m + [1, r-1])$.

Since $Q$ is reduced with respect to $p_0$, the polynomial $\tilde{p}_0(n)$ is relatively prime with $Q(n-j_m)$. Thus, there exist polynomials $A$ and $B$ such that

$$(22) \qquad \lambda_{j_m}(n) = A(n)\tilde{p}_0(n) + B(n)Q(n - j_m)^{s_{jm}}.$$

Then

$$\frac{A(n)\tilde{p}_0(n)}{Q(n - j_m)^{s_{jm}}} = \frac{\lambda_{j_m}(n)}{Q(n - j_m)^{s_{jm}}} - B(n),$$

so that

$$R - L^* \left( \frac{A(n)}{Q(n - j_m)^{s_{jm} + \text{ord}_{jm}(p_0)}} \right)$$

is equivalent to $R$ modulo $L^*(\mathbb{K}_Q(n))$ and with all its poles at zeros of $Q(n-j)$ with $j \in \tilde{J}$. This operation can be repeated a finite number of times until all poles are at zeros of $Q(n-j)$ with $j \ge 0$.
Similarly, let $j_M = \max(J)$ and $\mathrm{ord}_{j_M+r}(p_r)$ be the largest integer $\ell$ such that $Q(n-j_M)^\ell$ divides $p_r(n-r)$. Then

$$Q(n - j_M)^{s_{j_M}} L^* \left( \frac{1}{Q(n - j_M + r)^{s_{j_M} + \mathrm{ord}_{j_M+r}(p_r)}} \right) = \tilde{p}_r(n) \mod Q(n - j_M)^{s_{j_M}},$$

where $\tilde{p}_r(n)$ is the remainder in the Euclidean division of $p_r/Q(n-j_M)^{\mathrm{ord}_{j_M+r}(p_r)}$ by $Q(n-j_M)^{s_{j_M}}$. The poles of this rational function are at zeros of $Q(n-j)$ with $j \in \hat{J}' := J \setminus \{j_M\} \cup (j_M - [1, r-1])$. Again, since $Q$ is reduced with respect to $p_r$, the polynomial $\tilde{p}_r$ is relatively prime with $Q(n-j_M)$. Thus there exist two polynomials $A'$ and $B'$ such that

$$(23) \qquad \lambda_{j_M}(n) = A'(n)\tilde{p}_r(n) + B'(n)Q(n-j_M)^{s_{j_M}}$$

so that

$$R - L^* \left( \frac{A'(n+r)}{Q(n - j_M + r)^{s_{j_M} + \mathrm{ord}_{j_M+r}(p_r)}} \right)$$

is equivalent to $R$ modulo $L^*(\mathbb{K}_Q(n))$ and with all its poles at zeros of $Q(n-j)$ with $j \in \hat{J}'$. This operation can be repeated a finite number of times until all poles are at zeros of $Q(n-j)$ with $j \in [0, r-1]$.

Each step being $\mathbb{K}$-linear, so is the algorithm. $\square$

**Example 1.** Let

$$R = 8nx + \frac{n(8x^2 + x) - 16x - 1}{2(n-1)^2} - \frac{4(x-1)x^2}{n} - \frac{(4x^3 + 8x^2 - 31x - 32)n + 4x^3 - 31x - 32}{2(n+1)^2} + \frac{4x(x^2 - x - 8)}{n+2} + \frac{2x^3}{n+3}$$

and

$$L^* = x^2(n-2)S_n^{-3} - n(4n^2 - x^2 - 4n)S_n^{-2} + n(4n^2 - x^2 - 4n)S_n^{-1} - x^2(n+2).$$

The poles of $R$ are at $\{1, 0, -1, -2, -3\}$. We take $Q = n+1$ and follow the steps of the algorithm.

The pole at $-3$ is easy: from

$$L^* \left( \frac{1}{n+3} \right) = -x^2 + 4n + \frac{-x^2 + 8}{n+1} + \frac{2x^2 - 16}{n+2} + \frac{x^2}{n+3}$$

and the coefficient $2x^3$ of $(n+3)^{-1}$ in $R$, the algorithm performs the subtraction

$$R \leftarrow R - 2xL^* \left( \frac{1}{n+3} \right) = \frac{n(8x^2 + x) - 16x - 1}{2(n-1)^2} + \frac{4x^2}{n} - \frac{(8x^2 + x - 32)n + x - 32}{2(n+1)^2} - \frac{4x^2}{n+2}.$$

Next, the pole $-2$ is a simple root of the constant coefficient of $L^*$, leading to the computation of

$$L^* \left( \frac{1}{(n+2)^2} \right) = \frac{x^2(n-2)}{(n-1)^2} + \frac{x^2}{n} - \frac{n(x^2-4)-4}{(n+1)^2} - \frac{x^2}{n+2}$$

so that the pole is removed by

$$(24) \qquad R \leftarrow R - 4L^* \left( \frac{1}{(n+2)^2} \right) = \frac{x}{2(n-1)} - \frac{x}{2(n+1)}.$$

R now has all its poles in $\{-1, 0, 1\}$ and the weak reduction is finished.

4.3. **Strong reduction of the polar part.** By Corollary 1, the weak reduction produces rational functions all whose poles differ from those of $Q$ by an integer in $[0, r-1]$. The next step of the reduction is to subtract rational functions in $L^*(\mathbb{K}_Q(n))$ that have this property.

It turns out to be possible to focus on a finite-dimensional subspace of $L^*(\mathbb{K}_Q(n))$ thanks to the following.

**Lemma 2.** If $j < 0$, $s > \mathrm{ord}_j(p_0)$ and $\ell < s \deg(Q)$ or if $j \ge 0$, $s > \mathrm{ord}_{j+r}(p_r)$ and $\ell < s \deg(Q)$ then

$$\left[ L^* \left( \frac{n^\ell}{Q(n-j)^s} \right) \right]_Q = 0.$$
**Algorithm 3** Weak reduction of polynomials [·]∞

**Input** P and (σ,p) from Eq. (25)

**Output** a reduced form of P

a ← 0

while deg(P) ≥ σ do

if deg(P) − σ is a root of p then a ← a + lt(P); P ← P − lt(P)

else P ← P − (lc(P) / p(deg(P) − σ)) L $^{∗}$(n$^{deg}$ ( P ) − σ )

return a + P

Proof. Let j, s, ℓ be three integers that satisfy the first assumption. Then L $^{∗}$( n $^{ℓ}$/Q ( n − j ) $^{s}$) has a denomi- nator that is divisible by Q ( n − j ) with j < 0 by assumption. No smaller k is such that Q ( n − k ) divides the denominator. Thus the first pass through the first loop of the weak reduction subtracts L $^{∗}$( n $^{ℓ}$/Q ( n − j ) $^{s}$) to itself and reduces it to zero. When the second assumption is satisfied, then L $^{∗}$( n $^{ℓ}$/Q ( n − j ) $^{s}$) has a denominator that is divisible by Q ( n − ( j + r )) with j + r ≥ r by assumption. No larger k is such that Q ( n − k ) divides the denominator. Thus again, the second loop reduces that fraction to 0. □

**Corollary 1.** Let I$_{0}$ := { j ∈ Z < 0 | gcd( p$_{0}$ ( n ) ,Q ( n − j )) ̸ = 1 } and I$_{r}$ := { j ∈ Z ≥ 0 | gcd( p$_{r}$ ( n ) ,Q ( n − j )) ̸ = 1 } . The K -vector space [ L $^{∗}$( K $_{Q}$( n ))]$_{Q}$ is generated by the fractions

$$
\left[ L^* \left( \frac{n^\ell}{Q(n-j)^{s_j}} \right) \right]_Q, \quad \text{with} \quad \begin{cases} j \in I_0 \text{ and } 1 \le s_j \le \text{ord}_j(p_0) \text{ and } 0 \le \ell < s_j \deg Q \\ \text{or} \\ j \in I_r \text{ and } 1 \le s_j \le \text{ord}_{j+r}(p_r) \text{ and } 0 \le \ell < s_j \deg Q. \end{cases}
$$

Corollary 1 gives a generating family of the finite dimensional K-vector space [ L $^{∗}$( K $_{Q}$( n ))]$_{Q}$ . These rational functions can be written in the basis ( n $^{i}$/Q ( n − j ) $^{k}$)$_{i,k}$$_{∈}$$_{N}$ ,$_{j}$$_{∈}$$_{Z}$ and one can then compute an echelon basis of this finite-dimensional space. This precomputation step corresponds to the computation of the B$_{Q}$$_{i}$ ’s in Algorithm 4. The strong reduction of a rational function R ∈ K$_{Q}$ then consists in reducing [ R ]$_{Q}$ with this echelon basis. By this process, we obtain the following.

**Proposition 5.** *Strong reduction reduces every rational function R ∈ L $^{∗}$( K $_{Q}$( n )) to a polynomial in K [ n ].*

**Example 2.** With the same notation as in example 1, Corollary 1 shows that [ L $^{∗}$( K $_{n}$+1 ( n ))]$_{n}$$_{+1}$ is generated by

$$
\left[ L^* \left( \frac{1}{n-2} \right) \right]_{n+1} = 4n + \frac{x^2}{n+1} - \frac{x^2}{n-1}, \quad \left[ L^* \left( \frac{1}{n-1} \right) \right]_{n+1} = 4n + \frac{x^2}{n-1} - \frac{x^2}{n+1}.
$$

Thus the strong reduction of the rational function R from Eq. (24) is the polynomial

$$
R + L^{*}((n-2)^{-1})/(2x) = -2n/x,
$$

concluding the reduction.

4.4. Weak reduction of polynomials. The weak reduction of polynomials is a direct adaptation of the differential case [6]. The indicial polynomial of L* at infinity is the polynomial p ∈ K[s] defined by

$$
(25) \qquad L^*(n^s) = n^{s+\sigma}(p(s) + O(1/n)),
$$

with σ ∈ N . The ensuing weak reduction is presented in Algorithm 3.

**Example 3.** In Examples 1 and 2, the indicial equation at infinity is

$$
L^*(n^s) = n^{s+2}(8 + 4s + O(1/n)).
$$

The polynomial − 2 n/x found in Example 2 cannot be reduced further by weak reduction since its degree in n is smaller than 2.

The following properties are proved exactly as those for weak reduction at a pole.

**Lemma 3.** Algorithm [·]$_{∞}$ terminates and is K -linear. For all P ∈ K [ n ], there exists Q ∈ K [ n ] such that P − [ P ] = L $^{∗}$( Q ) . If s ∈ N is not a root of p , then [ L $^{∗}$( n $^{s}$)]$_{∞}$ = 0 .
**Algorithm 4 PrecomputeBases**

**Input** $Q_1, \dots, Q_v$ polynomials that occur in the shiftless decomposition of $p_0p_r$

**Output** The echelon bases $B_{Q_1}, \dots, B_{Q_v}, B_{pol}$

$$B_{pol} \leftarrow \{\}$$

**for** $i = 1$ **to** $v$ **do**

$$B_{Q_i} \leftarrow \text{Echelon}\left(\left\{\left[L^*\left(\frac{n^l}{Q(n-j)^s}\right)\right]_{Q_i} \mid l, j, s_j \text{ as in Cor.1}\right\}\right)$$

$$B_{pol} \leftarrow B_{pol} \cup (B_{Q_i} \cap \mathbb{K}[n])$$

$$B_{pol} \leftarrow \text{Echelon}(B_{pol} \cup \{[L^*(n^s)]_\infty \mid s \text{ integer root of } p\})$$

**return** $B_{Q_1}, \dots, B_{Q_v}, B_{pol}$

**Algorithm 5 Reduction of rational functions [·]**

**Input** $R$ and $B_{Q_1}, \dots, B_{Q_v}, B_{pol}$ computed by Algorithm 4

**Output** a reduced form of $R$ by $L^*(\mathbb{K}(n))$.

Decompose $R$ as $P_\infty + \sum_{i=1}^v R_i$ with $R_i \in \mathbb{K}_{Q_i}$

$$\text{for } i = 1 \text{ to } v \text{ do} \quad R_{Q_i} \leftarrow \text{StrongReduce}([R_i]_{Q_i}, B_{Q_i})$$

$$R \leftarrow P_\infty + \sum_{i=1}^s R_{Q_i}$$

Write $R = P + \tilde{R}$ with $P$ a polynomial and $\deg(\tilde{R}) < 0$

$$P \leftarrow \text{StrongReduce}([P]_\infty, B_{pol})$$

**return** $P + \tilde{R}$

4.5. **Strong reduction of polynomials.** The final step is to subtract polynomials in $L^*(\mathbb{K}(n))$. Here again, a finite number of generators can be obtained thanks to the following.

**Lemma 4.** Let $Q_1, \dots, Q_v$ be the polynomials that occur in a shiftless decomposition of $p_0p_r$ and let $P$ be a polynomial in $L^*(\mathbb{K}(n))$. Then

$$P \in E_{pol} \stackrel{\text{def}}{=} L^*(\mathbb{K}[n]) + \sum_{i=1}^{v} [L^*(\mathbb{K}_{Q_i}(n))]_{Q_i} \cap \mathbb{K}[n].$$

*Proof.* If $R \in \mathbb{K}(n)$ is such that $L^*(R)$ is a polynomial, then the poles of $R$ must be cancelled by the zeros of $p_0p_r$ or their shifts. It follows that $R$ decomposes as

$$R = R_\infty + \sum_{i=1}^{v} R_i$$

with $R_i \in \mathbb{K}_{Q_i}(n)$ and $R_\infty \in \mathbb{K}[n]$. Each $L^*(R_i)$ has to be a polynomial and thus invariant by $[\cdot]_{Q_i}$. This concludes the proof. $\square$

By Lemma 4, the vector space $[L^*(\mathbb{K}[n])]_\infty$ is generated by $\{[L^*(n^s)]_\infty \mid s \in \mathbb{N}$ and $p(s) = 0\}$ with $p$ the indicial polynomial of $L^*$ at infinity. Generators of each $[L^*(\mathbb{K}_{Q_i})]_{Q_i} \cap \mathbb{K}[n]$ are obtained from the echelon basis used in the strong reduction with respect to $Q_i$. This gives a finite set of generators for $[E_{pol}]_\infty$, which is easily transformed into a basis by a row echelon computation. Strong reduction consists in reducing modulo this basis. The following consequence is as in the polar case.

**Lemma 5.** The strong reduction of polynomials reduces every polynomial $P \in L^*(\mathbb{K}(n))$ to zero.

**Example 4.** Continuing Examples 1 to 3, the polynomial $p(s) = 8+4s$ has no positive integer root therefore $[L^*(\mathbb{K}[n])]_\infty = \{0\}$. A basis of $[L^*(\mathbb{K}_{n+1})]_{n+1} \cap \mathbb{K}[n]$ is $\{n\}$ according to example 2. Therefore $2n/x$ reduces to 0.

4.6. **Canonical Form.** Algorithm 4 and Algorithm 5 combine the previous algorithms to produce a canonical form.

**Theorem 2.** Algorithm 5 computes a canonical form.

*Proof.* Algorithm 5 is linear as every step is linear. By Proposition 5 and Lemma 5, $[L^*(\mathbb{K}(n))]$ reduces to 0 and $[R] - R \in L^*(\mathbb{K}(n))$ as only functions in this image were subtracted to $R$. $\square$
## 5. CERTIFICATES

Reduction-based creative telescoping algorithms allow to find a telescoper without having to compute an associated certificate. This has led to faster algorithms as certificates are known to be larger than telescopers [8]. This approach makes sense in the differential case when it is known in advance that the integral of a certificate over a cycle that avoids singularities is equal to zero. The framework is not as favorable for sums. Indeed, it is necessary to detect whether the certificate has poles in the range of summation and it is often unclear whether the certificate becomes 0 at the boundaries of the summation interval.

It is however possible to compute the certificates in a compact way during the execution of our algorithm, with almost no impact on the execution time. The idea is to make the computation and storage of certificates efficient by storing them as directed acyclic graphs (dags) rather than operators with normalized rational function coefficients. These dags have a number of internal nodes of the same order as the number of operations performed when computing the telescoper, so that their computation does not burden the complexity. They can then be evaluated at the endpoints of the range of summation, or expanded in Laurent series there.

### 5.1. Computation and structure of certificates.
Equations (17) and (18) show how one can compute telescopers without computing certificates. We now show how to compute the certificates simultaneously. These are obtained as sums of monomials in the $\partial_i$ multiplied by rational functions. These rational functions have denominators of the form $Q(n-j)^s$, with integers $j, s$ and polynomials $Q$ that occur either in the shiftless decomposition of $p_0p_r$ or in the denominator of the operator $A_1 \in \hat{\mathbb{K}}\langle S_n \rangle$ such that $1 = A_1\gamma$, or in the denominator of one of the operators $B_i$ from Proposition 3. Those sums are not reduced to a common denominator. They share many common coefficients and denominators that are efficiently compacted into dags by sharing common subexpressions (this is how Maple stores them by default).

The starting point is the cyclic vector $\gamma \in \mathbb{A}/\mathcal{I}$ and the operator $A_1 \in \hat{\mathbb{K}}\langle S_n \rangle$ such that $1 = A_1\gamma$. By Euclidean division by $\Delta_n$ on the left,

$$A_1 = R_0 + \Delta_n g_0.$$

In general, the certificate $G_\alpha$ in Eq. (17) is stored as an (unreduced) element $g_\alpha$ of $\mathbb{A}$ such that $G_\alpha = g_\alpha\gamma \mod \mathcal{I}$. The computation is incremental. It is initialized with $R_0$ and $g_0$ as above, corresponding to $\alpha = 0$. In many cases, the vector 1 is cyclic, so that one can take $\gamma = 1 = A_1 = R_0$ and $g_0 = 0$. Otherwise, the denominators of $R_0$ and $g_0$ are shifts of the denominator of $A_1$.

Next, Eq. (18) leads to the computation of the canonical form of the rational function $\phi_i(R_\beta)$. This computation is performed via a sequence of reductions which consist of subtractions of elements in $L^*(\mathbb{K}(n))$. Keeping track of these rational functions (without normalizing them) gives the canonical form $R_\alpha$ as

$$\phi_i(R_\beta) = R_\alpha - L^*(c_\alpha)$$

for some rational function $c_\alpha \in \hat{\mathbb{K}}$. Both $R_\alpha$ and $c_\alpha$ have denominators that are shifts of those of $R_\beta$ or of factors of $p_0p_r$. In view of Eq. (16), it follows that

$$\partial^{\alpha} = \partial_i \partial^{\beta} = R_{\alpha} \gamma + \Delta_n (\partial_i G_{\beta} + P_{B_i}(\sigma_i(R_{\beta}))\gamma + P_L(c_{\alpha})\gamma) \mod \mathcal{I}.$$

Thus, the certificate $G_\alpha$ is obtained as $g_\alpha\gamma$ with

$$(26) \qquad g_{\alpha} = \partial_i g_{\beta} + P_{B_i}(\sigma_i(R_{\beta})) + P_L(c_{\alpha}).$$

This proves the claim concerning the structure of the certificates and the factors of their denominators. In our implementation, $\partial_i$ is commuted with the coefficients of the certificate $g_\beta$ only. If desired, one can further use $\partial_i = B_i \mod \mathcal{I}$ so as to write $g_\alpha$ as an operator in $S_n$ only.

### 5.2. Evaluation of the certificates.
The output of Algorithm 1 is a set of elements $T$ of the telescoping ideal, which means that

$$T(F) = g(F)(n+1, x_1, \dots, x_m) - g(F)(n, x_1, \dots, x_m),$$

with $g$ a certificate as described above. Summing over $n$, the right-hand side telescopes and only the values of the certificate $g(F)$ at the endpoints are needed.

It is possible to prove that these evaluations are zero without any evaluation in two important cases. First, if the summand $F$ has finite support (e.g., binomial sums), then the sum of any certificate over $\mathbb{Z}$
will be zero provided it has no pole in the summation range. The second case is when one can prove that $R(n)\gamma(F), \dots, R(n)S_n^{r-1}\gamma(F)$ tend to zero as $n$ tends to $\pm\infty$ for any rational function $R \in \mathbb{K}[n]$ (as in the introductory example). Then again the sum of any certificate over $\mathbb{Z}$ will be zero provided it has no pole in the summation range.

**5.3. Integer pole detection.** By its very nature, the method of creative telescoping requires the certificate not to have poles in the range of summation, so that telescoping can occur. The structure of the certificates described above does not allow the efficient computation of its denominator exactly. However it is possible to compute a multiple of it by taking the least common multiple of the denominators of every rational function in the representation. This can be done efficiently by performing the computation on the dag representation of the certificates.

From this multiple of the denominator of the certificate, one can compute the set of roots that lie in the summation range; this amounts to computing the roots that differ from the endpoints of the summation range by an integer. If that set is not empty, then one can compute a Laurent series expansion of the certificate at any point to check whether it is a pole or not, again by exploiting the dag representation of the rational function coefficients.

**5.4. Examples.**

**5.4.1. Neumann's Addition Theorem for Bessel functions.** On input $S(x) = \sum_{n=1}^{\infty} J_n(x)^2$ where $J_n(x)$ is the Bessel function of the first kind, Algorithm 1 outputs the telescoper $\partial_x$ and a certificate $G$ in dag form that exhibits poles at $n \in \{-1, 0, 1, 2\}$. Our implementation produces the polynomial $n(n+1)(n-2)(n-1)$ containing this information. In this example, the certificate is small enough that it can easily be normalized and one gets its value as

$$-\frac{x}{4(n+1)}S_n^2 + \frac{n+1}{x}S_n - \frac{8n^2 - x^2 + 8n}{4x(n+1)},$$

showing that the poles at 1 and 2 vanish in the normalization. Without normalizing the certificate, one can still evaluate the series expansions of the certificate at those points to establish that it has a finite limit there, making the summation legitimate. At $n=1$, the evaluation is found to be

$$-J_0(x)J_1(x) = \frac{1}{2}(J_0(x))^2.$$

Thus we have proved

$$\frac{d}{dx}(\frac{1}{2}J_0^2(x) + J_1(x)^2 + J_2(x)^2 + \dots) = 0.$$

This shows that the sum is constant and the value is revealed by its value at 0, which follows from $J_k(0) = 0$ for $k > 0$ and $J_0(0) = 1$, so that in the end, we recover the classical identity [27, 10.23.3]

$$1 = J_0(x)^2 + 2 \sum_{k \ge 1} J_k(x)^2.$$

**5.4.2. Apéry's Sequence.** The classical sum

$$(27) \qquad A_n = \sum_{k=0}^{n} \binom{n}{k}^2 \binom{n+k}{k}^2,$$

used by Apéry in his proof of the irrationality of $\zeta(3)$, has telescoper

$$(28) \qquad \mathcal{T} := (n+2)^3 S_n^2 - (2n+3)(17n^2 + 51n + 39)S_n + (n+1)^3.$$

The singularities of the certificate obtained by our implementation are at $k \in \{n+1, n+2\}$. Indeed, once normalized, the certificate is found to be

$$C := \frac{4(3k - 4n - 8)k^4}{(k - n - 2)^2}S_n - \frac{4k^4(3k + 4n + 4)}{(k - n - 1)^2}.$$

Let $U_{n,k}$ denote the product of binomials in the sum. Summing $\mathcal{T}U_{n,k}$ from $k=0$ to $k=n+2$ gives $\mathcal{T}(A_n)$. If telescoping is legitimate, then the values of the endpoints are the values of $CU_{n,k}$ at $k=0$ and $k=n+3$, that are both easily checked to be 0. For this to allow to conclude that $\mathcal{T}(A_n) = 0$, it is then sufficient to check that $CU_{n,k}$ is not singular at $k=n+1$ and $k=n+2$, even though $C$ is. Indeed, a series expansion
<table><caption>TABLE 1. Timings. The notation (*) means that we could not check whether the telescopes returned by HF-FCT were minimal.</caption><thead><tr><th></th><th>HF-CT</th><th>HF-FCT</th><th>redctsum</th></tr></thead><tbody><tr><td>easy examples</td><td>6.7s</td><td>7s</td><td>0.9s</td></tr><tr><td>Eq. (30)</td><td>101s</td><td>49s</td><td>0.8s</td></tr><tr><td>Eq. (31)</td><td>52s</td><td>4s</td><td>1.4s</td></tr><tr><td>Eq. (32)</td><td>62s</td><td>1.7s</td><td>5.7s</td></tr><tr><td>Eq. (33)</td><td>4.9s</td><td>1.4s</td><td>10.3s</td></tr><tr><td>Eq. (34)</td><td>4.9s</td><td>1.4s</td><td>13.5s</td></tr><tr><td>Eq. (35)</td><td>1200s</td><td>13s</td><td>205s</td></tr><tr><td>Eq. (36)</td><td>&gt; 6h</td><td>108s</td><td>3338s</td></tr><tr><td>Eq. (37)</td><td>1703s</td><td>4.7s</td><td>580s</td></tr><tr><td>Eq. (38)</td><td>&gt; 1h</td><td>3.2s(*)</td><td>&gt; 1h</td></tr><tr><td>Eq. (39)</td><td>&gt; 1h</td><td>&gt; 1h</td><td>0.4s</td></tr></tbody></table>

of the evaluation of $CU_{n,k}$ at $k = n + 1$ and $k = n + 2$ is possible for our implementation, and finds that the sequence has a finite limit there, which concludes the proof that the telescoper Eq. (28) cancels the sum Eq. (27); see also [18] for more on these issues.

5.5. **A larger example.** The computation of telescoper and certificate for the sum

$$ (29) \qquad \sum_{n=0}^{\infty} \frac{(4n+1)(2n)!}{n!^2 2^{2n} \sqrt{x}} J_{2n+1/2}(x) P_{2n}(u) $$

takes less than 15 sec. with our current implementation (see Table 1). The telescopes are quite small:

$$ (1 - u^2)\partial_u + xu\partial_x, \quad (u^3 - u)\partial_u^2 + (1 + u^2)\partial_u - u^3 x^2. $$

In this example, not normalizing the certificates during their computation has a cost. The actual certificates, once reduced by the Gröbner basis of the annihilating ideal of the summand, are not very large. They are easily computed by Koutschan's program. Still, the corresponding dags are large. Nonetheless, it takes less than 1 sec. to compute a multiple of the denominators of the certificates and detect that they do not have integer roots. Evaluating the certificates at $n = 0$ using their dag representation takes less than 2 min. and proves that the telescopes cancel the sum in Eq. (29).

## 6. IMPLEMENTATION

This algorithm is implemented in Maple<fnref n="2" />. Table 1 gives a comparison of our code with Koutschan's heuristic (HF-FCT) and Chyzak's algorithm (HF-CT)<fnref n="3" />. They are both implemented in Koutschan's Holo-nomicFunctions package in Mathematica [25]. The column 'redctsum' corresponds to our algorithm.

These programs have been executed on a list of 21 easy examples that were compiled by Koutschan, as well as more difficult ones given in Eqs. (30) to (39) below. Eq. (30) comes from recent identities involving determinants [2], Eqs. (31), (32), (37) and (38) have been chosen because they looked natural to experiment with, Eq. (34) is a harder example found in Koutschan's list, Eq. (33) as well as Eq. (36) and it special case Eq. (35) come from the classical book of integral and series by Prudnikov et al. [28], and finally Eq. (39) is an example where Koutschan's heuristic does not stop as it does not guess correctly the form of the ansatz to use [11].

<footnote n="2">²The implementation is available at https://github.com/HBrochet/CreativeTelescoping.git, together with sessions of examples.</footnote>

<footnote n="3">³The code was run on a Intel Core i7-1265U with 32 GB of RAM.</footnote>
<table><caption>TABLE 2. Timings on the family $S_r$ from Eq. (40).</caption><thead><tr><th></th><th>HF-CT</th><th>HF-FCT</th><th>reduction</th></tr></thead><tbody><tr><td>S<sub>6</sub></td><td>11s</td><td>64s</td><td>0.4s</td></tr><tr><td>S<sub>7</sub></td><td>32s</td><td>331s</td><td>0.6s</td></tr><tr><td>S<sub>8</sub></td><td>106s</td><td>1044s</td><td>1.0s</td></tr><tr><td>S<sub>9</sub></td><td>325s</td><td>3341s</td><td>2.5s</td></tr><tr><td>S<sub>10</sub></td><td>1035s</td><td>&gt;1h</td><td>5.7s</td></tr></tbody></table>

$$ (30) \quad \sum_{j=1}^{n} \binom{m+x}{m-i+j} c_{n,j} \quad \text{where } c_{n,j} \text{ satisfies recurrences of order 2 [2, p. 6]} $$

$$ (31) \quad \sum_{n=0}^{\infty} C_n^{(k)}(x) C_n^{(k)}(y) \frac{u^n}{n!} $$

$$ (32) \quad \sum_{n=0}^{\infty} J_n(x) C_n^{(k)}(y) \frac{u^n}{n!} $$

$$ (33) \quad \sum_{k=0}^{\infty} (-1)^k (4k+1) J_{2k+1/2}(w) P_{2k}(z) $$

$$ (34) \quad \sum_{n=0}^{\infty} \frac{(4n+1)(2n)!}{n!^2 2^{2n} \sqrt{x}} J_{2n+1/2}(x) P_{2n}(u) $$

$$ (35) \quad \sum_{k=0}^{\infty} \frac{(b+3/2)_k}{(3/2)_k (b+1)_k} J_k^{(1/2,b)}(x) J_k^{(1/2,b)}(y) $$

$$ (36) \quad \sum_k \frac{(a+b+1)_k}{(a+1)_k (b+1)_k} J_k^{(a,b)}(x) J_k^{(a,b)}(y) $$

$$ (37) \quad \sum_{n=0}^{\infty} P_n(x) P_n(y) P_n(1/2) $$

$$ (38) \quad \sum_{n=0}^{\infty} P_n(x) P_n(y) P_n(z) $$

$$ (39) \quad \sum_{y} \frac{4x+2}{(45x+5y+10z+47)(45x+5y+10z+2)(63x-5y+2z+58)(63x-5y+2z-5)} $$

The family ($S_r$) is defined by [22]

$$ (40) \quad S_r = \sum_{k=0}^{n} \frac{(-1)^k (rn - (r-1)k)!(r!)^k}{(n-k)!rk!} $$

For any $r$, our algorithm produces a minimal telescoper of order $r$ and degree $r(r-1)/2$. The timings are reported in Table 2. It is unclear why the heuristic HF-FCT does not perform well on this family.

On most of these examples, the main part of the time of the computation is spent in the reductions in the call to CanonicalForm in Algorithm 1. For the two similar sums of Eqs. (33) and (34), almost half of the time is spent in Algorithm 4 performing the reductions needed to compute the bases for the strong reduction. This step is crucial to ensure that the minimal order elements in the telescoping ideal are found.

There are cases, like Eq. (36) and the family $S_r$, where the intermediate rational functions $R_\alpha$ in Eq. (17) become much larger than the telescopes found after linear algebra on them. In such situations, the direct, non-incremental approach taken by HF-CT and HF-FCT can be more efficient, by avoiding an unnecessarily large basis of rational functions.
REFERENCES

[1] Gert Almkvist and Doron Zeilberger. The method of differentiating under the integral sign. *Journal of Symbolic Computation*, 10:571–591, 1990.

[2] Tewodros Amdeberhan, Christoph Koutschan, and Doron Zeilberger. Yay for determinants! Technical Report 2023-18, RICAM Reports Series, 2023. Submitted for publication (preprint on arXiv:2307.01912).

[3] Louis C. Barrett and Forrest Dristy. Particular solutions for nonhomogeneous, linear, ordinary difference equations. *The American Mathematical Monthly*, 67(1):71–73, 1960.

[4] Alin Bostan, Shaoshi Chen, Frédéric Chyzak, and Ziming Li. Complexity of creative telescoping for bivariate rational functions. In *Proceedings of the 2010 International Symposium on Symbolic and Algebraic Computation (ISSAC 2010)*, pages 203–210. ACM Press, 2010.

[5] Alin Bostan, Shaoshi Chen, Frédéric Chyzak, Ziming Li, and Guoce Xin. Hermite reduction and creative telescoping for hyperexponential functions. In *Proceedings of the 38th International Symposium on Symbolic and Algebraic Computation*, ISSAC '13, pages 77–84, New York, NY, USA, 2013. ACM.

[6] Alin Bostan, Frédéric Chyzak, Pierre Lairez, and Bruno Salvy. Generalized Hermite reduction, creative telescoping and definite integration of D-finite functions. pages 95–102, 2018.

[7] Alin Bostan, Louis Dumont, and Bruno Salvy. Efficient algorithms for mixed creative telescoping. In *ISSAC'16—Proceedings of the 2016 ACM International Symposium on Symbolic and Algebraic Computation*, pages 127–134. ACM Press, 2016.

[8] Alin Bostan, Pierre Lairez, and Bruno Salvy. Creative telescoping for rational functions using the Griffiths-Dwork method. In Manuel Kauers, editor, *ISSAC '13: Proceedings of the 38th International Symposium on Symbolic and Algebraic Computation*, pages 93–100. ACM Press, 2013.

[9] Manuel Bronstein and Marko Petkovšek. An introduction to pseudo-linear algebra. *Theoretical Computer Science*, 157:3–33, 1996.

[10] Shaoshi Chen, Lixin Du, Manuel Kauers, and Rong-Hua Wang. Reduction-based creative telescoping for p-recursive sequences via integral bases. Technical Report 2311.05246, arXiv, 2023.

[11] Shaoshi Chen, Qing-Hu Hou, Hui Huang, George Labahn, and Rong-Hua Wang. Constructing minimal telescopes for rational functions in three discrete variables. *Advances in Applied Mathematics*, 141:Paper No. 102389, 31, 2022.

[12] Shaoshi Chen, Hui Huang, Manuel Kauers, and Ziming Li. A modified Abramov-Petkovšek reduction and creative telescoping for hypergeometric terms. In *ISSAC'15: Proceedings of the 2015 ACM on International Symposium on Symbolic and Algebraic Computation*, pages 117–124. ACM, 2015.

[13] Shaoshi Chen, Manuel Kauers, and Christoph Koutschan. Reduction-based creative telescoping for algebraic functions. In *Proceedings of the ACM on International Symposium on Symbolic and Algebraic Computation*, pages 175–182, New York, NY, USA, 2016. ACM.

[14] Shaoshi Chen, Mark van Hoeij, Manuel Kauers, and Christoph Koutschan. Reduction-based creative telescoping for fuchsian D-finite functions. *Journal of Symbolic Computation*, 85:108–127, 2018.

[15] R. C. Churchill and Jerald J. Kovacic. Cyclic vectors. In *Differential Algebra and Related Topics*, pages 191–218. World Scientific, 2002.

[16] Frédéric Chyzak. An extension of Zeilberger's fast algorithm to general holonomic functions. *Discrete Mathematics*, 217(1-3):115–134, 2000.

[17] Frédéric Chyzak. *The ABC of Creative Telescoping*. Mémoire d'habilitation à diriger des recherches, Université Paris-Sud, 2014.

[18] Frédéric Chyzak, Assia Mahboubi, Thomas Sibut-Pinote, and Enrico Tassi. A computer-algebra-based formal proof of the irrationality of $\zeta(3)$. In *Interactive Theorem Proving*, pages 160–176. Springer International Publishing, 2014.

[19] Frédéric Chyzak and Bruno Salvy. Non-commutative elimination in Ore algebras proves multivariate holonomic identities. *Journal of Symbolic Computation*, 26(2):187–227, August 1998.

[20] Lixin Du. On the existence of telescopes for p-recursive sequences. Technical Report 2311.06065, arXiv, 2023.

[21] J. Gerhard, M. Giesbrecht, A. Storjohann, and E. V. Zima. Shiftless decomposition and polynomial-time rational summation. In *Proceedings of the 2003 International Symposium on Symbolic and Algebraic Computation*, pages 119–126, New York, 2003. ACM.

[22] J. Gillis, B. Reznick, and D. Zeilberger. On elementary methods in positivity theory. *SIAM Journal on Mathematical Analysis*, 14(2):396–398, 1983.

[23] P. A. Hendricks and M. F. Singer. Solving difference equations in finite terms. *Journal of Symbolic Computation*, 27(3):239–259, 1999.

[24] N. Jacobson. Pseudo-linear transformations. *Annals of Mathematics. Second Series*, 38(2):484–507, 1937.

[25] Christoph Koutschan. *Advanced Applications of the Holonomic Systems Approach*. PhD thesis, RISC-Linz, 2009.

[26] Christoph Koutschan. A fast approach to creative telescoping. *Mathematics in Computer Science*, 4:259–266, 2010.

[27] Frank W. J. Olver, Daniel W. Lozier, Ronald F. Boisvert, and Charles W. Clark, editors. *NIST Handbook of Mathematical Functions*. Cambridge University Press, 2010.

[28] A. P. Prudnikov, Yu. A. Brychkov, and O. I. Marichev. *Integrals and Series. Volume 2: Special functions*. Gordon and Breach, 1986. 750 pages. First edition in Moscow, Nauka, 1983.

[29] Nobuki Takayama. An algorithm of constructing the integral of a module — an infinite dimensional analog of Gröbner basis. In *Symbolic and algebraic computation*, pages 206–211. ACM, 1990. Proceedings of ISSAC'90, Kyoto.
[30] Nobuki Takayama. Gröbner basis, integration and transcendental functions. In *Symbolic and algebraic computation*, pages 152–156. ACM, 1990. Proceedings ISSAC’90, Kyoto.

[31] Joris Van Der Hoeven. Constructing reductions for creative telescoping. Technical Report 014365877, Hal, 2017. Working paper or preprint.

[32] Joris van Der Hoeven. Creative telescoping using reductions. Technical Report 01773137, Hal, April 2018. Working paper or preprint.

[33] Joris van der Hoeven. Constructing reductions for creative telescoping: the general differentially finite case. *Applicable Algebra in Engineering, Communication and Computing*, 32(5):575–602, 2021.

[34] Joachim von zur Gathen and Jürgen Gerhard. *Modern computer algebra*. Cambridge University Press, New York, 3rd edition, 2013.

[35] Herbert S. Wilf and Doron Zeilberger. An algorithmic proof theory for hypergeometric (ordinary and "q") multisum/integral identities. *Inventiones Mathematicae*, 108:575–633, 1992.

[36] Doron Zeilberger. A fast algorithm for proving terminating hypergeometric identities. *Discrete Mathematics*, 80(2):207–211, 1990.

[37] Doron Zeilberger. A holonomic systems approach to special functions identities. *Journal of Computational and Applied Mathematics*, 32(3):321–368, 1990.

[38] Doron Zeilberger. The method of creative telescoping. *Journal of Symbolic Computation*, 11:195–204, 1991.