Exploiting Binary Floating-Point Representations for Constraint Propagation

Roberto Bagnara

BUGSENG srl and Dept. of Mathematics and Computer Science, University of Parma, Italy,

                          bagnara@cs.unipr.it, http://www.cs.unipr.it/~bagnara

Matthieu Carlier

INRIA Rennes Bretagne Atlantique, France

Roberta Gori

Dept. of Computer Science, University of Pisa, Italy,

                             gori@di.unipi.it, http://www.di.unipi.it/~gori

Arnaud Gotlieb

Certus Software V&V Center, SIMULA Research Laboratory, Norway,

                            arnaud@simula.no, http://simula.no/people/arnaud

Floating-point computations are quickly finding their way in the design of safety- and mission-critical systems, despite the fact that designing floating-point algorithms is significantly more difficult than designing integer algorithms. For this reason, verification and validation of floating-point computations is a hot research topic. An important verification technique, especially in some industrial sectors, is testing. However, generating test data for floating-point intensive programs proved to be a challenging problem. Existing approaches usually resort to random or search-based test data generation, but without symbolic reasoning it is almost impossible to generate test inputs that execute complex paths controlled by floating-point computations. Moreover, as constraint solvers over the reals or the rationals do not natively support the handling of rounding errors, the need arises for efficient constraint solvers over floating-point domains. In this paper, we present and fully justify improved algorithms for the propagation of arithmetic IEEE 754 binary floating-point constraints. The key point of these algorithms is a generalization of an idea by B. Marre and C. Michel that exploits a property of the representation of floating-point numbers.

Key words : software verification; testing; floating-point numbers; constraint solving

1. Introduction

During the last decade, the use of floating-point computations in the design of critical systems has become increasingly acceptable. Even in the civil and military avionics domain, which are among the most critical domains for software, floating-point numbers are now seen as a sufficiently-safe, faster and cheaper alternative to fixed-point arithmetic. To the point that, in modern avionics, floating-point is the norm rather than the exception (Burdy et al. 2012).

Acceptance of floating-point computations in the design of critical systems took a long time. In fact, rounding errors can cause subtle bugs which are often missed by non experts (Monniaux

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

2  00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS

2008), and can lead to catastrophic failures. For instance, during the first Persian Gulf War, the failure of a Patriot missile battery in Dhahran was traced to an accumulating rounding error in the continuous execution of tracking and guidance software: this failure prevented the interception of an Iraqi Scud that hit the barracks in Dhahran, Saudi Arabia, killing 28 US soldiers (Skeel 1992). A careful analysis of this failure revealed that, even though the rounding error obtained at each step of the floating-point computation was very small, the propagation during a long loop-iterating path could lead to dramatic imprecision.

Adoption of floating-point computations in critical systems involves the use of thorough unit testing procedures that are able to exercise complex chains of floating-point operations. In particular, a popular practice among software engineers in charge of the testing of floating-point-intensive computations consists in executing carefully chosen loop-iterating paths in programs. They usually pay more attention to the paths that are most likely to expose the system to unstable numerical computations.^{1} For critical systems, a complementary requirement is to demonstrate the infeasibility of selected paths, in order to convince a third-party certification authority that certain unsafe behaviors of the systems cannot be reached. As a consequence, software engineers face two difficult problems:

1. How to accurately predict the expected output of a given floating-point computation?^{2}

2. How to find a test input that is able to exercise a given path, the execution of which depends on the results of floating-point computations, or to guarantee that such a path is infeasible?

The first problem has been well addressed in the literature (Kuliamin 2010) through several techniques. Ammann and Knight (1988) report on a technique known as the data diversity approach, which uses multiple related program executions of a program to check their results. Metamorphic testing (Chan et al. 1998) generalizes this technique by using known numerical relations of the function implemented by a program to check the results of two or more executions. Goubault (2001) proposes using the abstract interpretation framework (Cousot and Cousot 1977) to estimate the deviation of the floating-point results with respect to an interpretation over the reals. Scott et al. (2007) propose using a probabilistic approach to estimate round-off error propagation. More recently, Tang et al. (2010) propose to exploit perturbation techniques to evaluate the stability of a numerical program. In addition to these approaches, it is possible to use a (partial) specification, a prototype or an old implementation in order to predict the results for a new implementation.

1 A computation can be called numerically stable if it can be proven not to magnify approximation errors. It can be called (potentially) unstable otherwise. 2 This is the well-known oracle problem (see Weyuker 1982).

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS  3

In contrast, the second problem received only little attention. Beyond the seminal work of Miller and Spooner (1976), who proposed to guide the search of floating-point inputs to execute a selected path, few approaches try to exactly reason over floating-point computations. The work of Miller and Spooner (1976) paved the way to the development of search-based test data generation techniques, which consist in searching test inputs by minimizing a cost function, evaluating the distance between the currently executed path and a targeted selected path (Korel 1990, Lakhotia et al. 2010a, McMinn 2004). Although these techniques enable quick and efficient coverage of testing criteria such as “all decisions,” they are unfortunately sensitive to the rounding errors incurred in the computation of the branch distance (Arcuri 2009). Moreover, search-based test data generation cannot be used to study path feasibility, i.e., to decide whether a possible execution path involving floating-point computations is feasible or not in the program. In addition, these techniques can be stuck in local minima without being able to provide a meaningful result (Arcuri 2009). An approach to tackle these problems combines program execution and symbolic reasoning (Godefroid et al. 2005), and requires solving constraints over floating-point numbers in order to generate test inputs that exercise a selected behavior of the program under test. However, solving floating-point constraints is hard and requires dedicated filtering algorithms (Michel 2002, Michel et al. 2001). According to our knowledge, this approach is currently implemented in four solvers only: ECLAIR^{3}, FPCS (Blanc et al. 2006), FPSE^{4} (Botella et al. 2006), and GATeL, a test data generator for Lustre programs (Marre and Blanc 2005). It is worth noticing that existing constraint solvers dedicated to continuous domains (such as, e.g., RealPaver (Granvilliers and Benhamou 2006), IBEX and Quimper (Chabert and Jaulin 2009) or ICOS (Lebbah 2009)) correctly handle real or rational computations, but they cannot preserve the solutions of constraints over floating-point computations in all cases (see Section 5 for more on this subject). “Surprising” properties of floating-point computations such as absorption and cancellation (Goldberg 1991) show that the rounding operations can severely compromise the preservation of the computation semantics between the reals and the floats.

Example 1. Consider the C functions f1 and f2:

