Merge PR #9867: Add primitive floats (binary64 floating-point numbers)

Ack-by: SkySkimmer
Reviewed-by: Zimmi48
Ack-by: ejgallego
Reviewed-by: maximedenes
Ack-by: proux01
Ack-by: silene
Ack-by: vbgl

7 files changed, 348 insertions, 24 deletions

diff --git a/kernel/byterun/coq_fix_code.c b/kernel/byterun/coq_fix_code.c
index 0865487c98..931b509f48 100644
--- a/kernel/byterun/coq_fix_code.c
+++ b/kernel/byterun/coq_fix_code.c
@@ -44,6 +44,7 @@ void init_arity () {
     arity[PUSHCONST0]=arity[PUSHCONST1]=arity[PUSHCONST2]=arity[PUSHCONST3]=
     arity[ACCUMULATE]=arity[STOP]=arity[MAKEPROD]= 
     arity[ADDINT63]=arity[SUBINT63]=arity[LTINT63]=arity[LEINT63]=
+    arity[LTFLOAT]=arity[LEFLOAT]=
     arity[ISINT]=arity[AREINT2]=0;
   /* instruction with one operand */
   arity[ACC]=arity[PUSHACC]=arity[POP]=arity[ENVACC]=arity[PUSHENVACC]=
@@ -63,7 +64,15 @@ void init_arity () {
     arity[CHECKLSLINT63]=arity[CHECKLSRINT63]=arity[CHECKADDMULDIVINT63]=
     arity[CHECKLSLINT63CONST1]=arity[CHECKLSRINT63CONST1]=
     arity[CHECKEQINT63]=arity[CHECKLTINT63]=arity[CHECKLEINT63]=
-    arity[CHECKCOMPAREINT63]=arity[CHECKHEAD0INT63]=arity[CHECKTAIL0INT63]=1;
+    arity[CHECKCOMPAREINT63]=arity[CHECKHEAD0INT63]=arity[CHECKTAIL0INT63]=
+    arity[CHECKEQFLOAT]=arity[CHECKLTFLOAT]=arity[CHECKLEFLOAT]=
+    arity[CHECKOPPFLOAT]=arity[CHECKABSFLOAT]=arity[CHECKCOMPAREFLOAT]=
+    arity[CHECKCLASSIFYFLOAT]=
+    arity[CHECKADDFLOAT]=arity[CHECKSUBFLOAT]=arity[CHECKMULFLOAT]=
+    arity[CHECKDIVFLOAT]=arity[CHECKSQRTFLOAT]=
+    arity[CHECKFLOATOFINT63]=arity[CHECKFLOATNORMFRMANTISSA]=
+    arity[CHECKFRSHIFTEXP]=arity[CHECKLDSHIFTEXP]=
+    arity[CHECKNEXTUPFLOAT]=arity[CHECKNEXTDOWNFLOAT]=1;
   /* instruction with two operands */
   arity[APPTERM]=arity[MAKEBLOCK]=arity[CLOSURE]=
   arity[PROJ]=2;
diff --git a/kernel/byterun/coq_float64.h b/kernel/byterun/coq_float64.h
new file mode 100644
index 0000000000..c41079c8ff
--- /dev/null
+++ b/kernel/byterun/coq_float64.h
@@ -0,0 +1,48 @@
+#ifndef _COQ_FLOAT64_
+#define _COQ_FLOAT64_
+
+#include <math.h>
+
+#define DECLARE_FREL(name, e)                                           \
+  int coq_##name(double x, double y) {                                  \
+    return e;                                                           \
+  }                                                                     \
+                                                                        \
+  value coq_##name##_byte(value x, value y) {                           \
+    return coq_##name(Double_val(x), Double_val(y));                    \
+  }
+
+#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_FREL(feq, x == y)
+DECLARE_FREL(flt, x < y)
+DECLARE_FREL(fle, x <= y)
+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))
+
+value coq_is_double(value x) {
+  return Val_long(Is_double(x));
+}
+
+#endif /* _COQ_FLOAT64_ */
diff --git a/kernel/byterun/coq_interp.c b/kernel/byterun/coq_interp.c
index 4b45608ae5..ca1308696c 100644
--- a/kernel/byterun/coq_interp.c
+++ b/kernel/byterun/coq_interp.c
@@ -17,11 +17,13 @@
 #include <signal.h>
 #include <stdint.h>
 #include <caml/memory.h>
