F R Éd

Éric Goualard

F Goualard

A How do you compute the midpoint of an interval?

Keywords: G.1.0 [General]: Computer arithmetic, Error analysis, Interval arithmetic General Terms: Algorithms, Reliability, Experimentation floating-point number, IEEE 754 standard, interval arithmetic, midpoint, rounding error

The algorithm that computes the midpoint of an interval with floating-point bounds requires some careful devising to handle all possible inputs correctly. We review several implementations from prominent C/C++ interval arithmetic packages and analyze their potential failure to deliver the expected results. We then show how to amend them to avoid common pitfalls. The results presented are also relevant to non-interval arithmetic computation such as the implementation of bisection methods. Enough background on IEEE 754 floating-point arithmetic is provided for this paper to serve as a practical introduction to the analysis of floating-point computation.

INTRODUCTION

In his 1966 report [START_REF] Forsythe | How do you solve a quadratic equation[END_REF]] "How do you solve a quadratic equation? ", Forsythe considers the seemingly simple problem of reliably solving a quadratic equation on a computer using floating-point arithmetic. Forsythe's goal is both to warn a large audience away from unstable classical textbook formulae as well as get them acquainted with the characteristics of pre-IEEE 754 standard floating-point arithmetic, a dual objective shared by his later paper "Pitfalls in Computation, or Why a Math Book isn't Enough" [START_REF] Forsythe | Pitfalls in computation, or why a math book isn't enough[END_REF]]. Following Forsythe's track, we consider here the problem of computing a good approximation to the midpoint between two floating-point numbers. We strive to provide both a reliable algorithm for midpoint computation and an introduction to floating-point computation according to the IEEE 754 standard [IEEE 1985] (1)

But, if a and b are two floating-point numbers from a set F of finite floating-point numbers, the sum a + b may not be a floating-point number itself, and we therefore have to take care of rounding it correctly to ensure that our floating-point implementation of m(I) does not violate the fundamental property:

m(I) ∈ I, (2)
viz., that the floating-point midpoint of an interval I should belong to I. This property is a prerequisite to use the midpoint operator in the dichotomic search process of some numerical algorithm, even though it is not sufficient to ensure that the search space is separated into sub-spaces of the same size2 . There are, however, problems that require more from a midpoint operator, viz., that it computes the floating-point number closest to the real midpoint of I3 :

∀v ∈ F : |c -v| |c -m(I)| , with c = a + b 2 , c ∈ R. (3)
In interval arithmetic [START_REF] Moore | Interval Analysis[END_REF][START_REF] Alefeld | Introduction to Interval Computations[END_REF][START_REF] Hayes | A lucid interval[END_REF]], some examples of such problems are:

-On average, the least overestimation in defining a centered form [START_REF] Moore | Interval Analysis[END_REF][START_REF] Ratschek | Centered forms[END_REF]] of a rational function is obtained by choosing a value closest to the midpoint as developing point [Ratschek and Rokne 1984, p. 41]; -For the Krawczyk operator [Neumaier 1990, p. 177], the midpoint is the value that gives the tightest result over all other points of the interval under consideration [Neumaier 1990, Theo. 5.1.9, p. 178].

Inclusion and precision are not, however, the only possible problems for floating-point implementations of the midpoint operator: what should be, for example, the midpoint of an empty interval? As we shall see in Section 2, the IEEE 754 standard comes to the rescue here by defining a special Not-a-Number value (NaN) that should be used for such situation; another difficulty lies in the fact that, as we will see in the next section, one or both bounds of an interval with floating-point bounds may be an infinite value "±∞". What should then be the midpoint of an interval of the form [a, +∞], [-∞, b], or even [-∞, +∞]? Though the midpoint is mathematically undefined, many applications (e.g., constraint solvers such as IBEX [START_REF] Chabert | The IBEX constraint solver[END_REF] or Realpaver [START_REF] Granvilliers | Algorithm 852: Realpaver: an interval solver using constraint satisfaction techniques[END_REF]) typically expect finite values and would behave incorrectly, should the midpoint operator return a NaN for these intervals. As a consequence, the current draft of the future IEEE P1788 standard on interval arithmetic defines the midpoint operator as follows:

         m(∅) is NaN, m([-∞, b]) = -realmax, m([a, +∞]) = realmax, m([-∞, +∞]) = 0, m([a, b]) = rndnr (a + b)/2 , for (a, b) ∈ F 2 , (4)
where rndnr x is defined in the next section as returning the floating-point value closest to the real number x, and realmax is the largest positive element of F.

For non-empty intervals with finite bounds, Equation (3) implies Equation (2). However, some implementations may choose to ensure only the latter and a relaxation of the former.

Due to its many uses (centered forms, Krawczyk operator, Newton operator [START_REF] Moore | Interval Analysis[END_REF]], computation of preconditioning matrices [START_REF] Kearfott | Preconditioners for the interval Gauss-Seidel method[END_REF]], representation of intervals in the How do you compute the midpoint of an interval? A:3 midpoint/radius format [Rump 1999a], . . . ), the midpoint operator is a staple of interval arithmetic libraries. It is, therefore, paramount that its floating-point implementation at least satisfies Equation (2). Accuracy, as stipulated by Equation (3), is also desirable, as seen in the examples above. Nevertheless, we will see in Section 3 that some formulae implemented in popular C/C++ interval libraries may not ensure even the containment requirement for some inputs. However, this study should not be considered as a report on the quality of these libraries since different implementations of the midpoint operator might suit specific needs, and some libraries were defined with no provision to support intervals with infinite bounds in the first place. We are only concerned here with the properties of the formulae these libraries implement.

In Section 3, we analyze the various formulae both theoretically and practically; contrary to most expositions, we consider the impact of both overflow and underflow on the accuracy and correctness of the formulae.

The error analysis conducted in this paper requires slightly more than a simple working knowledge of floating-point arithmetic as defined by the IEEE 754 standard. As a consequence, the basic facts on floating-point arithmetic required in Section 3 are presented in Section 2 for the sake of self-containedness.

It turns out that the study of the correct implementation of a floating-point midpoint operator may serve as a nice introduction to many important aspects of floating-point computation at large: the formulae studied are simple enough for their analysis to be easily understandable, while the set of problems raised is sufficiently broad in scope to be of general interest. We then hope that this paper will be valuable as both a technical presentation of reliable, accurate, and fast methods to compute a midpoint as well as an introduction to error analysis of floating-point formulae.

FLOATING-POINT ARITHMETIC IN A NUTSHELL

According to the IEEE 754 standard [IEEE 1985], a floating-point number ϕ is represented by a sign bit s, a significand m (where m is a bit string of the form "0.f " or "1.f ", with f the fractional part) and an integral exponent E:

ϕ = (-1) s × m × 2 E .
(5)

The IEEE 754 standard defines several formats varying in the number of bits l(f ) and l(E) allotted to the representation of f and E, the most prominent ones being single precision-l(E), l(f ) = 8, 23 -and double precision-l(E), l(f ) = 11, 52 . We will also use for pedagogical purposes an ad hoc IEEE 754 standard-compliant tiny precision formatl(E), l(f ) = 3, 3 . Wherever possible, the significand must be of the form "1.f " since it is the form that stores the largest number of significant figures for a given size of m:

ϕ = 0.01101 ×2 0 = 0.1101 ×2 -1 = 1.101 ×2 -2 .
Floating-point numbers with such a significand are called normal numbers. Such prevalence is given to normal numbers that the leading "1" is left implicit in the representation of an IEEE 754 number, and only the fractional part f is stored (see Figure 1).

