Claude-Pierre Jeannerod

Hervé Knochel

Christophe Monat

Guillaume Revy

Gilles Villard

A new binary floating-point division algorithm and its software implementation on the ST231 processor

Keywords: binary floating-point division, correct rounding, polynomial evaluation, code generation and validation, VLIW integer processor

This paper deals with the design and implementation of low latency software for binary floating-point division with correct rounding to nearest. The approach we present here targets a VLIW integer processor of the ST200 family, and is based on fast and accurate programs for evaluating some particular bivariate polynomials. We start by giving approximation and evaluation error conditions that are sufficient to ensure correct rounding. Then we describe the heuristics used to generate such evaluation programs, as well as those used to automatically validate their accuracy. Finally, we propose, for the binary32 format, a complete C implementation of the resulting division algorithm. With the ST200 compiler and compared to previous implementations, the speed-up observed with our approach is by a factor of almost 1.8.

Introduction

Although floating-point divisions are less frequent in applications than other basic arithmetic operations, reducing their latency is often an issue [START_REF] Oberman | Design issues in division and other floating-point operations[END_REF]. Since low latency implementations may typically be obtained by expressing and exploiting instruction parallelism, intrinsically parallel algorithms tend to be favored. Some examples are [START_REF] Piñeiro | High-speed doubleprecision computation of reciprocal, division, square root and inverse square root[END_REF][START_REF] Piso | Analysis of the impact of different methods for division/square root computation in the performance of a superscalar microprocessor[END_REF] for hardware, and [START_REF] Markstein | Software division and square root using Goldschmidt's algorithms[END_REF][START_REF] Raina | FLIP: a Floating-point Library for Integer Processors[END_REF] for software.

In this paper we focus on designing and implementing low latency floating-point division software for STMicroelectronics' ST231 processor. For this 32-bit VLIW integer processor, various correctly-rounded binary floating-point division algorithms have already been implemented [START_REF] Raina | FLIP: a Floating-point Library for Integer Processors[END_REF]. There, the lowest latency measured is of 48 cycles; it was obtained by exposing instruction-level parallelism (ILP) by means of a suitable combination of univariate polynomial evaluation and Goldschmidt's method, in a way similar to [START_REF] Piñeiro | High-speed doubleprecision computation of reciprocal, division, square root and inverse square root[END_REF]. However, we will see in this paper that much more ILP can in fact be exposed by relying exclusively on polynomial evaluation, leading to a latency of 27 cycles and thus a speed-up by a factor of almost 1.8. To get this result we shall take three main steps, which we summarize now.

As a first step, we extend to division the polynomialbased algorithm introduced in [START_REF] Jeannerod | Computing floating-point square roots via bivariate polynomial evaluation[END_REF] for square rooting. This extension provides a set of approximation and evaluation error conditions that will be proven to be sufficient to ensure correct rounding.

The second step consists in generating an efficient polynomial evaluation program, and in the automatic validation of its accuracy. On the one hand, the generation relies on several heuristics aimed at maximizing ILP exposure. On the other hand, checking that such a program is accurate enough turns out to be more involved than for square root, and will require a specific splitting strategy.

Third, besides the polynomial evaluation program, additional subroutines and associated C code sequences are proposed. Those include sign and exponent computation, rounding, and handling overflow, underflow and special values. This provides a complete division code (for rounding to nearest and without subnormal numbers), optimized for the ST231 processor yet portable.

The paper is organized as follows. After some preliminaries and notation in Section 2, the three contributions above will be described in Sections 3, 4, 5, respectively. Some experimental results obtained with the ST231 processor and the ST200 compiler will be reported in Section 6, and we shall conclude in Section 7.

Preliminaries and notation

Floating-point data and rounding. The floating-point data we shall consider are ±0, ±∞, quiet or signaling Not-a-Numbers (qNaN, sNaN), as well as normal binary floating-point numbers

x = (-1) sx • m x • 2 ex , (1)
with s x ∈ {0, 1}, m x = (1.m x,1 . . . m x,p-1 ) 2 and e x ∈ {e min , . . . , e max }. (In particular subnormal numbers [1, §3.3] will not be considered.) The precision p and extremal exponents e min and e max are assumed to be integers such that p ≥ 2 and e min = 1 -e max .

The rounding attribute chosen here is "to nearest even" (roundTiesToEven [1, §4.3.1]) and will be referred to as RN. Except for some special input (x, y) for which special values must be delivered as detailed in Section 5.4, the standard requires that RN(x/y) be returned whenever x and y are as in (1).

Correctly rounded division. That RN(x/y) essentially reduces to a correctly rounded ratio of (possibly scaled) significands can be recalled as follows. First, using RN(

-x) = -RN(x) gives RN(x/y) = (-1) sr • RN(|x/y|),
where s r is the XOR of s x and s y . Then, taking c = 1 if m x ≥ m y , and 0 otherwise, one has |x/y| = • 2 d , where

