Module FloatIntervalsBackward

Require Import Coqlib Utf8 ToString Util FastArith.
Require Import AdomLib.
Require Import Integers.
Require Import ZIntervals.
Require Import FloatLib.
Require Import FloatIntervals.

Definition RFF_lt (x:R) (y:full_float) :=
  if is_finite_FF y then (x < FF2R radix2 y)
  else y = F754_infinity false.

Definition RF_lt {prec emax} (x:R) (y:binary_float prec emax) :=
  RFF_lt x y.

Local Instance gammaR : gamma_op fitv R := fun i x =>
  let 'FITV mx Mx := i in
  RF_lt x Mx /\ RF_lt (-x) (Float.neg mx).

Definition next_ulp_ff prec emax (f:full_float) : full_float :=
  match f return _ with
  | F754_nan _ _ => f
  | F754_infinity true => binary_overflow prec emax mode_UP true
  | F754_finite true m e =>
    let m' := Pos.pred m in
    if Z.eqb e (3-emax-prec) then
      if Pos.eqb m 1 then F754_zero false
      else F754_finite true m' e
    else
      let minmant := (match 2 ^ (prec-1) with Zpos m => m | _ => xH end)%Z in
      if Pos.eqb m minmant then F754_finite true (Pos.pred_double minmant) (Z.pred e)
      else F754_finite true m' e
  | F754_zero _ => F754_finite false 1 (3-emax-prec)
  | F754_finite false m e =>
    let m' := Pos.succ m in
    let maxmant := (match 2 ^ prec with Zpos m => m | _ => xH end)%Z in
    if Pos.eqb m' maxmant then
      if Z.eqb e (emax-prec) then F754_infinity false
      else F754_finite false (Pos.div2 maxmant)%Z (Z.succ e)
    else
      F754_finite false m' e
  | F754_infinity false => f
  end.

Lemma next_ulp_ff_is_nan:
  forall prec emax f,
    is_nan_FF f = false ->
    is_nan_FF (next_ulp_ff prec emax f) = false.

Lemma pow2_pos :
  forall prec, (0 < prec)%Z ->
    exists n, (2 ^ (prec-1))%Z = Zpos n /\
              Z.pos (digits2_pos n) = prec /\
              Z.pos (digits2_pos (Pos.pred_double n)) = prec.

Lemma next_ulp_ff_valid_binary:
  forall prec emax
         (prec_gt_0_ : (0 < prec)%Z) (Hmax : (prec < emax)%Z)
         (f:binary_float prec emax),
    valid_binary prec emax (next_ulp_ff prec emax f) = true.

Definition next_ulp prec emax prec_gt_0_ Hmax
           (f:binary_float prec emax) : binary_float prec emax :=
  FF2B prec emax _ (next_ulp_ff_valid_binary _ _ prec_gt_0_ Hmax f).

Lemma next_ulp_is_nan:
  forall prec emax prec_gt_0_ Hmax f,
    is_nan f = false ->
    is_nan (next_ulp prec emax prec_gt_0_ Hmax f) = false.

Lemma next_ulp_round :
  forall prec emax prec_gt_0_ Hmax, (1 < prec)%Z ->
  forall (f:binary_float prec emax), is_finite f = true ->
  forall x mode,
    (round radix2 (FLT_exp (3 - emax - prec) prec) (round_mode mode) x <= f) ->
    RF_lt x (next_ulp _ _ prec_gt_0_ Hmax f).

Lemma next_ulp_round_overflow :
    forall prec emax prec_gt_0_ Hmax x mode (f:binary_float _ _), (1 < prec)%Z ->
    (if Rlt_bool (Rabs (round radix2 (FLT_exp (3 - emax - prec) prec) (round_mode mode) x))
                 (bpow radix2 emax)
     then (f:R) = round radix2 (FLT_exp (3 - emax - prec) prec) (round_mode mode) x /\
          is_finite f = true
     else (f:full_float) = binary_overflow prec emax mode (Rlt_bool x 0)) ->
  forall fM, Fleb f fM = true ->
  RF_lt x (next_ulp prec emax prec_gt_0_ Hmax fM).

Lemma next_ulp_lub:
  forall prec emax prec_gt_0_ Hmax (f1 f2:binary_float prec emax), (1 < prec)%Z ->
    is_nan f1 = false ->
    is_nan f2 = false ->
    Fleb f1 f2 = true \/
    Fleb (next_ulp _ _ prec_gt_0_ Hmax f2) f1 = true.

Definition backward_round (x:fitv) : fitv :=
  let 'FITV mx Mx := x in
  FITV (Float.neg (next_ulp 53 1024 eq_refl eq_refl (Float.neg mx)))
       (next_ulp 53 1024 eq_refl eq_refl Mx).

