Fuchsian holonomic sequences Joris van der Hoeven CNRS, École polytechnique, Institut Polytechnique de Paris Laboratoire d'informatique de l'École polytechnique (LIX, UMR 7161) 1, rue Honoré d'Estienne d'Orves Bâtiment Alan Turing, CS35003 91120 Palaiseau, France Email: vdhoeven@lix.polytechnique.fr

1.Introduction

1.1.Statement of the problems

Let be a subfield of . A sequence is said to be holonomic over if it satisfies a difference equation

(1.1)

where is a linear difference operator in with . (Note that some authors prefer the terminology D-finite or P-finite instead of holonomic.) Many interesting sequences from combinatorics, the analysis of algorithms, and number theory are holonomic [46, 13]. We say that is -holonomic if for all . We say that the equation (1.1) is non-degenerate if for all . In that case, the sequence is entirely determined by its first coefficients and is -holonomic if and only if .

Throughout this paper, we assume that is the field of algebraic numbers. Given a non-degenerate -holonomic sequence , one may raise several natural questions:

Q1

How to compute the asymptotic expansion of when ?

Q2

What kind of constants coefficients can occur in the asymptotic expansion of ?

Q3

How to compute terms of the sequence efficiently as a function of ?

Q4

How to decide whether or for large, all, or infinitely many ? (For this question, we assume that for all .)

These questions are related and can be further refined. For instance, it is natural to compute terms as elements of . However, if becomes large, then the bit-size of is generally at least proportional to . When that happens, it may be preferable to switch to a floating point representation. If we have an asymptotic expansion of with suitable error bounds, then we may exploit that to quickly compute floating point approximations of the for large . Similarly, if the dominant term of the expansion of is positive and provably dominates the other terms, then we may deduce that for large .

Assuming that the generating function

is analytic at the origin, it is well known that the asymptotic behavior of the sequence is closely related to the behavior of near its dominant singularities (i.e. the singularities of smallest absolute value). As will be recalled in section 2.1, the generating function is again holonomic in the sense that it satisfies a non-trivial linear differential equation with polynomial coefficients. Holonomic functions can be evaluated extremely efficiently and their singularities are well understood.

In this paper we will restrict our attention to the special case when at least the dominant singularities of are Fuchsian (see section 2.2 for a definition). In that case, the behavior of near its dominant singularities becomes much simpler and the evaluation of near these singularities even more efficient.

In their full generality, the questions Q1, Q3, and Q4 turn out to be very difficult, even in the Fuchsian case. Indeed, if is actually a rational function, then the last question Q4 is related to Skolem's problem, which asks whether for some . Now if is a rational function, then so is and for all . The hard cases for Skolem's problem occur when has several dominant singularities. One archetype example is

(1.2)

with . A variant of this problem is also relevant for the question Q3: if certain terms of this sequence (1.2) can become “absurdly small”, then the computation of floating point approximations for these terms may take much longer than expected, since we need to compute with a precision that exceeds the order of cancellation. We refer to [32] for more information on Skolem's problem and to [40] for some recent related progress in the context of hypergeometric sequences.

On the positive side, examples like (1.2) are fairly pathological, so it remains reasonable to hope answering our questions for most practical examples from combinatorics or the analysis of algorithms. One interesting concrete example was studied in [35]. The authors exploit the fact that the then open problem about the uniqueness of the Canham model for biomembranes reduces to proving the positivity of a certain holonomic sequence. They solve the latter problem using singularity analysis [13] and techniques for reliable numerical computations with holonomic functions [9, 20, 21, 37, 36]. In the present paper, among other things, we will extend this approach to more general holonomic sequences. Independently of this work, the implementations from [35] were further extended in [10]. Note that other approaches to automatically prove the positivity of sequences were proposed in [17, 31].

In fact, the main purpose of this paper is to provide the best possible answers to questions Q1–Q4 as long as we do not run into number theoretic trouble. We will also identify the precise nature of possible trouble, thereby clarifying which problems need to be overcome if we want to give even better answers. Most of our results rely on well known techniques. Our main contribution is therefore a detailed analysis of how to answer the questions Q1–Q4 as well as possible using these techniques.

1.2.Overview of our contributions

Let us briefly outline the structure of this paper. Section 2 contains reminders about holonomic functions, Fuchsian singularities, and holonomic constants.

In section 3, we start with questions Q1 and Q2. In order to obtain asymptotic expansions, we use Cauchy's classical contour integral formula for and deform the contour into a finite number of loop integrals around the smallest singularities of (section 3.1). Each of these loop integrals is a truncated Mellin-style integral, whose asymptotic expansion can be computed using classical formulas (section 3.2). (In the context of difference equations, note that some authors use the terminology “Pincherle transform” [42] instead of “Mellin transform”.) Adding up the contributions from each of the singularities, we obtain an asymptotic expansion for (Theorem 3.2). The coefficients of this asymptotic expansion can be computed explicitly and expressed in terms of (non-singular) holonomic constants and values of higher derivatives of at points in . Using reliable numeric techniques from [21, 35], one may also compute explicit bounds for the error of the asymptotic expansion (section 3.3).

