Put axioms on ldshiftexp and frshiftexp

Axioms on ldexp and frexp are replaced by proofs inside FloatLemmas.

The shift value has been increased to 2 * emax + prec because in ldexp
we want to be able to transform the smallest denormalized to the biggest
float value in one call.

8 files changed, 374 insertions, 44 deletions

diff --git a/kernel/byterun/coq_values.h b/kernel/byterun/coq_values.h
index 14f3f152bf..fa51b2d31f 100644
--- a/kernel/byterun/coq_values.h
+++ b/kernel/byterun/coq_values.h
@@ -45,6 +45,6 @@
 #define coq_tag_Some 1
 #define coq_None Val_int(0)
 
-#define FLOAT_EXP_SHIFT (1022 + 52)
+#define FLOAT_EXP_SHIFT (2101) /* 2*emax + prec */
 
 #endif /* _COQ_VALUES_ */
diff --git a/kernel/float64.ml b/kernel/float64.ml
index 0b22e07e82..a625c0f551 100644
--- a/kernel/float64.ml
+++ b/kernel/float64.ml
@@ -51,7 +51,7 @@ let normfr_mantissa f =
   if f >= 0.5 && f < 1. then Uint63.of_float (ldexp f prec)
   else Uint63.zero
 
-let eshift = 1022 + 52 (* minimum negative exponent + binary precision *)
+let eshift = 2101 (* 2*emax + prec *)
 
 (* When calling frexp on a nan or an infinity, the returned value inside
    the exponent is undefined.
diff --git a/theories/Floats/FloatAxioms.v b/theories/Floats/FloatAxioms.v
index 4d74edddab..d057d641da 100644
--- a/theories/Floats/FloatAxioms.v
+++ b/theories/Floats/FloatAxioms.v
@@ -1,8 +1,4 @@
-Require Import ZArith Int63 SpecFloat PrimFloat.
-
-(* Properties of the Binary64 IEEE 754 format *)
-Definition prec := 53%Z.
-Definition emax := 1024%Z.
+Require Import ZArith Int63 SpecFloat PrimFloat FloatOps.
 
 Notation valid_binary := (valid_binary prec emax).
 
@@ -12,32 +8,6 @@ Definition SF64sub := SFsub prec emax.
 Definition SF64div := SFdiv prec emax.
 Definition SF64sqrt := SFsqrt prec emax.
 
-Definition Prim2SF f :=
-  if is_nan f then S754_nan
-  else if is_zero f then S754_zero (get_sign f)
-       else if is_infinity f then S754_infinity (get_sign f)
-            else
-              let (r, exp) := frexp f in
-              let e := (exp - prec)%Z in
-              let (shr, e') := shr_fexp prec emax [| normfr_mantissa r |]%int63 e loc_Exact in
-              match shr_m shr with
-              | Zpos p => S754_finite (get_sign f) p e'
-              | Zneg _ | Z0 => S754_zero false (* must never occur *)
-              end.
-
-Definition SF2Prim ef :=
-  match ef with
-  | S754_nan => nan
-  | S754_zero false => zero
-  | S754_zero true => neg_zero
-  | S754_infinity false => infinity
-  | S754_infinity true => neg_infinity
-  | S754_finite s m e =>
-    let pm := of_int63 (of_Z (Zpos m)) in
-    let f := ldexp pm e in
-    if s then (-f)%float else f
-  end.
-
 Axiom Prim2SF_valid : forall x, valid_binary (Prim2SF x) = true.
 Axiom SF2Prim_Prim2SF : forall x, SF2Prim (Prim2SF x) = x.
 Axiom Prim2SF_SF2Prim : forall x, valid_binary x = true -> Prim2SF (SF2Prim x) = x.
@@ -62,5 +32,5 @@ 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 frexp_spec : forall f, let (m,e) := frexp f in (Prim2SF m, e) = SFfrexp prec emax (Prim2SF f).
-Axiom ldexp_spec : forall f e, Prim2SF (ldexp f e) = SFldexp prec emax (Prim2SF f) e.
+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)).
diff --git a/theories/Floats/FloatLemmas.v b/theories/Floats/FloatLemmas.v
new file mode 100644
index 0000000000..7b3714ab77
--- /dev/null
+++ b/theories/Floats/FloatLemmas.v
@@ -0,0 +1,325 @@
+Require Import ZArith Int63 SpecFloat PrimFloat FloatOps FloatAxioms.
+Require Import Psatz.
+
+Lemma shift_value : [|shift|]%int63 = (2*emax + prec)%Z.
+  reflexivity.
+Qed.
+
+Theorem frexp_spec : forall f, let (m,e) := frexp f in (Prim2SF m, e) = SFfrexp prec emax (Prim2SF f).
+  intro.
+  unfold frexp.
+  case_eq (frshiftexp f).
+  intros.
+  assert (H' := frshiftexp_spec f).
+  now rewrite H in H'.
+Qed.
+
+Theorem ldexp_spec : forall f e, Prim2SF (ldexp f e) = SFldexp prec emax (Prim2SF f) e.
+  intros.
+  unfold ldexp.
+  rewrite (ldshiftexp_spec f _).
+  assert (Hv := Prim2SF_valid f).
+  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).
+  {
+    intro.
+    rewrite Int63.add_spec.
+    rewrite of_Z_spec.
+    rewrite shift_value.
+    simpl.
+    unfold wB.
+    unfold size.
+    simpl.
+    unfold Z.pow_pos.
+    simpl.
+    set (n := Z.max (Z.min _ _) _).
+    set (wB := 9223372036854775808%Z).
+    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.
+  }
+  rewrite Hmod_elim.
+  clear Hmod_elim.
+  revert Hv.
+  unfold valid_binary, bounded, canonical_mantissa.
+  unfold fexp.
+  rewrite Bool.andb_true_iff.
+  intro H'.
+  destruct H' as (H1,H2).
+  apply Zeq_bool_eq in H1.
+  apply Z.max_case_strong.
+  apply Z.min_case_strong.
+  - reflexivity.
+  - intros He _.
+    destruct (Z.max_spec (Z.pos (digits2_pos m) + e0 - prec) emin) as [ (H, Hm) | (H, Hm) ].
+    + rewrite Hm in H1.
+      rewrite <- H1.
+      rewrite !Z.max_l by (revert He; unfold emax, emin, prec; lia).
+      replace (emin + _)%Z with emax by ring.
+      unfold shl_align.
+      rewrite <- H1 in H.
+      replace (Z.pos _ + _ - _ - _)%Z with (Z.pos (digits2_pos m) - prec)%Z by ring.
+      remember (Zpos _ - _)%Z as z'.
+      destruct z' ; [ lia | lia | ].
+      unfold binary_round_aux.
+      unfold shr_fexp.
+      unfold fexp.
+      unfold Zdigits2.
+      unfold shr_record_of_loc, shr.
+      rewrite !Z.max_l by (revert H He; unfold emax, emin, prec; lia).
+      replace (_ - _)%Z with (Z.pos (digits2_pos (shift_pos p m)) - prec)%Z by ring.
+      assert (Hs : (Z.pos (digits2_pos (shift_pos p m)) <= prec)%Z).
+      {
+        assert (H' : forall p p', digits2_pos (shift_pos p p') = (digits2_pos p' + p)%positive).
+        {
+          induction p0.
+          intro.
+          simpl.
+          rewrite IHp0.
+          rewrite IHp0.
+          lia.
+          intro.
+          simpl.
+          rewrite IHp0.
+          rewrite IHp0.
+          lia.
+          intro.
+          simpl.
+          lia.
+        }
+        rewrite H'.
+        lia.
+      }
+      replace (Z.pos (digits2_pos m) + (emin + e) - prec - (emin + e))%Z with (Z.neg p) by lia.
+      unfold shr_m, loc_of_shr_record.
+      unfold round_nearest_even.
+      remember (Z.pos (digits2_pos (shift_pos p m)) - prec)%Z as ds.
+      destruct ds.
+      * rewrite Z.max_l by (revert He; unfold emax, emin, prec; lia).
+        replace (_ - _)%Z with Z0 by lia.
+        replace (_ <=? _)%Z with false by (symmetry; rewrite Z.leb_gt; lia).
+        rewrite Z.max_l by (revert He; unfold emax, emin, prec; lia).
+        replace (_ - _)%Z with Z0 by lia.
+        rewrite Z.max_l by (revert He; unfold emax, emin, prec; lia).
+        replace (_ - _)%Z with Z0 by lia.
+        replace (_ <=? _)%Z with false by (symmetry; rewrite Z.leb_gt; lia).
+        reflexivity.
+      * exfalso; lia.
+      * rewrite Z.max_l by (revert He; unfold emax, emin, prec; lia).
+        replace (_ - _)%Z with (Zneg p0) by lia.
+        replace (_ <=? _)%Z with false by (symmetry; rewrite Z.leb_gt; lia).
+        rewrite Z.max_l by (revert He; unfold emax, emin, prec; lia).
+        replace (_ - _)%Z with (Zneg p0) by lia.
+        rewrite Z.max_l by (revert He; unfold emax, emin, prec; lia).
+        replace (_ - _)%Z with (Zneg p0) by lia.
+        replace (_ <=? _)%Z with false by (symmetry; rewrite Z.leb_gt; lia).
+        reflexivity.
+    + rewrite !Z.max_l by (revert H He; unfold emax, emin, prec; lia).
+      rewrite Hm in H1.
+      clear Hm.
+      replace (Zpos _ + _ - _)%Z with (e0 + (emax - emin))%Z by (rewrite <- H1 at 1; ring).
+      replace (Zpos _ + _ - _)%Z with (e0 + e)%Z by (rewrite <- H1 at 1; ring).
+      unfold shl_align.
+      replace (_ - _)%Z with Z0 by ring.
+      replace (e0 + e - _)%Z with Z0 by ring.
+      unfold binary_round_aux.
+      unfold shr_fexp.
+      unfold fexp.
+      unfold Zdigits2.
+      rewrite !Z.max_l by (revert H He; unfold emax, emin, prec; lia).
+      unfold shr_record_of_loc.
+      unfold shr.
+      unfold Zdigits2.
+      replace (Zpos _ + _ - _ - _)%Z with Z0 by lia.
+      unfold shr_m.
+      unfold loc_of_shr_record.
+      unfold round_nearest_even.
+      rewrite Z.max_l by (revert H He; unfold emax, emin, prec; lia).
+      replace (Zpos _ + _ - _ - _)%Z with Z0 by lia.
+      replace (_ <=? _)%Z with false by (symmetry; rewrite Z.leb_gt; lia).
+      replace (Zpos _ + _ - _ - _)%Z with Z0 by lia.
+      rewrite Z.max_l by (revert H He; unfold emax, emin, prec; lia).
+      replace (Zpos _ + _ - _ - _)%Z with Z0 by lia.
+      replace (_ <=? _)%Z with false by (symmetry; rewrite Z.leb_gt; lia).
+      reflexivity.
+  - rewrite Z.min_le_iff.
+    intro H.
+    destruct H as [ He | Habs ]; [ | revert Habs; now unfold emin, emax ].
+    unfold shl_align.
+    assert (Hprec : (Z.pos (digits2_pos m) <= prec)%Z).
+    {
+      destruct (Z.max_spec (Z.pos (digits2_pos m) + e0 - prec) emin) as [ (Hpi, Hpe) | (Hpi, Hpe) ]; rewrite Hpe in H1; lia.
+    }
+
+    assert (Hshr : forall p s, Zdigits2 (shr_m (iter_pos shr_1 p s)) = Z.max Z0 (Zdigits2 (shr_m s) - Z.pos p)%Z).
+    {
+      assert (Hshr1 : forall s, Zdigits2 (shr_m (shr_1 s)) = Z.max 0 (Zdigits2 (shr_m s) - 1)%Z).
+      {
+        intro.
+        destruct s0.
+        unfold shr_1.
+        destruct shr_m; try (simpl; lia).
+        - destruct p; unfold Zdigits2, shr_m, digits2_pos; lia.
+        - destruct p; unfold Zdigits2, shr_m, digits2_pos; lia.
+      }
+      induction p.
+      simpl.
+      intro.
+      do 2 rewrite IHp.
+      rewrite Hshr1.
+      lia.
+      intros.
+      simpl.
+      do 2 rewrite IHp.
+      lia.
+      apply Hshr1.
+    }
+
+    assert (Hd0 : forall z, Zdigits2 z = 0%Z -> z = 0%Z).
+    {
+      intro.
+      unfold Zdigits2.
+      now destruct z.
+    }
+
+    assert (Hshr_p0 : forall p0, (prec < Z.pos p0)%Z -> shr_m (iter_pos shr_1 p0 {| shr_m := Z.pos m; shr_r := false; shr_s := false |}) = Z0).
+    {
+      intros p0 Hp0.
+      apply Hd0.
+      rewrite Hshr.
+      rewrite Z.max_l; [ reflexivity | ].
+      unfold shr_m.
+      unfold Zdigits2.
+      lia.
+    }
+
+    assert (Hshr_p0_r : forall p0, (prec < Z.pos p0)%Z -> shr_r (iter_pos shr_1 p0 {| shr_m := Z.pos m; shr_r := false; shr_s := false |}) = false).
+    {
+      intros p0 Hp0.
+
+      assert (Hshr_p0m1 : shr_m (iter_pos shr_1 (p0-1) {| shr_m := Z.pos m; shr_r := false; shr_s := false |}) = Z0).
+      {
+        apply Hd0.
+        rewrite Hshr.
+        rewrite Z.max_l; [ reflexivity | ].
+        unfold shr_m.
+        unfold Zdigits2.
+        lia.
+      }
+
+      assert (Hiter_pos : forall A (f : A -> A) p e, iter_pos f (p + 1) e = f (iter_pos f p e)).
+      {
+        assert (Hiter_pos' : forall A (f : A -> A) p e, iter_pos f p (f e) = f (iter_pos f p e)).
+        {
+          intros A f'.
+          induction p.
+          intro e'.
+          simpl.
+          now do 2 rewrite IHp.
+          intro e'.
+          simpl.
+          now do 2 rewrite IHp.
+          intro e'.
+          now simpl.
+        }
+        intros A f'.
+        induction p.
+        intros.
+        simpl.
+        rewrite <- Pos.add_1_r.
+        do 2 rewrite IHp.
+        now do 3 rewrite Hiter_pos'.
+        intros.
+        simpl.
+        now do 2 rewrite Hiter_pos'.
+        intros.
+        now simpl.
+      }
+      replace p0 with (p0 - 1 + 1)%positive.
+      rewrite Hiter_pos.
+      unfold shr_1 at 1.
+      remember (iter_pos _ _ _) as shr_p0m1.
+      destruct shr_p0m1.
+      unfold SpecFloat.shr_m in Hshr_p0m1.
+      now rewrite Hshr_p0m1.
+      rewrite Pos.add_1_r.
+      rewrite Pos.sub_1_r.
+      apply Pos.succ_pred.
+      lia.
+    }
+
+    rewrite Z.leb_le in H2.
+
+    destruct (Z.max_spec (Z.pos (digits2_pos m) + (e0 + (emin - emax - 1)) - prec) emin) as [ (H, Hm) | (H, Hm) ].
+    + rewrite Hm.
+      replace (_ - _)%Z with (emax - e0 + 1)%Z by ring.
+      remember (emax - e0 + 1)%Z as z'.
+      destruct z'; [ exfalso; lia | | exfalso; lia ].
+      unfold binary_round_aux.
+      unfold shr_fexp, fexp.
+      unfold shr, shr_record_of_loc.
+      unfold Zdigits2.
+      rewrite Hm.
+      replace (_ - _)%Z with (Z.pos p) by (rewrite Heqz'; ring).
+      set (rne := round_nearest_even _ _).
+      assert (rne = 0%Z).
+      {
+        unfold rne.
+        unfold round_nearest_even.
+
+        assert (Hp0 : (prec < Z.pos p)%Z) by lia.
+
+        unfold loc_of_shr_record.
+        specialize (Hshr_p0_r _ Hp0).
+        specialize (Hshr_p0 _ Hp0).
+        revert Hshr_p0_r Hshr_p0.
+        set (shr_p0 := iter_pos shr_1 _ _).
+        destruct shr_p0.
+        unfold SpecFloat.shr_r, SpecFloat.shr_m.
+        intros Hshr_r Hshr_m.
+        rewrite Hshr_r, Hshr_m.
+        now destruct shr_s.
+      }
+
+      rewrite H0.
+      rewrite Z.max_r by (rewrite Heqz'; unfold prec; lia).
+      replace (_ - _)%Z with 0%Z by lia.
+      unfold shr_m.
+
+      rewrite Z.max_r by lia.
+      remember (emin - (e0 + e))%Z as eminmze.
+      destruct eminmze; [ exfalso; lia | | exfalso; lia ].
+
+      rewrite Z.max_r by lia.
+      rewrite <- Heqeminmze.
+
+      set (rne' := round_nearest_even _ _).
+      assert (Hrne'0 : rne' = 0%Z).
+      {
+        unfold rne'.
+        unfold round_nearest_even.
+
+        assert (Hp1 : (prec < Z.pos p0)%Z) by lia.
+
+        unfold loc_of_shr_record.
+        specialize (Hshr_p0_r _ Hp1).
+        specialize (Hshr_p0 _ Hp1).
+        revert Hshr_p0_r Hshr_p0.
+        set (shr_p1 := iter_pos shr_1 _ _).
+        destruct shr_p1.
+        unfold SpecFloat.shr_r, SpecFloat.shr_m.
+        intros Hshr_r Hshr_m.
+        rewrite Hshr_r, Hshr_m.
+        now destruct shr_s.
+      }
+
+      rewrite Hrne'0.
+      rewrite Z.max_r by (rewrite Heqeminmze; unfold prec; lia).
+      replace (_ - _)%Z with 0%Z by lia.
+      reflexivity.
+    + exfalso; lia.
+Qed.
diff --git a/theories/Floats/FloatOps.v b/theories/Floats/FloatOps.v
new file mode 100644
index 0000000000..8a3ec6c181
--- /dev/null
+++ b/theories/Floats/FloatOps.v
@@ -0,0 +1,40 @@
+Require Import ZArith Int63 SpecFloat PrimFloat.
+
+(* Properties of the Binary64 IEEE 754 format *)
+Definition prec := 53%Z.
+Definition emax := 1024%Z.
+Notation emin := (emin prec emax).
+
+Definition frexp f :=
+  let (m, se) := frshiftexp f in
+  (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).
+
+Definition Prim2SF f :=
+  if is_nan f then S754_nan
+  else if is_zero f then S754_zero (get_sign f)
+       else if is_infinity f then S754_infinity (get_sign f)
+            else
+              let (r, exp) := frexp f in
+              let e := (exp - prec)%Z in
+              let (shr, e') := shr_fexp prec emax [| normfr_mantissa r |]%int63 e loc_Exact in
+              match shr_m shr with
+              | Zpos p => S754_finite (get_sign f) p e'
+              | Zneg _ | Z0 => S754_zero false (* must never occur *)
+              end.
+
+Definition SF2Prim ef :=
+  match ef with
+  | S754_nan => nan
+  | S754_zero false => zero
+  | S754_zero true => neg_zero
+  | S754_infinity false => infinity
+  | S754_infinity true => neg_infinity
+  | S754_finite s m e =>
+    let pm := of_int63 (of_Z (Zpos m)) in
+    let f := ldexp pm e in
+    if s then (-f)%float else f
+  end.
diff --git a/theories/Floats/Floats.v b/theories/Floats/Floats.v
index 818de9ffb6..a1d49e87ee 100644
--- a/theories/Floats/Floats.v
+++ b/theories/Floats/Floats.v
@@ -1,3 +1,5 @@
 Require Export SpecFloat.
 Require Export PrimFloat.
+Require Export FloatOps.
 Require Export FloatAxioms.
+Require Export FloatLemmas.
diff --git a/theories/Floats/PrimFloat.v b/theories/Floats/PrimFloat.v
index 4cc0a9c4d5..21a1be8708 100644
--- a/theories/Floats/PrimFloat.v
+++ b/theories/Floats/PrimFloat.v
@@ -1,4 +1,3 @@
-Require Import ZArith.
 Require Import Int63.
 
 Primitive float := #float64_type.
@@ -41,16 +40,10 @@ Primitive of_int63 := #float64_of_int63.
 Primitive normfr_mantissa := #float64_normfr_mantissa.
 
 (* Exponent manipulation functions *)
-Definition shift := (1022 + 52)%int63.
+Definition shift := (2101)%int63. (* = 2*emax + prec *)
 Primitive frshiftexp := #float64_frshiftexp.
 Primitive ldshiftexp := #float64_ldshiftexp.
 
-Definition frexp f :=
-  let (m, se) := frshiftexp f in
-  (m, ([| se |] - [| shift |])%Z%int63).
-
-Definition ldexp f e := ldshiftexp f (of_Z e + shift).
-
 Local Open Scope float_scope.
 
 (* Special values *)
diff --git a/theories/Floats/SpecFloat.v b/theories/Floats/SpecFloat.v
index 60a70fc230..337c00b20f 100644
--- a/theories/Floats/SpecFloat.v
+++ b/theories/Floats/SpecFloat.v
@@ -332,6 +332,6 @@ Section FloatOps.
       else
         let d := (prec - Z.pos (digits2_pos mx))%Z in
         (S754_finite sx (shift_pos (Z.to_pos d) mx) (-prec), (ex+prec-d)%Z)
-    | _ => (f, emin%Z)
+    | _ => (f, (-2*emax-prec)%Z)
     end.
 End FloatOps.