Mioara Joldes
email: joldes@laas.fr

Jean-Michel Muller
email: jean-michel.muller@ens-lyon.fr

J.-M Muller

V Popescu

5res Tight and Rigorous Error Bounds for Basic Building Blocks of Double-Word Arithmetic

Keywords: CCS Concepts:, Mathematics of computing → Mathematical software, • Theory of computation → Design and analysis of algorithms, Numeric approximation algorithms, Floating-point arithmetic, double-word arithmetic, double-double arithmetic, error-free transforms

We analyze several classical basic building blocks of double-word arithmetic (frequently called "double-double arithmetic" in the literature): the addition of a double-word number and a floating-point number, the addition of two double-word numbers, the multiplication of a double-word number by a floating-point number, the multiplication of two double-word numbers, the division of a double-word number by a floating-point number, and the division of two double-word numbers. For multiplication and division we get better relative error bounds than the ones previously published. For addition of two double-word numbers, we show that the previously published bound was incorrect, and we provide a new relative error bound. We introduce new algorithms for division. We also give examples that illustrate the tightness of our bounds.

INTRODUCTION AND NOTATION

Some calculations require a precision significantly higher than the one offered by the binary64 (also known as "double-precision") format. A typical example is the evaluation of transcendental functions in binary64 arithmetic with correct rounding: If all intermediate calculations are done in the target precision, then it is very difficult to guarantee last-bit accuracy in the final result. For instance, the CRLibm library of correctly rounded elementary functions uses "double-double" or "triple-double" operations in critical parts [START_REF] Florent De Dinechin | Towards the post-ultimate libm[END_REF]. Double-double arithmetic has also been used with success in Basic Linear Algebra Subroutines (BLAS) [START_REF] Li | Design, implementation and testing of extended and mixed precision BLAS[END_REF]. Other examples where higher-precision arithmetic has been useful, mentioned by [START_REF] Briggs | The doubledouble library[END_REF] or [START_REF] Bailey | High-precision computation: Mathematical physics and dynamics[END_REF], are studies of dynamical systems, the calculation of two-loop integrals for radiative corrections in muon decay, experimental mathematics, supernova simulations, and studies of the fine structure constant of physics.

There exist very good arbitrary precision libraries, such as GNU-MPFR [START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with correct rounding[END_REF]). However, if one only needs calculations accurate within around 120 bits in a few critical parts of a numerical program, using such libraries will involve a significant penalty in terms of speed and memory consumption.

Although the binary128 format (frequently called "quad-precision") was specified by the IEEE 754-2008 Standard on Floating-Point Arithmetic, it is seldom implemented in hardware. To our knowledge, the only commercially significant platform that has supported binary128 in hardware for the last decade has been the IBM z Systems [START_REF] Lichtenau | Quad precision floating point on the IBM z13[END_REF]). Thus, one will be tempted to use "double-double" arithmetic at times. Furthermore, even if hardwired binary128 arithmetic becomes commonplace, there will be a need for "double-quad" operations for carefully implementing very accurate binary128 elementary functions. Hence, designing and analyzing algorithms for double-word arithmetic is of interest. Double-word arithmetic, called "double-double" in most of the literature, consists in representing a real number as the unevaluated sum of two floating-point numbers. In all existing implementations, the underlying floating-point format is the binary64 format of the IEEE 754 Standard on Floating-Point Arithmetic (IEEE Computer Society 2008;[START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF], commonly called "double-precision" (hence the name "double-double").

Double-word arithmetic is not similar to a conventional, IEEE 754-like, floating-point arithmetic with twice the precision. It lacks many nice properties such as Lemma 1.2 below, clearly defined roundings, and so on. Furthermore, many algorithms have been published without a proof, or with error bounds that are sometimes loose, sometimes fuzzy (the error is "less than a small integer times u 2 "), and sometimes unsure. Kahan qualifies double-double arithmetic as an "attractive nuisance except for the BLAS" and even compares it to an unfenced backyard swimming pool! He also mentions [START_REF] Kahan | Lecture Notes on the Status of IEEE-754[END_REF]) that it "undermines the incentive to provide quadruple precision correctly rounded." The purpose of this article is to provide a rigorous error analysis of some double-word algorithms and to introduce a few new algorithms. We cannot suppress all the drawbacks mentioned by Kahan: Clearly, having in hardware a "real" floating-point arithmetic with twice the precision would be a better option. And yet, if rigorously proven and reasonably tight error bounds are provided, then expert programmers can rely on double-word arithmetic for extending the precision of calculations in critical places where the available floating-point arithmetic does not suffice.

