1.Introduction

Statement of the problem and the main result

Let be a subfield of . A holonomic function over is a solution to a linear differential equation , where is a monic linear differential operator of order . Many classical special functions, such as , , , , , hypergeometric functions, Bessel functions, the Airy function, etc. are holonomic. Moreover, the class of holonomic functions is stable under many operations, such as addition, multiplication, differentiation, integration and postcomposition with algebraic functions.

In the sequel, and unless stated otherwise, we will assume that is the field of algebraic numbers. The only singularities of a holonomic function as above can occur at the poles of the rational functions ; let denote the finite set of these poles. We will say that has initial conditions in if for a certain non-singular point . In this paper, we are interested in the design of efficient algorithms for the numeric evaluation of such a function , with a particular focus on high precision and uniform efficiency as a function of the argument .

For a fixed non singular evaluation point, say , an efficient general purpose algorithm was first given by the Chudnovsky brothers [3]. More precisely, in the case when , they proved that an -bit approximation of can be computed in time . Here stands for a complexity bound for integer multiplication and it has recently been proved that one may take , where . The Chudnovsky–Chudnovsky algorithm was rediscovered in [7] and generalized to the case when is the field of algebraic numbers. An early precursor and further variants can be found in [2, 10, 6].

In order to design uniformly efficient evaluation algorithms, it is crucial to control the efficiency when approaches one of the singularities in . Actually, one first question concerns the computation of the limit of a holonomic function at a singularity if this limit exists. This was first done in [8] for so called regular singularities (achieving the same complexity bound as for non singular points), and in [9] for irregular singularities (in which case we showed that -bit approximations of limits can be computed in time ). We refer to [8, 9] for the definitions of the concepts of regular and irregular singularities.

The main aim of this paper is to achieve the same kind of complexity bounds uniformly in . Such bounds need to be stated with a lot of care. First of all, a holonomic function such as grows exponentially fast at infinity: given the -bit number , one needs bits to merely write down the closest integer approximation of . Using floating point approximations for both and does not help, since a similar explosion then occurs for the exponent. But we may hope for a good uniform complexity bound if we use fixed point approximations for and floating point approximations for .

Another complication is due to the number zero, which should be regarded as a singularity when using floating point representations: it is difficult to compute accurate floating point approximations for if is close to a zero of . Predicting the exact locations of zeros of holonomic functions is a notoriously difficult problem. Even the basic question to decide whether for admits no algorithmic answer for the moment. Nevertheless, the number is often the approximation of some other complex number with a precision of binary digits behind the dot. In that case, it is natural to consider the more general evaluation of on the ball with center and radius , and to require that admits no zeros on this ball.

We are almost in a position to state the main result of this paper. Let be the set of dyadic numbers. Given , we denote by the bitsize of . Given , we also denote . The set of floating point numbers is defined in the same way as , but the exponent of a floating point numbers is represented in binary notation, so that the bitsize of is now fiven by .

Let be the point at which we specified the initial conditions of . We define to be the open subset of of all points such that the straightline segment from to does not intersect . We take to be the unique solution of on that matches the prescribed initial conditions at . Let denote the set of zeros of . The main theorem of this paper is the following.

Proof strategy

As long as remains in a compact subset of in Theorem 1, the conclusion essentially follows from the existing complexity bounds in [3, 7]; using a refinement [11] of the complexity analysis from [7], one even obtains the stronger complexity bound for the evaluation of . Using the techniques from [8], these complexity bounds generalize to subsets of , where is a compact set that contains none of the irregular singularities of . If the point at infinity is a regular singularity, then the bound also applies on subsets for sufficiently large , modulo the change of coordinates .

The above discussion shows that the proof of Theorem 1 involves two main difficulties: controlling the complexity near irregular singularities and controlling the complexity of evaluating near zeros of . For the first task, we will adapt the technique of accelero-summation from [9]. For the second task, we rely on the idea that can never simultaneously become “smaller than expected”. A precise statement will be presented in Section 4; this statement can be regarded as a quantitative version of the well-known property that cannot vanish simultaneously unless vanishes itself.