Unfortunately, Theorem 3.2 is imperfect in the sense that some of the terms in the “asymptotic expansion” for may cancel out (in which case the bound for the error becomes larger than the expansion itself). This may actually happen in three different ways that will be analyzed in detail in section 4.

First of all, consider a holonomic function that converges on the closed unit disk . Its value at is a holonomic constant (and any non-singular holonomic constant can be obtained in this way). Now is also a holonomic function and the first term of the asymptotic expansion of is given by if and only if . This shows that we need a zero-test for holonomic constants in order to detect cancellations of terms in the asymptotic expansion of .

A second kind of cancellation may occur between distinct terms in the asymptotic expansion and only for certain values of . Consider for example

Then the dominant terms and cancel out for odd values of . We call this phenomenon “resonance” and the good news is that it can be entirely eliminated by considering a finite number of cases in which is replaced by a subsequence with and (see section 4.2).

A non-periodic variant of resonance is “quasi-resonance”. Consider for instance the sequence from (1.2). We say that this sequence is quasi-resonant if, for every and , there exist an infinite number of with . In fact, we conjecture that quasi-resonance implies resonance (Conjecture 4.5), but a proof seems currently out of reach.

Assuming a zero-test for holonomic constants and the absence of quasi-resonance, it is possible to automatically compute the asymptotic expansion of in a strong sense, without suffering from uncontrolled cancellations (Theorem 4.6). Under these hypotheses, we also show in section 5 that -bit floating point approximations for can be computed in smoothly linear time : see Theorem 5.2. This bound is uniform as long as .

In section 6 we turn to question Q4. Whenever we can compute an asymptotic expansion of for which the error is strictly smaller than for , the positivity of can be deduced from the positivity of for and the positivity of . However, for a sequence like (1.2) and , it can be hard decide whether for all : what we need here is an even more precise version of Conjecture 4.5. Nevertheless, this example is fairly pathological. For any and , we always have for all sufficiently large . Conversely, if , , , and are -linearly independent, then for infinitely many . In section 6, we will show that something similar holds in general, by relying on sequence counterparts of results from [19].

For our partial answers to questions Q1–Q4, we only need the dominant singularities of to be Fuchsian, except in the case of cancellations that also require the examination of subdominant singularities. In the last section 7, we finally mention a few interesting results that hold if is globally Fuchsian. From bases of local solutions of the differential equation for , it is then classical that we may construct a basis of solutions to the difference equation (1.1) using Mellin transforms based at the corresponding singularities [42, 38, 15]. We show that this theory can be made fully effective and also develop a difference counterpart for the concept of transition matrices from [21]. The Mellin transforms can still be applied in the case when is replaced with a complex variable such that is sufficiently large. This can be used for the construction and efficient evaluation of meromorphic solutions to the difference equation .

One obvious limitation of the present paper is that we only consider the case when the dominant singularities of are Fuchsian. The irregular case has been studied extensively from a theoretical point of view [39, 16, 11, 2, 28, 29, 30]. In a forthcoming work, we intend to investigate this case from a similar perspective as in this paper.

A minor restriction of this paper is that we assumed to be the field of algebraic numbers. This enables us to prove a softly optimal uniform complexity bound for the computation of a -bit floating point approximation of . For more general subfields of , the results in this paper still go through, using ideas from [9], but the uniform complexity bounds in section 5.3 have to be replaced by .

2.Preliminaries

2.1.Holonomic functions

A function is said to be holonomic over if it satisfies a differential equation

(2.1)

where is a linear differential operator in with . Without loss of generality, we may assume that we normalized so that . We say that is -holonomic if at a certain non-singular point .

It is not hard to see that a sequence is holonomic over if and only if its generating function is holonomic over . Indeed, using the dictionary

and modulo normalization, any non-zero operator can be rewritten as a non-zero operator and vice versa. Then (1.1) holds if and only if (2.1) holds.

Example 2.1. The harmonic numbers satisfy the difference equation

for all and the equation

for all . According to the above dictionary, this yields the equation

which can be rewritten as

Taking , we thus get .

Remark 2.2. In the above example, we multiplied the equation by to make it hold for all instead of all . In general, given a difference equation (1.1) that holds for all , we can transform it into a difference equation that holds for all through multiplication by .

From an analytic perspective, if is holomorphic at the origin, then we may retrieve the coefficients using the Cauchy integral

(2.2)

where we integrate over a circle with center and a sufficiently small radius.

2.2.Fuchsian singularities

