1 files changed, 435 insertions, 0 deletions
diff --git a/flocq/IEEE754/SpecFloatCompat.v b/flocq/IEEE754/SpecFloatCompat.v
new file mode 100644
index 00000000..e2ace4d5
--- /dev/null
+++ b/flocq/IEEE754/SpecFloatCompat.v
@@ -0,0 +1,435 @@
+(**
+This file is part of the Flocq formalization of floating-point
+arithmetic in Coq: http://flocq.gforge.inria.fr/
+
+Copyright (C) 2018-2019 Guillaume Bertholon
+#<br />#
+Copyright (C) 2018-2019 Érik Martin-Dorel
+#<br />#
+Copyright (C) 2018-2019 Pierre Roux
+
+This library is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 3 of the License, or (at your option) any later version.
+
+This library is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+COPYING file for more details.
+*)
+
+Require Import ZArith.
+
+(** ** Inductive specification of floating-point numbers
+
+Similar to [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 (Lt | Eq) => true
+    | _ => false
+    end.
+
+  Variant float_class : Set :=
+    | PNormal | NNormal | PSubn | NSubn | PZero | NZero | PInf | NInf | NaN.
+
+  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.