Creative telescoping using reductions

Joris van der Hoeven

Laboratoire d'informatique, UMR 7161 CNRS

Campus de l'École polytechnique

1, rue Honoré d'Estienne d'Orves

Bâtiment Alan Turing, CS35003

91120 Palaiseau

Email: vdhoeven@lix.polytechnique.fr

Draft version, November 4, 2023

Creative telescoping is a popular method for proving combinatorial identities and the computation of parametric integrals that involve special functions. Traditional implementations of this method admit an exponential bit complexity and it is an open problem under which conditions creative telescoping can be achieved in polynomial time. More efficient reduction-based algorithms were recently introduced in order to get a better grip on such complexity issues. Initially, reduction-based algorithms only applied to special cases such as rational, algebraic, or hyperexponential functions. More recently, constructions of reductions appeared for larger classes of Fuchsian D-finite and general differentially-finite functions.

In this paper, we show how to construct reductions for mixed differential-difference systems, where the difference operators are either shift operators or -difference operators. We recall how this yields an algorithm for creative telescoping and specify under which precise conditions on singularities this algorithm works. For creative telescoping of differential type, we next examine the complexity of our algorithms and prove a polynomial complexity bound. The algorithm for which this bound holds computes generators for a D-finite ideal of telescopers, but not necessarily a Gröbner basis for the ideal of all telescopers.

Keywords: creative telescoping, holonomic function, Hermite reduction, residues

A.C.M. subject classification:

I.1.2 Algebraic algorithms

A.M.S. subject classification: 33F10, 68W30

1.Introduction

1.1.Creative telescoping

The technique of creative telescoping is a powerful tool for proving and finding combinatorial identities and the computation of parametric integrals that involve special functions. The name was coined by van der Poorten [77] in relation with Apéry's irrationality proof of , but the first systematic algorithms were only developed about one decade later by Zeilberger and followers [90, 4, 91, 62, 85, 86]. Early precursors are [35, 89, 42].

Let us briefly recall the main idea behind creative telescoping on an extremely simple example: assume that we wish to prove the well known identity

Setting and , the defining relation of in Pascal's triangle yields

Summing this relation from until , the right hand side becomes a telescoping sum, and we obtain . Since , the recurrence relation proves .

More generally, for the evaluation of a sum of the form , the idea is to construct a relation of the form

(1.1)

where the left-hand side does not depend on and the right-hand telescopes when summing over . If we manage to construct such a relation, then the creative telescoping process provides us with a non-trivial linear recurrence relation for . Writing for the operator that sends to , the difference operator that annihilates in (1.1) is called a telescoper, and the function its corresponding certificate.

A similar idea works for integrals of the form instead of discrete sums. In that case, we need a relation of the form

(1.2)

instead of (1.1). Again, is called the telescoper and the certificate. Several other variants can be considered as well, such as -difference equations, multiple sums and integrals, etc.

After Zeilberger's initial impetus, the development of creative telescoping focussed on making implementations more efficient, while guaranteeing termination and widening the class of functions on which the theory can be applied. Whereas early implementations were particularly successful for first or low order recurrence relations (e.g. hypergeometric summation), efficient algorithms soon appeared for more general holonomic functions and D-finite ideals in Ore algebras [26, 29, 27, 63]. For a historical perspective on these advances and further references, we refer to [28].

One major application of creative telescoping is proving functional identities. In the purely differential setting, it should be noticed that alternative methods have been developed for larger classes of non-linear equations [32, 33, 83, 75, 46, 49].

1.2.From Hermite reduction to reduction-based telescoping

Despite the above successes, the efficiency of creative telescoping remained a major issue. Part of the problem was due to the fact that the construction of relations of the form (1.1) or (1.2) often relied on brute force computations rather than a better insight into their precise shape.

The next major step was the introduction of reduction-based algorithms as a remedy to these problems. The concept of a reduction goes back to the works of Ostrogradsky and Hermite on the integration of rational functions [74, 44]. Hermite reduction is nowadays a fundamental tool in symbolic integration and generalizations for algebraic, hyperexponential and hypergeometric functions were introduced with this kind of application in mind [88, 30, 2]. It is interesting to notice that the importance of Hermite reduction for integrals that depend on parameters was understood quite early [76]; see also [39].

Hermite reduction was first related to creative telescoping in [41]; the complexity perspective appeared in [13]. Initially, only quite restricted types of functions could be treated this way and it became a major challenge to develop reduction-based algorithms for similarly general classes of functions as for the other approaches. This triggered a lot of activity in recent years [24, 14, 67, 22, 34, 23, 54, 68], culminating with an algorithm for Fuchsian differential equations [21]. This last result made it plausible that a fully general algorithm in the differential case might be quite complex, since it might involve elaborate arguments in order to deal with irregular singularities. Fortunately this concern turned out to be unwarranted: last year we constructed suitable reductions for arbitrary differential equations in a purely algebraic way [50]. Our technique has recently been further improved in [16]. In the present paper, we continue its development and extend it to mixed differential-difference equations.

Before we proceed, it is useful to recall how reductions show up in relation with creative telescoping. Assume that we wish to find a linear differential equation over for the function , where . Hermite reduction in provides us with a unique decomposition of any as a sum

(1.3)

where admits only simple poles in and . When restricting ourselves to functions whose poles in belong to a fixed finite set (i.e. for some fixed square-free polynomial ), it follows that is “confined” to a finite dimensional -vector space . Indeed, for some algebraic extension of that contains all roots of , the reduction is always an -linear combination of the fractions , where runs over the set of roots of in . In our case, we simply take to be the square-free part of the denominator of .

Now the crucial observation is that , , , etc. all belong to the above finite dimensional vector space, whence there exists a non-trivial relation

(1.4)

with . By construction, there exists an with for . It is not hard to check that the reduction commutes with partial derivation with respect to , whence (1.4) implies

(1.5)

This new relation is of the desired form (1.2). Another advantage of this method with respect to previous ones is that the computation of certificates is possible, but only optional. Since certificates tend to be more voluminous than the telescopers themselves, this is very interesting from a complexity point of view.

In order to generalize the above reduction-based approach to functions that satisfy higher order differential equations, one essentially needs to define a confined reduction on a so called narrow -submodule of that contains . In the differential case, such narrow modules are still of the same form . In the difference case, the set of allowed poles are translates of a finite number of points and their orders need to be bounded in a way that will be made precise.

1.3.On the choice of our setup

The theory of creative telescoping has become very wide, with many variants that cover specific cases more or less efficiently. For the exposition in this paper, we therefore had to make several choices concerning the setup. Let us briefly comment on these choices.

Non-commutative operator algebras. It was recognized by Zeilberger that the framework of holonomic functions is particularly suitable for the development of creative telescoping [90]. The reason is that there exists a systematic elimination theory for the equations satisfied by such functions that can be regarded as a non-commutative counterpart of the theory of Gröbner bases.

Elimination methods for differential equations actually predate Buchberger's famous algorithm [19]. One line of development is due to Riquier–Janet–Thomas [79, 56, 87] and another one to Ritt–Kolchin [80, 61]. The most “effective” variant of the latter theory [81, 17] essentially covers both commutative Gröbner basis theory, its skew differential counterpart for Weyl algebras, and non-linear generalizations that specialize to triangular system type solving in the algebraic case. Yet another line of development is due to Ore, who introduced a systematic theory for non-commutative polynomials [72, 73]. A more direct inspiration for Zeilberger's work were developments in the theory of D-modules in the seventies [9, 10, 12] and the development of non-commutative Gröbner basis theory in the computer algebra community in the eighties and later [40, 84, 60, 71, 64, 26, 70].

In this paper, we have chosen to work over differential-difference algebras generated by a finite number of coordinates , derivations , shift operators , and -difference operators . All derivations and automorphisms , are assumed to commute pairwise. Our setting is not the most general possible, but it has a natural geometric interpretation, as well as a non-linear counterpart [45]. The commutation of the derivations is not really essential, except for the main derivation that will be eliminated. It should also be quite straightforward to replace the algebra in section 6.4 by more general algebras of solvable type [60, 64], Ore algebras [26], or G-algebras [70, 43]. On the other hand, when telescoping with respect to a difference operator it is essential that singularities are moved in a bijective way. It is therefore natural to restrict one's attention to homographies. Then, modulo a change of variables, homographies always operate like shift operators or -difference operators: see section 5.1. Other interesting generalizations would concern nested extensions by solutions of linear differential-difference equations as in [78, 82].

Scalar equations versus first order systems. An eternal dilemma when dealing with linear differential or difference equations is whether one should privilege scalar higher order equations or first order systems. Thanks to the cyclic vector lemma [3], both approaches are essentially equivalent, so the choice is mainly a matter of taste. The best solution is probably to work with the representation that was used for the input, whenever possible. Our personal taste is slightly in favour of scalar equations, although we opted for first order systems in [50]. In [16], the authors preferred to work with scalar equations instead. In the present paper, we mainly consider scalar equations, but also outline how to adapt the theory to first order systems. At the end of the day, it seems that both approaches are more or less equally diligent.

Algebraic extensions versus rational counterparts. From the early days of symbolic computation, algorithms that require computations in algebraic extensions have acquired the reputation of being slow. Starting with Trager's work on the integration of rational functions [88], it has therefore become common practice to develop “rational” counterparts for reduction algorithms that avoid explicit computations in algebraic extensions. Nevertheless, from a complexity point of view, it is not so clear that working in algebraic extensions is necessarily that bad: computing in an extension of degree is roughly times more expensive, but typically deals with conjugate roots in a single pass. Although the development of rational counterparts remains an interesting topic, we therefore do not regard such algorithms as intrinsically better. In this paper, we have chosen to systematically work in algebraic extensions whenever appropriate, which is somewhat simpler from a conceptual standpoint.

Characteristic zero. It will be convenient to assume that all fields considered in this paper are of characteristic zero, although most results can be generalized in a straightforward way to fields of (sufficiently large) positive characteristic.

1.4.Structure of the paper and outline of the main results

The main purpose of this paper is to present reduction-based algorithms for creative telescoping in the above setting of differential-shift--difference equations and to derive complexity bounds for several of these algorithms. In the purely differential case, our algorithms are completely systematic. In the presence of difference operators it is well known that indefinite summation and integration under the integral sign do not always preserve D-finiteness. Nevertheless, we provide a sufficient condition (explicit telescopability) for our algorithms to work, as well as a reduction-based algorithm that is able to check D-finiteness in a particular case. Let us briefly outline the structure of the paper and summarize our contributions.

Abstract reductions. For a gentle introduction to reduction-based creative telescoping, we refer to [34, Section 1.2]. In section 2, we start with a more abstract presentation of the basic theory that covers all cases that occurred in the literature so far. We also introduce the notion of “local reductions” and describe various related constructions on an abstract level.

The Lagrange identity. Let be a field. One major new idea in [16] is the consideration of reductions with respect to general differential operators instead of just the derivation as in (1.3). More precisely, a reduction with respect to is a -linear projection such that for all . Now in the case of a rational function , ordinary Hermite reduction provided us with a way to derive telescoping relations (1.5) from linear dependencies (1.4). In the case of a function that satisfies a linear differential equation , it is shown how to do something similar using reductions with respect to the adjoint operator of . The technical tool that makes this possible is the Lagrange identity [55, section 5.3]. We recall this formula in section 3 and generalize it to difference and matrix operators.

Differential reductions. In section 4, we show how to construct reductions for differential equations. This case was already covered in [50], but it is instructive to present it with the formalisms from sections 2 and 3 at hand. In sections 4.2 and 4.3, we start with a variant of the construction from [16], whereas section 4.5 explains the link with [50] and also clarifies how the Lagrange identity can be used to simplify the head/tail chopping process. We also show how new-style reductions with respect to the adjoint operator give rise to old-style reductions in the sense of [50] and vice versa. In section 4.6, we finally consider mixtures of both reduction styles. This makes it possible to extend creative telescoping to the resolution of more general linear differential or difference equations that involve parameters: see Remark 6.7.

Difference reductions. The analogue construction of reductions for difference operators is the subject of section 5. Whereas derivations do not introduce new singularities, but only aggravate the order of poles, difference operators move singularities, but do not increase their severity. The idea behind the construction of differential reductions is to diminish the order of the poles of the function to be reduced. Likewise, difference reductions are based on the idea to shift singularities back until they are confined in a finite set. Another difference with section 4 is that the construction is somewhat less canonical due to the absence of a privileged section of orbits of singularities (see Remark 5.7), but this is of no consequence for the application to creative telescoping. Apart from that, the theory from section 4 naturally adapts to the difference setting.

