FPScan: An Automated Constraint-Based Analyzer for Floating-Point Anomaly Detection Julien Bortolussi, Dorra Ben Khalifa, and Pierre-Loı̈c Garoche
arXiv:2609.07492v1 [cs.SE] 7 Sep 2026
Fédération ENAC ISAE-SUPAERO ONERA, Université de Toulouse, Toulouse, France {julien.bortolussi,dorra.ben-khalifa,pierre-loic.garoche}@enac.fr
Abstract—Writing error-free floating-point programs is a challenging task, especially for programmers who lack a strong background in numerical analysis and rounding-error propagation. State-of-the-art techniques typically aim to bound such errors using static or dynamic analysis. However, only a few tools explicitly address critical floating-point pitfalls such as absorption and catastrophic cancellation. These anomalies represent situations in which rounding errors are significantly amplified, causing the semantics of the finite-precision computation to deviate substantially from the real-number semantics. In this article, we present FPScan, a novel tool to formally define and detect both catastrophic cancellation and absorption in floating-point programs. Our approach starts with a custom static analyzer based on abstract interpretation to infer the order of magnitude of all program variables. This magnitude information is then used to build a set of first-order constraints that model error propagation and numerical precision within the program. Finally, we employ an off-the-shelf SMT solver to determine whether the program exhibits any of these critical numerical pitfalls. Experiments were conducted on FPBench, a well-known benchmark suite of floating-point programs, to evaluate the effectiveness of our tool. We also present a comparison with state-of-the-art tools regarding soundness and analysis time. Index Terms—Floating-point arithmetic, Static analysis, Catastrophic cancellation, Absorption, Constraint generation, SMT solver.
I. I NTRODUCTION Floating-point arithmetic is the standard approach for approximating real-number computations in modern computing systems [1]. By representing numbers in a scientific-notation format with a finite significand and a bounded exponent, it supports a large dynamic range while keeping relative errors small. However, this representation implies that floating-point numbers are not uniformly distributed over the real line. As a consequence, certain operations may introduce significant numerical errors. Two well-known examples are absorption and catastrophic cancellation. Absorption occurs when adding numbers with vastly different magnitudes, causing the smaller value to be lost due to rounding. Catastrophic cancellation arises when subtracting nearly equal numbers, amplifying previously accumulated rounding errors. In practice, many floating-point expressions are implemented by developers without specialized training in numerical analysis, which makes diagnosing such numerical issues difficult. As a result, developers often rely on iterative trialand-error modifications guided by test inputs, an approach that is both inefficient and unreliable [2].
Floating-point experts typically describe these phenomena as numerical pitfalls that explain the origins of many floatingpoint errors [2]–[4]. However, relatively little work has investigated how these pitfalls could be systematically exploited, either to assist developers in identifying problematic computations or to automatically generate more accurate expressions. In this article, we introduce FPScan1 , a tool for soundly and automatically detecting floating-point pitfalls. FPScan is based on the principle that pitfalls are closely related to the relative orders of magnitude of operands and the magnitude of the associated rounding errors. To exploit this property, FPScan first performs a static range analysis to infer interval bounds for program variables and intermediate expressions. These bounds are then used to approximate the orders of magnitude of values involved in computations. FPScan abstracts these orders of magnitude using integer variables and generates constraints that characterize conditions under which numerical pitfalls, such as absorption and catastrophic cancellation, may occur. The resulting constraints are checked using the Z3 SMT solver to determine whether such situations are detected within the inferred ranges. We evaluate FPScan using the FPBench benchmark suite [5] and compare it with FPChecker [6], [7], a state-of-the-art dynamic tool for detecting floatingpoint cancellations, and with a bitblasting-based approach. Our results show that FPScan achieves high precision while remaining computationally efficient. The three main contributions of this article are: 1) An efficient constraint-based method for detecting absorption and catastrophic cancellation in floating-point computations (including denormalized numbers) based on order-of-magnitude reasoning (sections II and III) . 2) FPScan, a tool based on static analysis that combines interval range analysis with SMT solving to automatically detect floating-point pitfalls (Section III) . 3) An experimental evaluation on the FPBench benchmark suite, including comparisons with FPChecker and a bitblasting-based approach (Section IV). It demonstrates the efficiency and accuracy of FPScan. II. BACKGROUND ON F LOATING -P OINT A RITHMETIC In this section, we review the fundamentals of floating-point arithmetic, with a focus on its key concepts and main sources of numerical errors. 1 Software artifact available at https://github.com/JBortolussi/FPScan
A. Floating-Point Numbers Floating-point arithmetic approximates real-number computations using a finite number of digits, making it suitable for computer systems. According to the IEEE 754 Standard [1], a floating-point number is represented as shown x = (−1)s × b0 .b1 b2 ...bp−1 × β e ,
1: 2: 3: 4:
Require x ∈ [1, 226 ] ∧ y ∈ [1, 10] z ←x+y result ← z − x Return result
Fig. 1. Running example program.
(1)
where s denotes the sign bit, b0 .b1 b2 ...bp−1 is the significand written in base β (here we consider β = 2), e is the exponent, and p is the precision. The base β, the precision p, and the bounds of the exponent emin and emax together define a floating-point format. The IEEE 754 Standard defines several floating-point formats that differ in precision and exponent range. The most common formats are binary32, binary64, and binary128, which provide 24, 53, and 113 bits of precision p, respectively. The floating-point representation of real numbers is not unique. For instance, with β = 10, the numbers 9.0 and 0.9 × 101 represent the same value. To address this ambiguity, floating-point numbers are normalized so that their significand is maximized, i.e., in this example, 9.0. According to Eq. (1), the leading digit of the significand (b0 ) is equal to one for all normalized numbers. Numbers that do not satisfy this property are called denormalized. The process of selecting which of the two nearest floatingpoint numbers represents a real value is called rounding. The IEEE 754 Standard [1] defines several rounding modes. In this work, we focus exclusively on the most common setting: the ↑p∼ rounding mode, which rounds a number to its nearest floating-point neighbour. B. Absorption and Catastrophic Cancellation in Floatingpoint Arithmetic Rounding errors arise in all floating-point operations and may accumulate, leading to significant loss of accuracy. Certain situations, known as floating-point pitfalls, can further amplify this effect. The specific focus of this work is to identify such pitfalls: absorption and catastrophic cancellation. While these phenomena are widely discussed in the literature, they are rarely given formal definitions. Let us first introduce them informally before giving a more precise definition in Section III-D. 1) Absorption: occurs when a large number is added to a small number, causing the contribution of the smaller number to be similar to a rounding error. Although only a regular rounding error is introduced in this situation, it can lead to substantial problems, as the result differs from what the programmer intended. Consider the program in Figure 1 with precision p = 24 (binary32), x = 226 , and y = 1. Line 1 computes z = ↑p∼ (226 + 1) = 226 . Although z is affected only by a rounding error, the information provided by y is lost. 2) Catastrophic Cancellation: occurs when subtracting two nearby numbers whose most significant bits coincide and cancel each other out making the result significantly smaller than the operands. Consequently, any error affecting one of the operands has a proportionally larger impact on the result.
Fig. 2. FPScan’s workflow.
Again, consider the program in Figure 1, with precision p = 24 (binary32), x = 226 , y = 1, and z = 226 . Line 2 computes result = ↑p∼ (226 −226 ) = 0 instead of 1. In this case, the subtraction between z and x leads to catastrophic cancellation, which reveals the pitfall caused by the absorption in Line 1. This large relative error in the result can then propagate significantly through subsequent computations. Unlike the definitions given in [7], [8], in this work we rely on the formalization of Unit in the First Place (ufp), Unit in the Last Place (ulp), and Number of Significant Bits (nsb). In the next section, we formally define absorption and catastrophic cancellation in terms of these quantities at the abstract-semantics level, expressed as the satisfiability of a constraint. III. C ONSTRAINT-BASED D ETECTION OF F LOATING -P OINT A BSORPTION AND C ATASTROPHIC C ANCELLATION In this section, we present the approach implemented in FPScan. The overall workflow of FPScan is illustrated in Figure 2. It first parses the floating-point program and performs a range analysis to infer bounds on variables and intermediate expressions. These bounds are then used to generate constraints characterizing the conditions under which absorption and catastrophic cancellation may occur. The resulting constraints are checked independently using an SMT solver to determine whether such situations are feasible. Floating-point pitfalls depend on the orders of magnitude of values and their errors. Therefore, both can be abstracted by their orders of magnitude, yielding integer constraints that are faster to solve. Notation: In this section, we use the following notations. Let p ∈ N and emin ∈ Z denote the precision and the minimum exponent of the floating-point format. Let Fp denote the set of floating-point numbers in this format. Let V be a set of variables and let Σ = (V → R × Fp ) be the set of environments, i.e., the set of machine states mapping each variable to its real and floating-point values. To simplify the notation, let σR and σFp be the projections of σ ∈ Σ onto
R and Fp , respectively. For all variables v ∈ V, we define ṽ = σFp (v). A. Quantifying Magnitude, Precision, and Error
Variables:
(xv , yv , zv , tv ) ∈ V4
Operators:
♢ ∈ {+, −, ×, ÷}
Functions:
We now introduce three quantities that characterize the order of magnitude of a floating-point number, its precision, and its associated error.
Comparisons: Statements:
f ∈ {cos, sin, atan, sqrt} ▷ ∈ {>, ≥, =, ̸=} s
:= | | | | | | |
Definition 1 (Unit in the First Place (ufp)). ufp : R x
−→ Z 7−→
min{i ∈ Z | 2i+1 > |x|} emin − p + 1
if x ̸= 0 else.
(2)
The ufp corresponds to the order of magnitude of the number: ∀x ∈ R, ∀n ∈ Z, ufp(x) ≤ n ⇐⇒ |x| < 2n+1 .
(3)
Definition 2 (Unit in the Last Place (ulp)). ulp : Fp x̃
−→ Z 7−→ max(ufp(x̃) − p + 1, emin − p + 1) .
(4)
The ulp corresponds to the exponent of the least significant bit and is valid for all floating-point formats. This definition requires special care when handling denormalized numbers. Let us note that the ulp is closely related to rounding errors [3], as shown in the following equation: p
∀x ∈ R, |x − ↑p∼ (x)| ≤ 2ulp(↑∼ (x))−1 .
s; s xv = n, n ∈ R xv = yv xv = yv ♢zv xv = f (yv ) if tv ▷ 0 then sthen else selse endif while tv ▷ 0 do s done noop
Fig. 3. FPScan language.
SJs1 ; s2 K(σ) := SJs2 K SJs1 K(σ) SJxv = nK(σ) := σ xv ← (n, ↑p∼ (n)) σR (xv )♢σR (yv ), SJzv = xv ♢yv K(σ) := σ zv ← p ↑∼ (σFp (xv )♢σFp (yv )) SJif tv ▷ 0 then sthen else selse endif K(σ) := if σR (tv ) ▷ 0 then SJsthen K else SJselse K SJxv = f (yv )K(σ) := σ xv ← (f (σR (yv )), ↑p∼ (f(σFp (yv )))) SJxv = yv K(σ) := σ [xv ← σ(yv )]
SJnoopK(σ) := σ
(5) Fig. 4. Concrete semantics of the language.
Definition 3 (Number of Significant Bits (nsb)). Let xv ∈ V a variable and σ ∈ Σ an environment. nsbσ : V −→ Z xv 7−→ ufp σFp (xv ) − ufp σR (xv ) − σFp (xv ) .
(6)
The nsb corresponds to the number of bits whose weight exceeds the rounding error. Consequently, it can take negative values when rounding error exceeds the magnitude of the value itself. B. Concrete Program Semantics The syntax of the language accepted by the FPScan parser is given in Figure 3. It supports basic arithmetic operators, mathematical functions, conditionals, and loops. In practice, FPScan takes a subset of FPCore [5] as input and translates it into this internal syntax. As highlighted by Titolo et al. [9], most static analysis tools support only statically bounded loops. FPScan is no exception and therefore only supports bounded loops, with the maximum number of iterations specified by the user. During a pre-processing stage, these loops are unrolled and replaced with a sequence of nested conditional statements that execute the loop body as long as the loop condition holds. The user-defined bound determines the depth of this conditional chain. Additionally, we assume that both real and floating-point executions follow the same control-flow path in conditional statements, i.e., test instabilities are not modeled. Handling unbounded loops and unstable tests is left for future work. The semantic function SJsK : Σ → Σ defines the semantics of a statement s by mapping a machine state to another
machine state according to the execution of s. Its definition is given in Figure 4. Here, σ ∈ Σ denotes a machine state; s1 and s2 denote statements; xv , yv , zv , and tv denote variables; n ∈ R denotes a real-valued constant; ▷ ∈ {>, ≥, =, ̸=} denotes a comparison operator; and f ∈ {cos, sin, atan, sqrt} denotes a supported function. In Eq. (7), we define CJsK : P(Σ) → P(Σ), the collecting semantic function that returns the set of reachable machine states resulting from executing a statement over a set of input states. Here, P(Σ) denotes the powerset of Σ. ∀X ∈ P(Σ), CJsK(X) = {SJsK(σ) | σ ∈ X}.
(7)
C. Approximated Constraint-based Program Semantics FPScan abstracts the semantic function C, defined in Eq. (7), using constraints over ufp, ulp, and nsb. Let L denote the set of the constraints. A constraint l : Σ → {true, false} is a boolean-valued function over a machine state. To formalize this abstraction, we introduce a concretization-based abstract interpretation [10]. The concretization function γ, defined in Eq. (8), maps each constraint to the set of machine states that satisfy it. γ: L l
−→ P(Σ) 7−→ {σ ∈ Σ|l(σ) = true}
(8)
The function C # JsK : L → L represents the abstract program semantics. To be considered correct, it must satisfy Eq. (9). In other words, it must not omit any reachable machine states. ∀l ∈ L, CJsK(γ(l)) ⊆ γ(C # JsK(l))
(9)
1) Range analysis: FPScan first performs an interval analysis to detect and reject programs that may produce overflows, NaNs, or divisions by zero. The inferred intervals, of the form [x, x], provide bounds on program variables and are then used to derive bounds on their ufp values using Eq. (3). To this end, we introduce the following three bounding constraints: UpperBoundJxv K = λσ. ufp(σFp (xv )) ≤ ufp max(|x|, |x|) , LowerBoundJxv K = λσ. ufp min(|x|, |x|) ≤ ufp(σFp (xv )), UpperBoundJxv K(σ) BoundJxv K = λσ. ∧ LowerBoundJxv K(σ) if 0 ∈ / [x, x].
2) Constraint generation: The constraints used to abstract program executions are derived from the concrete semantics of the language. For the remainder of this section, let xv , yv , zv and tv be program variables in V. Let s1 , s2 , sthen , and selse be statements in the language. We now define the constraints associated with each statement of the language. Since our constraints are expressed solely in terms of ufp, ulp, and nsb, they do not capture the sign of the operands. As a consequence, addition and subtraction are treated uniformly at the constraint level. Definition 4 (Addition/Subtraction). For all l ∈ L, the abstraction of the statement zv = xv ± yv is defined as C # Jzv = xv ± yv K(l) := λσ.
with
ufpe = max
max
+ 2,
.
Definition 5 (Multiplication). ∀ l ∈ L, the abstraction of the statement zv = xv × yv is defined as l(σ) ∧ BoundJzv K(σ) # C Jzv = xv × yv K(l) := λσ. ∧ nsbσ (zv ) ≥ ufp(z̃) − ufpe
ufp(x̃) + ufp(ỹ) − min nsbσ (xv ), nsbσ (yv ) + 3, ufpe = max ufp(x̃) + ufp(ỹ) − nsbσ (xv ) − nsbσ (yv ) + 3, . ulp(z̃)
Definition 6 (Division). ∀l ∈ L, the abstraction of the statement zv = xv ÷ yv is defined as
l(σ) ∧ BoundJzv K(σ) nsbσ (zv ) ≥ ufp(z̃) − ufpe C Jzv = xv ÷ yv K(l) := λσ. ∧ ∨ nsbσ (yv ) ≤ 1 #
with ufpe = max
ufp(x̃) − ufp(ỹ) − min(nsbσ (xv ), nsbσ (yv )) + 4 ulp(z̃) − 1
v
e
with ufp(x̃) − nsbσ (xv ) + 1 , ulp(z̃) − 1) . 2 Definition 9 (Assign Constant). Let n ∈ R be a constant, ∀l ∈ L, the abstraction of the statement xv := n is defined as " ufpe = max(
C # Jxv := nK(l) = λσ.
l(σ) ∧ BoundJxv K(σ) ∧ nsb(xv ) ≥ p
.
Definition 10 (Assign Variable). ∀l ∈ L, the abstraction of the statement xv := yv is defined as C # Jxv := yv K(l) = λσ.l(σ) ∧ xv = yv .
Definition 11 (Sequence). Let s1 and s2 denote two statements. The abstract semantics of the sequence is defined as ∀l ∈ L, C # Js1 ; s2 K(l) := C # Js2 K(C # Js1 K(l)) .
Definition 12 (If Then Else). Let ▷ be a comparison operator, then the abstract branching condition is defined as
Theorem 1 (Soundness of abstract semantics C # ). The abstract semantics C # , is a sound over-approximation of the concrete semantics C, i.e. satisfies Eq. (9).
!
ulp(z̃)
with
σ
∀l ∈ L, C # Jif tv ▷ 0 then sthen else selse endif K(l) := C # Jsthen K(l)(σ) ∨ C # Jselse K(l)(σ)
l(σ) ∧ BoundJzv K(σ) ∧ nsbσ (zv ) ≥ ufp(z̃) − ufpe
ufp(x̃) − nsbσ (xv ), ufp(ỹ) − nsbσ (yv )
Definition 8 (Sqrt). ∀l ∈ L, the abstraction of the statement √ zv = xv is defined as √ l(σ) ∧ BoundJzv K(σ) C # Jzv = xv K(l) := λσ. ∧ nsb (z ) ≥ ufp(z̃) − ufp
.
Definition 7 (Cos, Sin, and Atan). ∀l ∈ L, the abstraction of the statement zv = f (xv ) with f ∈ {cos, sin, atan} is defined as l(σ) ∧ BoundJzv K(σ) # C Jzv = f (xv )K(l) := λσ. ∧ nsbσ (zv ) ≥ ufp(z̃) − ufpe with ufp(x̃) − nsbσ (xv ) + 1 ufpe = max ulp(z̃) − 1 .
Proof. Soundness proof is performed by structural induction on the program constructs. Abstract math. functions soundness: Defs 4–8 constraints are all made of three parts. The first propagates previous constraints. The second bounds the ufp of variables using BoundJ.K. It is sound because the interval analysis producing the bound is. The third bounds the nsb. Let us develop the derivation of the bound on the nsb for the (Addition/Subtraction) constraint. The soundness proof of Defs 5–8 are provided in Appendix. Abstract addition soundness: Let εx , εy , and εz denote the floating-point errors: for v ∈ {x, y, z}, εv = v − ṽ ⇔ ṽ = v − εv ⇔ v = ṽ + εv . Let e+ be the rounding error introduced by the addition. We can derive z̃ = ↑p∼ (x̃ + ỹ) = x̃ + ỹ + e+ = (x − εx ) + (y − εy ) + e+ = x + y − (εx + εy − e+ ), characterizing εz as εx + εy − e+ . The error |εz | is bounded using the triangle inequality |εz | ≤ |εx | + |εy | + |e+ | while |εx | and |εy | are bounded by Eq. (3) and e+ by |e+ | ≤ 2ulp(z̃)−1 using Eq. (5). We then have |εz | < 2ufp(εx )+1 + 2ufp(εy )+1 + 2ulp(z̃)−1 < 2max(ufp(εx )+1,ufp(εy )+1,ulp(z̃)−1)+2 . Using Eq. (3) we obtain ufp(εz ) ≤ max(ufp(εx ) + 1, ufp(εy ) + 1, ulp(z̃) − 1) + 1. Abstract control flow soundness: Let us now focus on the statements of the language described in Defs 9–12. Def. 9 is sound because to round a real number to make it fit the floating-point format adds no more than a rounding error, hence nsb ≥ p. Def. 10 states that no additional error is introduced copying a variable and Def. 11 translates constraint
propagation. Def. 12 states the control flow may take either path of a conditional statement without any consideration for the condition. This is an overapproximation because any possible execution of the program necessarily follows one of the branches. Therefore, this abstraction is sound. D. A Taxonomy of Floating-Point Pitfalls In this section, we formally define absorption and catastrophic cancellation as constraints. These numerical issues arise specifically in the context of addition and subtraction operations. We recall that addition and subtraction are treated in the same way in our analysis. An addition or subtraction operation may exhibit two types of absorption: either the left-hand-side operand or the right-hand-side operand may be absorbed. Definition 13 (Absorption). ∀xv , yv ∈ V, xv is absorbed by yv in the computation xv ± yv if and only if ufp(x̃) ≤ ulp(ỹ) − 1 . Def. 13 defines the condition under which xv is absorbed by yv . Intuitively, this occurs when the magnitude of x (measured by its ufp) is smaller than the resolution of y, i.e., the magnitude of a rounding unit at y given by ulp(ỹ). Swapping the roles of xv and yv yields the symmetric condition corresponding to the absorption of yv by xv . We now formalize our definition of catastrophic cancellation. Definition 14 (Catastrophic Cancellation). Let xv , yv , zv ∈ V be three variables such that zv = xv ± yv . zv is affected by a catastrophic cancellation if and only if ufp(z̃) ≤ ufp(x̃) − nsbσ (xv ) ∨ ufp(z̃) ≤ ufp(ỹ) − nsbσ (yv ) . Def. 14 defines the conditions under which a catastrophic cancellation occurs during the addition or subtraction of xv and yv . Intuitively, this occurs when the magnitude of the result (measured by its ufp) is smaller than the magnitude of the error in either of the operands (measured by its nsb). Property 1 (Detection of same sign arguments). If the errors are sufficiently small, no spurious catastrophic cancellation is detected for same sign additions, or different sign subtraction. Proof. Let us consider zv = xv +yv with x̃ ≥ 0 and ỹ ≥ 0. By symmetry, let us focus on xv and assume nsbσ (xv ) > 0. Since ỹ ≥ 0, we have x̃ ≤ x̃ + ỹ. By monotonicity of rounding, x̃ = ↑p∼ (x̃) ≤ ↑p∼ (x̃ + ỹ) = z̃. Moreover, x̃ > |εx | and therefore z̃ > |εx |. Using Eq. (3) and the definition of nsb, 2ufp(z̃)+1 > 2ufp(x̃)−nsbσ (xv ) and thus ufp(z̃) + 1 > ufp(x̃) − nsbσ (xv ) which means that the formula used to characterize catastrophic cancellation is false. E. Detection of Floating-Point Pitfalls Using an SMT solver FPScan constructs a set of constraints to be checked using an SMT solver. Each constraint abstracts the set of reachable machine states exhibiting a specific numerical pitfall. One
File "..." => z = x + y: ABSORPTION: y absorbed by x File "..." => __result__ = z - x: CANCELLATION ABSORPTION: z absorbed by x Fig. 5. FPScan Output for Figure 1.
constraint is generated for each pitfall that may occur in the program. Each constraint consists of two components. The first characterizes the set of reachable machine states and is derived using C # . The second identifies machine states exhibiting a specific pitfall, as defined in Defs. 13 and 14. Again, we consider the program illustrated in Figure 1. Line 1 may be affected by three numerical pitfalls: absorption of x by y, absorption of y by x, and catastrophic cancellation. Each case is analyzed independently. We first consider the absorption of x by y. FPScan constructs a constraint characterizing the machine states reachable after Line 1 and satisfying the absorption condition. This constraint is obtained by combining Defs. 4 and 13. The resulting formula is submitted to the Z3 [11] SMT solver. Since the constraint is unsatisfiable, this case is proven impossible. The same procedure is applied to the absorption of y by x. In this case, the constraint is satisfiable, and FPScan reports a warning indicating that the pitfall may occur. The same analysis is repeated for catastrophic cancellation. Line 2 is treated in the same way. Note that the set of reachable machine states after Line 2 is constrained by both Line 1 and Line 2. The report produced by FPScan is shown in Figure 5. It was generated in less than 0.15 seconds. It contains two warnings on Line 2 and one warning on Line 1, which implies that the remaining three potential pitfalls are proven impossible. In particular, catastrophic cancellation cannot occur on Line 1. The report highlights the root causes of numerical inaccuracies in Figure 5. For large values of x, y is absorbed on Line 1, introducing a small error and, more importantly, a semantic gap that is later revealed by cancellation on Line 2. However, the absorption reported on Line 2 is spurious, as it would imply that z becomes larger than x, which is impossible. This issue arises from the non-relational nature of the constraints used to bound ufp. This is expected due to the sound yet incomplete nature of our analysis. We will see in the experimental evaluation that the conservativeness of the method is limited, giving very accurate results. IV. E XPERIMENTAL E VALUATION We now investigate the following research questions by evaluating FPScan on the FPBench benchmark suite [5]: • RQ1: To what extent can floating point pitfalls be detected using an abstraction based on ufp, ulp, and nsb? • RQ2: How does the runtime performance of our constraint solving approach compare to that of existing stateof-the-art tools?
FPScan: pitfalls FPScan: guaranteed impossible
15 10
7
5 5 0
4 4 1
1
2
3
18
FPChecker: guaranteed possible FPChecker: missed possible
14 14
10
8
1
2
1 1 1 1 1
2
1 1
2
4
6 3 3 3
4
8
6 2
1
6 6 6
7 7
3
na nt nt2
tes te t06 tesst06_sum btrian t0 _su s sp gle tes5_no ms44__sluine3 t0 n _ m tes5_nloin1___sum1 t03 nlin tes 2 tes _no 1__t2 test02 nlinr4 t01_su 2 _ m nosum83 n nlin veonlin1 pre x_rhu 2 da by_lst t hyorPrexy s int ec pot y ro_ int 4_e hy 32 ex ro_ xa po am ex mp t ple am le _ p sq mixele sq ua sq rt_add reR ua sq d oo reR root cat3Invoot3 rb a turonGlaid turbine s turbine1 kebine23 p ca ler0 f v1 flolouda0 floudas s ud 1 a do sus3 dopple m pp r1 m maatr rigdoppler2 tri ixD himidBoler3 xD et m dy eteermilbe 1 rm in au i a
Catastrophic Cancellations
20
Benchmarks
Fig. 6. Comparison of catastrophic cancellation detection by FPScan and FPChecker. The FPScan bar is split into cancellations detected by FPScan and cancellations proven impossible by FPScan, while the FPChecker bar is split into cancellations detected by FPChecker and occurring cancellations missed by FPChecker.
A. Experimental Setup
processor and 16 GB of memory.
To answer RQ1, we compare FPScan against a ground truth produced by a sound and complete bitblasting-based tool that we developed for this purpose. Bitblasting is a technique used by SMT solvers to reason about floating-point arithmetic by encoding floating-point values as bit vectors, possibly after dedicated optimizations [12]. Although sound and complete for decidable theories, this approach remains too computationally expensive for large-scale practical use. Our implementation relies on Z3 [11] and uses the following detection criteria: (i) catastrophic cancellation is reported when all bits of the result are lost, and (ii) absorption is reported when the least significant bit of one operand is smaller than the most significant bit of the other. This definition of catastrophic cancellation is stricter than the one implemented by FPScan. Consequently, some genuine cancellations may be incorrectly tagged as spurious by bitblasting. Although this does not undermine the overall quality of the evaluation, it may lead to an underestimation of FPScan’s completeness. We evaluate FPScan on a subset of the FPBench benchmark suite [5], which contains 130 numerical programs drawn from numerical analysis papers and textbooks; computations are performed in binary32 format. We retain 58 benchmarks whose preconditions are interval based and compatible with our tool. Among them, 41 are also compatible with bitblasting. We also compare FPScan with FPChecker [6], [7], a stateof-the-art dynamic tool for detecting floating-point catastrophic cancellations. FPChecker reports a cancellation when the number of significant bits in an addition or subtraction falls below 10. Although FPChecker also detects other numerical issues, such as overflow and underflow, it does not detect absorption. We therefore compare FPChecker and FPScan only with respect to cancellation detection. Since FPChecker requires concrete inputs, we sample each input interval using 1000 points. All experiments were conducted on a machine running Ubuntu 22.04.5, equipped with an Intel i7 1370P
B. RQ1: Soundness Validation and Completeness Evaluation. Table I reports the comparison between FPScan and the bitblasting-based ground truth. We classify all detected pitfalls into four categories: true positives (TP), corresponding to pitfalls reported as possible by FPScan that indeed occur; true negatives (TN), corresponding to pitfalls proven impossible by FPScan that do not occur; false positives (FP), corresponding to pitfalls reported as possible by FPScan that are actually impossible; and false negatives (FN), corresponding to pitfalls proven impossible by FPScan that actually occur. FPScan is sound by construction; consequently, the FN columns (highlighted in light gray) in Table I contain only zeros. This means that every pitfall proven impossible by FPScan is indeed impossible. FPScan is also effective at proving the absence of pitfalls: out of 665 checked cases, 482 (72%) are proven impossible. Among the 41 benchmark analyses, 36 contain no false positives. Overall, only 14 (8%) of the 183 reported pitfalls are false positives, indicating that the number of spurious reports is relatively limited. This incompleteness is a consequence of FPScan’s abstraction of machine states: it includes states that are not actually reachable. Consequently, FPScan may fail to prove that some non-occurring pitfalls cannot occur. These false positives are reported in the FP columns of Table I. Figure 6 shows a comparison of the detection capabilities of FPScan and FPChecker. FPChecker fails to detect any catastrophic cancellations on some programs, due to its sampling strategy, which is effective for certain programs but less suitable for others. Addressing this limitation by designing a sampling strategy that is robust across programs is left for future work. Additionally, Figure 6 graphically illustrates the accuracy of FPScan. For most programs, the blue portion of the left bar, representing the catastrophic cancellations detected by FPScan, closely matches the right bar, which represents all occurring catastrophic cancellations.
102
FPScan Bitblasting FPChecker
Time (s)
101
100
10 1
int
ro_ ex a tes intmplbespli t05 ro_ _m ne3 _no exa ixe nli nomp d n1 nl le __ in tes x_btest21 t te 0 flo y_x sqst05squ3_nouday ua _n are nli s reR on Ro n2 oo lin1 ot3 t3I __ h nva r4 pre ypot lid da hy 32 tes vteorPpr ot t02 rh ey u tes d _su lst test06 dooppl m8 t06_su d pp er1 _sums oppler3 ms4__ ler 4__sum2 su 2 ca m1 sec trianv10 4_enon gle x li fl am n2 tes oudaple t s rig01_s sum3 ca idBoum3 himrbo dy mi nGa1 lb s sm sq eau art roo sq Ro t turrt_adot ma turbined t marixD turbbine3 t e i 1 po rixDtermkeplne2 l a po rTo eter inaer0 lar Ca mi nt2 ToC rt flou na arthes da nt he ian s1 s _ aziian___y mu x NMtes SE t04 d sinth _ex _dq el e a m ta ins tan rigmpleom94 tan idB _3. eo ke ody1 us pl 2 Cu er sprren2t he re NM si kedelti6 SE trianeO ple a _ex ng rd r1 am le_ er3 ple tes _3. t 6
10 2
Benchmarks
Fig. 7.
Computation time (s) for FPScan, FPChecker (1000 point sampled) [6], [7], and bitblasting.
Figure 6 also highlights a potential synergy between FPScan and FPChecker. For all programs shown in bold, the catastrophic cancellations detected by FPScan and FPChecker exactly match. This is particularly interesting because, in these cases, catastrophic cancellations are divided into two categories: they are either proven impossible by FPScan or guaranteed possible by FPChecker through a counterexample. The two tools provide complementary information. C. RQ2 Execution Time Evaluation We measure the execution time of FPScan, FPChecker and bitblasting on the 58 benchmarks in the FPBench subset. The results are shown in Figure 7. Programs are ordered by increasing bitblasting execution time. Some programs have no gray bar because bitblasting either timed out or does not support some required operations. Due to the large variation in computation time, a logarithmic scale is used for the vertical axis. The execution time of FPChecker is proportional to the number of samples used in this experiment, namely 1000. For reference, its execution time is comparable to FPScan with 100 samples. Three programs exhibit significantly higher execution times, namely triangle, smartRoot, and instantaneousCurrent. This is due to the interaction between the iterative over-approximation of square root [13] and the use of rational numbers in the interval analysis. Despite issuing one SMT solver query per pitfall check, our method remains efficient. It is comparable to FPChecker and significantly outperforms bitblasting in terms of execution time. This efficiency stems from the simple structure of the constraint system, which involves integer variables, linear operations, min and max functions, and a limited number of conjunctions. Although scalability remains a concern, as each statement increases the constraint size and may trigger up to three SMT solver queries, the approach still demonstrates strong practical performance on the evaluated benchmarks. Future work will investigate strategies to improve
scalability and reduce the number of solver queries per statement. V. R ELATED W ORK A large body of work focuses on bounding rounding errors [14]–[19] or optimizing floating-point expressions [2], [20]–[22]. All these works focus on the characterization of the error, usually computing its upper bound. This is not the focus of this paper. Our primary objective is to formally prove the absence of floating-point pitfalls that can cause severe precision loss. This can help programmers identify weaknesses of floating-point expressions. A. ufp, ulp, and nsb constraint-based analysis Martel [23] gives definitions of ufp and ulp that translate order of magnitude of floating-point numbers and rounding errors into integer values. He uses them to compute errors of constants to propagate them in abstract domains. We rely on a similar approach by formalizing ufp and ulp as constraints and extending their definitions to handle denormalized numbers. Ben Khalifa et al. [21] provide a sound constraint-based precision-tuning tool. They leverage the definition of ufp and ulp by Martel [23] to establish constraints linking the precisions of all variables of a program. Solving these constraints yields a set of intermediate variable precisions that make the output meet a target precision. We rely on a similar approach: producing constraints to model program behavior, but our constraint set is completely different. We also develop a new constraint set including nsb and refining error propagation. In addition, rather than targeting the satisfiability of these constraints, we target their unsatisfiability to prove the absence of pitfalls. B. Detection of Floating-Point Pitfalls We can divide the state-of-the-art on the detection of floating-point pitfalls in two categories: dynamic and static methods. The former execute the programs, yielding unsound
TABLE I C OMPARISON OF FPS CAN WITH THE BITBLASTING BASED GROUND TRUTH .
0 triangle bspline3 0 test06 sums4 sum(x2) 12 0 test05 nonlin1 test2 test05 nonlin1 r4 0 test03 nonlin2 2 test02 sum8 0 4 test01 sum3 nonlin1 1 nonlin2 0 0 verhulst x by xy 0 predatorPrey 0 0 hypot(x2) sec4 example 0 intro example 1 0 intro example mixed smartRoot 0 sqrt add 0 sqroot 4 squareRoot3(x2) 6 carbonGas 1 1 turbine1 turbine2 1 1 turbine3 kepler0 5 3 cav10 floudas 2 floudas1 20 4 floudas3 sum 4 9 doppler(x3) 6 rigidBody1 himmilbeau 8 matrixDeterminant(x2) 20
Cancellation FPScan
Absorption FPScan
Program
FP TP FN TN FP TP FN TN 0 0 0 10 0 0 0 0 5 0 0 0 2 0 0 0 0 1 0 12 0 16 6 0 6 0 2 0 0 0 2 0 0 0 0 1 0 0 0 4 0 0 0 0 2 0 2 0 6 1 0 1 0 2 0 0 0 14 0 0 0 0 7 1 3 0 15 5 2 3 0 3 0 1 0 2 0 0 0 0 1 0 0 0 4 0 0 0 0 2 0 0 0 2 0 0 0 0 1 0 0 0 2 0 0 0 0 1 0 0 0 2 0 0 0 0 1 0 0 0 4 0 0 0 0 2 0 0 0 4 0 0 0 0 2 0 1 0 2 0 0 0 0 1 0 0 0 2 0 0 0 0 1 0 0 0 12 0 0 0 0 6 0 0 0 4 0 0 0 0 2 0 4 0 8 0 0 0 0 4 0 6 0 12 2 0 2 0 4 0 1 0 6 0 0 0 0 3 0 1 0 12 1 0 1 0 5 0 1 0 8 2 0 2 0 2 0 1 0 12 1 0 1 0 5 4 1 0 16 6 2 4 0 4 1 2 0 3 1 0 1 0 1 0 2 0 2 1 0 1 0 0 0 20 0 36 17 0 17 0 1 1 3 0 5 2 0 2 0 1 1 3 0 15 5 2 3 0 3 0 9 0 36 0 0 0 0 18 0 6 0 14 3 0 3 0 4 0 8 0 14 5 0 5 0 2 0 20 0 56 10 0 10 0 18
Programs annotated with “(x2)” or “(x3)” denote groups of two or three similar benchmarks. FPScan col. contains the number of detected pitfalls.
yet complete results. The latter reason about the programs yielding sound yet incomplete results. a) Dynamic Methods: Most existing approaches are dynamic. The program is executed on a set of inputs. These methods are often fast and scalable to large and complex programs. However, they are in practice computationally expensive, both in terms of execution time and memory consumption, especially when a large number of inputs must be explored. Moreover, their results are only valid for the explored inputs, with no guarantees over the rest of the input space. In this context, Lam et al. [8], Benz et al. [24], Jézéquel et al.[19], and Laguna et al. [6], [7] propose similar approaches. Their tools analyze programs to estimate the number of bits lost after each addition and subtraction. A catastrophic cancellation is reported when this number drops below a given threshold. Ben Khalifa et al. [25] rely on shadow execution at a higher precision to approximate floating-point errors. They distinguish two types of cancellation: catastrophic cancellation,
which occurs when the magnitude of the accumulated error becomes comparable to or larger than the magnitude of the computed result during an addition or subtraction, and benign cancellation, where only part of the significant bits are lost. We formalized their characterization of catastrophic cancellation as a constraint in our approach. Overall, these dynamic techniques target cancellations rather than absorptions. They achieve good performance, but are not sound unless the entire input space is explored, which is generally infeasible in practice. Our method uses a constraint-based abstraction of reachable machine states to soundly reason over all executions and prove the absence of both catastrophic cancellations and absorptions. b) Static Methods: Static methods infer the behavior of a program over all possible inputs by analyzing its code rather than executing it. The goal is to efficiently compute properties that hold for all program executions. To the best of our knowledge, only Lopes et al. [26] address the sound detection of both absorption and catastrophic cancellation. Their approach first performs a range analysis to bound all intermediate variables, and then applies decision rules to determine whether a numerical pitfall may occur. In their definition, an absorption occurs when the difference between the result of an operation and one of the operands is smaller than the rounding error. However, the paper does not discuss denormalized numbers, so it is difficult to assess whether the method extends to them soundly. Moreover, the treatment of loops is not described in sufficient detail to determine how pitfall conditions are propagated across iterations. Our approach addresses these issues by supporting denormalized values, soundly analysing loops up to a given unrolling bound, and taking rounding errors into account in the catastrophic cancellation detection criterion. However, we could not directly compare our approach with their tool because it is not freely available. VI. C ONCLUSION We presented FPScan, a novel tool and approach to formally detect both catastrophic cancellation and absorption in floating-point programs. Our approach combines a range analysis using abstract interpretation with an axiomatization as a set of logical constraints of the propagation of numerical errors in the program. Proving the absence of floating-point pitfalls amount to prove that a formula is unsatisfiable with an SMT solver. Experimental results show that FPScan is both precise and efficient, proving 72% of cases impossible with only 8% false positives, while outperforming bitblasting and comparing favorably with FPChecker. Future work will extend FPScan to mixed-precision programs. We also plan to rely on identified pitfalls to propose guaranteed repairs, for example, using precision tuning [21], expression transformations [2], and Taylor expansions [2]. R EFERENCES [1] “IEEE Standard for Floating-Point Arithmetic,” IEEE Std 754-2019 (Revision of IEEE 754-2008), pp. 1–84, 2019.
[2] P. Panchekha, A. Sanchez-Stern, J. R. Wilcox, and Z. Tatlock, “Automatically improving accuracy for floating point expressions,” in Proceedings of the 36th ACM SIGPLAN Conference on Programming Language Design and Implementation, ser. PLDI ’15. New York, NY, USA: Association for Computing Machinery, 2015, pp. 1–11. [3] D. Goldberg, “What every computer scientist should know about floating-point arithmetic,” ACM Computing Surveys, vol. 23, no. 1, pp. 5–48, Mar. 1991. [4] E. Misback, C. C. Chan, B. Saiki, E. Jun, Z. Tatlock, and P. Panchekha, “Odyssey: An interactive workbench for expert-driven floating-point expression rewriting,” in Proceedings of the 36th Annual ACM Symposium on User Interface Software and Technology, ser. UIST ’23. New York, NY, USA: Association for Computing Machinery, 2023. [Online]. Available: https://doi.org/10.1145/3586183.3606819 [5] N. Damouche, M. Martel, P. Panchekha, C. Qiu, A. Sanchez-Stern, and Z. Tatlock, “Toward a standard benchmark format and suite for floatingpoint analysis,” in International Workshop on Numerical Software Verification. Springer, 2016, pp. 63–77. [6] I. Laguna, “FPChecker: Detecting Floating-Point Exceptions in GPU Applications,” in 2019 34th IEEE/ACM International Conference on Automated Software Engineering (ASE), Nov. 2019, pp. 1126–1129, iSSN: 2643-1572. [7] I. Laguna, T. Tirpankar, X. Li, and G. Gopalakrishnan, “FPChecker: Floating-Point Exception Detection Tool and Benchmark for Parallel and Distributed HPC,” in 2022 IEEE International Symposium on Workload Characterization (IISWC), Nov. 2022, pp. 39–50. [8] M. O. Lam, J. K. Hollingsworth, and G. W. Stewart, “Dynamic floatingpoint cancellation detection,” Parallel Comput., vol. 39, no. 3, pp. 146– 155, 2013. [9] L. Titolo, M. M. Moscato, M. A. Feliú, P. Masci, and C. A. Muñoz, “Rigorous floating-point round-off error analysis in precisa 4.0,” in Formal Methods - 26th International Symposium, FM 2024, Milan, Italy, September 9-13, 2024, Proceedings, Part II, ser. Lecture Notes in Computer Science, A. Platzer, K. Y. Rozier, M. Pradella, and M. Rossi, Eds., vol. 14934. Springer, 2024, pp. 20–38. [Online]. Available: https://doi.org/10.1007/978-3-031-71177-0 2 [10] P. Cousot and R. Cousot, “Abstract interpretation: a unified lattice model for static analysis of programs by construction or approximation of fixpoints,” in Proceedings of the 4th ACM SIGACT-SIGPLAN symposium on Principles of programming languages, 1977, pp. 238–252. [11] L. De Moura and N. Bjørner, “Z3: An efficient SMT solver,” in International conference on Tools and Algorithms for the Construction and Analysis of Systems. Springer, 2008, pp. 337–340. [12] A. Brillout, D. Kroening, and T. Wahl, “Mixed abstractions for floatingpoint arithmetic,” in 2009 Formal Methods in Computer-Aided Design, 2009, pp. 69–76. [13] M. Daumas, D. Lester, and C. Muñoz, “Verified Real Number Calculations: A Library for Interval Arithmetic,” Aug. 2007, arXiv:0708.3721 [cs.MS]. [Online]. Available: http://arxiv.org/abs/0708.3721 [14] F. de Dinechin, C. Lauter, and G. Melquiond, “Certifying the floating-point implementation of an elementary function using gappa,” IEEE Trans. Comput., vol. 60, no. 2, p. 242–253, Feb. 2011. [Online]. Available: https://doi.org/10.1109/TC.2010.128 [15] M. Moscato, L. Titolo, A. Dutle, and C. A. Muñoz, “Automatic estimation of verified floating-point round-off errors via static analysis,” in Computer Safety, Reliability, and Security, S. Tonetta, E. Schoitsch, and F. Bitsch, Eds. Cham: Springer International Publishing, 2017, pp. 213–229. [16] E. Goubault, M. Martel, and S. Putot, “Asserting the precision of floating-point computations: A simple abstract interpreter,” in Programming Languages and Systems, D. Le Métayer, Ed. Berlin, Heidelberg: Springer Berlin Heidelberg, 2002, pp. 209–212. [17] B. Blanchet, P. Cousot, R. Cousot, J. Feret, L. Mauborgne, A. Miné, D. Monniaux, and X. Rival, “Design and implementation of a specialpurpose static program analyzer for safety-critical real-time embedded software,” The essence of computation: complexity, analysis, transformation, pp. 85–108, 2002. [18] S. Blazy, D. Bühler, and B. Yakobowski, “Structuring abstract interpreters through state and value abstractions,” in International Conference on Verification, Model Checking, and Abstract Interpretation. Springer, 2017, pp. 112–130. [19] F. Jézéquel and J.-M. Chesneaux, “CADNA: a library for estimating round-off error propagation,” Computer Physics Communications,
vol. 178, no. 12, pp. 933–955, Jun. 2008. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0010465508000775 [20] A. Ioualalen and M. Martel, “Sardana: an automatic tool for numerical accuracy optimization,” in SCAN: Scientific Computing, Computer Arithmetic and Validated Numerics, 2012, pp. 1–4. [21] D. Ben Khalifa, M. Martel, and A. Adjé, “Pop: A tuning assistant for mixed-precision floating-point computations,” in Formal Techniques for Safety-Critical Systems, O. Hasan and F. Mallet, Eds. Cham: Springer International Publishing, 2020, pp. 77–94. [22] C. Rubio-González, C. Nguyen, H. D. Nguyen, J. Demmel, W. Kahan, K. Sen, D. H. Bailey, C. Iancu, and D. Hough, “Precimonious: tuning assistant for floating-point precision,” in Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis, ser. SC ’13. New York, NY, USA: Association for Computing Machinery, 2013. [Online]. Available: https://doi.org/10.1145/2503210.2503296 [23] M. Martel, “Floating-point format inference in mixed-precision,” in NASA Formal Methods Symposium. Springer, 2017, pp. 230–246. [24] F. Benz, A. Hildebrandt, and S. Hack, “A dynamic program analysis to find floating-point accuracy problems,” in Proceedings of the 33rd ACM SIGPLAN Conference on Programming Language Design and Implementation, ser. PLDI ’12. New York, NY, USA: Association for Computing Machinery, 2012, p. 453–462. [Online]. Available: https://doi.org/10.1145/2254064.2254118 [25] D. Ben Khalifa, X. Li, I. Laguna, M. Martel, and G. Gopalakrishnan, “Toward Increasing Trust in Exascale Simulations,” in 2022 4th Annual Workshop on Extreme-scale Experiment-in-the-Loop Computing (XLOOP), Nov. 2022, pp. 26–31. [26] M. K. Lopes, R. B. França, C. M. Hirata, and L. A. V. Dias, “Method to detect floating-point absorption and cancellation phenomena in softwarecritical design models,” in 2018 Eighth Latin-American Symposium on Dependable Computing (LADC). IEEE, 2018, pp. 125–134.
A PPENDIX
B. Soundness of abstract multiplication (cf. Def. 5)
In this section, we provide the proofs of the theorems presented in Section III. We first recall the soundness proof for Addition/Subtraction. Note that the proofs are valid in the two edge cases of zero and denormalized numbers. This follows from the relations between rounding error and ulp, and between order of magnitude and ufp established by Eqs.(3) and (5), respectively. Both equations are valid for zero and denormalized numbers. A. Soundness of abstract addition/subtraction (cf. Def. 4)
z
z =x×y
(11)
εz
The error |εz | can be bounded using the triangle inequality as shown
z = (x̃ + εx ) × (ỹ + εy ) z = (x̃ × ỹ) + x̃εy + ỹεx + εx εy z = ↑p∼ (x̃ × ỹ) +(x̃εy + ỹεx + εx εy − e× ) | {z } | {z }
(12)
The errors |εx | and |εy | can be bounded using Eq. (3) and e+ is a simple rounding error; thus, |e+ | ≤ 2ulp(z̃)−1 using Eq. (5). Hence, Eq. (12) can be rewritten as Eq. 13. |εz | < 2ufp(εx )+1 + 2ufp(εy )+1 + 2ulp(z̃)−1 .
(13)
Eq. (13) is the sum of three powers of two. It can thus be bounded as shown in Eq. (14). |εz | < 2max(ufp(εx )+1,ufp(εy )+1,ulp(z̃)−1)+2 .
(17)
εz
z̃
The error |εz | can be bounded using triangle inequality: |εz | ≤ |x̃εy | + |ỹεx | + |εx εy | +|e× | . | {z } | {z } | {z } b2
b1
(18)
b3
The digits b1 , b2 , b3 and the error e× can be bounded independently. Equations (19), (21) and (22) can be obtained using equations (3) and (6). |x̃| < 2ufp(x̃)+1 .
|εz | ≤ |εx | + |εy | + |e+ | .
(16)
(10)
By applying Eq. (6), a bound on nsbσ (zv ) can be derived from a bound on ufp(εz ). This bound, in turn, can be obtained by bounding |εz | using Eq. (3). An explicit expression for εz is obtained through a series of manipulations, as demonstrated in Eq. (11). Here, we denote by e+ the rounding error introduced by the addition. z̃ = ↑p∼ (x̃ + ỹ) = x̃ + ỹ + e+ = (x − εx ) + (y − εy ) + e+ = x + y −(εx + εy − e+ ) | {z } | {z }
εx = x − x̃ ⇔ x̃ = x − εx ⇔ x = x̃ + εx ; εy = y − ỹ ⇔ ỹ = y − εy ⇔ y = ỹ + εy ; εz = z − z̃ ⇔ z̃ = z − εz ⇔ z = z̃ + εz .
Using the same reasoning as in the proof of Section A, the constraint will bound nsbσ (zv ). Eq. (17) shows how to bound εz . We denote by e× the rounding error introduced by the multiplication.
Let εx , εy , and εz denote the floating-point errors: εx = x − x̃ ⇔ x̃ = x − εx ⇔ x = x̃ + εx ; εy = y − ỹ ⇔ ỹ = y − εy ⇔ y = ỹ + εy ; εz = z − z̃ ⇔ z̃ = z − εz ⇔ z = z̃ + εz .
Let σ denotes a machine state and xv , yv , zv denote three variables. x, y and z denotes their respective real values, x = σR (xv ), y = σR (yv ) and z = σR (zv ). x̃, ỹ and z̃ denote their respective floating-point values, x̃ = σFp (xv ), ỹ = σFp (yv ) and z̃ = σFp (zv ).
|ỹ| < 2
ufp(ỹ)+1
(19)
.
(20)
ufp(ỹ)−nsbσ (yv )+1
.
(21)
ufp(x̃)−nsbσ (xv )+1
.
(22)
|εy | < 2
|εx | < 2
Together, equations (19) and (21) allow to bound b1 as shown in Equations (23). Eq. (24) is obtained using a similar reasoning. b1 = |x̃εy | < 2ufp(x̃)+ufp(ỹ)−nsbσ (yv )+2 . b2 = |ỹεx | < 2
ufp(x̃)+ufp(ỹ)−nsbσ (xv )+2
.
(14)
Eq. (25) is obtained using Equations (21) and (22).
We can deduce ufp(εz ) from Eq. (14) using Eq. (3) and thus obtain
b3 = |εx εy | < 2ufp(x̃)+ufp(ỹ)−nsbσ (xv )−nsbσ (yv )+2 .
(23) (24)
(25)
ufp(εz ) ≤ max(ufp(εx )+1, ufp(εy )+1, ulp(z̃)−1)+1 . (15)
Finally, we obtain Eq. (26) using Eq. (5) since e× is a rounding error: |e× | ≤ 2ulp(z̃)−1 . (26)
Eq. (15) together with Eq. (6) allows us to obtain the desired result.
Eq. (18) can be rewritten using equations (23) to (26) to finally obtain Eq. (27).
|εz | <
ufp(x̃)+ufp(ỹ)−nsbσ (yv )+2
2 +2ufp(x̃)+ufp(ỹ)−nsbσ (xv )+2 +2ufp(x̃)+ufp(ỹ)−nsbσ (xv )−nsbσ (yv )+2 +2ulp(z̃)−1
.
(27)
Since Eq. (27) is the sum of four terms, each a power of two, it can be bounded by a single power of two. Let M denote the maximum exponent, as defined in Eq. (28). ufp(x̃) + ufp(ỹ) − nsbσ (yv ) + 2 ufp(x̃) + ufp(ỹ) − nsbσ (xv ) + 2 = max ufp(x̃) + ufp(ỹ) − nsbσ (xv ) − nsbσ (yv ) + 2 ulp(z̃) − 1 ufp(x̃) + ufp(ỹ) − min(nsbσ (xv ), nsbσ (yv )) + 2 = max ufp(x̃) + ufp(ỹ) − nsbσ (xv ) − nsbσ (yv ) + 2 ulp(z̃) − 1
M
.
We obtain Eq. (34) using Eq. (5) since e÷ is a rounding error: |e× | ≤ 2ulp(z̃)−1 . . (34) The error |εz | can be bound using the triangle inequality as shown: ỹεx − x̃εy − e÷ | |εz | = | ỹ(ỹ + εy ) ỹεx − x̃εy ≤| | + |e÷ | ỹ(ỹ + εy ) ỹεx − x̃εy ≤| | + 2ulp(z)−1 . (35) ỹ(ỹ + εy ) To bound |εz |, one therefore needs to find an upper bound of |ỹεx − x̃εy | and a lower bound of |ỹ(ỹ + εy )|.
(28)
It is thus possible to establish a bound on |εz |. Given the maximum exponent M , this upper bound is attained when all four exponents are equal. Under such circumstances, two carry bits may arise, which yields |εz | < 2M +2 . Thus, using Eq. (3) we obtain ufp(εz ) ≤ M + 1 and then Eq. (29) replacing M by its value.
ufp(x̃) + ufp(ỹ) − min(nsbσ (xv ), nsbσ (yv )) + 3 . ufp(εz ) ≤ max ufp(x̃) + ufp(ỹ) − nsbσ (xv ) − nsbσ (yv ) + 3 ulp(z̃)
(29)
|x̃| < 2ufp(x̃)+1 . |ỹ| < 2
ufp(ỹ)+1
(36)
.
(37)
ufp(x̃)−nsbσ (xv )+1
.
(38)
ufp(ỹ)−nsbσ (yv )+1
.
(39)
|εx | < 2 |εy | < 2
Together, Equations (36), (37), (38), and (39) allow us to bound |ỹεx − x̃εy | as shown in Eq. (40). As Eq. (40) is made of two power of two, the upper bound can be refined as shown in Eq. (41).
Equations (29) and (6) allow to conclude. |ỹεx − x̃εy | < 2ufp(x̃)+ufp(ỹ)−nsbσ (xv )+2
■ C. Soundness of abstract division (cf. Def. 6) Let σ denotes a machine state and xv , yv , zv denote three variables. x, y and z denotes their respective real values, x = σR (xv ), y = σR (yv ) and z = σR (zv ). x̃, ỹ and z̃ denote their respective floating-point values, x̃ = σFp (xv ), ỹ = σFp (yv ) and z̃ = σFp (zv ) εx = x − x̃ ⇔ x̃ = x − εx ⇔ x = x̃ + εx ; εy = y − ỹ ⇔ ỹ = y − εy ⇔ y = ỹ + εy ; εz = z − z̃ ⇔ z̃ = z − εz ⇔ z = z̃ + εz .
(30)
+ 2ufp(x̃)+ufp(ỹ)−nsbσ (yv )+2 ufp(x̃)+ufp(ỹ)−min(nsb(xv ),nsbσ (yv ))+3
<2
(41)
|ỹ + εy | ≥ ||ỹ| − |εy || ≥ |ỹ| − |εy | > 2ufp(ỹ) − 2ufp(εy )+1 > 2ufp(ỹ)−1
(42)
Eq. (42) allows to conclude |ỹ(ỹ + εy )| > 22ufp(ỹ)−1 and then to get Eq. (43):
εz =
(31) x̃ + εx x̃ − − e÷ ỹ + εy ỹ ỹ(x̃ + εx ) − x̃(ỹ + εy ) = − e÷ ỹ(ỹ + εy ) ỹεx − x̃εy = − e÷ . ỹ(ỹ + εy )
.
Now, let us focus on a lower bound for |ỹ(ỹ + εy )|. The main issue is that ỹ + εy can be close to 0 for large errors. We limit the analysis to scenarios with small errors, i.e. when ufp(ỹ) > ufp(εy )+1 which can be simplified to nsbσ (yv ) > 1. Otherwise, no constraint is applied. To use Equations (37) and (39) allows to get Eq. (42)
Using the same reasoning as in the soundness proof of Section A, the constraint will bound nsbσ (zv ). Eq. (33) shows how to bound εz . We denote by e÷ the rounding error introduced by the division. x̃ x − ↑p∼ ( ) y ỹ x x̃ = − − e÷ y ỹ
(40)
|
2ufp(x̃)+ufp(ỹ)−min(nsbσ (xv ),nsbσ (yv ))+3 ỹεx − x̃εy |< ỹ(ỹ + εy ) 22ufp(ỹ)−1 < 2ufp(x̃)−ufp(ỹ)−min(nsbσ (xv ),nsbσ (yv ))+4 . (43)
=
Equations (35) and (43) can be put together to obtain Eq. (44): (32) (33)
|εz | < 2ufp(x̃)−ufp(ỹ)−min(nsbσ (xv ),nsbσ (yv ))+4 + 2ulp(z̃)−1 . (44)
Eq. (44) is the sum of two power of two and can thus be bounded as shown in Eq. 45: max
|εz | < 2
how to bound εz . We denote by e the rounding error introduced by the operation.
ufp(x̃) − ufp(ỹ) − min(nsbσ (xv ), nsbσ (yv )) + 4 +1 ulp(z̃) − 1
. (45)
Which finally allow to conclude with Eq. (46) If nsbσ (yv ) > 1 then ufp(x̃) − ufp(ỹ) − min(nsbσ (xv ), nsbσ (yv )) + 4 ufp(εz ) ≤ max ulp(z̃) − 1
(46) ■
D. Soundness of abstract Cos, Sin, and Atan (cf. Def. 7) Let εx , and εz denote the floating-point errors: εx = x − x̃ ⇔ x̃ = x − εx ⇔ x = x̃ + εx ; εz = z − z̃ ⇔ z̃ = z − εz ⇔ z = z̃ + εz .
(47)
Using the same reasoning as in the soundness proof of Section A, the constraint will bound nsbσ (zv ). Let f denote either sin, cos, or atan. Let z̃ = ↑p∼ (f (x̃)), z = f (x), and e a rounding error.
εz = z − z̃ √ √ = x − ↑p∼ ( x + εx ) √ √ = x − x − εx − e √ √ |εz | ≤ | x − x + εx | +|e| . {z } | ∆
p As stated by Moscato et al. [15], ∆ ≤ |εx |. This result can be obtain computing the derivative of ∆ with respect to x. Using Eq. (3) bounding the ufp together with Def. 3 defining the nsb, we obtain |εx | < 2ufp(x̃)−nsbσ (xv )+1
z
R
∆≤
p |εx |
|R| ≤ max|f ′ | |εx | ≤ |εx | < 2ufp(x)−nsbσ (x)+1 .
(49)
We denote by e a simple rounding error thus |e| ≤ 2ulp(z̃)−1 using Eq. (5). Equations (48) and (49) can be used to bound |εz |: |εz | = |R + e| ≤ |R| + |e| < 2ufp(x)−nsbσ (x)+1 + 2ulp(z)−1 < 2max(ufp(x)−nsbσ (x)+1, ulp(z)−1)+1 .
(50)
Equations (3) and (50) allow us to bound ufp(εz ) as shown ufp(εz ) ≤ max
ufp(x) − nsbσ (x) + 1 ulp(z) − 1
.
(51)
■ E. Soundness of abstract Sqrt (cf. Def. 8) Let σ denotes a machine state and xv , zv denote two variables. x, and z denotes their respective real values, x = σR (xv ), and z = σR (zv ). x̃ and z̃ denote their respective floating-point values, x̃ = σFp (xv ), and z̃ = σFp (zv ). εx = x − x̃ ⇔ x̃ = x − εx ⇔ x = x̃ + εx . εz = z − z̃ ⇔ z̃ = z − εz ⇔ z = z̃ + εz .
(52)
Using the same reasoning as in the soundness proof of Section A, the constraint will bound nsbσ (zv ). Eq. (53) shows
ufp(x̃)−nsbσ (xv )+1 2
.
(55)
Eq. (5) allows us to bound the rounding error e using ulp(z̃). |e| ≤ 2ulp(z̃)−1 .
εz
We denote by R the remainder that can be bounded using Taylor-Lagrange Inequality as shown in Eq. (49). Note that, since f is either cos, sin, or atan, max|f ′ | = 1.
(54)
that allows to bound ∆:
<2 z̃ = ↑p∼ (f (x̃)) = f (x̃) + e = f (x − εx ) + e = f (x) + (f (x − εx ) − f (x)) + e = f (x) − (−R − e) . (48) | {z } |{z} | {z }
(53)
(56)
Combining Eq. (55) and (56) leads to |εz | < 2
ufp(x̃)−nsbσ (xv )+1 2
+ 2ulp(z̃)−1
(57)
which is the sum of two power of two, hence: |εz | < 2max(
ufp(x̃)−nsbσ (xv )+1 ,ulp(z̃)−1))+1 2
(58)
Eq. (58) bounds |εz | using a power of two. Therefore, we use Eq. (3) to deduce a bound on ufp(εz ): ufp(x̃) − nsbσ (xv ) + 1 , ulp(z̃) − 1) (59) 2 Eq. (59) combined with Def. (3) allows us to conclude. ■ ufp(εz ) ≤ max(