Let us return to the evaluation of near an irregular singulary, say . At the origin, it is well-known that admits a basis of formal solution of the form

for , where , , , , and where for some . In [9], it is shown that the series are accelero-summable and that we can associate actual functions to them that are defined on sectors of the form

Moreover, a finite number of these sectors can be made to cover a punctured neighbourhood of the origin. One crucual step toward the design of an efficient evaluation algorithm for on such a sector is to deal with the special case when for some , which further reduces to the case when .

In the remainder of this paper, we will assume that the reader is familiar with [9] and the notations that we used there. For simplicity, we will also restrict to accelerations and Laplace transforms such that we integrate on the positive real axis. Using a change of variables for a suitable , this entails no loss of generality. More precisely, we assume that we are in the following situation. The function is the result

of an accelero-summation process with critical times , , and all integrals taken on the positive real axis. The accelero-sum is defined in some sector

for any with .

For any fixed , the accelero-summation process from [9] provides us with an algorithm to compute digits of in time . In Section 3.1, we will show that this complexity is uniform in , provided that for some computable constant . In other words: accelero-summation is a good numerical scheme under the condition that we really need a lot of digits. In [9], we also showed that the technique of “summation until the least term” [12] allows to compute digits of in time , provided that for some computable constant . This complexity bound is also uniform in .

The above uniform complexity bounds still leave a gap for precisions between and . In order to fill this gap, we introduce the technique of expedito-summation in Section 2. Roughly speaking, we perform the accelero-summation process until some critical time with and then “expedite” the process by directly taking a truncated Laplace transform with respect to . We will show in Section 3.3 that there exist computable constants such that expedito-summation until the critical time allows us compute digits of in time , uniformly in provided that .

Notational conventions

This paper should be regarded as a supplement to [9]. For this reason, and as we already stressed before, we will freely use concepts and notations from that paper. In this area it also frequently happens that there exist algorithms to explicitly compute various constants involved in error bounds, but that the precise values of these constants are irrelevant. In [9], we strived to make all error bounds as explicit as possible, but in this paper we will simply denote strictly positive constants of this kind by . In analysis, the habit to write for “some bounded function” is somewhat analoguous. For instance, given a real function and a constant , saying that

for all means that we can compute an explicit exponential bound for on the interval .

Acknowledgments. We wish to like Grégoire Lecerf for some helpful remarks.

2.Expedito-summation

Throughout this and the next section, we make the following assumptions:

• with is a formal solution to .

• is accelero-summable with critical times and .

• The holonomic equations satisfied by the Borel counterparts at the various critical times admit no singularities on the positive real axis.

• All acceleration integrals and the final Laplace transforms are performed on the positive real axis.

2.1.Introduction to expedito-summation

In [9], we provided a detailed analysis of two summation methods of . The usual accelero-summation process associates the accelero-sum to using

In the appendix, we also considered “summation up to the least term”: given , one may approximate by , where

Taking , we proved that

for all .

Summation up to the least term completely shortcuts the whole accelero-summation process. It provides approximations of a precision that correspond to stopping the accelero-summation process at the first singularity for the first critical time. It is natural to consider more general shortcuts, where we perform the usual accelero-summation process up till a given critical time and then “expedite” the remainer of the process by directly performing a truncated Laplace transform on for a suitable . More precisely, given and , we define

Here denotes the contour from to , turning around and then back from to .

As for summation to the least term, it is natural to chose such that is minimal. Since satisfies a bound of the form

at infinity, this means that we should take , where

Our main aim is to prove the error bound

for . When taking , this bound further simplifies to

2.2.The expedited approximation

The truncated Laplace transform

Let for and for , so that

Since we know how to compute a bound for on the contour , we may compute an explicit bound of the form

(1)