= 2m x /m y • 2 -c , d = e x -e y -1 + c. (2)
Since both m x and m y are in [1, 2), is in [1, 2) as well, and • 2 d will be called the normalized representation of |x/y|. Tighter enclosures of can in fact be given:

Property 1. If m x ≥ m y then ∈ [1, 2 -2 1-p ] else ∈ (1, 2 -2 1-p ). Proof. If m x ≥ m y then c = 1, and we deduce from 1 ≤ m x ≤ 2 -2 1-p and 0 < 1/m x ≤ 1/m y ≤ 1 that 1 ≤ ≤ 2 -2 1-p . If m x < m y then m x ≤ m y -2 1-p and thus ≤ 2 -2 2-p /m y . Hence, using m x ≥ 1 and 1/m y > 1/2, we obtain 1 < < 2 -2 1-p as desired.
These enclosures of will be used explicitly when handling underflow in Section 5.3. For now, we simply note that both of them give RN( )

∈ [1, 2 -2 1-p ], so that RN(|x/y|) = RN( ) • 2 d .
If e min ≤ d ≤ e max then we shall return the normal number RN(x/y) = (-1) sr • m r • 2 er , where

s r = s x ⊕ s y , m r = RN( ), e r = d. (3)
Else, d is either smaller than e min or greater than e max . Since subnormals are not supported, some special values will be returned in either case, as detailed in Section 5.3. Thus, to get RN(x/y) the main task consists in deducing from m x and m y the correctly rounded value RN( ) in (3), the sign s r and exponent e r being computable in parallel and at lower cost (see Section 5.1).

Correct rounding from one-sided approximations. Among the many methods known for getting RN( ) (see [5, §8.6] and the references therein), the one we focus on in this paper uses one-sided approximations: as shown in [5, p. 459], this method reduces the computation of RN( ) to that of an approximation v of such that

-2 -p < -v ≤ 0. ( 4
)
Here v is representable with, say, k bits while has in most cases an infinite binary expansion (1.

1 . . . p-1 p . . .) 2 .
Once such a v is known, correct rounding follows easily (see [5, p.460] and Section 5.2) and it remains to deduce from each possible pair (m x , m y ) a value v that satisfies [START_REF] Piso | Analysis of the impact of different methods for division/square root computation in the performance of a superscalar microprocessor[END_REF].

Before presenting in Sections 3 and 4 a novel approach for computing v in a certified and efficient way, we give next a brief description of the implementation context.

Software implementation on integer processors. Our implementation of division is for the binary32 format of [1]: k = 32, p = 24, e max = 127.

It will consist of a piece of C code using exclusively 32-bit integers, for input/output encoding as well as for intermediate variables.

Concerning the input data x, y the standard encoding into 32-bit unsigned integers X, Y is assumed [1, §3.4]. For example, the bit string X 31 . . . X 0 of X = 31 i=0 X i 2 i is such that X 31 = s x (sign bit of x), 7 i=0 X i+23 2 i = e x + 127 (biased exponent of x, for x normal), and X i = m x,23-i for i = 0, . . . , 22 (fraction bits of x). The output is encoded similarly and our division code eventually takes the form of a C function like the one below:

unsigned int binary32div( unsigned int X , unsigned int Y ) {...}
Concerning the operations on input or intermediate variables, we assume available the basic arithmetic and logical operators +, &, etc. as well as a max instruction and a multiplier mul defined for A, B ∈ {0, . . . , 2 32 -1} as mul(A, B) = A • B • 2 -32 . Here • denotes the usual floor function, and mul thus gives the 32 most significant bits of the exact product A • B.

Therefore, up to emulating max and mul as shown in Appendix A, all we need is a C compiler that implements 32-bit arithmetic. However some features of the ST231 processor and compiler have influenced the design and optimizations described in Sections 4 and 5. In particular, 4 instructions can be launched simultaneously, among which at most two mul instructions. Also, the latency of max and of the other basic operators is 1, while it is 3 for mul.

Sufficient conditions for correct rounding

This section gives sufficient conditions for (4) to hold. To do so we extend to division the bivariate polynomialbased framework introduced in [START_REF] Jeannerod | Computing floating-point square roots via bivariate polynomial evaluation[END_REF] for computing correctlyrounded floating-point square roots. A preliminary step is to restrict (4) to the following more symmetric, but only slightly stronger, condition

( + 2 -p-1 ) -v < 2 -p-1 .
(

) 5
Sufficient conditions to have (5) will be obtained by introducing a suitable bivariate polynomial approximant along with approximation and evaluation error bounds. We will rely on these conditions for designing our validation approach in Section 4.2, especially for providing hint values to the software tools Sollya and Gappa that we use.

Polynomial approximation

We start by interpreting the rational number + 2 -p-1 as the value of a suitable function of m x and m y . Defining F (s, t) = 2 -p-1 + s/(1 + t), we have that + 2 -p-1 is the exact value of F at the particular point

(s * , t * ) = (2m x • 2 -c , m y -1). (6)
When m x and m y vary then s * and t * range in the domains