Lemma backward_round_correct:
  forall x_r (x:float),
    (if Rlt_bool
          (Rabs (round radix2 (FLT_exp (3 - 1024 - 53) 53)
                       (round_mode mode_NE) x_r))
          (bpow radix2 1024)
     then (x:R) = round radix2 (FLT_exp (3 - 1024 - 53) 53) (round_mode mode_NE) x_r /\
          is_finite x = true
     else (x:full_float) = binary_overflow 53 1024 mode_NE (Rlt_bool x_r 0)) ->
  forall x_ab, x ∈ γ x_ab ->
  x_r ∈ γ (backward_round x_ab).

Definition backward_floatofz (res_ab:fitv+⊤) (x_ab:zitv+⊤) : zitv+⊤+⊥ :=
  match res_ab with
  | Just res_ab' =>
    let 'FITV mx Mx := backward_round res_ab' in
    let low : option Zfast :=
      match Zoffloat_DN (Float.neg mx), x_ab with
      | Some m, Just (ZITV m' _) => Some (Zfastmax (Zfastopp m) m')
      | Some m, All => Some (Zfastopp m)
      | None, Just (ZITV m' _) => Some (m':Zfast)
      | None, All => None
      end
    in
    let high : option Zfast :=
      match Zoffloat_DN Mx, x_ab with
      | Some M, Just (ZITV _ M') => Some (Zfastmin M M')
      | Some M, All => Some (M:Zfast)
      | None, Just (ZITV _ M') => Some M'
      | None, All => None
      end
    in
    match low, high with
    | Some m, Some M =>
      if Zfastleb m M then NotBot (Just (ZITV m M))
      else Bot
    | _, _ => NotBot All
    end
  | All => NotBot All
  end.

Lemma backward_floatofz_correct :
  forall x x_ab, x ∈ γ x_ab ->
  forall res_ab, BofZ 53 1024 (@eq_refl _ Lt) (@eq_refl _ Lt) x ∈ γ res_ab ->
    x ∈ γ (backward_floatofz res_ab x_ab).

Definition half_backward_zoffloat (a:Zfast)
           prec emax Hprec Hmax : binary_float prec emax :=
  if Zfastleb a F0 then
    next_ulp prec emax Hprec Hmax
      (binary_normalize prec emax Hprec Hmax mode_DN (Zfastsub a F1) 0 false)
  else
    binary_normalize prec emax Hprec Hmax mode_UP a 0 false.

Lemma half_backward_zoffloat_correct :
  forall prec emax Hprec Hmax (f:binary_float prec emax) i a, (1 < prec)%Z ->
    ZofB f = Some i ->
    (a <= i)%Z ->
    Fleb (half_backward_zoffloat a _ _ Hprec Hmax) f = true.

Definition backward_zoffloat (res_ab:zitv+⊤) : fitv+⊤ :=
  let new_x_ab :=
    match res_ab with
    | All =>
      FITV (FF2B 53 1024 (binary_overflow 53 1024 mode_UP true) eq_refl)
           (FF2B 53 1024 (binary_overflow 53 1024 mode_DN false) eq_refl)
    | Just (ZITV a b) =>
      FITV (half_backward_zoffloat a 53 1024 eq_refl eq_refl)
           (Float.neg (half_backward_zoffloat (Zfastopp b) 53 1024 eq_refl eq_refl))
    end
  in
  Just new_x_ab.

Transparent Float.cmp.

Lemma backward_zoffloat_correct :
  forall x res res_ab, res ∈ γ res_ab ->
    ZofB x = Some res ->
    x ∈ γ (backward_zoffloat res_ab).

Definition backward_le_l (y:fitv+⊤) : fitv :=
  let My := match y with
            | All => B754_infinity false
            | Just (FITV _ My) => My
            end in
  FITV (B754_infinity true) My.

Lemma backward_le_l_correct :
  forall x y y_ab, y ∈ γ y_ab ->
    Fleb x y = true ->
    x ∈ γ (backward_le_l y_ab).

Definition backward_eq_l (x y:fitv+⊤) (res:bool) : fitv+⊤+⊥ :=
  if res then
    Just (backward_le_l y) ⊓ forward_fneg (Just (backward_le_l (forward_fneg y)))
  else
    match y with
    | All => NotBot All
    | Just (FITV my My) =>
      if Fleb My my then
        match x with
        | All => NotBot All
        | Just (FITV mx Mx) =>
          if Float.cmp Ceq my mx then
            if Fleb Mx mx then Bot
            else NotBot (Just (FITV (next_ulp 53 1024 eq_refl eq_refl mx) (B754_infinity false)))
          else
            if Float.cmp Ceq my Mx then
              NotBot (Just (FITV (B754_infinity true)
                                 (Float.neg (next_ulp 53 1024 eq_refl eq_refl (Float.neg Mx)))))
            else NotBot All
        end
      else NotBot All
    end.

Lemma backward_eq_l_correct :
  forall x x_ab, x ∈ γ x_ab ->
  forall y y_ab, y ∈ γ y_ab ->
  forall res, Float.cmp Ceq x y = res ->
    x ∈ γ (backward_eq_l x_ab y_ab res).

Definition backward_lt_l (y:fitv+⊤) : fitv :=
  let My := match y with
            | All => B754_infinity false
            | Just (FITV _ My) => My
            end in
  FITV (B754_infinity true)
       (Float.neg (next_ulp 53 1024 eq_refl eq_refl (Float.neg My))).

Lemma backward_lt_l_correct :
  forall x y y_ab, y ∈ γ y_ab ->
    Float.cmp Clt x y = true ->
    x ∈ γ (backward_lt_l y_ab).

Definition backward_not_ge_l (x y:fitv+⊤) : fitv+⊤ :=
  match x, y with
  | Just _, Just _ => Just (backward_lt_l y)
  | _, _ => All
  end.

Lemma backward_not_ge_l_correct :
  forall x x_ab, x ∈ γ x_ab ->
  forall y y_ab, y ∈ γ y_ab ->
    Fleb y x = false ->
    x ∈ γ (backward_not_ge_l x_ab y_ab).

Definition backward_not_gt_l (x y:fitv+⊤) : fitv+⊤ :=
  match x, y with
  | Just _, Just _ => Just (backward_le_l y)
  | _, _ => All
  end.

Lemma backward_not_gt_l_correct :
  forall x x_ab, x ∈ γ x_ab ->
  forall y y_ab, y ∈ γ y_ab ->
    Float.cmp Cgt x y = false ->
    x ∈ γ (backward_not_gt_l x_ab y_ab).

Definition backward_cmpf_l (c:comparison) (x y:fitv+⊤) (res:bool) : fitv+⊤+⊥ :=
  match c, res with
  | Ceq, _ => backward_eq_l x y res
  | Cne, _ => backward_eq_l x y (negb res)
  | Cle, true =>
    NotBot (Just (backward_le_l y))
  | Cge, false =>
    NotBot (backward_not_ge_l x y)
  | Cge, true =>
    NotBot (forward_fneg (Just (backward_le_l (forward_fneg y))))
  | Cle, false =>
    NotBot (forward_fneg (backward_not_ge_l (forward_fneg x) (forward_fneg y)))
  | Clt, true =>
    NotBot (Just (backward_lt_l y))
  | Cgt, true =>
    NotBot (forward_fneg (Just (backward_lt_l (forward_fneg y))))
  | Cgt, false =>
    NotBot (backward_not_gt_l x y)
  | Clt, false =>
    NotBot (forward_fneg (backward_not_gt_l (forward_fneg x) (forward_fneg y)))
  end.

Lemma Bcompare_fneg:
  ∀ x y:float, Bcompare _ _ x y = Bcompare _ _ (Float.neg y) (Float.neg x).

Lemma backward_cmpf_l_correct:
  ∀ x x_ab, x ∈ γ x_ab ->
  ∀ y y_ab, y ∈ γ y_ab ->
  ∀ c res, Float.cmp c x y = res ->
    x ∈ γ (backward_cmpf_l c x_ab y_ab res).

Definition backward_absf (res_ab:fitv+⊤) : fitv+⊤+⊥ :=
  match res_ab with
  | All => NotBot All
  | Just (FITV _ _) =>
    do res_ab' <- res_ab ⊓ (Just (FITV (B754_zero false) (B754_infinity false)));
    forward_fneg res_ab' ⊔ res_ab'
  end.

Lemma backward_absf_correct :
  forall x res_ab, Float.abs x ∈ γ res_ab ->
    x ∈ γ (backward_absf res_ab).

Definition backward_addf_l (res y:fitv+⊤) : fitv+⊤+⊥ :=
  match y with
  | All => NotBot All
  | Just (FITV my My) =>
    match res with
    | All => NotBot All
    | Just res =>
      let 'FITV mres Mres := backward_round res in
      let m:= Bminus 53 1024 (@eq_refl _ Lt) (@eq_refl _ Lt) Float.binop_pl mode_UP mres My in
      let M:= Bminus 53 1024 (@eq_refl _ Lt) (@eq_refl _ Lt) Float.binop_pl mode_DN Mres my in
      if is_nan m || is_nan M then NotBot All
      else NotBot (Just (FITV m M))
    end
  end.

Lemma backward_addf_l_correct:
  ∀ x y y_ab, y ∈ γ y_ab ->
  ∀ res_ab, (Float.add x y)%Z ∈ γ res_ab ->
    x ∈ γ (backward_addf_l res_ab y_ab).

Definition backward_addf_r := backward_addf_l.

Lemma backward_addf_r_correct:
  ∀ x x_ab, x ∈ γ x_ab ->
  ∀ y res_ab, (Float.add x y)%Z ∈ γ res_ab ->
    y ∈ γ (backward_addf_r res_ab x_ab).

Definition backward_subf_l (res y:fitv+⊤) : fitv+⊤+⊥ := backward_addf_l res (forward_fneg y).

Lemma backward_subf_l_correct:
  ∀ x y y_ab, y ∈ γ y_ab ->
  ∀ res_ab, (Float.sub x y)%Z ∈ γ res_ab ->
    x ∈ γ (backward_subf_l res_ab y_ab).

Definition backward_subf_r (res x:fitv+⊤) : fitv+⊤+⊥ :=
  fmap forward_fneg (backward_addf_l res x).

Lemma backward_subf_r_correct:
  ∀ x x_ab, x ∈ γ x_ab ->
  ∀ y res_ab, (Float.sub x y)%Z ∈ γ res_ab ->
    y ∈ γ (backward_subf_r res_ab x_ab).

Definition binary32offloat_DN (f: float) : binary32 :=
  Bconv 53 1024 24 128 eq_refl eq_refl Float.to_single_pl mode_DN f.

Lemma binary32offloat_DN_Fleb :
  ∀ f1 f2,
    Fleb (Float.of_single f1) f2 = true ->
    Fleb f1 (binary32offloat_DN f2) = true.

Definition backward_round_single (x:fitv) : fitv :=
  let 'FITV mx Mx := x in
  FITV (Float.neg (Float.of_single
                     (next_ulp 24 128 eq_refl eq_refl (binary32offloat_DN (Float.neg mx)))))
       (Float.of_single (next_ulp 24 128 eq_refl eq_refl (binary32offloat_DN Mx))).

Lemma backward_round_single_correct :
  forall x_r (x:float32),
    (if Rlt_bool
          (Rabs (round radix2 (FLT_exp (3 - 128 - 24) 24)
                       (round_mode mode_NE) x_r))
          (bpow radix2 128)
     then (x:R) = round radix2 (FLT_exp (3 - 128 - 24) 24) (round_mode mode_NE) x_r /\
          is_finite x = true
     else (x:full_float) = binary_overflow 24 128 mode_NE (Rlt_bool x_r 0)) ->
  forall x_ab, (Float.of_single x) ∈ γ x_ab ->
  gammaR (backward_round_single x_ab) x_r.

Definition backward_singleoffloat (res_ab:fitv+⊤) : fitv+⊤+⊥ :=
  match res_ab with
  | All => NotBot All
  | Just res_ab => NotBot (Just (backward_round_single res_ab))
  end.

Lemma backward_singleoffloat_correct :
  forall x res_ab, Float.of_single (Float.to_single x) ∈ γ res_ab ->
    x ∈ γ (backward_singleoffloat res_ab).

Definition backward_singleofz (res_ab:fitv+⊤) (x_ab:zitv+⊤) : zitv+⊤+⊥ :=
  match res_ab with
  | All => NotBot x_ab
  | Just res_ab =>
    let 'FITV mx Mx := backward_round_single res_ab in
    let low :=
      match Zoffloat_DN (Float.neg mx), x_ab with
      | Some m, Just (ZITV m' _) => Some (Zfastmax (Zfastopp m) m')
      | Some m, All => Some (Zfastopp m)
      | None, Just (ZITV m' _) => Some (m')
      | None, All => None
      end
    in
    let high :=
      match Zoffloat_DN Mx, x_ab with
      | Some M, Just (ZITV _ M') => Some (Zfastmin M M')
      | Some M, All => Some (M:Zfast)
      | None, Just (ZITV _ M') => Some (M')
      | None, All => None
      end
    in
    match low, high with
    | Some m, Some M =>
      if Zfastleb m M then NotBot (Just (ZITV m M))
      else Bot
    | _, _ => NotBot All
    end
  end.

Lemma backward_singleofz_correct :
  forall x x_ab, x ∈ γ x_ab ->
  forall res_ab, Float.of_single (BofZ 24 128 (@eq_refl _ Lt) (@eq_refl _ Lt) x) ∈ γ res_ab ->
    x ∈ γ (backward_singleofz res_ab x_ab).