Fix ldshiftexp

* Fix the implementations and add tests
* Change shift from int63 to Z (was always used as a Z)
* Update FloatLemmas.v accordingly

Co-authored-by: Erik Martin-Dorel <erik.martin-dorel@irit.fr>

13 files changed, 60 insertions, 46 deletions

diff --git a/kernel/byterun/coq_interp.c b/kernel/byterun/coq_interp.c
index c21aeecb16..ca1308696c 100644
--- a/kernel/byterun/coq_interp.c
+++ b/kernel/byterun/coq_interp.c
@@ -1711,7 +1711,8 @@ value coq_interprete
         print_instr("CHECKLDSHIFTEXP");
         CheckPrimArgs(Is_double(accu) && Is_uint63(sp[0]), apply2);
         Swap_accu_sp;
-        Int_of_uint63(accu);
+        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;
       }
diff --git a/kernel/byterun/coq_uint63_emul.h b/kernel/byterun/coq_uint63_emul.h
index e09803ae2d..143a6d098c 100644
--- a/kernel/byterun/coq_uint63_emul.h
+++ b/kernel/byterun/coq_uint63_emul.h
@@ -169,5 +169,5 @@ DECLARE_UNOP(of_float)
 DECLARE_UNOP(of_int)
 #define Uint63_of_int(x) CALL_UNOP(of_int, x)
 
-DECLARE_UNOP(to_int_saturate)
-#define Int_of_uint63(x) CALL_PREDICATE(accu, to_int_saturate, 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 0ab374194e..5be7587091 100644
--- a/kernel/byterun/coq_uint63_native.h
+++ b/kernel/byterun/coq_uint63_native.h
@@ -155,10 +155,9 @@ value coq_uint63_to_float_byte(value x) {
 
 #define Uint63_of_int(x) (accu = (x))
 
-#define Int_of_uint63(x) do{ \
-  uint64_t int_of_uint63__ = uint63_of_value(x); \
-  if ((int_of_uint63__ & 0xFFFFFFFF80000000L) == 0) \
-    accu = (int)int_of_uint63__; \
+#define Uint63_to_int_min(n, m) do { \
+  if (uint63_lt((n),(m))) \
+    accu = (n); \
   else \
-    accu = 0x7FFFFFFF; \
+    accu = (m); \
   }while(0)
diff --git a/kernel/float64.ml b/kernel/float64.ml
index 86b14c5cd2..07fb25734b 100644
--- a/kernel/float64.ml
+++ b/kernel/float64.ml
@@ -117,7 +117,7 @@ let frshiftexp f =
     m, Uint63.of_int (e + eshift)
 [@@ocaml.inline always]
 
-let ldshiftexp f e = ldexp f (snd (Uint63.to_int2 e) - eshift)
+let ldshiftexp f e = ldexp f (Uint63.to_int_min e (2 * eshift) - eshift)
 [@@ocaml.inline always]
 
 external next_up : float -> float = "coq_next_up_byte" "coq_next_up"
diff --git a/kernel/uint63.mli b/kernel/uint63.mli
index 7ed3d415e4..e0bf44da35 100644
--- a/kernel/uint63.mli
+++ b/kernel/uint63.mli
@@ -19,8 +19,10 @@ val to_int2 : t -> int * int (* msb, lsb *)
 val of_int64 : Int64.t -> t
 (*
 val of_uint : int -> t
-*)
-val to_int_saturate : t -> int (* maxuint31 in case of overflow *)
+ *)
+(** [int_min n m] returns the minimum of [n] and [m],
+    [m] must be in [0, 2^30-1]. *)
+val to_int_min : t -> int -> int
 
       (* conversion to float *)
 val of_float : float -> t
diff --git a/kernel/uint63_31.ml b/kernel/uint63_31.ml
index a5646e87c3..e38389ca13 100644
--- a/kernel/uint63_31.ml
+++ b/kernel/uint63_31.ml
@@ -27,9 +27,8 @@ let of_int i = Int64.of_int i
 let to_int2 i = (Int64.to_int (Int64.shift_right_logical i 31), Int64.to_int i)
 let of_int64 i = i
 
-let to_int_saturate i =
-  let r = if Int64.(equal (logand i maxuint31) i) then i else maxuint31 in
-  Int64.to_int r
+let to_int_min n m =
+  if Int64.(compare n (of_int m)) < 0 then Int64.to_int n else m
 
 let of_float f = mask63 (Int64.of_float f)
 let to_float = Int64.to_float
@@ -225,4 +224,4 @@ let () =
   Callback.register "uint63 of_float" of_float;
   Callback.register "uint63 to_float" to_float;
   Callback.register "uint63 of_int" of_int;