Consider a holonomic function that satisfies (2.1). The only possible singularities of are located at the roots of or at infinity. Modulo a change of variables with or (for a singularity at infinity), the study of near such singularities can be reduced to the case of a singularity at the origin. If and if (modulo multiplication by a power of ) we can rewrite the equation (2.1) as an equation

with and , then we say that is regular-singular or Fuchsian at . If this is the case at all singularities (modulo the above changes of variables), then we say that is Fuchsian . Sometimes, we will also apply this terminology to solutions of the equation or to the sequence instead of . If is holomorphic at the origin and the non-zero singularities of smallest absolute value of are all Fuchsian

)

If is Fuchsian at the origin, then it is well known [14] that (2.1) has a fundamental system of local solutions of the form

where , , , and denotes the ring of convergent power series in . Moreover, there exists a unique such system of solutions (up to a permutation of indices) with the property that the coefficient of in vanishes whenever . We call it the canonical system of solutions at . We call with a fundamental monomial for .

If is Fuchsian at a point , we thus have a corresponding canonical system of solutions with

If is Fuchsian at infinity, then we also have a corresponding canonical system of solutions with

Let be the row vector with entries . Given a local solution to at , there exists a unique column vector such that . We call the initial condition or generalized value of at . This definition still makes sense at non-singular points, in which case is simply the column vector with entries .

2.3.Holonomic constants

For a precise answer to question Q2, it is important to first introduce various relevant classes of “holonomic constants” that can be obtained as values of holonomic functions. We will follow [26, section 4.4 and appendix B], while restricting us to non-singular and regular-singular holonomic constants.

Let and denote the sets of monic whose coefficients are respectively defined on and . Let be the set of such that is at worst regular-singular at . We also define to be the set of monic operators whose coefficients are defined on and such that is at worst regular-singular at .

We define , , and to be the sets of solutions to an equation , where , , or , respectively, and such that exists. Then we define

We call the set of holonomic constants and the set of Fuchsian holonomic constants. Each of these three sets actually forms a ring [26, Proposition B.1]. Moreover, these three rings turn out to be closely related (see also [12]):

Theorem 2.3. [26, Theorem B.5] We have

where

Now consider a holonomic sequence whose generating function belongs to . So for some . Assume that has a Fuchsian singularity at and let be the canonical system of local solutions at . In particular, there exist unique with . In fact, we have (see [26, Proposition B.3]) and we can compute using the algorithms from [21]. Here “computing” is understood in the following sense: for any and , we can compute an approximation of with . Note that this does not imply the existence of a zero-test for the constants . The zero-test problem for holonomic constants will be discussed in section 4.4 below.

3.Asymptotic expansions

3.1.Decomposing Cauchy contour integrals into Mellin integrals

The traditional method to determine the asymptotics of a sequence whose generating function is holomorphic at the origin is based on the Cauchy contour integral (2.2):

(3.1)

If is holonomic, then has only a finite number of singularities , which are all in . For some and , assume that for and for . Then we may deform the contour from (3.1) into a contour that consists of truncated Hankel contours and residual arcs on the circle with center and radius : see Figure 3.1. Then (3.1) becomes:

(3.2)

We may chose the truncated Hankel contours to depart radially from the origin toward infinity. In the degenerate case when the arguments of certain singularities coincide, we turn the contours clockwise by a fixed sufficiently small angle: see Figure 3.1 .

Figure 3.1. Deformation of a Cauchy contour into a contour of radius that avoids a finite number of singularities. Since and are aligned with the origin, we slightly modified the directions of the corresponding truncated Hankel contours and to avoid collisions.

Integrals of the form

(3.3)

are called truncated Mellin integrals

,

where . Note that an upper bound for can be computed efficiently using the algorithms from [20].

Most of the remainder of this section is dedicated to determining the asymptotics of integrals of the form (3.3) in the case when is Fuchsian at . The integrals for which is smallest typically dominate the other ones, but cancellations may sometimes occur, in which case we will also need to examine the subdominant singularities. Let us consider one simple example of this phenomenon.

Example 3.1. The rational function

from the introduction is certainly holonomic, with singularities at , , and . We have

where we note that for odd values of . In other words, the asymptotics of depends on the parity of :

			(n∈2*ℕ)
			(n∈2*ℕ+1)

3.2.Elementary Mellin integrals

Let us first study the very special case when

with , , and . Explicit formulas for the asymptotics of the Taylor coefficients are well known in this case, but it is convenient to recall the details of this computation. Modulo a change of variables , we may assume assume without loss of generality that .

In this special case, for , the contour integral (3.1) can rewritten into the full Mellin integral

where is a Hankel contour that starts at , then turns clockwise around , and finally returns to . Setting , this formula becomes

where is a Hankel contour that starts at , then turns clockwise around , and finally returns to . We regard as an element of , where denotes the set of polynomials in of degree in :

with (Kronecker delta) and for all . Then

(3.4)

Let

