Add primitive floats to 'vm_compute'

* This commit add float instructions to the VM, their encoding in bytecode
  and the interpretation of primitive float values after the reduction.

* The flag '-std=c99' could be added to the C compiler flags to ensure
  that float computation strictly follows the norm (ie. i387 80-bits
  format is not used as an optimization).

  Actually, we use '-fexcess-precision=standard' instead of '-std=c99'
  because the latter would disable GNU asm used in the VM.

4 files changed, 158 insertions, 24 deletions

diff --git a/kernel/byterun/coq_fix_code.c b/kernel/byterun/coq_fix_code.c
index 0865487c98..bca2cc3bd9 100644
--- a/kernel/byterun/coq_fix_code.c
+++ b/kernel/byterun/coq_fix_code.c
@@ -63,7 +63,12 @@ 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[CHECKOPPFLOAT]=arity[CHECKABSFLOAT]=arity[CHECKCOMPAREFLOAT]=
+    arity[CHECKADDFLOAT]=arity[CHECKSUBFLOAT]=arity[CHECKMULFLOAT]=
+    arity[CHECKDIVFLOAT]=arity[CHECKSQRTFLOAT]=
+    arity[CHECKFLOATOFINT63]=arity[CHECKFLOATNORMFRMANTISSA]=
+    arity[CHECKFRSHIFTEXP]=arity[CHECKLDSHIFTEXP]=1;
   /* instruction with two operands */
   arity[APPTERM]=arity[MAKEBLOCK]=arity[CLOSURE]=
   arity[PROJ]=2;
diff --git a/kernel/byterun/coq_interp.c b/kernel/byterun/coq_interp.c
index 4b45608ae5..55b973dcdb 100644
--- a/kernel/byterun/coq_interp.c
+++ b/kernel/byterun/coq_interp.c
@@ -17,6 +17,7 @@
 #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"
@@ -167,38 +168,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 +1530,128 @@ 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 (CHECKCOMPAREFLOAT) {
+        double x, y;
+        print_instr("CHECKCOMPAREFLOAT");
+        CheckFloat2();
+        x = Double_val(accu);
+        y = Double_val(*sp++);
+        if(x < y) {
+          Alloc_small(accu, 1, coq_tag_Some);
+          Field(accu, 0) = coq_Lt;
+        }
+        else if(x > y) {
+          Alloc_small(accu, 1, coq_tag_Some);
+          Field(accu, 0) = coq_Gt;
+        }
+        else if(x == y) {
+          Alloc_small(accu, 1, coq_tag_Some);
+          Field(accu, 0) = coq_Eq;
+        }
+        else { // nan value
+          accu = coq_None;
+        }
+        Next;
+      }
+
+      Instruct (CHECKADDFLOAT) {
+        print_instr("CHECKADDFLOAT");
+        CheckFloat2();
+        Coq_copy_double(Double_val(accu) + Double_val(*sp++));
+        Next;
+      }
+
+      Instruct (CHECKSUBFLOAT) {
+        print_instr("CHECKSUBFLOAT");
+        CheckFloat2();
+        Coq_copy_double(Double_val(accu) - Double_val(*sp++));
+        Next;
+      }
+
+      Instruct (CHECKMULFLOAT) {
+        print_instr("CHECKMULFLOAT");
+        CheckFloat2();
+        Coq_copy_double(Double_val(accu) * Double_val(*sp++));
+        Next;
+      }
+
+      Instruct (CHECKDIVFLOAT) {
+        print_instr("CHECKDIVFLOAT");
+        CheckFloat2();
+        Coq_copy_double(Double_val(accu) / Double_val(*sp++));
+        Next;
+      }
+
+      Instruct (CHECKSQRTFLOAT) {
+        print_instr("CHECKSQRTFLOAT");
+        CheckFloat1();
+        Coq_copy_double(sqrt(Double_val(accu)));
+        Next;
+      }
+
+      Instruct (CHECKFLOATOFINT63) {
+        print_instr("CHECKFLOATOFINT63");
+        CheckInt1();
+        Coq_copy_double(uint63_to_double(accu));
+        Next;
+      }
+
+      Instruct (CHECKFLOATNORMFRMANTISSA) {
+        double f;
+        print_instr("CHECKFLOATNORMFRMANTISSA");
+        CheckFloat1();
+        f = fabs(Double_val(accu));
+        if (f >= 0.5 && f < 1) {
+          accu = uint63_of_double(ldexp(f, DBL_MANT_DIG));
+        }
+        else {
+          accu = Val_int(0);
+        }
+        Next;
+      }
+
+      Instruct (CHECKFRSHIFTEXP) {
+        int exp;
+        double f;
+        print_instr("CHECKFRSHIFTEXP");
+        CheckFloat1();
+        f = frexp(Double_val(accu), &exp);
+        if (fpclassify(f) == FP_NORMAL) {
+          exp += FLOAT_EXP_SHIFT;
+        }
+        else {
+          exp = 0;
+        }
+        Coq_copy_double(f);
+        *--sp = accu;
+        Alloc_small(accu, 2, coq_tag_pair);
+        Field(accu, 0) = *sp++;
+        Field(accu, 1) = Val_int(exp);
+        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));
+        Next;
+      }
+
 /* Debugging and machine control */
 
       Instruct(STOP){
diff --git a/kernel/byterun/coq_uint63_native.h b/kernel/byterun/coq_uint63_native.h
index 9fbd3f83d8..a14ed5c262 100644
--- a/kernel/byterun/coq_uint63_native.h
+++ b/kernel/byterun/coq_uint63_native.h
@@ -138,3 +138,6 @@ 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))
diff --git a/kernel/byterun/coq_values.h b/kernel/byterun/coq_values.h
index 0cf6ccf532..14f3f152bf 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,9 @@
 #define coq_Eq Val_int(0)
 #define coq_Lt Val_int(1)
 #define coq_Gt Val_int(2)
+#define coq_tag_Some 1
+#define coq_None Val_int(0)
+
+#define FLOAT_EXP_SHIFT (1022 + 52)
 
 #endif /* _COQ_VALUES_ */