1 float f 1 ( float x ) {  float f 2 ( float x ) {
2 float y = 1 . 0 e12F ;  float y = 1 . 0 e12F ;
3 i f ( x < 10000.0F)  i f ( x > 0 . 0F)
4 z = x + y ; z = x + y ;
5 i f ( z > y )  i f ( z = y )
6  …  …

3 http://bugseng.com/products/eclair
4 http://www.irisa.fr/celtique/carlier/fpse.html

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

4  00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS

For both functions, let us ask the question whether the paths traversing lines 2-3-4-5-6 are feasible. The condition that must be satisfied in order for a certain path to be traversed is called path condition. For f1, the path conditions x < 10000.0 and x + 1.0e12 > 1.0e12, which on the reals are equivalent to x ∈ (0, 10000) whereas on the floats they have no solution. Conversely, for f2 the path conditions are x > 0.0 and x + 1.0e12 = 1.0e12, which have no solutions on the reals but are satisfied by all IEEE 754 single precision floating-point numbers in the range (0, 32767.99 ⋯).

A promising approach to improve the filtering capabilities of constraints over floating-point variables consists in using some peculiar numerical properties of floating-point numbers. For linear constraints, this led to a relaxation technique where floating-point numbers and constraints are converted into constraints over the reals by using linear programming approaches (Belaid et al. 2012). For interval-based consistency approaches, Marre and Michel (2010) identified a property of the representation of floating-point numbers and proposed to exploit it in filtering algorithms for addition and subtraction constraints. Carlier and Gotlieb (2011) proposed a reformulation of the Marre-Michel property in terms of “filtering by maximum ULP” (Units in the Last Place) that is generalizable to multiplication and division constraints.

Bagnara et al. (2013) addressed the question of whether the Marre-Michel property can be useful for the automatic solution of realistic test input generation problems: they sketched (without proofs) a reformulation and correction of the filtering algorithm proposed in (Marre and Michel 2010), along with a uniform framework that generalizes the property identified by Marre and Michel to the case of multiplication and division. Most importantly, (Bagnara et al. 2013) presented the implementation of filtering by maximum ULP in FPSE and some of its critical design choices, and an experimental evaluation on constraint systems that have been extracted from programs engaging into intensive floating-point computations. These results show that the Marre-Michel property and its generalization are effective, practical properties for solving constraints over the floats with an acceptable overhead. The experiments reported in (Bagnara et al. 2013) showed that improvement of filtering procedures with these techniques brings speedups of the overall constraint solving process that can be substantial (we have observed up to an order of magnitude); in the cases where such techniques do not allow significant extra-pruning, the slowdowns are always very modest (up to a few percent on the overall solution time).

The present paper is, on the one hand, the theoretical counterpart of (Bagnara et al. 2013) in that all the results are thoroughly proved; on the other hand, this paper generalizes and extends (Bagnara et al. 2013) as far as the handling of subnormals and floating-point division are concerned. More precisely, the contributions of the paper are:

1. a uniform framework for filtering by maximum ULP is thoroughly defined and justified;

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS  5

2. the framework is general enough to encompass all floating-point arithmetic operations and subnormals (the latter are not treated in (Bagnara et al. 2013));

3. a second indirect projection by maximum ULP for division (not present in any previous work);

4. all algorithms only use floating-point machine arithmetic operations on the same formats used by the analyzed computations.

The plan of the paper is as follows. Next section presents the IEEE 754 standard of binary floating-point numbers and introduces the notions and notations used throughout the paper. Section 3 recalls the basic principles of interval-based consistency techniques over floating-point variables and constraints. Section 4 presents our generalization of the Marre-Michel property along with a precise definition and motivation of all the required algorithms. Section 5 discusses related work. Section 6 concludes. The most technical proofs are available in the Online Supplement to this paper (Bagnara et al. 2015).

2. Preliminaries

In this section we recall some preliminary concepts and introduce the used notation.

2.1. IEEE 754 This section recalls the arithmetic model specified by the IEEE 754 standard for binary floatingpoint arithmetic (IEEE Computer Society 2008). Note that, although the IEEE 754 standard also specifies formats and methods for decimal floating-point arithmetic, in this paper we only deal with binary floating-point arithmetic.

IEEE 754 binary floating-point formats are uniquely identified by quantities: p ∈ ℕ, the number of significant digits (precision); e_{max} ∈ ℕ, the maximum exponent; −e_{min} ∈ ℕ, the minimum exponent.^{5} The single precision format has p = 24 and e_{max} = 127, the double precision format has p = 53 and e_{max} = 1023 (IEEE 754 also defines extended precision formats). A finite, non-zero IEEE 754 floating-point number z has the form (−1)^{s}b_{1}.m × 2^{e} where s is the sign bit, b_{1} is the hidden bit, m is the (p − 1)-bit significand and the exponent e is also denoted by e_{z} or exp(z). Hence the number is positive when s = 0 and negative when s = 1. b_{1} is termed “hidden bit” because in the binary interchange format encodings it is not explicitly represented, its value being encoded in the exponent (IEEE Computer Society 2008).

Each format defines several classes of numbers: normal numbers, subnormal numbers, signed zeros, infinities and NaNs (Not a Number). The smallest positive normal floating-point number is f_{min}^{nor} = 1.0 ⋯ 0 × 2^{e_{min}} = 2^{e_{min}} and the largest is f_{max} = 1.1 ⋯ 1 × 2^{e_{max}} = 2^{e_{max}}(2 − 2^{1−p}); normal numbers have the hidden bit b_{1} = 1. The non-zero floating-point numbers whose absolute value

5 Note that, although the IEEE 754 formats have e_{min} = 1 − e_{max}, we never use this property and decided to keep the extra-generality, which might be useful to accommodate other formats.

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

6  00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS

is less than 2^{e_{min}} are called subnormals: they always have exponent equal to e_{min} and fewer than p significant digits as their hidden bit is b_{1} = 0. Every finite floating-point number is an integral multiple of the smallest subnormal f_{min} = 0.0 ⋯ 01 × 2^{e_{min}} = 2^{e_{min}+1−p}. There are two infinities, denoted by +∞ and −∞, and two signed zeros, denoted by +0 and −0: they allow some algebraic properties to be maintained (Goldberg 1991).^{6} NaNs are used to represent the results of invalid computations such as a division of two infinities or a subtraction of infinities with the same sign: they allow the program execution to continue without being halted by an exception.

IEEE 754 defines five rounding directions: toward negative infinity (roundTowardNegative or, briefly, down), toward positive infinity (roundTowardPositive, a.k.a. up), toward zero (roundTowardZero, a.k.a. chop) and toward the nearest representable value (a.k.a. near); the latter comes in two flavors that depend on different tie-break rules for numbers exactly halfway between two representable numbers: roundTiesToEven (a.k.a. tail-to-even) or roundTiesToAway (a.k.a. tail-to-away) in which values with even significand or values away from zero are preferred, respectively. This paper is only concerned with roundTiesToEven, which is, by far, the most widely used. The roundTiesToEven value of a real number x will be denoted by [x]_{n}.

The most important requirement of IEEE 754 arithmetic is the accuracy of floating-point computations: add, subtract, multiply, divide, square root, remainder, conversion and comparison operations must deliver to their destination the exact result rounded as per the rounding mode in effect and the format of the destination. It is said that these operations are “correctly rounded.”

The accuracy requirement of IEEE 754 can still surprise the average programmer: for example the single precision, round-to-nearest addition of 999999995904 and 10000 (both numbers can be exactly represented) gives 999999995904, i.e., the second operand is absorbed. The maximum error committed by representing a real number with a floating-point number under some rounding mode can be expressed in terms of the function ulp: ℝ → ℝ (Muller 2005). Its value on 1.0 is about 10^{−7} for the single precision format.

2.2. Notation The set of real numbers is denoted by ℝ while 𝔽_{p,e_{max}} denotes a subset of the binary floating-point numbers, defined from a given IEEE 754 format, which includes the infinities −∞ and +∞, the signed zeros +0 and −0, but neither subnormal numbers nor NaNs. Subnormals are introduced in the set 𝔽_{p,e_{max}}^{sub} = 𝔽_{p,e_{max}} ∪ { (−1)^{s}0.m × 2^{e_{min}} | s ∈ {0, 1}, m ≠ 0 }. In some cases, the exposition can be much simplified by allowing the e_{max} of 𝔽_{p,e_{max}} to be ∞, i.e., by considering an idealized set of floats where the exponent is unbounded. Among the advantages is the fact that subnormals in

6 Examples of such properties are \sqrt{1/z} = 1/\sqrt{z} and 1/(1/x) = x for x = ±∞.

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS  7

                              z = x ⊕ y                        z = x ⊖ y

\overline{z} = \overline{x} ⊕ \overline{y},  (direct)

                   \underline{z} = \underline{x} ⊕ \underline{y}

\overline{x} = mid(\overline{z}, \overline{z}^{+}) ⊖ \underline{y} (1^{st} indirect)
\underline{x} = mid(\underline{z}, \underline{z}^{−}) ⊖ \overline{y}
\overline{y} = mid(\overline{z}, \overline{z}^{+}) ⊖ \underline{x} (2^{nd} indirect)
\underline{y} = mid(\underline{z}, \underline{z}^{−}) ⊖ \overline{x}

\overline{z} = \overline{x} ⊖ \underline{y},  (direct)

                                                   \underline{z} = \underline{x} ⊖ \overline{y}

\overline{x} = mid(\overline{z}, \overline{z}^{+}) ⊕ \overline{y} (1^{st} indirect)
\underline{x} = mid(\underline{z}, \underline{z}^{−}) ⊕ \underline{y}
\overline{y} = \overline{x} ⊖ mid(\underline{z}, \underline{z}^{−}) (2^{nd} indirect)
\underline{y} = \underline{x} ⊖ mid(\overline{z}, \overline{z}^{+})

Figure 1 Formulas for direct/indirect projections of addition/subtraction

𝔽_{p,e_{max}}^{sub} can be represented as normal floating-point numbers in 𝔽_{p,∞}. Given a set of floating-point numbers 𝔽, 𝔽^{+} denotes the “non-negative” subset of 𝔽, i.e., with s = 0.

For a finite, non-zero floating-point number x, we will write even(x) (resp., odd(x)) to signify that the least significant digit of x’s significand is 0 (resp., 1).

When the format is clear from the context, a real decimal constant (such as 10^{12}) denotes the corresponding roundTiesToEven floating-point value (i.e., 999999995904 for 10^{12}).

Henceforth, for x ∈ ℝ, x^{+} (resp., x^{−}) denotes the smallest (resp., greatest) floating-point number strictly greater (resp., smaller) than x with respect to the considered IEEE 754 format. Of course, we have f_{max}^{+} = +∞ and (−f_{max})^{−} = −∞.

Binary arithmetic operations over the floats will be denoted by ⊕, ⊖, ⊗ and ⊘, corresponding to +, −, · and / over the reals, respectively. According to IEEE 754, they are defined, under roundTiesToEven, by

x ⊕ y = [x + y]_{n},  x ⊖ y = [x − y]_{n},
x ⊗ y = [x · y]_{n},  x ⊘ y = [x/y]_{n}.

As IEEE 754 floating-point numbers are closed under negation, we denote the negation of x ∈ 𝔽_{p,e_{max}}^{sub} simply by −x. Note that negation is a bijection. The symbol ⊙ denotes any of ⊕, ⊖, ⊗ or ⊘. A floating-point variable x is associated with an interval of possible floating-point values; we will write x ∈ [\underline{x}, \overline{x}], where \underline{x} and \overline{x} denote the smallest and greatest value of the interval, \underline{x} ≤ \overline{x} and either \underline{x} ≠ +0 or \overline{x} ≠ −0. Note that [+0, −0] is not an interval, whereas [−0, +0] is the interval denoting the set of floating-point numbers {−0, +0}.

3. Background on Constraint Solving over Floating-Point Variables

In this section, we briefly recall the basic principles of interval-based consistency techniques over floating-point variables and constraints.

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

8  00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS

3.1. Interval-based Consistency on Arithmetic Constraints Program analysis usually starts with the generation of an intermediate code representation in a form called three-address code (TAC). In this form, complex arithmetic expressions and assignments are decomposed into sequences of assignment instructions of the form

                        result := operand_{1} operator operand_{2}.

A further refinement consists in the computation of the static single assignment form (SSA) whereby, labeling each assigned variable with a fresh name, assignments can be considered as if they were equality constraints. For example, the TAC form of the floating-point assignment z := z∗z+z is t := z ∗ z; z := t + z, which in SSA form becomes t_{1} := z_{1} ∗ z_{1}; z_{2} := t_{1} + z_{1}, which, in turn, can be regarded as the conjunction of the constraints t_{1} = z_{1} ⊗ z_{1} and z_{2} = t_{1} ⊕ z_{1}.

In an interval-based consistency approach to constraint solving over the floats, constraints are used to iteratively narrow the intervals associated with each variable: this process is called filtering. A projection is a function that, given a constraint and the intervals associated with two of the variables occurring in it, computes a possibly refined interval^{7} for the third variable (the projection is said to be over the third variable). Taking z_{2} = t_{1} ⊕ z_{1} as an example, the projection over z_{2} is called direct projection (it goes in the same sense of the TAC assignment it comes from), while the projections over t_{1} and z_{1} are called indirect projections.^{8} Note that, for constraint propagation, both direct and indirect projections are applied in order to refine the intervals for t_{1}, z_{1} and z_{2}. In this paper we propose new filtering algorithms for improving indirect projections.

A projection is called optimal if the interval constraints it infers are as tight as possible, that is, if both bounds of the inferred intervals are attainable (and thus cannot be pruned).

Figure 1 gives non-optimal projections for addition and subtraction. For finite x, y ∈ 𝔽_{p,e_{max}}, mid(x, y) denotes the number that is exactly halfway between x and y; note that either mid(x, y) ∈ 𝔽_{p,e_{max}} or mid(x, y) ∈ 𝔽_{p+1,e_{max}}. Non-optimal projections for multiplication and division can be found in (Botella et al. 2006, Michel 2002). Optimal projections are known for monotonic functions over one argument (Michel 2002), but they are generally not available for other functions. Note, however, that optimality is not required in an interval-based consistency approach to constraint solving, as filtering is just used to remove some, not necessarily all, inconsistent values.

7 That is, tighter than the original interval. 8 Note that direct and indirect projections are idempotent, as their inputs and outputs do not intersect. Consider z = x ⊕ y: direct projection propagates information on x and y onto z, and doing it twice in a row would not enable any further inference. Likewise for indirect projections, which propagate information on z onto x and y.

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS  9

zα α^{+} α^{++} β β^{+} zΔ 2Δ 2Δ 0

Figure 2 An illustration of the Marre-Michel property: the segment z, if it has to represent the difference between two floats, cannot be moved past α

3.2. The Marre-Michel Property Marre and Michel (2010) published an idea to improve the filtering of the indirect projections for addition and subtraction. This is based on a property of the distribution of floating-point numbers among the reals: the greater a float, the greater the distance between it and its immediate successor. More precisely, for a given float x with exponent e_{x}, if x^{+} − x = Δ, then for y of exponent e_{x} + 1 we have y^{+} − y = 2Δ.

Proposition 1. (Marre and Michel 2010, Proposition 1) Let z ∈ 𝔽_{p,∞} be such that 0 < z < +∞; let also

z = 1.b_{2} ⋯ b_{i} \overbrace{0 ⋯ 0}^{k} × 2^{e_{z}},  with b_{i} = 1;
α = \overbrace{1.1 ⋯ 1}^{p} × 2^{e_{z}+k},  with k = p − i;
β = α ⊕ z.

Then, for each x, y ∈ 𝔽_{p,∞}, z = x ⊖ y implies that x ≤ β and y ≤ α. Moreover, β ⊖ α = β − α = z.

This property, which can be generalized to subnormals, can intuitively be explained on Figure 2 as follows. Let z ∈ 𝔽_{p,∞} be a strictly positive constant such that z = x ⊖ y, where x, y ∈ 𝔽_{p,∞} are unknown. The Marre-Michel property says that y cannot be greater than α. In fact, α is carefully positioned so that α^{++} − α^{+} = 2(α^{+} − α), e_{α} + 1 = e_{β} and z = β − α; if we take y = α^{+} we need x > β if we want z = x − y; however, the smallest element of 𝔽_{p,∞} that is greater than β, β^{+}, is 2Δ away from β, i.e., too much. Going further with y does not help: if we take y ≥ α^{+}, then y − α is an odd multiple of Δ (one Δ step from α to α^{+}, all the subsequent steps being even multiples of Δ), whereas for each x ≥ β, x − β is an even multiple of Δ. Hence, if y > α, |z − (x − y)| ≥ Δ = 2^{e_{z}+1−i}. However, since k ≠ p − 1, z^{+} − z = z − z^{−} = 2^{e_{z}+1−p} ≤ Δ. The last inequality, which holds because p ≥ i, implies z ≠ x ⊖ y. A similar reasoning allows one to see that x cannot be greater than β independently from the value of y. In order to improve the filtering of the addition/subtraction projectors, Marre and Michel (2010) presented an algorithm to maximize the values of α and β over an interval. As we will see, that algorithm is not correct for some inputs. In Section 4.5, the main ideas behind the work presented in (Marre and Michel 2010) will be revisited, corrected and discussed.

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

10  00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS

4. Filtering by Maximum ULP

In this section we first informally present, by means of worked numerical examples, the techniques that are precisely defined later. We then reformulate the Marre-Michel property so as to generalize it to subnormals and to multiplication and division operators. The filtering algorithms that result from this generalization are collectively called filtering by maximum ULP.

4.1. Motivating Example Consider the IEEE 754 single-precision constraint z = x ⊕ y with initial intervals z ∈ [−∞, +∞], x ∈ [−1.0 × 2^{50}, 1.0 × 2^{50}] and y ∈ [−1.0 × 2^{30}, 1.0 × 2^{30}]. Forward projection gives

z ∈ [−1. \overbrace{0 ⋯ 0}^{19} 1 × 2^{50}, 1. \overbrace{0 ⋯ 0}^{19} 1 × 2^{50}],

which is optimal, as both bounds are attainable. Suppose now the interval for z is further restricted to z ∈ [1.0, 2.0] due to, say, a constraint from an if-then-else in the program or another indirect projection.

With the classical indirect projection we obtain x, y ∈ [−1.0 × 2^{30}, 1.0 × 2^{30}], which, however, is not optimal. For example, pick x = 1.0 × 2^{30}: for y = −1.0 × 2^{30} we have x ⊕ y = 0 and x ⊕ y^{+} = 64. By monotonicity of ⊕, for no y ∈ [−1.0 × 2^{30}, 1.0 × 2^{30}] we can have x ⊕ y ∈ [1.0, 2.0].

With our indirect projection, fully explained later, we obtain, from z ∈ [1.0, 2.0], the much tighter intervals x, y ∈ [−1.1 ⋯ 1 × 2^{24}, 1.0 × 2^{25}]. These are actually optimal as −1.1 ⋯ 1 × 2^{24} ⊕ 1.0 × 2^{25} = 1.0 × 2^{25} ⊕ −1.1 ⋯ 1 × 2^{24} = 2.0. This example shows that filtering by maximum ULP can be stronger than classical interval-consistency based filtering. However, the opposite phenomenon is also possible. Consider again z = x ⊕ y with z ∈ [1.0, 2.0]. Suppose now the constraints for x and y are x ∈ [1.0, 5.0] and y ∈ [−f_{max}, f_{max}]. As we have seen, our indirect projection gives y ∈ [−1.1 ⋯ 1 × 2^{24}, 1.0 × 2^{25}]; in contrast, the classical indirect projection exploits the available information on x to obtain y ∈ [−4, 1]. Indeed, classical and maximum ULP filtering for addition and subtraction are orthogonal: both should be applied in order to obtain precise results.

For an example on multiplication, consider the IEEE 754 single-precision constraint z = x⊗y with initial intervals z ∈ [1.0×2^{−50}, 1.0×2^{−30}] and x, y ∈ [−∞, +∞]. In this case, classical projections do not allow pruning the intervals. However, take x = 1.1×2^{119}: for y = 0 we have x⊗y = 0 and x⊗y^{+} = 1.1 × 2^{−30}. By monotonicity of ⊗, for no y ∈ [−∞, +∞] we can have x ⊗ y ∈ [1.0 × 2^{−50}, 1.0 × 2^{−30}].

On the same example, x, y ∈ [−1.0 ⋯ 0 × 2^{119}, 1.0 ⋯ 0 × 2^{119}] are the constraints inferred by our indirect projection. These are optimal because 1.0 × 2^{−30} = −1.0 ⋯ 0 × 2^{119} ⊗ −1.0 ⋯ 0 × 2^{149} = 1.0 ⋯ 0 × 2^{119} ⊗ 1.0 ⋯ 0 × 2^{149}. As is the case for addition, classical indirect projection can be more precise. Consider again z = x ⊗ y with z ∈ [1.0 × 2^{−50}, 1.0 × 2^{−30}], x ∈ [2.0, 4.0] and y ∈ [−f_{max}, f_{max}]. Classical indirect projection infers y ∈ [1.0 × 2^{−52}, 1.0 × 2^{−31}] by exploiting the information on x.

z^{−} zz^{+} z^{−} zz^{+} z^{−} − Δ_{z^{−}}^{−} /2 z^{−} + Δ_{z^{−}}^{+} /2 z − Δ_{z}^{−} /2 z + Δ_{z}^{+} /2 z^{+} − Δ_{z^{+}}^{−} /2 z^{+} + Δ_{z^{+}}^{+} /2 ℝ 𝔽

Figure 3 Rounding of real numbers in the neighborhood of an even floating-point number z under roundTiesToEven

4.2. Round-To-Nearest Tail-To-Even We now formally define the roundTiesToEven rounding mode. To do that, we first introduce two functions: Δ_{z}^{+} and Δ_{z}^{−} give the distance between z^{+} and z and the distance between z and z^{−}.

Definition 1. The partial functions Δ^{−}: 𝔽_{p,e_{max}}^{sub} ↣ ℝ and Δ^{+} : 𝔽_{p,e_{max}}^{sub} ↣ ℝ are defined as follows, for each finite z ∈ 𝔽_{p,e_{max}}^{sub}:

Δ_{z}^{+} = \begin{cases}2^{1−p+e_{max}}, & if z = f_{max}; \\ f_{min}, & if z = +0 or z = −0; \\ z^{+} − z, & otherwise; \end{cases}

Δ_{z}^{−} = \begin{cases}2^{1−p+e_{max}} & if z = −f_{max}; \\ f_{min}, & if z = +0 or z = −0; \\ z − z^{−}, & otherwise. \end{cases}

Note the special cases when z = ±0: since both +0 and −0 represent the real number 0, the distance between z^{+} = f_{min} and z = ±0 is f_{min}. We can now define the function [·]_{n} that captures roundTiesToEven.

Definition 2. For x ∈ ℝ, [x]_{n} is defined as follows:

[x]_{n} = \begin{cases}+0, & if 0 ≤ x ≤ Δ_{0}^{+} /2; \\ −0, & if − Δ_{0}^{−} /2 ≤ x < 0; \\ z, & if z ∈ 𝔽_{p,e_{max}}^{sub} \ \{−∞, +∞\} and either even(z) and \\ & z − Δ_{z}^{−} /2 ≤ x ≤ z + Δ_{z}^{+} /2, or odd(z) and \\ & z − Δ_{z}^{−} /2 < x < z + Δ_{z}^{+} /2; \\ +∞, & if x ≥ f_{max} + Δ_{f_{max}}^{+} /2; \\ −∞, & if x ≤ −f_{max} − Δ_{−f_{max}}^{−} /2. \end{cases}

Figure 3 illustrates the roundTiesToEven rounding mode; if z is even, each real number between z − Δ_{z}^{−} /2 and z + Δ_{z}^{+} /2, including extremes, is rounded to the same floating-point number z. As z is even, z^{−} is odd, and each real number between z^{−} − Δ_{z^{−}}^{−} /2 and z^{−} + Δ_{z^{−}}^{+} /2, excluding extremes, is rounded to z^{−}. Similarly for z^{+}. Note that point z − Δ_{z}^{−} /2 coincides with z^{−} + Δ_{z^{−}}^{+} /2 and z + Δ_{z}^{+} /2 coincides with z^{+} − Δ_{z^{+}}^{−} /2.

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

12  00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS

All rounding modes are monotonic; in particular, for each x, y ∈ ℝ, x ≤ y implies [x]_{n} ≤ [y]_{n}. Moreover, the chop and near rounding modes are symmetric, i.e., the value after rounding does not depend on the sign: for each x ∈ ℝ, [x]_{n} = −[−x]_{n}.

4.3. Upper Bound It is worth pointing out that, while arithmetic operations on reals are strictly monotone, that is if x+y = z then x_{1} +y > z for any x_{1} > x, in floating-point arithmetic, operations are just monotone. If x + y = z then we may still have x_{1} + y = z for some (or many) x_{1} > x since addition over the floats is absorbing. Therefore, for determining the greatest (or the smallest) x_{1} satisfying x_{1} +y = z and correctly filter intervals over the floats, we need to introduce an appropriate, duly justified, framework.

For each IEEE 754 floating-point operation ⊙ ∈ {⊕, ⊖, ⊗, ⊘}, in later sections we will define the sets 𝔽_{⊙} ⊆ 𝔽_{p,e_{max}}^{sub} and 𝔽̄_{⊙} ⊆ 𝔽_{p,∞}. Then we will define functions δ̄_{⊙}: 𝔽_{⊙} → 𝔽̄_{⊙} (see Definition 3 in Section 4.5 for ⊕ and, consequently, ⊖, Definition 5 in Section 4.6 for ⊗, and Definition 6 in Section 4.7 for ⊘) that satisfy the following property, for each z ∈ 𝔽_{⊙} \ {−0, +0, −∞}:

δ̄_{⊙}(z) = max{ v ∈ 𝔽̄_{⊙} | ∃y ∈ 𝔽̄_{⊙} . v ⊙ y = z }.  (1)

In words, δ̄_{⊙}(z) is the greatest float in 𝔽̄_{⊙} that can be the left operand of ⊙ to obtain z.

Verifying that a function δ̄_{⊙} satisfies (1) is equivalent to proving that it satisfies the following properties, for each z ∈ 𝔽_{⊙} \ {−0, +0, −∞}:

δ̄_{⊙}(z) ∈ { v ∈ 𝔽̄_{⊙} | ∃y ∈ 𝔽̄_{⊙} . v ⊙ y = z };  (2)
∀z′ ∈ 𝔽̄_{⊙}: z′ > δ̄_{⊙}(z) ⟹ z′ ∉ { v ∈ 𝔽̄_{⊙} | ∃y ∈ 𝔽̄_{⊙} . v ⊙ y = z }.  (3)

Property (3) means δ̄_{⊙}(z) is a correct upper bound for the possible values of x, whereas (2) implies that δ̄_{⊙}(z) is the most precise upper bound we could choose.

Note that we may have 𝔽̄_{⊙} ⊈ 𝔽_{p,e_{max}}: property (1) refers to an idealized set of floating-point numbers with unbounded exponents.

Since we are interested in finding the upper bound of δ̄_{⊙}(z) for z ∈ [\underline{z}, \overline{z}], we need the following

Proposition 2. Let w, v_{1}, …, v_{n} ∈ 𝔽_{⊙} \ {−0, +0, −∞} be such that, for each i = 1, …, n, δ̄_{⊙}(w) ≥ δ̄_{⊙}(v_{i}). Then, for each w′ ∈ 𝔽̄_{⊙} with w′ > δ̄_{⊙}(w) and each z ∈ 𝔽_{⊙} \ {−0, +0, −∞}, we have that w′ ∉ { v ∈ 𝔽̄_{⊙} | ∃y ∈ 𝔽̄_{⊙} . v ⊙ y = z }.

Proof. Follows directly from (1).

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS  13

Let z = x ⊙ y be a floating-point constraint where −0, +0, −∞ ∉ [\underline{z}, \overline{z}] and let w ∈ [\underline{z}, \overline{z}] be such that δ̄_{⊙}(w) ≥ δ̄_{⊙}(v) for each v ∈ [\underline{z}, \overline{z}]: then no element of x that is greater than δ̄_{⊙}(w) can participate to a solution of the constraint.

Dually, in order to refine the upper bound of y subject to z = x ⊙ y, it is possible to define a function δ̄′_{⊙} satisfying the following property, for each z ∈ 𝔽_{⊙} \ {−0, +0, −∞}:

δ̄′_{⊙}(z) = max{ v ∈ 𝔽̄_{⊙} | ∃x ∈ 𝔽̄_{⊙} . x ⊙ v = z }.  (4)

Due to (4), a result analogous to the one of Proposition 2 holds for δ̄′_{⊙}, which allows refining the interval for y. Note, though, that when ⊙ is commutative (i.e., it is ⊕ or ⊗), δ̄_{⊙} = δ̄′_{⊙}.

4.4. Lower bound For computing the lower bound, we will introduce functions δ̄_{⊙}: 𝔽_{⊙} → 𝔽̄_{⊙} (defined in terms of the corresponding δ̄_{⊙} functions in Section 4.5 for ⊕ and ⊖, in Section 4.6 for ⊗, and in Section 4.7 for ⊘) satisfying the following property, for each z ∈ 𝔽_{⊙} \ {−0, +0, +∞}:

δ̄_{⊙}(z) = min{ v ∈ 𝔽̄_{⊙} | ∃y ∈ 𝔽̄_{⊙} . v ⊙ y = z }.  (5)

This property entails a result similar to Proposition 2: given constraint z = x ⊙ y where −0, +0, +∞ ∉ [\underline{z}, \overline{z}] and w ∈ [\underline{z}, \overline{z}] such that δ̄_{⊙}(w) ≤ δ̄_{⊙}(v) for each v ∈ [\underline{z}, \overline{z}], the float δ̄_{⊙}(w) is a possibly refined lower bound for x.

In a dual way, in order to refine the lower bound of y subject to z = x⊙y, we will define functions δ̄′_{⊙} satisfying, for each z ∈ 𝔽_{⊙} \ {−0, +0, +∞}:

δ̄′_{⊙}(z) = min{ v ∈ 𝔽̄_{⊙} | ∃x ∈ 𝔽̄_{⊙} . x ⊙ v = z }.  (6)

Property (6) ensures that, under z = x ⊙ y where −0, +0, +∞ ∉ [\underline{z}, \overline{z}], if w ∈ [\underline{z}, \overline{z}] is such that δ̄′_{⊙}(w) ≤ δ̄′_{⊙}(v) for each v ∈ [\underline{z}, \overline{z}], then the float δ̄′_{⊙}(w) is a possibly refined lower bound for y.

Again, when ⊙ is commutative δ̄_{⊙} = δ̄′_{⊙}.

4.5. Filtering by Maximum ULP on Addition/Subtraction In this section we introduce the functions δ̄_{⊕}, δ̄_{⊕}, δ̄_{⊖}, δ̄′_{⊖}, δ̄_{⊖} and δ̄′_{⊖}. Note that, since ⊕ is commutative, we have δ̄′_{⊕} = δ̄_{⊕} and δ̄′_{⊕} = δ̄_{⊕}. Moreover, the function δ̄_{⊕}: 𝔽_{⊕} → 𝔽̄_{⊕} can be defined in terms of the function δ̄_{⊕} as follows: for each z ∈ 𝔽_{⊕} \ {−0, +0, +∞}, δ̄_{⊕}(z) = − δ̄_{⊕}(−z). We see that, if δ̄_{⊕} satisfies Property (1), then δ̄_{⊕} satisfies Property (5). Again, since ⊕ is commutative, δ̄′_{⊕} = δ̄_{⊕}.

The first step for defining δ̄_{⊕} consists in extending Proposition 1 in order to explicitly handle subnormal numbers. Such extension was already sketched by Marre and Michel (2010): here we fully describe it and prove its correctness. Subnormals, which in 𝔽_{p,e_{max}}^{sub} are represented by numbers

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

14  00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS

having the hidden bit b_{1} = 0 and exponent e_{min}, can be represented in 𝔽_{p,∞} by numbers with b_{1} = 1 and exponent strictly smaller than e_{min}. Namely, the element of 𝔽_{p,e_{max}}^{sub}

0.0 ⋯ 01b_{j+1} ⋯ b_{p} × 2^{e_{min}}

can be represented in 𝔽_{p,∞} by the (normal) float

1.b_{j+1} ⋯ b_{p} \overbrace{0 ⋯ 0}^{j−1} × 2^{e_{min}−(j−1)}.

Based on this observation we can state the following

Proposition 3. Let z ∈ 𝔽_{p,e_{min}}^{sub} be such that 0 < z < f_{min}^{nor}; define also

z = 0.0 ⋯ 01b_{j+1} ⋯ b_{i} \overbrace{0 ⋯ 0}^{k} × 2^{e_{min}},  with b_{i} = 1;
α = \overbrace{1.1 ⋯ 1}^{p} × 2^{e_{min}+k},  with k = p − i;
β = α ⊕ z.

Then, for each x, y ∈ 𝔽_{p,e_{max}}^{sub}, z = x ⊖ y implies that x ≤ β and y ≤ α. Moreover, β ⊖ α = β − α = z.

Proof. The subnormal z is represented in 𝔽_{p,∞} by the normal float

ẑ = 1.b_{j+1} ⋯ b_{i} \overbrace{0 ⋯ 0}^{k} \overbrace{0 ⋯ 0}^{j−1} × 2^{e_{min}−(j−1)} = 1.b_{j+1} ⋯ b_{i} \overbrace{0 ⋯ 0}^{k+j−1} × 2^{e_{min}−(j−1)}.

We can apply Proposition 1 to ẑ and obtain α = 1.1 ⋯ 1 × 2^{e_{min}−(j−1)+k+j−1} = 1.1 ⋯ 1 × 2^{e_{min}+k}. Moreover, Proposition 1 assures that

β = α ⊕ 1.b_{j+1} ⋯ b_{i} \overbrace{0 ⋯ 0}^{k+j−1} × 2^{e_{min}−(j−1)}

is such that, for each x, y ∈ 𝔽_{p,∞}, z = x⊖y implies x ≤ β and y ≤ α and β ⊖α = β −α = z. Since each number in 𝔽_{p,e_{max}}^{sub} has an equivalent representation in 𝔽_{p,∞}, we only need to prove that β = α ⊕ z, which holds, since

β = α ⊕ 1.b_{j+1} ⋯ b_{i} \overbrace{0 ⋯ 0}^{k+j−1} × 2^{e_{min}−(j−1)}

= α ⊕ 0.0 ⋯ 01b_{j+1} ⋯ b_{i}\underbrace{0 ⋯ 0}_{k} × 2^{e_{min}}

= α ⊕ z.

□

Using Propositions 1 and 3, we formally define the function δ̄_{⊕} as follows.

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS  15

Definition 3. Let 𝔽_{⊕} = 𝔽_{p,e_{max}}^{sub}, 𝔽̄_{⊕} = 𝔽_{p,∞}^{+}, and z ∈ 𝔽_{⊕} be such that |z| = b_{1}.b_{2} ⋯ b_{i}0 ⋯ 0 × 2^{e_{z}}, with b_{i} = 1. Similarly to Propositions 1 and 3, let k = p − i, α = 1.1 ⋯ 1 × 2^{e_{z}+k} and β = α ⊕ |z|. Then δ̄_{⊕}: 𝔽_{⊕} → 𝔽̄_{⊕} is defined, for each z ∈ 𝔽_{⊕}, by

δ̄_{⊕}(z) = \begin{cases}+∞, & if z = −∞ or z = +∞; \\ α, & if −∞ < z < 0; \\ +0, & if z = −0 or z = +0; \\ β, & if 0 < z < +∞. \end{cases}

Theorem 1. Function δ̄_{⊕} is well-defined and satisfies (2) and (3).

Proof. We first show that δ̄_{⊕}(z) is well-defined, i.e., that it is a total function from 𝔽_{p,e_{max}}^{sub} to 𝔽_{p,∞}^{+}. To this aim note that α and β are always non-negative normal floating-point numbers belonging to 𝔽_{p,∞}, and that δ̄_{⊕}(z) is defined for each z ∈ 𝔽_{p,e_{max}}^{sub}. Secondly, let us consider the following cases:

z = +∞: for each y ≠ −∞ we have +∞ ⊕ y = +∞; thus, as δ̄_{⊕}(z) = +∞, (2) holds and (3) vacuously holds.

f_{min}^{nor} ≤ z < +∞: we can apply Proposition 1 to obtain z = β ⊖α. Then note that β ⊖α = [β −α]_{n} = [β + −α]_{n} = β ⊕ −α. Hence, β ⊕ −α = z. Thus, δ̄_{⊕}(z) ⊕ −α = β ⊕ −α = z and (2) is satisfied with y = −α. For proving (3), first note that β > −α since β > 0 and α > 0. Moreover, by Proposition 1, we know that there does not exist an x ∈ 𝔽_{p,∞} with x > β such that there exists y ∈ 𝔽_{p,∞} that satisfies x ⊖ y = z. Since x ⊖ y = x ⊕ −y we can conclude that, for each z′ > β = δ̄_{⊕}(z), it does not exist y′ ∈ 𝔽_{p,∞} such that z′ ⊕ y′ = z. Hence also (3) holds.

0 < z < f_{min}^{nor}: by applying Proposition 3 instead of Proposition 1 we can reason exactly as in the previous case.

−∞ < z ≤ −f_{min}^{nor}: since 0 < −z < +∞ we can apply Proposition 1 to −z and obtain β ⊖ α = −z and thus −(β ⊖ α) = z. As [·]_{n} is a symmetric rounding mode, we have −(β ⊖ α) = −[β − α]_{n} = [α − β]_{n} = α ⊕ −β = z. Thus, δ̄_{⊕}(z) ⊕ −β = α ⊕ −β = z and (2) is satisfied with y = −β. For proving (3), first note that α > −β since α > 0 and β > 0. Moreover, by Proposition 1, we know that there does not exist an y ∈ 𝔽_{p,∞} with y > α such that there exists x ∈ 𝔽_{p,∞} that satisfies x⊖y = −z. Since x ⊖ y = −z is equivalent to y ⊕ −x = z, we can conclude that, for each z′ > α = δ̄_{⊕}(z), it does not exist y′ ∈ 𝔽_{p,∞} such that z′ ⊕ y′ = z. Therefore, also in this case, (3) holds.

−f_{min}^{nor} < z < 0: by applying Proposition 3 instead of Proposition 1 we can reason exactly as in the previous case. □

As we have already observed, since ⊕ is commutative we have δ̄′_{⊕} = δ̄_{⊕}, that is, the same function δ̄_{⊕} is used to filter both x and y with respect to z = x ⊕ y.

We now need algorithms to maximize δ̄_{⊕} and minimize δ̄_{⊕} over an interval of floating-point values. Since the two problems are dual to each other, we will focus on the maximization of δ̄_{⊕}.

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

16  00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS

As δ̄_{⊕} is not monotonic, a nontrivial analysis of its range over an interval is required. When the interval contains only finite, nonzero and positive (resp., negative) values, the range of δ̄_{⊕} has a simple shape. We are thus brought to consider an interval [\underline{z}, \overline{z}] such that \underline{z}∉ \{−∞, −0, +0\}, \overline{z}∉ \{−0, +0, +∞\}, and \underline{z} and \overline{z} have the same sign. We will now revisit, correct and extend to subnormal floating-point numbers the algorithm originally proposed by Marre and Michel (2010) to maximize δ̄_{⊕} over [\underline{z}, \overline{z}].

The idea presented in (Marre and Michel 2010) is the following. When dealing with an interval [\underline{z}, \overline{z}] with \underline{z} > 0, α (and thus β and, therefore, our δ̄_{⊕}) grows (i) with the exponent and (ii) with the number of successive 0 bits to the right of the significand, i.e., k in Propositions 1 and 3 and in Definition 3. Thus, maximizing these two criteria allows one to maximize α over the interval.

Definition 4. Let z be a variable over 𝔽_{p,e_{max}}^{sub}. If we have 0 < \underline{z} < \overline{z} < +∞, then μ_{⊕}(z) ∈ [\underline{z}, \overline{z}] is given by:

1. μ_{⊕}(z) = 1.0 ⋯ 0 × 2^{e\overline{z}}, if e\underline{z} ≠ e\overline{z};

2. μ_{⊕}(z) = b_{1}.b_{2} ⋯ b_{i−1}a0 ⋯ 0 × 2^{e\overline{z}}, if e\underline{z} = e\overline{z}, where, for some b_{i} ≠ b′_{i}:

\underline{z} = b_{1}.b_{2} ⋯ b_{i−1}b_{i} ⋯ × 2^{e\overline{z}};
\overline{z} = b_{1}.b_{2} ⋯ b_{i−1}b′_{i} ⋯ × 2^{e\overline{z}};
a = \begin{cases}0, & if b_{1}.b_{2} ⋯ b_{i−1}0 ⋯ 0 × 2^{e\overline{z}} = \underline{z}; \\ 1, & otherwise. \end{cases}

If 0 < \underline{z} = \overline{z} < +∞, then μ_{⊕}(z) = \underline{z}. If −∞ < \underline{z} ≤ \overline{z} < 0, then μ_{⊕}(z) ∈ [\underline{z}, \overline{z}] is simply defined by μ_{⊕}(z) = − μ_{⊕}(w) where w ∈ [−\overline{z}, −\underline{z}]. We leave μ_{⊕}(z) undefined otherwise. Note that Definition 4 cannot be usefully extended to intervals containing zeros or infinities, as no interesting bounds can be derived for x and y in such cases. Consider, for example, the constraint x ⊕ y = z with z = +0: for each x ∈ [−f_{max}, +f_{max}] we have x ⊕ −x = +0. Hence, when the interval of z contains zeros or infinities, only the classical filtering (Botella et al. 2006, Michel 2002) is applied.

Theorem 2. Let z be over 𝔽_{p,e_{max}}^{sub} with \underline{z}∉ \{−∞, −0, +0\} and \overline{z}∉ \{−0, +0, +∞\} having the same sign. Then, for each z ∈ [\underline{z}, \overline{z}], δ̄_{⊕}(z) ≤ δ̄_{⊕}(μ_{⊕}(z)).

Proof. Without loss of generality, assume \underline{z} > 0. If \underline{z} = \overline{z} the result holds. Let us now assume \underline{z} < \overline{z}. We start proving that α and β of Definition 3 computed over μ_{⊕}(z) are greater than or equal to the α’s and β’s computed over any other value in [\underline{z}, \overline{z}].

We first prove that μ_{⊕}(z) maximizes α. For z ∈ [\underline{z}, \overline{z}] we have

α = 1.1 ⋯ 1 × 2^{e_{z}+k},

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS  17

where k is the number of successive 0’s to the right of the significand of z. Let us consider the maximum exponent of the values in z, which is e\overline{z}. Among the values in [\underline{z}, \overline{z}] with such an exponent, we want to select the one with the highest number of successive zeros to the right of the significand. Since \underline{z} > 0, the maximum value for α would be attained by the float 1.0 ⋯ 0 × 2^{e\overline{z}}, if this belongs to [\underline{z}, \overline{z}]. This happens in three cases:

1. e\underline{z} ≠ e\overline{z} and μ_{⊕}(z) = 1.0 ⋯ 0 × 2^{e\overline{z}}, by the first case of Definition 4.

2. e\underline{z} = e\overline{z} and \underline{z} = 1.0 ⋯ 0 × 2^{e\overline{z}}; in this case we have, again, μ_{⊕}(z) = 1.0 ⋯ 0 × 2^{e\overline{z}}, so defined by the second case of Definition 4; in fact, for some i ∈ \{2, …, p − 1\} that depends on \overline{z}, we have

\overline{z} = 1.b_{2} ⋯ b_{i−1}10 ⋯ 0 × 2^{e\overline{z}},
\underline{z} = 1.b_{2} ⋯ b_{i−1}00 ⋯ 0 × 2^{e\overline{z}}

with b_{2} = ⋯ = b_{i−1} = 0, and the algorithm gives 1.b_{2} ⋯ b_{i−1}a0 ⋯ 0×2^{e\overline{z}} with a = 0, i.e., 1.0 ⋯ 0×2^{e\overline{z}}.

3. e\underline{z} = e\overline{z}, \underline{z} = 0.b_{2} ⋯ b_{p} × 2^{e_{min}} and \overline{z} = 1.b′_{2} ⋯ b′_{p} × 2^{e_{min}}; thus we have, μ_{⊕}(z) = 1.0 ⋯ 0 × 2^{e_{min}}, once again by the second case of Definition 4 where i = 1, hence μ_{⊕}(z) = a.0 ⋯ 0 × 2^{e_{min}}. Moreover, since \underline{z} > 0, necessarily \underline{z} ≠ 0.0 ⋯ 0 × 2^{e_{min}} and we must have a = 1. We are now left with the case when 1.0 ⋯ 0 × 2^{e\overline{z}} ∉ [\underline{z}, \overline{z}]. This occurs when e\underline{z} = e\overline{z} but either \underline{z} > 1.0 ⋯ 0 × 2^{e\overline{z}} or \overline{z} < 1.0 ⋯ 0 × 2^{e\overline{z}}. In both cases, all the floats in [\underline{z}, \overline{z}] have the same exponent and the same most significant bit (b_{1}). Therefore, in order to maximize α, we need to choose among them the one with the greatest number of successive zeros to the right of the significand. The first step is to find the index of the most significant significand bit where \underline{z} and \overline{z} differ: since \underline{z} < \overline{z}, such an index must exist. Let then

\underline{z} = b_{1}.b_{2} ⋯ b_{i−1}b_{i} ⋯ × 2^{e\overline{z}},
\overline{z} = b_{1}.b_{2} ⋯ b_{i−1}b′_{i} ⋯ × 2^{e\overline{z}},

where b_{i} = 0 and b′_{i} = 1 for some i > 1. The significand maximizing α is b_{1}.b_{2} ⋯ b_{i−1}0 ⋯ 0. Indeed, any float having a significand with a larger number of consecutive zeros to the right does not belong to [\underline{z}, \overline{z}]. However, it is not always the case that b_{1}.b_{2} ⋯ b_{i−1}0 ⋯ 0 × 2^{e\overline{z}} belongs to [\underline{z}, \overline{z}]: we must have

\underline{z} = b_{1}.b_{2} ⋯ b_{i−1}b_{i}0 ⋯ 0 × 2^{e\overline{z}}.  (7)

If (7) is true, then the second case of Definition 4 gives

μ_{⊕}(z) = b_{1}.b_{2} ⋯ b_{i−1}a0 ⋯ × 2^{e\overline{z}},  with a = 0,

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

18  00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS

which is indeed equal to \underline{z}. On the other hand, if (7) is false, then no float with significand b_{1}.b_{2} ⋯ b_{i−1}00 ⋯ 0 belongs to [\underline{z}, \overline{z}], hence the significand maximizing α is necessarily the one with one less zero to the right, i.e., b_{1}.b_{2} ⋯ b_{i−1}10 ⋯ 0, which is guaranteed to belong to [\underline{z}, \overline{z}]. This is consistent with the second case of Definition 4, which gives

μ_{⊕}(z) = b_{1}.b_{2} ⋯ b_{i−1}a0 ⋯ 0 × 2^{e\overline{z}},  with a = 1.

We have proved that Definition 4 gives a float μ_{⊕}(z) that maximizes the value α. We now prove that μ_{⊕}(z) also maximizes the value of β. By Propositions 1 and 3 and Definition 3, β = α⊕z. Note that μ_{⊕}(z) maximizes α; however, since β also depends on z, we have to prove that no z ∈ [\underline{z}, \overline{z}] such that z > μ_{⊕}(z) results into a greater β. Observe first that, by construction, μ_{⊕}(z) has the maximum exponent in [\underline{z}, \overline{z}]. Therefore any z > μ_{⊕}(z) in [\underline{z}, \overline{z}] must have a larger significand. Assume that μ_{⊕}(z) = b_{1}.b_{2} ⋯ b_{j}0 ⋯ 0×2^{e\overline{z}} with b_{j} = 1 for some j ∈ \{1, …, p\}. The exponent of the corresponding α is e\overline{z} + p − j. Suppose now there exists z > μ_{⊕}(z) in [\underline{z}, \overline{z}] with a larger significand: this must have the form b_{1}.b_{2} ⋯ b_{ℓ}0 ⋯ 0 × 2^{e\overline{z}} with b_{ℓ} = 1 and j < ℓ ≤ p. The exponent of the corresponding α is e\overline{z} + p − ℓ, which is smaller than the α computed for μ_{⊕}(z) by at least one unit. Hence, we can conclude that b_{1}.b_{2} ⋯ b_{j}0 ⋯ 0 × 2^{e\overline{z}} + 1.1 ⋯ 1 × 2^{e\overline{z}+p−j} > b_{1}.b_{2} ⋯ b_{l}0 ⋯ 0 × 2^{e\overline{z}} + 1.1 ⋯ 1 × 2^{e\overline{z}+p−ℓ}, since ℓ > j. This shows that the float μ_{⊕}(z) also maximizes the value of β. We have proved that Definition 4 gives a float μ_{⊕}(z) that maximizes the value of both α and β over z. Since Definition 3 defines δ̄_{⊕}(z) = α for −∞ < z < 0 and δ̄_{⊕}(z) = β for 0 < z < +∞, we can conclude that, for each z ∈ [\underline{z}, \overline{z}], δ̄_{⊕}(z) ≤ δ̄_{⊕}(μ_{⊕}(z)). □

As we have already pointed out, the algorithm of Definition 4, if restricted to normal numbers, is similar to the algorithm presented in (Marre and Michel 2010). There is an important difference, though, in the case when \underline{z} = b_{1}.b_{2} ⋯ b_{i−1}b_{i}0 ⋯ 0 × 2^{e\overline{z}}, \overline{z} = b_{1}.b_{2} ⋯ b_{i−1}b′_{i} ⋯ × 2^{e\overline{z}} and \underline{z} > 0. In this case the algorithm of Marre and Michel (2010) erroneously returns b_{1}.b_{2} ⋯ b_{i−1}10 ⋯ 0×2^{e\overline{z}} instead of the value that maximizes α, i.e., \underline{z}, which is correctly computed by our algorithm.

For efficiency reasons, filtering by maximum ULP might be applied only when δ̄_{⊕}(μ_{⊕}(z)) ≤ f_{max} so as to avoid the use of wider floating-point formats.

In order to define δ̄_{⊖}, δ̄′_{⊖}, δ̄_{⊖} and δ̄′_{⊖}, we can use the following observation. Since x⊖y = [x−y]_{n} = [x + −y]_{n} = x ⊕ −y, the constraints z = x ⊖ y and z = x ⊕ −y are equivalent. Thus we have δ̄_{⊖} = δ̄_{⊕} and δ̄_{⊖} = δ̄_{⊕}, while δ̄′_{⊖} = − δ̄_{⊕} and δ̄′_{⊖} = − δ̄_{⊕} since, if −y ∈ [δ̄_{⊕}(z), δ̄_{⊕}(z)], then y ∈ [− δ̄_{⊕}(z), − δ̄_{⊕}(z)]. Moreover, since μ_{⊕}(z) maximizes δ̄_{⊕} and minimizes δ̄_{⊕} over an interval of floating-point values z, μ_{⊕}(z) can be used as well to maximize δ̄′_{⊖} and minimize δ̄′_{⊖} on z.

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS  19

4.6. Filtering by Maximum ULP on Multiplication For filtering multiplication constraints of the form z = x ⊗ y (and similarly for division), we cannot rely on the same maximum ULP property identified by Marre and Michel upon which the treatment of addition and subtraction rests. This is because the ULP property of z is only loosely related to the ULP property of x and y when they are being multiplied. Our generalized property, instead, covers also multiplication (and division, as we will see in Section 4.7). As indicated in (1) and (5), we have to the determine the maximum and minimum values for x satisfying z = x ⊗ y.

Consider a strictly positive constant z ∈ 𝔽_{p,e_{max}} and two unknowns x, y ∈ 𝔽_{p,e_{max}}^{sub} such that z = x ⊗ y. If z ≤ f_{max}/f_{min}, there exists a greatest float x_{m} ∈ 𝔽_{p,e_{max}}^{sub} such that there exists y ∈ 𝔽_{p,e_{max}}^{sub} satisfying z = x_{m} ⊗ y. More precisely, x_{m} must satisfy z = x_{m} ⊗ f_{min} and it turns out that we can take x_{m} = z ⊘ f_{min}. Since, for z ≤ f_{max}/f_{min}, division of z by f_{min} = 2^{e_{min}+1−p} amounts to an exponent shifting, we have that 𝔽_{p,e_{max}}^{sub} ∋ x_{m} = z/f_{min}. Moreover, we have that x_{m} = z/f_{min} is the greatest float such that z = x_{m} ⊗ f_{min}.^{9}

On the other hand, there is no other float y < f_{min} such that z = x ⊗ y, since y must be greater than +0, for otherwise x ⊗ y would not be strictly positive. However, for no y ∈ 𝔽_{p,e_{max}}^{sub} we have +0 < y < f_{min}. Therefore, the greatest value x_{m} such that z = x_{m} ⊗ f_{min} is the greatest value for x that can satisfy z = x ⊗ y for some y ∈ 𝔽_{p,e_{max}}^{sub}.

When dealing with subnormal floating-point numbers a similar argument applies. In fact, also in this case there exists a greatest float x_{m} ∈ 𝔽_{p,e_{max}}^{sub} satisfying z = x_{m} ⊗ y for some y ∈ 𝔽_{p,e_{max}}^{sub}. As before, such x_{m} must satisfy z = x_{m} ⊗ f_{min}. However, it turns out that, when z is subnormal, there may exist values for x_{m} greater than z/f_{min} that still satisfy z = x_{m} ⊗ f_{min}. This is because the distance between subnormal numbers, being fixed to f_{min}, does not depend on z.

Based on the previous reasoning, we can define δ̄_{⊗} and δ̄_{⊗}.

Definition 5. Let 𝔽_{⊗} = { z ∈ 𝔽_{p,e_{max}}^{sub} | |z|/f_{min} ≤ f_{max} } and 𝔽̄_{⊗} = 𝔽_{p,e_{max}}. Then δ̄_{⊗}: 𝔽_{⊗} → 𝔽̄_{⊗} is defined, for each z ∈ 𝔽_{⊗}, by

δ̄_{⊗}(z) = \begin{cases}|z| ⊘ f_{min}, & if |z| ≥ f_{min}^{nor}; \\ (|z| ⊘ f_{min}) ⊕ 2^{−1}, & if 0 < |z| < f_{min}^{nor} and even(z); \\ ((|z| ⊘ f_{min}) ⊕ 2^{−1})^{−}, & if 0 < |z| < f_{min}^{nor} and odd(z). \end{cases}

Theorem 3. Function δ̄_{⊗} is well-defined and satisfies (2) and (3).

Proof. Given in the Online Supplement (Bagnara et al. 2015).

A monotonicity property of δ̄_{⊗} simplifies the identification an element of the interval z that maximizes the value of δ̄_{⊗} over z.

9 See the proof of forthcoming Theorem 3 available in the Online Supplement (Bagnara et al. 2015).

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

20  00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS

Proposition 4. Let z ∈ 𝔽_{⊗} be nonzero. If z > 0, then δ̄_{⊗}(z^{+}) ≥ δ̄_{⊗}(z); on the other hand, if z < 0, then δ̄_{⊗}(z^{−}) ≥ δ̄_{⊗}(z).

Proof. Given in the Online Supplement (Bagnara et al. 2015).

Since ⊗ is commutative, δ̄′_{⊗} = δ̄_{⊗}, and the same bounds can be used to filter both x and y in the constraint z = x ⊗ y.

The function δ̄_{⊗}: 𝔽_{⊗} → 𝔽̄_{⊗} is defined dually: for each z ∈ 𝔽_{⊗} \ {−0, +0}, δ̄_{⊗}(z) = − δ̄_{⊗}(z). We can see that properties (2) and (3) of δ̄_{⊗} entail property (5) of δ̄_{⊗}. Again, since ⊗ is commutative we have δ̄′_{⊗} = δ̄_{⊗}.

Thanks to Proposition 4 we know that the value M ∈ [\underline{z}, \overline{z}] that maximizes δ̄_{⊗} is the one with the greatest absolute value, i.e., M = max\{|\underline{z}|, |\overline{z}|\}. Since δ̄_{⊗} is defined as − δ̄_{⊗}(z), the value that minimizes δ̄_{⊗} is again M. Hence, if [\underline{z}, \overline{z}] does not contain zeros, δ̄_{⊗}(M) (resp., δ̄_{⊗}(M)) is an upper bound (resp., a lower bound) of x with respect to the constraint z = x ⊗ y.

The restriction to intervals z not containing zeros is justified by the fact that, e.g., if z = 0 then z = x ⊗ y holds with x = f_{max} and y = 0, hence, in this case, no useful filtering can be applied to x. The same thing happens when max\{|\underline{z}|, |\overline{z}|\}/f_{min} > f_{max}. Moreover, whenever the interval of y does not contain zeros, filtering by maximum ULP for multiplication, in order to refine x, is subsumed by the standard indirect projection, which, in this case, can usefully exploit the information on y. In contrast, when the interval of y does contain zeros, our filter is able to derive bounds that cannot be obtained with the standard indirect projection, which, in this case, does not allow any refinement of the interval. Thus, for multiplication (and, as we will see, for division as well), the standard indirect projection and filtering by maximum ULP are mutually exclusive: one applies when the other cannot derive anything useful. Commenting on a previous version of the present paper, Claude Michel observed that one could modify the standard indirect projections with interval splitting so that indirect projections are always applied to source intervals not containing zeros. This idea rests on the observation that, for z = x ⊙ y with ⊙ ∈ {⊗, ⊘}, when the interval of z is a subset of the finite non zero floats neither x nor y do have any support for ±0 and ±∞. For multiplication, ordinary standard indirect projection would be modified as follows, assuming that z is positive and that we want to apply the standard indirect projection to z and y in order to refine x (the other cases being similar):

• we apply the ordinary standard indirect projection to z and y ∩ [−f_{max}, −f_{min}], intersecting the resulting interval with [−f_{max}, −f_{min}];

• we apply the ordinary standard indirect projection to z and y ∩ [f_{min}, f_{max}], intersecting the resulting interval with [f_{min}, f_{max}];

• finally, we use the convex union of the two intervals so computed to refine x.

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS  21

We believe that, when the applied ordinary (i.e., non-splitting) standard indirect projection is as precise as the one specified by Michel (2002), the refining interval computed for x by the modified procedure is very precise: it coincides with the result of the ordinary standard indirect projection, when 0 ∉ y and thus filtering by maximum ULP is not applicable, or it coincides with the result of filtering by maximum ULP, when 0 ∈ y and therefore the ordinary standard indirect projection would not help.^{10} This approach has the advantage to be applicable to any rounding mode. On the other hand the standard indirect projections specified in (Michel 2002) require working on rationals or on larger floating-point formats, whereas one of our aims is to always work with machine floating-point numbers of the same size of those used in the analyzed computation.

Example 2. Consider the IEEE 754 single-precision constraint z = x ⊗ y with z subnormal, z ∈ [−0.00000000000000010001001 × 2^{−126}, −0.00000000000010000000000 × 2^{−126}], and x and y unconstrained, x, y ∈ [−∞, +∞]. Our indirect projection infers the constraints x, y ∈ [−1.00000000001 × 2^{10}, 1.00000000001 × 2^{10}], while classical inverse projections do not allow pruning the intervals for x and y, no matter what they are.

4.7. Filtering by Maximum ULP on Division We now define filtering by maximum ULP for floating-point constraints of the form z = x ⊘ y. We begin defining the first indirect projection. We will then tackle the problem of defining the second indirect projection, which, as we will see, is significantly more involved than the first one: the solution we propose is new to this paper.

4.7.1. The First Indirect Projection A role similar to the one of f_{min} in the definition of filtering by maximum ULP on multiplication is played by f_{max} in the definition of the first indirect projection for division.

Definition 6. Let us define the sets 𝔽′_{⊘} = { z ∈ 𝔽_{p,e_{max}}^{sub} | |z| ⊗ f_{max} ≤ f_{max} } and 𝔽̄′_{⊘} = 𝔽_{p,e_{max}}. Let also q = 1 − p + e_{min} + e_{max}.^{11} Then δ̄_{⊘}: 𝔽′_{⊘} → 𝔽̄′_{⊘} is defined, for each z ∈ 𝔽′_{⊘}, by

δ̄_{⊘}(z) = \begin{cases}|z| ⊗ f_{max}, & if f_{min}^{nor} ≤ |z| ≤ 1; \\ (|z| ⊗ f_{max}) ⊕ 2^{q}, & if 0 ≤ |z| < f_{min}^{nor} \\ & ∧ (|z| ≠ 1 × 2^{e_{z}} ∨ e_{z} = e_{min} − 1); \\ ((|z| ⊗ f_{max}) ⊕ 2^{q})^{−}, & otherwise. \end{cases}

Observe that we have |z| ⊗ f_{max} ≤ f_{max} if and only if |z| ≤ 1. In fact, for z = 1^{+} = 1 + 2^{1−p}, we obtain

|z| ⊗ f_{max} = (1 + 2^{1−p}) ⊗ f_{max}

10 We are indebted to Claude Michel for this observation. 11 In the very common case where e_{min} = 1 − e_{max} we have q = 2 − p.

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

22  00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS
= [(1 + 2^{1−p})f_{max}]_{n}
= [f_{max} + (2 − 2^{1−p})2^{e_{max}+1−p}]_{n}
= +∞,  (8)

where (8) holds by Definition 2, since (2 − 2^{1−p})2^{e_{max}+1−p} > Δ_{f_{max}}^{+} /2 = 2^{e_{max}−p}. By monotonicity of ⊗ we can conclude that z ∈ 𝔽′_{⊘} if and only if |z| ≤ 1.

Theorem 4. δ̄_{⊘} is well-defined and satisfies (2) and (3).

Proof. Given in the Online Supplement (Bagnara et al. 2015).

The function δ̄_{⊘} is defined, for each z ∈ 𝔽′_{⊘}, by δ̄_{⊘} = − δ̄_{⊘}(z).

As for multiplication, a monotonicity property of δ̄_{⊘} enables quickly identifying the value of z that maximizes the function.

Proposition 5. Let z ∈ 𝔽_{⊘} be nonzero. If z > 0, then δ̄_{⊘}(z^{+}) ≥ δ̄_{⊘}(z); on the other hand, if z < 0, then δ̄_{⊘}(z^{−}) ≥ δ̄_{⊘}(z).

Proof. Given in the Online Supplement (Bagnara et al. 2015).

By monotonicity, the value M ∈ [\underline{z}, \overline{z}] that maximizes δ̄_{⊘} is the one that has the greatest absolute value, i.e., M = max\{|\underline{z}|, |\overline{z}|\}. Since δ̄_{⊘} is defined as − δ̄_{⊘}(z), M is also the value that minimizes δ̄_{⊘}. Hence, if [\underline{z}, \overline{z}] does not contain zeros, δ̄_{⊘}(M) (resp., δ̄_{⊘}(M)) is an upper bound (resp. a lower bound) of x with respect to the constraint z = x⊘y. The restriction to intervals not containing zeros is justified by the fact that, e.g., if z = 0 then z = x⊘y holds with x = f_{max} and y = ∞; hence, in this case, no useful filtering can be applied to x. The same happens when max\{|\underline{z}|, |\overline{z}|\} ⊗ f_{max} > f_{max}. In addition, the same phenomenon we saw for multiplication manifests itself here: whenever the interval of the variable y does not contain infinities, filtering by maximum ULP for division in order to refine x is subsumed by the standard indirect projection. On the other hand, when the interval of y does contain infinities, the standard indirect projection gives nothing whereas filtering by maximum ULP provides nontrivial bounds. Thus, the standard indirect projection and filtering by maximum ULP for division are mutually exclusive: one applies when the other cannot derive anything useful. And, just as for multiplication, if using rationals or extended floating-point formats is an option, then a pruning variant (one that cuts off infinities) of the indirect projection specified in (Michel 2002) will be equally precise.

Example 3. Consider the IEEE 754 single-precision constraint z = x ⊘ y with initial intervals z ∈ [−1.0 × 2^{−110}, −1.0 × 2^{−121}] and x, y ∈ [−∞, +∞]. We have

δ̄_{⊘}(1.0 × 2^{−110}) = 1.0 × 2^{−110} · 1.1 ⋯ 1 × 2^{127}
= 1.1 ⋯ 1 × 2^{17},

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS  23
δ̄_{⊘}(1.0 × 2^{−110}) = −1.0 × 2^{−110} · 1.1 ⋯ 1 × 2^{127}
= −1.1 ⋯ 1 × 2^{17}.

Filtering by maximum ULP improves upon classical filtering, which would not restrict any interval, with x ∈ [−1.1 … 1 × 2^{17}, 1.1 … 1 × 2^{17}].

For an example involving subnormals, consider z = x ⊘ y with initial interval for z equal to [0.00000000000000000000001×2^{−126}, 0.01×2^{−126}] and x, y ∈ [−∞, +∞]: our algorithm produces the constraint x ∈ [−1.00000000000000000000001 × 2^{−46}, 1.00000000000000000000001 × 2^{−46}] whereas classical filtering is unable to infer anything on x.

4.7.2. The Second Indirect Projection The discussion in Section 4.7.1 shows that, for |z| ≤ 1, we have δ̄′_{⊘}(z) = f_{max}. We thus need to study δ̄′_{⊘}(z) for |z| > 1. It turns out that, due to rounding, the restriction of δ̄′_{⊘} over that subdomain is not a simple function. Given z ∈ 𝔽_{p,e_{max}}^{sub}, δ̄′_{⊘}(z) is the maximum y such that x ⊘ y = z. Note that, in order to maximize y, x must be maximized as well. A qualitative reasoning on the reals tells us that, since f_{max}/(f_{max}/z) = z, y should be roughly equal to f_{max}/|z|. Indeed, it can be proved that, for |z| > 1, f_{max} ⊘ (f_{max} ⊘ |z|) is equal to z, z^{−} or z^{+} depending on the value of z. This allows the determination of a rather small upper bound to the values that z may take, which is ultimately our goal for filtering y values. To this aim we define the function δ̃′_{⊘}.

Definition 7. The function δ̃′_{⊘}: 𝔽_{p,e_{max}}^{sub} → 𝔽_{p,e_{max}}^{+} is defined, for each z ∈ 𝔽_{p,e_{max}}^{sub}, as follows:

δ̃′_{⊘}(z) = \begin{cases}f_{max} ⊘ |z|^{− −}, & if 1^{+} < |z| ≤ f_{max}; \\ f_{max}, & otherwise. \end{cases}

It turns out that δ̃′_{⊘}(z) satisfies the dual of Property (3), i.e., it is a correct upper bound, while it does not satisfy the dual of Property (2), i.e., smaller correct upper bounds might exist.

Theorem 5. Let 𝔽′′_{⊘} = 𝔽_{p,e_{max}}^{sub} and 𝔽̄′′_{⊘} = 𝔽_{p,e_{max}}^{+}. Let δ̄′_{⊘}: 𝔽′′_{⊘} → 𝔽̄′′_{⊘} be a function satisfying (4). Then, for 0 < |z| ≤ 1^{+} or z = +∞, δ̄′_{⊘}(z) ≤ δ̃′_{⊘}(z); moreover, for 1^{+} < |z| ≤ f_{max}, δ̄′_{⊘}(z) < δ̃′_{⊘}(z).

Proof. Given in the Online Supplement (Bagnara et al. 2015).

Dually, a lower bound for the function δ̄′_{⊘} can be obtained by means of the function \underline{δ̃}′_{⊘}, defined by \underline{δ̃}′_{⊘}(z) = − δ̃′_{⊘}(z).

The value N ∈ [\underline{z}, \overline{z}] that maximizes δ̃′_{⊘} is the one that has the smallest absolute value, i.e., N = min\{|\underline{z}|, |\overline{z}|\}. Since \underline{δ̃}′_{⊘} is defined as − δ̃′_{⊘}(z), N is also the value that minimizes \underline{δ̃}′_{⊘}. Thus, if [\underline{z}, \overline{z}] does not contain zeros, δ̃′_{⊘}(N) (resp., \underline{δ̃}′_{⊘}(N)) is an upper bound (resp. a lower bound) for x with respect to the constraint z = x ⊘ y. The restriction to intervals not containing zeros is justified by the fact that if, e.g., z = 0, then the equality z = x⊘y holds with y = ∞ for each x such

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

24  00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS

that 0 ≤ x ≤ f_{max}. Hence, as in the case of the first projection, no useful filtering can be applied to y. Analogously to the case of the filter for the first projection, this filter is useful whenever the interval of x contains infinities. In this case, in fact, it is able to derive useful bounds for y where the standard indirect projection does not allow any refinement of the interval. Just as is the case for multiplication and the first indirect projection of division, the standard indirect projection and filtering by maximum ULP are mutually exclusive: one applies when the other cannot derive anything useful.

Note that, only for this projection, we have chosen to compute a (very small) upper bound that, in general, is not the least upper bound. We did so in order to trade precision for efficiency: this way we have an algorithm that only uses floating-point machine arithmetic operations on the same format used by the analyzed constraint z = x ⊘ y. When using rationals or larger floating-point formats is an option, a pruning variant (as in the previous case, one that cuts off infinities) of a second indirect projection satisfying the precision constraints set forth in (Michel 2002) may result in extra precision at a higher computational cost.

Example 4. Consider the IEEE 754 single-precision division constraint z = x ⊘ y with initial intervals z ∈ [1.0 ⋯ 010 × 2^{110}, 1.0 × 2^{121}] and x, y ∈ [−∞, +∞]. We have

δ̃′_{⊘}(1.0 ⋯ 01 × 2^{110}) = 1.1 ⋯ 1 × 2^{127} ⊘ ((1.0 ⋯ 01 × 2^{110})^{−})^{−}
= 1.1 ⋯ 1 × 2^{127} ⊘ 1.1 ⋯ 1 × 2^{109}
= 1.0 × 2^{18},
\underline{δ̃}′_{⊘}(1.0 ⋯ 01 × 2^{110}) = −1.1 ⋯ 1 × 2^{127} ⊘ ((1.0 ⋯ 01 × 2^{110})^{−})^{−}
= −1.0 × 2^{18}.

Filtering by maximum ULP improves upon classical filtering, which gives nothing, with the constraint y ∈ [−1.0 × 2^{18}, 1.0 × 2^{18}].

4.8. Synthesis Table 1 provides a compact presentation of filtering by maximum ULP.

Table 1 Filtering by maximum ULP synopsis

Constraint  x ⊆ ·  y ⊆ ·  Condition(s)
z = x ⊕ y, 0 < z ≤ f_{max}  [ δ̄_{⊕}(ζ ), δ̄_{⊕}(ζ )] [ δ̄_{⊕}(ζ ), δ̄_{⊕}(ζ )] ζ = μ_{⊕}(z), −f_{max} ≤ δ̄_{⊕}(ζ), δ̄_{⊕}(ζ) ≤ f_{max}
z = x ⊕ y, −f_{max} ≤ z < 0  [− δ̄_{⊕}(ζ′), − δ̄_{⊕}(ζ′)] [− δ̄_{⊕}(ζ′), − δ̄_{⊕}(ζ′)] ζ′ = μ_{⊕}(−z), −f_{max} ≤ δ̄_{⊕}(ζ′), δ̄_{⊕}(ζ′) ≤ f_{max}
z = x ⊖ y, 0 < z ≤ f_{max}  [ δ̄_{⊕}(ζ ), δ̄_{⊕}(ζ )] [− δ̄_{⊕}(ζ ), − δ̄_{⊕}(ζ )] ζ = μ_{⊕}(z), −f_{max} ≤ δ̄_{⊕}(ζ), δ̄_{⊕}(ζ) ≤ f_{max}
z = x ⊖ y, −f_{max} ≤ z < 0  [− δ̄_{⊕}(ζ′), − δ̄_{⊕}(ζ′)] [ δ̄_{⊕}(ζ′), δ̄_{⊕}(ζ′)] ζ′ = μ_{⊕}(−z), −f_{max} ≤ δ̄_{⊕}(ζ′), δ̄_{⊕}(ζ′) ≤ f_{max}
z = x ⊗ y, |z| ≤ 2^{2−p}(2 − 2^{1−p}) [ δ̄_{⊗}(m), δ̄_{⊗}(m)] [ δ̄_{⊗}(m), δ̄_{⊗}(m)] m = max\{|\underline{z}|, |\overline{z}|\}
z = x ⊘ y, 0 < |z| ≤ 1  [ δ̄_{⊘}(m), δ̄_{⊘}(m)] [−f_{max}, + f_{max}] m = max\{|\underline{z}|, |\overline{z}|\}
z = x ⊘ y, 1 < |z| ≤ f_{max}  [ \underline{δ̃}′_{⊘}(n), δ̃′_{⊘}(n)] n = min\{|\underline{z}|, |\overline{z}|\}

δ̄_{⊕}(z) = \begin{cases}β, & if 0 < z < +∞, \\ α, & if −∞ < z < 0; \end{cases}  δ̄_{⊕}(z) = − δ̄_{⊕}(−z);

δ̄_{⊗}(z) = \begin{cases}|z| ⊘ f_{min}, & if z ≥ f_{min}^{nor}; \\ (|z| ⊘ f_{min}) ⊕ 2^{−1}, & if 0 < z < f_{min}^{nor} and even(z); \\ ((|z| ⊘ f_{min}) ⊕ 2^{−1})^{−}, & if 0 < z < f_{min}^{nor} and odd(z); \end{cases}  δ̄_{⊗}(z) = − δ̄_{⊗}(z);

δ̄_{⊘}(z) = \begin{cases}|z| ⊗ f_{max}, & if f_{min}^{nor} ≤ |z| ≤ 1; \\ (|z| ⊗ f_{max}) ⊕ 2^{q},^{(∗)} & if 0 ≤ |z| < f_{min}^{nor} ∧ (|z| ≠ 1 × 2^{e_{z}} ∨ e_{z} = e_{min} − 1); \\ ((|z| ⊗ f_{max}) ⊕ 2^{q})^{−}, & otherwise; \end{cases}  δ̄_{⊘}(z) = − δ̄_{⊘}(z);

δ̃′_{⊘}(z) = \begin{cases}f_{max} ⊘ |z|^{− −}, & if 1^{+} < |z| ≤ f_{max}; \\ f_{max}, & otherwise; \end{cases}  \underline{δ̃}′_{⊘}(z) = − δ̃′_{⊘}(z);

(∗) q = 1 − p + e_{min} + e_{max}.

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

26  00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS

5. Discussion

This work is part of a long-term research effort concerning the correct, precise and efficient handling of floating-point constraints (Belaid 2013, Belaid et al. 2012, Botella et al. 2006, Carlier and Gotlieb 2011, Marre and Michel 2010, Michel 2002, Michel et al. 2001) for software verification purposes.

Restricting the attention to test data generation other authors have considered using searchbased techniques with a specific notion of distance in their fitness function (Lakhotia et al. 2010a,b). For instance, search-based tools like AUSTIN and FloPSy can generate a test input for a specific path by evaluating the path covered by some current input with respect to a targeted path in the program. However, they cannot solve the constraints of path conditions, since: 1) they cannot determine unsatisfiability when the path is infeasible, and 2) they can fail to find a test input while the set of constraints is satisfiable (Bagnara et al. 2013).