for on a small positive sector near zero.

Borel transforms of at other critical times

For , we define

where

We may also represent as the analytic Borel transform of with respect to . Using the bound (1), this allows us to compute a bound

(2)

for .

The difference between and

Let . For , we also define

For , and setting , we thus have

for . One major topic of this section will be to compute bounds at the origin for and .

Behaviour of the at infinity

Let . Combining the bound (2) with the superexponential bound for , as provided by the accelero-summation process, we may compute a bound

(3)

for . Notice that we may also compute a bound

(4)

for and .

Majorants for specific accelerates and Laplace transforms

The following bounds will be useful for proving precise error estimates for the and . The proofs are a routine application of the saddle point technique.

Lemma 2. Let and be critical times with . Then

(5)

for all . If , then

(6)

for all .

2.3.The first acceleration

Lemma 3. Let . We can compute such that, for , we have

Proof. We have

Using (3) and (4), it follows that

on an interval for some computable .

2.4.Subsequent accelerations

Lemma 4. For each , we can compute a constant , together with a bound

for .

Proof. We will prove the lemma by induction over . We have

where

If , then Lemma 3 yields a bound

for . For , the induction hypothesis yields the bound

for . Now, in view of (6), we may compute a suffiently small such that

for all . If , then a similar computation yields

for all , modulo a decrease of if necessary. Putting these bounds together, we obtain

for . Using (3) and (4), we may also compute a bound

for , modulo a further decrease of if necessary, and where we used the fact that . Combining the bounds for and , the result follows.

2.5.The final Laplace transform

Lemma 5. For any aperture , we can compute a and a bound

for all .

Proof. We have

where

Using Lemma 4 and (5), we can compute and a bound

for . Using (2) and the exponential bound for as provided by the accelero-summation process, we may also compute a bound

for . Modulo a further increase of if necessary, this allows us to compute a bound

for . Adding up the bounds for and , the results follows.

Corollary 6. For any aperture , and assuming that

we can compute a and a bound

for all .

Proof. This directly follows from the fact that .

3.Uniform complexity on local sectors

3.1.Uniform complexity of accelero-summation

Proposition 7. Let be a fixed aperture and let . Then we can compute a constant with the following property: given and with

and , we can compute an approximation with in time , where the complexity bound holds uniformly in under the above conditions.

Proof. Recall that we may compute an exponential bound

for at infinity. For and , this yields a bound

We now wish to compute by approximating the truncated Laplace integral

(7)

with precision , i.e. and .

Let us first consider the case when the bitsize of is bounded by . Under the assumption that , we observe that . This implies that we can chose the contour to use a circle of fixed radius around the origin (which does not depend on ). We next evaluate (7) using the algorithm from [9, Section 6]. Our hypothesis that implies that the primitive of satisfies a holonomic equation of size , uniformly in . Consequently, it can be checked that the complexity bound from [9] holds uniformly in . This means that the required -approximation of can be computed in time , uniformly in .

For general , we approximate in two steps. Let be the growth rate of the linear differential equation satisfied by at the origin. In [9, Theorem 5.2], we showed that in the sector , we have the following bound for the transition matrix on a straightline path in :

For , it follows that

(8)

Now let be an approximation of with and . By what precedes, we may compute -approximations of in time , uniformly in . Using the usual “bitburst” algorithm from [3, 7, 11], together with (8), it follows that we may compute a -approximation of using an additional time of , uniformly in . Adding up these complexity bounds, the result follows.

3.2.Uniform complexity of summation until the least term

Proposition 8. Let . Then we may compute a constants such that

for all , and with

and . Moreover, if and , then we can compute an approximation with in time , where the complexity bound holds uniformly in under the above conditions.

Proof. Direct consequence of [9, Theorem A.1].

3.3.Uniform complexity of expedito-summation

Proposition 9. Let be a fixed aperture and let , where . Then we may compute a constant such that

for all and with

and , where

