<\body> >||<\author-address> Mathématiques, CNRS Bâtiment 425 Université Paris-Sud 91405 Orsay Cedex France ||>> <\abstract> Until now, the area of symbolic computation has mainly focused on the manipulation of algebraic expressions. Based on earlier, theoretical work, the author has started to develop a systematic library for mathematically correct computations with more analytic objects, like complex numbers and analytic functions. While implementing the library, we found that several of our theoretical ideas had to be further improved or adapted. In this paper, we report on the current implementation, we present several new results and suggest directions for future improvements. Although the field of symbolic computation has given rise to several softwares for mathematically correct computations with algebraic expressions, similar tools for analytic computations are still somewhat inexistent. Of course, a large amount of software for numerical analysis does exist, but the user generally has to make several error estimates by hand in order to guarantee the applicability of the method being used. There are also several systems for interval arithmetic, but the vast majority of them works only for fixed precisions. Finally, several systems have been developed for certified arbitrary precision computations with polynomial systems. However, such systems cannot cope with transcendental functions or differential equations. The central concept of a systematic theory for certified computational analysis is the notion of an . Such a number \> is given by an which takes \\=\*2>> with \0> on input and which produces an >-approximation> \\> for with -x\|\\>. One defines effective complex numbers in a similar way. Effective real and complex numbers are a bit tricky to manipulate: although they can easily be added, multiplied, , there exists no test for deciding whether an effective real number is identically zero. Some interesting questions from a complexity point of view are also raised: if we want to compute an >-approximation of +x>, how to determine +\=\> so that the computation of >-approximations of the > is most efficient? Concrete approaches and implementations for computations with effective real numbers have been proposed, often independently, by several authors . A first step in these approaches is often to implement some interval arithmetic. As an optional second step, one may then provide a class for which a real number is given by an which, given \\>>, computes a closed interval\x> with endpoints in >, of radius>\>. In this paper, we report on the implementation of a C++ class for effective real numbers in the library . This implementation is based on, but it also contains some new ideas. In section , we start by quickly reviewing interval arithmetic and the computationally more efficient variant of ``ball arithmetic'' (see also ). We also try to establish a more precise semantics for this kind of arithmetic in the multi-precision context and discuss the use of interval classes as parameters of template classes such as complex numbers, matrices or polynomials. Our implementation relies on the and libraries. In section , we mainly review previous work: equivalent definitions of effective real numbers, representation by dags and the different techniques for priori> and error estimations. We also state improved versions of the ``global approximation problem'', which serve as a base for the complexity analysis of our library. We finally correct an error which slipped into . In section , we describe the current implementation, which is based on the sole use of error estimates. In section, we prove that our implementation is optimal up to a linear overhead in the input size() and a logarithmic overhead in the time complexity(). It is interesting to compare these results with previous theoretical work on the complexity of interval arithmetic . In the last section, we indicate how to use error estimates in a more efficient way than in . In a forthcoming paper, we hope to work out the details and remove the linear overhead() in the input size, at least under certain additional assumptions. Since real numbers cannot be stored in finite memory, a first approach to certified computations with real numbers is to compute with intervals instead of numbers. For instance, when using interval arithmetic, a real number like > would be approximated by an interval like =[|\>,|\>]=[3.141592,3.141593]>. Evaluation of a real function at a point \=[|\>,|\>]> then corresponds to finding an interval =[|\>,|\>]> with \> for all \>. When all functions under consideration are Lipschitz continuous and when all computations are performed using a sufficiently high precision, this technique provides arbitrarily good approximations for the real numbers we want to compute. Our current interval library is based on MPFR, which implements a generalization of the IEEE754 standard. More precisely, let {32,64,128}> be the word precision of the computer under consideration. Given a bit precision 2> and a maximal exponent 2> fixed by the user, the MPFR library implements arithmetic on -bit floating point numbers and exponents in the range m,\,m-1,m]> with exact rounding. We recall that the IEEE754 norm (and thus MPFR) includes special numbers 0>, \> and (not a number). Assuming that has been fixed once and for all, we will denote by > the set of -bit numbers. Several representations are possible for the computation with -bit intervals: In this representation, an -bit interval =[|\>,|\>]> is determined by its end-points |\>\|\>\\\{NaN}>. We also allow for the exceptional value =NaIn=[NaN,NaN]>. Since the MPFR library provides exact rounding, it is particularly simple to base an interval library on it, when using this representation: whenever a function is monotonic on a certain range, it suffices to consider the values at the end-points with opposite rounding modes. However, every high precision evaluation of an interval function requires two high precision evaluations of floating point numbers. In this representation, an -bit interval \NaIn> is represented by a ball =\(c>,r>)> with center >\\\{NaN}> and radius >\\>>. If W> and \NaIn>, then we also require that the ball (c>,r>)> is in the sense that >\2*\|c>\|>. The underlying idea is that the endpoints of a high precision interval are usually identical apart from a few bits at the end, whence it is more efficient to only store the common part and the difference. As a consequence, one high precision operation on balls reduces to one high precision operation on floating point numbers and several operations on low precision numbers. However, the ball representation has essentially less expressive power. For instance, it is impossible to represent the interval 0,\\]>. Also, positivity is not naturally preserved by addition if . This may be problematic in the case we want to compute quantities like ,\)=+\>>. <\remark> The normality condition admits several variants. For instance, given \\{\0,\\,NaN}>, we define the > as the exponent of minus . A ball =\(c>,r>)> may then be called normal if >>\\>>> instead of >\2*\|c>\|>. Alternatively, one may require >\{0,\,2-1}*2>>>>. In , we currently use the endpoint representation for low precision numbers W> and the ball representation for high precision numbers W>. Modulo some additional overhead, this combines the advantages of both representations while removing their drawbacks. Let > denote the set of -bit intervals for the hybrid representation. If W>, then it should be noticed that the set > is not stable under the usual arithmetic operations, due to the phenomenon of . Indeed, the sum =\+\> of ,\\\> is typically computed using >=c>+>c>> and >=r> +> r>+>\>, where > is a small bound for the rounding error. However, if \\\>, then> is not necessarily normal. Nevertheless, the set =k\l>\> is stable under all usual arithmetic operations. Indeed, any ball or interval > may be normalized by replacing >> by a lower bit approximation. More precisely, consider an abnormal interval |\>,|\>]> with >\\> and >\\>. Given \>, let |\>=2*\2k>*|\>\> and |\>=2*\2k>*|\>\>, so that >,r>\2*\>. Let \> be minimal such that >\\>>, >\\> and >\2>\|c>\|> for some \{W+1,\,l}>. Then =norm(\)> is called the of >. If no such exists, then )> is defined to be the smallest interval with endpoints in >, which contains >. It can be shown that \>>> and )>\(1+2W>)*r>>. <\remark> The normalization procedure is not very canonical and it admits several variants (in particular, it has to be adapted whenever we change the definition of normality). Nevertheless, all good normalization procedures share the property that )>\(1+B*2W>)*r>> for some small fixed constant \>. Dually, it may occur that the result of an operation can be given with more precision than the argument. For instance, if =\(1,0.5)*2>, then > can be computed with a precision of about binary digits. Similarly, > can be computed with a precision of about binary digits. We call this the phenomenon of . The point here is that it is not necessarily desirable to compute the results of > and > with the maximal possible precision. Taking into account the phenomena of precision loss and gain, we propose the following ``ideal'' semantics for operations on intervals. First of all, a default precision for computations is fixed by the user. Now assume that we wish to evaluate an -ary function|^>\|^>> at intervals ,\,\\\>, where |^>=\\{\\,NaN}>. Then there exists a smallest interval > with |\>,|\>\\>, which satisfies either <\itemize> ,\,x)\\> for all \\,\,x\\>. =NaIn> and ,\,x)=NaN> for some \\,\,x\\>. Whenever 0,\\}\\\\>, then =\(\,\,\)> is taken to be the smallest interval of > which contains>. Otherwise, we take =norm(\)>. <\remark> For some purposes one may use the alternative convention for exceptions that > is the smallest interval which contains all non exceptional values. Since the normalization procedure is somewhat arbitrary (see remark ), the ideal semantics may be loosened a little bit for implementation purposes. Instead of requiring the optimal return value =\(c>,r>)>, we rather propose to content oneself with a return value |~>=\(c|~>>,r|~>>)> with |~>>-c>\|\(1+B*2l>)*\|c>\|> and |~>>-r>\|\(1+B*2l>)*\|c>\|>, for some fixed small constant \>, in the case when > has precision W>. This remark also applies to the underlying MPFR layer: exact rounding is not really necessary for our purpose. It would be sufficient to have a ``looser'' rounding mode, guaranteeing results up to times the last bit. A tempting way to implement the complex analogue of interval arithmetic in C++ is to use the template class from the standard library. Unfortunately, this approach leads to a lot of overestimation for the bounding rectangles. In order to see this, consider the complex rectangle >+\>*\> and the sequence =x>, =x*a>. Because multiplication with > ``turns'' the bounding rectangle, the error > is roughly multiplied by > at each step. In other words, we loose one bit of precision every two steps. The above phenomenon can be reduced in two ways. First of all, one may use a better algorithm for computing>, like repeated squaring. In the case of complex numbers though, the best solution is to systematically use complex ball representations. However, standardization of the operations requires more effort. Indeed, given an operation on balls ,\,\> of precision , it can be non-trivial to design an algorithm which computes a ball \f(\,\,\)> of almost minimal radius (up to W>>). The precision loss phenomenon is encountered more generally when combining interval arithmetic with template types. The best remedy is again to modify the algorithms and/or data types in a way that the errors in the data are all of a similar order of magnitude. For instance, when computing a monodromy matrix as a product *\*\> of connection matrices, it is best to compute this product by dichotomy *\*\k/2\>)*(\k/2\+1>*\*\)>. Similarly, when computing the product of two truncated power series, it is good to first perform a change of variables z/\> which makes the errors in the coefficients of and of the same order of magnitude For users of computer algebra systems, it would be convenient to provide a data type for real numbers which can be used in a similar way as the types of rational numbers, polynomials, etc. Since interval arithmetic already provides a way to perform certified computations with ``approximate real numbers'', this additional level of abstraction should mainly be thought of as a convenient interface. However, due to the fact that real numbers can only be represented by infinite structures like Cauchy sequences, their manipulation needs more care. Also, the actual implementation of a library of functions on effective real numbers raises several interesting computational complexity issues. In this section, we review some previous work on this matter. Let =\*2>> denote the set of dyadic numbers. Given \> and \\>>, we recall from the introduction that an >-approximation> of is a dyadic number \\> with -x\|\\>. We say that is , if it admits an , which takes \\>> on input and which returns an >-approximation> for . The of such an approximation algorithm is the time it takes to compute a n>>-approximation for , when \>. We denote by > the set of effective real numbers. The above definition admits numerous variants . For instance, instead of require an approximation algorithm, one may require the existence of an algorithm which associates a closed interval =[|\>,|\>]> with end-points in > to each \>, such that \\\\> and \> r>=0> (an interval \x> with end-points in > will also be called an >-bounding interval> for ). Similarly, one may require the existence of an effective and rapidly converging Cauchy sequence \\;n\x> for which there exists a number \>> with -x\|\M*2n>> for all . All these definitions have in common that an effective real number is determined by an algorithm which provides more and more precise approximations of on demand. In an object oriented language like , this can be implemented by providing an abstract representation class with a purely virtual method which corresponds to this approximation algorithm. The class is implemented as a pointer to . Since effective real numbers should be thought of as algorithms, the zero-test problem in > can be reduced to the halting problem for Turing machines. Consequently, there exist no algorithms for the basic relations >, >>, >>, >>, >> and >> on >. Given an open domain > of )>, a real function \\> is said to be if there exists an algorithm >> which takes an approximation algorithm >=(>,\,>)> for ,\,x)\\> on input and which produces an approximation algorithm >> for ,\,x)>. Here we understand that >=>(>)> approximates the same number , if >> is another approximation algorithm for . Most common operations, like >, >, >>, , , , >, >, , can easily shown to be effective. On the other hand, without any of the operations for comparison, it seems more difficult to implement functions like \x\>. In fact, it turns out that effective real functions are necessarily continuous . A concrete library for computations with effective real numbers consists of a finite number of functions like , , >, >, >>, , etc. Given inputs ,\,x> of type , such an operation should produce a new instance ,\,x)> of . Usually, the representation class for in particular contains members for ,\,x>, which can then be used in the method which implements the approximation algorithm for . For instance, a very simple implementation of addition might look as follows: <\small> <\code> class add_real_rep: public real_rep { \ \ real x, y; public: \ \ add_real_rep (const real& x2, const real& y2): \ \ \ \ x (x2), y (y2) {} \ \ dyadic approximate (const dyadic& err) { \ \ \ \ return x-\approximate (err/2) + y-\approximate (err/2); \ \ } }; When implementing a library of effective real functions ,f,\> in this way, we notice in particular that any effective real number computed by the library reproduces the expression by which it was computed in memory. Such effective real numbers may therefore be modeled faithfully by rooted dags (directed acyclic graphs) , whose nodes are labeled by ,f,\>. More generally, finite sets of effective real numbers can be modeled by general dags of this type. Figure shows an example of such a dag, together with some parameters for measuring its complexity. Since storing entire computations in memory may require a lot of space, the bulk of computations should not be done on the effective real numbers themselves, but rather in their approximation methods. In particular, should not be thought of as some kind of improved type, which can be plugged into existing numerical algorithms: the class rather provides a user-friendly high-level interface, for which new algorithms need to be developed. <\with|par-mode|center> <\small-figure> |gr-frame|>|gr-geometry||gr-grid||gr-grid-old||1>|gr-edit-grid-aspect|||>|gr-edit-grid||gr-edit-grid-old||1>|gr-text-at-halign|center||>>|>|>>||>|>|>>|>>|>>||>||>||>|>>|>>||>|||>|||>||>||>>|||>>> <|small-figure> Example of a dag with roots. The dag has ( the total number of nodes) and 3 ( the longest path from a root to a leaf). The of the dag corresponds to the sum of the sizes of the trees obtained by ``copying'' each of the roots. In our example, the weight is . Finally, the of the dag is defined to be the maximum number of ancestors of a leaf. In our example, the ancestrality is . Let ,f,\> be a library of effective real functions as in the previous section, based on a corresponding library ,\,\> of functions on intervals. In order to study the efficiency of our library, it is important to have a good model for the computational complexity. In this section, we will describe a static and a dynamic version of the , which are two attempts to capture the computational complexity issues in a precise manner. In its static version, the input of the global approximation problem consists of <\itemize> A dag, whose nodes are labeled by ,f,\> Achallenge >\\>\{\\}> for each node \\>. Denoting by >> the function associated to the node > and by ,\,\\|>> its children, we may recursively associate a real value >=f>(x>,\,x\|>>)> to >. On output, we require for each node\G> an interval >> with endpoints in >, such that <\itemize> >\\>> and >>\\>>. For certain >\\>,\,\\|>>\\\|>>>, we have <\equation*> \>=\>(\>,\,\\|>>). Notice that the second condition implies in particular that >\\>(\>,\,\\|>>)>. The dynamic version of the global approximation problem consists of a sequence of static global approximation problems for a fixed labeled dag , when we incrementally add challenges >> for nodes >. More precisely, we are given <\itemize> A dag , whose nodes are labeled by ,f,\> A finite sequence ,\),\,(\,\)> of pairs ,\)\G\(\>\{\})>. On output, we require for each {1,\,k}> a solution to the -th static global approximation problem, which consists of the labeled dag with challenges >=min {\:j\i,\=\}>. Here we understand that a solution at the stage may be presented as the set of changes the solution at stage . Let us explain why we think that the dynamic global approximation problem models the complexity of the library in an adequate way. For this, consider a computation by the library. The set of all effective real numbers constructed during the computation forms a labeled dag. The successive calls of the approximation methods of these numbers naturally correspond to the sequence ,\),\,(\,\)>. It is reasonable to assume that the library itself does not construct any real numbers, all nodes of correspond to explicit creations of real numbers by the user. Indeed, if new numbers are created from inside an approximation method, then all computations which are done with these numbers can be seen as parts of the approximation method, so they should not be taken into account during the complexity analysis. Similarly, if the constructor of a number ,\,x)> induces the construction of other real numbers, then ,\,x)> may be expressed in terms of more basic real functions, so we may consider as a function outside our library. Now assume that another, possibly better library were used for the same computation. It is reasonable to assume that the corresponding dag > and challenges ,\),\,(a>,\>)> coincide with the previous ones. Indeed, even though it might happen that the first and second library return different bounding intervals >> and >> for a given challenge ,\)>, the libraries cannot know what the user wants to do with the result. Hence, for a fair comparison between the first and second library, we should assume that the user does not take advantage out of possible differences between >> and >>. This reduces to assuming that =G>, =k> and ,\)=(a,\)> for all . Finally, it is reasonable to assume that all actual approximations >> of the >> are done using a fixed interval library ,\,\>. This means for instance that the second library has no better algorithms for multiplication, exponentiation, than the first one. When putting all our ``reasonable assumptions'' together, the time of the computation which was spent in the library now corresponds to the time which was required to solve the corresponding dynamic global approximation problem. error estimates> Let us now consider the problem of obtaining an >-approximation for the result of an operation f(x,\,x)>. For simplicity, we will focus on the case of addition +x>. In this and the next section, we briefly recall several strategies, which are discussed in more detail in. In the case of error estimates, the tolerance> is distributed over >and>. In other words, we first determine > and > with +\\\>, next compute >-approximations for the >, and finally add the results. The systematic choice of =\=\/2> can be very inefficient: in the case of badly balanced trees like in figure (this occurs in practice when evaluating polynomials using Horner's rule), this requires the approximation of > with a much lower tolerance than > (/2> versus /2>). This problem can be removed by balancing the error according to the weights > and > of > and > (, by taking =\*w/(w+w)>). For ``non degenerate'' cases of the global approximation problem for a dag of weight and size , it can be shown this technique requires tolerances which are never worse than times the optimal ones. Unfortunately, while implementing the algorithms from, it turned out that is often of the same order as and therefore far from good enough. This is for instance the case when the expressions are obtained some iterative process or as the coefficient of a lazy power series. For this reason, we have currently abandoned the use of error estimates in our implementation. However, this situation is quite unsatisfactory, since this technique is still most efficient in many cases. We will come back to this problem in section . <\with|par-mode|center> | > || > ||>| > >>>|A badly balanced tree.> error estimates and relaxed evaluation> A second strategy consists of computing error estimates : if we want to compute an >-approximation for +x>, we start with the computation of a bounding interval for at precision . As long as the obtained result is not sufficiently precise, we keep doubling the precision and repeating the same computation. As explained in, this strategy can be optimized in two ways. First of all, the strategy may be carried out locally, by storing a ``best available approximation'' (together with the corresponding precision) for each instance of . Indeed, when increasing the precision for the computation of , sufficiently precise approximations for > and > might already be known, in which case their recomputation is unnecessary. Secondly, instead of doubling the precision at each step, it is better to double the expected computation time. For instance, consider the computation of , where has time complexity>\*n>> ( admits an >\*n>+T(n+O(1))> approximation algorithm, whenever admits a approximation algorithm). Evaluate at successive precisions >,2>,\,2>>, where \*log n\> and is the smallest precision at which the evaluation yields a sufficiently precise result. Then the total computation time >\+2*\+\+2*\\2*\> never exceeds *n>\2*\> by a factor more than (see also ). Unfortunately, an error slipped into , because the successive recursive approximations of may not be sufficiently precise in order to allow for evaluations of at successive precisions >,2>,\,2>>. For instance, if is given by an algorithm of exponential time complexity >, then successive approximations of will only yield one more bit at every step. This error can be repaired up to a logarithmic factor in two ways. First of all, we notice that the error only concerns the cumulative cost of the successive reevaluation of . In section, we will prove that the total cost of reevaluating nodes of the dag remains good. Secondly, it is possible to adapt the technique of relaxed formal power series to real numbers. Roughly speaking, this approach relies on the recursive decomposition of a ``relaxed mantissa'' of length into a fixed part > of length \l/2> and a relaxed remainder > (so that +x>). Given an operation , we then compute )> and (x)> at precision > and obtain a formula for the relaxed decomposition +y>, since > is a truncation of )> and =f(x)-y+f(x)*x>. As soon as the precision of > exceeds2>, we take a new value for > and recompute )> and (x)> at a doubled precision. Working out the details of this construction shows that most common real functions can be evaluated in a relaxed way with the same complexity as usual, multiplied by an overhead. However, the relaxed strategy accounts for a lot of additional implementation work and no noticeable improvement with respect to the global bound () which will be proved in section. Therefore, it is mainly interesting from a theoretical point of view. > and > Inside , dyadic numbers in > are represented using generalized floating point numbers in>, where is bounded by a precision of the order of > or >. Effective real numbers (of type ) are implemented as pointers to an abstract representation class with a virtual method for the computation of >-bounding intervals. Usually, such a number is of the form ,\,x)>, where is an effective real function and ,\,x> are other effective real numbers. The number is concretely represented by an instance of a class which derives from and with fields corresponding to ,\,x>. The current implementation is based on the technique of error bounds from section with the two optimizations mentioned there: remembering the best currently available approximations for each real number and doubling computations times instead of precisions. These strategies are reflected as follows in the data type: <\small> <\code> class real_rep { protected: \ \ double \ \ cost; \ \ interval best; \ \ real_rep (): cost (1.0) { compute (); } \ \ virtual int as_precision (double cost); \ \ virtual interval compute (); public: \ \ interval improve (double new_cost); \ \ interval approximate (const dyadic& err); }; The field corresponds to the best currently available bounding interval for . The value of is recomputed several times by the purely virtual method at increasing costs, the last one of which is stored in . More precisely, is recomputed as a function of approximations ,\,\> of ,\,x> at the same costs. When these approximations are sufficiently precise, then the cost of the computation of will be more or less equal to . Otherwise, the actual computation may take less time (see the discussion at the end of section). The costs are normalized (we start with ) and doubled at each iteration. The purely virtual method is used to convert an intended cost to the corresponding intended precision. The user interface is given by the routines and . The first one computes an approximation of at intended cost : <\small> <\code> interval real_rep::improve (double new_cost) { \ \ if (new_cost \= cost) return best; \ \ cost= max (new_cost, 2.0 * cost); \ \ set_precision (as_precision (cost)); \ \ best= compute (); \ \ restore_precision (); \ \ return best; } The method returns an >-bounding interval > for as a function of >: <\small> <\code> interval real_rep::approximate (const dyadic& eps) { \ \ while (radius (best) \= eps) \ \ \ \ (void) improve (2.0 * cost); \ \ return best; } <\remark> In practice, the method also avoids the call of if the new precision associated to is equal to the old one. This may indeed happen if the cost of the operation increases more than linearly as a function of the bit precision. > Let us illustrate the mechanism from the previous section in the case of exponentiation. The exponential of a number is represented by an instance of <\small> <\code> class exp_real_rep: public real_rep { \ \ real x; \ \ int as_precision (double cost); \ \ interval compute (); public: \ \ exp_real_rep (const real& x2): x (x2) {} }; The computation of bits of takes a time proportional to > for small values of and a more or less linear time for large values of . Therefore, a simple implementation of would be as follows: <\small> <\code> int exp_real_rep::as_precision (double cost) { \ \ if (cost \= 256.0) return (int) sqrt (cost); \ \ return min (MAX_PREC, (int) cost/16.0); } Of course, this is a very rough approximation of the real time complexity of . For the theoretical bounds in the next sections, better approximations are required. In practice however, a simple implementation like the above one is quite adequate. If necessary, one may implement a more precise algorithm, based on benchmarks. One may also gradually increase precisions and use a timer. The actual approximation of is done using the overloaded function on intervals: <\small> <\code> interval exp_real_rep::compute () { \ \ return exp (x-\improve (cost)); } In the case of functions with arity more than one it is often possible to avoid unnecessarily precise computations of one of the arguments, when the approximations of the other argument are far less precise. For instance, in the case of addition, may be implemented as follows: <\small> <\code> interval add_real_rep::compute () { \ \ dyadic eps= pow (2.0, -BITS_IN_WORD) \ \ if (radius (y-\best) \ eps * radius (x-\best)) { \ \ \ \ (void) x-\improve (cost); \ \ \ \ while (y-\cost \ cost && \ \ \ \ \ \ \ \ \ \ \ radius (y-\best) \= eps * radius (x-\best)) \ \ \ \ \ \ (void) y-\improve (2.0 * y-\cost); \ \ } \ \ else \ \ else return x-\improve (cost) + y-\improve (cost); } Notice that our implementation, as outlined in the previous section, naturally solves both the static and the dynamic versions of the global approximation problem: we first construct the dag and then either compute an >>-approximation for each >>, or successive >-approximations for each > (,k>). In this section, we examine the efficiency of this approach. Since >> is approximated several times during our algorithm, let us first study the difference between the total computation time and the time taken by the final and most precise approximations of the >>.\ For each node >, let ,0>,\,t,p>>> be the successive timings for the approximation of >>. We will also denote by ,0>\\\T,p>>> the intended computation times and precisions. By construction, we have ,i>\T,i>> for all and >=1,\,T,p>>=2>>>. For each >, let >=t,0>+\+t,p>>\T>=T,0>+\+T,p>>> and >=t,p>>>. We define \G>t>>, \G>T>> and =\G>t>>. We already warned against the possibility that ,i>\T,i>>. Nevertheless, we necessarily have ,i>=T,i>> if > is a leaf. Also, any operation >> of cost ,i>> triggers an operation of cost ,i>> for one of the children of >. By induction, it follows that there exists at least one leaf ,i>> descending from> which spends a time ,i>,i>=T,i>,i>=T,i>>. Hence, denoting by the ancestrality of the dag and by > its subset of leaves, we have <\equation> T\\G,i>T,i>,i>*\a*\\,i>T,i>=a*\\,i>t,i>\a*t. We also have 2*T> since ,1>+\+T,p>>\2*T,p>>> for all>. Consequently, <\equation> >*T\*t\2*T. The bound () is sharp in the case when the dag has only one leaf> and ) the computation of an digit approximation of >> requires exponential time; ) all other operations can be performed in linear or polynomial time. A similar situation occurs when cancellations occur during the computation of >>, in which case the computation of>> at many bits precision still produces a -bit result. A variant of (), which is usually better, is obtained as follows. Since the precision of the result of an operation on intervals increases with the precision of the arguments, and similarly for the computation times, we have ,1>\\\t,p>>>. Let> be a node (which can be assumed to be a leaf by what precedes) for which >> is maximal. Then <\equation*> t=\G,i>t,i>\\G>p>*t>\p*t. It follows that <\equation> t\t\(log t)*t, since T>=log t>\log t>. Let us now compare > with the computation time> for an optimal solution to the global approximation problem. In fact, it suffices to compare with an optimal solution for the static version: in the dynamic case, we consider the last static global approximation problem. Denote by >> the computation time at each node > for a fixed optimal solution, so that =\G>t>>. If is the size of the dag, then we claim that <\equation> t\t\2*s*t. For simplicity, we will assume that >\p>> whenever > is a descendant from >. This is no longer the case if we apply the optimization from the end of section, but the reasoning can probably be adapted to that case. Assume for contradiction that () does not hold. Now consider the moment during the execution at which we first call with a maximal cost > for some node >. At that point, the current cost >> of each of the descendants> of> is >=2=T>/2>. When such a > is a leaf, it follows that >=t>/2\t/(2*s)\t\t>>. By structural induction over the descendants > of >, it follows that >\t>> and the best available ( optimal) approximation >> ( >>) for >> satisfies >>\r>>\\>>. In particular >>\\>>. On the other hand, the first call of with a maximal cost > was necessarily triggered by , whence >>\\>>. This contradiction proves our claim. Up to a constant factor, the bound () is sharp. Indeed, consider the case of a multiplication *\*x> of numbers which are all zero. When gradually increasing the precisions for the computation of ,\,x>, it can happen that one of the > produces bounding intervals> whose radii quickly converge to zero, contrary to each of the other>. In that case, the time spent on improving each of the >(i>) is a waste, whence we loose a factor with respect to the optimal solution. On the other hand, without additional knowledge about the functions>, it is impossible to design a deterministic procedure for choosing the most efficient index. In this sense, our current solution is still optimal. However, under additional monotonicity hypotheses on the cost functions, efficient indices can be found, by taking into account the ``cost per digit''. error estimates> Although the approach from the previous section has the advantage of never being extremely bad and rather easy to implement on top of an existing layer for interval arithmetic, there are even simple cases in which the factor in the bound () is not necessary: in the truncated power series evaluation +a/2+\+a/2> with \|\1> for all , the computation of a n>>-approximation of induces the computation of -bit approximations of each of the >. If s>, this means that we spend a time n*s> instead of n>. In order to remedy to this problem, we suggest to improve the balanced estimate technique from and cleverly recombine it with the current approach. In this section, we will briefly sketch how this could be done. The results are based on joint ideas with Kreinovich>, which we plan to work out in a forthcoming paper. Let us start by isolating those situations in which error estimates should be efficient. Consider a labeled dag so that >> admits an initial interval approximation >\x>> for each \G>. Assume also that for each node > and each child> of >, we have an interval ,i>> with x>/\ x>)(\>,\,\\|>>)\\,i>>. If ,i>\|\\>, then we say that (together with the >>) is a dag. If, in addition, we have ,i>\\(c,i>>,\*\|c,i>>\|)> for some \\1> and all ,i>, then we say that is >-rigid>. A typical obstruction to the Lipschitz property occurs in dags like >. Similarly, a dag like0> is typically Lipschitz, but not rigid. Given a Lipschitz dag, a variant of automatic differentiation provides us with bounds for the error in >> in terms of the errors in the >>, where > ranges over the leaves below >. If is >-rigid, and especially when \2W>>, then these bounds actually become very sharp. For instance, given a rooted Lipschitz dag and a challenge > at the root >, one may compute a sufficient precision for obtaining an >-approximation of >> as follows. Let >> be the error in >> at each leaf>, when computing with precision . We have >=r>*2l>> for some>> which depends on >. Then we recursively estimate the error >> at each node > by <\equation*> \>=\||\>,1>*\>+\+\||\>,\|\\|>*\\|>>. This provides us with a bound of the form >=r>*2l>> for the error at the root >. We may thus take log r>/\\>. The approach can be further improved using similar ideas as in the implementation of addition in section . Instead of working with a fixed precision , a better idea is to compute the contribution >=\\>/\\>> of the error >> at each leaf to the error >> at >. This problem is dual to the problem of automatic differentiation, since it requires us to look at the opposite dag >> of , which is obtained by inverting the direction of the edges. Indeed, if ,\,\]>> denote all parents of a node >, and if > is the ,j>>-th child of > for each , then we take <\equation*> \>=\||\>,i,1>>*\>+\+\||\>]>,i],1>>*\]>>. Together with the initial condition >=1>, this allows us to compute >> for all leafs >. In order to compute >> with error >\\>, we may now balance > over the leaves > according to the >>. More precisely, we compute an >=\/(p*\>)>-approximation of each >>, where is the number of leaves, and recompute all other nodes using interval arithmetic. As an additional optimization, one may try to balance according to the computational complexities of the >>. The above strategy is a bit trickier to implement in an incremental way. Indeed, in the dynamic global approximation problem, we ask for >-approximations at different nodes > and, since good previous approximations may already be present, it is not always necessary to compute the complete dag below >. A solution to this problem is to keep track of the ``creation date'' of each node > and to compute the >> from the top down to the leaves, while first considering nodes with the latest creation date ( by means of a heap). Whenever the computed >> is so small that the current error >=r>>> at > contributes only marginally to the error >\\> at the top ( >*\>\\*2W>>), then it is not necessary to consider the descendants of >. <\bibliography|bib|abbrv|all.bib> <\bib-list|10> G.Alefeld and J.Herzberger. . Academic Press, 1983. J.Blanck. General purpose exact real arithmetic. Technical Report CSR 21-200, Luleå University of Technology, Sweden, 2002. . J.Blanck, V.Brattka, and P.Hertling, editors. , volume 2064 of Springer, 2001. A.Edalat and P.Sünderhauf. A domain-theoretic approach to real number computation. , 210:73--98, 1998. A.Gaganov. Computational complexity of the range of the polynomial in several variables. , pages 418--425, 1985. T.Granlund et al. GMP, the GNU multiple precision arithmetic library. , 1991--2006. M.Grimmer, K.Petras, and N.Revol. Multiple precision interval packages: Comparing different approaches. Technical Report RR 2003-32, LIP, École Normale Supérieure de Lyon, 2003. G.Hanrot, V.Lefèvre, K.Ryde, and P.Zimmermann. MPFR, a C library for multiple-precision floating-point computations with exact rounding. , 2000--2006. V.Kreinovich. For interval computations, if absolute accuracy is NP-hard, then so is relative accuracy+optimization. Technical Report UTEP-CS-99-45, UTEP-CS, 1999. V.Kreinovich and S.Rump. Towards optimal use of multi-precision arithmetic: a remark. Technical Report UTEP-CS-06-01, UTEP-CS, 2006. B.Lambov. The RealLib project. , 2001--2006. V.Ménissier-Morain. Arbitrary precision real arithmetic: design and algorithms. N.Müller. iRRAM, exact arithmetic in C++. , 2000--2006. R.O'Connor. A monadic, functional implementation of real numbers. Technical report, Institute for Computing and Information Science, Radboud University Nijmegen, 2005. N.Revol. MPFI, a multiple precision interval arithmetic library. , 2001--2006. S.Rump. Fast and parallel inteval arithmetic. , 39(3):534--554, 1999. A.Turing. On computable numbers, with an application to the Entscheidungsproblem. , 2(42):230--265, 1936. J.vander Hoeven. GMPX, a C-extension library for gmp. , 1999. No longer maintained. J.vander Hoeven. Relax, but don't be too lazy. , 34:479--542, 2002. J.vander Hoeven. Computations with effective real numbers. , 351:52--60, 2006. J.vander Hoeven et al. Mmxlib: the standard library for Mathemagix, 2002--2006. . K.Weihrauch. . Springer-Verlag, Berlin/Heidelberg, 2000. <\initial> <\collection> <\references> <\collection> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > <\auxiliary> <\collection> <\associate|bib> Tur36 Wei00 ES98 MM96 BBH01 Bl02 IRRAM REALLIB Con05 vdH:effreal AH83 Rump99 GPR03 MPFI vdH:gmpx vdH:mml vdH:effreal AH83 Rump99 Bl02 GMP MPFR vdH:effreal Ga85 Krei99 vdH:effreal MPFR vdH:relax Wei00 Wei00 vdH:effreal vdH:effreal vdH:effreal KR06 vdH:effreal vdH:effreal <\associate|figure> <\tuple|normal> Example of a dag with |2> roots. The dag has |size> |8> ( the total number of nodes) and |depth> 3 ( the longest path from a root to a leaf). The |weight> of the dag corresponds to the sum of the sizes of the trees obtained by ``copying'' each of the roots. In our example, the weight is |13>. Finally, the |ancestrality> of the dag is defined to be the maximum number of ancestors of a leaf. In our example, the ancestrality is |5>. > A badly balanced tree.|> <\associate|toc> |math-font-series||1Introduction> |.>>>>|> |math-font-series||2Interval arithmetic> |.>>>>|> |2.1Representation issues |.>>>>|> > |Endpoint representation |.>>>>|> > |Ball representation |.>>>>|> > |Hybrid representation |.>>>>|> > |2.2Precision changes and semantics |.>>>>|> > |2.3Complex numbers and other template types |.>>>>|> > |math-font-series||3Effective real numbers> |.>>>>|> |3.1Definitions and theoretical properties |.>>>>|> > |3.2Dag models |.>>>>|> > |3.3The global approximation problems |.>>>>|> > |3.4|A priori> error estimates |.>>>>|> > |3.5|A posteriori> error estimates and relaxed evaluation |.>>>>|> > |math-font-series||4Effective real numbers in |Mmxlib>> |.>>>>|> |4.1The classes |language||real> and |language||real_rep> |.>>>>|> > |4.2Examples of derived classes of |language||real_rep> |.>>>>|> > |math-font-series||5Complexity analysis> |.>>>>|> |5.1Total complexity versus final complexity |.>>>>|> > |5.2Final complexity versus optimal complexity |.>>>>|> > |math-font-series||6Back to |a priori> error estimates> |.>>>>|> |6.1Rigid dags |.>>>>|> > |6.2Backward error bounds |.>>>>|> > |math-font-series||Bibliography> |.>>>>|>