Stef Graillat

Ping Tak

Peter Tang

Christoph Lauter

Shin ' Ichi

A Efficient calculations of faithfully rounded l 2 -norms of n-vectors

Keywords: Design, Algorithms, Performance Floating-point arithmetic, error-free transformations, faithful rounding, 2-norm, underflow, overflow

In this paper, we present an efficient algorithm to compute the faithful rounding of the l 2 -norm of a floatingpoint vector. This means that the result is accurate to within one bit of the underlying floating-point type. This algorithm does not generate overflows or underflows spuriously, but does so when the final result indeed calls for such a numerical exception to be raised. Moreover, the algorithm is well suited for parallel implementation and vectorization. The implementation runs up to 3 times faster than the netlib version on current processors.

INTRODUCTION

Computing the l 2 -norm x 2 = n j=1 x 2 j of a vector x = [x 1 , x 2 , . . . , x n ] T is prevalent in scientific and engineering applications. This operation is part of the first (lowest) level of the Basic Linear Algebra Subroutine BLAS1. The simplicity of the formula n j=1 x 2 j is misleading. Summing the squares can cause unwarranted (spurious) overflows or underflows in many instances when in fact x 2 is well within the normal range of the working-precision floating-point arithmetic. Common implementations such as the public version of LAPACK [START_REF] Anderson | LAPACK Users' guide[END_REF] released by netlib 1 essentially compute the l 2 -norm as x × x/ x 2 where x is max j |x j |. That implementation requires n divisions in total, which is significantly more expensive than the naïve formula would suggest. Spurious exceptions aside, accuracy is also an issue with implementations that rely on accumulation of squares in working-precision arithmetic. In the worst-case scenario, the last log 2 (n) bits of the binary floating-point result could be corrupted. Equivalently the last log 10 (n) digits of the result, when displayed in decimal, could be corrupted. Improving the accuracy of l 2 -norm computation enhances the qualities of the larger computational tasks relying on it. Moreover, highly accurate l 2 -norm improves the chances of obtaining reproducible numerical results should the l 2 -norm computation be done in parallel, with threads or vector-floating-point (SIMD) instructions2 .

In this paper, we present a new division-free l 2 -norm algorithm that is amenable to straightforward parallel implementations. We prove that the algorithm always returns a faithfully rounded result, to be defined rigorously later. For now (and informally), this means that the result is accurate to within one bit of the underlying floating-point type. The algorithm also reports overflow and underflow faithfully, a property to be defined later. Loosely speaking, these exceptions are triggered only when the true value x 2 calls for the event. On current processors (with AVX extensions), our implementation runs at least as fast as the netlib version, and up to three times faster when the IEEE754 fused-multiply-add instruction is available.

There are two main features of our algorithm. First, it accumulates the squares, x 2 j , using a pair of floating-point numbers, providing essentially double the underlying floating-point precision. Our technique is similar to the addition operator in the double-double library [START_REF] Li | Design, implementation and testing of extended and mixed precision BLAS[END_REF] but at almost twice the speed by exploiting the nonnegative nature of sum-of-squares. We provide a rigorous analysis of the accuracy properties of our accumulation process. Second, we eliminate all spurious exceptions by scaling the input data but without using division. The technique is a "binning" method and is similar to the one proposed in [START_REF] Blue | A portable Fortran program to find the Euclidean norm of a vector[END_REF]]: the vector elements x j are grouped, i.e. binned, into small, medium and large inputs, such that, after appropriate scaling up or down, their squares in each bin can be computed without spurious overflow or underflow. However, we improve the binning technique in that [START_REF] Blue | A portable Fortran program to find the Euclidean norm of a vector[END_REF] requires three bins and we use only two bins here, resulting in economy of registers usage and performance improvement. The accumulation and binning are amenable to straightforward parallel implementations. Our reference implementation uses data parallelism through SIMD instructions, but adding thread parallelism is straightforward. Our technique does not incur any memory overhead; in particular each vector element x j is read only once, the same way it would be in other approaches.

We organize the rest of the paper as follows. Section 2 defines the key technical terms to make our paper self contained. We formulate our main problem and state our objectives. We present and establish the main theorem in Section 3. The essence is that if x 2 2 is computed to enough accuracy as a floating-point pair as "leading"-plus-"correction," then a standard IEEE-conforming square root on the "leading" part yields a faithfully-rounded l 2 -norm. Section 4 presents a serial and a parallel l 2 -norm algorithm without binning. In the absence of exceptions, we prove the numerical properties of these two algorithms that will lay the foundation for the actual binned algorithms that provide the spurious-exception-free property. Section 5 presents the binned versions and the proofs of the faithful rounding and spurious-exception-free nature. Section 6 shows numerical and performance test results to corroborate with the theoretical analysis presented.

Efficient calculations of faithfully rounded l2-norms of n-vectors A:3

BACKGROUND

Throughout this paper we consider a specific IEEE 754 binary floating-point type. Let F denote the entire set of finite values in this type, identified by three parameters ε, e min , and e max :

F = {±2 e+1 m ε | m ∈ N, e min ≤ e ≤ e max , 0 ≤ mε < 1}.
For example, (ε, e min , e max ) is (2 -24 , -126, 127) for the binary32 format, and (2 -53 , -1022, 1023) for the binary64 format. Denote the smallest and largest positive normalized numbers by F small = 2 emin and F large = 2 emax+1 (1 -ε). We assume ε ≤ 2 -24 throughout this paper.

Closely related to F is the set F where no upper limit of the exponent is imposed:

F = {±2 e+1 m ε | m ∈ N, e min ≤ e, 0 ≤ mε < 1}. It is clear that F ∩ [0, 2 emax+1 ) = F ∩ [0, 2 emax+1 ).
In particular, numbers x ∈ F where 0 < |x| < 2 emin are also denormalized.

For any real number α ∈ R, •(α) is the IEEE round-to-nearest-even function that maps any finite real number to F , • : R → F . The rounding •(α) of a real number α is the number in F that is closest to α, with a tie broken by choosing •(α) to have an even mantissa (least-significant bit being zero). A crucial property of •(α) is best stated in terms of the ulp (units of last place) function defined as follows. For all α ∈ R,

ulp(α) = 2 e+1 ε if |α| ∈ [2 e , 2 e+1
), e ≥ e min , 2 emin+1 ε otherwise.

Note that this definition of ulp reflects the floating-point arithmetic properties in the denormalized range. For any real number α and integer k such that both |α| and |2 k α| ≥ F small , then ulp(2 k α) = 2 k ulp(α), and ulp(α) ≤ 2|α|ε where equality ulp(α) = 2|α|ε holds if and only if log 2 (|α|) is an integer. It is easy to see that •(α) is at worst half an unit of last place away from α: For a real number α, |α| ∈ [2 e , 2 e+1 ], then |•(α) -α| ≤ ulp(α)/2. Note that this holds even for e < e min .

The clipping function clip maps floating-point numbers in

F to F ∪ {-∞, +∞}: For all x ∈ F , clip(x) = x for |x| < 2 emax+1 , x |x| ∞ for |x| ≥ 2 emax+1 .
The four basic IEEE arithmetic operations ⊕, , ⊗, can be defined in terms of the rounding and clipping functions:

a b = clip(•(a * b))
for any * ∈ {+, -, ×, /} and a, b ∈ F (with b = 0 if the operation " * " is division). Similarly, the IEEE square root function sqrt is defined as sqrt

(x) = •( √ x) for all x ∈ F ∩ [0, ∞). The clipping function is not needed here as √ x ≤ F large .
For any α ∈ R, we define α's faithful set of floating-point numbers ♦(α) as follows. ♦(α) is the singleton {α} if α ∈ F . Otherwise, ♦(α) is the set of two numbers in F that are closest to α from below and above. In other words,

♦(α) = max y {y ∈ F | y ≤ α}, min y {y ∈ F | y ≥ α} .
Clipping the set ♦(α) is by definition the set consisting of the clipping of each of the elements in ♦(α):

clip(♦(α)) = {clip(a)|a ∈ ♦(α)}.
This paper presents a parallelizable and division-free algorithm AccuNrm2 such that for all vectors x of practical length, AccuNrm2(x) ∈ clip(♦( x 2 )) where x = [x 1 , x 2 , . . . , x n ] T , x j ∈ F. We call this numerical property a faithful rounding of x 2 . Note that our definition of faithful rounding is slightly stronger than that in [Muller et al. 2010], differing with it only in the case when x 2 ≥ 2 emax+1 . In this case, our definition requires that +∞ be returned, while the definition in [Muller et al. 2010] allows either +∞ or F large to be a faithful rounding.

In addition to computing a faithfully rounded numerical value, AccuNrm2(x) also reports overflow and underflow faithfully in the following sense. Given an implementation G(x) of a function g(x), we say G(x) reports

(1) overflow faithfully if:

-G(x) never reports an overflow when |g(x)| ≤ F large .

-G(x) always reports an overflow when |g(x)| ≥ 2 emax+1 . (2) underflow faithfully:

-

When |g(x)| ≥ F small , G(x) never reports underflow. -When 0 < |g(x)| ≤ F small -2 emin+1 ε, then -|G(x)| ≤ F small -2 emin+1 ε always, -if underflow is unmasked, G(x) reports underflow.
-if underflow is masked and g(x) ∈ F, then G(x) = g(x) and does not report underflow.

We emphasize that achieving faithful rounding is non-trivial. One can show that

