From cfde041aad06fdbf61b16d880637c2e94ca2b078 Mon Sep 17 00:00:00 2001 From: FelixKrayer Date: Mon, 25 Jul 2022 17:40:30 +0200 Subject: [PATCH] Add math_funeval option and better abstractions using c-stubs --- src/cdomains/floatDomain.ml | 104 ++++++++++++++++--- src/cdomains/floatOps/floatOps.ml | 22 ++++ src/cdomains/floatOps/floatOps.mli | 6 ++ src/cdomains/floatOps/stubs.c | 12 +++ src/common/util/options.schema.json | 7 ++ tests/regression/57-floats/22-math-funeval.c | 67 ++++++++++++ 6 files changed, 205 insertions(+), 13 deletions(-) create mode 100644 tests/regression/57-floats/22-math-funeval.c diff --git a/src/cdomains/floatDomain.ml b/src/cdomains/floatDomain.ml index 39d3744401a..0a61a74c9a0 100644 --- a/src/cdomains/floatDomain.ml +++ b/src/cdomains/floatDomain.ml @@ -618,8 +618,74 @@ module FloatIntervalImpl(Float_t : CFloatType) = struct | (l, _) when l > Float_t.zero -> false_zero_IInt () | _ -> unknown_IInt () (**any interval containing zero has to fall in this case, because we do not distinguish between 0. and -0. *) - (**This Constant overapproximates pi to use as bounds for the return values of trigonometric functions *) + (**These constants over-/underpproximate pi *) let overapprox_pi = 3.1416 + let underapprox_pi = 3.1415 + + let trig_helper l h k = (** computes dist = distance, l'' = (l/(k*pi) % 10) - floor(l/(k*pi)), h'' = (h/(k*pi)) - floor(h/(k*pi))*) + let ft_over_kpi = (Float_t.mul Up (Float_t.of_float Up k) (Float_t.of_float Up overapprox_pi)) in + let ft_under_kpi = (Float_t.mul Down (Float_t.of_float Down k) (Float_t.of_float Down underapprox_pi)) in + let l' = + if l >= Float_t.zero then (Float_t.div Down l ft_over_kpi) + else (Float_t.div Up l ft_under_kpi) + in + let h' = + if h >= Float_t.zero then (Float_t.div Up h ft_under_kpi) + else (Float_t.div Down h ft_over_kpi) + in + let dist = (Float_t.sub Up h' l') in + let l'' = + if l' >= Float_t.zero then + Float_t.sub Down l' (Float_t.of_float Down (Big_int_Z.float_of_big_int (Float_t.to_big_int l'))) + else + Float_t.sub Down l' (Float_t.of_float Down (Big_int_Z.float_of_big_int (IntOps.BigIntOps.sub (Float_t.to_big_int l') (Big_int_Z.big_int_of_int 1)))) + in + let h'' = + if h' >= Float_t.zero then + Float_t.sub Up h' (Float_t.of_float Up (Big_int_Z.float_of_big_int (Float_t.to_big_int h'))) + else + Float_t.sub Up h' (Float_t.of_float Up (Big_int_Z.float_of_big_int (IntOps.BigIntOps.sub (Float_t.to_big_int h') (Big_int_Z.big_int_of_int 1)))) + in + (dist, l'', h'') + + let eval_cos_cfun l h = + let (dist, l'', h'') = trig_helper l h 2. in + if Messages.tracing then Messages.trace "CstubsTrig" "sin: dist %s; l'' %s; h'' %s\n" (Float_t.to_string dist) (Float_t.to_string l'') (Float_t.to_string h''); + if dist <= (Float_t.of_float Down 0.5) && (h'' <= Float_t.of_float Up 0.5) && (h'' >= l'') then + Interval (Float_t.cos Down h, Float_t.cos Up l) + else if dist <= (Float_t.of_float Down 0.5) && (l'' >= Float_t.of_float Down 0.5) && (h'' >= l'') then + Interval (Float_t.cos Down l, Float_t.cos Up h) + else if dist <= (Float_t.of_float Down 1.) && (l'' <= h'') then + Interval (Float_t.of_float Down (-.1.), max (Float_t.cos Up l) (Float_t.cos Up h)) + else if dist <= (Float_t.of_float Down 1.) && (l'' >= Float_t.of_float Down 0.5) && (h'' <= Float_t.of_float Down 0.5) then + Interval (min (Float_t.cos Down l) (Float_t.cos Down h), Float_t.of_float Up 1.) + else + of_interval (-. 1., 1.) + + let eval_sin_cfun l h = + let (dist, l'', h'') = trig_helper l h 2. in + if Messages.tracing then Messages.trace "CstubsTrig" "cos: dist %s; l'' %s; h'' %s\n" (Float_t.to_string dist) (Float_t.to_string l'') (Float_t.to_string h''); + (* phase shift -0.25 -> same phase as cos *) + let l'' = if l'' <= Float_t.of_float Down 0.25 then Float_t.add Up l'' (Float_t.of_float Up 0.75) else Float_t.sub Down l'' (Float_t.of_float Up 0.25) in + let h'' = if h'' <= Float_t.of_float Down 0.25 then Float_t.add Up h'' (Float_t.of_float Up 0.75) else Float_t.sub Down h'' (Float_t.of_float Up 0.25) in + if dist <= (Float_t.of_float Down 0.5) && (h'' <= Float_t.of_float Up 0.5) && (h'' >= l'') then + Interval (Float_t.sin Down h, Float_t.sin Up l) + else if dist <= (Float_t.of_float Down 0.5) && (l'' >= Float_t.of_float Down 0.5) && (h'' >= l'') then + Interval (Float_t.sin Down l, Float_t.sin Up h) + else if dist <= (Float_t.of_float Down 1.) && (l'' <= h'') then + Interval (Float_t.of_float Down (-.1.), max (Float_t.sin Up l) (Float_t.sin Up h)) + else if dist <= (Float_t.of_float Down 1.) && (l'' >= Float_t.of_float Down 0.5) && (h'' <= Float_t.of_float Down 0.5) then + Interval (min (Float_t.sin Down l) (Float_t.sin Down h), Float_t.of_float Up 1.) + else + of_interval (-. 1., 1.) + + let eval_tan_cfun l h = + let (dist, l'', h'') = trig_helper l h 1. in + if Messages.tracing then Messages.trace "CstubsTrig" "tan: dist %s; l'' %s; h'' %s\n" (Float_t.to_string dist) (Float_t.to_string l'') (Float_t.to_string h''); + if dist <= (Float_t.of_float Down 1.) && Bool.not ((l'' <= Float_t.of_float Up 0.5) && (h'' >= Float_t.of_float Down 0.5)) then + Interval (Float_t.tan Down l, Float_t.tan Up h) + else + top () let eval_fabs = function | (l, h) when l > Float_t.zero -> Interval (l, h) @@ -644,33 +710,45 @@ module FloatIntervalImpl(Float_t : CFloatType) = struct let eval_acos = function | (l, h) when l = h && l = Float_t.of_float Nearest 1. -> of_const 0. (*acos(1) = 0*) - | (l, h) -> - if l < (Float_t.of_float Down (-.1.)) || h > (Float_t.of_float Up 1.) then - Messages.warn ~category:Messages.Category.Float "Domain error might occur: acos argument might be outside of [-1., 1.]"; - of_interval (0., (overapprox_pi)) (**could be more exact *) + | (l, h) when l < (Float_t.of_float Down (-.1.)) || h > (Float_t.of_float Up 1.) -> + Messages.warn ~category:Messages.Category.Float "Domain error might occur: acos argument might be outside of [-1., 1.]"; + of_interval (0., (overapprox_pi)) + | (l, h) when GobConfig.get_bool "ana.float.math_funeval" -> + norm @@ Interval (Float_t.acos Down h, Float_t.acos Up l) (** acos is monotonic decreasing in [-1, 1]*) + | _ -> of_interval (0., (overapprox_pi)) let eval_asin = function | (l, h) when l = h && l = Float_t.zero -> of_const 0. (*asin(0) = 0*) - | (l, h) -> - if l < (Float_t.of_float Down (-.1.)) || h > (Float_t.of_float Up 1.) then - Messages.warn ~category:Messages.Category.Float "Domain error might occur: asin argument might be outside of [-1., 1.]"; - div (of_interval ((-. overapprox_pi), overapprox_pi)) (of_const 2.) (**could be more exact *) + | (l, h) when l < (Float_t.of_float Down (-.1.)) || h > (Float_t.of_float Up 1.) -> + Messages.warn ~category:Messages.Category.Float "Domain error might occur: asin argument might be outside of [-1., 1.]"; + div (of_interval ((-. overapprox_pi), overapprox_pi)) (of_const 2.) + | (l, h) when GobConfig.get_bool "ana.float.math_funeval" -> + norm @@ Interval (Float_t.asin Down l, Float_t.asin Up h) (** asin is monotonic increasing in [-1, 1]*) + | _ -> div (of_interval ((-. overapprox_pi), overapprox_pi)) (of_const 2.) let eval_atan = function | (l, h) when l = h && l = Float_t.zero -> of_const 0. (*atan(0) = 0*) - | _ -> div (of_interval ((-. overapprox_pi), overapprox_pi)) (of_const 2.) (**could be more exact *) + | (l, h) when GobConfig.get_bool "ana.float.math_funeval" -> + norm @@ Interval (Float_t.atan Down l, Float_t.atan Up h) (** atan is monotonic increasing*) + | _ -> div (of_interval ((-. overapprox_pi), overapprox_pi)) (of_const 2.) let eval_cos = function | (l, h) when l = h && l = Float_t.zero -> of_const 1. (*cos(0) = 1*) - | _ -> of_interval (-. 1., 1.) (**could be exact for intervals where l=h, or even for Interval intervals *) + | (l, h) when GobConfig.get_bool "ana.float.math_funeval" -> + norm @@ eval_cos_cfun l h + | _ -> of_interval (-. 1., 1.) let eval_sin = function | (l, h) when l = h && l = Float_t.zero -> of_const 0. (*sin(0) = 0*) - | _ -> of_interval (-. 1., 1.) (**could be exact for intervals where l=h, or even for some intervals *) + | (l, h) when GobConfig.get_bool "ana.float.math_funeval" -> + norm @@ eval_sin_cfun l h + | _ -> of_interval (-. 1., 1.) let eval_tan = function | (l, h) when l = h && l = Float_t.zero -> of_const 0. (*tan(0) = 0*) - | _ -> top () (**could be exact for intervals where l=h, or even for some intervals *) + | (l, h) when GobConfig.get_bool "ana.float.math_funeval" -> + norm @@ eval_tan_cfun l h + | _ -> top () let eval_sqrt = function | (l, h) when l = Float_t.zero && h = Float_t.zero -> of_const 0. diff --git a/src/cdomains/floatOps/floatOps.ml b/src/cdomains/floatOps/floatOps.ml index a4e39d930ef..b73e7b424e4 100644 --- a/src/cdomains/floatOps/floatOps.ml +++ b/src/cdomains/floatOps/floatOps.ml @@ -36,6 +36,14 @@ module type CFloatType = sig val mul: round_mode -> t -> t -> t val div: round_mode -> t -> t -> t val sqrt: round_mode -> t -> t + + val acos: round_mode -> t -> t + val asin: round_mode -> t -> t + val atan: round_mode -> t -> t + val cos: round_mode -> t -> t + val sin: round_mode -> t -> t + val tan: round_mode -> t -> t + val atof: round_mode -> string -> t end @@ -77,6 +85,13 @@ module CDouble = struct external div: round_mode -> t -> t -> t = "div_double" external sqrt: round_mode -> t -> t = "sqrt_double" + external acos: round_mode -> t -> t = "acos_double" + external asin: round_mode -> t -> t = "asin_double" + external atan: round_mode -> t -> t = "atan_double" + external cos: round_mode -> t -> t = "cos_double" + external sin: round_mode -> t -> t = "sin_double" + external tan: round_mode -> t -> t = "tan_double" + external atof: round_mode -> string -> t = "atof_double" end @@ -111,6 +126,13 @@ module CFloat = struct external div: round_mode -> t -> t -> t = "div_float" external sqrt: round_mode -> t -> t = "sqrt_float" + external acos: round_mode -> t -> t = "acos_float" + external asin: round_mode -> t -> t = "asin_float" + external atan: round_mode -> t -> t = "atan_float" + external cos: round_mode -> t -> t = "cos_float" + external sin: round_mode -> t -> t = "sin_float" + external tan: round_mode -> t -> t = "tan_float" + external atof: round_mode -> string -> t = "atof_float" let of_float mode x = add mode zero x diff --git a/src/cdomains/floatOps/floatOps.mli b/src/cdomains/floatOps/floatOps.mli index cf24f75ed5d..75d04818224 100644 --- a/src/cdomains/floatOps/floatOps.mli +++ b/src/cdomains/floatOps/floatOps.mli @@ -39,6 +39,12 @@ module type CFloatType = sig val mul: round_mode -> t -> t -> t val div: round_mode -> t -> t -> t val sqrt: round_mode -> t -> t + val acos: round_mode -> t -> t + val asin: round_mode -> t -> t + val atan: round_mode -> t -> t + val cos: round_mode -> t -> t + val sin: round_mode -> t -> t + val tan: round_mode -> t -> t val atof: round_mode -> string -> t end diff --git a/src/cdomains/floatOps/stubs.c b/src/cdomains/floatOps/stubs.c index 50e4a2fb313..4489ef71e16 100644 --- a/src/cdomains/floatOps/stubs.c +++ b/src/cdomains/floatOps/stubs.c @@ -49,6 +49,18 @@ static void change_round_mode(int mode) UNARY_OP(sqrt, double, sqrt); UNARY_OP(sqrt, float, sqrtf); +UNARY_OP(acos, double, acos); +UNARY_OP(acos, float, acosf); +UNARY_OP(asin, double, asin); +UNARY_OP(asin, float, asinf); +UNARY_OP(atan, double, atan); +UNARY_OP(atan, float, atanf); +UNARY_OP(cos, double, cos); +UNARY_OP(cos, float, cosf); +UNARY_OP(sin, double, sin); +UNARY_OP(sin, float, sinf); +UNARY_OP(tan, double, tan); +UNARY_OP(tan, float, tanf); #define BINARY_OP(name, type, op) \ CAMLprim value name##_##type(value mode, value x, value y) \ diff --git a/src/common/util/options.schema.json b/src/common/util/options.schema.json index 7c921c4f53d..a3145d54690 100644 --- a/src/common/util/options.schema.json +++ b/src/common/util/options.schema.json @@ -463,6 +463,13 @@ "Use FloatDomain: (float * float) option.", "type": "boolean", "default": false + }, + "math_funeval": { + "title": "ana.float.math_funeval", + "description": + "Allow evaluation of functions from math.h with c-stubs. Only sound, if the same math.h implementation is used for the programm that is also used in Goblint", + "type": "boolean", + "default": false } }, "additionalProperties": false diff --git a/tests/regression/57-floats/22-math-funeval.c b/tests/regression/57-floats/22-math-funeval.c new file mode 100644 index 00000000000..51a9454fa8a --- /dev/null +++ b/tests/regression/57-floats/22-math-funeval.c @@ -0,0 +1,67 @@ +// PARAM: --enable ana.float.interval --enable ana.float.math_funeval +#include +#include +#include + +int main() +{ + int x; + double s1, s2, s3, c1, c2, sc1, sc2, t1, t2; + //s1: 0.5pi < [2.0, 7.5] < 2.5pi + //s2: 1.5pi < [5.0, 10.5] < 3.5pi + //s3: -0.5pi < [-1.0, -1.0] < 0.5pi + //c1: 0pi < [0.2, 6.1] < 2pi + //c2: pi < [3.3, 9.0] < 3pi + //sc1: 0.5pi < [2.0, 3.0] < pi + //sc2: 4.5pi < [14.5, 15.5] < 5pi + //t1: 0pi-0.1 <= [-0.1, -0.1] <= 0pi+0.1 + //t2: 6pi-0.1 < [18.8, 18.9] < 6pi+0.1 + if (x) { + s1 = 2.0; + s2 = 5.0; + s3 = -1.0; + c1 = 0.2; + c2 = 3.3; + sc1 = 2.0; + sc2 = 14.5; + t1 = -0.1; + t2 = 18.8; + } + else { + s1 = 7.5; + s2 = 10.5; + s3 = 1.0; + c1 = 6.1; + c2 = 9.0; + sc1 = 3.0; + sc2 = 15.5; + t1 = 0.1; + t2 = 18.9; + } + + //acos(x) + assert(1.4 < acos(0.1) && acos(0.1) < 1.5); // SUCCESS + + //asin(x) + assert(0.6 < asin(0.6) && asin(0.6) < 0.7); // SUCCESS + + //atan(x) + assert(0.7 < atan(1.) && atan(1.) < 0.8); // SUCCESS + + //cos(x) + assert(-1. <= cos(c1) && cos(c1) < 0.99); // SUCCESS + assert(-0.99 < cos(c2) && cos(c2) <= 1.0); // SUCCESS + assert(-0.99 < cos(sc1) && cos(sc1) < 0.); // SUCCESS + assert(-0.99 < cos(sc2) && cos(sc2) < 0.); // SUCCESS + + //sin(x) + assert(-1. <= sin(s1) && sin(s1) < 0.99); // SUCCESS + assert(-0.99 < sin(s2) && sin(s2) <= 1.); // SUCCESS + assert(-0.99 < sin(s3) && sin(s3) <= 0.99); // SUCCESS + assert(0. < sin(sc1) && sin(sc1) < 0.99); // SUCCESS + assert(0. < sin(sc2) && sin(sc2) < 0.99); // SUCCESS + + //tan(x) + assert(-0.11 < tan(t1) && tan(t1) < 0.11 ); // SUCCESS + assert(-0.1 < tan(t2) && tan(t2) < 0.1 ); // SUCCESS +}