Make primitive float work on x86_32

Flag -fexcess-precision=standard is not enough on x86_32
where -msse2 -mfpmath=sse is required (-msse is not enough)
to avoid double rounding issues in the VM.

Most floating-point operation are now implemented in C because OCaml
is suffering double rounding issues on x86_32 with 80 bits extended
precision registers used for floating-point values, causing double
rounding making floating-point arithmetic incorrect with respect to
its specification.

Add a runtime test for double roundings.

12 files changed, 139 insertions, 57 deletions

diff --git a/Makefile.build b/Makefile.build
index f9286c9635..90b3408d79 100644
--- a/Makefile.build
+++ b/Makefile.build
@@ -568,13 +568,13 @@ $(FAKEIDEBYTE): $(FAKEIDECMO) | $(IDETOPBYTE)
 
 VOTOURCMO:=clib/cObj.cmo kernel/uint63.cmo kernel/float64.cmo checker/analyze.cmo checker/values.cmo checker/votour.cmo
 
-bin/votour: $(call bestobj, $(VOTOURCMO))
+bin/votour: $(call bestobj, $(VOTOURCMO)) $(LIBCOQRUN)
 	$(SHOW)'OCAMLBEST -o $@'
-	$(HIDE)$(call bestocaml, -I checker)
+	$(HIDE)$(call bestocaml, -I checker -cclib -lcoqrun) $(VMBYTEFLAGS)
 
-bin/votour.byte: $(VOTOURCMO)
+bin/votour.byte: $(VOTOURCMO) $(LIBCOQRUN)
 	$(SHOW)'OCAMLC -o $@'
-	$(HIDE)$(call ocamlbyte, -I checker)
+	$(HIDE)$(call ocamlbyte, -I checker -cclib -lcoqrun) $(VMBYTEFLAGS)
 
 ###########################################################################
 # Csdp to micromega special targets
diff --git a/configure.ml b/configure.ml
index 7f678246dd..4adde7f956 100644
--- a/configure.ml
+++ b/configure.ml
@@ -456,7 +456,7 @@ let coq_bin_annot_flag = if !prefs.bin_annot then "-bin-annot" else ""
 let coq_safe_string = "-safe-string"
 let coq_strict_sequence = "-strict-sequence"
 
-let cflags = "-Wall -Wno-unused -g -O2 -fexcess-precision=standard"
+let cflags = "-Wall -Wno-unused -g -O2 -fexcess-precision=standard -msse2 -mfpmath=sse"
 
 (** * Architecture *)
 
diff --git a/kernel/byterun/coq_float64.h b/kernel/byterun/coq_float64.h
new file mode 100644
index 0000000000..6814c31642
--- /dev/null
+++ b/kernel/byterun/coq_float64.h
@@ -0,0 +1,32 @@
+#ifndef _COQ_FLOAT64_
+#define _COQ_FLOAT64_
+
+#include <math.h>
+
+#define DECLARE_FBINOP(name, e)                                         \
+  double coq_##name(double x, double y) {                               \
+    return e;                                                           \
+  }                                                                     \
+                                                                        \
+  value coq_##name##_byte(value x, value y) {                           \
+    return caml_copy_double(coq_##name(Double_val(x), Double_val(y)));  \
+  }
+
+#define DECLARE_FUNOP(name, e)                                          \
+  double coq_##name(double x) {                                         \
+    return e;                                                           \
+  }                                                                     \
+                                                                        \
+  value coq_##name##_byte(value x) {                                    \
+    return caml_copy_double(coq_##name(Double_val(x)));                 \
+  }
+
+DECLARE_FBINOP(fmul, x * y)
+DECLARE_FBINOP(fadd, x + y)
+DECLARE_FBINOP(fsub, x - y)
+DECLARE_FBINOP(fdiv, x / y)
+DECLARE_FUNOP(fsqrt, sqrt(x))
+DECLARE_FUNOP(next_up, nextafter(x, INFINITY))
+DECLARE_FUNOP(next_down, nextafter(x, -INFINITY))
+
+#endif /* _COQ_FLOAT64_ */
diff --git a/kernel/byterun/coq_interp.c b/kernel/byterun/coq_interp.c
index 06042bb753..e67325eb9a 100644
--- a/kernel/byterun/coq_interp.c
+++ b/kernel/byterun/coq_interp.c
@@ -23,6 +23,7 @@
 #include "coq_fix_code.h"
 #include "coq_memory.h" 
 #include "coq_values.h" 
+#include "coq_float64.h"
 
 #ifdef ARCH_SIXTYFOUR
 #include "coq_uint63_native.h"
