From a85a6e7b7b148793ad7ab1bc3ca93ec4dc981cae Mon Sep 17 00:00:00 2001
From: Maxime Jacquemin <maxime.jacquemin@cea.fr>
Date: Wed, 24 Jan 2024 18:23:36 +0100
Subject: [PATCH] [Integer] power_int_positive_int returns an option
This function relies on Z.power_int_positive_int, which can raise the
exception Invalid_argument when dealing with an exponent too big to
produce a valid result in GMP. This exception is now catched and thus
the function produce an option.
Modifications at callsites have been taken care of, and the file
floating_point.ml has been made more human readable.
---
src/libraries/stdlib/integer.ml | 4 +-
src/libraries/stdlib/integer.mli | 2 +-
src/libraries/utils/floating_point.ml | 555 ++++++++----------
.../eva/engine/subdivided_evaluation.ml | 5 +-
4 files changed, 254 insertions(+), 312 deletions(-)
diff --git a/src/libraries/stdlib/integer.ml b/src/libraries/stdlib/integer.ml
index 47a589ee970..3c50b502862 100644
--- a/src/libraries/stdlib/integer.ml
+++ b/src/libraries/stdlib/integer.ml
@@ -38,7 +38,9 @@ let two_power n =
else
two_power_of_int k
-let power_int_positive_int = Big_int_Z.power_int_positive_int
+let power_int_positive_int n e =
+ try Some (Big_int_Z.power_int_positive_int n e)
+ with Invalid_argument _ -> None
let popcount = Z.popcount
diff --git a/src/libraries/stdlib/integer.mli b/src/libraries/stdlib/integer.mli
index fd3e19bf0d6..8647dbcbd7b 100644
--- a/src/libraries/stdlib/integer.mli
+++ b/src/libraries/stdlib/integer.mli
@@ -203,7 +203,7 @@ val two_power : t -> t
val two_power_of_int : int -> t
(** Computes [2^n] *)
-val power_int_positive_int: int -> int -> t
+val power_int_positive_int: int -> int -> t option
(** Exponentiation *)
val extract_bits : start:t -> stop:t -> t -> t
diff --git a/src/libraries/utils/floating_point.ml b/src/libraries/utils/floating_point.ml
index 736b3a549b2..cf038e2fac3 100644
--- a/src/libraries/utils/floating_point.ml
+++ b/src/libraries/utils/floating_point.ml
@@ -20,6 +20,8 @@
(* *)
(**************************************************************************)
+
+
type c_rounding_mode =
FE_ToNearest | FE_Upward | FE_Downward | FE_TowardZero
@@ -29,253 +31,235 @@ let string_of_c_rounding_mode = function
| FE_Downward -> "FE_DOWNWARD"
| FE_TowardZero -> "FE_TOWARDZERO"
+
+
+
external set_round_downward: unit -> unit = "set_round_downward" [@@noalloc]
external set_round_upward: unit -> unit = "set_round_upward" [@@noalloc]
external set_round_nearest_even: unit -> unit = "set_round_nearest_even" [@@noalloc]
external set_round_toward_zero : unit -> unit = "set_round_toward_zero" [@@noalloc]
external get_rounding_mode: unit -> c_rounding_mode = "get_rounding_mode" [@@noalloc]
external set_rounding_mode: c_rounding_mode -> unit = "set_rounding_mode" [@@noalloc]
-
external round_to_single_precision_float: float -> float = "round_to_float"
-external sys_single_precision_of_string: string -> float =
- "single_precision_of_string"
+external sys_single_precision_of_string: string -> float = "single_precision_of_string"
(* TODO two functions above: declare "float",
must have separate version for bytecode, see OCaml manual *)
let max_single_precision_float = Int32.float_of_bits 0x7f7fffffl
let most_negative_single_precision_float = -. max_single_precision_float
-type parsed_float = {
- f_nearest : float ;
- f_lower : float ;
- f_upper : float ;
-}
+
+
+type parsed_float = { f_nearest : float ; f_lower : float ; f_upper : float }
+
+let zero = { f_lower = 0.0 ; f_nearest = 0.0 ; f_upper = 0.0 }
let inf ~man_size ~max_exp =
- let biggest_not_inf = ldexp (2.0 -. ldexp 1.0 (~- man_size)) max_exp in
- {
- f_lower = biggest_not_inf ;
- f_nearest = infinity ;
- f_upper = infinity ;
- }
+ let f_lower = ldexp (2.0 -. ldexp 1.0 (~- man_size)) max_exp in
+ { f_lower ; f_nearest = infinity ; f_upper = infinity }
+
+let neg f =
+ let f_lower = -. f.f_upper in
+ let f_nearest = -. f.f_nearest in
+ let f_upper = -. f.f_lower in
+ { f_lower ; f_nearest ; f_upper }
(* [s = num * 2^exp / den] hold *)
let make_float ~num ~den ~exp ~man_size ~min_exp ~max_exp =
- assert (Integer.gt num Integer.zero);
- assert (Integer.gt den Integer.zero);
-(*
- Format.printf "make_float: num den exp:@\n%a@\n@\n%a@\n@\n%d@.min_exp:%d max_exp:%d@."
- Datatype.Integer.pretty num Datatype.Integer.pretty den exp min_exp max_exp;
-*)
- let size_bi = Integer.of_int man_size in
- let ssize_bi = Integer.of_int (succ man_size) in
- let min_exp = min_exp - man_size in
-
- let den = ref den in
- let exp = ref exp in
- while
- Integer.ge num (Integer.shift_left !den ssize_bi)
- || !exp < min_exp
- do
- den := Integer.shift_left !den Integer.one;
- incr exp
- done;
- let den = !den in
- let shifted_den = Integer.shift_left den size_bi in
- let num = ref num in
- while
- Integer.lt !num shifted_den && !exp > min_exp
- do
- num := Integer.shift_left !num Integer.one;
- decr exp
- done;
- let num = !num in
- let exp = !exp in
-(*
- Format.printf "make_float2: num den exp:@\n%a@\n@\n%a@\n@\n%d@."
- Datatype.Integer.pretty num Datatype.Integer.pretty den exp;
-*)
- if exp > max_exp - man_size then inf ~man_size ~max_exp
- else
- let man,rem = Integer.e_div_rem num den in
- let rem2 = (* twice the remainder *)
- Integer.shift_left rem Integer.one
- in
- let man = Integer.to_int64_exn man in
- (* Format.printf "pre-round: num den man rem:@\n%a@\n@\n%a@\n@\n%Ld@\n@\n%a@."
- Datatype.Integer.pretty num Datatype.Integer.pretty den
- man Datatype.Integer.pretty rem; *)
- let lowb = ldexp (Int64.to_float man) exp in
- if Integer.is_zero rem2 then {
- f_lower = lowb ;
- f_nearest = lowb ;
- f_upper = lowb ;
- } else
- let upb = ldexp (Int64.to_float (Int64.succ man)) exp in
- if Integer.lt rem2 den ||
- (Integer.equal rem2 den && (Int64.logand man Int64.one) = 0L)
- then {
- f_lower = lowb ;
- f_nearest = lowb ;
- f_upper = upb ;
- }
- else {
- f_lower = lowb ;
- f_nearest = upb ;
- f_upper = upb ;
- }
+ assert Integer.(ge num zero && gt den zero) ;
+ if not (Integer.is_zero num) then
+ let size_bi = Integer.of_int man_size in
+ let ssize_bi = Integer.of_int (succ man_size) in
+ let min_exp = min_exp - man_size in
+ let num = ref num and den = ref den and exp = ref exp in
+ while Integer.(ge !num (shift_left !den ssize_bi)) || !exp < min_exp do
+ den := Integer.shift_left !den Integer.one ;
+ incr exp
+ done ;
+ let shifted_den = Integer.shift_left !den size_bi in
+ while Integer.lt !num shifted_den && !exp > min_exp do
+ num := Integer.shift_left !num Integer.one ;
+ decr exp
+ done ;
+ if !exp <= max_exp - man_size then
+ let man, rem = Integer.e_div_rem !num !den in
+ let rem2 = Integer.shift_left rem Integer.one in
+ let man = Integer.to_int64_exn man in
+ let lowb = ldexp Int64.(to_float man) !exp in
+ let upb = ldexp Int64.(to_float (succ man)) !exp in
+ let rem2_zero = Integer.is_zero rem2 in
+ let rem2_lt_den = Integer.lt rem2 !den in
+ let rem2_is_den = Integer.equal rem2 !den in
+ let last_zero = Int64.(logand man one) = 0L in
+ let near_down = rem2_zero || rem2_lt_den || (rem2_is_den && last_zero) in
+ let f_lower = lowb in
+ let f_upper = if rem2_zero then lowb else upb in
+ let f_nearest = if near_down then lowb else upb in
+ { f_lower ; f_nearest ; f_upper }
+ else inf ~man_size ~max_exp
+ else zero
+
+
let reg_exp = "[eE][+]?\\(-?[0-9]+\\)"
let reg_dot = "[.]"
let reg_numopt = "\\([0-9]*\\)"
let reg_num = "\\([0-9]+\\)"
-let numdotfrac = Str.regexp (reg_numopt ^ reg_dot ^ reg_numopt)
-let numdotfracexp = Str.regexp (reg_numopt ^ reg_dot ^ reg_numopt ^ reg_exp)
-let numexp = Str.regexp (reg_num ^ reg_exp)
+let num_dot_frac = Str.regexp (reg_numopt ^ reg_dot ^ reg_numopt)
+let num_dot_frac_exp = Str.regexp (reg_numopt ^ reg_dot ^ reg_numopt ^ reg_exp)
+let num_exp = Str.regexp (reg_num ^ reg_exp)
+
+
+
+let float_of_string_opt s =
+ try Some (float_of_string s)
+ with Failure _ -> None
+
+let sys_single_precision_of_string_opt s =
+ try Some (sys_single_precision_of_string s)
+ with Failure _ -> None
+
+let is_hexadecimal s =
+ String.length s >= 2 && s.[0] = '0' && (Char.uppercase_ascii s.[1] = 'X')
+
+
exception Shortcut of parsed_float
-let zero = { f_lower = 0.0 ; f_nearest = 0.0 ; f_upper = 0.0 }
+let match_exp ~man_size ~min_exp ~max_exp group s =
+ let s = Str.matched_group group s in
+ match int_of_string_opt s with
+ | Some n -> n
+ | None when s.[0] = '-' ->
+ let f_upper = ldexp 1.0 (min_exp - man_size) in
+ raise (Shortcut { f_lower = 0.0 ; f_nearest = 0.0 ; f_upper })
+ | None -> raise (Shortcut (inf ~man_size ~max_exp))
(* [man_size] is the size of the mantissa, [min_exp] the frontier exponent
between normalized and denormalized numbers *)
-let parse_float ~man_size ~min_exp ~max_exp s =
- (* Format.printf "parse: %s@." s; *)
- let match_exp group =
- let s = Str.matched_group group s in
- try
- int_of_string s
- with Failure _ ->
- (* Format.printf "Error in exponent: %s@." s; *)
- if s.[0] = '-'
- then raise (Shortcut {
- f_lower = 0.0 ;
- f_nearest = 0.0 ;
- f_upper = ldexp 1.0 (min_exp - man_size) ;
- })
- else raise (Shortcut (inf ~man_size ~max_exp))
- in
- try
- (* At the end of the function, [s = num * 2^exp / den] *)
- let num, den, exp =
- if Str.string_match numdotfracexp s 0
- then
- let n = Str.matched_group 1 s in
- let frac = Str.matched_group 2 s in
- let len_frac = String.length frac in
- let num = Integer.of_string (n ^ frac) in
- let den = Integer.power_int_positive_int 5 len_frac in
- if Integer.is_zero num then raise (Shortcut zero);
- let exp10 = match_exp 3
- in
- if exp10 >= 0
- then
- Integer.mul num (Integer.power_int_positive_int 5 exp10),
- den,
- exp10 - len_frac
- else
- num,
- Integer.mul den (Integer.power_int_positive_int 5 (~- exp10)),
- exp10 - len_frac
- else if Str.string_match numdotfrac s 0
- then
- let n = Str.matched_group 1 s in
- let frac = Str.matched_group 2 s in
- let len_frac = String.length frac in
- Integer.of_string (n ^ frac),
- Integer.power_int_positive_int 5 len_frac,
- ~- len_frac
- else if Str.string_match numexp s 0
- then
- let n = Str.matched_group 1 s in
- let num = Integer.of_string n in
- if Integer.is_zero num then raise (Shortcut zero);
- let exp10 = match_exp 2 in
- if exp10 >= 0
- then
- Integer.mul num (Integer.power_int_positive_int 5 exp10),
- Integer.one,
- exp10
- else
- num,
- (Integer.power_int_positive_int 5 (~- exp10)),
- exp10
- else (Format.printf "Could not parse floating point number %S@." s;
- assert false)
- in
- if Integer.is_zero num
- then zero
- else
- make_float ~num ~den ~exp ~man_size ~min_exp ~max_exp
- with Shortcut r -> r
-
-let is_hex s =
- let l = String.length s in
- l >= 2 && s.[0] = '0' && (s.[1] = 'x' || s.[1] = 'X')
-
-let opp_parse_float f =
- { f_lower = -. f.f_upper ; f_nearest = -. f.f_nearest ; f_upper = -. f.f_lower }
+let parse_positive_float_with_shortcut ~man_size ~min_exp ~max_exp s =
+ let open Option.Operators in
+ let match_exp group = match_exp ~man_size ~min_exp ~max_exp group s in
+ let return = make_float ~man_size ~min_exp ~max_exp in
+ (* At the end of the function, [s = num * 2^exp / den] *)
+ if Str.string_match num_dot_frac_exp s 0 then
+ let n = Str.matched_group 1 s in
+ let frac = Str.matched_group 2 s in
+ let len_frac = String.length frac in
+ let num = Integer.of_string (n ^ frac) in
+ let* den = Integer.power_int_positive_int 5 len_frac in
+ if Integer.is_zero num then raise (Shortcut zero) ;
+ let exp10 = match_exp 3 in
+ let+ pow5 = Integer.power_int_positive_int 5 (abs exp10) in
+ let num = if exp10 >= 0 then Integer.mul num pow5 else num in
+ let den = if exp10 >= 0 then den else Integer.mul den pow5 in
+ let exp = exp10 - len_frac in
+ return ~num ~den ~exp
+ else if Str.string_match num_dot_frac s 0 then
+ let n = Str.matched_group 1 s in
+ let frac = Str.matched_group 2 s in
+ let len_frac = String.length frac in
+ let num = Integer.of_string (n ^ frac) in
+ let+ den = Integer.power_int_positive_int 5 len_frac in
+ return ~num ~den ~exp:(~- len_frac)
+ else if Str.string_match num_exp s 0 then
+ let n = Str.matched_group 1 s in
+ let num = Integer.of_string n in
+ if Integer.is_zero num then raise (Shortcut zero) ;
+ let exp10 = match_exp 2 in
+ let+ pow5 = Integer.power_int_positive_int 5 (abs exp10) in
+ let num = if exp10 >= 0 then Integer.mul num pow5 else num in
+ let den = if exp10 >= 0 then Integer.one else pow5 in
+ return ~num ~den ~exp:exp10
+ else None
+
+let parse_positive_float ~man_size ~min_exp ~max_exp s =
+ try parse_positive_float_with_shortcut ~man_size ~min_exp ~max_exp s
+ with Shortcut r -> Some r
+
+
let rec single_precision_of_string s =
+ let open Option.Operators in
if s.[0] = '-' then
- opp_parse_float (single_precision_of_string (String.sub s 1 (String.length s - 1)))
- else if is_hex s
- then
- try
- let f = sys_single_precision_of_string s in
- { f_lower = f ; f_nearest = f ; f_upper = f }
- with Failure _ ->
- Kernel.fatal "could not parse single-precision float string: %s" s
- else (* decimal *)
- parse_float ~man_size:23 ~min_exp:(-126) ~max_exp:127 s
-
-(* May raise Failure("float_of_string"). *)
+ let s = String.(sub s 1 (length s - 1)) in
+ let+ f = single_precision_of_string s in
+ neg f
+ else if is_hexadecimal s then
+ let+ f = sys_single_precision_of_string_opt s in
+ { f_lower = f ; f_nearest = f ; f_upper = f }
+ else parse_positive_float ~man_size:23 ~min_exp:(-126) ~max_exp:127 s
+
let rec double_precision_of_string s =
+ let open Option.Operators in
if s.[0] = '-' then
- opp_parse_float (double_precision_of_string (String.sub s 1 (String.length s - 1)))
- else if is_hex s
- then
- let f = float_of_string s in
+ let s = String.(sub s 1 (length s - 1)) in
+ let+ f = double_precision_of_string s in
+ neg f
+ else if is_hexadecimal s then
+ let+ f = float_of_string_opt s in
{ f_lower = f ; f_nearest = f ; f_upper = f }
- else (* decimal *)
- parse_float ~man_size:52 ~min_exp:(-1022) ~max_exp:1023 s
+ else parse_positive_float ~man_size:52 ~min_exp:(-1022) ~max_exp:1023 s
-let parse_kind kind string =
- match kind with
+let parse_fkind string = function
| Cil_types.FFloat -> single_precision_of_string string
- | Cil_types.FDouble
- | Cil_types.FLongDouble -> double_precision_of_string string
+ | Cil_types.(FDouble | FLongDouble) -> double_precision_of_string string
+
+
+
+let fkind_of_char = function
+ | 'F' -> Cil_types.FFloat, true
+ | 'D' -> Cil_types.FDouble, true
+ | 'L' -> Cil_types.FLongDouble, true
+ | _ -> Cil_types.FDouble, false
+
+let suffix_of_fkind = function
+ | Cil_types.FFloat -> 'F'
+ | Cil_types.FDouble -> 'D'
+ | Cil_types.FLongDouble -> 'L'
+
+let pretty_fkind fmt = function
+ | Cil_types.FFloat -> Format.fprintf fmt "single precision"
+ | Cil_types.FDouble -> Format.fprintf fmt "double precision"
+ | Cil_types.FLongDouble -> Format.fprintf fmt "long double precision"
+
+
+
+let cannot_be_parsed string fkind =
+ Kernel.abort ~current:true
+ "The string %s cannot be parsed as a %a floating-point constant"
+ string pretty_fkind fkind
+
+let empty_string () =
+ Kernel.abort ~current:true
+ "Parsing an empty string as a floating-point constant."
let parse string =
let l = String.length string - 1 in
- if l < 0
- then Kernel.fatal ~current:true
- "Parsing an empty string as a floating-point constant."
- else
+ if l >= 0 then
let last = Char.uppercase_ascii string.[l] in
- let suffix, kind =
- match last with
- | 'F' -> true, Cil_types.FFloat
- | 'D' -> true, Cil_types.FDouble
- | 'L' -> true, Cil_types.FLongDouble
- | _ -> false, Cil_types.FDouble
- in
+ let fkind, suffix = fkind_of_char last in
let baseint = if suffix then String.sub string 0 l else string in
- try
- let basefloat = parse_kind kind baseint in
- kind, basefloat
- with Failure _ -> (* should never happen, suffix already stripped *)
- Kernel.fatal ~current:true
- "Unexpected error parsing floating-point constant: %s." string
+ match parse_fkind baseint fkind with
+ | Some result -> fkind, result
+ | None -> cannot_be_parsed string fkind
+ else empty_string ()
+
+
let has_suffix fk lit =
- let s = match fk with
- | Cil_types.FFloat -> 'F'
- | Cil_types.FDouble -> 'D'
- | Cil_types.FLongDouble -> 'L' in
let ln = String.length lit in
- ln > 0 && Char.uppercase_ascii lit.[ln-1] = s
+ let suffix = suffix_of_fkind fk in
+ ln > 0 && Char.uppercase_ascii lit.[ln - 1] = suffix
+
+let is_not_finite f =
+ match classify_float f with
+ | FP_normal | FP_subnormal | FP_zero -> false
+ | FP_infinite | FP_nan -> true
+
+let is_not_integer s =
+ String.(contains s '.' || contains s 'e' || contains s 'E')
let pretty_normal ~use_hex fmt f =
let double_norm = Int64.shift_left 1L 52 in
@@ -286,127 +270,79 @@ let pretty_normal ~use_hex fmt f =
let exp = Int64.to_int (Int64.shift_right_logical i 52) in
let man = Int64.logand i double_mask in
let s = if s then "-" else "" in
- if exp = 2047
- then begin
- if man = 0L
- then
- Format.fprintf fmt "%sinf" s
- else
- Format.fprintf fmt "NaN"
- end
+ if exp = 2047 then
+ Format.(if man = 0L then fprintf fmt "%sinf" s else fprintf fmt "NaN")
else
- let firstdigit, exp =
- if exp <> 0
- then 1, (exp - 1023)
- else 0, if f = 0. then 0 else -1022
- in
- if not use_hex
- then begin
- let firstdigit, man, exp =
- if 0 < exp && exp <= 12
- then begin
- Int64.to_int
- (Int64.shift_right_logical
- (Int64.logor man double_norm)
- (52 - exp)),
- Int64.logand (Int64.shift_left man exp) double_mask,
- 0
- end
- else firstdigit, man, exp
- in
- let d =
- Int64.float_of_bits
- (Int64.logor 0x3ff0000000000000L man)
- in
- let d, re =
- if d >= 1.5
- then d -. 1.5, 5000000000000000L
- else d -. 1.0, 0L
- in
- let d = d *. 1e16 in
+ let firstdigit = if exp <> 0 then 1 else 0 in
+ let exp = if exp <> 0 then exp - 1023 else if f = 0. then 0 else -1022 in
+ if not use_hex then
+ let in_bound = 0 < exp && exp <= 12 in
+ let doubled_man = Int64.logor man double_norm in
+ let shifted = Int64.shift_right_logical doubled_man (52 - exp) in
+ let firstdigit = if in_bound then Int64.to_int shifted else firstdigit in
+ let doubled = Int64.(logand (shift_left man exp) double_mask) in
+ let man = if in_bound then doubled else man in
+ let exp = if in_bound then 0 else exp in
+ let d = Int64.(float_of_bits (logor 0x3ff0000000000000L man)) in
+ let re = if d >= 1.5 then 5000000000000000L else 0L in
+ let shift = if d >= 1.5 then 1.5 else 1.0 in
+ let d = (d -. shift) *. 1e16 in
let decdigits = Int64.add re (Int64.of_float d) in
if exp = 0 || (firstdigit = 0 && decdigits = 0L && exp = -1022)
- then
- Format.fprintf fmt "%s%d.%016Ld"
- s
- firstdigit
- decdigits
- else
- Format.fprintf fmt "%s%d.%016Ld*2^%d"
- s
- firstdigit
- decdigits
- exp
- end
- else
- Format.fprintf fmt "%s0x%d.%013Lxp%d"
- s
- firstdigit
- man
- exp
+ then Format.fprintf fmt "%s%d.%016Ld" s firstdigit decdigits
+ else Format.fprintf fmt "%s%d.%016Ld*2^%d" s firstdigit decdigits exp
+ else Format.fprintf fmt "%s0x%d.%013Lxp%d" s firstdigit man exp
+
+let ensure_round_nearest_even () =
+ if get_rounding_mode () <> FE_ToNearest then
+ let mode = string_of_c_rounding_mode (get_rounding_mode ()) in
+ let () = Kernel.failure "pretty: rounding mode (%s) <> FE_TONEAREST" mode in
+ set_round_nearest_even ()
let pretty fmt f =
- let use_hex = Kernel.FloatHex.get() in
+ let use_hex = Kernel.FloatHex.get () in
(* should always arrive here with nearest_even *)
- if get_rounding_mode () <> FE_ToNearest then begin
- Kernel.failure "pretty: rounding mode (%s) <> FE_TONEAREST"
- (string_of_c_rounding_mode (get_rounding_mode ()));
- set_round_nearest_even();
- end;
- if use_hex || (Kernel.FloatNormal.get ())
- then
- pretty_normal ~use_hex fmt f
- else begin
+ ensure_round_nearest_even () ;
+ if not (use_hex || Kernel.FloatNormal.get ()) then
let r = Format.sprintf "%.*g" 12 f in
- if (String.contains r '.' || String.contains r 'e' ||
- String.contains r 'E')
- || (match classify_float f with
- | FP_normal | FP_subnormal | FP_zero -> false
- | FP_infinite | FP_nan -> true)
- then Format.pp_print_string fmt r
- else Format.fprintf fmt "%s." r
- end
+ let dot = if is_not_integer r || is_not_finite f then "" else "." in
+ Format.fprintf fmt "%s%s" r dot
+ else pretty_normal ~use_hex fmt f
+
type sign = Neg | Pos
exception Float_Non_representable_as_Int64 of sign
+let min_64_float = -9.22337203685477581e+18
+let max_64_float = 9.22337203685477478e+18
+
(* If the argument [x] is not in the range [min_64_float, 2*max_64_float],
raise Float_Non_representable_as_Int64. This is the most reasonable as
a floating-point number may represent an exponentially large integer. *)
-let truncate_to_integer =
- let min_64_float = -9.22337203685477581e+18
- (* Int64.to_float (-0x8000000000000000L) *)
- in
- let max_64_float = 9.22337203685477478e+18
- (* let open Int64 in
- float_of_bits (pred (bits_of_float (to_float max_int))) *)
- in
- fun x ->
- let max_64_float = Fun.id max_64_float in
- if x < min_64_float
- then raise (Float_Non_representable_as_Int64 Neg);
- if x > (max_64_float +. max_64_float)
- then raise (Float_Non_representable_as_Int64 Pos);
- if x <= max_64_float then
- Integer.of_int64 (Int64.of_float x)
- else
- Integer.add
- (Integer.of_int64 (Int64.of_float (x +. min_64_float)))
- (Integer.two_power_of_int 63)
+let truncate_to_integer x =
+ if x < min_64_float then raise (Float_Non_representable_as_Int64 Neg) ;
+ if x > 2. *. max_64_float then raise (Float_Non_representable_as_Int64 Pos) ;
+ let convert x = Integer.of_int64 (Int64.of_float x) in
+ let shift n = Integer.(add (two_power_of_int 63)) n in
+ if x <= max_64_float then convert x else convert (x +. min_64_float) |> shift
let bits_of_max_double =
Integer.of_int64 (Int64.bits_of_float max_float)
+
let bits_of_most_negative_double =
Integer.of_int64 (Int64.bits_of_float (-. max_float))
(** See e.g. http://www.h-schmidt.net/FloatConverter/IEEE754.html *)
let bits_of_max_float = Integer.of_int64 0x7F7FFFFFL
let bits_of_most_negative_float =
- let v = Int64.of_int32 0xFF7FFFFFl in(* cast to int32 to get negative value *)
+ (* cast to int32 to get negative value *)
+ let v = Int64.of_int32 0xFF7FFFFFl in
Integer.of_int64 v
+
+
external fround: float -> float = "c_round"
external trunc: float -> float = "c_trunc"
@@ -427,6 +363,7 @@ external atanf: float -> float = "c_atanf"
external atan2f: float -> float -> float = "c_atan2f"
+
(** C math-like functions *)
let isnan f =
@@ -447,27 +384,31 @@ let neg_min_single_precision_denormal = -. min_single_precision_denormal
(* auxiliary functions for nextafter/nextafterf *)
let min_denormal_float ~is_f32 =
if is_f32 then min_single_precision_denormal else min_denormal
+
let nextafter_aux ~is_f32 fincr fdecr x y =
- if x = y (* includes cases "(0.0, -0.0) => -0.0" and its symmetric *)
- then y
+ let sign = if x < y then 1.0 else -.1.0 in
+ if x = y then y
else if isnan x || isnan y then nan
- else if x = 0.0 (* or -0.0 *) then
- if x < y then min_denormal_float ~is_f32 else -. (min_denormal_float ~is_f32)
- (* the following conditions might be simpler if we had unsigned ints
- (uint32/uint64) *)
- else if x = neg_infinity (* && y = neg_infinity *) then fdecr x
- else if (x < y && x > 0.0) || (x > y && x < 0.0) then fincr x else fdecr x
+ else if x = 0.0 then sign *. min_denormal_float ~is_f32
+ else if x = neg_infinity then fdecr x
+ else if (x < y && x > 0.0) || (x > y && x < 0.0) then fincr x
+ else fdecr x
+
let incr_f64 f =
- Int64.float_of_bits (Int64.succ (Int64.bits_of_float f))
+ if f = neg_infinity then -. max_float
+ else Int64.(float_of_bits (succ (bits_of_float f)))
+
let decr_f64 f =
if f = infinity then max_float
- else Int64.float_of_bits (Int64.pred (Int64.bits_of_float f))
+ else Int64.(float_of_bits (pred (bits_of_float f)))
+
let incr_f32 f =
if f = neg_infinity then most_negative_single_precision_float
- else Int32.float_of_bits (Int32.succ (Int32.bits_of_float f))
+ else Int32.(float_of_bits (succ (bits_of_float f)))
+
let decr_f32 f =
if f = infinity then max_single_precision_float
- else Int32.float_of_bits (Int32.pred (Int32.bits_of_float f))
+ else Int32.(float_of_bits (pred (bits_of_float f)))
let nextafter x y =
nextafter_aux ~is_f32:false incr_f64 decr_f64 x y
diff --git a/src/plugins/eva/engine/subdivided_evaluation.ml b/src/plugins/eva/engine/subdivided_evaluation.ml
index 60949a0ef0a..0e0c9d130e8 100644
--- a/src/plugins/eva/engine/subdivided_evaluation.ml
+++ b/src/plugins/eva/engine/subdivided_evaluation.ml
@@ -767,10 +767,9 @@ module Make
subvalues that are all evaluated. Limits the number of splits to
keep the number of evaluations linear on [nb]. *)
let subdivnb =
+ let pow n e = Integer.power_int_positive_int n e |> Option.get in
if nb > 3
- then
- let pow = Integer.power_int_positive_int in
- (subdivnb * nb) / (Integer.to_int_exn (pow 2 (nb - 1)))
+ then (subdivnb * nb) / (Integer.to_int_exn (pow 2 (nb - 1)))
else subdivnb
in
Self.result ~current:true ~once:true ~dkey
--
GitLab