Throughout this article, we assume a radix-2, precision-p floating-point (FP) arithmetic system, with unlimited exponent range and correct rounding. This means that our results will apply to "real-world" binary floating-point arithmetic, such as the one specified by the IEEE 754-2008Standard (IEEE Computer Society 2008;[START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF], provided that and overflow do not occur.

The notation RN(t ) stands for t rounded to the nearest FP number, ties-to-even. For instance, RN(c • d ) is the result of the FP multiplication c×d, assuming round-to-nearest rounding mode. The number ulp(x ), for x 0 is 2 log 2 |x | -p+1 , and u = 2 -p = 1 2 ulp(1) denotes the roundoff error unit. We will frequently use the three following, classical lemmas. Lemma 1.1. Let t ∈ R. If |t | ≤ 2 k , where k is an integer, then

|RN(t ) -t | ≤ u 2 • 2 k .
Lemma 1.2 (Sterbenz Lemma (Sterbenz 1974)). Let x and y be two positive FP numbers. If

x 2 ≤ y ≤ 2x, then x -y is a floating-point number, so RN(x -y) = x -y.
Lemma 1.3. If t ∈ R, then there exist ϵ 1 and ϵ 2 , both of absolute value less than or equal to u, such that

RN(t ) = t • (1 + ϵ 1 ) = t 1 + ϵ 2 .
The algorithms analyzed in this article use as basic blocks Algorithms 1, 2, and 3 below. They have been coined as "error-free transforms" by [START_REF] Rump | Transformations and ill-conditioned problems[END_REF].

Algorithms 1 and 2, introduced by [START_REF] Møller | Quasi double-precision in floating-point addition[END_REF], [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF], and [START_REF] Knuth | The Art of Computer Programming[END_REF] make it possible to compute both the result and the rounding error of a FP addition. We will choose between them depending on the information that we have on the input numbers.

ALGORITHM 1: -Fast2Sum(a, b). The Fast2Sum algorithm [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF].

s ← RN(a + b) z ← RN(s -a) t ← RN(b -z)
If a = 0 or b = 0, or if the floating-point exponents e a and e b satisfy e a ≥ e b , then s + t = a + b. Hence, t is the error of the FP addition s ← RN(a + b). In practice, condition "e a ≥ e b " may be hard to check. However, if |a| ≥ |b |, then that condition is satisfied.

ALGORITHM 2: -2Sum(a, b). The 2Sum algorithm [START_REF] Knuth | The Art of Computer Programming[END_REF][START_REF] Møller | Quasi double-precision in floating-point addition[END_REF].

s ← RN(a + b) a ← RN(s -b) b ← RN(s -a ) δ a ← RN(a -a ) δ b ← RN(b -b ) t ← RN(δ a + δ b )
Algorithm 2 gives the same results as Algorithm 1, but without any requirement on the exponents of a and b. It uses six FP operations for computing the result (instead of three for Algorithm 1), but on modern processors, comparing the absolute values of a and b and swapping them if needed before calling Algorithm 1 will in general be more time consuming than directly calling Algorithm 2. Hence, in general, Algorithm 1 is to be used only if we have preliminary information on the respective orders of magnitude of a and b. However, in all the algorithms presented below, a call to 2Sum can be replaced by a test and a call to Fast2Sum without changing the error bounds.

Let a and b be two FP numbers, with exponents e a and e b , respectively. Define π = RN(ab). The number ρ = abπ is a FP number. When the exponent range is not unbounded, this holds provided that e a + e b ≥ e min + p -1, where e min is the minimum exponent of the underlying FP format. See [START_REF] Nievergelt | Scalar fused multiply-add instructions produce floating-point matrix arithmetic provably accurate to the penultimate digit[END_REF] for a proof. The first algorithm introduced for computing π and ρ is in [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF] and [START_REF] Boldo | Pitfalls of a full floating-point proof: Example on the formal proof of the veltkamp/dekker algorithms[END_REF]. It requires 17 FP operations. When an FMA instruction is available, Algorithm Fast2Mult (Algorithm 3 below), mentioned by [START_REF] Kahan | Lecture Notes on the Status of IEEE-754[END_REF], only requires 2 FP operations for computing the same values.

ALGORITHM 3: -Fast2Mult(a, b). The Fast2Mult algorithm (see, for instance, [START_REF] Kahan | Lecture Notes on the Status of IEEE-754[END_REF], [START_REF] Nievergelt | Scalar fused multiply-add instructions produce floating-point matrix arithmetic provably accurate to the penultimate digit[END_REF][START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]). It requires the availability of a fused multiply-add (FMA) instruction for computing RN(abπ ).

π ← RN(a • b) ρ ← RN(a • b -π )
In the following, we will denote 2Prod an algorithm that computes π and ρ. It can be either Dekker's algorithm or Algorithm 3. However, when we count the number of floating-point operations required by the various algorithms presented in this article (in Table 1), we assume that Algorithm 3 is used. [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF] was the first to suggest using algorithms similar to Algorithm 1 and the equivalent (without FMA) of Agorithm 3 to manipulate numbers represented as unevaluated sums of two FP numbers. He called such numbers doublelength numbers. Dekker presented algorithms for adding, multiplying, and dividing double-word numbers. His addition and multiplication algorithms are very similar (in fact, mathematically equivalent) to Algorithms 5 and 10, analyzed below. His division algorithm was quite different (and less accurate) than the algorithms considered in this article. [START_REF] Linnainmaa | Software for doubled-precision floating-point computations[END_REF] suggested similar algorithms, assuming that an underlying extended precision format is available. We will not assume that hypothesis here.

Libraries that offer double-word arithmetic (with binary64 as the underlying floating-point format) have been written by Bailey [START_REF] Hida | C++/Fortran-90 double-double and quad-double package[END_REF] and [START_REF] Briggs | The doubledouble library[END_REF]. Briggs no longer maintains his library. Fairly recent functions for double-word arithmetic are included in the QD ("quad-double") library by Hida, Li, andBailey (Hida et al. 2012, 2001).

In Definition 1.4, we formally introduce the concept of double-word representation.

Definition 1.4. A double-word number x is the unevaluated sum x h + x of two floating-point numbers x h and x such that

x h = RN(x ).
The sequel of the article is organized as follows: Section 2 deals with the sum of a double-word number and a floating-point number; Section 3 is devoted to the sum of two double-word numbers; in Section 4, we consider the product of a double-word number by a floating-point number; in Section 5, we consider the product of two double-word numbers; Section 6 deals with the division of a double-word number by a floating-point number; and Section 7 is devoted to the division of two double-word numbers. All algorithms considered in this article return their results as a double-word number. We summarize our results in Table 1, in the Conclusion section.

ADDITION OF A DOUBLE-WORD NUMBER AND A FLOATING-POINT NUMBER

The algorithm implemented in the QD library [START_REF] Hida | C++/Fortran-90 double-double and quad-double package[END_REF] for adding a double-word number and a floating-point number is Algorithm 4 below.

Algorithm 4, or variants of it, implicitly appears in many "compensated summation" algorithms. Compensated summation algorithms aim at accurately computing the sum of several FP numbers. Most such algorithms implicitely represent, at intermediate steps of the summation, the sum of all input numbers accumulated so far as a double-word number. For instance the first two lines of Algorithm 4 constitute the internal loop of Rump, Ogita and Oishi's "cascaded summation" algorithm [START_REF] Ogita | Accurate sum and dot product[END_REF].

ALGORITHM 4: -DWPlusFP(x h , x , y). Algorithm for computing (x h , x ) + y in binary, precision-p, floating-point arithmetic, implemented in the QD library. The number x = (x h , x ) is a double-word number (i.e., it satisfies Definition 1.4).

1: (s h , s ) ← 2Sum(x h , y) 2: v ← RN(x + s ) 3: (z h , z ) ← Fast2Sum(s h , v) 4: return (z h , z )
To prove the correctness and bound the error of Algorithm 4 (and Algorithm 6 below), we will need the following lemma. That lemma is an immediate consequence of Property (2.16) in [START_REF] Rump | Accurate floating-point summation part I: Faithful rounding[END_REF]. (1)

Proof. First, the case x h + y = 0 is trivial, since s h = s = 0 and the computation is errorless. Now, without loss of generality, we can assume |x h | ≥ |y|. If this is not the case, since x h and y play a symmetrical role in the algorithm, then we can exchange them in our proof: We add the double word number (y, x ) and the floating-point number x h . 1 We also assume that x h is positive (otherwise we change the sign of all the operands) and that 1 ≤ x h ≤ 2 -2u (otherwise we scale the operands by a power of 2).

Define ϵ as the error committed at step 2, that is,

ϵ = v -(x + s ). 1. If -x h < y ≤ -x h /2, then Sterbenz Lemma implies s h = x h + y and s = 0. It follows that v = x . Lemma 2.1 implies |s h | ≥ 1 2 ulp(x h ), which implies |s h | ≥ |x |.
Hence Algorithm Fast2Sum introduces no error at line 3 of the algorithm, so z h + z = s h + v = x + y exactly.

If

-x h /2 < y ≤ x h , then 1 2 ≤ x h 2 < x h + y ≤ 2x h , so s h ≥ 1/2.
|ϵ | 2 -u ≤ 2u 2 2 -u < 2u 2 .
Notice that the bound (1) is very sharp. In fact, it is asymptotically optimal. This is shown by the following example:

x h = 1, x = (2 p -1) • 2 -2p , and y = -1 2 (1 -2 -p ), for which the computed sum is 1 2 + 3 • 2 -p-1 and the exact sum is 1 2 + 3 • 2 -p-1 -2 -2p
, resulting in a relative error,

2u 2 1 + 3u -2u 2 ≈ 2u 2 -6u 3 .
In the binary64 format (p = 53), this generic example gives an error,

1.99999999999999933 • • • × 2 -106 .

ADDITION OF TWO DOUBLE-WORD NUMBERS

Algorithm 5 below was first given by [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF], under the name of add2, with a slightly different presentation. Dekker did not use the 2Sum algorithm: Instead of Line 1, there was a comparison of |x h | and |y h | followed by a possible swap of x and y and a call to Fast2Sum. However, from a mathematical point of view, Dekker's algorithm and Algorithm 5 are equivalent: They always return the same result. This algorithm was then implemented by Bailey in the QD library [START_REF] Hida | C++/Fortran-90 double-double and quad-double package[END_REF] under the name of "sloppy addition."

ALGORITHM 5: -SloppyDWPlusDW(x h , x , y h , y ). "Sloppy" calculation of (x h , x ) + (y h , y ) in binary, precision-p, floating-point arithmetic. 1: (s h , s ) ← 2Sum(x h , y h ) 2: v ← RN(x + y ) 3: w ← RN(s + v) 4: (z h , z ) ← Fast2Sum(s h , w ) 5: return (z h , z )
Dekker proved an error bound on the order of (|x | + |y|) • 4u 2 . Notice the absolute values: When x and y do not have the same sign, there is no proof that the relative error is bounded. Indeed, the relative error can be so large that the obtained result has no significance at all. Consider, for instance, the case

x h = 1 + 2 -p+3 , x = -2 -p , y h = -1 -6 • 2 -p , and y = -2 -p + 2 -2p
. It leads to a computed value of the sum equal to zero, whereas the exact value is 2 -2p : The relative error is equal to 1. This is why the use of Algorithm 5 should be restricted to special cases such as, for instance, when we know in advance that the operands will have the same sign. When accurate computations are required, it is much more advisable to use the following algorithm, presented by [START_REF] Li | Design, Implementation and Testing of Extended and Mixed Precision BLAS[END_REF][START_REF] Li | Design, implementation and testing of extended and mixed precision BLAS[END_REF] and implemented in the QD library under the name of "IEEE addition." [START_REF] Li | Design, Implementation and Testing of Extended and Mixed Precision BLAS[END_REF][START_REF] Li | Design, implementation and testing of extended and mixed precision BLAS[END_REF] claim that in binary64 arithmetic (p = 53) the relative error of Algorithm 6 is upper bounded by 2 • 2 -106 . This bound is incorrect, as shown by the following example: If

ALGORITHM 6: -AccurateDWPlusDW(x h , x , y h , y ). Calculation of (x h , x ) + (y h , y ) in binary, precision-p, floating-point arithmetic. 1: (s h , s ) ← 2Sum(x h , y h ) 2: (t h , t ) ← 2Sum(x , y ) 3: c ← RN(s + t h ) 4: (v h , v ) ← Fast2Sum(s h , c) 5: w ← RN(t + v ) 6: (z h , z ) ← Fast2Sum(v h , w ) 7: return (z h , z )
x h = 9007199254740991, x = -9007199254740991/2 54 , y h = -9007199254740987/2, and y = -9007199254740991/2 56 , (2) then the relative error of Algorithm 6 is 2.24999999999999956 • • • × 2 -106 .
Note that this example is somehow "generic": In precision-p FP arithmetic, the choice

x h = 2 p -1, x = -(2 p -1) • 2 -p-1 , y h = -(2 p -5)/2, and y = -(2 p -1) • 2 -p-3
leads to a relative error that is asymptotically equivalent (as p goes to infinity) to 2.25u 2 . Now let us try to find a relative error bound. We are going to show the following result.

Theorem 3.1. If p ≥ 3, then the relative error of Algorithm 6 (AccurateDWPlusDW) is bounded by

3u 2 1 -4u = 3u 2 + 12u 3 + 48u 4 + • • • , (3)
which is less than 3u 2 + 13u 3 as soon as p ≥ 6.

Note that the conditions on p (p ≥ 3 for the bound (3) to hold, p ≥ 6 for the simplified bound 3u 2 + 13u 3 ) are satisfied in all practical cases.

Proof. First, we exclude the straightforward case in which one of the operands is zero. We can also quickly proceed with the case x h +y h = 0: The returned result is 2Sum(x , y ), which is equal to x + y, that is, the computation is errorless. Now, without loss of generality, we assume 1

≤ x h < 2, x ≥ |y| (which implies x h ≥ |y h |), and x h + y h nonzero. Notice that 1 ≤ x h < 2 implies 1 ≤ x h ≤ 2 -2u, since x h is a FP number.
Define ϵ 1 as the error committed at Line 3 of the algorithm:

ϵ 1 = c -(s + t h ) (4)
and ϵ 2 as the error committed at Line 5:

ϵ 2 = w -(t + v ).
(5)

1. If -x h < y h ≤ -x h /2.
Sterbenz Lemma, applied to the first line of the algorithm, implies

s h = x h + y h , s = 0, and c = RN(t h ) = t h . Define σ = 2 if y h ≤ -1, 1 if -1 < y h ≤ -x h /2. We have -x h < y h ≤ (1 -σ ) + x h 2 (σ -2), so 0 ≤ x h + y h ≤ 1 + σ • ( x h 2 -1) ≤ 1 -σu. Also, since x h is a multiple of 2u and y h is a multiple of σu, s h = x h + y h is a multiple of σu. Since s h is nonzero, we finally obtain σu ≤ s h ≤ 1 -σu. ( 6
)
We have |x | ≤ u and |y | ≤ σ 2 u, so

|t h | ≤ 1 + σ 2 u and |t | ≤ u 2 . ( 7
)
From Equation ( 6), we deduce that the floating-point exponent of s h is at least -p + σ -1. From Equation ( 7), the floating-point exponent of c = t h is at most -p + σ -1. Therefore, the Fast2Sum algorithm introduces no error at line 4 of the algorithm, which implies

v h + v = s h + c = s h + t h = x + y -t .
Equations ( 6) and ( 7) imply

|s h + t h | ≤ 1 + 1 - σ 2 u ≤ 1 + u 2 , so |v h | ≤ 1 and |v | ≤ u 2 .
From the bounds on |t | and |v |, we obtain:

|ϵ 2 | ≤ 1 2 ulp(t + v ) ≤ 1 2 ulp u 2 + u 2 = u 2 2 (8)
and

|ϵ 2 | ≤ 1 2 ulp 1 2 ulp(x + y ) + 1 2 ulp (x + y) + 1 2 ulp(x + y ) . (9)
Lemma 2.1 and |s h | ≥ σu imply that either

s h + t h = 0, or |v h | = |RN(s h + c)| = |RN(s h + t h )| ≥ σu 2 . If s h + t h = 0, then v h = v =
0 and the sequel of the proof is straightforward. Therefore, in the following, we assume |v h | ≥ σu 2 . Now,

• If |v h | = σu 2 , then |v + t | ≤ u |v h | + u 2 = σu 3 + u 2 , which implies |w | = |RN(t + v )| ≤ σu 2 = |v h |; • If |v h | > σu 2 , then, since v h is a FP number, |v h | is larger than or equal to the FP number immediately above σu 2 , which is σ (1 + 2u)u 2 . Hence |v h | ≥ σu 2 /(1 -u), so |v h | ≥ u • |v h | + σu 2 ≥ |v | + |t |. So, |w | = |RN(t + v )| ≤ |v h |.
Therefore, in all cases, Fast2Sum introduces no error at line 6 of the algorithm, and we have

z h + z = v h + w = x + y + ϵ 2 . ( 10
)
Directly using Equation ( 10) and the bound u 2 /2 on |ϵ 2 | to get a relative error bound would result in a large bound, because x + y may be small. However, when x + y is very small, some simplification occurs thanks to Sterbenz Lemma. First,

x h + y h is a nonzero multiple of σu. Hence, since |x + y | ≤ (1 + σ 2 )u, we have |x + y | ≤ 3 2 (x h + y h ).
Let us now consider the two possible cases:

• If -3 2 (x h + y h ) ≤ x + y ≤ -1 2 (x h + y h ), which implies -3 2 s h ≤ t h ≤ -1 2 s h , then
Sterbenz lemma applies to the floating-point addition of s h and c = t h . Therefore line 4 of the algorithm results in v h = s h and v = 0. An immediate consequence is ϵ 2 = 0, so z h + z = v h + w = x + y: the computation of x + y is errorless;

• If -1 2 (x h + y h ) < x + y ≤ 3 2 (x h + y h ), then 5 2 (x + y ) ≤ 3 2 (x h + y h + x + y ) = 3 2 (x + y), and -1 2 (x + y) < 1 2 (x + y ). Hence, |x + y | < |x + y|, so ulp(x + y ) ≤ ulp(x + y).
Combined with Equation ( 9), this gives

|ϵ 2 | ≤ 1 2 ulp 3 2 ulp(x + y) ≤ 2 -p ulp(x + y) ≤ 2 • 2 -2p • (x + y). 2. If -x h /2 < y h ≤ x h Notice that we have x h /2 < x h + y h ≤ 2x h , so x h /2 ≤ s h ≤ 2x h . Also notice that we have |x | ≤ u. • If 1 2 < x h + y h ≤ 2 -4u. Define σ = 1 if x h + y h ≤ 1 -2u, 2 if 1 -2u < x h + y h ≤ 2 -4u. We have σ 2 (1 -2u) ≤ s h ≤ σ (1 -2u) and |s | ≤ σ 2 u. ( 11
) When σ = 1, we necessarily have -x h /2 < y h < 0, so |y | ≤ u/2. And when σ = 2, |y h | ≤ x h ≤ 2 -2u implies |y | ≤ u. Hence we always have |y | ≤ σ 2 u. This implies |x + y | ≤ (1 + σ /2)u, therefore |t h | ≤ 1 + σ 2 u and |t | ≤ u 2 . (12) Now, |s + t h | ≤ (1 + σ )u, so |c | ≤ (1 + σ )u and |ϵ 1 | ≤ σu 2 . ( 13
)
Since s h ≥ 1/2 and |c | ≤ 3u, if p ≥ 3, then Algorithm Fast2Sum introduces no error at line 4 of the algorithm, that is,

v h + v = s h + c. Therefore |v h + v | = |s h + c | ≤ σ (1 -2u) + (1 + σ )u ≤ σ . This implies |v h | ≤ σ and |v | ≤ σ 2 u. ( 14
) Thus |t + v | ≤ u 2 + σ 2 u, so |w | ≤ σ 2 u + u 2 and |ϵ 2 | ≤ σ 2 u 2 . ( 15
)
From Equations ( 11) and ( 13), we deduce

s h + c ≥ σ 2 -u (2σ + 1), so |v h | ≥ σ 2 -u (2σ + 1). If p ≥ 3, then |v h | ≥ |w |, so Algorithm Fast2Sum introduces no error at line 6 of the algorithm, that is, z h + z = v h + w. Therefore, z h + z = x + y + η, with |η| = |ϵ 1 + ϵ 2 | ≤ 3σ 2 u 2 . Since x + y ≥ (x h -u) + (y h -u/2) > 1 2 -3 2 u if σ = 1, 1 -4u if σ = 2, the relative error |η|/(x + y) is upper bounded by 3u 2 1 -4u . • If 2 -4u < x h + y h ≤ 2x h , then 2 -4u ≤ s h ≤ RN(2x h ) = 2x h ≤ 4 -4u
v h + v = s h + c ≤ 4 -4u + 4u = 4, so v h ≤ 4 and |v | ≤ 2u. Thus, |t + v | ≤ 2u + u 2 . Hence, either |t + v | < 2u and |ϵ 2 | ≤ 1 2 ulp(t + v ) ≤ u 2 , or 2u ≤ t + v ≤ 2u + u 2 , in which case w = RN(t + v ) = 2u and |ϵ 2 | ≤ u 2 . In all cases, |ϵ 2 | ≤ u 2 . Also, s h ≥ 2 -4u and |c | ≤ 4u imply v h ≥ 2 -8u, and |t + v | ≤ 2u + u 2 implies |w | ≤ 2u.
Hence if p ≥ 3, then Algorithm Fast2Sum introduces no error at line 6 of the algorithm. All this gives

z h + z = v h + w = x + y + η, with |η| = |ϵ 1 + ϵ 2 | ≤ 3u 2 . Since x + y ≥ (x h -u) + (y h -u) > 2 -6u, the relative error |η|/(x + y) is upper bounded by 3u 2 2 -6u ,
The largest bound obtained in the various cases we have analyzed is 3u 2 1 -4u .

Elementary calculus shows that for u ∈ [0, 1/64] (i.e., p ≥ 6) this is always less than 3u 2 + 13u 3 .

The bound (3) is probably not optimal. The largest relative error we have obtain through many tests is around 2.25 × 2 -2p = 2.25u 2 . An example is the input values given in Equation (2), for which, with p = 53 (binary64 arithmetic), we obtain a relative error equal to 2.24999999999999956

• • • × 2 -106 .

MULTIPLICATION OF A DOUBLE-WORD NUMBER BY A

FLOATING-POINT NUMBER We first consider the following algorithm, suggested by [START_REF] Li | Design, Implementation and Testing of Extended and Mixed Precision BLAS[END_REF]: [START_REF] Li | Design, Implementation and Testing of Extended and Mixed Precision BLAS[END_REF][START_REF] Li | Design, implementation and testing of extended and mixed precision BLAS[END_REF] (with more detail in the technical report [START_REF] Li | Design, Implementation and Testing of Extended and Mixed Precision BLAS[END_REF], which is a preliminary version of the journal article [START_REF] Li | Design, implementation and testing of extended and mixed precision BLAS[END_REF]), give a relative error bound 4 • 2 -106 for Algorithm 7 when the underlying floating-point arithmetic is binary64 (i.e., p = 53). Below, we prove an improved sharp relative error bound, even in the more general context of precision-p arithmetic. More precisely, we have the following. Theorem 4.1. If p ≥ 4, then the relative error of Algorithm 7 (DWTimesFP1) is bounded by

ALGORITHM 7: -DWTimesFP1(x h , x , y). Calculation of (x h , x ) × y in binary, precision-p, floating- point arithmetic. 1: (c h , c 1 ) ← 2Prod(x h , y) 2: c 2 ← RN(x • y) 3: (t h , t 1 ) ← Fast2Sum(c h , c 2 ) 4: t 2 ← RN(t 1 + c 1 ) 5: (z h , z ) ← Fast2Sum(t h , t 2 ) 6: return (z h , z )
3 2 u 2 + 4u 3 .
Proof. One easily notices that if x = 0, or y = 0, or y is a power of 2, the obtained result is exact. Therefore, without loss of generality, we can assume 1 ≤ x h ≤ 2 -2u and 1 + 2u ≤ y ≤ 2 -2u. This gives 1

+ 2u ≤ x h y ≤ 4 -8u + 4u 2 , so 1 + 2u ≤ c h ≤ 4 -8u (16)
and

|c 1 | ≤ 1 2 ulp(4 -8u) = 2u. ( 17
)
From |x | ≤ u and y ≤ 2 -2u, we deduce

|c 2 | ≤ 2u -2u 2 , ( 18
) so ϵ 1 = x y -c 2 satisfies |ϵ 1 | ≤ u 2 .
From Equations ( 16) and ( 18), we deduce that Algorithm Fast2Sum introduces no error at line 3 of the algorithm, that is,

t h + t 1 = c h + c 2 . Also, we deduce that 1 = RN(1 + 2u 2 ) ≤ t h ≤ RN(4 -6u -2u 2 ) = 4 -8u (19) and |t 1 | ≤ 1 2 ulp(4 -8u) = 2u. ( 20
)
From Equations ( 17) and ( 20), we obtain

|t 2 | ≤ RN(4u) = 4u, (21)
and we find that

ϵ 2 = t 2 -(t 1 + c 1 ) satisfies |ϵ 2 | ≤ 2u 2 . Define ϵ = ϵ 2 -ϵ 1 .
Using Equations ( 19) and ( 21), we deduce that Algorithm Fast2Sum introduces no error at line 5 of the algorithm. Therefore,

z h + z = t h + t 2 = t h + t 1 + c 1 + ϵ 2 = c h + c 2 + c 1 + ϵ 2 = x h y + x y -ϵ 1 + ϵ 2 = xy + ϵ. ( 22
)
Hence the absolute error of Algorithm 7 is

|ϵ | ≤ |ϵ 1 | + |ϵ 2 | ≤ 3u 2 .
Let us now consider two possible cases:

1. If x h y ≥ 2, then xy ≥ x h (1 -u)y ≥ 2 -2u. This leads to a relative error |ϵ/(xy)| bounded by 3u 2 2 -2u = 3 2 u 2 + 3 2 u 3 + 3 2 u 4 + • • • . ( 23
2 2 -4u -3u 2 = 3 2 u 2 + 3u 3 + 33 4 u 4 + • • • . ( 24
) If t h < 2, then |t 1 | ≤ u, |t 2 | ≤ 2u, and |ϵ 2 | ≤ u 2 .
Hence, a first upper bound on |ϵ | is 2u 2 . However, some refinement is possible.

• First, if |c 2 | < u, then |ϵ 1 | ≤ u 2 /2, which implies |ϵ | ≤ 3u 2 /2. • Second, if |c 2 | ≥ u, then c 2 is a multiple of ulp(u) = 2u 2
, so t 1 is a multiple of 2u 2 . Also, since x h and y are multiple of 2u, x h y is a multiple of 4u 2 , so c 1 is a multiple of 4u 2 . Hence, t 1 + c 1 is a multiple of 2u 2 of absolute value less than or equal to 2u. This implies that t 1 + c 1 is a FP number, hence RN(t 1 + c 1 ) = t 1 + c 1 and ϵ 2 = 0.

Therefore, when t h < 2, |ϵ | is upper bounded by 3u 2 /2 so the relative error |ϵ/(xy)| is bounded by

3 2 u 2 (1 -u)(1 + 2u) ≤ 3 2 u 2 . ( 25
)
The largest of the three bounds (Equations ( 23), ( 24), and ( 25)) is the second one. It is less than 3 2 u 2 + 4u 3 as soon as u ≤ 1/16. This proves the theorem. The bound given by Theorem 4.1 is very sharp. For instance, in binary32 arithmetic (p = 24), with x h = 8388609, x = 4095/8192, and y = 8389633, the relative error of Algorithm 7 is 1.4993282

• • • × 2 -48 .
In Bailey's QD library [START_REF] Hida | C++/Fortran-90 double-double and quad-double package[END_REF] as well as in Briggs' library [START_REF] Briggs | The doubledouble library[END_REF], another algorithm (Algorithm 8 below) is suggested for multiplying a double-word number by a floatingpoint number.

ALGORITHM 8: -DWTimesFP2(x h , x , y). Algorithm for computing (x h , x ) × y in binary, precision-p, floating-point arithmetic, implemented in the QD library.

1: (c h , c 1 ) ← 2Prod(x h , y) 2: c 2 ← RN(x • y) 3: c 3 ← RN(c 1 + c 2 ) 4: (z h , z ) ← Fast2Sum(c h , c 3 ) 5: return (z h , z )
Algorithm 8 is faster than Algorithm 7 (we save one call to Fast2Sum), but it is less accurate: There are input values for which the error attained using Algorithm 8 is larger than the bound given by Theorem 4.1. An example with p = 53 is x h = 4525788557405064, x = 8595672275350437/2 54 , and y = 5085664955107621, for which the relative error is 2.517

• • • × 2 -106 .
Hence, the relative error bound we are going to prove for Algorithm 8 is necessarily larger than the one we had for Algorithm 7. More precisely, we have the following. Theorem 4.2. If p ≥ 3, then the relative error of Algorithm 8 (DWTimesFP2) is less than or equal to 3u 2 .

Proof. The proof is very similar to (in fact, simpler than) the proof of Theorem 4.1. Without loss of generality, we can assume 1 ≤ x h ≤ 2 -2u and 1 ≤ y ≤ 2 -2u. Since the analysis of the case y = 1 is straightforward, we even assume 1 + 2u ≤ y ≤ 2 -2u. This gives 1

+ 2u ≤ x h y ≤ 4 -8u + 4u 2 , so 1 + 2u ≤ c h ≤ 4 -8u and |c 1 | ≤ 2u. From |x | ≤ u and y ≤ 2 -2u we deduce |c 2 | ≤ 2u - 2u 2 , so ϵ 1 = x y -c 2 satisfies |ϵ 1 | ≤ u 2 . Now, |c 1 + c 2 | ≤ 4u -2u 2 , hence |c 3 | ≤ 4u, and c 3 = c 1 + c 2 + ϵ 2 , with |ϵ 2 | ≤ 2u 2 .
From |c 3 | ≤ 4u and c h ≥ 1 + 2u we deduce that Algorithm Fast2Sum introduces no error at line 4 of the algorithm.

Hence,

z h + z = c h + c 3 = xy -ϵ 1 + ϵ 2 , and | -ϵ 1 + ϵ 2 | ≤ 3u 2 . Since xy ≥ (x h -u)y ≥ (1 -u)(1 + 2u
) ≥ 1, we deduce that the relative error of Algorithm 8 is less than 3u 2 .

If an FMA instruction is available, then we can improve Algorithm 8 by merging lines 2 and 3 of the algorithm and obtain Algorithm 9: ALGORITHM 9: -DWTimesFP3(x h , x , y). Algorithm for computing (x h , x ) × y in binary, precision-p, floating-point arithmetic, assuming an FMA instruction is available.

1: (c h , c 1 ) ← 2Prod(x h , y) 2: c 3 ← RN(c 1 + x y) 3: (z h , z ) ← Fast2Sum(c h , c 3 ) 4: return (z h , z )
This results in a better error bound as follows.

Theorem 4.3. If p ≥ 3, then the relative error of Algorithm 9 (DWTimesFP3) is less than or equal to 2u 2 .

The proof is very similar to the proof of Theorem 4.2, so we omit it. The bound provided by Theorem 4.3 is sharp. For instance, in binary64 arithmetic (p = 53), we attain error 1.984 • • • × 2 -106 for x h = 4505619370757448, x = -9003265529542491/2 54 , and y = 4511413997183120.

MULTIPLICATION OF TWO DOUBLE-WORD NUMBERS

Algorithm 10 below was first suggested by Dekker (under the name mul2 in [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF]). It has been implemented in the QD library [START_REF] Hida | C++/Fortran-90 double-double and quad-double package[END_REF]) and in Briggs' library [START_REF] Briggs | The doubledouble library[END_REF] for multiplying two double-word numbers.

Dekker proved a relative error bound 11u 2 . We are going to show the following:

ALGORITHM 10: -DWTimesDW1(x h , x , y h , y ). Algorithm for computing (x h , x ) × (y h , y ) in binary, precision-p, floating-point arithmetic, implemented in the QD library.

1: (c h , c 1 ) ← 2Prod(x h , y h ) 2: t 1 ← RN(x h • y ) 3: t 2 ← RN(x • y h ) 4: c 2 ← RN(t 1 + t 2 ) 5: c 3 ← RN(c 1 + c 2 ) 6: (z h , z ) ← Fast2Sum(c h , c 3 ) 7: return (z h , z )
Theorem 5.1. If p ≥ 4, then the relative error of Algorithm 10 (DWTimesDW1) is less than or equal to 7u 2 /(1 + u) 2 < 7u 2 .

In the proof of Theorem 5.1, we will use the following lemma: Lemma 5.2. Let a and b be two positive real numbers.

If ab ≤ 2, a ≥ 1 and b ≥ 1, then a + b ≤ 2 √ 2.
The proof of the lemma is straightforward calculus. Let us focus on the proof of Theorem 5.1.

Proof. Without loss of generality, we assume that 1 ≤ x h ≤ 2 -2u and 1 ≤ y h ≤ 2 -2u. We have x h y h < 4, and

c h + c 1 = x h y h ,
with |c 1 | ≤ 2u. We also have

t 1 = x h y + ϵ 1 , with |x h y | ≤ 2u -2u 2 ,
z h + z = c h + c 3 = (x h y h -c 1 ) + c 1 + c 2 + ϵ 4 = x h y h + t 1 + t 2 + ϵ 3 + ϵ 4 = x h y h + x h y + x y h + ϵ 1 + ϵ 2 + ϵ 3 + ϵ 4 = xy -x y + ϵ 1 + ϵ 2 + ϵ 3 + ϵ 4 = xy + η, ( 26
)
with |η| ≤ u 2 + |ϵ 1 + ϵ 2 + ϵ 3 + ϵ 4 | ≤ 9u 2 .
Let us consider the following cases.

• If x h y h > 2, then, since x ≥ x h (1 -u) and y ≥ y h (1 -u), the relative error is bounded by 9u 2 2(1 -u) 2 . ( 27
)
• If x h y h ≤ 2, then |c 1 | ≤ u. Furthermore, Lemma 5.2 implies

x h + y h ≤ 2 √ 2. ( 28
)
We have

|t 1 | = |RN(x h y )| ≤ RN(x h u) = x h u,
and, similarly, |t 2 | ≤ y h u, so, using Equation ( 28)

|t 1 + t 2 | ≤ x h u + y h u ≤ 2 √ 2u.
Therefore, c 2 now satisfies

|c 2 | ≤ |t 1 + t 2 | + |ϵ 3 | ≤ 2 √ 2u + 2u 2 .
We now deduce

|c 1 + c 2 | ≤ u • (2 √ 2 + 1) + 2u 2 ≤ 4u
(as soon as u ≤ 1/16, that is, p ≥ 4). Therefore, |ϵ 4 | ≤ 2u 2 . In Equation [START_REF] Linnainmaa | Software for doubled-precision floating-point computations[END_REF], this results in |η| ≤ 7u 2 instead of 9u 2 . Notice that if x h = 1 or y h = 1, then either ϵ 1 = 0 or ϵ 2 = 0, which results in a significantly smaller bound for |η|. So we can assume that x h ≥ 1 + 2u (hence, x > 1 + u) and y h ≥ 1 + 2u (hence, y > 1 + u). Therefore the relative error is bounded by

7u 2 (1 + u) 2 < 7u 2 . ( 29
)
If p ≥ 4, then the bound of Equation ( 27) is less than the bound of Equation ( 29). This proves the theorem.

The bound 7u 2 provided by Theorem 5.1 is probably too pessimistic. The largest relative error we have encountered in our tests was 4.9916 × 2 -106 , obtained for p = 53, x h = 4508231565242345, x = -9007199254524053/2 54 , y h = 4504969740576150, and y = -4503599627273753/2 53 . In binary32 arithmetic (p = 24), the largest error obtained in our tests was 4.947 × 2 -48 for x h = 8399376, x = 16763823/2 25 , y h = 8414932, and y = 16756961/2 25 . Now, if a fused multiply-add instruction (FMA) is available, then we can slightly improve Algorithm 10, both in terms of speed and accuracy, by merging lines 3 and 4. Consider Algorithm 11.

ALGORITHM 11: -DWTimesDW2(x h , x , y h , y ). Algorithm for computing (x h , x ) × (y h , y ) in binary, precision-p, floating-point arithmetic, assuming an FMA instruction is available.

1: (c h , c 1 ) ← 2Prod(x h , y h ) 2: t ← RN(x h • y ) 3: c 2 ← RN(t + x y h ) 4: c 3 ← RN(c 1 + c 2 ) 5: (z h , z ) ← Fast2Sum(c h , c 3 ) 6: return (z h , z )
We have the following.

Theorem 5.3. If p ≥ 5, then the relative error of Algorithm 11 ( DWTimesDW2) is less than or equal to

6u 2 + 1 2 u 3 (1 + u) 2 < 6u 2 .
The proof is very similar to (in fact, simpler than) the proof of Theorem 5.1, and follows the same structure, so we omit it.

We do not know if the bound given by Theorem 5.3 is sharp. The largest relative error we have encountered during our intensive tests was, for binary64 (p = 53), 4.9433 × 2 -106 , obtained for x h = 4515802244422058, x = -2189678420952711/2 52 , y h = 4503988428047019, and y = -2248477851812015/2 52 . In binary32 arithmetic (p = 24), the largest error obtained in our tests was 4.936 × 2 -48 , for x h = 8404039, x = -8284843/2 24 , y h = 8409182, and y = -4193899/2 23 .

The accuracy of the multiplication of two double-word numbers can be improved even more by also computing the partial product x y . This gives Algorithm 12 below.

ALGORITHM 12: -DWTimesDW3(x h , x , y h , y ). Algorithm for computing (x h , x ) × (y h , y ) in binary, precision-p, floating-point arithmetic, assuming an FMA instruction is available.

1: (c h , c 1 ) ← 2Prod(x h , y h ) 2: t 0 ← RN(x • y ) 3: t 1 ← RN(x h • y + t 0 ) 4: c 2 ← RN(t 1 + x • y h ) 5: c 3 ← RN(c 1 + c 2 ) 6: (z h , z ) ← Fast2Sum(c h , c 3 ) 7: return (z h , z )
We have the following.

Theorem 5.4. If p ≥ 4, then the relative error of Algorithm 12 (DWTimesDW3) is less than or equal to

5u 2 + 1 2 u 3 (1 + u) 2 < 5u 2 .
The proof is very similar to the proof of Theorem 5.1, and follows the same structure, so we omit it. We do not know if the bound given by Theorem 5.4 is sharp. The largest relative error we have encountered n intensive tests was (for p = 53) 3.936 × 2 -106 , obtained for x h = 4510026974538724, x = 4232862152422029/2 53 , y h = 4511576932111935, and y = 2250098448199619/2 52 .

DIVISION OF A DOUBLE-WORD NUMBER BY A FLOATING-POINT NUMBER

Before presenting algorithms for dividing a double-word number by a floating-point number, let us recall a classical result (see, for instance, [START_REF] Boldo | Representable correcting terms for possibly underflowing floating point operations[END_REF]), easy to prove, and common knowledge among the designers of Newton-Raphson-based division algorithms. Property 6.1 (Theorem 4 in [START_REF] Boldo | Representable correcting terms for possibly underflowing floating point operations[END_REF] with the additional assumption that we have an unbounded exponent range). If x and y are FP numbers with y 0, and if t = RN(x/y), then ytx is a FP number.

The algorithm suggested by [START_REF] Li | Design, Implementation and Testing of Extended and Mixed Precision BLAS[END_REF] for dividing a double-word number by a floatingpoint number is similar to Algorithm 13 below.

ALGORITHM 13: -DWDivFP1(x h , x , y). Calculation of (x h , x ) ÷ y in binary, precision-p, floating-point arithmetic.

1: t h ← RN(x h /y) 2: (π h , π ) ← 2Prod(t h , y) 3: (δ h , δ ) ← 2Sum(x h , -π h ) 4: δ ← RN(x -π ) 5: δ ← RN(δ + δ ) 6: δ ← RN(δ h + δ ) 7: t ← RN(δ /y) 8: (z h , z ) ← Fast2Sum(t h , t ) 9: return (z h , z )
Algorithm 13 can be simplified. We have

t h = (x h /y)(1 + ϵ 0 ) and π h = t h y(1 + ϵ 1 ), with |ϵ 0 |, |ϵ 1 | ≤ u. Hence, (1 -u) 2 x h ≤ π h ≤ (1 + u) 2 x h .
Therefore, as soon as p ≥ 2 (i.e., u ≤ 1/4), π h is within a factor 2 from x h . Sterbenz Lemma (Lemma 1.2) therefore implies that x hπ h is a floating-point number. As a consequence, we always have δ h = x hπ h , δ = 0, line 3 of the algorithm can be replaced by a simple subtraction, and we always have δ = δ = RN(xπ ). Therefore, the significantly simpler Algorithm 14 always returns the same result as Algorithm 13.

The authors of [START_REF] Li | Design, Implementation and Testing of Extended and Mixed Precision BLAS[END_REF] claim that their binary64 (i.e., p = 53) implementation of Algorithm 13 has a relative error bounded by 4 • 2 -106 . It is possible to show a slightly better bound, namely (7/2) • 2 -p . The proof can be found in the Ph.D. dissertation of one of the authors of this article [START_REF] Popescu | Towards Fast and Certified Multiple-precision Libraries[END_REF]). We will not detail it here, since we will suggest a slightly more accurate algorithm,2 Algorithm 15, obtained by modifying lines 4 and 5 of Algorithm 14.