-  Callback.register "uint63 to_int_saturate" to_int_saturate
+  Callback.register "uint63 to_int_min" to_int_min
diff --git a/kernel/uint63_63.ml b/kernel/uint63_63.ml
index 1776f904d6..85b44528a7 100644
--- a/kernel/uint63_63.ml
+++ b/kernel/uint63_63.ml
@@ -27,8 +27,6 @@ let to_int2 i = (0,i)
 
 let of_int64 _i = assert false
 
-let to_int_saturate i = i
-
 let of_float = int_of_float
 
 external to_float : int -> (float [@unboxed])
@@ -104,6 +102,10 @@ let le (x : int) (y : int) =
   (x lxor 0x4000000000000000) <= (y lxor 0x4000000000000000)
 [@@ocaml.inline always]
 
+let to_int_min n m =
+  if lt n m then n else m
+[@@ocaml.inline always]
+
     (* division of two numbers by one *)
 (* precondition: xh < y *)
 (* outputs: q, r s.t. x = q * y + r, r < y *)
diff --git a/test-suite/primitive/float/frexp.v b/test-suite/primitive/float/frexp.v
index f13d5cebf6..2a600429b1 100644
--- a/test-suite/primitive/float/frexp.v
+++ b/test-suite/primitive/float/frexp.v
@@ -1,5 +1,3 @@
-(* TODO add tests for ldexp (particularly with overflow with 31 and 63 bits integers) *)
-
 Require Import ZArith Floats.
 
 Definition denorm := Eval compute in ldexp one (-1074)%Z.
diff --git a/test-suite/primitive/float/ldexp.v b/test-suite/primitive/float/ldexp.v
new file mode 100644
index 0000000000..a725deeeca
--- /dev/null
+++ b/test-suite/primitive/float/ldexp.v
@@ -0,0 +1,21 @@
+Require Import ZArith Int63 Floats.
+
+Check (eq_refl : ldexp one 9223372036854773807%Z = infinity).
+Check (eq_refl infinity <: ldexp one 9223372036854773807%Z = infinity).
+Check (eq_refl infinity <<: ldexp one 9223372036854773807%Z = infinity).
+
+Check (eq_refl : ldshiftexp one 9223372036854775807 = infinity).
+Check (eq_refl infinity <: ldshiftexp one 9223372036854775807 = infinity).
+Check (eq_refl infinity <<: ldshiftexp one 9223372036854775807 = infinity).
+
+Check (eq_refl : ldexp one (-2102) = 0%float).
+Check (eq_refl 0%float <: ldexp one (-2102) = 0%float).
+Check (eq_refl 0%float <<: ldexp one (-2102) = 0%float).
+
+Check (eq_refl : ldexp one (-3) = 0.125%float).
+Check (eq_refl 0.125%float <: ldexp one (-3) = 0.125%float).
+Check (eq_refl 0.125%float <<: ldexp one (-3) = 0.125%float).
+
+Check (eq_refl : ldexp one 3 = 8%float).
+Check (eq_refl 8%float <: ldexp one 3 = 8%float).
+Check (eq_refl 8%float <<: ldexp one 3 = 8%float).
diff --git a/theories/Floats/FloatAxioms.v b/theories/Floats/FloatAxioms.v
index dc637e50a6..8ca64aac42 100644
--- a/theories/Floats/FloatAxioms.v
+++ b/theories/Floats/FloatAxioms.v
@@ -51,8 +51,8 @@ Axiom sqrt_spec : forall x, Prim2SF (sqrt x) = SF64sqrt (Prim2SF x).
 Axiom of_int63_spec : forall n, Prim2SF (of_int63 n) = binary_normalize prec emax (to_Z n) 0%Z false.
 Axiom normfr_mantissa_spec : forall f, to_Z (normfr_mantissa f) = Z.of_N (SFnormfr_mantissa prec (Prim2SF f)).
 
-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 frshiftexp_spec : forall f, let (m,e) := frshiftexp f in (Prim2SF m, ((to_Z e) - 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) - 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/FloatLemmas.v b/theories/Floats/FloatLemmas.v
index 4e1f14610d..81cb7120e0 100644
--- a/theories/Floats/FloatLemmas.v
+++ b/theories/Floats/FloatLemmas.v
@@ -3,7 +3,7 @@ Require Import Psatz.
 
 (** * Support results involving frexp and ldexp *)
 
-Lemma shift_value : [|shift|]%int63 = (2*emax + prec)%Z.
+Lemma shift_value : shift = (2*emax + prec)%Z.
   reflexivity.
 Qed.
 
@@ -24,23 +24,15 @@ Theorem ldexp_spec : forall f e, Prim2SF (ldexp f e) = SFldexp prec emax (Prim2S
   destruct (Prim2SF f); auto.
   unfold SFldexp.
   unfold binary_round.
-  assert (Hmod_elim :  forall e, ([| of_Z (Z.max (Z.min e (emax - emin)) (emin - emax - 1)) + shift |]%int63 - [|shift|]%int63 = Z.max (Z.min e (emax - emin)) (emin - emax - 1))%Z).
+  assert (Hmod_elim :  forall e, ([| of_Z (Z.max (Z.min e (emax - emin)) (emin - emax - 1) + shift)|]%int63 - shift = Z.max (Z.min e (emax - emin)) (emin - emax - 1))%Z).
   {
-    intro.
-    rewrite Int63.add_spec.
-    rewrite of_Z_spec.
-    rewrite shift_value.
-    simpl.
-    unfold wB.
-    unfold size.
-    simpl.
-    unfold Z.pow_pos.
-    simpl.
+    intro e1.
+    rewrite of_Z_spec, shift_value.
+    unfold wB, size; simpl.
+    unfold Z.pow_pos; simpl.
     set (n := Z.max (Z.min _ _) _).
-    set (wB := 9223372036854775808%Z).
+    set (wB := 9223372036854775808%Z). (* Z.pow_pos 2 63 *)
     assert (-2099 <= n <= 2098)%Z by (unfold n; lia).
-    rewrite Z.add_mod_idemp_l by (unfold wB; lia).
-    destruct H as (H1, H2).
     rewrite Z.mod_small by (unfold wB; lia).
     now rewrite Z.add_simpl_r.
   }
@@ -79,17 +71,17 @@ Theorem ldexp_spec : forall f e, Prim2SF (ldexp f e) = SFldexp prec emax (Prim2S
         assert (H' : forall p p', digits2_pos (shift_pos p p') = (digits2_pos p' + p)%positive).
         {
           induction p0.
-          intro.
+          intro p'.
           simpl.
           rewrite IHp0.
           rewrite IHp0.
           lia.
-          intro.
+          intro p'.
           simpl.
           rewrite IHp0.
           rewrite IHp0.
           lia.
-          intro.
+          intro p'.
           simpl.
           lia.
         }
@@ -161,7 +153,7 @@ Theorem ldexp_spec : forall f e, Prim2SF (ldexp f e) = SFldexp prec emax (Prim2S
     {
       assert (Hshr1 : forall s, Zdigits2 (shr_m (shr_1 s)) = Z.max 0 (Zdigits2 (shr_m s) - 1)%Z).
       {
-        intro.
+        intro s0.
         destruct s0.
         unfold shr_1.
         destruct shr_m; try (simpl; lia).
@@ -170,7 +162,7 @@ Theorem ldexp_spec : forall f e, Prim2SF (ldexp f e) = SFldexp prec emax (Prim2S
       }
       induction p.
       simpl.
-      intro.
+      intro s0.
       do 2 rewrite IHp.
       rewrite Hshr1.
       lia.
@@ -183,7 +175,7 @@ Theorem ldexp_spec : forall f e, Prim2SF (ldexp f e) = SFldexp prec emax (Prim2S
 
     assert (Hd0 : forall z, Zdigits2 z = 0%Z -> z = 0%Z).
     {
-      intro.
+      intro z.
       unfold Zdigits2.
       now destruct z.
     }
diff --git a/theories/Floats/FloatOps.v b/theories/Floats/FloatOps.v
index 5ffbfc7215..f0d3bcced9 100644
--- a/theories/Floats/FloatOps.v
+++ b/theories/Floats/FloatOps.v
@@ -6,13 +6,15 @@ Definition prec := 53%Z.
 Definition emax := 1024%Z.
 Notation emin := (emin prec emax).
 
+Definition shift := 2101%Z. (** [= 2*emax + prec] *)
+
 Definition frexp f :=
   let (m, se) := frshiftexp f in
-  (m, ([| se |] - [| shift |])%Z%int63).
+  (m, ([| se |] - shift)%Z%int63).
 
 Definition ldexp f e :=
   let e' := Z.max (Z.min e (emax - emin)) (emin - emax - 1) in
-  ldshiftexp f (of_Z e' + shift).
+  ldshiftexp f (of_Z (e' + shift)).
 
 Definition ulp f := ldexp one (fexp prec emax (snd (frexp f))).
 
diff --git a/theories/Floats/PrimFloat.v b/theories/Floats/PrimFloat.v
index 880252c2b9..bc1727469d 100644
--- a/theories/Floats/PrimFloat.v
+++ b/theories/Floats/PrimFloat.v
@@ -73,8 +73,6 @@ The sign bit is always ignored. *)
 Primitive normfr_mantissa := #float64_normfr_mantissa.
 
 (** ** Exponent manipulation functions *)
-Definition shift := 2101%int63. (** [= 2*emax + prec] *)
-
 (** [frshiftexp]: convert a float to fractional part in $[0.5, 1.)$#[0.5, 1.)#
 and integer part. *)
 Primitive frshiftexp := #float64_frshiftexp.