Tom Hubrecht

Claude-Pierre Jeannerod

Jean-Michel Muller

Useful applications of correctly-rounded operators of the form ab + cd + e

Keywords: Floating-point arithmetic, correct rounding, fused operators, dot-product instructions, complex arithmetic, double-word arithmetic, rounding error analysis

We show that the availability of fused arithmetic operators that evaluate expressions of the form ab + cd (FD2 instruction) or ab + cd + e (FD2A instruction) in floatingpoint arithmetic with one final rounding only would significantly facilitate many calculations that are hard to perform with high accuracy at small cost using only the traditional operations +, -, ×, ÷, √ , and fused multiply-add (FMA).

I. INTRODUCTION

A. Fused dot-product operators

The recent years have seen the emergence of new arithmetic instructions, that allow to perform fused "exact" dotproduct (or dot-product-add) calculations for low-dimensional vectors: the arithmetic operations involved are fused, with one final rounding only. Most frequently, these instructions correspond to some types of dot products, such as ab + cd or ab + cd + e or longer [START_REF] Tao | Correctly rounded architectures for floating-point multi-operand addition and dot-product computation[END_REF], [START_REF] Brunie | Modified fused multiply and add for exact low precision product accumulation[END_REF], [START_REF] Bertaccini | MiniFloat-NN and ExSdotp: An ISA extension and a modular open hardware unit for low-precision training on RISC-V cores[END_REF], [START_REF] Desrentes | Exact fused dot product add operators[END_REF]. Typical applications of such instructions are the implementation of ML algorithms, but the FFT is also a target [START_REF] Dupont De Dinechin | In-Place Multicore SIMD Fast Fourier Transforms[END_REF]. The predecessor of these instructions is the fused multiply-add (FMA), which evaluates expressions of the form ab + c with only one final rounding. The first widely available processor with an FMA was the IBM POWER [START_REF] Montoye | Design of the IBM RISC System/6000 floating-point execution unit[END_REF]. In addition to speeding up (and generally making more accurate) the computation of dot products and the evaluation of polynomials, the availability of an FMA instruction allows efficient software implementation of correctly rounded division and square root [START_REF] Markstein | Computation of elementary functions on the IBM RISC System/6000 processor[END_REF], [START_REF] Cornea-Hasegan | Correctness proofs outline for Newton-Raphson based floating-point divide and square root algorithms[END_REF], and makes it easy to compute the error of a floating-point multiplication, paving the way for efficient double-word multiplication and better complex arithmetic [START_REF] Jeannerod | Further analysis of Kahan's algorithm for the accurate computation of 2 × 2 determinants[END_REF]. The goal of this paper is to show that these new fused instructions can bring to the field of numerical computing a progress comparable to the progress brought by the FMA 30 years ago.

So far, the fused exact dot-product instructions presented in the literature are mainly targeted at low-precision arithmetic. Furthermore, they frequently use different formats at the same time: typically, the products are accumulated in a format larger than the format of the input operands. Brunie [START_REF] Brunie | Modified fused multiply and add for exact low precision product accumulation[END_REF] considers 16-bit inputs and an 80-bit accumulator, Bertaccini et al. [START_REF] Bertaccini | MiniFloat-NN and ExSdotp: An ISA extension and a modular open hardware unit for low-precision training on RISC-V cores[END_REF] present such a sum-of-dot-product operator for NN training and inference, which takes four inputs a, b, c, d in a w-bit format (with w = 8 or 16) and one input e in a 2w-bit format, and outputs ab + cd + e rounded to the 2w-bit format. Desrentes et al. [START_REF] Desrentes | Exact fused dot product add operators[END_REF] investigate the general case of 2n + 1 input operands with n ≥ 2. Kalray's MPPA architectures offer an FDMDA 1 (fused dual-multiply dual-add) instruction, that computes ab+cd+e correctly rounded in binary32 arithmetic.

Lauter [START_REF] Lauter | An efficient software implementation of correctly rounded operations extending FMA: a + b + c and a × b + c × d[END_REF] describes software implementations of the fused, correctly-rounded operations a + b + c and ab + cd, and lists a few applications: accurate complex multiplication, faster compensated algorithms, simplification of the last rounding step in correctly rounded elementary functions, and accurate computation of discriminants of quadratic polynomials.

Several other fused sum or dot-product operators have been introduced in the literature, especially for ML applications; see for example [START_REF] Sohn | A fused floating-point four-term dot product unit[END_REF], [START_REF] Hickmann | Experimental analysis of matrix multiplication functional units[END_REF], [START_REF] Hickmann | Intel Nervana Neural Network Processor-T (NNP-T) fused floating point many-term dot product[END_REF]. They are of little use for our study, either because they do not provide correct rounding or because the way they round the results is not fully specified.

B. Our contributions

In this paper, we consider only homogeneous operators (that is, operators whose input and output operands are all represented in the same format) that evaluate expressions of the form ab+cd or ab+cd+e with correct rounding to nearest. We show that these correctly-rounded fused operators allow for simple and accurate implementation of complex arithmetic (section II); that they simplify the evaluation of some rounding errors, either exactly (error-free transforms) or with correct rounding, paving the way for efficient compensated algorithms (section III); that they allow for more accurate computations of products of three or four numbers, with an application to the accurate computation of integer powers (section IV); that they simplify the evaluation of the sign of the discriminant of some quadratic and cubic equations (section V) as well as the implementation of correctly-rounded multiplication by constants such as π (section VI); that they allow for implementing sums and inner products with fewer instructions and smaller error bounds than the classical algorithms (section VII); and that they allow for evaluating polynomials as accurately as the classical and compensated FMA-based Horner schemes but with significantly fewer instructions (section VIII). Finally, we show that such operators allow for simple and efficient implementations of double-word arithmetic (section IX).

For each of these eight application domains, the possible gains are quantified precisely and systematically via detailed operation counts and/or the derivation of rigorous, worst-case accuracy bounds (several of them being shown to be tight).

C. Notation and background