+#include <math.h>
 #include "coq_gc.h"
 #include "coq_instruct.h"
 #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"
@@ -167,38 +169,34 @@ if (sp - num_args < coq_stack_threshold) {                                     \
 #endif
 #endif
 
-#define CheckInt1() do{                            \
-    if (Is_uint63(accu)) pc++;			   \
+#define CheckPrimArgs(cond, apply_lbl) do{         \
+    if (cond) pc++;                                \
     else{					   \
       *--sp=accu;				   \
       accu = Field(coq_global_data, *pc++);	   \
-      goto apply1;				   \
-    }						   \
-  }while(0)
-
-#define CheckInt2() do{    			   \
-    if (Is_uint63(accu) && Is_uint63(sp[0])) pc++;	   \
-    else{					   \
-      *--sp=accu;				   \
-      accu = Field(coq_global_data, *pc++);	   \
-      goto apply2;				   \
+      goto apply_lbl;                              \
     }						   \
   }while(0)
 
-
-
-#define CheckInt3() do{						      \
-    if (Is_uint63(accu) && Is_uint63(sp[0]) && Is_uint63(sp[1]) ) pc++;     \
-    else{							      \
-      *--sp=accu;						      \
-      accu = Field(coq_global_data, *pc++);			      \
-      goto apply3;						      \
-    }								      \
-  }while(0)
+#define CheckInt1() CheckPrimArgs(Is_uint63(accu), apply1)
+#define CheckInt2() CheckPrimArgs(Is_uint63(accu) && Is_uint63(sp[0]), apply2)
+#define CheckInt3() CheckPrimArgs(Is_uint63(accu) && Is_uint63(sp[0]) \
+                                                  && Is_uint63(sp[1]), apply3)
+#define CheckFloat1() CheckPrimArgs(Is_double(accu), apply1)
+#define CheckFloat2() CheckPrimArgs(Is_double(accu) && Is_double(sp[0]), apply2)
 
 #define AllocCarry(cond) Alloc_small(accu, 1, (cond)? coq_tag_C1 : coq_tag_C0)
 #define AllocPair() Alloc_small(accu, 2, coq_tag_pair)
 
+/* Beware: we cannot use caml_copy_double here as it doesn't use
+   Alloc_small, hence doesn't protect the stack via
+   Setup_for_gc/Restore_after_gc. */
+#define Coq_copy_double(val) do{                   \
+  double Coq_copy_double_f__ = (val);              \
+  Alloc_small(accu, Double_wosize, Double_tag);    \
+  Store_double_val(accu, Coq_copy_double_f__);     \
+  }while(0);
+
 #define Swap_accu_sp do{                        \
     value swap_accu_sp_tmp__ = accu;            \
     accu = *sp;                                 \
@@ -1533,6 +1531,206 @@ value coq_interprete
       }
 
 
