Truncation Bounds for Differentially Finite Series

We describe a exible symbolic-numeric algorithm for computing bounds on the tails of series solutions of linear di erential equations with polynomial coe cients. Such bounds are useful in rigorous numerics, in particular in rigorous versions of the Taylor method of numerical integration of ODEs and related algorithms. The focus of this work is on obtaining tight bounds in practice at an acceptable computational cost, even for equations of high order with coe cients of large degree. Our algorithm fully covers the case of generalized series expansions at regular singular points. We provide a complete implementation in SageMath and use it to validate the method in practice. Table of contents

In order for the method to yield arbitrarily tight enclosures, the bound B(ζ, N ) should at least tend to zero as N grows.
In the present work, we are interested in the case where the coefficients u n of the series (1.1) are generated by a linear recurrence with polynomial coefficients. In other words, there exist rational functions b 1 , . . . , b s over some subfield K ⊆ C such that, for large enough n, the coefficient u n is given by (1.3) u n = b 1 (n)u n−1 + · · · + b s (n)u n−s .
Functions with this property are called differentially finite (D-finite) or holonomic. Under this assumption, computing the coefficients and the partial sum given sufficiently many initial values is not a problem in principle (though one needs to bound the round-off error if the computation is done in approximate arithmetic, and numerical instability issues may occur). The question we consider here is the following: assuming that the series (1.1) is given by means of the differential equation (1.4) and an appropriate set of initial values, how can we compute bounds of the form (1.2) on its remainders?
Theoretical answers to this question have long been known. In fact, textbook proofs of existence theorems for solutions of linear ODEs implicitly contain error bounds whose computation often can be made algorithmic, even in the nonlinear case. More specifically, the basic idea of the very method we use here can be traced back to Cauchy's 1842 memoir [Cau42] on the "calcul des limites", or method of majorants as it is now called. (See Henrici [Hen86, Section 9.2] or Hille [Hil97, Sections 2.4-2.6] for modern presentations of the classical method of majorants.) Yet, the basic method is not enough, by itself, to produce tight bounds in practice at a reasonable computational cost.

Contribution
The goal of the present article is to describe a practical and versatile algorithm for these bound computations. Our method fully covers the case of expansions at regular singular points of the ODE (see Section 5.1 for reminders), including expansions in generalized series with non-integer powers and logarithms. (It does not apply in the irregular case, even to individual solutions that happen to be convergent.) It is designed to be easy to implement in a symbolic-numeric setting using interval arithmetic, and to be applicable to fast evaluation algorithms based on binary splitting and related techniques [CC90]. We provide a complete implementation in the SageMath computer algebra system, and use it to assess the tightness of the bounds with experiments on non-trivial examples.
To see more precisely what the algorithm does, consider first the following inequality that one might write to bound the remainder of the exponential series with pen and paper: We are aiming for bounds of a similar shape: like this one, our bounds decompose into a factor that mainly depends on the analytic behavior of the function of interest (actually, in our case, of the general solution of the ODE), multiplied by one that is related to the first few neglected terms of the series. Also like this one, they are parameterized by the truncation order and the evaluation point; however, we focus on the evaluation algorithm giving the bound as a function of N and ζ rather than expressing it as a formula.
Internally, what the algorithm actually computes for fixed N is a majorant series (Definition 3.2 below) of the remainder u N : (z). Given such a series, one readily deduces bounds on u N : (ζ) for a given ζ, but also on related quantities like remainders of the derivative u (z) or higher-order remainders u N : (z), N N . This last feature means that the method provides both "a priori" and "a posteriori" bounds: the knowledge of the coefficients u 0 , . . . , u N −1 (for large enough N ) can be used to bound u N : (ζ) for any N N , and these bounds already tend to zero as N increases, but taking N closer to N leads to tighter bounds, as illustrated on Figure 8.1.
For the method to be of any interest, the majorant series we are computing need to be significantly simpler than the series we are to bound. We seek hyperexponential majorant series, that is, majorant series of the form ∞ n=0û n z n = exp z 0 a(w) dw, a(z) ∈ R(z), or, equivalently, series that satisfy linear ODEs of the form (1.4) but of order r = 1. The algorithm to compute these hyperexponential majorant series is based on a novel combination of two classical ideas: (1) estimating the error in the solution of linear equations using residuals, (2) bounding the solutions of ODEs with analytic coefficients by the method of majorants. An analogy with an elementary situation might be helpful to see how this works. Consider the linear system Ax = b, for some invertible matrix A ∈ C r×r . Suppose that we have computed an approximationx of the exact solution x and we want to bound the approximation error x − x . Suppose, additionally, that we are given an M such that A −1 M . Writing all we need to conclude is an approximation of the residualb − b, which is easy to obtain. In our setting, the residuals can be computed from a small number of terms of the coefficient sequence of the series, while the method of majorants provides the required "bound on the inverse" of the differential operator.
Continuing with the analogy, the more effort we spend computing M , the tighter we can make the final bound. Nevertheless, even if the inequality A −1 M is loose, as soon as the residual is computed accurately enough, the bound (1.6) overestimates the actual error by a constant factor only asx tends to x for fixed A. In the differential case as well, tighter bounds come at a higher cost even for fixed N . Our algorithm contains additional parameters that can be adjusted to influence the trade-off between computational cost and accuracy. Thus, the same general algorithm can be used to obtain simple but sufficient bounds in "easy" cases or, with more effort, better bounds in "hard" cases where the simpler one are unusable.

Outline
The rest of this article is organized as follows. Section 2 compares the approach taken here with existing work and discusses some applications. Section 3 presents the notation used in the sequel while recalling a few classical results. Section 4 sketches our algorithm and the inequalities on which it relies, under some simplifying assumptions. The reader only interested in the general idea can stop there. The following sections are concerned with technical details and proofs related to the bound computation algorithm. In Section 5, we prove a majorization theorem that covers expansions at regular singular points, laying the foundation for the detailed description of our main algorithm in Section 6. We then consider the practical aspects of computing certain intermediate bounds on sequences defined by rational functions in Section 7, and the derivation of various kinds of concrete numerical bounds from the output of the main algorithm in Section 8. Finally, Section 9 describes experiments that illustrate the quality of these bounds and the effect of the tuning parameters.

Related work
While articles directly comparable to the present one are relatively few, similar questions appear naturally as sub-problems of various computational tasks. Somewhat arbitrarily, we group previous work by context in two categories: the evaluation of classical functions, and the numerical solution of ordinary differential equations.

Differentially finite functions as special functions
An important application of the summation of power series where tail bounds are needed is the rigorous multiple-precision evaluation of special functions [BZ10]. For a fixed function, it is usually not too hard to derive good ad hoc bounds [MPF ].
Bounds that cover wider classes of functions become more complicated as the number of parameters increases, yet adequate tail bounds are available in the literature and used in practice for common special functions depending on parameters. For example, Du and Yap [DY05, Section 3] or Johansson [Joh16,Section 4.1] present bounds that cover general hypergeometric functions.
From this perspective, our goal is to describe an algorithm for bound computations that applies to the "general differentially finite function" of a complex variable, viewed as a special function parameterized by the complete coefficient list of the differential equation (1.4). Van der Hoeven [vdH99, Section 2] [vdH01, Section 2.4] [vdH03, Section 3.5] already gave several such algorithms, some based, like ours, on the method of majorants, and mentioned (in a different context) a variant of the idea of combining majorants with residuals in order to obtain tighter bounds [vdH01,Section 5.3]. Like the theoretical algorithms read between the lines of existence proofs mentioned in the Introduction, van der Hoeven's algorithms lack details and suffer from overestimation issues that make them unsatisfactory in practice. Our method can be seen as a refinement of these ideas yielding practical bounds.
The present author already considered related tightness issues in an article with B. Salvy [MS10] on asymptotically tight bounds on linear recurrence sequences, with some extensions (including a few of the ideas developed in the present paper) in the author's doctoral thesis [Mez11]. While the focus of that article was not on tails of power series, the remainder bounds on tails of power series that came as a corollary seemed to perform well in experiments with simple special functions. Unfortunately, they later proved too pessimistic and costly to compute with more complicated equations. The setting of the present work is different in that we do not insist on asymptotic tightness, nor on producing human-readable formulae, and focus instead on making the bounds easy to evaluate numerically and as tight as possible for realistic values of the parameters.

Interval methods for ODEs
Directly summing the series (1.1) to compute u(ζ) only works when ζ lies within its disk of convergence D. When that is not the case, an option is to evaluate the function by "numerical analytic continuation", that is, to first use (1.1) to compute the Taylor expansion of u at some point ζ 1 ∈ D closer to the boundary, then use that series to compute the expansion of u at a point ζ 2 , possibly outside of D, and so on until we reach ζ. Because u satisfies the differential equation (1.4), it is enough at each step i to compute the first r derivatives of u at ζ i+1 using the Taylor expansion at ζ i , and the rest of the expansion follows from the differential equation.
We can also take a slightly different perspective and view the evaluation algorithm as a numerical ODE solver. Seen from this angle, what we just outlined is the numerical solution of (1.4) by a Taylor method. Taylor methods are among the oldest numerical methods for ODEs. While often considered too costly for machineprecision computations, they remain the methods of choice for high precision and for interval computations [BRAB11,NJC99]. In addition, they can be made particularly efficient in the case of ODEs of the form (1.4) [CC90, Sections 2-3], and the efficient variants generalize to deal with singular problems [vdH07].
Rigorous Taylor methods (not limited to the linear case) have attracted a lot of interest from the field of Interval Analysis, starting with Moore's 1962 PhD thesis [Moo62, Section 6]. Since a Taylor method reduces the solution of an ODE to the summation of series expansions of (derivatives of) solutions, the problem of bounding the remainders occurs naturally. Instructive overviews of the ideas used in this context can be found in the literature [Rih94,NJC99,NJN07]. A popular approach is to start with a coarse enclosure of the solution and its first derivatives over the whole time step, obtained by a fixed-point argument based for instance on Picard's iteration. One then deduces an enclosure of the derivative of order N using the ODE, and a bound for the remainder of order N by Taylor's formula. If, instead of values on a grid, one computes a polynomial approximation of the solution over some domain, it is also possible to apply a fixed-point argument directly to a neighborhood of this approximation in function space, by evaluating an integral form of the differential equation in polynomial interval arithmetic (using Taylor models for instance). The way our algorithm uses residuals bears some conceptual similarity to this last approach, but our method does not require an explicit polynomial approximation of the solution as input.
Still in the context of Taylor methods, but closer in spirit to van der Hoeven's approach or ours, Warne et al. [WWS + 06] develop explicit majorants for solutions of nonlinear differential equations. Neher [Neh99,Neh01] considers the case of linear ODEs with analytic coefficients and bounds the coefficients of the local series expansions of the solution by geometric sequences, starting from bounds of the same type on the coefficients of the equation. In the special case of ordinary points, our method might be described as a blend of Neher's approach, simplified by the language of majorant series, with elements of the "continuous approximation" strategy.

