Module FloatLib

Require Import Coqlib Utf8 Util Axioms ToString.
Require Import ZArith.
Require Import Integers.
Require Import Reals.
Require Export Fappli_IEEE.
Require Export Fappli_IEEE_bits.
Require Export Fappli_IEEE_extra.
Require Export Fappli_double_round.
Require Export Fcore_FLT.
Require Export Fcore_float_prop.
Require Export Fcore_defs.
Require Export Fcore_Raux.
Require Export Fcore_Zaux.
Require Export Fcore_generic_fmt.
Require Export Fcore_rnd_ne.
Require Export Fcore_digits.
Require Export Floats.

Open Scope R.

Arguments is_finite {prec emax} f.
Arguments is_finite_strict {prec emax} f.
Arguments is_nan {prec emax} f.
Arguments B2R {prec emax} f.
Arguments Bsign {prec emax} x.
Arguments ZofB {prec emax} f.
Arguments B2FF {prec emax} x.
Arguments B754_zero {prec emax} _.
Arguments B754_infinity {prec emax} _.
Arguments B754_nan {prec emax} _ _.
Arguments B754_finite {prec emax} _ m e _.

Identity Coercion f2b64 : float >-> binary64.
Identity Coercion f2b32 : float32 >-> binary32.
Identity Coercion b642bf : binary64 >-> binary_float.
Identity Coercion b322bf : binary32 >-> binary_float.
Coercion B2R : binary_float >-> R.
Coercion Z2R : Z >-> R.
Coercion B2FF : binary_float >-> full_float.

Ltac Rlt_bool_case H :=
  match type of H with
    | context [Rlt_bool ?x ?y] =>
      destruct (Rlt_bool_spec x y)
  end.

Ltac Rle_bool_case H :=
  match type of H with
    | context [Rle_bool ?x ?y] =>
      destruct (Rle_bool_spec x y)
  end.

Ltac check_nodup H :=
  let typ := type of H in
  match goal with
    | [H':typ |- _] =>
      match H' with
        | H => fail 1
        | _ => fail 2
      end
    | _ => idtac
  end.

Ltac use_float_spec :=
  let H := fresh in
  match goal with
    | [_ : context [binary_normalize ?prec ?emax ?Hprec ?Hemax ?mode ?m ?e ?s] |- _] =>
      assert (H:=binary_normalize_correct prec emax Hprec Hemax mode m e s);
      check_nodup H
    | [|- context [binary_normalize ?prec ?emax ?Hprec ?Hemax ?mode ?m ?e ?s]] =>
      assert (H:=binary_normalize_correct prec emax Hprec Hemax mode m e s);
      check_nodup H
    | [_ : context [Bconv ?prec1 ?emax1 ?prec2 ?emax2 ?Hprec2 ?Hemax2 ?pl ?mode ?x] |- _] =>
      assert (H:=Bconv_correct prec1 emax1 prec2 emax2 Hprec2 Hemax2 pl mode x);
      check_nodup H
    | [|- context [Bconv ?prec1 ?emax1 ?prec2 ?emax2 ?Hprec2 ?Hemax2 ?pl ?mode ?x]] =>
      assert (H:=Bconv_correct prec1 emax1 prec2 emax2 Hprec2 Hemax2 pl mode x);
      check_nodup H
    | [_ : context [BofZ ?prec ?emax ?Hprec ?Hemax ?x] |- _] =>
      assert (H:=BofZ_correct prec emax Hprec Hemax x);
      check_nodup H
    | [|- context [BofZ ?prec ?emax ?Hprec ?Hemax ?x]] =>
      assert (H:=BofZ_correct prec emax Hprec Hemax x);
      check_nodup H
    | [_ : context [Bplus ?prec ?emax ?Hprec ?Hemax ?pl ?mode ?x ?y] |- _] =>
      assert (H:=Bplus_correct prec emax Hprec Hemax pl mode x y);
      check_nodup H
    | [|-context [Bplus ?prec ?emax ?Hprec ?Hemax ?pl ?mode ?x ?y]] =>
      assert (H:=Bplus_correct prec emax Hprec Hemax pl mode x y);
      check_nodup H
    | [_ : context [Bminus ?prec ?emax ?Hprec ?Hemax ?pl ?mode ?x ?y] |- _] =>
      assert (H:=Bminus_correct prec emax Hprec Hemax pl mode x y);
      check_nodup H
    | [|-context [Bminus ?prec ?emax ?Hprec ?Hemax ?pl ?mode ?x ?y]] =>
      assert (H:=Bminus_correct prec emax Hprec Hemax pl mode x y);
      check_nodup H
    | [_ : context [Bmult ?prec ?emax ?Hprec ?Hemax ?pl ?mode ?x ?y] |- _] =>
      assert (H:=Bmult_correct prec emax Hprec Hemax pl mode x y);
      check_nodup H
    | [|-context [Bmult ?prec ?emax ?Hprec ?Hemax ?pl ?mode ?x ?y]] =>
      assert (H:=Bmult_correct prec emax Hprec Hemax pl mode x y);
      check_nodup H
    | [_ : context [Bdiv ?prec ?emax ?Hprec ?Hemax ?pl ?mode ?x ?y] |- _] =>
      assert (H:=Bdiv_correct prec emax Hprec Hemax pl mode x y);
      check_nodup H
    | [|-context [Bdiv ?prec ?emax ?Hprec ?Hemax ?pl ?mode ?x ?y]] =>
      assert (H:=Bdiv_correct prec emax Hprec Hemax pl mode x y);
      check_nodup H
  end.

Ltac use_float_specs := repeat use_float_spec.

Transparent Float32.add Float.add Float32.sub Float.sub Float32.mul Float.mul
            Float32.div Float.div Float.to_single Float.of_single Float.to_int
            Float32.to_int Float.to_intu Float32.to_intu Float.to_long
            Float.to_longu Float32.to_long Float32.to_longu
            Float.neg Float.cmp.
Transparent Float.transform_quiet_pl Float.reduce_pl Float.to_single_pl.
Transparent Archi.choose_binop_pl_32 Archi.choose_binop_pl_64
            Archi.default_pl_32 Archi.default_pl_64
            Archi.float_of_single_preserves_sNaN.

Lemma prec53_gt_0: Fcore_FLX.Prec_gt_0 53.

Hint Resolve fexp_correct valid_rnd_round_mode prec53_gt_0.

Lemma Bsign_false_ge_0:
  ∀ prec emax (f:binary_float prec emax), Bsign f = false -> 0 <= f.

Lemma Bsign_true_le_0:
    ∀ prec emax (f:binary_float prec emax), Bsign f = true -> f <= 0.

Definition Fleb {prec emax} f1 f2 :=
  match Bcompare prec emax f1 f2 with
  | Some (Eq|Lt) => true
  | Some Gt | None => false
  end.

Lemma Fleb_Rle:
  ∀ prec emax (f1 f2:binary_float prec emax),
    is_finite f1 = true -> is_finite f2 = true ->
    Fleb f1 f2 = Rle_bool f1 f2.

Lemma Bcompare_nan:
  ∀ prec emax f1 f2,
    is_nan f1 = false -> is_nan f2 = false ->
    Bcompare prec emax f1 f2 <> None.

Lemma Fleb_total:
  ∀ prec emax (f1 f2:binary_float prec emax),
    is_nan f1 = false -> is_nan f2 = false ->
    Fleb f1 f2 = true \/ Fleb f2 f1 = true.

Lemma Bcompare_refl:
  ∀ prec emax f,
    is_nan f = false ->
    Bcompare prec emax f f = Some Eq.

Lemma Fleb_refl:
  ∀ prec emax (f:binary_float prec emax), is_nan f = false -> Fleb f f = true.

Hint Resolve Fleb_refl.

Lemma Bcompare_trans:
  ∀ prec emax f1 f2 f3 c,
    is_nan f2 = false ->
    Bcompare prec emax f1 f2 = c ->
    Bcompare prec emax f2 f3 = c ->
    Bcompare prec emax f1 f3 = c.

Lemma Bcompare_Eq_trans:
  ∀ prec emax f1 f2 f3 c,
    Bcompare prec emax f1 f2 = Some Eq ->
    Bcompare prec emax f2 f3 = c ->
    Bcompare prec emax f1 f3 = c.

Lemma Bcompare_trans_Eq:
  ∀ prec emax f1 f2 f3 c,
    Bcompare prec emax f1 f2 = c ->
    Bcompare prec emax f2 f3 = Some Eq ->
    Bcompare prec emax f1 f3 = c.

Lemma Fleb_nan_l:
  ∀ prec emax (f1 f2:binary_float prec emax),
    Fleb f1 f2 = true -> is_nan f1 = false.

Hint Resolve Fleb_nan_l.

Lemma Fleb_nan_r:
  ∀ prec emax (f1 f2:binary_float prec emax),
    Fleb f1 f2 = true -> is_nan f2 = false.

Hint Resolve Fleb_nan_r.

Lemma Fleb_trans:
  ∀ prec emax (f1 f2 f3:binary_float prec emax),
    Fleb f1 f2 = true ->
    Fleb f2 f3 = true ->
    Fleb f1 f3 = true.

Lemma fneg_decr:
  forall a b, Fleb (Float.neg b) (Float.neg a) = Fleb a b.

Lemma binary_normalize_is_nan:
  forall prec emax H1 H2 m e s mode,
    is_nan (binary_normalize prec emax H1 H2 mode m e s) = false.

Lemma floatofsingle_exact:
  forall f, (Float.of_single f:R) = f
            /\ (is_finite f = true -> is_finite (Float.of_single f) = true).

Lemma floatofsingle_nan:
  ∀ f, is_nan f = false ->
       is_nan (Float.of_single f) = false.

Lemma zoffloat_incr:
  forall prec emax (f1 f2:binary_float prec emax) z1 z2,
    Fleb f1 f2 = true ->
    ZofB f1 = Some z1 ->
    ZofB f2 = Some z2 ->
    (z1 <= z2)%Z.

Definition Zoffloat_DN (f:float): option Z :=
  match f with
  | B754_finite s m (Zpos e) _ => Some (cond_Zopp s (Zpos m) * Zpower_pos radix2 e)%Z
  | B754_finite s m 0 _ => Some (cond_Zopp s (Zpos m))
  | B754_finite s m (Zneg e) _ => Some (cond_Zopp s (Zpos m) / Zpower_pos radix2 e)%Z
  | B754_zero _ => Some 0%Z
  | _ => None
  end.

Lemma Zoffloat_DN_correct:
  forall f : float,
    match Zoffloat_DN f with
    | Some n =>
      is_finite f = true /\ (n:R) = round radix2 (Fcore_FIX.FIX_exp 0) (round_mode mode_DN) f
    | None => is_finite f = false
    end.

Lemma fadd_commut:
  forall a b pl,
    Bplus 53 1024 (@eq_refl _ Lt) (@eq_refl _ Lt) pl mode_NE a b =
    Bplus 53 1024 (@eq_refl _ Lt) (@eq_refl _ Lt) (fun a b => pl b a) mode_NE b a.

Lemma Bconv_Bopp:
  ∀ prec emax
    prec' emax' Hprec' Hemax'
    conv_pl opp_pl opp_pl'
    x,
    (∀ s pl, (let '(s, pl) := opp_pl s pl in conv_pl s pl) =
             (let '(s, pl) := conv_pl s pl in opp_pl' s pl)) ->
    Bopp prec' emax' opp_pl' (Bconv prec emax prec' emax' Hprec' Hemax' conv_pl mode_NE x) =
    Bconv prec emax prec' emax' Hprec' Hemax' conv_pl mode_NE (Bopp prec emax opp_pl x).

Lemma Bconv_Babs:
  ∀ prec emax
    prec' emax' Hprec' Hemax'
    conv_pl abs_pl abs_pl'
    x,
    (∀ s pl, (let '(s, pl) := abs_pl s pl in conv_pl s pl) =
             (let '(s, pl) := conv_pl s pl in abs_pl' s pl)) ->
    Babs prec' emax' abs_pl' (Bconv prec emax prec' emax' Hprec' Hemax' conv_pl mode_NE x) =
    Bconv prec emax prec' emax' Hprec' Hemax' conv_pl mode_NE (Babs prec emax abs_pl x).

Lemma Bconv_Bplus:
  ∀ prec emax Hprec Hemax
    prec' emax' Hprec' Hemax'
    conv_pl conv_pl' plus_pl plus_pl'
    x y,
    (2*prec + 1 <= prec')%Z -> (emax < emax')%Z ->

    is_nan x = false -> is_nan y = false ->

    (∀ x y, is_nan x = false -> is_nan y = false ->
     ∀ x' y', is_nan x' = false -> is_nan y' = false ->
     (let '(b, pl) := plus_pl' x' y' in conv_pl b pl) = plus_pl x y) ->

    Bplus prec emax Hprec Hemax plus_pl mode_NE x y =
    Bconv prec' emax' prec emax Hprec Hemax conv_pl mode_NE
      (Bplus prec' emax' Hprec' Hemax' plus_pl' mode_NE
        (Bconv prec emax prec' emax' Hprec' Hemax' conv_pl' mode_NE x)
        (Bconv prec emax prec' emax' Hprec' Hemax' conv_pl' mode_NE y)).

Lemma Bconv_Bmult:
  ∀ prec emax Hprec Hemax
    prec' emax' Hprec' Hemax'
    conv_pl conv_pl' mult_pl mult_pl'
    x y,
    (2*prec <= prec')%Z -> (2*emax <= emax')%Z ->

    is_nan x = false -> is_nan y = false ->

    (∀ x y, is_nan x = false -> is_nan y = false ->
     ∀ x' y', is_nan x' = false -> is_nan y' = false ->
     (let '(b, pl) := mult_pl' x' y' in conv_pl b pl) = mult_pl x y) ->

    Bmult prec emax Hprec Hemax mult_pl mode_NE x y =
    Bconv prec' emax' prec emax Hprec Hemax conv_pl mode_NE
      (Bmult prec' emax' Hprec' Hemax' mult_pl' mode_NE
        (Bconv prec emax prec' emax' Hprec' Hemax' conv_pl' mode_NE x)
        (Bconv prec emax prec' emax' Hprec' Hemax' conv_pl' mode_NE y)).

Lemma Bconv_Bdiv:
  ∀ prec emax Hprec Hemax
    prec' emax' Hprec' Hemax'
    conv_pl conv_pl' div_pl div_pl'
    x y,
    (2 * prec <= prec')%Z -> (2 + emax + 2 * prec <= emax' + prec')%Z ->
    (-3 + 2 * emax + prec <= emax')%Z ->

    is_nan x = false -> is_nan y = false ->

    (∀ x y, is_nan x = false -> is_nan y = false ->
     ∀ x' y', is_nan x' = false -> is_nan y' = false ->
     (let '(b, pl) := div_pl' x' y' in conv_pl b pl) = div_pl x y) ->

    Bdiv prec emax Hprec Hemax div_pl mode_NE x y =
    Bconv prec' emax' prec emax Hprec Hemax conv_pl mode_NE
      (Bdiv prec' emax' Hprec' Hemax' div_pl' mode_NE
        (Bconv prec emax prec' emax' Hprec' Hemax' conv_pl' mode_NE x)
        (Bconv prec emax prec' emax' Hprec' Hemax' conv_pl' mode_NE y)).

Lemma transform_quiet_32_64:
  ∀ n,
    Float32.transform_quiet_pl n =
    Float.reduce_pl (Float.transform_quiet_pl (Float.expand_pl n)).

Lemma plus_double_round:
  ∀ (x y:float32),
    Float32.add x y = Float.to_single (Float.add (Float.of_single x) (Float.of_single y)).

Lemma sub_double_round:
  ∀ (x y:float32),
    Float32.sub x y = Float.to_single (Float.sub (Float.of_single x) (Float.of_single y)).

Lemma mult_double_round:
  ∀ (x y:float32),
    Float32.mul x y = Float.to_single (Float.mul (Float.of_single x) (Float.of_single y)).

Lemma div_double_round:
  ∀ (x y:float32),
    Float32.div x y = Float.to_single (Float.div (Float.of_single x) (Float.of_single y)).

Instance FloatToString : ToString float :=
  fun f =>
    match f with
    | Fappli_IEEE.B754_zero sg => if sg then "-0.0" else "+0.0"
    | Fappli_IEEE.B754_infinity sg => if sg then "-∞" else "+∞"
    | Fappli_IEEE.B754_nan sg pl =>
      (if sg then "NAN(-," else "NAN(+,") ++ (to_string (proj1_sig pl)) ++ ")"
    | Fappli_IEEE.B754_finite sg m e _ =>
      if Z.eqb e (-1074) then
        (if sg then "-" else "+") ++ (to_string m) ++ ".p" ++ (to_string (e))
      else
        (if sg then "-." else "+.") ++ (to_string m) ++ "p" ++ (to_string (e+56)%Z)
    end.

Instance Float32ToString : ToString float32 :=
  fun f =>
    match f with
    | Fappli_IEEE.B754_zero sg => if sg then "-0.0" else "+0.0"
    | Fappli_IEEE.B754_infinity sg => if sg then "-∞" else "+∞"
    | Fappli_IEEE.B754_nan sg pl =>
      (if sg then "NAN(-," else "NAN(+,") ++ (to_string (proj1_sig pl)) ++ ")"
    | Fappli_IEEE.B754_finite sg m e _ =>
      if Z.eqb e (-149) then
        (if sg then "-" else "+") ++ (to_string m) ++ ".p" ++ (to_string e)
      else
        (if sg then "-." else "+.") ++ (to_string m) ++ "p" ++ (to_string (e+24)%Z)
    end.

Definition one : float :=
  BofZ 53 1024 eq_refl eq_refl 1.