From the classical formula [45, Section 12.22]

we deduce

whence

			(3.5)

When plugging this identity into (3.4), we obtain an asymptotic expansion

where can be computed explicitly. Note that the series that underlies this expansion usually diverges. It will be useful to introduce

3.3.Effective local error bounds

Assume now that has a Fuchsian singularity at and consider a local solution of the form

where , . We recall that general local solutions to are linear combinations of local solutions of this special form. For some sufficiently small, consider

where is a Hankel contour from around and then back to . Our aim is to determine both the asymptotics of and an effective bound for the remainder.

As in section 3.2, we first reduce to the case when and then perform the change of variables . Let and , so that

Let be a truncation order and let

be the truncation of modulo , with . Using (3.5), we can explicitly compute

as an element of .

Given , let us now study the difference

We have

Moreover, assuming that , we have , , and for , so

			(i⩽Re κ)
			(i⩾Re κ)

Since is a linear combination of functions with , this allows us to compute an explicit bound

for a suitable constant .

Since is Fuchsian at , by taking sufficiently small, we can use the algorithms from [21] to compute a bound

for all with and . Given , we may further transform this into a bound

If , then

			(3.6)

Altogether, we may compute constants and polynomials such that

for all . Such a bound can also be proved without the assumption , by doing partial integrations before estimating the left-hand side of (3.6). We omit the details, since the case when suffices for the proof of our main Theorem 3.2 below.

To conclude this section, let us finally consider a solution to with for all . Let be the fundamental system of local solutions to with

for . Using the algorithms from [21], we may compute the unique constants such that . Let and let be minimal integers with for . We have shown above how to compute polynomials and constants such that

for all . Setting and , it follows that

for all .

3.4.Asymptotic expansions with effective error bounds

Let us now return to the general setting from subsection 3.1: for some and , we assume that for and for . We also assume that each of the singularities with is Fuchsian.

In the previous subsection, we have seen how to compute asymptotic expansions with effective error bounds for the truncated Hankel integrals

for . Assuming that , the truncated Hankel contour consists of and two straight stretches between and . Using the algorithms from [20], we may compute a bound for on these two stretches. Then we have

Combining this with the bounds for the residual integrals, this gives a first answer to questions Q1 and Q2:

Theorem 3.2. Let , , , , , and be as in the text above. For any , we can compute an asymptotic expansion for in together with bounds and such that

(3.7)

for all .

Remark 3.3. In principle, the error at the right hand side of (3.7) can be replaced by a single term of the form , where is smallest among . However, from a numerical point of view, some of the and may be very small (or even zero), in which case it may be preferable to use the error bound from the theorem. We will return to this issue in section 4.4 below.

4.Periodic or quasi-periodic cancellations

4.1.Classical results for rational functions

In Example 3.1, we have seen that the contributions of more than one dominant singularity may cancel out in a periodic fashion. Is it possible to predict when this phenomenon occurs? As demonstrated by Example 3.1, this is already an interesting question in the case when is a rational function. In that case, exact cancellations can actually be predicted, as we will recall now. In section 4.2, we will deal more generally with approximate cancellations as in Example 3.1.

So consider a sequence whose generating function is is rational. Then the classical Skolem-Mahler-Lech theorem tells us that the zero-set is ultimately periodic. This was first proved by Skolem in the case when , next by Mahler for , and finally by Lech for general fields of characteristic zero:

Theorem 4.1. [44, 34, 33] Let be a field of characteristic zero and let be a sequence whose generating function is rational. Then there exists a period and finite sets and such that

(4.1)

The periodic part in the above decomposition is actually computable [4], whereas the computability of the exceptional part is currently an open problem [41]. Let be the singularities of and let . We say that is resonant if for some . Setting

(4.2)

this is the case if and only if . Berstel and Mignotte proved the following:

It is instructive to detail the computability statements. Let and . Consider for which with . Since is a subfield of , we must have . Using that for [43, Theorem 15], it follows that . This allows us to compute as the smallest with . We conclude that is computable.

Now for every , the generating function of satisfies

(4.3)

where . This allows us in particular to test whether and to compute the set . If needed, we may also deduce the minimal period and corresponding for which (4.1) holds.

4.2.Removing resonance

Let us now consider an arbitrary holonomic sequence whose generating function is convergent at the origin. Let be the non-zero singularities of . We define as in (4.2) and say that is resonant if . Note that the arguments at the end of the previous subsection still allow us to compute . As in the algebraic case, the coefficients of a resonant holonomic sequence can vanish in a periodic manner, or become disproportionally small on a regular basis. This problem can be removed through the consideration of subsequences and case separation, as in Example 3.1:

Proof. The generating function of satisfies (4.3). Using standard closure properties, it follows that is holonomic. The singularities of the right-hand side of (4.3) are contained in the set . Therefore the singularities of are contained in the set . So it suffices to show that whenever . Assume for contradiction that , but for some . Without loss of generality, we may assume that is minimal with this property. But then we have for some with , which implies by (4.2): a contradiction.