Creative telescoping. After the formal introduction of the setting of DD-operator algebras in section 6.1, we next consider the application to creative telescoping. For suitable “telescopable” D-finite ideals of a DD-operator algebra , the theory from sections 4 and 5 allows us to define a computable confined reduction on a narrow submodule of . In section 6.4 we show how this can be used to construct telescopers. The algorithm is a straightforward generalization of the technique from section 1.2 and only relies on linear algebra. As an end-result, it produces a finite set of generators for a D-finite ideal of telescopers; in the favourable case when the reduction is “normal”, these generators actually form a Gröbner basis of the ideal of all telescopers. In the context of Gröbner basis, the algorithm is reminiscent of the FGML algorithm [37] and a variant of the one used in [16].

D-finiteness tests for telescoping ideals. In the purely differential case, it is a well-known fact from the theory of D-modules (reproved in section 7) that all D-finite ideals are telescopable, so the algorithm from section 6.4 systematically works. In the presence of difference operators, the telescopability of a D-finite ideals depends on the nature of its singularities. In section 7, we examine such singularities in more detail and derive sufficient conditions for the telescopability of a D-finite ideal with a given basis for the quotient module (a particular was treated before in [16, section 4.4]). Sometimes these “explicit telescopability” conditions are only met after a change of basis. For differential reductions, we present an algorithm to compute such a base change if it exists. Even when the explicit telescopability conditions do not hold, it sometimes still happens that the telescoping ideal is D-finite; in the last subsection 7.4, we present an algorithm that decides whether this is the case and that computes a Gröbner basis of the telescoping ideal if it is D-finite.

Complexity of rational function arithmetic. The last two sections are devoted to analyzing the complexity of reduction-based algorithms for creative telescoping. For this, we first recall basic complexity results for computations with rational functions in several variables, including partial fraction decomposition. When using the common dense representation, these complexity bounds are quite pessimistic, yet always polynomial if the number of variables is fixed. When allowing for probabilistic algorithms of Las Vegas type, it is more efficient to represent rational functions by the straight line programs (SLPs) [20], and fast dense univariate polynomial arithmetic can be used at every evaluation point.

Complexity of creative telescoping. In the last section we finally turn to the cost of creative telescoping itself. For this, we use a separate analysis for the complexity of the reduction process and the complexity of the FGLM-type algorithm for computing linear relations. The latter algorithm relies on standard linear algebra and its complexity analysis is straightforward in terms of the dimension of the vector space of reduced functions. We next turn to the reduction process and restrict our attention to differential reductions. In this case, we prove a polynomial complexity bound, as well as a polynomial bound for as a function of the input size. Altogether, this yields a polynomial time algorithm for creative telescoping, as promised in [50]. The complexity exponent is again extremely bad when using dense representations, but essentially cubic per evaluation point in the SLP model.

Remark 1.1. It is well known that it is easy to construct operators in of small bit size that admit polynomial solutions with a large number of terms. For instance, the operator admits as a solution. For the design of polynomial time algorithms for operators in it is therefore crucial to avoid the explicit computation of such solutions. For a similar reason, telescopers of minimal order (e.g. of the form in (1.5) with minimal ) can admit very large coefficients. The telescopers computed by our polynomial time algorithm in section 9 are therefore not necessarily minimal. The above phenomenon also implies that the algorithms at the end of section 7 (for base changes to reach an explicitly telescopable form and creative telescoping for not explicitly telescopable ideals) do not run in polynomial time, in general.

1.5.Pending issues

Let us briefly list a few problems that have not been tackled in the present draft by lack of time. We intend to address them in a future version or in a follow-up paper.

Base changes for explicit telescopability. The base change algorithm to remove apparent singularities from the end of section 7.3 has been presented in the differential case only. We are still looking for existing work on this topic in the literature and a similar algorithm for the difference case.

Testing D-finiteness. The general algorithm for creative telescoping from section 7.4 (in absence of explicit telescopability) has also been presented in the differential case only. The difference analogue seems to require a bit more work, but it should be feasible to work out an algorithm along similar lines.

Fast arithmetic for rational functions. Several complexity bounds from section 8 are not very sharp or a bit sketchy. We intend to improve this section in the next version.

Complexity bounds for the difference case. The final complexity bounds in sections 9.2 and 9.3 have only been presented in the differential case (which still allows for difference operators among the ). If and in Proposition 7.7 do not grow too fast, then it should be possible to prove polynomial complexity bounds in the difference case along similar lines (this in particular covers the case when for ). If or becomes large, then the defining equations of tend to become large as well due to Proposition 7.2. The hope is that a polynomial complexity bound can somehow be derived from this.

Acknowledgments. We are grateful to Grégoire Lecerf for a few helpful thoughts concerning section 8.

2.Abstract reductions

2.1.Definitions and basic properties

Let be a field and consider a linear map on some -vector space . A reduction with respect to is a linear map that satisfies:

R1

for all ;

R2

for all .

Such a reduction is said to be confined if

R3

is finite;

and normal if

R4

.

The condition R2 stipulates that is a projection. It will be convenient to also introduce the complementary projection involved in R1, called the remainder:

A normal reduction with respect to can be regarded as a projection of onto a subvector space of that is isomorphic to the cokernel of .

Proposition 2.1. If is a bijective -linear map, then a (normal, confined) reduction with respect to is also a (normal, confined) reduction with respect to and vice versa.

Proof. This follows from the observation that R1 and R4 only involve .

Remark 2.2. Sometimes, the linear map is only defined on a -vector space that contains . Our definitions easily extend to this case module the replacement of by in the conditions R1 and R4.

Proposition 2.3. Let be an algebraic extension of and consider a reduction with respect to the natural lift of . If for all automorphisms over and , then for all .

Proof. For any , we have for all automorphisms of , which implies .

2.2.Local reductions

Let be a projection and . Given a reduction with respect to a linear map , we say that the reduction is local for the projection if

(2.1)

for all .

Proposition 2.4. If is a local reduction with respect to for , then is a reduction with respect to .

Proof. Given , there exists a with , whence . Since , we also have . We conclude that .

Proposition 2.5. If is a reduction with respect to and a linear mapping with , then we define a local reduction with respect to for by

Proof. Given , there exists a with and . Furthermore, , whence and . We finally have .

2.3.Gluing local reductions together

Two projections are said to be orthogonal if . A family of projections is said to be orthogonal if its members are pairwise orthogonal. Such a family is said to be summable if is finite for all . In that case, the sum with is well defined for all . Moreover, is a projection that is orthogonal to each and any can be decomposed canonically as

Two local reductions for are said to be independent if and are orthogonal and for all , we have

(2.2)

Such reductions necessarily commute: we have , and in a similar way. A family of reductions for projections is said to be independent if the are pairwise independent and if is summable.

Proposition 2.6. Let be a family of independent local reductions with respect to and for the projections . Let and . Then we may define a reduction with respect to by

(2.3)

This reduction is local for . For all , we have

(2.4)
(2.5)

where are the indices with .

Proof. For all and , we get from (2.1) that , whence . It follows that , which proves (2.4). We also have , whence , and . It follows that . Since , this shows that is a reduction with respect to . We also have , so is local for .

Let us prove (2.5) by induction over . To simplify notations, assume that for all , and write for . For , we have nothing to prove. For , the result follows from (2.1). If , then let and . Notice that for . The induction hypothesis yields , whence

We conclude by induction.

2.4.Composition of reductions

Proposition 2.7. The composition of pairwise commuting reductions with respect to is again a reduction with respect to .

Proof. Let us first assume that . For each , we have . We also have . The general case follows using an easy induction on .

Let be two local reductions with respect to for two orthogonal projections . We have seen that Proposition 2.6 applies whenever (2.2) holds. If we only have for all , then the following still holds:

Proposition 2.8. If for all , then is a reduction with respect to .

Proof. For each , we again have . Furthermore, setting with , we have

It follows that .

2.5.Normalization of confined reductions

Consider a confined reduction with respect to . Then is a finite dimensional vector space and a finite dimensional subvector space of that we call the space of exceptional functions for . Notice that . Given any projection with and , let be defined by

(2.6)

Proposition 2.9. The relation (2.6) defines a normal confined reduction with respect to .

Proof. For any , we have

Secondly, for any , there exists a with , whence . For any , we thus obtain

Of course, we have . Finally, for any , the construction ensures that whence .

This proposition also admits a local analogue. Assume that we are given a local reduction with respect to for the projection . We say that is locally confined if is finite dimensional, and locally normal if . Assume that is locally confined and decompose with . Let a linear mapping such that is a projection of on (such a map is easily constructed from a basis for ). Then we define

(2.7)

Proposition 2.10. The relation (2.7) defines a local reduction with respect to that is locally normal and confined for .

Proof. Given , let . Then we have . Letting be such that , we also have , whence , since and coincide on . From , we conclude that , whence is a reduction with respect to .

The reduction is local for , since . By construction, . For and , we finally have .

3.Lagrange's identity and generalizations

3.1.The original differential case

Let be a differential field for the derivation and let denote the skew ring of linear differential operators over . For instance, one may take and . Recall that the adjoint of a differential operator is defined by

Given two indeterminates and , we write for the set of -linear combinations of products with . The following identity is due to Lagrange [66]:

Proposition 3.1. For any of order , there exists a

called the “bilinear concomitant”, such that

(3.1)

More specifically, we have

(3.2)

Proof. Let us first prove the existence of by induction on . If , then , and we may take . If , then we write with of order and . By the induction hypothesis, there exists a with

It follows that

We conclude by taking

From this relation, it follows by induction over that

which can be rewritten as (3.2).

3.2.The difference case

Assume now that is a difference field for the automorphism and let denote the corresponding ring of skew difference operators. For instance, one may take and or with . We recall that the adjoint of a difference operator is defined by

We also define the finite difference operator associated to by .

Proposition 3.2. For any of order , there exists a

such that

(3.3)

More specifically, we have

(3.4)

Proof. We again prove the existence of by induction on . If , then , and we may take . If , then we write with of order and . By the induction hypothesis, there exists a with

It follows that

We conclude by taking

From this relation, it follows by induction over that

which can be rewritten as (3.4).

3.3.Matrix versions

Let again be a differential field for and consider a not necessarily commutative differential algebra over . One important example of such an algebra is the algebra of matrices . Proposition 3.1 naturally generalizes to the case when admits coefficients in , modulo a few precautions: in this setting, it is important to distinguish between usual operators that operate on the left

and their adjoints that operate on the right:

Similarly, a bilinear operator acts via

where are the coefficients of . Similar precautions apply to the difference setting, and we have:

Proposition 3.3. Let be a differential algebra over a differential field for . For any of order , there exists a

such that

More specifically, we have

Proposition 3.4. Let be a difference algebra over a difference field for . For any of order , there exists a

such that

More specifically, we have

3.4.Twisting

Let be a differential field for and assume that is the logarithmic derivative of another element . In fact, we may allow to live in some abstract extension field of . Then we have for all and we call the twist of by . More generally, for any , we define the twist . One may check that and for all . Using twisting, it is possible to replace the operator at the right hand side of (3.1) by any first order differential operator:

Proposition 3.5. Let be a differential field for . For any and , there exists a

such that

Proof. Let be a formal solution of in some extension field of . Let and . Then Proposition 3.1 provides us with such that

Let be the twist of in by with . Then we conclude by taking .

This proposition again admits a matrix generalization: for , it suffices to replace and by matrices with invertible and . We also need to assume that commutes with each coefficient of .

One interesting special case occurs as follows: let , let be monic of order and assume that we wish to obtain a generalization of (3.1) with instead of in the right hand side. Then we lift into an operator by sending each coefficient to and take

We next apply Proposition 3.5 to obtain a with

This equation can be considered as the “desired” generalization of (3.1) with instead of in the right hand side.

As usual, Proposition 3.5 also admits a difference version. Let be a differential field for and let be such that . The twist of an operator by is again defined as , so that .

Proposition 3.6. Let be a differential field for . For any and , there exists a

such that

3.5.Another twisted variant

There are many more variants and generalizations of Lagrange's identity. Let us describe one more variant in the differential case. Let be a possibly non-commutative -algebra, where is a field.

Proposition 3.7. Let be a monic first order operator in and let be a monic operator of order . Then there exist differential operators such that