The exponent E is a signed integer stored as a biased exponent e = E + bias, with bias = 2 l(E)-1 -1. The biased exponent e is a non-negative integer that ranges from e min = 0 to e max = 2 l(E) -1. However, for the representation of normal numbers, E only ranges from E min = (e minbias) + 1 to E max = e maxbias -1 because the smallest and largest values are reserved for special purposes (see below). As an example of what precedes, the bias for the tiny format is equal to 3, e ranges from 0 to 7, and E ranges from -2 to +3. Consider the binary number ρ = 1.0011. It cannot be represented as a tiny floating-point number since its fractional part has four bits, and the tiny format has room for only three. It therefore has to be rounded to the floating-point number fl ρ according to one of four rounding directions4 (see Figure 2 Note the use of angles " " instead of more conventional parentheses for the rounding operators. They are used to express the fact that each individual value and/or operation comprising the expression is individually rounded according to the leading operator. For example:

fl ρ 1 + ρ 2 ≡ fl(fl(ρ 1 ) + fl(ρ 2 )), ∀(ρ 1 , ρ 2 ) ∈ R 2 .
When rounding to nearest, if ρ is equidistant from two consecutive floating-point numbers, it is rounded to the one whose rightmost bit of the fractional part is zero (the "even" number).

It is not possible to represent 0 as a normal number. Additionally, consider the number ϕ:

ϕ = 0.00011 × 2 0 .
To store it as a tiny floating-point normal number requires shifting the leftmost "1" of the fractional part to the left of the radix point:

ϕ = 1.100 × 2 -4 .
How do you compute the midpoint of an interval? A:5 However, doing so requires an exponent smaller than E min . It is nevertheless possible to represent ϕ, provided we accept to store it with a "0" to the left of the radix point:

ϕ = 0.011 × 2 -2 .
Numbers with a significand of the form "0.f " (with f = 0) are called subnormal numbers. Their introduction is necessary to reduce the large gap that would otherwise occur around 0 (compare Figure 3(a) and Figure 3(b)). The value "0" is also represented with the form "0.f " (with f = 0) but it is considered by the IEEE 754 standard as neither normal nor subnormal. To signal that a number is a subnormal number, the biased exponent e stored is set to the reserved value 0, even though the unbiased exponent E is E min (otherwise, it would not be possible to distinguish between a normal number and a subnormal number whose unbiased exponent is E min ). Exceptional steps must be taken to handle subnormal numbers correctly since their leading bit is "0", not "1". This has far reaching consequences in terms of performance, as we will see at the end of Section 3.

Figure 4 shows the repeated division by two of a tiny number across the normal/subnormal divide. As seen in that figure, dividing an IEEE 754 floating-point number by two is an error-free operation if the result is not a subnormal number. Otherwise, the rightmost bit of the fractional part is shifted out, and if it is non-zero, the operation is not error-free. The values λ and µ are, respectively, the smallest positive normal and the smallest positive subnormal floating-point numbers.

s e f 0 011 000 0 010 000 0 001 000 0 000 100 0 000 010 0 000 001 0 000 000

1.000 × 2 0 = 1 1.000 × 2 -1 = 0 .5 1.000 × 2 -2 = 0 .25 (λ) 0.100 × 2 -2 = 0 .125 0.010 × 2 -2 = 0 .0625 0.001 × 2 -2 = 0 .03125 (µ) 0.000 × 2 -2 = 0 ÷2 ÷2 ÷2 ÷2 ÷2 ÷2
Normal numbers Subnormal numbers Fig. 4. Repeated division by two from 1.0 to 0.0 in the tiny format (rounding to nearest/even).

The condition whereby an operation leads to the production of a subnormal number or an inexact zero is called underflow. On the other side of the spectrum, an operation may produce a number that is too large to be represented in the floating-point format on hand, potentially leading to an overflow. In that case, depending on the rounding policy in effect, the number is replaced either by an infinity of the right sign-which is a special floatingpoint number whose biased exponent e is set to e max and whose fractional part is set to zero-, or by the largest (in absolute value) finite floating-point number with the right sign5 . In the rest of this paper, we denote by F the set of normal and subnormal floating-point numbers and by F = F ∪ {-∞, +∞} its affine extension, the size of the underlying format (tiny or double precision, mainly) being unambiguously drawn from the context.

To ensure non-stop computation even in the face of a meaningless operation, the IEEE 754 standard defines the outcome of all operations, the undefined ones generating a NaN (Not a Number ), which is a floating-point datum whose biased exponent is set to e max and whose fractional part is any value different from zero. We will have, e.g.:

√ -1 = NaN, ∞ -∞ = NaN.
The NaNs are supposed to be unordered, and any test in which they take part is false6 . As a consequence, the right way to test whether an interval [a, b] is empty is to check whether ¬(a b) is true since that form returns a correct result even if either a or b is a NaN.

To sum up, the interpretation of a bit string representing a floating-point number depends on the values of e and f as follows: Figure 5 presents the distribution of all tiny numbers from F on the real line. The farther from 0, the larger the gap from a floating-point number to the next. More specifically, as Figure 6 shows, the difference between two consecutive floating-point numbers doubles every 2 l(e) numbers in the normal range, while it is constant throughout the whole subnormal range.

         e = 0, f = 0 : ϕ = (-1) s × 0.0 e = 0, f = 0 : ϕ = (-1) s × (0.f ) × 2 1-bias 0 < e < e max : ϕ = (-1) s × (1.f ) × 2 e-bias e = e max f = 0 : ϕ = (-1) s × ∞ e = e max f = 0 : ϕ = NaN. 0 0.
-0.9375 -1 -1.125 -1.25 -1.375 -1.5 -1.625 -1.75 -1.875 -2 -2.25 -2.5 -2.75 -3 -3.25 -3.5 -3.75 -4 -4.5 -5 -5.5 -6 -6.5 -7 -7.5 -8 -9 -10 -11 -12 -13 -14 -15 0
Rounding is order-preserving (monotonicity of rounding):

∀(ρ 1 , ρ 2 ) ∈ R 2 : ρ 1 ρ 2 =⇒ fl ρ 1 fl ρ 2
for the same instantiation of fl in both occurrences. The IEEE 754 standard also mandates that real constants and the result of some operations (addition, subtraction, multiplication, and division, among others) be correctly rounded : depending on the current rounding direction, the floating-point number used to represent a real value must be one of the two floating-point numbers surrounding it. Disregarding overflow-for which error analysis is not very useful-, a simple error model may be derived from this property: let ρ be a positive real number 7 , ϕ l the largest floating-point number smaller or equal to ρ and ϕ r the smallest floating-point number greater or equal to ρ. Let also rndnr ρ be the correctly s e f 0 010 010 0 010 001 0 010 000 0 001 111 . . . 0 001 001 0 001 000 . . . 0 000 001 0 000 000 rounded representation to the nearest-even of ρ by a floating-point number (rndnr ρ is either ϕ l or ϕ r ). If ρ is a floating-point number, we have:

1.010 • 2 -1 1.001 • 2 -1 1.000 • 2 -1 1.111 • 2 -2 1.001 • 2 -2 1.000 • 2 -2 0.001 • 2 -2 0.000 • 2 -2 2 -3 × 2 -1 2 -3 × 2 -1 2 -3 × 2 -2 2 -3 × 2 -2 2 -3 × 2 -2
rndnr ρ = ρ. (6)
Otherwise, we have:

|rndnr ρ -ρ| ϕ r -ϕ l 2 . (7)
We may consider two cases depending on whether ϕ l is normal or not:

-Case ϕ l normal. The real ρ may be expressed as m ρ × 2 E , with 1 m ρ < 2. Then, ϕ l can be put into the form m ϕ l × 2 E , with 1 m ϕ l < 2. If we call ε M the machine epsilon corresponding to the value 2 -l(f ) of the last bit in the fractional part8 of a floating-point number, we have ϕ r = ϕ l + ε M × 2 E . From Equation (7), we get:

|rndnr ρ -ρ| ε M 2 × 2 E .
We have9 :

|rndnr ρ -ρ| ρ (ε M /2) × 2 E m ρ × 2 E .
We know 1 m ρ < 2. Then:

|rndnr ρ -ρ| ε M 2 ρ
or, alternatively:

rndnr ρ = ρ(1 + δ), with |δ| ε M 2 . ( 8
)
-Case ϕ l subnormal. The distance between any subnormal number to the next floatingpoint number is constant and equal to µ = ε M × 2 Emin . From Equation ( 7), we then have:

|rndnr ρ -ρ| µ 2 ,
which may be expressed as:

rndnr ρ = ρ + η, with |η| µ 2 . ( 9
)
We may unify Equation ( 6), Equation ( 8), and Equation ( 9) as:

rndnr ρ = ρ(1 + δ) + η, with δη = 0, |δ| ε M 2 , |η| µ 2 , (10)
where either one of δ and η has to be null (δ = 0 in case of underflow, and η = 0 otherwise). This error model is valid for all correctly rounded operations and then:

∀(ϕ 1 , ϕ 2 ) ∈ F 2 , ∀ ∈ {+, -, ×, ÷} : rndnr ϕ 1 ϕ 2 = (ϕ 1 ϕ 2 )(1 + δ) + η, with δη = 0, |δ| ε M 2 |η| µ 2 . ( 11
)
For the other rounding directions, the same error model can be used, except that the bounds on δ and η are now larger:

|δ| < ε M and |η| < µ.

Table I. Characteristics of floating-point formats.

Format l(E) l(f ) E min Emax ε † M λ µ tiny 3 3 -2 +3 2 -3 2 -2 2 -5 single 8 23 -126 +127 2 -23 2 -126 2 -149 double 11 52 -1022 +1023 2 -52 2 -1022 2 -1074 †) machine epsilon, equal to 2 -l(f ) ) realmin (smallest positive normal number), equal to 2 E min ) subrealmin (smallest positive subnormal number), equal to ε M λ
Table I summarizes the values of the various constants encountered for the tiny, single, and double precision floating-point formats.

