More precise abstractions of trigonometric functions using c-stubs

src/cdomain/value/cdomains/floatDomain.ml

@@ -5,6 +5,23 @@ open FloatOps
 
 exception ArithmeticOnFloatBot of string
 
+(** Define records that hold mutable variables representing different Configuration values.
+  * These values are used to keep track of whether or not the corresponding Config values are en-/disabled  *)
+type ana_float_config_values = {
+  mutable evaluate_math_functions : bool option;
+}
+
+let ana_float_config: ana_float_config_values = {
+  evaluate_math_functions = None;
+}
+
+let get_evaluate_math_functions () =
+  if ana_float_config.evaluate_math_functions = None then
+    ana_float_config.evaluate_math_functions <- Some (GobConfig.get_bool "ana.float.evaluate_math_functions");
+  Option.get ana_float_config.evaluate_math_functions
+
+let reset_lazy () = ana_float_config.evaluate_math_functions <- None
+
 module type FloatArith = sig
   type t
 
@@ -273,12 +290,12 @@ module FloatIntervalImpl(Float_t : CFloatType) = struct
     | _ -> Bot
 
   (** [widen x y] assumes [leq x y]. Solvers guarantee this by calling [widen old (join old new)]. *)
-  let widen v1 v2 = (**TODO: support 'threshold_widening' option *)
+  let widen v1 v2 = (* TODO: support 'threshold_widening' option *)
     match v1, v2 with
     | Top, _ | _, Top -> Top
     | Bot, v | v, Bot -> v
     | Interval (l1, h1), Interval (l2, h2) ->
-      (**If we widen and we know that neither interval contains +-inf or nan, it is ok to widen only to +-max_float,
+      (* If we widen and we know that neither interval contains +-inf or nan, it is ok to widen only to +-max_float,
          because a widening with +-inf/nan will always result in the case above -> Top *)
       let low = if l1 <= l2 then l1 else Float_t.lower_bound in
       let high = if h1 >= h2 then h1 else Float_t.upper_bound in
@@ -289,7 +306,7 @@ module FloatIntervalImpl(Float_t : CFloatType) = struct
     | _ -> Top
 
   let narrow v1 v2 =
-    match v1, v2 with (**we cannot distinguish between the lower bound beeing -inf or the upper bound beeing inf. Also there is nan *)
+    match v1, v2 with (* we cannot distinguish between the lower bound beeing -inf or the upper bound beeing inf. Also there is nan *)
     | Bot, _ | _, Bot -> Bot
     | Top, _ -> v2
     | Interval (l1, h1), Interval (l2, h2) ->
@@ -623,15 +640,79 @@ module FloatIntervalImpl(Float_t : CFloatType) = struct
       else
         unknown_IInt ()
 
-  (**it seems strange not to return an explicit 1 for negative numbers, but in c99 signbit is defined as: *)
-  (**<<The signbit macro returns a nonzero value if and only if the sign of its argument value is negative.>> *)
+  (**it seems strange not to return an explicit 1 for negative numbers, but in c99 signbit is defined as:
+   **<<The signbit macro returns a nonzero value if and only if the sign of its argument value is negative.>> *)
   let eval_signbit = function
     | (_, h) when h < Float_t.zero -> true_nonZero_IInt ()
     | (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 *)
-  let overapprox_pi = 3.1416
+    | _ -> unknown_IInt () (* any interval containing zero has to fall in this case, because we do not distinguish between 0. and -0. *)
+
+  (** returns the max of the given mfun once computed with rounding mode Up and once with rounding mode down*)
+  let safe_mathfun_up mfun f = max (mfun Down f) (mfun Up f)
+  (** returns the min of the given mfun once computed with rounding mode Up and once with rounding mode down*)
+  let safe_mathfun_down mfun f = min (mfun Down f) (mfun Up f)
+
+  (** This function does two things: 
+   ** 1. projects l and h onto the interval [0, k*pi] (for k = 2 this is the phase length of sin/cos, for k = 1 it is the phase length of tan)
+   ** 2. compresses/transforms the interval [0, k*pi] to the interval [0, 1] to ease further computations
+   ** i.e. the function computes dist = distance, l'' = (l/(k*pi)) - floor(l/(k*pi)), h'' = (h/(k*pi)) - floor(h/(k*pi))*)
+  let project_and_compress l h k =
+    let ft_over_kpi = (Float_t.mul Up (Float_t.of_float Up k) Float_t.pi) in
+    let ft_under_kpi = (Float_t.mul Down (Float_t.of_float Down k) Float_t.pi) in
+    let l' =
+      if l >= Float_t.zero then (Float_t.div Down l ft_over_kpi)
+      else (Float_t.div Down 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 Up 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 Up (Z.to_float (Float_t.to_big_int l')))
+      else
+        Float_t.sub Down l' (Float_t.of_float Up (Z.to_float (Z.pred (Float_t.to_big_int l'))))
+    in
+    let h'' =
+      if h' >= Float_t.zero then
+        Float_t.sub Up h' (Float_t.of_float Down (Z.to_float (Float_t.to_big_int h')))
+      else
+        Float_t.sub Up h' (Float_t.of_float Down (Z.to_float (Z.pred (Float_t.to_big_int h'))))
+    in
+    (dist, l'', h'')
+
+  let eval_cos_cfun l h =
+    let (dist, l'', h'') = project_and_compress 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'');
+    if (dist <= Float_t.of_float Down 0.5) && (h'' <= Float_t.of_float Down 0.5) && (l'' <= h'') then
+      (* case: monotonic decreasing interval*)
+      Interval (safe_mathfun_down Float_t.cos h, safe_mathfun_up Float_t.cos l)
+    else if (dist <= Float_t.of_float Down 0.5) && (l'' >= Float_t.of_float Up 0.5) && (l'' <= h'') then
+      (* case: monotonic increasing interval*)
+      Interval (safe_mathfun_down Float_t.cos l, safe_mathfun_up Float_t.cos h)
+    else if (dist <= Float_t.of_float Down 1.) && (l'' <= h'') then
+      (* case: contains at most one minimum*)
+      Interval (Float_t.of_float Down (-.1.), max (safe_mathfun_up Float_t.cos l) (safe_mathfun_up Float_t.cos h))
+    else if (dist <= Float_t.of_float Down 1.) && (l'' >= Float_t.of_float Up 0.5) && (h'' <= Float_t.of_float Down 0.5) then
+      (* case: contains at most one maximum*)
+      Interval (min (safe_mathfun_down Float_t.cos l) (safe_mathfun_down Float_t.cos h), Float_t.of_float Up 1.)
+    else
+      of_interval (-. 1., 1.)
+
+  let eval_sin_cfun l h =
+    let lcos = Float_t.sub Down l (Float_t.div Up Float_t.pi (Float_t.of_float Down 2.0)) in
+    let hcos = Float_t.sub Up h (Float_t.div Down Float_t.pi (Float_t.of_float Up 2.0)) in
+    eval_cos_cfun lcos hcos
+
+  let eval_tan_cfun l h =
+    let (dist, l'', h'') = project_and_compress 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 Up 0.5))) then
+      (* case: monotonic increasing interval*)
+      Interval (safe_mathfun_down Float_t.tan l, safe_mathfun_up Float_t.tan h)
+    else
+      top ()
 
   let eval_fabs = function
     | (l, h) when l > Float_t.zero -> Interval (l, h)
@@ -656,39 +737,51 @@ 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.]";
+      Interval (Float_t.of_float Down 0., Float_t.pi)
+    | (l, h) when get_evaluate_math_functions () ->
+      norm @@ Interval (safe_mathfun_down Float_t.acos h, safe_mathfun_up Float_t.acos l) (* acos is monotonic decreasing in [-1, 1]*)
+    | _ -> Interval (Float_t.of_float Down 0., Float_t.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 (Interval (Float_t.neg Float_t.pi, Float_t.pi)) (of_const 2.)
+    | (l, h) when get_evaluate_math_functions () ->
+      norm @@ Interval (safe_mathfun_down Float_t.asin l, safe_mathfun_up Float_t.asin h) (* asin is monotonic increasing in [-1, 1]*)
+    | _ -> div (Interval (Float_t.neg Float_t.pi, Float_t.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 get_evaluate_math_functions () ->
+      norm @@ Interval (safe_mathfun_down Float_t.atan l, safe_mathfun_up Float_t.atan h) (* atan is monotonic increasing*)
+    | _ -> div (Interval (Float_t.neg Float_t.pi, Float_t.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 get_evaluate_math_functions () ->
+      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 get_evaluate_math_functions () ->
+      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 get_evaluate_math_functions () ->
+      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.
-    | (l, h) when l >= Float_t.zero ->
-      let low = Float_t.sqrt Down l in
-      let high = Float_t.sqrt Up h in
+    | (l, h) when l >= Float_t.zero && get_evaluate_math_functions () ->
+      let low = safe_mathfun_down Float_t.sqrt l in
+      let high = safe_mathfun_up Float_t.sqrt h in
       Interval (low, high)
     | _ -> top ()
 

src/cdomain/value/cdomains/floatDomain.mli

@@ -4,6 +4,8 @@ open GoblintCil
 
 exception ArithmeticOnFloatBot of string
 
+val reset_lazy: unit -> unit
+
 module type FloatArith = sig
   type t
 

src/common/cdomains/floatOps/floatOps.ml

@@ -13,6 +13,7 @@ module type CFloatType = sig
   val upper_bound: t
   val lower_bound: t
   val smallest : t
+  val pi : t
 
   val of_float: round_mode -> float -> t
   val to_float: t -> float option
@@ -37,6 +38,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
 
@@ -56,6 +65,7 @@ module CDouble = struct
   let upper_bound = Float.max_float
   let lower_bound = -. Float.max_float
   let smallest = Float.min_float
+  let pi = Float.pi
 
   let of_float _ x = x
   let to_float x = Some x
@@ -78,6 +88,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
 
@@ -89,10 +106,12 @@ module CFloat = struct
 
   external upper': unit -> float = "max_float"
   external smallest': unit -> float = "smallest_float"
+  external pi': unit -> float = "pi_float"
 
   let upper_bound = upper' ()
   let lower_bound = -. upper_bound
   let smallest = smallest' ()
+  let pi = pi' ()
 
   let to_float x = Some x
   let to_big_int = big_int_of_float
@@ -112,6 +131,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

src/common/cdomains/floatOps/floatOps.mli

@@ -14,6 +14,7 @@ module type CFloatType = sig
   val upper_bound: t
   val lower_bound: t
   val smallest : t
+  val pi : t
 
   val of_float: round_mode -> float -> t
   val to_float: t -> float option
@@ -39,6 +40,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
 

src/common/cdomains/floatOps/stubs.c

@@ -1,3 +1,4 @@
+#define _GNU_SOURCE // necessary for M_PI to be defined
 #include <stdio.h>
 #include <math.h>
 #include <float.h>
@@ -54,6 +55,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)   \
@@ -109,6 +122,8 @@ CAMLprim value atof_float(value mode, value str)
     // No need to use CAMLreturn because we don't use CAMLparam.
 }
 
+// These are only given for floats as these operations involve no rounding and their OCaml implementation (Float module) can be used
+
 CAMLprim value max_float(value unit)
 {
     // No need to use CAMLparam to keep unit as GC root,
@@ -124,3 +139,8 @@ CAMLprim value smallest_float(value unit)
     return caml_copy_double(FLT_MIN);
     // No need to use CAMLreturn because we don't use CAMLparam.
 }
+
+CAMLprim value pi_float(value unit)
+{
+    return caml_copy_double(M_PI);
+}

src/config/options.schema.json

@@ -463,6 +463,13 @@
                 "Use FloatDomain: (float * float) option.",
               "type": "boolean",
               "default": false
+            },
+            "evaluate_math_functions": {
+              "title": "ana.float.evaluate_math_functions",
+              "description":
+                "Allow a more precise evaluation of some functions from math.h. Evaluation of functions may differ depending on the implementation of math.h. Caution: For some implementations of functions it is not guaranteed that they behave monotonic where they mathematically should, thus possibly leading to unsoundness.",
+              "type": "boolean",
+              "default": false
             }
           },
           "additionalProperties": false

src/util/server.ml

@@ -284,6 +284,7 @@ let analyze ?(reset=false) (s: t) =
   InvariantCil.reset_lazy ();
   WideningThresholds.reset_lazy ();
   IntDomain.reset_lazy ();
+  FloatDomain.reset_lazy ();
   StringDomain.reset_lazy ();
   PrecisionUtil.reset_lazy ();
   ApronDomain.reset_lazy ();

tests/regression/39-signed-overflows/07-abs-sqrt.c

@@ -1,4 +1,4 @@
-//PARAM: --enable ana.int.interval --enable ana.float.interval --set ana.activated[+] tmpSpecial
+//PARAM: --enable ana.int.interval --enable ana.float.interval --enable ana.float.evaluate_math_functions --set ana.activated[+] tmpSpecial
 #include<math.h>
 int main() {
     int data;

tests/regression/39-signed-overflows/09-labs-sqrt.c

@@ -1,4 +1,4 @@
-//PARAM: --enable ana.int.interval --enable ana.float.interval --set ana.activated[+] tmpSpecial
+//PARAM: --enable ana.int.interval --enable ana.float.interval --enable ana.float.evaluate_math_functions --set ana.activated[+] tmpSpecial
 #include<math.h>
 int main() {
     int data;