ALGORITHM 14: -DWDivFP2(x h , x , y). Calculation of (x h , x ) ÷ y in binary, precision-p, floating-point arithmetic.

1: t h ← RN(x h /y) 2: (π h , π ) ← 2Prod(t h , y) 3: δ h ← RN(x h -π h ) = x h -π h (exact operation) 4: δ ← RN(x -π ) 5: δ ← RN(δ h + δ ) 6: t ← RN(δ /y) 7: (z h , z ) ← Fast2Sum(t h , t ) 8: return (z h , z )
ALGORITHM 15: -DWDivFP3(x h , x , y). Calculation of (x h , x ) ÷ y in binary, precision-p, floating-point arithmetic.

1: t h ← RN(x h /y) 2: (π h , π ) ← 2Prod(t h , y) 3: δ h ← RN(x h -π h ) = x h -π h (exact operation) 4: δ t ← RN(δ h -π ) = δ h -π (exact operation) 5: δ ← RN(δ t + x ) 6: t ← RN(δ /y) 7: (z h , z ) ← Fast2Sum(t h , t ) 8: return (z h , z )
Theorem 6.2. If p ≥ 4, then the relative error of Algorithm 15 (DWDivFP3) is bounded by 3u 2 .

The bound is sharp: In practice, the largest relative errors we have found in calculations were slightly less than 3u 2 . For instance, for p = 53, relative error 2.95157083 • • • × 2 -106 is attained for x h = 4588860379563012, x = -4474949195791253/2 53 , and y = 4578284000230917.