The error model presented in Equation ( 10) and Equation ( 11) is pessimistic by nature as it only gives an upper bound for the error. In some instances, it is possible to refine it to more accurate bounds. For example, we have for any finite floating-point number ϕ:

rndnr ϕ 2 = ϕ 2 + η with η ∈ - µ 2 , 0, µ 2 . ( 12
)
The same holds when rounding toward zero; for upward and downward rounding, the possible values for η are reduced to, respectively, {0, µ/2} and {-µ/2, 0}. The proof of Equation ( 12) is easy: in the absence of underflow, the division by 2 is errorfree since it amounts to decrementing the biased exponent by one; in case of underflow, the division by 2 is realized by shifting the fractional part one bit to the right (see Figure 4), which can only introduce an error of ±2

-(l(f )+1) × 2 Emin = ± ε M 2 × 2 Emin = ± µ 2 if the discarded bit is a "1" and no error if the bit is a "0".
We may also enjoy an error-free addition/subtraction for certain operands, as stated by the following three theorems, which will be used in the next section. Beforehand, note that we may consider for all practical purposes that, for ϕ 1 , ϕ 2 and ϕ 3 some floating-point numbers, fl ϕ 1 + ϕ 2 = fl ϕ 1 -(-ϕ 2 ) and fl ϕ 1ϕ 2 = fl ϕ 1 + (-ϕ 2 ) since taking the opposite of a floating-point number only involves flipping its sign bit, which is error-free.

Theorem 2.1 (Hauser's theorem 3.4.1 [START_REF] Hauser | Handling floating-point exceptions in numeric programs[END_REF]). If ϕ 1 and ϕ 2 are floatingpoint numbers, and if fl ϕ 1 + ϕ 2 underflows, then fl

ϕ 1 + ϕ 2 = ϕ 1 + ϕ 2 .
In effect, Theorem 2.1 means that η is always equal to 0 in Equation ( 11) when " " is an addition or a subtraction.

Theorem 2.2 (Hauser's theorem 3.4.1a [START_REF] Hauser | Handling floating-point exceptions in numeric programs[END_REF]). If ϕ 1 and ϕ 2 are floating-point numbers such that

|ϕ 1 + ϕ 2 | < 2λ, then fl ϕ 1 + ϕ 2 = ϕ 1 + ϕ 2 .
Theorem 2.3 (Hauser's theorem 3.4.2 [Hauser 1996] 10 ). If ϕ 1 and ϕ 2 are floating-point numbers such that

1/2 ϕ 1 /ϕ 2 2, then fl ϕ 1 -ϕ 2 = ϕ 1 -ϕ 2 .
Many other results to narrow the error bounds for some expressions may be found in the writings of Hauser [START_REF] Hauser | Handling floating-point exceptions in numeric programs[END_REF]], Goldberg [START_REF] Goldberg | What every computer scientist should know about floating-point arithmetic[END_REF]] or Higham [START_REF] Higham | Accuracy and Stability of Numerical Algorithms Second Ed[END_REF]], to name a few, but the ones presented here will be enough for our purpose.

Hauser's theorems give sufficient conditions to have error-free additions and subtractions. On the other hand, these operations may generate large errors through absorption and cancellation.

Absorption. Consider the two tiny floating-point numbers a = 1.000 × 2 3 and b = 1.000 × 2 -1 . To add them, we have first to align their exponents. This is done by scaling the operand with the smallest exponent:

a 1.000 ×2 3 b +0.0001 ×2 3 1.0001 ×2 3
The addition is performed with one extra bit for the significand of all values involved, the guard bit, which is necessary to ensure the validity of the error model of Equation ( 11). However, the result loses this guard bit when stored in the floating-point format on hand; if we round to the nearest-even, we then obtain a + b = 1.000 × 2 3 = a. The absolute error is:

|rndnr a + b -(a + b)| = 2 -1 ,
while the relative error is:

|rndnr a + b -(a + b)| |a + b| = 2 -1 2 3 + 2 -1 = 1 17 .
As expected by the error model in Equation ( 11), the relative error is bounded by ε M /2 = 2 -4 = 1/16. A loss of significance occurs whenever the difference in magnitude between the two operands of the addition leads to a shifting of some significant bits of one of the operands to the right of the rightmost bit of its fractional part. Absorption will occur when all the significant bits of one of the operands are shifted out.

Cancellation. Consider two nearly equal tiny numbers a = 1.001 × 2 3 and b = 1.000 × 2 3 and suppose they are the rounded approximations of earlier computation, respectively a = 1.00011 × 2 3 and b = 1.00010 × 2 3 . If we perform the subtraction ab, we get:

a 1.001 ×2 3 b -1.000 ×2 3 0.001 ×2 3
Even though the subtraction performed is error-free, we get the result "1" while a -b = 2 -2 . All correct figures from a and b cancel each other out and only the wrong figures remain and are magnified in the process-the more so for large a and b. As a rule of thumb, we should therefore avoid subtracting quantities spoiled by previous errors.

A:10 F. Goualard

THE FLOATING-POINT MIDPOINT OPERATOR

As said at the beginning of Section 1, a naive floating-point implementation of Equation ( 1) must consider the rounding of ρ = a + b, for a and b two floating-point numbers11 . A floating-point implementation of Equation ( 1) may fatally fail in two different ways with respect to the definition given by Equation ( 4):

Pitfall 1. For the interval [-∞, +∞], a + b returns a not-a-number (NaN), as defined by the IEEE 754 standard for the addition involving two infinities of differing signs. Mathematically, this does not violate the definition of midpoint, since the midpoint of [-∞, +∞] is undefined. From a practical point of view, however, it violates Property (2) since a NaN would not be considered to be between -∞ and +∞. As said previously, for many current applications, it is better to return a finite value. Otherwise, the user must take proper care to check for the occurrence of a NaN each time some midpoint is computed, even for non-empty intervals; Pitfall 2. If a and b have the same sign and their sum overflows, a + b will be an infinity, which divided by two will still give an infinity12 .

We have presented in Section 1 the definition chosen for the midpoint operator. In the rest of this paper, we then require the implementation of midpoint to comply with Equation (2) and Equation (3), and to also have the properties:

m(I) is NaN ⇐⇒ I = ∅.
That is, the midpoint operator should return a NaN only for empty intervals. And:

∀I = [a, b] : a = -b =⇒ m(I) = 0.
That is, the midpoint of a symmetric interval should be 0 (the interval [-∞, +∞] being considered symmetric).

The most straightforward implementation of the midpoint operator-such as used in, e.g., FI LIB [START_REF] Hofschuster | Fi lib, eine schnelle und portable funktionsbibliothek für reelle argumente und reelle intervalle im ieee-double-format[END_REF]]-has only one equation rounded to the nearest:

m([a, b]) = rndnr a + b 2 . ( 13
)
Evidently, it falls into both pitfalls highlighted above, and will return an infinite value in case a + b overflows, and a NaN for the interval [-∞, +∞].

The C++ Boost Interval Arithmetic Library [START_REF] Brönnimann | The design of the Boost interval arithmetic library[END_REF]] uses a slightly more complex implementation of the midpoint operator, viz.:

m([a, b]) =    -realmax if a = -∞; realmax if b = +∞; rndnr a+b 2 otherwise. (Boost rev84)
The cases are not mutually exclusive; however, only one rule is used (the first one that applies, in the order given), noting that the actual implementation adds a zeroth case for empty interval inputs, for which it returns a NaN.

The error analysis for the third rule of Equation (Boost rev84)-also valid, obviously, for FI LIB-is fairly simple with the model of Equation ( 11):

rndnr a + b 2 = (a + b)(1 + δ 1 ) + η 1 2 (1 + δ 2 ) + η 2 .
Thanks to Theorem 2.1, we have η 1 = 0. We also have δ 2 = 0, since a division by two is error-free in the absence of underflow. Hence the simplified expression:

rndnr a + b 2 = (a + b) 2 (1 + δ 1 ) + η 2 , (14)
with |δ 1 | ε M /2 and η 2 ∈ {-µ 2 , 0, µ 2 } (see Equation ( 12)). In the absence of overflow/underflow, we then have a correctly rounded result:

rndnr a + b 2 = a + b 2 (1 + δ 1 ),
which is the best we can expect.

Thanks to the first rule, the Boost library will not return a NaN for m([-∞, +∞]) but -realmax. This result does not violate Equation ( 2), and Equation ( 3) is not relevant in the context of intervals with an infinite radius. It might, however, appear odd to the unsuspecting user to discover that the midpoint of a non-empty symmetric interval may not be 0.

That the Boost library may, in some instances, not fall prey to Pitfall 2 may come as a surprise; in fact, the occurrence of overflows depends on minute details of the actual compilation of the expression rndnr (a + b)/2 : for example, on a computer with an ix86based processor offering at least the SSE2 instruction set [Intel 2007, Chap. 11], there are essentially two ways to compute with double precision floating-point numbers:

(1) With the FPU instruction set and 80-bits registers stack [Intel 2007, Chap. 8]: computation is performed internally with an extended precision format having 64 bits of significand and 15 bits of exponent. It is still possible to restrict the size of the significand to 53 bits in order to emulate double precision13 but this feature is rendered almost useless by the fact that it is not possible to restrict the size of the exponent; (2) With the SSE2 instruction set and the SSE2 registers: computation is performed exclusively in double precision.

If a program that uses the Boost library is compiled without any particular option on an ix86-based processor, the expression rndnr (a + b)/2 might be computed as follows:

(1) The double precision variables a and b will be loaded into the FPU floating-point registers and promoted to the 80 bits extended format;

(2) The expression rndnr (a + b)/2 will be computed entirely with extended precision in the FPU registers, possibly incurring a first rounding in the extended format;

(3) The result of the expression will be demoted back into a double precision variable, possibly incurring a second rounding in the double format.

In this situation, the overflow cannot happen since rndnr a + b is always representable by a finite number in extended precision; once divided by 2, the result is again representable in double precision. Additionnally, the two consecutive roundings may invalidate our error analysis since, even in the absence of underflow and overflow, the resulting error may then get larger than half a unit in the last place in rounding to nearest. Now, suppose the program is compiled for an ix86 processor with a compiler option that instructs it to emit code to store back floating-point values into double precision variables after each operation 14 (this behavior has an obvious adverse effect on performance. It is useful, however, to promote repeatability of the computation across different platforms, some of which may not offer extended precision registers of the same size, if at all) 15 . In that case, the expression rndnr a+b 2 will be computed as follows:

(1) The double precision variables a and b will be loaded into the FPU floating-point registers and promoted to the 80 bits extended format;

(2) The expression rndnr a + b will be computed with extended precision in the FPU registers;

(3) The value of rndnr a + b will be stored back as a double precision floating-point number, at which point an overflow may occur; (4) The value of rndnr a + b will be loaded into an extended register, promoted to 80 bits, and then divided by 2; (5) The result of the division will be stored in a double precision variable with some possible rounding.

In this situation, the Boost implementation of the midpoint operator may fail to comply with its requirements, and will return an infinite value as the midpoint of any interval with finite bounds whose absolute value of the sum is greater than realmax. The same problem would arise on architectures that do not rely internally on extended precision floating-point registers.

A different way to achieve the same failure on an ix86 processor relies on the possibility for modern compilers to transparently-or upon explicit request by the programmer-use the SIMD facilities offered by the SSE2 instruction set: some floating-point expressions may be computed in double precision SSE2 registers instead of the FPU extended precision registers. Once again, rndnr a + b may overflow, as it will be computed in double precision only. 16 Such possibilities for unexpected failures are compounded by the advent of processors featuring new floating-point operations-such as fused multiply-add (fma), defined by the most recent version of the IEEE 754 standard for floating-point arithmetic [IEEE 2008]whose rounding characteristics may invalidate the error analysis performed in this paper.

To avoid the overflow problem that affects the Boost implementation, Intlab V5.5 [Rump 1999b] uses the following formula: The first rule protects explicitly against the creation of a NaN for the interval [-∞, +∞], while the second rule avoids an overflow by first dividing the operands by two before adding the results, instead of the other way round.

The error analysis is as follows:

rndup a 2 + b 2 = a 2 (1 + δ 1 ) + η 1 + b 2 (1 + δ 2 ) + η 2 (1 + δ 3 ) + η 3 .
Using Theorem 2.1 and Equation ( 12), we may simplify the expression as:

rndup a 2 + b 2 = a + b 2 + η 1 + η 2 (1 + δ 3 ), ( 15
)
15 However, due to the unavoidably larger range of the exponent, perfect repeatability cannot be guaranteed for all cases.

16 This transparent use of SIMD registers instead of FPU ones leads to another problem: the assembler instructions that set the rounding direction for the FPU and SSE2 instructions are separate; therefore, the rounding direction set by the programmer for some FPU computation will not apply if, unbeknownst to him, the compiler chooses to use SSE2 instructions instead. However, this problem disappears if the programmer can make use of recent enough versions of high level functions such as the C99 fesetround() and Microsoft's controlfp(), which set the rounding direction for both FPU and SSE2 intructions.

with η 1 ∈ {0, µ/2}, η 2 ∈ {0, µ/2}, and 0 δ 3 < ε M , since results are rounded upward now.

In the absence of underflow/overflow, Equation (15) simplifies to:

rndup a 2 + b 2 = a + b 2 (1 + δ 3 ),
meaning that the result is correctly rounded (though not to the nearest). Note that, in practice, the quality of this formula is not as good as Boost's since the bound on δ 3 in Equation ( 15) is twice the one on δ 1 in Equation ( 14). With this algorithm, the midpoint of an interval with one infinite bound (but not both) is this infinity (for example: m([-∞, 3]) = -∞ and m([3, +∞]) = +∞). In addition, a symmetric interval with finite bounds may have a midpoint different from 0: the midpoint of the interval [-µ, µ] is computed as µ because the expressions rndup 0.5 × -µ and rndup 0.5 × µ are rounded, respectively, to 0 and µ. For the interval [-µ, µ], this has the unfortunate effect of violating Equation (3) since 0 is a representable floating-point number closer to the midpoint of the interval.

Additionally and more importantly, Intlab V5.5 may also violate Equation ( 2). Consider the reduction of Equation ( 15) for subnormal inputs:

rndup a 2 + b 2 = a + b 2 + η 1 + η 2 .
We have η 1 + η 2 ∈ {0, µ/2, µ}, which means that in some cases, the result may not be correctly rounded. Take for example, the interval [µ, µ]. Since we have:

rndup 0.5 × µ = µ, it follows: m([µ, µ]) = 2µ / ∈ [µ, µ].
BIAS 2.0.8 [START_REF] Knüppel | PROFIL/BIAS-a fast interval library[END_REF]] follows a different path and computes the midpoint as follows:

m([a, b]) = rndup a + b -a 2 .
(BIAS 2.0.8) Despite using a subtraction in place of an addition, this algorithm has the same flaws as a naive implementation of midpoint: if ba overflows (say, because a and b are huge, with differing signs), the result will be infinite; additionally

, m([-∞, +∞]) is a NaN because b -a = +∞ --∞ = +∞ and a + (b -a) = -∞ + ∞ = NaN.
The error analysis is also interesting:

rndup a + b -a 2 = a + (b -a)(1 + δ 1 ) + η 1 2 (1 + δ 2 ) + η 2 (1 + δ 3 ) + η 3 ,
which reduces to:

rndup a + b -a 2 = a + (b -a) 2 (1 + δ 1 ) + η 2 (1 + δ 3 ), ( 16
)
with η 2 ∈ {0, µ/2}, 0 δ 1 < ε M and 0 δ 3 < ε M . In the absence of overflow/underflow, we get:

rndup a + b -a 2 = a + b 2 + δ 1 b -a 2 + δ 3 a + b 2 + δ 1 δ 3 b -a 2 .
In theory, BIAS's formula is then even worse than Intlab V5.5's in the absence of underflow.

On the other hand, the reduction of Equation ( 16) for the case of subnormal inputs shows that it is not affected by Intlab V5.5's problem:

rndup a + b -a 2 = a + b 2 + η 2 .
For subnormal inputs, BIAS computes a correctly rounded result.

CXSC 2.4.0 [START_REF] Hammer | C++ Toolbox for Verified Scientific Computing I: Basic Numerical Problems[END_REF]] offers two functions to compute the midpoint: a "fast" one, Mid(I), and an accurate one, mid(I). The fast one computes the midpoint by a formula akin to both BIAS's and Intlab's, except for the rounding, which is now downward: The error analysis for the "fast" midpoint operator of CXSC gives:

m([a, b]) = rnddn a + b 2 - a 2 . ( CXSC
rnddn a + b 2 - a 2 = a + b 2 (1 + δ 1 ) + η 1 - a 2 (1 + δ 2 ) + η 2 (1 + δ 3 ) + η 3 (1 + δ 4 ) + η 4 ,
which simplifies to:

rnddn a + ( b 2 - a 2 ) = a + b -a 2 + η 1 -η 2 (1 + δ 3 ) (1 + δ 4 ), ( 17
)
with η 1 ∈ {-µ/2, 0}, η 2 ∈ {-µ/2, 0}, -ε M < δ 3 0, and -ε M < δ 4 0.

In the absence of overflow/underflow, CXSC's formula is as bad as BIAS's:

rnddn a + ( b 2 - a 2 ) = a + b 2 + δ 3 b -a 2 + δ 4 a + b 2 + δ 3 δ 4 b -a 2 .
For subnormal inputs, we have:

rnddn a + ( b 2 - a 2 ) = a + b 2 + η 1 -η 2 , with η 1 -η 2 ∈ {-µ/2, 0, µ/2}.
Filib++ 2.0 [START_REF] Lerch | Filib++, a fast interval library supporting containment computations[END_REF]] uses the same formula as Intlab V5.5 rounded differently, with additional rules to avoid some of the pitfalls previously encountered:

m([a, b]) =    a if a = b, 0 if a = -b, rndnr a 2 + b 2 otherwise.
(Filib++ 2.0)

The error analysis for Filib++'s formula leads to the same expression as for Intlab V5.5, except for the bounds on the errors:

rndnr a 2 + b 2 = a + b 2 + η 1 + η 2 (1 + δ 3 ), ( 18
)
with η 1 ∈ {-µ/2, 0, µ/2}, η 2 ∈ {-µ/2, 0, µ/2}, and

|δ 3 | ε M /2.
As for Intlab V5.5, the result is correctly rounded whenever no overflow/underflow occurs. On the other hand, for subnormal inputs, we still get:

rndnr a 2 + b 2 = a + b 2 + η 1 + η 2 ,
with η 1 + η 2 ∈ {-µ, -µ/2, 0, µ/2, µ}, which means that the computed midpoint may be one of the floating-point numbers that surround the correctly rounded midpoint. For a degenerate interval reduced to one point, the first rule protects against computing a midpoint outside the interval, as was the case with Intlab V.5.5; obviously, the formula ensures that the midpoint of an interval with a width greater or equal to 2µ is included in it. For an interval whose bounds a and b are consecutive subnormal floating-point numbers (hence, with a width precisely equal to µ), it suffices to notice that exactly one of them has a rightmost bit of the fractional part equal to "1" and the other one has the corresponding bit equal to "0". Consequently, η 1 + η 2 ∈ {-µ/2, 0, µ/2} and the inclusion is, once again, ensured.

The second rule in Equation (Filib++ 2.0) is more general than the first one of Intlab V5.5: it protects against the computation of a NaN for the interval [-∞, +∞] as well as ensuring that symmetric intervals in general have 0 as their midpoint (another "flaw" of the Intlab V5.5 formula). The expression in the third rule cannot overflow for finite a and b, which makes the whole formula immune to the overflow problem. On the other hand, the midpoint of any interval with an infinite bound (but not both) is infinite (e.g.: m([3, +∞]) = +∞). This result does not violate Equation (2) and Equation (3); it may nevertheless lead to serious trouble when the returned value is used in, say, a bisection algorithm or for a centered form without proper care.

The recent Version 6 of the Intlab package has a different algorithm from Version 5.5 to compute the midpoint of an interval:

m([a, b]) =        0 if a = -∞ and b = +∞, b if a = -∞, a if b = +∞, rndnr a + ( b 2 -a 2 ) otherwise. (Intlab V6)
Once again, these rules are not mutually exclusive, but the first that applies is executed to the exclusion of the others. The first rule explicitly avoids the NaN problem for the interval [-∞, +∞]; the second and third rules ensure that the midpoint of an interval with one infinite bound (and one only) is finite and equal to the other-finite-bound. This is one possible choice among many to resolve the problem, though different from Equation (4).

Apart from the rounding direction, the formula in the fourth rule is the same as CXSC's formula, and consequently, an error analysis along the same lines may be made:

rndnr a + ( b 2 - a 2 ) = a + ( b 2 + η 1 ) -( a 2 + η 2 ) (1 + δ 3 ) (1 + δ 4 ), (19)
with

η 1 ∈ {-µ/2, 0, µ/2}, η 2 ∈ {-µ/2, 0, µ/2}, |δ 3 | ε M /2, and |δ 4 | ε M /2.
For subnormal inputs, it may violate Equation (3), as evident from the error formula:

rndnr a + ( b 2 - a 2 ) = a + b 2 + η 1 -η 2 ,
with η 1η 2 ∈ {-µ, -µ/2, 0, µ/2, µ}. For example, with Intlab V6's formula, we have:

m([-µ, µ]) = -µ.
On the other hand, for finite a and b, the fourth rule is immune to overflow, and it ensures that Equation ( 2) is satisfied, as proven by Prop. 3.1. Proposition 3.1. Given a and b two finite floating-point numbers, with a b, we have:

a rndnr a + b 2 - a 2 b.
Proof. Since a b, we have rndnr a/2 rndnr b/2 , thanks to rounding monotonicity. This is true even if an underflow occurs. As a consequence, rndnr b/2a/2 0, and then rndnr a + (b/2a/2) a, again by monotonicity of rounding. For finite a and b, rndnr b/2a/2 cannot overflow.

Let us now prove that rndnr a + (b/2a/2) b. We first prove that rndnr b/2a/2 ba, restricting the discussion to the case a = b since the proof is trivial otherwise. We have:

rndnr b 2 - a 2 = b 2 + η 1 - a 2 + η 2 (1 + δ 3 ) . ( 20
)
Using Equation (20), let us now find the conditions such that rndnr b/2a/2 > ba:

rndnr b 2 - a 2 > b -a ⇐⇒ (b -a) 1 + δ 3 2 + (η 1 -η 2 )(1 + δ 3 ) > b -a ⇐⇒ b -a < 2(η 1 -η 2 ) 1 + δ 3 1 -δ 3 .
We have