Recently, Borges et al. (2012) combined a search-based test data generation engine with the RealPaver (Granvilliers and Benhamou 2006) interval constraint solver, which is well-known in the Constraint Programming community. Even though constraint solvers over continuous domains (e.g., RealPaver (Granvilliers and Benhamou 2006), Quimper (Chabert and Jaulin 2009) or ICOS (Lebbah 2009)) and the work described in the present paper are based on similar principles, the treatment of intervals is completely different. While our approach preserves all the solutions over the floats, it is not at all concerned with solutions over the reals. In contrast, RealPaver preserves solutions over the reals by making the appropriate choices in the rounding modes used for computing the interval bounds, but RealPaver can lose solutions over the floats. For instance, a constraint like (x > 0.0 ∧ x ⊕ 10000.0 ≤ 10000.0) is shown to be unsatisfiable on the reals by RealPaver, while it is satisfied by many IEEE 754 floating-point values of single or double precision format for x (Botella et al. 2006). Note that RealPaver has recently been used to tackle test input generation in presence of transcendental functions (Borges et al. 2012), but this approach, as mentioned by the authors of the cited paper, is neither correct nor complete due to the error rounding of floating-point computations.

The CBMC model checker (Clarke et al. 2004) supports floating-point arithmetic using a bit-precise floating-point decision procedure based on propositional encoding. According to (D’Silva et al. 2012) “CBMC translates the floating-point arithmetic to large propositional circuits which are hard for SAT solvers.” The technique presented in this paper is orthogonal to decision procedures over floating-point computations, such as those used in CBMC (Clarke et al. 2004) or CDFL (D’Silva et al. 2012). Of course, implementing the filtering procedures suggested here would require dedicated bitwise encodings, but this would enable to perform more constraint-based reasoning over these computations.

