CONFIDENTIAL Precise approximation of Floating-Point Computations - - PowerPoint PPT Presentation
CONFIDENTIAL Precise approximation of Floating-Point Computations for C/C++ Software Using the Mathematical Libraries Roberto Bagnara (BUGSENG, UNIPR) Michele Chiari (UNIPR, BUGSENG, POLIMI) Roberta Gori (UNIPI), Abramo Bagnara (BUGSENG)
Precise approximation of Floating-Point Computations for C/C++ Software Using the Mathematical Libraries
Roberto Bagnara (BUGSENG, UNIPR) Michele Chiari (UNIPR, BUGSENG, POLIMI) Roberta Gori (UNIPI), Abramo Bagnara (BUGSENG)
http://bugseng.com
IMDEA, Madrid, Spain — November 22nd, 2019
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Introduction
Floating-Point Numbers in Critical Systems
The use of floating-point computations for the implementation of critical systems is perceived as increasingly acceptable Even in modern avionics floating-point numbers are increasingly used instead of fixed-point arithmetic Three main goals of floating-point program verification:
1
Ensure absence of undefined or undesirable behaviors: infinities, NaNs, overflows in conversion to integers. . .
2
Identifying sources of round-off errors
3
Prove that the program implements a certain numerical computation up to some specified error bound Our work is concerned with the first and the second goal
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Introduction
The Issue
The IEEE 754 Standard for Floating-Point Arithmetic offers, to some extent, a predictable, reliable and well-defined numerical environment Compliance with significant parts of it is increasingly spread among software and hardware platforms The state of the art of software testing and verification techniques allows to tackle the most basic operations The same cannot be said about the mathematical library functions
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Introduction
C/C++ libm Functions
C and C++ offer a floating-point mathematical library, accessible via the <math.h>/<cmath> header files, with functions such as exp, log, sin, cos, tan, sinh, cosh, . . . Problem: they come with very limited guarantees on the computed results “Therefore many of the functions in the math library have
which is exposed by one of the existing tests in the test
the actual maximum error (or at least a ballpark figure) but this is often not achieved due to the large search space.” (excerpted from the GNU libc documentation)
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Introduction
A Real-World Example
1 #include <math.h> 2 #include <stdint.h> 3 4 // Omitted macros CScal, CAdd, CSub... 10 #define CI(z) { float tmp = z.re; z.re = - z.im; z.im = tmp; } 11 #define CExp(z) { float e = exp(z.re); z.re = e*cos(z.im); \ 12 z.im = e*sin(z.im); } 13 #define CSin(z) { CI(z); struct complex _z = {-z.re, -z.im}; \ 14 float e = exp(z.re); float cos_z_im = cos(z.im); z.re = e*cos_z_im; \ 15 float sin_z_im = sin(z.im); z.im = e*sin_z_im; _z.re = cos_z_im/e; \ 16 _z.im = -sin_z_im/e; CSub(_z, z); CScal(-0.5, z); CI(z); } 17 18 static inline float isometric_latitude(float phi, float e) { 19 return log(tan(M_PI_4 + phi / 2.0)) 20
21 } 22 23 static inline float isometric_latitude0(float phi) { 24 return log(tan(M_PI_4 + phi / 2.0)); 25 }
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Introduction
A Real-World Example (cont’d)
27 void latlong_utm_of(float phi, float lambda, uint8_t utm_zone) { 28 float lambda_c = LambdaOfUtmZone(utm_zone); 29 float ll = isometric_latitude(phi, E); 30 float dl = lambda - lambda_c; 31 float phi_ = asin(sin(dl) / cosh(ll)); 32 float ll_ = isometric_latitude0(phi_); 33 float lambda_ = atan(sinh(ll) / cos(dl)); 34 struct complex z_ = { lambda_, ll_ }; 35 CScal(serie_coeff_proj_mercator[0], z_); 36 uint8_t k; 37 for(k = 1; k < 3; k++) { 38 struct complex z = { lambda_, ll_ }; 39 CScal(2*k, z); 40 CSin(z); 41 CScal(serie_coeff_proj_mercator[k], z); 42 CAdd(z, z_); 43 } 44 CScal(N, z_); 45 latlong_utm_x = XS + z_.im; 46 latlong_utm_y = z_.re; 47 }
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Introduction
Questions
What do the mathematical functions used in the code compute? More specifically:
1
Can infinities and NaNs be generated?
2
Can sin, cos, and tan be invoked on ill-conditioned arguments?
3
If anomalies of any one of these kinds are possible, which inputs may cause them? Ill-conditioned arguments: if x is an IEEE 754 single- precision number and x ≥ 223, then the smallest single- precision range containing [x, x + 2π) contains no more than three floating-point numbers
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Introduction
Can We Answer such Questions Statically?
Lack of specification for the libm functions would seem to prevent any attempt at formal verification of programs that use them Without a specification for the functions, no conclusion can be drawn statically about the behavior of the program We propose an alternative to surrender: If we do not have a proper specification for the functions, we might have a partial specification At the very least we have their implementation, for which a partial specification can sometimes be extracted When even a partial specification is unavailable, we can still detect flaws in the program and facilitate its partial verification with testing
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Solving Constraints on libm Functions
libm constraints
Let f : F → F be a libm mathematical function: we are interested in solving constraints of the form y = f (x), We associate an interval domain to each variable: x : [xl, xu] ∈ IF y : [yl, yu] ∈ IF Then we iteratively refine the intervals by means of projections: Direct Projection y ← x Inverse Projection y → x
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Solving Constraints on libm Functions
Projection Functions
Direct Projection Determine [y ′
l , y ′ u] ⊂ [yl, yu] such that
∀x ∈ [xl, xu] : f (x) ∈ [yl, yu] = ⇒ f (x) ∈ [y ′
l , y ′ u]
Inverse Projection Determine [x′
l , x′ u] ⊂ [xl, xu] such that
∀x ∈ [xl, xu] : f (x) ∈ [yl, yu] = ⇒ x ∈ [x′
l , x′ u]
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Solving Constraints on libm Functions
Exact Projections
Exact projection functions are known in literature: they require correctly rounded implementations of f and f −1 No implementation of libm we know of, among those actually used in the embedded software industry, is correctly rounded The IEEE 754 Standard only recommends correct rounding The C/C++ standards give no prescriptions The POSIX Standard requires correctness for special values, e.g.: log(±0.0) = −∞ −1.0 ≤ sin(x) ≤ +1.0 log(−1) = NaN asin(−2.0) = NaN GNU libc provides a “ballpark figure” of the maximum error
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Solving Constraints on libm Functions
Quasi-Monotonicity To the Rescue
However, one important property of the functions is almost preserved: monotonicity Functions are almost piece-wise monotonic, with a few irregularities, that we call glitches They are usually few, shallow and narrow
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Solving Constraints on libm Functions
Glitches
G1 G2
width(G1) = 5, depth(G1) = 4, width(G2) = depth(G2) = 2.
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Solving Constraints on libm Functions
Glitch Data
Our approach only requires safe approximations of: ng ∈ N: the total number of glitches dM ∈ N: the maximal depth of glitches (in ULPs) wM ∈ N: the maximal width of glitches (in ULPs) α ∈ F: the point where the first glitch begins ω ∈ F: the point where the last glitch ends
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Solving Constraints on libm Functions
Glitch Data (cont’d)
function Dmin DM near up down zero ng dM wM ng dM wM ng dM wM ng dM wM acosf −1 1 acoshf 1 ∞ 1 1 2 1 1 2 asinf −1 1 asinhf −∞ ∞ 2 1 2 2 1 2 atanf −∞ ∞ 1 1 108 atanhf −1 1 2 1 2 2 1 2 cbrtf −∞ ∞ 106 1 2 106 1 2 106 1 2 106 1 2 coshf −∞ ∞ 456 1 2 462 1 2 442 1 2 448 1 2 erff −∞ ∞ expf −∞ ∞ exp10f −∞ ∞ exp2f −∞ ∞ 2 1 2 expm1f −∞ ∞ lgammaf 2 ∞ logf ∞ log10f ∞ log1pf −1 ∞ 1 1 2 1 1 2 log2f ∞ sinhf −∞ ∞ sqrtf ∞ tanhf −∞ ∞ 1 1 2 2 1 3 tgammaf 2 ∞ 104 2 3 104 4 4 104 2 3 104 3 4 cosf −223 223 104 1 3 104 1 3 104 1 3 104 1 3 sinf −223 223 tanf −223 223 Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Solving Constraints on libm Functions
Direct Projection Outside Glitch Region
If [xl, xu] ∩ [α, ω] = ∅, we can take [y ′
l , y ′ u] ≡
xu f (xl) f (xu)
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Solving Constraints on libm Functions
Direct Projection with Glitches
There might be a glitch right after xl, where the function is lower than f (xl): if xl ≻ α, take yl = preddM(f (xl)) if xl ≺ α, take yl = min
somewhere before xu: if xu ≻ ω, take yu = f (xu) if xu ≺ ω, take yu = succdM(f (xu))
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Solving Constraints on libm Functions
Direct Projection with Glitches (cont’d)
xl xu α f (xl) f (xu) yl yu preddM(f (α)) succdM(f (xu))
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Solving Constraints on libm Functions
Inverse Projection
We separately solve yl ≤ f
[xl, xu] An adapted dichotomic search algorithm is used rough inverse f i is used to find a starting point: e.g., exp is used as a rough inverse of log and the other way around Glitch data are used to mitigate the quasi-monotonicity of the functions, so that the bisection process can continue even in the presence of glitches The more precise the glitch data, the quicker is overall convergence Crucial point: logarithmic complexity is maintained
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Solving Constraints on libm Functions
Trigonometric Functions
sin, cos and tan change their monotonicity periodically in multiples of π
2
How to distinguish between monotonicity changes and glitches? If we can recognize and generate multiples of π
2, we can split
interval [xl, xu] in quasi-monotonic branches, and locally apply the described projection algorithms Range reduction algorithms can be used for this purpose.
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Solving Constraints on libm Functions
Trigonometric Functions
y = sinf(x)
2
2
2
2
2
2
2
2
Precise approximation of Floating-Point Computations
Solving Constraints on libm Functions
Range Reduction
Given x ∈ F, compute k ∈ Z and r ∈ F with sufficient precision such that x ≈ k π 2 + r We actually just need to compute x 2
π and k π 2 :
for the single precision (32 bit) format, 53 precision bits are sufficient for the double precision (64 bit) format, at least 110 bits are needed
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Solving Constraints on libm Functions
Constraint Solving Example
Take this line:
31 float phi_ = asin(sin(dl) / cosh(ll));
Convert it to SSA:
1 float phi_; double z1, z2, z3, z4, z5, z6; 2 z1 = (double) dl; 3 z2 = sin(z1); 4 z3 = (double) ll; 5 z4 = cosh(z3); 6 z5 = z2 / z4; 7 z6 = asin(z5); 8 phi_ = (float) z6;
Suppose we want to prove z6 cannot be NaN: we try to solve z5 > 1 with dl, ll ∈ {+∞, −∞, nan}.
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Solving Constraints on libm Functions
Constraint Solving Example (cont’d)
We try to solve z5 > 1 with dl, ll ∈ {+∞, −∞, nan}:
z5>1
− − − → z5 ∈ [1.0000000000000002, 1.798 · 10308]
i6
− → z4 ∈ [−1.798 · 10308, 1.798 · 10308]
d5
− → z4 ∈ [1, 1.798 · 10308]
i4
− → z3 ∈ [−710.5, 710.5]
d3
− → z2 ∈ [−1, 1]
d6
− → z5 ∈ ∅ = [−1, 1] ∩ [1.0000000000000002, 1.798 · 10308]
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Implementation and Results
Implementation
The described projection algorithms have been integrated into the ECLAIR platform ECLAIR’s test-data generation and symbolic model checking module has been run, internally, on real-world safety critical code The implementation also contains efficient and precise projection algorithms for all floating-point and integer operations It uses multi-intervals (intervals are used in the papers and in this presentation for simplicity) and it exploits aliasing Work is still needed to turn this functionality in a real product
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Implementation and Results
Analysis of Paparazzi UAV Code
1 #include <math.h> 2 #include <stdint.h> 3 4 // Omitted macros CScal, CAdd, CSub... 10 #define CI(z) { float tmp = z.re; z.re = - z.im; z.im = tmp; } 11 #define CExp(z) { float e = exp(z.re); z.re = e*cos(z.im); \ 12 z.im = e*sin(z.im); } 13 #define CSin(z) { CI(z); struct complex _z = {-z.re, -z.im}; \ 14 float e = exp(z.re); float cos_z_im = cos(z.im); z.re = e*cos_z_im; \ 15 float sin_z_im = sin(z.im); z.im = e*sin_z_im; _z.re = cos_z_im/e; \ 16 _z.im = -sin_z_im/e; CSub(_z, z); CScal(-0.5, z); CI(z); } 17 18 static inline float isometric_latitude(float phi, float e) { 19 return log(tan(M_PI_4 + phi / 2.0)) 20
21 } 22 23 static inline float isometric_latitude0(float phi) { 24 return log(tan(M_PI_4 + phi / 2.0)); 25 }
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Implementation and Results
Analysis of Paparazzi UAV Code (cont’d)
27 void latlong_utm_of(float phi, float lambda, uint8_t utm_zone) { 28 float lambda_c = LambdaOfUtmZone(utm_zone); 29 float ll = isometric_latitude(phi, E); 30 float dl = lambda - lambda_c; 31 float phi_ = asin(sin(dl) / cosh(ll)); 32 float ll_ = isometric_latitude0(phi_); 33 float lambda_ = atan(sinh(ll) / cos(dl)); 34 struct complex z_ = { lambda_, ll_ }; 35 CScal(serie_coeff_proj_mercator[0], z_); 36 uint8_t k; 37 for(k = 1; k < 3; k++) { 38 struct complex z = { lambda_, ll_ }; 39 CScal(2*k, z); 40 CSin(z); 41 CScal(serie_coeff_proj_mercator[k], z); 42 CAdd(z, z_); 43 } 44 CScal(N, z_); 45 latlong_utm_x = XS + z_.im; 46 latlong_utm_y = z_.re; 47 }
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Implementation and Results
Analysis of Paparazzi UAV Code (cont’d)
We use ECLAIR to find finite-to-infinite and numeric-to-NaN transitions, and invocation of trigonometric functions with ill-conditioned arguments Ill-conditioned trig. arguments (test cases): line 19 (−0x86487f.p − 18F, +0.0F, 1), where −0x86487f.p − 18F ≈ −33.570797, line 20 (−0x8a3ae7.p − 19F, +0.0F, 1), where −0x8a3ae7.p − 19F ≈ −17.278761, line 31 (+0.0F, −0x98b6c1.p − 19F, 1), where −0x98b6c1.p − 19F ≈ −19.089235.
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Implementation and Results
Analysis of Paparazzi UAV Code (cont’d)
We use ECLAIR to find finite-to-infinite and numeric-to-NaN transitions, and invocation of trigonometric functions with ill-conditioned arguments numeric-to-NaN (test cases): line 19 (+0xc90fdb.p − 23F, +0.0F, 1), where +0xc90fdb.p − 23F ≈ 1.570796, line 24 (−0x8a9483.p − 53F, +0xcfc398.p − 23F, 1), where −0x8a9483.p − 53F ≈ −1.0083 · 10−9, and +0xcfc398.p − 23F ≈ −17.278761.
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Implementation and Results
Analysis of Paparazzi UAV Code (cont’d)
Taking into account the caller context, we add
1 assert(-M_PI_2 <= phi && phi <= M_PI_2); 2 assert(-M_PI <= lambda && lambda <= M_PI); 3 assert(1 <= utm_zone && utm_zone <= 60);
Now ECLAIR finds that if gps lat = −899 999 991 gps lon = −1 800 000 000 a numeric-to-NaN transition is triggered at line 19.
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Implementation and Results
Analysis of Paparazzi UAV Code (cont’d)
If we add
18 assert( -899999990 <= gps_lat && gps_lat <= 899999990); 19 assert(-1800000000 <= gps_lon && gps_lon <= 1800000000);
ECLAIR finds no more anomalies If the assumptions on the input and on glitch data are correct, none of the 154 potential run-time anomalies can occur E.g., no division by 0 can occur at line 33; 4 integer overflows, 4 inexact conversions, 10 ill-conditioned trig. arguments, 70 finite-to-infinity and 66-numeric-to-NaN transitions can be excluded
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Implementation and Results
Analysis of GNU Scientific Library (GSL) Code
This function computes the scaled regular modified cylindrical Bessel function of first order, exp(−|x|)I1(x)
1 int gsl_sf_bessel_i1_scaled_e(const double x, gsl_sf_result * result) 2 { 3 double ax = fabs(x); 4 5 /* CHECK_POINTER(result) */ 6 7 if(x == 0.0) { 8 result->val = 0.0; 9 result->err = 0.0; 10 return GSL_SUCCESS; 11 } 12 else if(ax < 3.0*GSL_DBL_MIN) { 13 UNDERFLOW_ERROR(result); 14 }
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Implementation and Results
h, s, g = hard, soft, gradual underflow; a = absorption; i = finite-to-infinity; n = numeric-to-NaN.
15 else if(ax < 0.25) { 16 const double eax = exp(-ax); 17 const double y = x*hx; 18 const double c1 = 1.0/10.0; 19 const double c2 = 1.0/280.0; 20 const double c3 = 1.0/15120.0; 21 const double c4 = 1.0/1330560.0; 22 const double c5 = 1.0/172972800.0; 23 const double sum = 24 1.0 +a y*sg(c1 +a y*sg(c2 +a y*sg(c3 +a y*sg(c4 +a y*sgc5)))); 25 result->val = eax * x/3.0 * sum; 26 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); 27 return GSL_SUCCESS; 28 } 29 else { 30 double ex = exphs(-2.0*iax); 31 result->val = 0.5 * (ax*(1.0+aex) -a (1.0-aex)) /n (ax*iax); 32 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); 33 if(x < 0.0) result->val = -result->val; 34 return GSL_SUCCESS; 35 } 36 }
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Implementation and Results
Analysis GSL Code (cont’d)
1 numeric-to-NaN transition 2 finite-to-infinite transitions 2 hard underflows 5 gradual underflows 6 soft underflows line exception input 17 hard underflow x = -0x1.8p-1021 ≈ −6.6752 × 10−308 30 finite-to-infinity x = +/-0x1p+1023 ≈ ±8.9885 × 10307 31 finite-to-infinity x = +/-0x1p+1023 ≈ ±8.9885 × 10307 31 numeric-to-NaN x = −∞
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Conclusion
Conclusion
“Without a specification, a system cannot be right or wrong, it can
This was the current state of affairs for C/C++ software using the math.h/cmath functions We can now properly verify the absence of run-time anomalies for such software: as far as we know, this was completely out of reach Verification in the strong sense is only feasible modulo the possibility of conservatively bounding the size of glitches When this is not possible, our approach can still detect many high-severity program issues via the automatic generation of test input data
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Conclusion
The End
info@bugseng.com roberto.bagnara@bugseng.com
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Conclusion
SMT Solvers
Some of the most popular tools for the verification of numerical programs Source code ⇒ Verification Conditions (VCs) ⇒ SMT Solver Constraints encoded with, e.g., SMT-LIB Decision procedures based on bit-vectors, Abstract CDCL, etc. No direct support for libm functions (other than sqrt)
Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations