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.

13 files changed, 211 insertions, 28 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_ */
diff --git a/kernel/cemitcodes.ml b/kernel/cemitcodes.ml
index 181211d237..82dd7bd85d 100644
--- a/kernel/cemitcodes.ml
+++ b/kernel/cemitcodes.ml
@@ -234,7 +234,18 @@ let check_prim_op = function
   | Int63lt        -> opCHECKLTINT63
   | Int63le        -> opCHECKLEINT63
   | Int63compare   -> opCHECKCOMPAREINT63
-  | _ -> 0 (* TODO: BERTHOLON add float64 operations *)
+  | Float64opp     -> opCHECKOPPFLOAT
+  | Float64abs     -> opCHECKABSFLOAT
+  | Float64compare -> opCHECKCOMPAREFLOAT
+  | Float64add     -> opCHECKADDFLOAT
+  | Float64sub     -> opCHECKSUBFLOAT
+  | Float64mul     -> opCHECKMULFLOAT
+  | Float64div     -> opCHECKDIVFLOAT
+  | Float64sqrt    -> opCHECKSQRTFLOAT
+  | Float64ofInt63 -> opCHECKFLOATOFINT63
+  | Float64normfr_mantissa -> opCHECKFLOATNORMFRMANTISSA
+  | Float64frshiftexp -> opCHECKFRSHIFTEXP
+  | Float64ldshiftexp -> opCHECKLDSHIFTEXP
 
 let emit_instr env = function
   | Klabel lbl -> define_label env lbl
diff --git a/kernel/csymtable.ml b/kernel/csymtable.ml
index 6c9e73b50d..cbffdc731e 100644
--- a/kernel/csymtable.ml
+++ b/kernel/csymtable.ml
@@ -184,7 +184,16 @@ and eval_to_patch env (buff,pl,fv) =
     | Reloc_proj_name p -> slot_for_proj_name p
   in
   let tc = patch buff pl slots in
-  let vm_env = Array.map (slot_for_fv env) fv in
+  let vm_env =
+    (* Beware, this may look like a call to [Array.map], but it's not.
+       Calling [Array.map f] when the first argument returned by [f]
+       is a float would lead to [vm_env] being an unboxed Double_array
+       (Tag_val = Double_array_tag) whereas eval_tcode expects a
+       regular array (Tag_val = 0).
+       See test-suite/primitive/float/coq_env_double_array.v
+       for an actual instance. *)
+    let a = Array.make (Array.length fv) crazy_val in
+    Array.iteri (fun i v -> a.(i) <- slot_for_fv env v) fv; a in
   eval_tcode tc (get_atom_rel ()) (vm_global global_data.glob_val) vm_env
 
 and val_of_constr env c =
diff --git a/kernel/float64.ml b/kernel/float64.ml
index e74fd2e9f1..0b22e07e82 100644
--- a/kernel/float64.ml
+++ b/kernel/float64.ml
@@ -19,6 +19,8 @@ let is_nan f = f <> f
 let to_string f = if is_nan f then "nan" else string_of_float f
 let of_string = float_of_string
 
+let of_float f = f
+
 let opp = ( ~-. )
 let abs = abs_float
 
diff --git a/kernel/float64.mli b/kernel/float64.mli
index fd84f9e61d..7ced535dc0 100644
--- a/kernel/float64.mli
+++ b/kernel/float64.mli
@@ -19,6 +19,8 @@ val is_nan : t -> bool
 val to_string : t -> string
 val of_string : string -> t
 
+val of_float : float -> t
+
 val opp : t -> t
 val abs : t -> t
 
diff --git a/kernel/genOpcodeFiles.ml b/kernel/genOpcodeFiles.ml
index a8a4ffce9c..7deffd030b 100644
--- a/kernel/genOpcodeFiles.ml
+++ b/kernel/genOpcodeFiles.ml
@@ -137,6 +137,18 @@ let opcodes =
     "CHECKTAIL0INT63";
     "ISINT";
     "AREINT2";
+    "CHECKOPPFLOAT";
+    "CHECKABSFLOAT";
+    "CHECKCOMPAREFLOAT";
+    "CHECKADDFLOAT";
+    "CHECKSUBFLOAT";
+    "CHECKMULFLOAT";
+    "CHECKDIVFLOAT";
+    "CHECKSQRTFLOAT";
+    "CHECKFLOATOFINT63";
+    "CHECKFLOATNORMFRMANTISSA";
+    "CHECKFRSHIFTEXP";
+    "CHECKLDSHIFTEXP";
     "STOP"
   |]
 
diff --git a/kernel/vconv.ml b/kernel/vconv.ml
index 414c443c4e..5d36ad54a2 100644
--- a/kernel/vconv.ml
+++ b/kernel/vconv.ml
@@ -73,6 +73,9 @@ and conv_whd env pb k whd1 whd2 cu =
       else raise NotConvertible
   | Vint64 i1, Vint64 i2 ->
     if Int64.equal i1 i2 then cu else raise NotConvertible