We assume radix-2, precision-p, floating-point (FP) arithmetic, with an unbounded exponent range (which means that the algorithms and properties presented below apply to an arithmetic conforming to the IEEE 754 standard [START_REF]IEEE, IEEE Standard for Floating-Point Arithmetic (IEEE Std[END_REF] as long as underflows and overflows do not occur).

We write F to denote the set of FP numbers. We also assume that the rounding function, denoted RN, is round to nearest with ties to even, which is the default in IEEE 754 arithmetic. The associated unit roundoff is then u = 2 -p .

Given a, b, c, d, e ∈ F, we assume that the following two instructions are available for performing fused dot-products of order two (FD2), possibly followed by an addition (FD2A):

• FD2(a, b, c, d) returns RN(ab + cd); • FD2A(a, b, c, d, e) returns RN(ab + cd + e).
The necessary background in FP arithmetic and error analysis can be found for example in [START_REF] Boldo | Floatingpoint arithmetic[END_REF]. In particular, recall that for x ∈ R, ufp(x) = 2 ⌊log2 |x|⌋ if x ̸ = 0 (and zero otherwise), that ulp(x) = 2u • ufp(x), and that the absolute and relative errors due to rounding to nearest are bounded, respectively, as |RN(x) -x| ≤ 1 2 ulp(x) and as RN(x) = x(1 + δ), |δ| ≤ u/(1 + u). Recall also that if x ∈ F, then it is an integral multiple of its ulp: x ∈ ulp(x)Z.

II. COMPLEX ARITHMETIC

For complex FP numbers a + ib and c + id, addition and subtraction can be implemented perfectly well with basic FP arithmetic, as RN(a ± c) + i RN(b ± d). However, for complex multiplication and division, even FMA-based implementations can suffer heavy cancellation in the real or imaginary part of the result, and so an FD2 instruction will be most useful here to get highly-accurate results.

Specifically, for complex products (a + ib)(c + id) with a, b, c, d ∈ F, two FD2 calls yield the correctly-rounded real and imaginary parts RN(ac -bd) and RN(ad + bc). When c+id = a-ib (complex multiplication by complex conjugate), this ensures that the imaginary part of the result is zero; with an FMA, this algebraic property can be lost, since RN(-ab + RN(ab)) ̸ = 0 unless ab ∈ F.

For complex division, we have (a + ib)/(c + id) = R + iI with R = ac+bd c 2 +d 2 and I = bc-ad c 2 +d 2 , so that three FD2 calls and two FP divisions suffice to get approximations R and I having high relative accuracy. More precisely, R = RN(RN(ac + bd)/RN(c 2 + d 2 )), where the relative error of FD2 is at most u/(1 + u) and the one of FP division is at most u -2u 2 [START_REF] Jeannerod | On relative errors of floating-point operations: optimal bounds and applications[END_REF]; the same holds for I, and it is then easily checked that this implies

R = R(1 + δ R ), I = I(1 + δ I ), |δ R |, |δ I | ≤ 3u.
Another immediate consequence of the availability of an FD2 instruction is a more accurate complex absolute value √ a 2 + b 2 (i.e., hypotenuse): with or without an FMA, the relative error is known to be tightly bounded by 2u [START_REF] Jeannerod | On relative errors of floating-point operations: optimal bounds and applications[END_REF], [START_REF] Jeannerod | The classical relative error bounds for computing √ a 2 + b 2 and c/ √ a 2 + b 2 in binary floatingpoint arithmetic are asymptotically optimal[END_REF]; with an FD2, we now compute h := RN RN(a 2 + b 2 ) . By applying twice the relative error bound u/(1 + u), we deduce that

h = a 2 + b 2 (1 + δ h ), |δ h | < 3 2 u. Furthermore, this new bound is tight: taking a = 1 and b = RN( √ 3u(1 -u)) gives |δ h | > 1 -1/ √ 1 + 3u -12u 2 for p ≥ 3, which is 3u/2 -O(u 2 ) as u → 0.
Note also that having an FD2A (and not just an FD2) gives directly a complex fused multiply-add operation, which is for example key for the FFT [START_REF] Desrentes | Exact fused dot product add operators[END_REF]: since (a + ib)(c + id) + e + if = ac -bd + e + i(ad + bc + f ), two FD2A calls suffice to get the correctly-rounded values of the real and imaginary parts of this complex number.

Finally, another advantage of these FD2(A)-based solutions is to allow to avoid spurious underflows and overflows, either completely (for multiplication and complex FMA) or partially (for division and absolute value).

III. ERROR-FREE TRANSFORMS

It is well known that under mild conditions the error of an FP multiplication and the error of a (rounded to nearest) FP addition are FP numbers, and that these errors can be computed using the 2Sum, Fast2Sum, and 2MultFMA algorithms; see for example [15, §4.2]. These algorithms are called errorfree transformations (EFT) in [START_REF] Ogita | Accurate sum and dot product[END_REF].

The 2Sum algorithm requires 6 FP additions and, assuming some kind of ordering on the input, the Fast2Sum algorithm needs only 3 FP additions. Without an FMA, computing the exact error of an FP product can be done using 17 FP operations [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF], but with an FMA this becomes much easier and cheaper: if one computes P h = RN(ab) and P ℓ = RN(ab -P h ), then P h + P ℓ = ab exactly. 2 In other words, the FMA makes it possible to implement the EFT of multiplication with only two operations.

Similarly, as noted in [START_REF] Ogita | Accurate sum and dot product[END_REF], the availability of an ADD3 operation, which for a, b, c ∈ F returns RN(a + b + c), would greatly simplify the computation of the error of an FP addition: if S h = RN(a + b) and S ℓ = RN(a + b -S h ), then S h + S ℓ = a + b exactly.

Since the FMA and ADD3 are special cases of the FD2 and FD2A operations, respectively, having an FD2 makes it easy to obtain the error of an FP multiplication, and having an FD2A makes it easy to obtain the error of an FP addition: in each case, two operations now suffice to implement the EFT.

What about the error of the FMA ab + c -RN(ab + c)? As shown in [START_REF] Boldo | Some functions computable with a Fused-Mac[END_REF], [START_REF]Exact and approximated error of the FMA[END_REF], this error is in general not an FP number but can always be expressed as the exact sum e 1 + e 2 of two FP numbers such that e 1 = RN(e 1 + e 2 ); furthermore, the pair (e 1 , e 2 ) can be produced using only basic FP operations, but the cost is high: Algorithm ErrFmac from [START_REF] Boldo | Some functions computable with a Fused-Mac[END_REF] uses 2 FMA, 1 multiplication, and 17 additions, for a total of 20 basic FP operations. Interestingly, with an ADD3 (and thus a fortiori with an FD2A), this operation count drops to 10: by inspecting Algorithm ErrFmac and its proof in [START_REF] Boldo | Some functions computable with a Fused-Mac[END_REF], one sees that 2 FMA, 1 multiplication, 3 additions, and 4 ADD3 are sufficient.

If one needs only the correctly-rounded value of the error of the FMA, that is,

e 1 = RN(ab + c -RN(ab + c)), (1)
then Algorithm ErrFmaNearest from [START_REF]Exact and approximated error of the FMA[END_REF] can be used: it performs 18 basic FP operations and, again, it is easily seen that with an ADD3 this operation count is halved too, dropping to 9. However, in this case, having an FD2A makes an even bigger difference, since e 1 in (1) can now be obtained using only 2 operations: one FMA followed by one FD2A. We will see in section VIII how this can be used to design a new, FMA/FD2A-based, compensated Horner algorithm, which avoids using the expensive EFT of the FMA without sacrificing for accuracy.

Similarly, only one call to the FD2A makes it possible to obtain the correct rounding of the error of an FD2 (and thus the correct rounding of the errors in a complex product):

P h = RN(ab + cd), P ℓ = RN(ab + cd -P h ). (2)
As for the FMA, we will have P h + P ℓ ̸ = ab + cd in general. However, the theorem below shows that the scheme in (2) does define an EFT for ab + cd provided some conditions on the ratio ab/cd are satisfied.

Theorem 1. Given a, b, c, d ∈ F, let P h and P ℓ be as in ( 2). If ab/2 ≤ -cd ≤ 2ab, then ab + cd = P h + P ℓ .

Proof. We will use two basic properties: for x, y, z, t ∈ R,

(i) if |x| ≤ 2 k then |RN(x) -x| ≤ 2 -p-1 • 2 k ; (ii) if |xy| ≤ |zt| then ulp(x)ulp(y) ≤ 2ulp(z)ulp(t).
Since ab and cd play a symmetric role, we can assume that -cd = |cd| ≤ ab, so 0

≤ ab + cd ≤ |cd| ≤ 2 2p ulp(c)ulp(d). Using (i) we deduce that |ab + cd -P h | ≤ 2 p-1 ulp(c)ulp(d). On the other hand, ab + cd ∈ µZ with µ := min{ulp(a)ulp(b), ulp(c)ulp(d)}. Since |cd| ≤ ab, it follows from (ii) that µ ≥ 1 2 ulp(c)ulp(d) =: 2 e ,
so that ab+cd and P h are in 2 e Z. Hence |ab + cd -P h | is an integral multiple of 2 e that is less than or equal to 2 p+e . This implies ab+cd-P h ∈ F and thus

P ℓ = RN(ab + cd -P h ) = ab + cd -P h .
We conclude this section by noting that even when the exact sum P h + P ℓ differs from the exact value ab + cd, its relative error is extremely small, as we always have the following bound, which is asymptotically optimal:

P h + P ℓ = (ab + cd)(1 + δ), |δ| ≤ 1 2 u 2 .
This bound and its asymptotic optimality follow as a special case of the more general property below (which will also be useful in order to establish Property 2 in section IV, and when designing double-word arithmetic algorithms in section IX).

Property 1. Let x ∈ R, x = RN(x)
, and e = RN(x -x).

Then the relative error of the exact sum x + e is at most 1 2 u 2 . Furthermore, this bound is asymptotically optimal.

Proof. The bound is clear for x = 0, so let x ̸ = 0. Writing e = x -x, we have |e| ≤ u ufp(x) =: B. This bound B is an integral power of 2, so

| x + e -x| = |RN(e) -e| ≤ u 2 B = 1 2 u 2 ufp(x).
Using |x| ≥ ufp(x) then leads to the relative error bound. Asymtpotic optimality can be shown by considering

x = 1 + u + 1 2 u 2 , for which x = 1 + 2u, e ∈ {-u, -u + u 2 }, and |( x + e)/x -1| = u 2 /2 1+u+ 1 2 u 2 ∼ 1
2 u 2 as u → 0. IV. PRODUCTS OF THREE OR FOUR FP NUMBERS Note first that given a, b, c, d ∈ F, one can compute the correctly-rounded value of abc + d very easily with an FD2A: compute p h := RN(ab), p ℓ := RN(ab -p h ), and q h := RN(p h c + p ℓ c + d); since p h + p ℓ = ab, we see immediately that q h satisfies q h = RN(abc + d).

(

) 3
This has several interesting properties. A first one (detailed in section V) is to ease the correct computation of the sign of the discriminant of some cubic equations. Evaluating expressions of the form RN(abc + d) can also be useful for elementary function libraries. For example, the cosine of a small argument x is typically approximated by 1 + x 2 P (x) with P a low-degree polynomial.

Of course, (3) leads in particular to correctly-rounded products of three FP numbers. But what if we want to multiply four FP numbers together? Interestingly, the FD2A can again be used easily in order to get an approximation to abcd with relative error at most u, as follows: by applying the method underlying (3), compute q h := RN(abc) and q ℓ := RN(abc -q h ), and then deduce r := RN(q h d + q ℓ d).

(4)

Note that since abc -q h is not always in F, the exact sum q h +q ℓ can differ from abc, and so r may not be the correctlyrounded value of abcd. Such a situation occurs for example when (a, b, c, d) = (1 + 2u, 1 + 2u, 1 -u, 1): in this case, one can check that for u ≤ 2 -3 , q h = 1 + 2u and r = 1 + 4u. Despite this, the relative accuracy of r in ( 4) is essentially best possible, as the next property shows.

Property 2. For r as in (4), | r -abcd| ≤ u|abcd|.

Proof. By Property 1, q h +q ℓ = abc(1+δ) with |δ| ≤ 1 2 u 2 ; on the other hand, r = (q

h d + q ℓ d)(1 + δ ′ ) with |δ ′ | ≤ u/(1 + u). Hence r = abcd(1 + δ ′′ ) with δ ′′ = (1 + δ)(1 + δ ′ ) -1, which satisfies |δ ′′ | ≤ u for u ≤ 1/4.
For comparison, when evaluating abcd in the most usual way with three FP multiplications, the associated relative error bound is about 3u. Longer products, of the form x 1 x 2 • • • x n with n ≥ 5, will obviously benefit from the above lemma in about the same way, with the usual bound ≈ (n -1)u (see for example [START_REF] Rump | Improved error bounds for floating-point products and Horner's scheme[END_REF]) being replaced by a smaller one

≈ n 4 u + ( n 4 -1)u = nu/2.

V. QUADRATIC AND CUBIC EQUATIONS

When solving quadratic or cubic equations in FP arithmetic, a first key step is to get the sign of their discriminant right, as this tells about the nature of the solutions and in particular the number of real solutions.

For quadratic equations ax 2 +bx+c = 0 with a, b, c ∈ F, the discriminant is ∆ = b 2 -4ac and, as already noted in [START_REF] Lauter | An efficient software implementation of correctly rounded operations extending FMA: a + b + c and a × b + c × d[END_REF], one FD2 instruction suffices to produce RN(∆) and thus sign(∆) as well. (This holds at least for binary arithmetic, since 4a ∈ F is obtained exactly as RN(4a).)

By contrast, other ways to evaluate ∆, such as computing ∆ +× := RN(RN(b 2 ) -RN(4ac)), ∆ FMA 1 := RN(RN(b 2 ) -4ac), and ∆ FMA 2 := RN(b 2 -RN(4ac)) can be so wrong that even the sign is not obtained correctly [START_REF] Kahan | Lecture notes on the status of IEEE-754[END_REF]. For example, one has ∆ +× = 0 and ∆ ̸ = 0 when (a, b, c) = (1/4 -u/2, 1, 1 + 2u). One can also have ∆ FMA 1 < 0 or ∆ FMA 2 < 0, while ∆ = 0, which leads to conclude wrongly that there is no real solution; an example for

∆ FMA 1 is given by (a, b, c) = (1/4 -u/4, 1 -u, 1 -u).
An FD2 instruction can also be immediately useful when solving so-called depressed cubic equations x3 + bx + c = 0 with b, c ∈ F. Here the discriminant is ∆ = -4b 3 -27c 2 and its sign can be computed wrongly by direct implementations, with or without FMA. For example, in binary32 arithmetic with b = -14628329 and c = 21534767104, we have

∆ FMA < 0 < ∆ with ∆ FMA := RN(-RN(4b • RN(b 2 )) -c • RN(27c)).
With an FD2 instruction, however, one can use the solution presented in section IV to easily obtain the correctly-rounded products δ 1 = RN(4b 3 ) and δ 2 = 2 ), and then compute ∆ FD2 = RN(-δ 1 -δ 2 ). By monotonicity of rounding, this approach ensures that ∆ FD2 > 0 implies ∆ > 0 (and that there are three distinct real roots), and that ∆ FD2 < 0 implies ∆ < 0 (and that there is one real root). Hence only the case ∆ FD2 = 0 will require more work in order to compute sign(∆) correctly.

VI. CORRECTLY-ROUNDED MULTIPLICATION BY REAL

CONSTANTS

In [START_REF] Brisebarre | Correctly rounded multiplication by arbitrary precision constants[END_REF], the authors consider the problem of performing the correctly rounded (to nearest) multiplication by a constant C / ∈ F. Typical useful values of C are π, 1/π, ln(2), √ 2, etc. Assuming that C h = RN(C) and C ℓ = RN(C -C h ) are precomputed, for an input x ∈ F they assume that y * = C • x is approximated by

y 1 = RN(C h • x + RN(C ℓ • x)), (5)
whose computation requires an FP multiplication followed by an FMA. There are values of C and x for which y ̸ = RN(y * ). However, it is shown in [START_REF] Brisebarre | Correctly rounded multiplication by arbitrary precision constants[END_REF] that for a given C, a simple algorithm can check if y 1 = RN(y * ) for all x ∈ F, and it is observed that this is indeed the case for many useful values of C. The underlying reasoning is as follows:

• First, there are ≤ 4 values of the significand of x for which y 1 and y * may belong to different binades. One easily checks if y 1 = RN(y * ) for these values.

• Then, for all other numbers x and with ϵ := |C -

(C h + C ℓ )|, one can bound |y 1 -y * | by 1 2 ulp(y * ) + η 1 , η 1 := 1 2 ulp(C ℓ x) + ϵ|x|. (6)
• Finally, [START_REF] Montoye | Design of the IBM RISC System/6000 floating-point execution unit[END_REF] implies that if y 1 ̸ = RN(y * ), then Cx is within distance ≤ η 1 from the middle of two consecutive FP numbers. This in turn implies that C is very close to a rational number of the form (2A + 1)/(2X) with A and X two integers between 2 p-1 and 2 p , which can be checked using the continued-fraction expansion of C.

Method 3 in [START_REF] Brisebarre | Correctly rounded multiplication by arbitrary precision constants[END_REF] implements this strategy, and either tells that there is no x such that y 1 ̸ = RN(y * ), or gives all values x (up to a multiplication by a power of 2) for which y 1 ̸ = RN(y * ).

With an FD2 instruction we can change the way the approximation to y * is computed and, instead of (5), use

y 2 = RN(C h • x + C ℓ • x).
We easily obtain that there are at most 2 values of the significand of x for which y and y * may belong to different binades (as previously, one checks these values separately), and that for all other values we have [START_REF] Brisebarre | Correctly rounded multiplication by arbitrary precision constants[END_REF] tells us there is no x such that y 1 ̸ = RN(y * ), then there is no x such that y 2 ̸ = RN(y * ) either. If Method 3 yields x such that y 1 ̸ = RN(y * ), then we must check if y 2 = RN(y * ). Using this strategy we can show that with π, 1/π, ln(2), 1/ ln(2), ln [START_REF] Lauter | An efficient software implementation of correctly rounded operations extending FMA: a + b + c and a × b + c × d[END_REF], 1/ ln [START_REF] Lauter | An efficient software implementation of correctly rounded operations extending FMA: a + b + c and a × b + c × d[END_REF], in the binary32, binary64, and binary128 formats we always have y 2 = RN(y * ).

|y 2 -y * | ≤ 1 2 ulp(y * ) + η 2 , η 2 := ϵ|x|. Hence, if y 2 ̸ = RN(y * ), then Cx is within distance ≤ η 2 from the middle of two consecutive FP numbers. As η 2 ≤ η 1 , if Method 3 in

VII. INNER PRODUCTS AND SUMS

Another natural application of the FD2A is the computation of inner products with both fewer instructions and fewer rounding errors than when just using an FMA. Given two floating-point vectors x, y ∈ F n , recall from [START_REF] Lange | Error estimates for the summation of real numbers with application to floating-point summation[END_REF] that a sequence of n FMAs suffices to produce an approximation r FMA to the inner product x T y such that 3

| r FMA -x T y| ≤ nu|x| T |y|. ( 7
)
This bound is asymptotically optimal: for fixed n, there exist x, y ∈ F n such that both sides of ( 7) are asymptotically equivalent as u → 0. For example, with x T = [1+2u, 1, . . . , 1] and y T = [1 -u, u -u 2 , . . . , u -u 2 ], we have r FMA = 1 and

x T y = 1 + nu -(n + 1)u 2 .
With an FD2A this scheme is easily adapted by considering the following recursion: let r 1 := RN(x 1 y 1 + x 2 y 2 ) and, for i = 2, . . . , ⌈n/2⌉, let

r i := RN( r i-1 + x 2i-1 y 2i-1 + x 2i y 2i ),
with x n+1 = y n+1 = 0 if n is odd. As expected, the output r FD2A := r ⌈n/2⌉ has an error bound about twice smaller than the one in [START_REF] Markstein | Computation of elementary functions on the IBM RISC System/6000 processor[END_REF]. The next theorem makes this claim precise. Theorem 2. Given x, y ∈ F n , a sequence of ⌈n/2⌉ FD2A instructions suffices to produce an approximation r FD2A ∈ F to the inner product x T y such that

| r FD2A -x T y| ≤ mu|x| T |y|, m := ⌈n/2⌉.
Furthermore, this bound is asymptotically optimal.

Proof. Let the vectors z, ε, δ, s ∈ R m be defined by

z i = x 2i-1 y 2i-1 + x 2i y 2i , s i = r i , s 1 = z 1 + δ 1 = z 1 (1 + ε 1 ) and, for i ≥ 2, s i = z i + s i-1 + δ i = (z i + s i-1 )(1 + ε i ). For i ≥ 2, we thus have s i = RN(s i-1 + z i ) and |δ i | = |RN(s i-1 +z i )-(s i-1 +z i )| = min f ∈F |f -(s i-1 +z i )|; since s i-1 ∈ F, taking in particular f = s i-1 implies that |δ i | ≤ |z i |.
Hence the assumption needed to apply [25, Lemma 2.1] to ε, δ, s) is satisfied and since