•( •(σ)) ∈ ♦( x 2 )
, where σ = j x 2 j = x T x is the exact sum of squares (or inner product). This says the correctly rounded square root of the correctly rounded sum of squares is a faithfully rounded l 2 -norm. Nevertheless, computing the correctly rounded sums of squares is very expensive [START_REF] Ogita | Accurate Sum And Dot Product[END_REF]Rump et al. 2008a;2008b]. On the other hand, examples exist where •( √ S) / ∈ ♦( x 2 ) for some S ∈ ♦(σ). That is, computing the σ to only slightly worse than the correctly rounded sum of squares can cause unfaithful rounding. Our algorithm in essence computes a floating-point number S ≈ σ very accurately and yet efficiently. Theorem 3.3 gives a condition on the accuracy of S that guarantees •( √ S) to be a faithful rounding of x 2 . The fact •( •(σ)) ∈ ♦( x 2 ) alluded to earlier follows trivially from Theorem 3.3 as well.

The fundamental task is that of computing σ = j x 2 j accurately. As we will see later, it is possible to transform this computation into a sum without loss of information (no rounding error). The problem is now transformed into the accurate computation of a sum. There is an abundant literature about floating-point summation (see [Higham 2002, chap. 4], [START_REF] Knuth | The Art of Computer Programming[END_REF][START_REF] Ogita | Accurate Sum And Dot Product[END_REF]Rump et al. 2008a;2008b;[START_REF] Rump | Ultimately fast accurate summation[END_REF][START_REF] Zhu | Correct Rounding and a Hybrid Approach to Exact Floating-Point Summation[END_REF]2010] and references therein). For our purpose, we only need an accurate summation algorithm whose precision is doubled because the entries are nonnegative numbers. For such a precision, a straightforward adaptation of the algorithm Sum2 [START_REF] Ogita | Accurate Sum And Dot Product[END_REF] is very efficient since it requires 8(n -1) floating-point operations (flops) where n is the size of the vector. The resulting sum has a relative error no more than on the order of n 2 ε 2 . Another choice is to use the double-double arithmetic presented in [START_REF] Li | Design, implementation and testing of extended and mixed precision BLAS[END_REF]. The resulting sum is much more accurate, having a relative error no more than 2nε 2 /(1 -2nε 2 ). The cost, however, is 20(n -1) flops. The difference between the two algorithms comes from the "renormalization steps" that are present in the double-double library. As will be shown later, an error bound in the order of nε 2 as opposed to n 2 ε 2 is crucial if faithful rounding is to be guaranteed for a general vector length n. We therefore devise an algorithm that is faster than, but of comparable accuracy to, the addition operator in the double-double library. Instead of performing two renormalization steps in the addition of two double-double numbers, we perform only one renormalization step in a careful manner. The resulting error bound is 3nε 2 /(1 -3nε 3 ), at a cost of 11(n -1) flops, which is almost twice as fast as 20(n -1) flops.

A naïve computation of the l 2 -norm can cause spurious overflow and underflow. The netlib library addresses the overflow issue but not that of underflow. Spurious underflow can cause significant performance degradation. Blue's work [START_REF] Blue | A portable Fortran program to find the Euclidean norm of a vector[END_REF]] uses three bins to eliminate spurious overflows and underflows. Accuracies of the netlib and Blue algorithms are comparable, where close to log 10 (n) digits can be corrupted in the worst case. In summary, our algorithm is accurate to within one binary bit, free from spurious over/underflows, and run faster than the netlib version.

We state without proofs the following elementary facts about floating-point arithmetic. These facts and notations will be used freely in the subsequent sections.

-Let a, b ∈ F . The absolute error bound

|•(a op b) -(a op b)| ≤ ulp(a op b)/2
holds for all op ∈ {+, -, ×, /} (excluding division by zero). The relative error bound

•(a op b) = (a op b)(1 + δ), |δ| ≤ ε holds for op ∈ {+, -}. For op ∈ {×, /}, this relative error bound holds if F small ≤ |a op b| (excluding division by zero). -Rounding to nearest is monotonic: Given real numbers α, β ∈ R, α ≤ β implies •(α) ≤ •(β). -For a, b ∈ F, if S = •(a + b) ∈ F, then a + b -S ∈ F.
In particular the value s = a + b -S satisfies the relationship S + s = a + b and •(S + s) = S. Furthermore, S and s can be computed from a and b by a sequence of instructions involving ⊕ and (see for example [START_REF] Knuth | The Art of Computer Programming[END_REF]] Thm B, page 236, and [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF]). We encapsulate these facts by the two functions TwoSum(a, b) and FastTwoSum(a, b). They deliver S and s where S + s = a + b and •(S + s) = S. The former handles general a and b and requires 6 floating-point operations; the latter requires only 3 floating-point operations, but relies on the assumption that |a| ≥ |b|. Mapping a, b to S, s is usually called an error-free transformation as a + b = S + s exactly. Our algorithm and its implementation use both functions.

-Similar to error-free transformations for sum, we have error-free transformations for product. For a, b ∈ F, if |ab| ≥ F small /ε, and

P = •(a × b) ∈ F, then a × b -P ∈ F.
In particular the value p = a × b -P satisfies the relationship P + p = a × b and •(P + p) = P . P and p can be computed from a and b with the IEEE754-2008 fused multiply-add (FMA) instruction: P := a⊗b and p := •(a×b-P ). Alternatively, one can use a sequence of ⊕, and ⊗ instructions as outlined in [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF]]. We denote the exact product function that returns P and p by TwoProd(a, b). Our implementation uses either the FMA-based sequence or the sequence given in [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF]], depending on whether the FMA instruction is supported on the available hardware.

MAIN THEOREM

The goal is to compute

x 2 = √ σ, σ = n j=1
x 2 j , faithfully. Our algorithm is built on the core case when σ is within the "normal" range of [F small , F large ]. We first compute an ac-curate floating-point approximation, S, to σ. The IEEE square root (correctly rounded) •( √ S) is returned as the final result. This section shows that if S is accurate to within a specific threshold, the final result is faithful:

•( √ S) ∈ ♦( x 2 ) = clip(♦( x 2 )).
LEMMA 3.1. Let α, α be two real numbers in the interval [2 e , 2 e+1 ] for some integer e. If |α -α| < ulp(2 e )/2, then •(α ) ∈ ♦(α).

PROOF. Let a = •(α ). We first note that a and α are close to one another:

|a -α| ≤ |a -α | + |α -α| ≤ ulp(2 e )/2 + |α -α|, because a = •(α ) < ulp(2 e )/2 + ulp(2 e )/2, by assumption, |a -α| < ulp(2 e ).
(1)

To prove a ∈ ♦(α), we analyze all of the three possibilities of a = α, a < α and a > α.

The case of a = α is trivial because a ∈ F and therefore a ∈ {a} = ♦(a) = ♦(α).

Consider the case of a < α. This means that a < 2 e+1 . Hence a ∈ [2 e , 2 e+1 ) if e ≥ e min , and a ∈ [0, 2 emin ) if e < e min . Regardless, the next floating-point number in F that is bigger than a is a + ulp(2 e ). We have

min{y ∈ F |y > a} = a + ulp(2 e ) > a + |a -α|, by (1) = a + (α -a), because a < α = α.
Therefore, while a < α, the next floating-point number above a is strictly bigger than α. In other words, a = max{y ∈ F |y ≤ α}. This says that a ∈ ♦(α) by definition.

Consider now the final case: a > α. Since a > 2 e , we either have a ∈ (2 e , 2 e+1 ] if e ≥ e min , or a ∈ (0, 2 emin ] if e < e min . In either case, the next number in F that is smaller than a is a -ulp(2 e ).

max{y

∈ F |y < a} = a -ulp(2 e ) < a -|a -α|, by (1) = a -(a -α), because a > α = α.
Therefore, while a > α, the next floating-point number below a is strictly less than α.

In other words, a = min{y ∈ F |y ≥ α}. This says once again that a ∈ ♦(α). The proof is now complete.

LEMMA 3.2. Let σ ∈ [F small , F large ] be a real number in the interval [2 e , 2 e+1 ). In particular e min ≤ e ≤ e max . Let S, s ∈ F be such that

•(S + s) = S. If |(S + s) -σ| < σε/2, then S ∈ [2 e , 2 e+1 ] ∩ [0, F large ] and |s| ≤ 2 e ε.
PROOF. We will first establish the fact that S

∈ [2 e , 2 e+1 ] ∩ [0, F large ]. By assumption σ(1 -ε/2) < S + s < σ(1 + ε/2). Since σ ∈ [2 e , 2 e+1
), e ≥ e min , we have

2 e -ulp(2 e )/4 = 2 e (1 -ε/2) < S + s < 2 e+1 (1 + ε/2) = 2 e+1 + ulp(2 e+1 )/4. Consequently •(S + s) ∈ [2 e , 2 e+1
]. Moreover, S + s < σ + σε/2 implies also that S + s < F large + ulp(2 emax )/2 which leads to •(S + s) ≤ F large . But S = •(S + s) by assumption, and thus we have established that S ∈ [2 e , 2 e+1 ] ∩ [0, F large ].

Turning now to s, there are only two possibilities:

S + s ∈ [2 e , 2 e+1 ] or S + s / ∈ [2 e , 2 e+1 ]. If S + s ∈ [2 e , 2 e+1 ], then |s| = |•(S + s) -(S + s)| ≤ ulp(2 e )/2 = 2 e ε. Consider now S + s / ∈ [2 e , 2 e+1
]. But since S = •(S + s) and that we have established previously that •(S + s) ∈ [2 e , 2 e+1 ], we are left with only two cases:

(1) S = 2 e and s < 0: In this case, the inequality 2 e -ulp(2 e )/4 < S + s implies s > -ulp(2 e )/4 = -2 e ε/2. Thus |s| < 2 e ε/2. (2) S = 2 e+1 and s > 0: In this case, the inequality S + s < 2 e+1 + ulp(2 e+1 )/4 implies that s < ulp(2 e+1 )/4 = 2 e ε. Thus |s| < 2 e ε.

To review, |s| ≤