Moreover, if and , then we can compute an approximation with in time , where the complexity bound holds uniformly in under the above conditions.

Proof. Our hypothesis on implies that

By Corollary 6, it follows that for all with and , we have

For any suitable point close to the origin and , we have shown in [9] how to compute decimal digits of in time . This provides us with the required initial conditions for the analytic continuation of the integrant of the truncated Laplace integral

In a similar way as in the proof of Proposition 7, we may therefore approximate to precision in time , where the complexity bound is uniform in under our conditions.

3.4.The combined local strategy

Putting Propositions 7, 8 and 9 together, we obtain:

Proof. Let be as in Propositions 8, 9 and 7. For any with and , at least one of the following three statements holds:

1. We have .

2. We have for some .

3. We have .

In these cases we respectively apply Proposition 7, 9 or 8 in order to obtain the desired result.

4.Globally efficient evaluation

4.1.Local analysis of cancellations

Assume that is singular at the origin. Then for some , there exists a basis of formal solutions of the form

(9)

for , where , , , and where for some . Moreover, each belongs to the subset of accelero-summable series.

For each fixed accelero-summation scheme, there exist , and such that the and give rise to analytic functions and on the sector . A sector for which this happens is said to be admissible. Moreover, there exist a finite number of admissible sectors with whose interiors cover a small neighbourhood of . We will call this an admissible cover.

Let be one of the sectors in an admissible cover and let and denote the accelero-sums of and on this sector. For each , let . Let denote the subset of all such that

More generally, given a permutation of , let denote the subset of all with

Clearly, .

Let be a non zero solution to on and let be the column vectors with entries . Although can vanish on due to cancellations among the terms and , the vector cannot vanish unless . We will now prove a stronger version of this observation by showing that the sup-norm of cannot become much smaller than .

Theorem 11. There exist constants and such that

for all .

Proof. Without loss of generality, we may assume that . For each , let be the Wronskian matrix

We may decompose

where

and where the entries of are in for some that only depends on the degrees of . It follows that

where by the linear independence of . Now and the entries of are all elements of . It follows that there exists a constant such that for all .

Now consider our fixed linear combination and let

where , so that and . Also let be the column vector with entries . For the sup-norm on vectors, the above discussion shows that

For some fixed constant , this means that

(10)

There also exist constants and such that for all and ,

(11)

Now we may partition into subsets as follows. By induction over , we define to be the subset of all such that

where we understand that . If , then it follows that

Still for , the relation (10) also implies the existence of an such that

Using (11), it follows that

whence

We conclude that for and , using our assumption that .

Remark 12. It is plausible that a bound for can be stated in terms of and the degrees of . We have not pursued this line of thought any further since any constant will do for our purposes.

4.2.Existence of zeros on disks

Consider the power series expansion of at . For each , let be the vector with entries . Theorem 11 provides us with a uniform lower bound for in terms of . We also have the following upper bound for the remaining coefficients.

Lemma 13. There exist constants , and such that

for all with .

Proof. Since is holonomic, there exists a matrix with coefficients in such that

Consequently, there exists a uniform majorant equation for of the form

for suitable constants and , and where denotes the matrix whose coefficients are all one. Taking to be the vector with entries , it follows that is the vector with entries . By construction .

Proof. Let . We may factor with . Let be such that for all . Then we have for all with , whence . By Rouché's theorem, it follows that and admit the same number of zeros in . Hence admits at least one zero inside .

Lemma 15. There exist positive constants , and such that

for all and with .

Proof. Let and be as in Theorem 11 and , and as in Lemma 13. Take . We thus have and for all . Let . Then it follows that and for . Assuming that , we also obtain . Decreasing if necessary, we may arrange ourselves so that . Consequently, satisfies for . We now conclude by Lemma 14.

4.3.Global uniform complexity bounds

