# Reduction-Based Creative Telescoping for Fuchsian D-finite Functions

Continuing a series of articles in the past few years on creative telescoping using reductions, we adapt Trager's Hermite reduction for algebraic functions to fuchsian D-finite functions and develop a reduction-based creative telescoping algorithm for this class of functions, thereby generalizing our recent reduction-based algorithm for algebraic functions, presented at ISSAC 2016.

## Authors

• 17 publications
• 8 publications
• 23 publications
• 15 publications
02/01/2016

### Reduction-Based Creative Telescoping for Algebraic Functions

Continuing a series of articles in the past few years on creative telesc...
02/12/2021

### Lazy Hermite Reduction and Creative Telescoping for Algebraic Functions

Bronstein's lazy Hermite reduction is a symbolic integration technique t...
01/26/2017

### Bounds for Substituting Algebraic Functions into D-finite Functions

It is well known that the composition of a D-finite function with an alg...
10/24/2017

### Definite Sums of Hypergeometric Terms and Limits of P-Recursive Sequences

The ubiquity of the class of D-finite functions and P-recursive sequence...
10/19/2017

### A Fast and Generic GPU-Based Parallel Reduction Implementation

Reduction operations are extensively employed in many computational prob...
02/01/2021

### A Brief Account of Klein's Icosahedral Extensions

We present an alternative relatively easy way to understand and determin...
02/07/2020

### Integral P-Recursive Sequences

In an earlier paper, the notion of integrality known from algebraic numb...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The classical question in symbolic integration is whether the integral of a given function can be written in “closed form”. In its most restricted form, the question is whether for a given function  belonging to some domain there exists another function , also belonging to , such that . For example, if is the field of rational functions, then for we can find , while for no suitable exists. When no exists in , there are several other questions we may ask. One possibility is to ask whether there is some extension  of such that in there exists some with . For example, in the case of elementary functions, Liouville’s principle restricts the possible extensions , and there are algorithms which construct such extensions whenever possible. Another possibility is to ask whether for some modification of  there exists a such that . Creative telescoping leads to a question of this type. Here we are dealing with domains  containing functions in several variables, say and , and the question is whether there is a linear differential operator , nonzero and free of , such that there exists a with , where denotes the derivative of with respect to . Typically,  itself has the form for some operator (which may be zero and need not be free of ). In this case, we call a telescoper for , and a certificate for .

Creative telescoping is the backbone of definite integration, because implies, for instance, . A telescoper for  thus gives rise to an annihilating operator for the definite integral .

###### Example 1 ((Manin, 1958)).

The algebraic function

 f(x,t)=1√x(x−1)(x−t)

does not admit an elementary integral with respect to . However, we have for

 P=4(t−1)tD2t+4(2t−1)Dt+1,Q=2x(x−1)t−x.

This implies

 P⋅∫10f(x,t)dx=[2x(x−1)t−xf(x,t)]x=1x=0,

so the integral satisfies the differential equation

 4(t−1)tF′′(t)+4(2t−1)F′(t)+F(t)=0.

In the common case when the right-hand side collapses to zero, we say that the integral has “natural boundaries”. Readers not familiar with creative telescoping are referred to the literature (Petkovšek et al., 1996; Zeilberger, 1990a, 1991, b; Koepf, 1998; Kauers and Paule, 2011) for additional motivation, theory, algorithms, implementations, and applications. There are several ways to find telescopers for a given . In recent years, an approach has become popular which has the feature that it can find a telescoper without also constructing the corresponding certificate. This is interesting because certificates tend to be much larger than telescopers, and in some applications, for instance when an integral has natural boundaries, only the telescoper is of interest. This approach was first formulated for rational functions by Bostan et al. (2010) and later generalized to rational functions in several variables (Bostan et al., 2013b; Lairez, 2016), to hyperexponential functions (Bostan et al., 2013a) and, for the shift case, to hypergeometric terms (Chen et al., 2015; Huang, 2016) and binomial sums (Bostan et al., 2016). At ISSAC’16, three of the present authors have given a version for algebraic functions (Chen et al., 2016). In the present article, we extend this algorithm to fuchsian D-finite functions.

The basic principle of the general approach is as follows. Assume that the -constants form a field and that

is a vector space over the field of