+  | Vfloat64 f1, Vfloat64 f2 ->
+    if Float64.(equal (of_float f1) (of_float f2)) then cu
+    else raise NotConvertible
   | Vatom_stk(a1,stk1), Vatom_stk(a2,stk2) ->
       conv_atom env pb k a1 stk1 a2 stk2 cu
   | Vfun _, _ | _, Vfun _ ->
@@ -80,7 +83,7 @@ and conv_whd env pb k whd1 whd2 cu =
       conv_val env CONV (k+1) (apply_whd k whd1) (apply_whd k whd2) cu
 
   | Vprod _, _ | Vfix _, _ | Vcofix _, _  | Vconstr_const _, _ | Vint64 _, _
-  | Vconstr_block _, _ | Vatom_stk _, _ -> raise NotConvertible
+  | Vfloat64 _, _ | Vconstr_block _, _ | Vatom_stk _, _ -> raise NotConvertible
 
 
 and conv_atom env pb k a1 stk1 a2 stk2 cu =
diff --git a/kernel/vm.ml b/kernel/vm.ml
index 319a26d824..5f08720f77 100644
--- a/kernel/vm.ml
+++ b/kernel/vm.ml
@@ -169,7 +169,8 @@ let rec apply_stack a stk v =
 let apply_whd k whd =
   let v = val_of_rel k in
   match whd with
-  | Vprod _ | Vconstr_const _ | Vconstr_block _ | Vint64 _ -> assert false
+  | Vprod _ | Vconstr_const _ | Vconstr_block _ | Vint64 _ | Vfloat64 _ ->
+     assert false
   | Vfun f -> reduce_fun k f
   | Vfix(f, None) -> 
       push_ra stop;
diff --git a/kernel/vmvalues.ml b/kernel/vmvalues.ml
index fe3c76c960..5acdd964b1 100644
--- a/kernel/vmvalues.ml
+++ b/kernel/vmvalues.ml
@@ -76,6 +76,8 @@ let rec eq_structured_values v1 v2 =
        Int.equal (Obj.size o1) (Obj.size o2)
     then if Int.equal t1 Obj.custom_tag
       then Int64.equal (Obj.magic v1 : int64) (Obj.magic v2 : int64)
+    else if Int.equal t1 Obj.double_tag
+      then Float64.(equal (of_float (Obj.magic v1)) (of_float (Obj.magic v2)))
     else begin
       assert (t1 <= Obj.last_non_constant_constructor_tag &&
               t2 <= Obj.last_non_constant_constructor_tag);
@@ -289,6 +291,7 @@ type whd =
   | Vconstr_const of int
   | Vconstr_block of vblock
   | Vint64 of int64
+  | Vfloat64 of float
   | Vatom_stk of atom * stack
   | Vuniv_level of Univ.Level.t
 
@@ -320,6 +323,7 @@ let uni_lvl_val (v : values) : Univ.Level.t =
         | Vconstr_const _i -> str "Vconstr_const"
         | Vconstr_block _b -> str "Vconstr_block"
         | Vint64 _ -> str "Vint64"
+        | Vfloat64 _ -> str "Vfloat64"
         | Vatom_stk (_a,_stk) -> str "Vatom_stk"
         | Vuniv_level _ -> assert false
       in
@@ -379,6 +383,8 @@ let rec whd_accu a stk =
       end
   | i when Int.equal i Obj.custom_tag ->
     Vint64 (Obj.magic i)
+  | i when Int.equal i Obj.double_tag ->
+    Vfloat64 (Obj.magic i)
   | tg ->
     CErrors.anomaly
       Pp.(strbrk "Failed to parse VM value. Tag = " ++ int tg ++ str ".")
@@ -408,6 +414,7 @@ let whd_val : values -> whd =
            | 3 -> Vatom_stk(Aid(RelKey(int_tcode (fun_code o) 1)), [])
            | _ -> CErrors.anomaly ~label:"Vm.whd " (Pp.str "kind_of_closure does not work."))
         else if Int.equal tag Obj.custom_tag then Vint64 (Obj.magic v)
+        else if Int.equal tag Obj.double_tag then Vfloat64 (Obj.magic v)
         else
            Vconstr_block(Obj.obj o)
 
@@ -681,6 +688,7 @@ and pr_whd w =
   | Vconstr_const i -> str "Vconstr_const(" ++ int i ++ str ")"
   | Vconstr_block _b -> str "Vconstr_block"
   | Vint64 i -> i |> Format.sprintf "Vint64(%LiL)" |> str
+  | Vfloat64 f -> str "Vfloat64(" ++ str (Float64.(to_string (of_float f))) ++ str ")"
   | Vatom_stk (a,stk) -> str "Vatom_stk(" ++ pr_atom a ++ str ", " ++ pr_stack stk ++ str ")"
   | Vuniv_level _ -> assert false)
 and pr_stack stk =
diff --git a/kernel/vmvalues.mli b/kernel/vmvalues.mli
index 1e40801be0..9c24006ff0 100644
--- a/kernel/vmvalues.mli
+++ b/kernel/vmvalues.mli
@@ -128,6 +128,7 @@ type whd =
   | Vconstr_const of int
   | Vconstr_block of vblock
   | Vint64 of int64
+  | Vfloat64 of float
   | Vatom_stk of atom * stack
   | Vuniv_level of Univ.Level.t