Bagnara et al.: Exploiting Binary Floating-Point Representations for Constraint Propagation

00(0), pp. 000–000, \begin{matrix}c \\ ◯ \end{matrix} 0000 INFORMS  27

6. Conclusion

This paper concerns constraint solving over binary floating-point numbers. Interval-based consistency techniques are very effective for the solution of such numerical constraints, provided precise and efficient filtering algorithms are available. We reformulated and corrected the filtering algorithm proposed by Marre and Michel (2010) for addition and subtraction. We proposed a uniform framework that generalizes the property identified by Marre and Michel to the case of multiplication and division. We also revised, corrected and extended our initial ideas, sketched in Carlier and Gotlieb (2011), to subnormals and to the effective treatment of floating-point division. All algorithms have been proved correct. In order to gain further confidence on the algorithms, we have exhaustively tested a first prototype, symbolic implementation on floating-point numbers with a small number of bits (e.g., p = 6 and e_{max} = 3). The implementation working on IEEE 754 formats was also tested with a variety of methodologies with the help of test-suites like the one by the IBM Labs in Haifa (2008).

An important objective of this work has been to allow maximum efficiency by defining all algorithms in terms of IEEE 754 elementary operations on the same formats as the ones of the filtered constraints. Indeed, the computational cost of filtering by maximum ULP as defined in the present paper and properly implemented is negligible. In fact, all the filters defined in this paper can be directly translated into constant-time algorithms (as IEEE 754 formats have fixed size) based on IEEE 754 elementary operations and simple bitwise manipulations. Moreover, for multiplication and division, precise conditions are given in order to decide whether standard filtering or our filtering by maximum ULP is applied: it is one or the other, never both of them. As shown in (Bagnara et al. 2013), the improvement of filtering procedures with these techniques brings significant speedups of the overall constraint solving process, with only occasional, negligible slowdowns. Note that the choice of different heuristics concerning the selection of constraints and variables to subject to filtering and the labeling strategy has a much more dramatic effect on solution time, even though the positive or negative effects of such heuristics change wildly from one analyzed program to the other. Filtering by maximum ULP contributes to reducing this variability. To understand this, consider the elementary constraint z = x ⊙ y: if x and y are subject to labeling before z, then filtering by maximum ULP will not help. However, z might be labeled before x or y: this can happen under any labeling heuristic and constitutes a performance bottleneck. In the latter case, filtering by maximum ULP may contribute to a much improved pruning of the domains of x and y and remove the bottleneck.