The regular singular case
To the best of our knowledge, the present work is the first that directly applies to logarithmic solutions at regular singular points (with an early version of some of the ideas in the author's dissertation [Mez11, Section 6.4]). The other methods mentioned above are limited to series solutions, often at ordinary points only, with the exception that van der Hoeven [vdH01, Section 3.3.2] sketches an adaptation to the case of logarithmic solutions of the earliest of his algorithms for plain power series. The same author [vdH03, Section 5.2] considers majorants for solutions of linear differential systems at regular singular points, which provides an indirect way of handling logarithmic solutions.

Formal power series
For any commutative ring R, we denote by R [[z]] the ring of formal power series over R in the variable z, and by R((z)) the ring of formal Laurent series. Given f ∈ R((z)), we denote by f n or [z n ]f (z) the coefficient of z n in f (z). As in Equation (1.2), we also write f n: (z) for the remainder i n f i z i . More generally, f ν is the coefficient of z ν in a generalized series ν∈Λ f ν z ν indexed by Λ ⊂ C. In this context, Latin letters denote integers, while Greek coefficient indices can be arbitrary complex numbers. We sometimes identify expressions representing analytic functions with their series expansions at the origin. Thus, for instance, [z n ](1 − αz) −1 is the coefficient of z n in the Taylor expansion of (1 − αz) −1 at 0, that is, α n .
We also consider polynomials in log(z) with formal power (or Laurent) series coefficients. The space R((z))[log z] of such expressions embeds in R[log z]((z)), and for f ∈ R((z))[log z], the notations f n and f n: refer to the coefficients of f viewed as an element of R[log z]((z)).

Differential equations and recurrences
We assume some familiarity with the classical theory of linear differential equations with analytic coefficients, as described for instance in Henrici's or Hille's books [Hen77,Hil97], and with the basic properties of D-finite formal power series, for which we refer the reader to Stanley [Sta99, Section 6.4] or Kauers and Paule [KP11]. In particular, we will freely use the properties summarized below.
We usually write linear differential equations a r (z)y (r) (z) + · · · + a 0 y(z) = q(z) in operator form, i.e. as L · y = q where L is a linear operator acting on some space of functions of interest. When the coefficients a k are formal Laurent series, such an operator can be written as L = P (z, D z ), that is, as a polynomial P (X, Y ) ∈ C((X))[Y ] evaluated on the operators z : y → (w → wy(w)) that multiplies a function by the identity function and the standard differentiation operator D z : y → y . Note that the operators z and D z do not commute; one can bring the coefficients a k (z) "to the right" of D k z using the relation D z z = zD z + 1 deduced from Leibniz' rule. By abuse of notation, we do not always make the difference between the operator P (z, D z ) and the polynomial P (X, Y ).
Similarly, recurrences b s (n)u n+s + · · · + b 0 (n)u n = v n with polynomial coefficients can be written P (n, S n ) · u = v where n : (u k ) → (ku k ) is the operator that multiplies a sequence by its index, S n : (u n ) → (u n+1 ) is the shift operator, and . The corresponding commutation rule reads S n n = (n + 1)S n . We typically consider bi-infinite sequences (indexed by Z or λ + Z, λ ∈ C), so that we also have at our disposal the backward shift operator S −1 n : (u n ) → (u n−1 ). In the differential case, when working in a neighborhood of the origin, it proves convenient to write differential operators in terms of the Euler operator defined by θ = zD z instead of the standard derivation D z . Any operator L = a r (z)D z + · · · + a 0 (z) ∈ C((z))[D z ] can be written L = L(z, θ) with L =ã r (X)Y r + · · · +ã 0 (X) ∈ C((X))[Y ]. Conversely, an operator L ∈ C[z][D z ] can be rewritten as an element of C[z][θ] at the price of multiplying it on the left by a power of z. This change does not affect the space of solutions of the equation L · y = 0. One can check that the quotient a r (z)/ã(z) of the leading coefficients with respect to D z and to θ is a power of z. "Moving" the coefficients of either form "to the right" of the differentiation operators leaves the leading coefficient unchanged.
It is classical that the coefficients of a linear ODE with power series coefficients obey a recurrence relation obtained by "changing θ to n and z −1 to S n " in the equation. The precise result is as follows, see Proposition 5.1 below for a sketch of a proof in a more general setting.
. A formal power series y ∈ C[[z]] satisfies the differential equation L(z, θ) · y = 0 if and only if its coefficient sequence (y n ), extended with zeroes for negative n, satisfies L(S −1 n , n) · (y n ) n∈Z = 0. Note that in general, the recurrence operator L(S −1 n , n) has infinite order, but since y n vanishes for n < 0, each coefficient of the sequence L(S −1 n , n) · (y n ) only involves a finite number of y n . When L(z, θ) has polynomial coefficients, the corresponding recurrence has finite order.

Majorant series
The language of majorant series offers a flexible framework to express bounds on series solutions of differential equations.
if it is a majorant series of f 0 , . . . , f k . In both cases, we write f (z) f (z).
Series denoted with a hat in this article always have non-negative coefficients, and f is typically some kind of bound on f , but not necessarily a majorant series in the sense of the above definition.
The following properties are classical and easy to check (see, e.g., Hille [Hil97, Section 2.4]). Note the assumption that f and g are plain power series: while convenient to state some results, the extension to f ∈ C[[z]][log z] does not share all the nice properties of majorant series in the usual sense.
be such that f f and g ĝ. The following assertions hold: (c) f n: f n: for n ∈ N, Additionally, the disk of convergenceD off is contained in that of f , and when g 0 ∈D, we have f (g(z)) f (ĝ(z)). In particular, |f (ζ)| is bounded byf (|ζ|) for all ζ ∈D.

A sketch of the method
In this section, we outline how our bounds are obtained, in a simplified setting and without giving detailed algorithms. Our simplified setting covers all solutions at ordinary (non-singular) points of the differential equation, where Cauchy's theorem applies and all solutions are analytic. The contents of this section may be enough for the reader only interested in the general idea of the algorithm. Sections 5 to 7 essentially repeat the same algorithm in the general case, while adding more detail and introducing refinements that help obtain tighter bounds.

Truncated solutions and residuals
We start with a linear differential equation (1.4), with rational coefficients and right-hand side. As observed in Section 3.2, such an equation can always be brought into the form without changing its solution space. The unusual choice of writing the coefficients p k to the right of θ k is not essential, but will prove convenient later. We also assume without loss of generality that the polynomials p 0 , . . . , p r are globally coprime. Most importantly, we assume that the leading coefficient p r does not vanish at 0. This is the case when the origin is an ordinary point of (1.4) (1) , i.e., when none of the fractions a k (z)/a r (z), 0 k < r, has a pole at 0. Let u(z) be a solution of P (z, θ) · u(z) = 0, and consider the truncatioñ u n z n of the series u(z) to some large enough order N 1. Our goal is to obtain an explicit majorant series of the remainder u(z) −ũ(z). This remainder satisfies where q(z), the residual associated to the approximate solutionũ(z), is an explicit, computable polynomial. Becauseũ is a truncation of an exact series solution, the residual q(z) has "small" support: more precisely, it can be checked to be of the form q(z) = q N z N + · · · + q N +s z N +s where s = deg z P (z, θ).

The majorant equation
Setting (4.1) (1) Assuming p r (0) = 0 also allows for a regular singular point at the origin, cf. Section 5.1. However, we restrict ourselves in this section to power series solutions at regular singular points, leaving out a number of technicalities related to non-integer exponents and logarithms.
we get an equation with rational function coefficients in which we immediately rewrite L(z, θ) as by expanding p r (z) −1 in power series and rearranging. Crucially, since p r (0) = 0, the polynomial Q 0 has degree r, while the degree of Q j for j 1 is at most r − 1. By the correspondence of Proposition 3.1, the coefficient sequence (y n ) n∈Z of y(z) then satisfies When n is large enough, Q 0 (n) does not vanish and y n is given recursively by Because of the degree constraints noted above, n q n /Q 0 (n) and n Q j (n)/Q 0 (n) (for fixed j) remain bounded as n grows. Assume for a moment that we have at our disposal boundsq n andâ j such that n q n Q 0 (n) q n < ∞, n n 0 , (4.4) n Q j (n) Q 0 (n) â j < ∞, n n 0 , j 1, (4.5) withâ j = O(α j ) for some α as j → ∞.

Majorant series for the remainders
The general solution of (4.8) reads for an arbitrary constant c ∈ C. We have just seen that anyŷ(z) of this form such that (4.7) holds is a majorant series for y(z). It remains to choose the parametersq,â, and c in such a way that (4.4), (4.5), and (4.7) hold. Recall that, in our context, y(z) = p r (z) (ũ(z) − u(z)) is a power series of valuation at least N 1. Provided that N is large enough, we can assume n 0 = N and take c = 0. Equation (4.7) is then automatically satisfied. Since at most s of the constraints (4.4) are nontrivial, it is not hard to compute explicit valuesq n fulfilling these constraints. Knowing h(z), we can even chooseq(z) as a polynomial multiple (2) of h(z), so that the integral in the expression ofŷ(z) reduces to an explicit polynomialĝ(z).
The case of (4.5) is slightly more involved, as, in appearance, it entails an infinite number of inequalities. However, the polynomials Q j stem from the power series expansions of the finitely many rational functions p k /p r of (4.2). By bounding the numerator and denominator of L(z, θ) separately, we can find a single rational functionâ(z) that satisfies all these inequalities; moreover, we can arrange thatâ(z) has radius of convergence arbitrarily close to that of p r (z) −1 .
A similar (but simpler) computation yields a polynomialp(z) such thatp(z) −1 is a majorant series of p r (z) −1 with the same condition on the radius of convergence. Putting everything together, we obtain an explicit majorant seriesû(z) for the error u(z) −ũ(z), in the form As noted in Proposition 3.3, we have in particular |u(ζ) −ũ(ζ)| û(|ζ|) for all ζ in the disk of convergence ofû(z). Observe that, in (4.10),ĝ(z) can be chosen of valuation N , so thatû(z) itself has valuation N . Moreover, its coefficients can be taken of roughly the same order of magnitude as those of the first neglected terms of the series u(z). The radius of convergence ofâ(z) andp(z) −1 can be made arbitrarily close to the distance from the (2) Since the power series h(z) has nonnegative coefficients and starts with h(0) = 1, replacing any polynomialq(z) computed from (4.4) alone byq(z) h(z) will do. origin to the smallest nonzero singularity of the differential equation (1.4) (or even equal to it, but manipulating this quantity exactly may be costly). In the absence of apparent singularities, this distance is the local radius of convergence of a generic solution of the differential equation. Thus, there is reason to hope thatû(|ζ|) does not stray too far from |u(ζ) −ũ(ζ)| as N and ζ vary.

Majorant equations: the general case
The key step of the reasoning outlined above is the construction in Section 4.2 of a "majorant equation" of the inhomogeneous equation (4.2). Our goal is now to extend this construction to the general case of solutions at regular singular points. To make them applicable to variants of the algorithms adapted to other situations, the results are stated in slightly more general form than actually needed in the sequel.

Regular singular points
Let us start with some reminders on regular singular points. Consider an operator where a 0 , . . . , a r are analytic functions with no common zero. Recall that ζ ∈ C is called a singular point of the operator D when a r (ζ) = 0, and an ordinary point otherwise.
At an ordinary point, by Cauchy's theorem, the equation D ·y = 0 admits r linearly independent analytic solutions. If ζ is a singular point, analytic solutions defined on domains with ζ on their boundary are in general not analytic at ζ. Regular singular points are a class of singular points that share many similarities with ordinary points, starting with the existence of convergent (generalized) series solutions.
A singular point ζ of D is called regular singular when, for all k, its multiplicity as a pole of a k /a r is at most r − k. We say that an operator D is regular at a point ζ (or that ζ is a regular point of D) if ζ is either an ordinary point or a regular singular point of D. An operator written in the formã r (z)θ r + · · · +ã 0 (z) where at least one ofã 0 , . . . ,ã r has a nonzero constant coefficient is regular at the origin if and only ifã r (0) = 0. This is equivalent to saying that the univariate polynomial Q 0 (θ) =ã r (0)θ r + · · · +ã 0 (0), called the indicial polynomial of the operator, has degree r.
When ζ ∈ C is a regular point of D, the equation D · y = 0 admits r linearly independent solutions of the form , where the f i are analytic at ζ. The possible λ are exactly the roots of Q 0 . Formally at least, it is often convenient to rewrite (5.1) as a single series We call expressions of this form logarithmic series, and refer to Henrici [Hen86, Section 9.5] for a rigorous presentation of their algebraic manipulation.

Recurrences on the coefficients of solutions
The coefficients of the logarithmic series expansions of solutions of linear ODEs at regular points are related by systems of recurrence relations whose existence and explicit description go back to Frobenius [Fro73] and Heffter [Hef94, Chapter VIII].
Here we recall these results in a compact formalism [Mez10,Mez11] close to that of Poole [Poo36] which makes them appear as direct generalizations of Proposition 3.1.
where it is understood that y λ+n,k = 0 for n 0 or k K, and similarly for q(z).
where S ν and S k are the shift operators on C (λ+Z)×N mapping a double sequence (u ν,k ) respectively to (u ν+1,k ) and (u ν,k+1 ).
Proof sketch. -By calculating the action of the operators z and θ on logarithmic monomials z ν log(z) k /k!, one checks that the coefficient sequence of The result follows. Observe that, restricted to k k 0 where k 0 is the largest index for which the u ν,k and q ν,k are not all zero, (5.2) reduces to the univariate recurrence L(S −1 ν , ν)·(y ν,k 0 ) ν = (q ν,k 0 ) ν . In particular, if y(z) and q(z) are plain power series (k 0 = 0, λ = 0), we get the formula of Proposition 3.1.
In the general case, we can put the "implicit" equation (5.2) into "solved" form as follows. Consider again an operator , and, similar to what we did in Section 4.2, write Note that Q 0 coincides with the indicial polynomial introduced above. Let µ(ν) denote the multiplicity of a complex number ν as a root of Q 0 . To simplify the notation somewhat, we limit ourselves here to right-hand sides of the form Q 0 (θ)·q(z).
and define y ν,k and q ν,k as in Proposition 5.1. Let κ(n) 0 be such that q λ+n,k = 0 for all k κ(n). Let τ (n) 0 be such that y λ,k = 0 for k τ (0), and Then, the coefficients of y(z) satisfy y λ+n,k = 0 for all k τ (n), and are given by k+t for all n ∈ Z and k ∈ N. (Both sums can be extended to range from the indicated lower index to infinity.) Conversely, any solution (y λ+n,k ) of (5.4) with y λ+n, In general, to satisfy (5.3), one can take Proof. -Consider the equation which results from Proposition 5.1. By extracting the sequence term of index ν on both sides, we get a relation between sequences y ν−j = (y ν−j,k ) k∈N and q ν = (q ν,k ) k∈N of a single index k: Let us first prove by induction that for all n, the sequence y ν where ν = λ + n vanishes under the action of S τ (n) k . This is the case by assumption for n 0. Now assume that this holds true for all n < n.
to (5.6) and using the induction hypothesis givesQ 0 (S k )S denote the canonical representative of its inverse, and write Considering again (5.6) and applying A ν (S k ), we get We have seen that S Finally, for j 1, the sequence S τ (n−j) k · y ν−j is zero, and A ν (X)Q j (ν + X) coincides at least to the order τ (n)−µ(ν) τ (n−j) with the series X µ(ν) Q j (ν +X)/Q 0 (ν +X). Equation (5.4) then follows by extracting the coefficient of index k.
Remark 5.3. -While the recurrence (5.4) giving the whole sequence y ν "at once" is useful in the context of bound computations, specializing (5.6) without "inverting" Q 0 yields an alternative expression of y ν,µ(ν)+k , as a function of the y ν−j , j 1, and the y ν,k , k + 1 k < µ(ν), which can be better adapted to the iterative computation of the y ν,k .
Thus E is a set of cardinality r. As already mentioned, it is well known that the homogeneous equation L(z, θ)·y = 0 has an r-dimensional vector space of convergent solutions spanned by elements of z λ C[[z]][log z] for Q 0 (λ) = 0. Corollary 5.2 suggests a specific choice of basis, giving rise to a dual notion of "initial values".
Proof sketch. -According to Corollary 5.2, the coefficient y ν,k of a solution y ∈ z λ C[[z]][log z] is determined by the y ν−j,k , j 1, as soon as k µ(n). In contrast, (5.4) imposes no constraint on the y ν,k with (ν, k) ∈ E, and it is not too hard to see that taking y ν 0 ,k 0 = 1 for a certain (ν 0 , k 0 ) ∈ E and y ν,k = 0 for (ν, k) ∈ E\{(ν 0 , k 0 )} determines a solution (y ν,k ) ν∈ν 0 +Z,k∈N such that y ν,k = 0 for k τ (ν). The collection of the logarithmic series associated to these sequences for all choices of (ν 0 , k 0 ) is a basis of the solution space of L(z, θ).
For each (ν 0 , k 0 ) ∈ E, Corollary 5.2 applied with λ = ν 0 , κ ≡ 0 and τ (0) = k 0 provides a bound of the form (5.5) on the degree with respect to log(z) of the coefficients of the corresponding element of this basis. One gets a uniform bound valid for all coefficients of all solutions y ∈ z λ 0 C[[z]][log z] by taking λ = λ 0 − 1, κ ≡ 0, τ (0) = 0, and letting n tend to infinity.

The majorant equation
Equipped with these preliminaries, it is now easy to extend the method of Section 4.2 to logarithmic series solutions. The next proposition, in some sense, reduces the problem of bounding the solutions of arbitrary regular ODEs with polynomial coefficients to the case of first-order equations.
Aiming for first-order majorant equations means that the solutions are plain power series without logarithms, and makes it possible to write them down explicitly, while leaving us with enough freedom to match the analytic behavior of solutions of the original equation. A small restriction is that first-order equations cannot capture the asymptotics of entire series of non-integral order. As further discussed in Section 8.3, this limitation is usually of little consequence in practice, and it is possible to circumvent the issue if necessary. The method can also be adapted to other forms of majorant equations.
We still consider an operator L(z, θ) = j 0 Q j (θ)z j and a solution y(z) of ] for a fixed λ ∈ C, and define µ(ν), E, y ν,k and q ν,k as above.
] be a solution of Assume, additionally, that the following assertions hold for a certain n 0 1: (i) for all j 1 and all n n 0 , where τ (n) is a non-decreasing sequence such that y λ+n,k = 0 for k τ (n), (ii) for all n n 0 and all k 0, n |q λ+n,µ(λ+n)+k | q n , (iii) |y λ+n,k | ŷ n for 0 n < n 0 and for all k 0, (iv) |y λ+n,k | ŷ n for n n 0 and 0 k < µ(λ + n). We then have the bound |y λ+n,k | ŷ n for all n, k 0. In other words,ŷ(z) is a majorant series (in the extended sense) of z −λ y(z).
Proof. -We prove the result by induction on n. The inequality holds for n < n 0 by assumption (iii). Fix n n 0 (in particular, n > 0), and assume that |y λ+n ,k | ŷ n for all n < n and k 0. Denote m = µ(λ + n). For k < m, the required inequality is given by assumption (iv). Now consider y λ+n,m+k for some k 0. By Corollary 5.2 applied to the equation L(z, θ) · y = Q 0 (θ) · q, we have It follows that and hence, using assumptions (i)-(ii) and the induction hypothesis, Extracting the coefficient of z n on both sides of the differential equation (5.8) yields so that (5.10) becomes |y λ+n,m+k | ŷ n . This completes the induction step and the proof.
The general solution of (5.8) readŝ for an arbitrary constant c. Observe thatĥ(z) has valuation zero, and henceŷ(z) either has valuation zero as well (if c = 0) or has the same valuation asq(z) (if c = 0). Conditions (i) and (ii) in Proposition 5.5 ensure that the solutions of (5.8) can be used to control those of the equation L(z, θ) · y = q. For the proposition to be applicable,â(z) andq(z) should be chosen so that these conditions hold.
Conditions (iii) and (iv) depend on the specific solution y we are interested in. They express that the bound holds "initially", and most importantly for the generalized initial values described in Corollary 5.4 (which are not determined by the "previous terms" of the series). In particular, condition (iii) becomes trivial if n 0 val(z −λ y(z)), where val(f ) is the valuation of f viewed as an element of C[log z] [[z]]. As for condition (iv), it vanishes as soon as n 0 is larger than all the zeros n of Q 0 (λ + n), or even than the zeros of Q 0 (λ + n) corresponding to nonzero generalized initial values in y(z).
Of special interest is the situation where both (iii) and (iv) are automatically satisfied. This happens in particular when the valuation of z −λ y(z) is larger than that of any h(z) ∈ C((z)) such that z λ h(z) satisfies the homogeneous part L(z, θ)·(z λ h(z)) = 0 of the differential equation on y (i.e., larger than max({ν − λ : (ν, k) ∈ E} ∩ N)), and n 0 is chosen between these two values.
Even in the general case, assumingâ 1 > 0, it is always possible to adjust c or increase selected coefficients ofq(z) so that (iii) and (iv) hold. Therefore, for any choice of n 0 , the other parameters can be selected in such a way that Proposition 5.5 yields a nontrivial bound.