2 e ε if S + s ∈ [2 e , 2 e+1 ] and |s| < 2 e ε if S + s / ∈ [2 e , 2 e+1
]. Thus |s| ≤ 2 e ε always and the proof is complete.

THEOREM 3.3. Let σ ∈ [F small , F large ] be a real number and S, s ∈ F where •(S +s) = S. If |(S + s) -σ| < εσ/8, then •( √ S) ∈ ♦( √ σ).
PROOF. We establish this theorem by showing √ σ and √ S satisfy the conditions for α and α in Lemma 3.1. The assumption

|(S + s) -σ| < εσ/8 implies S + s > (1 -ε/8)σ. •(S + s) = S and that S + s is obviously positive imply S ≥ (S + s)(1 -ε). Thus S > (1 -ε)(1 -ε/8)σ, > (1 -3ε/2)σ, because ε ≤ 2 -24 √ S > (1 -3ε/2) √ σ, because 1 -3ε/2 > 1 -3ε/2, √ S + √ σ > 2(1 -3ε/4) √ σ. (2)
We derive an upper bound of | √ S -√ σ| as follows.

| √ S - √ σ| = |S -σ| √ S + √ σ , < |S -σ| 2 √ σ(1 -3ε/4) , by (2) < |S -σ| 2 √ σ (1 + ε), because (1 -3ε/4) -1 < (1 + ε) < εσ/8 + |s| 2 √ σ (1 + ε), by assumption, | √ S - √ σ| < ε 16 √ σ + |s| 2 √ σ (1 + ε). (3)
There is a unique interval of the form (2 e ε). We accomplish this by using Equation (3) and the simple case analysis tabulated below:

[2 2e , 2 2e+2 ) that contains σ. If σ is in the "left" half: σ ∈ [2 2e , 2 2e+1 ), then 2e ≥ e min and Lemma 3.2 shows that S ∈ [2 2e , 2 2e+1 ] and |s| ≤ 2 2e ε. If σ is in the "right" half: σ ∈ [2 2e+1 , 2 2e+2 ),
Case of ε √ σ 2 e ε |s| √ σ2 e+1 ε √ S- √ σ 2 e ε < ε √ σ 2 e+4 ε + |s| √ σ2 e+1 ε (1 + ε) σ ∈ [2 2e , 2 2e+1 ) < √ 2 ≤ 2 2e ε 2 2e+1 ε = 1 2 < √ 2 16 + 1 2 (1 + ε) < 1 σ ∈ [2 2e+1 , 2 2e+2 ) < 2 ≤ 2 2e+1 ε 2 2e+1 √ 2ε = 1 √ 2 < 1 8 + 1 √ 2 (1 + ε) < 1 Clearly, | √ S - √ σ| < 2 e ε always and •( √ S) ∈ ♦( √ σ) as claimed. 4. FAITHFUL 2-NORM -CORE CASE Let x = [x 1 , x 2 , . . . ,
x n ] T be the vector in question. The core algorithm considers the case of normal range in the sense that σ = x T x is safely away from overflow and that the least significant bit of each of the x 2 j does not underflow. More formally, we consider vectors x in safe range, defined as (1) σ ≤ F large /2, and (2) for all j, x j is either 0, or in the range F small /ε 2 ≤ x 2 j ≤ F large . We use a data type, double-FP, consisting of two floating-point numbers. We denote them, for example, by A = [A, a]. A is the double-FP variable, and the actual pairs of floating-point numbers are A and a. They have the characteristics •(A + a) = A, which means that a is a "tail" part to add extra precision to A. The mathematical sum of the two components represents a value of at least twice the precision of the underlying floating-point number. Here is the core algorithm for computing the sum of squares j x 2 j . We make the crucial distinction between an assignment operation ":=" in an algorithm and the mathematical equality sign "=". shows that the relative error is bounded by ∆ n-1 (δ) where ∆ (δ) = δ/(1 -δ) and δ is the relative error bound on the underlying addition operation, which is the SumNonNeg function. Theorem 4.1 shows that this δ is 3ε 2 . From that we can deduce the length limit of x within which |(S + s) -σ| < εσ/8.

/ A = [A, a], B = [B, b] nonnegative: A + a, B + b ≥ 0 H := TwoSum(A, B) // H = [H, h], H + h = A + B exactly c := a ⊕ b // c = a + b + δc d := h ⊕ c // d = h + c + δ d . S := FastTwoSum(H, d) // S = [S, s], S + s = H + d exactly //
If F small /ε 2 ≤ σ ≤ F large /2, then |(S + s) -σ| ≤ 3ε 2 σ.
PROOF. The theorem clearly holds if one of α or β is zero. Moreover, it is clear that SumNonNeg is insensitive to the order of its two input arguments. It suffices therefore to consider α ≥ β > 0. There are only two rounding errors in the entire function:

δ c = c -(a + b) = •(a + b) -(a + b) and δ d = d -(h + c) = •(h + c) -(h + c). More precisely, S + s = H + d and H + d = H + h + a + b + δ c + δ d = σ + δ c + δ d . Thus |(S + s) -σ| ≤ |δ c | + |δ d |. The rest of the proof establishes the fact |δ c | + |δ d | ≤ 3ε 2 σ.
We use this fact heavily: For any real number Because σ ≥ F small /ε 2 , ulp(ulp(σ/2)) = ulp(σ)ε ≤ 2σε 2 . Hence

γ ∈ R, |•(γ) -γ| ≤ ulp(γ)/2. Let α = 2 eα (1 + f α ), β = 2 e β (1 + f β ), and σ = 2 eσ (1 + f σ ) where 0 ≤ f α , f β , f σ < 1. σ ≥ α ≥ β =⇒ ulp(σ) ≥ ulp(α) ≥ ulp(β). Because |a| ≤ ulp(α)/2, |a| ≤ ulp(σ)/2. Similarly, |b| ≤ ulp(σ)/2. Therefore, A + B = σ -(a + b) ≤ σ + ulp(σ), which implies ulp(A + B) ≤ 2ulp(σ).
|δ c | ≤ σε 2 . ( 4
)
We now show that |δ d | ≤ 2σε 2 . Consider the case of e α ≥ e β + 1. We have |a| ≤ ulp(σ)/2 and |b| ≤ ulp(σ)/4. Therefore,

|δ d | = | • (h + c) -(h + c)| ≤ 1 2 ulp(h + c) ≤ 1 2 ulp(ulp(A + B)/2 + |c|), |δ d | ≤ 1 2 ulp(ulp(σ) + |c|). (5
|c| = | • (a + b)| ≤ |a + b|(1 + ε) ≤ 3 4 (1 + ε)ulp(σ).
Equation 5 implies

|δ d | ≤ 1 2 ulp(ulp(σ) + 3 4 (1 + ε)ulp(σ)) = 1 2 ulp(ulp(σ)), |δ d | ≤ 2σε 2 . ( 6
)
Consider the case of e α = e β . In this situation, we must have e σ = e α + 1 and ulp(α) = ulp(β) = ulp(σ)/2. As a result, |a| + |b| ≤ ulp(σ)/2 and |c| ≤ (1 + ε)ulp(σ)/2, so

|δ d | ≤ 1 2 ulp(ulp(σ) + 1 2 (1 + ε)ulp(σ)) = 1 2 ulp(ulp(σ)), |δ d | ≤ 2σε 2 . ( 7
)
Equations 4, 6 and 7 together show that |δ c |+|δ d | ≤ 3σε 2 , and the theorem is proved.

THEOREM 4.2. Let n be the length of a vector x in safe range and σ denote j x 2 j . Let SumOfSquares(x) return the result [S, s]. Then

|(S + s) -σ| ≤ ∆ n-1 (3ε 2 )σ where ∆ (δ) = δ/(1 -δ). In particular, if the length n satisfies n < ((24 + ε)ε) -1 , then |(S + s) -σ| < εσ/8.
PROOF. From Theorem 4.1 and standard error bound on adding n nonnegative floating-point types [START_REF] Nicholas | Accuracy and stability of numerical algorithms[END_REF]] with an addition operation of relative error bounded by 3ε 2 ,

|(S + s) -σ| ≤ ∆ n-1 (3ε 2 )σ.
Because ∆ (δ) is an increasing function in in the range 0 ≤ < 1/δ, for n < ((24 + 3ε)ε) -1 ,

∆ n-1 (3ε 2 ) < ∆ n (3ε 2 ) < ∆ L (3ε 2 ) = ε 8 ,
where L = ((24 + 3ε)ε) -1 .

That maximum vector length bound L = ((24 + 3ε)ε) -1 corresponds to L = 699050 for IEEE754 binary32 and L ≤ 3.76 • 10 14 for IEEE754 binary64.

We parallelize SumOfSquares in an obvious manner: partition the input vector x to τ subvectors of roughly equal length. Perform the sum of squares on each subvector in parallel. This parallelism can be realized either at the thread level or data level, the latter using SIMD vector instructions such as SSE or AVX. The partial sums of squares are then accumulated in a serial manner.

function SumOfSquaresP(x) // Parallel SumOfSquares Partition x into τ portions, x (t) , t = 1, 2, . . . , τ // length of each x (t) is no more than m = n/τ . S (t) := SumOfSquares(x (t) ), t = 1, 2, . . . , τ . // In parallel, each S (t) = [S (t) , s (t) ] is a double-FP. S := [0, 0]; S := SumNonNeg(S, S (t) ), t = 1, 2, . . . , τ . // In serial, summing the τ partial sums of squares // S = [S, s] at this point; S + s ≈ n j x 2 j . return S end SumOfSquaresP THEOREM 4.3. Let n be the length of x and S = [S, s] be the result of SumOfSquaresP(x) with τ portions and m = n/τ . Then