S = [1, 2 -2 1-p ] ∪ [2, 4 -2 3-p ] and T = [0, 1 -2 1-p ].
(In particular, the second interval for S comes from the fact that c = 0 implies m x ≤ 2 -2 2-p .) Then, following [START_REF] Jeannerod | Computing floating-point square roots via bivariate polynomial evaluation[END_REF], the next step is to approximate F over S × T by a suitable bivariate polynomial P . Since F is linear with respect to s, one can reduce to univariate approximation by taking

P (s, t) = 2 -p-1 + s • a(t), (7)
with a(t) a polynomial approximating 1/(1 + t) over T . If we define

α(a) = max t∈T |1/(1 + t) -a(t)| , (8)
then the approximation error for F at (s * , t * ) satisfies

|F (s * , t * ) -P (s * , t * )| ≤ 4 -2 3-p α(a).
Since from (5) the overall error must be less than

2 -p-1 , we take α(a) < 2 -p-1 / 4 -2 3-p (9)
as a target approximation error bound for computing the approximant a.

The approximant construction is done using the software environment Sollya 1 , from the function 1/(1 + t), the domain T , and the error bound [START_REF] Kahan | A brief tutorial on gradual underflow[END_REF]. For reducing the final cost we choose a polynomial of smallest degree δ that satisfies the above constraints. For p = 24, the Sollya function guessdegree leads to δ = 10. Then, with the additional constraint that the coefficients should have no more than k = 32 fraction bits, we deduce a suitable polynomial a(t) from the remez function and by truncation. One obtains

a(t) = 10 i=0 a i t i , with each |a i | defined by a 32-bit un- signed integer A i such that |a i | = A i • 2 -32 . In our case, it turns out that 1 > |a 0 | > • • • > |a 10 |
and that the signs alternate, so that

a i = (-1) i • A i • 2 -32 ∈ (-1, 1), 0 ≤ i ≤ 10. (10)
Finally, a certified bound on the actual approximation error is obtained using Sollya's infnorm function:

α(a) ≤ 3 • 2 -29 ≈ 2 -27.41 < 2 -25 4 -2 -21 ≈ 2 -27 . (11)

Subdomain-based error conditions

Once a is available, for establishing [START_REF] Ercegovac | Digital Arithmetic[END_REF] we need to consider the rounding errors of the evaluation of P given by [START_REF] Hauser | The SoftFloat and TestFloat Packages[END_REF]. In [START_REF] Jeannerod | Computing floating-point square roots via bivariate polynomial evaluation[END_REF], for the design of a correctly-rounded square root function in a similar context, errors have been bounded by a single value, on the whole domain S × T . Here we extend this approach since the room left by [START_REF] Markstein | Computation of elementary functions on the IBM RISC System/6000 processor[END_REF] between the actual approximation error and 2 -p-1 , is not sufficient for taking into account rounding errors. (We refer to the experimental data in Section 4.2.) The validation of our division algorithm will require to see the interval T as a union of n subintervals: T = n i=1 T (i) and, accordingly, the approximation error α(a) of (8) will then split up into

α (i) (a) = max t∈T (i) |1/(1 + t) -a(t)| , 1 ≤ i ≤ n. (12)
The value v in (5) will result from the evaluation of P by a finite precision program P that produces, for each subdomain S × T (i) , the rounding error

ρ (i) (P) = max (s,t)∈S×T (i) |P (s, t) -P(s, t)| . ( 13
)
Let i be such that (s * , t * ) belongs to S × T (i) . Then, from ( 12) and ( 13), the overall error is eventually bounded as ( + 2 -p-1 ) -v ≤ 4 -23-p α (i) (a) + ρ (i) (P), and we arrive at the following sufficient conditions for (5): Property 2. If the approximation and rounding errors satisfy, for

1 ≤ i ≤ n, (4 -2 3-p ) α (i) (a) + ρ (i) (P) < 2 -p-1 (14)
then ( 5) holds.

Evaluation code generation and validation

Here we present our strategy to produce automatically an efficient evaluation program P along with an a posteriori certificate on its accuracy (in the sense of Property 2).

For p = 24, the program P will implement a finite precision evaluation of P (s * , t * ) = 2 -25 + s * • a(t * ). Here a is the polynomial obtained in Section 3.1 and whose coefficients a i satisfy [START_REF] Lauter | Arrondi correct de fonctions mathématiques -fonctions univariées et bivariées[END_REF]. The program will in fact store only the 32-bit integers A i , the coefficient signs being handled through an appropriate choice of arithmetic operators (+ or -) made when generating P. Concerning the evaluation point (s * , t * ) defined in [START_REF] Harrison | The computation of transcendental functions on the IA-64 architecture[END_REF], it can be encoded by a pair (S, T ) of 32-bit unsigned integers such that

s * = S • 2 -30 , t * = T • 2 -32 . (15)
The computation of S and T from X and Y can be implemented in C as follows: The above code has been written with ILP exposure in mind. For example, c is computed by comparing the trailing significands m x -1 and m y -1, rather than m x and m y . On ST231, this piece of code will typically take 3 cycles, the delay between S and T being of 2 cycles.

