More precise abstractions of trigonometric functions using c-stubs

src/cdomain/value/cdomains/floatDomain.ml

@@ -618,8 +618,87 @@
     | (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
+
+  (** 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 it 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.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 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 (IntOps.BigIntOps.sub (Float_t.to_big_int l') (Z.of_int 1))))
+    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 (IntOps.BigIntOps.sub (Float_t.to_big_int h') (Z.of_int 1))))
+    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" "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 Down 0.5) && (h'' >= l'') then
+      (** case: monotonic decreasing interval*)
+      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 Up 0.5) && (h'' >= l'') then
+      (** case: monotonic increasing interval*)
+      Interval (Float_t.cos Down l, Float_t.cos Up 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 (Float_t.cos Up l) (Float_t.cos Up 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 (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'') = 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'');
+    (* phase shift -0.25 -> same phase as cos *)
+    let l'' = if l'' <= Float_t.of_float Down 0.25 then Float_t.add Down l'' (Float_t.of_float Down 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 Up h'' (Float_t.of_float Down 0.25) in
+    if (dist <= Float_t.of_float Down 0.5) && (h'' <= Float_t.of_float Down 0.5) && (h'' >= l'') then
+      (** case: monotonic decreasing interval*)
+      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 Up 0.5) && (h'' >= l'') then
+      (** case: monotonic increasing interval*)
+      Interval (Float_t.sin Down l, Float_t.sin Up 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 (Float_t.sin Up l) (Float_t.sin Up 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 (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'') = 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 (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 +723,33 @@
 
   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) -> norm @@ Interval (Float_t.acos Down h, Float_t.acos Up l) (** acos is monotonic decreasing in [-1, 1]*)
 
   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) -> norm @@ Interval (Float_t.asin Down l, Float_t.asin Up h) (** asin is monotonic increasing in [-1, 1]*)
 
   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) -> norm @@ Interval (Float_t.atan Down l, Float_t.atan Up h) (** atan is monotonic increasing*)
 
   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) -> norm @@ eval_cos_cfun l h
 
   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) -> norm @@ eval_sin_cfun l h
 
   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) -> norm @@ eval_tan_cfun l h
 
   let eval_sqrt = function
     | (l, h) when l = Float_t.zero && h = Float_t.zero -> of_const 0.

src/common/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

src/common/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
 

src/common/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)   \

tests/regression/57-floats/22-math-funeval.c

@@ -0,0 +1,67 @@
+// PARAM: --enable ana.float.interval
+#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
+}