|(S + s) -σ| ≤ ∆ m+τ (3ε 2 )σ.
In particular,

|(S + s) -σ| ≤ ∆ n-1 (3ε 2 )σ whenever m + τ ≤ n -1.
PROOF. We document the two stages of errors with the following notations. Let σ (t) = (x (t) ) T x (t) , t = 1, 2, . . . , τ , and σ = τ t=1 σ (t) denote the exact partial and exact complete inner products, respectively. Let

σ (t) = S (t) + s (t) , approximate partials 1 ≤ t ≤ τ , σ = τ t=1 σ (t) ,
exact sum of approximate partials, σ = S + s, approximate sum of approximate partials.

For the computed partials, where we are summing no more than m = n/τ double-FP types, we have

| σ (t) -σ (t) | ≤ ∆ m-1 (3ε 2 ) σ (t) , t = 1, 2, . . . , τ,
and t) .

τ t=1 ( σ (t) -σ (t) ) ≤ ∆ m-1 (3ε 2 ) τ j=1 σ (
This implies

| σ -σ| ≤ ∆ m-1 (3ε 2 ) σ, (8)
σ ≤ (1 + ∆ m-1 (3ε 2 )) σ. (9)
Similarly,

| σ -σ| ≤ ∆ τ -1 (3ε 2 ) σ. (10)
Combining Equations 8 through 10,

| σ -σ| σ ≤ ∆ τ -1 (3ε 2 ) ( σ/σ) + ∆ m (3ε 2 ) ≤ ∆ τ -1 (3ε 2 ) (1 + ∆ m-1 (3ε 2 )) + ∆ m (3ε 2 ) ≤ ∆ τ (3ε 2 ) + ∆ m (3ε 2 ) (11) ≤ ∆ m+τ (3ε 2 ). (12)
Equation 11 follows from its preceding line as long as 3ε 2 ≤ (n + m + 2) -1 , and Equation 12 follows from Equation 11because ∆ (3ε 2 ) + ∆ (3ε 2 ) ≤ ∆ + (3ε 2 ). Finally, whenever m + τ ≤ n -1, ∆ m+τ (3ε 2 ) ≤ ∆ n-1 (3ε 2 ). The proof is now complete.

We remark that the number of threads τ is typically much smaller than the vector length n. In a common scenario, τ equals the number of cores which is in the order of 10 or so. Moreover, threading is beneficial only when there is enough work per thread, implying m = n/τ is in the order of 100 or more. Thus, m + τ ≈ n/τ , implying m + τ ≤ n -1. This says that a parallel sum of squares is in general more accurate then the serial version. In particular, as long as n < ((24 + ε)ε) -1 and [S, s] is obtained from with SumOfSquares or SumOfSquaresP, an IEEE conforming square root evaluation sqrt(S) produces a faithfully rounded x 2 for x whose elements fall in the core range discussed here.

FAITHFUL 2-NORM -GENERAL CASE

That the simple accumulation of σ = n j=1 x 2 j is susceptible to spurious exceptions can be illustrated by the simple example of n = 8, x j = •(2 F large ) for j ≤ 4 and x j = •( √ F small /2), j > 4. While x 2 is approximately 4 F large , overflows in computing x 2 j for j ≤ 4 leads to a computed σ of +∞, rendering the final computed l 2 -norm completely wrong. This is why general-purpose software such as LAPACK's public release essentially computes x j (x j /x) 2 where x = max j |x j |. While this strategy resolves the spurious overflow problem satisfactorily, spurious underflows can still be triggered. If underflow is masked, spurious underflow is harmless to the final numerical results. Nevertheless, these spurious underflows may significantly degrade performance on computing platforms that handle underflow via a trapping mechanism.

We present here an algorithm that returns a value in clip(♦( x 2 )) and reports overflows and underflows faithfully as defined in Section 2. Let σ = j x 2 j denote the sum of squares. Our algorithm first computes Z, a faithful rounding of a scaled l 2 -norm

x 2 = γ -m/2 x 2 , that is Z ∈ ♦( x 2 ).
The factor γ -m/2 is chosen so that γ m is an even power of 2, γ -m/2 ∈ F, and x 2 ∈ [F small , F large ]. We shall describe a way to choose γ and to compute m below. The final result is returned, naturally, as γ m/2 ⊗ Z. Theorem 5.1 establishes rigorously that not only does the numerical value γ m/2 ⊗ Z ∈ clip(♦( x 2 )), but the multiplication also reports overflow and underflow faithfully. The remaining of this section then focuses on the computation of Z ∈ ♦( x 2 ). THEOREM 5.1. Let ζ ∈ [F small , F large ]∪{0} be a real number and Z be a floating-point number where Z ∈ ♦(ζ). Let t be an integer where 2 t ∈ F. Then the IEEE multiplication 2 t ⊗ Z satisfies 2 t ⊗ Z ∈ clip(♦(2 t ζ)) and 2 t ⊗ Z reports overflow and underflow faithfully as an implementation of 2 t ζ.

PROOF. 2 t ⊗ Z ∈ clip(♦(2 t ζ)) follows easily if •(2 t Z) ∈ ♦(2 t ζ), which is what we will prove. The case of ζ = 0 is trivial as ♦(ζ) = {0}, implying Z = 0 as well. Obviously •(2 t Z) ∈ ♦(2 t ζ). It therefore suffices to consider ζ ∈ [F small , F large ]. There is a unique integer e, e ≥ e min , such that ζ ∈ [2 e , 2 e+1 ). Z ∈ ♦(ζ) implies that Z ∈ [2 e , 2 e+1
] and |ζ -Z| < ulp(2 e ). There are only two possibilities: t + e ≥ e min and t + e < e min . If t + e ≥ e min , we have

•(2 t Z) = 2 t Z ∈ {2 t α | α ∈ ♦(ζ)} = ♦(2 t ζ).
If t + e < e min , both 2 t ζ and 2 t Z lie inside [2 t+e , 2 t+e+1 ] and

|2 t ζ -2 t Z| < 2 t ulp(2 e ) ≤ ulp(2 t+e )/2. By virtue of Lemma 3.1, •(2 t Z) ∈ ♦(2 t ζ). Now that W = •(2 t Z) ∈ ♦(2 t ζ), it is clear that given any Y ∈ F , 2 t ζ ≤ Y implies W ≤ Y , and 2 t ζ ≥ Y implies W ≥ Y . Therefore, if 2 t ζ ∈ [F small , F large ], 2 t ⊗ Z will not report overflow or underflow. If 2 t ζ ≥ 2 emax+1
, 2 t ⊗ Z will definitely report an overflow. Finally, observe that if 2 t ζ ∈ F for some t < 0, then ζ ∈ F and thus Z ∈ ♦(ζ) = {ζ}, implying that Z = ζ and 2 t Z = 2 t ζ. Consequently, whenever 2 t ζ ≤ F small -2 emin+1 ε, 2 t ⊗ Z will report underflow faithfully as defined in Section 2. This completes the proof.

We turn now to the computation of a scaled l 2 -norm. A known strategy [START_REF] Blue | A portable Fortran program to find the Euclidean norm of a vector[END_REF]] partitions the input data into three bins: one for small inputs, one for large and one for "medium" inputs. The data in each bin are scaled by a common, statically chosen scale factor so that sums of squares of elements in each bin incur no spurious underflow or overflow exceptions. The final result is constructed by appropriate combination of the partial sums-of-squares from the bins. We follow the same approach but with two enhancements. First, while we will explain our approach using three bins, we will point out later that our implementation keeps only two bins of data at any given time. This is important as the bins are in practice kept in the scarce SIMD vector registers3 . Second, our combination of the binned partial sums of squares are done in a way that guarantees a faithfully rounded scaled l 2 -norm.

The basic idea is to divide the entire input range of |x j | into three "equal" subranges, where sums of squares of data in the middle (interior) subrange does not generate exceptions. Data in the two exterior ranges are scaled into the interior subrange. Now the specifics. Define the even integer E by

E = min{e | 3e ≥ e max -e min -log 2 (ε), e is even }.
From E, we define a scale factor γ = 2 -E and use the following notations.

E hi = e max + 1 -E, β hi = 2 E hi , E lo = e max + 1 -2E, β lo = 2 E lo .
In particular, γβ hi = β lo . We tabulate the specific values for binary32 and binary64 here.

E E hi E lo β hi β lo binary32 94 34 -60 2 -26 = ε/4 binary64 700 324 -376 2 -52 = 2ε
Given the input vector x = [x 1 , x 2 , . . . , x n ] T , the three bins are

A = { γx j | |x j | ≥ β hi }, B = { x j | β lo ≤ |x j | < β hi }, C = { x j /γ | |x j | < β lo }.
By design

β lo ≤ | x j | < β hi for x j ∈ A ∪ B ∪ C.
Denote the partial, scaled, sums-of-squares as

σ A = xj ∈A x 2 j , σ B = xj ∈B x 2 j , and σ C = xj ∈C x 2 j .
Clearly, the "bin-sums" are in the range

σ A , σ B , σ C ∈ {0} ∪ [β 2 lo , nβ 2 hi ) ⊆ {0} ∪ [β 2 lo , β 2 hi /ε), (13)
by assuming n ≤ 1/ε. Furthermore,

σ = j x 2 j = γ -2 σ A + σ B + γ 2 σ C . (14)
A straightforward implementation can collect all three bins, and invoke the SumOfSquares function (or the parallel version) on each of the bins, followed by some appropriate combination method to arrive at the final result. Our implementation in fact only keeps and processes two bins. The observation is that A and C are never needed simultaneously. If A is nonempty, then Equation 13 and the fact that

γβ hi = β lo show that γ 2 σ C /(γ -2 σ A ) ≤ γ 2 /ε ε 2 .
Neglecting γ 2 σ C altogether incurs a relative error (much) less than ε 2 . A similar estimate shows that keeping σ B is nevertheless necessary in order not to loose too much accuracy.