+      Instruct (CHECKOPPFLOAT) {
+        print_instr("CHECKOPPFLOAT");
+        CheckFloat1();
+        Coq_copy_double(-Double_val(accu));
+        Next;
+      }
+
+      Instruct (CHECKABSFLOAT) {
+        print_instr("CHECKABSFLOAT");
+        CheckFloat1();
+        Coq_copy_double(fabs(Double_val(accu)));
+        Next;
+      }
+
+      Instruct (CHECKEQFLOAT) {
+        print_instr("CHECKEQFLOAT");
+        CheckFloat2();
+        accu = coq_feq(Double_val(accu), Double_val(*sp++)) ? coq_true : coq_false;
+        Next;
+      }
+
+      Instruct (CHECKLTFLOAT) {
+        print_instr("CHECKLTFLOAT");
+        CheckFloat2();
+      }
+      Instruct (LTFLOAT) {
+        print_instr("LTFLOAT");
+        accu = coq_flt(Double_val(accu), Double_val(*sp++)) ? coq_true : coq_false;
+        Next;
+      }
+
+      Instruct (CHECKLEFLOAT) {
+        print_instr("CHECKLEFLOAT");
+        CheckFloat2();
+      }
+      Instruct (LEFLOAT) {
+        print_instr("LEFLOAT");
+        accu = coq_fle(Double_val(accu), Double_val(*sp++)) ? coq_true : coq_false;
+        Next;
+      }
+
+      Instruct (CHECKCOMPAREFLOAT) {
+        double x, y;
+        print_instr("CHECKCOMPAREFLOAT");
+        CheckFloat2();
+        x = Double_val(accu);
+        y = Double_val(*sp++);
+        if(x < y) {
+          accu = coq_FLt;
+        }
+        else if(x > y) {
+          accu = coq_FGt;
+        }
+        else if(x == y) {
+          accu = coq_FEq;
+        }
+        else { // nan value
+          accu = coq_FNotComparable;
+        }
+        Next;
+      }
+
+      Instruct (CHECKCLASSIFYFLOAT) {
+        double x;
+        print_instr("CHECKCLASSIFYFLOAT");
+        CheckFloat1();
+        x = Double_val(accu);
+        switch (fpclassify(x)) {
+        case FP_NORMAL:
+          accu = signbit(x) ? coq_NNormal : coq_PNormal;
+          break;
+        case FP_SUBNORMAL:
+          accu = signbit(x) ? coq_NSubn : coq_PSubn;
+          break;
+        case FP_ZERO:
+          accu = signbit(x) ? coq_NZero : coq_PZero;
+          break;
+        case FP_INFINITE:
+          accu = signbit(x) ? coq_NInf : coq_PInf;
+          break;
+        default:  /* FP_NAN */
+          accu = coq_NaN;
+          break;
+        }
+        Next;
+      }
+
+      Instruct (CHECKADDFLOAT) {
+        print_instr("CHECKADDFLOAT");
+        CheckFloat2();
+        Coq_copy_double(coq_fadd(Double_val(accu), Double_val(*sp++)));
+        Next;
+      }
+
+      Instruct (CHECKSUBFLOAT) {
+        print_instr("CHECKSUBFLOAT");
+        CheckFloat2();
+        Coq_copy_double(coq_fsub(Double_val(accu), Double_val(*sp++)));
+        Next;
+      }
+
+      Instruct (CHECKMULFLOAT) {
+        print_instr("CHECKMULFLOAT");
+        CheckFloat2();
+        Coq_copy_double(coq_fmul(Double_val(accu), Double_val(*sp++)));
+        Next;
+      }
+
+      Instruct (CHECKDIVFLOAT) {
+        print_instr("CHECKDIVFLOAT");
+        CheckFloat2();
+        Coq_copy_double(coq_fdiv(Double_val(accu), Double_val(*sp++)));
+        Next;
+      }
+
+      Instruct (CHECKSQRTFLOAT) {
+        print_instr("CHECKSQRTFLOAT");
+        CheckFloat1();
+        Coq_copy_double(coq_fsqrt(Double_val(accu)));
+        Next;
+      }
+
+      Instruct (CHECKFLOATOFINT63) {
+        print_instr("CHECKFLOATOFINT63");
+        CheckInt1();
+        Uint63_to_double(accu);
+        Next;
+      }
+
+      Instruct (CHECKFLOATNORMFRMANTISSA) {
+        double f;
+        print_instr("CHECKFLOATNORMFRMANTISSA");
+        CheckFloat1();
+        f = fabs(Double_val(accu));
+        if (f >= 0.5 && f < 1) {
+          Uint63_of_double(ldexp(f, DBL_MANT_DIG));
+        }
+        else {
+          Uint63_of_int(Val_int(0));
+        }
+        Next;
+      }
+
+      Instruct (CHECKFRSHIFTEXP) {
+        int exp;
+        double f;
+        print_instr("CHECKFRSHIFTEXP");
+        CheckFloat1();
+        /* frexp(infinity) incorrectly returns nan on mingw */
+#if defined(__MINGW32__) || defined(__MINGW64__)
+        if (fpclassify(Double_val(accu)) == FP_INFINITE) {
+          f = Double_val(accu);
+        } else
+#endif
+        f = frexp(Double_val(accu), &exp);
+        if (fpclassify(f) == FP_NORMAL) {
+          exp += FLOAT_EXP_SHIFT;
+        }
+        else {
+          exp = 0;
+        }
+        Coq_copy_double(f);
+        *--sp = accu;
+#ifdef ARCH_SIXTYFOUR
+        Alloc_small(accu, 2, coq_tag_pair);
+        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);
+        Swap_accu_sp;
+        Uint63_to_int_min(accu, Val_int(2 * FLOAT_EXP_SHIFT));
+        accu = Int_val(accu);
+        Coq_copy_double(ldexp(Double_val(*sp++), accu - FLOAT_EXP_SHIFT));
+        Next;
+      }
+
+      Instruct (CHECKNEXTUPFLOAT) {
+        print_instr("CHECKNEXTUPFLOAT");
+        CheckFloat1();
+        Coq_copy_double(coq_next_up(Double_val(accu)));
+        Next;
+      }
+
+      Instruct (CHECKNEXTDOWNFLOAT) {
+        print_instr("CHECKNEXTDOWNFLOAT");
+        CheckFloat1();
+        Coq_copy_double(coq_next_down(Double_val(accu)));
+        Next;
+      }
+
 /* Debugging and machine control */
 
       Instruct(STOP){
diff --git a/kernel/byterun/coq_uint63_emul.h b/kernel/byterun/coq_uint63_emul.h
index 528cc6fc1f..143a6d098c 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_BINOP(to_int_min)
+#define Uint63_to_int_min(n, m) CALL_BINOP(to_int_min, n, m)
diff --git a/kernel/byterun/coq_uint63_native.h b/kernel/byterun/coq_uint63_native.h
index 9fbd3f83d8..5be7587091 100644
--- a/kernel/byterun/coq_uint63_native.h
+++ b/kernel/byterun/coq_uint63_native.h
@@ -138,3 +138,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(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 Uint63_to_int_min(n, m) do { \
+  if (uint63_lt((n),(m))) \
+    accu = (n); \
+  else \
+    accu = (m); \
+  }while(0)
diff --git a/kernel/byterun/coq_values.h b/kernel/byterun/coq_values.h
index 0cf6ccf532..b027673ac7 100644
--- a/kernel/byterun/coq_values.h
+++ b/kernel/byterun/coq_values.h
@@ -14,6 +14,8 @@
 #include <caml/alloc.h>
 #include <caml/mlvalues.h>
 
+#include <float.h>
+
 #define Default_tag 0
 #define Accu_tag 0
 
@@ -29,8 +31,9 @@
 /* Les blocs accumulate */
 #define Is_accu(v) (Is_block(v) && (Tag_val(v) == Accu_tag))
 #define IS_EVALUATED_COFIX(v) (Is_accu(v) && Is_block(Field(v,1)) && (Tag_val(Field(v,1)) == ATOM_COFIXEVALUATED_TAG))
+#define Is_double(v) (Tag_val(v) == Double_tag)
 
-/* */
+/* coq values for primitive operations */
 #define coq_tag_C1 2
 #define coq_tag_C0 1
 #define coq_tag_pair 1
@@ -39,5 +42,20 @@
 #define coq_Eq Val_int(0)
 #define coq_Lt Val_int(1)
 #define coq_Gt Val_int(2)
+#define coq_FEq Val_int(0)
+#define coq_FLt Val_int(1)
+#define coq_FGt Val_int(2)
+#define coq_FNotComparable Val_int(3)
+#define coq_PNormal Val_int(0)
+#define coq_NNormal Val_int(1)
+#define coq_PSubn Val_int(2)
+#define coq_NSubn Val_int(3)
+#define coq_PZero Val_int(4)
+#define coq_NZero Val_int(5)
+#define coq_PInf Val_int(6)
+#define coq_NInf Val_int(7)
+#define coq_NaN Val_int(8)
+
+#define FLOAT_EXP_SHIFT (2101) /* 2*emax + prec */
 
 #endif /* _COQ_VALUES_ */
diff --git a/kernel/byterun/dune b/kernel/byterun/dune
index 20bdf28e54..d0145176ea 100644
--- a/kernel/byterun/dune
+++ b/kernel/byterun/dune
@@ -1,3 +1,16 @@
+; Dune doesn't use configure's output, but it is still necessary for
+; some Coq files to work; will be fixed in the future.
+(rule
+ (targets dune.c_flags)
+ (mode fallback)
+ (deps %{project_root}/configure.ml %{project_root}/dev/ocamldebug-coq.run (env_var COQ_CONFIGURE_PREFIX))
+ (action (chdir %{project_root} (run %{ocaml} configure.ml -no-ask -native-compiler no))))
+
+(env
+ (dev     (c_flags (:include dune.c_flags)))
+ (release (c_flags (:include dune.c_flags)))
+ (ireport (c_flags (:include dune.c_flags))))
+
 (library
  (name byterun)
  (synopsis "Coq's Kernel Abstract Reduction Machine [C implementation]")