The above proposition has an operator counterpart. Let and consider a monic operator of order with . Let now be the singularities of (i.e. the zeros of the dominators of its coefficients). We define as in (4.2) and say that is resonant if .

Let us show how to compute annihilators for the generating functions . First of all, for each , we observe that for the monic operator . It follows that is an annihilator for . For each , we next observe that , so is still a monic annihilator of . Furthermore, since , we must have . Setting and , we finally note that . This allows us to rewrite as a monic operator in with for all .

We note that the sets of non-zero singularities of and (i.e. the zeros of their denominators) both coincide with the set from the proof of Proposition 4.3. Consequently, the set of non-zero singularities of is . In other words, is non-resonant.

Remark 4.4. We may regard the transformation that maps to as a differential counterpart of the Graeffe transform [27]. Indeed, if is a rational function, then the denominator of must be the -fold Graeffe transform of for each , up to constant multiples.

4.3.Quasi-resonance

Even for a non-resonant holonomic sequence, the contributions of the dominant singularities to its asymptotics may occasionally cancel out. For instance, theoretically speaking, the sequence from (1.2) might occasionally become small, although this would necessarily happen in a non-periodic manner.

Let us make this unlikely phenomenon more precise in the case when is dominant-Fuchsian. Let be the dominant singularities for this sequence (we assume that ). We say that is quasi-resonant if for any constants and , there exist infinitely many with

(4.4)

In fact, we conjecture that this never happens:

Although Conjecture 4.5 is far beyond the current state of knowledge in number theory, let us provide meager evidence why we believe that it might hold. The conjecture is already interesting in the case when is a rational function. Even more specifically, assume that

where are such that , but . If , then we clearly have for all . Assume that and that , , and are -linearly independent. Then Baker's theorem [1] implies the existence of a (computable) constant such that

for all but a finite number of . Taking exponentials, it follows that

and

for all but a finite number of . If with , then Baker's theorem implies the existence of a constant with

for all but a finite number of . Using a similar reasoning as above, we deduce that our conjecture again holds in this particular case.

Of course, our general conjecture is far more ambitious and we have no convincing further evidence for it yet. It may be interesting to consider the slightly weaker conjecture for which (4.4) is replaced by

(4.5)

for some fixed constant .

4.4.Asymptotic expansions

Using Proposition 4.3, we can generalize the technique from Example 3.1 and reduce the determination of the asymptotic expansion of a holonomic sequence to the non-resonant case, modulo a finite number of case separations.

In order to detect periodic cancellations we still need a way to decide whether a particular solution to the equation in section 3.1 is actually analytic at a given regular singularity of . Now we may write as a -linear combination of the canonical basis of local solutions at . Let be the subset of indices for which is singular at . Then is analytic at if and only if for every .

In other words, if we have a way to decide whether a given Fuchsian holonomic constant is zero, then can determine the actual dominant singularity of . By Theorem 2.3, it actually suffices to be able to decide whether a holonomic constant in is zero. We will denote by Hol an oracle to do this. Here we assume an exact representation for elements as as the value at of a holonomic function in that is given via a vanishing operator in and a finite number of initial conditions in .

It is interesting to point out that if we can determine the dominant singularity of a holonomic function that is analytic at the origin, then we also have an algorithm to test whether holonomic constants in are zero. Indeed, given , let be such that . Then is the dominant singularity of the holonomic function if and only if .

The above considerations lead to the following variant of Theorem 3.2:

Theorem 4.6. Let , , , , , and be as in section 2.3. Assume that we have an oracle Hol, that is non-resonant, and that is singular at for some . Modulo re-ordering indices, assume that , , and for . Then, for any , we can compute an asymptotic expansion for in together with and such that

(4.6)

for all . Moreover, if Conjecture 4.5 holds, then we may require in addition that

for all sufficiently large .

5.Uniformly fast evaluation

In this section, we investigate question Q3. We assume that satisfies the difference equation (1.1). We also assume that its generating function is convergent at the origin and that it satisfies a holonomic equation (2.1) that is Fuchsian at the origin. Our goal is to compute for large using a precision of bits, with a good uniform complexity in both and .

5.1.Preliminaries

Before we proceed, let us briefly recall some notations and basic facts about fast arithmetic. We define to be the time that is needed to multiply two -bit integers and we make that customary assumption that is non-decreasing. It was shown recently in [18] that we may take .

Approximate computations with real and complex numbers can be done using either fixed point or floating point arithmetic. Let be the set of dyadic numbers. A -bit fixed point approximation of a complex number is a number with . A -bit floating point approximation of is a number with , , , , and either or .