z 1 + • • • + z m = x T y and |z 1 | + • • • + |z m | ≤ |x| T |y|, we deduce from this lemma that | r m -x T y| ≤ mu|x| T |y|.
The asymptotic optimality of this bound can be seen by considering x T = [1+2u, 0, 1, 0, 1, . . .], and y T = [1 -u, 0, u -u 2 , 0, u -u 2 , . . .], for which we have r m = 1 and

x T y = 1 + mu -(m + 1)u 2 .
Alternatively, one could approximate x T y by first computing the m values z i := RN(x 2i-1 y 2i-1 + x 2i y 2i ) and then, if m ≥ 2, adding them together by means of a ternary summation tree having ⌊m/2⌋ inner nodes. This would use m calls to FD2 followed by ⌊m/2⌋ calls to ADD3, and thus m + ⌊m/2⌋ ≈ 3n/4 operations. Although this operation count is larger than m by about n/4, the error bound associated with this second approach will be smaller than the one in Theorem 2 by at least about n/4, and thus at least about twice smaller. Specifically, classical analysis as in [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF]Ch. 4] would show that the constant mu obtained so far can be replaced by (h + 1)u + O(u 2 ), with h denoting the height of the underlying ternary tree. Since this tree has ⌊m/2⌋ inner nodes, its height h ranges from ⌈log 3 ⌊m/2⌋⌉ to ⌊m/2⌋. Hence