The main algorithm
Let us now go back to homogeneous equations with polynomial coefficients, and specialize the previous results to develop an algorithm that computes bounds on remainders of logarithmic series solutions in this case. The description of the algorithm is split into two parts. The first part, Algorithm 6.1, does not depend on the particular solution or truncation order. Roughly speaking, it aims at satisfying condition (i) of Proposition 5.5, which yields what one might call a "pseudo-inverse bound" on the differential operator. The second part, Algorithm 6.11, uses the result of the first one to bound a specific remainder of a specific solution.

Setting
For the whole of this section, P (z, θ) ∈ K[z][θ] denotes a differential operator, assumed to be regular at the origin, over a fixed number field K ⊂ C. We consider a solution of the logarithmic series expansion of u.
Our goal is to give an algorithm to bound the error u −ũ, based on Proposition 5.5. Note that a solution of P (z, θ) can in general be a linear combination of solutions of the form (6.1) with λ belonging to different cosets in C/Z. We limit ourselves here to a single λ, and will if necessary bound the remainders corresponding to exponents in other cosets λ + Z separately.

"Pseudo-inverse bounds" on differential operators
Algorithm 6.1 below corresponds to the first part of the main bound computation algorithm.
Compared with the outline in Section 4, and besides supporting regular singular points, this version of the algorithm provides more freedom in the choice of the rational functionâ(z), via the input parameter . The method of Section 4 corresponds to = 1. Increasing leads to tighter bounds, at the price of increased computation time. Taking ≈ s/2 often gives good results in practice for evaluations far from the border of the disk of convergence of p r (z) −1 (but see Section 9.4 for more on that). Except for step 3, the algorithm can be run in interval arithmetic and will return finite bounds if all operations are performed with sufficient precision. Besides, the "exact" data needed at step 3, that is the roots of Q 0 in λ + Z and their multiplicities, will typically be known beforehand by the choice of λ.
Algorithm 6.1. -Bound differential operator Input: (To do so, expand the fractions p i /p r in power series up to the order and set If working in approximate arithmetic, force the coefficient of θ r in Q 0 to the exact value 1 and the coefficients of θ r in Q 2 , . . . , Q −1 , U 0 , . . . , U s−1 to 0.) 2. Compute lower bounds 0 < ρ i |ζ i | on the roots ζ 1 , ζ 2 , . . . of p r (z), and a lower bound c on the absolute value of its leading coefficient. Seť where m i is the multiplicity of the root ζ i .