-constants. Assume further that there is some -linear map such that for every there exists a with . Such a map is called a reduction. For example, in Hermite reduction (Hermite, 1872) decomposes any into with such that is a proper rational function (i.e., the numerator degree is smaller than the denominator degree) with a squarefree denominator. In this case, we can take . In order to find a telescoper, we can compute , , , …, until we find that they are linearly dependent over . Once we find a relation , then, by linearity, , and then, by definition of , there exists a such that . In other words, is a telescoper.

There are two ways to guarantee that this method terminates. The first requires that we already know for other reasons that a telescoper exists. The idea is then to show that the reduction has the property that when is such that there exists a with , then . If this is the case and is a telescoper for , then is integrable in , so , and by linearity , …, are linearly dependent over . This means that the method won’t miss any telescoper. In particular, this argument has the nice feature that we are guaranteed to find a telescoper of smallest possible order . This approach was taken by Chen et al. (2015). The second way consists in showing that the -vector space generated by has finite dimension. This approach was taken by Bostan et al. (2010, 2013a). It has the nice additional feature that every bound for the dimension of this vector space gives rise to a bound for the order of the telescoper. In particular, it implies the existence of a telescoper.

In this paper, we generalize Trager’s Hermite reduction for algebraic functions to fuchsian D-finite functions (Section 4), which yields a reduction-based creative telescoping algorithm via the first approach (Section 5). In Section 6, we introduce a second reduction, called polynomial reduction, which together with the Hermite reduction gives rise to a reduction-based creative telescoping algorithm via the second approach. It also delivers a new bound for the order of the telescoper of a fuchsian D-finite function, and in particular an independent proof for its existence.

## 2 Fuchsian D-finite Functions

Throughout the paper, let be a field of characteristic zero. We consider linear differential operators with belonging to some ring containing . Typical choices for will be or . When , we say that is the order of .

Let be a differential ring and write for its derivation. We write for the algebra consisting of all linear differential operators, together with the usual addition and the unique non-commutative multiplication satisfying for all . We shall assume throughout that . The algebra acts on a differential -module  via

 (ℓ0+ℓ1∂x+⋯+ℓn∂nx)⋅f=ℓ0f+ℓ1f′+⋯+ℓnf(n).

An element is called a solution of an operator if .

By we denote some algebraically closed field containing  (not necessarily the smallest). An operator of order  is called fuchsian at a point if it admits linearly independent solutions in

 ¯C[[[x−a]]]:=⋃ν∈C(x−a)ν¯C[[x−a]][log(x−a)].

It is called fuchsian at if it admits linearly independent solutions in

 ¯C[[[x−1]]]:=⋃ν∈Cx−ν¯C[[x−1]][log(x)].

It is simply called fuchsian if it is fuchsian at all . Note that the exponents are restricted to , not to the larger field . For simplicity, the dependence on is not reflected in the notation.

Examples for fuchsian operators are operators that have a basis of algebraic function solutions, the Gauss hypergeometric differential operator, or the operator , whose solutions are and . However, the class of fuchsian D-finite functions considered in this paper is not as rich as it may seem at first glance, because we require that the operators under consideration should be fuchsian at all points including infinity. Functions such as , , , Bessel functions, etc. are only fuchsian at all finite points but not at infinity.

For a fixed fuchsian operator , we will consider the left -module , where denotes the left ideal generated by  in . Then is a solution of , because we have in . We can say that consists of all the “functions” which can be obtained from a “generic” solution  of by applying some operator to it. When is a field, then is an -vector space of dimension , generated by .

It is instructive to compare this setup to the situation for algebraic functions. Comparing to an algebraic function field (when is a field), our operator plays the role of the minimal polynomial . In the algebraic case,  is a formal solution of the equation , similar as is a formal solution of . Besides these formal solutions there are, for each fixed , exactly different Puiseux series solutions of at places above . They correspond in the differential setting to the series solutions of in , which generate a -vector space of dimension .

The exponents of an element at a point are the values such that one of the series in (or , respectively) associated to has (or ) as initial term. For an element , let be the minimal order of an operator with . We say that is an ordinary point of if the set of exponents of at is and the solutions at do not involve logarithms. There can be at most finitely many non-ordinary points; these are called the singular points. The defect of at , denoted , is defined as the sum of the exponents of at minus . Then the Fuchs relation (Schlesinger, 1895; Ince, 1926) says that we have

 ∑a∈¯C∪{∞}defecta(f)=nf(1−nf)

