add sqrt handling in float domain

src/analyses/base.ml

@@ -2356,6 +2356,7 @@ struct
           | Isunordered (x,y) -> Int(ID.cast_to IInt (apply_binary FDouble FD.unordered x y))
           | Fmax (fd, x ,y) -> Float (apply_binary fd FD.fmax x y)
           | Fmin (fd, x ,y) -> Float (apply_binary fd FD.fmin x y)
+          | Sqrt (fk, x) -> Float (apply_unary fk FD.sqrt x)
         end
       in
       begin match lv with

src/analyses/libraryDesc.ml

@@ -38,7 +38,8 @@ type math =
   | Atan2 of (CilType.Fkind.t * Basetype.CilExp.t * Basetype.CilExp.t)
   | Cos of (CilType.Fkind.t * Basetype.CilExp.t)
   | Sin of (CilType.Fkind.t * Basetype.CilExp.t)
-  | Tan of (CilType.Fkind.t * Basetype.CilExp.t) [@@deriving eq, ord, hash]
+  | Tan of (CilType.Fkind.t * Basetype.CilExp.t)
+  | Sqrt of (CilType.Fkind.t * Basetype.CilExp.t) [@@deriving eq, ord, hash]
 
 (** Type of special function, or {!Unknown}. *)
 (* Use inline record if not single {!Cil.exp} argument. *)
@@ -170,6 +171,7 @@ module MathPrintable = struct
     | Cos (fk, exp) -> Pretty.dprintf "(%a )cos(%a)" d_fkind fk d_exp exp
     | Sin (fk, exp) -> Pretty.dprintf "(%a )sin(%a)" d_fkind fk d_exp exp
     | Tan (fk, exp) -> Pretty.dprintf "(%a )tan(%a)" d_fkind fk d_exp exp