We are now in a position to prove our main theorem. We start with proving the uniform bound on “super-admissible” sectors near singularities. Here the sector is said to be super-admissible if we may take in Lemma 15, as well as in the analoguous statement on for each permutation of . Given and with , we will say that is an -approximation of .

Proof. Let and denote the accelero-sums of and on . By Theorem 10, we may compute -approximations of the evaluations in time , uniformly for . In particular, the constants with can be evaluated with a precision of bits in time .

For a given , we first determine a permutation such that . Modulo a permutation of the basis elements , we may assume without loss of generality that . In order to evaluate at , we perform tentative evaluations at increasing bit precisions until the desired approximation with a relative precision of bits is found. For the tentative evaluations, we proceed as follows:

• We compute -approximations of .

• We compute -approximations of .

• Summing up, we obtain a -approximation of .

If the -approximation of has a relative precision of at least bits, then we obtain using one final multiplication with a floating point approximation of . If has a smaller relative precision, then we set and keep iterating.

Now whenever both and , Lemma 15 implies that . In other words, the iteration will stop whenever . Since , this happens for . Since we double at every iteration, the total running time is dominated by the running time of the last tentative evaluation at precision . The most expensive step of this tentative evaluation is the computation of the -approximations of . By Theorem 10, this can be done in time , uniformly in .

Proof of Theorem 1. Let be one of the singularities and let be an admissible ball cover in the neighbourhood of . For each admissible sector and each connected component of (there are at most two such connected components), we also arbitrarily pick a point in . We may compute -approximations for in time . These values may be used as initial conditions for on .

For sufficiently close to , we use the following algorithm for the evaluation of . Among the sectors that contain , we pick the one for which is maximal. In particular, for some fixed constant . Let be the connected component of of that contains . We now evaluate using the algorithm from Lemma 16, by using the initial conditions for at . Applying Lemma 16 on each of the sectors , we obtain a constant such that can be approximated with a relative precision of bits in time , uniformly in , provided that .

Considering the change of variables , one may prove in a similar way that, for some sufficiently large , we can approximate with a relative precision of bits in time , uniformly in , provided that .

Let . The complement is a compact set that contains none of the singularities of . Using the complexity bounds from [7], it follows that a -approximation for can be computed in time , uniformly in . Now admits only a finite number of zeros on and each zero has multiplicity at most . Considering the local power series expansions around any of these zeros , we observe that for some computable contant and sufficiently close to . Provided that , this implies that we can also compute an approximation for with a relative precision of bits in time , uniformly for .

5.Further thoughts and challenges

There are several directions in which the results of this paper can be extended or made more precise.

More general constants

In our main Theorem 1, we assumed that is the field of algebraic numbers. Following the Chudnovsky's [3], and using the baby-step-giant-step technique, one may replace with more general effective subfield of whose elements can be approximated fast. More precisely, if for any constants in we can compute a -approximation of in time , then Theorem 1 still holds, but one should replace the uniform complexity bound by .

Riemann surfaces

In this paper, we used for the domain of our holonomic function . Of course, is really defined on the covering space of which is a Riemann surface. Points on this Riemann surface can be represented by broken line paths as in [7, 8, 9]. By using a suitable size function for broken line paths with vertices in , one may extend Theorem 1 to the evaluation of at points above on this Riemann surface.

Fast approximation of zeros

Given a sufficiently good approximation of a zero of of multiplicity (we must have ), we may use Newton method's to compute a better approximation . Since the evaluations of and can be done with good uniform complexity, this should make it possible to compute a -approximation of in time , uniformly in under suitable conditions. It would be a useful contribution to prove a more precise statement of this kind.

Ball evaluations

In this paper, we assumed that the points where we evaluate are exactly known. An interesting question concerns the efficient computation of high quality ball lifts of . In that case, the evaluation point is replaced by an explicit ball with and , and the evaluation should be a similar ball with the property that and contains two points with distance at least . It would be worthwhile to extend Theorem 1 to this kind of arithmetic.

Multi-summation

Bibliography