(h + 1)u ≲ m/2 • u ≈ n/4 • u.
Another way to reduce even further the error bound ( 7) is via compensation: besides the FMA-based recursion

r i := RN( r i-1 + x i y i ), i = 1, . . . , n,
which yields r FMA := r n , one computes correctly-rounded approximations to the errors e i := r i-1 + x i y i -r i by means of n calls to FD2A:

e i := RN( r i-1 + x i y i -r i ), i = 1, . . . , n.
(Of course, we set r 0 := 0 and so an FMA suffices to produce e 1 = x 1 y 1 -r 1 ∈ F.) These n values e i are then added together by means of ⌊n/2⌋ calls to ADD3 into a single value e. Finally, this value, which approximates the exact overall error e 1 + • • • + e n = x T y -r FMA , is used to produce the refined approximation r COMP := RN( r FMA + e).

Theorem 3. Under the assumptions of the previous theorem, a sequence of n FMA instructions and ⌊3n/2⌋ FD2A instructions suffices to produce r COMP ∈ F such that

| r COMP -x T y| ≤ u|x T y| + λ COMP |x| T |y| with λ COMP = 1 2 (n 2 + 3n)u 2 + 1 4 (n 3 + n 2 )u 3 . Proof.
| e -e| ≤ n 2 v|e 1 | + n 2 (v + v 2 ) + v n i=2 |e i |. (9)
Second, to bound |e i |, let us define r i = j≤i x i y i and r i = j≤i |x i y i |. Since the r i 's are produced recursively by the FMA, we deduce from [25, Lemma 2.1] that | r i -r i | ≤ iv r i . Consequently, for i ≥ 1 and with r 0 = 0,

|e i | ≤ v| r i-1 + x i y i | ≤ v r i + (i -1)v r i-1 . Using r i-1 ≤ r i ≤ r n = r for i ≤ n, we obtain |e i | ≤ (v + (i -1)v 2 )
r and, therefore,

|e 1 | ≤ v r and n i=2 |e i | ≤ nv 1 + n-1 2 v r. (10)
The claimed error bound then follows from the inequalities ( 8), ( 9), [START_REF] Lauter | An efficient software implementation of correctly rounded operations extending FMA: a + b + c and a × b + c × d[END_REF]. Furthermore, this was achieved using n + 1 FMA, n -1 FD2A, and ⌊n/2⌋ ADD3.

What is interesting with Theorem 3 is that λ COMP ∼ 1 2 n 2 u 2 , which is about half of the classical one from [START_REF] Ogita | Accurate sum and dot product[END_REF], and that this smaller bound is obtained using only about 3n/2 FD2A besides the initial FMA-based recursion.

Note also that λ COMP can be reduced even further by using a more balanced ternary tree when adding up the n approximate values e i 's (and thus without increasing the operation count): a ternary tree of height h leads to λ COMP ∼ (h + 1)nu 2 , so that with the minimal height h = ⌈log 3 ⌊n/2⌋⌉ we arrive at λ COMP ∼ 0.91n ln n u 2 .

