Skip to content

Commit

Permalink
Add math_funeval option and better abstractions using c-stubs
Browse files Browse the repository at this point in the history
  • Loading branch information
FelixKrayer committed Nov 29, 2023
1 parent 7fa7bfd commit cfde041
Show file tree
Hide file tree
Showing 6 changed files with 205 additions and 13 deletions.
104 changes: 91 additions & 13 deletions src/cdomains/floatDomain.ml
Original file line number Diff line number Diff line change
Expand Up @@ -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')))

Check warning

Code scanning / Semgrep

Semgrep Finding: semgrep.big_int_z-to_float Warning

use Z instead

Check warning

Code scanning / Semgrep

Semgrep Finding: semgrep.big_int_z Warning

use Z instead
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))))

Check warning

Code scanning / Semgrep

Semgrep Finding: semgrep.big_int_z-to_float Warning

use Z instead

Check warning

Code scanning / Semgrep

Semgrep Finding: semgrep.big_int_z-of_int Warning

use Z instead

Check warning

Code scanning / Semgrep

Semgrep Finding: semgrep.big_int_z Warning

use Z instead

Check warning

Code scanning / Semgrep

Semgrep Finding: semgrep.big_int_z Warning

use Z instead
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')))

Check warning

Code scanning / Semgrep

Semgrep Finding: semgrep.big_int_z-to_float Warning

use Z instead

Check warning

Code scanning / Semgrep

Semgrep Finding: semgrep.big_int_z Warning

use Z instead
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))))

Check warning

Code scanning / Semgrep

Semgrep Finding: semgrep.big_int_z-to_float Warning

use Z instead

Check warning

Code scanning / Semgrep

Semgrep Finding: semgrep.big_int_z-of_int Warning

use Z instead

Check warning

Code scanning / Semgrep

Semgrep Finding: semgrep.big_int_z Warning

use Z instead

Check warning

Code scanning / Semgrep

Semgrep Finding: semgrep.big_int_z Warning

use Z instead
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)
Expand All @@ -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.
Expand Down
22 changes: 22 additions & 0 deletions src/cdomains/floatOps/floatOps.ml
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down
6 changes: 6 additions & 0 deletions src/cdomains/floatOps/floatOps.mli
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
12 changes: 12 additions & 0 deletions src/cdomains/floatOps/stubs.c
Original file line number Diff line number Diff line change
Expand Up @@ -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) \
Expand Down
7 changes: 7 additions & 0 deletions src/common/util/options.schema.json
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
67 changes: 67 additions & 0 deletions tests/regression/57-floats/22-math-funeval.c
Original file line number Diff line number Diff line change
@@ -0,0 +1,67 @@
// PARAM: --enable ana.float.interval --enable ana.float.math_funeval
#include <assert.h>
#include <math.h>
#include <float.h>

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
}

0 comments on commit cfde041

Please sign in to comment.