Proof. Writing for the derivative of a differential operator , and using the expansion

we notice that

Taking for , the result follows.

We typically read the above identity as

(3.5)

where plays a similar role as in (3.1), where plays the role of , and where the derivation got replaced by a general monic operator .

4.Differential reductions

4.1.Partial fraction decomposition and Hermite reduction

Let be a field with algebraic closure . Any rational function admits a partial fraction decomposition

(4.1)

where ,

and, for each ,

Given , we will write for the projection and we call the polar part of at . The number stands for the order of the pole of . We also write for the projection and call the constant part of .

The decomposition (4.1) has the advantage of being symmetric with respect to all points in . Often, the polar part at infinity is combined with the constant part, in which case (4.1) becomes

(4.2)

where . More generally, for subsets , it is convenient to denote

We also denote and write for the projection . Notice that forms a ring for any .

The Hermite reduction of with respect to the ordinary differentiation with respect to is defined by

It is not hard to check that this indeed defines a normal reduction on with respect to . Its restriction to determines a normal reduction on with respect to , by Proposition 2.3.

4.2.Local reduction with respect to differential operators

Let us now consider an arbitrary non-zero linear differential operator instead of . For any , it is well-known [55] that there exists an indicial polynomial at and a shift at with

for all . Similarly, at infinity there exist and with

Given , let us now construct a local reduction at . Given , let and be such that

Setting , we define by induction over , as follows:

(4.3)

In a similar way, we may construct a local reduction at infinity. Given , let and be such that

Setting , we define by induction over :

(4.4)

Notice that for any .

Proposition 4.1. Let and . The mapping is a reduction with respect to such that:

  1. The reduction is local for .

  2. We have .

  3. The reductions and are independent for any with .

  4. If , then .

Proof. We adopt the notations from above. Whenever , we let if and if . Whenever , we also denote if and if .

By induction on , let us first show that there exists a with and ; notice that this in particular implies ii, as well as . If , then we may take . If and , then there exists a with and , whence . If and , then there exists a with and . Hence , so we may take with .

Again by induction on , let us next show that . If , then we have . If and , then implies . Consequently, . If and , then . This completes the proof that is a reduction with respect to .

Decomposing with , let us next show by induction on that . If , then we have and . If and , then we may decompose and with , after which . If and , then implies . It follows that . This shows i.

Now let be such that . The projections and are clearly orthogonal and it follows from ii that and ; this shows iii. The last fact iv also follows from ii.

Proposition 4.2. For each , let

(α∈𝕂¯)
(α=∞)

Then is included in the vector space spanned by . We also have for and .

Proof. Adopt the same notations as in the proof of Proposition 4.1. Given , let us prove by induction on that . This is clear if . If and , then the induction hypothesis yields , whence . If and , then the induction hypothesis directly yields . This completes the proof that . If , then we also have and . Since admits at most roots, the cardinality bounds follow.

Remark 4.3. Proposition 4.2 shows that is locally confined for all . This makes it possible to define a locally normal reduction using Proposition 2.10. In order to make this construction fully effective, we still have to determine the space . Let

(α∈𝕂¯)
(α=∞)

We claim that . Assume for contradiction that there exists a function with and chose such that is minimal with this property. Let as in the proof of Proposition 4.1. Since we must have , it follows that . Therefore provides a counterexample with , contracting our minimality hypothesis.

Having proved our claim, one may use linear algebra to determine elements with and such that form a basis for . Setting and , we define to be the “distinguished” projection of on such that for all . In order to apply Proposition 2.10, we finally construct the mapping by sending each to and then extend to by setting .

Since can be arbitrarily large, we notice that the determination of can be very expensive. In particular, there exists no polynomial time algorithm to compute them. We also notice that the construction of is symmetric in the sense that for all automorphisms of over .

4.3.Global reduction with respect to differential operators

From Propositions 4.1 and 2.6, it follows that we can glue the local reductions at the finite points together into a reduction , defined by

(4.5)

Notice that Proposition 2.3 implies that whenever . Using Propositions 4.1-iv and 2.8, we may next define the global reduction with respect to by

(4.6)

Notice that whenever .

Remark 4.4. Applying the above construction to the local reductions from Remark 4.3 instead of , we claim that we obtain a normal reduction. Indeed, given , we get for all , whence and . Since is locally normal, we conclude that . We notice that is again symmetric in the sense that . It is also easily checked that is locally confined for each projection .

It is interesting to consider the restriction of the global reduction to certain differential subrings of and . More precisely, given a subset of poles, the set forms a differential -subalgebra of . From Proposition 4.1-ii, it follows that is stable under reduction. If is the set of zeros of a monic separable polynomial , then we observe that coincides with . In that case, Proposition 2.3 implies that is stable under reduction as well. We notice that is finitely generated by as a -submodule of ; for this reason we call it a narrow submodule of .

Theorem 4.5. Let be of order and let be the global reduction with respect to . Assume that is the set of zeros of a monic separable polynomial of degree . Then the restriction of to is a confined reduction with respect to such that

Proof. Given and with , we notice that , whence . Given with , it follows that . In other words, is stable under reduction and so is , since the reduction commutes with all automorphisms of over .

Let us now turn to the dimension bound. Since is -linear, it suffices to show that . Let be defined as in Proposition 4.2 for all . Then the proposition implies that

contains at most elements. Now let , so that . We have

For any , we have by Proposition 4.1-ii, whence . Since by Proposition 4.2, it follows that , whence

Now by Proposition 4.1-ii and by Proposition 4.2. We conclude that , with .

Remark 4.6. A bit of flexibility is possible regarding the choice of in the definitions of the local reductions . Taking for as in [16] may lead to slightly more confined reductions in Theorem 4.5. Conversely, by taking somewhat larger, one may ensure that for all in Proposition 4.1-ii, so that and commute for all in Proposition 4.1-iii. This makes it possible to treat the singularity at infinity in a more symmetric way and replace definition (4.6) by

Remark 4.7. Since the reduction form Remark 4.4 is both normal and locally confined for each projection with , its restriction to yields a normal confined reduction.

4.4.Back to reductions with respect to

Assume still that is a non-zero differential operator. Let be its order, let be its dominant coefficient, and let be an arbitrary monic separable polynomial such that is contained in

Given formal indeterminates , the set

(4.7)

admits a natural -module structure for the derivation with

By construction, is a formal solution to in .

Let us show how to construct a confined reduction with respect to on . We start with the confined reduction with respect provided by Theorem 4.5. Now consider

with and let us define by induction on . If , then we take

By construction, we have , whence there exists a such that . Using Proposition 3.1, it follows that

Assume now that . Then we define

We notice that

It follows that

In fact, an easy induction on shows that

(4.8)

In particular, the constructed reduction is confined with :

Theorem 4.8. Let be of order and let be a monic separable polynomial of degree such that (4.7) admits the structure of a -module. Then (4.8) defines a confined reduction with respect to and we have

4.5.First order systems

Let us now turn to first order systems and study analogues for the theory and results from the previous subsections. This means that the operator from section 4.4 is now a first order operator in and we need to adapt the theory of sections 4.2 and 4.3 to the adjoint operator ; recall that operates at the right on matrices in . By Proposition 2.1, we also recall that the construction of a reduction with respect to is equivalent to the construction of a reduction with respect to , for any invertible matrix . In particular, for the construction of a local reduction at , it suffices to consider operators of the form with .

As the analogue for the operator from sections 4.2 and 4.3, we will therefore start with an operator and show how to construct the local reduction at with respect to . It will be convenient to first put into an even simpler form through multiplication by a suitable invertible matrix . Let . In analogy with [50], we say that an invertible matrix is a tail chopper for at if we either have and , or , , and for some :

T1

If , then and for the matrix formed by the first rows of .

T2

If , then and for the matrix formed by the last rows of .

Tail choppers can be computed using the following algorithm:

Algorithm 4.1

Input:

Output: a tail chopper for at

Let ,

Let ,

while and do

  • Let .

  • Using row sweeping, we determine an invertible matrix such that the first rows of and are the same, the first rows of have rank , and the last rows of are zero.

  • Set and next multiply the last rows of with .

  • Set and let and be the matrices formed by the first and last rows of .

return

Proposition 4.9. Algorithm 4.1 terminates and computes a tail chopper for at .

Proof. At every iteration of this loop, we notice that increases by one (as long as ) and it is easily checked that T1 and are loop invariants.

Remark 4.10. The above algorithm for the computation of tail choppers is somewhat reminiscent of Abramov's EG-elimination method [1].

Given a tail chopper as above, let be the denominator of (so that ), and consider the operator . We have

(4.9)

We consider as the “preconditioned” version of with respect to which we will now construct the local reduction . From (4.9) it follows that the shift of at is simply and its indicial polynomial is given by

By construction, is invertible as a matrix in , whence there exist at most values of where . We are now in a position to formulate a counterpart of (4.3).

Given , let and be such that

We define by induction over , as follows:

(4.10)

Although this definition works, it is slightly suboptimal in the sense that values of for which are “skipped” altogether even though the system might actually be solvable. Given an arbitrary scalar matrix , let denote a pseudo-inverse such that for any . Also let denote a matrix with for all . We may now replace (4.10) by

(4.11)

This “more confined” definition ensures that instead of . The local reduction at infinity can be defined in a similar fashion, via the change of variables .

Given an arbitrary first order operator , this completes our construction of the local reductions and with respect to . Proposition 4.1 and the theory from section 4.3 can be adapted mutatis mutandis. In particular, and are indeed local reductions for and that can be glued together into a global reduction with respect to . Assuming that the mappings and were taken to commute with all automorphisms of over , we again have for all . The analogue of Theorem 4.5 is:

Theorem 4.11. Let be of order one with invertible and let be the global reduction with respect to . Assume that is the set of zeros of a monic separable polynomial of degree . Then the restriction of the reduction to is a confined reduction with respect to such that

As to the counterpart for section 4.4, assume that with . Let be an arbitrary monic separable polynomial such that , where . Given a formal column vector with entries , the set

(4.12)

admits the natural structure of a -module for the derivation with . Given a row vector , we define the reduction on by

where is the confined reduction on with respect to as constructed above. By construction, there exists a with , whence

The further results from section 4.4 now adapt mutatis mutandis.

Remark 4.12. Starting from a reduction for with , one may also define a reduction for : given , there exists a with , and we take . Then there exists a with . Using Lagrange's identity the other way around, we obtain , whence . This shows that is indeed a reduction with respect to .

4.6.Reductions on D-modules with respect to general operators

We have seen how to construct reductions on with respect to general operators and how to construct reductions on D-modules with respect to . It is natural to ask whether we can construct reductions on D-modules with respect to more general operators.

Let us first show how to do this in the setting of first order systems from the previous subsection. More precisely, assume that , with and are as in the previous subsection. Given another matrix , our aim is to construct a reduction on with respect to . The idea is to use an explicit variant of Proposition 3.5. Given a matrix , we notice that

(4.13)

Now we may reinterpret as a first order differential operator of dimension on , for which we may construct a confined reduction using the theory from the previous subsection. Given , this means that there exists a with . Setting

(4.14)

it follows from (4.13) and the defining equation of that

This shows that (4.14) indeed defines a reduction with respect to on .

The above construction admits several variants. Let us briefly sketch what happens if the first order matrix operator is replaced by a monic operator . This time, we rather rely on Proposition 3.7 and more precisely on formula (3.5). Given a row matrix and using the fact that , we have

The idea is now to construct a confined reduction with respect to on . This can be done by generalizing the theory of sections 4.2-4.4 to allow for matrix coefficients, along the lines of section 4.5. We next define the confined reduction on as usual by . Given , there then exists a with , whence

We conclude that is a confined reduction with respect to on .

5.Difference reductions

5.1.The field as a difference field

Let be a field and consider a birational map from the projective line over into itself. Such an map is called a homography and is necessarily of the form

where the matrix is invertible. Given and , we will denote .

Given a homography , the field becomes a difference field for the automorphism that postcomposes with :

Notice that

for all , , and . More generally, if was already equipped with an automorphism , then we may extend into an automorphism on by taking

where acts coefficientwise on rational functions . Inversely, we must have for automorphisms of this kind, and we call the homography associated to .

Two important special cases are the shift operators and the -difference operator :

Notice that is a fixed point for for both of these examples.

