Skip to content
Snippets Groups Projects
Commit a85a6e7b authored by Maxime Jacquemin's avatar Maxime Jacquemin
Browse files

[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.
parent 6ee3fd7a
No related branches found
No related tags found
No related merge requests found
...@@ -38,7 +38,9 @@ let two_power n = ...@@ -38,7 +38,9 @@ let two_power n =
else else
two_power_of_int k 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 let popcount = Z.popcount
......
...@@ -203,7 +203,7 @@ val two_power : t -> t ...@@ -203,7 +203,7 @@ val two_power : t -> t
val two_power_of_int : int -> t val two_power_of_int : int -> t
(** Computes [2^n] *) (** Computes [2^n] *)
val power_int_positive_int: int -> int -> t val power_int_positive_int: int -> int -> t option
(** Exponentiation *) (** Exponentiation *)
val extract_bits : start:t -> stop:t -> t -> t val extract_bits : start:t -> stop:t -> t -> t
......
...@@ -20,6 +20,8 @@ ...@@ -20,6 +20,8 @@
(* *) (* *)
(**************************************************************************) (**************************************************************************)
type c_rounding_mode = type c_rounding_mode =
FE_ToNearest | FE_Upward | FE_Downward | FE_TowardZero FE_ToNearest | FE_Upward | FE_Downward | FE_TowardZero
...@@ -29,253 +31,235 @@ let string_of_c_rounding_mode = function ...@@ -29,253 +31,235 @@ let string_of_c_rounding_mode = function
| FE_Downward -> "FE_DOWNWARD" | FE_Downward -> "FE_DOWNWARD"
| FE_TowardZero -> "FE_TOWARDZERO" | FE_TowardZero -> "FE_TOWARDZERO"
external set_round_downward: unit -> unit = "set_round_downward" [@@noalloc] external set_round_downward: unit -> unit = "set_round_downward" [@@noalloc]
external set_round_upward: unit -> unit = "set_round_upward" [@@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_nearest_even: unit -> unit = "set_round_nearest_even" [@@noalloc]
external set_round_toward_zero : unit -> unit = "set_round_toward_zero" [@@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 get_rounding_mode: unit -> c_rounding_mode = "get_rounding_mode" [@@noalloc]
external set_rounding_mode: c_rounding_mode -> unit = "set_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 round_to_single_precision_float: float -> float = "round_to_float"
external sys_single_precision_of_string: string -> float = external sys_single_precision_of_string: string -> float = "single_precision_of_string"
"single_precision_of_string"
(* TODO two functions above: declare "float", (* TODO two functions above: declare "float",
must have separate version for bytecode, see OCaml manual *) must have separate version for bytecode, see OCaml manual *)
let max_single_precision_float = Int32.float_of_bits 0x7f7fffffl let max_single_precision_float = Int32.float_of_bits 0x7f7fffffl
let most_negative_single_precision_float = -. max_single_precision_float let most_negative_single_precision_float = -. max_single_precision_float
type parsed_float = {
f_nearest : float ;
f_lower : float ; type parsed_float = { f_nearest : float ; f_lower : float ; f_upper : 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 inf ~man_size ~max_exp =
let biggest_not_inf = ldexp (2.0 -. ldexp 1.0 (~- man_size)) max_exp in let f_lower = ldexp (2.0 -. ldexp 1.0 (~- man_size)) max_exp in
{ { f_lower ; f_nearest = infinity ; f_upper = infinity }
f_lower = biggest_not_inf ;
f_nearest = infinity ; let neg f =
f_upper = infinity ; 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 *) (* [s = num * 2^exp / den] hold *)
let make_float ~num ~den ~exp ~man_size ~min_exp ~max_exp = let make_float ~num ~den ~exp ~man_size ~min_exp ~max_exp =
assert (Integer.gt num Integer.zero); assert Integer.(ge num zero && gt den zero) ;
assert (Integer.gt den Integer.zero); if not (Integer.is_zero num) then
(* let size_bi = Integer.of_int man_size in
Format.printf "make_float: num den exp:@\n%a@\n@\n%a@\n@\n%d@.min_exp:%d max_exp:%d@." let ssize_bi = Integer.of_int (succ man_size) in
Datatype.Integer.pretty num Datatype.Integer.pretty den exp min_exp max_exp; let min_exp = min_exp - man_size in
*) let num = ref num and den = ref den and exp = ref exp in
let size_bi = Integer.of_int man_size in while Integer.(ge !num (shift_left !den ssize_bi)) || !exp < min_exp do
let ssize_bi = Integer.of_int (succ man_size) in den := Integer.shift_left !den Integer.one ;
let min_exp = min_exp - man_size in incr exp
done ;
let den = ref den in let shifted_den = Integer.shift_left !den size_bi in
let exp = ref exp in while Integer.lt !num shifted_den && !exp > min_exp do
while num := Integer.shift_left !num Integer.one ;
Integer.ge num (Integer.shift_left !den ssize_bi) decr exp
|| !exp < min_exp done ;
do if !exp <= max_exp - man_size then
den := Integer.shift_left !den Integer.one; let man, rem = Integer.e_div_rem !num !den in
incr exp let rem2 = Integer.shift_left rem Integer.one in
done; let man = Integer.to_int64_exn man in
let den = !den in let lowb = ldexp Int64.(to_float man) !exp in
let shifted_den = Integer.shift_left den size_bi in let upb = ldexp Int64.(to_float (succ man)) !exp in
let num = ref num in let rem2_zero = Integer.is_zero rem2 in
while let rem2_lt_den = Integer.lt rem2 !den in
Integer.lt !num shifted_den && !exp > min_exp let rem2_is_den = Integer.equal rem2 !den in
do let last_zero = Int64.(logand man one) = 0L in
num := Integer.shift_left !num Integer.one; let near_down = rem2_zero || rem2_lt_den || (rem2_is_den && last_zero) in
decr exp let f_lower = lowb in
done; let f_upper = if rem2_zero then lowb else upb in
let num = !num in let f_nearest = if near_down then lowb else upb in
let exp = !exp in { f_lower ; f_nearest ; f_upper }
(* else inf ~man_size ~max_exp
Format.printf "make_float2: num den exp:@\n%a@\n@\n%a@\n@\n%d@." else zero
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 ;
}
let reg_exp = "[eE][+]?\\(-?[0-9]+\\)" let reg_exp = "[eE][+]?\\(-?[0-9]+\\)"
let reg_dot = "[.]" let reg_dot = "[.]"
let reg_numopt = "\\([0-9]*\\)" let reg_numopt = "\\([0-9]*\\)"
let reg_num = "\\([0-9]+\\)" let reg_num = "\\([0-9]+\\)"
let numdotfrac = Str.regexp (reg_numopt ^ reg_dot ^ reg_numopt) let num_dot_frac = Str.regexp (reg_numopt ^ reg_dot ^ reg_numopt)
let numdotfracexp = Str.regexp (reg_numopt ^ reg_dot ^ reg_numopt ^ reg_exp) let num_dot_frac_exp = Str.regexp (reg_numopt ^ reg_dot ^ reg_numopt ^ reg_exp)
let numexp = Str.regexp (reg_num ^ 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 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 (* [man_size] is the size of the mantissa, [min_exp] the frontier exponent
between normalized and denormalized numbers *) between normalized and denormalized numbers *)
let parse_float ~man_size ~min_exp ~max_exp s = let parse_positive_float_with_shortcut ~man_size ~min_exp ~max_exp s =
(* Format.printf "parse: %s@." s; *) let open Option.Operators in
let match_exp group = let match_exp group = match_exp ~man_size ~min_exp ~max_exp group s in
let s = Str.matched_group group s in let return = make_float ~man_size ~min_exp ~max_exp in
try (* At the end of the function, [s = num * 2^exp / den] *)
int_of_string s if Str.string_match num_dot_frac_exp s 0 then
with Failure _ -> let n = Str.matched_group 1 s in
(* Format.printf "Error in exponent: %s@." s; *) let frac = Str.matched_group 2 s in
if s.[0] = '-' let len_frac = String.length frac in
then raise (Shortcut { let num = Integer.of_string (n ^ frac) in
f_lower = 0.0 ; let* den = Integer.power_int_positive_int 5 len_frac in
f_nearest = 0.0 ; if Integer.is_zero num then raise (Shortcut zero) ;
f_upper = ldexp 1.0 (min_exp - man_size) ; let exp10 = match_exp 3 in
}) let+ pow5 = Integer.power_int_positive_int 5 (abs exp10) in
else raise (Shortcut (inf ~man_size ~max_exp)) let num = if exp10 >= 0 then Integer.mul num pow5 else num in
in let den = if exp10 >= 0 then den else Integer.mul den pow5 in
try let exp = exp10 - len_frac in
(* At the end of the function, [s = num * 2^exp / den] *) return ~num ~den ~exp
let num, den, exp = else if Str.string_match num_dot_frac s 0 then
if Str.string_match numdotfracexp s 0 let n = Str.matched_group 1 s in
then let frac = Str.matched_group 2 s in
let n = Str.matched_group 1 s in let len_frac = String.length frac in
let frac = Str.matched_group 2 s in let num = Integer.of_string (n ^ frac) in
let len_frac = String.length frac in let+ den = Integer.power_int_positive_int 5 len_frac in
let num = Integer.of_string (n ^ frac) in return ~num ~den ~exp:(~- len_frac)
let den = Integer.power_int_positive_int 5 len_frac in else if Str.string_match num_exp s 0 then
if Integer.is_zero num then raise (Shortcut zero); let n = Str.matched_group 1 s in
let exp10 = match_exp 3 let num = Integer.of_string n in
in if Integer.is_zero num then raise (Shortcut zero) ;
if exp10 >= 0 let exp10 = match_exp 2 in
then let+ pow5 = Integer.power_int_positive_int 5 (abs exp10) in
Integer.mul num (Integer.power_int_positive_int 5 exp10), let num = if exp10 >= 0 then Integer.mul num pow5 else num in
den, let den = if exp10 >= 0 then Integer.one else pow5 in
exp10 - len_frac return ~num ~den ~exp:exp10
else else None
num,
Integer.mul den (Integer.power_int_positive_int 5 (~- exp10)), let parse_positive_float ~man_size ~min_exp ~max_exp s =
exp10 - len_frac try parse_positive_float_with_shortcut ~man_size ~min_exp ~max_exp s
else if Str.string_match numdotfrac s 0 with Shortcut r -> Some r
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 rec single_precision_of_string s = let rec single_precision_of_string s =
let open Option.Operators in
if s.[0] = '-' then if s.[0] = '-' then
opp_parse_float (single_precision_of_string (String.sub s 1 (String.length s - 1))) let s = String.(sub s 1 (length s - 1)) in
else if is_hex s let+ f = single_precision_of_string s in
then neg f
try else if is_hexadecimal s then
let f = sys_single_precision_of_string s in let+ f = sys_single_precision_of_string_opt s in
{ f_lower = f ; f_nearest = f ; f_upper = f } { f_lower = f ; f_nearest = f ; f_upper = f }
with Failure _ -> else parse_positive_float ~man_size:23 ~min_exp:(-126) ~max_exp:127 s
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 rec double_precision_of_string s = let rec double_precision_of_string s =
let open Option.Operators in
if s.[0] = '-' then if s.[0] = '-' then
opp_parse_float (double_precision_of_string (String.sub s 1 (String.length s - 1))) let s = String.(sub s 1 (length s - 1)) in
else if is_hex s let+ f = double_precision_of_string s in
then neg f
let f = float_of_string s in else if is_hexadecimal s then
let+ f = float_of_string_opt s in
{ f_lower = f ; f_nearest = f ; f_upper = f } { f_lower = f ; f_nearest = f ; f_upper = f }
else (* decimal *) else parse_positive_float ~man_size:52 ~min_exp:(-1022) ~max_exp:1023 s
parse_float ~man_size:52 ~min_exp:(-1022) ~max_exp:1023 s
let parse_kind kind string = let parse_fkind string = function
match kind with
| Cil_types.FFloat -> single_precision_of_string string | Cil_types.FFloat -> single_precision_of_string string
| Cil_types.FDouble | Cil_types.(FDouble | FLongDouble) -> double_precision_of_string string
| Cil_types.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 parse string =
let l = String.length string - 1 in let l = String.length string - 1 in
if l < 0 if l >= 0 then
then Kernel.fatal ~current:true
"Parsing an empty string as a floating-point constant."
else
let last = Char.uppercase_ascii string.[l] in let last = Char.uppercase_ascii string.[l] in
let suffix, kind = let fkind, suffix = fkind_of_char last in
match last with
| 'F' -> true, Cil_types.FFloat
| 'D' -> true, Cil_types.FDouble
| 'L' -> true, Cil_types.FLongDouble
| _ -> false, Cil_types.FDouble
in
let baseint = if suffix then String.sub string 0 l else string in let baseint = if suffix then String.sub string 0 l else string in
try match parse_fkind baseint fkind with
let basefloat = parse_kind kind baseint in | Some result -> fkind, result
kind, basefloat | None -> cannot_be_parsed string fkind
with Failure _ -> (* should never happen, suffix already stripped *) else empty_string ()
Kernel.fatal ~current:true
"Unexpected error parsing floating-point constant: %s." string
let has_suffix fk lit = 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 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 pretty_normal ~use_hex fmt f =
let double_norm = Int64.shift_left 1L 52 in let double_norm = Int64.shift_left 1L 52 in
...@@ -286,127 +270,79 @@ let pretty_normal ~use_hex fmt f = ...@@ -286,127 +270,79 @@ let pretty_normal ~use_hex fmt f =
let exp = Int64.to_int (Int64.shift_right_logical i 52) in let exp = Int64.to_int (Int64.shift_right_logical i 52) in
let man = Int64.logand i double_mask in let man = Int64.logand i double_mask in
let s = if s then "-" else "" in let s = if s then "-" else "" in
if exp = 2047 if exp = 2047 then
then begin Format.(if man = 0L then fprintf fmt "%sinf" s else fprintf fmt "NaN")
if man = 0L
then
Format.fprintf fmt "%sinf" s
else
Format.fprintf fmt "NaN"
end
else else
let firstdigit, exp = let firstdigit = if exp <> 0 then 1 else 0 in
if exp <> 0 let exp = if exp <> 0 then exp - 1023 else if f = 0. then 0 else -1022 in
then 1, (exp - 1023) if not use_hex then
else 0, if f = 0. then 0 else -1022 let in_bound = 0 < exp && exp <= 12 in
in let doubled_man = Int64.logor man double_norm in
if not use_hex let shifted = Int64.shift_right_logical doubled_man (52 - exp) in
then begin let firstdigit = if in_bound then Int64.to_int shifted else firstdigit in
let firstdigit, man, exp = let doubled = Int64.(logand (shift_left man exp) double_mask) in
if 0 < exp && exp <= 12 let man = if in_bound then doubled else man in
then begin let exp = if in_bound then 0 else exp in
Int64.to_int let d = Int64.(float_of_bits (logor 0x3ff0000000000000L man)) in
(Int64.shift_right_logical let re = if d >= 1.5 then 5000000000000000L else 0L in
(Int64.logor man double_norm) let shift = if d >= 1.5 then 1.5 else 1.0 in
(52 - exp)), let d = (d -. shift) *. 1e16 in
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 decdigits = Int64.add re (Int64.of_float d) in let decdigits = Int64.add re (Int64.of_float d) in
if exp = 0 || (firstdigit = 0 && decdigits = 0L && exp = -1022) if exp = 0 || (firstdigit = 0 && decdigits = 0L && exp = -1022)
then then Format.fprintf fmt "%s%d.%016Ld" s firstdigit decdigits
Format.fprintf fmt "%s%d.%016Ld" else Format.fprintf fmt "%s%d.%016Ld*2^%d" s firstdigit decdigits exp
s else Format.fprintf fmt "%s0x%d.%013Lxp%d" s firstdigit man exp
firstdigit
decdigits let ensure_round_nearest_even () =
else if get_rounding_mode () <> FE_ToNearest then
Format.fprintf fmt "%s%d.%016Ld*2^%d" let mode = string_of_c_rounding_mode (get_rounding_mode ()) in
s let () = Kernel.failure "pretty: rounding mode (%s) <> FE_TONEAREST" mode in
firstdigit set_round_nearest_even ()
decdigits
exp
end
else
Format.fprintf fmt "%s0x%d.%013Lxp%d"
s
firstdigit
man
exp
let pretty fmt f = 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 *) (* should always arrive here with nearest_even *)
if get_rounding_mode () <> FE_ToNearest then begin ensure_round_nearest_even () ;
Kernel.failure "pretty: rounding mode (%s) <> FE_TONEAREST" if not (use_hex || Kernel.FloatNormal.get ()) then
(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
let r = Format.sprintf "%.*g" 12 f in let r = Format.sprintf "%.*g" 12 f in
if (String.contains r '.' || String.contains r 'e' || let dot = if is_not_integer r || is_not_finite f then "" else "." in
String.contains r 'E') Format.fprintf fmt "%s%s" r dot
|| (match classify_float f with else pretty_normal ~use_hex fmt f
| 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
type sign = Neg | Pos type sign = Neg | Pos
exception Float_Non_representable_as_Int64 of sign 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], (* 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 raise Float_Non_representable_as_Int64. This is the most reasonable as
a floating-point number may represent an exponentially large integer. *) a floating-point number may represent an exponentially large integer. *)
let truncate_to_integer = let truncate_to_integer x =
let min_64_float = -9.22337203685477581e+18 if x < min_64_float then raise (Float_Non_representable_as_Int64 Neg) ;
(* Int64.to_float (-0x8000000000000000L) *) if x > 2. *. max_64_float then raise (Float_Non_representable_as_Int64 Pos) ;
in let convert x = Integer.of_int64 (Int64.of_float x) in
let max_64_float = 9.22337203685477478e+18 let shift n = Integer.(add (two_power_of_int 63)) n in
(* let open Int64 in if x <= max_64_float then convert x else convert (x +. min_64_float) |> shift
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 bits_of_max_double = let bits_of_max_double =
Integer.of_int64 (Int64.bits_of_float max_float) Integer.of_int64 (Int64.bits_of_float max_float)
let bits_of_most_negative_double = let bits_of_most_negative_double =
Integer.of_int64 (Int64.bits_of_float (-. max_float)) Integer.of_int64 (Int64.bits_of_float (-. max_float))
(** See e.g. http://www.h-schmidt.net/FloatConverter/IEEE754.html *) (** See e.g. http://www.h-schmidt.net/FloatConverter/IEEE754.html *)
let bits_of_max_float = Integer.of_int64 0x7F7FFFFFL let bits_of_max_float = Integer.of_int64 0x7F7FFFFFL
let bits_of_most_negative_float = 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 Integer.of_int64 v
external fround: float -> float = "c_round" external fround: float -> float = "c_round"
external trunc: float -> float = "c_trunc" external trunc: float -> float = "c_trunc"
...@@ -427,6 +363,7 @@ external atanf: float -> float = "c_atanf" ...@@ -427,6 +363,7 @@ external atanf: float -> float = "c_atanf"
external atan2f: float -> float -> float = "c_atan2f" external atan2f: float -> float -> float = "c_atan2f"
(** C math-like functions *) (** C math-like functions *)
let isnan f = let isnan f =
...@@ -447,27 +384,31 @@ let neg_min_single_precision_denormal = -. min_single_precision_denormal ...@@ -447,27 +384,31 @@ let neg_min_single_precision_denormal = -. min_single_precision_denormal
(* auxiliary functions for nextafter/nextafterf *) (* auxiliary functions for nextafter/nextafterf *)
let min_denormal_float ~is_f32 = let min_denormal_float ~is_f32 =
if is_f32 then min_single_precision_denormal else min_denormal if is_f32 then min_single_precision_denormal else min_denormal
let nextafter_aux ~is_f32 fincr fdecr x y = let nextafter_aux ~is_f32 fincr fdecr x y =
if x = y (* includes cases "(0.0, -0.0) => -0.0" and its symmetric *) let sign = if x < y then 1.0 else -.1.0 in
then y if x = y then y
else if isnan x || isnan y then nan else if isnan x || isnan y then nan
else if x = 0.0 (* or -0.0 *) then else if x = 0.0 then sign *. min_denormal_float ~is_f32
if x < y then min_denormal_float ~is_f32 else -. (min_denormal_float ~is_f32) else if x = neg_infinity then fdecr x
(* the following conditions might be simpler if we had unsigned ints else if (x < y && x > 0.0) || (x > y && x < 0.0) then fincr x
(uint32/uint64) *) else fdecr x
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
let incr_f64 f = 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 = let decr_f64 f =
if f = infinity then max_float 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 = let incr_f32 f =
if f = neg_infinity then most_negative_single_precision_float 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 = let decr_f32 f =
if f = infinity then max_single_precision_float 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 = let nextafter x y =
nextafter_aux ~is_f32:false incr_f64 decr_f64 x y nextafter_aux ~is_f32:false incr_f64 decr_f64 x y
......
...@@ -767,10 +767,9 @@ module Make ...@@ -767,10 +767,9 @@ module Make
subvalues that are all evaluated. Limits the number of splits to subvalues that are all evaluated. Limits the number of splits to
keep the number of evaluations linear on [nb]. *) keep the number of evaluations linear on [nb]. *)
let subdivnb = let subdivnb =
let pow n e = Integer.power_int_positive_int n e |> Option.get in
if nb > 3 if nb > 3
then then (subdivnb * nb) / (Integer.to_int_exn (pow 2 (nb - 1)))
let pow = Integer.power_int_positive_int in
(subdivnb * nb) / (Integer.to_int_exn (pow 2 (nb - 1)))
else subdivnb else subdivnb
in in
Self.result ~current:true ~once:true ~dkey Self.result ~current:true ~once:true ~dkey
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment