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.

6 files changed, 172 insertions, 25 deletions

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_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             *)
+(* <O___,, * (see version control and CREDITS file for authors & dates) *)
+(*   \VV/  **************************************************************)
+(*    //   *    This file is distributed under the terms of the         *)
+(*         *     GNU Lesser General Public License Version 2.1          *)
+(*         *     (see LICENSE file for the text of the license)         *)
+(************************************************************************)
+
+include Float64_common
+
+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]
+
+(*** 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             *)
+(* <O___,, * (see version control and CREDITS file for authors & dates) *)
+(*   \VV/  **************************************************************)
+(*    //   *    This file is distributed under the terms of the         *)
+(*         *     GNU Lesser General Public License Version 2.1          *)
+(*         *     (see LICENSE file for the text of the license)         *)
+(************************************************************************)
+
+include Float64_common
+
+let mul (x : float) (y : float) : float = x *. y
+[@@ocaml.inline always]
+
+let add (x : float) (y : float) : float = x +. y
+[@@ocaml.inline always]
+
+let sub (x : float) (y : float) : float = x -. y
+[@@ocaml.inline always]
+
+let div (x : float) (y : float) : float = x /. y
+[@@ocaml.inline always]
+
+let sqrt (x : float) : float = sqrt x
+[@@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.ml b/kernel/float64_common.ml
index 76005a3dc6..2991a20b49 100644
--- a/kernel/float64.ml
+++ b/kernel/float64_common.ml
@@ -88,21 +88,6 @@ let classify x =
   | 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]
 
@@ -157,12 +142,3 @@ let total_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_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             *)
+(* <O___,, * (see version control and CREDITS file for authors & dates) *)
+(*   \VV/  **************************************************************)
+(*    //   *    This file is distributed under the terms of the         *)
+(*         *     GNU Lesser General Public License Version 2.1          *)
+(*         *     (see LICENSE file for the text of the license)         *)
+(************************************************************************)
+
+(** [t] is currently implemented by OCaml's [float] type.
+
+Beware: NaNs have a sign and a payload, while they should be
+indistinguishable from Coq's perspective. *)
+type t = float
+
+(** Test functions for special values to avoid calling [classify] *)
+val is_nan : t -> 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