In fact, we claim that the general case can essentially be reduced to one of these two special cases via a change of variables . Indeed, such a change of variables transforms difference equations in and into difference equations in and . Assuming that we are allowed to extend with the roots of the characteristic polynomial of , we may first put in Jordan normal form and then normalize it through division by . After this, we have for some or for some .

For what follows, we will always assume that has infinite order. In the case of a shift operator , this means that we should have ; for a -difference operator , the number should not be a root of unity.

5.2.Difference modules

One particularity of automorphisms of as above with respect to the ordinary differentiation is that the sets are no longer stable under . On the other hand, the set is stable under , where and . It will be convenient to expand somewhat more on this observation and introduce a few more notations in addition to those from section 4.1.

Let us first introduce truncated analogues of the sets with : given , we define

The definitions extend to the case when by taking . Given and a map , we also denote

Each of these vector spaces also come with projections similar to the and that were already defined before. For example, stands for truncation at order of . For , it is also convenient to introduce .

Given , we write if there exists an such that . This clearly defines an equivalence relation. We also write if there exists an such that . Notice that is a partial ordering on by our assumption that has infinite order. Given , we introduce the subsets , and of by

For , we understand that , , etc. We notice that is a -module, whereas is a -module. Since is a fixed point for , then we notice that is even a -module.

More generally, let be an increasing mapping in the sense that for all . Then is a -module. If , then the set is actually a -module. We say that is symmetric if it is stable under all automorphisms of over . Similarly, is symmetric if is symmetric and commutes with all automorphisms of over . In that case, the intersection is a -module and a -module, whenever .

Assume that and consider a -submodule of with . It is not hard to check that such a module is necessarily of the form for some increasing . We say that is narrow if it is a finitely generated as a -module. More generally, a finitely generated -module of of the form is said to be narrow. Such a narrow submodule of is always the intersection of with a narrow submodule of for some symmetric .

Consider a narrow submodule that is finitely generated by . Without loss of generality, we may replace each by the finite set and remove any elements such that for some . This means that we may assume that is finite and for some function . Now consider the map defined by

Then coincides with the -module generated by , so that .

5.3.Local and semi-local reduction with respect to shift operators

Let us now move our attention to a shift operator , where with . We assume that with both and . We took an operator with respect to rather than since, for the applications that we have in mind (see Proposition 3.2), the operator will actually be the adjoint of an operator in .

Since , the action of at infinity can be approximated to any order by a differential operator in . In particular, there still exist an indicial polynomial and a shift (in fact, , since admits coefficients in ) such that

Consequently, we may define the local reduction at infinity with respect to in a similar way as in subsection 4.2.

Our next aim is to define a “semi-local” reduction at . Since admits a singularity at if and only if admits a singularity at for all , we also need to take into account the behaviour at each of these shifted singularities. For this reason, the reduction will only be “semi-local”.

Given , we define the -span of to be the largest integer such that . If no such integer exists, then we set . Assuming that , consider

If , then the orders of and at coincide,

and . If , then it follows that

This allows us to define by induction on as follows:

(5.1)

We call the downward semi-local reduction with respect to at .

Proposition 5.1. Let and . The mapping is a reduction with respect to such that:

  1. The reduction is local for .

  2. We have .

  3. The reductions and are independent for any with .

  4. If , then .

Proof. If , then we have already noticed that is a reduction with respect to such that i and ii are satisfied. If , then the proof follows the same scheme as in Proposition 4.1; for the sake of completeness, we give it in full. Throughout the proof, we denote and .

By induction on , let us first show that there exists a with ; notice that this in particular implies ii and also . If , then we may take . If and , then there exists a with , whence . If and , then there exists a with , where . Hence , so we may take .

Again by induction on , let us next show that . If , we have . If and , then implies . It follows that . If and , then we have . This completes the proof that is a reduction with respect to .

Decomposing with , let us next show by induction on that . If , then we have and . If and , then we may decompose and with , after which . If and , then implies . It follows that . This shows i.

Now let be such that . Then the projections and are clearly orthogonal and it follows from ii that and ; this shows iii. The next fact iv also follows from ii.

Proposition 5.2. Let and .

  1. We have , where

    (5.2)
  2. For some , we have , , and .

  3. We have .

  4. We have .

Proof. Assume the notations from the previous proposition. Let us prove i by induction on . If , then . If and , then ; since the induction hypothesis implies , we obtain . If and , then the induction hypothesis yields .

Let us next prove ii, again by induction over . If , then we may take . If and , then let be such that , , and . Then we have , , and . If and , then let be such that , , and . We have . Since , we have , , and . It thus suffices to take .

With as above, we have , which implies iii. Similarly, yields iv.

5.4.Aligned downward reduction with respect to shift operators

Let us now turn to the problem of combining the semi-local reductions where for some subset . Given a general subset , the set admits the partitioning

The notion of -span can also be generalized:

In this subsection we start with the case when is aligned in the sense that with . In that case, we notice that

(i=1,…,t-1)

Although the construction from the previous subsection provides us with a local reduction for the projection , this reduction is sometimes suboptimal due to the fact that it may take a long time to reduce functions for which is small. For this reason, we introduce an alternative semi-local reduction by

We call the downward reduction with respect to for . It is “less confined” than but relies on the reduction in order to reduce the polar parts for the singularities in .

Proposition 5.3. Let be aligned with and . Then the mapping is a reduction with respect to that satisfies:

  1. The reduction is local for .

  2. We have .

  3. The reductions and are independent for any aligned with .

  4. We have .

Proof. Let us denote and for , so that . Also let with and for .

Let us first show that is a reduction with respect to . We clearly have . By Propositions 5.1-ii and 5.2-iii, we notice that for , whence and . Then equals .

The reduction is also local for since . The relation shows that . The properties iii and iv are shown in a similar way as for Proposition 5.1.

Proposition 5.4. Let be aligned with and .

  1. We have , where , using the notation from (5.2).

  2. For some , we have and ().

  3. We have .

  4. We have ().

Proof. With the same notations as in the previous proof, Proposition 5.2-i implies for , whence . This proves i.

Proposition 5.2-ii implies the existence of with , , and , for . In particular, this yields . Now consider with . For , we conclude that .

As to iii, we have , thanks to Proposition 5.2-iii. Proposition 5.2-iv similarly implies iv, since for .

Corollary 5.5. Let be aligned with , let , and consider the narrow -module . Then is stable under .

Proof. Let and decompose it as as above. From Proposition 5.3-ii, it follows that . Notice that if and only if for . Assuming that that , it follows from Proposition 5.4-iv that for , whence .

Remark 5.6. In analogy with Remark 4.3, one may define a local normalization of . This time, it follows from Proposition 5.4 that the space is spanned by . Again the computation of a basis of can be expensive and there exists no general polynomial time algorithm for doing so.

5.5.General downward reduction with respect to shift operators

Let us now assume that is such that is finite for every . The equivalence relation restricted to then leads to a natural partitioning

(5.3)

where each component is aligned, and whenever . This allows us to define

By Propositions 5.3-i, 5.3-iii and 2.6, we get that is a local reduction for with respect to . From Propositions 5.3-iv and 2.8, we deduce that is a reduction with respect to . We call and the downward reductions with respect to for and . If is symmetric, then and , by Proposition 2.3.

Remark 5.7. When taking to contain a section of with respect to the equivalence relation , we ensure that . Nevertheless, each can only contain “half” of the set , since is assumed to be finite. This means that cannot cover , so cannot truly qualify as a “global reduction”. In order to obtain a global reduction, one needs to combine downward and “upward” reductions (constructed from downward reductions with respect to the operator ). Even then, the obtained “global” reduction crucially depends on the choice of . In fact, this kind of reductions will not be useful for our application to creative telescoping; we rather need reductions on suitable narrow submodules.

In order to obtain a counterpart for Theorem 4.5, let us now turn our attention to restrictions of downward reductions to suitable narrow submodules of .

Theorem 5.8. Let be of order . Consider a narrow -submodule

of , where is a symmetric finite set and for some symmetric function . Let , , and . Then the restriction of to is a confined reduction with respect to with

Given , there exists a with .

Proof. Let be partitioned as in (5.3). Let us first show that the narrow -submodule

of is stable under downward reduction . The module is clearly stable under the projection for any subset . Since is local for , we also have stability under the reduction . Now given and , we have with and by Corollary 5.5. For any , it follows that and . Since commutes with all automorphisms of over , we also get the stability of under .

Given , let us next prove the existence of some with . For each , setting , Proposition 5.4-ii yields an element with and for all , whence . At infinity, there also exists a with , and . Taking , it follows that and

If and admits coefficients in an extension field of of degree , then we may replace by , since and for every automorphism of over . Notice that this shows in particular that is a reduction on .

Let us finally show that

(5.4)

which also implies the bound on . For each , let be as in (5.2) and notice that satisfies

Given , let

We also define ,

and . Now given , Proposition 5.3 implies for all , where . It follows that . We also have , whence . Since , we conclude that (5.4) indeed holds.

Remark 5.9. When using the locally normal reductions from Remark 5.6 instead of the reductions and similarly for the local reduction at infinity, it can be checked that the resulting reduction in Theorem 5.8 is both confined and normal.

5.6.Back to reductions with respect to

Assume now that is a difference operator of order in with . Consider a narrow -submodule of , where is a symmetric finite set, a symmetric function, and . Assume that for every zero of of multiplicity , we have

(5.5)

where we understand that if .

Given formal indeterminates , we claim that the set

(5.6)

admits the natural structure of a difference module over for the shift operator with

Indeed, given , this forces us to take

and it suffices to check that . Now for any root of of multiplicity , we have , whence , thanks to our assumption (5.5). We conclude that , as desired. Notice that is a formal solution to the equation in .

Let us now show how to construct a confined reduction with respect to on . Let denote the restriction of the directed reduction at with respect to as constructed in the previous subsection. Now consider

with and let us show how to define by induction on . If , then we take

By construction, we have , whence there exists a such that . By Theorem 5.8, we actually have . Using Proposition 3.2, it follows that

Assume now that . Then we define

We notice that

By what precedes, it follows that

In fact, an easy induction on shows that

(5.7)

In particular, the constructed reduction is confined with :

Theorem 5.10. Let be of order with and consider , , and as above. Let , , and . Then is a confined reduction such that

Remark 5.11. The smallest set and map that satisfy our hypotheses are given by and . Assume that and are taken this way. If for all , then it follows that for all , whence . However, in the other extreme case when is aligned, we only have for all and the growth of can be quadratic in . This larger growth is due to the fact that we required functions of small span to be reducible fast in section 5.4. If we drop this requirement, then linear growth can be restored by changing the definition (5.1) to

(5.8)

where and . Compromises between both definitions are also possible, which should ultimately make it possible to systematically gain a factor with respect to some of the complexity bounds later in this paper. However, the details are technical, so we reserve them for a future work.

5.7.The case of -difference operators

The theory of the previous section naturally adapts to -difference operators (where is not a root of unity), except that both zero and infinity now need to be treated apart.

So let us start with the construction of the local reductions and for . This time, the behaviour of at zero and at infinity is given by

where and (in fact, we have and ). Since is not a root of unity, we notice that admits at most one solution for any .

Now consider . Let and be such that

Setting , we define by induction over :

(5.9)

Similarly, assume now that and are such that

Setting , we define by induction over :

(5.10)

Notice that for any .

The local reductions at other points are constructed in a similar way as in subsection 5.3. It is not hard to check that the with are indeed local reductions that satisfy similar properties as in Proposition 5.1. Given a subset , this allows us to glue them together into a reduction in a similar way as in subsections 5.4 and 5.5. Mutatis mutandis, this leads to the following analogue of Theorem 5.8:

Theorem 5.12. Let be of order . Consider a narrow -submodule

of , where is a symmetric finite set and for some symmetric function . Let , , and . Then the restriction of to is a confined reduction with respect to with

Given , there exists a with .

The results from subsection 5.6 also naturally adapt to the -difference setting.

5.8.First order systems

Let us now outline how to adapt the theory of this section to first order systems. We assume that for ; the case of -difference operators can be treated in a similar way. Concerning the analogue of section 5.3, let us start with an operator , where is a monic separable polynomial and where is an invertible matrix. Recall that operates at the right on matrices in .

The local reduction at infinity is constructed in a similar way as for differential operators. The semi-local reductions at other points are somewhat easier to construct in the sense that there is no need for any auxiliary tail choppers. Given , we define by induction over :

(5.11)

As in section 4.5, this definition can be further optimized through the use of pseudo-inverses:

(5.12)