Concerning the output of P, note first that by Property 1, v in (4) must satisfy v ∈ [1, 2 -2 -p ). Moreover, because of the formats introduced in ( 10) and ( 15) and of the fact that P is essentially a fixed-point code, the output of P will be a 32-bit unsigned integer V representing v as

v = V • 2 -30 . ( 16
)
(How to implement correct rounding using this format will be detailed in Section 5.2.)

Evaluation program generation

Once the integers S, T, A 0 , A 1 , . . . , A 10 are available, we determine automatically an efficient program P for evaluating V by means of 32-bit additions/subtractions and mul instructions. By efficient program we mean a way of parenthesizing the arithmetic expression P (s * , t * ) that reduces the execution latency as well as the number of mul instructions. This depends on a set of heuristics, still under development, but that already enabled to produce an evaluation program suitable for floating-point division on ST231.

Evaluation tree generation. Following [START_REF] Harrison | The computation of transcendental functions on the IA-64 architecture[END_REF], we generate evaluation trees whose height is kept low through ILP exposure, thus avoiding highly sequential schemes like Horner's. To do so, we proceed in two substeps. Given the polynomial degree δ, the delay between S and T , and some latencies for addition/subtraction and mul, we first heuristically compute a target latency τ for P, assuming unbounded parallelism. This hint for the latency is a priori feasible, still reasonably low. Then we generate automatically a set of evaluation trees with height no more than this target latency. An example of such a tree is given in Figure 1, whose height corresponds to the target latency τ = 14. This latency2 is more than three times lower than the latency of (3 + 1) × 11 = 44 cycles that would result from computing P (s * , t * ) with Horner's rule. Our approach is heuristic in the sense that if no evaluation tree is found that satisfies the target latency τ , we increase τ and restart the process. However, the practice has shown that the number of evaluation trees found given the target τ is usually extremely large, already for degrees much lower than the degree δ = 10 we have here. Consequently, we have implemented several filters in order to reduce significantly this number during the generation, and thus to speed up the whole process.

Arithmetic operator choice. A first filter consists in restricting to evaluation trees for which all intermediate values are positive. This restriction allows to work with the full precision k = 32, instead of loosing one bit because of signed arithmetic. Such special trees can be found by choosing the addition/subtraction operators appropriately, for example considering a 0 -(-a 1 t) rather than a 0 + a 1 t, for a 1 is negative. If the sign of one of the intermediate values computed by the tree changes when the input (S, T ) varies, then that evaluation tree is rejected. This first filter is implemented using the MPFI 3 library for interval arithmetic. (An interval containing zero is interpreted as a sign change for the corresponding variable.)

Scheduling verification. A second filter consists in checking if a given evaluation tree can be scheduled without latency increase on a simplified model of the real target architecture. In our case of division on ST231, the evaluation tree has a latency of 14 cycles in theory (that is, assuming unbounded parallelism) as well as in practice. In fact, the possible scheduling displayed in Evaluation code. Finally, the first evaluation tree that passes both filters is chosen, and the corresponding evaluation program is printed out as a piece of C code. In our case, the obtained evaluation program is presented in Listing 1.

unsigned int r0 = mul( T , 0xffffe7d7 ); unsigned int r1 = 0xffffffe8 -r0; unsigned int r2 = mul( S , r1 ); unsigned int r3 = 0x00000020 + r2; unsigned int r4 = mul( T , T ); unsigned int r5 = mul( S , r4 ); unsigned int r6 = mul( T , 0xffbad86f ); unsigned int r7 = 0xfffbece7 -r6; unsigned int r8 = mul( r5 , r7 ); unsigned int r9 = r3 + r8; unsigned int r10 = mul( r4 , r5 ); unsigned int r11 = mul( T , 0xf3672b51 ); unsigned int r12 = 0xfd9d3a3e -r11; unsigned int r13 = mul( T , 0x9a3c4390 ); unsigned int r14 = 0xd4d2ce9b -r13; unsigned int r15 = mul( r4 , r14 ); unsigned int r16 = r12 + r15; unsigned int r17 = mul( r10 , r16 ); unsigned int r18 = r9 + r17; unsigned int r19 = mul( r4 , r4 ); unsigned int r20 = mul( T , 0x1bba92b3 ); unsigned int r21 = 0x525a1a8b -r20; unsigned int r22 = mul( r4 , 0x0452b1bf ); unsigned int r23 = r21 + r22;

3 http://gforge.inria.fr/projects/mpfi/ unsigned int r24 = mul( r19 , r23 ); unsigned int r25 = mul( r10 , r24 ); unsigned int V = r18 + r25;

Listing 1. Generated evaluation program.

Evaluation program validation

We have validated the above evaluation code using Gappa4 . Gappa is a software tool intended to help verifying and formally proving properties on numerical programs. It manipulates logical formulas involving the inclusion of expressions in intervals, and allows to bound rounding errors in a certified way. The first validation step is to check that no overflow can occur. This is easily done with Gappa by verifying that no intermediate value can be greater than 2 32 -1.