Briefly speaking, our implementation starts the binning process by keeping the interior (middle) bin and one exterior bin. When the first element belonging to A appears, the existing exterior bin is replaced with that element, and elements from C is never collected from that point onwards. Let us mention that the logic required to decide whether the first element belonging to A has already appeared or not can be implemented using nothing but masks, not requiring branches. Denote the two actually maintained bins by U and V, then

σ 2 = γ k ( σ U + γ 2 σ V ),
where k = -2 if U corresponds to A, and k = 0 if U corresponds to B. In either case, the "two-bin" sum of squares σ 2 satisfies |σ 2 -σ| ≤ ε 2 σ, and

σ 2 ≤ (1 + ε 2 )σ. (15)
To compute σ 2 , the SumOfSquares function is applied to each of the two resulting bins, yielding two double-FP variables U = [U, u] and V = [V, v]. Both U + u and V + v approximate their targets with high relative accuracies:

|(U + u) -σ U | σ U , |(V + v) -σ V | σ V ≤ ∆ n-1 (3ε 2 ),
where ∆ (δ) = δ/(1 -δ). Using n instead of n -1 for simplicity, we have

γ k [(U + u) + γ 2 (V + v)] -σ 2 ≤ ∆ n (3ε 2 )σ 2 , (16)
and

γ k [(U + u) + γ 2 (V + v)] ≤ (1 + ∆ n (3ε 2 ))σ 2 . ( 17
) We handle γ k [(U + u) + γ 2 (V + v)] as follows. For nonzero U and V , β 2 lo ≤ U + u, V + v < β 2 hi /ε (as we assume n ≤ 1/ε). If U ≥ β 2 lo /ε 3 , γ 2 (V + v) < γ 2 β 2 hi /ε = β 2 lo /ε ≤ ε 2 U ≤ ε 2 (1 + ε)(U + u). Similarly, if V ≤ β 2 lo ε 2 /γ 2 = β 2 hi ε 2 , γ 2 (V + v) ≤ γ 2 (1 + ε)V < β 2 lo ε 2 (1 + ε) ≤ ε 2 (1 + ε)(U + u).
Stage 1 and Stage 2 errors have already been discussed in Equations 15 through 17.

Putting these together,

|γ m (S + s) -σ|/σ ≤ γ m (S + s) -γ k [(U + u) + γ 2 (V + v)] /σ + γ k [(U + u) + γ 2 (V + v)] -σ 2 /σ + |σ 2 -σ|/σ, ≤ 5ε 2 (1 + ∆ n (3ε 2 ))(1 + ε 2 ) + ∆ n (3ε 2 )(1 + ε 2 ) + ε 2 , ≤ ∆ n (3ε 2 ) + 7ε 2 .
Consequently,

|(S + s) -σ| ≤ (∆ n (3ε 2 ) + 7ε 2 ) σ, ≤ (∆ n (3ε 2 ) + 9ε 2 ) σ, ≤ ∆ n+3 (3ε 2 ) σ. (21)
From Theorems 4.2 and 4.3, n + 3 < ((24 + 3ε)ε) -1 implies

|(S + s) -σ| < ε σ/8, a condition that guarantees, by Theorem 3.3, that •( √ S) ∈ ♦( √ σ).
AccuNrm2 below is the straightforward synthesis of the previous discussions. Theorem 5.3 that follows is a formal statement that summarizes the technical results of this paper.

function AccuNrm2(x) // general faithful l 2 -norm (m, S) := SumOfSquaresBins(x)
// m is an integer in the range [-2, 2] and γ m (S + s) ≈ j x 2 j // By design, γ m is an even power of 2. Z := sqrt(S) return γ m/2 ⊗ Z // Value in clip(♦( x 2 )) and reports overflow/underflow faithfully (Theorem 5.3) end AccuNrm2 THEOREM 5.3. Let x be a vector of length n. If n < L with L = ((24 + 3ε)ε) -1 -3, then AccuNrm2(x) ∈ clip(♦( x 2 )) and reports overflow and underflow faithfully.

PROOF. This is a direct consequence of Theorems 5.1 and 5.2.

The bound L on the vector length n induced by 5.3 translates as follows for IEEE754 binary32 and binary64:

Vector length bound n < L binary32 L = 699047 binary64 L = 3.75299968947538 • 10 14

IMPLEMENTATION AND TESTING

The complete set of codes, together with testing and performance measurement auxiliary sources, is available at under an open source license.

We implemented and tested our faithfully rounded, division-free l 2 -norm with faithful reporting of underflow and overflow. The implementation referred as FaithfulNorm closely follows the algorithmic description given in the previous sections.

We used IEEE 754 binary64 as working-precision and we restricted ourselves to a SIMD environment, targeting in particular Intel SSE/AVX units, with or without support for the IEEE754 fused-multiply-and-add (FMA) instruction. Recent versions of SSE and all versions of AVX support IEEE 754 binary64 precision. The rationale for the choice of a SIMD environment is twofold: to use an environment most similar to the existing codes we compare our algorithm to, and to maximize our performance in a typical processor without the use of threads.

To achieve high performance on modern pipelined floating-point units, it is important to avoid branching (when possible) as well as avoid the use of expensive operations such as floating-point division. By design, our l 2 -norm algorithm is division free. We are also able to make our inner loop branching free based on three observations. First, the SSE/AVX units offer comparison instructions that return their results as masks of all-ones or all-zeros. Second, logical bit-and a floating-point variable with all-ones leave the variable unchanged while bit-and with all-zeros turn it into a floating-point value of zero. As a matter of course, multiplying floating-point zeros and accumulating them is innocuous. Third, discarding the C bin when the first A is found and maintaining in these registers the U bin from that point onward (cf. Section 5) just means maintaining a binary flag, which can also be implemented as a bit-mask.

We shall repeat that our implementation has no memory overhead nor memory access overhead: each input vector element x j is read only once and the intermediate values (accumulators etc.) are kept in registers.

We compared the implementation of our faithfully rounded l 2 -norm with implementations for other approaches with respect to both accuracy and performance. To do so, we implemented a naïve l 2 -norm, called NaiveNorm, that plainly uses working precision for squaring the x i and accumulating these squares, without any underflow and overflow avoidance. We further implemented the algorithm found in netlib [START_REF] Anderson | LAPACK Users' guide[END_REF]; we call this implementation NetlibNorm. Finally, we implemented another faithfully rounded l 2 -norm using the arbitrary precision library MPFR [MPFR ]. This implementation is simply based on an exact accumulation of the squares x 2 i in an accumulator that provides enough precision: using an MPFR variable with precision p = 2 (e max -e min -log 2 ε) + log 2 n is just enough. We refer to it as MPFRNorm.

We performed testing on a 4-core Intel Core i7 at 2.67 GHz with 4Gb of RAM and on a 8-core Intel Xeon E3-1275 v3 at 3.50 GHz with 32Gb of RAM. All implementations were written in C -using builtins for access to SIMD instructions-and compiled using gcc version 4.8 and options -std=c99 -O3 -march=native. Timings are given cycles per vector element, obtained using the Read-Time-Step-Counter instruction with serialization, subtracting off the measured overhead for a call to an empty function and dividing by the number of elements.

We used pseudo-random floating-point input vectors in our tests. These pseudorandom values were constructed as follows: we separately generated a uniformly distributed exponent value in range and a uniformly distributed significand for that chosen exponent value. We then constructed a floating-point value out of this exponent and significand value. When generating values for a subdomain We did accuracy testing with test vectors of various lengths n and input types. The testing results are resumed in Table I.

First, we considered input vectors chosen such that the final l 2 -norm result is a normal floating-point number. Second, we tested the algorithms on input vectors for which the final result gradually underflows. Third, we performed testing on vectors with inputs around 1.0, i.e. where underflow or overflow avoidance is not necessary. Finally, we constructed input vectors with x 1 = 1 and subsequent x j s are chosen to be much smaller than 1 but with the contrived property that •(1 + x 2 j ) produces a positive absolute error very close to the "half-ulp" bound of ε. This test case is admittedly artificial, but nonetheless demonstrates a near worse-case error scenario for NaiveNorm and NetlibNorm.

Testing shows that both the NaiveNorm and NetlibNorm implementations fail to provide faithfully rounded results. It demonstrates also that both our FaithfulNorm as well as the MPFRNorm algorithm do yield faithfully rounded results.

In cases when the final l 2 -norm does not overflow, our algorithm FaithfulNorm returns a result with an error well below 1ulp. As expected, the maximum error does not vary with vector length, whereas it does for the netlib l 2 -norm.

Accuracy testing also shows that with random inputs, netlib l 2 -norm can result in hundreds of ulps of error when the vectors lengths n get to be ten million or more. In deliberately constructed inputs, errors as large as about 0.25 n ulp can be observed.

Turning now to performance testing, Tables II, III, IV and V summarize our observations. We used vectors of floating-point numbers in various domains of interest. We tested for vectors for which the final l 2 -norm results gradually underflow, overflow or stay in the range of normal floating-point numbers. Measurements were done for varying vector lengths. Starting with some minimal vector length (a couple of dozen ele- ments), vector length had no influence on computation time per element; we therefore report only the numbers obtained for vectors of length 10 6 .

These performance results speak in favor of our FaithfulNorm implementation. For cases when the final result is a normal floating-point number, our implementation is up to 3 times faster than the netlib implementation. As already explained, this is due to several factors: avoidance of spurious underflow, no use of expensive divisions and an algorithm that is branch-free in the inner loop. In particular the relative cost of the divisions and branches in the netlib implementation can be seen in the performance data: on Intel Core i7, which does not yet implement the recent AVX extensions, our algorithm is up to 4.5 times faster than netlib whereas on Intel Xeon E3-1275, the same SSE codes run equally fast. However, on Intel Xeon E3-1275, AVX and FMA are available, allowing our l 2 -norm to be up to 3 times faster than netlib. It also worth mentioning that spurious underflow in netlib hurts netlib performance on some processors -such as the Intel Core i7-but does not on others, such as the Intel Xeon E3-1285.

