Hoeven Ball

Hal-00432152v3

Joris Van Der Hoeven

Ball arithmetic *

Keywords: Ball arithmetic, interval arithmetic, reliable computing, computable analysis A.M.S. subject classification: 65G20, 03F60, 65F99, 37-04

Introduction

Motivation

Computer algebra systems are widely used today in order to perform mathematically correct computations with objects of algebraic or combinatorial nature. It is tempting to develop a similar system for analysis.

Typical problems of this nature are the numerical resolution of a system of equations or the integration of a dynamical system with known initial conditions. In a "computer analysis" system, it should be possible to solve such problems in an exact way. In particular, it should be possible to compute arbitrarily close approximations for all numbers occurring in the input and the output. Computable analysis [Wei00, BB85, Abe80, Tur36] provides a suitable logical framework for computing with this kind of "computable numbers".

Assume that we require an approximation y ˜∈ Q of a number y = f (x) ∈ R with |y ˜y | ε, where ε ∈ Q > = {a ∈ Q: a > 0} is provided by the user. A typical strategy in order to compute such an approximation is to evaluate y p ≈ f (x) numerically at a fixed bit precision p, while computing a rigorous upper bound δ p for the error |y pf (x)|. As soon as δ p ε, then we are done. Otherwise, we keep repeating the same procedure at higher and higher precisions. Interval analysis constitutes a classical theory for the computation of such rigorous upper bounds [Moo66, AH83, Neu90, JKDW01, Kul08, MKC09, Rum10]. The website [Kreb] provides a wealth of information about the interval analysis community.

The preceding discussion imposes some minimal requirements on the system we are striving for: it should contain at least three layers for multiple precision arithmetic, interval analysis and computable analysis. It is natural to add another layer, or at least an interface, for classical numerical analysis. Taken on its own, each of these layers is well understood. However, it is also well known that multiple precision operations are at least one order of magnitude slower than hardware operations on single or double precision numbers. More generally, the naive combination of known techniques leads to highly inefficient algorithms, as we will see in section 5 on the hand of a problem as simple as matrix multiplication.

In order to design a reasonably efficient system for computer analysis, we therefore need to rethink many of the algorithms from the individual layers in a more integrated way. Moreover, in a reasonably complete system, this rethinking should be performed in a very systematic way for the most basic mathematical objects: besides operations on numbers, we need efficient, numerically stable and reliable algorithms for standard operations on matrices, polynomials, series and analytic functions.

In this paper, we will focus on the single layer of interval analysis (with the short exception of section 5). The first part of the paper provides a survey of well known techniques, but from our perspective of usefulness for the design of a computer analysis system. In particular, it turns out that ball enclosures are usually more convenient than interval enclosures in our context. This is partly due to the fact that multiple precision enclosures are usually very narrow. In the second part of the paper, we provide a more detailed study of basic operations on non scalar objects, such as matrices and polynomials. To the best of our knowledge, several of the algorithms in the second part are new.

Although ball arithmetic is often referred to as "midpoint-radius interval analysis", we prefer our naming for three main reasons. First, since we will be using balls all around, it is nice to have a short name. Secondly, complex balls are not intervals. More generally, there is a considerable latitude concerning the choice of a center type and a radius type. Finally, balls are standard in mathematics, where they correspond to classical δ-ε calculus.

Existing software

There are several specialized systems and libraries which contain work into the direction of system for computer analysis. A fairly complete list of software for interval computations is maintained at [Krea]. For instance, Taylor models have been used with success for the validated long term integration of dynamical systems [START_REF] Berz | Cosy infinity version 8 reference manual[END_REF][START_REF] Makino | Remainder differential algebras and their applications[END_REF][START_REF] Makino | Suppression of the wrapping effect by Taylor model-based validated integrators[END_REF]. A fairly complete Matlab interval arithmetic library for linear algebra and polynomial computations is Intlab [START_REF] Rump | INTLAB -INTerval LABoratory[END_REF][START_REF] Rump | Fast and parallel interval arithmetic[END_REF]. Another historical interval library, which continues to be developed is [ea67].

There exist several libraries for multiple precision arithmetic with correct or at least well specified rounding [START_REF] Hanrot | MPFR, a C library for multipleprecision floating-point computations with exact rounding[END_REF][START_REF] Haible | CLN,a Class Library for Numbers[END_REF], as well as a library for multiple precision interval arithmetic [START_REF] Revol | MPFI, a multiple precision interval arithmetic library[END_REF], and libraries for computable real numbers [M0, [START_REF] Lambov | Reallib: An efficient implementation of exact real arithmetic[END_REF][START_REF] Balluchi | Ariadne: a framework for reachability analysis of hybrid automata[END_REF]. There are also libraries for very specific problems. For instance, the MPSolve library [START_REF] Bini | Design, analysis, and implementation of a multiprecision polynomial rootfinder[END_REF] allows for the certified computation of the roots of multiple precision polynomials. Similarly, PHCpack [START_REF] Verschelde | PHCpack: A general-purpose solver for polynomial systems by homotopy continuation[END_REF] and Bertini [START_REF] Bates | Bertini: Software for numerical algebraic geometry[END_REF] are systems for polynomial system solving via numerical homotopy methods.

Unfortunately, the mere agglomeration of existing software is insufficient for building an efficient general purpose system for computer analysis. Besides licensing issues and technical interfacing problems, one of the main reasons is that the reliable computation community has a traditional focus on limited precision computations. Even though this situation is evolving, we have already mentioned the fact that the mere replacement of machine arithmetic by arbitrary precision arithmetic induces a huge performance penalty. Software in which this efficiency problem is addressed (libraries for computable real numbers, MPSolve, etc.) usually focusses on numbers or a very specific problem.

Mmx] use "analyziz"

Mmx] bit_precision := 128; Mmx] z == series (ball 0.0, ball 1.0); Mmx] B == exp (exp z -1) 1.0000000000000000000000000000000000000 + 1.00000000000000000000000000000000\ 0000 z + 1.000000000000000000000000000000000000 z 2 + 0.8333333333333333333333333333\ 3333333 z 3 + 0.625000000000000000000000000000000000 z 4 + 0.433333333333333333333333\ 33333333333 z 5 + 0.28194444444444444444444444444444444 z 6 + 0.174007936507936507936\ 507936507936508 z 7 + 0.102678571428571428571428571428571429 z 8 + 0. These considerations have motivated the development of Mathemagix [vdHLM+02], a new system for computer analysis, which is distributed under a free licence and does not depending on proprietary software such as Matlab. On the one hand, Mathemagix provides a wide variety of numerical types, such as machine doubles, multiple precision numbers, complex numbers, tangent numbers, intervals, balls, computable numbers, but also integers, rationals, algebraic numbers, modular integers, etc. On the other hand, we implemented standard mathematical containers, such as matrices, polynomials and power series. We systematically attempt to make these implementations as efficient as possible for the various kinds of numerical types. For instance, the algorithm for matrix multiplication will be very different for, say, rational matrices and complex ball matrices.

The Mathemagix system has been designed in four layers, with multiple precision arithmetic and classical numeric methods at the bottom, and the certification of numerical results and computable analysis at the top. Although we briefly discuss this general architecture in section 5, this paper is mainly devoted to the automatic computation of error bounds via ball arithmetic. We will outline both the mathematical and implementation techniques which are used in order to make this certification layer as efficient as possible. With the exception of some ideas in section 7.5, all algorithms in this paper have been implemented in our system.

Although we will focus on ball arithmetic, it should be noticed that high efficiency of the bottom layers is crucial in order to achieve high performance. Actually, one important design principle is to delay certification methods to the very end, so that the bulk of the execution time is spent on the actual numerical computations and not on their certification. Now there are many books on numerical analysis at limited precision [START_REF] Press | Numerical recipes, the art of scientific computing[END_REF][START_REF] Stoer | Introduction to Numerical Analysis[END_REF]. The systematic development of multiple precision algorithms which are both efficient and numerically stable is still a major challenge. Since switching from machine precision to multiple precision involves a huge overhead, another challenge is to rely as much as possible on algorithms which first compute an approximation at low precision and then adjust the result at high precision. Inside Mathemagix, we are helped by the fact that most current packages are C++ template libraries, which can be instantiated for numbers at single, double, quadruple and multiple precision (see [START_REF] Bailey | QD, a C/C++/Fortran library for double-double and quad-double arithmetic[END_REF][START_REF] Bailey | Algorithms for quad-double precision floating point arithmetic[END_REF] for a double-double and quadruple-double precision library).

In figure 1, we have illustrated the new algorithms available in Mathemagix on a simple example from combinatorics: the reliable approximation of Bell numbers B n in limited precision (128 bits in our case). Although one might use an asymptotic for-mula for this particularly well known example, we have preferred to compute B n from the exponential generating series B(z) = n∈N B n z n /n! = e e z -1 . Indeed, this very general technique can also be applied in many other cases and has applications to the random generation of combinatorial objects [START_REF] Flajolet | A calculus for the random generation of labelled combinatorial structures[END_REF]. In figure 1, each floating point number in the output is certified to be correct modulo an error of the order of its last digit.

Overview

Ball arithmetic, also known as midpoint-radius interval arithmetic, provides a systematic tool for the automatic computation of error bounds. In section 3, we provide a short introduction to this topic and describe its relation to computable analysis. In particular, we recall how to compute with approximate real and complex numbers. Although ball arithmetic is really a variant of interval arithmetic, we will give several arguments in section 4 why we actually prefer balls over intervals for most of our applications. Sections 2, 3 and 4 mainly provide a survey of well known results in the reliable computing community. Only our perspective might be slightly less classical. Some of the considerations in section 3.5 also seem to be less well known.