Before proving Theorem 6.2, let us prove the following Lemma.

Lemma 6.3. Assume a radix-2, precision-p, FP arithmetic. Let a and b be FP numbers between 1 and 2. Let u = 2 -p . The distance between RN(a/b) and a/b is less than

u -2u 2 /b if a/b ≥ 1; u/2 -u 2 /b otherwise.
Proof. It suffices to estimate the smallest possible distance between a/b and a "midpoint" (i.e., a number exactly halfway between two consecutive FP numbers

). Let a = M a • 2 -p+1 , b = M b • 2 -p+1 , with 2 p-1 ≤ M a , M b ≤ 2 p -1.
• If a/b ≥ 1, then a midpoint μ between 1 and 2 has the form

(2M μ + 1)/2 p , with 2 p-1 ≤ M μ ≤ 2 p -1. We have a b -μ = 2 p M a -M b (2Mμ + 1) 2 p M b .
The numerator, 2 p M a -M b (2Mμ +1), of that fraction cannot be zero: since 2Mμ +1 is odd, having 2 p M a = M b (2Mμ + 1) would require M b to be a multiple of 2 p , which is impossible since M b ≤ 2 p -1. Hence that numerator has absolute value at least 1. Hence

a b -μ ≥ 1 2 p M b = 2u 2 b .
• If a/b < 1, then the proof is similar. The only change is that a midpoint is of the form

