From stdpp Require Import prelude ssreflect.
From Flocq.IEEE754 Require Import
  Binary BinarySingleNaN (mode_NE, mode_DN, mode_UP) Bits.
From stdpp Require Import options.

Global Arguments B754_zero {_ _}.
Global Arguments B754_infinity {_ _}.
Global Arguments B754_nan {_ _}.
Global Arguments B754_finite {_ _}.

(** The bit representation of floats is not observable in Nix, and it appears
that only negative NaNs are ever produced. So we setup the Flocq floats in
the way that it always produces [-NaN], i.e., [indef_nan]. *)
Definition float := binary64.

Variant round_mode := Floor | Ceil | NearestEven.
Global Instance round_mode_eq_dec : EqDecision round_mode.
Proof. solve_decision. Defined.

Module Float.
  Definition prec : Z := 53.
  Definition emax : Z := 1024.

  Lemma Hprec : (0 < 53)%Z.
  Proof. done. Qed.
  Lemma Hprec_emax : (53 < 1024)%Z.
  Proof. done. Qed.

  Global Instance inhabited : Inhabited float := populate (B754_zero false).

  Global Instance eq_dec : EqDecision float.
  Proof.
    refine (λ f1 f2,
      match f1, f2 with
      | B754_zero s1, B754_zero s2 => cast_if (decide (s1 = s2))
      | B754_infinity s1, B754_infinity s2 => cast_if (decide (s1 = s2))
      | B754_nan s1 pl1 _, B754_nan s2 pl2 _ =>
          cast_if_and (decide (s1 = s2)) (decide (pl1 = pl2))
      | B754_finite s1 m1 e1 _, B754_finite s2 m2 e2 _ =>
          cast_if_and3 (decide (s1 = s2)) (decide (m1 = m2)) (decide (e1 = e2))
      | _, _ => right _
      end); abstract naive_solver (f_equal; apply (proof_irrel _)).
  Defined.

  Definition of_Z (x : Z) : float :=
    binary_normalize prec emax (refl_equal _) (refl_equal _) mode_NE x 0 false.

  Definition to_Z (f : float) : option Z :=
    match f with
    | B754_zero _ => Some 0
    | B754_finite s m e _ => Some (Zaux.cond_Zopp s (Zpos m) ≪ e)
    | _ => None
    end%Z.

  (** QNaN Floating-Point Indefinite; see Table 4-3. Floating-Point Number and
  NaN Encodings. *)
  Definition indef_nan : { f | is_nan prec emax f = true } :=
    @B754_nan prec emax true (2^51) (refl_equal _) ↾ eq_refl _.

  Definition to_flocq_round_mode (m : round_mode) : BinarySingleNaN.mode :=
    match m with Floor => mode_DN | Ceil => mode_UP | NearestEven => mode_NE end.

  Definition round (m : round_mode) : float → float :=
    Bnearbyint prec emax (refl_equal _) (λ _, indef_nan) (to_flocq_round_mode m).

  (* For add: not [mode_DN]; otherwise [+0.0 + -0.0] would yield [-0.0], but
  [inf / (+0.0 + -0.0)] yields [inf] in C++, so this cannot be the case. *)
  (* C++ compiles floating point addition to the x86 ADDSD instruction. Looking
  at the Intel x86 Software Developer's Manual, it seems that the default
  rounding mode on x86 is Round to Nearest (even); see table 4-8. (In §4.8.4.) *)
  Definition add : float → float → float :=
    Bplus _ _ Hprec Hprec_emax (λ _ _, indef_nan) mode_NE.
  Definition sub : float → float → float :=
    Bminus _ _ Hprec Hprec_emax (λ _ _, indef_nan) mode_NE.
  Definition mul : float → float → float :=
    Bmult _ _ Hprec Hprec_emax (λ _ _, indef_nan) mode_NE.
  Definition div : float → float → float :=
    Bdiv _ _ Hprec Hprec_emax (λ _ _, indef_nan) mode_NE.

  Definition eqb (f1 f2 : float) : bool :=
    bool_decide (b64_compare f1 f2 = Some Eq).

  Definition ltb (f1 f2 : float) : bool :=
    bool_decide (b64_compare f1 f2 = Some Lt).
End Float.