Future work includes coupling filtering by maximum ULP with sophisticated implementations of classical filtering based on multi-intervals and with dynamic linear relaxation algorithms (Denmat et al. 2007) using linear relaxation formulas such as the ones proposed by Belaid et al.

(2012). Another extension, by far more ambitious, concerns the correct handling of transcendental functions (i.e., sin, cos, exp, …): as IEEE 754 only provides recommendations rather than formal requirements for these functions, solutions will be dependent on the particular implementation and/or be imprecise; in other words, generated test inputs will not be applicable to other implementations and/or may fail to exercise the program paths they were supposed to traverse.

References

Ammann, P. E., J. C. Knight. 1988. Data diversity: An approach to software fault tolerance. IEEE Trans-

actions on Computers 37 418–425. Arcuri, A. 2009. Theoretical analysis of local search in software testing. O. Watanabe, T. Zeugmann, eds.,

Proceedings of the 5th International Symposium on Stochastic Algorithms: Foundations and Applications (SAGA 2009), Lecture Notes in Computer Science, vol. 5792. Springer-Verlag, Berlin, Sapporo, Japan, 156–168.

Bagnara, R., M. Carlier, R. Gori, A. Gotlieb. 2013. Symbolic path-oriented test data generation for floatingpoint programs. Proceedings of the 6th IEEE International Conference on Software Testing, Verification and Validation. IEEE Press, Luxembourg City, Luxembourg.