We shall however mention that our algorithm does have lower performance than the netlib in two cases. First, for inputs when the vector length is (very) short -typically less than a dozen elements. In this case our algorithm has a much higher static overhead due to the elaborate computations needed in the reduction of the bins and square root. In future work we shall address this problem with a call-out to a specialized l 2norm for very small vectors. Second, the netlib norm can be faster on some processors that have a faster floating-point division instruction (with respect to multiplication) and that do not suffer a performance impact due to branching nor on subnormal handling. This effect can already be measured on recent Intel Xeons. But we point out that these processors come equipped with the FMA instructions which FaithfulNorm can exploit. Table V shows FaithfulNorm regains the speed advantage when it uses the FMA instructions on these processors appropriately.

CONCLUSIONS

In this paper, we presented an efficient algorithm to compute the faithful rounding of the l 2 -norm of a floating-point vector. While our algorithm is very accurate, it is also faster than previous algorithms like the one of netlib that gives no information about the accuracy of the result. Moreover, our algorithm avoids spurious overflow and underflow. It is also suitable for parallel implementations. We have hitherto focused our implementation on vector-parallelism using SIMD instructions. Implementation and testing of our algorithm in a threaded environment as well as formulating the algorithm in terms of auto-vectorizable, auto-parallelizable code is left to future work.

Disclaimers

Software and workloads used in performance tests may have been optimized for performance only on Intel microprocessors. Performance tests, such as SYSmark and MobileMark, are measured using specific computer systems, components, software, operations and functions. Any change to any of those factors may cause the results to vary. You should consult other information and performance tests to assist you in fully evaluating your contemplated purchases, including the performance of that product when combined with other products. Configurations: Refer to table 1. For more information go to http://www.intel.com/performance Intel's compilers may or may not optimize to the same degree for non-Intel microprocessors for optimizations that are not unique to Intel microprocessors. These optimizations include SSE2, SSE3, and SSE3 instruction sets and other optimizations. Intel does not guarantee the availability, functionality, or effectiveness of any optimization on microprocessors not manufactured by Intel. Microprocessor-dependent optimizations in this product are intended for use with Intel microprocessors. Certain optimizations not specific to Intel microarchitecture are reserved for Intel microprocessors. Please refer to the applicable product User and Reference Guides for more information regarding the specific instruction sets covered by this notice. Notice revision #20110804

function

  SumOfSquares(x) // Accurate accumulation S := [0, 0] for j = 1, 2, . . . , n do: P := TwoProd(x j , x j ) // P = [P, p], P + p = x 2 j exactly S := SumNonNeg(S, P) return S end SumOfSquares function SumNonNeg(A, B) // [A, a] + [B, b] // error bound of this operation 3ε 2 (Theorem 4.1) /

  THEOREM 4.1. Let S = [S, s] be the result from applying SumNonNeg on nonnegatives A = [A, a] and B = [B, b]. Let α = A+a ≥ 0, β = B +b ≥ 0 denote the exact input values, and σ = α + β denote the exact sum.

  We note that |δ c | = |•(a+b)-(a+b)| ≤ ulp(a+b)/2 and |a+b| ≤ ulp(σ). But |a+b| = ulp(σ) only when |a| = |b| = ulp(σ)/2, which implies a + b is representable exactly in F and δ c = 0. When |a + b| < ulp(σ), we have ulp(a + b) ≤ ulp(ulp(σ/2)). Hence using basic properties of the ulp function stated in Section 2,

)

  To complete the estimate on |δ d |, we analyze |c|. There are only two possibilities: either e α ≥ e β + 1 or e α = e β . We show that each situation leads to |δ d | ≤ 2σε 2 .

  [a; b] ⊆ F, we per-formed that random-generation process for an exponent range completely covering the possible exponents of floating-point values in [a; b], discarding all generated floatingpoint values that were outside of [a; b].

  then 2e+1 ≥ e min and Lemma 3.2 shows that S ∈ [2 2e+1 , 2 2e+2 ] and |s| ≤ 2 2e+1 ε.

	Summarizing, both [2 e , 2 e+1 ] with e ≥ e min . By virtue of Lemma 3.1, we can establish the fact that •( √ σ and √ S are in √ S) ∈ ♦( √ σ) provided we can show | √ S -√ s| < ulp(2 e )/2. Since e ≥ e min , ulp(2 e ) = 2 e+1 ε, this is equivalent to establishing | √ S -√ σ| <

  see Section 2 for TwoSum and FastTwoSum. return S end SumNonNeg Theorem 3.3 guarantees that if [S, s] returned by SumOfSquares satisfies |(S + s) -σ| < εσ/8, then •( √ S) is a faithful rounding of √ σ = x 2 . Since SumOfSquares is summing n double-FP type, standard error analysis (see for example Chapter 3 of [Higham 2002])

Table I .

 I Maximum error in ulps observed for various domains and vector lengths n, plain SSE implementation

		vectors with	vectors for which	vectors with	vectors with chosen
		normal results	results underflow	entries around 1.0	"half-ulp" entries
		n = 10 3	n = 10 7	n = 10 3	n = 10 7	n = 10 3	n = 10 7	n = 10 3	n = 10 7
	NaiveNorm	∞	∞	8.84 • 10 12	5.46 • 10 10	7.73	861	250	2.50 • 10 6
	NetlibNorm	2.01	524	0.496	0.698	7.58	609	250	2.50 • 10 6
	MPFRNorm	0.494	0.481	0.490	0.498	0.468	0.497	0.0749	0.484
	FaithfulNorm	0.620	0.628	0.497	0.499	0.605	0.701	0.0749	0.484

Table II .

 II Computation time in cycles per vector element, plain SSE version on Intel Core i7

		vectors with	vectors for	vectors with	vectors for	vectors
		normal	which results entries around	which results provoking spuri-
		results	underflow	1.0	overflow	ous underflow
						in NetlibNorm
	NaiveNorm	47	137	3.48	46.8	128
	NetlibNorm	156	472	19.1	156	274
	MPFRNorm	1080	2670	818	1090	1660
	FaithfulNorm	34.2	289	25.3	34.2	62.2

Table III .

 III Computation time in cycles per vector element, plain SSE version on Intel Xeon E3-1275

		vectors with	vectors for	vectors with	vectors for	vectors
		normal	which results entries around which results provoking spuri-
		results	underflow	1.0	overflow	ous underflow
						in NetlibNorm
	NaiveNorm	4.95	4.75	4.72	4.70	4.52
	NetlibNorm	21.9	158	12.8	21.1	21.8
	MPFRNorm	810	1160	536	803	717
	FaithfulNorm	21.5	87.3	21.8	21.7	20.3
	Table IV. Computation time in cycles per vector element, AVX version w/o FMA on Intel Xeon E3-1275
		vectors with	vectors for	vectors with	vectors for	vectors
		normal	which results entries around which results provoking spuri-
		results	underflow	1.0	overflow	ous underflow
						in NetlibNorm
	NaiveNorm	4.85	4.61	4.68	4.86	4.52
	NetlibNorm	21.1	157	13.3	21.6	21.8
	MPFRNorm	795	1250	552	765	720
	FaithfulNorm	12	50.7	12.5	12.6	14.8

Table V .

 V Computation time in cycles per vector element, AVX version using FMA on Intel Xeon E3-1275

		vectors with	vectors for	vectors with	vectors for	vectors
		normal	which results entries around which results provoking spuri-
		results	underflow	1.0	overflow	ous underflow
						in NetlibNorm
	NaiveNorm	4.52	4.52	4.52	4.52	4.52
	NetlibNorm	20.5	151	12.6	20.5	22
	MPFRNorm	722	1110	481	723	770
	FaithfulNorm	6.94	42.3	6.94	6.94	10.4

ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY.

Clearly, a l 2 -norm routine that returns the correctly rounded floating-point result is always reproducible.ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY.

It is not possible to go down to one bin only as long as the bin boundaries and scaling factors are chosen statically to avoid division operations.ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY.

http://www.christoph-lauter.org/faithfulnorm.tgz ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY.

Part of this work was done while Prof. ShinStef Graillat was supported by a CNRS/JSPS "Exchange Scientist" grant.

Thus if U ≥ β 2 lo /ε 3 or V ≤ β 2 hi ε 2 , the error by dropping the V term is in the order of ε 2 :

If U < β 2 lo /ε 3 and V > β 2 hi ε 2 , then neither PROOF. We group the total errors incurred in computing σ as γ m σ into three stages. In Stage 1, the value σ, which is exactly represented in terms of γ and the three binsums (Equation 14), is approximated by σ 2 = γ k ( σ U + γ 2 σ V ). In Stage 2, the two binsums σ U and σ V are approximated by the double-FP [U, u] and [V, v]. Finally, in Stage 3, γ k ((U + u) + γ 2 (V + v)) is approximated as γ m (S + s) by possibly dropping v or both V and v and the use of SumNonNeg.

Consider the Stage-3 error. There are three possible points of exit in the procedure SumOfSquaresBins. The first point of exit corresponds to a zero Stage-3 error as there is actually only at one nonempty bin. The second point of exit corresponds to Stage-3 error bounded by ε 2 (1 + ε) as given by Equation 18. If the last point of exit is taken, Stage-3 error consists of one part that is due to a single application of SumNonNeg, which is bounded by 3ε 2 (Theorem 4.1), and one due to possibly dropping the v term, which is bounded by ε 2 (1 + ε) (Equation 18). Thus we bound the error in Stage 3 conservatively by 5ε 2 :