Denoting
use Algorithm 7.4 and Remark 7.5 to get boundsQ 1 , . . . , Remark 6.2. -A simple strategy for step 2 (cf. Grégoire [Gré12]) is to compute a single tight lower bound ρ on the smallest root of p r by the method based on Graeffe iterations of Davenport and Mignotte [DM90], and takep(z) = c (ρ − z) deg p . This method is fast and typically yields good bounds in simple cases like standard special functions. Unsurprisingly, though, it can also lead to large overestimations for operators with many singularities. At the other extreme, one can call a (rigorous) numerical root finder to obtain arbitrarily tight bounds. While, on paper, this option complexifies the algorithm a lot, very good implementations of complex root-finding exist, so that it works well in practice. It would also likely be feasible to extend the method of Davenport and Mignotte to get better estimates of larger roots. Proposition 6.3. -Given an operator P (z, θ) = p r (z)θ r + · · · , regular at 0, Algorithm 6.1 returns a pair consisting of a polynomialp ∈ R[z] such that p r (z) −1 p(z) −1 , and a rational functionâ(z) withâ(0) = 0 that satisfies condition (i) of Proposition 5.5 for any solution Proof. -At step 1, the fractions p k /p r can be expanded in power series since p r (0) is not zero. The operator Q(z, θ) defined by taking [θ k ]Q j = [z k ](p k /p r ) for 0 j < then agrees with P (z, θ)p r (z) −1 up to the order O(z ). In particular, it is monic with respect to θ, so that [θ r ]Q 0 = 1 and deg Q j r − 1 for j 1. For the same reason, the difference P (z, θ) − Q(z, θ) p r is of the form U (z, θ)z , with deg U j r − 1 for all j since 1. The operator Q(z, θ) p r has degree (with respect to z) at most − 1 + s, and P (z, θ) has degree at most s, hence the degree of their difference is strictly less than + s, and the coefficients we compute correctly represent U (z, θ).
The polynomials Q 0 , . . . , Q −1 computed here coincide with those defined before by the series expansion L(z, θ) = j Q j (θ)z j . Additionally, the remaining coefficients of the expansion are related to U (z, θ) by After step 2, we have for all i p(z) −1 . By Corollary 5.2, setting τ (n) = n k=0 µ(λ + n) at step 3 guarantees that u λ+n,k = y λ+n,k = 0 for all k τ (n). When the origin is an ordinary point of P (z, θ), though, this generic choice only gives τ (n) = r for n r − 1, but we can take τ (n) = 1 since all solutions of P (z, Let us now show that, when the bounds from steps 2 and 4 hold, the coefficientsâ j of the rational seriesâ (z) =Q(z) + z Û (z) p(z) satisfy condition (i) of Proposition 5.5, that is, F (Q j , n) â j for all j 1 and n n 0 . The first − 1 inequalities of step 4 say that this holds true for the initial coefficientsâ 1 =Q 1 , . . . ,â =Q . The last s inequalities translate into (6.5) According to (6.4), we have for fixed n ∞ j=0 Q +j (λ + n + X) Using the bound p r (z) −1 p(z) −1 from step 2 (and Proposition 3.3), it follows that for all t 0. By summing over 0 t < τ (n) and multiplying by n, then applying (6.5), we get Remark 6.4. -In our implementation, the code corresponding to Algorithm 6.1 actually does not take n 0 as input. Instead, the boundsQ j andÛ j of step 4 are functions of n 0 , and the algorithm returns an object (called a DiffOpBound) that can then be evaluated at a given n 0 to get the pair (p,â). The DiffOpBound objects serve to share part of the computation when one needs to bound several tails of the same series-for example, to keep adding terms to the sum until the bound becomes small enough-or tails of several solutions of the same equation.
Similarly, when Algorithm 6.1 is run on the same operator with increasing values of , almost all of the previous computations can be reused from one run to the next. We use that in the implementation to provide refinable bounds: DiffOpBound objects start off as relatively coarse and inexpensive bounds, and offer a method to increase if the bounds turn out to be too pessimistic. The switch from the "fast" strategy for step 2 discussed in Remark 6.2 to the more accurate one works in the same way.
Remark 6.5. -We also experimented with a variant of Algorithm 6.1 that performs a partial fraction decomposition with respect to z of the operator instead of splitting it into a truncated series and a rational remainder term, but did not manage to obtain better bounds using this approach.
It is not hard to compute a majorant seriesû(z) of u(z) from the output of this algorithm, as illustrated by the following corollary (not used in the sequel).
Corollary 6.6. -Let (p,â) be the output of Algorithm 6.1 called on the input (L, λ, , n 0 ). For any solution u(z) of P (z, θ) · u = 0 lying in is a majorant series (in the extended sense of Definition 3.2) of u(z).

ANNALES HENRI LEBESGUE
Proof. -Apply Proposition 5.5 with L(z, θ) = P (z, θ)p r (z) −1 and q(z) = 0. By the previous proposition, Algorithm 6.1 provides a seriesâ(z) satisfying the first condition. The same is true a fortiori ofâ(z) + z. In addition, the correspondinĝ y(z) = exp z 0 (1+w −1â (w)) dw has strictly positive coefficients. The second condition is trivial for homogeneous equations. We choose c large enough that |y λ+n,k | cŷ n for n < n 0 or k < µ(λ + n), fulfilling the last two constraints. Proposition 5.5 then implies y(z) cŷ(z), and hence u(z) = y(z)/p r (z) û(z). The remainders ofû(z) are majorants of the corresponding remainders of u(z) and, whenû converges at |ζ|, provide us with a sequence of bounds on the error |u(ζ)−ũ(ζ)| that tends to zero as the truncation order increases. While the remainders ofû(z) may not admit nice closed-form expressions, it is not hard to extract reasonably good bounds on their values algorithmically, as discussed in Mezzarobba and Salvy [MS10] or in Section 8.2 below, thus solving our problem in principle. We can do better, though, by taking into account the residual associated to u(z).

Normalized residuals
Like with Corollary 5.2 and Proposition 5.5, instead of working with the residual P (z, θ) ·ũ(z), it is convenient to introduce the normalized residual q(z) defined by with the additional condition that q ν,k = 0 for (ν, k) ∈ E. That the normalized residual is well defined follows from the following lemma. When 0 is an ordinary point and N r, or more generally whenũ ∈ C[[z]] and Q 0 (n) = 0 for all n N , the normalized residual is simply the residual with the coefficient of z n divided by Q 0 (n), cf. (4.4).
Lemma 6.7. -Given f ∈ z λ C[[z]][log z] of degree at most K with respect to log(z), the differential equation Q 0 (θ) · q(z) = f (z) has a unique logarithmic series solution Proof. -This is a consequence of Proposition 5.1, applied to the operator Q 0 (θ) (which has the same associated set E as P (z, θ)) and the right-hand side f (z). More precisely, the recurrence (5.2) for ν ∈ λ+Z reduces in this case to Q 0 (ν +S k )·(q ν,k ) = (f ν,k ), that is, For all ν, this equation has a (unique) solution (q ν,k ) k with q ν,k = 0 outside of the interval The corresponding logarithmic series q(z) ∈ z λ C[[z]][log z] satisfies Q 0 (θ) · q = f . By Corollary 5.4 applied to the homogeneous equation Q 0 (θ) · q = 0, it is the only solution of Q 0 (θ) · q = f with q ν,k = 0 for all (ν, k) ∈ E.
Now consider the special case of f (z) = P (z, θ) ·ũ(z) used to define the normalized residual.
Lemma 6.8. -The logarithmic series q(z) of Lemma 6.7 corresponding to f (z) = P (z, θ) ·ũ(z) is of the form where s = deg z P (z, θ) and K is the largest power of log(z) appearing inũ(z).
Proof. -With the convention that the exponents are ordered by real part, call z-degree of a logarithmic series the maximum exponent of z that appears, z-valuation the minimum one, and log-degree the largest power of log(z) involved. The application of θ to a term cz ν log(z) k yields a (possibly empty) sum of terms of the same z-degree and at most the same log-degree. Since Q 0 (θ) · q = L(z, θ) ·ũ, whereũ has z-degree less than λ + N , the z-degree of Q 0 (θ) · q is less than λ + N + s. Since we also have Q 0 (θ) · q = −L(z, θ) · (u −ũ), where u −ũ is of z-valuation λ + N or more, Q 0 (θ) · q has z-valuation at least λ + N . The same bounds apply to q(z) due to (6.6) and the requirement that q ν,k = 0 for (ν, k) ∈ E. By a similar reasoning, the log-degree of Q 0 (θ) · q is at most K. By Lemma 6.7, this implies that the q ν,k can only be nonzero for µ(ν) k µ(ν) + K.
When computing the series expansion of the solution u(z) using the recurrence from Proposition 3.1 or its generalization (5.4), the state one needs to maintain from one iteration to the next is a vector (u λ+n−1 , . . . , u λ+n−s ) where, in general, each entry may be a polynomial in log(z). The following algorithm details how to compute the normalized residual at a given iteration from these coefficients only. For simplicity, we limit ourselves to the generic case where neither the truncation index nor its next few shifts by integers are roots of the indicial polynomial. Note that the most expensive part is typically the evaluation at ν + j of the polynomial coefficients of the recurrence and their derivatives (step 2), but these coefficients can be recycled from the iterative computation of the u λ,k . Additionally, fast algorithms dedicated to the evaluation of polynomials at regularly spaced points [Ger04,BS05] are applicable. 1. Write P (S −1 n , n) = b 0 (n) + b 1 (n)S −1 n + · · · + b s (n)S −s n . Let c = p r (0) where p r (z) is the leading coefficient of P (z, θ). Let K = max j max {k : u j,k = 0} (and K = 1 if s = 1).
Proposition 6.10. -In the notation of Algorithm 6.9, let u j,k = u λ+N −j,k be coefficients as in (6.1) of a solution u(z) of P (z, θ). Assume that the indicial polynomial Q 0 of P (z, θ) does not vanish at any of the points λ+N, . . . , λ+N +s−1. Then, Algorithm 6.9 returns the coefficients q j,k = q λ+N +j,k , 0 j < s, k 0, of the normalized residual (6.8) associated to the truncationũ(z) of u(z) defined by (6.2).
Proof. -The algorithm amounts to the computation of the coefficients (v ν,k ) of v(z) = P (z, θ) · u(z) by application of P (S −1 n , ν + S k ) to (u ν,k ), interleaved with the solution of a triangular linear system that expresses the relation Q(ν + S k ) · (q ν,k ) = (v ν,k ). The table v is filled with v j,k = v λ+N +j,k , using the relation with ν = λ + N − j, and noting that u ν−i,· = 0 when j = i − j 0. The factor c accounts for the discrepancy between the monic indicial polynomial Q 0 used in the definition of the normalized residual and the leading coefficient of the recurrence: they are related by b 0 (n) = c Q 0 (n), so that the equation

Bounds on tails of differentially finite series
In terms of L(z, θ) and y(z), the normalized residual we just defined satisfies L(z, θ) · y = Q 0 (θ) · q. Given suitable bounds on the coefficients of L, Proposition 5.5 applies and provides a bound on y(z). If (p,â) is the pair computed by Algorithm 6.1, it follows that the truncation error satisfies whereq(z) is a power series satisfying condition (ii) of Proposition 5.5. The obvious choice is to takeq(z) as a polynomial with the same support as z −λ q(z). However, choosingq(z) as a polynomial multiple ofĥ(z) instead makes the integral explicitly computable. Combined with a parameter choice that render the last two conditions of Proposition 5.5 trivial (cf. the discussion at the end of Section 5.3), this leads to the following algorithm. 1. Rewrite L in the form P (z, θ) with θ on the left: compute coefficients p 0 , . . . , p r in K[z] such that z ρ L = θ r p r (z) + · · · + θp 1 (z) + p 0 where ρ ∈ Z is such that p r (0) = 0. Let s = max k deg p k .
6. Return the symbolic expression In our implementation, the seriesĥ(z) (actuallyĥ(z)/p(z)) andû(z) are represented by objects of type "hyperexponential majorant" that encode formal power series of the form where the polynomial G(z) is represented in factored form and H(z) is an unevaluated sum of rational functions, also with factored denominators. With the algorithm of this paper, the sum reduces to two terms (a polynomial part related toQ(z) in Algorithm 6.1, and a rational part of denominator p r (z)), but variants like that of Remark 6.5 can introduce additional terms. This representation makes it easy to extract numerical bounds from the majorant series, as discussed in Section 8 below.
Proposition 6.12. -Consider the generalized initial values (u ν,k ) (ν,k)∈E associated to u, and assume that N is at least equal to the largest index of a nonzero initial value, that is, Then Algorithm 6.11 returns an expression representing a seriesû(z) Proof. - Step 1 of the algorithm consists in rewriting the equation L · u = 0 into an equivalent one of a special form. We have seen in Section 3.2 that this can always be done, and in Section 5.1 that, due to the assumption that 0 is a regular point, the resulting p r does not vanish at 0. By Lemmas 6.7 and 6.8, the normalized residual q(z) is well defined and the polynomials f [k] have degree bounded by s. Additionally, the rational function returned by Algorithm 6.1 satisfiesâ(0) = 0 (by Proposition 6.3), so thatĥ(z) is analytic at the origin. Thus, at step 5, the computation of g(z) amounts to that of a Taylor expansion ofâ(z) (a rational function with no pole at the origin) followed by routine operations on truncated formal power series. Therefore, the algorithm makes sense and runs without error.
To see that the result satisfies the claim, let us check that Proposition 5.5 is applicable toâ(z) andq(z) = zĝ (z)ĥ(z), whereĝ(z) is the polynomial computed at step 5. By Proposition 6.3, the seriesâ(z) satisfies condition (i). By definition of g(z), we have g (z)ĥ(z) = z N −1f (z) + O(z N +s ). Sincef (z) has degree s, the first N + s coefficients of g (z)ĥ(z) are exactly the coefficients of the polynomial z N −1f (z). While these coefficients are non-negative, this may not be true of the remaining coefficients of g (z)ĥ(z). However, replacing g(z) byĝ(z) yields a productĝ (z)ĥ(z) with and hence z Nf (z) q(z). In combination with the inequalities resulting from steps 2 and 4, this implies and condition (ii) is satisfied. Because N n 0 and y(z) = p r (z)y(z), we have y λ+n,k = 0 for all n < n 0 , hence condition (iii) holds. Similarly, condition (iv) holds due to the assumption (6.9). Therefore, the relation y(z) ŷ(z) holds for any solution of zŷ (z) =â(z)ŷ(z) + q(z), and in particular for It follows that u(z) −ũ(z) = p r (z) −1 y(z) p(z)ŷ(z).
Remark 6.13. -Algorithm 6.11 is but one way to compute remainder bounds based on Proposition 5.5. It admits many variants that use the additional flexibility of the framework of the previous section.
(1) If N is replaced by n 0 in (6.9), the proof of Proposition 6.12 applies verbatim to any N n 0 instead of only N = n 0 . One can thus compute majorants corresponding to multiple truncation orders without running Algorithm 6.1 again, or even specializing again the majorant sequence of Remark 6.4.
(2) It may also happen that (6.9) fails to hold: typically, if the indicial polynomial Q 0 has a root at λ + n 1 for some very large n 1 , one may want to truncate the series expansion of u(z) at an order N n 1 even if one of the generalized initial values at λ + n 1 is nonzero. One can modify Algorithm 6.11 with a more general choice ofŷ(z), as discussed at the end of Section 5.3, so as to cover this case.
(3) There is no need to run the algorithm several times to bound the tails of several solutions of the same equation (corresponding to the same λ) truncated at the same order: it is enough to compute each of the corresponding normalized residuals, and modify step 4 to take them all into account. In particular, using the fact that derivatives of majorants are majorants of derivatives, bounding the remainder of a fundamental matrix of the equation at an ordinary point only requires a single call to Algorithm 6.11. (4) It would be possible to choose the quantity τ (n) needed at step 4 of Algorithm 6.1 in a slightly tighter way, based on (5.5) and the observed degree of the normalized residual after step 2 of Algorithm 6.11. One would then recover the special case of ordinary points hard-coded in Algorithm 6.1. However, proceeding this way complicates the reuse of computations when n 0 varies. (5) Step 5 of Algorithm 6.11 is optional: takingĝ(z) = z 0 w N −1 f (w) dw also yields a valid, if coarser, bound. Indeed, we havef (z) f (z)ĥ(z) sincê h 0 = 1, meaning that we can replacef (z) byf (z)ĥ(z) without contradicting the inequality from step 4, but then the integral at step 5 reduces to

Bounds on rational sequences
The main algorithm presented in the previous section crucially relies on bounds on quantities of the form sup n n 0 |f (n)|, where f (n) is a rational function with complex coefficients. While it is not hard to come up with such bounds (isolating the poles and local extrema of f yields optimal bounds), their computation can easily become costly in practice. This section describes an approach that we found to be a good trade-off between speed and quality for an implementation based on interval arithmetic.

The generic case
When n is large (larger than any of the poles of the rational function f of interest, at least), a simple and effective way to bound |f (n)| for n n 0 is to perform an interval evaluation.
More precisely, if f (n) = p(n)/q(n), write f (x −1 ) = x δp (x)/q(x), wherep,q are the reciprocal polynomials of p, q, and δ 0. The evaluation of this expression on x = [0, 1/n 0 ] in interval arithmetic yields a bound on f (n) valid for all n n 0 . Moreover, this bound converges to the same limit as |f (n)| as n tends to infinity. Note that the rewriting step is essential, as illustrated by the example of f (n) = (n − 1)/n, whose naive interval evaluation on [10, ∞] gives [0, ∞], to be compared with [0.9, 1] after rewriting the expression.
A big advantage of this approach is that it generalizes naturally to the simultaneous computation of bounds on several derivatives of f (n), as required by step 4 of Algorithm 6.1. The generalization is formalized as Algorithm 7.1 below. Given f = p/q ∈ C(n) and T 1, we denote (this is not exactly the same notation as in Algorithm 6.1). This definition already covers the bounds on rational sequences needed in the main algorithm when n 0 is larger than the roots of the indicial polynomial. Algorithm 7.1 also accepts a parameter Z that can make it "ignore" some of the poles of f , and will be useful when we turn to the remaining (non-generic) cases. In this algorithm and the next one, boldface letters stand for intervals or polynomials with real or complex interval coefficients. All interval operations are extended in a natural way to handle intervals containing ∞. An instruction like "compute g(ε) = ϕ(ε) + O(ε T )" means "compute a polynomial g of degree at most T − 1 whose coefficients are interval enclosures of the first T Taylor coefficients of ϕ(ε)". The computation reduces to routine operations on truncated power series.
Before proving that this algorithm works as stated, let us discuss the step marked as optional. Without this step, the algorithm is a direct generalization of the method for T = 1 sketched above. It returns an infinite bound as soon as q has a real root in [n 0 , ∞), and only gives satisfying results when n 0 is sufficiently larger than the largest real root. Note that this may be enough in the context of solutions of differential equations at ordinary points, since the only denominator that arises is then n(n − 1) · · · (n − r + 1), where r is the order of the equation. In the general case, the next lemma offers a simple way to mitigate the issue.
Lemma 7.2. -Let q ∈ C[n] be a monic polynomial of degree d, and let n 0 1.
Then, for all n n 0 , we have where the product is over the multi-set of roots α of q with Re α > 0.
Note that this lower bound can be somewhat expensive to compute compared to the rest of the algorithm. For this reason, our implementation actually decides whether to run the optional step based on the accuracy of the intervalq(0) at the beginning of step 5.
In practice, there is no need for the check that Re α > 0 to be exact, since false positives can only decrease the result. Also, rough enclosures of the roots of q are sufficient: one can add more terms to the minimum in the definition of b α (n) if the enclosure of α does not uniquely determine π(α). Similarly, given Z ⊂ N, one obtains a lower bound valid for n n 0 with n ∈ Z by replacing π(α) by the adjacent integers when α ∈ Z.
Proposition 7.3. -Algorithm 7.4 returns a quantity M ∈ R + ∪ {∞} such that F [T ] (p/q, n) M for all n n 0 with n ∈ Z. The version that includes the optional step returns a finite bound as soon as Z contains all the roots of q in N n 0 (provided that the working precision for interval operations is large enough).
Proof. -When n 0 = 0, the algorithm returns ∞. Assume n 0 1. Fix n n 0 such that q(n) = 0, and let x = n −1 . The quantitiesp andq defined in the algorithm are polynomials in ε (since d = deg q > deg p) with complex coefficients. Letting they satisfy .
Since q(n) = 0, one can compute the first T Taylor coefficients of nf (n + ε) by evaluating the right-hand side of this expression in C[[ε]] and truncating the intermediate results to the order T after each operation. Without step 5, this is exactly what the algorithm does (in interval arithmetic) to compute s(ε). Since the interval x contains n −1 , we have n [ε t ]f (n + ε) ∈ s t for 0 t < T . Now consider the case where we include step 5. Suppose additionally that n ∈ Z, and let c = n dq (n −1 ) = n −d q(n). Note that c = 0. After step 5, the modified polynomialq(ε) satisfies Indeed, the relation for t = 0 reduces to 1 = 1, while for t > 0, it follows from the fact that |c| ρ and hence c −1 ∈ γ. Therefore, we have c n [ε t ]f (n + ε) ∈ s t after step 6.
In both cases, we conclude that n [ε t ]f (n + ε) ∈ γ (|s 0 | + |s 1 | + · · · + |s T −1 |) for all n ∈ N n 0 \Z such that q(n) = 0. If q(n) vanishes for some n ∈ N n 0 \Z, then the algorithm returns infinity (either because the constant coefficient ofq contains zero or because γ is unbounded), hence the bound holds for all n ∈ N n 0 \Z. Otherwise, we can take ρ > 0, and the result is finite.

The general case
In its general form, Algorithm 6.1 requires bounds for n n 0 on sequences similar to F [T ] (f, n), but with a number of summands T that varies with n and a special handling of poles of the denominator. Both modifications to the formula can be summed up by introducing as an additional parameter the sequence of multiplicities of the poles that we wish to treat in a special way. Thus, given a rational function f = p/q as before and a sequence m(n) with finitely many nonzero terms, let us define where τ (n) = n k=0 m(k). The quantities that we have to bound at step 4 of Algorithm 6.1 are of this form (3) . The following algorithm computes the associated bounds.
Algorithm 7.4. -Bound rational sequence, general case Input: Polynomials p and q, a set Z, and an integer n 0 ∈ N\Z as in Algorithm 7.1.
A "multiplicity" function m : N → N with m(n) = 0 for n ∈ Z. Output: A bound M ∈ R + ∪ {∞} such that the function F : N → R + ∪ {∞} defined by (7.2) satisfies F (n) M for all n n 0 .

1.3.
If b is larger than the maximum of the S(k) defined so far, set S(n) = b.
2. Find the smallest n n 0 on which S is defined. Set M exn = S(n).

Return max(M exn , M gen ).
(3) There is no harm in also replacing µ( · ) by zero in (6.3) when we take τ (n) = 1 because the origin is an ordinary point: doing so only makes the first few terms of the bound infinite.
Remark 7.5. -Algorithms 7.1 and 7.4 immediately extend to vectors of rational functions with the same denominator. A large part of the computation can be shared between the entries.
Proposition 7.6. -Given parameters that define a sequence F (n) of the form (7.2) and an integer n 0 , Algorithm 7.4 returns a bound M such that sup n n 0 If, for every integer n n 0 , the multiplicity of n as a zero of q is at most m(n) (and the working precision is large enough), then M is finite.
Proof. -First observe that, due to step 5.3, the values S(n) defined at step 1 are decreasing (in the sense that S(m) < S(n) when both are defined and m > n). We extend S to a staircase function defined on N ∪ {∞} by setting the undefined S(n) to the smallest values that make S non-increasing (see Figure 7.1).
Let us prove by induction that This is true for n > max Z. Then fix n ∈ Z and assume (7.3) holds for larger indices.
Step 5.1 of the algorithm ensures that S(n) F (n). For n < k < min(N >n ∩ Z), we have F (k) = F [τ (n)] (p/q, k), hence S(n) F (k) by the correction of Algorithm 7.1 (Proposition 7.3). Finally, by the induction hypothesis, we also have S(n) F (k) for larger k.
When n 0 ∈ Z, the inequality sup n n 0 F (n) M holds by (7.3). When n 0 ∈ Z, as with steps 15.2-15.3, we have M S(n 0 ) S(n 0 ) S(n) for n n 0 , where n 0 = min(N n 0 ∩ Z), and M F [τ (n 0 )] (p/q, n) for n 0 n < n 0 . It follows that M F (n) for all n n 0 .

Numerical bounds
The algorithms developed at this stage bound the tails of differentially finite series by hyperexponential majorant series. We have seen that a seriesû(z) with u(z) û(z) encodes bounds |u (j) (ζ)| û (j) (|ζ|) on the values of u and all its derivatives. Yet, in order to use the majorants in a concrete setting, we still have to explain how to effectively evaluate theû (j) (|ζ|). Let us first consider the details of this operation, and then outline two other ways-appropriate for different settings-of deriving numerical tail bounds from the majorant series.

Values
When j is large, it is essential for performance to compute the values of the derivatives by working with truncated power series (or another similar compact representation) rather than by naive symbolic differentiation. The process is laid out in Algorithm 8.1 below, which takes as input a hyperexponential majorant represented in the form discussed after Algorithm 6.11. Despite the heavy notation, the algorithm is mostly straightforward, the only subtlety being that we bound sub-expressions of the form (f /g) by ( f ) /g to avoid computing integrals of rational functions. As in the previous section, boldface letters stand for intervals or polynomials with interval coefficients, and "compute ϕ(ε) + O(ε m )" means "compute a Taylor expansion with interval coefficients of ϕ(ε), truncated to the order m, using power series operations".
Proposition 8.2. -In Algorithm 8.1, suppose that for all i, the seriesû i (z) anď v i (z) −1 have non-negative coefficients. Then the values M j returned by the algorithm are bounds for the derivatives ofû(z) at x, with 0 û (j) (x) M j , 0 j < m.
Proof. -Steps 1 to 3 compute the Taylor expansion to the order m at z = x of

An integration by parts shows that
where the left-hand side and both terms on the right-hand side are series with non-negative coefficients, so that we have

ANNALES HENRI LEBESGUE
It follows that the coefficients of S(ε) are enclosures of upper bounds on the corresponding Taylor coefficients ofû(|ζ| + ε).
When u(z) is a plain power series and u n: (z) û(z), the above algorithm yields bounds on |(u n: ) (j) (ζ)|, that is, on values of the successive derivatives of the tail. This is what one typically needs in applications, yet, since (û n: ) (z) = Nû n z n−1 + (û ) n: (z), the same quantities also bound the tails of the derivatives.
If u(z) = z λ k f k (z) log(z) k /k! is a logarithmic series, though, we only get bounds on derivatives of tails of the components f k (z). To deduce bounds that apply to tails of the generalized series expansion of u(z) itself, it is usually best to form the expansion of f k (|ζ| + ε) with respect to ε and compute (|ζ| + ε) λ k f k (|ζ| + ε) log(|ζ| + ε) k /k! in power series arithmetic, similar to what Algorithm 8.1 does.

"A priori" bounds
A small drawback of the technique discussed above is that it requires knowing the last few coefficients u n−1 , . . . , u n−s just before the truncation point in order to bound the tail u n: (ζ) of a series u(z). Thus, the results are in a sense a posteriori bounds, that cannot be used to decide with certainty where to truncate the series before even starting the computation of the coefficients. (It is of course possible to "guess" a plausible truncation order based on asymptotic considerations, or to compute the required coefficients at low precision before running the full-precision computation.) If however, given a majorantû(z) of u n 0 : obtained by the previous methods for some n 0 , we can get sufficiently tight bounds on the values of the remaindersû n 1 : (|ζ|) for n 1 > n 0 , then we have a way of bounding the higher-order remainders u n 1 : (ζ) without computing all the coefficients up to n 1 . We now adapt to the setting of this paper a technique of Mezzarobba and Salvy [MS10] for doing so, based on the classical saddle point method in asymptotics. The idea is that for suitably chosen ρ, the bound (8.1) below is relatively tight.
Lemma 8.3. -Letv(z) be a power series with non-negative coefficients. Fix real numbers ρ x > 0 within the disk of convergence ofv. For all n ∈ N, the series expansion at x of the remainder of order n ofv(z) is bounded as In particular, we havev n: Since ρ x, we havev n: (xz) x ρ nv n: (ρz) x ρ nv (ρz).
The result follows by substituting 1 + x −1 ε for z.
For simplicity, let us focus on seriesû(z) of the shape returned by Algorithm 6.11, viz., and assume thatâ(z) andb(z) are not both constant (4) . Additionally, we limit ourselves in the analysis to the bound (x/ρ) nv (ρ) on the value of the remainder, leaving to the reader the case of derivatives.
Let ρ * be the pole ofâ(z)b(z) closest to the origin, with ρ * = ∞ if bothâ(z) andb(z) are polynomials. Asâ(z) andb(z) have non-negative coefficients, ρ * is real (or infinite) and positive. Besides, in our setting, the denominators ofâ(z) and b(z) divide the polynomialp(z) computed by Algorithm 6.1, hence ρ * can be taken arbitrarily close to the modulus of the singularity of the differential equation closest to the origin. The next lemma shows that, in (8.1), it is possible to choose ρ = ρ n close to ρ * , so that (x/ρ n ) n decreases fast, without lettingv(ρ n ) grow too large. In particular, the resulting sequence of bounds tends to zero at least exponentially fast, with the same exponential rate x/ρ * as the coefficients ofv n in the case ρ * < ∞.
In both cases, the other factor ofv(ρ n ) satisfiesb(ρ n ) = n O(1) . Finally, the prefactor (x/ρ n ) n is bounded as The result follows by combining these estimates. While, asymptotically, these formulas yield tight bounds for any fixed c, the actual values of (x/ρ n ) nv (ρ n ) for finite n are quite sensitive to that of ρ n , and the bounds implied by the previous lemma are typically very poor for moderate n. We could, of course, derive more precise asymptotic estimates of the optimal choice of ρ n as a function of n, but in practice it is best to minimize F n (ρ) = (x/ρ) nv (ρ) numerically, not necessarily in a rigorous way. Note that, for large n, the function F n has a unique local minimum (5) on [0, ρ * ), so that simple numerical methods work well for this purpose.
To sum up, we can use the following algorithm to compute bounds onû n: (x) and its first m derivatives given a seriesû(z) of the form (8.2): first call Algorithm 8.1 onû(z) and x. If n n 0 , return its result. Otherwise, define ρ * as above, and search for a value ρ ∈ [x, ρ * ) that minimizes log((x/ρ) n f (ρ)), where f (ρ) is the value computed by Algorithm 8.1 called with m = 1.
Compare with the bounds computed at the first step, and return the better one. (This last step is not strictly necessary, since (8.1) reduces to the trivial boundv n: (x + ε) v(x + ε) when ρ = x, but it can help if the numerical computation of ρ is inaccurate.) Proposition 8.5. -If u(z) û(z) and n n 0 , the algorithm outlined above yields valid bounds on u n: (ζ), . . . , (u n: ) (m−1) (ζ) for all ζ such that |ζ| x. Provided that the numerical method used to choose ρ is accurate enough, the bound on u n: (ζ) tends to zero at least exponentially fast as n grows. Figure 8.1 shows the results of a simple implementation of this strategy. While the sub-exponential overhead promised by Lemma 8.4 can be observed in practice, at least for the simpler of the equations, we also see that these "a priori" bounds are much more pessimistic than the "a posteriori" ones in slightly more complicated cases.

Remainders of Laplace transforms
Another downside of our framework is that hyperexponential majorant series are not expressive enough to capture the remainder asymptotics of solutions of arbitrary differential equations with polynomial coefficients. For example, the Airy function (5) Indeed, its logarithmic derivative reads F n (ρ)/F n (ρ) = ρ −1 (G(ρ)−n) where G(ρ) = ρb (ρ)/b(ρ)+ a(ρ). Since the numerator ofb(z) has non-negative coefficients, G(ρ) is a rational function without poles on [0, ρ * ), and hence a continuous function with finitely many local extrema. Additionally, it satisfies G(0) = 0 and G(ρ) → +∞ as ρ → ρ * . Therefore, for large n, the equation G(ρ) = n has exactly one solution. satisfies Ai(z) = exp(σ ϕ z 3/2 + O(log z)) as z → ∞ in a generic direction ϕ, and [z n ] Ai(z) = O(n! −2/3 ), but the best majorants our method can express are of the formû(z) =ĝ(z) exp(σz 3 ),ĝ(z) ∈ Q[z], with u n of the order of n! −1/3 . More generally, for any rational κ 0, it is not hard to find an equation whose solutions are entire functions u(z) of order κ or higher, and hence have tails u n: (ζ) bounded by n! −1/κ e O(n) (ζ fixed). In contrast, hyperexponential entire functions are necessarily of integer order, with tails that decrease like n! 1/d e O(n) where d ∈ N >0 . Besides, even when the fastest-growing solution of the differential equation we are considering is of integer order, there is no guarantee that the hyperexponential majorant will respect that order.
This overestimation is not too much of an issue for the "a posteriori" error bounds discussed in the main part of this paper because, even if the choice ofĥ(z) in Section 6.4 is not as tight as it could, the factor derived from the residual has the "correct" asymptotics with respect to n. For the "a priori" bounds of the previous subsection this argument does not apply, and they can be far from matching the actual speed of convergence of the solutions, even in non-degenerate cases.
A possible remedy (essentially a streamlined version of a similar idea proposed by Mezzarobba and Salvy [MS10], to which we refer for more information) is to re-scale the solutions of the original equation by a generalized Laplace transform Lemma 8.6. -Letû(z) = n 0ûn z n ∈ R + [[z]], and letv(z) = n 0 φ(n)û n z n for integers p ∈ Z, q ∈ N >0 chosen so thatv(z) has a nonzero radius of convergence. Let ρ > 0 lie within the disk of convergence ofv(z). Then, for all x and n such that 0 x ρ (n + 1) p/q , we have the bound Proof. -Writê Since the function t → Γ(t+a)/Γ(t) is non-decreasing for fixed a, and using Gautschi's inequality [OLBC10a, Section 5.6.4], we have Γ(1 + n/q) Γ(1 + (n + k)/q) The result follows. In practice, when applying this lemma, x is given, n must be large enough that x ρ * (n + 1) p/q holds, and ρ is to be chosen as a function of n as in the previous subsection. For large x, Lemma 8.6 starts being applicable when n ≈ x p/q , which typically matches the position where the terms of the series u(z) "start converging" (cf. Section 9.2).

Implementation
We have implemented the algorithms of this paper in ore_algebra [KJJ15, Mez16], a library for working with Ore polynomials in the SageMath (Sage) computer algebra system. Among the features of Sage and the many external libraries on which it is based, our code relies in an essential way on Arb [Joh17], which provides all necessary basic operations on intervals and truncated power series with interval coefficients. The implementation also makes use of PARI for complex root finding, and of ore_algebra itself for basic arithmetic with differential and recurrence operators. The operations on rational numbers and exact polynomials mainly come from GMP/MPIR, Flint and Singular.
The ore_algebra package is available from https://github.com/mkauers/ore_ algebra/ under the GNU General Public License (version 2 or later). The code for computing error bounds can be found in the file src/ore_algebra/analytic/ bounds.py of the source tree. We performed the experiments described below using the git revision dd88e8d5 under Sage 8.2.
Except when otherwise noted, we run the full version of the algorithm described in this paper, that is, the combination of Algorithms 6.1, 6.11, 7.4, and 8.1, including their optional steps, and using PARI's root finder at step 2 of Algorithm 6.1. Numerical coefficients are represented by real or complex intervals almost everywhere, with a fixed precision of 53 bits.
In examples involving generalized series expansions at singularities, the displayed truncation error and error bound correspond to the complete logarithmic series, including the non-analytic factors.
Remark 9.1. -Testing the implementation for correctness is a significant issue, as the final truncation bounds are likely to be valid even if the code is wrong, due to overestimation. Without a complete formal proof of the implementation, it is very hard to be certain that it is fully correct. Nevertheless, the fact that our bounds are close to optimal in simple cases is of great help in catching bugs. We further limit the risks of missing issues hidden by overestimations by testing not only the final bounds, but also various intermediate results. In our experience, plausible but incorrect changes to the bound computation algorithm tend to be caught by the test suite. (z 2 + 1) u (z) + 2z u (z) = 0, u (z) + 2z u (z) = 0, u (z) − z u(z) = 0, 4z 2 u (z) − (z 2 − 8z + 11) u(z) = 0. For each of these equations, we consider the series expansion at the origin of one or several solutions (respectively the exponential, the arctangent, the error function, the Airy functions Ai and Bi, and certain linear combination of Whittaker functions of parameters κ = 2, µ = √ 3), evaluated at different points ζ of their disk of convergence. The first four equations are ordinary at the origin. The last one is regular singular, with irrational exponents 1/2 ± √ 3. The local expansion of M κ,µ lies in z 1/2− √ 3 C[[z]], and the constant (6) c is chosen so that W κ,µ − cM κ,µ ∈ z 1/2+

Elementary and special functions
We plot the truncation error |u n: (ζ)| and the corresponding bound on a logarithmic scale, as a function of the number of terms n, until they become smaller than 10 −10 . All bounds are computed using Algorithm 6.1 with = 3. The left column corresponds to evaluations at points ζ where the series u(ζ) converges nicely-halfway (6) c = π/(sin(απ)Γ(α)Γ(− √ 3 − 3/2)) where α = 1 + 2 √ 3. from the circle of convergence in the case of arctan(z), and at points where the magnitude of the terms of u(ζ) start decreasing early on in the case of entire functions. For such simple equations, the experiments confirm that the bounds are very tight, in spite of the very low accuracy target and the moderate value of . Though a bit larger, the overestimation remains moderate in the two other columns. We will come back to the case of arctan(7/8) later. We can observe that when the terms exhibit a "hump" before the series really starts converging, our bounds tend to overestimate the true error by a factor comparable to the height of the hump. A heuristic explanation is that, in the majorant seriesû(z) =p(z) −1ĝ (z)ĥ(z) that we compute for u n: (z), the coefficients ofĝ(z) are roughly of the order of |u n |, whileĥ(z) is not too far from j 0 |u j | z j . Thus, the value of the majorant series at x = |ζ| > 0 is about j 0 |u j | x j ≈ max j 0 (|u j | x j ) times larger than |u n: (ζ)| ≈ |u n | x n . This overestimation has a limited impact in applications, because one typically tries to avoid humps in the coefficients in the first place-for example, in a Taylor method, by adjusting the step size. Nevertheless, it would be interesting to find a way of avoiding it.

"Real-world" examples
Next, we consider a family of larger examples borrowed from an application. Other families of "interesting" operators, as well as tools to produce plots similar to the ones in this paper, are available in the ore_algebra source tree. Figure 9.2 shows remainder bounds for solutions at the origin of the differential equations for lattice Green functions of face-centered hypercubic lattices discovered (heuristically or rigorously) by several authors over the last decade [Gut09,Bro09,Kou13,ZHM15,HKMZ16]. In particular, the case d = 4 corresponds to the operator of Example 1.1, but without the shift from the neighborhood of 0 to that of 1/2. The other operators can be found on Ch. Koutschan's web page (7) ; here we limit ourselves to collecting some statistics about their size: The value at z = 1 of the lattice Green function gives access to the return probability of the random walk on the lattice, and, to the best of our knowledge, going through the differential equation is the only known way of computing that probability to high precision in good complexity.
All these operators have regular singularities with integer exponents at the origin. For example, the monic indicial polynomial for d = 9 is n 9 (n − 1) 5 (n − 2) 3 (n − 3). In each case, we select small initial values at random and evaluate the logarithmic series solution characterized by these initial values about halfway from the boundary of its circle of convergence. We run Algorithm 6.1 with = s/2 + 2, in accordance with the heuristic suggested on page 117. It is especially noticeable here that the left endpoint of the upper curve is located significantly to the right of that of the other curve, even when no part of the curve is clipped out of the plot. This is because our implementation can return infinite bounds, and does so at least until it has passed all initial values.
The fast convergence of the power series as d increases can be explained by the presence of apparent singularities close to the origin, which leads us to select evaluation points much smaller than the actual radius of convergence of the solutions. Thus, roughly speaking, a rigorous order-adaptive Taylor method using our algorithm would perform many more steps to evaluate the lattice Green function at z = 1 for d = 9 than for d = 4, but sum fewer terms at each step. Although our bounds tend to infinity when the magnitude of the evaluation point approaches that of the singular point of the equation closest to the origin-apparent or not-, we can see that in these examples, they adapt nicely to the unexpectedly fast convergence, thanks to the use of residuals. Having to perform many small steps is nevertheless sub-optimal, which raises the question of adapting the algorithm to handle apparent singularities in a way that allows for larger step sizes.
We can also observe that overestimation becomes significant for the larger of these operators. Further increasing the tuning parameter does not seem to help. Other methods may hence be needed for the really efficient numerical solution of equations with coefficients of very large degree at comparatively low precision.
It would be hard to draw any conclusion from a detailed analysis of the running time of our Python implementation. Still, anecdotal data confirms that the bounds of this paper are usable in practice as part of symbolic-numeric algorithms. For example, the simple variable-order interval Taylor method implemented in ore_algebra can compute an enclosure with radius less than 10 −60 of the value at z = 1 of the lattice Green function for d = 4 in less than a second on a modest laptop. The computation of the remainder bounds accounts for an estimated 5 to 10% of the total running time. Of course, the solver does not compute a complete tail bound at every iteration, but starts with a heuristic convergence check, and reuses various intermediate results when the rigorous check does not succeed on the first try, as discussed in particular in Remark 6.4. The evaluation point is a simple rational halfway from the circle of convergence, and we select the initial values at random. (d) An operator of order 7 and degree 8 of the form ((20i − 1)z 8 + · · · + 178371756i + 577700775)D 7 z + · · · + ((1 − i)z 8 + · · · − 42079467i − 23600391) obtained in a similar way, also evaluated at a point halfway from the circle of convergence and with random initial values. We can observe that going from = 1 to = 2 already brings very significant improvements, and that larger values of are useful for nontrivial examples. Clearly, the heuristic ≈ s/2 suggested earlier is but a crude rule of thumb. An implementation that adaptively increases starting from = 2 based on a heuristic measure of tightness is likely to perform better. These conclusions are consistent with experiments on various other examples. We also see again that operators with a few dozens of monomials and large (but not huge) coefficients can be challenging in terms of overestimation, at least at low precision.

Effort-tightness trade-off
A second parameter involved in the trade-off between speed and accuracy is the optional step of Algorithm 7.1, whose effect already was illustrated on Figure 7.1. This optional step can bring huge benefits on some inputs but makes no difference in the most common cases, while being relatively expensive. Hence, it makes sense to only run it when the previous step does not suffice.
to save time by doing with a single shared lower bound on the singularities, as suggested in Remark 6.2, only makes sense for very simple operators.

Comparison with existing software
A systematic comparison with the related methods discussed in Section 2 would be meaningless, as they were designed for different goals, widely differ in scope, often (even more than ours) involve tuning parameters with a crucial influence on the quality of the results, and usually are not implemented. Nevertheless, it makes sense to check that, in simple cases where several methods apply, our bounds are not worse than what could be computed by other approaches.
We limit ourselves to a very simple example borrowed from Neher [Neh03, Example 3], namely the problem of bounding the remainders of the Taylor expansion of u(z) = cos(z)/(z 2 + 101), solution of (9.1) ((z 2 + 101)D 2 z + 4zD z + (z 2 + 103)) · u = 0, u(0) = 1/101, u (0) = 0. Table 9.1 lists bounds on u n: (ζ) for several values of ζ and n computed by the following methods: NumGfun: The author's own Maple package NumGfun [Mez10], whose limitations prompted the present work. We use the function bound_diffeq_tail() with no particular options or configuration, and substitute the appropriate ζ and n in the result. vdH2001: A method described by van der Hoeven [vdH01, Section 2], here applied manually with the help of non-rigorous floating-point computations. While the description leaves out a number of algorithmic details related to the computation of bounds on rational functions, the optimal choice of parameters (λ = 101 −1/2 , M 0 = 103/101, M 1 = 0) is unambiguous here thanks to the special shape of our example. Like the "a priori" bounds of Section 8.2 and for similar reasons, this method is highly sensitive to the choice of the remaining parameter µ, which we select by numerical optimization to minimize the bound. The six test cases of Table 9.1 respectively use µ = 0.139, µ = 0.121, µ = 0.139, µ = 0.121, µ = 0.105, and µ = 0.105.

vdH2003:
A related method also suggested by van der Hoeven [vdH03, Section 3.5], no longer involving µ. Again, it is possible for our particular example to choose optimal parameter values α = 101 −1/2 , M = 103/101, C = 1/101. We assume that the remainder of the majorant series obtained as output can be bounded as tightly as needed. ACETAF: We also list bounds obtained by Neher [Neh03] using version 2.71 of ACETAF [EN03]. ACETAF is mainly designed for computing bounds on the coefficients and remainders of series coefficients of linear analytic ODEs before applying the method of majorants. Therefore, unlike the other techniques we are comparing, it does not start from the equation (9.1) but from the closedform expression of u(z). Both the quality of the bounds it returns and the running time vary widely depending on the values of several tuning parameters. We only list the best of the three bounds reported in Neher's article for each instance. ore_algebra: Finally, we run the implementation in ore_algebra of our algorithm, with = 2, using the simplified method of Remark 6.2 for bounding the denominator. Only ACETAF approaches the accuracy of the bounds of the present paper in some cases. In addition, the alternative methods considered here tend to degrade rapidly, in quality, running time, or both, as the size of the operator increases, so that the gap would likely be wider for realistic problems.