|δ 3 | ε M /2, |η 1 | µ/2, and |η 2 | µ/2. As soon as l(f ) > 1, (1 + δ 3 )/(1 -δ 3 ) is strictly smaller than 2.
Hence, ba has to be less than 4µ for rndnr b/2a/2 to be greater than ba.

We then deduce that rndnr b/2a/2 is safely strictly less than λ for any floating-point format of practical use (just use ba < 4µ in Equation (20)), that is rndnr b/2a/2 underflows. In that case, δ 3 = 0, thanks again to Theorem 2.1, and then we have:

rndnr b 2 - a 2 > b -a ⇐⇒ b -a 2 + η 1 -η 2 > b -a ⇐⇒ b -a < 2(η 1 -η 2 ).
With the known bounds on η 1 and η 2 , we now deduce that ba must be strictly less than 2µ for rndnr b/2a/2 to be greater than b-a. But b-a must be a strictly positive floating-point number (we have b > a, and b-a < λ implies that Theorem 2.1 applies). Consequently, b-a must be equal to µ, which means that b and a are two consecutive floating-point numbers.

As a result, either a or b has an "even" binary significand (rightmost bit equal to 0), which means that either η 1 = 0 or η 2 = 0. We may then refine our condition:

rndnr b 2 - a 2 > b -a ⇐⇒ b -a 2 + η > b -a with |η| µ/2 ⇐⇒ b -a < µ.
This condition is not possible for a = b since there is no positive floating-point number smaller than µ. Consequently:

rndnr b 2 - a 2 b -a ∀(a, b) ∈ F 2 .
We then have:

a + rndnr b 2 - a 2 b.
And, by monotonicity of rounding:

rndnr a + ( b 2 - a 2 ) b.
The interval arithmetic library Gaol 4.0RC [START_REF] Goualard | GAOL 4.0RC: Not Just Another Interval Arithmetic Library 5.0[END_REF]] implements the midpoint operator by taking advantage of Hauser's theorems:

m([a, b]) =        0 if a = -b, -realmax if a = -∞, realmax if b = +∞, rndnr (a -a 2 ) + b 2 otherwise.
(Gaol 4.0RC)

As before, only the first rule whose condition is true is considered. Case 4 rewrites the expression "a + (b/2a/2)" from Intlab V6 as "(aa/2) + b/2", a small algebraic manipulation indeed, but one that can have a significant impact on the accuracy of the whole expression. In effect, the error analysis gives:

rndnr a - a 2 + b 2 = a - a 2 (1 + δ 1 ) + η 1 (1 + δ 2 ) + η 2 + b 2 (1 + δ 3 ) + η 3 (1 + δ 4 ) + η 4 . (21)
Lemma 3.2. For any finite floating-point number a ∈ F, we have:

rndnr a - a 2 = a - a 2 + η , with η ∈ - µ 2 , 0, µ 2 .
Proof. There are two cases to consider: let us first suppose no underflow occurs when halving a. Then, we have rndnr a/2 = a/2 since halving is error-free in the absence of underflow. We may then use Hauser's Theorem 2.3 to state that rndnr aa/2 = aa/2. Now, suppose an underflow occurs when halving a; then, rndnr a/2 = a/2 + η. This event is only possible if a/2 is strictly smaller than λ (by definition), and this can happen only if a is strictly smaller than 2λ. But then, rndnr a -(a/2 + η) is also strictly smaller than 2λ, and we may then use Theorem 2.2 to state that rndnr aa/2 = a -(a/2 + η).

Using Lemma 3.2 and the property on halving given by Equation ( 12), we are able to simplify Equation ( 21) to:

rndnr a - a 2 + b 2 = a + b 2 + η 3 -η 1 (1 + δ 4 ), (22)
with

η 1 ∈ {-µ/2, 0, µ/2}, η 3 ∈ {-µ/2, 0, µ/2}, |δ 4 | ε M /2.
It is then straightforward to show that containment is always ensured: no overflow can occur in computing (aa/2) + b/2 and, in the absence of underflow, a correctly rounded midpoint is computed; on the other hand, in case of underflow, the error is bounded by µ and not µ/2. The benefit of using "aa/2" instead of "a/2" directly is apparent when we compare Equation ( 22) with Equation ( 18): when a = b, we have η 1 + η 2 = ±µ in the worst case for Equation (18), and η 3η 1 = 0 in all cases for Equation (22) (errors in halving a and b cancel out when a = b).

The formulae investigated so far all have defects with respect to the definition in Equation (4). Table II synthesizes the error formulae in three cases:

-GC. The general case, where an underflow may or may not occur at any step of the computation (for the experiments, overflows are allowed in the general case as well, although the formulae in Table II are then no longer relevant); -No OF/SN. The case where no overflow nor any underflow occurs at any step of the computation, and neither infinities nor subnormal numbers are used as bounds of the test intervals; -SN. The case where all quantities manipulated and produced (counting the initial inputs) are subnormal numbers.

(1 + δ 1 ) + η 2 a+b 2 (1 + δ 1 ) a+b 2 + η 2 rndnr (a + b)/2 |η 2 | ∈ {0, µ/2}, |δ 1 | ε M /2 Intlab V5.5 ( a 2 + η 1 ) + ( b 2 + η 2 ) (1 + δ 3 ) a+b 2 (1 + δ 3 ) a+b 2 + η 1 + η 2 rndup a/2 + b/2 η 1 ∈ {0, µ/2}, η 2 ∈ {0, µ/2}, 0 δ 3 < ε M BIAS a + (b-a) 2 (1 + δ 1 ) + η 2 (1 + δ 3 ) a + (b-a) 2 (1 + δ 1 ) (1 + δ 3 ) a+b 2 + η 2 rndup a + (b -a)/2 η 2 ∈ {0, µ/2}, 0 δ 1 < ε M , 0 δ 3 < ε M CXSC a + ( b 2 + η 1 ) -( a 2 + η 2 ) (1 + δ 3 ) (1 + δ 4 ) a + (b-a) 2 (1 + δ 3 ) (1 + δ 4 ) a+b 2 + η 1 - η 2 rnddn a + (b/2 -a/2) η 1 ∈ {-µ/2, 0}, η 2 ∈ {-µ/2, 0}, -ε M < δ 3 0,-ε M < δ 4 0 Filib++ ( a 2 + η 1 ) + ( b 2 + η 2 ) (1 + δ 3 ) a+b 2 (1 + δ 3 ) a+b 2 + η 1 + η 2 rndnr a/2 + b/2 |η 1 | ∈ {0, µ/2}, |η 2 | ∈ {0, µ/2}, |δ 3 | ε M /2 Intlab V6 a + ( b 2 + η 1 ) -( a 2 + η 2 ) (1 + δ 3 ) (1 + δ 4 ) a + (b-a) 2 (1 + δ 3 ) (1 + δ 4 ) a+b 2 + η 1 - η 2 rndnr a + (b/2 -a/2) |η 1 | ∈ {0, µ/2}, |η 2 | ∈ {0, µ/2}, |δ 3 | ε M /2, |δ 4 | ε M /2 Gaol 4.0RC a+b 2 + η 3 -η 1 (1 + δ 4 ) a+b 2 (1 + δ 4 ) a+b 2 + η 3 - η 1 rndnr (a -a/2) + b/2 |η 1 | ∈ {0, µ/2}, |η 3 | ∈ {0, µ/2}, |δ 4 | ε M /2
For the general case (GC), Table II presents worst-case error bounds that may not be reachable in practice. On close inspection, we can even refine these bounds. Consider for example the formula for Filib++:

rndnr a/2 + b/2 = ( a 2 + η 1 ) + ( b 2 + η 2 ) (1 + δ 3 ).
It is not possible for η 1 , η 2 , and δ 3 to be simultaneously non-null: according to Theorem 2.2, for δ 3 to be different from 0, we must have:

( a 2 + η 1 ) + ( b 2 + η 2 ) 2λ. ( 23
)
But, for η 1 and η 2 to be both different from 0, we must have |a| < 2λ and |b| < 2λ, and then, |( a 2 + η 1 )| and |( b 2 + η 2 )| are each at most equal to λ. As a consequence, the only case where Equation ( 23) is satisfied is when |( a 2 + η 1 ) + ( b 2 + η 2 )| is precisely equal to 2λ, for which case δ 3 = 0.

How do you compute the midpoint of an interval?

A:19

On the other hand, for the two other columns of the table, it is possible to construct examples for which all δs and ηs are non-null simultaneously, and the worst-case error bound is reached.

With respect to error bounds, the most promising expressions are:

rndnr (a + b)/2 (Boost), which has the best error bounds for all situations, but may return an infinite value if a + b overflows; -rndnr a/2 + b/2 (Filib++), which may return an infinite value for semi-unbounded intervals (unbounded intervals being treated separately); -rndnr (aa/2) + b/2 (Gaol 4.0RC).

The algorithms that use the first two expressions can be amended to remove the associated flaws:

-Boost's algorithm (Equation (Boost rev84)) can be supplemented by some a posteriori test to check whether the output overflowed; if it is the case, the midpoint is recomputed using a different formula for which an overflow cannot occur (e.g., rndnr a/2 + b/2 ); -Amending Filib++'s algorithm is easy by adding cases as follows:

m([a, b]) =              a if a = b, 0 if a = -b, -realmax if a = -∞, realmax if b = +∞, rndnr a 2 + b 2 otherwise. ( 24
)
From Theory to Practice

So far, we have only theoretically analyzed the various formulae to compute the midpoint of an interval [a, b]. According to Table II, we should be able to separate the various methods into four classes for the case where no underflow nor any overflow occurs:

-Boost, Filib++ and Gaol, which deliver a correctly rounded result to the unit roundoff (half the epsilon of the machine); -Intlab V5.5, which delivers a correctly rounded result to the epsilon of the machine; -Intlab V6, which delivers a result with an error proportional to the unit roundoff; -Bias and CXSC, which deliver a result with an error that may get twice as large as Intlab V6's for the same inputs. This is indeed verified by the experiments, as shown in Table III, where the first column "Disc" indicates the number of floating-point numbers between the correctly rounded midpoint and the midpoint computed by the method tested. The last row "Fail" in the table indicates the number of times a midpoint computed was not included in the interval input. The table reports the results for the computation of the midpoint of ten million non-empty random intervals.

The comparison of all formulae when subnormal numbers are produced for each operation shows less drastic differences (see Table IV): the best we can compute in that case is (a + b)/2 + η, with |η| ∈ {0, µ/2}, which is exactly what is computed by Boost. BIAS exhibits the same error bound, but the computed result may differ by one unit in the last place from Boost's, due to the upward rounding. Intlab V5.5 shows a discrepancy of one unit in the last place almost half of the time: the error is η 1 + η 2 with η 1 and η 2 both positive or null, which means that they can never cancel each other out. On the other hand, CXSC offers a good accuracy because the error term is η 1η 2 with both η 1 and η 2 being of the same sign. Lastly, Filib++, Intlab V6 and Gaol are of similar accuracy, with Intlab V6 and Gaol being rigorously equivalent, as Table II suggests already.

If we allow both underflows and overflows (see Table V), we now have five classes: -Boost, which always returns a correctly rounded result, except when overflows make it fail to ensure inclusion of the midpoint; -Gaol and Filib++, which never fail to meet the requirement of inclusion while never departing by more than one unit in the last place from the correctly rounded result; -Intlab V5.5, which delivers half of the time a result departing from the correctly rounded result by one ulp. Although it is not apparent from Table V, recall that it may fail to ensure containment in some rare cases (see page 13); -Intlab V6, whose results may significantly depart from the correctly rounded one, although it always meets the inclusion requirement; -BIAS and CXSC, whose results may significantly depart from the correctly rounded one, and which may not meet the inclusion requirement.

Discrepancies for intervals with at least one infinite bound are not reported in Table V.

Figure 7 compares the various methods performance-wise for two different platforms (Intel Xeon E5520 on Ubuntu 10.04.4, and Intel i7 "Sandy Bridge" on Ubuntu 11.04). All the methods have been implemented in C++ in the same environment. Times are given in seconds for the computation of the midpoint of fifty million non-empty random intervals, with the same three categories (GC, No OF/SN, SN) as in Table II. The C++ compiler used is GCC 4.8.0 with the options "-ffloat-store -O0 -mpc64 -frounding-math". Note the difference in time between the first two categories ("GC" and "No OF/SN") and the last one ("SN") for both architectures: subnormal numbers have to be handled as exceptional values, which slows the computation down dramatically. However, Intel Xeon and Intel i7 "Sandy Bridge" differ in their handling of subnormal numbers, the latter implementing part of its support in hardware rather than as microcode exceptions [Fog 2012, p. 102], which drastically reduces the overhead incurred in using subnormals. 3): returning a value different from the correctly rounded center of the interval in the absence of overflow.

B o o s t

Table VI summarizes the potential shortcomings of the midpoint implementation for the packages we have considered. Constraint (3) is violated by almost all packages. Pending further investigation, SYM and IR may be considered benign errors because they should not compromise the reliability of algorithms based on packages that do exhibit them. All other problems may put reliability in jeopardy for applications that expect a definition of the midpoint operator matching Equation (4).

The corrected versions of Boost's and Filib++'s formulae (Boost ∞ and Equation ( 24), respectively) are the algorithms with the fewest number of "defects". However, it seems easier to avoid the symmetry problem of Boost ∞ by returning 0 when the input is [-∞, +∞] than to correct Equation (24)'s rounding problem in presence of subnormal numbers.

Performance-wise, Boost ∞ is also not an unreasonable choice if we keep in mind that slightly faster methods may lead to incorrect results. Experiments also show that adding a test for the interval [-∞, +∞] as input to return 0 does not incur any significant overhead, and gives us an algorithm that matches the requirements of Equation ( 4). The resulting algorithm to compute the midpoint would then be the one given in pseudo-code in Table VII.

CONCLUSION

The IEEE Standards Association is currently supporting the definition of Standard P1788 [IEEE 1788[IEEE 2013] ] for interval arithmetic. There is an ongoing effort to rigorously define the properties that functions such as the midpoint operator should exhibit. As the numerous exchanges on the IEEE P1788 mailing list have shown, the consensus on these properties is not trivial to reach, as different properties may be needed for different applications, which might give rise eventually to several specialized midpoint operators17 . Motion 37, passed in November 2012, defined after some heated debate the midpoint operator as stated in Equation (4).

  1 . Given two real numbers a and b from R, with a b, the midpoint m(I) of the closed interval I = [a, b] = {x ∈ R | a x b} can easily be obtained with the straightforward formula: m([a, b]) = a + b 2 .

Fig. 1 .

 1 Fig. 1. Binary representation as a seven bit string of a tiny floating-point number.

Fig. 2 .

 2 Fig. 2. Rounding a real number according to the IEEE 754 standard.

Fig. 3 .

 3 Fig. 3. The tiny floating-point format with and without subnormals (focus on 0).

Fig. 5 .

 5 Fig. 5. IEEE 754 tiny floating-point (normal and subnormal) numbers.

Fig. 6 .

 6 Fig. 6. Gap between consecutive floating-point numbers in the tiny format.

  2.4.0) That formula is still subject to the NaN problem whenever [a, b]=[-∞, +∞]. The accurate version uses extended precision, and then, contrary to the fast version, it always returns the correctly rounded value for the midpoint, except for the interval [-∞, +∞]. For that interval, CXSC aborts abruptly with an exception.

Fig. 7 .

 7 Fig. 7. Performances of nine midpoint() implementations for two CPU architectures.

Table II .

 II Synthesis of worst-case error bounds.

		Error
	Method	GC	No OF/SN	SN
	Boost	a+b 2