Bagnara, R., M. Carlier, R. Gori, A. Gotlieb. 2015. Online Supplement to “Exploiting Binary Floating-Point Representations for Constraint Filtering”. Available at

https://www.informs.org/Pubs/IJOC/Online-Supplements.

Belaid, M. S. 2013. Résolution de contraintes sur les flottants dédié à la vérification de programmes. Thèse pour obtenir le titre de “Docteur en Sciences”, École doctorale STIC, Université de Nice — Sophia Antipolis, Nice, France.

Belaid, M. S., C. Michel, M. Rueher. 2012. Boosting local consistency algorithms over floating-point numbers. M. Milano, ed., Proceedings of the 18th International Conference on Principles and Practice of Constraint Programming, Lecture Notes in Computer Science, vol. 7514. Springer-Verlag, Berlin, Québec City, Canada, 127–140.

Blanc, B., F. Bouquet, A. Gotlieb, B. Jeannet, T. Jeron, B. Legeard, B. Marre, C. Michel, M. Rueher. 2006. The V3F project. Proceedings of the 1st Workshop on Constraints in Software Testing, Verification and Analysis (CSTVA’06). Nantes, France.

Borges, M., M. d’Amorim, S. Anand, D. Bushnell, C. S. Pasareanu. 2012. Symbolic execution with interval solving and meta-heuristic search. Proceedings of the 5th IEEE International Conference on Software Testing, Verification and Validation. IEEE Computer Society, Montreal, Canada, 111–120.