(2M μ + 1)/2 p+1 .
In a recent article (see Jeannerod and Rump (2016, Table 1)), a similar bound is given for floatingpoint division. It could be used instead of Lemma 6.3, but we included the lemma for completeness. Let us now prove Theorem 6.2.

Proof. The case where y is a power of 2 is straightforward, so we omit it. Without loss of generality, we assume 1

≤ x h ≤ 2 -2u, so |x | ≤ u; and 1 + 2u ≤ y ≤ 2 -2u. Therefore, we have 1 2 -2u ≤ x h y ≤ 2 -2u 1 + 2u . ( 30
)
The quotient 1/(2 -2u) is always larger than 1/2 + u/2, and, as soon as p ≥ 4,

(2 -2u)/(1 + 2u) is less than 2 -5u. Therefore, 1 2 + u ≤ t h = RN x h y ≤ 2 -6u. ( 31
)
We have already proved, when discussing Algorithm 13, that

δ h = x h -π h . An immediate con- sequence is δ h -π = (x h -π h ) -(t h y -π h ) = x h -t h y. Property 6.1 implies that x h -t h y is an FP number. Hence δ t = RN(δ h -π ) = δ h -π = x h -t h y. We immediately deduce δ = RN(δ h -π + x ) = RN(x h + x -π h -π ) = RN(x -t h y).
Define ϵ 1 and ϵ 2 as the errors committed at lines 5 and 6 of the algorithm, more precisely,

ϵ 1 = δ -(x -t h y), ϵ 2 = t -δ/y. Lemma 6.3 implies |t h -x h /y| ≤ u -2u 2 /y, hence |t h y -x h | ≤ uy -2u 2 and hence |t h y -x | ≤ u (y + 1) -2u 2 . An immediate consequence of this is |t h y -x | < 3u, so |ϵ 1 | ≤ 2u 2 . Also δ y ≤ x y -t h + ϵ 1 y ≤ x h y -t h + x h y - x y + ϵ 1 y ≤ u - 2u 2 y + u y + 2u 2 y = u + u y < 2u.
= ϵ 1 x + ϵ 2 • y x . ( 32
)
Let us now consider two possible cases.

(1) If x ≥ y, then from Equation (32) we immediately deduce that the relative error of Algorithm 15 is bounded by

|ϵ 1 | + |ϵ 2 | ≤ 3u 2 .
(2) If x < y, which implies x h ≤ y and t h ≤ 1. The case x h = y is easily handled. It leads to t h = 1, π h = x h , π = 0, δ = x , and z h + z = t h + t = x/y + η, with |η| ≤ u |x |/y ≤ u 2 x/y. We can now focus on the case x h < y. Lemma 6.3 now implies

|t h -x h /y| ≤ u 2 - u 2 /y, so |t h y -x h | ≤ u 2 • y -u 2 ≤ u -2u 2 . Therefore |t h y -x | ≤ 2u -2u 2 , which implies |ϵ 1 | ≤ u 2 .
From Equation (32) we deduce that the relative error of Algorithm 15 is bounded by u 2 + 2u 2 = 3u 2 .

DIVISION OF TWO DOUBLE-WORD NUMBERS

The algorithm implemented in the QD library for dividing two double-word numbers is the following.

ALGORITHM 16: -DWDivDW1(x h , x , y h , y ). Calculation of (x h , x ) ÷ (y h , y ) in binary, precision-p, floating-point arithmetic.

1: t h ← RN(x h /y h ) 2: (r h , r l ) ← DWTimesFP1(y h , y , t h ) {approximation to (y h + y ) • t h using Alg. 7} 3: (π h , π ) ← 2Sum(x h , -r h ) 4: δ h ← RN(π -r ) 5: δ ← RN(δ h + x ) 6: δ ← RN(π h + δ ) 7: t ← RN(δ /y h ) 8: (z h , z ) ← Fast2Sum(t h , t ) 9: return (z h , z )
Let us quickly analyze the beginning of Algorithm 16. This will lead us to suggest another algorithm, faster yet mathematically equivalent as soon as p ≥ 3. Without loss of generality, we assume x h > 0 and y h > 0. Define ϵ x and ϵ y such that x h = x (1 + ϵ x ) and y h = y/(1 + ϵ y ). These two numbers ϵ x and ϵ y have an absolute value less than or equal to u. We have

t h = x h y h (1 + ϵ 0 ), with |ϵ 0 | ≤ u, ( 33
)
and, from Theorem 4.1,

r h + r = t h y(1 + η), with|η| ≤ 3 2 u 2 + 4u 3 . ( 34
)
There exists |ϵ 1 | ≤ u such that r h = (r h + r )(1 + ϵ 1 ). This can be rewritten r = -ϵ 1 (r h + r ), so, using Equation ( 34), r = -ϵ 1 t h y(1 + η). We finally obtain

r h = t h y h (1 + ϵ y )(1 + ϵ 1 )(1 + η) = x h (1 + ϵ y )(1 + ϵ 0 )(1 + ϵ 1 )(1 + η), ( 35
) so (1 -u) 3 (1 -2u 2 )x h ≤ r h ≤ (1 + u) 3 3 2 u 2 + 4u 3 x h ,
from which we deduce

