Claude-Pierre Jeannerod

A radix-independent error analysis of the Cornea-Harrison-Tang method

Assuming floating-point arithmetic with a fused multiply-add operation and rounding to nearest, the Cornea-Harrison-Tang method aims to evaluate expressions of the form ab + cd with high relative accuracy. In this paper we provide a rounding error analysis of this method, which unlike previous studies is not restricted to binary floating-point arithmetic but holds for any radix β. We show first that an asymptotically optimal bound on the relative error of this method is 2βu+2u 2 β-2u 2 = 2u+ 2 β u 2 +O(u 3 ), where u = 1 2 β 1-p is the unit roundoff in radix β and precision p. Then we show that the possibility of removing the O(u 2 ) term from this bound is governed by the radix parity and the tie-breaking strategy used for rounding: if β is odd or rounding is to nearest even, then the simpler bound 2u is obtained, while if β is even and rounding is to nearest away, then there exist floating-point inputs a, b, c, d that lead to a relative error larger than 2u + 2 β u 2 -4u 3 . All these results hold provided underflows and overflows do not occur and under some mild assumptions on β and p satisfied by IEEE 754-2008 formats.

Introduction

Given four floating-point numbers a, b, c, d the Cornea-Harrison-Tang method [1, p. 273] aims to evaluate x = ab + cd efficiently and accurately using the fused multiply-add operation. Writing RN to denote rounding to nearest, this method can be described as follows:

algorithm CHT(a, b, c, d) p 1 := RN(ab); p 2 := RN(cd); e 1 := RN(ab -p 1 ); e 2 := RN(cd -p 2 ); // these two operations are exact. r := RN(p 1 + p 2 ); e := RN(e 1 + e 2 ); x := RN(r + e); return x

One key feature of this algorithm is its use of the fused multiply-add operation to compute the rounding errors of the two multiplications exactly in the absence of underflow and overflow, so that e 1 = ab -p 1 and e 2 = cd -p 2 . The rounded sum e of these error terms is then added to the (possibly highly inaccurate) rounded sum r of the two products in order to obtain an approximation x having a tiny relative error. Another attractive feature of this method is its symmetry, which ensures that ab + cd and cd + ab are approximated by the same quantity and thus makes it straightforward to provide implementations of complex floating-point multiplication that preserve commutativity.

The accuracy of the CHT algorithm has been studied extensively in radix 2: assuming p-bit floating-point numbers and an unbounded exponent range, Cornea, Harrison, and Tang [1, pp. 273-275] showed that the relative error | x -x|/|x| is always in O(u) with u = 2 -p the unit roundoff; this result was refined recently by Muller [START_REF] Muller | On the error of computing ab + cd using Cornea, Harrison and Tang's method[END_REF], who derived the upper bound 2u + 7u 2 + 6u 3 and found that | x -x|/|x| can be as large as 2u -7u 2 + O(u 3 ) for some values of a, b, c, d. In other words, in radix 2 the relative error of algorithm CHT is bounded by 2u + O(u 2 ), and this bound is asymptotically optimal in the sense that there are inputs for which the ratio error/(error bound) tends to one as u tends to zero.

These results raise two questions, however, which we answer in this paper:

1. Does the bound 2u + O(u 2 ) hold beyond radix β = 2, that is, for β > 2 and u equal to 1 2 β 1-p ? 2. Is it possible to remove the quadratic term O(u 2 ) and thus to bound the relative error simply by 2u?

The first question is natural since the IEEE 754-2008 standard [START_REF]IEEE Standard for Floating-Point Arithmetic, IEEE Standard[END_REF] specifies floating-point arithmetic not only for radix 2, but also for radix 10. Furthermore, although the techniques developed in [START_REF] Muller | On the error of computing ab + cd using Cornea, Harrison and Tang's method[END_REF] for β = 2 extend to β > 2, the resulting bound on | x-x|/|x| would be larger than 3β+4 β+4 u and thus larger than 2.2u when β 6.

The second question is motivated by the rounding error analysis of another method for evaluating ab + cd with a fused multiply-add, namely Kahan's algorithm [2, p. 60]. Kahan's algorithm computes only one product and its error term (say, p 1 and e 1 ), then handles the other product directly by using the fused multiply-add operation r = RN(p 1 + cd), and finally returns RN(r + e 1 ). This approach thus saves three floating-point operations compared with algorithm CHT and, as shown in [START_REF] Jeannerod | Further analysis of Kahan's algorithm for the accurate computation of 2 × 2 determinants[END_REF], it admits 2u as an asymptotically optimal bound on its relative error. However, this comes at the price of a lack of symmetry, and it is therefore important to understand whether this simple O(u 2 )-free bound 2u can still be achieved in the symmetrized version provided by algorithm CHT.

Main results. Our first contribution is to answer the first question above positively, by proving that the bound 2u + O(u 2 ) holds for p 6; this extends in particular the result of [START_REF] Muller | On the error of computing ab + cd using Cornea, Harrison and Tang's method[END_REF] to the practical case β = 10. Our second contribution is to show that, perhaps surprisingly, the answer to the second question depends on the parity of β and the way RN breaks ties: in some cases (say, when β is odd or ties are rounded to even), the bound 2u + O(u 2 ) can be replaced by 2u, while in other cases the O(u 2 ) term cannot be removed.

More precisely, we shall work under the following customary assumptions (all of which are implicitly or explicitly used for the analyses in radix 2 given in [START_REF] Cornea | Scientific Computing on Itanium R -based Systems[END_REF][START_REF] Muller | On the error of computing ab + cd using Cornea, Harrison and Tang's method[END_REF]), and establish Theorems 1 and 2 below. Here and hereafter a, b, c, d are taken from a set F of finite floating-point numbers in base β and precision p. We assume that β 2 and p 2, and that the exponent range of F is unbounded, so

F = {0} ∪ {S • β e : S, e ∈ Z, β p-1 |S| < β p }.
We also assume that the exact result of every operation on some element(s) of F is rounded back to F using a round-to-nearest function RN satisfying the following properties: for all t ∈ R,

|RN(t) -t| = min s∈F |s -t|; RN(-t) = -RN(t); RN(β i t) = β i RN(t) for all i ∈ Z.
The last two properties say that the way of breaking ties is independent of the sign and magnitude of the number being rounded. This assumption is in particular satisfied by roundTiesToEven and roundTiesToAway, the two specifications of rounding to nearest given in the IEEE 754-2008 standard: when t is a midpoint, that is, a number halfway between two consecutive elements of F, then roundTiesToEven requires that the significand S of RN(t) is an even integer, while roundTiesToAway requires that |S| is maximal. For example, writing

u = 1 2 β 1-p
for the unit roundoff associated with F and RN, the midpoint 1 + u is rounded down to 1 ∈ F by roundTiesToEven, and up to 1+2u ∈ F by roundTiesToAway.

We can now state our main results more formally:

Theorem 1. If β p-1 24 then the value x computed by algorithm CHT satisfies

x = x(1 + ), | | 2u if β is odd or RN(1 + u) = 1, 2βu+2u 2 β-2u 2 otherwise.
Furthermore, these bounds on the relative error | | are asymptotically optimal.

Since 2βu+2u 2 β-2u 2 = 2u + 2 β u 2 + O(u 3
), this first result shows that the relative error of algorithm CHT is always bounded by 2u + O(u 2 ) and that the leading constant 2u is best possible as u tends to zero. The next result shows that when β is even and RN is so that the midpoint 1 + u is rounded up to 1 + 2u, then the term O(u 2 ) cannot, in general, be removed. Theorem 2. Assume β is even and RN(1 + u) = 1 + 2u. If β = 2 and 2 p + 1 is not a Fermat prime, or if β = 2, then there exist a, b, c, d in F for which the value x computed by algorithm CHT satisfies

x = x(1 + ), | | > 2u + 2 β u 2 -4u 3 .
Consequences for IEEE arithmetic. When β = 2 Theorem 2 excludes values of p such that 2 p + 1 is a Fermat prime, that is, a prime number of the form 2 2 q +1 with q ∈ N. However, this is not a restriction in practice, since 2 p +1 is known to be composite for any of the binary interchange formats specified by the IEEE 754-2008 standard; see [START_REF] Jeannerod | On relative errors of floating-point operations: optimal bounds and applications[END_REF]. Similarly, it is easily checked that the assumption β p-1 24 used in Theorem 1 is satisfied for all formats. Third, roundTiesToEven and roundTiesToAway imply RN(1 + u) = 1 and RN(1 + u) = 1 + 2u, respectively. Therefore, in the specific context of IEEE arithmetic Theorems 1 and 2 lead to the following conclusion:

Corollary 1. Assume floating-point arithmetic as specified by the IEEE 754-2008 standard, with radix β and unit roundoff u. Then, in the absence of underflow and overflow, algorithm CHT has a relative error bounded by 2u when RN is roundTiesToEven, and by 2βu+2u 2 β-2u 2 = 2u + O(u 2 ) when RN is roundTiesTo-Away. Furthermore, for roundTiesToAway, the O(u 2 ) term cannot be removed, since there exist floating-point numbers a, b, c, d leading to a relative error larger than 2u + 2 β u 2 -4u 3 .

Outline of the paper and additional background. The rest of this paper is devoted to the proof of Theorems 1 and 2. We begin in Section 2 by presenting the three main tools used to establish the upper bounds in Theorem 1. Then we give in Section 3 an outline of our proof of that theorem, showing that it mainly consists of analyzing separately several cases which depend on the features of the exact or rounded values of the products ab and cd; overall, this case analysis leads to about ten different error bounds, whose detailed proofs are postponed to the appendix for readability. Finally, the lower bound given in Theorem 2 is established independently in Section 4.

A useful tool for our analyses will be the unit in the first place function [START_REF] Rump | Accurate floatingpoint summation, Part I: Faithful rounding[END_REF], denoted by ufp and defined for t ∈ R by

ufp(t) = 0 if t = 0, β log β |t| if t = 0.
By definition of F, RN, and u we have the classical relations

|RN(t) -t| u ufp(t) u min{|t|, |RN(t)|} for all t ∈ R,
which lead in particular to the standard models

RN(t) = t(1 + ) with | | u and RN(t) = t/(1 + ) with | | u.
For | |, the upper bound u can in fact be replaced by the slightly smaller quantity

u 1 := u 1 + u ,
giving the following refined model:

RN(t) = t(1 + ), | | u 1 . (1)
This refined bound appears for example in [7, p. 232] and its attainability was noted in [START_REF] Jeannerod | On relative errors of floating-point operations: optimal bounds and applications[END_REF] along with the attainability of the bound | | u.

When using | | u 1 instead of | | u, the bound 2u + 7u 2 + O(u 3
) obtained in [START_REF] Muller | On the error of computing ab + cd using Cornea, Harrison and Tang's method[END_REF] immediately becomes 2u + 5u 2 + O(u 3 ). However, this sharper bound is not enough for our purposes, since it assumes β = 2 and has a O(u 2 ) term regardless of the tie-breaking strategy.

2 Tools for the proof of Theorem 1

Inequalities resulting from the refined model

Applying the refined model ( 1) to the operations in algorithm CHT, we deduce that there exist rational numbers 1 , 2 , 3 , 4 , 5 such that

p 1 = ab(1 + 1 ), | 1 | u 1 , p 2 = cd(1 + 2 ), | 2 | u 1 , e = (e 1 + e 2 )(1 + 3 ), | 3 | u 1 , r = (p 1 + p 2 )(1 + 4 ), | 4 | u 1 , x = (r + e)(1 + 5 ), | 5 | u 1 .
Recalling that e 1 = ab -p 1 and e 2 = cd -p 2 , we have

x = p 1 + p 2 + e 1 + e 2
and, therefore,

x = x(1 + 4 )(1 + 5 ) + (e 1 + e 2 )( 3 -4 )(1 + 5 ). (2)
On the other hand, the definition of 1 and 2 implies that e 1 = -1 ab and e 2 = -2 cd.

Hence, for x nonzero and with

K = |ab| + |cd| |ab + cd| , (3)
we arrive at the following inequalities:

| x -x| |x| | 4 + 5 + 4 5 | + max{| 1 |, | 2 |} • | 3 -4 |(1 + 5 ) • K | 4 + 5 + 4 5 | + 2u 2 1 (1 + u 1 ) • K. (4)

Range constraints resulting from large relative errors

In addition to the usual unit roundoff u and to the quantity u 1 = u 1+u , the following generalization will prove very useful in the sequel:

u k := u 1 + ku , k ∈ R 0 .
In particular, for fixed k we see that u k = u -ku 2 + O(u 3 ) as u tends to zero.

Having defined u k , we can now state the following two properties, which indicate the range constraints implied by large enough relative errors. Property 1 says that if rounding a real number t yields a relative error that is larger than u k , then |t| is necessarily close enough to its ufp. This immediate property is then refined by Property 2, which exploits further the sign of the relative error in order to confine |t| to unions of about k/2 intervals of width O(u 2 ) • ufp(t).

Property 1. Let k ∈ R >1 . Then for t ∈ R =0 we have the following implication:

|RN(t) -t| |t| > u k ⇒ 1 < |t| ufp(t) < 1 + ku.
Proof. The lower bound on |t|/ufp(t) follows from the fact that t ∈ F, and the upper bound follows from u k |t| < |RN(t) -t| u ufp(t).

Property 2. Given k ∈ R >1 , let = (k -1)/2 and define the half-open intervals

I j = 1 + (2j -1)u, 1 + 2ju 1 + u k , j = 1, . . . , , and
I j = 1 + 2ju 1 -u k , 1 + (2j + 1)u , j = 0, . . . , -1.
Then for t ∈ R =0 we have the following implications:

(i) RN(t) -t t > u k ⇒ |t| ufp(t) ∈ I 1 ∪ I 2 ∪ • • • ∪ I ; (ii) RN(t) -t t < -u k ⇒ |t| ufp(t) ∈ I 0 ∪ I 1 ∪ • • • ∪ I -1 .
Proof. We can assume t > 0 and ufp(t) = 1, so that 1 t < β and |RN(t) -t| u. To prove (i), note first that since RN(t) is in F and larger than t, it has the form RN(t) = 1 + 2ju for some integer j 1. The assumption

RN(t)-t t > u k is thus equivalent to t < 1 + 2ju 1 + u k . ( 5a
)
In addition, |RN(t) -t| u implies RN(t) -t u, that is,

1 + (2j -1)u t. ( 5b
)
Hence the expression for interval I j follows from the inequalities [START_REF] Jeannerod | Further analysis of Kahan's algorithm for the accurate computation of 2 × 2 determinants[END_REF]. Recalling from Property 1 that t < 1 + ku, we deduce from (5b) that the integer j and the real number k satisfy 2j -1 < k, that is, j (k -1)/2 = . This concludes the proof of (i).

Let us now prove (ii). From 1 RN(t) ∈ F it follows that RN(t) = 1 + 2ju for some integer j 0. The assumption

RN(t)-t t < -u k is then equivalent to 1 + 2ju 1 -u k < t. (6a)
On the other hand, |RN(t) -t| u implies -u RN(t) -t, that is,

t 1 + (2j + 1)u. ( 6b
)
From the inequalities (6) we deduce the definition of interval I j . Finally, using again Property 1 we have t < 1 + ku, which together with (6a) and

u k = u/(1 + ku) leads to 1 + 2ju < (1 + ku)(1 -u k ) = 1 + (k -1)u.
The latter inequality is equivalent to j -1, which concludes the proof.

In practice, when analyzing the CHT algorithm we shall avoid using the unwieldy rational functions involved in the right endpoint of I j and the left endpoint of I j . Instead, it will be enough to consider the following simpler intervals I j and I j , defined by degree-2 polynomials in u: for j = 1, . . . , ,

I j := α j , α j + j,k , α j := 1 + (2j -1)u, j,k := (k -2j + 1)u 2
and, for j = 0, . . . , -1,

I j := α j -j,k , α j , α j := 1 + (2j + 1)u, j,k := (k -2j -1)u 2 .
Since is defined in Property 2 as = (k -1)/2 , it is easily checked that

I j ⊂ I j , I j-1 ⊂ I j-1 , j = 1, . . . , .

Properties of floating-point products

Algorithm CHT is built upon the fact that for an unbounded exponent range, the rounding error of the product of two elements of F is always itself in F. The next two properties show that if this error is nonzero then its ulp cannot be too small. Those properties (which are straightforward extensions to radix β of those used in [START_REF] Muller | On the error of computing ab + cd using Cornea, Harrison and Tang's method[END_REF] for radix two) will be useful when dealing with the case where p 1 and -p 2 are so close to each other that p 1 + p 2 is computed exactly.

min{β i-2p+2 , β i-2p+1 } = β i-2p+1 .
The property above can be refined in the sense that either the error is an integer multiple of a larger quantity, or the product admits a smaller upper bound:

Property 4. Let i ∈ Z and a, b ∈ F be such that β i |ab| < β i+1 . Then • either ab -RN(ab) is an integer multiple of β i-2p+2 , • or |ab| (1 -2u β ) 2 β i+1 . Proof.
Using the same notation and reasoning as in the proof of Property 3, if i = e a +e b then ab-RN(ab) is an integer multiple of

β i-2p+2 . Else i = e a +e b +1, but then |ab| = |AB|β 1-2p+i (β p -1) 2 β 1-2p+i = (1 -2u β ) 2 β i+1 .
3 Proof outline for Theorem 1

The first step of the proof of the bounds in Theorem 1 consists of restricting the input set using symmetry arguments. Since CHT(a, b, c, d) = CHT(c, d, a, b), we can exchange ab and cd to ensure |cd| |ab|. On the other hand, using RN(-t) = -RN(t) gives CHT(-a, b, -c, d) = -CHT(a, b, c, d), so that we can also restrict further to ab 0 and eventually assume that |cd| ab.

As a second step, we consider the situation where either ab or cd or x is zero. In this case, by propagating the equality RN(0) = 0 within the CHT algorithm, we see that x = RN(x) and, recalling the refined model in (1),

x = x(1 + ), | | u 1 . (7)
The third and main step of the proof is the analysis of the remaining cases, namely when a, b, c, d are such that ab > 0 and cd = 0 and -ab < cd ab.

To derive the bounds in Theorem 1 for inputs as in [START_REF] Muller | On the error of computing ab + cd using Cornea, Harrison and Tang's method[END_REF], we shall analyze separately several subcases and, using the tools described in Section 2, obtain the nine bounds from ( 9) to (17) below. The rest of this section only gives an overview of those subcases and the associated bounds, the detailed proofs being deferred to Appendix.

Analysis when ab and cd have the same sign

In this case, the number K = (|ab| + |cd|)/|ab + cd| introduced in (3) satisfies K = 1. If min{ 4 , 5 } u 3 then, using (4), we easily obtain

| x -x| |x| < 2u -u 2 + 4u 3 ; (9)
else, using further Property 1, we deduce that

| x -x| |x| < 2u -u 2 + u 4 for β p-1 3. ( 10
)
3.2 Analysis when ab and cd have opposite signs

3.2.1 When 1 2 p 1 |p 2 |
In this case the sum p 1 + p 2 is computed exactly thanks to Sterbenz's theorem [START_REF] Sterbenz | Floating-Point Computation[END_REF], and a direct consequence of this will be that if K 1/u 1 then the bound in [START_REF] Sterbenz | Floating-Point Computation[END_REF] still applies here. If K > 1/u 1 , then ufp(cd) turns out to be either ufp(ab) or β -1 ufp(ab), and we can handle those two cases separately: in the first case, by applying Property 3 to both ab and cd, we can show that e 1 + e 2 is a floating-point number and then deduce that

| x -x| |x| u 1 ; (11)
in the second case, applying Property 3 to ab and Property 4 to cd leads either to the same bound as in (11) or to the bound

| x -x| |x| 3 2 u. (12)
3.2.2 When 1 2 p 1 > |p 2 | In this case, noting that K is at most about 3, we shall consider separately four cases defined by the pair ( 4 , 5 ). The first two subcases can be handled using this bound on K together with (4): if 4 and 5 have opposite signs, then

| x -x| |x| u + 7u 2 for β p-1 4; (13) if min{| 4 |, | 5 |} u 7 , then | x -x| |x| < 2u for β p-1 24. ( 14
)
The remaining two subcases, which correspond to min{| 4 |, | 5 |} > u 7 with 4 and 5 either both positive or both negative, turn out to be more involved and their analysis relies on either part (i) or part (ii) of Property 2 with k = 7. When 4 > u 7 and 5 > u 7 , we show the following:

• if β is odd or RN(1 + u) = 1, then | x -x| |x| < 2u for β p-1 10; (15) • else | x -x| |x| 2βu + 2u 2 β -2u 2 for β p-1 10. ( 16
)
When 4 < -u 7 and 5 < -u 7 , we show that

| x -x| |x| < 2u -u 2 for β p-1 10. ( 17
)
The two error bounds in Theorem 1 then follow directly from the intermediate bounds ( 7) and (9-17) shown above: these ten bounds require no more than β p-1 24; furthermore, all of them are less than 2u, except the one in ( 16), which has the form 2u + 2 β u 2 + O(u 3 ). Finally, we know from [4, Corollary 4.2] that if β p-1 11, then there exist a, b, c, d ∈ F for which the CHT algorithm returns x such that

| x -x| |x| > 2u -8u 1.5 -6u 2 ;
this proves the asymptotic optimality of the two error bounds in Theorem 1.

Proof of Theorem 2

The assumption on β and p implies that there exist a, b ∈ F such that Define further

ab = 1 + u; (18)
c = u + 2u 2 and d = -1 + β -1 β • 2u. Recalling that u = 1 2 β 1-p , we have c = C • β 1-2p with C = 1 2 β p + β 2 and d = -D • β -p with D = β p -β + 1.
Since C and D are integers such that β p-1 C, D < β p , we deduce that c and d are in F. In addition,

cd = -u + 2 β u 2 -4(1 - 1 β )u 3 , (19)
which implies

u < |cd| < u + 2 β u 2 = u + u ufp(u)
and thus p 2 = -u.

which gives

| x -x| |x| u 1 + u 3 + u 1 u 3 + 2u 2 1 + 2u 3 1 < 2u -u 2 + 4u 3
and thus proves [START_REF] Rump | Accurate floatingpoint summation, Part I: Faithful rounding[END_REF].

A.2 Case where 4 > u 3 and 5 > u 3

Since both ab and cd are positive and since the exponent range of F is unbounded, we have p 1 + p 2 > 0. Thus, applying Property 1 with k = 3 gives

β i < p 1 + p 2 < (1 + 3u)β i , β i = ufp(p 1 + p 2 ).
Since 4 is positive, rounding to nearest coincides here with rounding up, and the rounded sum r = RN(p 1 + p 2 ) must be

r = (1 + 2u)β i .
Let us now bound |e|. We have |e|

(1 + u 1 ) |e 1 | + |e 2 | (u 1 + u 2 1 )
x. Furthermore, it follows from r = ab(1

+ 1 )(1 + 4 ) + cd(1 + 2 )(1 + 4 ) that (1 -u 1 ) 2 x r (1 + u 1 ) 2 x. ( 20
)
Using the lower bound in (20) thus leads to

|e| u 1 + u 2 1 (1 -u 1 ) 2 r = (1 + 2u)ur < 2uβ i for β p-1 3.
Consequently,

β i < r + e < (1 + 4u)β i β i+1 .
Now, the assumption 5 > u 3 implies that when rounding r + e to nearest then rounding up occurs and, by Property 1, that r + e must be less than (1 + 3u)β i . In other words, (1 + u)β i r + e < (1 + 2u)β i and x = (1 + 2u)β i .

Therefore, x = r and, using (20), we conclude that

| x -x| |x| = |r -x| |x| 2u 1 + u 2 1 < 2u -u 2 + u 4 for β p-1 3. ( 21
)
This proves [START_REF] Sterbenz | Floating-Point Computation[END_REF] and concludes the analysis of the case where the two products ab and cd have the same sign.

B Analysis when the two products have opposite signs

When ab and cd have opposite signs, the assumption in ( 8) can be rewritten

-ab < cd < 0 < ab. ( 22
)
For K as in (3), this implies

K = ab + |cd| ab -|cd| (23)
and, rounding being monotonic,

-p 1 p 2 < 0 < p 1 .
Note that neither p 1 nor p 2 can be zero (because the exponent range of F is unbounded) and that p 1 + p 2 0.

B.1 When

1 2 p 1 |p 2 |
In this case, the two floating-point numbers p 1 and -p 2 satisfy 1 2 p 1 -p 2 p 1 , so that for any radix β the sum p 1 + p 2 is computed exactly by Sterbenz's theorem [10, p. 138] (see also [2, p. 45]). Hence 4 = 0 and, using (2),

| x -x| |x| | 5 | + |(e 1 + e 2 ) 3 | |x| (1 + 5 ) ( 24
)
u 1 + u 2 1 (1 + u 1 )K.
If K 1/u 1 then we deduce immediately from the latter bound that the relative error on x is bounded by 2u 1 + u 2 1 and thus as in (21). Hence the rest of this section is devoted to handling the case

K > 1 u 1 = 1 u + 1. ( 25
)
(This case of a huge value of K does occur, for example when a = 1 + 2u, b = 1 -u, c = 1, and d = -1.) From ( 22), (23), and (25) we deduce that

1 β ab 1 1 + 2u ab < |cd| < ab. Consequently, if ab ∈ [β i , β i+1 ) with i ∈ Z, then |cd| is either in [β i-1 , β i ) or in [β i , β i+1
). This yields two sub-cases which we handle separately.

B.1.1 Case β i |cd| < ab < β i+1

In this case, ufp(ab) = ufp(cd) = β i and, using |e 1 | u ufp(ab) and |e 2 | u ufp(cd), we deduce that

|e 1 + e 2 | 2uβ i = β i-p+1 .
On the other hand, Property 3 implies the existence of integers

E 1 , E 2 such that e 1 = E 1 • β i-2p+1 , e 2 = E 2 • β i-2p+1
.

Hence We now have ufp(ab) = β i and ufp(cd) = β i-1 , so that

e 1 + e 2 = (E 1 + E 2 )β i-
|e 1 + e 2 | uβ i + uβ i-1 < 2uβ i . (26)
Property 3 still gives

e 1 = E 1 • β i-2p+1
for some integer E 1 , and by applying Property 4 to the product cd we have either

e 2 = E 2 • β i-2p+1
for some integer E 2 , or

|cd| 1 - 2u β 2 β i .
We handle these two situations independently as follows.

If

e 2 = E 2 • β i-2p+1 , then |e 1 + e 2 | = |E 1 + E 2 |β i-2p+1 . Using (26) gives |E 1 + E 2 | < β p ,
from which we deduce e 1 + e 2 ∈ F, that is, 3 = 0. We conclude as in Section B.1.

1 that | x -x|/|x| u 1 . Assume now that |cd| 1 -2u β 2 β i . Then |x| = ab -|cd| β i -1 - 2u β 2 β i .
On the other hand, the strict inequality in (26) implies ufp(e 1 + e 2 ) β i-p and, since RN(e 1 + e 2 ) = (e 1 + e 2 )(1 + 3 ), we obtain

|(e 1 + e 2 ) 3 | uβ i-p . Applying (24) thus gives | x -x|/|x| u 1 + (1 + u 1 )ϕ with ϕ := |(e 1 + e 2 ) 3 | |x| uβ -p 4u β -4u 2 β 2 . B.2.
3 Case where 4 > u 7 and 5 > u 7

In this case, we will derive the bounds in ( 15) and ( 16), depending on radix parity and the tie-breaking strategy of rounding to nearest. Defining the condition (C) :

β is odd or RN(1 + u) = 1,
our goal in this section is thus to show that for β p-1 10,

| x -x| |x|    2u -η for some η > 0 if (C) holds, 2βu+2u 2 β-2u 2 = 2u + 2 β u 2 + O(u 3 ) otherwise. (29)
To establish (29), we shall apply Property 2 in order to obtain suitable ranges for the exact sum p 1 + p 2 from which we can then deduce some values for x together with some ranges for x, and eventually some bounds on | x -x|/|x|.

Preliminaries. Since 4 is nonzero, p 1 + p 2 is not in F (and thus nonzero as well), so there exists an integer i such that β i < p 1 + p 2 < β i+1 . In order to simplify the expressions used in the sequel, we shall assume that i = 0 (which is possible up to a scaling by an integer power of the base β and because the exponent range of F is unbounded). Therefore,

1 < p 1 + p 2 < β.
Since 4 > u 7 , applying part (i) of Property 2 with k = 7 gives

p 1 + p 2 ∈ [1 + u, 1 + u + 6u 2 ) =: I 1 ∪ [1 + 3u, 1 + 3u + 4u 2 ) =: I 2 ∪ [1 + 5u, 1 + 5u + 2u 2 ) =: I 3 . Furthermore, 4 > 0 implies that r = RN(p 1 + p 2 ) satisfies r > p 1 + p 2 .
From this strict inequality and the range of p 1 + p 2 shown above, we deduce that r ∈ {1 + 2u, 1 + 4u, 1 + 6u}.

Furthermore, the right endpoint of I 3 leads to the following bound on |e|: since

1 2 p 1 > |p 2 | by assumption, we have p 1 + |p 2 | < 3(p 1 + p 2 ) and thus |e| (1 + u 1 ) |e 1 | + |e 2 | u(1 + u 1 ) p 1 + |p 2 | < 3u(1 + u 1 )(p 1 + p 2 ) < 4u when β p-1 10. (31)
From ( 30) and (31) we deduce that 1 -2u < r + e < 1 + 10u. Since β p-1 10, this implies r + e ∈ [β -1 , 1) ∪ [1, β). On the other hand, if 1 -2u < r + e 1 then the relative error | 5 | is at most uβ -1 /(1 -2u), which contradicts the assumption 5 > u 7 . Thus, overall, we must have

1 < r + e < β,
and, using 5 > 0, the rounded value x = RN(r + e) must be such that

x > r + e.

Finally, applying Property 2 (i) with k = 7, we deduce from 5 > u 7 that

r + e ∈ I 1 ∪ I 2 ∪ I 3 . ( 32
)
Analysis depending on whether p 1 + p 2 belongs to I 1 , I 2 , or I 3 . We now consider each of these three cases in turn in order to deduce the possible values for x together with the corresponding intervals for x.

If p 1 + p 2 ∈ I 1 then r = 1 + 2u and, using (32), we deduce that e ∈ [-u, -u + 6u 2 ) =:

I (1) 1 ∪ [u, u + 4u 2 ) =: I (2) 1 ∪ [3u, 3u + 2u 2 ) =: I (3) 1
.

These intervals are valid no matter what the radix β and the tie-breaking strategy of RN. If in addition β is odd or RN rounds 1 + u down to 1, as is the case when condition (C) holds, then we have further -u < e and -u < e 1 + e 2 (33) (see Appendix C for a detailed proof); thus, in this special case one can in particular replace I

1 \{-u} = (-u, -u + 6u 2 ).

Then, for each of the four intervals I we deduce the value of x and a range for e 1 +e 2 and for x, as shown in Table 1 below. The value of x follows immediately from rounding the sum 1 + 2u + e up to the nearest floating-point number. Let us now bound e 1 +e 2 from below when e ∈ I , we use the strict lower bound already mentioned in (33). The remaining lower bounds and upper bounds for e 1 + e 2 are all deduced from the fact that e 1 + e 2 = e(1 + δ) with |δ| u. Finally, since x = p 1 + p 2 + e 1 + e 2 , the range of x is obtained

1 + 2u [-u -2 β u 2 , -u + 7u 2 -6u 3 ) [1 -2 β u 2 , 1 + 13u 2 -6u 3 ) I (1,C) 1 1 + 2u (-u, -u + 7u 2 -6u 3 ) (1, 1 + 13u 2 -6u 3 ) I (2) 1 1 + 4u [u -u 2 , u + 5u 2 + 4u 3 ) [1 + 2u -u 2 , 1 + 2u + 11u 2 + 4u 3 ) I (3) 1 1 + 6u [3u -3u 2 , 3u + 5u 2 + 2u 3 ) [1 + 4u -3u 2 , 1 + 4u + 11u 2 + 2u 3 )
simply by adding the range Applying (32) with r = 1 + 4u and using the fact that -3u + 6u 2 -5 2 u(1 + u 1 ) for β p-1 10, we finally obtain e ∈ [-u, -u + 4u 2 ) =:

I 1 = [1 + u, 1 + u + 6u
I (1) 2 ∪ [u, u + 2u 2 ) =: I (2) 2 .
Then we proceed in the same way as in the previous case and deduce for the intervals I

1 + 4u [-u -u 2 , -u + 5u 2 -4u 3 ) [1 + 2u -u 2 , 1 + 2u + 9u 2 -4u 3 ) I (2) 2 1 + 6u [u -u 2 , u + 3u 2 + 2u 3 ) [1 + 4u -u 2 , 1 + 4u + 7u 2 + 2u 3 ) If p 1 + p 2 ∈ I 3 then r = 1 + 6u
and, using (32) and recalling from (31) that e cannot be smaller than -4u, we deduce that for β p-1 10

e ∈ [-3u, -3u + 4u 2 ) =:

I (1) 3 ∪ [-u, -u + 2u 2 ) =: I (2) 3
.

Proceeding as in the two previous cases, we then obtain the values of x and the ranges of e 1 + e 2 and x shown in Table 3.

1 + 4u [-3u -3u 2 , -3u + 7u 2 -4u 3 ) [1 + 2u -3u 2 , 1 + 2u + 9u 2 -4u 3 ) I (2) 3 1 + 6u [-u -u 2 , -u + 3u 2 -2u 3 ) [1 + 4u -u 2 , 1 + 4u + 5u 2 -2u 3 )
Conclusion. For β p-1 10, the second and fourth columns of Tables 1-3 lead to x x > 0 and to the following relative error bounds:

| x -x| |x| = x x -1                                            2βu+2u 2 β-2u 2 = 2u + 2 β u 2 + O(u 3 ) if e ∈ I (1) 1
, 2u -η for some η > 0 if e ∈ I (1,C) 1 , 2u+u 2 1+2u-u 2 = 2u -3u 2 + O(u 3 ) if e ∈ I (2)
1 ,

2u+3u 2 1+4u-3u 2 = 2u -5u 2 + O(u 3 ) if e ∈ I (3) 1 , 2u+u 2 1+2u-u 2 = 2u -3u 2 + O(u 3 ) if e ∈ I (1)
2 ,

2u+u 2 1+4u-u 2 = 2u -7u 2 + O(u 3 ) if e ∈ I (2) 2 , 2u+3u 2 1+2u-3u 2 = 2u -u 2 + O(u 3 ) if e ∈ I (1) 3
, 2u+u 2 1+4u-u 2 = 2u -7u 2 + O(u 3 ) if e ∈ I (2) 3 .
From these eight cases and for β p-1 10, it is easily deduced that when discarding the interval I

(1)
1 , the error is always less than 2u, while it is at most . This shows (29) and, therefore, concludes the analysis of the case where 4 > u 7 and 5 > u 7 .

B.2.4

Case where 4 < -u 7 and 5 < -u 7

In this case our goal is to show (17), that is, | x -x|/|x| < 2u -u 2 for β p-1 10. This bound shall be obtained in the same way as in Section B.2.3, but since it is less than 2u independently of the tie-breaking strategy of RN, the analysis will be slightly simpler.

Preliminaries. We can assume as before that 1 < p 1 + p 2 < β and, applying Property 2 (ii) with k = 7, we deduce from 4 < -u 7 that (35)

Analysis depending on whether p 1 + p 2 belongs to I 0 , I 1 , or I 2 . As in the previous section, we will now consider each of these three cases in turn in order to deduce values for x and intervals for x.

If p 1 + p 2 ∈ I 0 then r = 1 and, using (35) together with the fact that e is less than 4u, we deduce that e ∈ (u -6u 2 , u] =: I and in each case the value of x and the ranges of e 1 + e 2 and x are as shown in Table 5.

If p 1 + p 2 ∈ I 2 then r = 1 + 4u. and for each of these three intervals, the corresponding information about x, e 1 + e 2 , and x appears in Table 6.

Property 3 .

 3 Let i ∈ Z and a, b ∈ F be such that β i |ab| < β i+1 . Then ab -RN(ab) is an integer multiple of β i-2p+1 . Proof. Writing a = Aβ ea-p+1 and b = Bβ e b -p+1 with A, B two integers such that β p-1 |A|, |B| β p -1, we have ab = ABβ ea+e b -2p+2 . Hence ab, RN(ab), and ab -RN(ab) are integer multiples of β ea+e b -2p+2 . Now, i is either e a + e b or e a + e b + 1, so ab -RN(ab) is always an integer multiple of

see [ 6 , 1 = 1 + 2u and e 1 =

 6111 Theorem 3.2]. Hence, using the assumption that RN(1 + u) = 1 + 2u, we have p -u.

-2 β u 2 .

 2 in this case |e| u, we have ufp(e) ufp(u) = β -p ; using |e 1 +e 2 -e| u ufp(e) then gives e 1 + e 2 e -u ufp(e) -u When e ∈ I (1,C) 1

  collected in Table2: rounding 1+4u+e up to the nearest floating-point number gives the values of x, applying e 1 + e 2 = e(1 + δ) with |δ| u gives the ranges of e 1 + e 2 , and adding the ranges of e 1 + e 2 to the one of p 1 + p 2 leads to the ranges of x.

p 1 +=: I 1 ∪=: I 2 .

 112 p 2 ∈ (1 + u -6u 2 , 1 + u] =: I 0 ∪ (1 + 3u -4u 2 , 1 + 3u] (1 + 5u -2u 2 , 1 + 5u]Since 4 < 0, we have the strict inequality r < p 1 + p 2 and then, no matter what the tie-breaking strategy of rounding to nearest, r ∈ {1, 1 + 2u, 1 + 4u}. Since β p-1 10 and since the right endpoint of I 2 is not larger than the one of I 3 from Section B.2.3, the bound |e| < 4u established in (31) still holds. From β p-1 10 and 5 < -u 7 it then follows that 1 < r + e < β and x < r + e and, using again Property 2 (ii) with k = 7, that r + e ∈ I 0 ∪ I 1 ∪ I 2 .

.

  Proceeding as in Section B.2.3 we can set up the table below. If p 1 + p 2 ∈ I 1 then r = 1 + 2u. Consequently, (35) implies e ∈ (-u -6u 2 , -u]

Using ( 35 ) 2 ∪

 352 , we deduce e ∈ (-3u -6u 2 , -3u] (-u -4u 2 , -u]

  2p+1 and the integer E 1 + E 2 must satisfy |E 1 + E 2 | β p . This means that e 1 + e 2 ∈ F or, equivalently, 3 = 0. It then follows from (24) that | x -x|/|x| | 5 | u 1 , which proves (11).

B.1.2 Case β i-1 |cd| < β i ab < β i+1

Table 1 :

 1 Ranges or values of e, x, e 1 + e 2 , x in the case p 1 + p 2 ∈ I 1 .

		e	x	e 1 + e 2	x
	I	(1) 1

  2 ) of p 1 + p 2 to the range of e 1 + e 2 just computed.

	If p 1 + p 2 ∈ I 2 then
				r = 1 + 4u.
	Furthermore, since 0 < -p 2 < 1 2 p 1 , we also have the lower bound
	-	5 2	u e 1 + e 2 .	(34)
	(See Appendix C for a detailed proof.) Recalling that e = (e 1 + e 2 )(1 + 3 ) with
	| 3 | u 1 , we deduce that
	-	5 2	u(1 + u 1 ) e.

Table 2 :

 2 Ranges or values of e, x, e 1 + e 2 , x in the case p 1 + p 2 ∈ I 2 .

		e	x	e 1 + e 2	x
	I	(1) 2

Table 3 :

 3 Ranges or values of e, x, e 1 + e 2 , x in the case p 1 + p 2 ∈ I 3 .

		e	x	e 1 + e 2	x
	I	(1) 3

Table 4 :

 4 Ranges or values of e, x, e 1 + e 2 , x in the casep 1 + p 2 ∈ I 0 . -7u 2 + 6u 3 , u + u 2 ] (1 + 2u -13u 2 + 6u 3 , 1 + 2u + u 2 ] + 4u 3 , 3u + 3u 2 ] (1 + 4u -13u 2 + 4u 3 , 1 + 4u + 3u 2 ]

	e	x	e 1 + e 2	x
	I (u I (1) 0 1 (2) 0 1 + 2u (3u -7u 2

Table 5 :

 5 Ranges or values of e, x, e 1 + e 2 , x in the casep 1 + p 2 ∈ I 1 . -7u 2 -6u 3 , -u + u 2 ] (1 + 2u -11u 2 -6u 3 , 1 + 2u + u 2 ] -5u 2 + 4u 3 , u + u 2 ] (1 + 4u -9u 2 + 4u 3 , 1 + 4u + u 2 ] 5u 2 + 2u 3 , 3u + 3u 2 ] (1 + 6u -9u 2 + 2u 3 , 1 + 6u + 3u 2 ]

	e	x	e 1 + e 2	x
	I (-u I (1) 1 1 (2) 1 1 + 2u (u I (3) 1 1 + 4u (3u -

Table 6 :

 6 Ranges or values of e, x, e 1 + e 2 , x in the casep 1 + p 2 ∈ I 2 . -6u 3 , -3u + 3u 2 ] (1 + 2u -11u 2 -6u 3 , 1 + 2u + 3u 2 ] -5u 2 -4u 3 , -u + u 2 ] (1 + 4u -7u 2 -4u 3 , 1 + 4u + u 2 ] -3u 2 + 2u 3 , u + u 2 ] (1 + 6u -5u 2 + 2u 3 , 1 + 6u + u 2 ]

	e	x	e 1 + e 2	x
	I (-3u -9u 2 I (1) 2 1 (2) 2 1 + 2u (-u I (3) 2 1 + 4u (u

Acknowledgments

This work was supported in part by the French National Research Agency, under grant ANR-13-INSE-0007 (MetaLibm project). I am also grateful to the referees for carefully reading the manuscript and offering helpful suggestions.

Consequently, p 1 + p 2 = 1 + u, which rounds to r = 1 + 2u.

On the other hand, noticing that e 1 = p 2 we obtain e = RN(p 2 + e 2 ) = p 2 , that is, e = -u.

Hence r + e = 1 + u, which rounds to

Finally, we deduce from ( 18) and ( 19) that

which is such that 0 < x < x. Thus, overall, | x-x| |x| = x x -1 = 2u β+u+(2-2β)u 2 β-2u 2 +(4β-4)u 3 , and one can check that the latter quantity is larger than 2u + 2 β u 2 -4u 3 . This concludes the proof of Theorem 2.

APPENDIX A Analysis when the two products have the same sign

When the products ab and cd have the same sign, the assumption in [START_REF] Muller | On the error of computing ab + cd using Cornea, Harrison and Tang's method[END_REF] implies that they are both positive and that K in ( 3) is equal to one. Using (4), we deduce that

Using the obvious inequality

), which is slightly larger than 2u for any radix and tie-breaking rule. Instead, we consider two sub-cases separately, as follows.

A.1 Case where 4 u 3 or 5 u 3

In this case

1 , we have

which is the bound claimed in (12).

B.2 When

We begin by noting that in this case the ratio K = ab+|cd| ab-|cd| is at most about 3: since (1 + u 1 )ab p 1 and |p 2 | (1 -u 1 )|cd|, the assumption 1 2

and, since

Combining ( 27) and ( 4), we obtain

For the same reason as in Section A, applying to (28) the straightforward inequality

is not enough for our purpose, since the resulting relative error bound would then have the form 2u + 5u 2 + O(u 3 ). In order to achieve the sharper bounds claimed in Theorem 1, we shall refine further this analysis by examining separately four cases defined by the pair ( 4 , 5 ).

B.2.1 Case where 4 and 5 have opposite signs

In this case | 4 + 5 + 4 5 | u 1 and it follows from (28) that

It is then easily checked that this bound is at most u + 7u 2 when u 1/8 or, equivalently, when β p-1 4, which proves (13).

B.2.2 Case where

In this case

, and applying (28) then leads to

< 2u for β p-1 24, which proves (14).

Conclusion.

For β p-1 10, the second and fourth columns of Tables 4-6 imply x x and thus the following relative error bounds:

0 ,

1 ,

1 ,

1 ,

2 ,

2 ,

Finally, it is easily checked that all these bounds are less than 2u -u 2 , which proves (17) and finishes the case where 4 < -u 7 and 5 < -u 7 .

C Proofs of (33) and (34)

We begin with the following lemma, which will be used to prove (33) .

Lemma 1. If β is even then u ∈ F. If β is odd then u is a midpoint for F and its expansion in radix β has the form

Proof. By definition, u = β 2 • β -p . If β is even then β/2 is an integer less than β p , which implies u ∈ F. Assume now that β is odd. In this case

which gives the announced expansion in radix β. It remains to check that u is a midpoint for F. From the radix-β expansion of u, we deduce that the two consecutive elements

The associated midpoint for F is thus f1+f2

Proof of (33). This amounts to checking that -u < e and -u < e 1 + e 2 when (i) r = 1 + 2u;

(ii) -u e ∈ F;

(iii) x := RN(r + e) > r + e;

(iv) condition (C) holds, that is,

If β is odd, then Lemma 1 implies -u ∈ F and it follows from (ii) that -u < e. If RN(1 + u) = 1, then e = -u, for otherwise (i) and (iii) would yield RN(1 + u) > 1 + u, a contradiction. Hence, we have in both cases -u < e.

Let us now check that -u < e 1 + e 2 . If β is odd, then -u < e with e = RN(e 1 + e 2 ) and, by Lemma 1, -u is a midpoint for F. The definition of RN thus implies that -u e 1 + e 2 . Now, since e 1 , e 2 are in F and since, on the other hand, -u has infinitely many radix-β digits, we must have -u = e 1 + e 2 . If β is even then Lemma 1 implies -u ∈ F, so that by the monotonicity of RN the strict inequality -u < RN(e 1 + e 2 ) shown above leads to -u < e 1 + e 2 .

Proof of (34). This amounts to checking that - 5 2 u e 1 + e 2 when (i)

Applying (i) gives 1 + 3u -p 2 p 1 < 1 + 3u + 4u 2 -p 2 and then, using (ii) and u 1/4, we obtain 1 + 3u < p 1 < 2 + 6u + 8u