Add next_{up,down} primitive float functions

15 files changed, 187 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
diff --git a/test-suite/primitive/float/next_up_down.v b/test-suite/primitive/float/next_up_down.v
new file mode 100644
index 0000000000..95f5248bf7
--- /dev/null
+++ b/test-suite/primitive/float/next_up_down.v
@@ -0,0 +1,89 @@
+Require Import ZArith Int63 Floats.
+
+Open Scope float_scope.
+
+Definition f0 := zero.
+Definition f1 := neg_zero.
+Definition f2 := Eval compute in ldexp one 0.
+Definition f3 := Eval compute in -f1.
+(* smallest positive float *)
+Definition f4 := Eval compute in ldexp one (-1074).
+Definition f5 := Eval compute in -f3.
+Definition f6 := infinity.
+Definition f7 := neg_infinity.
+Definition f8 := Eval compute in ldexp one (-1).
+Definition f9 := Eval compute in -f8.
+Definition f10 := Eval compute in of_int63 42.
+Definition f11 := Eval compute in -f10.
+(* max float *)
+Definition f12 := Eval compute in ldexp (of_int63 9007199254740991) 1024.
+Definition f13 := Eval compute in -f12.
+(* smallest positive normalized float *)
+Definition f14 := Eval compute in ldexp one (-1022).
+Definition f15 := Eval compute in -f14.
+
+Check (eq_refl : Prim2SF (next_up f0) = SF64succ (Prim2SF f0)).
+Check (eq_refl : Prim2SF (next_down f0) = SF64pred (Prim2SF f0)).
+Check (eq_refl : Prim2SF (next_up f1) = SF64succ (Prim2SF f1)).
+Check (eq_refl : Prim2SF (next_down f1) = SF64pred (Prim2SF f1)).
+Check (eq_refl : Prim2SF (next_up f2) = SF64succ (Prim2SF f2)).
+Check (eq_refl : Prim2SF (next_down f2) = SF64pred (Prim2SF f2)).
+Check (eq_refl : Prim2SF (next_up f3) = SF64succ (Prim2SF f3)).
+Check (eq_refl : Prim2SF (next_down f3) = SF64pred (Prim2SF f3)).
+Check (eq_refl : Prim2SF (next_up f4) = SF64succ (Prim2SF f4)).
+Check (eq_refl : Prim2SF (next_down f4) = SF64pred (Prim2SF f4)).
+Check (eq_refl : Prim2SF (next_up f5) = SF64succ (Prim2SF f5)).
+Check (eq_refl : Prim2SF (next_down f5) = SF64pred (Prim2SF f5)).
+Check (eq_refl : Prim2SF (next_up f6) = SF64succ (Prim2SF f6)).
+Check (eq_refl : Prim2SF (next_down f6) = SF64pred (Prim2SF f6)).
+Check (eq_refl : Prim2SF (next_up f7) = SF64succ (Prim2SF f7)).
+Check (eq_refl : Prim2SF (next_down f7) = SF64pred (Prim2SF f7)).
+Check (eq_refl : Prim2SF (next_up f8) = SF64succ (Prim2SF f8)).
+Check (eq_refl : Prim2SF (next_down f8) = SF64pred (Prim2SF f8)).
+Check (eq_refl : Prim2SF (next_up f9) = SF64succ (Prim2SF f9)).
+Check (eq_refl : Prim2SF (next_down f9) = SF64pred (Prim2SF f9)).
+Check (eq_refl : Prim2SF (next_up f10) = SF64succ (Prim2SF f10)).
+Check (eq_refl : Prim2SF (next_down f10) = SF64pred (Prim2SF f10)).
+Check (eq_refl : Prim2SF (next_up f11) = SF64succ (Prim2SF f11)).
+Check (eq_refl : Prim2SF (next_down f11) = SF64pred (Prim2SF f11)).
+Check (eq_refl : Prim2SF (next_up f12) = SF64succ (Prim2SF f12)).
+Check (eq_refl : Prim2SF (next_down f12) = SF64pred (Prim2SF f12)).
+Check (eq_refl : Prim2SF (next_up f13) = SF64succ (Prim2SF f13)).
+Check (eq_refl : Prim2SF (next_down f13) = SF64pred (Prim2SF f13)).
+Check (eq_refl : Prim2SF (next_up f14) = SF64succ (Prim2SF f14)).
+Check (eq_refl : Prim2SF (next_down f14) = SF64pred (Prim2SF f14)).
+Check (eq_refl : Prim2SF (next_up f15) = SF64succ (Prim2SF f15)).
+Check (eq_refl : Prim2SF (next_down f15) = SF64pred (Prim2SF f15)).
+
+Check (eq_refl (SF64succ (Prim2SF f0)) <: Prim2SF (next_up f0) = SF64succ (Prim2SF f0)).
+Check (eq_refl (SF64pred (Prim2SF f0)) <: Prim2SF (next_down f0) = SF64pred (Prim2SF f0)).
+Check (eq_refl (SF64succ (Prim2SF f1)) <: Prim2SF (next_up f1) = SF64succ (Prim2SF f1)).
+Check (eq_refl (SF64pred (Prim2SF f1)) <: Prim2SF (next_down f1) = SF64pred (Prim2SF f1)).
+Check (eq_refl (SF64succ (Prim2SF f2)) <: Prim2SF (next_up f2) = SF64succ (Prim2SF f2)).
+Check (eq_refl (SF64pred (Prim2SF f2)) <: Prim2SF (next_down f2) = SF64pred (Prim2SF f2)).
+Check (eq_refl (SF64succ (Prim2SF f3)) <: Prim2SF (next_up f3) = SF64succ (Prim2SF f3)).
+Check (eq_refl (SF64pred (Prim2SF f3)) <: Prim2SF (next_down f3) = SF64pred (Prim2SF f3)).
+Check (eq_refl (SF64succ (Prim2SF f4)) <: Prim2SF (next_up f4) = SF64succ (Prim2SF f4)).
+Check (eq_refl (SF64pred (Prim2SF f4)) <: Prim2SF (next_down f4) = SF64pred (Prim2SF f4)).
+Check (eq_refl (SF64succ (Prim2SF f5)) <: Prim2SF (next_up f5) = SF64succ (Prim2SF f5)).
+Check (eq_refl (SF64pred (Prim2SF f5)) <: Prim2SF (next_down f5) = SF64pred (Prim2SF f5)).
+Check (eq_refl (SF64succ (Prim2SF f6)) <: Prim2SF (next_up f6) = SF64succ (Prim2SF f6)).
+Check (eq_refl (SF64pred (Prim2SF f6)) <: Prim2SF (next_down f6) = SF64pred (Prim2SF f6)).
+Check (eq_refl (SF64succ (Prim2SF f7)) <: Prim2SF (next_up f7) = SF64succ (Prim2SF f7)).
+Check (eq_refl (SF64pred (Prim2SF f7)) <: Prim2SF (next_down f7) = SF64pred (Prim2SF f7)).
+Check (eq_refl (SF64succ (Prim2SF f8)) <: Prim2SF (next_up f8) = SF64succ (Prim2SF f8)).
+Check (eq_refl (SF64pred (Prim2SF f8)) <: Prim2SF (next_down f8) = SF64pred (Prim2SF f8)).
+Check (eq_refl (SF64succ (Prim2SF f9)) <: Prim2SF (next_up f9) = SF64succ (Prim2SF f9)).
+Check (eq_refl (SF64pred (Prim2SF f9)) <: Prim2SF (next_down f9) = SF64pred (Prim2SF f9)).
+Check (eq_refl (SF64succ (Prim2SF f10)) <: Prim2SF (next_up f10) = SF64succ (Prim2SF f10)).
+Check (eq_refl (SF64pred (Prim2SF f10)) <: Prim2SF (next_down f10) = SF64pred (Prim2SF f10)).
+Check (eq_refl (SF64succ (Prim2SF f11)) <: Prim2SF (next_up f11) = SF64succ (Prim2SF f11)).
+Check (eq_refl (SF64pred (Prim2SF f11)) <: Prim2SF (next_down f11) = SF64pred (Prim2SF f11)).
+Check (eq_refl (SF64succ (Prim2SF f12)) <: Prim2SF (next_up f12) = SF64succ (Prim2SF f12)).
+Check (eq_refl (SF64pred (Prim2SF f12)) <: Prim2SF (next_down f12) = SF64pred (Prim2SF f12)).
+Check (eq_refl (SF64succ (Prim2SF f13)) <: Prim2SF (next_up f13) = SF64succ (Prim2SF f13)).
+Check (eq_refl (SF64pred (Prim2SF f13)) <: Prim2SF (next_down f13) = SF64pred (Prim2SF f13)).
+Check (eq_refl (SF64succ (Prim2SF f14)) <: Prim2SF (next_up f14) = SF64succ (Prim2SF f14)).
+Check (eq_refl (SF64pred (Prim2SF f14)) <: Prim2SF (next_down f14) = SF64pred (Prim2SF f14)).
+Check (eq_refl (SF64succ (Prim2SF f15)) <: Prim2SF (next_up f15) = SF64succ (Prim2SF f15)).
+Check (eq_refl (SF64pred (Prim2SF f15)) <: Prim2SF (next_down f15) = SF64pred (Prim2SF f15)).
diff --git a/theories/Floats/FloatAxioms.v b/theories/Floats/FloatAxioms.v
index d78e3192e7..142883171e 100644
--- a/theories/Floats/FloatAxioms.v
+++ b/theories/Floats/FloatAxioms.v
@@ -8,6 +8,8 @@ Definition SF64add := SFadd prec emax.
 Definition SF64sub := SFsub prec emax.
 Definition SF64div := SFdiv prec emax.
 Definition SF64sqrt := SFsqrt prec emax.
+Definition SF64succ := SFsucc prec emax.
+Definition SF64pred := SFpred prec emax.
 
 Axiom Prim2SF_valid : forall x, valid_binary (Prim2SF x) = true.
 Axiom SF2Prim_Prim2SF : forall x, SF2Prim (Prim2SF x) = x.
@@ -45,3 +47,6 @@ Axiom normfr_mantissa_spec : forall f, to_Z (normfr_mantissa f) = Z.of_N (SFnorm
 
 Axiom frshiftexp_spec : forall f, let (m,e) := frshiftexp f in (Prim2SF m, ((to_Z e) - (to_Z shift))%Z) = SFfrexp prec emax (Prim2SF f).
 Axiom ldshiftexp_spec : forall f e, Prim2SF (ldshiftexp f e) = SFldexp prec emax (Prim2SF f) ((to_Z e) - (to_Z shift)).
+
+Axiom next_up_spec : forall x, Prim2SF (next_up x) = SF64succ (Prim2SF x).
+Axiom next_down_spec : forall x, Prim2SF (next_down x) = SF64pred (Prim2SF x).
diff --git a/theories/Floats/FloatOps.v b/theories/Floats/FloatOps.v
index 8a3ec6c181..6cc7cb0568 100644
--- a/theories/Floats/FloatOps.v
+++ b/theories/Floats/FloatOps.v
@@ -13,6 +13,8 @@ Definition ldexp f e :=
   let e' := Z.max (Z.min e (emax - emin)) (emin - emax - 1) in
   ldshiftexp f (of_Z e' + shift).
 
+Definition ulp f := ldexp one (fexp prec emax (snd (frexp f))).
+
 Definition Prim2SF f :=
   if is_nan f then S754_nan
   else if is_zero f then S754_zero (get_sign f)
diff --git a/theories/Floats/PrimFloat.v b/theories/Floats/PrimFloat.v
index bc5c49d085..b84965a11a 100644
--- a/theories/Floats/PrimFloat.v
+++ b/theories/Floats/PrimFloat.v
@@ -50,6 +50,9 @@ Definition shift := (2101)%int63. (* = 2*emax + prec *)
 Primitive frshiftexp := #float64_frshiftexp.
 Primitive ldshiftexp := #float64_ldshiftexp.
 
+Primitive next_up := #float64_next_up.
+Primitive next_down := #float64_next_down.
+
 Local Open Scope float_scope.
 
 (* Special values *)
diff --git a/theories/Floats/SpecFloat.v b/theories/Floats/SpecFloat.v
index fc26ba8cab..900739812a 100644
--- a/theories/Floats/SpecFloat.v
+++ b/theories/Floats/SpecFloat.v
@@ -349,4 +349,35 @@ Section FloatOps.
         (S754_finite sx (shift_pos (Z.to_pos d) mx) (-prec), (ex+prec-d)%Z)
     | _ => (f, (-2*emax-prec)%Z)
     end.
+
+  Definition SFone := binary_round false 1 0.
+
+  Definition SFulp x := SFldexp SFone (fexp (snd (SFfrexp x))).
+
+  Definition SFpred_pos x :=
+    match x with
+    | S754_finite _ mx _ =>
+      let d :=
+        if (mx~0 =? shift_pos (Z.to_pos prec) 1)%positive then
+          SFldexp SFone (fexp (snd (SFfrexp x) - 1))
+        else
+          SFulp x in
+      SFsub x d
+    | _ => x
+    end.
+
+  Definition SFmax_float :=
+    S754_finite false (shift_pos (Z.to_pos prec) 1 - 1) (emax - prec).
+
+  Definition SFsucc x :=
+    match x with
+    | S754_zero _ => SFldexp SFone emin
+    | S754_infinity false => x
+    | S754_infinity true => SFopp SFmax_float
+    | S754_nan => x
+    | S754_finite false _ _ => SFadd x (SFulp x)
+    | S754_finite true _ _ => SFopp (SFpred_pos (SFopp x))
+    end.
+
+  Definition SFpred f := SFopp (SFsucc (SFopp f)).
 End FloatOps.
diff --git a/vernac/g_vernac.mlg b/vernac/g_vernac.mlg
index aeae86e1f7..635b9b2c88 100644
--- a/vernac/g_vernac.mlg
+++ b/vernac/g_vernac.mlg
@@ -286,6 +286,8 @@ GRAMMAR EXTEND Gram
        | "#float64_normfr_mantissa" -> { CPrimitives.Float64normfr_mantissa }
        | "#float64_frshiftexp" -> { CPrimitives.Float64frshiftexp }
        | "#float64_ldshiftexp" -> { CPrimitives.Float64ldshiftexp }
+       | "#float64_next_up" -> { CPrimitives.Float64next_up }
+       | "#float64_next_down" -> { CPrimitives.Float64next_down }
        ] ]
     ;