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

Valentina Popescu
email: valentina.popescu@ens-lyon.fr

Ping Tak

Peter Tang
email: peter.tang@intel.com

J.-M Muller

A new multiplication algorithm for extended precision using floating-point expansions

Keywords: floating-point arithmetic, floating-point expansions, high precision arithmetic, multiple-precision arithmetic, multiplication !

Some important computational problems must use a floating-point (FP) precision several times higher than the hardwareimplemented available one. These computations critically rely on software libraries for high-precision FP arithmetic. The representation of a high-precision data type crucially influences the corresponding arithmetic algorithms. Recent work showed that algorithms for FP expansions, that is, a representation based on unevaluated sum of standard FP types, benefit from various high-performance support for native FP, such as low latency, high throughput, vectorization, threading, etc. Bailey's QD library and its corresponding Graphics Processing Unit (GPU) version, GQD, are such examples. Despite using native FP arithmetic as the key operations, QD and GQD algorithms are focused on double-double or quad-double representations and do not generalize efficiently or naturally to a flexible number of components in the FP expansion. In this paper, we introduce a new multiplication algorithm for FP expansion with flexible precision, up to the order of tens of FP elements in mind. The main feature consists in the partial products being accumulated in a special designed data structure that has the regularity of a fixed-point representation while allowing the computation to be naturally carried out using native FP types. This allows us to easily avoid unnecessary computation and to present rigorous accuracy analysis transparently. The algorithm, its correctness and accuracy proofs and some performance comparisons with existing libraries are all contributions of this paper.

INTRODUCTION

M ANY numerical problems in dynamical systems or planetary orbit dynamics require higher precisions than the standard double-precision (now called binary64 [START_REF]IEEE Computer Society, IEEE Standard for Floating-Point Arithmetic. IEEE Standard[END_REF]). Examples include, to list a few, long-term stability analysis of the solar system [START_REF] Laskar | Existence of collisional trajectories of Mercury, Mars and Venus with the Earth[END_REF], finding sinks in the Henon Map [START_REF] Joldes | Searching for sinks for the hénon map using a multipleprecision gpu arithmetic library[END_REF] and iterating the Lorenz attractor [START_REF] Abad | Computing periodic orbits with arbitrary precision[END_REF]. Since quadruple or higher precision is rarely implemented in hardware, these calculations rely on software-emulated higher-precision libraries, also called arbitrary-precision libraries.

A crucial design point of these libraries is the way higher-precision numbers are represented as this will influence the choice of arithmetic algorithms and their subsequent performance on various computer architectures, including highly parallel accelerators such as Graphics Processing Units (GPUs). There are mainly two ways of representing numbers in higher precision. The first one is the multiple-digit representation: each number is represented by one exponent, followed by a sequence of possibly high-radix digits. This representation gives a fixed-point flavor as the radix positions of the digits follow a very regular pattern. This characteristic greatly facilitates rigorous accuracy analysis on the associated arithmetic algorithms. The GNU package MPFR [START_REF] Fousse | MPFR: A Multiple-Precision Binary Floating-Point Library with Correct Rounding[END_REF] is an open-source C library that uses the multiple-digit representation. Besides arbitrary precision, it also provides correct rounding for each atomic operation. One drawback, however, is that full exploitation of highlyoptimized native floating-point (FP) instructions can be hard to achieve. The more memory intensive nature of the algorithms render them not a natural fit for GPUs, to which MPFR has yet to be ported.

The second way is the multiple-term representation in which each number is expressed as the unevaluated sum of several standard FP numbers. This sum is usually called a FP expansion. Because each term of this expansion has its own exponent, the term's positions can be quite irregular. Bailey's library QD [START_REF] Hida | Algorithms for quad-double precision floating-point arithmetic[END_REF] uses this approach and supports double-double and quad-double computations, i.e. numbers are represented as the unevaluated sum of 2 or 4 standard double-precision FP numbers. The FP centric nature allows QD to take great advantage of optimized native FP infrastructure; and the GPU version, GQD, is already available. Nevertheless, there are some drawbacks of FP expansions as implemented in QD. These algorithms do not straightforwardly generalize to an arbitrary number of terms, for example, a multiplication of two n-term FP expansion can easily lead to O(n 3 ) complexity. Moreover, their nature made rigorous error analysis unnatural, that is why the operations implemented in QD do not come with any form of guaranteed error bounds.

In this article we focus on multiplication of FP expansions. Our algorithm accumulates the partial products of a multiplication using a fixed-point structure that is nevertheless FP friendly. The overall algorithm can be viewed as a generalized version of the paper-andpencil method for long multiplication. The accumulation is done using a method by Rump, first presented in [START_REF] Rump | Ultimately Fast Accurate Summation[END_REF]. The regularity of this algorithm allows us to provide a thorough error analysis and tight error bound, while its FP friendliness make it a natural fit for GPUs.

The outline of the paper is the following: after recalling some basic notions about FP expansions in Section 2, we detail the new algorithm for multiplying FP expansions in Section 3. Here we also include a correctness proof, the computation of the specific error bound (Section 3.1) and the complexity analysis (Section 3.2). Finally, in Section 4 we assess the performance of our algorithm in terms of performance comparing to other existing libraries and we draw the conclusions.

FLOATING-POINT EXPANSIONS

A normal binary precision-p floating-point (FP) number is a number of the form

x = M x • 2 ex-p+1 , with 2 p-1 ≤ |M x | ≤ 2 p -1.
The integer e x is called the exponent of x, and M x •2 -p+1 is called the significand of x. We denote accordingly to Goldberg's definition ulp

(x) = 2 ex-p+1 [8, Chap. 2].
A natural extension of the notion of double-double or quad-double is the notion of floating-point expansion.

Definition 2.1.

A FP expansion u with n terms is the unevaluated sum of n FP numbers u 0 , u 1 , . . . , u n-1 , in which all nonzero terms are ordered by magnitude (i.e., if v is the sequence obtained by removing all zeros in the sequence u, and if sequence v contains m terms,

|v i | ≥ |v i+1 |, for all 0 ≤ i < m -1).
Arithmetics on FP expansions have been introduced by Priest [START_REF] Priest | Algorithms for arbitrary precision floating point arithmetic[END_REF], and later on by Shewchuk [START_REF] Shewchuk | Adaptive precision floating-point arithmetic and fast robust geometric predicates[END_REF].

To make sure that such an expansion carries significantly more information than only one FP number, it is required that the u i 's do not "overlap". The notion of (non-)overlapping varies depending on the authors. We give here different definitions, using the above introduced notation. The first two, by Priest and Bailey, have already been presented in the literature, but the third one introduces a new, weaker condition, that allows for a more relaxed handling of the FP expansions.

We specify first that even if a FP expansion may contain interleaving zeros, all the definitions that follow apply only to the non-zero terms of the expansion (i.e., the array v in Definition 2.1). Definition 2.2. Assuming x and y are normal numbers with representations M x • 2 ex-p+1 and M y • 2 ey-p+1 (with 2 p-1 ≤ |M x |, |M y | ≤ 2 p -1), they are P-nonoverlapping (that is, nonoverlapping according to Priest's definition [START_REF] Priest | On properties of floating-point arithmetics: Numerical stability and the cost of accurate computations[END_REF]) if |e y -e x | ≥ p. Definition 2.3. An expansion is P-nonoverlapping (that is, nonoverlapping according to Priest's definition [START_REF] Priest | On properties of floating-point arithmetics: Numerical stability and the cost of accurate computations[END_REF]) if all its components are mutually P-nonoverlapping.

