Merge PR #9867: Add primitive floats (binary64 floating-point numbers)

Ack-by: SkySkimmer
Reviewed-by: Zimmi48
Ack-by: ejgallego
Reviewed-by: maximedenes
Ack-by: proux01
Ack-by: silene
Ack-by: vbgl

7 files changed, 978 insertions, 0 deletions

diff --git a/theories/Floats/FloatAxioms.v b/theories/Floats/FloatAxioms.v
new file mode 100644
index 0000000000..8ca64aac42
--- /dev/null
+++ b/theories/Floats/FloatAxioms.v
@@ -0,0 +1,58 @@
+Require Import ZArith Int63 SpecFloat PrimFloat FloatOps.
+
+(** * Properties of the primitive operators for the Binary64 format *)
+
+Notation valid_binary := (valid_binary prec emax).
+
+Definition SF64classify := SFclassify prec.
+Definition SF64mul := SFmul prec emax.
+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.
+Axiom Prim2SF_SF2Prim : forall x, valid_binary x = true -> Prim2SF (SF2Prim x) = x.
+
+Theorem Prim2SF_inj : forall x y, Prim2SF x = Prim2SF y -> x = y.
+  intros. rewrite <- SF2Prim_Prim2SF. symmetry. rewrite <- SF2Prim_Prim2SF. now rewrite H.
+Qed.
+
+Theorem SF2Prim_inj : forall x y, SF2Prim x = SF2Prim y -> valid_binary x = true -> valid_binary y = true -> x = y.
+  intros. rewrite <- Prim2SF_SF2Prim by assumption. symmetry. rewrite <- Prim2SF_SF2Prim by assumption. rewrite H. reflexivity.
+Qed.
+
+Axiom opp_spec : forall x, Prim2SF (-x)%float = SFopp (Prim2SF x).
+Axiom abs_spec : forall x, Prim2SF (abs x) = SFabs (Prim2SF x).
+
+Axiom eqb_spec : forall x y, (x == y)%float = SFeqb (Prim2SF x) (Prim2SF y).
+Axiom ltb_spec : forall x y, (x < y)%float = SFltb (Prim2SF x) (Prim2SF y).
+Axiom leb_spec : forall x y, (x <= y)%float = SFleb (Prim2SF x) (Prim2SF y).
+
+Definition flatten_cmp_opt c :=
+  match c with
+  | None => FNotComparable
+  | Some Eq => FEq
+  | Some Lt => FLt
+  | Some Gt => FGt
+  end.
+Axiom compare_spec : forall x y, (x ?= y)%float = flatten_cmp_opt (SFcompare (Prim2SF x) (Prim2SF y)).
+
+Axiom classify_spec : forall x, classify x = SF64classify (Prim2SF x).
+Axiom mul_spec : forall x y, Prim2SF (x * y)%float = SF64mul (Prim2SF x) (Prim2SF y).
+Axiom add_spec : forall x y, Prim2SF (x + y)%float = SF64add (Prim2SF x) (Prim2SF y).
+Axiom sub_spec : forall x y, Prim2SF (x - y)%float = SF64sub (Prim2SF x) (Prim2SF y).
+Axiom div_spec : forall x y, Prim2SF (x / y)%float = SF64div (Prim2SF x) (Prim2SF y).
+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) - 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/FloatClass.v b/theories/Floats/FloatClass.v
new file mode 100644
index 0000000000..627cb648f9
--- /dev/null
+++ b/theories/Floats/FloatClass.v
@@ -0,0 +1,2 @@
+Variant float_class : Set :=
+  | PNormal | NNormal | PSubn | NSubn | PZero | NZero | PInf | NInf | NaN.
diff --git a/theories/Floats/FloatLemmas.v b/theories/Floats/FloatLemmas.v
new file mode 100644
index 0000000000..81cb7120e0
--- /dev/null
+++ b/theories/Floats/FloatLemmas.v
@@ -0,0 +1,319 @@
+Require Import ZArith Int63 SpecFloat PrimFloat FloatOps FloatAxioms.
+Require Import Psatz.
+
+(** * Support results involving frexp and ldexp *)
+
+Lemma shift_value : shift = (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 = Z.max (Z.min e (emax - emin)) (emin - emax - 1))%Z).
+  {
+    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). (* Z.pow_pos 2 63 *)
+    assert (-2099 <= n <= 2098)%Z by (unfold n; lia).
+    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 p'.
+          simpl.
+          rewrite IHp0.
+          rewrite IHp0.
+          lia.
+          intro p'.
+          simpl.
+          rewrite IHp0.
+          rewrite IHp0.
+          lia.
+          intro p'.
+          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 s0.
+        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 s0.
+      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 z.
+      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..f0d3bcced9
--- /dev/null
+++ b/theories/Floats/FloatOps.v
@@ -0,0 +1,48 @@
+Require Import ZArith Int63 SpecFloat PrimFloat.
+
+(** * Derived operations and mapping between primitive [float]s and [spec_float]s *)
+
+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).
+
+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))).
+
+(** [Prim2SF] is an injective function that will be useful to express
+the properties of the implemented Binary64 format (see [FloatAxioms]).
+*)
+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
new file mode 100644
index 0000000000..700c69b99d
--- /dev/null
+++ b/theories/Floats/Floats.v
@@ -0,0 +1,17 @@
+(** The Floats library is split in 6 theories:
+- FloatClass: define the [float_class] inductive
+- PrimFloat: define the floating-point values and operators as kernel primitives
+- SpecFloat: specify the floating-point operators with binary integers
+- FloatOps: define conversion functions between [spec_float] and [float]
+- FloatAxioms: state properties of the primitive operators w.r.t. [spec_float]
+- FloatLemmas: prove a few results involving frexp and ldexp
+
+For a brief overview of the Floats library,
+see {{https://coq.inria.fr/distrib/current/refman/language/coq-library.html#floats-library}} *)
+
+Require Export FloatClass.
+Require Export PrimFloat.
+Require Export SpecFloat.
+Require Export FloatOps.
+Require Export FloatAxioms.
+Require Export FloatLemmas.
diff --git a/theories/Floats/PrimFloat.v b/theories/Floats/PrimFloat.v
new file mode 100644
index 0000000000..bc1727469d
--- /dev/null
+++ b/theories/Floats/PrimFloat.v
@@ -0,0 +1,118 @@
+Require Import Int63 FloatClass.
+
+(** * Definition of the interface for primitive floating-point arithmetic
+
+This interface provides processor operators for the Binary64 format of the
+IEEE 754-2008 standard. *)
+
+(** ** Type definition for the co-domain of [compare] *)
+Variant float_comparison : Set := FEq | FLt | FGt | FNotComparable.
+
+Register float_comparison as kernel.ind_f_cmp.
+
+Register float_class as kernel.ind_f_class.
+
+(** ** The main type *)
+(** [float]: primitive type for Binary64 floating-point numbers. *)
+Primitive float := #float64_type.
+
+(** ** Syntax support *)
+Declare Scope float_scope.
+Delimit Scope float_scope with float.
+Bind Scope float_scope with float.
+
+Declare ML Module "float_syntax_plugin".
+
+(** ** Floating-point operators *)
+Primitive classify := #float64_classify.
+
+Primitive abs := #float64_abs.
+
+Primitive sqrt := #float64_sqrt.
+
+Primitive opp := #float64_opp.
+Notation "- x" := (opp x) : float_scope.
+
+Primitive eqb := #float64_eq.
+Notation "x == y" := (eqb x y) (at level 70, no associativity) : float_scope.
+
+Primitive ltb := #float64_lt.
+Notation "x < y" := (ltb x y) (at level 70, no associativity) : float_scope.
+
+Primitive leb := #float64_le.
+Notation "x <= y" := (leb x y) (at level 70, no associativity) : float_scope.
+
+Primitive compare := #float64_compare.
+Notation "x ?= y" := (compare x y) (at level 70, no associativity) : float_scope.
+
+Primitive mul := #float64_mul.
+Notation "x * y" := (mul x y) : float_scope.
+
+Primitive add := #float64_add.
+Notation "x + y" := (add x y) : float_scope.
+
+Primitive sub := #float64_sub.
+Notation "x - y" := (sub x y) : float_scope.
+
+Primitive div := #float64_div.
+Notation "x / y" := (div x y) : float_scope.
+
+(** ** Conversions *)
+
+(** [of_int63]: convert a primitive integer into a float value.
+    The value is rounded if need be. *)
+Primitive of_int63 := #float64_of_int63.
+
+(** Specification of [normfr_mantissa]:
+- If the input is a float value with an absolute value inside $[0.5, 1.)$#[0.5, 1.)#;
+- Then return its mantissa as a primitive integer.
+  The mantissa will be a 53-bit integer with its most significant bit set to 1;
+- Else return zero.
+
+The sign bit is always ignored. *)
+Primitive normfr_mantissa := #float64_normfr_mantissa.
+
+(** ** Exponent manipulation functions *)
+(** [frshiftexp]: convert a float to fractional part in $[0.5, 1.)$#[0.5, 1.)#
+and integer part. *)
+Primitive frshiftexp := #float64_frshiftexp.
+
+(** [ldshiftexp]: multiply a float by an integral power of 2. *)
+Primitive ldshiftexp := #float64_ldshiftexp.
+
+(** ** Predecesor/Successor functions *)
+
+(** [next_up]: return the next float towards positive infinity. *)
+Primitive next_up := #float64_next_up.
+
+(** [next_down]: return the next float towards negative infinity. *)
+Primitive next_down := #float64_next_down.
+
+(** ** Special values (needed for pretty-printing) *)
+Definition infinity := Eval compute in div (of_int63 1) (of_int63 0).
+Definition neg_infinity := Eval compute in opp infinity.
+Definition nan := Eval compute in div (of_int63 0) (of_int63 0).
+
+Register infinity as num.float.infinity.
+Register neg_infinity as num.float.neg_infinity.
+Register nan as num.float.nan.
+
+(** ** Other special values *)
+Definition one := Eval compute in (of_int63 1).
+Definition zero := Eval compute in (of_int63 0).
+Definition neg_zero := Eval compute in (-zero)%float.
+Definition two := Eval compute in (of_int63 2).
+
+(** ** Predicates and helper functions *)
+Definition is_nan f := negb (f == f)%float.
+
+Definition is_zero f := (f == zero)%float. (* note: 0 == -0 with floats *)
+
+Definition is_infinity f := (abs f == infinity)%float.
+
+Definition is_finite (x : float) := negb (is_nan x || is_infinity x).
+
+(** [get_sign]: return [true] for [-] sign, [false] for [+] sign. *)
+Definition get_sign f :=
+  let f := if is_zero f then (one / f)%float else f in
+  (f < zero)%float.
diff --git a/theories/Floats/SpecFloat.v b/theories/Floats/SpecFloat.v
new file mode 100644
index 0000000000..fd0aa5e075
--- /dev/null
+++ b/theories/Floats/SpecFloat.v
@@ -0,0 +1,416 @@
+Require Import ZArith FloatClass.
+
+(** * Specification of floating-point arithmetic
+
+This specification is mostly borrowed from the [IEEE754.Binary] module
+of the Flocq library (see {{http://flocq.gforge.inria.fr/}}) *)
+
+(** ** Inductive specification of floating-point numbers
+
+Similar to [Flocq.IEEE754.Binary.full_float], but with no NaN payload. *)
+Variant spec_float :=
+  | S754_zero (s : bool)
+  | S754_infinity (s : bool)
+  | S754_nan
+  | S754_finite (s : bool) (m : positive) (e : Z).
+
+(** ** Parameterized definitions
+
+[prec] is the number of bits of the mantissa including the implicit one;
+[emax] is the exponent of the infinities.
+
+For instance, Binary64 is defined by [prec = 53] and [emax = 1024]. *)
+Section FloatOps.
+  Variable prec emax : Z.
+
+  Definition emin := (3-emax-prec)%Z.
+  Definition fexp e := Z.max (e - prec) emin.
+
+  Section Zdigits2.
+    Fixpoint digits2_pos (n : positive) : positive :=
+      match n with
+      | xH => xH
+      | xO p => Pos.succ (digits2_pos p)
+      | xI p => Pos.succ (digits2_pos p)
+      end.
+
+    Definition Zdigits2 n :=
+      match n with
+      | Z0 => n
+      | Zpos p => Zpos (digits2_pos p)
+      | Zneg p => Zpos (digits2_pos p)
+      end.
+  End Zdigits2.
+
+  Section ValidBinary.
+    Definition canonical_mantissa m e :=
+      Zeq_bool (fexp (Zpos (digits2_pos m) + e)) e.
+
+    Definition bounded m e :=
+      andb (canonical_mantissa m e) (Zle_bool e (emax - prec)).
+
+    Definition valid_binary x :=
+      match x with
+      | S754_finite _ m e => bounded m e
+      | _ => true
+      end.
+  End ValidBinary.
+
+  Section Iter.
+    Context {A : Type}.
+    Variable (f : A -> A).
+
+    Fixpoint iter_pos (n : positive) (x : A) {struct n} : A :=
+      match n with
+      | xI n' => iter_pos n' (iter_pos n' (f x))
+      | xO n' => iter_pos n' (iter_pos n' x)
+      | xH => f x
+      end.
+  End Iter.
+
+  Section Rounding.
+    Inductive location := loc_Exact | loc_Inexact : comparison -> location.
+
+    Record shr_record := { shr_m : Z ; shr_r : bool ; shr_s : bool }.
+
+    Definition shr_1 mrs :=
+      let '(Build_shr_record m r s) := mrs in
+      let s := orb r s in
+      match m with
+      | Z0 => Build_shr_record Z0 false s
+      | Zpos xH => Build_shr_record Z0 true s
+      | Zpos (xO p) => Build_shr_record (Zpos p) false s
+      | Zpos (xI p) => Build_shr_record (Zpos p) true s
+      | Zneg xH => Build_shr_record Z0 true s
+      | Zneg (xO p) => Build_shr_record (Zneg p) false s
+      | Zneg (xI p) => Build_shr_record (Zneg p) true s
+      end.
+
+    Definition loc_of_shr_record mrs :=
+      match mrs with
+      | Build_shr_record _ false false => loc_Exact
+      | Build_shr_record _ false true => loc_Inexact Lt
+      | Build_shr_record _ true false => loc_Inexact Eq
+      | Build_shr_record _ true true => loc_Inexact Gt
+      end.
+
+    Definition shr_record_of_loc m l :=
+      match l with
+      | loc_Exact => Build_shr_record m false false
+      | loc_Inexact Lt => Build_shr_record m false true
+      | loc_Inexact Eq => Build_shr_record m true false
+      | loc_Inexact Gt => Build_shr_record m true true
+      end.
+
+    Definition shr mrs e n :=
+      match n with
+      | Zpos p => (iter_pos shr_1 p mrs, (e + n)%Z)
+      | _ => (mrs, e)
+      end.
+
+    Definition shr_fexp m e l :=
+      shr (shr_record_of_loc m l) e (fexp (Zdigits2 m + e) - e).
+
+    Definition round_nearest_even mx lx :=
+      match lx with
+      | loc_Exact => mx
+      | loc_Inexact Lt => mx
+      | loc_Inexact Eq => if Z.even mx then mx else (mx + 1)%Z
+      | loc_Inexact Gt => (mx + 1)%Z
+      end.
+
+    Definition binary_round_aux sx mx ex lx :=
+      let '(mrs', e') := shr_fexp mx ex lx in
+      let '(mrs'', e'') := shr_fexp (round_nearest_even (shr_m mrs') (loc_of_shr_record mrs')) e' loc_Exact in
+      match shr_m mrs'' with
+      | Z0 => S754_zero sx
+      | Zpos m => if Zle_bool e'' (emax - prec) then S754_finite sx m e'' else S754_infinity sx
+      | _ => S754_nan
+      end.
+
+    Definition shl_align mx ex ex' :=
+      match (ex' - ex)%Z with
+      | Zneg d => (shift_pos d mx, ex')
+      | _ => (mx, ex)
+      end.
+
+    Definition binary_round sx mx ex :=
+      let '(mz, ez) := shl_align mx ex (fexp (Zpos (digits2_pos mx) + ex))in
+      binary_round_aux sx (Zpos mz) ez loc_Exact.
+
+    Definition binary_normalize m e szero :=
+      match m with
+      | Z0 => S754_zero szero
+      | Zpos m => binary_round false m e
+      | Zneg m => binary_round true m e
+      end.
+  End Rounding.
+
+  (** ** Define operations *)
+
+  Definition SFopp x :=
+    match x with
+    | S754_nan => S754_nan
+    | S754_infinity sx => S754_infinity (negb sx)
+    | S754_finite sx mx ex => S754_finite (negb sx) mx ex
+    | S754_zero sx => S754_zero (negb sx)
+    end.
+
+  Definition SFabs x :=
+    match x with
+    | S754_nan => S754_nan
+    | S754_infinity sx => S754_infinity false
+    | S754_finite sx mx ex => S754_finite false mx ex
+    | S754_zero sx => S754_zero false
+    end.
+
+  Definition SFcompare f1 f2 :=
+    match f1, f2 with
+    | S754_nan , _ | _, S754_nan => None
+    | S754_infinity s1, S754_infinity s2 =>
+      Some match s1, s2 with
+      | true, true => Eq
+      | false, false => Eq
+      | true, false => Lt
+      | false, true => Gt
+      end
+    | S754_infinity s, _ => Some (if s then Lt else Gt)
+    | _, S754_infinity s => Some (if s then Gt else Lt)
+    | S754_finite s _ _, S754_zero _ => Some (if s then Lt else Gt)
+    | S754_zero _, S754_finite s _ _ => Some (if s then Gt else Lt)
+    | S754_zero _, S754_zero _ => Some Eq
+    | S754_finite s1 m1 e1, S754_finite s2 m2 e2 =>
+      Some match s1, s2 with
+      | true, false => Lt
+      | false, true => Gt
+      | false, false =>
+        match Z.compare e1 e2 with
+        | Lt => Lt
+        | Gt => Gt
+        | Eq => Pcompare m1 m2 Eq
+        end
+      | true, true =>
+        match Z.compare e1 e2 with
+        | Lt => Gt
+        | Gt => Lt
+        | Eq => CompOpp (Pcompare m1 m2 Eq)
+        end
+      end
+    end.
+
+  Definition SFeqb f1 f2 :=
+    match SFcompare f1 f2 with
+    | Some Eq => true
+    | _ => false
+    end.
+
+  Definition SFltb f1 f2 :=
+    match SFcompare f1 f2 with
+    | Some Lt => true
+    | _ => false
+    end.
+
+  Definition SFleb f1 f2 :=
+    match SFcompare f1 f2 with
+    | Some Le => true
+    | _ => false
+    end.
+
+  Definition SFclassify f :=
+    match f with
+    | S754_nan => NaN
+    | S754_infinity false => PInf
+    | S754_infinity true => NInf
+    | S754_zero false => NZero
+    | S754_zero true => PZero
+    | S754_finite false m _ =>
+      if (digits2_pos m =? Z.to_pos prec)%positive then PNormal
+      else PSubn
+    | S754_finite true m _ =>
+      if (digits2_pos m =? Z.to_pos prec)%positive then NNormal
+      else NSubn
+    end.
+
+  Definition SFmul x y :=
+    match x, y with
+    | S754_nan, _ | _, S754_nan => S754_nan
+    | S754_infinity sx, S754_infinity sy => S754_infinity (xorb sx sy)
+    | S754_infinity sx, S754_finite sy _ _ => S754_infinity (xorb sx sy)
+    | S754_finite sx _ _, S754_infinity sy => S754_infinity (xorb sx sy)
+    | S754_infinity _, S754_zero _ => S754_nan
+    | S754_zero _, S754_infinity _ => S754_nan
+    | S754_finite sx _ _, S754_zero sy => S754_zero (xorb sx sy)
+    | S754_zero sx, S754_finite sy _ _ => S754_zero (xorb sx sy)
+    | S754_zero sx, S754_zero sy => S754_zero (xorb sx sy)
+    | S754_finite sx mx ex, S754_finite sy my ey =>
+      binary_round_aux (xorb sx sy) (Zpos (mx * my)) (ex + ey) loc_Exact
+    end.
+
+  Definition cond_Zopp (b : bool) m := if b then Z.opp m else m.
+
+  Definition SFadd x y :=
+    match x, y with
+    | S754_nan, _ | _, S754_nan => S754_nan
+    | S754_infinity sx, S754_infinity sy =>
+      if Bool.eqb sx sy then x else S754_nan
+    | S754_infinity _, _ => x
+    | _, S754_infinity _ => y
+    | S754_zero sx, S754_zero sy =>
+      if Bool.eqb sx sy then x else
+      S754_zero false
+    | S754_zero _, _ => y
+    | _, S754_zero _ => x
+    | S754_finite sx mx ex, S754_finite sy my ey =>
+      let ez := Z.min ex ey in
+      binary_normalize (Zplus (cond_Zopp sx (Zpos (fst (shl_align mx ex ez)))) (cond_Zopp sy (Zpos (fst (shl_align my ey ez)))))
+        ez false
+    end.
+
+  Definition SFsub x y :=
+    match x, y with
+    | S754_nan, _ | _, S754_nan => S754_nan
+    | S754_infinity sx, S754_infinity sy =>
+      if Bool.eqb sx (negb sy) then x else S754_nan
+    | S754_infinity _, _ => x
+    | _, S754_infinity sy => S754_infinity (negb sy)
+    | S754_zero sx, S754_zero sy =>
+      if Bool.eqb sx (negb sy) then x else
+      S754_zero false
+    | S754_zero _, S754_finite sy my ey => S754_finite (negb sy) my ey
+    | _, S754_zero _ => x
+    | S754_finite sx mx ex, S754_finite sy my ey =>
+      let ez := Z.min ex ey in
+      binary_normalize (Zminus (cond_Zopp sx (Zpos (fst (shl_align mx ex ez)))) (cond_Zopp sy (Zpos (fst (shl_align my ey ez)))))
+        ez false
+    end.
+
+  Definition new_location_even nb_steps k :=
+    if Zeq_bool k 0 then loc_Exact
+    else loc_Inexact (Z.compare (2 * k) nb_steps).
+
+  Definition new_location_odd nb_steps k :=
+    if Zeq_bool k 0 then loc_Exact
+    else
+      loc_Inexact
+      match Z.compare (2 * k + 1) nb_steps with
+      | Lt => Lt
+      | Eq => Lt
+      | Gt => Gt
+      end.
+
+  Definition new_location nb_steps :=
+    if Z.even nb_steps then new_location_even nb_steps else new_location_odd nb_steps.
+
+  Definition SFdiv_core_binary m1 e1 m2 e2 :=
+    let d1 := Zdigits2 m1 in
+    let d2 := Zdigits2 m2 in
+    let e' := Z.min (fexp (d1 + e1 - (d2 + e2))) (e1 - e2) in
+    let s := (e1 - e2 - e')%Z in
+    let m' :=
+      match s with
+      | Zpos _ => Z.shiftl m1 s
+      | Z0 => m1
+      | Zneg _ => Z0
+      end in
+    let '(q, r) := Z.div_eucl m' m2 in
+    (q, e', new_location m2 r).
+
+  Definition SFdiv x y :=
+    match x, y with
+    | S754_nan, _ | _, S754_nan => S754_nan
+    | S754_infinity sx, S754_infinity sy => S754_nan
+    | S754_infinity sx, S754_finite sy _ _ => S754_infinity (xorb sx sy)
+    | S754_finite sx _ _, S754_infinity sy => S754_zero (xorb sx sy)
+    | S754_infinity sx, S754_zero sy => S754_infinity (xorb sx sy)
+    | S754_zero sx, S754_infinity sy => S754_zero (xorb sx sy)
+    | S754_finite sx _ _, S754_zero sy => S754_infinity (xorb sx sy)
+    | S754_zero sx, S754_finite sy _ _ => S754_zero (xorb sx sy)
+    | S754_zero sx, S754_zero sy => S754_nan
+    | S754_finite sx mx ex, S754_finite sy my ey =>
+      let '(mz, ez, lz) := SFdiv_core_binary (Zpos mx) ex (Zpos my) ey in
+      binary_round_aux (xorb sx sy) mz ez lz
+    end.
+
+  Definition SFsqrt_core_binary m e :=
+    let d := Zdigits2 m in
+    let e' := Z.min (fexp (Z.div2 (d + e + 1))) (Z.div2 e) in
+    let s := (e - 2 * e')%Z in
+    let m' :=
+      match s with
+      | Zpos p => Z.shiftl m s
+      | Z0 => m
+      | Zneg _ => Z0
+      end in
+    let (q, r) := Z.sqrtrem m' in
+    let l :=
+      if Zeq_bool r 0 then loc_Exact
+      else loc_Inexact (if Zle_bool r q then Lt else Gt) in
+    (q, e', l).
+
+  Definition SFsqrt x :=
+    match x with
+    | S754_nan => S754_nan
+    | S754_infinity false => x
+    | S754_infinity true => S754_nan
+    | S754_finite true _ _ => S754_nan
+    | S754_zero _ => x
+    | S754_finite sx mx ex =>
+      let '(mz, ez, lz) := SFsqrt_core_binary (Zpos mx) ex in
+      binary_round_aux false mz ez lz
+    end.
+
+  Definition SFnormfr_mantissa f :=
+    match f with
+    | S754_finite _ mx ex =>
+      if Z.eqb ex (-prec) then Npos mx else 0%N
+    | _ => 0%N
+    end.
+
+  Definition SFldexp f e :=
+    match f with
+    | S754_finite sx mx ex => binary_round sx mx (ex+e)
+    | _ => f
+    end.
+
+  Definition SFfrexp f :=
+    match f with
+    | S754_finite sx mx ex =>
+      if (Z.to_pos prec <=? digits2_pos mx)%positive then
+        (S754_finite sx mx (-prec), (ex+prec)%Z)
+      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, (-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.