The remainder of the theory and result from sections 5.3, 5.4, 5.5, and 5.6 can be adapted mutatis mutandis. In particular, for every symmetric finite set , we obtain a reduction that satisfies the following analogue of Theorem 5.8:

Theorem 5.13. Let be of order be of order one. Assume that is invertible and that for some monic . Consider a narrow -submodule

of , where is a symmetric finite set and for some symmetric function . Let , , and . Then the restriction of to is a confined reduction with respect to with

Given , there exists a with .

6.Creative telescoping

6.1.D-finite ideals of DD-operator algebras

Let be a field that is equipped with pairwise commuting operators , together with derivations and automorphisms such that for each we have either and or and . Given such a field , we will say that the skew polynomial ring

is a DD-operator algebra. In practice, we usually have for some constant field , and each is given by either one of the formulas

Modulo a homographic change of variables (and quadratic algebraic extensions), we have seen in section 5.1 that we may further reduce to the case when each is either a shift or a -difference operator. DD-operator algebras of this special type will be called standard.

A (left) ideal of is said to be D-finite if is finite dimensional as a vector space over . Given an -module and a “function” , we say that is D-finite if its annihilating ideal

is D-finite. We say that an ideal of is reflexive if for every and , we have .

It is well known [40, 60, 71, 64, 26, 70] that the theory of Gröbner bases generalizes to skew polynomial rings such as . As usual, this assumes a total ordering on the set of “monomials”

that refines the divisibility relation on . Let denote this divisibility relation. Given a finite subset of , a Gröbner basis for the left ideal generated by this subset can be computed using a non-commutative version of Buchberger's algorithm. The initial segment of reduced monomials for under the Gröbner stairs for admits the usual property that forms a basis for the vector space . Given , the reduction of with respect to yields with .

6.2.Matrix representations of operators

Let be a reflexive D-finite ideal of . Consider any basis of and let be the column vector with entries (it will be convenient to call such a column vector a basis of as well). Given an operator , the action of on gives rise to a matrix with

(6.1)

In particular, for each there is a matrix with

(6.2)

If are such that , then we may compute the -th column of by reducing with respect to and writing the result as a linear combination of .

Proposition 6.1. If is an automorphism, then is invertible.

Proof. Assume for contradiction that there exists a non-zero row vector with . Then we would get . Writing for the column matrix with entries , we thus obtain and , which yields a non-trivial relation between the .

Proposition 6.2. For any with , we have:

i=∂i,θj=∂j)
ii,θj=∂j)
ii,θjj)

Proof. If and , then (6.2) yields

Since both derivations commute and the entries of form a basis of , this yields the first relation. The two other relations are proved in a similar way.

Given , we say that is a cyclic vector for if form a basis for . If and contains an element with , then it is well known that such a cyclic vector always exists and can be computed [3, 25, 15]. These results can naturally be adapted to the case when is a shift operator or a -difference operator such that contains an element such that is infinite. We will say that acts non trivially on if we are in either one of these situations. In the terminology of [6, section 5.3], this implies that is simple, whence [6, Corollary 5.3.6] yields a way to compute cyclic vectors. Given a cyclic vector , the left ideal of of operators that annihilate is principal and its unique monic generator has order . In terms of matrices, these properties may be restated as follows:

Proposition 6.3. Let be such that acts non trivially on . Then there exists a basis of with entries such that

The operator is the unique monic operator of order with .

6.3.Narrowness

Let us now consider a reflexive D-finite ideal of a DD-operator algebra

with one additional operator . We abbreviate , , , and assume that either

where is the constant field of . Moreover, if , then we require that is either a shift operator with or a -difference operator , where is not a root of unity. As usual, we set , where we understand that if . We also write if and if .

Recall that a -submodule of is said to be narrow in if it is stable under and finitely generated as a -module. Given such a module , we have shown in section 5.2 that there exists a finite symmetric set and a symmetric such that . Let be a fixed basis of . A -submodule of is said to be narrow in if it is of the form for a narrow submodule of .

Given a narrow submodule of as above and , we have seen in sections 4.5 and 5.8 how to construct a computable confined reduction with respect to the first order operator as well as a corresponding computable confined reduction with respect to . Alternatively, one may rely on the constructions of reductions with respect to scalar linear differential and difference operators , as we will detail now.

Let be the entries of the basis . We say that is cyclic if for . Such a basis always exists by Proposition 6.3 and it comes with a unique monic operator of order with . Multiplying with the lcm of the denominators of its coefficients, we obtain a new operator with and .

Assume that and that we are given a narrow submodule in as above with respect to a cyclic basis . Then we have shown in sections 4.2, 4.3 and 4.4 how to construct a computable confined reduction with respect to as well as a corresponding computable confined reduction with respect to .

Assume now that and that we are given a narrow submodule in as above with respect to a cyclic basis . Assume also that for all , where denotes the set of fixed points of . Then we have shown in section 5 how to construct a computable confined reduction with respect to as well as a corresponding computable confined reduction with respect to , where . For , we notice that for some -vector space of dimension at most . Given with and , we may write with and for . By induction on , we may then define a confined reduction on by . We have .

Assume that there exists a -submodule of that comes with a (computable) confined reduction with respect to . Then we say that is telescopable in for the chosen basis , and we call the associated reduction. By what precedes, this is the case whenever there exists a -submodule of that is narrow in .

6.4.Creative telescoping

With the notations from the previous subsection, given a function in some -module, the set of telescopers for with respect to is defined by

From this definition it is immediately apparent that forms an ideal of . Given a telescoper , an element such that is called a certificate for .

Assume now that is telescopable in , with associated reduction . We say that a function is primitive if for each , there exists an with . We write for the set of such functions. Given and an operator with , we notice that implies the existence of some with , whence and . In particular, given an operator with and , we have

(6.3)

Since lives in the finite dimensional -vector space , this yields a way to determine telescopers through linear algebra, simply by searching -linear relations between the where runs over . For instance, taking , one obtains a telescoper in , which generalizes the introductory example (1.4). Doing this for each shows in particular that is D-finite.

In our DD-algebra setting, it is also natural to perform the linear algebra incrementally following the term order on , as in the FGLM algorithm [37]. This leads to the following algorithm to compute a D-finite ideal contained in .

Algorithm 6.1

Input: a telescopable D-finite ideal with associated reduction and

Output: the reduced Gröbner basis for a D-finite ideal included in

monomials remaining to be treated
Gröbner basis being constructed
monomials under stairs

while do

Let be minimal for and set

if is not reducible with respect to then

Compute

(If , then we may also take )

if , then

else

return the reduced Gröbner basis for the ideal generated by

Theorem 6.4. Algorithm 6.1 is correct, terminates, and returns the Gröbner basis of a D-finite ideal with

(6.4)

Moreover, if is normal, then .

Proof. Throughout the algorithm, we observe that the elements with are -linearly independent. Consequently, never exceeds , which ensures termination of the algorithm. By (6.3), it is also clear that we only insert telescopers in to . Let us now prove the bound (6.4).

Throughout the algorithm, it is easily verified that contains only monomials that are smaller than . It follows that the leading monomial of is , when inserting a new element into . This in turn implies that, at the start of each iteration of the main loop, the set contains precisely those monomials below that cannot be reduced with respect to . At the end of the loop, this means that consists exactly of those monomials that cannot be reduced with respect to . Consequently, .

Assume finally that is a normal reduction. Then it is straightforward to verify that (6.3) becomes an equivalence. In particular, given a telescoper of the form , we get that , whence for all , since the elements with are linearly independent. In other words, , whence . At the end of the main loop, this also means that the S-polynomials of any two elements in necessarily reduce to zero; in other words, is already a reduced Gröbner basis, so we may directly return it after the main loop.

Remark 6.5. The important property of is that . If , then it follows that we may also take instead of , since

where and .

Remark 6.6. Reduction-based algorithms for creative telescoping do not require the computation of certificates. Nevertheless, it is not hard to modify the algorithms such that certificates are computed along with the telescopers themselves.

Remark 6.7. It is possible to generalize the theory of this section to the case when is replaced by a more general operator that commutes with . In the differential case, it essentially suffices to replace the confined reduction with respect to by a confined reduction with respect to , as constructed in section 4.6. In the difference case, such reductions can be constructed along similar lines.

7.D-finiteness tests for telescoping ideals

7.1.Location of singularities

Adopt the same notations as in section 6.3. Let be a fixed basis of with entries . Given an operator , we define its order at in and with respect to the basis by

and its (finite) set of singularities in by

Thanks to Proposition 6.2, it turns out that is closely related to for all .

Proposition 7.1. Assume that and let . Then

i=∂i)
ii)

Proof. Assuming for contradiction that and , we have , whence Proposition 6.2 implies

Similarly, if and , then we have , whence

Proposition 7.2. Assume that and let . Let

ii)
ii)

Then for any with , there exist integers and with

Proof. Since , we notice that

Assume that and . Then Proposition 6.2 implies

as well as

whence . Assuming that , it follows by induction on , we get that implies . Since is finite, we must have for some . If , it follows similarly that implies , whence the existence of some with .

Let us next consider the case when and . Then Proposition 6.2 implies

whence and . We conclude in a similar way as above.

7.2.Orders at singularities

Given a set of operators , we define

Let be the monoid generated by the with and . Similarly, let denote the multiplicative monoid generated by the such that .

Proposition 7.3. Assume that and let with be such that for with . Then we have

Proof. Let be sufficiently large such that for all with . This implies that . Let . Since and is invertible, we get

(7.1)

On the other hand, there are constants with

Since for , we have in the above formula. By our choice of , we also have , whence

Combined with (7.1), this yields

Since , we have and similarly .

Proposition 7.4. Assume that . Let and with . Then

Proof. Let be sufficiently large such that . Since and is invertible, we have

Since , we also have

The combination of both formulas yields

We conclude in a similar way as in the proof of Proposition 7.3.

7.3.Explicit telescopability

Proposition 7.5. Assume that . If is finite, then is telescopable in .

Proof. Assume that is finite and let , where . Notice that is stable under both and . Given , Proposition 7.1 implies that . Given , this means that

i=∂i)
ii)

This proves that the -module is narrow, as required.

Example 7.6. If is a standard DD-operator algebra with , then is finite if and only if does not depend on for every with . In particular, we may take in the proof.

Proposition 7.7. Assume that and let . Assume that the following holds for every :

Then is telescopable in .

Proof. We define the function by

Since is finite, we may also define by

By construction, this function is increasing. We also notice that

is finite, with , and the restrictions and of to and satisfy . It follows that is a narrow -module of . In particular, the -module comes with a reduction with respect to .

We claim that is stable under . The stability under is clear. Let , and with , so that . Then for any , we have

In other words, , which completes the proof of our claim.

Given an iterated derivative and , we next observe that there exist integers with

where we notice that for all . For any , it follows that

By Proposition 7.3, there are only a finite number of points where and at these points. In other words, setting , there exists a finite dimensional -vector space with . We extend to a confined reduction with respect to by setting for all and . Notice that is still stable under , since commutes with .

Example 7.8. Consider the case when is a standard DD-operator algebra with . If and , then it follows that the denominator of as a rational function in should be a power product of linear forms with , and . If this is not the case, then one may try replacing the operator by its inverse , for which one might have more luck. If either or is a -difference operator, and the other one is a shift operator, then the denominator of should be a polynomial in . If both and are -difference operators, then the denominator of should be a power product of terms with , , and .

We say that is explicitly telescopable in with respect to the basis if it satisfies the sufficient conditions from Propositions 7.5 or 7.7. These conditions present the advantage that they are formulated exclusively in terms of the singularities of . For a fixed basis , they are also easy to check. It remains to investigate how to compute a basis with respect to which is explicitly telescopable, whenever such a basis exists.

We will restrict our attention to the case when ; an algorithm can probably be developed for the difference case as well, but we will leave this for a future work. Given a singularity , let us show how to compute a basis with whenever such a basis exists. We outline the algorithm over , but it is not hard to adapt it to work over . The computed basis is no longer cyclic; with additional effort, it can be shown that there also exist cyclic bases with .

We start with the computation of a basis of formal local transseries solutions of the system at , given in the form of a fundamental matrix . If there exist local solutions that are not Laurent series, then the singularity clearly cannot be removable. Otherwise, the truncation of at a sufficiently high order gives rise to a basis such that admits no singularity at .