For proving the correct rounding inequality (5), we use Property 2 for providing hints to Gappa. With the polynomial a computed in Section 3.1, and no subdivision of the domain T (that is, n = 1), the approximation error satisfies

α(a) ≤ θ 0 = 3 • 2 -29 ≈ 2 -27.41 . ( 17
)
According to [START_REF] Melquiond | De l'arithmétique d'intervalles à la certification de programmes[END_REF], we then take η 0 = 2 -25 -(4-2 -21 )•θ 0 , and use the hint

ρ(P) < η 0 ≈ 2 -26.99 (18)
as a sufficient condition on the rounding error during the evaluation of the program P. This approach succeeds in the case m x ≥ m y where, with the help of Gappa, we have checked that (18) is true. Note that the strict inequality is checked with Gappa through an inequality ρ(P) ≤ η 0for a small enough, positive .

In the case m x < m y , with s * ∈ [2, 4 -2 -21 ] and t * ∈ [2 -23 , 1 -2 -23 ], we need to split up the domain T . Indeed, we find points t * at which (18) is not satisfied. For example, for t * = 0.97490441799163818359375, with t * in a small enough interval T (i) (see the last row of Table 2), we can check that ρ(P) < ρ (i) (P) < η 4 ≈ 2 -26.77 . Nevertheless, considering the approximation error on the same interval, we can establish [START_REF] Melquiond | De l'arithmétique d'intervalles à la certification de programmes[END_REF] since, in compensation, it can be proven that α (i) (a) ≤ θ 4 ≈ 2 -27.49 < α(a).

For validating the division program, and finding adequate subintervals T (i) 's, we have implemented a splitting of T by dichotomy. This process has split the definition domain up into n = 36127 subintervals. For each of them, using Sollya, we first bound the approximation error α (i) (a). This gives α (i) (a) ≤ θ for some rational number θ. Then the sufficient condition ( 14) is established by checking that ρ (i) (P) < η = 2 -25 -(4 -2 -21 ) • θ. The dichotomy is illustrated by Table 2 where "no" in the last column means that the interval in the second column has to be split up. The bounds on α (•) (a) and ρ (•) (P) (where the • stands for the index of the interval in the subdivision) that are discovered with Sollya and Gappa are given in the third and fourth columns. For the exact values of the θ j 's and η l 's we refer to Appendix B.

Depth Subintervals α (•) (a) ≤ ρ (•) (P) < Does (14) hold? 1 I1,1 = [2 -23 , 1 -2 -23 ] θ1 ≈ 2 -27.41 η1 ≈ 2 -26.99 no 2 I2,1 = [2 -23 , 0.5 -2 -23 ] θ2 ≈ 2 -27.41 η2 ≈ 2 -26.99 yes I2,2 = [0.5, 1 -2 -23 ] θ1 ≈ 2 -27.41 η1 ≈ 2 -26.99 no • • • j Ij,1 = [2 -23 , 0.5 -2 -23 ] θ2 ≈ 2 -

Implementation of a complete division code

So far we have presented an efficient way of deducing from m x and m y an integer [START_REF] Piso | Analysis of the impact of different methods for division/square root computation in the performance of a superscalar microprocessor[END_REF]. To obtain a complete code for binary floating-point division it remains to compute the sign and the exponent of the result, to deduce from v a correctlyrounded significand, and to pack those three fields according to the standard encoding of the binary32 format. If either x or y is a special operand then a special value must be returned, as prescribed in [1], which requires specific handling. The following sections detail algorithms for each of those tasks along with some C codes well-suited for a target like the ST231 processor. (All the variables are unsigned int, except for D which is int.)

V such that v = V • 2 -30 = (1.v 1 . . . v 30 ) 2 satisfies
From Section 5.1 to Section 5.3, both x and y are supposed to be normal numbers, while in Section 5.4 at least one of them is special, that is, ±0, ±∞, qNaN or sNaN.

Sign and exponent computation

The sign s r of the result r is trivially obtained by taking the XOR of the sign bits of X and Y :

Sr = ( X ˆY ) & 0x80000000;
The result exponent e r will be obtained by first computing the integer

D = d + e max -1. ( 19
)
If e min ≤ d ≤ e max then r is normal and the bits of its biased exponent E r = D + 1 will be deduced by adding the correctly-rounded significand to D • 2 p-1 (see Section 5.2); otherwise, since subnormals are not supported, r must take particular values like ±0, ±2 emin or ±∞, and we will see in Section 5.3 how such situations can be detected directly by inspecting the integer D.

Let us now compute D. Since x and y are normal, their biased exponents are E x = e x + e max and E y = e y + e max . Therefore, using [START_REF] Raina | FLIP: a Floating-point Library for Integer Processors[END_REF] together with the definition of d in (2),