It is well known [8] that -bit approximations of bounded exponentials and logarithms can be computed in time , both for fixed point and floating point representations. Here a bounded exponential is a number with and for some fixed constant . Similarly, a bounded logarithm is a number with . For non-zero , this allows us to convert between -bit floating point approximations of and -bit fixed point approximations of .

We also observe that, given with , we can compute a -bit approximation of in time . Indeed, it suffices to compute and with and then use the fact that .

We will finally need the following result that was essentially proved in [21].

Proof. We approximate using the technique from [21, section 3]. The matrices from (3.10) and (3.12) in that section have size . Since , it suffices to evaluate for in order to obtain -bit approximations for the values . Using binary splitting, this requires bit operations. Since -bit approximations of can be computed using bit operations for all , this allows us to compute a -bit fixed point approximation of in time .

5.2.Exact computation of the terms of a holonomic sequence

If , then the most efficient strategy is to compute exactly as an element of and convert the result into a fixed point or floating point approximation. In fact, as in [20, 21], it suffices to compute in the algebraic number field generated by the coefficients of and the initial conditions with . Let us recall from [9, 21] how to compute using binary splitting.

For each , let be the column vector with entries . If , then

where

(5.1)

More generally, for and , we have

and can be computed efficiently using binary splitting:

Since is Fuchsian at the origin, the bit-size of the entries of is bounded by . Using a classical complexity analysis [9, 21], it follows that can be computed in time . In particular, we may compute in time . Converting the result into -bit fixed point or floating point notation can be done in time .

5.3.Fast computation of floating point approximations

Having dealt with the case when in the previous subsection, let us now assume that and . If has a dominant singularity at , then typically grows like , in first approximation. Consequently, a -bit fixed point approximation of typically requires bits if . In particular, unless , then it is hopeless to compute such approximations in smoothly linear time in . From now on, we will focus on the computation of a -bit floating point approximation of . Under the assumption that is not quasi-resonant, we will show that such an approximation can be computed in uniform smoothly linear time.

So assume that is not quasi-resonant and let be the singularities of . Assume that are the dominant singularities of with . We compute using

where is an axial truncated Hankel contour around until and is a circular arc around from to (or to if ), for . Let be a temporarily increased working precision with and . We will specify later. For some sufficiently small with and any , we take

Since is Fuchsian at , we may compute positive constants , , and such that

for all (uniformly, for all ). It follows that

For , we have

where

We note that the integrand satisfies a holonomic equation whose total size is bounded by . Consequently, the same holds for . Since is analytic at and bounded by for , the uniformity assumptions of Proposition 5.1 are satisfied. It follows that we may compute a -bit fixed point approximation of using bit operations. We may also compute a -bit fixed point approximation of using bit operations. Altogether, this allows us to compute a number with

using bit operations. By construction, it follows that

provided that . Since is not quasi-resonant, there also exist constants , , and such that

for all . In order to obtain -bit floating point approximations of and then , it suffices to chose in such a way that

Using that , this is certainly the case if we take

Note that we have indeed have for this choice of . In fact, we even have as soon as .

In combination with the results from section 5.2, we have proved the following:

Remark 5.4. Note that the operator with may have singularities with , as long as the particular solution remains analytic at . Assuming that is Fuchsian and that we have an oracle Hol, we may verify whether this is the case by checking that the coefficients of the singular canonical solutions in all vanish. Note that our theorem and its corollary only claim the existence of an efficient algorithm to compute ; for this, we do not need the oracle Hol.

Remark 5.5. In practice, for the fast evaluation of a Fuchsian holonomic sequence , we first make it non-resonant using the pre-treatment from section 4.2, then apply the algorithm from the previous subsection for the evaluation of , while falling back on the slower algorithm from section 5.2 whenever we detect massive cancellation. In particular, this mixed strategy takes care of exceptional values of for which vanishes. Whenever this algorithm does not run in time when , we note that would actually provide an explicit counterexample to Conjecture 4.5.

Remark 5.6. If we replace (4.4) by (4.5) for some fixed constant , then the theorem and its corollary still hold, but the complexity bound becomes .

6.Positivity testing

In this section, we study question Q4. We assume that is a holonomic sequence whose generating function is convergent at the origin. Modulo the pre-treatment from section 4.2, we may assume without loss of generality that is non-resonant. We also assume that is dominant-Fuchsian. Our aim is to decide whether or for all or for all sufficiently large .

6.1.A density theorem for sequences

Our positivity test will rely on a way to compute limsups and liminfs of certain oscillating sequences. For this, we will adapt results from [6, 19]. Recall that a Hardy field is a field of germs of differentiable real functions at infinity that is closed under differentiation [7]. In particular, and are Hardy fields. The following is a direct consequence of [19, Theorem 3].

Theorem 6.1. Let be real functions whose germs at infinity belong to a Hardy field and such that , i.e. . For each , let be -linearly independent numbers in . Consider the function