and


Christoph Lauter
email: christoph.lauter@lip6.fr

Easing Development of Precision-Sensitive Applications with a Beyond-Quad-Precision Library

The IEEE754 Standard offers essentially two binary floating-point formats, binary32 and binary64, natively supported by current hardware. Whenever these precisions do not suffice, developers are restricted to arbitrary precision libraries, such as MPFR.

These libraries however leave a gap in the mid-precision range (64 to 512 bits). Their very nature -as object code to be linked inprevents modern compilers from inlining the code or optimizing it e.g. with loop unrolling.

We propose the libwidefloat software meant to fill this gap. It offers precisions from 64 through 512 bits. It supports all basic operations (+, -, ×, /, √ , FMA, comparisons etc.) It is fully implemented in header files for automatic optimization.

I. INTRODUCTION

The IEEE754 Standard for Floating-Point (FP) Arithmetic defines several binary and decimal FP formats [1]. The basic ones, such as binary32 (single precision) and binary64 (double precision), are natively supported in hardware by current general-purpose processors. The IEEE754 binary128 (quad precision) format is supported in the current FP environment mainly through software emulation. The additional extended formats defined by IEEE754, with precisions beyond quadprecision, have no practical realization [1].

While binary32 and binary64 gear the majority of numerical software, some applications need support for precisions beyond quad-precision, i.e. with more than 113 bits of significand length. Examples for such applications include astronomy [START_REF] Laskar | A long-term numerical solution for the insolation quantities of the earth[END_REF], [START_REF] Laskar | High precision astronomical solution for paleoclimate studies[END_REF], [START_REF] Farrés Basiana | High precision symplectic integrators for the solar system[END_REF], computational physics [START_REF] Wittig | Rigorous high precision interval arithmetic in cosy infinity[END_REF] or experimental mathematics [START_REF] Bailey | High-precision computation and mathematical physics[END_REF]. In some cases the use of extended precision might be overcome by subtle FP accuracy analysis and FP tricks [START_REF] Graillat | Contribution à l'amélioration de la précision et à la validation des algorithmes numériques[END_REF], [START_REF] Boldo | Formal verification of tricky numerical computations[END_REF]. However, some developers may simply prefer to speed up development with extended precision, used in some compute kernels, even if this means breaking a butterfly on a wheel [START_REF] Thomas | A general-purpose method for faithfully rounded floatingpoint function approximation in fpgas[END_REF].

As we will discuss in more detail below, such applications typically resort to arbitrary precision FP libraries, such as MPFR 1 or MPFUN 2 . These libraries are properly designed and optimized: they offer very good, if not optimal, performance in the precision range where arbitrary precision overhead becomes negligible [START_REF] Ghazi | Why and how to use arbitrary precision[END_REF]. This range typically starts at precisions beyond a couple of hundreds to a thousand bits. This leaves a gap in the mid-precision range between quad-precision and a couple of hundred bits.

In this article, we propose the libwidefloat software designed to filled this gap. The package offers FP types with precisions between 64 and 512 bits, to be chosen at application compile time in 32 bit steps. Our FP types come with a 32 bit wide exponent, which essentially eliminates the need for support of subnormal numbers. They come with support for IEEE754-like infinities, Not-A-Numbers, rounding modes and FP flags [1]. The latter two features can be deactivated (resp. set to default round-to-nearest mode) at application compile time to allow for enhanced performance.

This article is organized as follows: in Section II, we first present existing solutions to extended FP precision, analyzing the causes of their shortcomings in the mid-precision range. We then give an overview of modern compiler techniques that we will leverage to overcome these shortcomings (Section III). In Section IV, we present our libwidefloat library, detailing the data-types and operations it supports. We briefly describe the algorithms used in the library in Section V, give experimental results in Section VI and conclude with Section VII.

II. EXISTING EXTENDED PRECISION LIBRARIES AND THEIR SHORTCOMINGS

We have claimed that the existing arbitrary precision libraries, such as MPFR, are very well optimized for higher precisions (starting at about 1000 bits) but hindered by considerable overhead in the low-and mid-precision range. This overhead is due to several factors.

First, the very nature of arbitrary precision means that the length of the limb arrays that hold the FP significands of the different formats is not statically known. As a consequence, all variables having a type defined by one of these arbitrary precision libraries need not only be declared but also be initialized, using a call to a certain initialization function [START_REF] Fousse | Mpfr: A multiple-precision binary floating-point library with correct rounding[END_REF]. This function actually attributes a precision setting to the arbitrary precision variable and calls a memory allocation function such as malloc to provide space for the FP variable's significand. All accesses to a program's arbitrary precision variables' significands are in consequence indirect, requiring pointers to be followed onto the stack. Together with the dynamic nature of the significand length, this mostly prevents compilers from performing any reasonable optimization involving the arbitrary precision FP significands. In addition, for certain applications the requirement to allocate the arbitrary precision data on the memory heap and not hold them on the stack induces important overhead, as memory allocation may have significant cost with respect to computations for smaller precisions.

Second, the fact that existing arbitrary precision software comes packaged as software libraries adds to their overhead. While operations on IEEE754 binary32 and binary64 datatypes are "inlined" as basic operations into the code, each and every use of an arbitrary precision operation involves a function call with all its cost in terms of register usage, variable spills, instruction cache misses etc. In addition, the correspondence between function names and the code behind them gets lost to optimization. Typically, when calling the arbitrary precision library, the compiler does not know of the functions' use cases and when compiling the application code, it does not anything about the library functions but their name. Due to the use of assembly in the underlying libraries such as GMP 3 , the arbitrary precision libraries even push this issue further to a point where the compiler can almost no optimization, even when using Link-Time-Optimization [START_REF] Glek | Optimizing real world applications with GCC link time optimization[END_REF], [START_REF] Fousse | Mpfr: A multiple-precision binary floating-point library with correct rounding[END_REF].

III. LEVERAGING COMPILER OPTIMIZATION TECHNIQUES

In order to take the observations described in the previous Section II into account, the following requirements to a software package for mid-precision FP arithmetic seem sensible.

First, the piece of software should not be packaged as a library to be linked in but should come in the form of code to be directly integrated into the application code. With modern compilers, such as gcc, this can be achieved by condensing the mid-precision FP library into header files that contain the functions provided by the software package declared static and inline. This way, the compiler 4 will be able to inline the FP functions into the application, propagating constants as far as possible, eliminating all unnecessary parts (such as NaN handling where no NaNs are possible) and optimizing loops by unrolling them.

Second, the software package should not use any microoptimizations resorting to inline assembly or linking with objects written in assembly. While it is tempting for performance optimization at a function level, the use of assembly prevents compilers from propagating the information about the program it needs for optimization through the assembly sequence [START_REF] Glek | Optimizing real world applications with GCC link time optimization[END_REF]. In most cases, modern compilers have sufficient support for doubled-precision integer arithmetic, in particular full 64 × 64 bit multiplication, through builtins (such as the __uint128_t type in gcc 5 ), so that resorting to assembly is no longer necessary.

Finally, all wide integer and FP data-types to be declared by a FP software package optimized for the mid-precision range should be statically sized, in order to be held directly on the stack. In addition, the data should be properly aligned in order to avoid cache misses. As a matter of course, the endianess of wide integers, represented as limbs of a statically sized array should match the system's natural endianess, allowing the compiler to make use of up-converting and down-converting memory move operations [START_REF]Intel Corporation, Intel 64 and IA-32 Architectures Software Developer's Manual[END_REF]. In order to help the compiler, all accesses to these data-types involving pointers (or indexing in an array) should be done in a manner that allows for optimization. For instance, functions taking references on such data-types should be declared in a way that makes the compiler understand6 that the pointers are not aliased (using the restrict keyword) and constant (using const) [START_REF]Programming Languages[END_REF], [START_REF]Information technology -Programming languages -C[END_REF].

As a matter of course, the last requirement of using statically sized types in a FP mid-precision software package supporting more than one precision comes at the cost of a combinatorial explosion of functions to be provided. Indeed, if the ease of mixing FP types of different precisions that is provided by arbitrary precision libraries is to be maintained, where an arbitrary precision library provides one function for, say, addition, a software package with statically sized types has to provide several functions, one per combination of input and output types that is possible and sensible. In our work on libwidefloat, we have made extensive use of macroprogramming and code generation to overcome this issue of combinatorial explosion.

IV. THE libwidefloat LIBRARY

The libwidefloat software was implemented with the shortcomings of the existing arbitrary precision libraries and modern compiler technology in mind. The software is implemented with two header files only, using modern C99/C11 [START_REF]Programming Languages[END_REF], [START_REF]Information technology -Programming languages -C[END_REF] but no assembly in order not to break the compiler's opportunities to optimize the code. All functions are declared in a way that allows them to be inlined and then optimized in their use context by the compiler.

The libwidefloat header starts by declaring statically sized data-types for libwidefloat FP values with precisions between 32 and 512 bits, by 32 bit increments. All these types support the basic IEEE754 FP data, such as positive and negative real numbers, signed zeros, signed infinities and Not-A-Number (NaN) data. The significand length is adapted in size according to the type's precision; the exponent width is kept at 32 bits for all precisions. Since this exponent width is pretty wide, we chose not to provide support for subnormal numbers in the libwidefloat types. The NaNs are signed but do not provide support for payload. They are all quiet NaNs; signaling NaNs are not supported.

Even though the libwidefloat types appear to be IEEE754-compliant from a programmer's perspective, their encoding is an ad-hoc one and does not respect the provisions of the IEEE754 Standard on extended precision types [1].

