
In a series of previous articles, we have given efficient algorithms for the evaluation of holonomic functions over the algebraic numbers and for the computation of their limits at singularities. The focus of these articles was mainly on the efficient evaluation at a fixed point. In the present note, we will show that there exist uniformly efficient algorithms for evaluating holonomic functions. The main technical difficulty is to maintain uniform efficiency near irregular singularities. We will introduce a variant of acceleratosummation for this purpose that we call “expeditosummation”.

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 nonsingular 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.
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 accelerosummation 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 wellknown 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 wellknown 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 accelerosummable 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 accelerosummation process with critical times , , and all integrals taken on the positive real axis. The accelerosum is defined in some sector
for any with .
For any fixed , the accelerosummation 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: accelerosummation 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 expeditosummation in Section 2. Roughly speaking, we perform the accelerosummation 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 expeditosummation until the critical time allows us compute digits of in time , uniformly in provided that .
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
Throughout this and the next section, we make the following assumptions:
with is a formal solution to .
is accelerosummable 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.
In [9], we provided a detailed analysis of two summation methods of . The usual accelerosummation process associates the accelerosum 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 accelerosummation process. It provides approximations of a precision that correspond to stopping the accelerosummation process at the first singularity for the first critical time. It is natural to consider more general shortcuts, where we perform the usual accelerosummation 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
Since we know how to compute a bound for on the contour , we may compute an explicit bound of the form
for on a small positive sector near zero.
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
for .
For , and setting , we thus have
for . One major topic of this section will be to compute bounds at the origin for and .
for . Notice that we may also compute a bound
for and .
for all . If , then
for all .
Proof. We have
Using (3) and (4), it follows that
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 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 accelerosummation process, we may also compute a bound
for . Modulo a further increase of if necessary, this allows us to compute a bound
we can compute a and a bound
for all .
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
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
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.
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
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:
We have .
We have for some .
We have .
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 accelerosummable series.
For each fixed accelerosummation 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 accelerosums 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 supnorm of cannot become much smaller than .
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 supnorm on vectors, the above discussion shows that
For some fixed constant , this means that
There also exist constants and such that for all and ,
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
Remark
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.
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 all and with .
We are now in a position to prove our main theorem. We start with proving the uniform bound on “superadmissible” sectors near singularities. Here the sector is said to be superadmissible 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 accelerosums 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.
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 .
There are several directions in which the results of this paper can be extended or made more precise.
[1] W. Balser. From divergent power series to analytic functions. Theory and application of multisummable power series, volume 1582 of Lect. Notes in Math. SpringerVerlag, Berlin, 1994.
[2] R. P. Brent. The complexity of multiprecision arithmetic. In Complexity of Computational problem solving, 1976.
[3] D. V. Chudnovsky and G. V. Chudnovsky. Computer algebra in the service of mathematical physics and number theory (Computers in mathematics, Stanford, CA, 1986). In Lect. Notes in Pure and Applied Math., volume 125, pages 109–232, NewYork, 1990. Dekker.
[4] J. Écalle. L'accélération des fonctions résurgentes (survol). Unpublished manuscript, 1987.
[5] J. Écalle. Introduction aux fonctions analysables et preuve constructive de la conjecture de Dulac. Hermann, collection: Actualités mathématiques, 1992.
[6] B. Haible and T. Papanikolaou. Fast multipleprecision evaluation of elementary functions. Technical Report TI7/97, Universität Darmstadt, 1997.
[7] J. van der Hoeven. Fast evaluation of holonomic functions. TCS, 210:199–215, 1999.
[8] J. van der Hoeven. Fast evaluation of holonomic functions near and in singularities. JSC, 31:717–743, 2001.
[9] J. van der Hoeven. Efficient accelerosummation of holonomic functions. JSC, 42(4):389–428, 2007.
[10] E. A. Karatsuba. Fast evaluation of Bessel functions. Integral Transforms and Special Functions, 1(4):269–276, 1993.
[11] Marc Mezzarobba. Autour de l'évaluation numérique des fonctions Dfinies. PhD thesis, École polytechnique, 2011.
[12] H. Poincaré. Sur les intégrales irrégulières des équations linéaires. Acta Math., 8:295–344, 1886.