Botella, B., A. Gotlieb, C. Michel. 2006. Symbolic execution of floating-point computations. Software Testing,

Verification and Reliability 16 97–121. Burdy, L., J.-L. Dufour, T. Lecomte. 2012. The B method takes up floating-point numbers. Proceedings of

the 6th International Conference & Exhibition on Embedded Real Time Software and Systems (ERTS

2012). Toulouse, France. Available at http://www.erts2012.org/Site/0P2RUC89/5C-2.pdf.

Carlier, M., A. Gotlieb. 2011. Filtering by ULP maximum. Proceedings of the 23rd IEEE International Conference on Tools with Artificial Intelligence (ICTAI 2011). IEEE Computer Society, Boca Raton, Florida, USA, 209–214.

Chabert, G., L. Jaulin. 2009. Contractor programming. Artificial Intelligence 173 1079–1100. Chan, F. T., T. Y. Chen, S. C. Cheung, M. F. Lau, S. M. Yiu. 1998. Application of metamorphic testing

in numerical analysis. Proceedings of the IASTED International Conference on Software Engineering (SE’98). ACTA Press, Las Vegas, Nevada, USA, 191–197.

Clarke, E. M., D. Kroening, F. Lerda. 2004. A tool for checking ANSI-C programs. K. Jensen, A. Podelski, eds., Tools and Algorithms for the Construction and Analysis of Systems, Proceedings of the 10th International Conference (TACAS 2004), Lecture Notes in Computer Science, vol. 2988. Springer, Barcelona, Spain, 168–176.