from into the torus of dimension . Then is dense in for any .

What we really need is a counterpart of this theorem for sequences:

Proof. For , let

We also define

Let and let be any constant with . If is sufficiently large, then for all . By the theorem, the image is dense in . Given and , we may thus find an such that

where . Let be an integer with minimal distance to . Since , we have in particular , whence . It follows that

whence

This shows that is indeed dense in .

We will also need the following counterpart of [19, Theorem 5]:

Theorem 6.3. Let be real functions with . For each , let . Consider the function

from into the torus of dimension . Let

be a Laurent polynomial that takes only real values on . Then and are both computable numbers in .

Proof. Using the rewriting techniques from the proof of [19, Theorem 5], we first reduce the general case to the case when are -linearly independent for . (Note that there exists an algorithm to find -linear dependencies between algebraic numbers, so we do not need the general oracle to find -linear dependencies between exp-log constants.)

We next reduce to the case when are -linearly independent. Assume on the contrary that with and . Without loss of generality, we may assume that . Taking for , we may rewrite and as Laurent polynomials in , while performing the corresponding corresponding substitutions in . We repeat this procedure with in the role of until are -linearly independent.

After the above reductions, we are in a position to apply Corollary 6.2. This yields

Now and can be computed using classical algorithms from effective real algebraic geometry [3].

6.2.Positivity testing

Let be the dominant singularities of . By Theorem 4.6, we may compute constants , , , , , , and such that

(6.1)

for all and for all . Setting

we may rewrite (6.1) as

(6.2)

Now observe that can be interpreted as a polynomial

Since for all , this allows us to apply Theorem 6.3 and compute and . If , then (6.2) yields for all . If , then for infinitely many . For any with , the relation (6.2) then yields .

The only remaining case is when . In lucky cases, we may look at the next subdominant term of the asymptotic expansion of and prove the positivity of in a similar way as above. However, in general, the positivity of can be hard to decide. In fact, a general decision procedure would allow us to answer difficult questions about diophantine approximability. For instance, given a real algebraic number , the positivity of the holonomic sequence

is related to a rate of diophantine approximability of by of the form

We do not know of a general decision procedure for this kind of inequalities.

7.Fuchsian holonomic sequences

Let be the generating function of a holonomic sequence . So far, we have mainly been interested in the case when is convergent at the origin and dominant-Fuchsian (possibly modulo a reduction to the non-resonant case as in section 4.2). This is indeed sufficient for obtaining information about through its asymptotic properties. In this section, we assume that is globally Fuchsian and study nice additional properties that hold in this case.

For simplicity, we assume that satisfies a non-degenerate difference equation (1.1). We normalize to make it divisible by and we let be the corresponding differential operator with , as constructed in section 2.1. We let be the singularities of .

7.1.Full Mellin integrals

First of all, for Fuchsian , the truncated Mellin integrals from (3.3) tend to a limit

(7.1)

when tends to infinity and is sufficiently large. This leads to the exact representation

(7.2)

We call (7.1) a full Mellin integral that is based at .

Remark 7.1. We may take to be the contours from infinity to a point close to , which next performs a complete circular turn around , and then goes back to infinity along the same direction where it came from. Then we have

(7.3)

where denotes the circle around and denotes the difference between the analytic continuations of that get around on the left and right, respectively. Note that is a solution of the same differential equation as , i.e. . The formula (7.3) is convenient for machine computations due to the fact that we only have a single stretch going to infinity.

More generally, for an analytic function on and with at infinity, we define

Then we note that

where the second relation is proved using integration by parts. Consequently, the sequence satisfies the same recurrence relation (1.1) as .

7.2.Canonical solutions via Mellin transforms

We observe that (7.1) only depends on the behavior of at the singularity . Expressing in terms of the canonical basis of local solutions to at ,

and setting

it follows that

(7.4)

Here we note that whenever is analytic at . Moreover, we recall that

for some and . If or , then the sequence has a non-zero formal transseries in as its asymptotic expansion, by the formulas from subsection 3.2. We will sometimes identify with this transseries.

Let be the collection of all non-zero with and . Since the dominant monomials of are pairwise distinct when considering them as transseries, these sequences are -linearly independent. From (7.2) and (7.4), we also know that any sequence solution of (1.1) can be written as a -linear combination of . Conversely, we noted at the end of subsection 7.1 that each is actually a solution of (1.1). This shows that forms basis of the solution space of (1.1) in . Since we assumed to be non-degenerate, this space has dimension , so . We denote by the row vector with entries .

7.3.Canonical formal solutions at infinity

We have shown that (1.1) has a basis of formal transseries solutions

In fact, it is well known [5] that we may compute a canonical system of formal solutions in , similar to the ones that we saw in the differential case in subsection 2.2. Let us briefly describe how to do this.