@@ -1593,42 +1594,42 @@ value coq_interprete
       Instruct (CHECKADDFLOAT) {
         print_instr("CHECKADDFLOAT");
         CheckFloat2();
-        Coq_copy_double(Double_val(accu) + Double_val(*sp++));
+        Coq_copy_double(coq_fadd(Double_val(accu), Double_val(*sp++)));
         Next;
       }
 
       Instruct (CHECKSUBFLOAT) {
         print_instr("CHECKSUBFLOAT");
         CheckFloat2();
-        Coq_copy_double(Double_val(accu) - Double_val(*sp++));
+        Coq_copy_double(coq_fsub(Double_val(accu), Double_val(*sp++)));
         Next;
       }
 
       Instruct (CHECKMULFLOAT) {
         print_instr("CHECKMULFLOAT");
         CheckFloat2();
-        Coq_copy_double(Double_val(accu) * Double_val(*sp++));
+        Coq_copy_double(coq_fmul(Double_val(accu), Double_val(*sp++)));
         Next;
       }
 
       Instruct (CHECKDIVFLOAT) {
         print_instr("CHECKDIVFLOAT");
         CheckFloat2();
-        Coq_copy_double(Double_val(accu) / Double_val(*sp++));
+        Coq_copy_double(coq_fdiv(Double_val(accu), Double_val(*sp++)));
         Next;
       }
 
       Instruct (CHECKSQRTFLOAT) {
         print_instr("CHECKSQRTFLOAT");
         CheckFloat1();
-        Coq_copy_double(sqrt(Double_val(accu)));
+        Coq_copy_double(coq_fsqrt(Double_val(accu)));
         Next;
       }
 
       Instruct (CHECKFLOATOFINT63) {
         print_instr("CHECKFLOATOFINT63");
         CheckInt1();
-        Coq_copy_double(uint63_to_double(accu));
+        Uint63_to_double(accu);
         Next;
       }
 
@@ -1638,10 +1639,10 @@ value coq_interprete
         CheckFloat1();
         f = fabs(Double_val(accu));
         if (f >= 0.5 && f < 1) {
-          accu = uint63_of_double(ldexp(f, DBL_MANT_DIG));
+          Uint63_of_double(ldexp(f, DBL_MANT_DIG));
         }
         else {
-          accu = Val_int(0);
+          Uint63_of_int(Val_int(0));
         }
         Next;
       }
@@ -1660,31 +1661,39 @@ value coq_interprete
         }
         Coq_copy_double(f);
         *--sp = accu;
+#ifdef ARCH_SIXTYFOUR
         Alloc_small(accu, 2, coq_tag_pair);
-        Field(accu, 0) = *sp++;
         Field(accu, 1) = Val_int(exp);
+#else
+        Uint63_of_int(Val_int(exp));
+        *--sp = accu;
+        Alloc_small(accu, 2, coq_tag_pair);
+        Field(accu, 1) = *sp++;
+#endif
+        Field(accu, 0) = *sp++;
         Next;
       }
 
       Instruct (CHECKLDSHIFTEXP) {
         print_instr("CHECKLDSHIFTEXP");
         CheckPrimArgs(Is_double(accu) && Is_uint63(sp[0]), apply2);
-        Coq_copy_double(ldexp(Double_val(accu),
-                              uint63_of_value(*sp++) - FLOAT_EXP_SHIFT));
+        Swap_accu_sp;
+        Int_of_uint63(accu);
+        Coq_copy_double(ldexp(Double_val(*sp++), accu - FLOAT_EXP_SHIFT));
         Next;
       }
 
       Instruct (CHECKNEXTUPFLOAT) {
         print_instr("CHECKNEXTUPFLOAT");
         CheckFloat1();
-        Coq_copy_double(nextafter(Double_val(accu), INFINITY));
+        Coq_copy_double(coq_next_up(Double_val(accu)));
         Next;
       }
 
       Instruct (CHECKNEXTDOWNFLOAT) {
         print_instr("CHECKNEXTDOWNFLOAT");
         CheckFloat1();
-        Coq_copy_double(nextafter(Double_val(accu), -INFINITY));
+        Coq_copy_double(coq_next_down(Double_val(accu)));
         Next;
       }
 
diff --git a/kernel/byterun/coq_uint63_emul.h b/kernel/byterun/coq_uint63_emul.h
index 528cc6fc1f..e09803ae2d 100644
--- a/kernel/byterun/coq_uint63_emul.h
+++ b/kernel/byterun/coq_uint63_emul.h
@@ -156,3 +156,18 @@ DECLARE_BINOP(mulc_ml)
     *(h) = Field(uint63_return_value__, 0); \
     accu = Field(uint63_return_value__, 1); \
   }while(0)