Table III .

 III Comparison of the quality of the midpoint() implementations -No OF/SN.

	Disc.	Boost Intlab V5.5	BIAS	CXSC Filib++ Intlab V6 Gaol 4.0RC
	0 10000000	5008451 2508795 2508307 10000000 9864663	10000000
	1		4991549 7486334 7486883		132983
	2			3132	3028		1178
	3			460	467
	4			585	569		603
	6			103	139
	8			287	291		294
	12			18	34
	16			126	125		127
	24			7	3
	32			73	74		74
	48			1	2
	64			36	36		36
	96			1
	128			25	25		25
	256			8	8		8
	512			4	4		4
	1024			4	4		4
	16384			1	1		1
	Fail	0	0	0	0	0	0	0

Table IV .

 IV Comparison of the quality of the midpoint() implementations -SN.

	Disc.	Boost	Intlab V5.5	BIAS	CXSC Filib++	Intlab V6 Gaol 4.0RC
	0	10000000	5001983 7502568 7500491	6251593	6250468	6250468
	1		4998017 2497432 2499509	3748407	3749532	3749532
	Fail	0	0	0	0	0	0	0

Table V .

 V Comparison of the quality of the midpoint() implementations -GC.

	Disc.	Boost	Intlab V5.5	BIAS	CXSC Filib++	Intlab V6 H-method
	0	9990382	5002986 2506554 2506186	9990379	9855400	9990381
	1		4987404 7478969 7479251	11	132539	9
	2			3104	3144		1226
	3			428	468
	4			601	612		632
	6			105	116
	8			277	275		283
	12			31	20
	16			160	161		162
	24			6	8
	32			68	70		70
	48			2	1
	64			46	46		46
	96			2
	128			15	15		15
	256			6	6		6
	512			5	5		5
	1024			3	3		3
	8192			2	2		2
	32768			1	1		1
	65536	8		5
	Fail	8	0	4847	4842	0	0	0

  Boost formula where Filib++'s formula is used in case of overflow ( ) Filib++ formula with tests added to avoid returning an infinite value.

	10.62 GC No OF/SN SN	11.56	14.0	13.62	17.35	25.09	18.69	23.34	25.65 Intel Xeon
	3.73 3.84	4.01 4.04	3.47 3.38	3.64 3.32	3.04 3.51	3.47 3.15	4.09 3.96	3.85 3.79	4.28 4.26
	3.5 3.68 6.01	3.4 3.34 6.28	2.86 3.21 7.83	2.97 2.8 5.66	3.11 3.06 7.49	2.98 3.26 8.07	3.47 3.39 9.13	3.77 3.52 8.67	3.37 3.79 8.78 Intel i7
	( †)

Table VI .

 VI Midpoint formulae shortcomings. NaN. Returning a NaN for a non-empty interval (e.g., for m([-∞, +∞])). Inf. Returning an infinite midpoint (e.g., in case of overflow in the addition a + b, or if a bound is infinite). NC. Failing to meet the containment requirement (e.g., m([a, b]) / ∈ [a, b] in case of underflow). SYM. Violation of Equation (3): not returning 0 for a non-empty symmetric interval. IR. Violation of Equation (

	Package	NaN	Inf	NC	SYM	IR
	Boost rev84 (P. 10)
	Boost∞
	Intlab V5.5 (P. 12)
	BIAS 2.0.8 (P. 13)
	CXSC 2.4.0 (Mid) (P. 14)
	Filib++ 2.0 (P. 14)
	Intlab V6 (P. 15)
	Gaol 4.0RC (P. 17)
	Equation (24) (P. 19)

A new version of the original standard from 1985 was approved and released in

[IEEE 2008]. Most prominently, it considers both binary and decimal representations of floating-point numbers, while the 1985 version only considered a binary representation. In this paper, we only consider the features standardized by the 1985 version-by far the most used to this day-for the sake of simplicity.

Where the size w(I) of an interval I = [a, b] with floating-point bounds is defined as w(I) = ba.

If the real midpoint is midway between two floating-point numbers, there exist two possible values for m(I) that satisfy Equation (3).ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY.

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

The actual direction chosen may depend on the settings of the Floating-Point Unit at the time, or alternatively, on the machine instruction used, for some architectures.ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY.

For example, if the sum a + b overflows with a + b > 0, we have rnddn a + b = realmax and rndup a + b = +∞.

Surprisingly enough to the unsuspecting, even the equality test turns false for a NaN, so much so that the statement x = x is an easy way to check whether x is a NaN.

The negative case may be handled analogously.ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY.

Alternatively, ε M may be defined as the difference between 1.0 and the next larger floating-point number.

Here, ρ is different from 0 since it is not a floating-point number at this point.ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY.

The better known Sterbenz's theorem[START_REF] Sterbenz | Floating-point Computation[END_REF]] states the same result in a restricted form where underflows are forbidden[START_REF] Higham | Accuracy and Stability of Numerical Algorithms Second Ed[END_REF]].ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY.

In the rest of this paper, we consider the double precision format, even though all results are applicable to other binary formats defined by the IEEE 754 standard[IEEE 1985] as well.

See below for a proviso on this, though.ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY.

Since version 4.3.0, the GNU C Compiler gcc offers the option -mpc64 to do that easily.

This is the -ffloat-store option of gcc, or the -mp option for the Intel C++ compiler.ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY.

Some applications might even benefit from different versions of the midpoint operator. A dichotomic process, for example, might take advantage of using a geometric mean for wildly unbalanced intervals, while resorting to the "regular" definition otherwise.ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY.

ACKNOWLEDGMENTS

The origin of the work can be traced back to some email exchanges with Prof. Arnold Neumaier; Prof. Siegfried Rump read a preliminary version of this paper and gave very valuable advice and feedback, pointing out in particular an error in the first version of the proof of Prop. 3.1. Dr. Alexandre Goldsztejn

Part of the work presented here was funded by IFCPAR/CEFIPRA Project 4502-1.

A:23

Table VII. Algorithm to compute the midpoint of an interval according to the specifications of the draft of the future IEEE 1788 Standard on Interval Arithmetic.

def midpoint ([ a , b ]):

if [a , b ] In particular, not everyone agrees on what the midpoint of unbounded or semi-bounded intervals should be. For semi-bounded intervals, Arnold Neumaier [Neumaier 2010] argues for returning the other, finite, bound of the interval, his justification being that this is a sensible choice when using a midpoint in the context of centered forms [START_REF] Ratschek | Computer Methods for the Range of Functions[END_REF]. On the other hand, that choice is not so useful in the context of a dichotomic exploration algorithm (Neumaier's stance on the subject is that we should use a different, specialized operator for this kind of algorithm).

For exploration algorithms, a simple choice seems to be the one advocated by Motion 37, viz. to return ±realmax. Obviously, any finite value from the interval could be used too. In that situation, Zuras and Hayes [START_REF] Zuras | Midpoint and unbounded intervals[END_REF] suggest to return the midpoint of the interval obtained by replacing in the original interval the infinite bound γ by sign(γ) × realmax; like the choice advocated by Motion 37, it preserves the monotonicity properties:

A downside of that approach is that it is more computationally demanding.

Even after having set the definition of the midpoint operator, there is still a lot of work to do to ensure the reliability of the actual implementation. It is not up to the IEEE P1788 standard to prescribe an implementation, as that would restrict the possibilities of implementers to take advantage of particular features of some platforms. Even after its release and widespread acceptance, there might then still be some room for some implementations to get it wrong, an easy thing to do, as this paper has shown.

A:24 F. Goualard made insightful comments on an early draft of the paper. I am indebted to the ACM TOMS reviewers of this paper, whose thoroughtness, expertise and dedication helped to greatly improve its contents and presentation. All remaining errors are, naturally, my own.


est destinée au dépôt

  Model calculation of the magnetic induction generated by magnetic nanoparticles flowing into a microfluidic system. Performance analysis of the detection Kamel Fodil, Matthieu Denoual, Christophe Dolabdjian, A. Treizebre, Vincent Senez To cite this version: Kamel Fodil, Matthieu Denoual, Christophe Dolabdjian, A. Treizebre, Vincent Senez. Model calculation of the magnetic induction generated by magnetic nanoparticles flowing into a microfluidic system. Performance analysis of the detection. IEEE Transactions on Magnetics, 2014, 50, pp.1-8. 10.1109/TMAG.2013.2291546. hal-00980048