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.

4 files changed, 92 insertions, 15 deletions

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)