Unfortunately, may introduce new singularities, so this algorithm needs to be tweaked a little more. We replace with a new approximation at a order that can be written as a product of diagonal matrices with entries of the form with , invertible constant matrices in , and upper/lower triangular matrices of the form . It is not hard to compute such an approximation. By construction, its inverse admits a factorization of the same form, whence it introduces no new singularities.

Even small equations such as may admit solutions with a high valuation . It follows that the above algorithm for changing bases may be very expensive and its cost is generally not polynomial in the input size.

7.4.Testing D-finiteness

If is explicitly telescopable in with respect to the chosen basis and associated reduction , then the telescoping ideal is certainly D-finite for any . Moreover, Algorithm 6.1 allows us to efficiently compute a Gröbner basis for a D-finite subideal of ; with more computational effort, one may even compute a Gröbner basis for itself.

If is not explicitly telescopable in , even after a suitable change of basis, then telescoping ideals are generally not D-finite, although they exceptionally might be. In this subsection, we will show how to modify Algorithm 6.1 in order to test D-finiteness of and to compute a Gröbner basis for if so. We will restrict our attention to the case when , but a similar approach is likely to work in the difference case as well. It will also be convenient to assume that for any and with , we have for all . This is the case for standard DD-operator algebras.

Assume that we have fixed a cyclic basis of with entries and let be the vanishing operator of in . Denote , , , and . Then is a -submodule of that contains , but that is not necessarily narrow in .

Let be the normal global reduction on with respect to as constructed in Remark 4.7. The construction from section 4.4 still works and leads to a corresponding normal reduction on . The restrictions and of these reductions to and are also normal. We call the associated normal reduction to (and for the chosen basis ). For any subset , we notice that . In particular, for all .

Lemma 7.9. Let and . Let be the entries of the first row of . Then , where

and .

Proof. Let be the row vector with entries . Then implies . The result now follows from (4.8).

Proposition 7.10. Let be the set of all such that there exists an with and . Let with be such that . Then the -vector space generated by all with has infinite dimension.

Proof. Let be such that and let be such that and . Modulo the replacement of by , we may assume without loss of generality that for all . Recall that are pairwise distinct.

Consider the sequence with and for all . Let us show by induction on that and for all . These assertions hold by assumption for , so assume that . By Lemma 7.9, we have

For each , the assumption that therefore implies . Assume for contradiction that . Let be such that

Then we have

and Lemma 7.9 implies

Since and , it follows that . From , , and the normality of , we also get that . We conclude that , which contradicts the induction hypothesis.

We claim that are -linearly independent. Indeed, given a non-trivial relation with , we would obtain the contradiction

We conclude by noticing that .

We are now in a position to state the adaptation of Algorithm 6.1 to the case when the ideal is not necessarily telescopable.

Algorithm 7.1

Input: a D-finite ideal with associated normal reduction and

Output: the reduced Gröbner basis for if is D-finite and otherwise

, ,

while do

Let be minimal for and set

if is not reducible with respect to then

if then let

else decompose and let

Let be such that and let be as in Proposition 7.10

if then return

if , then

else ,

return

Theorem 7.11. Algorithm 7.1 is correct and terminates.

Proof. Let us first show that is finite. Given , we either have and , or for some . In the latter case, our assumption implies the existence of a maximal with . We claim that there also exists a maximal with . Our claim clearly implies that is finite, since we are left with a finite number of possible for each .

Assume for contradiction that the claim does not hold. Then , whence there exists a with . For each , there exists an with . For some and , it follows that . Since , we have and acts on as . In particular, there exist minimal and maximal with . For , we conclude that there cannot exist an with .

Having shown that is finite, we either hit some with , in which case the correctness follows from Proposition 7.10, or all computations take place in the narrow submodule of , and the correctness and termination are proved in a similar way as for Theorem 6.4.

8.Complexity of rational function arithmetic

8.1.Complexity of arithmetic in

Consider a standard DD-operator algebra with , as in section 6.1. We will use the algebraic complexity model: running times are measured in terms of the required number of field operations in and space complexity in terms of the required number of coefficients in .

We use a dense representation for polynomials in . A polynomial of total degree thus requires space . From now on, we assume that the dimension is a fixed constant, whereas may become large. Given of total degrees and , it is well known for that their product can be computed in quasi-linear time ; see [53] for a particularly efficient algorithm. Partial derivatives and -differences can clearly be computed in linear time . It is also well-known [11] that partial shifts reduce to coefficientwise univariate multiplications, so they can be computed in time .

The deterministic computation of the greatest common divisor of two polynomials is more expensive. The best currently known algorithms first select a principal variable, say , and reinterpret and as univariate polynomials in . They next recursively factor out the content and then compute the univariate gcd using a division-free algorithm for the computation of subresultants such as Berkowitz' algorithm [8]; see also [69]. This leads to the complexity bound with . When using randomized algorithms of Las Vegas type, one may take instead; since , this means that gcds can be computed in quasi-linear time (more precisely: after a generic linear change of variables, both and can be assumed to be monic in , and using evaluation/interpolation at points, one reduces to fast univariate gcds [18]).

Rational functions in are represented as fractions with , , and monic with respect to a suitable term ordering. We write for the degree of such a rational function. Given and , we notice that . One may compute in non simplified form in quasi-linear time , whereas the simplification of the resulting fraction requires further operations. Of course, it is possible to delay the simplification of fractions in many cases.

When computing with vectors of matrices with entries in , we always assume that denominators have been factored out. In other words, a matrix is represented as a matrix divided by a monic polynomial such that the gcd of and all entries of equals one (here is understood to be monic with respect to a suitable monomial ordering on ). We also define . Simplifying a general fraction can be done in time . Ignoring the cost of the final simplification, two matrices can be added in time and multiplied in time if . Here denotes the “exponent of matrix multiplication”, i.e. two matrices in can be multiplied using operations in .

When allowing for probabilistic algorithms (of Las Vegas type), it is interesting to consider the alternative straight line program (SLP) representation for polynomials and rational functions [20]. In this model, the length of the SLP that evaluates is a natural mesure for the complexity of . Using sparse interpolation [7, 58, 57, 51, 5], the actual coefficients of a polynomial of degree bounded by can be retrieved from its SLP representation in expected time . Similarly, the dense representation of a rational function can be retrieved in expected time ; see [59].

Algorithms that rely on SLP representations and sparse interpolation are particularly interesting when intermediate expressions swell in size, but the end-result is small in comparison. It turns out that the sizes of SLPs for creative telescoping remain reasonably small, which makes this approach relevant in our context. For this reason, we will also analyze the complexity of our algorithms in this model.

8.2.Complexity of arithmetic in

Now consider a standard DD-operator algebra extension of as in section 6.3. Polynomials in and rational fractions in can be represented in a similar way as in the previous subsection, but for our complexity analysis, it is useful to distinguish between the degree in and the total degree in . As in the case of vectors and matrices, we also assume that the denominators of polynomials in have been factored out.

With the above conventions, the monic gcd of two polynomials can be computed in time and we have

(8.1)

In fact, for any monic polynomial that divides , we even have

(8.2)

It is also possible to compute the corresponding cofactors with

still with the same complexity. These cofactors satisfy similar degree bounds

(8.3)
(8.4)

in addition to the usual bounds and in .

Similarly, given two polynomials with monic , the Euclidean division of by yields a relation

in which the quotient and the remainder satisfy the degree bounds

(8.5)
(8.6)

Setting , the division itself can be done in time .

8.3.Complexity of arithmetic in

The computation of the local reductions in sections 4 and 5 requires us to compute with roots of the leading coefficient of our differential or difference operator . If the factorization of is known, then this comes down to computations over the extension field , where is the monic irreducible factor of with .

Let us first examine the cost of arithmetic in such an extension field and denote . Elements in are represented by their pre-images with and . It is natural to define the total degree of in by . Given , it follows from (8.6) that

where we recall from (8.2) that . Similarly, the inverse of a non-zero satisfies

thanks to (8.3). Since additions and subtraction do not require any Euclidean divisions, we also have

The computations of and can respectively be performed in time and , whereas can be computed in time . One may take in the first two complexity bounds if we only require unsimplified results.

For simplicity, we will assume that the irreducible factorization of is known in what follows. In order to be complete, we still have to discuss the cost of computing such a factorization. Since this task can be rather expensive, it is actually better to avoid it by relying on the strategy of “dynamic evaluation” [31].

Starting with the square-free part of , the idea is to directly work in the extension as if it were a field. Whenever we hit a non-trivial zero-divisor during the computation of an inverse in , we simply abort all computations and resume with and in the role of . This can happen at most times, so all complexity bounds that will be proved in the sequel have to be multiplied by in order to take into account the cost of the factorization of . In fact, the overhead can be reduced at the expense of technical complications; since we merely want to establish polynomial complexity bounds, we can spare ourselves this effort.

Computations with algebraic functions in are often more efficient in the SLP model. In that case, at every evaluation point , we work in the algebraic extension of instead of . Here we used angular brackets for the evaluation of at in order to avoid confusion with evaluation in . Now SLPs over this algebraic extension of can be reinterpreted as SLPs over : additions/subtractions in translate into additions/subtractions in , whereas each multiplication or division in give rise to ring operations in . When implementing dynamic evaluation, one also must be careful to “take the same branch of the computation tree” at each evaluation point: whenever we hit a non-trivial zero-divisor of , we store the corresponding SLP for , and first relaunch the computation with in the role of (for all evaluation points), and then once more with in the role of (again for all evaluation points).

8.4.Partial fraction decomposition

Consider a fraction with and . Assume that is monic and that its factorization over is known, say . Let us analyze the cost of computing the partial fraction decomposition of .

The partial fraction decomposition can actually be computed over as well as over its algebraic closure. When working over , the partial fraction decomposition can be written at choice in one of the following two forms

where and . It is well known [65] that partial fraction decompositions of these kinds can be computed in quasi-linear time in terms of operations over .

In this paper, we have chosen to systematically work over the algebraic closure of . Now given a root of and any automorphism of over , we notice that . It therefore suffices to compute for only one root in each conjugacy class. We define the partial fraction decomposition of over to consist of and the collection of coefficients , where runs over and one element in each conjugacy class among the set of roots of , and where . The “tangling” and “untangling” morphisms from [52, Section 4] allow for conversions between partial fraction decompositions over and in quasi-linear time. It will be convenient to write

for the maximal degree in of a coefficient in the partial fraction decomposition of over .

We have the following deterministic (but rather crude) bounds for the cost of partial fraction decomposition and the resulting degree swell in .

Proposition 8.1. Let be as above. Then

If the factorization of is known, then the partial fraction decomposition of over can be computed in time

Proof. From (8.5) and the fact that we may compute as the quotient of the Euclidean division of by , we first notice that

Given a root of of multiplicity , let us next prove that

for . Let be the monic minimal polynomial of . We first notice that for all , whence for all and similarly for . Let and be the power series expansions of and at with . It follows that and are bounded by for all , and . Now consider the power series

in . By what precedes, we have and we notice that for . Since and , we conclude that .

The quotient of the Euclidean division of and can be computed in time . The naive expansion of the series and until order takes time and the division of the results with precision can be performed in time . Doing this for one in each conjugacy class yields the announced complexity bound.

Partial fraction decomposition becomes much faster when carried out directly over . In the SLP model, we have the following complexity bound:

Proposition 8.2. Let be as above and assume that the factorization of into irreducibles is known. Given a joint SLP for the evaluation of , , and at points , there exists an SLP for the evaluation of the partial fraction decomposition of of length at most .

Proof. Direct consequence of [65] and the quasi-linear time algorithm for the “untangling” morphism from [52, Sections 4.3 and 4.4].

Since the local reductions really operate on partial fraction decompositions, it is also interesting to study the inverse operation that recovers a rational function from its partial fraction decomposition over .

Proposition 8.3. Let the partial fraction decomposition of over be given and define . Then

and the standard representation can be computed in time

Proof. Let be one the poles in the partial fraction decomposition of and let be its monic minimal polynomial. Denote , , and . Consider the polar part in the partial fraction decomposition of with and . This components gives rise to a corresponding component

in the partial fraction decomposition of over . Now , whence and . The traces of can be computed using Newton–Girard's formulas and for all . Since is linear, it follows that

Summing over all poles, the bound for follows. The time complexity bound is proved in a similar way.

Proposition 8.4. Let and assume that we are given an SLP of length for the evaluation of the partial fraction decomposition of over at points (this includes the evaluation of the roots of the denominator of ). Then there exists an SLP for the evaluation of the numerator and the denominator of of length at most .