A slightly stronger sense of nonoverlapping was introduced by Hida, Li and Bailey [START_REF] Hida | Algorithms for quad-double precision floating-point arithmetic[END_REF]: Definition 2.4. An expansion u 0 , u 1 , . . . , u n-1 is Bnonoverlapping (that is, nonoverlapping according to Bailey's definition [START_REF] Hida | Algorithms for quad-double precision floating-point arithmetic[END_REF]) if for all 0 < i < n, we have |u i | ≤ 1 2 ulp(u i-1 ). Remark 2.5. Note that for P-nonoverlapping expansions we have |u i | ≤ 2 p -1 2 p ulp(u i-1 ). Intuitively, the stronger the sense of the nonoverlapping definition, the more difficult it is to obtain, implying extra manipulation of the FP expansions. In order to save operations we chose to use a slightly weaker sense of nonoverlapping, referred to as ulp-nonoverlapping , that we define in what follows.

Definition 2.6. An expansion u 0 , u 1 , . . . , u n-1 is ulpnonoverlapping if for all 0 < i < n, we have |u i | ≤ ulp(u i-1 ).

In other words, the components are either Pnonoverlapping or they overlap by one bit, in which case the second component is a power of two.

Depending on the nonoverlapping type of an expansion, when using standard FP formats for representation, the exponent range forces a constraint on the number of terms. The largest expansion can be obtained when the largest term is close to overflow and the smallest is close to underflow. We remark that, when using any of the above nonoverlapping definitions, for the two most common FP formats, the constraints are:

• for double-precision (exponent range [-1022, 1023]) the maximum expansion size is 39; • for single-precision (exponent range [-126, 127]) the maximum is 12. The majority of algorithms performing arithmetic operations on FP expansions are based on the socalled error-free transforms (such as the algorithms 2Sum, Fast2Sum, Dekker's product and 2MultFMA presented for instance in [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]), that make it possible to compute both the result and the error of a FP addition or multiplication. This implies that each such error-free transform applied to two FP numbers, returns still two FP numbers. So, when adding two expansions x and y with n and m terms, respectively, the exact result is going to have at most n + m terms. Similarly, for multiplication, the exact product is going to have at most 2nm terms [START_REF] Priest | Algorithms for arbitrary precision floating point arithmetic[END_REF]. A potential problem appears when subsequent computations are done using this results; the size of the exact result is going to increase more and more. To avoid this, some "truncation" methods (both on-the-fly or a-posteriori) may be used to compute only an approximation of the exact result. Also, so-called (re-)normalization algorithms are used to render the result nonoverlapping, which implies also a potential reduction in the number of components.

In what follows we will present a new multiplication algorithm that allows for computing the "truncated" product of two FP expansions. We also provide a full correctness proof and error analysis for it.

MULTIPLICATION

ALGORITHM FOR FLOATING-POINT EXPANSIONS

In Algorithm 1 we present the multiplication algorithm. Throughout this paper we assume that no underflow / overflow occurs during the calculations. We consider two ulp-nonoverlapping expansions x and y, with n and m terms, respectively, and we compute the r most significant components of the expansion π = π 0 , π 1 , π 2 , . . . that represents the product x • y.

Algorithm 1 Algorithm for multiplication of FP expansions. We denote by p, the precision of the used FP format and by b, the size of the bins. The algorithms called are going to be detailed later on in this paper.

Input: ulp-nonoverlapping FP expansions x = x 0 + . . . + x n-1 ; y = y 0 + . . . + y m-1 . Output: ulp-nonoverlapping FP expansion π = π 0 + . . . + π r-1 . 1: e ← e x0 + e y0 2: for i ← 0 to r • p/b + 1 do 3: B i ← 1.5 • 2 e-(i+1
B i ← B i -1.5 • 2 e-(i+1)b+p-1 23: end for 24: π[0 : r -1] ← Renormalize(B[0 : r • p/b + 1]) 25: return FP expansion π = π 0 + . . . + π r-1 .
The way we compute the partial products is based on term-times-expansion products, x i • y, and it follows the paper-and-pencil technique. In order to gain performance we "truncate" the computations by discarding the partial products that have an order of magnitude less than π r , meaning that we compute the approximate result based on the first r+1 k=1 k partial products. The products with the same order of magnitude as π r are intended as an extra correction term, that is why we compute them using only standard FP multiplication.

The main idea of the algorithm is to accumulate numbers of size at most b in "containers", referred to as bins, that are FP variables whose least significant bit (LSB) has a fixed weight. This allows for errorless accumulation provided that we do not add more than 2 c numbers to one bin, where b + c = p -1.

Even though in this article we use a general notation, our implementation is done for standard FP formats. When using double-precision (p = 53) we defined bins of size b = 45, that allows for c = 7 bits of carry to happen; this means that we can add 128 numbers that satisfy the above condition to each bin and the result is still going to be exact. For single-precision (p = 24) we chose bins with b = 18, which implies c = 5, allowing us to add up to 32 numbers to one bin. In both cases, these values also satisfy 3b > 2p, a property that we are going to use later on. We suggest the use of the latter one only if is the only standard FP format available or if, by any reason, is much faster than the first one.

The number of allocated bins is computed as r•p b + 2 and the LSB of each bin is set according to the starting exponent, e = e x0 + e y0 , at a distance of b bits. We start the algorithm by initializing each bin B i with the value

1.5 • 2 e-(i+1)•b+p-1 .
After the initialization step is done we start the actual computations. For each partial product computed using 2MultFMA we get the pair (P, E) (line 7 of Algorithm 1) and, using the formula (e -e xi -e yj )/b , we determine the corresponding bins in which we have to accumulate it. We know that one pair of FP numbers can "fall" into at most four bins (3b > 2p), and we can deduce three different cases (see Fig. 1):

• 2:2 case, in which both P and E fall into two bins each; • 2:3 case, in which P falls into two bins and E into three;

• 3:2 case, in which P falls into three bins and E into two. These cases are dealt with in Algorithm 2 that accumulates the partial products. Apart from the (P, E) pair and the bins array, B, the algorithm also receives two integer parameters. The sh value represents the first corresponding bin for the pair and , computed as e-e xi -e yj -sh•b, the number of leading bits. This value gives the difference between the LSB of B sh-1 and the sum of the exponents of x i and y j of the corresponding (P, E) pair.

We determine which one of the three cases apply depending on the value. So, if 2c + 1 < b -≤ b we are in the 2:2 case; if c < b -≤ 2c + 1 we need to consider the 2:3 case; and if 0 < b -≤ c the 3:2 case applies. Remark 3.1. For simplicity, we consider that the extra error correction partial products, the ones that are computed using standard FP multiplication, are dealt with the same way, using the pair (P, 0). This is not the case in our implementation, where we try to save operations by accumulating only the P term.

Algorithm 2 Accumulate(P, E, B, sh, ).

Input: FP numbers P, E; FP array B;

Integers sh and . Output:

FP array B. 1: if < b -2c -1 then 2:
(B sh , B sh+1 ) ← Deposit(P )

3:

// that is, (B sh , P ) ← Fast2Sum(B sh , P ), and As we stated before, all the bins are initialized with a constant value depending on their LSB, that is going to be subtracted before the renormalization step. This type of addition was first used by Rump (in [START_REF] Rump | Ultimately Fast Accurate Summation[END_REF]) for adding the elements of an array of FP numbers. In his paper he proved that the result is correct. For the sake of completeness we also give here a short correctness proof.

Proof of correctness for Algorithm 2

In what follows we will prove that when accumulating the numbers as we do in Algorithm 2 no rounding errors can occur and the result is exact.

We consider all the values that fall into the same bins B sh , B sh+1 as an array of FP numbers x 1 , . . . , x n that satisfy |x i | < 2 e+b , where 2 e is the LSB of B sh and b is the size of the bins. The lower part of each x i is denoted by x l i and represents the part that will be accumulated into B sh+1 . Proposition 3.2. Let x 1 , . . . , x n an array of precision-p FP numbers that satisfy |x i | < 2 e+b , for all 0 < i ≤ n, where b < p-1. We initialize a FP container with s 0 = 1.5•2 e+p-1 , we compute s 1 = RN(s 0 + x 1 ); . . . ; s i = RN(s i-1 + x i ); . . . ; s n = RN(s n-1 + x n ); and we return the value RN(s n -s 0 ). For each x i we also compute the lower part x i = RN(RN(s i-1 -s i ) + x i ). When using this method, no significant informations is lost in the process, though no rounding can occur, provided that n ≤ 2 p-b-2 -1 and that no underflow / overflow occurs during the calculations.

Proof: We first prove by induction that the following statement holds:

1.5 • 2 e+p-1 -i • 2 e+b ≤ s i ≤ 1.5 • 2 e+p-1 + i • 2 e+b . (1)
It is easy to see that it holds for i = 0. Now we assume that is true for i and we try to prove it for i + 1. We deduce:

1.5•2 e+p-1 -(i+1)2 e+b ≤ s i +x i+1 ≤ 1.5•2 e+p-1 +(i+1)2 e+b .
Hence, since rounding is a monotonic function:

RN(1.5 • 2 e+p-1 -(i + 1)2 e+b ) ≤ s i+1 ≤ RN(1.5 • 2 e+p-1 + (i + 1)2 e+b ).
The value 1.5 • 2 e+p-1 -(i + 1)2 e+b is an exact FP number because it is a multiple of 2 e and it is less than 2 e+p in absolute value, provided that i ≤ 2 p-b-1 • 3.5 -1, which holds in all practical cases (with our parameters: i ≤ 447 for double-precision -p = 53 and b = 45and i ≤ 111 for single-precision -p = 24 and b = 18-). The same holds for 1.5 • 2 e+p-1 + (i + 1)2 e+b provided that i ≤ 2 p-b-2 -1, which also holds in all practical cases (with our parameters: i ≤ 63 for double-precision and i ≤ 15 for single-precision).

Furthermore, we have

(s i , x i ) = Fast2Sum(s i-1 , x i ), therefore ∀i, s i + x i = s i-1 + x i ,
such that, by induction,

s i + x i + x i-1 + . . . + x 1 = s 0 + x 1 + x 2 + . . . + x i , which implies (s n -s 0 ) + x n + x n-1 + . . . + x 1 = x 1 + x 2 + . . . + x n .
From (1), we easily find that s n and s 0 are within a factor 2 (in practice, much less), such that (from Sterbenz' lemma), their difference is exactly computed: RN(s n -s 0 ) = s n -s 0 . We therefore conclude that

RN(s n -s 0 ) + x n + x n-1 + . . . + x 1 = x 1 + x 2 + . . . + x n .
The last part of Algorithm 1 consists in applying a renormalization algorithm (see Algorithm 3) on the bins array, in order to render the result ulp-nonoverlapping . This is a chain of Fast2Sum starting from the two most significant components and propagating the errors. If however, the error after a Fast2Sum block is zero, then we propagate the sum (this is shown in Figure 2).

Algorithm 3 Renormalize(x).

Input: FP array x = (x 0 , . . . , x n-1 ); Output: ulp-nonoverlapping FP expansion r = r 0 + . . . + r k-1 . 1: ε ← x 0 2: j ← 0 3: i ← 1 4: while i < n and j < k do 5: i ← i + 1 12: end while 13: if ε = 0 and j < k then 14: r j ← ε 15: end if 16: return FP expansion r = r 0 + . . . + r k-1 . Fig. 2: Renormalize with n terms. Each Fast2Sum box performs an error-free transform; the sum is outputted downwards and the error to the right. If the error is zero, the sum is propagated to the right, otherwise the error is propagated and the sum is outputted.

(r j , ε) ← Fast2Sum(ε, x i ) 6: if ε = 0 then 7: ε ← r j 8:

Proof of correctness for Algorithm 3

In what follows we will prove the correctness of Algorithm 3 by using the notations seen in Figure 2.

Proposition 3.3. Let an input array x = (x 0 , . . . , x n-1 ) that satisfies:

x i is multiple of 2 e0-ib = 2 ei and |x i+1 | < 2 ei+c+1 , for all 0 ≤ i < n -1, where b + c = p -1 (c << b). Also |x 0 | > |x 1 | and |x 0 | < 2 e0+p+1
. After applying Algorithm 3, the output array r = (r 0 , . . . , r k-1 ), with 0 ≤ k ≤ n-1 is an ulp-nonoverlapping expansion, provided that no underflow / overflow occurs during the calculations.

Proof: The case when x contains 1 or 2 elements is trivial. Consider now at least 3 elements. By the input type we know that:

x 0 = X 0 • 2 e0 , x 1 = X 1 • 2 e1 with e 1 = e 0 -b.
Hence, r 0 and ε 1 are both multiples of 2 e0-b . Two possible cases may occur:

(i) ε 1 = 0. If we choose to propagate directly ε 1 , then r 1 = x 2 and ε 2 = 0. This implies, by induction, that r i = x i+1 , ∀i ≥ 1. So, directly propagating the error poses a problem, since the whole remaining chain of Fast2Sum is executed without any change between the input and the output. So, as shown in Algorithm 3, line 7, when ε i = 0 we propagate the sum r j , so ε i ← r j .

(ii) ε 1 = 0. By definition of Fast2Sum, we have

|ε 1 | ≤ 1 2 ulp(r 0 ). We also have |x 0 | > |x 1 |, so |r 0 | = |RN(x 0 + x 1 )| ≤ 2 |x 0 |. Hence: |ε 1 | < 2 e0+2 .
We now prove by induction the following statement: at each step i ≥ 1 of the loop in Algorithm 3, both r i-1 and ε i are multiples of 2 ei and ε i = 0 or |ε i | < 2 ei+b+2 , meaning that ε i fits in at most b+1 bits. We proved above that for i = 1 this holds. Suppose now it holds for i and prove it for i + 1.

At this step we have

ε i + x i+1 = r i + ε i+1 . Since ε i is a multiple of 2 ei and x i+1 is a multiple of 2 ei+1
, with e i+1 = e i -b, then both r i and ε i+1 are multiples of 2 ei+1 . Two cases may occur: -if |ε i | < 2 ei+c+1 then ε i + x i+1 is a FP number, which implies r i = ε i +x i+1 exactly and ε i+1 = 0, in which case we propagate r i < 2 ei+c+2 . -if |ε i | > 2 ei+c we have |r i | ≤ 2 |ε i | and we get (by definition of Fast2Sum):

|ε i+1 | ≤ 1 2 ulp(r i ) < 2 -p • 2 • 2 ei+b+2 < 2 ei-c+2 (2)
This condition is even stronger than what we were trying to prove, so the induction holds. Finally, we prove the relation between r i-1 and r i . If

ε i = 0, we propagate r i-1 , i.e. ε i ← r i-1 . Otherwise |r i | = |RN(ε i + x i+1 )| ≤ 2 |ε i | and since ε i ≤ 1 2 ulp(r i-1 ), then |r i | ≤ ulp(r i-1
) and the proposition is proven. Remark 3.4. After using Algorithm 3, the result π = π 0 + . . . + π r-1 (of Algorithm 1) cannot have interleaving zeros; zeros may appear only at the end of the expansion.

Error analysis for multiplication

The multiplication algorithm explained in the above section returns only an approximation of the exact result, but the maximum error that can occur is bounded. We give here the error bound that we computed and we prove its correctness. Proposition 3.5. Let x and y be two ulp-nonoverlapping FP expansions, with n, and m terms, respectively. If, for computing the product xy we use Algorithm 1, the result π, with r terms, satisfies:

|xy -π| ≤ |x 0 y 0 | 2 -(p-1)r [1 + (r + 1)2 -p + 2 -(p-1) -2 -(p-1) (1-2 -(p-1) ) 2 + m+n-r-2 1-2 -(p-1)
], provided that no underflow / overflow occurs during the calculations.

For simplicity, we will first prove an intermediate proposition, and only after that we proceed to the proof on the final bound.

x i y j ≤ |x 0 y 0 | 2 -(p-1)(r+1) -2 -(p-1) (1-2 -(p-1) ) 2 + m+n-r-2 1-2 -(p-1)
, provided that no underflow / overflow occurs during the calculations.

Proof: From the definition of the ulp-nonoverlapping expansion we have x 1 ≤ ulp(x 0 ) ≤ 2 -p+1 |x 0 | and, by induction, we get x i ≤ 2 -pi+i |x 0 |, for all 0 < i < n. The same goes for y.

The discarded partial products satisfy:

m+n-2 k=r+1 i+j=k a i b j ≤ m+n-2 k=r+1 i+j=k 2 -p(i+j)+i+j |x 0 y 0 | ≤ |x 0 y 0 | m+n-2 k=r+1 i+j=k 2 -(p-1)k ≤ |x 0 y 0 | m+n-2 k=r+1 (m + n -1 -k)2 -(p-1)k ≤ |x 0 y 0 | m+n-r-3 k =0 (m + n -k -r -2)2 -(p-1)(k +r+1) ≤ |x 0 y 0 | 2 -(p-1)(r+1) × m+n-r-3 k =0 (m + n -r -2 -k )2 -(p-1)k . (3)
We define the function

φ(e) = ∞ k=0 (m + n -r -2 -k)e k
that satisfies:

φ(e) = ∞ k=0 -ke k + ∞ k=0 (m + n -r -2)e k = -e ∞ k=1 ke k-1 + (m + n -r -2) ∞ k=0 e k = -e d de ∞ k=1 e k + (m + n -r -2) 1 1 -e = -e d de 1 1 -e + m + n -r -2 1 -e = -e (1 -e) 2 + m + n -r -2 1 -e .
When applying function φ(2 -(p-1) ) in equation ( 3) we get: 1) , which concludes the proof.

m+n-2 k=r+1 i+j=k x i y j ≤ |x 0 y 0 | 2 -(p-1)(r+1) • • -2 -(p-1) (1 -2 -(p-1) ) 2 + m + n -r -2 1 -2 -(p-
We are now able to prove Prop. 3.5.

Proof: (of Prop. 3.5) When using Algorithm 1 we "truncate" the result by discarding the partial products with an order of magnitude less than π r . From Prop. 3.6 we know that this causes a maximum error that is less or equal to 1) .

|x 0 y 0 | 2 -(p-1)(r+1) -2 -(p-1) (1 -2 -(p-1) ) 2 + m + n -r -2 1 -2 -(p-
(4) Also, in the algorithm we do not use error-free transforms for computing the last r + 1 partial products, the ones with the same order of magnitude as π r , for which i + j = r. We know that |x i y j | ≤ 2 -(p-1)(i+j) |x 0 y 0 |, from where |x i y j -RN(x i y j )| ≤ 2 -(p-1)r |x 0 y 0 | 2 -p . This implies that the maximum error caused doing this is less or equal to:

(r + 1) |x 0 y 0 | 2 -(p-1)r • 2 -p . (5)
Apart from these two possible errors we also need to account for the error caused by the renormalization step. We already showed that Algorithm 3 returns an ulpnonoverlapping expansion, in which case the maximum error is less or equal to ulp(π r-1 ). This implies that is less or equal to:

|x 0 y 0 | 2 -(p-1)r . (6)
To get the final error bound we have to add the bounds on all the possible errors that can occur, Eq. ( 4), ( 5) and ( 6), and we get:

|x 0 y 0 |2 -(p-1)r [1 + (r + 1)2 -p + + 2 -(p-1) -2 -(p-1) (1 -2 -(p-1) ) 2 + m + n -r -2 1 -2 -(p-1) ]
This concludes our proof.

Complexity analysis for multiplication

As presented before our algorithm has the advantage of computing "truncated" expansions on the fly, by taking into account only the significant terms of the input expansions.

We present here the worst case FP operation count of our algorithm, by taking ( [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]) 3 FP operations for Fast2Sum and 2 for 2MultFMA. For the sake of simplicity, we will consider that the input expansions and the result have the same number of terms, say k. Proposition 3.7. When using Algorithm 1 for computing the product of two FP expansions with k terms, we perform Proof: During a preprocessing step we need to get the exponents of each term of the input expansions. We do this using the math.h library function, frexp, that uses only one FP operation, which we call 2k times.

The first step of the algorithm consists in initializing the bins; as mentioned in Sec. 3, we allocate binN r = k•p b + 2 bins. For this we use twice the math.h library function, ldexp (that can take up to 34 FLOPS, depending on the exponent's size), and after that we perform binN r -1 FP multiplications.

During the main loop of the algorithm we compute k i=1 i partial products using 2MultFMA, which we accumulate into the bins using 3 Fast2Sum calls and 2 FP additions. During this same loop we also compute k -1 correction terms using simple FP multiplication, which we can accumulate using only 2 Fast2Sum calls and one FP addition.

In the last part of the algorithm we first unbias the bins, by subtracting their initial values from each one of them, followed by the renormalization of the result. This accounts for a total of binN r FP subtractions followed by binN r -1 calls to Fast2Sum and comparisons.

When adding all the operations accounted for above, we conclude that Algorithm 1 requires a total of 13 2 k 2 +

COMPARISON AND DISCUSSION

In We chose to use it for comparison because it is more general than the QD/GQD one, which is limited to double-double and quad-double precision, plus, from our knowledge, is the only algorithm in the literature of FP expansions provided with a complete correctness proof and error analysis. This algorithm is not tuned for obtaining "truncated" results on the fly, meaning that it computes the result fully and only then truncates, which allows for a straightforward error analysis, but makes it more costly. In practice however, Baileys QD/GQD and MPFR libraries are very often used, so we present the performance results comparing to these. We do this for completeness reasons, even though the comparison is not completely fair since we target different things. To be exact: (i) MPFR uses a different representation for extended precision and is limited to CPU use; (ii) QD and GQD libraries are limited to double-double and quad-double precisions, plus the algorithm used in the implementation of quad-double are not provided with any correctness proof and there is no specific error bound given.

This algorithm is part of the CAMPARY (CudA Multiple Precision ARithmetic librarY) software that is available at: http://homepages.laas.fr/mmjoldes/campary/. Our implementation was done using CUDA C -an extension of the C language developed by NVIDIA [START_REF][END_REF] for their GPUs. Algorithm 1 is suitable for GPU use, since all basic operations (+, -, * , /, √ ) conform to the IEEE 754-2008 standard for FP arithmetic for singleand double-precision, support for the four rounding modes is provided and dynamic rounding mode change is supported without any penalties. The fma instruction is supported for all devices with Compute Capability at least 2.0, which allows us to use the 2MultFMA algorithm.

In the implementation we use templates for the number of terms in the expansions, meaning that we allow static generation of any input-output precision combinations (e.g. add a double-double with a quad-double and store the result on triple-double). All the functions are defined using host device specifiers, which allows for the library to be used on both CPU and GPU. We would like to stress here the fact that the algorithm is not itself parallelized; our intention is to provide a multiple precision alternative for problems that are suitable for GPU use.

In Table 2 we give some CPU performance measurements obtained on an Intel(R) Xeon(R) CPU E5-2609 v3 @ 1.90GHz computer. The values are given in MFlops/s and by one Flop we understand one operation using extended precision. Table 3 contains the performance values obtained on a GPU, in comparison to the GQD library. The tests were performed on a Tesla K40c GPU, using CUDA 7.5 software architecture, running on a single thread of execution. From the above mentioned tables we deduce that for small size expansions (2 up to 4 terms) the algorithm is too complex to be efficient. In these cases we suggest the use of a different algorithm, presented shortly in [START_REF] Joldes | Arithmetic algorithms for extended precision using floating-point expansions[END_REF] (Algorithm 5 and Figure 4), that takes advantage of the processor's pipeline by avoiding branches. That algorithm is a generalized version of the QD's multiplication combined with a proven renormalization algorithm. It allows for a tight error bound when small expansion results are needed (just slightly larger than the one obtained when using Algorithm 1), but it becomes exponentially worst as the size of the result increases, which limits its practical interest to expansions with at most 4 terms.

)b+p- 1 4: end for 5 : 6 :PB

 156 for i ← 0 to min(n -1, r) do for j ← 0 to min(m -1, r -1 -i) do 7: (P, E) ← 2MultFMA(x i , y j ) 8: ← e -e xi -e yj 9: sh ← /b 10: ← -sh • b 11: B ← Accumulate(P, E, B, sh, ) ← x i • y j 15: ← e -e xi -e yj ← Accumulate(P, 0., B, sh, ) 19: end if 20: end for 21: for i ← 0 to r • p/b + 1 do 22:

Fig. 1 :

 1 Fig. 1: Cases for accumulating the partial products into the bins.

4 :// B sh+1 ← B sh+1 + P 5 :

 45 (B sh+1 , B sh+2 ) ← Deposit(E) 6: else if < b -c then 7: (B sh , B sh+1 ) ← Deposit(P ) 8: (B sh+1 , E) ← Fast2Sum(B sh+1 , E) 9:(B sh+2 , B sh+3 ) ← Deposit(E) 10: else 11:(B sh , P ) ← Fast2Sum(B sh , P )12:(B sh+1 , B sh+2 ) ← Deposit(P ) 13: (B sh+2 , B sh+3 ) ← Deposit(E) 14: end if 15: return FP array B.

Proposition 3 . 6 .

 36 Let x and y be two ulp-nonoverlapping FP expansions, with n, and m terms, respectively. If, when computing the product xy we "truncate" the operations by adding only the first r+1 k=1 k partial productsproducts, where r is the required size of the final result, then the generated error satisfies: xy -r k=0 i+j=k

33 2 k

 2 + 6( k•p b + 2) + 55 FP operations.

Table 1

 1

we give effective values of the worst case FP operation count for our algorithm (see Sec. 3.2) vs. Priest's multiplication algorithm (

[START_REF] Priest | Algorithms for arbitrary precision floating point arithmetic[END_REF]

), which, when multiplying two P-nonoverlapping expansions with n and respectively m terms, it requires 81 mn 2 + 747 nm + 2 m -233 n FP operations.

TABLE 1 :

 1 Effective values of the worst case FP operation count for Algorithm 1 and Priest's multiplication algorithm[START_REF] Priest | Algorithms for arbitrary precision floating point arithmetic[END_REF] when the input and output expansions are of size r.

	r	2	4	8	16
	Algorithm 1	138	261	669	2103
	Priest's mul. ([9])	3174 16212 87432 519312

TABLE 2 :

 2 Performance in MFlops/s for multiplying two FP expansions on the CPU; d x and d y represent the number of terms in the input expansions and d r is the size of the computed result. * precision not supported

	dx, dy, dr	Algorithm 1	QD	MPFR
	2, 2, 2	11.69	99.16	18.64
	1, 2, 2	14.96	104.17	19.85
	3, 3, 3	6.97	*	12.1
	2, 3, 3	8.62	*	13.69
	4, 4, 4	4.5	5.87	10.64
	1, 4, 4	8.88	15.11	14.1
	2, 4, 4	6.38	9.49	13.44
	8, 8, 8	1.5	*	6.8
	4, 8, 8	2.04	*	9.15
	16, 16, 16	0.42	*	2.55

TABLE 3 :

 3 Performance in MFlops/s for multiplying two FP exansions on the GPU; d x and d y represent the number of terms in the input expansions and d r is the size of the computed result. * precision not suported

	dx, dy, dr	Algorithm 1	QD
	2, 2, 2	0.027	0.1043
	1, 2, 2	0.365	0.1071
	3, 3, 3	0.0149	*
	2, 3, 3	0.0186	*
	4, 4, 4	0.0103	0.0174
	1, 4, 4	0.0215	0.0281
	2, 4, 4	0.0142	*
	8, 8, 8	0.0034	*
	4, 8, 8	0.0048	*
	16, 16, 16	0.001	*

ACKNOWLEDGMENTS

The Authors would like to thank Région Rh ône-Alpes and ANR FastRelax Project for the grants that supports this activity.