+
+DECLARE_UNOP(to_float)
+#define Uint63_to_double(x) CALL_UNOP(to_float, x)
+
+DECLARE_UNOP(of_float)
+#define Uint63_of_double(f) do{ \
+  Coq_copy_double(f); \
+  CALL_UNOP(of_float, accu); \
+  }while(0)
+
+DECLARE_UNOP(of_int)
+#define Uint63_of_int(x) CALL_UNOP(of_int, x)
+
+DECLARE_UNOP(to_int_saturate)
+#define Int_of_uint63(x) CALL_PREDICATE(accu, to_int_saturate, x)
diff --git a/kernel/byterun/coq_uint63_native.h b/kernel/byterun/coq_uint63_native.h
index a14ed5c262..0ab374194e 100644
--- a/kernel/byterun/coq_uint63_native.h
+++ b/kernel/byterun/coq_uint63_native.h
@@ -139,5 +139,26 @@ value uint63_div21(value xh, value xl, value y, value* ql) {
 }
 #define Uint63_div21(xh, xl, y, q) (accu = uint63_div21(xh, xl, y, q))
 
-#define uint63_to_double(val) ((double) uint63_of_value(val))
-#define uint63_of_double(f) (Val_long((long int) f))
+#define Uint63_to_double(x) Coq_copy_double((double) uint63_of_value(x))
+
+double coq_uint63_to_float(value x) {
+  return (double) uint63_of_value(x);
+}
+
+value coq_uint63_to_float_byte(value x) {
+  return caml_copy_double(coq_uint63_to_float(x));
+}
+
+#define Uint63_of_double(f) do{ \
+  accu = Val_long((uint64_t)(f)); \
+  }while(0)
+
+#define Uint63_of_int(x) (accu = (x))
+
+#define Int_of_uint63(x) do{ \
+  uint64_t int_of_uint63__ = uint63_of_value(x); \
+  if ((int_of_uint63__ & 0xFFFFFFFF80000000L) == 0) \
+    accu = (int)int_of_uint63__; \
+  else \
+    accu = 0x7FFFFFFF; \
+  }while(0)
diff --git a/kernel/float64.ml b/kernel/float64.ml
index 5160d32892..55ad472ea9 100644
--- a/kernel/float64.ml
+++ b/kernel/float64.ml
@@ -70,20 +70,20 @@ let classify x =
   | FP_nan -> NaN
 [@@ocaml.inline always]
 
-let mul x y = x *. y
-[@@ocaml.inline always]
+external mul : float -> float -> float = "coq_fmul_byte" "coq_fmul"
+[@@unboxed] [@@noalloc]
 
-let add x y = x +. y
-[@@ocaml.inline always]
+external add : float -> float -> float = "coq_fadd_byte" "coq_fadd"
+[@@unboxed] [@@noalloc]
 
-let sub x y = x -. y
-[@@ocaml.inline always]
+external sub : float -> float -> float = "coq_fsub_byte" "coq_fsub"
+[@@unboxed] [@@noalloc]
 
-let div x y = x /. y
-[@@ocaml.inline always]
+external div : float -> float -> float = "coq_fdiv_byte" "coq_fdiv"
+[@@unboxed] [@@noalloc]
 
-let sqrt x = sqrt x
-[@@ocaml.inline always]
+external sqrt : float -> float = "coq_fsqrt_byte" "coq_fsqrt"
+[@@unboxed] [@@noalloc]
 
 let of_int63 x = Uint63.to_float x
 [@@ocaml.inline always]
@@ -111,25 +111,11 @@ let frshiftexp f =
 let ldshiftexp f e = ldexp f (snd (Uint63.to_int2 e) - eshift)
 [@@ocaml.inline always]
 
-let eta_float = ldexp 1. (-1074) (* smallest positive float (subnormal) *)
+external next_up : float -> float = "coq_next_up_byte" "coq_next_up"
+[@@unboxed] [@@noalloc]
 
-let next_up f =
-  match classify_float f with
-  | FP_nan -> f
-  | FP_infinite -> if 0. < f then f else -.max_float
-  | FP_zero | FP_subnormal ->
-     let f = f +. eta_float in
-     if f = 0. then -0. else f (* or next_down may return -0. *)
-  | FP_normal ->
-     let f, e = frexp f in
-     if 0. < f || f <> -0.5 || e = -1021 then
-       ldexp (f +. epsilon_float /. 2.) e
-     else
-       ldexp (-0.5 +. epsilon_float /. 4.) e
-[@@ocaml.inline always]
-
-let next_down f = -.(next_up (-.f))
-[@@ocaml.inline always]
+external next_down : float -> float = "coq_next_down_byte" "coq_next_down"
+[@@unboxed] [@@noalloc]
 
 let equal f1 f2 =
   match classify_float f1 with