D = E x -E y + e max -2 + c. ( 20
)
When k = 32, p = 24, and e max = 127, an implementation of ( 20) that exposes ILP is: However, we tend to prefer the version involving absX and absY, for these absolute values will be reused when computing special values in Section 5.4.

Rounding and packing

Given v that approximates from above as in (4), the correctly rounded significand m r = RN( ) can be deduced essentially as in [START_REF] Jeannerod | Computing floating-point square roots via bivariate polynomial evaluation[END_REF] and [5, p. 459] (see also [START_REF] Schwarz | Rounding for quadratically converging algorithms for division and square root[END_REF][START_REF] Oberman | Fast IEEE rounding for division by functional iteration[END_REF] and [17, §3.3]):

defining w = (1.v 1 . . . v p-1 v p ) 2
as the truncated value of v to p fraction bits, one has m r = w truncated to p-1 fraction bits if w ≥ , w+2 -p truncated to p-1 fraction bits if w < .

To check this, note first that by truncation to p fraction bits, 0 ≤ v -w < 2 -p , which together with (4) leads to

| -w| < 2 -p . (21)
Then it remains to verify for v p = 0 and for v p = 1, that the above definition of m r always gives RN( ), using in particular the classic fact that cannot be exactly halfway between two consecutive normal floating-point numbers (see for example [11, Lemma 1] or [3, p. 229]). To implement the above definition of m r one has to compute w and to decide whether w ≥ . Since v = V • 2 2-k and the bit string of V is 01v 1 . . . v k-2 , truncating v to p fraction bits means zeroing out the last k -p -2 bits of V . For (k, p) = (32, 24) the line below sets to zero the last six bits of the bit string 01v 1 . . . v 24 v 25 v 26 v 27 v 28 v 29 v 30 of V :

W = V & 0xFFFFFFC0;
The bit string of W is thus 01v 1 . . . v 24 000000 and w = W • 2 -30 . To know if w ≥ or not, let us introduce the integer M y such that m y = M y • 2 1-k . With k = 32 and S as in [START_REF] Oberman | Fast IEEE rounding for division by functional iteration[END_REF], we have the following characterization: Property 3. The condition w ≥ is true if and only if the condition mul(W,My) >= (S >> 1) is true.

Proof. From w = W • 2 -30 , m y = M y • 2 -31 , and s * = S • 2 -30 , we deduce that w ≥ = s * /m y is equivalent to W • M y • 2 -32 ≥ S/2
. Now, by definition of mul as a floor function, we have

W • M y • 2 -32 -1 < mul(W, M y ) ≤ W • M y • 2 -32 . (22)
Using the two inequalities in (22) together with the fact that mul(W, M y ) and S/2 are nonnegative integers representable with 32 bits, one can check that w ≥ if and only if mul(W, M y ) ≥ S/2. In C, the latter condition reads mul(W,My) >= (S >> 1).

To deduce the value of

m r = (1.m r,1 . . . m r,p-1 ) 2 let M r be the integer such that m r = M r •2 1-p . It then follows from w = W • 2 2-k that M r = W • 2 -(k-p-1) if w ≥ , (W + 2 k-p-2 ) • 2 -(k-p-1) if w < .
Thus, for (k, p) = (32, 24), M r is obtained by shifting either W or W + 2 6 to the right by seven positions.

Let us now pack the bits of m r with the bits of s r and E r . Since in Section 5.1 we have in fact computed D = E r -1, removing the leading 1 in m r = (1.m r,1 . . . m r,p-1 ) 2 can be avoided: the k-bit integer that encodes the result r is

s r • 2 k-1 + D • 2 p-1 + M r . (23)
In particular D ∈ {0, . . . , 2e max -1} implies that D • 2 p-1 + M r must be less than 2 k and thus never propagates a carry to the sign bit. Since we expect the computation of s r and D to be much cheaper than that of M r , one may first take the bitwise OR of s r •2 k-1 and D•2 p-1 , and then only add M r . For (k, p) = (32, 24) and p = 24 this gives the following code sequence:

Overflow and underflow detection

The value of the integer D = E r -1 in [START_REF] Schwarz | Rounding for quadratically converging algorithms for division and square root[END_REF] can be used to detect overflow and underflow. Indeed, the normality assumption on x and y implies that both E x and E y lie in the normal (biased) exponent range {1, . . . , 2e max }. Since c is either 0 or 1, it follows from (20) that D ranges from -e max -1 to 3e max -2. This range contains {0, . . . , 2e max -1}, for which the result is a normal number; it also contains two extremal ranges, for which either overflow or underflow occurs.

Let us start with overflow. If D ≥ 2e max then E r = e r + e max gives e r ≥ e max + 1. Therefore |x/y| ≥ 2 emax+1 and, by definition of the rounding operator RN(•) in [1], we have RN(x/y) = (-1) sr ∞.

For e max = 127 this case can be implemented as:

if( D >= 0xFE ) return Sr | 0x7F800000;
Now for underflow. If D < 0 then e r ≤ e min -1. Using < 2 then gives |x/y| < 2 emin . Following [START_REF] Hauser | The SoftFloat and TestFloat Packages[END_REF] we round |x/y| to nearest-even as if subnormals were supported. This rounded value can be either 2 emin (the smallest positive normal number) or a subnormal number. By definition of RN(•) the first case occurs if and only if