In sections 6 and 7, we turn our attention to more complicated mathematical objects, such as matrices, polynomials, power series and analytic functions. Although some of the techniques are still well known (such as Hansen's method for the inversion of matrices), many of the algorithms presented in these sections seem to be new or expand on earlier work by the same author [vdH07b]. Usually, there is a trade-off between efficiency and the quality of the obtained error bounds. One may often start with a very efficient algorithm which only computes rough bounds, or bounds which are only good in favourable wellconditioned cases. If the obtained bounds are not good enough, then we may switch to a more expensive and higher quality algorithm.

The main algorithmic challenge in the area of ball arithmetic is to reduce the overhead of the bound computations as much as possible with respect to the principal numeric computation. In favourable cases, this overhead indeed becomes negligible. For instance, for high working precisions p, the centers of real and complex balls are stored with the full precision, but we only need a small precision for the radii. Consequently, the cost of elementary operations (+, -, •, exp , etc.) on balls is dominated by the cost of the corresponding operations on their centers. Similarly, crude error bounds for products of large matrices can be obtained quickly with respect to the actual product computation, by using norm bounds for the rows and columns of the multiplicands; see (22).

Another algorithmic challenge is to use fast algorithms for the actual numerical computations on the centers of the balls. In particular, it is important to use existing high performance libraries, such as Blas, Linpack, etc., whenever this is possible [START_REF] Rump | Fast and parallel interval arithmetic[END_REF]. Similarly, we should systematically rely on asymptotically efficient algorithms for basic arithmetic, such as fast integer and polynomial multiplication [KO63, [START_REF] Cooley | An algorithm for the machine calculation of complex Fourier series[END_REF][START_REF] Schönhage | Schnelle Multiplikation grosser Zahlen[END_REF]. There are several techniques to achieve this goal:

1. The representations of objects should be chosen with care. For instance, should we rather work with ball matrices or matricial balls (see sections 6.4 and 7.1)?

2. If the result of our computation satisfies an equation, then we may first solve the equation numerically and only perform the error analysis at the end. In the case of matrix inversion, this method is known as Hansen's method; see section 6.2. It should also be noticed that Newton type iterations which are used in this kind of methods are also useful for doubling the accuracy in multiple precision contexts.

3. When considering a computation as a tree or as a dag, then the error tends to increase with the depth of the tree. If possible, algorithms should be designed so as to keep this depth small. Examples will be given in sections 6.4 and 7.5. Notice that this kind of algorithms are usually also more suitable for parallelization. When combining the above approaches to the series of problems considered in this paper, we are usually able to achieve a constant overhead for sharp error bound computations. In favourable cases, the overhead becomes negligible. For particularly ill conditioned problems, we need a logarithmic overhead. It remains an open question whether we have been lucky or whether this is a general pattern.

In this paper, we will be easygoing on what is meant by "sharp error bound". Regarding algorithms as functions f : R k → R l ; x y = f (x), an error δ x ∈ R k in the input automatically gives rise to an error δ y ≈ J f (x) δ x in the output. When performing our computations with bit precision p, we have to consider that the input error δ x is as least of the order of 2 -p |x| ∈ (R ) k (where

|x| i = |x i | for all i). Now given an error bound |δ x | ε x for the input, an error bound |δ y | ε y is considered to be sharp if ε y ≈ max |δ x | ε x |J f (x) δ x |.
More generally, condition numbers provide a similar device for measuring the quality of error bounds. A detailed investigation of the quality of the algorithms presented in this paper remains to be carried out. Notice also that the computation of optimal error bounds is usually NP-hard [START_REF] Kreinovich | Computational Complexity and Feasibility of Data Processing and Interval Computations[END_REF].

Notations and classical facts

Floating point numbers and the IEEE 754 norm

We will denote by

D = Z 2 Z = {m 2 e : m, e ∈ Z}
the set of dyadic numbers. Given fixed bit precisions p ∈ N and q ∈ N for the unsigned mantissa and signed exponent, we also denote by D p,q the corresponding set of floating point numbers. Numbers in D p,q can be stored in p + q + 1 bit space (the extra bit being used for the sign) and correspond to "ordinary numbers" in the IEEE 754 norm [START_REF]ANSI-IEEE Standard[END_REF][START_REF] Muller | Elementary Functions, Algorithms and Implementation[END_REF]. For double precision numbers, we have p = 51 and q = 12. For multiple precision numbers, we usually take q to be the size of a machine word, say q = 64, and denote D p = D p,64 .

The IEEE 754 norm also defines several special numbers. The most important ones are the infinities ±∞ and "not a number" NaN, which corresponds to the result of an invalid operation, such as -2 √ . We will denote

D p,q * = D p,q ∪ {-∞, +∞, NaN}.
Less important details concern the specification of signed zeros ±0 and several possible types of NaNs. For more details, we refer to [START_REF]ANSI-IEEE Standard[END_REF].

The other main feature of the IEEE 754 is that it specifies how to round results of operations which cannot be performed exactly. There are four basic rounding modes ↑ (upwards), ↓ (downwards), (nearest) and 0 (towards zero). Assume that we have an operation f :

U → R with U ⊆ R k . Given x ∈ U ∩ D p,q
k and y = f (x), we define f ↑ (x) to be the smallest number y ˜in D p,q * , such that y ˜ y. More generally, we will use the notation (f (x) + y) ↑ = f ↑ (x) + ↑ y to indicate that all operations are performed by rounding upwards. The other rounding modes are specified in a similar way. One major advantage of IEEE 754 arithmetic is that it completely specifies how basic operations are done, making numerical programs behave exactly in the same way on different architectures. Most current microprocessors implement the IEEE 754 norm for single and double precision numbers. A conforming multiple precision implementation also exists [START_REF] Hanrot | MPFR, a C library for multipleprecision floating-point computations with exact rounding[END_REF].

Classical complexity results

Using Schönhage-Strassen multiplication [START_REF] Schönhage | Schnelle Multiplikation grosser Zahlen[END_REF], it is classical that the product of two p-bit integers can be computed in time I(n) = O(n log n log log n). A recent algorithm by Fürer [START_REF] Fürer | Faster integer multiplication[END_REF] further improves this bound to I(n) = O(n log n 2 O(log * n) ), where log * denotes the iterator of the logarithm (we have log * n = log * log n + 1). Other algorithms, such as the division of two n-bit integers, have a similar complexity O(I(n)). Given k prime numbers p 1 < < p k and a number 0 q < p 1 p k , the computation of q mod p i for i = 1, , k can be done in time O(I(n) log n) using binary splitting [GG02, Theorem 10.25], where n = log (p 1 p k ). It is also possible to reconstruct q from the remainders q mod p i in time O(I(n) log n).

Let K be an effective ring, in the sense that we have algorithms for performing the ring operations in K. If K admits 2 p -th roots of unity for 2 p n, then it is classical [START_REF] Cooley | An algorithm for the machine calculation of complex Fourier series[END_REF] that the product of two polynomials P , Q ∈ K[z] with deg P Q < n can be computed using O(n log n) ring operations in K, using the fast Fourier transform. For general rings, this product can be computed using M(n) = O(n log n log log n) operations [START_REF] Cantor | On fast multiplication of polynomials over arbitrary algebras[END_REF], using a variant of Schönhage-Strassen multiplication.

A formal power series

f = f 0 + f 1 z + f 2 z 2 + ∈ K[[z]
] is said to be computable, if there exists an algorithm which takes n ∈ N on input and which computes f n ∈ K. We denote by K[[z]] com the set of computable power series. Many simple operations on formal power series, such as multiplication, division, exponentiation, etc., as well as the resolution of algebraic differential equations can be done modulo O(z n ) using a similar time complexity O(M(n)) as polynomial multiplication [START_REF] Brent | Fast algorithms for manipulating formal power series[END_REF][START_REF] Bostan | Fast computation of power series solutions of systems of differential equation[END_REF][START_REF] Van Der Hoeven | Newton's method and FFT trading[END_REF]. Alternative relaxed multiplication algorithms of time complexity R(n) = O(M(n)) are given in [vdH02, vdH07a]. In this case, the coefficients (f g) n are computed gradually and (f g) n is output as soon as f 0 , , f n and g 0 , , g n are known. This strategy is often most efficient for the resolution of implicit equations.

Let K m×n denote the set of m × n matrices with entries in K. Given two n × n matrices M , N ∈ K n×n , the naive algorithm for computing M N requires O(n ω ) ring operations with ω = 3. The exponent ω has been reduced by Strassen [Str69] to ω = log 7/log 2. Using more sophisticated techniques, ω can be further reduced to ω < 2.376 [START_REF] Pan | How to multiply matrices faster[END_REF][START_REF] Coppersmith | Matrix multiplication via arithmetic progressions[END_REF]. However, current machine implementations of matrix multiplication are far from reaching this asymptotic complexity.

Notations

Intervals and balls. Given x < y in a totally ordered set R, we will denote by [x, y] the closed interval

[x, y] = {z ∈ R: x z y }. (1)
Assume now that R is a totally ordered field and consider a normed vector space V over R. Given c ∈ V and r ∈ R = {r ∈ R: r 0}, we will denote by

B(c, r) = {x ∈ V : x -c r} (2)
the closed ball with center c and radius r. Notice that we will never work with open balls in what follows. We denote the set of all balls of the form (2) by B(V , R). Given X ∈ B(V , R), we will denote by X c and X r the center resp. radius of x.

Linear algebra. We will use the standard Euclidean norm

x = |x 1 | 2 + + |x n | 2 n
as the default norm on R n and C n . Occasionally, we will also consider the max-norm

x ∞ = max {|x 1 |, , |x n |}.
For matrices M ∈ R m×n or M ∈ C m×n , the default operator norm is defined by

M = max x =1 M x . (3)
Occasionally, we will also consider the max-norm

M ∞ = max i,j |M i,j |, which satisfies M ∞ M n M ∞ .
Polynomials and series. Given a polynomial P ∈ R[x] or P ∈ C[x], we define its maxnorm by

P ∞ = max |P i |. For power series f ∈ R[[z]] or f ∈ C[[z]
] which converge on the compact disk B(0, r), we define the norm

f r = max |z| r |f (z)|. (4)
After rescaling f (z) f (r z), we will usually work with the norm • = • 1 , which should not be confused with the L 1 norm on the coefficients of f .

Generalized radii.

In some cases, it is useful to generalize the concept of a ball and allow for radii in partially ordered rings. For instance, given two vectors X c ∈ R n and X r ∈ (R ) n , we define X = B(X c , X r ) by

X = {x ∈ R n |∀1 i n, |x i -(X c ) i | (X r ) i }.
We denote by B(R n , R n ) the set of all such "balls". The sets

B(R m×n , R m×n ), B(R[x], R[x]), B(R[[x]], R[[x]]
), etc. are defined in a similar componentwise way. Given x ∈ R n , it will also be convenient to denote by |x| the vector with

|x| i = |x i |. For matrices M ∈ R m×n , polynomials P ∈ R[x]

Balls and computable real numbers

Ball arithmetic

Let V be a normed vector space over R and recall that B(V, R) stands for the set of all balls with centers in V and radii in R. Given an operation ϕ: V k → V, the operation is said to lift to an operation

f • : B(V, R) k → B(V, R), if we have {f (x 1 , , x k ): x 1 ∈ X 1 , , x k ∈ X k } ⊆ f • (X 1 , , X k ),
for all X 1 , , X k ∈ B(V, R). For instance, both the addition +: V 2 → V and subtraction -:

V 2 → V admits lifts B(x, r) + • B(y, s) = B(x + y, r + s) (5) B(x, r) -• B(y, s) = B(x -y, r + s). (6)
Similarly, if V is a normed algebra, then the multiplication lifts to

B(x, r) • • B(y, s) = B(x y, (|x| + r) s + r |y |) (7)
The lifts + • , -• , • • , etc. are also said to be the ball arithmetic counterparts of addition, subtractions, multiplication, etc. Although ball arithmetic is really a variant of interval arithmetic, the use of ball enclosures has several advantages, which will be discussed in section 4. Both ball arithmetic and interval arithmetic provide a systematic way for the computation of error bounds when the input of a numerical operation is known only approximately. The bounds are usually not sharp. For instance, consider the mathematical function f : x ∈ R xx which evaluates to zero everywhere, even if x is only approximately known. Taking x = B(2, 0.001), we have x -• x = B(0, 0.002) B(0, 0). This phenomenon is known as overestimation. In general, ball algorithms have to be designed carefully so as to limit overestimation.

In the above definitions, R can be replaced by a subfield R ⊆ R, V by an R-vector space V ⊆ V, and the domain of f by an open subset of V k . If V and R are effective in the sense that we have algorithms for the additions, subtractions and multiplications in V and R, then basic ball arithmetic in B(V , R) is again effective. If we are working with finite precision floating point numbers in D p,q rather than a genuine effective subfield R, we will now show how to adapt the formulas (5), ( 6) and (7) in order to take into account rounding errors; it may also be necessary to allow for an infinite radius in this case.

Finite precision ball arithmetic

Let us detail what needs to be changed when using IEEE conform finite precision arithmetic, say V = R = D p,q . We will denote

B = B(D, D) B p,q = B(D p,q , D p,q ) B[i] p,q = B(D p,q [i], D p,q ).
When working with multiple precision numbers, it usually suffices to use low precision numbers for the radius type. Recalling that D p = D p,64 , we will therefore denote

B p = B(D p , D 64 ) B[i] p = B(D p [i], D 64 ).
We will write ǫ = ǫ p = 2 -p for the machine accuracy and η = η p,q = 2 2-2 q for the smallest normal positive number in D p,q (we might also take η = 2 2-2 q -p to be the smallest strictly positive subnormal number, but this choice deteriorates performance on architectures where computations with subnormal numbers are penalized).

Given an operation f : R k → R as in section 3.1, together with balls

X i = B(x i , r i ), it is natural to compute the center y of B(y, s) = f • (X 1 , , X s )
by rounding to the nearest:

y = f (x 1 , , x k ). (8)
One interesting point is that the committed error

δ = |y -f (x 1 , , x k )|
does not really depend on the operation f itself: we have the universal upper bound

δ ∆(y) ∆(y) 6 (|y | + ↑ η) • ↑ ǫ. (9)
It would be useful if this adjustment function ∆ were present in the hardware.

For the computation of the radius s, it now suffices to use the sum of ∆(y) and the theoretical bound formulas for the infinite precision case. For instance,

B(x, r) + • B(y, s) = B ∆ (x + y, r + ↑ s) (10) B(x, r) -• B(y, s) = B ∆ (x -y, r + ↑ s) (11) B(x, r) • • B(y, s) = B ∆ (x • y, [(|x| + r) s + r |y |] ↑ ), (12)
where B ∆ stands for the "adjusted constructor"

B ∆ (y, s) = B(y, s + ↑ ∆(y)).
The approach readily generalizes to other "normed vector spaces" V over D p,q , as soon as one has a suitable rounded arithmetic in V and a suitable adjustment function ∆ attached to it. Notice that ∆(y) = ∞, if y = ∞ or y is the largest representable real number in D p,q . Consequently, the finite precision ball computations naturally take place in domains of the form B p,q * = B(D p,q * , D p,q * ) rather than B p,q . Of course, balls with infinite radius carry no useful information about the result. In order to ease the reading, we will assume the absence of overflows in what follows, and concentrate on computations with ordinary numbers in B p,q . We will only consider infinities if they are used in an essential way during the computation.

Similarly, if we want ball arithmetic to be a natural extension of the IEEE 754 norm, then we need an equivalent of NaN. One approach consists of introducing a NaB (not a ball) object, which could be represented by B(NaN, NaN). A ball function f • returns NaB if f returns NaN for one selection of members of the input balls. For instance, sqrt • (B(1, 3)) = NaB. An alternative approach consists of the attachment of an additional flag to each ball object, which signals a possible invalid outcome. Following this convention, sqrt • (B(1, 3)) yields B(1, 1) NaN .

Implementation details

Using the formulas from the previous section, it is relatively straightforward to implement ball arithmetic as a C++ template library, as we have done in the Mathemagix system. However, in the case of multiple precision arithmetic this is far from optimal. Let us discuss several possible optimizations:

1. Multiple precision libraries such as Mpfr [HLRZ00] suffer from a huge overhead when it comes to moderate (e.g. quadruple) precision computations. Since the radii are always stored in low precision, it is recommended to inline all computations on the radii. In the case of multiplication, this divides the number of function calls by four.

2. When computing with complex numbers z ∈ D p [i], one may again save several function calls. Moreover, it is possible to regard z as an element of

Z[i] 2 Z rather than (Z 2 Z )[i],
i.e. use a single exponent for both the real and imaginary parts of z. This optimization reduces the time spent on exponent computations and mantissa normalizations.

3. Consider a ball B(c, r) ∈ B p and recall that c ∈ D p , r ∈ D 64 . If |c| < 2 p-64 r, then the ⌊log 2 r + p -log 2 |c| -64⌋ least significant binary digits of c are of little interest. Hence, we may replace c by its closest approximation in D p ′, with p ′ = ⌈log 2 |c| + 64 -log 2 r⌉, and reduce the working precision to p ′ . Modulo slightly more work, it is also possible to share the exponents of the center and the radius.

4. If we don't need large exponents for our multiple precision numbers, then it is possible to use machine doubles D 51,12 as our radius type and further reduce the overhead of bound computations.

When combining the above optimizations, it can be hoped that multiple precision ball arithmetic can be implemented almost as efficiently as standard multiple precision arithmetic. However, this requires a significantly higher implementation effort.

Remark 1. It should be noticed that each of the optimizations described above implies some loss of accuracy. For instance, if z is close to the real axis in the second optimization, then the accuracy of the imaginary part will deteriorate significantly with respect to the usual representation. This kind of examples where quality is traded against efficiency are encountered frequently in the area of reliable computing and we will see more of them in sections 6 and 7. Whether we are willing to accept a lower quality for higher efficiency usually depends on our specific subarea or problem. For instance, our second optimization is usually acceptable in a complex number library, since errors for complex numbers are usually bounded in norm. Furthermore, if the computed result is too inaccurate, then it is always possible to repeat the computation using a more expensive algorithm. [START_REF] Grzegorczyk | On the definitions of computable real continuous functions[END_REF][START_REF] Weihrauch | Computable analysis[END_REF] if there exists an approximation algorithm which takes an absolute tolerance ε ∈ D > on input and which returns an ε-approximation of x. We denote by R com the set of computable real numbers. We may regard R com as a data type, whose instances are represented by approximation algorithms (this is also known as the Markov representation [Wei00, Section 9.6]).

Computable numbers

Given x ∈ R, x ˜∈ D and ε ∈ D > = {ε ∈ D: ε > 0}, we say that x ˜is an ε-approximation of x if |x ˜-x| ε. A real number x ∈ R is said to be computable [Tur36,
In practice, it is more convenient to work with so called ball approximation algorithms: a real number is computable if and only if it admits a ball approximation algorithm, which takes a working precision p ∈ N on input and returns a ball approximation B(c p , r p ) ∈ B with x ∈ B(c p , r p ) and lim p→∞ r p = 0. Indeed, assume that we have a ball approximation algorithm. In order to obtain an ε-approximation, it suffices to compute ball approximations at precisions p = 32, 64, 128, which double at every step, until r p ε. Conversely, given an approximation algorithm x ˇ: D > → D of x ∈ R com , we obtain a ball approximation algorithm p B(c p , r p ) by taking r p = 2 -p and c p = x ˇ(r p ).

Given x, y ∈ R com with ball approximation algorithms x ˇ, y ˇ: N B, we may compute ball approximation algorithms for x + y, xy, x y simply by taking

(x + y )(p) = x ˇ(p) + • y ˇ(p) (x -y )(p) = x ˇ(p) -• y ˇ(p) (x y )(p) = x ˇ(p) • • y ˇ(p).
More generally, assuming a good library for ball arithmetic, it is usually easy to write a wrapper library with the corresponding operations on computable numbers.

From the efficiency point of view, it is also convenient to work with ball approximations. Usually, the radius r p satisfies log 2 r p = -p + o(1),

or at least log 2 r p = -c p + o(1),
for some c ∈ Q > . In that case, doubling the working precision p until a sufficiently good approximation is found is quite efficient. An even better strategy is to double the "expected running time" at every step [START_REF] Van Der Hoeven | Computations with effective real numbers[END_REF][START_REF] Kreinovich | Towards optimal use of multi-precision arithmetic: a remark[END_REF]. Yet another approach will be described in section 3.6 below.

The concept of computable real numbers readily generalizes to more general normed vector spaces. Let V be a normed vector space and let V dig be an effective subset of digital points in V, i.e. the elements of V dig admit a computer representation. For instance, if V = R, then we take V dig = D. Similarly, if V = R m×n is the set of real m × n matrices, with one of the matrix norms from section 2.3, then it is natural to take V dig = (R dig ) m×n = D m×n . A point x ∈ V is said to be computable, if it admits an approximation algorithm which takes ε ∈ D > on input, and returns an ε-approximation x ˜∈ V dig of x (satisfying x ˜x ε, as above).

Asymmetric computability

A real number x is said to be left computable if there exists an algorithm for computing an increasing sequence x ˇ: N D; n x ˇn with lim n→∞ x ˇn = x (and x ˇis called a left approximation algorithm). Similarly, x is said to be right computable if -x is left computable. A real number is computable if and only if it is both left and right computable. Left computable but non computable numbers occur frequently in practice and correspond to "computable lower bounds" (see also [START_REF] Weihrauch | Computable analysis[END_REF][START_REF] Van Der Hoeven | On effective analytic continuation[END_REF]).

We will denote by R lcom and R rcom the data types of left and right computable real numbers. It is convenient to specify and implement algorithms in computable analysis in terms of these data types, whenever appropriate [vdH07b]. For instance, we have computable functions

+: R lcom × R lcom → R lcom -: R lcom × R rcom → R lcom
More generally, given a subset S ⊆ R, we say that x ∈ S is left computable in S if there exists a left approximation algorithm x ˇ: N → S for x. We will denote by S lcom and S rcom the data types of left and computable numbers in S, and define S com = S lcom ∩ S rcom .

Identifying the type of boolean numbers T with {0, 1}, we have T lcom = T rcom = T as sets, but not as data types. For instance, it is well known that equality is non computable for computable real numbers [START_REF] Turing | On computable numbers, with an application to the Entscheidungsproblem[END_REF]. Nevertheless, equality is "ultimately computable" in the sense that there exists a computable function

=: R com × R com → T rcom .
Indeed, given x, y ∈ R com with ball approximation algorithms x ˇand y ˇ, we may take

(x = y ) n = 1 if (x ˇ0 ∩ y ˇ0 ∩ ∩ x ˇn ∩ y ˇn) ∅ 0 otherwise
Similarly, the ordering relation is ultimately computable. This asymmetric point of view on equality testing also suggest a semantics for the relations =, , etc. on balls. For instance, given balls x, y ∈ B, it is natural to take

x = y x ∩ y ∅ x y x ∩ y = ∅ x y ∃a ∈ x, b ∈ y, a b
These definitions are interesting if balls are really used as successive approximations of a real number. An alternative application of ball arithmetic is for modeling non-deterministic computations: the ball models a set of possible values and we are interested in the set of possible outcomes of an algorithm. In that case, the natural return type of a relation on balls becomes a "boolean ball". In the area of interval analysis, this second interpretation is more common [START_REF]IEEE interval standard working group[END_REF].

Remark 2. We notice that the notions of computability and asymmetric computability do not say anything about the speed of convergence. In particular, it is usually impossible to give useful complexity bounds for algorithms which are based on these mere concepts.

In the case of asymmetric computability, there even do not exist any recursive complexity bounds, in general.

Lipschitz ball arithmetic

Given a computable function f : R com → R com , x ∈ R com and ε ∈ D > , let us return to the problem of efficiently computing an approximation y ˜∈ D of y = f (x) with |y ˜y | < ε. In section 3.4, we suggested to compute ball approximations of y at precisions which double at every step, until a sufficiently precise approximation is found. This computation involves an implementation f • : B → B of f on the level of balls, which satisfies

{f (X): x ∈ X } ⊆ f • (X), (13)
for every ball X ∈ B. In practice, f is often differentiable, with f ′ (x) 0. In that case, given a ball approximation X of x, the computed ball approximation Y = f • (X) of y typically has a radius

Y r ∼ |f ′ (x)| X r , (14)
for X r → 0. This should make it possible to directly predict a sufficient precision at which Y r ε. The problem is that (14) needs to be replaced by a more reliable relation. This can be done on the level of ball arithmetic itself, by replacing the usual condition (13) by

B(f (X c ), M X r ) ⊆ f • (X) (15)
M = sup x∈X |f ′ (x)|.
Similarly, the multiplication of balls is carried out using

B(x, r) • • B(y, s) = B(x y, (|x| + r) s + (|y | + s) r). (16)
instead of (7). A variant of this kind of "Lipschitz ball arithmetic" has been implemented in [M0]. Although a constant factor is gained for high precision computations at regular points x, the efficiency deteriorates near singularities (i.e. the computation of 0 √ ).

Balls versus intervals

In the area of reliable computation, interval arithmetic has for long been privileged with respect to ball arithmetic. Indeed, balls are often regarded as a more or less exotic variant of intervals, based on an alternative midpoint-radius representation. Historically, interval arithmetic is also preferred in computer science because it is easy to implement if floating point operations are performed with correct rounding. Since most modern microprocessors implement the IEEE 754 norm, this point of view is well supported by hardware. Not less historically, the situation in mathematics is inverse: whereas intervals are the standard in computer science, balls are the standard in mathematics, since they correspond to the traditional ε-δ-calculus. Even in the area of interval analysis, one usually resorts (at least implicitly) to balls for more complex computations, such as the inversion of a matrix [START_REF] Hansen | Interval arithmetic in matrix computations, part ii[END_REF][START_REF] Moore | Interval Analysis[END_REF]. Indeed, balls are more convenient when computing error bounds using perturbation techniques. Also, we have a great deal of flexibility concerning the choice of a norm. For instance, a vectorial ball is not necessarily a Cartesian product of one dimensional balls.

In this section, we will give a more detailed account on the respective advantages and disadvantages of interval and ball arithmetic.

Standardization

One advantage of interval arithmetic is that the IEEE 754 norm suggests a natural and standard implementation. Indeed, let f be a real function which is increasing on some interval I. Then the natural interval lift f of f is given by

f ([l, h]) = [f ↓ (l), f ↑ (h)].
This implementation has the property that f ([l, h]) is the smallest interval with endpoints in D p,q ∪ {±∞}, which satisfies

{f (x): x ∈ [l, h]} ⊆ f ([l, h]).
For not necessarily increasing functions f this property can still be used as a requirement for the "standard" implementation of f . For instance, this leads to the following implementation of the cosine function on intervals:

cos ([l, h]) =              [cos ↓ l, cos ↑ h] if ⌊l/2 p⌋ = ⌊h/2 p⌋ ∈ 2 Z -1 [cos ↓ h, cos ↑ l] if ⌊l/2 p⌋ = ⌊h/2 p⌋ ∈ 2 Z [min (cos ↓ l, cos ↓ h), 1] if ⌊l/2 p⌋ = ⌊h/2 p⌋-1 ∈ 2 Z -1 [-1, max (cos ↑ l, cos ↑ h)] if ⌊l/2 p⌋ = ⌊h/2 p⌋-1 ∈ 2 Z [-1, 1] if ⌊l/2 p⌋ < ⌊h/2 p⌋-1
Such a standard implementation of interval arithmetic has the convenient property that programs will execute in the same way on any platform which conforms to the IEEE 754 standard.

By analogy with the above approach for standardized interval arithmetic, we may standardize the ball implementation f • of f by taking

f • (B(c, r)) = B(f (c), s),
where the radius s is smallest in D p,q ∪ {+∞} with

{f (x): x ∈ B(c, r)} ⊆ B(f (c), s).
Unfortunately, the computation of such an optimal s is not always straightforward. In particular, the formulas (10), ( 11) and (12) do not necessarily realize this tightest bound. In practice, it might therefore be better to achieve standardization by fixing once and for all the formulas by which ball operations are performed. Of course, more experience with ball arithmetic is required before this can happen.

Practical efficiency

The respective efficiencies of interval and ball arithmetic depend on the precision at which we are computing. For high precisions and most applications, ball arithmetic has the advantage that we can still perform computations on the radius at single precision. By contrast, interval arithmetic requires full precision for operations on both end-points. This makes ball arithmetic twice as efficient at high precisions.

When working at machine precision, the efficiencies of both approaches essentially depend on the hardware. A priori, interval arithmetic is better supported by current computers, since most of them respect the IEEE 754 norm [START_REF]ANSI-IEEE Standard[END_REF], whereas the function ∆ from (9) usually has to be implemented by hand. However, changing the rounding mode is often highly expensive (over hundred cycles). Therefore, additional gymnastics may be required in order to always work with respect to a fixed rounding mode. For instance, if ↑ is our current rounding mode, then we may take [START_REF] Gudenberg | Interval arithmetic on multimedia architectures[END_REF][START_REF] Lambov | Interval arithmetic using sse-2[END_REF] x

+ ↓ y = -((-x) + ↑ (-y)),
since the operation x -x is always exact (i.e. does not depend on the rounding mode). As a consequence, interval arithmetic becomes slightly more expensive. By contrast, when releasing the condition that centers of balls are computed using rounding to the nearest, we may replace (8) by

y = f ↑ (x 1 , , x k ) (17)
and (9) by

∆(y) 6 (|y | + ↑ η) • ↑ (2 ǫ).
Hence, ball arithmetic already allows us to work with respect to a fixed rounding mode. Of course, using (17) instead of (8) does require to rethink the way ball arithmetic should be standardized. An alternative technique for avoiding changes in rounding mode exists when performing operations on compound types, such as vectors or matrices. For instance, when adding two vectors, we may first add all lower bounds while rounding downwards and next add the upper bounds while rounding upwards. Unfortunately, this strategy becomes more problematic in the case of multiplication, because different rounding modes may be needed depending on the signs of the multiplicands. As a consequence, matrix operations tend to require many conditional parts of code when using interval arithmetic, with increased probability of breaking the processor pipeline. On the contrary, ball arithmetic highly benefits from parallel architecture and it is easy to implement ball arithmetic for matrices on top of existing libraries: see [START_REF] Rump | Fast and parallel interval arithmetic[END_REF] and section 6 below.

Quality

Besides the efficiency of ball and interval arithmetic for basic operations, it is also important to investigate the quality of the resulting bounds. Indeed, there are usually differences between the sets which are representable by balls and by intervals. For instance, when using the extended IEEE 754 arithmetic with infinities, then it is possible to represent [1, ∞] as an interval, but not as a ball.

These differences get more important when dealing with complex numbers or compound types, such as matrices. For instance, when using interval arithmetic for reliable computations with complex numbers, it is natural to enclose complex numbers by rectangles X + Y i, where X and Y are intervals. For instance, the complex number z = 1 + i may be enclosed by

z ∈ [1 -ε, 1 + ε] + [1 -ε, 1 + ε] i,
for some small number ε. When using ball arithmetic, we would rather enclose z by

z ∈ B(1 + i, ε).
Now consider the computation of u = z 2 . The computed rectangular and ball enclosures are given by

u ∈ [-2 ε, 2 ε] + [2 -2 ε, 2 + 2 ε] i + o(ε) u ∈ B(2, 2 √ ε) + o(ε).
Consequently, ball arithmetic yields a much better bound, which is due to the fact that multiplication by 1 + i turns the rectangular enclosure by 45 degrees, leading to an overestimation by a factor 2 √ when re-enclosing the result by a horizontal rectangle (see figure 2). The above computation of B(1 + i, ε) 2 is one of the simplest instances of the wrapping effect [START_REF] Moore | Interval Analysis[END_REF]. This example highlights another advantage of ball arithmetic: we have a certain amount of flexibility regarding the choice of the radius type. By choosing a simple radius type, we do not only reduce the wrapping effect, but also improve the efficiency: when computing with complex balls in B[i], we only need to bound one radius instead of two for every basic operation. More precisely, we replace (8) and (9) by

y = f (x 1 , , x k ) = (Re f ) (x 1 , , x k ) + (Im f ) (x 1 , , x k ) i ∆(y) 6 (abs ↑ (y) + ↑ η) • ↑ ( 2 √ ǫ). (18)
Of course, in some very specific situations, interval representations might still be preferred. For instance, consider the rectangle

Z = [a, b] + [-ε, ε] i.
Then any ball enclosure of Z is necessarily very rough. However, when developing a general purpose library for reliable complex numbers, this situation is very exotic (see also remark 1): as soon as we will be multiplying truly complex numbers with non trivial real and imaginary parts, the potential benefit of interval enclosures will be lost as a consequence of the wrapping effect. The choice of a real radius type completely eliminates the wrapping effect in the case of complex balls, due to the fact that |u z | = |u| |z | for any u, z ∈ C. In analogy with the complex case, we may prefer to compute the square of the matrix

M = 0 1 1 0 (19)
in B(D 2×2 , D) rather than B(D, D) 2×2 , while using the operator norm (3) for matrices. Unfortunately, for general 2 × 2 matrices, some overestimation may occur, due to the fact that we only have M N M N for M , N ∈ R 2×2 . Similarly, is not uncommon that the entries of the matrices are of different orders of magnitude, in which case a single bound on the norm leads to a deterioration of the accuracy. Finally, generalized norms may be harder to compute, even though a rough bound often suffices (e.g. replacing abs ↑ (y) by |Re y | + ↑ |Im y | in (18)).

In the area of reliable integration of dynamical systems, various techniques have been developed [Moo66, Loh88, GS88, Neu93, MB96] in order to reduce the wrapping effect induced by matrix computations. Many of these techniques rely on the use of suitable enclosures: linear transformations of higher dimensional rectangles, non linear transformations under Taylor models, ellipsoids, zonotopes, etc. Another, more algorithmic, technique for reducing the wrapping effect will be discussed in sections 6.4 and 7.5 below.

Benefits and cost of correct rounding

Implementations of interval arithmetic often rely on floating point arithmetic with correct rounding. One may question how good correct rounding actually is in order to achieve reliability. One major benefit is that it provides a simple and elegant way to specify what a mathematical function precisely does at limited precision. In particular, it allows numerical programs to execute exactly in the same way on many different hardware architectures.

On the other hand, correct rounding does have a certain cost. Although the cost is limited for field operations and elementary functions [START_REF] Muller | Elementary Functions, Algorithms and Implementation[END_REF], the cost increases for more complex special functions, especially if one seeks for numerical methods with a constant operation count. For arbitrary computable functions on R com , correct rounding even becomes impossible. Another disadvantage is that correct rounding is lost as soon we perform more than one operation: in general, g ↑ • f ↑ and (g • f ) ↑ do not coincide.

In the case of ball arithmetic, we only require an upper bound for the error, not necessarily the best possible representable one. In principle, this is just as reliable and usually more economic. In a similar way, when computing with intervals, we may require faithful rounding instead of correct rounding. That is, we require guaranteed lower and upper bounds for the interval, but not necessarily the best possible ones. Of course, faithful rounding is not unique, so we will lose the benefit of standardization (see section 4.1). Now in an ideal safe world, the development of numerical codes goes hand in hand with the systematic development of routines which compute the corresponding error bounds. In such a world, correct rounding becomes superfluous, since correctness is no longer ensured at the micro-level of hardware available functions, but rather at the top-level, via mathematical proof.

Development of new reliable algorithms

When conceiving a new reliable algorithm, it is often necessary to carry out a lot of tedious bound computations by hand. In the process of software development, not only execution time, but also development time can therefore be a major concern. Depending on the algorithm or application area, as well as personal taste, error bounds may be easier to obtain when using either interval arithmetic or ball arithmetic.

If we are essentially performing computations with real functions on regions where they are monotonic, then interval enclosures are often preferred, since they are both easy to obtain, efficient and tight. Another typical area where interval arithmetic is more convenient is the resolution of a system of equations using dichotomy. Indeed, it is easier to cut an n-dimensional block [a 1 , b 1 ] × × [a n , b n ] into 2 n smaller blocks than to perform a similar operation on balls.

For many other applications, ball representations are more convenient. Indeed, error bounds are usually obtained by perturbation methods. For any mathematical proof where error bounds are explicitly computed in this way, it is generally easy to derive a certified algorithm based on ball arithmetic. We will see several illustrations of this principle in the sections below.

Whenever possible, it is also useful to develop algorithms in such a way that the underlying arithmetic (e.g. interval or ball representation) can be chosen by the user. For instance, in the C++ libraries of the Mathemagix system, many reliable algorithms are templated by a type which can be either be taken to be an interval type or a ball type. Usually, the important feature of the template type is the ability to compute centers, radii, upper or lower bounds, etc., rather than a specific way of doing so.

The numerical hierarchy

Computable analysis provides a good high-level framework for the automatic and certified resolution of analytic problems. The user states the problem in a formal language and specifies a required absolute or relative precision. The program should return a numerical result which is certified to meet the requirement on the precision. A simple example is to compute an ε-approximation for p, for a given ε ∈ D > .

The Mathemagix system [vdHLM+02] aims the implementation of efficient algorithms for the certified resolution of numerical problems. Our ultimate goal is that these algorithms become as efficient as more classical numerical methods, which are usually non certified and only operate at limited precision. A naive approach is to systematically work with computable real numbers. Although this approach is convenient for theoretical purposes in the area of computable analysis, the computation with functions instead of ordinary floating point numbers is highly inefficient.

In order to address this efficiency problem, the Mathemagix libraries for basic arithmetic on analytic objects (real numbers, matrices, polynomials, etc.) are subdivided into four layers of the so called numerical hierarchy (see figure 3). We will illustrate this decomposition on the problem of multiplying two n × n computable real matrices. The numerical hierarchy turns out to be a convenient framework for more complex problems as well, such as the analytic continuation of the solution to a dynamical system. As a matter of fact, the framework incites the developer to restate the original problem at the different levels, which is generally a good starting point for designing an efficient solution.

Analytic problem

Robust problem formal algorithm

Numerical problem interval or ball arithmetic numerical method

Approximate solution

Certified solution

Analytic solution fast arithmetic

Arithmetic problem Exact solution

Mathematical level.

On the mathematical top level, we are given two computable real n × n matrices A, B ∈ (R com ) n×n and an absolute error ε ∈ D > . The aim is to compute ε-approximations for all entries of the product C = A B.

The simplest approach to this problem is to use a generic formal matrix multiplication algorithm, using the fact that R com is an effective ring. However, as stressed above, instances of R com are really functions, so that ring operations in R com are quite expensive. Instead, when working at precision p, we may first compute ball approximations for all the entries of A and B, after which we form two approximation matrices A ˜, B ˜∈ B n×n . The multiplication problem then reduces to the problem of multiplying two matrices in B n×n . This approach has the advantage that O(n 3 ) operations on "functions" in R com are replaced by a single multiplication in B n×n .

Reliable level. The aim of this layer is to implement efficient algorithms on balls. Whereas the actual numerical computation is delegated to the numerical level below, the reliable level should be able to perform the corresponding error analysis automatically.

When operating on non-scalar objects, such as matrices of balls, it is often efficient to rewrite the objects first. For instance, when working in fixed point arithmetic (i.e. all entries of A and B admit similar orders of magnitude), a matrix in B n×n may also be considered as a matricial ball in B(D n×n , D) for the matrix norm • ∞ . We multiply two such balls using the formula

B(x, r) • • B(y, s) = B(x y, [n (( x ∞ + r) s + r y ∞ )] ↑ ),
which is a corrected version of (7), taking into account that the matrix norm • ∞ only satisfies x y ∞ n x ∞ y ∞ . Whereas the naive multiplication in B n×n involves 4 n 3 multiplications in D, the new method reduces this number to n 3 + 2 n 2 : one "expensive" multiplication in D n×n and two scalar multiplications of matrices by numbers. This type of tricks will be discussed in more detail in section 6.4 below.

Essentially, the new method is based on the isomorphism

B(D, D) n×n D B(D n×n , D).
A similar isomorphism exists on the mathematical level:

(R com ) n×n D (R n×n ) com
As a variant, we directly may use the latter isomorphism at the top level, after which ball approximations of elements of (R n×n ) com are already in B(D n×n , D).

Numerical level.

Being able to perform an automatic error analysis, the actual numerical computations are done at the numerical level. In our example, we should implement an efficient algorithm to multiply two matrices in D n×n . In single or double precision, we may usually rely on highly efficient numerical libraries (Blas, Lapack, etc.). In higher precisions, new implementations are often necessary: even though there are efficient libraries for multiple precision floating point numbers, these libraries usually give rise to a huge overhead. For instance, when using Mpfr [HLRZ00] at double precision, the overhead with respect to machine doubles is usually comprised between 10 and 100. When working with matrices with multiple precision floating point entries, this overhead can be greatly reduced by putting the entries under a common exponent using the isomorphism

D n×n = (Z 2 Z ) n×n D Z n×n 2 Z .
This reduces the matrix multiplication problem for floating point numbers to a purely arithmetic problem. Of course, this method becomes numerically unstable when the exponents differ wildly; in that case, row preconditioning of the first multiplicand and column preconditioning of the second multiplicand usually helps.

We notice that a matrix multiplication algorithm has been proposed in [START_REF] Ozaki | Tight and efficient enclosure of matrix multiplication by using optimized blas[END_REF], which is essentially based on a similar principle, but with difference that floating point numbers are cut into smaller floating point numbers instead of integers. This approach can be more efficient when an efficient BLAS library is available and if the precision is not too high. Arithmetic level. After the above succession of reductions, we generally end up with an arithmetic problem such as the multiplication of two n × n matrices in Z n×n . The efficient resolution of this problem for all possible n and integer bit lengths p is again non-trivial.

Indeed, libraries such as Gmp [Gra91] do implement Schönhage-Strassen's algorithm [START_REF] Schönhage | Schnelle Multiplikation grosser Zahlen[END_REF] algorithm for integer multiplication. However, the corresponding naive algorithm for the multiplication of matrices has a time complexity O(I(p) n 3 ), which is far from optimal for large values of n.

Indeed, for large n, it is better to use multi-modular methods. For instance, choosing sufficiently many small primes q 1 , , q k < 2 32 (or 2 64 ) with q 1 q k > 2 n 4 p , the multiplication of the two integer matrices can be reduced to k multiplications of matrices in (Z/q i Z) n×n . Recall that a 2 p-bit number can be reduced modulo all the q i and reconstructed from these reductions in time O(I(p) log p). The improved matrix multiplication algorithm therefore admits a time complexity O(p n 3 + n 2 I(p) log p) and has been implemented in Linbox [DGG+02b, DGG+02a], Mathemagix and several other systems. FFT-based methods achieve similar practical time complexities O(p n 3 + n 2 I(p)) when n and p are of the same order of magnitude.

Reliable linear algebra

In this section, we will start the study of ball arithmetic for non numeric types, such as matrices. We will examine the complexity of common operations, such as matrix multiplication, linear system solving, and the computation of eigenvectors. Ideally, the certified variants of these operations are only slightly more expensive than the non certified versions. As we will see, this objective can sometimes be met indeed. In general however, there is a trade-off between the efficiency of the certification and its quality, i.e. the sharpness of the obtained bound. As we will see, the overhead of bound computations also tends to diminish for increasing bit precisions p.

Matrix multiplication

Let us first consider the multiplication of two n × n double precision matrices M , N ∈ B 51,12 n×n .

Naive strategy. The simplest multiplication strategy is to compute M N using the naive symbolic formula

(M N ) i,k = j =1 n M i,j N j ,k . (20)
Although this strategy is efficient for very small n, it has the disadvantage that we cannot profit from high performance Blas libraries which might be available on the computer.

Revisited naive strategy.

Reconsidering M and N as balls with matricial radii M , N ∈ B(D 51,12 n×n , D 51,12 n×n ),

we may compute M N using

M N = B([M c N c ] , [|M c | N r + M r (|N c | + N r ) + n ǫ |M c | |N c |] ↑ ), (21)
where |M | is given by |M | i,j = |M i,j |. A similar approach was first proposed in [START_REF] Rump | Fast and parallel interval arithmetic[END_REF].

Notice that the additional term n ǫ

|M c | • ↑ |N c | replaces ∆(M N ).
This extra product is really required: the computation of (M c N c ) i,j may involve cancellations, which prevent a bound for the rounding errors to be read off from the end-result. The formula (21) does assume that the underlying Blas library computes M c N c using the naive formula (20) and correct rounding, with the possibility to compute the sums in any suitable order. Less naive schemes, such as Strassen multiplication [START_REF] Strassen | Gaussian elimination is not optimal[END_REF], may give rise to additional rounding errors.

Fast strategy. The above naive strategies admit the disadvantage that they require four non certified n × n matrix products in D 51,12 n×n . If M and N are well-conditioned, then the following formula may be used instead:

M N = B(M c N c , R) (22) R i,k = [ (M i,• ) c (N •,k ) r + (M i,• ) r (N •,k ) c + n ǫ (M i,• ) c (N i,• ) c ] ↑ , ( 23
)
where M i,• and N •,k stand for the i-th row of M and the k-th column of N . Since the O(n) norms can be computed using only O(n 2 ) operations, the cost of the bound computation is asymptotically negligible with respect to the O(n 3 ) cost of the multiplication M c N c .

Hybrid strategy. For large n × n matrices, chances increase that M or N gets badly conditioned, in which case the quality of the error bound (23) decreases. Nevertheless, we may use a compromise between the naive and the fast strategies: fix a not too small constant K, such as K ≈ 16, and rewrite M and N as ⌈ n K ⌉ × ⌈ n K ⌉ matrices whose entries are K × K matrices. Now multiply M and N using the revisited naive strategy, but use the fast strategy on each of the K × K block coefficients. Being able to choose K, the user has an explicit control over the trade-off between the efficiency of matrix multiplication and the quality of the computed bounds.

As an additional, but important observation, we notice that the user often has the means to perform an "a posteri quality check". Starting with a fast but low quality bound computation, we may then check whether the computed bound is suitable. If not, then we recompute a better bound using a more expensive algorithm.

High precision multiplication.

Assume now that we are using multiple precision arithmetic, say M , N ∈ B p n×n . Computing M N using (21) requires one expensive multiplication in D p n×n and three cheap multiplications in D 64 n×n . For large p, the bound computation therefore induces no noticeable overhead.

Matrix inversion

Assume that we want to invert an n × n ball matrix M ∈ B p n×n . This a typical situation where the naive application of a symbolic algorithm (such as gaussian elimination or LRdecomposition) may lead to overestimation. An efficient and high quality method for the inversion of M is called Hansen's method [START_REF] Hansen | Interval arithmetic in matrix computations, part ii[END_REF][START_REF] Moore | Interval Analysis[END_REF]. The main idea is to first compute the inverse of M c using a standard numerical algorithm. Only at the end, we estimate the error using a perturbative analysis. The same technique can be used for many other problems.

More precisely, we start by computing an approximation N c ∈ D p n×n of (M c ) -1 . Putting N = B(N c , 0), we next compute the product M N using ball arithmetic. This should yield a matrix of the form 1 -E, where E is small. If E is indeed small, say E 2 -p/2 , then

(1 -E) -1 -1 -E E 2 1 -E . ( 24
)
Denoting by Ω n the n × n matrix whose entries are all B(0, 1), we may thus take

(1 -E) -1 6 1 + E + E 2 1 -E Ω n .
Having inverted 1 -E, we may finally take M -1 = N (1 -E) -1 . Notice that the computation of E can be quite expensive. It is therefore recommended to replace the check E 2 -p/2 by a cheaper check, such as n E ∞ 2 -p/2 . Unfortunately, the matrix E is not always small, even if M is nicely invertible. For instance, starting with a matrix M of the form

M = J n,K =     1 K 1 K 1     ,
with K large, we have

M -1 =        1 -K K 2 (-K) n-1 1 -K 1 K 2 -K 1       
.

Computing (M c ) -1 using bit precision p, this typically leads to

E ≈ K n-1 2 -p .
In such cases, we rather reduce the problem of inverting 1 -E to the problem of inverting 1 -E 2 , using the formula

(1 -E) -1 = (1 + E) (1 -E 2 ) -1 . ( 25
)
More precisely, applying this trick recursively, we compute E 2 , E 4 , E 8 , until E 2 k becomes small (say E 2 k 2 -p/2 ) and use the formula

(1 -E) -1 = (1 + E) (1 + E 2 ) (1 + E 2 k -1 ) (1 -E 2 k ) -1 (26) (1 -E 2 k ) -1 6 1 + E 2 k + E 2 k 2 1 -E 2 k Ω n .
We may always stop the algorithm for 2 > n, since (J n,K -1) n = 0. We may also stop the algorithm whenever E 2 k 1 and E 2 k E 2 k -1 , since M usually fails to be invertible in that case.

In general, the above algorithm requires O(log m) ball n × n matrix multiplications, where m is the size of the largest block of the kind J K ,m in the Jordan decomposition of M . The improved quality therefore requires an additional O(log n) overhead in the worst case. Nevertheless, for a fixed matrix M and p → ∞, the norm E will eventually become sufficiently small (i.e. 2 -p/2 ) for (24) to apply. Again, the complexity thus tends to improve for high precisions. An interesting question is whether we can avoid the ball matrix multiplication M N altogether, if p gets really large. Theoretically, this can be achieved by using a symbolic algorithm such as Gaussian elimination or LR-decomposition using ball arithmetic. Indeed, even though the overestimation is important, it does not depend on the precision p. Therefore, we have (M -1 ) r = O(2 -p ) and the cost of the bound computations becomes negligible for large p.

Eigenproblems

Let us now consider the problem of computing the eigenvectors of a ball matrix M ∈ B p n×n , assuming for simplicity that the corresponding eigenvalues are non zero and pairwise distinct. We adopt a similar strategy as in the case of matrix inversion. Using a standard numerical method, we first compute a diagonal matrix Λ ∈ D p n×n and an invertible transformation matrix T ∈ D p n×n , such that

T -1 M c T ≈ Λ. ( 27
)
The main challenge is to find reliable error bounds for this computation. Again, we will use the technique of small perturbations. The equation ( 27) being a bit more subtle than M c N c ≈ 1, this requires more work than in the case of matrix inversion. In fact, we start by giving a numerical method for the iterative improvement of an approximate solution.

A variant of the same method will then provide the required bounds. Again, this idea can often be used for other problems. The results of this section are work in common with B. Mourrain and Ph. Trebuchet; in [vdHMT], an even more general method is given, which also deals with the case of multiple eigenvalues.

Given M ˜close to M c , we have to find T ˜close to T and Λ ˜close to Λ, such that

T ˜-1 M ˜T ˜= Λ ˜. ( 28
)
Putting

T ˜= T (1 + E) Λ ˜= Λ (1 + ∆) N = T -1 M ˜T H = N -Λ,
this yields the equation

(1 + E) -1 N (1 + E) = Λ (1 + ∆).
Expansion with respect to E yields

1 1 + E N (1 + E) = N (1 + E) - E 1 + E N (1 + E) = N + N E -E N + E 2 1 + E N - E 1 + E N E = N + [Λ, E] + [H , E] + E 2 1 + E N - E 1 + E N E (29) = N + [Λ, E] + O([H , E]) + O(E 2 ).
Forgetting about the non-linear terms, the equation

[E , Λ] + Λ ∆ = H admits a unique solution E i,j =    H i, j Λ j , j -Λ i,i if j i 0 otherwise ∆ i,j =    H i,i Λ i,i if i = j 0 otherwise Setting κ = κ(Λ) = max max 1 |Λ j,j -Λ i,i | : 1 i < j n , max 1 Λ i,i : 1 i n , it follows that max { E , ∆ } κ n √ H . ( 30
)
Setting 29) also implies

T ′ = T (1 + E), Λ ′ = Λ (1 + ∆) and H ′ = (T ′ ) -1 M ˜T ′ -Λ ′ , the relation (
H ′ = [H , E] + E 2 1 + E (Λ + H) - E 1 + E (Λ + H) E.
Under the additional condition E 1 2 , it follows that

H ′ 3 H E + 4 E 2 Λ . ( 31
)
For sufficiently small H, we claim that iteration of the mapping Φ: (T , Λ) (T ′ , Λ ′ ) converges to a solution of (28).

Let us denote (T (k) ,

Λ (k) ) = Φ k (T , Λ), H (k) = (T (k) ) -1 M ˜T (k) -Λ (k) and let (E (k) , ∆ (k) ) be such that T (k+1) = T (k) (1 + E (k) ) and Λ (k+1) = Λ (k) (1 + ∆ (k) ). Assume that H 1 80 n κ 2 Λ (32)
and let us prove by induction over k that

H (k) 2 -k H (33) Λ (k) 2 Λ (34) max { E (k) , ∆ (k) } 2 n √ κ H (k) (35) 1 40 n √ κ Λ 2 k
This is clear for k = 0, so assume that k > 0. In a similar way as (31), we have

H (k) 3 H (k-1) E (k -1) + 4 E (k -1) 2 Λ (k-1) . ( 36
)
Using the induction hypotheses and κ Λ 1, it follows that

H (k) (3 + 16 n √ κ Λ ) E (k -1) H (k-1)
1 2 H (k-1) , which proves (33). Now let Σ (k) be such that

Λ (k) = Λ (1 + ∆ (0) ) (1 + ∆ (k -1) ) = Λ (1 + Σ (k) ).
From (35), it follows that

Σ (k) 1 4 κ Λ .
On the one hand, this implies (34). On the other hand, it follows that κ(Λ (k) ) 2 κ, whence (30) generalizes to (35). This completes the induction and the linear convergence of H (k) to zero. In fact, the combination of ( 35) and (36) show that we even have quadratic convergence.

Let us now return to the original bound computation problem. We start with the computation of H = T -1 M T -Λ using ball arithmetic. If the condition (32) is met (using the most pessimistic rounding), the preceding discussion shows that for every M ˜∈ M (in the sense that M ˜i,j ∈ M i,j for all i, j), the equation (28) admits a solution of the form

T ˜= T (1 + Υ ˜) = T (1 + E (0) ) (1 + E (1) ) Λ ˜= Λ (1 + Σ ˜) = Λ (1 + ∆ (0) ) (1 + ∆ (1) ) , with max { E (k) , ∆ (k) } 2 1-k n √ κ H , for all k. It follows that max { Υ ˜ , Σ ˜ } η 6 6 n √ κ H .
We conclude that

T ˜∈ T (1 + η Ω n ) Λ ˜∈ B(1, η) Λ.
We may thus return (T (1 + η Ω n ), B(1, η) Λ) as the solution to the original eigenproblem associated to the ball matrix M . The reliable bound computation essentially reduces to the computation of three matrix products and one matrix inversion. At low precisions, the numerical computation of the eigenvectors is far more expensive in practice, so the overhead of the bound computation is essentially negligible. At higher precisions p, the iteration (T , Λ) Φ(T , Λ) actually provides an efficient way to double the precision of a numerical solution to the eigenproblem at precision p/2. In particular, even if the condition (32) is not met initially, then it usually can be enforced after a few iterations and modulo a slight increase of the precision. For fixed M and p → ∞, it also follows that the numerical eigenproblem essentially reduces to a few matrix products. The certification of the end-result requires a few more products, which induces a constant overhead. By performing a more refined error analysis, it is probably possible to make the cost certification negligible, although we did not investigate this issue in detail.

Matricial balls versus ball matrices

In section 4.3, we have already seen that matricial balls in B(D n×n , D) often provide higher quality error bounds than ball matrices in B(D, D) n×n or essentially equivalent variants in B(D n×n , D n×n ). However, ball arithmetic in B(D n×n , D) relies on the possibility to quickly compute a sharp upper bound for the operator norm M of a matrix M ∈ D n×n . Unfortunately, we do not know of any really efficient algorithm for doing this.

One expensive approach is to compute a reliable singular value decomposition of M , since M coincides with the largest singular value. Unfortunately, this usually boils down to the resolution of the eigenproblem associated to M * M , with a few possible improvements (for instance, the dependency of the singular values on the coefficients of M is less violent than in the case of a general eigenproblem).

Since we only need the largest singular value, a faster approach is to reduce the computation of M to the computation of M * M , using the formula

M = M * M .
Applying this formula k times and using a naive bound at the end, we obtain

M n (M * M ) 2 k -1 ∞ 2 k
. This bound has an accuracy of ≈k -O(log 2 n) bits. Since M * M is symmetric, the k -1 repeated squarings of M * M only correspond to about k 2 matrix multiplications. Notice also that it is wise to renormalize matrices before squaring them, so as to avoid overflows and underflows.

The approach can be speeded up further by alternating steps of tridiagonalization and squaring. Indeed, for a symmetric tridiagonal matrix D, the computation of D 2 and its tridiagonalization only take O(n) steps instead of O(n 3 ) for a full matrix product. After a few k = O(n 2 ) steps of this kind, one obtains a good approximation µ of D . One may complete the algorithm by applying an algorithm with quadratic convergence for finding the smallest eigenvalues of Dµ. In the lucky case when D has an isolated maximal eigenvalue, a certification of this method will provide sharp upper bounds for D in reasonable time.

Even after the above improvements, the computation of sharp upper bounds for M remains quite more expensive than ordinary matrix multiplication. For this reason, it is probably wise to avoid ball arithmetic in B(D n×n , D) except if there are good reasons to expect that the improved quality is really useful for the application in mind.

Moreover, when using ball arithmetic in B(D n×n , D n×n ), it is often possible to improve algorithms in ways to reduce overestimation. When interpreting a complete computation as a dag, this can be achieved by minimizing the depth of the dag, i.e. by using an algorithm which is better suited for parallelization. Let us illustrate this idea for the computation of the k-th power M k of a matrix M . When using Horner's method (multiply the identity matrix k times by M ), we typically observe an overestimation of O(k) bits (as for the example ( 19)). If we use binary powering, based on the rule

M k = M ⌊k/2⌋ M ⌈k/2⌉ , (37)
then the precision loss drops down to O(log k) bits. We will encounter a less obvious application of the same idea in section 7.5 below.

Reliable power series arithmetic

Ball series versus serial balls

There are two typical applications of power series

f ∈ R[[z]] or f ∈ C[[z]
] with certified error bounds. When f occurs as a generating function in a counting problem or random object generator, then we are interested in the computation of the coefficients f n for large n, together with reliable error bounds. A natural solution is to systematically work with computable power series with ball coefficients in B p [[z]] com . For many applications, we notice that p is fixed, whereas n ≫ p may become very large.

The second typical application is when

f ∈ C[[z]]
is the local expansion of an analytic function on a disk B(0, ρ) and we wish to evaluate f at a point z with |z | < ρ. The geometric decrease of |f n z n | implies that we will need only n = O(p) coefficients of the series. In order to bound the remaining error using Cauchy's formula, we do not only need bounds for the individual coefficients f n , but also for the norm f ρ defined in (4). Hence, it is more natural to work with serial balls in B(D p [i][[z]] com , D 64 ), while using the • ρ norm. Modulo a rescaling f (z) f (ρ z), it will be convenient to enforce ρ = 1. In order to compute sharp upper bounds ⌈⌈f ⌉⌉ for f = f 1 , it will also be convenient to have an algorithm which computes bounds ⌈⌈f n; ⌉⌉ for the tails

f n; = f n z n + f n+1 z n+1 + .
Compared to the computation of the corresponding head

f ;n = f 0 + + f n-1 z n-1 ,
we will show in section 7.4 that the computation of such a tail bound is quite cheap.

Again the question arises how to represent f ;n in a reliable way. We may either store a global upper bound for the error, so that f ;n ∈ B(D p [i][z], D 64 ), or compute individual bounds for the errors, so that f

;n ∈ B[i] p [z]. If our aim is to evaluate f at a point z with |z | ≈ 1, then both representations P ∈ B[i] p [z] and P ˆ= B(P c , P r ∞ ) ∈ B(D p [i][z], D 64 )
give rise to evaluations P (z), P ˆ(z) ∈ B[i] p with equally precise error bounds. Since the manipulation of global error bounds is more efficient, the corresponding representation should therefore be preferred in this case. In the multivariate case, one has the additional benefit that "small" coefficients f i z i (e.g. |f i | ǫ p f ) can simply be replaced by a global error B(0, |f |), thereby increasing the sparsity of f . On the other hand, individual error bounds admit the advantage that rescaling f (z) f (λ z) is cheap. If we suddenly find out that f is actually convergent on a larger disk and want to evaluate f at a point z with |z | > 1, then we will not have to recompute the necessary error bounds for f ;n from scratch.

Serial ball representations similar to what has been described above are frequently used in the area of Taylor models [START_REF] Makino | Remainder differential algebras and their applications[END_REF][START_REF] Makino | Suppression of the wrapping effect by Taylor model-based validated integrators[END_REF] for the validated long term integration of dynamical systems. In the case of Taylor models, there is an additional twist: given a dynamical system of dimension d, we not only compute a series expansion with respect to the time t, but also with respect to small perturbations ε 1 , , ε d of the initial conditions. In particular, we systematically work with power series in several variables. Although such computations are more expensive, the extra information may be used in order to increase the sharpness of the computed bounds. A possible alternative is to compute the expansions in ε 1 , , ε d only up to the first order and to use binary splitting for the multiplication of the resulting Jacobian matrices on the integration path. This approach will be detailed in section 7.5.

Reliable multiplication of series and polynomials

In order to study the reliable multiplication of series f and g, let us start with the case when f , g ∈ B(D p [[z]] com , D 64 ), using the sup-norm on the unit disk. We may take

h = f g = B(f c g c , [ f c g r + f r ( g c + g r ) + δ] ↑ ),
where h c = f c g c stands for a δ-approximation of f c g c . Since h c is really a numeric algorithm for the computation of its coefficients, the difficulty resides in the fact that δ has to be chosen once and for all, in such a way that the bound h cf c g c will be respected at the limit. A reasonable choice is δ = ǫ p f c g c . We next distribute this error over the infinity of coefficients: picking some α < 1, each coefficient (h c ) n is taken to be an [(1α) α n δ]-approximation of f c g c . Of course, these computation may require a larger working precision than p. Nevertheless, f and g are actually convergent on a slightly larger disk B(0, ρ). Picking α = 1/ρ, the required increase of the working precision remains modest.

Let us now turn our attention to the multiplication of two computable series f , g ∈ B p [[z]] com with ball coefficients. Except for naive power series multiplication, based on the formula (f g) n = i+ j=n f i g j , most other multiplication algorithms (whether relaxed or not) use polynomial multiplication as a subalgorithm. We are thus left with the problem of finding an efficient and high quality algorithm for multiplying two polynomials P , Q ∈ B p [z] of degrees < n. In order to simplify the reading, we will assume that P 0 , P n-1 , Q 0 , Q n-1 are all non-zero.

As in the case of matrix multiplication, there are various approaches with different qualities, efficiencies and aptitudes to profit from already available fast polynomial arithmetic in D p [z]. Again, the naive O(n 2 ) approach provides almost optimal numerical stability and qualities for the error bounds. However, this approach is both slow from an asymptotic point of view and unable to rely on existant multiplication algorithms in D p [z].

If the coefficients of P and Q are all of the same order of magnitude, then we may simply convert P and Q into polynomial balls in B(D p [z], D 64 ) for the norm • ∞ and use the following crude formula for their multiplication:

P Q = B(P c Q c , [n P c ∞ Q r ∞ + n P r ∞ ( Q c ∞ + Q r ∞ ) + δ] ↑ ), (38)
where P c Q c stands for a δ-approximation of P c Q c . In other words, we may use any efficient multiplication algorithm in D p [z] for the approximation of P c Q c , provided that we have a means to compute a certified bound δ for the error.

In our application where P and Q correspond to ranges of coefficients in the series f and g, we usually have P i ≈ P 0 ρ f -i and Q i ≈ Q 0 ρ g -i for the convergence radii ρ f and ρ g of P and Q. In order to use (38), it is therefore important to scale P (z) P (z/ρ) and Q(z) Q(z/ρ) for a suitable ρ. If we are really interested in the evaluation of f g at points z in a disk B(0, r), then we may directly take ρ = r. In we are rather interested in the coefficients of f g, then ρ = min (ρ f , ρ g ) is the natural choice. However, since ρ f and ρ g are usually unknown, we first have to compute suitable approximations for them, based on the available coefficients of P and Q. A good heuristic approach is to determine indices i < n/2 j such that

P k P i j -i P j P i k-i
, for all k, and to take

ρ f ≈ P i P i 1 j -i
as the approximation for ρ f . Recall that the numerical Newton polygon of N P is the convex hull of all points (i, log |P i |λ) ∈ R 2 with λ 0. Consider the edge of N P whose projection on the x-axis contains n/2. Then i and j are precisely the extremities of this projection, so they can be computed in linear time.

For large precisions n = O(p), the scaling algorithm is both very efficient and almost of an optimal quality. For small p and large n, there may be some precision loss which depends on the nature of the smallest singularities of f and g. Nevertheless, for many singularities, such as algebraic singularities, the precision loss is limited to O(log n) bits. For a more detailed discussion, we refer to [vdH02, Section 6.2].

In other applications, where P and Q are not obtained from formal power series, it is usually insufficient to scale using a single factor ρ. This is already the case when multiplying P = 1 + z and Q = 1 + δ z for small δ ≪ ǫ p , since the error bound for (P Q) 2 = δ exceeds ǫ p . One possible remedy is to "precondition" P and Q according to their numerical Newton polygons and use the fact that N PQ is close to the Minkowski sum N P + N Q .

More precisely, for each i, let (N P ) i > 0 denote the number such that (i, log (N P ) i ) lies on one of the edges of N P . Then (N P ) i is of the same order of magnitude as P i , except for indices for which P i is small "by accident". Now consider the preconditioned relative error ν P = max i (P i ) r (N P ) i and similarly for Q. Then

[(P Q) i ] r n (ν P + ν Q + ν P ν Q ) (N P + N Q ) i , (39)
if P Q is computed using infinite precision ball arithmetic. Assuming a numerically stable multiplication algorithm for the centers, as proposed in [vdH08], and incorporating the corresponding bound for the additional errors into the right hand side of (39), we thus obtain an efficient way to compute an upper bound for (P Q) r . Notice that the numerical Newton polygon N P has a close relation to the orders of magnitudes of the roots of P . Even though the error bounds for some "accidentally small" coefficients (P Q) i may be bad for the above method, the error bounds have a good quality if we require them to remain valid for small perturbations of the roots of P Q.

Reliable series division and exponentiation

The algorithms from the previous section can in particular be used for the relaxed multiplication of two ball power series f , g ∈ B p [[z]] com . In particular, using the implicit equation

g = 1 + z f g, (40)
this yields a way to compute g = (1z f ) -1 . Unfortunately, the direct application of this method leads to massive overestimation. For instance, the computed error bounds for the inverses of

g = 1 1 -3 z -2 z 2 h = 1 1 -3 z + 2 z 2 coincide.
Indeed, even when using a naive relaxed multiplication, the coefficients of g and h are computed using the recurrences

g n = 3 g n-1 + 2 g n-2 h n = 3 h n-1 -2 h n-2 ,
but the error bound ε n for g n and h n is computed using the same recurrence

ε n = 3 ε n-1 + 2 ε n-2 ,
starting with ε 0 ≈ ǫ p . For n → ∞, it follows that g n ∼ λ α -n and ε n ∼ λ ǫ p α -n for some λ, where α ≈ 0.281 is the smallest root of 1 -3 α -2 α 2 . Hence the error bound for g n is sharp. On the other hand, h n ∼ µ β n for some µ, where β = 1 2 is the smallest root of 1 -3 β + 2 β 2 . Hence, the error bound ε n is n log 2 (β/α) bits too pessimistic in the case of h. The remedy is similar to what we did in the case of matrix inversion. We first introduce the series

ϕ = 1 1 -z f c ∈ D p [[z]] com ψ = [ϕ (z f -1)] 1; z ∈ B p [[z]] com
and next compute g using

g = ϕ 1 -z ψ ,
where 1z ψ is inverted using the formula (40). This approach has the advantage of being compatible with relaxed power series expansion and it yields high quality error bounds. Another typical operation on power series is exponentiation. Using relaxed multiplication, we may compute the exponential g of an infinitesimal power series f ∈ z B p [[z]] com using the implicit equation

g = 1 + f ′ g, (41)
where stands for "distinguished integration" in the sense that ( h) 0 = 0 for all h. Again, one might fear that this method leads to massive overestimation. As a matter of fact, it usually does not. Indeed, assume for simplicity that f r = 0, so that f ∈ D p [[z]]. Recall that |f | denotes the power series with |f | n = |f n |. Roughly speaking, the error bound for g n , when computed using the formula (41) will be the coefficient ε n of the power series ε = exp (ǫ p |f |). Since |f | has the same radius of convergence as f , it directly follows that the bit precision loss is sublinear o(n). Actually, the dominant singularity of |f | often has the same nature as the dominant singularity of f . In that case, the computed error bounds usually become very sharp. The observation generalizes to the resolution of linear differential equations, by taking a square matrix of power series for f .

Automatic tail bounds

In the previous sections, we have seen various reliable algorithms for the computation of the expansion f ;n of a power series f ∈ C[[z]] up till a given order n. Such expansions are either regarded as ball polynomials in

B[i] p [z] or as polynomial balls in B(D p [i][z], D 64 ).
Assuming that f is convergent on the closed unit disk B(0, 1), it remains to be shown how to compute tail bounds ⌈⌈f n; ⌉⌉. We will follow [vdH07b]. Given a ball polynomial P , then we notice that reasonably sharp upper and lower bounds ⌈⌈P ⌉⌉ and ⌊⌊P ⌋⌋ for P can be obtained efficiently by evaluating P at O(deg P ) primitive roots of unity using fast Fourier transforms [vdH07b, Section 6.2]. We will assume that the series f is either an explicit series, the result of an operation on other power series or the solution of an implicit equation. Polynomials are the most important examples of explicit series. Assuming that f is a polynomial of degree < d, we may simply take

⌈⌈f n; ⌉⌉ = ⌈⌈z -n f n;d ⌉⌉ if n < d 0 otherwise , (42)
where

f n;d = f n z n + + f d-1 z d-1 .
For simple operations on power series, we may use the following bounds:

⌈⌈(f + g) ;n ⌉⌉ = ⌈⌈f ;n ⌉⌉ + ⌈⌈g ;n ⌉⌉ (
⌈⌈( f ) n; ⌉⌉ = 1 n + 1 ⌈⌈f n; ⌉⌉, (46)
where

⌈⌈(f ;n g ;n ) n; ⌉⌉ k=0 n-1 |f k | l=n-k n-1 |g l |
can be computed in time O(n). Due to possible overestimation, division has to be treated with care. Given an infinitesimal power series ε and

f = 1 1 -ε , so that f n; = 1 + ε f ;n -f ;n 1 -ε , we take ⌈⌈f n; ⌉⌉ = ⌈⌈(ε f ;n ) n; ⌉⌉ ⌊⌊1 -ε ;n + B(0, ⌈⌈ε n; ⌉⌉)⌋⌋ , (47)
where ⌈⌈(ε f ;n ) n; ⌉⌉ is computed using (45). Let Φ(f ) be an expression which is constructed from f and polynomials, using the ring operations and distinguished integration (we exclude division in order to keep the discussion simple). Assume that each coefficient Φ(f ) n only depends on the previous coefficients f 0 , , f n-1 of f and not on f n , f n-1 , . Then the sequence 0, Φ(0), Φ 2 (0), tends to the unique solution f of the implicit equation

f = Φ(f ), (48)
Moreover, Φ k (0) n = f n for any k > n. In (40) and (41), we have already seen two examples of implicit equations of this kind.

The following technique may be used for the computation of tail bounds ⌈⌈f n; ⌉⌉. Given c ∈ D and assuming that ⌈⌈f n; ⌉⌉ c, we may use the rules (42-47) in order to compute a "conditional tail bound" ⌈⌈Φ(f ) n; |c⌉⌉ for Φ(f ) n; . Mathematically speaking, this bound has the property that for any power series g with g ;n = f ;n and g n; c, we have

Φ(g) n; ⌈⌈Φ(g) n; |c⌉⌉. If ⌈⌈Φ(g) n; |c⌉⌉ c, then it follows that [Φ (k) (f ;n )] ;n c for all k. Given any disk B(0, ρ) with ρ < 1, it follows that Φ (k) (f ;n ) -f ρ 2 c ρ k /(1 -ρ) for any k n, since Φ (k) (f ;n ) -f = O(z k ).
In other words, Φ (k) (f ;n ) uniformly converges to f on B(0, ρ). Therefore, f n; ρ c for any ρ < 1. Modulo application of the method on a disk of radius slightly larger than one and scaling back, this yields a way to compute a certified tail bound ⌈⌈f n; ⌉⌉.

Let us now consider the mapping ϕ: c ⌈⌈Φ(f ) n; |c⌉⌉ and assume that Φ involves no divisions. When computing ⌈⌈Φ(f ) n; |c⌉⌉ using infinite precision and the rules (42-45), we notice that ϕ is a real analytic function, whose power series expansion contains only positive coefficients. Finding the smallest c with ϕ(c) c thus reduces to finding the smallest fixed point c fix of ϕ, if such a fixed point exists. We may use the secant method:

c 0 6 0 c 1 6 ϕ(c 0 ) c k+2 6 c k + ϕ(c k ) -c k c k+1 -ϕ(c k+1 ) + ϕ(c k ) -c k (c k+1 -c k )
If c k+1 < c k for some k or if k exceeds a given threshold K, then the method fails and we set ⌈⌈f n; ⌉⌉ = +∞. Otherwise, c k converges quadratically to c fix . As soon as |c k+1 /c k -1| < ε, for some given ε > 0, we check whether ϕ(c ˜fix ) c ˜fix for c ˜fix = 2 c k+1c k , in which case we stop. The resulting c ˜fix is an approximation of c fix with relative accuracy ε > 0.

Assuming that Ψ(f ) ;n has been computed for every subexpression of Φ, we notice that the computation of ⌈⌈Φ(f ) n; |c⌉⌉ only requires O(n s) floating point operations, where s is the size of Φ as an expression. More generally, the evaluation of ⌈⌈Φ(f ) n; |c i ⌉⌉ for k different hypotheses c 1 , , c k only requires O(n s k) operations, since the heads Ψ(f ) ;n do not need to be recomputed. Our algorithm for the computation of c ˜fix therefore requires at most O(n s K) floating point operations. Taking K = o(R(n)/n), it follows that the cost of the tail bound computation remains negligible with respect to the series expansion itself.

The approach generalizes to the case when f is a vector or matrix of power series, modulo a more involved method for the fixed point computation [vdH07b, Section 6.3]. If f is indeed convergent on B(0, 1), then it can also be shown that ϕ indeed admits a fixed point if n is sufficiently large.

Long term integration of dynamical systems

Let us now consider a dynamical system

f = f (0) + Φ(f ), (49)
where f is a vector of d unknown complex power series, f (0) ∈ B[i] p d , and Φ an expression built up from the entries of f and polynomials in D p [i][z] using the ring operations. Given t ∈ D[i] \ {0}, denote by f z ×t the scaled power series f (t z) and by f z+t the shifted power series f (t + z), assuming that f converges at t. Also denote by Φ z ×t and Φ z+t the expressions which are obtained when replacing each polynomial P in Φ by P z ×t resp. P z+t . Then f z ×t and f z+t satisfy

f z ×t = f (0) + t Φ z ×t (f z ×t ) (50)
f z+t = f (t) + Φ z+t (f z+t ) (51)
Since ( 49) is a particular case of an implicit equation of the form (48), we have an algorithm for the computation of tail bounds ⌈⌈f n; ⌉⌉ ∈ (D 64 ∪ {∞}) d for f on B(0, 1). Modulo rescaling (50), we may also compute tail bounds ⌈⌈f n; ⌉⌉ ρ on more general disks B(0, ρ).

Assume that we want to integrate (49) along the real axis, as far as possible, and performing all computations with an approximate precision of p bits. Our first task is to find a suitable initial step size t and evaluate f at t. Since we require a relative precision of approximately p bits, we roughly want to take t ≈ ρ f /2, where ρ f is the radius of convergence of f , and evaluate the power series f up to order n ≈ p. We thus start by the numerical computation of f n; and the estimation ρ ˜f of ρ f , based on the numerical Newton polygon of f n; . Setting t 6 ρ ˜f/2, we next compute a tail bound ⌈⌈f n; ⌉⌉ t . This bound is considered to be of an acceptable quality if ⌈⌈f n; ⌉⌉ t ǫ p ⌈⌈f ;n ⌉⌉ t .

In the contrary case, we keep setting t 6 t/2 and recomputing ⌈⌈f n; ⌉⌉ t , until we find a bound of acceptable quality. It can be shown [vdH07b] that this process ultimately terminates, when p is sufficiently large. We thus obtain an appropriate initial step size t which allows us to compute a δ-approximation of f (t) with δ = 2 ǫ p ⌈⌈f ;n ⌉⌉ t + ↑ (f ;n (t)) r .

In principle, we may now perform the translation z t + z and repeat the analytic continuation process using the equation (51). Unfortunately, this approach leads to massive overestimation. Indeed, if the initial condition f (0) is given with relative precision η, then the relative precisions of the computed coefficients f k are typically a non trivial factor times larger than η, as well as the next "initial condition" f (t) at t. Usually, we therefore lose O(1) bits of precision at every step.

One remedy is the following. Let ∆ = ∆ 0,t be the analytic continuation operator which takes an initial condition c ∈ C d on input and returns the evaluation f (t) ∈ C d of the unique solution to (49) with f (0) = c. Now instead of computing f ;n using a ball initial condition f (0) ∈ B[i] p d , we rather use its center f (0) c as our initial condition and compute ∆(f (0) c ) using ball arithmetic. In order to obtain a reliable error bound for f (t), we also compute the Jacobian J ∆ of ∆ at f (0) c using ball arithmetic, and take

f (t) = ∆(f (0) c ) + B(0, [J ∆ (f (0) c ) f (0) r ] ↑ ).
The Jacobian can either be computed using the technique of automatic differentiation [START_REF] Baur | The complexity of partial derivatives[END_REF] or using series with coefficients in a jet space of order one.

With this approach ∆(f (0) c ) is computed with an almost optimal p-bit accuracy, and J ∆ (f (0) c ) with an accuracy which is slightly worse than the accuracy of the initial condition. In the lucky case when J ∆ (f (0) c ) is almost diagonal, the accuracy of f (t) will therefore be approximately the same as the accuracy of f (0). However, if J ∆ (f (0) c ) is not diagonal, such as in (19), then the multiplication J ∆ (f (0) c ) f (0) r may lead to overestimation. This case may already occur for simple linear differential equations such as

f g = 1 0 + 0 -1 1 0 f g (52)
and the risk is again that we lose O(1) bits of precision at every step. Let us describe a method for limiting the harm of this manifestation of the wrapping effect. Consider the analytic continuations of f at successive points 0 = t 0 < t 1 < < t k and denote J i = J ∆ t i-1 ,t i (f (t i-1 ) c ), i = 1, , k.

Instead of computing the error at t i due to perturbation of the initial condition f (0) using the formula

[f (t) r ] per(0) = J k (J k-1 ( J 1 f (0) r )),

we may rather use the formula

[f (t) r ] per(0) = (J k J k-1 J 1 ) f (0) r ,
where matrix chain products are computed using a variant of binary powering (37):

J j J i = (J j J ⌊(i+j)/2⌋+1 ) (J ⌊(i+ j)/2⌋ J i ).

In order to be complete, we also have to take into account the additional error δ i , which occurs during the computation of ∆ t i-1 ,t i (f (t i-1 ) c ). Setting δ 0 = f (0) r , we thus take

f (t i ) r = δ i + J i δ i-1 + + (J k J 1 ) δ 0 . (53)
When using this algorithm, at least in the case of simple systems such as (52), the precision loss at step k will be limited to O(log k) bits. Notice that we benefit from the fact that the Jacobian matrices remain accurate as long as the initial conditions remain accurate. Although we have reduced the wrapping effect, the asymptotic complexity of the above algorithm is cubic or at least quadratic in k when evaluating (53) using a binary splitting version of Horner's method. Let us now describe the final method which requires only O(k log k) matrix multiplications up till step k. For 0 i < j, we define δ i,j = δ j -1 + J j -1 δ j -2 + + (J j -1 J i+1 ) δ i . J i,j = J j -1 J i , where J 0 = 1. At stage i = 2 e 1 + + 2 e l with e 1 > > e l , we assume that δ [j] 6 δ [i;j] 6 δ 2 e 1 + +2 e j -1 ,2 e 1 + +2 e j J [j] 6 J [i;j] 6 J 2 e 1 + +2 e j -1 ,2 e 1 + +2 e j are stored in memory for all 1 j l. From these values, we may compute

f (t i-1 ) r = δ [l] + J [l] (δ [l-1] + J [3] (δ [2] + J [2] δ [1] ) )
using only O(l) = O(log i) matrix multiplications. Now assume that we go to the next stage i + 1. If m is maximal such that 2 m divides i + 1, then i + 1 = 2 e 1 + + 2 e l-m + 2 m . Consequently, δ [i+1;j] = δ [i;j] and J [i+1;j] = J [i;j] for j < l ′ 6 lm + 1 and

δ [i+1;l ′ ] = δ i + J i (δ [i;l] + J [i;l ′ +2] (δ [i;l ′ +1] + J [i;l ′ +1] δ [i;l ′ ] ) ) J [i+1;l ′ ] = J i J [i;l] J [i;l ′ ] .
Hence, the updated lists δ [i+1;j] and J [i+1;j] can be computed using O(m) = O(log i) matrix multiplications. Furthermore, we only need to store O(log i) auxiliary matrices at step i.

Discussion

There is a vast literature on validated integration of dynamical systems and reduction of the wrapping effect [Moo65, Moo66, Nic85, Loh88, GS88, Neu93, K98, Loh01, Neu02, MB04]. We refer to [START_REF] Lohner | On the ubiquity of the wrapping effect in the computation of error bounds[END_REF] for a nice review. Let us briefly discuss the different existing approaches.

The wrapping effect was noticed in the early days of interval arithmetic [START_REF] Moore | Automatic local coordinate transformations to reduce the growth of error bounds in interval computation of solutions to ordinary differential equation[END_REF] and local coordinate transformations were proposed as a remedy. The idea is to work as much as possible with respect to coordinates in which all errors are parallel to axes. Hence, instead of considering blocks x ∈ B p d , we rather work with parallelepipeds x = c + T B(0, r), with c ∈ D p d , T ∈ D p d×d and r ∈ D p d . A natural choice for T after k steps is T = J k J 1 , but more elaborate choices may be preferred [START_REF] Lohner | Einschließung der Lösung gewöhnlicher Anfangs-und Randwertaufgaben und Anwendugen[END_REF]. Other geometric shapes for the enclosures have been advanced in the literature, such as ellipsoids [START_REF] Neumaier | The wrapping effect, ellipsoid arithmetic, stability and confedence regions[END_REF], which are also invariant under linear transformations, and zonotopes [K98]. However, as long as we integrate (49) using a straightforward iterative method, and even if we achieve a small average loss ε ≪ 1 of the bit precision at a single step, the precision loss after k steps will be of the form k ε.

The idea of using dichotomic algorithms in order to reduce the wrapping effect was first described in the case of linear differential equations [START_REF] Gambill | Logarithmic reduction of the wrapping effect with application to ordinary differential equations[END_REF]. The previous section shows how to adapt that technique to non linear equations. We notice that the method may very well be combined with other geometric shapes for the enclosures: this will help to reduce the precision loss to (log k) ε instead of k ε in the above discussion.

Notice that there are two common misunderstandings about the dichotomic method. Contrary to what we have shown above (see also [START_REF] Gambill | Logarithmic reduction of the wrapping effect with application to ordinary differential equations[END_REF] in the linear case), it is sometimes believed that we need to keep all matrices J 1 , , J k in memory and that we have to reevaluate products J j J i over and over again. Secondly, one should not confuse elimination and reduction of the wrapping effect: if M is the 2 × 2 rotation matrix by a 30 • angle, then none of its repeated squarings M 2 l will be the identity matrix, so every squaring (M 2 l ) 2 will involve a wrapping. Even though we have not eliminated the wrapping effect, we did achieve to reduce the number of wrappings to l instead of 2 l .

Taylor models are another approach for the validated long term integration of dynamical systems [EKW84, EMOK91, MB96, MB04]. The idea is to rely on higher k-th order expansions of the continuation operator ∆. This allows for real algebraic enclosures which are determined by polynomials of degree k. Such enclosures are a priori more precise for non linear problems. However, the method requires us to work in order k jet spaces in d variables; the mere storage of such a jet involves d + k d coefficients. From a theoretical point of view it is also not established that Taylor methods eliminate the wrapping effect by themselves. Nevertheless, Taylor models can be combined with any of the above methods and the non linear enclosures seem to be more adapted in certain cases. For a detailed critical analysis, we refer to [START_REF] Neumaier | Taylor forms -use and limits[END_REF].

Let us finally investigate how sharp good error bounds actually may be. If ρ f (0) = f (0) r / f (0) c denotes the relative precision of the initial condition at the start and ρ f (t) = f (t) r / f (t) c the relative precision of the final result, then it is classical that ρ f (t) κ(J ∆ (f (0))) ρ f (0) , κ(J ∆ (f (0))) = J ∆ (f (0)) J ∆ (f (0)) -1 , where κ(J ∆ (f (0))) is the condition number of the matrix J ∆ (f (0)). We propose to define the condition number κ(Φ, f (0), 0, t) for the integration problem (49) on [0, t] by K = κ(Φ, f (0), 0, t) = max 0 t 1 t 2 t κ(J ∆ t 1 ,t 2 (f (t 1 ))).

Indeed, without using any particular mathematical properties of Φ or f , we somehow have to go through the whole interval [0, t]. Of course, it may happen that Φ is indeed special. For instance, if Φ = M f for a matrix M with a special triangular or diagonal form, then special arguments may be used to improve error bounds and more dedicated condition numbers will have to be introduced.

However, in general, we suspect that a precision loss of log K cannot be avoided. If K gets large, the only real way to achieve long term stability is to take p > 2 log K (say), whence the interest of efficient multiple precision and high order ball algorithms. It seems to us that the parallelepiped method should manage to keep the precision loss bounded by log K, for p > 2 log K and ε ≈ ǫ p . The algorithm from section 7.5 (without coordinate transformations) only achieves a log k log K in the worst case, although a log K bound is probably obtained in many cases of interest. We plan to carry out a more detailed analysis once we will have finished a first efficient multiple precision implementation.

Figure 1 .

 1 Figure 1. A sample Mathemagix session for the reliable computation of the Bell number B 1000 .

  and power series f ∈ R[[x]], we define |M |, |P | and |f | similarly.

Figure 2 .

 2 Figure 2. Illustration of the computation of z 2 using interval and ball arithmetic, for z = 1 + i.

Figure 3 .

 3 Figure 3. The numerical hierarchy.

Joris van der Hoeven

Ball arithmetic

* . This work has been supported by the ANR-09-JCJC-0098-01 MaGiX project, as well as a Digiteo 2009-36HD grant and Région Ile-de-France.