@@ -153,3 +139,13 @@ 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 then
+    failwith "Detected double rounding due to the use of intermediate \
+              80 bits floating-point representation. Use of Float is \
+              thus unsafe."
diff --git a/kernel/float64.mli b/kernel/float64.mli
index acc3a556ab..78dc1a7bd7 100644
--- a/kernel/float64.mli
+++ b/kernel/float64.mli
@@ -46,19 +46,14 @@ val classify : t -> float_class
 [@@ocaml.inline always]
 
 val mul : t -> t -> t
-[@@ocaml.inline always]
 
 val add : t -> t -> t
-[@@ocaml.inline always]
 
 val sub : t -> t -> t
-[@@ocaml.inline always]
 
 val div : t -> t -> t
-[@@ocaml.inline always]
 
 val sqrt : t -> t
-[@@ocaml.inline always]
 
 (** Link with integers *)
 val of_int63 : Uint63.t -> t
@@ -77,10 +72,8 @@ val ldshiftexp : t -> Uint63.t -> t
 [@@ocaml.inline always]
 
 val next_up : t -> t
-[@@ocaml.inline always]
 
 val next_down : t -> t
-[@@ocaml.inline always]
 
 (** Return true if two floats are equal.
  * All NaN values are considered equal. *)
diff --git a/kernel/uint63.mli b/kernel/uint63.mli
index c7d1e36451..7ed3d415e4 100644
--- a/kernel/uint63.mli
+++ b/kernel/uint63.mli
@@ -20,6 +20,7 @@ val of_int64 : Int64.t -> t
 (*
 val of_uint : int -> t
 *)
+val to_int_saturate : t -> int (* maxuint31 in case of overflow *)
 
       (* conversion to float *)
 val of_float : float -> t
diff --git a/kernel/uint63_31.ml b/kernel/uint63_31.ml
index 76d768e20a..a5646e87c3 100644
--- a/kernel/uint63_31.ml
+++ b/kernel/uint63_31.ml
@@ -27,6 +27,10 @@ let of_int i = Int64.of_int i
 let to_int2 i = (Int64.to_int (Int64.shift_right_logical i 31), Int64.to_int i)
 let of_int64 i = i
 
+let to_int_saturate i =
+  let r = if Int64.(equal (logand i maxuint31) i) then i else maxuint31 in
+  Int64.to_int r
+
 let of_float f = mask63 (Int64.of_float f)
 let to_float = Int64.to_float
 
@@ -217,4 +221,8 @@ let () =
   Callback.register "uint63 one" one;
   Callback.register "uint63 sub" sub;
   Callback.register "uint63 subcarry" subcarry;
-  Callback.register "uint63 tail0" tail0
+  Callback.register "uint63 tail0" tail0;
+  Callback.register "uint63 of_float" of_float;
+  Callback.register "uint63 to_float" to_float;
+  Callback.register "uint63 of_int" of_int;
+  Callback.register "uint63 to_int_saturate" to_int_saturate
diff --git a/kernel/uint63_63.ml b/kernel/uint63_63.ml
index 4c9377b628..1776f904d6 100644
--- a/kernel/uint63_63.ml
+++ b/kernel/uint63_63.ml
@@ -27,8 +27,13 @@ let to_int2 i = (0,i)
 
 let of_int64 _i = assert false
 
+let to_int_saturate i = i
+
 let of_float = int_of_float
-let to_float i = Int64.to_float (to_uint64 i)
+
+external to_float : int -> (float [@unboxed])
+  = "coq_uint63_to_float_byte" "coq_uint63_to_float"
+[@@noalloc]
 
 let hash i = i
 [@@ocaml.inline always]
diff --git a/test-suite/primitive/float/frexp.v b/test-suite/primitive/float/frexp.v
index 2a600429b1..f13d5cebf6 100644
--- a/test-suite/primitive/float/frexp.v
+++ b/test-suite/primitive/float/frexp.v
@@ -1,3 +1,5 @@
+(* TODO add tests for ldexp (particularly with overflow with 31 and 63 bits integers) *)
+
 Require Import ZArith Floats.
 
 Definition denorm := Eval compute in ldexp one (-1074)%Z.