Add Floats to standard library

All supported floating point operations are defined on specification floats.
Then we register the primitive type and functions, and add conversion
functions to and from the specification type.
Finally we put axioms to state that primitive operations behave exactly
the same as specification operations.

CREDITS: Most of the code inside SpecFloat is adapted from the Flocq
library.

NOTE: For the moment this code will not compile if native compilation is
enabled in the configuration phase. This will be resolved later when
native_compute will be supported by primitive floats.

So please use option "-native-compiler no" in ./configure currently.

4 files changed, 501 insertions, 0 deletions

diff --git a/theories/Floats/FloatAxioms.v b/theories/Floats/FloatAxioms.v
new file mode 100644
index 0000000000..4d74edddab
--- /dev/null
+++ b/theories/Floats/FloatAxioms.v
@@ -0,0 +1,66 @@
+Require Import ZArith Int63 SpecFloat PrimFloat.
+
+(* Properties of the Binary64 IEEE 754 format *)
+Definition prec := 53%Z.
+Definition emax := 1024%Z.
+
+Notation valid_binary := (valid_binary prec emax).
+
+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 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.
+
+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 compare_spec : forall x y, (x ?= y)%float = SFcompare (Prim2SF x) (Prim2SF y).
+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 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.
diff --git a/theories/Floats/Floats.v b/theories/Floats/Floats.v
new file mode 100644
index 0000000000..818de9ffb6
--- /dev/null
+++ b/theories/Floats/Floats.v
@@ -0,0 +1,3 @@
+Require Export SpecFloat.
+Require Export PrimFloat.
+Require Export FloatAxioms.
diff --git a/theories/Floats/PrimFloat.v b/theories/Floats/PrimFloat.v
new file mode 100644
index 0000000000..4cc0a9c4d5
--- /dev/null
+++ b/theories/Floats/PrimFloat.v
@@ -0,0 +1,95 @@
+Require Import ZArith.
+Require Import Int63.
+
+Primitive float := #float64_type.
+
+Declare Scope float_scope.
+Delimit Scope float_scope with float.
+Bind Scope float_scope with float.
+
+Primitive opp := #float64_opp.
+Notation "- x" := (opp x) : float_scope.
+
+Primitive abs := #float64_abs.
+
+Register option as kernel.ind_option.
+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.
+
+Primitive sqrt := #float64_sqrt.
+
+(* Convert a primitive integer into a float value.
+ * The value is rounded if necessary. *)
+Primitive of_int63 := #float64_of_int63.
+
+(* If the input is a float value with an absolute value
+ * inside [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 *)
+Definition shift := (1022 + 52)%int63.
+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 *)
+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).
+
+Definition one := Eval compute in (of_int63 1).
+Definition zero := Eval compute in (of_int63 0).
+Definition neg_zero := Eval compute in (-zero).
+Definition two := Eval compute in (of_int63 2).
+
+Definition is_nan f :=
+  match f ?= f with
+  | None => true
+  | _ => false
+  end.
+
+Definition is_zero f :=
+  match f ?= zero with
+  | Some Eq => true (* note: 0 == -0 with floats *)
+  | _ => false
+  end.
+
+Definition is_infinity f :=
+  match f ?= infinity with
+  | Some Eq => true
+  | Some Lt => match f ?= neg_infinity with
+           | Some Eq => true
+           | _ => false
+           end
+  | _ => false
+  end.
+
+Definition get_sign f := (* + => false, - => true *)
+  let f := if is_zero f then one / f else f in
+  match f ?= zero with
+  | Some Gt => false
+  | _ => true
+  end.
+
+Definition is_finite (x : float) := negb (is_nan x || is_infinity x).
diff --git a/theories/Floats/SpecFloat.v b/theories/Floats/SpecFloat.v
new file mode 100644
index 0000000000..60a70fc230
--- /dev/null
+++ b/theories/Floats/SpecFloat.v
@@ -0,0 +1,337 @@
+Require Import ZArith.
+
+Variant spec_float :=
+  | S754_zero (s : bool)
+  | S754_infinity (s : bool)
+  | S754_nan
+  | S754_finite (s : bool) (m : positive) (e : Z).
+
+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 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, emin%Z)
+    end.
+End FloatOps.