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) http://bugseng.com IMDEA, Madrid, Spain — November 22 nd , 2019 Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Introduction Floating-Point Numbers in Critical Systems CONFIDENTIAL 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: Ensure absence of undefined or undesirable behaviors: 1 infinities, NaNs, overflows in conversion to integers. . . Identifying sources of round-off errors 2 Prove that the program implements a certain numerical 3 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 CONFIDENTIAL 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 CONFIDENTIAL 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 errors. The table lists the maximum error for each function which is exposed by one of the existing tests in the test suite. The table tries to cover as much as possible and list 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 CONFIDENTIAL 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 - e / 2.0 * log((1.0 + e * sin(phi)) / (1.0 - e * sin(phi))); 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) CONFIDENTIAL 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 CONFIDENTIAL What do the mathematical functions used in the code compute? More specifically: Can infinities and NaNs be generated? 1 Can sin , cos , and tan be invoked on ill-conditioned 2 arguments? If anomalies of any one of these kinds are possible, which 3 inputs may cause them? Ill-conditioned arguments: if x is an IEEE 754 single- precision number and x ≥ 2 23 , 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? CONFIDENTIAL 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 CONFIDENTIAL 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 : [ x l , x u ] ∈ I F y : [ y l , y u ] ∈ I F Then we iteratively refine the intervals by means of projections: Direct Projection Inverse Projection y ← x y → x Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Solving Constraints on libm Functions Projection Functions CONFIDENTIAL Direct Projection Determine [ y ′ l , y ′ u ] ⊂ [ y l , y u ] such that ∀ x ∈ [ x l , x u ] : f ( x ) ∈ [ y l , y u ] = ⇒ f ( x ) ∈ [ y ′ l , y ′ u ] Inverse Projection Determine [ x ′ l , x ′ u ] ⊂ [ x l , x u ] such that ⇒ x ∈ [ x ′ l , x ′ ∀ x ∈ [ x l , x u ] : f ( x ) ∈ [ y l , y u ] = u ] Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Solving Constraints on libm Functions Exact Projections CONFIDENTIAL 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 CONFIDENTIAL 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 CONFIDENTIAL G 1 G 2 width ( G 1 ) = 5, depth ( G 1 ) = 4, width ( G 2 ) = depth ( G 2 ) = 2. Roberto Bagnara et al., BUGSENG & UNIPR, Italy Precise approximation of Floating-Point Computations
Solving Constraints on libm Functions Glitch Data CONFIDENTIAL Our approach only requires safe approximations of: n g ∈ N : the total number of glitches d M ∈ N : the maximal depth of glitches (in ULPs) w M ∈ 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
Recommend
More recommend
