Add next_{up,down} primitive float functions

9 files changed, 55 insertions, 2 deletions

diff --git a/kernel/byterun/coq_fix_code.c b/kernel/byterun/coq_fix_code.c
index fb39ca8358..3fe77afc2d 100644
--- a/kernel/byterun/coq_fix_code.c
+++ b/kernel/byterun/coq_fix_code.c
@@ -69,7 +69,8 @@ void init_arity () {
     arity[CHECKADDFLOAT]=arity[CHECKSUBFLOAT]=arity[CHECKMULFLOAT]=
     arity[CHECKDIVFLOAT]=arity[CHECKSQRTFLOAT]=
     arity[CHECKFLOATOFINT63]=arity[CHECKFLOATNORMFRMANTISSA]=
-    arity[CHECKFRSHIFTEXP]=arity[CHECKLDSHIFTEXP]=1;
+    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_interp.c b/kernel/byterun/coq_interp.c
index b862480fda..06042bb753 100644
--- a/kernel/byterun/coq_interp.c
+++ b/kernel/byterun/coq_interp.c
@@ -1674,6 +1674,20 @@ value coq_interprete
         Next;
       }
 
+      Instruct (CHECKNEXTUPFLOAT) {
+        print_instr("CHECKNEXTUPFLOAT");
+        CheckFloat1();
+        Coq_copy_double(nextafter(Double_val(accu), INFINITY));
+        Next;
+      }
+
+      Instruct (CHECKNEXTDOWNFLOAT) {
+        print_instr("CHECKNEXTDOWNFLOAT");
+        CheckFloat1();
+        Coq_copy_double(nextafter(Double_val(accu), -INFINITY));
+        Next;
+      }
+
 /* Debugging and machine control */
 
       Instruct(STOP){
diff --git a/kernel/cPrimitives.ml b/kernel/cPrimitives.ml
index 02a5351ccf..342cc29a22 100644
--- a/kernel/cPrimitives.ml
+++ b/kernel/cPrimitives.ml
@@ -46,6 +46,8 @@ type t =
   | Float64normfr_mantissa
   | Float64frshiftexp
   | Float64ldshiftexp
+  | Float64next_up
+  | Float64next_down
 
 let equal (p1 : t) (p2 : t) =
   p1 == p2
@@ -88,6 +90,8 @@ let hash = function
   | Float64normfr_mantissa -> 35
   | Float64frshiftexp -> 36
   | Float64ldshiftexp -> 37
+  | Float64next_up -> 38
+  | Float64next_down -> 39
 
 (* Should match names in nativevalues.ml *)
 let to_string = function
@@ -128,6 +132,8 @@ let to_string = function
   | Float64normfr_mantissa -> "normfr_mantissa"
   | Float64frshiftexp -> "frshiftexp"
   | Float64ldshiftexp -> "ldshiftexp"
+  | Float64next_up    -> "next_up"
+  | Float64next_down  -> "next_down"
 
 type prim_type =
   | PT_int63
@@ -165,7 +171,8 @@ let types =
   | Int63div21 ->
      [int_ty; int_ty; int_ty; PITT_ind (PIT_pair, (PT_int63, PT_int63))]
   | Int63addMulDiv -> [int_ty; int_ty; int_ty; int_ty]
-  | Float64opp | Float64abs | Float64sqrt -> [float_ty; float_ty]
+  | Float64opp | Float64abs | Float64sqrt
+  | Float64next_up | Float64next_down -> [float_ty; float_ty]
   | Float64ofInt63 -> [int_ty; float_ty]
   | Float64normfr_mantissa -> [float_ty; int_ty]
   | Float64frshiftexp -> [float_ty; PITT_ind (PIT_pair, (PT_float64, PT_int63))]
diff --git a/kernel/cPrimitives.mli b/kernel/cPrimitives.mli
index af95f6c6d7..3cb210233d 100644
--- a/kernel/cPrimitives.mli
+++ b/kernel/cPrimitives.mli
@@ -46,6 +46,8 @@ type t =
   | Float64normfr_mantissa
   | Float64frshiftexp
   | Float64ldshiftexp
+  | Float64next_up
+  | Float64next_down
 
 val equal : t -> t -> bool
 
diff --git a/kernel/cemitcodes.ml b/kernel/cemitcodes.ml
index 535034d8fa..908f84293c 100644
--- a/kernel/cemitcodes.ml
+++ b/kernel/cemitcodes.ml
@@ -247,6 +247,8 @@ let check_prim_op = function
   | Float64normfr_mantissa -> opCHECKFLOATNORMFRMANTISSA
   | Float64frshiftexp -> opCHECKFRSHIFTEXP
   | Float64ldshiftexp -> opCHECKLDSHIFTEXP
+  | Float64next_up    -> opCHECKNEXTUPFLOAT
+  | Float64next_down  -> opCHECKNEXTDOWNFLOAT
 
 let emit_instr env = function
   | Klabel lbl -> define_label env lbl
diff --git a/kernel/float64.ml b/kernel/float64.ml
index 7b54fd0c4b..351661f44d 100644
--- a/kernel/float64.ml
+++ b/kernel/float64.ml
@@ -76,6 +76,24 @@ let frshiftexp f =
 
 let ldshiftexp f e = ldexp f (snd (Uint63.to_int2 e) - eshift)
 
+let eta_float = ldexp 1. (-1074) (* smallest positive float (subnormal) *)
+
+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
+
+let next_down f = -.(next_up (-.f))
+
 let equal f1 f2 =
   match classify_float f1 with
   | FP_normal | FP_subnormal | FP_infinite -> (f1 = f2)
diff --git a/kernel/float64.mli b/kernel/float64.mli
index 1ad980a691..580004126d 100644
--- a/kernel/float64.mli
+++ b/kernel/float64.mli
@@ -49,6 +49,9 @@ val normfr_mantissa : t -> Uint63.t
 val frshiftexp : t -> t * Uint63.t (* float remainder, shifted exponent *)
 val ldshiftexp : t -> Uint63.t -> t
 
+val next_up : t -> t
+val next_down : t -> t
+
 (** Return true if two floats are equal.
  * All NaN values are considered equal. *)
 val equal : t -> t -> bool
diff --git a/kernel/genOpcodeFiles.ml b/kernel/genOpcodeFiles.ml
index 045a1e361d..52b7a822e3 100644
--- a/kernel/genOpcodeFiles.ml
+++ b/kernel/genOpcodeFiles.ml
@@ -150,6 +150,8 @@ let opcodes =
     "CHECKFLOATNORMFRMANTISSA";
     "CHECKFRSHIFTEXP";
     "CHECKLDSHIFTEXP";
+    "CHECKNEXTUPFLOAT";
+    "CHECKNEXTDOWNFLOAT";
     "STOP"
   |]
 
diff --git a/kernel/primred.ml b/kernel/primred.ml
index cfe6c8effe..2766793005 100644
--- a/kernel/primred.ml
+++ b/kernel/primred.ml
@@ -302,6 +302,10 @@ struct
       let f = get_float evd args 0 in
       let e = get_int evd args 1 in
       E.mkFloat env (Float64.ldshiftexp f e)
+    | Float64next_up ->
+      let f = get_float1 evd args in E.mkFloat env (Float64.next_up f)
+    | Float64next_down ->
+      let f = get_float1 evd args in E.mkFloat env (Float64.next_down f)
 
   let red_prim env evd p args =
     try