Note finally that similar gains can of course be obtained for the special case of summation, where y T = [1, 1, . . . , 1], using exclusively the ADD3 instruction.

VIII. POLYNOMIAL EVALUATION

To be able to rapidly evaluate a low-degree polynomial on a specific architecture, taking into account the number of available operators, the depth of the arithmetic pipelines, etc., one can conduct an exhaustive or heuristic-based search among the various possible evaluation schemes [START_REF] Revy | Implementation of binary floating-point arithmetic on embedded integer processors : polynomial evaluation-based algorithms and certified code generation[END_REF], [START_REF] Cornea | Scientific Computing on Itanium-Based Systems[END_REF]. Having an FD2A instruction increases the number of possible schemes, potentially leading to a better optimum. We illustrate this by focusing on the most classical way (Horner's scheme with an FMA) and detailing how it can be adapted to exploit the FD2A when computing a(x) = n i=0 a i x i given x, a i ∈ F. Recall first that with an FMA, Horner's recursion is given by r 0 := a n and

r i := RN( r i-1 x + a n-i ), i = 1, . . . , n, (11)
and that this sequence of n FMA instructions yields an approximation r FMA := r n ∈ F such that

| r FMA -a(x)| ≤ γ n a(|x|), (12)
where γ n := nu/(1 -nu) if nu < 1 and a(x) := i |a i |x i ; see [29, §6.1]. The bound ( 12) is essentially best possible for fixed n and as u → 0, since taking in particular x = a n = 1 and

a n-i = u -u 2 for i ≤ n gives r FMA = 1 and a(x) = a(|x|) = 1 + nu -nu 2 .
With an FD2A, the recursion in ( 11) can be adapted immediately by computing first y := RN(x 2 ) and then

r i := RN( r i-1 y + a n-2i+1 x + a n-2i ), i = 1 : ⌈n/2⌉.
(Here, r 0 = a n and, if n is odd, a -1 = 0.) This scheme costs about n/2 FD2A instructions. In terms of accuracy, if x 2 is computed exactly, then the factor γ n ≈ nu in ( 12) can be decreased down to γ ⌈n/2⌉ ≈ 1 2 nu; unfortunately, in the general case where y ̸ = x 2 , the rounding error due to this initial squaring operation must be taken into account, bringing the factor γ n back and thus leading to the same a priori bound [START_REF] Hickmann | Experimental analysis of matrix multiplication functional units[END_REF] as for Horner with FMA. This is in contrast with what we have seen for inner products in section VII and Theorem 2.

There is, however, a more interesting way of exploiting the FD2A instruction in the context of Horner evaluation: the ability of the FD2A to recover the correctly-rounded value of the error of an FMA makes it possible to design a very simple compensated version of (11), as follows. Set r 0 := a n and e := 0; then compute for i from n

r i := RN( r i-1 x + a n-i ); // FMA e i := RN( r i-1 x + a n-i -r i );
// FD2A e := RN( ex + e i );

// FMA and finally return r COMP := RN( r n + e). Since e is initialized to zero, its final value is obtained using n -1 FMA and is equal to the approximation to n i=1 e i x n-i produced by Horner with FMA. The result r COMP is thus obtained by means of 2n -1 FMA, n FD2A, and 1 addition, that is, 3n instructions.

An accuracy bound is given by the following theorem. Interestingly, this bound matches the best one from [29, §6.3.2], where it needed 19n FP instructions: n calls to an EFT for the FMA, n + 1 additions, and n -1 FMA. (Remark that the EFT used in [START_REF] Louvet | Algorithmes compensés en arithmétique flottante : Précision, validation, performances[END_REF] is ErrFmac [START_REF] Boldo | Some functions computable with a Fused-Mac[END_REF] without the last Fast2Sum, so its cost is of 17 FP instructions instead of 20.) Theorem 4. Given a 0 , a 1 , . . . , a n , x ∈ F, the compensated Horner algorithm described above yields an approximation

r COMP ∈ F to a(x) = n i=0 a i x i such that | r COMP -a(x)| ≤ u|a(x)| + γ 2 n a(|x|).
Proof. Let r = a(x), r = a(|x|), and v = u/(1 + u). For 1 ≤ i ≤ n, define also r i = i j=0 a n-j x i-j and let e i be such that r i-1 x + a n-i = r i + e i . Multiplying both sides of this equality by x n-i and then summing from i = 1 to n yields after simplification r = r n +e, where e = n i=1 e i x n-i . Since

r COMP = ( r n + e)(1 + δ) with |δ| ≤ v, it follows that | r COMP -r| ≤ v|r| + (1 + v)| e -e|. (13)
Let us now bound | e -e|. Since e is obtained by evaluating | n i=1 e i x n-i | by Horner with FMA, it follows from [START_REF] Hickmann | Experimental analysis of matrix multiplication functional units[END_REF] and

| e i -e i | ≤ v|e i | that | e -e| ≤ ((1 + v)γ n-1 + v) e, (14)
where e = n i=1 |e i ||x| n-i . Now, |e i | ≤ u| r i | with r i obtained by evaluating r i by Horner with FMA. Hence [START_REF] Hickmann | Experimental analysis of matrix multiplication functional units[END_REF] implies

| r i | ≤ |r i | + γ i i j=0 |a n-j ||x| i-j , so that |e i ||x| n-i ≤ u(1 + γ n ) r.
Therefore, e ≤ nu(1 + γ n ) r and the conclusion follows from combining this bound with the ones in ( 13) and [START_REF]IEEE, IEEE Standard for Floating-Point Arithmetic (IEEE Std[END_REF].

IX. DOUBLE-WORD ARITHMETIC

Double-word arithmetic (often referred to as double-double in the literature) deals with real numbers x represented as the unevaluated sums of two FP numbers x h and x ℓ such that x h = RN(x h + x ℓ ) = RN(x). Several algorithms for manipulating double-word (DW) numbers have been introduced since the pioneering work of Dekker [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF], in particular in [START_REF] Li | Design, implementation and testing of extended and mixed precision BLAS[END_REF], [START_REF] Hida | C++/Fortran-90 double-double and quad-double package, release 2.3.17[END_REF], [START_REF] Joldes | Tight and rigorous error bounds for basic building blocks of double-word arithmetic[END_REF], [START_REF] Popescu | Towards fast and certified multiple-precision librairies[END_REF], [START_REF] Lefèvre | Accurate calculation of Euclidean norms using double-word arithmetic[END_REF].

The availability of an FD2A instruction makes it possible to design new algorithms for DW arithmetic, that are cheaper in terms of numbers of instructions, often more accurate, and simpler to state and analyze than the traditional algorithms mentioned above. A key building-block will be the following implementation of the EFT for addition, which appears in [START_REF] Ogita | Accurate sum and dot product[END_REF]. This algorithm performs only two operations, instead of six for 2Sum. It does not require a condition such as |a| ≥ |b| as Fast2Sum does. These features will, of course, benefit to all the DW algorithms described below.

A. Sum of a DW number and an FP number

: (s, f ) ← TwoSum Add3(a h , b) 2: e ← RN(f + a ℓ ) 3: (c h , c ℓ ) ← TwoSum Add3(s, e)
This algorithm is deduced from Algorithm DWPlusFP in [START_REF] Joldes | Tight and rigorous error bounds for basic building blocks of double-word arithmetic[END_REF], using TwoSum Add3 instead of 2Sum and Fast2 Sum. It returns the same output in 5 FP instructions instead of 10, and one can proceed as in [32, §2] to check that an asymptotically-optimal bound on its relative error is 2u. The analysis is simpler as it avoids having to check the usual FastTwoSum condition |s| ≥ |e| used in [START_REF] Joldes | Tight and rigorous error bounds for basic building blocks of double-word arithmetic[END_REF].

B. Sum of two DW numbers

The authors of [START_REF] Joldes | Tight and rigorous error bounds for basic building blocks of double-word arithmetic[END_REF] analyze algorithm AccurateDWPlusDW of Hida, Li and Bailey's library [START_REF] Hida | C++/Fortran-90 double-double and quad-double package, release 2.3.17[END_REF]. They show that it has an asymptotically-optimal relative error bound of 3u 2 /(1 -4u). If an FD2A operation is available, one can use that algorithm with each 2Sum and Fast2Sum replaced by a TwoSum Add3, which saves 10 operations without changing the error. However, the availability of an FD2A operation allows one to design an algorithm that is both cheaper and more accurate:

Algorithm 3 : AccDWPlusDW FD2A(a h , a ℓ , b h , b ℓ ). Com- putes a DW approximation z to a + b. 1: (s h , s ℓ ) ← TwoSum Add3(a h , b h ) 2: (t h , t ℓ ) ← TwoSum Add3(a ℓ , b ℓ ) 3: (c h , c ℓ ) ← TwoSum Add3(s ℓ , t h ) 4: (v h , v ℓ ) ← TwoSum Add3(s h , c h ) 5: w ← RN(t ℓ + v ℓ + c ℓ ) 6: (z h , z ℓ ) ← TwoSum Add3(v h , w)
This algorithm uses 5 FP additions and 6 FD2A operations (instead of 20 additions for AccurateDWPlusDW), and its relative error is bounded by 2u 2 (the proof is similar to the proof of AccurateDWPlusDW in [START_REF] Joldes | Tight and rigorous error bounds for basic building blocks of double-word arithmetic[END_REF]).

C. Product of a DW number by an FP number

Without an FD2A, multiplying a DW number by an FP number is a rather complex operation (for example, the most accurate of the three algorithms analyzed in [START_REF] Joldes | Tight and rigorous error bounds for basic building blocks of double-word arithmetic[END_REF] takes 10 FP operations and has a tight error bound of 3 2 u 2 + 4u 3 ). With an FD2A, this becomes much simpler and more accurate: Algorithm 4 : DWTimesFP FD2A(a h , a ℓ , b). Computes a DW approximation z to ab.

1: s ← RN(a h • b + a ℓ • b) 2: e ← RN(a h • b + a ℓ • b -s) 3: (z h , z ℓ ) ← TwoSum Add3(s, e)
It follows from Property 1 that the relative error of this algorithm is at most 1 2 u 2 . This bound is asymptotically optimal: with the values a h = 1 + 2u, a ℓ = 3 2 u 2 , and b = 1 -u, we obtain s = 1 and e = u -u 2 , which leads to a relative error equal to

1 2 u 2 -3 2 u 3 1+u-1 2 u 2 -3 2 u 3 ∼ 1 2 u 2 as u → 0.
As e = RN(ab -s), in most applications the last line of the algorithm can be omitted. That line is however necessary if we wish to always output a DW, as we may have s ̸ = RN(s + e) (an example is with a

h = 1 + 2u, a ℓ = 5 2 u 2 , b = 1 -u).

D. Product of two DW numbers

Three algorithms for computing the product of two DW numbers are analyzed in [START_REF] Muller | Formalization of double-word arithmetic, and comments on "Tight and rigorous error bounds for basic building blocks of double-word arithmetic[END_REF]. The most accurate of them has an error bound of 4u 2 and requires 9 FP operations. The availability of an FD2A operation would improve the situation, both in terms of accuracy and number of operations. Consider the following algorithm, which uses 5 FP operations:

Algorithm 5 : DWTimesDW FD2A(a h , a ℓ , b h , b ℓ ). Com- putes a DW approximation z to ab. 1: (c h , e) ← TwoProd(a h , b h ) 2: c ℓ ← RN(a ℓ • b h + a h • b ℓ + e) 3: (z h , z ℓ ) ← TwoSum Add3(c h , c ℓ )
This algorithm has relative error ≤ 3u 2 . This can be shown as follows. The absolute error is |ϵ 1 +ϵ 2 | with ϵ 1 = RN(a ℓ b h + a h b ℓ + e) -(a ℓ b h + a h b ℓ + e) the rounding error occurring when computing c ℓ , and ϵ 2 the error due to a ℓ b ℓ being ignored. Without loss of generality we assume If u ≤ 1 16 (which always holds in practice), B 1 ≥ B 2 , giving the final bound 3u 2 . Furthermore, this bound is asymptotically optimal. To see this, let a h = b h = 1 + 2u and a ℓ = b ℓ = u -u 2 . One can check that c h = 1 + 4u, e = 4u 2 , and c ℓ = 2u + 4u 2 , so that the relative error of c h + c ℓ (and thus of z h + z ℓ as well) is 3u 2 -6u 3 +u 4 1+6u+7u 2 -6u 3 +u 4 ∼ 3u 2 as u → 0. A more accurate calculation is possible at a higher cost (8 FP operations), using the following algorithm. The relative error of this algorithm is bounded by

1 ≤ a h , b h ≤ 2 -2u. This implies that |a ℓ |, |b ℓ | ≤ u, so that |ϵ 2 | ≤ u 2 . • If a h b h ≤ 2,
4u 2 +3u 3 2(1-u) 2 ∼ 2u 2 + 11 2 u 3 .
The proof is similar to the one of the previous algorithm.

E. Division of a DW number by a DW number

This algorithm has the same general structure as DW-DivDW3 in [START_REF] Joldes | Tight and rigorous error bounds for basic building blocks of double-word arithmetic[END_REF], whose relative error bound is 9.8u 2 . Its relative error is ≤ 7.8u 2 when p ≥ 13. The proof follows the steps of the proof of [32, Theorem 7.2] (using the error bounds on DWTimesFP FD2A, DWPlusFP FD2A, and DW-TimesDW FD2A given in this paper). The authors of [START_REF] Lefèvre | Accurate calculation of Euclidean norms using double-word arithmetic[END_REF] analyze and formally prove a DW square root algorithm (that already appears in [START_REF] Rump | Faithfully rounded floating-point computation[END_REF], without the final renormalization). The availability of an FD2A operator makes it possible to fuse two steps of this algorithm, and to replace its final renormalization by a call to TwoSum Add3. This gives Algorithm 8 : DWSQRT FD2A(a h , a ℓ ). Computes a DW approximation z to √ a.

1: if a h = 0 then (z h , z ℓ ) ← (0, 0) 2: else 3: s h ← RN( √ a h )

4: ρ ← RN(a h -s 2 h + a ℓ ) 5: s ℓ ← RN(ρ/(2 • s h )) 6: (z h , z ℓ ) ← TwoSum Add3(s h , s ℓ )
The relative error bound is the same as for the original algorithm, namely, (25/8)u 2 = 3.125u 2 . The proof is the same as the one given in [START_REF] Lefèvre | Accurate calculation of Euclidean norms using double-word arithmetic[END_REF]Theorem 3.6].

X. CONCLUDING REMARKS

We have shown that many basic building blocks of numerical computing (complex arithmetic, error-free transforms, inner products, polynomial evaluation, double-word arithmetic, etc.) could greatly benefit from the availability of FD2 and FD2A operations. This would allow for more accurate calculations, via algorithms that use fewer FP instructions and are simpler to state and analyze. An additional benefit would be a somehow simpler management of spurious overflow issues.

  For readability define here r = x T y, r = |x| T |y|, r = r FMA , e = n i=1 e i , and v = u/(1+u). Since r COMP = RN( r+ e) and r + e = r, we have | r COMP -r| = |( r + e)(1 + δ) -r| = |rδ + ( e -e)(1 + δ)|, |δ| ≤ v. Hence | r COMP -r| ≤ v|r| + (1 + v)| e -e| (8) and it remains to bound | e-e|. First, | e-e| ≤ | e-n i=1 e i |+ n i=2 | e i -e i | since e 1 = e 1 . On the one hand, | e i -e i | ≤ v|e i | for all i. On the other hand, adding the e i 's recursively by means of ⌊n/2⌋ ADD3, we deduce from [25, Lemma 2.1] that | e -n i=1 e i | ≤ ⌊n/2⌋v n i=1 | e i |. Hence

Algorithm 1 :

 1 TwoSum Add3(a, b).

1 :

 1 s ← RN(a + b) 2: e ← RN(a + b -s) // a + b = e + s (EFT)

Algorithm 2 :

 2 DWPlusFP FD2A(a h , a ℓ , b). Computes a DW approximation c to a + b.

1

 1

Algorithm 6 :

 6 DWTimesDW2 FD2A(a h , a ℓ , b h , b ℓ ). Computes a DW approximation z to ab.

1 :

 1 (c, e) ← TwoProd(a h , b h ) 2: v h ← RN(a ℓ • b h + a h • b ℓ ) 3: v ℓ ← RN(a ℓ • b h + a h • b ℓ -v h ) 4: w ← RN(v ℓ + a ℓ b ℓ ) 5: s ← RN(e + v h + w) 6: (z h , z ℓ ) ← TwoSum Add3(c, s)

Algorithm 7 :

 7 DWDivDW FD2A(a h , a ℓ , b h , b ℓ ). Computes a DW approximation z to a/b.

  then |e| ≤ u and, from[START_REF] Muller | Formalization of double-word arithmetic, and comments on "Tight and rigorous error bounds for basic building blocks of double-word arithmetic[END_REF] Lemma 2.11],a h +b h ≤ 3-2u. Thus |a ℓ b h + a h b ℓ + e| ≤ u(a h +b h )+ u < 4u, so that |ϵ 1 | ≤ 2u 2 . Therefore, |(c h + c ℓ ) -ab| ≤ 3u 2 . As ab ≥ 1, the relative error is ≤ B 1 = 3u 2 . • If a h b h > 2, then |e| ≤ 2u, so that |a ℓ b h + a h b ℓ + e| ≤ 2u(2 -2u) + 2u < 6u and |ϵ 1 | ≤ 4u 2 . Hence |ϵ 1 + ϵ 2 | ≤ 5u 2 . Since ab ≥ a h (1 -u) • b h (1 -u), the relative error is ≤ B 2 = 5u 2 /(2(1 -u) 2 ).

  1: t h ← RN(1/b h ) 2: r h ← RN(1 -t h • b h ) 3: r l ← -RN(t h • b ℓ ) 4: (e h , e ℓ ) ← TwoSum Add3(r h , r ℓ ) 5: (δ h , δ ℓ ) ← DWTimesFP FD2A(e h , e ℓ , t h ) 6: (m h , m ℓ ) ← DWPlusFP FD2A(δ h , δ ℓ , t h ) 7: (z h , z ℓ ) ← DWTimesDW FD2A(a h , a ℓ , m h , m ℓ )F. Square root of a DW number

See https://patents.google.com/patent/US20220222073A1/en

We recall that we have assumed an unbounded exponent range. If this is not the case, then there is a condition on the exponents.

Note that instead of nu in[START_REF] Markstein | Computation of elementary functions on the IBM RISC System/6000 processor[END_REF], a term of the form (1 + u) n -1 = nu + O(u 2 ) can be obtained classically, following[START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF] Ch. 3]; however, the fact that the term O(u 2 ) can be removed is a consequence of[START_REF] Lange | Error estimates for the summation of real numbers with application to floating-point summation[END_REF] Lemma 2.1].