for all . This relation is the counterpart of the Bézout relation in the algebraic case. Note that when is an ordinary point, then , but does in general not imply that is an ordinary point.

In the context of creative telescoping, we let be some algebraically closed field containing the rational function field , and we use instead of . Integration will always be with respect to , but besides the derivation there is now also the derivation with respect to . The notation will always refer to the derivative of with respect to , not with respect to . In addition to the operator algebra , we consider the operator algebra , in which commute with each other (although they need not commute with elements of ).

The action of on or is extended to by letting act coefficient-wise. We further assume that the action of on is extended to an action of on  in a way that is compatible with the action of on series domains. This means that when is a solution of  and is an element of , so that is an element of , then we want to have , where the in refers to the action of on and the three other dots refer to the action on .

If is such that and is such that , then the annihilator of in contains and . We therefore have , which is the usual definition of D-finiteness in the case of several variables (Zeilberger, 1990b; Chyzak and Salvy, 1998; Koutschan, 2009; Kauers, 2015). Pathological situations, where we also have but does not have a basis of the form , are not considered in this paper.

## 3 Integral Bases

Trager’s Hermite reduction for algebraic functions rests on the notion of integral bases. The notion of integral bases has been generalized to D-finite functions (Kauers and Koutschan, 2015), and an algorithm was also given there for computing such bases. We recall here the relevant definitions and properties.

Although the elements of a generalized series ring are formal objects, the series notation suggests certain analogies with complex functions. For simplicity, let us assume throughout that . Terms or are called integral if . A series in or is called integral if it only contains integral terms. A non-integral series is said to have a pole at the reference point. Note that in this terminology also has a pole at , while does not; this convention differs slightly from the default setting of (Kauers and Koutschan, 2015, Ex. 2), but can be achieved by defining the function  (Kauers and Koutschan, 2015, Def. 1) accordingly. Note also that the terminology only refers to but not to .

Integrality at is not preserved by differentiation, but if is integral at , then so is . On the other hand, integrality at infinity is preserved by differentiation; we even have the stronger property that when is integral at infinity, then not only but also is integral at infinity.

Let be some field with . Let be a fuchsian operator. An element is called (locally) integral at if for every solution of in or , respectively, the series is integral.  is called (globally) integral if it is locally integral at every (“at all finite places”).

For an element to have a “pole” at means that is not locally integral at ; to have a “double pole” at means that (or if ) is not integral; to have a “double root” at means that (or if ) is integral, and so on.

The set of all globally integral elements forms a -submodule of . A basis of this module is called an integral basis for . Kauers and Koutschan (2015) proposed an algorithm which computes an integral basis for a given . This algorithm is a generalization of van Hoeij’s algorithm (van Hoeij, 1994) for computing integral bases of algebraic function fields (Trager, 1984; Rybowicz, 1991).

For a fixed , let be the ring of rational functions with , and write for the ring of all rational functions with . Then the set of all which are locally integral at some fixed forms a -module. A basis of this module is called a local integral basis at for . The algorithm given by Kauers and Koutschan (2015) for computing (global) integral bases computes local integral bases at finite points as an intermediate step. By an analogous algorithm, it is also possible to compute a local integral basis at infinity.

An integral basis is always also a -vector space basis of . A key feature of integral bases is that they make poles explicit. Writing an element as a linear combination for some , we have that has a pole at if and only if at least one of the has a pole there.

###### Lemma 2.

Let be a fuchsian operator and let be a local integral basis of at . Let and be such that . Then is integral at if and only if each is integral at .

###### Proof.

The direction “” is obvious. To show “”, suppose that is integral at . Then there exist such that . Thus , and then for all , because is a vector space basis of . As elements of , the are integral at , and hence also all the are integral at .

The lemma says in particular that poles of the in a linear combination cannot cancel each other.

###### Lemma 3.

Let be a fuchsian operator and let be an integral basis of . Let and be such that

 eω′i=n∑j=1mi,jωj

for and . Then is squarefree.

###### Proof.

Let be a root of . We show that is not a multiple root. Since is integral, it is in particular locally integral at . Therefore is locally integral at . Since is an integral basis, it follows that for all . Because of , no factor of can be canceled by all the . Therefore the factor can appear in only once.

###### Lemma 4.

Let be a fuchsian operator and let be a local integral basis at infinity of . Let and be defined as in Lemma 3. Then for all .

