From 6fe8c44ff828ef4ec89b49ada634ce87639f384f Mon Sep 17 00:00:00 2001 From: Pierre Roux Date: Tue, 6 Oct 2020 16:52:03 +0200 Subject: Use OCaml floating-point operations on 64 bits arch C functions were used for floating-point arithmetic operations, by fear of double rounding that could happen on old x87 on 32 bits architecture. This commit uses OCaml floating-point operations on 64 bits architectures. The following snippet is made 17% faster by this commit. From Coq Require Import Int63 BinPos PrimFloat. Definition foo n := let eps := sub (next_up one) one in Pos.iter (fun x => add x eps) two n. Time Eval native_compute in foo 1000000000. --- kernel/dune | 7 +- kernel/float64.ml | 168 ---------------------------------------------- kernel/float64_31.ml | 35 ++++++++++ kernel/float64_63.ml | 35 ++++++++++ kernel/float64_common.ml | 144 +++++++++++++++++++++++++++++++++++++++ kernel/float64_common.mli | 95 ++++++++++++++++++++++++++ kernel/kernel.mllib | 1 + 7 files changed, 316 insertions(+), 169 deletions(-) delete mode 100644 kernel/float64.ml create mode 100644 kernel/float64_31.ml create mode 100644 kernel/float64_63.ml create mode 100644 kernel/float64_common.ml create mode 100644 kernel/float64_common.mli (limited to 'kernel') diff --git a/kernel/dune b/kernel/dune index ce6fdc03df..bd663974da 100644 --- a/kernel/dune +++ b/kernel/dune @@ -3,7 +3,7 @@ (synopsis "The Coq Kernel") (public_name coq.kernel) (wrapped false) - (modules (:standard \ genOpcodeFiles uint63_31 uint63_63)) + (modules (:standard \ genOpcodeFiles uint63_31 uint63_63 float64_31 float64_63)) (libraries lib byterun dynlink)) (executable @@ -19,6 +19,11 @@ (deps (:gen-file uint63_%{ocaml-config:int_size}.ml)) (action (copy# %{gen-file} %{targets}))) +(rule + (targets float64.ml) + (deps (:gen-file float64_%{ocaml-config:int_size}.ml)) + (action (copy# %{gen-file} %{targets}))) + (documentation (package coq)) diff --git a/kernel/float64.ml b/kernel/float64.ml deleted file mode 100644 index 76005a3dc6..0000000000 --- a/kernel/float64.ml +++ /dev/null @@ -1,168 +0,0 @@ -(************************************************************************) -(* * The Coq Proof Assistant / The Coq Development Team *) -(* v * Copyright INRIA, CNRS and contributors *) -(* f - as comparison on floats rather than the polymorphic OCaml comparison - which is much slower. *) -let is_nan (f : float) = f <> f -let is_infinity f = f = infinity -let is_neg_infinity f = f = neg_infinity - -(* OCaml gives a sign to nan values which should not be displayed as - all NaNs are considered equal here. - OCaml prints infinities as "inf" (resp. "-inf") - but we want "infinity" (resp. "neg_infinity"). *) -let to_string_raw fmt f = - if is_nan f then "nan" - else if is_infinity f then "infinity" - else if is_neg_infinity f then "neg_infinity" - else Printf.sprintf fmt f - -let to_hex_string = to_string_raw "%h" - -(* Printing a binary64 float in 17 decimal places and parsing it again - will yield the same float. *) -let to_string = to_string_raw "%.17g" - -let of_string = float_of_string - -(* Compiles a float to OCaml code *) -let compile f = - Printf.sprintf "Float64.of_float (%s)" (to_hex_string f) - -let of_float f = f - -let sign f = copysign 1. f < 0. - -let opp = ( ~-. ) -let abs = abs_float - -type float_comparison = FEq | FLt | FGt | FNotComparable - -(* See above comment on [is_nan] for the [float] type annotations. *) -let eq (x : float) (y : float) = x = y -[@@ocaml.inline always] - -let lt (x : float) (y : float) = x < y -[@@ocaml.inline always] - -let le (x : float) (y : float) = x <= y -[@@ocaml.inline always] - -(* inspired by lib/util.ml; see also #10471 *) -let pervasives_compare (x : float) (y : float) = compare x y - -let compare (x : float) (y : float) = - if x < y then FLt - else - ( - if x > y then FGt - else - ( - if x = y then FEq - else FNotComparable (* NaN case *) - ) - ) -[@@ocaml.inline always] - -type float_class = - | PNormal | NNormal | PSubn | NSubn | PZero | NZero | PInf | NInf | NaN - -let classify x = - match classify_float x with - | FP_normal -> if 0. < x then PNormal else NNormal - | FP_subnormal -> if 0. < x then PSubn else NSubn - | FP_zero -> if 0. < 1. /. x then PZero else NZero - | FP_infinite -> if 0. < x then PInf else NInf - | FP_nan -> NaN -[@@ocaml.inline always] - -external mul : float -> float -> float = "coq_fmul_byte" "coq_fmul" -[@@unboxed] [@@noalloc] - -external add : float -> float -> float = "coq_fadd_byte" "coq_fadd" -[@@unboxed] [@@noalloc] - -external sub : float -> float -> float = "coq_fsub_byte" "coq_fsub" -[@@unboxed] [@@noalloc] - -external div : float -> float -> float = "coq_fdiv_byte" "coq_fdiv" -[@@unboxed] [@@noalloc] - -external sqrt : float -> float = "coq_fsqrt_byte" "coq_fsqrt" -[@@unboxed] [@@noalloc] - -let of_int63 x = Uint63.to_float x -[@@ocaml.inline always] - -let prec = 53 -let normfr_mantissa f = - let f = abs f in - if f >= 0.5 && f < 1. then Uint63.of_float (ldexp f prec) - else Uint63.zero -[@@ocaml.inline always] - -let eshift = 2101 (* 2*emax + prec *) - -(* When calling frexp on a nan or an infinity, the returned value inside - the exponent is undefined. - Therefore we must always set it to a fixed value (here 0). *) -let frshiftexp f = - match classify_float f with - | FP_zero | FP_infinite | FP_nan -> (f, Uint63.zero) - | FP_normal | FP_subnormal -> - let (m, e) = frexp f in - m, Uint63.of_int (e + eshift) -[@@ocaml.inline always] - -let ldshiftexp f e = ldexp f (Uint63.to_int_min e (2 * eshift) - eshift) -[@@ocaml.inline always] - -external next_up : float -> float = "coq_next_up_byte" "coq_next_up" -[@@unboxed] [@@noalloc] - -external next_down : float -> float = "coq_next_down_byte" "coq_next_down" -[@@unboxed] [@@noalloc] - -let equal f1 f2 = - match classify_float f1 with - | FP_normal | FP_subnormal | FP_infinite -> (f1 = f2) - | FP_nan -> is_nan f2 - | FP_zero -> f1 = f2 && 1. /. f1 = 1. /. f2 (* OCaml consider 0. = -0. *) -[@@ocaml.inline always] - -let hash = - (* Hashtbl.hash already considers all NaNs as equal, - cf. https://github.com/ocaml/ocaml/commit/aea227fdebe0b5361fd3e1d0aaa42cf929052269 - and http://caml.inria.fr/pub/docs/manual-ocaml/libref/Hashtbl.html *) - Hashtbl.hash - -let total_compare f1 f2 = - (* pervasives_compare considers all NaNs as equal, which is fine here, - but also considers -0. and +0. as equal *) - if f1 = 0. && f2 = 0. then pervasives_compare (1. /. f1) (1. /. f2) - else pervasives_compare f1 f2 - -let is_float64 t = - Obj.tag t = Obj.double_tag -[@@ocaml.inline always] - -(*** Test at runtime that no harmful double rounding seems to - be performed with an intermediate 80 bits representation (x87). *) -let () = - let b = ldexp 1. 53 in - let s = add 1. (ldexp 1. (-52)) in - if add b s <= b || add b 1. <> b || ldexp 1. (-1074) <= 0. then - failwith "Detected non IEEE-754 compliant architecture (or wrong \ - rounding mode). Use of Float is thus unsafe." diff --git a/kernel/float64_31.ml b/kernel/float64_31.ml new file mode 100644 index 0000000000..09b28e6cf0 --- /dev/null +++ b/kernel/float64_31.ml @@ -0,0 +1,35 @@ +(************************************************************************) +(* * The Coq Proof Assistant / The Coq Development Team *) +(* v * Copyright INRIA, CNRS and contributors *) +(* float -> float = "coq_fmul_byte" "coq_fmul" +[@@unboxed] [@@noalloc] + +external add : float -> float -> float = "coq_fadd_byte" "coq_fadd" +[@@unboxed] [@@noalloc] + +external sub : float -> float -> float = "coq_fsub_byte" "coq_fsub" +[@@unboxed] [@@noalloc] + +external div : float -> float -> float = "coq_fdiv_byte" "coq_fdiv" +[@@unboxed] [@@noalloc] + +external sqrt : float -> float = "coq_fsqrt_byte" "coq_fsqrt" +[@@unboxed] [@@noalloc] + +(*** Test at runtime that no harmful double rounding seems to + be performed with an intermediate 80 bits representation (x87). *) +let () = + let b = ldexp 1. 53 in + let s = add 1. (ldexp 1. (-52)) in + if add b s <= b || add b 1. <> b || ldexp 1. (-1074) <= 0. then + failwith "Detected non IEEE-754 compliant architecture (or wrong \ + rounding mode). Use of Float is thus unsafe." diff --git a/kernel/float64_63.ml b/kernel/float64_63.ml new file mode 100644 index 0000000000..0025531cb1 --- /dev/null +++ b/kernel/float64_63.ml @@ -0,0 +1,35 @@ +(************************************************************************) +(* * The Coq Proof Assistant / The Coq Development Team *) +(* v * Copyright INRIA, CNRS and contributors *) +(* b || ldexp 1. (-1074) <= 0. then + failwith "Detected non IEEE-754 compliant architecture (or wrong \ + rounding mode). Use of Float is thus unsafe." diff --git a/kernel/float64_common.ml b/kernel/float64_common.ml new file mode 100644 index 0000000000..2991a20b49 --- /dev/null +++ b/kernel/float64_common.ml @@ -0,0 +1,144 @@ +(************************************************************************) +(* * The Coq Proof Assistant / The Coq Development Team *) +(* v * Copyright INRIA, CNRS and contributors *) +(* f + as comparison on floats rather than the polymorphic OCaml comparison + which is much slower. *) +let is_nan (f : float) = f <> f +let is_infinity f = f = infinity +let is_neg_infinity f = f = neg_infinity + +(* OCaml gives a sign to nan values which should not be displayed as + all NaNs are considered equal here. + OCaml prints infinities as "inf" (resp. "-inf") + but we want "infinity" (resp. "neg_infinity"). *) +let to_string_raw fmt f = + if is_nan f then "nan" + else if is_infinity f then "infinity" + else if is_neg_infinity f then "neg_infinity" + else Printf.sprintf fmt f + +let to_hex_string = to_string_raw "%h" + +(* Printing a binary64 float in 17 decimal places and parsing it again + will yield the same float. *) +let to_string = to_string_raw "%.17g" + +let of_string = float_of_string + +(* Compiles a float to OCaml code *) +let compile f = + Printf.sprintf "Float64.of_float (%s)" (to_hex_string f) + +let of_float f = f + +let sign f = copysign 1. f < 0. + +let opp = ( ~-. ) +let abs = abs_float + +type float_comparison = FEq | FLt | FGt | FNotComparable + +(* See above comment on [is_nan] for the [float] type annotations. *) +let eq (x : float) (y : float) = x = y +[@@ocaml.inline always] + +let lt (x : float) (y : float) = x < y +[@@ocaml.inline always] + +let le (x : float) (y : float) = x <= y +[@@ocaml.inline always] + +(* inspired by lib/util.ml; see also #10471 *) +let pervasives_compare (x : float) (y : float) = compare x y + +let compare (x : float) (y : float) = + if x < y then FLt + else + ( + if x > y then FGt + else + ( + if x = y then FEq + else FNotComparable (* NaN case *) + ) + ) +[@@ocaml.inline always] + +type float_class = + | PNormal | NNormal | PSubn | NSubn | PZero | NZero | PInf | NInf | NaN + +let classify x = + match classify_float x with + | FP_normal -> if 0. < x then PNormal else NNormal + | FP_subnormal -> if 0. < x then PSubn else NSubn + | FP_zero -> if 0. < 1. /. x then PZero else NZero + | FP_infinite -> if 0. < x then PInf else NInf + | FP_nan -> NaN +[@@ocaml.inline always] + +let of_int63 x = Uint63.to_float x +[@@ocaml.inline always] + +let prec = 53 +let normfr_mantissa f = + let f = abs f in + if f >= 0.5 && f < 1. then Uint63.of_float (ldexp f prec) + else Uint63.zero +[@@ocaml.inline always] + +let eshift = 2101 (* 2*emax + prec *) + +(* When calling frexp on a nan or an infinity, the returned value inside + the exponent is undefined. + Therefore we must always set it to a fixed value (here 0). *) +let frshiftexp f = + match classify_float f with + | FP_zero | FP_infinite | FP_nan -> (f, Uint63.zero) + | FP_normal | FP_subnormal -> + let (m, e) = frexp f in + m, Uint63.of_int (e + eshift) +[@@ocaml.inline always] + +let ldshiftexp f e = ldexp f (Uint63.to_int_min e (2 * eshift) - eshift) +[@@ocaml.inline always] + +external next_up : float -> float = "coq_next_up_byte" "coq_next_up" +[@@unboxed] [@@noalloc] + +external next_down : float -> float = "coq_next_down_byte" "coq_next_down" +[@@unboxed] [@@noalloc] + +let equal f1 f2 = + match classify_float f1 with + | FP_normal | FP_subnormal | FP_infinite -> (f1 = f2) + | FP_nan -> is_nan f2 + | FP_zero -> f1 = f2 && 1. /. f1 = 1. /. f2 (* OCaml consider 0. = -0. *) +[@@ocaml.inline always] + +let hash = + (* Hashtbl.hash already considers all NaNs as equal, + cf. https://github.com/ocaml/ocaml/commit/aea227fdebe0b5361fd3e1d0aaa42cf929052269 + and http://caml.inria.fr/pub/docs/manual-ocaml/libref/Hashtbl.html *) + Hashtbl.hash + +let total_compare f1 f2 = + (* pervasives_compare considers all NaNs as equal, which is fine here, + but also considers -0. and +0. as equal *) + if f1 = 0. && f2 = 0. then pervasives_compare (1. /. f1) (1. /. f2) + else pervasives_compare f1 f2 + +let is_float64 t = + Obj.tag t = Obj.double_tag +[@@ocaml.inline always] diff --git a/kernel/float64_common.mli b/kernel/float64_common.mli new file mode 100644 index 0000000000..4fb1c114a5 --- /dev/null +++ b/kernel/float64_common.mli @@ -0,0 +1,95 @@ +(************************************************************************) +(* * The Coq Proof Assistant / The Coq Development Team *) +(* v * Copyright INRIA, CNRS and contributors *) +(* bool +val is_infinity : t -> bool +val is_neg_infinity : t -> bool + +val of_string : string -> t + +(** Print a float exactly as an hexadecimal value (exact decimal + * printing would be possible but sometimes requires more than 700 + * digits). *) +val to_hex_string : t -> string + +(** Print a float as a decimal value. The printing is not exact (the + * real value printed is not always the given floating-point value), + * however printing is precise enough that forall float [f], + * [of_string (to_decimal_string f) = f]. *) +val to_string : t -> string + +val compile : t -> string + +val of_float : float -> t + +(** Return [true] for "-", [false] for "+". *) +val sign : t -> bool + +val opp : t -> t +val abs : t -> t + +type float_comparison = FEq | FLt | FGt | FNotComparable + +val eq : t -> t -> bool + +val lt : t -> t -> bool + +val le : t -> t -> bool + +(** The IEEE 754 float comparison. + * NotComparable is returned if there is a NaN in the arguments *) +val compare : t -> t -> float_comparison +[@@ocaml.inline always] + +type float_class = + | PNormal | NNormal | PSubn | NSubn | PZero | NZero | PInf | NInf | NaN + +val classify : t -> float_class +[@@ocaml.inline always] + +(** Link with integers *) +val of_int63 : Uint63.t -> t +[@@ocaml.inline always] + +val normfr_mantissa : t -> Uint63.t +[@@ocaml.inline always] + +(** Shifted exponent extraction *) +val eshift : int + +val frshiftexp : t -> t * Uint63.t (* float remainder, shifted exponent *) +[@@ocaml.inline always] + +val ldshiftexp : t -> Uint63.t -> t +[@@ocaml.inline always] + +val next_up : t -> t + +val next_down : t -> t + +(** Return true if two floats are equal. + * All NaN values are considered equal. *) +val equal : t -> t -> bool +[@@ocaml.inline always] + +val hash : t -> int + +(** Total order relation over float values. Behaves like [Pervasives.compare].*) +val total_compare : t -> t -> int + +val is_float64 : Obj.t -> bool +[@@ocaml.inline always] diff --git a/kernel/kernel.mllib b/kernel/kernel.mllib index d4d7150222..5b2a7bd9c2 100644 --- a/kernel/kernel.mllib +++ b/kernel/kernel.mllib @@ -2,6 +2,7 @@ Names TransparentState Uint63 Parray +Float64_common Float64 Univ UGraph -- cgit v1.2.3