|x h -r h | ≤ 3u + 9 2 u 2 + 19 2 u 3 + 33 2 u 4 + 27 2 u 5 + 4u 6 • x h , which implies |x h -r h | ≤ (3u + 6u 2 ) • x h ( 36
) as soon as p ≥ 3. One easily checks that for p ≥ 3 (i.e., u ≤ 1/8), 3u + 6u 2 is less than 1/2. Hence, from Sterbenz Lemma (Lemma 1.2), the number x hr h is a floating-point number. Therefore the number π obtained at line 3 of Algorithm 16 is always 0, and that line can be replaced by a simple, errorless, subtraction. This gives π h = x hr h , and δ h = -r . Hence, without changing the final result, we can replace Algorithm 16 by the simpler Algorithm 17, below.

ALGORITHM 17: -DWDivDW2(x h , x , y h , y ). Calculation of (x h , x ) ÷ (y h , y ) in binary, precision-p, floating-point arithmetic: Improved version of Algorithm 16. Useless operations have been removed. The result is exactly the same.

1: t h ← RN(x h /y h ) 2: (r h , r l ) ← DWTimesFP1(y h , y , t h ) {approximation to (y h + y ) • t h using Algorithm 7} 3: π h ← RN(x h -r h ) = x h -r h (exact operation) 4: δ ← RN(x -r ) 5: δ ← RN(π h + δ ) 6: t ← RN(δ /y h ) 7: (z h , z ) ← Fast2Sum(t h , t ) 8: return (z h , z )
If an FMA instruction is available, then Algorithm 9 can be used at line 2 instead of Algorithm 7 without changing much the error bound provided by Theorem 7.1 below. We have Theorem 7.1. If p ≥ 7, then the relative error of Algorithms 16 (DWDivDW1) and 17 is upper bounded by 15u 2 + 56u 3 .

Proof. For reasons of symmetry, we can assume that x and y are positive. We will use the results of Equations ( 33) to (36) obtained when analyzing the beginning of Algorithm 16. Assume p ≥ 7. We have δ

= (x -r )(1 + ϵ 2 ), with |ϵ 2 | ≤ u. We have |x | ≤ u • x h and |r | ≤ u • r h , so |x -r | ≤ |x | + |r | ≤ u • x h + u • r h ≤ u • x h + u • ((r h -x h ) + x h ) ≤ u • x h + u • (|r h -x h | + x h ) .
Therefore, using Equation (36),

|x -r | ≤ u • x h + u • (3u + 6u 2 )x h + x h , which gives |x -r | ≤ (2u + 3u 2 + 6u 3 ) • x h . (37)
We have δ

= (π h + δ )(1 + ϵ 3 ), with |ϵ 3 | ≤ u, so δ = x h -r h + x -r + (x -r )(ϵ 2 + ϵ 3 + ϵ 2 ϵ 3 ) + (x h -r h ) • ϵ 3 , = x -(r h + r ) + α • x h ,
with (using Equation ( 36) and ( 37))

|α | ≤ (2u + 3u 2 + 6u 3 )(2u + u 2 ) + (3u + 6u 2 )u ≤ 7u 2 + 15u 3 (38) as soon as p ≥ 4. Hence δ = x -t h y(1 + η) + αx h , so δ y h = x -t h y y • y y h - ηt h y y h + α x h y h . (39)
The number xt h y is equal to x ht h y h + xt h y . From Equation ( 33), x ht h y h is equal to -x h ϵ 0 . Also, |x | is less than or equal to ux h , and

|t h y | ≤ |ut h y h | ≤ u (1 + u)x h . Hence, |x -t h y| ≤ x h • (u + u + u (1 + u)) = x h • (3u + u 2 ). (40)
From Equation (39), we deduce

δ y h = x -t h y y • (1 + ϵ y ) -ηt h (1 + ϵ y ) + α x h y h , = x -t h y y + β, (41)
with

|β | = ϵ y • x -t h y y -t h (1 + ϵ y )η + x h y h ≤ u (3u + u 2 ) x h y + (1 + u)(2u 2 ) x h y h + (7u 2 + 15u 3 ) x h y h ≤ u (3u + u 2 )(1 + u) x y + (1 + u) 3 (2u 2 ) x y + (7u 2 + 15u 3 )(1 + u) 2 x y = (12u 2 + 39u 3 + 44u 4 + 17u 5 ) • x y . (42)
Hence,

t = RN δ y h = δ y h (1 + ϵ 4 ) with |ϵ 4 | ≤ u, = x -t h y y + β (1 + ϵ 4 ) = x -t h y y + γ , (43)
with

|γ | = x -t h y y ϵ 4 + β + ϵ 4 β ≤ x h y (3u + u 2 )u + β + βu ≤ (3u + u 2 )u (1 + u) x y + β + βu = (15u 2 + 55u 3 + 84u 4 + 61u 5 + 17u 6 ) • x y . (44)
Hence,

t h + t = x y + γ .
Since we straightforwardly have

t h ≥ x y • (1 -u) 3 , (45)
we deduce

|t | ≤ x y • (15u 2 + 55u 3 + 84u 4 + 61u 5 + 17u 6 ) + (3u -3u 2 + u 3 ) . (46)
From Equations ( 45) and ( 46) we easily deduce that as soon as p ≥ 4 (i.e., u ≤ 1/16), t h is larger than |t |, so Algorithm Fast2Sum introduces no error at line 7 of the algorithm. Therefore,

z h + z = t h + t = x y + γ ,
so the relative error of Algorithm 17 (and Algorithm 16) is upper bounded by 15u 2 + 55u 3 + 84u 4 + 61u 5 + 17u 6 , which is less than 15u 2 + 56u 3 as soon as p ≥ 7 (i.e., u ≤ 128), which always holds in practice.

The bound provided by Theorem 7.1 is almost certainly not optimal. However, during our intensive tests, we have encountered cases for which the relative error, although significantly less than the bound 15u 2 + 56u 3 of Theorem 7.1, remains of a similar order of magnitudethat is, more than half the bound. For instance, for p = 53, relative error 8.465 • • • × 2 -106 is attained for x h = 4503607118141812, x = 4493737176494969/2 53 , y h = 4503600552333684, and y = -562937972998161/2 50 .

If an FMA instruction is available, then one can design a more accurate algorithm. What makes it work is Property 6.1, applied in the special case x = 1.

ALGORITHM 18: -DWDivDW3(x h , x , y h , y ). Calculation of (x h , x ) ÷ (y h , y ) in binary, precision-p, floating-point arithmetic: more accurate algorithm that requires the availability of an FMA instruction Let us now consider y 2 (t h -1/y) . That term is less than

1: t h ← RN(1/y h ) 2: r h ← (1 -y h t h ) = 1 -y h t h (exact operation) 3: r ← -RN(y • t h ) 4: (e h , e ) ← Fast2Sum(r h , r ) 5: (δ h , δ ) ← DWTimesFP3(e h ,
y 2 t h - 1 y h + y -y h yy h 2 ,
which is less than

y 2 u 2 1 2 + 1 y(y -u) 2 .
The largest value of

y 2 1 2 + 1 y(y -u) 2
for 1 ≤ y < 2 is always attained for y = 1, so as soon as p ≥ 6 (i.e., u ≤ 1/64), we have

y 2 t h - 1 y ≤ 1 2 + 1 1 -1 64 2 u 2 = 36481 15876 u 2 ≤ 2.298u 2 .
Hence, from Equation ( 53), we obtain This relative error bound is certainly a large overestimate, since we cumulate in its calculation the overestimates of the errors of Algorithms 9, 4, and 12. In practice, Algorithm 18 is rather accurate: the largest relative error found so far in our tests for p = 53, is 5.922 • • • × 2 -106 , obtained for x h = 4528288502329187, x = 1125391118633487/2 51 , y h = 4522593432466394, and y = -9006008290016505/2 54 . For each algorithm, we give the previously known bound (when we are aware of it, and when the algorithm already existed), the bound we have proved, the largest relative error observed in our fairly intensive tests, and the number of floating-point operations required by the algorithm.

(m h + m ) - 1 y ≤ 1 y • 2.

CONCLUSION

We have proven relative error bounds for several basic building blocks of double-word arithmetic, suggested a new algorithm for multiplying two double-word numbers, suggested an improvement of the algorithms used in the QD library for dividing a double-word number by a floating-point number and for dividing two double-word numbers. We have also suggested a new algorithm for dividing two double-word numbers when an FMA instruction is available. Table 1 summarizes the obtained results. For the functions for which an error bound was already published, we always obtain a significantly smaller bound, except in one case, for which the previously known bound turned out to be slightly incorrect. Our results make it possible to have more trust in double-word arithmetic. They also allow us to give some recommendations in what follows.

• For adding two double-word numbers, never use Algorithm 5, unless you are certain that both operands have the same sign. Double-word numbers can be added very accurately using the (unfortunately more expensive) Algorithm 6. • For multiplying a double-word number by a floating-point number, Algorithm 8 is less accurate, yet slightly faster, than Algorithm 7. Hence one cannot say that one is really better than the other one. Choose between them depending on whether you mainly need speed or accuracy. If an FMA instruction is available, then Algorithm 9 is a good candidate.

• For multiplying two double-word numbers, if an FMA instruction is available, then Algorithm 12 is to be favored. It is more accurate both from a theoretical (better error bound) and from a practical (smaller observed errors in our intensive testings) points of view. • There is no point in using Algorithm 13 for dividing a double-word number by a floatingpoint number: Algorithm 15 is faster and has a better error bound. • There is no point in using Algorithm 16 for dividing two double-word numbers:

Algorithm 17, presented in this article, always returns the same result and is faster. If an FMA instruction is available, depending whether the priority is speed or accuracy, then one might prefer Algorithm 18. It is almost certainly significantly more accurate (although we have no full proof of that: We can just say that our bounds are smaller, as well as the observed errors); however, it is slower.

  Lemma 2.1 (seeProperty (2.16) in[START_REF] Rump | Accurate floating-point summation part I: Faithful rounding[END_REF]). Let a and b be FP numbers, and let s = RN(a + b). If s 0, then |s | l.o.g., assume |a| ≥ |b |, so ulp(a) ≥ ulp(b). The number |a + b | is the distance between a and -b. Hence, since a -b (otherwise s would be 0), |a + b | is larger than or equal to the distance between a and the FP number nearest a, which is larger than or equal to 1 2 ulp(a). Therefore |RN(a + b)| = RN(|a + b |) ≥ RN( 1 2 ulp(a)) = 1 2 ulp(a). Let us now turn to the analysis of Algorithm 4. We have the following. Theorem 2.2. The relative error (z h + z ) -(x + y) x + y of Algorithm 4 (DWPlusFP) is bounded by 2 • u 2 .

  Since |x + s | ≤ 3u (see the two cases considered below), we have |v | ≤ 3u, so s h > |v |: Algorithm Fast2Sum introduces no error at line 3 of the algorithm. Therefore, z h + z = s h + v = x + y + ϵ, and the relative error of algorithm 4 is |ϵ/(x + y)|.• If x h + y ≤ 2, then |s | ≤ u, so |x + s | ≤ 2u, hence |ϵ | ≤ u 2 . If x + y ≥ 1/2,then this immediately implies that the relative error is less than 2u 2 . Now, if x + y < 1/2, then x h + y < 1/2 + u, so x h /2 < 1/2 + u. The only solution compatible with the range of x h is x h = 1. In such a case x = x h + x ≥ 1u/2, hence we must have y < -1/2 + u/2, and, hence, since

	y is a floating-point number, we must have y ≤ -1/2, which is not compatible with the
	assumption -x h /2 < y.
	• If x h + y > 2, then |s | ≤ 2u, so |x + s | ≤ 3u, hence |ϵ | ≤ 2u 2 and the relative error |ϵ |/|x +
	y| of the calculation is bounded by

  and |s | ≤ 2u. We have t h + t = x + y , with |x + y | ≤ 2u, hence |t h | ≤ 2u, and |t | ≤ u 2 . Now, |s + t h | ≤ 4u, so |c | ≤ 4u, and |ϵ 1 | ≤ 2u 2 . Since s h ≥ 2 -4u and |c | ≤ 4u, if p ≥ 3, then Algorithm Fast2Sum introduces no error at line 4 of the algorithm. Therefore,

  which implies |c h | ≤ 2, then we easily improve on some of the previously obtained bounds. We have |c 1 | ≤ u, andt h ≤ RN(2 + 2u -2u 2 ) = 2.The case t h = 2 is easily handled: Equation (22) implies xy = t h + t 2ϵ ≥ 2 -4u -3u 2 , and the relative error |ϵ/(xy)| is bounded by 3u

)

2. If x h y < 2,

  so |t 1 | ≤ 2u -2u 2 and |ϵ 1 | ≤ u 2 ; andt 2 = x y h + ϵ 2 , with |x y h | ≤ 2u -2u 2 , so |t 2 | ≤ 2u -2u 2 and |ϵ 2 | ≤ u 2 . Now, we have c 2 = t 1 + t 2 + ϵ 3 , with |t 1 + t 2 | ≤ 4u -4u 2 , which implies |c 2 | ≤ 4u -4u 2 and |ϵ 3 | ≤ 2u 2 .We finally obtainc 3 = c 1 + c 2 + ϵ 4 ,and, from |c 1 + c 2 | ≤ 6u -4u 2 , we deduce |c 3 | ≤ 6u and |ϵ 4 | ≤ 4u 2 . Since c h ≥ 1, Algorithm Fast2Sum introduces no error at line 6 of the algorithm. Therefore,

  e , t h ) {Approximation to (e h + e ) • t h with relative error ≤ 2u 2 using Algorithm 9} 6:(m h , m ) ← DWPlusFP(δ h , δ , t h ){Approximation to δ h + δ + t h with relative error ≤ 2u 2 using Algorithm 4} 7: (z h , z ) ← DWTimesDW2(x h , x , m h , m ){Approximation to (x h + x )(m h + m ) with relative error≤ 5u 2 using Algorithm 12}

8: return (z h , z ) ACM Transactions on Mathematical Software, Vol. 44, No. 2, Article 15res. Publication date: October 2017.

  + 9u 3 . Now, from Theorem 5.4, we have|z h + zx (m h + m )| ≤ 5u 2 |x (m h + m )|

								298u 2 + t h	5 2	u 2 + 9u 3 ,
	which implies
	x (m h + m ) -	x y	≤	x y	• 2.298u 2 + xt h	5 2	u 2 + 9u 3 .
								x y	≤	x y	• φ(u),	(54)
	with φ(u) = 2.298u 2 + (1 + u) 2 5 2 u 2 ≤ 5u 2 x y	+ 5u 2 x y	-x (m h + m )	(55)
								≤	x y	5u 2 + 5u 2 φ(u) .
	Combining Equations (54) and (55), we finally obtain
	z h + z -	x y	≤	x y	(5u 2 + φ(u) + 5u 2 φ(u))
			≤	x y		9.798u 2 + 14u 3 + 44.49u 4 + 79u 5 + 102.5u 6 + 45u 7
			≤ 9.8u 2 x y	as soon as p ≥ 13.

Notice that |t

h | ≤ (1 + u)/y h ≤ (1 + u) 2 /y, so x (m h + m ) -

Table 1 .

 1 Summary of the Results Presented in This Paper

					Largest
			Previously		relative error
			known		observed in	of FP
	Operation	Algorithm	bound	Our bound experiments	ops
	DW + FP	Algorithm 4	?	2u 2	2u 2 -6u 3	10
	DW + DW Algorithm 5	N/A	N/A	1	11
		Algorithm 6	2u 2 (incorrect) 3u 2 + 13u 3	2.25u 2	20
	DW × FP	Algorithm 7	4u 2	3 2 u 2 + 4u 3	1.5u 2	10
		Algorithm 8	?	3u 2	2.517u 2	7
		Algorithm 9	N/A	2u 2	1.984u 2	6
	DW × DW Algorithm 10 11u 2	7u 2	4.9916u 2	9
		Algorithm 11 N/A	6u 2	4.9433u 2	8
		Algorithm 12 N/A	5u 2	3.936u 2	9
	DW ÷ FP	Algorithm 13 4u 2	3.5u 2	2.95u 2	16
		Algorithm 14 N/A	3.5u 2	2.95u 2	10
		Algorithm 15 N/A	3u 2	2.95u 2	10
	DW ÷ DW Algorithm 16 ?	15u 2 + 56u 3 8.465u 2	24
		Algorithm 17 N/A	15u 2 + 56u 3 8.465u 2	18
		Algorithm 18 N/A	9.8u 2	5.922u 2	31

ACM Transactions on Mathematical Software, Vol. 44, No. 2, Article 15res. Publication date: October 2017.

(y, x ) may not be a double-word number, according to Definition 1.4, in the case x = 1

ulp(y ) = 1 2 ulp(x h ). However, one easily checks that in that case the algorithm returns an exact result.ACM Transactions on Mathematical Software, Vol. 44, No. 2, Article 15res. Publication date: October 2017.

The improvement was suggested by one of the anonymous reviewers. We are very grateful for that.ACM Transactions on Mathematical Software, Vol. 44, No. 2, Article 15res. Publication date: October 2017.

ACKNOWLEDGMENTS

We thank the anonymous reviewers for their very detailed and helpful comments and suggestions.

This work was supported by the FastRelax (ANR-14-CE25-0018-01) project of the French National Agency for Research (ANR).

We have the following. Theorem 7.2. As soon as p ≥ 13, and if y 0, the relative error of Algorithm 18 (DWDivDW3) is bounded by 9.8u 2 .

Proof. Roughly speaking, Algorithm 18 first approximates 1/y by t h = RN(1/y h ), then improves that approximation to 1/y by performing one step of Newton-Raphson iteration, and then multiplies the obtained approximation (m h , m ) by x.

Without loss of generality, we assume 1 ≤ y h ≤ 2 -2u, so 1/2 ≤ t h ≤ 1. We have

and (from Property 6.1)

Also, since

and from Theorem 2.2 we have

Combining Equations ( 49) and ( 50), we obtain

Therefore,

so, using Equation ( 47) and the bounds on η and α, (53)