(1 -2 -p ) • 2 emin ≤ |x/y| < 2 emin , (24)
and we shall return the normal number (-1) sr 2 emin . In the second case, since we do not support subnormals, we shall return (-1) sr 0. Using [START_REF] Schwarz | Rounding for quadratically converging algorithms for division and square root[END_REF] with c = 1 shows that the condition e xe y = -e max is equivalent to D = -1. With M x and M y computed as in Sections 4 and 5.2, Property 4 thus suggests the following code for handling underflow:

if( D < 0 ) { if( D == -1 && Mx == 0xFFFFFF00
&& My == 0x80000000 ) return Sr | 0x00800000; else return Sr; }

Computing special values

Assume now that (x, y) is a special input, that is, x or y is ±0, ±∞, qNaN, or sNaN. For each possible case the standard [1] requires that a special value be returned. These special values follow from those in Table 3

-point datum x 0 +0 [2 p-1 , 2 k-1 -2 p-1 ) positive normal number 2 k-1 -2 p-1 +∞ (2 k-1 -2 p-1 , 2 k-1 -2 p-2 ) sNaN [2 k-1 -2 p-2 , 2 k-1
) qNaN Table 4. Floating-point data encoded by X.

which are X and Y modulo 2 k-1 (and are computed in Section 5.1 for k = 32), to filter out all special input (x, y). Indeed, x or y is in {±0, ±∞, qNaN, sNaN} if and only if

absX or absY is in {0} ∪ {2 k-1 -2 p-1 , . . . , 2 k-1 -1}, which is equivalent to max(absXm1, absYm1) ≥ 2 k-1 -2 p-1 -1, (25)
with absXm1 = (absX -1) mod 2 k and absYm1 = (absY -1) mod 2 k . Assume now that (x, y) has been filtered out by (25). Table 3 then requires that a qNaN be returned if absX equals absY, or if max(absX, absY) > 2 k-1 -2 p-1 . Otherwise, if x is finite (that is, absX < 2 k-1 -2 p-1 ) and y is nonzero (that is, absY is nonzero) then ±0 is returned, else ±∞ is returned. Hence the code below when (k, p) = (32, 24): Here absXm1 and absYm1 are unsigned int, so that addition is done modulo 2 32 . Note also that thanks to the use of the Max variable the payload of the qNaN result is one the input payloads, as recommended in [1, §6.2.3].

Experimental results

We report here some experiments done with the complete division code that follows immediately from the C codes of Sections 4 and 5.

Validation of the complete division code. Although our code has not been tested exhaustively, it has been validated using two classical test programs for IEEE binary floatingpoint arithmetic, and especially for division: the Extremal Rounding Tests Set [START_REF] Matula | Extremal rounding test sets[END_REF] and the TestFloat package [START_REF] Hauser | The SoftFloat and TestFloat Packages[END_REF]. This was done by compiling the code with Gcc under Linux and by using the software implementations of max and mul given in Appendix A.

Results obtained with the ST200 VLIW compiler. The same code has then been compiled with the ST200 VLIW compiler (version 2-0-0-A 20061215) based on the Open64 compiler technology, re-targeted to the ST200 family.

This compiler can produce an annotated assembly output, with specific scheduling annotations giving the date at which each bundle (a group of up to 4 instructions) is scheduled. In addition, one peculiarity of the ST200 architecture is that it provides a partial predication support in the form of conditional selection instructions, that are used by the compiler [START_REF] Bruel | If-conversion SSA framework for partially predicated VLIW architectures[END_REF] to generate branch-free code for all the subroutines described in this article.

Using optimization level -O3, the assembly produced in our case indeed consists of fully if-converted straight-line code in which the number of instructions as well as the latency is the same regardless of the nature of the input (normal numbers or special values). The table below collects these two values together with the number of instructions per cycle (IPC) and the code size: Since the ST231 has four issues, the IPC value indicates clearly the parallel nature of our approach.

Also, since one cycle is necessary for the last two instructions to select and return the result, the other 85 instructions must have been scheduled on the 4 issues during the first 26 cycles. This value of 26 cycles is close to the ideal value of 22 = 85/4 .

Finally, in the same context and for the same problem (that is, software implementation of binary32 floating-point division, correctly-rounded to nearest even and without subnormal numbers), the previously fastest implementation had a latency of 48 cycles [19, p. 104]. Thus, the speed-up brought by our approach is by a factor of 1.78.

Conclusion and future work

In this paper, we have focused on the binary32 format and on a particular 32-bit architecture. However, the approach we have presented in Sections 4 and 5 for generating/validating polynomial evaluation codes and for designing complete division codes, should in principle be usable for other values of the key parameters k, p, and e max . This could be of interest when optimizing division codes for other processors and/or other floating-point formats (like binary64, binary128, or smaller formats used in computer graphics).