###### Proof.

Since every is locally integral at infinity, so is every . Since is an integral basis at infinity, it follows that for all . This means that for all , and therefore , as claimed.

A -vector space basis of is called normal at if there exist such that is a local integral basis at . Trager (1984) shows for the case of algebraic function fields how to construct an integral basis which is normal at infinity from a given integral basis and a given local integral basis at infinity. The same procedure also applies in the present situation. It works as follows.

Let be a global integral basis and be a local integral basis at infinity. Let be such that

 ωi=n∑j=1mi,jνj.

For each , let be the largest integer such that has no pole at infinity for any . Then each is locally integral at infinity. Let be the matrix obtained by evaluating at infinity (this is possible by the choice of ). If is invertible, then the form a local integral basis at infinity and we are done. Otherwise, there exists a nonzero vector with . Among the indices with choose one where is minimal, and then replace by . Note that the resulting basis is still global integral. Repeating the process, it can be checked that the value of strictly increases in each iteration. According to the following lemma, the sum is bounded, so the procedure must terminate after a finite number of iterations.

###### Lemma 5.

Let , , and be as above. Let be the number of points where at least one of the does not have distinct exponents in  (i.e., counts the finite singular points of  that are not “apparent” singularities). Then

 τ1+⋯+τn≤12n(n−1)(N−1).
###### Proof.

We show that when is such that is locally integral at infinity, then , for every . Let be the minimal order of an operator with . By the Fuchs relation we have

 ∑a∈¯C∪{∞}defecta(ωi)=ni(1−ni),

hence

 defect∞(ωi)=ni(1−ni)−∑a∈¯Cdefecta(ωi).

When all exponents of the series associated to at form a subset of  of size , then . At all other points , of which there are at most

by assumption, we still have the estimate

, because is integral at all finite places. It follows that

 defect∞(ωi)≤ni(1−ni)+12ni(ni−1)N=12ni(ni−1)(N−2).

Next, for every we have . Moreover, if is such that is integral at infinity, then we must have , i.e.,

 defect∞(ωi)−τni≥−12ni(ni−1),

and hence,

 τ ≤1ni(12ni(ni−1)+defect∞(ωi)) ≤12(ni−1)+12(ni−1)(N−2)=12(ni−1)(N−1)≤12(n−1)(N−1)

as claimed.

Although normality is a somewhat weaker condition on a basis than integrality, it also excludes the possibility that poles in the terms of a linear combination of basis elements can cancel:

###### Lemma 6.

Let be a fuchsian operator and let be a basis of which is normal at some . Let for some . Then has a pole at if and only if there is some such that has a pole at .

###### Proof.

Let be such that is a local integral basis at . By and by Lemma 2,  is integral at  iff all are integral at .

We will mostly be using bases that are integral at every point in except one. For the case of algebraic functions, the reason is that the only algebraic functions which are integral at all finite places and also at infinity are the constant functions (Chevalley’s theorem (Chevalley, 1951, p. 9, Cor. 3)). The results of Chen et al. (2016) depend heavily on this fact. There is no analogous result for fuchsian D-finite functions: such functions may be integral at all finite places and also at infinity without being constant. It is easy to construct examples using the Papperitz symbol.

###### Example 7.

The operator has three singular points . Infinity is an ordinary point of . At all three singularities, there is one local solution starting with exponent  and another starting with exponent , so all the solutions are integral everywhere according to our standard definition of integrality of generalized series.

Fortunately, we can still be sure that there are not too many such functions.

###### Lemma 8.

Let for some fuchsian operator , let be a global integral basis of  which is normal at infinity, and let be such that is a local integral basis at infinity. Denote by the set of all which are integral at all finite places and at infinity. Then is a -vector space of finite dimension, and is a basis of .

###### Proof.

It is clear that is closed under taking -linear combinations, so it is clearly a vector space. We show that is a basis.

Every is by definition integral at all finite places and because of also integral at infinity. Therefore , and therefore generates a subspace of .

Conversely, let be arbitrary. Then is in particular integral at all finite places, and since is a global integral basis we can write for some polynomials . If we had for some , then would not be integral at infinity (by definition of the numbers ), and then, because is normal at infinity, Lemma 6 implies that would not be integral at infinity. But is in and therefore integral at all points, including infinity. It follows that for all , and therefore is a -linear combination of elements of .