+    | Sqrt (fk, exp) -> Pretty.dprintf "(%a )sqrt(%a)" d_fkind fk d_exp exp
 
   include Printable.SimplePretty (
     struct

src/analyses/libraryFunctions.ml

@@ -944,9 +944,9 @@ let math_descs_list: (string * LibraryDesc.t) list = LibraryDsl.[
     ("scalbln", unknown [drop "arg" []; drop "exp" []]);
     ("scalblnf", unknown [drop "arg" []; drop "exp" []]);
     ("scalblnl", unknown [drop "arg" []; drop "exp" []]);
-    ("sqrt", unknown [drop "x" []]);
-    ("sqrtf", unknown [drop "x" []]);
-    ("sqrtl", unknown [drop "x" []]);
+    ("sqrt", special [__ "x" []] @@ fun x -> Math { fun_args = (Sqrt (FDouble, x)) });
+    ("sqrtf", special [__ "x" []] @@ fun x -> Math { fun_args = (Sqrt (FFloat, x)) });
+    ("sqrtl", special [__ "x" []] @@ fun x -> Math { fun_args = (Sqrt (FLongDouble, x)) });
     ("tgamma", unknown [drop "x" []]);
     ("tgammaf", unknown [drop "x" []]);
     ("tgammal", unknown [drop "x" []]);

src/cdomains/floatDomain.ml

@@ -40,6 +40,8 @@ module type FloatArith = sig
   (** sin(x) *)
   val tan : t -> t
   (** tan(x) *)
+  val sqrt : t -> t
+  (** sqrt(x) *)
 
   (** {inversions of unary functions}*)
   val inv_ceil : ?asPreciseAsConcrete:bool -> t -> t
@@ -670,6 +672,14 @@ module FloatIntervalImpl(Float_t : CFloatType) = struct
     | (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 *)
 
+  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
+      Interval (low, high)
+    | _ -> top ()
+
   let eval_inv_ceil ?(asPreciseAsConcrete=false) = function
     | (l, h) ->
       if (Float_t.sub Up (Float_t.ceil l) (Float_t.sub Down (Float_t.ceil l) (Float_t.of_float Nearest 1.0)) = (Float_t.of_float Nearest 1.0)) then (
@@ -784,6 +794,7 @@ module FloatIntervalImpl(Float_t : CFloatType) = struct
   let cos = eval_unop eval_cos
   let sin = eval_unop eval_sin
   let tan = eval_unop eval_tan
+  let sqrt = eval_unop eval_sqrt
 
   let inv_ceil ?(asPreciseAsConcrete=false) = eval_unop ~warn:false (eval_inv_ceil ~asPreciseAsConcrete:asPreciseAsConcrete)
   let inv_floor ?(asPreciseAsConcrete=false) = eval_unop ~warn:false (eval_inv_floor ~asPreciseAsConcrete:asPreciseAsConcrete)
@@ -899,6 +910,7 @@ module FloatIntervalImplLifted = struct
   let cos = lift (F1.cos, F2.cos)
   let sin = lift (F1.sin, F2.sin)
   let tan = lift (F1.tan, F2.tan)
+  let sqrt = lift (F1.sqrt, F2.sqrt)
 
   let inv_ceil ?(asPreciseAsConcrete=BoolDomain.MustBool.top ()) = function
     | F32 a -> F32 (F1.inv_ceil ~asPreciseAsConcrete:true a)
@@ -1159,6 +1171,8 @@ module FloatDomTupleImpl = struct
     map { f1= (fun (type a) (module F : FloatDomain with type t = a) -> F.sin); }
   let tan =
     map { f1= (fun (type a) (module F : FloatDomain with type t = a) -> F.tan); }
+  let sqrt =
+    map { f1= (fun (type a) (module F : FloatDomain with type t = a) -> F.sqrt); }
 
   (*"asPreciseAsConcrete" has no meaning here*)
   let inv_ceil ?(asPreciseAsConcrete=BoolDomain.MustBool.top ()) =

src/cdomains/floatDomain.mli

@@ -57,6 +57,9 @@ module type FloatArith = sig
   val tan : t -> t
   (** tan(x) *)
 
+  val sqrt : t -> t
+  (** sqrt(x) *)
+
   (** {inversions of unary functions}*)
   val inv_ceil : ?asPreciseAsConcrete:bool -> t -> t
   (** (inv_ceil z -> x) if (z = ceil(x)) *)

src/cdomains/floatOps/floatOps.ml

@@ -35,6 +35,7 @@ module type CFloatType = sig
   val sub: round_mode -> t -> t -> t
   val mul: round_mode -> t -> t -> t
   val div: round_mode -> t -> t -> t
+  val sqrt: round_mode -> t -> t
   val atof: round_mode -> string -> t
 end
 
@@ -74,6 +75,7 @@ module CDouble = struct
   external sub: round_mode -> t -> t -> t = "sub_double"
   external mul: round_mode -> t -> t -> t = "mul_double"
   external div: round_mode -> t -> t -> t = "div_double"
+  external sqrt: round_mode -> t -> t = "sqrt_double"
 
   external atof: round_mode -> string -> t = "atof_double"
 end
@@ -107,6 +109,7 @@ module CFloat = struct
   external sub: round_mode -> t -> t -> t = "sub_float"
   external mul: round_mode -> t -> t -> t = "mul_float"
   external div: round_mode -> t -> t -> t = "div_float"
+  external sqrt: round_mode -> t -> t = "sqrt_float"
 
   external atof: round_mode -> string -> t = "atof_float"
 

src/cdomains/floatOps/floatOps.mli

@@ -38,6 +38,7 @@ module type CFloatType = sig
   val sub: round_mode -> t -> t -> t
   val mul: round_mode -> t -> t -> t
   val div: round_mode -> t -> t -> t
+  val sqrt: round_mode -> t -> t
   val atof: round_mode -> string -> t
 end
 

src/cdomains/floatOps/stubs.c

@@ -36,6 +36,20 @@ static void change_round_mode(int mode)
     }
 }
 
+#define UNARY_OP(name, type, op)                                 \
+    CAMLprim value name##_##type(value mode, value x)            \
+    {                                                            \
+        int old_roundingmode = fegetround();                     \
+        change_round_mode(Int_val(mode));                        \
+        volatile type r, x1 = Double_val(x);                     \
+        r = op(x1);                                              \
+        fesetround(old_roundingmode);                            \
+        return caml_copy_double(r);                              \
+    }
+
+UNARY_OP(sqrt, double, sqrt);
+UNARY_OP(sqrt, float, sqrtf);
+
 #define BINARY_OP(name, type, op)                                \
     CAMLprim value name##_##type(value mode, value x, value y)   \
     {                                                            \