This work also suggests two research directions, which we are currently investigating. First, gradual underflow is important [START_REF] Kahan | A brief tutorial on gradual underflow[END_REF] and our division code should definitely be extended in order to support subnormal numbers while keeping the latency as low as possible; similarly, all the rounding attributes prescribed in [1] should be implemented. Second, in some important special cases such as inversion (instead of division) or faithfully rounded results (instead of correctly rounded results), further optimizations should be easy to implement and validate by using the generation/validation tools and the design method that we have presented here.

  Tx = X << 9; T = Y << 9; X8 = X << 8; c = Tx >= T; Mx = X8 | 0x80000000; S = Mx >> c;

Figure 1 .

 1 Figure 1. Generated evaluation tree.

  absX = X & 0x7FFFFFFF; absY = Y & 0x7FFFFFFF; Ex = absX >> 23; Ey = absY >> 23; int D = (Ex + 125) -(Ey -c); Note that the biased exponent values E x and E y have been obtained by first computing the absolute values of x and y (by setting the sign bits to zero), and then shifting right by p -1 = 23 positions. Of course, there are other ways of extracting these biased exponent values, such as, for example for E x , Ex = ( X >> 23 ) & 0xFF; or Ex = ( X << 1 ) >> 24;

  My = (Y << 8) | 0x80000000; if( mul(W,My) >= (S >> 1) ) return (Sr | (D << 23)) + (W >> 7); else return (Sr | (D << 23)) + ((W + 0x40) >> 7);

Property 4 .

 4 One has (24) if and only if e x -e y = -e max and m x = 2 -2 1-p and m y = 1.Proof. Assume that (24) holds. Then it follows from |x/y| = • 2 d with ∈ [1, 2) that d must be equal to e min -1 = -e max , and that ≥ 2 -2 1-p . Property 1 implies further that = 2 -2 1-p and m x ≥ m y . Using the definition of in (2) then gives m x /m y = 2 -2 1-p and, since 1 ≤ m x , m y ≤ 2 -2 1-p , we conclude that the only possible value for (m x , m y ) is (2 -2 1-p , 1). Applying (2) with (c, d) = (1, -e max ) gives furthermore e x -e y = -e max . Conversely, one may check that if (m x , m y ) = (2 -2 1-p , 1) and e x -e y = -e max then |x/y| = (1 -2 -p ) • 2 emin .

Table 1 . Feasible scheduling on ST231.

 1 Table 1 uses only 3 out of the 4 issues offered by the ST231. For now, this second filter is implemented using a naive list scheduling algorithm.

		Issue 1	Issue 2	Issue 3	Issue 4
	Cycle 0	r0	r4
	Cycle 1	r6	r13
	Cycle 2	r11	r20
	Cycle 3	r1	r5	r22
	Cycle 4	r2	r14	r19
	Cycle 5	r12	r15	r21
	Cycle 6	r7	r10	r23
	Cycle 7	r3	r8	r24
	Cycle 8	r16
	Cycle 9	r17
	Cycle 10	r9	r25
	Cycle 11
	Cycle 12	r18
	Cycle 13	V

Table 2 . Splitting steps.

 2

	27.41	η2 ≈ 2 -26.99	yes

  by adjoining the correct sign, using x/y = (-1) sr |x|/|y|. (The standard does not specify the sign of a NaN result [1, §6.3].)

		|x|/|y|	+0	normal	|y|	+∞	NaN
		+0	qNaN	+0		+0	qNaN
	|x|	normal +∞	+∞ +∞	RN(|x|/|y|) +∞	+0 qNaN	qNaN qNaN
		NaN	qNaN	qNaN		qNaN	qNaN

Table 3 . Special values for |x|

 3 /|y|.

	The standard binary encoding [1, §3.4] implies in par-
	ticular

Table 4 .

 4 This table allows to reuse absX and absY,

	Value or range of integer X	Floating

http://sollya.gforge.inria.fr/ and[START_REF] Lauter | Arrondi correct de fonctions mathématiques -fonctions univariées et bivariées[END_REF].

On the tree, it corresponds to computing r 4 , and then r 22 up to V in

+ 3 + 1 + 3 + 3 + 1 = 14 cycles.

http://lipforge.ens-lyon.fr/www/gappa/ and[START_REF] Melquiond | De l'arithmétique d'intervalles à la certification de programmes[END_REF].

Acknowledgements

This work was supported by "Pôle de compétitivité mondial" Minalogic and by the ANR project EVA-Flo.

A. Implementing the max and mul instructions A C99 implementation of max and mul is as follows: The ST200 compiler is able to generate a single instruction when compiling the max or the mul function: it generates respectively the maxu and mul64hu ST200 instructions. This is interesting, since there is no need to write code using specific instrinsic functions, while reaching best efficiency.

B. Exact values of error bounds

C. Useful hexadecimal constants

For convenience' sake the table below displays the decimal value and/or the bit string of each of the hexadecimal constants appearing in some code sequences of Section 4 and Section 5.