We have seen that for any singularity of , there exists at least one formal solution of the form . Modulo a transformation , we may assume without loss of generality that . We next replace by in , where . After multiplication by a suitable power of , this yields an operator

with . Whenever we have a formal solution with , then this implies that is divisible by . Inversely, if is a root of multiplicity of in , then for any , there exists a unique formal solution to (1.1). Indeed, writing with , we have

which yields the recurrence relation

(7.5)

for the computation of the coefficients . The solution is unique when requiring that and that is divisible by whenever is a root of multiplicity of in .

Let be the collection of all formal solutions of (1.1) in that we obtain in the way that we just described, with . Since the dominant monomials of these solutions are pairwise distinct, this again forms a fundamental system of solutions; we call it the canonical system of solutions of (1.1) at infinity and we denote by the row vector with entries . Since and are both fundamental systems of solutions, there exists a matrix with

(7.6)

Since has coefficients in and has coefficients in , the matrix must actually be in .

For machine computations, the recurrence relation (7.5) is not very efficient if we want to compute a large number of terms. In that case it is better to rewrite the original equation (1.1) with respect to , which transforms the shift operator into . Compositions of a power series with can be computed efficiently in a relaxed manner using the algorithm from [22, section 3.4.2]. The equation (1.1) is not necessarily “recursive”, so it is not always possible to directly solve it using the techniques from [22]. Nevertheless, it can always be rewritten as a recursive equation using the algorithms from [24]. Altogether, this allows us to compute the first coefficients of the canonical solutions in time , which is softly optimal in the bit-size of the result.

Remark 7.2. The equation (7.6) can be used to map formal transseries solutions of (1.1) to actual holonomic sequences. It is interesting to note that this association actually preserves all difference ring operations, in a similar way as accelero-summation in the differential setting [23].

7.4.Transition matrices

Let us briefly recall the concept of a transition matrix, which forms an important ingredient for the efficient evaluation of holonomic functions in [9, 20, 21]. We will use the notations and from section 2.2 for canonical systems of local solutions and generalized values at .

Given a non-singular path between two non-singular points , the analytic continuation of the canonical solutions at can be expressed as linear combinations of the canonical solutions at . In other words, there exists a matrix with

We call the transition matrix along the path . We naturally have the relation

for composed paths. In terms of generalized values of a solution to (2.1), we also obtain

These notions extend to the case when the paths start and/or end at singular points, modulo the precaution that we specify the angles that are used to approach the singularities (in order to determine the branch of the logarithm).

7.5.Transition matrices for sequences

Let us now study the analogue of transition matrices for sequences. Given and , assume that (1.1) has a unique solution with and for . Then we will denote this solution by . We denote by the row vector with entries if these solutions are all defined and call it the canonical system of solutions at . Given a general solution to (1.1), we call the column vector with entries the generalized value of at , so that

By definition, the satisfy a recurrence relation

where was defined in (5.1). More generally, for and , we have

Setting , this relation extends to the case when . Dually, we also have

We call the transition matrix between and . We have seen in section 5.2 how to compute in time using binary splitting.

7.6.Transition matrices at infinity

In section 7.2, we introduced the canonical system of solutions of (1.1) at infinity. Given a solution of (1.1), this leads to the corresponding notion of generalized value at infinity with

(7.7)

For any , the matrix with

is called the transition matrix between and , whenever it exists.

Using a combination of the techniques so far, we may compute the transition matrix as follows. Consider one of the canonical solutions of (1.1) at and the corresponding power series solution of (2.1) at the origin. Given one of the singularities of , we may use the algorithms from [21] to compute the transition matrix for and then re-express as a -linear combination of the canonical solutions at . The collection of these relations yields as a -linear combination of the canonical solutions . Using the methods from sections 7.2 and 7.3, we finally obtain as a -linear combination of . Doing this for each with , this yields the transition matrix . For fixed and large , we may compute -approximations of the entries of in time , using the algorithms from [23].

7.7.Analytic solutions to the difference equation

Our main focus in this paper is on the study of sequence solutions to holonomic equations (1.1). Yet, it is interesting to note that the definition of Mellin transforms of the canonical solutions generalizes to complex numbers for which decreases sufficiently fast on :

(7.8)

Applying this to the theory from section 7.2, this yields a fundamental system of analytic solutions to the difference equation

that extends our fundamental system of sequence solutions. For , we also note that the integrand of (7.8) is still holonomic over and Fuchsian. Consequently, we may evaluate it with a precision of bits in time , using the algorithms from [21]. For general , the evaluation can be done in time using the baby-step-giant-step technique from [9]. It is also possible to extend the uniform complexity analysis from section 5 to this setting, but additional care is needed for the treatment of non-real arguments . We intend to carry out the detailed analysis in a forthcoming paper.

Acknowledgments. We are grateful to Marc Mezzarobba and Ruiwen Dong for some references and helpful discussions.

Bibliography