###### Corollary 9.

With the notation of Lemmas 5 and 8, we have

 dimK(V)≤n(12(n−1)(N−1)+1).
###### Proof.

Note that the exponents  in Lemma 8 are the same as in Lemma 5; in the proof of the latter it was shown that . The estimate on the dimension of  then follows immediately.

## 4 Hermite Reduction

Hermite reduction was first introduced by Hermite (Hermite, 1872) for rational functions. This reduction was later extended to elementary functions by Risch (Risch, 1969, 1970; Geddes et al., 1992; Bronstein, 1998b, 2005), to algebraic functions by Trager (Trager, 1984; Geddes et al., 1992; Bronstein, 1998b), and to hyperexponential functions (Bostan et al., 2013a). These generalizations are the key step in many integration algorithms, including the telescoping algorithm presented in this paper. It turns out that the Hermite reduction for fuchsian D-finite functions is literally the same as Trager’s reduction for algebraic functions.

We start with a technical lemma, which is needed later to ensure that the Hermite reduction always works. The analogous statement for algebraic functions and its proof can be found in (Trager, 1984, pp. 46–47); Trager’s proof for the algebraic case directly carries over and is reproduced here only for the convenience of the reader.

Throughout this section, let be a fuchsian operator of order  and let .

###### Lemma 10.

Let be a squarefree polynomial and let be a basis of  that is locally integral at all roots of . For some integer we define ; then is a local integral basis at each root of .

###### Proof.

By expanding one sees that the themselves are integral at all roots of . Let now be an arbitrary but fixed root of . We have to show that each that is integral at  can be expressed as a linear combination of the with coefficients in . To the contrary, assume that there exists an integral element  that requires in the denominator of some coefficient, i.e.,

 f=1vn∑i=1ciψiwith ci∈¯C(x)a and ci(a)≠0 for some i

(here we use the fact that is squarefree). Further let , which is obviously integral. Then also their sum

must be integral. Since is an integral basis at , there exists for each a series solution of  such that involves a term with and . Let now be an index such that ; this implies that appears in . Using the fact that the form a local integral basis, it follows by Lemma 2 that is also present in where . Let now be the dominant term of , i.e., among all terms with minimal the one with the largest exponent . From

 ((x−a)μ−1∂x(x−a)1−μ)⋅T=(1−μ+α)(x−a)α−1log(x−a)β+β(x−a)α−1log(x−a)β−1 (1)

it follows that is the dominant term of ; here we use the assumption that , because for and the coefficient in (1) is zero. This calculation reveals that is not integral at , which contradicts our assumption on the integrality of . Hence is a local integral basis at .

Let be an integral basis for . Further let () be such that and as in Lemma 3. For describing the Hermite reduction we fix an integrand and represent it in the integral basis, i.e., with . The purpose is to find such that and with and denoting the squarefree part of . As differentiating the can introduce denominators, namely the factors of , it is convenient to consider those denominators from the very beginning on, which means that we shall assume . Note that can then be nontrivial.

We now execute one step of the Hermite reduction, where the multiplicity of some nontrivial squarefree factor  of is reduced. Let be such that ; it follows that and . We want to find such that

 n∑i=1fiuvμωi=(n∑i=1givμ−1ωi)′+n∑i=1hiuvμ−1ωi. (2)

By a repeated application of such reduction steps one can decompose any as where the denominators of the coefficients of are squarefree and the coefficients of are proper rational functions.

In order to determine the unknown polynomials in (2), clearing the denominator yields

 n∑i=1fiωi=n∑i=1(uvg′iωi+uvμgi(v1−μωi)′+vhiωi), (3)

and then this equation is reduced modulo :

 n∑i=1fiωi=n∑i=1giuvμ(v1−μωi)′modv. (4)

By Lemma 10 and from it follows that the elements form a local integral basis at each root of , which implies that the coefficients are uniquely determined modulo .

By Lemma 3 the polynomial  is squarefree and therefore ; hence we can write for some . By rewriting the derivatives of the in terms of the integral basis, Equation (4) turns into

 n∑i=1fiωi =n∑i=1gi(uvω′i−(μ−1)uv′ωi)modv =n∑i=1gi(wn∑j=1mi,jωj−(μ−1)uv′ωi)modv.

Comparing coefficients w.r.t.  yields a system of linear equations over for the unknown functions