The libwidefloat software has provision for all five IEEE754 exception flags (inexact, invalid, overflow, underflow and division-by-zero) as well as for all four IEEE754 rounding modes (round-to-nearest-ties-to-even, round-down, round-up, round-towards-zero) [1]. The software can be configured at compile time in a way that allows this FP state, made of the flags and the rounding-mode, to be either

• held in a global variable,

• passed as a handle to any of the libwidefloat functions or • omitted, assuming round-to-nearest as the rounding-mode and setting no flags. The code is written in a way such that the compiler is able to eliminate all unnecessary computation steps that are not required when no flags are to be set and the rounding mode is statically set to round-to-nearest.

Currently, libwidefloat includes 1303 functions covering the following types of operations:

• Conversion between all combination of libwidefloat types. The up-conversions are exact, the downconversions are correctly rounded. • Conversion from and to the IEEE754 binary32 and binary64 formats [1]. These conversions are correctly rounded and signal all required exceptions, in particular inexact when rounding is necessary. These exceptions are signaled in and the rounding mode is taken from the FP environment of the type the conversion converts to. For instance, when converting from IEEE754 to libwidefloat, rounding is done according to the libwidefloat rounding mode and not the IEEE754 one.

• Conversion from and to signed and unsigned integer types. These conversions are provided in all variants (signaling inexact or not, taking the rounding mode from an immediate or from the FP state) that are defined by the IEEE754 Standard. • Rounding to integral as defined by the IEEE754 Standard [1]. • Signaling and quiet comparisons, returning the comparison result as a enumerated type [1]. • Homogeneous-type absolute value and negation operations. The heterogeneous variants, where the operands' type differs from the one of the result can be provided through composition with a conversion operation [1]. • Operations to access the libwidefloat FP environment with respect to flags and rounding modes. • And finally, addition, Subtraction, Multiplication, Division, Fused-Multiply-And-Add (FMA) and Square Root with correct rounding. All these operations are provided in a homogeneous variant, where all operands and the result have the same type and in all heterogeneous variants, where the result is of a smaller type, such that all correctly rounded heterogeneous operations required by IEEE754 can be provided through composition with the conversion operations [1].

V. ARITHMETICAL ALGORITHMS USED IN libwidefloat

There is no real novelty in the algorithms used in libwidefloat. In order to manage the combinatorial explosion due to support for the 16 formats covering the precisions from 32 to 512 bits in homogeneous and heterogeneous operations, the final rounding step towards any of the supported formats has been factored into final rounding functions. These final rounding functions receive an intermediate result that has undergone already some rounding with respect to the mathematically exact result.

We shall hence explain the techniques used for intermediate and final rounding first, before detailing the algorithms behind multiplication, addition, FMA, division and square root.

Let us just note that libwidefloat also contains conversion functions from and to integer and the IEEE754 binary32 and binary64 formats. These conversions have been coded in ad-hoc ways, which are not described inhere.

A. Intermediate and Final Rounding

In order to provide correct rounding •(y) = •(f (x)) for a function f , an intermediate result ỹ must be computed in such a way that •(ỹ) = •(y) [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]. For all operations supported by libwidefloat, computing such an approximation ỹ is pretty easy: it is either possible to compute the exactly representable result f (x) (multiplication), or to compute a good approximation together with a flag indicating inexactness. Indeed, for division and square root, the backward error can be exactly computed and compared to zero to yield the inexactness predicate. For addition and FMA, it is just necessary to control the part of the significand that gets discarded during the alignment shift.

However, when computing the approximation ỹ out of the exact result or a better approximation assorted with an inexact indication, pure truncation is not enough. Such a rounding scheme would suffer from the so-called doublerounding issue [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF], [START_REF] Boldo | When double rounding is odd[END_REF].

In order to overcome the double-rounding issue, libwidefloat uses rounding to odd for the intermediate results ỹ. This particular rounding mode is due to Boldo and Melquiond [START_REF] Boldo | When double rounding is odd[END_REF] and can be though of as a intermediate rounding longer than final precision where the last few additional bits represent the classical round, guard and sticky bits [START_REF] Ercegovac | Digital arithmetic[END_REF].

B. Multiplication

As already stated, multiplication in libwidefloat is based on an exact multiplication of the input significands, considered to be integer. At the current state of the library development, naïve schoolbook multiplication using O(n 2 ) machine multiplications is used. The partial products are first represented in a carry-save format and then reduced [START_REF] Defour | Software carry-save for fast multipleprecision algorithms[END_REF].

The use of Karatsuba or Toom-Cook Multiplication [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF] is left for future work but might prove necessary.

C. Addition, Subtraction and Fused-Multiply-Add

When implemented in software, Floating-Point Addition and Subtraction reduce to the same algorithm. The FMA operation can also share this part of code; the only difference is in whether one of the arguments to addition comes directly from one of the operands or from an intermediate exact product.

The core algorithm for addition and subtraction used in libwidefloat is a classical one [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF], [START_REF] Ercegovac | Digital arithmetic[END_REF] and very simple: after reordering of the FP operands for exponent, the significand of the lesser operand is shifted right for alignment, unless this shift would be larger than the output precision. The aligned significands are then added resp. subtracted. A normalization step then ensures correct construction of a FP significand in the case of elimination for subtractions. In software, due to overhead in branching and code bloat, it is sensible to run that normalization step in any case and not to distinguish between a near-path and a far-path.

D. Division and Square Root

The algorithms for division and square root are both based on an approximation to 1/ √ z in libwidefloat. We compute division x/y as x/( √ y) 2 and square root √ x as x × 1/ √

x. An approximation to 1/ √ x can easily be computed using the Newton-Raphson iteration

y n+1 = y n • 1 2 • 3 -x • y 2 n
seeded with a small table for y 0 [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]. In libwidefloat we use a 6 bit approximation for y 0 .

In an extended-precision software environment, it is useful not to perform the Newton iteration on software-simulated FP variables but purely on integer variables, using a fixed-point representation. We use this approach in libwidefloat. However, we chose to generate that code sequence instead of writing it manually in order not to have to determine all alignment shifts in that point by hand.

VI. EXPERIMENTAL RESULTS

Our current implementation of libwidefloat provides 1303 user-available functions. We have tested it for performance against MPFR 3.1.2-p3 based on GMP 6.0.0. We used gcc 4.9.2 on a 64bit Linux system featuring a 4 core Intel i7-5500U running at 2.4 GHz. The optimization level was set to -O3.

There is no one-to-one relation between the functions provided by libwidefloat and those provided by MPFR. We therefore do not give performance results per function. We rather wrote two benchmark applications and compiled them against both libwidefloat and MPFR. The first benchmark application consists of a LU-decomposition, triangular system solving and residue computation. The second one uses secant's method to approximate the zero of a high-degree polynomial.

All benchmark codes were run in all precisions supported by libwidefloat and compared to MPFR, which was run at the same precisions. While we focused on performance testing, we also compared the actual numerical results of libwidefloat and MPFR. On all benchmark codes and all precisions, the results were bitwise identical between both libraries, as expected, given that both libraries provide correct rounding.

We ran our LU-decomposition, triangular system solving and residue benchmark for matrix sizes 1000 × 1000. The matrices and respective right-hand sides were filled with random IEEE754 binary64 values between -107 and 10 matrices was included in the performance testing. Matrix and vector-allocation for MPFR required calling the MPFR initialization procedures for each element of the matrix resp. vector; libwidefloat was not burdened by this requirement.

The results of the LU-decomposition benchmark are given in Figure 1. They show that libwidefloat performs well better than MPFR for small precisions but that starting with precisions around 224 bits, the use of Karatsuba-multiplication helps in keeping MPFR's timings low while libwidefloat exhibits quadratic behavior. We shall address this issue in our library in future releases.

In our second benchmark, we used secant's method to approximate a zero of a polynomial of degree 150 up to the accuracy that was achievable within the given compute precision (which we made vary for the test.) The polynomial was input as an expression tree and evaluated using a recursive algorithm. This example is typical for numerical Computer Algebra Systems, such as Sollya 7 or Sage 8 .

As shown with Figure 2, for that second expressionevaluation benchmark, libwidefloat typically performed better with respect to MPFR. However, our library is still hindered by the use of a naïve multiplication algorithm for precisions higher than 320 bits.

VII. CONCLUSION AND FUTURE WORK

The libwidefloat software fills the gap between IEEE754 binary64 (or binary128 where available) and arbitrary precision libraries, which involve a certain amount of overhead in the mid-precision range. Our library leverages modern compilers' optimization techniques, such as constant propagation, code inlining and dead code elimination, to achieve this goal.

At the current state, libwidefloat supports all basic operations, like e.g. addition, FMA, square root or comparisons. Development is ongoing: we are working on a C++ wrapper to allow for an easy drop-in-replacement of the IEEE754 binary64 format by libwidefloat formats. Further, we plan to use GNU's printf/scanf extension techniques to allow for easy input/output operations with libwidefloat. The development of the elementary intrinsic functions like the exponential, the logarithm or the trigonometric operations is also part of future work.

Experimentally, we see that libwidefloat's performance still needs some improvement in the range between 256 and 512 bits. The current implementation of the library lacks sufficient performance due to the use of a naïve (schoolbook) multiplication algorithm. We shall extend the library to use the Karatsuba or Toom-Cook multiplication algorithms in that higher precision range.

7

 7

Figure 2 .

 2 Figure 1. Experimental results: Linear Algebra

http://www.mpfr.org/

http://www.davidhbailey.com/dhbsoftware/

see https://gmplib.org/

see e.g. https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html

see e.g. https://gcc.gnu.org/onlinedocs/gcc/ 005f 005fint128.html

see e.g. https://gcc.gnu.org/onlinedocs/gcc/Restricted-Pointers.html

see http://sollya.gforge.inria.fr/

see http://www.sagemath.org/