From 34dedf4f613e618499432f0d1fb7aa2dfb9cf2ac Mon Sep 17 00:00:00 2001 From: Maxime Jacquemin <maxime.jacquemin@cea.fr> Date: Thu, 21 Dec 2017 12:22:29 +0100 Subject: [PATCH] [Eva - Numerors] Propagate module name and split the main file --- .../numerors_domain.ml} | 0 .../numerors_domain.mli} | 0 .../value/values/numerors/numerors_float.ml | 184 ++++++++ .../value/values/numerors/numerors_float.mli | 53 +++ .../value/values/numerors/numerors_utils.ml | 77 ++++ .../value/values/numerors/numerors_utils.mli | 43 ++ .../numerors_value.ko.ml} | 0 .../numerors_value.ml} | 414 +++++++----------- .../numerors_value.mli} | 0 .../numerors_value.ok.ml} | 414 +++++++----------- 10 files changed, 679 insertions(+), 506 deletions(-) rename src/plugins/value/domains/{errors_domain.ml => numerors/numerors_domain.ml} (100%) rename src/plugins/value/domains/{errors_domain.mli => numerors/numerors_domain.mli} (100%) create mode 100644 src/plugins/value/values/numerors/numerors_float.ml create mode 100644 src/plugins/value/values/numerors/numerors_float.mli create mode 100644 src/plugins/value/values/numerors/numerors_utils.ml create mode 100644 src/plugins/value/values/numerors/numerors_utils.mli rename src/plugins/value/values/{errors_value.ko.ml => numerors/numerors_value.ko.ml} (100%) rename src/plugins/value/values/{errors_value.ml => numerors/numerors_value.ml} (69%) rename src/plugins/value/values/{errors_value.mli => numerors/numerors_value.mli} (100%) rename src/plugins/value/values/{errors_value.ok.ml => numerors/numerors_value.ok.ml} (69%) diff --git a/src/plugins/value/domains/errors_domain.ml b/src/plugins/value/domains/numerors/numerors_domain.ml similarity index 100% rename from src/plugins/value/domains/errors_domain.ml rename to src/plugins/value/domains/numerors/numerors_domain.ml diff --git a/src/plugins/value/domains/errors_domain.mli b/src/plugins/value/domains/numerors/numerors_domain.mli similarity index 100% rename from src/plugins/value/domains/errors_domain.mli rename to src/plugins/value/domains/numerors/numerors_domain.mli diff --git a/src/plugins/value/values/numerors/numerors_float.ml b/src/plugins/value/values/numerors/numerors_float.ml new file mode 100644 index 00000000000..31d282d3e86 --- /dev/null +++ b/src/plugins/value/values/numerors/numerors_float.ml @@ -0,0 +1,184 @@ +open Numerors_utils + + +(*----------------------------------------------------------------------------- + * Float representation. The MPFR numbers are always positives. + *---------------------------------------------------------------------------*) +type t = Sign.t * num +and num = Nan | Inf | Float of Mpfrf.t + + +(*----------------------------------------------------------------------------- + * Pretty printer + *---------------------------------------------------------------------------*) +let pretty fmt (s, n) = + let pp_num fmt = function + | Nan -> Format.fprintf fmt "Nan" + | Inf -> Format.fprintf fmt "Inf" + | Float f -> Format.fprintf fmt "%a" Mpfrf.print f + in Format.fprintf fmt "%a%a" Sign.pretty s pp_num n + + +(*----------------------------------------------------------------------------- + * Internal functions + *---------------------------------------------------------------------------*) +(* Set asked precision and call the function *) +let set_precision p f = Mpfr.set_default_prec (Precisions.get p) ; f + +(* Rebuild the correct MPFR with the sign *) +let rebuild ?rnd s n = + let r = match rnd with Some r -> r | None -> Mpfr.Near in + match s with Sign.Positive -> n | Sign.Negative -> Mpfrf.neg n r + + +(*----------------------------------------------------------------------------- + * Methods to get informations on floats + *---------------------------------------------------------------------------*) +let sign_of (s, _) = s +let precision_of (_, n) = + match n with + | Nan | Inf -> None + | Float f -> Precisions.of_int (Mpfr.get_prec (Mpfrf.to_mpfr f)) + + +(*----------------------------------------------------------------------------- + * Constructors + *---------------------------------------------------------------------------*) +let nan = Sign.Positive, Nan +let mpfr_zero p = set_precision p Mpfrf.of_float 0.0 Mpfr.Zero +let pos_zero p = Sign.Positive, Float (mpfr_zero p) +let neg_zero p = Sign.Negative, Float (mpfr_zero p) +let pos_inf = Sign.Positive, Inf +let neg_inf = Sign.Negative, Inf + +let of_mpfr f = + let s = Sign.of_mpfr f in + let n = + if Mpfrf.inf_p f then Inf + else if Mpfrf.nan_p f then Nan + else Float (Mpfrf.abs f Mpfr.Near) + in s, n + +let of_int ?(rnd = Mpfr.Near) p i = + of_mpfr (set_precision p Mpfrf.of_int i rnd) + +let of_float ?(rnd = Mpfr.Near) p f = + of_mpfr (set_precision p Mpfrf.of_float f rnd) + +let of_string ?(rnd = Mpfr.Near) p s = + let s = List.hd (String.split_on_char 'f' s) in + of_mpfr (set_precision p Mpfrf.of_string s rnd) + + +(*----------------------------------------------------------------------------- + * Comparison methods + *---------------------------------------------------------------------------*) +let compare (sx, nx) (sy, ny) = + match nx, ny with + | Nan, Nan -> 0 | Nan, _ -> 1 | _, Nan -> -1 + | Inf, _ -> if Sign.eq sx Sign.Negative then -1 else 1 + | _, Inf -> if Sign.eq sy Sign.Positive then -1 else 1 + | Float fx, Float fy -> Mpfrf.cmp (rebuild sx fx) (rebuild sy fy) +let eq a b = compare a b = 0 +let le a b = compare a b <= 0 +let lt a b = compare a b < 0 +let ge a b = compare a b >= 0 +let gt a b = compare a b > 0 + + +(*----------------------------------------------------------------------------- + * Methods to check properties on floats + *---------------------------------------------------------------------------*) +let is_nan (_, n) = match n with Nan -> true | _ -> false +let is_inf (_, n) = match n with Inf -> true | _ -> false +let is_pos (s, _) = Sign.eq s Sign.Positive +let is_neg (s, _) = Sign.eq s Sign.Negative +let is_a_zero f = + match precision_of f with + | Some p -> eq f (pos_zero p) + | None -> false +let is_pos_zero f = is_pos f && is_a_zero f +let is_neg_zero f = is_neg f && is_a_zero f +let is_strictly_pos f = is_pos f && not (is_a_zero f) +let is_strictly_neg f = is_neg f && not (is_a_zero f) + + +(*----------------------------------------------------------------------------- + * Cast a float to the asked precision + *---------------------------------------------------------------------------*) +let cast ?(rnd = Mpfr.Near) ~prec (s, n) = + match n with + | Float x -> + let x' = Mpfrf.to_mpfr x in + let r = Mpfr.init () in + let old = Mpfr.get_default_rounding_mode () in + Mpfr.set_default_rounding_mode rnd ; + Mpfr.set_prec r (Precisions.get prec) ; + let _ = Mpfr.set r x' Mpfr.Near in + Mpfr.set_default_rounding_mode old ; + s, Float (Mpfrf.of_mpfr r) + | _ -> (s, n) + + +(*----------------------------------------------------------------------------- + * Return the second argument with the sign of the first + *---------------------------------------------------------------------------*) +let apply_sign (s, _) (_, n) = (s, n) + +(*----------------------------------------------------------------------------- + * Arithmetics + *----------------------------------------------------------------------------- + * The IEEE-754 norm does not specify the propagation of the sign bit with + * Nan operands. For the addition and the substraction, we set it to Positive. + * For the multiplication and the division, we set it to the exclusive or + * between the signs of the operands. + *---------------------------------------------------------------------------*) +let normal_op rnd p op (sa, fa) (sb, fb) = + let x, y = rebuild sa fa, rebuild sb fb in + of_mpfr (set_precision p op x y rnd) + +let add ?(rnd = Mpfr.Near) ~prec (sa, na) (sb, nb) = + match na, nb with + | Nan, _ | _, Nan -> Sign.Positive, Nan + | Inf, Inf when Sign.eq sa sb -> sa, Inf + | Inf, Inf -> Sign.Positive, Nan + | Inf, Float _ -> sa, Inf + | Float _, Inf -> sb, Inf + | Float fa, Float fb -> normal_op rnd prec Mpfrf.add (sa, fa) (sb, fb) + +let sub ?(rnd = Mpfr.Near) ~prec (sa, na) (sb, nb) = + match na, nb with + | Nan, _ | _, Nan -> Sign.Positive, Nan + | Inf, Inf when not (Sign.eq sa sb) -> sa, Inf + | Inf, Inf -> Sign.Positive, Nan + | Inf, Float _ -> sa, Inf + | Float _, Inf -> Sign.neg sb, Inf + | Float fa, Float fb -> normal_op rnd prec Mpfrf.sub (sa, fa) (sb, fb) + +let mul ?(rnd = Mpfr.Near) ~prec (sa, na) (sb, nb) = + match na, nb with + | Nan, _ | _, Nan -> Sign.mul sa sb, Nan + | Inf, _ | _, Inf -> Sign.mul sa sb, Inf + | Float fa, Float fb -> normal_op rnd prec Mpfrf.mul (sa, fa) (sb, fb) + +let div ?(rnd = Mpfr.Near) ~prec (sa, na) (sb, nb) = + match na, nb with + | Nan, _ | _, Nan -> Sign.mul sa sb, Nan + | Inf, Inf -> Sign.mul sa sb, Nan + | Inf, Float _ -> Sign.mul sa sb, Inf + | Float _, Inf -> Sign.mul sa sb, Float (mpfr_zero prec) + | Float fa, Float fb -> normal_op rnd prec Mpfrf.div (sa, fa) (sb, fb) + +let neg (sa, na) = Sign.neg sa , na +let abs (_ , na) = Sign.Positive, na + +let sqrt ?(rnd = Mpfr.Near) ~prec (sa, na) = + match na with + | Inf when Sign.eq sa Sign.Positive -> sa, Inf + | Float x when Sign.eq sa Sign.Positive -> + of_mpfr (set_precision prec Mpfrf.sqrt x rnd) + | _ -> sa, Nan + +(* Min and Max *) +let min x y = if compare x y <= 0 then x else y +let max x y = if compare x y <= 0 then y else x diff --git a/src/plugins/value/values/numerors/numerors_float.mli b/src/plugins/value/values/numerors/numerors_float.mli new file mode 100644 index 00000000000..23c90c7241e --- /dev/null +++ b/src/plugins/value/values/numerors/numerors_float.mli @@ -0,0 +1,53 @@ +open Numerors_utils + +type t + +val pretty : Format.formatter -> t -> unit + +val sign_of : t -> Sign.t +val precision_of : t -> Precisions.t option + +val nan : t +val pos_inf : t +val neg_inf : t +val pos_zero : Precisions.t -> t +val neg_zero : Precisions.t -> t + +val of_mpfr : Mpfrf.t -> t +val of_int : ?rnd:Mpfr.round -> Precisions.t -> int -> t +val of_float : ?rnd:Mpfr.round -> Precisions.t -> float -> t +val of_string : ?rnd:Mpfr.round -> Precisions.t -> string -> t + +val compare : t -> t -> int +val eq : t -> t -> bool +val le : t -> t -> bool +val lt : t -> t -> bool +val ge : t -> t -> bool +val gt : t -> t -> bool + +val is_nan : t -> bool +val is_inf : t -> bool +val is_pos : t -> bool +val is_neg : t -> bool +val is_a_zero : t -> bool +val is_pos_zero : t -> bool +val is_neg_zero : t -> bool +val is_strictly_pos : t -> bool +val is_strictly_neg : t -> bool + +val cast : ?rnd:Mpfr.round -> prec:Precisions.t -> t -> t + +val apply_sign : t -> t -> t + +val add : ?rnd:Mpfr.round -> prec:Precisions.t -> t -> t -> t +val sub : ?rnd:Mpfr.round -> prec:Precisions.t -> t -> t -> t +val mul : ?rnd:Mpfr.round -> prec:Precisions.t -> t -> t -> t +val div : ?rnd:Mpfr.round -> prec:Precisions.t -> t -> t -> t + +val neg : t -> t +val abs : t -> t + +val sqrt : ?rnd:Mpfr.round -> prec:Precisions.t -> t -> t + +val min : t -> t -> t +val max : t -> t -> t diff --git a/src/plugins/value/values/numerors/numerors_utils.ml b/src/plugins/value/values/numerors/numerors_utils.ml new file mode 100644 index 00000000000..e96224f95fb --- /dev/null +++ b/src/plugins/value/values/numerors/numerors_utils.ml @@ -0,0 +1,77 @@ +open Cil_types + +(*----------------------------------------------------------------------------- + * Module describing the different precisions that will be manipulated + *---------------------------------------------------------------------------*) +module Precisions = struct + + type t = Simple | Double | Long_Double | Real + + let rp : int = Value_parameters.RealSignificandBits.get () + + let pretty fmt = function + | Simple -> Format.fprintf fmt "Simple" + | Double -> Format.fprintf fmt "Double" + | Long_Double -> Format.fprintf fmt "Long Double" + | Real -> Format.fprintf fmt "Real" + + let of_fkind = function + | FFloat -> Simple + | FDouble -> Double + | FLongDouble -> Long_Double + + let get = function + | Simple -> 23 | Double -> 52 + | Long_Double -> 64 | Real -> rp + + let of_int i = + if i = get Simple then Some Simple + else if i = get Double then Some Double + else if i = get Long_Double then Some Long_Double + else if i = get Real then Some Real + else None + + let compare a b = Pervasives.compare (get a) (get b) + let eq a b = compare a b = 0 + let le a b = compare a b <= 0 + let lt a b = compare a b < 0 + let ge a b = compare a b >= 0 + let gt a b = compare a b > 0 + + let max a b = if compare a b <= 0 then b else a + let min a b = if compare a b <= 0 then a else b + +end + + +(*----------------------------------------------------------------------------- + * Sign type for infinites + *---------------------------------------------------------------------------*) +module Sign = struct + + type t = Positive | Negative + + let pretty fmt = function + | Positive -> Format.fprintf fmt "+" + | Negative -> Format.fprintf fmt "-" + + let of_int i = if i < 0 then Negative else Positive + let of_mpfr f = of_int (Mpfrf.sgn f) + + let compare a b = + match a, b with + | Positive, Positive | Negative, Negative -> 0 + | Positive, Negative -> 1 + | Negative, Positive -> -1 + let eq a b = compare a b = 0 + let le a b = compare a b <= 0 + let lt a b = compare a b < 0 + let ge a b = compare a b >= 0 + let gt a b = compare a b > 0 + + let neg s = if s = Positive then Negative else Positive + let mul a b = if eq a b then Positive else Negative + +end + + diff --git a/src/plugins/value/values/numerors/numerors_utils.mli b/src/plugins/value/values/numerors/numerors_utils.mli new file mode 100644 index 00000000000..2179797f12e --- /dev/null +++ b/src/plugins/value/values/numerors/numerors_utils.mli @@ -0,0 +1,43 @@ +module Precisions : sig + + type t = Simple | Double | Long_Double | Real + + val pretty : Format.formatter -> t -> unit + + val of_fkind : Cil_types.fkind -> t + val of_int : int -> t option + val get : t -> int + + val compare : t -> t -> int + val eq : t -> t -> bool + val le : t -> t -> bool + val lt : t -> t -> bool + val ge : t -> t -> bool + val gt : t -> t -> bool + + val max : t -> t -> t + val min : t -> t -> t + +end + + +module Sign : sig + + type t = Positive | Negative + + val pretty : Format.formatter -> t -> unit + + val of_int : int -> t + val of_mpfr : Mpfrf.t -> t + + val compare : t -> t -> int + val eq : t -> t -> bool + val le : t -> t -> bool + val lt : t -> t -> bool + val ge : t -> t -> bool + val gt : t -> t -> bool + + val neg : t -> t + val mul : t -> t -> t + +end diff --git a/src/plugins/value/values/errors_value.ko.ml b/src/plugins/value/values/numerors/numerors_value.ko.ml similarity index 100% rename from src/plugins/value/values/errors_value.ko.ml rename to src/plugins/value/values/numerors/numerors_value.ko.ml diff --git a/src/plugins/value/values/errors_value.ml b/src/plugins/value/values/numerors/numerors_value.ml similarity index 69% rename from src/plugins/value/values/errors_value.ml rename to src/plugins/value/values/numerors/numerors_value.ml index ab58c843bbd..8d080b16d7e 100755 --- a/src/plugins/value/values/errors_value.ml +++ b/src/plugins/value/values/numerors/numerors_value.ml @@ -22,207 +22,13 @@ open Cil_types open Eval - - -(*----------------------------------------------------------------------------- - * Module describing the different precisions that will be manipulated - *---------------------------------------------------------------------------*) -module Precisions = struct - type t = Simple | Double | Long_Double | Real - let rp = Value_parameters.RealSignificandBits.get () - let of_fkind = function - | FFloat -> Simple - | FDouble -> Double - | FLongDouble -> Long_Double - let get = function - | Simple -> 24 | Double -> 53 - | Long_Double -> 64 | Real -> rp - let pretty fmt = function - | Simple -> Format.fprintf fmt "Simple" - | Double -> Format.fprintf fmt "Double" - | Long_Double -> Format.fprintf fmt "Long Double" - | Real -> Format.fprintf fmt "Real" - let compare a b = Pervasives.compare (get a) (get b) - let le a b = compare a b <= 0 - let gt a b = compare a b > 0 - let equal a b = compare a b = 0 - let max a b = if compare a b <= 0 then b else a - let min a b = if compare a b <= 0 then a else b -end - - -(*----------------------------------------------------------------------------- - * Sign type for infinites - *---------------------------------------------------------------------------*) -module Sign = struct - type t = Positive | Negative - let pretty fmt = function - | Positive -> Format.fprintf fmt "+" - | Negative -> Format.fprintf fmt "-" - let compare a b = - match a, b with - | Positive, Positive | Negative, Negative -> 0 - | Positive, Negative -> 1 - | Negative, Positive -> -1 - let equal a b = compare a b = 0 - let of_int i = if i < 0 then Negative else Positive - let of_mpfr f = of_int (Mpfrf.sgn f) - let neg s = if s = Positive then Negative else Positive - let mul a b = if equal a b then Positive else Negative -end +open Numerors_utils (*----------------------------------------------------------------------------- * Module representing numbers with a specified precision *---------------------------------------------------------------------------*) -module F = struct - - (* Numbers representation. Float contains only positives MPFR numbers *) - type t = Sign.t * num - and num = Nan | Inf | Float of Mpfrf.t - - (* Getters on the representation *) - (* let precision_of (_, _, p) = p *) - let sign_of (s, _) = s - - (* Check if the number is Nan *) - let is_nan (_, n) = match n with Nan -> true | _ -> false - - (* Check if the number is Inf *) - let is_inf (_, n) = match n with Inf -> true | _ -> false - - (* Pretty printers *) - let pretty fmt (s, n) = - Format.fprintf fmt "%a%a" Sign.pretty s - ( fun fmt -> function - | Nan -> Format.fprintf fmt "Nan" - | Inf -> Format.fprintf fmt "Inf" - | Float f -> Format.fprintf fmt "%a" Mpfrf.print f - ) n - - (* Rebuild the correct MPFR with the sign *) - let rebuild ?rnd s n = - let r = match rnd with Some r -> r | None -> Mpfr.Near in - match s with Sign.Positive -> n | Sign.Negative -> Mpfrf.neg n r - - (* Comparisons. For elements of the same sign, Float < Inf < Nan. *) - let compare (sx, nx) (sy, ny) = - match nx, ny with - | Nan, Nan -> 0 | Nan, _ -> 1 | _, Nan -> -1 - | Inf, _ -> if Sign.equal sx Sign.Negative then -1 else 1 - | _, Inf -> if Sign.equal sy Sign.Positive then -1 else 1 - | Float fx, Float fy -> Mpfrf.cmp (rebuild sx fx) (rebuild sy fy) - - let le x y = compare x y <= 0 - let equal x y = compare x y = 0 - - (* Create a new element containing the value of the input but with - * the asked precision *) - let create_new_prec ~rnd ~prec (s, n) = - match n with - | Float x -> - let x' = Mpfrf.to_mpfr x in - let r = Mpfr.init () in - let old = Mpfr.get_default_rounding_mode () in - Mpfr.set_default_rounding_mode rnd ; - Mpfr.set_prec r (Precisions.get prec) ; - let _ = Mpfr.set r x' Mpfr.Near in - Mpfr.set_default_rounding_mode old ; - s, Float (Mpfrf.of_mpfr r) - | _ -> (s, n) - - (* Set asked precision and call the function *) - let set_precision p f = Mpfr.set_default_prec (Precisions.get p) ; f - - (* Create a new t element from a Mpfrf.t and a Precisions.t. - * Check Nan and Infinite. *) - let create f = - let s = Sign.of_mpfr f in - let n = - if Mpfrf.inf_p f then Inf - else if Mpfrf.nan_p f then Nan - else Float (Mpfrf.abs f Mpfr.Near) - in s, n - - (* Constructors *) - let mpfr_zero p = set_precision p Mpfrf.of_float 0.0 Mpfr.Zero - let pos_zero p = Sign.Positive, Float (mpfr_zero p) - let neg_zero p = Sign.Negative, Float (mpfr_zero p) - let pos_inf = Sign.Positive, Inf - let neg_inf = Sign.Negative, Inf - - let of_int ?(rnd = Mpfr.Near) p i = - create (set_precision p Mpfrf.of_int i rnd) - - let of_float ?(rnd = Mpfr.Near) p f = - create (set_precision p Mpfrf.of_float f rnd) - - let of_string ?(rnd = Mpfr.Near) p s = - let s = List.hd (String.split_on_char 'f' s) in - create (set_precision p Mpfrf.of_string s rnd) - - (* Internal function to choose the precision to use in calculations *) - let internal_prec prec pa pb = - match prec with - | None -> Precisions.max pa pb - | Some p -> p - - (* Arithmetics. The IEEE-754 norm does not specify the propagation - * of the sign bit with Nan operands. For the addition and the - * substraction, we set it to Positive. For the multiplication - * and the division, we set it to the exclusive or between the - * signs of the operands. *) - let normal_op rnd p op (sa, fa) (sb, fb) = - let x, y = rebuild sa fa, rebuild sb fb in - create (set_precision p op x y rnd) - - let add ?(rnd = Mpfr.Near) ~prec (sa, na) (sb, nb) = - match na, nb with - | Nan, _ | _, Nan -> Sign.Positive, Nan - | Inf, Inf when Sign.equal sa sb -> sa, Inf - | Inf, Inf -> Sign.Positive, Nan - | Inf, Float _ -> sa, Inf - | Float _, Inf -> sb, Inf - | Float fa, Float fb -> normal_op rnd prec Mpfrf.add (sa, fa) (sb, fb) - - let sub ?(rnd = Mpfr.Near) ~prec (sa, na) (sb, nb) = - match na, nb with - | Nan, _ | _, Nan -> Sign.Positive, Nan - | Inf, Inf when not (Sign.equal sa sb) -> sa, Inf - | Inf, Inf -> Sign.Positive, Nan - | Inf, Float _ -> sa, Inf - | Float _, Inf -> Sign.neg sb, Inf - | Float fa, Float fb -> normal_op rnd prec Mpfrf.sub (sa, fa) (sb, fb) - - let mul ?(rnd = Mpfr.Near) ~prec (sa, na) (sb, nb) = - match na, nb with - | Nan, _ | _, Nan -> Sign.mul sa sb, Nan - | Inf, _ | _, Inf -> Sign.mul sa sb, Inf - | Float fa, Float fb -> normal_op rnd prec Mpfrf.mul (sa, fa) (sb, fb) - - let div ?(rnd = Mpfr.Near) ~prec (sa, na) (sb, nb) = - match na, nb with - | Nan, _ | _, Nan -> Sign.mul sa sb, Nan - | Inf, Inf -> Sign.mul sa sb, Nan - | Inf, Float _ -> Sign.mul sa sb, Inf - | Float _, Inf -> Sign.mul sa sb, Float (mpfr_zero prec) - | Float fa, Float fb -> normal_op rnd prec Mpfrf.div (sa, fa) (sb, fb) - - let neg (sa, na) = Sign.neg sa , na - let abs (_ , na) = Sign.Positive, na - - let sqrt ?(rnd = Mpfr.Near) ~prec (sa, na) = - match na with - | Inf when Sign.equal sa Sign.Positive -> sa, Inf - | Float x when Sign.equal sa Sign.Positive -> - create (set_precision prec Mpfrf.sqrt x rnd) - | _ -> sa, Nan - - (* Min and Max *) - let min x y = if compare x y <= 0 then x else y - let max x y = if compare x y <= 0 then y else x - -end +module F = Numerors_float (*----------------------------------------------------------------------------- @@ -258,8 +64,8 @@ module I = struct if Precisions.le p prec then Mpfr.Near, Mpfr.Near else Mpfr.Up, Mpfr.Down in - let x = F.create_new_prec ~rnd:down ~prec:prec x in - let y = F.create_new_prec ~rnd:up ~prec:prec y in + let x = F.cast ~rnd:down ~prec:prec x in + let y = F.cast ~rnd:up ~prec:prec y in let r = make ~nan:n x y in r, prec @@ -298,7 +104,7 @@ module I = struct (* Partial order, two intervals with different precisions * are not comparable. *) let is_included (a, pa) (b, pb) = - if Precisions.equal pa pb then + if Precisions.eq pa pb then match a, b with | I (x, y, n), I (x', y', n') -> F.le x x' && F.le y y' && (not n || n') @@ -408,10 +214,10 @@ module I = struct let add r = results := r :: !results in let treat_nan x y = nan := true ; - if F.is_inf x && not (F.equal b1 e1) then - add (F.sign_of x, F.Float (Mpfrf.of_int 1 Mpfr.Near)) ; - if F.is_inf y && not (F.equal b2 e2) then - add (F.sign_of y, F.Float (Mpfrf.of_int 1 Mpfr.Near)) ; + if F.is_inf x && not (F.eq b1 e1) then + add (F.apply_sign x (F.of_int prec 1)) ; + if F.is_inf y && not (F.eq b2 e2) then + add (F.apply_sign y (F.of_int prec 1)) ; in let op rnd x y = let r = op ~rnd:rnd ~prec:p x y in @@ -467,17 +273,28 @@ module I = struct let r = if has_ninf then join ~prec:prec (neg_inf ~prec:prec) r else r in if nan then add_nan r else r + let abs ~prec a = + let a, _ = cast ~prec:prec a in + let r = a >>- fun (x, y, n) -> + if contains_pos_zero ~prec:prec (a, prec) then + make ~nan:n (F.pos_zero prec) (F.max (F.abs x) (F.abs y)) + else if F.compare y (F.pos_zero prec)< 0 then + make ~nan:n (F.neg y) (F.neg x) + else a + in r, prec + (* Min and Max *) let max_abs ~prec a = - let a, _ = cast ~prec:prec a in + let a, _ = abs ~prec:prec a in match a with - | Nan -> Sign.Positive, F.Nan - | I (x, y, _) -> F.max (F.abs x) (F.abs y) + | Nan | I (_, _, true) -> F.nan + | I (x, y, _) -> if F.compare x y < 0 then y else x let min_abs ~prec a = - let a, _ = cast ~prec:prec a in + let a, _ = abs ~prec:prec a in match a with - | Nan -> Sign.Positive, F.Nan - | I (x, y, _) -> F.min (F.abs x) (F.abs y) + | Nan | I (_, _, true) -> F.nan + | I (x, y, _) -> if F.compare x y < 0 then x else y + end @@ -486,10 +303,17 @@ end * Constants definitions *---------------------------------------------------------------------------*) module Cst = struct - let e = F.of_float Precisions.Real (2.0 ** (~-.53.0)) - let e_interval = I.make (F.neg e) e, Precisions.Real - let one_interval = I.of_floats ~prec:Precisions.Real (1.0, 1.0) - let std_error = I.add ~prec:Precisions.Real one_interval e_interval + let machine_epsilon p = + let m = Pervasives.float_of_int (Precisions.get p) in + F.of_float Precisions.Real (2.0 ** (~-.m)) + let eps_interval p = + let e = machine_epsilon p in + I.make (F.neg e) e, Precisions.Real + let one_interval = + I.of_floats ~prec:Precisions.Real (1.0, 1.0) + let std_error p = + let eps_interval = eps_interval p in + I.add ~prec:Precisions.Real one_interval eps_interval end @@ -511,7 +335,7 @@ type err = } let get_prec (_, pa) (_, pb) = - if Precisions.equal pa pb then pa else assert false + if Precisions.eq pa pb then pa else assert false (* Lattice Structure. Normally, we should implement a Lattice for each * C float type, because it change the precision of the approx field. @@ -564,6 +388,7 @@ let narrow x y = `Value { exact ; approx ; abs_err ; rel_err } | _, _, _, _ -> `Bottom +(* let apply_constraints exact app_opt in_abs in_rel = match app_opt with | None -> { exact ; approx = None ; abs_err = in_abs ; rel_err = in_rel } @@ -571,41 +396,58 @@ let apply_constraints exact app_opt in_abs in_rel = let extract t = function `Value v -> v | `Bottom -> t in let try_narrow x y = extract x (I.narrow ~prec:Precisions.Real x y) in let abs_from_vals = I.sub ~prec:Precisions.Real approx exact in + let abs_from_rel = in_abs in + (* let abs_from_rel = I.mul ~prec:Precisions.Real in_rel exact in + *) let rel_from_vals = I.div ~prec:Precisions.Real abs_from_vals exact in - let rel_from_abs = I.div ~prec:Precisions.Real in_abs exact in + let rel_from_abs = in_rel in + (* + let rel_from_abs = I.div ~prec:Precisions.Real in_abs exact in + *) let rel_err = try_narrow (try_narrow in_rel rel_from_vals) rel_from_abs in let abs_err = try_narrow (try_narrow in_abs abs_from_vals) abs_from_rel in + (* + let rel_err = try_narrow in_rel rel_from_vals in + let abs_err = try_narrow in_abs abs_from_vals in + *) (* Debug purpose, will disappear. *) - let is_better x y = I.compare x y != 0 && I.is_included x y in - if is_better abs_err in_abs || is_better rel_err in_rel then - Format.printf - "Constraints applications :@.@[<v>\ - Input Exact : %a@.Input Approx : %a@.\ - Input Abs : %a@.Input Rel : %a@.\ - Abs From Vals : %a@.Abs From Rel : %a@.\ - Rel From Vals : %a@.Rel From Abs : %a@.\ - Output Abs : %a@.Output Rel : %a@.\ - @]@." - I.pretty exact I.pretty approx - I.pretty in_abs I.pretty in_rel - I.pretty abs_from_vals I.pretty abs_from_rel - I.pretty rel_from_vals I.pretty rel_from_abs - I.pretty abs_err I.pretty rel_err ; + (* + Format.printf + "Constraints applications :@.@[<v>\ + Input Exact : %a@.Input Approx : %a@.\ + Input Abs : %a@.Input Rel : %a@.\ + Abs From Vals : %a@.Abs From Rel : %a@.\ + Rel From Vals : %a@.Rel From Abs : %a@.\ + Output Abs : %a@.Output Rel : %a@.\ + @]@." + I.pretty exact I.pretty approx + I.pretty in_abs I.pretty in_rel + I.pretty abs_from_vals I.pretty abs_from_rel + I.pretty rel_from_vals I.pretty rel_from_abs + I.pretty abs_err I.pretty rel_err ; + *) (* End of debugging part. *) { exact ; approx = Some approx ; abs_err ; rel_err } +*) let reduce fval t = match Fval.min_and_max fval, t.approx with | (Some (min, max), false), Some app -> let p = I.precision_of app in let itv = I.of_floats p (Fval.F.to_float min, Fval.F.to_float max) in + I.narrow ~prec:p app itv >>-: fun app -> { t with approx = Some app } + (* I.narrow ~prec:p app itv >>-: fun approx -> apply_constraints t.exact (Some approx) t.abs_err t.rel_err + *) | (Some (min, max), false), None -> let itv = Fval.F.to_float min, Fval.F.to_float max in let itv = I.of_floats Precisions.Double itv in + `Value { t with approx = Some itv } + (* `Value (apply_constraints t.exact (Some itv) t.abs_err t.rel_err) + *) | _ -> `Value t (* Associated Datatype *) @@ -730,7 +572,7 @@ module Approx_Sem = let handle_prec op x y = let px = I.precision_of x in let py = I.precision_of y in - if Precisions.equal px py then op ~prec:px x y + if Precisions.eq px py then op ~prec:px x y else assert false let add x y = match x.approx, y.approx with @@ -759,29 +601,33 @@ module Abs_Err_Sem : Semantic = let sqrt v = match Approx_Sem.sqrt v with | None -> I.top ~prec:Precisions.Real - | Some approx -> Reals.sub approx (Exact_Sem.sqrt v) + | Some approx -> Reals.sub approx (Exact_Sem.sqrt v) let add x y = (x, y) >>-% fun x_approx y_approx -> + let e_interval = Cst.eps_interval (I.precision_of x_approx) in let abs_err = Reals.add x.abs_err y.abs_err in let approx = Reals.add x_approx y_approx in - let t = Reals.mul approx Cst.e_interval in + let t = Reals.mul approx e_interval in Reals.add abs_err t let sub x y = (x, y) >>-% fun x_approx y_approx -> + let e_interval = Cst.eps_interval (I.precision_of x_approx) in let abs_err = Reals.sub x.abs_err y.abs_err in let approx = Reals.sub x_approx y_approx in - let t = Reals.mul approx Cst.e_interval in + let t = Reals.mul approx e_interval in Reals.add abs_err t let mul x y = (x, y) >>-% fun x_approx y_approx -> + let e_interval = Cst.eps_interval (I.precision_of x_approx) in let approx = Reals.mul x_approx y_approx in - let mul = Reals.mul approx Cst.e_interval in + let mul = Reals.mul approx e_interval in let xey = Reals.mul x.exact y.abs_err in let yex = Reals.mul y.exact x.abs_err in let exy = Reals.mul x.abs_err y.abs_err in Reals.add (Reals.add mul (Reals.add xey yex)) exy let div x y = (x, y) >>-% fun x_approx y_approx -> + let e_interval = Cst.eps_interval (I.precision_of x_approx) in let xy = Reals.div x.exact y.exact in let eax = x.abs_err in let eay = Reals.mul xy y.abs_err in - let phi = Reals.mul x_approx Cst.e_interval in + let phi = Reals.mul x_approx e_interval in Reals.div (Reals.add (Reals.sub eax eay) phi) y_approx end @@ -790,32 +636,71 @@ module Abs_Err_Sem : Semantic = * during the forward operations *) module Rel_Err_Sem : Semantic = struct + let get_prec v = match v.approx with + | None -> Precisions.Real + | Some t -> I.precision_of t let neg v = Reals.neg v.rel_err let sqrt v = + let std_error = Cst.std_error (get_prec v) in let s = Reals.sqrt (Reals.add v.rel_err Cst.one_interval) in - Reals.sub (Reals.mul s Cst.std_error) Cst.one_interval + Reals.sub (Reals.mul s std_error) Cst.one_interval + let add x y = - let d = F.max (Reals.max_abs x.rel_err) (Reals.max_abs y.rel_err) in - let num = R.add (Reals.max_abs x.exact) (Reals.max_abs y.exact) in - let den = Reals.min_abs (Reals.add x.exact y.exact) in - let eps = R.add (F.of_float Precisions.Real 1.0) Cst.e in - let v = R.add (R.mul (R.mul (R.div num den) d) eps) Cst.e in - I.make (R.neg v) v, Precisions.Real + let r1 = (x, y) >>-% fun x_approx y_approx -> + let e = Cst.eps_interval (get_prec x) in + let div = Reals.add x.exact y.exact in + let a = Reals.mul x.exact x.rel_err in + let b = Reals.mul y.exact y.rel_err in + let f = Reals.mul (Reals.add x_approx y_approx) e in + Reals.div (Reals.add (Reals.add a b) f) div + in + let r2 = + let e = Cst.machine_epsilon (get_prec x) in + let d = F.max (Reals.max_abs x.rel_err) (Reals.max_abs y.rel_err) in + let num = R.add (Reals.max_abs x.exact) (Reals.max_abs y.exact) in + let den = Reals.min_abs (Reals.add x.exact y.exact) in + let eps = R.add (F.of_float Precisions.Real 1.0) e in + let v = R.add (R.mul (R.mul (R.div num den) d) eps) e in + I.make (R.neg v) v, Precisions.Real + in + Format.printf "Add results :@[<v>%a@ %a@]@." I.pretty r1 I.pretty r2 ; + match I.narrow ~prec:Precisions.Real r1 r2 with + | `Bottom -> assert false + | `Value v -> v + let sub x y = - let d = F.max (Reals.max_abs x.rel_err) (Reals.max_abs y.rel_err) in - let num = R.add (Reals.max_abs x.exact) (Reals.max_abs y.exact) in - let den = Reals.min_abs (Reals.sub x.exact y.exact) in - let eps = R.add (F.of_float Precisions.Real 1.0) Cst.e in - let v = R.add (R.mul (R.mul (R.div num den) d) eps) Cst.e in - I.make (R.neg v) v, Precisions.Real + let r1 = (x, y) >>-% fun x_approx y_approx -> + let e = Cst.eps_interval (get_prec x) in + let div = Reals.sub x.exact y.exact in + let a = Reals.mul x.exact x.rel_err in + let b = Reals.mul y.exact y.rel_err in + let f = Reals.mul (Reals.sub x_approx y_approx) e in + Reals.div (Reals.add (Reals.sub a b) f) div + in + let r2 = + let e = Cst.machine_epsilon (get_prec x) in + let d = F.max (Reals.max_abs x.rel_err) (Reals.max_abs y.rel_err) in + let num = R.add (Reals.max_abs x.exact) (Reals.max_abs y.exact) in + let den = Reals.min_abs (Reals.sub x.exact y.exact) in + let eps = R.add (F.of_float Precisions.Real 1.0) e in + let v = R.add (R.mul (R.mul (R.div num den) d) eps) e in + I.make (R.neg v) v, Precisions.Real + in + Format.printf "Sub results :@[<v>%a@ %a@]@." I.pretty r1 I.pretty r2 ; + match I.narrow ~prec:Precisions.Real r1 r2 with + | `Bottom -> assert false + | `Value v -> v + let mul x y = + let std_error = Cst.std_error (get_prec x) in let ex = Reals.add Cst.one_interval x.rel_err in let ey = Reals.add Cst.one_interval y.rel_err in - Reals.sub (Reals.mul (Reals.mul ex ey) Cst.std_error) Cst.one_interval + Reals.sub (Reals.mul (Reals.mul ex ey) std_error) Cst.one_interval let div x y = + let std_error = Cst.std_error (get_prec x) in let ex = Reals.add Cst.one_interval x.rel_err in let ey = Reals.add Cst.one_interval y.rel_err in - Reals.sub (Reals.mul (Reals.div ex ey) Cst.std_error) Cst.one_interval + Reals.sub (Reals.mul (Reals.div ex ey) std_error) Cst.one_interval end @@ -832,9 +717,16 @@ let constant _ = function | Some s -> I.of_strings Precisions.Real (s, s) | None -> I.of_floats Precisions.Real (r, r) in - let top_real = I.top Precisions.Real in let approx = Some (I.of_floats prec (r, r)) in - apply_constraints exact approx top_real top_real + let eps = Cst.eps_interval Precisions.Real in + let abs_err = Reals.mul eps exact in + let rel_err = eps in + { exact ; approx ; abs_err ; rel_err } + (* + let top_real = I.top Precisions.Real in + let approx = I.of_floats prec (r, r) in + apply_constraints exact (Some approx) top_real top_real + *) | _ -> top (* Forward operations *) @@ -850,20 +742,36 @@ let forward_binop ~context:_ _ op x y = match op with | PlusA -> let exact , approx = Exact_Sem.add x y, Approx_Sem.add x y in + let abs_err, rel_err = Abs_Err_Sem.add x y, Rel_Err_Sem.add x y in + return { exact ; approx ; abs_err ; rel_err } + (* let prg_abs, prg_rel = Abs_Err_Sem.add x y, Rel_Err_Sem.add x y in return (apply_constraints exact approx prg_abs prg_rel) + *) | MinusA -> let exact , approx = Exact_Sem.sub x y, Approx_Sem.sub x y in + let abs_err, rel_err = Abs_Err_Sem.sub x y, Rel_Err_Sem.sub x y in + return { exact ; approx ; abs_err ; rel_err } + (* let prg_abs, prg_rel = Abs_Err_Sem.sub x y, Rel_Err_Sem.sub x y in return (apply_constraints exact approx prg_abs prg_rel) + *) | Mult -> let exact , approx = Exact_Sem.mul x y, Approx_Sem.mul x y in + let abs_err, rel_err = Abs_Err_Sem.mul x y, Rel_Err_Sem.mul x y in + return { exact ; approx ; abs_err ; rel_err } + (* let prg_abs, prg_rel = Abs_Err_Sem.mul x y, Rel_Err_Sem.mul x y in return (apply_constraints exact approx prg_abs prg_rel) + *) | Div -> let exact , approx = Exact_Sem.div x y, Approx_Sem.div x y in + let abs_err, rel_err = Abs_Err_Sem.div x y, Rel_Err_Sem.div x y in + return { exact ; approx ; abs_err ; rel_err } + (* let prg_abs, prg_rel = Abs_Err_Sem.div x y, Rel_Err_Sem.div x y in return (apply_constraints exact approx prg_abs prg_rel) + *) | _ -> return top (* Backward operations *) diff --git a/src/plugins/value/values/errors_value.mli b/src/plugins/value/values/numerors/numerors_value.mli similarity index 100% rename from src/plugins/value/values/errors_value.mli rename to src/plugins/value/values/numerors/numerors_value.mli diff --git a/src/plugins/value/values/errors_value.ok.ml b/src/plugins/value/values/numerors/numerors_value.ok.ml similarity index 69% rename from src/plugins/value/values/errors_value.ok.ml rename to src/plugins/value/values/numerors/numerors_value.ok.ml index ab58c843bbd..8d080b16d7e 100755 --- a/src/plugins/value/values/errors_value.ok.ml +++ b/src/plugins/value/values/numerors/numerors_value.ok.ml @@ -22,207 +22,13 @@ open Cil_types open Eval - - -(*----------------------------------------------------------------------------- - * Module describing the different precisions that will be manipulated - *---------------------------------------------------------------------------*) -module Precisions = struct - type t = Simple | Double | Long_Double | Real - let rp = Value_parameters.RealSignificandBits.get () - let of_fkind = function - | FFloat -> Simple - | FDouble -> Double - | FLongDouble -> Long_Double - let get = function - | Simple -> 24 | Double -> 53 - | Long_Double -> 64 | Real -> rp - let pretty fmt = function - | Simple -> Format.fprintf fmt "Simple" - | Double -> Format.fprintf fmt "Double" - | Long_Double -> Format.fprintf fmt "Long Double" - | Real -> Format.fprintf fmt "Real" - let compare a b = Pervasives.compare (get a) (get b) - let le a b = compare a b <= 0 - let gt a b = compare a b > 0 - let equal a b = compare a b = 0 - let max a b = if compare a b <= 0 then b else a - let min a b = if compare a b <= 0 then a else b -end - - -(*----------------------------------------------------------------------------- - * Sign type for infinites - *---------------------------------------------------------------------------*) -module Sign = struct - type t = Positive | Negative - let pretty fmt = function - | Positive -> Format.fprintf fmt "+" - | Negative -> Format.fprintf fmt "-" - let compare a b = - match a, b with - | Positive, Positive | Negative, Negative -> 0 - | Positive, Negative -> 1 - | Negative, Positive -> -1 - let equal a b = compare a b = 0 - let of_int i = if i < 0 then Negative else Positive - let of_mpfr f = of_int (Mpfrf.sgn f) - let neg s = if s = Positive then Negative else Positive - let mul a b = if equal a b then Positive else Negative -end +open Numerors_utils (*----------------------------------------------------------------------------- * Module representing numbers with a specified precision *---------------------------------------------------------------------------*) -module F = struct - - (* Numbers representation. Float contains only positives MPFR numbers *) - type t = Sign.t * num - and num = Nan | Inf | Float of Mpfrf.t - - (* Getters on the representation *) - (* let precision_of (_, _, p) = p *) - let sign_of (s, _) = s - - (* Check if the number is Nan *) - let is_nan (_, n) = match n with Nan -> true | _ -> false - - (* Check if the number is Inf *) - let is_inf (_, n) = match n with Inf -> true | _ -> false - - (* Pretty printers *) - let pretty fmt (s, n) = - Format.fprintf fmt "%a%a" Sign.pretty s - ( fun fmt -> function - | Nan -> Format.fprintf fmt "Nan" - | Inf -> Format.fprintf fmt "Inf" - | Float f -> Format.fprintf fmt "%a" Mpfrf.print f - ) n - - (* Rebuild the correct MPFR with the sign *) - let rebuild ?rnd s n = - let r = match rnd with Some r -> r | None -> Mpfr.Near in - match s with Sign.Positive -> n | Sign.Negative -> Mpfrf.neg n r - - (* Comparisons. For elements of the same sign, Float < Inf < Nan. *) - let compare (sx, nx) (sy, ny) = - match nx, ny with - | Nan, Nan -> 0 | Nan, _ -> 1 | _, Nan -> -1 - | Inf, _ -> if Sign.equal sx Sign.Negative then -1 else 1 - | _, Inf -> if Sign.equal sy Sign.Positive then -1 else 1 - | Float fx, Float fy -> Mpfrf.cmp (rebuild sx fx) (rebuild sy fy) - - let le x y = compare x y <= 0 - let equal x y = compare x y = 0 - - (* Create a new element containing the value of the input but with - * the asked precision *) - let create_new_prec ~rnd ~prec (s, n) = - match n with - | Float x -> - let x' = Mpfrf.to_mpfr x in - let r = Mpfr.init () in - let old = Mpfr.get_default_rounding_mode () in - Mpfr.set_default_rounding_mode rnd ; - Mpfr.set_prec r (Precisions.get prec) ; - let _ = Mpfr.set r x' Mpfr.Near in - Mpfr.set_default_rounding_mode old ; - s, Float (Mpfrf.of_mpfr r) - | _ -> (s, n) - - (* Set asked precision and call the function *) - let set_precision p f = Mpfr.set_default_prec (Precisions.get p) ; f - - (* Create a new t element from a Mpfrf.t and a Precisions.t. - * Check Nan and Infinite. *) - let create f = - let s = Sign.of_mpfr f in - let n = - if Mpfrf.inf_p f then Inf - else if Mpfrf.nan_p f then Nan - else Float (Mpfrf.abs f Mpfr.Near) - in s, n - - (* Constructors *) - let mpfr_zero p = set_precision p Mpfrf.of_float 0.0 Mpfr.Zero - let pos_zero p = Sign.Positive, Float (mpfr_zero p) - let neg_zero p = Sign.Negative, Float (mpfr_zero p) - let pos_inf = Sign.Positive, Inf - let neg_inf = Sign.Negative, Inf - - let of_int ?(rnd = Mpfr.Near) p i = - create (set_precision p Mpfrf.of_int i rnd) - - let of_float ?(rnd = Mpfr.Near) p f = - create (set_precision p Mpfrf.of_float f rnd) - - let of_string ?(rnd = Mpfr.Near) p s = - let s = List.hd (String.split_on_char 'f' s) in - create (set_precision p Mpfrf.of_string s rnd) - - (* Internal function to choose the precision to use in calculations *) - let internal_prec prec pa pb = - match prec with - | None -> Precisions.max pa pb - | Some p -> p - - (* Arithmetics. The IEEE-754 norm does not specify the propagation - * of the sign bit with Nan operands. For the addition and the - * substraction, we set it to Positive. For the multiplication - * and the division, we set it to the exclusive or between the - * signs of the operands. *) - let normal_op rnd p op (sa, fa) (sb, fb) = - let x, y = rebuild sa fa, rebuild sb fb in - create (set_precision p op x y rnd) - - let add ?(rnd = Mpfr.Near) ~prec (sa, na) (sb, nb) = - match na, nb with - | Nan, _ | _, Nan -> Sign.Positive, Nan - | Inf, Inf when Sign.equal sa sb -> sa, Inf - | Inf, Inf -> Sign.Positive, Nan - | Inf, Float _ -> sa, Inf - | Float _, Inf -> sb, Inf - | Float fa, Float fb -> normal_op rnd prec Mpfrf.add (sa, fa) (sb, fb) - - let sub ?(rnd = Mpfr.Near) ~prec (sa, na) (sb, nb) = - match na, nb with - | Nan, _ | _, Nan -> Sign.Positive, Nan - | Inf, Inf when not (Sign.equal sa sb) -> sa, Inf - | Inf, Inf -> Sign.Positive, Nan - | Inf, Float _ -> sa, Inf - | Float _, Inf -> Sign.neg sb, Inf - | Float fa, Float fb -> normal_op rnd prec Mpfrf.sub (sa, fa) (sb, fb) - - let mul ?(rnd = Mpfr.Near) ~prec (sa, na) (sb, nb) = - match na, nb with - | Nan, _ | _, Nan -> Sign.mul sa sb, Nan - | Inf, _ | _, Inf -> Sign.mul sa sb, Inf - | Float fa, Float fb -> normal_op rnd prec Mpfrf.mul (sa, fa) (sb, fb) - - let div ?(rnd = Mpfr.Near) ~prec (sa, na) (sb, nb) = - match na, nb with - | Nan, _ | _, Nan -> Sign.mul sa sb, Nan - | Inf, Inf -> Sign.mul sa sb, Nan - | Inf, Float _ -> Sign.mul sa sb, Inf - | Float _, Inf -> Sign.mul sa sb, Float (mpfr_zero prec) - | Float fa, Float fb -> normal_op rnd prec Mpfrf.div (sa, fa) (sb, fb) - - let neg (sa, na) = Sign.neg sa , na - let abs (_ , na) = Sign.Positive, na - - let sqrt ?(rnd = Mpfr.Near) ~prec (sa, na) = - match na with - | Inf when Sign.equal sa Sign.Positive -> sa, Inf - | Float x when Sign.equal sa Sign.Positive -> - create (set_precision prec Mpfrf.sqrt x rnd) - | _ -> sa, Nan - - (* Min and Max *) - let min x y = if compare x y <= 0 then x else y - let max x y = if compare x y <= 0 then y else x - -end +module F = Numerors_float (*----------------------------------------------------------------------------- @@ -258,8 +64,8 @@ module I = struct if Precisions.le p prec then Mpfr.Near, Mpfr.Near else Mpfr.Up, Mpfr.Down in - let x = F.create_new_prec ~rnd:down ~prec:prec x in - let y = F.create_new_prec ~rnd:up ~prec:prec y in + let x = F.cast ~rnd:down ~prec:prec x in + let y = F.cast ~rnd:up ~prec:prec y in let r = make ~nan:n x y in r, prec @@ -298,7 +104,7 @@ module I = struct (* Partial order, two intervals with different precisions * are not comparable. *) let is_included (a, pa) (b, pb) = - if Precisions.equal pa pb then + if Precisions.eq pa pb then match a, b with | I (x, y, n), I (x', y', n') -> F.le x x' && F.le y y' && (not n || n') @@ -408,10 +214,10 @@ module I = struct let add r = results := r :: !results in let treat_nan x y = nan := true ; - if F.is_inf x && not (F.equal b1 e1) then - add (F.sign_of x, F.Float (Mpfrf.of_int 1 Mpfr.Near)) ; - if F.is_inf y && not (F.equal b2 e2) then - add (F.sign_of y, F.Float (Mpfrf.of_int 1 Mpfr.Near)) ; + if F.is_inf x && not (F.eq b1 e1) then + add (F.apply_sign x (F.of_int prec 1)) ; + if F.is_inf y && not (F.eq b2 e2) then + add (F.apply_sign y (F.of_int prec 1)) ; in let op rnd x y = let r = op ~rnd:rnd ~prec:p x y in @@ -467,17 +273,28 @@ module I = struct let r = if has_ninf then join ~prec:prec (neg_inf ~prec:prec) r else r in if nan then add_nan r else r + let abs ~prec a = + let a, _ = cast ~prec:prec a in + let r = a >>- fun (x, y, n) -> + if contains_pos_zero ~prec:prec (a, prec) then + make ~nan:n (F.pos_zero prec) (F.max (F.abs x) (F.abs y)) + else if F.compare y (F.pos_zero prec)< 0 then + make ~nan:n (F.neg y) (F.neg x) + else a + in r, prec + (* Min and Max *) let max_abs ~prec a = - let a, _ = cast ~prec:prec a in + let a, _ = abs ~prec:prec a in match a with - | Nan -> Sign.Positive, F.Nan - | I (x, y, _) -> F.max (F.abs x) (F.abs y) + | Nan | I (_, _, true) -> F.nan + | I (x, y, _) -> if F.compare x y < 0 then y else x let min_abs ~prec a = - let a, _ = cast ~prec:prec a in + let a, _ = abs ~prec:prec a in match a with - | Nan -> Sign.Positive, F.Nan - | I (x, y, _) -> F.min (F.abs x) (F.abs y) + | Nan | I (_, _, true) -> F.nan + | I (x, y, _) -> if F.compare x y < 0 then x else y + end @@ -486,10 +303,17 @@ end * Constants definitions *---------------------------------------------------------------------------*) module Cst = struct - let e = F.of_float Precisions.Real (2.0 ** (~-.53.0)) - let e_interval = I.make (F.neg e) e, Precisions.Real - let one_interval = I.of_floats ~prec:Precisions.Real (1.0, 1.0) - let std_error = I.add ~prec:Precisions.Real one_interval e_interval + let machine_epsilon p = + let m = Pervasives.float_of_int (Precisions.get p) in + F.of_float Precisions.Real (2.0 ** (~-.m)) + let eps_interval p = + let e = machine_epsilon p in + I.make (F.neg e) e, Precisions.Real + let one_interval = + I.of_floats ~prec:Precisions.Real (1.0, 1.0) + let std_error p = + let eps_interval = eps_interval p in + I.add ~prec:Precisions.Real one_interval eps_interval end @@ -511,7 +335,7 @@ type err = } let get_prec (_, pa) (_, pb) = - if Precisions.equal pa pb then pa else assert false + if Precisions.eq pa pb then pa else assert false (* Lattice Structure. Normally, we should implement a Lattice for each * C float type, because it change the precision of the approx field. @@ -564,6 +388,7 @@ let narrow x y = `Value { exact ; approx ; abs_err ; rel_err } | _, _, _, _ -> `Bottom +(* let apply_constraints exact app_opt in_abs in_rel = match app_opt with | None -> { exact ; approx = None ; abs_err = in_abs ; rel_err = in_rel } @@ -571,41 +396,58 @@ let apply_constraints exact app_opt in_abs in_rel = let extract t = function `Value v -> v | `Bottom -> t in let try_narrow x y = extract x (I.narrow ~prec:Precisions.Real x y) in let abs_from_vals = I.sub ~prec:Precisions.Real approx exact in + let abs_from_rel = in_abs in + (* let abs_from_rel = I.mul ~prec:Precisions.Real in_rel exact in + *) let rel_from_vals = I.div ~prec:Precisions.Real abs_from_vals exact in - let rel_from_abs = I.div ~prec:Precisions.Real in_abs exact in + let rel_from_abs = in_rel in + (* + let rel_from_abs = I.div ~prec:Precisions.Real in_abs exact in + *) let rel_err = try_narrow (try_narrow in_rel rel_from_vals) rel_from_abs in let abs_err = try_narrow (try_narrow in_abs abs_from_vals) abs_from_rel in + (* + let rel_err = try_narrow in_rel rel_from_vals in + let abs_err = try_narrow in_abs abs_from_vals in + *) (* Debug purpose, will disappear. *) - let is_better x y = I.compare x y != 0 && I.is_included x y in - if is_better abs_err in_abs || is_better rel_err in_rel then - Format.printf - "Constraints applications :@.@[<v>\ - Input Exact : %a@.Input Approx : %a@.\ - Input Abs : %a@.Input Rel : %a@.\ - Abs From Vals : %a@.Abs From Rel : %a@.\ - Rel From Vals : %a@.Rel From Abs : %a@.\ - Output Abs : %a@.Output Rel : %a@.\ - @]@." - I.pretty exact I.pretty approx - I.pretty in_abs I.pretty in_rel - I.pretty abs_from_vals I.pretty abs_from_rel - I.pretty rel_from_vals I.pretty rel_from_abs - I.pretty abs_err I.pretty rel_err ; + (* + Format.printf + "Constraints applications :@.@[<v>\ + Input Exact : %a@.Input Approx : %a@.\ + Input Abs : %a@.Input Rel : %a@.\ + Abs From Vals : %a@.Abs From Rel : %a@.\ + Rel From Vals : %a@.Rel From Abs : %a@.\ + Output Abs : %a@.Output Rel : %a@.\ + @]@." + I.pretty exact I.pretty approx + I.pretty in_abs I.pretty in_rel + I.pretty abs_from_vals I.pretty abs_from_rel + I.pretty rel_from_vals I.pretty rel_from_abs + I.pretty abs_err I.pretty rel_err ; + *) (* End of debugging part. *) { exact ; approx = Some approx ; abs_err ; rel_err } +*) let reduce fval t = match Fval.min_and_max fval, t.approx with | (Some (min, max), false), Some app -> let p = I.precision_of app in let itv = I.of_floats p (Fval.F.to_float min, Fval.F.to_float max) in + I.narrow ~prec:p app itv >>-: fun app -> { t with approx = Some app } + (* I.narrow ~prec:p app itv >>-: fun approx -> apply_constraints t.exact (Some approx) t.abs_err t.rel_err + *) | (Some (min, max), false), None -> let itv = Fval.F.to_float min, Fval.F.to_float max in let itv = I.of_floats Precisions.Double itv in + `Value { t with approx = Some itv } + (* `Value (apply_constraints t.exact (Some itv) t.abs_err t.rel_err) + *) | _ -> `Value t (* Associated Datatype *) @@ -730,7 +572,7 @@ module Approx_Sem = let handle_prec op x y = let px = I.precision_of x in let py = I.precision_of y in - if Precisions.equal px py then op ~prec:px x y + if Precisions.eq px py then op ~prec:px x y else assert false let add x y = match x.approx, y.approx with @@ -759,29 +601,33 @@ module Abs_Err_Sem : Semantic = let sqrt v = match Approx_Sem.sqrt v with | None -> I.top ~prec:Precisions.Real - | Some approx -> Reals.sub approx (Exact_Sem.sqrt v) + | Some approx -> Reals.sub approx (Exact_Sem.sqrt v) let add x y = (x, y) >>-% fun x_approx y_approx -> + let e_interval = Cst.eps_interval (I.precision_of x_approx) in let abs_err = Reals.add x.abs_err y.abs_err in let approx = Reals.add x_approx y_approx in - let t = Reals.mul approx Cst.e_interval in + let t = Reals.mul approx e_interval in Reals.add abs_err t let sub x y = (x, y) >>-% fun x_approx y_approx -> + let e_interval = Cst.eps_interval (I.precision_of x_approx) in let abs_err = Reals.sub x.abs_err y.abs_err in let approx = Reals.sub x_approx y_approx in - let t = Reals.mul approx Cst.e_interval in + let t = Reals.mul approx e_interval in Reals.add abs_err t let mul x y = (x, y) >>-% fun x_approx y_approx -> + let e_interval = Cst.eps_interval (I.precision_of x_approx) in let approx = Reals.mul x_approx y_approx in - let mul = Reals.mul approx Cst.e_interval in + let mul = Reals.mul approx e_interval in let xey = Reals.mul x.exact y.abs_err in let yex = Reals.mul y.exact x.abs_err in let exy = Reals.mul x.abs_err y.abs_err in Reals.add (Reals.add mul (Reals.add xey yex)) exy let div x y = (x, y) >>-% fun x_approx y_approx -> + let e_interval = Cst.eps_interval (I.precision_of x_approx) in let xy = Reals.div x.exact y.exact in let eax = x.abs_err in let eay = Reals.mul xy y.abs_err in - let phi = Reals.mul x_approx Cst.e_interval in + let phi = Reals.mul x_approx e_interval in Reals.div (Reals.add (Reals.sub eax eay) phi) y_approx end @@ -790,32 +636,71 @@ module Abs_Err_Sem : Semantic = * during the forward operations *) module Rel_Err_Sem : Semantic = struct + let get_prec v = match v.approx with + | None -> Precisions.Real + | Some t -> I.precision_of t let neg v = Reals.neg v.rel_err let sqrt v = + let std_error = Cst.std_error (get_prec v) in let s = Reals.sqrt (Reals.add v.rel_err Cst.one_interval) in - Reals.sub (Reals.mul s Cst.std_error) Cst.one_interval + Reals.sub (Reals.mul s std_error) Cst.one_interval + let add x y = - let d = F.max (Reals.max_abs x.rel_err) (Reals.max_abs y.rel_err) in - let num = R.add (Reals.max_abs x.exact) (Reals.max_abs y.exact) in - let den = Reals.min_abs (Reals.add x.exact y.exact) in - let eps = R.add (F.of_float Precisions.Real 1.0) Cst.e in - let v = R.add (R.mul (R.mul (R.div num den) d) eps) Cst.e in - I.make (R.neg v) v, Precisions.Real + let r1 = (x, y) >>-% fun x_approx y_approx -> + let e = Cst.eps_interval (get_prec x) in + let div = Reals.add x.exact y.exact in + let a = Reals.mul x.exact x.rel_err in + let b = Reals.mul y.exact y.rel_err in + let f = Reals.mul (Reals.add x_approx y_approx) e in + Reals.div (Reals.add (Reals.add a b) f) div + in + let r2 = + let e = Cst.machine_epsilon (get_prec x) in + let d = F.max (Reals.max_abs x.rel_err) (Reals.max_abs y.rel_err) in + let num = R.add (Reals.max_abs x.exact) (Reals.max_abs y.exact) in + let den = Reals.min_abs (Reals.add x.exact y.exact) in + let eps = R.add (F.of_float Precisions.Real 1.0) e in + let v = R.add (R.mul (R.mul (R.div num den) d) eps) e in + I.make (R.neg v) v, Precisions.Real + in + Format.printf "Add results :@[<v>%a@ %a@]@." I.pretty r1 I.pretty r2 ; + match I.narrow ~prec:Precisions.Real r1 r2 with + | `Bottom -> assert false + | `Value v -> v + let sub x y = - let d = F.max (Reals.max_abs x.rel_err) (Reals.max_abs y.rel_err) in - let num = R.add (Reals.max_abs x.exact) (Reals.max_abs y.exact) in - let den = Reals.min_abs (Reals.sub x.exact y.exact) in - let eps = R.add (F.of_float Precisions.Real 1.0) Cst.e in - let v = R.add (R.mul (R.mul (R.div num den) d) eps) Cst.e in - I.make (R.neg v) v, Precisions.Real + let r1 = (x, y) >>-% fun x_approx y_approx -> + let e = Cst.eps_interval (get_prec x) in + let div = Reals.sub x.exact y.exact in + let a = Reals.mul x.exact x.rel_err in + let b = Reals.mul y.exact y.rel_err in + let f = Reals.mul (Reals.sub x_approx y_approx) e in + Reals.div (Reals.add (Reals.sub a b) f) div + in + let r2 = + let e = Cst.machine_epsilon (get_prec x) in + let d = F.max (Reals.max_abs x.rel_err) (Reals.max_abs y.rel_err) in + let num = R.add (Reals.max_abs x.exact) (Reals.max_abs y.exact) in + let den = Reals.min_abs (Reals.sub x.exact y.exact) in + let eps = R.add (F.of_float Precisions.Real 1.0) e in + let v = R.add (R.mul (R.mul (R.div num den) d) eps) e in + I.make (R.neg v) v, Precisions.Real + in + Format.printf "Sub results :@[<v>%a@ %a@]@." I.pretty r1 I.pretty r2 ; + match I.narrow ~prec:Precisions.Real r1 r2 with + | `Bottom -> assert false + | `Value v -> v + let mul x y = + let std_error = Cst.std_error (get_prec x) in let ex = Reals.add Cst.one_interval x.rel_err in let ey = Reals.add Cst.one_interval y.rel_err in - Reals.sub (Reals.mul (Reals.mul ex ey) Cst.std_error) Cst.one_interval + Reals.sub (Reals.mul (Reals.mul ex ey) std_error) Cst.one_interval let div x y = + let std_error = Cst.std_error (get_prec x) in let ex = Reals.add Cst.one_interval x.rel_err in let ey = Reals.add Cst.one_interval y.rel_err in - Reals.sub (Reals.mul (Reals.div ex ey) Cst.std_error) Cst.one_interval + Reals.sub (Reals.mul (Reals.div ex ey) std_error) Cst.one_interval end @@ -832,9 +717,16 @@ let constant _ = function | Some s -> I.of_strings Precisions.Real (s, s) | None -> I.of_floats Precisions.Real (r, r) in - let top_real = I.top Precisions.Real in let approx = Some (I.of_floats prec (r, r)) in - apply_constraints exact approx top_real top_real + let eps = Cst.eps_interval Precisions.Real in + let abs_err = Reals.mul eps exact in + let rel_err = eps in + { exact ; approx ; abs_err ; rel_err } + (* + let top_real = I.top Precisions.Real in + let approx = I.of_floats prec (r, r) in + apply_constraints exact (Some approx) top_real top_real + *) | _ -> top (* Forward operations *) @@ -850,20 +742,36 @@ let forward_binop ~context:_ _ op x y = match op with | PlusA -> let exact , approx = Exact_Sem.add x y, Approx_Sem.add x y in + let abs_err, rel_err = Abs_Err_Sem.add x y, Rel_Err_Sem.add x y in + return { exact ; approx ; abs_err ; rel_err } + (* let prg_abs, prg_rel = Abs_Err_Sem.add x y, Rel_Err_Sem.add x y in return (apply_constraints exact approx prg_abs prg_rel) + *) | MinusA -> let exact , approx = Exact_Sem.sub x y, Approx_Sem.sub x y in + let abs_err, rel_err = Abs_Err_Sem.sub x y, Rel_Err_Sem.sub x y in + return { exact ; approx ; abs_err ; rel_err } + (* let prg_abs, prg_rel = Abs_Err_Sem.sub x y, Rel_Err_Sem.sub x y in return (apply_constraints exact approx prg_abs prg_rel) + *) | Mult -> let exact , approx = Exact_Sem.mul x y, Approx_Sem.mul x y in + let abs_err, rel_err = Abs_Err_Sem.mul x y, Rel_Err_Sem.mul x y in + return { exact ; approx ; abs_err ; rel_err } + (* let prg_abs, prg_rel = Abs_Err_Sem.mul x y, Rel_Err_Sem.mul x y in return (apply_constraints exact approx prg_abs prg_rel) + *) | Div -> let exact , approx = Exact_Sem.div x y, Approx_Sem.div x y in + let abs_err, rel_err = Abs_Err_Sem.div x y, Rel_Err_Sem.div x y in + return { exact ; approx ; abs_err ; rel_err } + (* let prg_abs, prg_rel = Abs_Err_Sem.div x y, Rel_Err_Sem.div x y in return (apply_constraints exact approx prg_abs prg_rel) + *) | _ -> return top (* Backward operations *) -- GitLab