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/float64_63.ml | 35 +++++++++++++++++++++++++++++++++++
 1 file changed, 35 insertions(+)
 create mode 100644 kernel/float64_63.ml

(limited to 'kernel/float64_63.ml')

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."
-- 
cgit v1.2.3