Proof. We first convert the partial fraction decomposition over back into a partial fraction decomposition over ; using the tangling algorithm from [52, Section 4.5], this can be done in quasi-linear time. We next recover the numerator and denominator using traditional rational function arithmetic, again in quasi-linear time.

9.Complexity of creative telescoping

9.1.Complexity of creative telescoping

Let be a DD-algebra as in the previous section. For all our complexity bounds, we recall that is assumed to be constant. Consider a telescopable D-finite ideal with associated reduction . Let with entries be a basis for the finite dimensional vector space . In Algorithm 6.1, it often occurs that around iterations are necessary before the main loop terminates. In that case, it may be profitable to use the alternative way to compute the and to further optimize the reduction process by precomputing for all basis elements . For each , this leads to a matrix such that

Given , we then get . We call the reduction matrices for and . With these optimizations, Algorithm 6.1 becomes:

Algorithm 9.1

Input: reduction matrices for and with

Output: D-finite generators of a D-finite ideal included in

, ,

while do

Let be minimal for and set

if is not reducible with respect to then

if then

else decompose and let

if , then

else ,

return

Remark 9.1. Finding sharp bounds for the complexity of the Gröbner basis computation in the final step of Algorithm 6.1 is a problem that is somewhat independent from the cost of creative telescoping itself. For this reason, we removed this step here. The returned set still admits the property that the set of monomials that are not divisible by a leading monomial of an element of is finite; for this reason we call a D-finite set of generators.

Theorem 9.2. Algorithm 9.1 is correct and terminates. Assuming that for and , its running time is bounded by

Proof. Since each is obtained as the product of at most matrices of total degree , we have for all encountered in the algorithm. The number of iterations of the main loop is always bounded by . It follows that the computation of the can be done in time . The linear algebra of testing the existence of a relation (and computing such a relation if it exists) reduces to the computation of the kernel of a matrix in of total degree in and rank , where . This can be done in time through evaluation-interpolation of the numerator at points, using the fact that the determinants of all minors of this numerator have total degree . Multiplying with the maximal number of iterations , the main result follows.

Theorem 9.3. Given an SLP for the joint evaluation of the reduction matrices and , Algorithm 9.1 gives rise to an SLP for the joint evaluation of the polynomials in , whose length is bounded by

Proof. The evaluation of the input SLPs at a point takes operations and yields scalar matrices and a scalar vector . Applying the algorithm to instead of produces instead of . This time, each iteration amounts to one matrix-vector product of cost and adding one row to the matrix of the ; we incrementally put this matrix in echelon normal form in order to determine the relations , again with cost . Since the algorithm uses at most iterations, the complexity bound follows.

Remark 9.4. The cost of the linear algebra in Theorem 9.3 can be further lowered if the matrices are sparse or in the case when we are using a lexicographic admissible ordering. For details, see [38] and references therein.

It may happen that some of the basis elements are actually superfluous for the computation, in which case it is possible to replace the reduction matrices by smaller “quasi-reduction” matrices. More precisely, assume that the vector space spanned by the basis elements is stable under the mappings for . Then each is a block matrix

where . If , then we call quasi-reduction matrices for and . Given with , it is not hard to check that Algorithm 9.1 still applies to instead of , and that the above theorems generalize to this case.

9.2.Deterministic complexity of differential reduction

In this subsection, we assume that and that our D-finite ideal is explicitly telescopable for some cyclic basis . We write for the associated reduction. There also exists an operator of order than annihilates the first entry of . We may write , where and is the monic annihilator of the zero-set of . We assume that the irreducible factors of are all known.

In order to apply Theorem 9.2, we need to bound the dimension of the space of reduced functions and to analyze the cost of computing the reduction matrices . Let us start with the dimension bound.

Lemma 9.5. Under the above assumptions, the reduction can be constructed in such a way that

Proof. We may take in the proof of Proposition 7.5 (see Example 7.6). It follows that . When applying Theorem 4.5 to , we thus have , whence the complexity bound, since .

Given a vector with entries , let us now analyze the cost to compute the reduction . By construction, we recall that , where and denotes the reduction with respect to the adjoint operator . It thus suffices to analyze the cost of the reduction .

Lemma 9.6. Given , we have

and can be computed in time

Proof. Consider one of the poles with , let be its monic minimal polynomial, and . Given , written as a Laurent polynomial, let denote the maximal degree in of one of its coefficients. Let us show by induction on that

This is clear for . Otherwise, consider the expansion for . If , then we get

If and , then the induction hypothesis yields

By construction, we have

Putting all these bounds together, we conclude that .

Let us now bound the quantity . Rewriting in terms of and amounts to a Taylor shift from which it is also possible to read off the indicial polynomial . This yields

Since divides , we also have and . It follows that

For , a similar reasoning at infinity shows that

Let us now turn our attention to a general . Proposition 8.1 implies that . From what precedes, we have at every pole with , whence

From Proposition 8.3, we also have

Combined with the previous bounds, this yields

The time complexity bound is proved in a similar way.

Lemma 9.7. Let

Then there exist quasi-reduction matrices for and with

and we may compute them in time

Proof. Denote . For each , let be minimal such that and for . Since admits at most roots and , we notice that . Consider the -subvector space of and the corresponding -subvector space . The way we defined ensures that we actually have , whence and for . This means that is stable under the mappings . We construct as the quasi-reduction matrices with respect to a basis of .

The basis elements can all be taken to be of the form or , where is the minimal polynomial of some and . For each basis element , it then follows that . Writing for the column vector with entries , we also have and for . Similarly, we have , , and . Applying Lemma 9.6 to , it follows that . Regrouping the coefficients of the partial fraction decomposition of by conjugate roots, we directly obtain the expression of in terms of the basis elements . Using Proposition 8.3, it follows that .

We need to compute for and , which means that we need applications of Lemma 9.6, in view of Lemma 9.5. The complexity bound then follows from Lemma 9.6 and Proposition 8.3.

Theorem 9.8. Assume that and that we are given an explicitly telescopable D-finite ideal for some cyclic basis with entries . Then there exists a polynomial time algorithm to compute a D-finite set of generators of a D-finite ideal that is contained in , as a function of the matrices and the operator associated to .

Proof. In view of the discussion at the end of section 8.3, this is a direct consequence of Lemmas 9.5 and 9.7, together with the straightforward generalization of Theorem 9.2 to the case of quasi-reduction matrices.

Remark 9.9. We did not put a lot of effort in lowering the exponents in the complexity bounds in Lemmas 9.6 and 9.7. One idea that should allow for significant improvements is to avoid conversions between the default representation of rational functions and partial fraction decompositions: it should be possible to carry out most computations directly for the partial fraction representation; see [36] for a similar line of thought in a different context.

9.3.SLP complexity of differential reduction

With the same assumptions as in the previous subsection, let us now carry out the complexity analysis in the SLP model. Before stating the counterpart of Lemma 9.6, let us first show that the local reductions from section 4.2 can essentially be computed in quasi-linear time.

Lemma 9.10. Let with monic minimal polynomial , and . Given a Laurent polynomial , we may compute as in (4.3) using

arithmetic operations in . Similarly, given a Laurent polynomial , we may compute as in (4.4) using

arithmetic operations in .

Proof. The operator can be rewritten as an operator in using a Taylor shift; this can be done in time . Rewritten in this way, the operator naturally operates on Laurent polynomials in . We may also decompose with and, for all ,

This decomposition can be done in linear time since both and are of the form for suitable sets . Now consider the linear operators , on that act on monomials by

Then we observe that the evaluation of (4.3) at can be rewritten as

The second equation is “recursive”, which allows us to compute its solution in quasi-linear time in terms of the size of the equation using relaxed evaluation [48, 47]. This directly implies the first complexity bound. The bound at infinity is proved in a similar way.

Lemma 9.11. Given an SLP for the joint evaluation of , the irreducible factors of and at points , we can compute an SLP for the evaluation of , whose length is bounded by

Proof. Denote and let be an evaluation point. By Proposition 8.2, we can compute the rational fraction decomposition of over in quasi-linear time . For each , let be a root of . Lemma 9.10 allows us to compute in time . Using the fact that

it follows that all for can be computed in time . Using Proposition 8.4, these values allow us to compute in quasi-linear time . Again by Lemma 9.10, we obtain using more operations in .

Lemma 9.12. Let . Given an SLP for the joint evaluation of the coefficients of , the irreducible factors of , and the matrices at points , we can compute an SLP for the evaluation of quasi-reduction matrices for and , whose length is bounded by

Proof. Similarly as in the proof of Lemma 9.7, we apply Lemma 9.11 to for and . We again have , whence each individual can be treated in time . Since there are values of to consider, the conclusion follows.

Theorem 9.13. Assume that and that we are given an explicitly telescopable D-finite ideal for some cyclic basis with entries . Let . Given an SLP for the joint evaluation of the coefficients of , the irreducible factors of , and the matrices at points , we can compute an SLP for the evaluation of a D-finite set of generators of a D-finite ideal that is contained in , whose length is bounded by

Proof. This is a direct consequence of Lemmas 9.5 and 9.12, together with the straightforward generalization of Theorem 9.3 to the case of quasi-reduction matrices.

Remark 9.14. The bound in the proof of Lemma 9.12 is actually quite pessimistic, since it corresponds to highly unlucky values of the roots of the indicial polynomials. In practice, the better bound is likely to hold. In that case, the time complexity per evaluation point in Theorem 9.13 drops to .

Bibliography

[1]

S. A. Abramov. Eg-eliminations. Journal of Difference Equations and Applications, 5(4–5):393–433, 1999.

[2]

S. A. Abramov and M. Petrovšek. Rational normal forms and minimal decompositions of hypergeometric terms. JSC, 33(5):521–543, 2002.

[3]

K. Adjamagbo. Sur l'effectivité du lemme du vecteur cyclique. C. R. Acad. Sci. Paris Sér. I, 306(13):543–546, 1988.

[4]

G. Almkvist and D. Zeilberger. The method of differentiating under the integral sign. Journal of Symbolic Computation, 10(6):571–591, 1990.

[5]

A. Arnold, M. Giesbrecht, and D. S. Roche. Faster sparse multivariate polynomial interpolation of straight-line programs. Technical Report http://arxiv.org/abs/1412.4088, Arxiv, 2015. To appear in JSC.

[6]

M. Aschenbrenner, L. van den Dries, and J. van der Hoeven. Asymptotic Differential Algebra and Model Theory of Transseries. Number 195 in Annals of Mathematics studies. Princeton University Press, 2017.

[7]

M. Ben-Or and P. Tiwari. A deterministic algorithm for sparse multivariate polynomial interpolation. In STOC '88: Proceedings of the twentieth annual ACM symposium on Theory of computing, pages 301–309. New York, NY, USA, 1988. ACM Press.

[8]

S. J. Berkowitz. On computing the determinant in small parallel time using a small number of processors. Inform. Process. Lett., 18:147–140, 1984.

[9]

I. N. Bernshtein. Modules over a ring of differential operators, study of the fundamental solutions of equations with constant coefficients. Functional Anal. Appl., 5(2):89–101, 1971.

[10]

I. N. Bernshtein. The analytic continuation of generalized functions with respect to a parameter. Functional Anal. Appl., 6(4):273–285, 1972.

[11]

D. Bini and V. Y. Pan. Polynomial and matrix computations. Vol. 1. Birkhäuser Boston Inc., Boston, MA, 1994. Fundamental algorithms.

[12]

J. E. Björk. Rings of differential operators. North-Holland, New York, 1979.

[13]

A. Bostan, F. Chen, S. Chyzak, and Z. Li. Complexity of creative telescoping for bivariate rational functions. In Proc. ISSAC '12, pages 203–210. New York, NY, USA, 2010. ACM.

[14]

A. Bostan, S. Chen, F. Chyzak, Z. Li, and G. Xin. Hermite reduction and creative telescoping for hyperexponential functions. In Proc. ISSAC'13, pages 77–84. ACM, 2013.

[15]

A. Bostan, F. Chyzak, and É. de Panafieu. Complexity estimates for two uncoupling algorithms. In Proc. ISSAC 2013, ISSAC '13, pages 85–92. New York, NY, USA, 2013. ACM.

[16]

A. Bostan, F. Chyzak, P. Lairez, and B. Salvy. Generalized Hermite reduction, creative telescoping and definite integration of differentially finite functions. Technical Report, ArXiv, 2018. http://arxiv.org/abs/1805.03445.

[17]

F. Boulier. Étude et implantation de quelques algorithmes en algèbre différentielle. PhD thesis, University of Lille I, 1994.

[18]

R. P. Brent, F. G. Gustavson, and D. Y. Y. Yun. Fast solution of Toeplitz systems of equations and computation of Padé approximants. J. Algorithms, 1(3):259–295, 1980.

[19]

B. Buchberger. Ein Algorithmus zum auffinden der Basiselemente des Restklassenringes nach einem null-dimensionalen Polynomideal. PhD thesis, University of Innsbruck, 1965.

[20]

Bürgisser, P. and Clausen, M. and Shokrollahi, M. A. Algebraic complexity theory, volume 315 of Grundlehren der Mathematischen Wissenschaften. Springer-Verlag, 1997.

[21]

S. Chen, M. van Hoeij, M. Kauers, and C. Koutschan. Reduction-based creative telescoping for Fuchsian D-finite functions. Technical Report, ArXiv, 2016. http://arxiv.org/abs/1611.07421.

[22]

S. Chen, H. Huang, M. Kauers, and Z. Li. A modified Abramov-Petkovsˇek reduction and creative telescoping for hypergeometric terms. In Proc. ISSAC'15, pages 117–124. ACM, 2015.

[23]

S. Chen, M. Kauers, and C. Koutschan. Reduction-based creative telescoping for algebraic functions. In Proc. ISSAC '16, pages 175–182. New York, NY, USA, 2016. ACM.

[24]

S. Chen, M. Kauers, and M. Singer. Telescopers for rational and algebraic functions via residues. In Proc. ISSAC'12, pages 130–137. ACM, 2012.

[25]

R. C. Churchill and J. J. Kovacic. Differential Algebra and Related Topics, chapter Cyclic vectors, pages 191–218. World Scientific, 2002.

[26]

F. Chyzak. Fonctions holonomes en calcul formel. PhD thesis, École polytechnique, France, 1998.

[27]

F. Chyzak. An extension of Zeilberger's fast algorithm to general holonomic functions. Discrete Math., 217(1–3):115–134, 2000.

[28]

F. Chyzak. The ABC of Creative Telescoping — Algorithms, Bounds, Complexity. Habilitation, École polytechnique, 2014.

[29]

F. Chyzak and B. Salvy. Non–commutative elimination in Ore algebras proves multivariate identities. JSC, 26(2):187–227, 1998.

[30]

J. H. Davenport. The Risch differential equation problem. SIAM J. Comput., 15(4):903–918, 1986.

[31]

J. Della Dora, C. Dicrescenzo, and D. Duval. A new method for computing in algebraic number fields. In G. Goos and J. Hartmanis, editors, Eurocal'85 (2), volume 174 of Lect. Notes in Comp. Science, pages 321–326. Springer, 1985.

[32]

J. Denef and L. Lipshitz. Power series solutions of algebraic differential equations. Math. Ann., 267:213–238, 1984.

[33]

J. Denef and L. Lipshitz. Decision problems for differential equations. The Journ. of Symb. Logic, 54(3):941–950, 1989.

[34]

L. Dumont. Efficient algorithms for the symbolic computation of some contour integrals depending on one parameter. PhD thesis, École Polytechnique, 2016.

[35]

M. C. Fasenmyer. Some generalized hypergeometric polynomials. PhD thesis, Univ. of Michigan, 1945.

[36]

R. J. Fateman. Rational function computing with poles and residues. http://www.cs.berkeley.edu/~fateman/papers/poles.pdf, 2010.

[37]

J. C. Faugère, P. Gianni, D. Lazard, and T. Mora. Efficient computation of zero-dimensional Gröbner bases by change of ordering. JSC, 16(4):329–344, 1993.

[38]

J.-C. Faugère, P. Gaudry, L. Huot, and G. Renault. Sub-cubic Change of Ordering for Gröbner Basis: A Probabilistic Approach. In Proc. ISSAC 2014, pages 170–177. Kobe, Japon, jul 2014. ACM.

[39]

L. Fuchs. Die Periodicitaätsmoduln der hyperelliptischen Integrale als Functionen eines Parameters aufgefasst. J. Reine Angew. Math., 71:91–127, 1870.

[40]

A. Galligo. Some algorithmic questions on ideals of differential operators. In Proc. EUROCAL '85, volume 204 of Lecture Notes in Computer Science, pages 413–421. Springer-Verlag, 1985.

[41]

K. Geddes, H. Le, and Z. Li. Differential rational normal forms and a reduction algorithm for hyperexponential functions. In Proc. ISSAC'04, pages 183–190. ACM, 2004.

[42]

R. W. Gosper, Jr. Decision procedure for indefinite hypergeometric summation. Proc. Nat. Acad. Sci. U.S.A., 75(1):40–42, 1978.

[43]

G.-M. Greuel, V. Levandovskyy, O. Motsak, and H. Schönemann. Plural. a singular 3-1 subsystem for computations with non-commutative polynomial algebras. http://www.singular.uni-kl.de, 2010.

[44]

C. Hermite. Sur l'intégration des fractions rationnelles. Ann. Sci. École Norm. Sup. Série 2, 1:215–218, 1972.

[45]

J. van der Hoeven. Differential and mixed differential-difference equations from the effective viewpoint. Technical Report LIX/RR/96/11, LIX, École polytechnique, France, 1996.

[46]

J. van der Hoeven. A new zero-test for formal power series. In Teo Mora, editor, Proc. ISSAC '02, pages 117–122. Lille, France, July 2002.

[47]

J. van der Hoeven. Relax, but don't be too lazy. JSC, 34:479–542, 2002.

[48]

J. van der Hoeven. New algorithms for relaxed multiplication. JSC, 42(8):792–802, 2007.

[49]

J. van der Hoeven. Computing with D-algebraic power series. Technical Report, HAL, 2014. http://hal.archives-ouvertes.fr/hal-00979367, accepted for publication in AAECC.

[50]

J. van der Hoeven. Constructing reductions for creative telescoping. Technical Report, HAL, 2017. http://hal.archives-ouvertes.fr/hal-01435877.

[51]

J. van der Hoeven and G. Lecerf. Sparse polynomial interpolation in practice. ACM Commun. Comput. Algebra, 48(3/4):187–191, 2015.

[52]

J. van der Hoeven and G. Lecerf. Composition modulo powers of polynomials. In Proc. ISSAC '17, pages 445–452. New York, NY, USA, 2017. ACM.

[53]

J. van der Hoeven and É. Schost. Multi-point evaluation in higher dimensions. AAECC, 24(1):37–52, 2013.

[54]

H. Huang. New bounds for hypergeometric creative telescoping. In Proc. ISSAC'16, pages 279–286. ACM, 2016.

[55]

E. L. Ince. Ordinary differential equations. Longmans, Green and Co., 1926. Reprinted by Dover in 1944 and 1956.

[56]

M. Janet. Sur les systèmes d'équations aux dérivées partielles. PhD thesis, Faculté des sciences de Paris, 1920. Thèses françaises de l'entre-deux-guerres. Gauthiers-Villars.

[57]

M. Javadi and M. Monagan. Parallel sparse polynomial interpolation over finite fields. In Proceedings of PASCO 2010, pages 160–168. ACM Press, 2010.

[58]

E. Kaltofen, Y. N. Lakshman, and J.-M. Wiley. Modular rational sparse multivariate polynomial interpolation. In ISSAC '90: Proceedings of the international symposium on Symbolic and algebraic computation, pages 135–139. New York, NY, USA, 1990. ACM Press.

[59]

E. Kaltofen and B. M. Trager. Computing with polynomials given by black boxes for their evaluations: greatest common divisors, factorization, separation of numerators and denominators. JSC, 9(3):301–320, 1990.

[60]

A. Kandri-Rody and V. Weispfenning. Non-commutative Gröbner bases in algebras of solvable type. JSC, 9:1–26, 1990.

[61]

E. R. Kolchin. Differential algebra and algebraic groups. Academic Press, New York, 1973.

[62]

T. H. Koornwinder. On Zeilberger's algorithm and its q-analogue. In Proceedings of the Seventh Spanish Symposium on Orthogonal Polynomials and Applications (VII SPOA), volume 48, pages 91–111. Granada (1991), 1993.

[63]

C. Koutschan. Advanced Applications of the Holonomic Systems Approach. PhD thesis, RISC-Linz, 2009.

[64]

H. Kredel. Solvable polynomial rings. PhD thesis, Univ. Passau, 1993. Published by Shaker.

[65]

H. T. Kung and D. M. Tong. Fast algorithms for partial fraction decomposition. SIAM J. Comput., 6:582–593, 1977.

[66]

J.-L. Lagrange. Miscellanea taurinensia iii. p. 179.

[67]

P. Lairez. Periods of rational integrals: algorithms and applications. PhD thesis, École polytechnique, Nov 2014.

[68]

P. Lairez. Computing periods of rational integrals. Math. Comp., 85(300):1719–1752, 2016.

[69]

G. Lecerf. On the complexity of the Lickteig-Roy subresultant algorithm. Technical Report, CNRS & École polyechnique, 2017. https://hal.archives-ouvertes.fr/hal-01450869.

[70]

V. Levandovskyy. Non-commutative Computer Algebra for polynomial algebras: Gröbner bases, applications and implementation. PhD thesis, Univ. Kaiserslautern, 2005.

[71]

T. Mora. Groebner bases in non-commutative algebras. In Proc. ISSAC'88, number 358 in LNCS, pages 150–161. 1989.

[72]

O. Ore. Theorie der linearen Differentialgleichungen. J. Reine und Angew. Math., 167:221–234, 1932.

[73]

O. Ore. Theory of non-commutative polynomials. Ann. of Math., 34(3):480–508, 1933.

[74]

M. Ostrogradsky. De l'integration des fractions rationelles. Bull. de la Classe Physico- Mathématique de l'Académie Imperiale des Sciences de Saint-Petersburg, IV:147–168, 1845.

[75]

A. Péladan-Germa. Tests effectifs de nullité dans des extensions d'anneaux différentiels. PhD thesis, Gage, École Polytechnique, Palaiseau, France, 1997.

[76]

É. Picard. Sur les intégrales doubles de fonctions rationnelles dont tous les résidus sont nuls. Bull. Sci. Math. (Série 2), 26:143–152, 1902.

[77]

A. van der Poorten. A proof that Euler missed: Apéry's proof of the irrationality of . Math. Intelligencer, 1(4):195–203, 1979.

[78]

C. G. Raab. Definite Integration in Differential Fields. PhD thesis, Johannes Kepler Universität Linz, 2012.

[79]

Ch. Riquier. Les systèmes d'équations aux dérivées partielles. Gauthier-Villars, 1910.

[80]

J. F. Ritt. Differential algebra. Amer. Math. Soc., New York, 1950.

[81]

A. Rosenfeld. Specializations in differential algebra. Trans. Amer. Math. Soc., 90:394–407, 1959.

[82]

C. Schneider. Simplifying Multiple Sums in Difference Fields. In C. Schneider and J. Bluemlein, editors, Computer Algebra in Quantum Field Theory: Integration, Summation and Special Functions, Texts and Monographs in Symbolic Computation, pages 325–360. Springer, 2013.

[83]

J. Shackell. A differential-equations approach to functional equivalence. In Proc. ISSAC '89, pages 7–10. Portland, Oregon, A.C.M., New York, 1989. ACM Press.

[84]

N. Takayama. Gröbner basis and the problem of contiguous relations. Japan J. Appl. Math., 6:147–160, 1989.

[85]

N. Takayama. Kan: a system for computation in algebraic analysis. http://www.math.kobe-u.ac.jp/KAN/, 1991.

[86]

N. Takayama. An approach to the zero recognition problem by Buchberger algorithm. JSC, 14:265–282, 1992.

[87]

J. M. Thomas. Differential Systems, volume XXI. American Mathematical Society Colloquium Publications, New York, 1937.

[88]

B. M. Trager. Integration of Algebraic Functions. PhD thesis, MIT, 1984.

[89]

P. Verbaeten. Rekursiebetrekkingen voor lineaire hypergeometrische funkties. PhD thesis, K. U. Leuven, 1976.

[90]

D. Zeilberger. A holonomic systems approach to special functions identities. Journal of Comp. and Appl. Math., 32:321–368, 1990.

[91]

D. Zeilberger. The method of creative telescoping. JSC, 11(3):195–204, 1991.