Cousot, P., R. Cousot. 1977. Abstract interpretation: A unified lattice model for static analysis of programs by construction or approximation of fixpoints. Proceedings of the Fourth Annual ACM Symposium on Principles of Programming Languages. ACM Press, Los Angeles, CA, USA, 238–252.

Denmat, T., A. Gotlieb, M. Ducassé. 2007. Improving constraint-based testing with dynamic linear relaxations. Proceedings of the 18th IEEE International Symposium on Software Reliability (ISSRE 2007). IEEE Computer Society, Trollhättan, Sweden, 181–190.

D’Silva, V., L. Haller, D. Kroening, M. Tautschnig. 2012. Numeric bounds analysis with conflict-driven learning. C. Flanagan, B. König, eds., Tools and Algorithms for the Construction and Analysis of Systems, Proceedings of the 18th International Conference (TACAS 2012), Lecture Notes in Computer Science, vol. 7214. Springer, Tallinn, Estonia, 48–63.

Godefroid, P., N. Klarlund, K. Sen. 2005. DART: Directed automated random testing. V. Sarkar, M. W. Hall, eds., Proceedings of the ACM SIGPLAN 2005 Conference on Programming Language Design and Implementation (PLDI 2005). ACM, Chicago, IL, USA, 213–223.

Goldberg, D. 1991. What every computer scientist should know about floating-point arithmetic. ACM

Computing Surveys 23 5–48. Goubault, E. 2001. Static analyses of the precision of floating-point operations. P. Cousot, ed., Static

Analysis: 8th International Symposium, SAS 2001 , Lecture Notes in Computer Science, vol. 2126. Springer-Verlag, Berlin, Paris, France, 234–259.

Granvilliers, L., F. Benhamou. 2006. Algorithm 852: RealPaver: An interval solver using constraint satisfac-

tion techniques. ACM Transactions on Mathematical Software 32 138–156.

IBM Labs in Haifa, FPgen Team. 2008.  Floating-point test-suite for IEEE.  Available at

https://www.research.ibm.com/haifa/projects/verification/fpgen/papers/ieee-test-suite-v2.pdf.

Version 1.02.

IEEE Computer Society. 2008. IEEE Standard for Floating-Point Arithmetic. The Institute of Electrical

and Electronics Engineers, Inc., IEEE Std 754-2008 (revision of IEEE Std 754-1985) ed. Korel, B. 1990. Automated software test data generation. IEEE Transactions on Software Engineering 16

870–879. Kuliamin, V. V. 2010. Standardization and testing of mathematical functions. A. Pnueli, I. Virbitskaite,

A. Voronkov, eds., Perspectives of Systems Informatics, Revised Papers from the 7th International Andrei Ershov Memorial Conference (PSI 2009), Novosibirsk, Russia, June 15–19, 2009 , Lecture Notes in Computer Science, vol. 5947. Springer-Verlag, Berlin, 257–268.

Lakhotia, K., M. Harman, H. Gross. 2010a. AUSTIN: A tool for search based software testing for the C language and its evaluation on deployed automotive systems. Proceedings of the 2nd International Symposium on Search Based Software Engineering (SSBSE ’10). IEEE Computer Society, 101–110.

Lakhotia, K., N. Tillmann, M. Harman, J. De Halleux. 2010b. FloPSy: Search-based floating point constraint solving for symbolic execution. Proceedings of the 22nd IFIP WG 6.1 International Conference on Testing Software and Systems. Springer-Verlag, Berlin, Heidelberg, Natal, Brazil, 142–157.

Lebbah, Y. 2009. ICOS: a branch and bound based solver for rigorous global optimization. Optimization

Methods and Software 24 709–726. Marre, B., B. Blanc. 2005. Test selection strategies for Lustre descriptions in GATeL. Y. Gurevich, A. K.

Petrenko, A. Kossatchev, eds., Proceedings of the Workshop on Model Based Testing (MBT 2004), Electronic Notes in Theoretical Computer Science, vol. 111. Elsevier Science Publishers B. V., Barcelona, Spain, 93–111.

Marre, B., C. Michel. 2010. Improving the floating point addition and subtraction constraints. D. Cohen, ed., Proceedings of the 16th International Conference on Principles and Practice of Constraint Programming (CP 2010), Lecture Notes in Computer Science, vol. 6308. Springer, St. Andrews, Scotland, UK, 360– 367.

McMinn, P. 2004. Search-based software test data generation: A survey. Software Testing, Verification and

Reliability 14 105–156. Michel, C. 2002. Exact projection functions for floating point number constraints. Proceedings of the 7th

International Symposium on Artificial Intelligence and Mathematics. Fort Lauderdale, FL, USA. Michel, C., M. Rueher, Y. Lebbah. 2001. Solving constraints over floating-point numbers. T. Walsh, ed.,

Proceedings of the 7th International Conference on Principles and Practice of Constraint Programming (CP 2001), Lecture Notes in Computer Science, vol. 2239. Springer-Verlag, Berlin, Paphos, Cyprus, 524–538.

Miller, W., D. L. Spooner. 1976. Automatic generation of floating-point test data. IEEE Transactions on Software Engineering 2 223–226.

Monniaux, D. 2008. The pitfalls of verifying floating-point computations. ACM Transactions on Program-

ming Languages and Systems 30. Muller, J.-M. 2005. On the definition of ulp(x). Rapport de recherche 5504, INRIA. Scott, N. S., F. Jézéquel, C. Denis, J.-M. Chesneaux. 2007. Numerical ‘health check’ for scientific codes: The

CADNA approach. Computer Physics Communications 176 507–521. Skeel, R. 1992. Roundoff error and the Patriot missile. SIAM News 25 11. Tang, E., E. T. Barr, X. Li, Z. Su. 2010. Perturbing numerical calculations for statistical analysis of floating-

point program (in)stability. P. Tonella, A. Orso, eds., Proceedings of the 19th International Symposium on Software Testing and Analysis (ISSTA 2010). ACM, Trento, Italy, 131–142.

Weyuker, E. J. 1982. On testing non-testable programs. The Computer Journal 25 465–470.

Acknowledgments

We are grateful to Abramo Bagnara (BUGSENG srl, Italy) for the many fruitful discussions we had on the subject of this paper, and to Paul Zimmermann (INRIA Lorraine, France) for the help he gave us proving a crucial result. We are also indebted to Claude Michel for several constructive remarks that allowed us to improve the paper. Finally, we wish to express our gratitude to the anonymous reviewers for the many useful suggestions they contributed.