Library Float.Expansions.EFast2Sum

Require Export Fast2Sum.

Section EFast.
Variable b : Fbound.
Variable precision : nat.

Let radix := 2%Z.

Let radixMoreThanOne : (1 < radix)%Z.
unfold radix in |- *; red in |- *; simpl in |- *; auto.
Qed.

Let radixMoreThanZERO := Zlt_1_O _ (Zlt_le_weak _ _ radixMoreThanOne).
Hint Resolve radixMoreThanZERO radixMoreThanOne: zarith.

Coercion Local FtoRradix := FtoR radix.
Hypothesis precisionGreaterThanOne : 1 < precision.
Hypothesis pGivesBound : Zpos (vNum b) = Zpower_nat radix precision.
Variable Iplus : float -> float -> float.
Hypothesis
  IplusCorrect :
    forall p q : float,
    Fbounded b p -> Fbounded b q -> Closest b radix (p + q) (Iplus p q).
Hypothesis
  IplusComp :
    forall p q r s : float,
    Fbounded b p ->
    Fbounded b q ->
    Fbounded b r ->
    Fbounded b s -> p = r :>R -> q = s :>R -> Iplus p q = Iplus r s :>R.
Hypothesis IplusSym : forall p q : float, Iplus p q = Iplus q p.
Hypothesis
  IplusOp : forall p q : float, Fopp (Iplus p q) = Iplus (Fopp p) (Fopp q).
Variable Iminus : float -> float -> float.
Hypothesis IminusPlus : forall p q : float, Iminus p q = Iplus p (Fopp q).

Theorem IminusComp :
 forall p q r s : float,
 Fbounded b p ->
 Fbounded b q ->
 Fbounded b r ->
 Fbounded b s -> p = r :>R -> q = s :>R -> Iminus p q = Iminus r s :>R.
intros p q r s H' H'0 H'1 H'2 H'3 H'4.
repeat rewrite IminusPlus.
apply IplusComp; auto; try apply oppBounded; auto.
unfold FtoRradix in |- *; repeat rewrite Fopp_correct; auto.
unfold FtoRradix in H'4; rewrite H'4; auto.
Qed.

Theorem EvenBound :
 forall (p : float) (m : Z),
 Even m ->
 (Zpred (Zpower_nat radix precision) <= m)%Z ->
 (m <= Zpower_nat radix (S precision) - radix)%Z ->
 Fbounded b p -> exists q : float, Fbounded b q /\ q = Float m (Fexp p) :>R.
intros p m H' H'0 H'1 H'2.
cut (0 < radix)%Z; [ intros Hr | unfold radix in |- *; auto with zarith ].
repeat split; simpl in |- *; auto with float zarith.
case H'; intros am2 H'3.
exists (Float am2 (Zsucc (Fexp p))); split.
repeat split; simpl in |- *; auto with float zarith.
apply Zmult_gt_0_lt_reg_r with (p := radix); auto with zarith.
rewrite Zmult_comm.
pattern radix at 1 in |- *; rewrite <- (Zabs_eq radix); auto with zarith.
rewrite <- Zabs_Zmult.
unfold radix in |- *; rewrite <- H'3.
rewrite Zabs_eq; fold radix in |- *; auto.
apply Zle_lt_trans with (1 := H'1).
rewrite pGivesBound; replace (S precision) with (precision + 1);
 [ idtac | ring ].
rewrite Zpower_nat_is_exp; auto; rewrite Zpower_nat_1; unfold Zminus in |- *;
 auto with zarith.
apply Zle_trans with (2 := H'0); auto with zarith arith.
unfold FtoRradix, FtoR in |- *; simpl in |- *; rewrite powerRZ_Zs;
 auto with real.
rewrite H'3; rewrite Rmult_IZR; simpl in |- *.
ring.
Qed.

Theorem ExtMDekker1 :
 forall p q : float,
 Fbounded b p ->
 Fbounded b q ->
 (Fexp q <= Fexp p)%Z ->
 (0 <= p)%R -> Iminus (Iplus p q) p = (Iplus p q - p)%R :>R.
intros p q H' H'0 H'1 H'2.
cut (0 <= Fnum p)%Z;
 [ intros Z1 | apply LeR0Fnum with (radix := radix); auto ].
case (Rle_or_lt (Rabs q) (Rabs p)); intros Rl0.
apply MDekker with (b := b) (precision := precision); auto.
rewrite Rabs_right in Rl0; auto with real.
case (LessExpBound _ TwoMoreThanOne b p (Fabs q)); auto.
apply absFBounded; auto.
unfold FtoRradix in |- *; rewrite Fabs_correct; auto; apply Rlt_le; auto.
replace (Fexp (Fabs q)) with (Fexp q);
 [ idtac | case q; simpl in |- *; auto ].
intros m E; elim E; fold radix in |- *; intros H'4 H'5; clear E.
cut (0 <= Fnum (Float m (Fexp q)))%Z;
 [ simpl in |- *; intros Z2
 | apply LeR0Fnum with (radix := radix); auto with zarith;
    unfold FtoRradix in H'4; rewrite H'4 ]; auto.
unfold FtoRradix in |- *; rewrite <- H'4; auto.
apply
 trans_eq
  with
    (y := FtoRradix (Iminus (Iplus (Float m (Fexp q)) q) (Float m (Fexp q)))).
apply IminusComp; auto.
apply IplusBounded; auto.
apply IplusBounded; auto.
fold FtoRradix in |- *.
rewrite (IplusComp p q (Float m (Fexp q)) q); auto.
case (Rle_or_lt 0 q); intros Rl1.
cut (0 <= Fnum q)%Z;
 [ intros Z3 | apply LeR0Fnum with (radix := radix); auto ].
rewrite Rabs_right in Rl0; auto with real.
case (Zle_or_lt (Zpos (vNum b)) (m + Fnum q)); intros Rl2.
2: apply (MDekkerAux3 b); auto.
2: unfold Fplus in |- *; simpl in |- *.
2: repeat rewrite Zmin_n_n; repeat rewrite <- Zminus_diag_reverse;
    simpl in |- *.
2: rewrite Zpower_nat_O; repeat rewrite Zmult_1_r.
2: repeat split; simpl in |- *; auto with float zarith.
2: rewrite Zabs_eq; auto with zarith.
cut (m + Fnum q <= Zpower_nat radix (S precision) - radix)%Z;
 [ intros Rl4 | idtac ].
2: replace (Zpower_nat radix (S precision) - radix)%Z with
    (Zpred (Zpower_nat radix precision) + Zpred (Zpower_nat radix precision))%Z;
    auto with zarith.
2: apply Zplus_le_compat; auto.
2: apply Zle_Zpred; rewrite <- (Zabs_eq m); auto with float zarith.
2: rewrite <- pGivesBound; case H'5; simpl in |- *; auto.
2: apply Zle_Zpred; rewrite <- (Zabs_eq (Fnum q)); auto with float zarith.
2: rewrite <- pGivesBound; case H'0; auto.
2: replace (S precision) with (1 + precision);
    [ idtac | simpl in |- *; auto ].
2: rewrite Zpower_nat_is_exp; auto with zarith; rewrite Zpower_nat_1.
2: replace radix with (1 + 1)%Z;
    [ idtac | unfold radix in |- *; simpl in |- *; auto ].
2: unfold Zpred in |- *; replace (-1)%Z with (- (1))%Z;
    [ ring | simpl in |- *; auto ].
case (OddEvenDec (m + Fnum q)); intros E0.
case
 (RoundedModeRep b radix precision)
  with
    (P := Closest b radix)
    (p := Fplus radix (Float m (Fexp q)) q)
    (q := Iplus (Float m (Fexp q)) q); auto.
apply ClosestRoundedModeP with (precision := precision); auto.
rewrite Fplus_correct; auto.
replace (Fexp (Fplus radix (Float m (Fexp q)) q)) with (Fexp q).
simpl in |- *; intros m' Hm'.
2: unfold Fplus in |- *; simpl in |- *.
2: repeat rewrite Zmin_n_n; repeat rewrite <- Zminus_diag_reverse;
    simpl in |- *; auto.
cut (m' - m <= Zsucc (Fnum q))%Z; [ intros Rl5 | idtac ].
case (Zle_lt_or_eq _ _ Rl5); clear Rl5; intros Rl5.
cut (Zpred (Fnum q) <= m' - m)%Z; [ intros Rl6 | idtac ].
unfold FtoRradix in |- *; rewrite Hm'.
unfold FtoRradix in |- *; rewrite <- Fminus_correct; auto.
apply sym_eq; apply (ClosestIdem b radix); auto.
unfold Fminus, Fplus in |- *; simpl in |- *.
repeat rewrite Zmin_n_n; repeat rewrite <- Zminus_diag_reverse; simpl in |- *.
rewrite Zpower_nat_O; repeat rewrite Zmult_1_r.
repeat split; simpl in |- *; auto with float zarith.
apply Zlt_Zabs_intro.
apply Zlt_le_trans with (Zpred 0); auto with zarith.
replace (Zpred 0) with (- (1))%Z; auto with zarith.
apply Zlt_Zopp; apply (vNumbMoreThanOne radix) with (precision := precision);
 auto with zarith.
apply Zle_trans with (2 := Rl6); auto with zarith.
apply Zle_lt_trans with (Fnum q); auto with zarith.
rewrite <- (Zabs_eq (Fnum q)); auto; case H'0; auto.
rewrite Fminus_correct; auto.
rewrite <- Hm'; auto.
apply (IminusCorrect b Iplus); auto.
apply IplusBounded; auto.
apply Zplus_le_reg_l with (p := m); auto.
replace (m + (m' - m))%Z with m'; auto with zarith.
apply le_IZR.
apply Rmult_le_reg_l with (r := powerRZ radix (Fexp q));
 auto with real zarith.
repeat rewrite (Rmult_comm (powerRZ radix (Fexp q))).
change (Float (m + Zpred (Fnum q)) (Fexp q) <= Float m' (Fexp q))%R in |- *.
replace (m + Zpred (Fnum q))%Z with (Zpred (m + Fnum q)).
2: unfold Zpred in |- *; ring.
case (EvenBound (Float m (Fexp q)) (Zpred (m + Fnum q))); auto.
apply OddSEvenInv; auto.
rewrite <- Zsucc_pred; auto.
case (Zle_lt_or_eq 0 (m + Fnum q)).
replace 0%Z with (0 + 0)%Z; auto with zarith.
intros H'3; apply Zlt_succ_le.
cut (forall r : Z, Zsucc (Zpred r) = r);
 [ intros Er; rewrite Er; auto
 | intros r; unfold Zsucc, Zpred in |- *; ring ].
rewrite <- pGivesBound; auto with zarith.
intros H'3; Contradict E0; rewrite <- H'3; simpl in |- *.
apply EvenNOdd; apply EvenO.
apply Zle_trans with (2 := Rl4); auto with zarith.
simpl in |- *; intros x H'3; elim H'3; intros H'6 H'7; rewrite <- H'7;
 clear H'3.
unfold FtoRradix in |- *; rewrite <- Hm'.
apply (ClosestMonotone b radix x (Float m (Fexp q) + q)%R); auto.
unfold FtoRradix in |- *; rewrite <- Fplus_correct; auto.
unfold FtoRradix in H'7; rewrite H'7.
unfold Fplus in |- *.
repeat rewrite Zmin_n_n; repeat rewrite <- Zminus_diag_reverse; simpl in |- *;
 auto.
rewrite Zpower_nat_O; repeat rewrite Zmult_1_r.
unfold FtoR in |- *; simpl in |- *.
repeat rewrite (fun x : Z => Rmult_comm x (powerRZ 2 (Fexp q))).
apply Rmult_lt_compat_l; auto with real arith.
apply Rlt_IZR; auto.
apply Zlt_pred; auto.
apply (RoundedModeProjectorIdem b radix (Closest b radix)); auto.
apply ClosestRoundedModeP with (precision := precision); auto.
case H'0; rewrite Zabs_eq; auto with zarith.
intros Hq Hq'; case (Zlt_next _ _ Hq); intros H'3.
2: cut
    (Closest b radix
       (FtoR radix (Iplus (Float m (Fexp q)) q) -
        FtoR radix (Float m (Fexp q)))
       (Iminus (Iplus (Float m (Fexp q)) q) (Float m (Fexp q)))).
2: unfold FtoRradix in |- *; rewrite Hm'.
2: rewrite <- Fminus_correct; auto.
2: unfold Fminus, Fplus in |- *; simpl in |- *; auto.
2: repeat rewrite Zmin_n_n; repeat rewrite <- Zminus_diag_reverse;
    simpl in |- *; auto.
2: rewrite Zpower_nat_O; repeat rewrite Zmult_1_r.
2: replace (m' + - m)%Z with (m' - m)%Z; [ idtac | ring ].
2: intros H'6; apply sym_eq; apply (ClosestIdem b radix); auto.
2: repeat split; simpl in |- *; auto with float zarith.
2: rewrite Rl5; rewrite Zabs_eq; auto with zarith.
2: apply (IminusCorrect b Iplus); auto.
2: apply IplusBounded; auto.
cut
 (Closest b radix
    (FtoR radix (Iplus (Float m (Fexp q)) q) - FtoR radix (Float m (Fexp q)))
    (Iminus (Iplus (Float m (Fexp q)) q) (Float m (Fexp q)))).
unfold FtoRradix in |- *; rewrite Hm'.
rewrite <- Fminus_correct; auto.
unfold Fminus, Fplus in |- *; simpl in |- *; auto.
repeat rewrite Zmin_n_n; repeat rewrite <- Zminus_diag_reverse; simpl in |- *;
 auto.
rewrite Zpower_nat_O; repeat rewrite Zmult_1_r.
replace (m' + - m)%Z with (m' - m)%Z; [ idtac | ring ].
rewrite Rl5.
rewrite <- H'3.
replace (FtoR radix (Float (Zpos (vNum b)) (Fexp q))) with
 (FtoR radix (Float (nNormMin radix precision) (Zsucc (Fexp q)))).
intros H'6; apply sym_eq; apply (ClosestIdem b radix); auto.
repeat split; simpl in |- *; auto with float arith.
rewrite Zabs_eq; auto with zarith arith.
apply ZltNormMinVnum; auto with zarith.
unfold nNormMin in |- *; auto with zarith.
apply Zle_trans with (Fexp q); auto with float zarith.
rewrite (PosNormMin radix b precision); unfold nNormMin in |- *;
 auto with zarith.
rewrite Zmult_comm; unfold FtoR in |- *; simpl in |- *;
 rewrite powerRZ_Zs; auto with real.
rewrite Rmult_IZR; auto with real.
apply (IminusCorrect b Iplus); auto.
apply IplusBounded; auto.
apply Zplus_le_reg_l with (p := m); auto.
replace (m + (m' - m))%Z with m'; auto with zarith.
apply le_IZR.
apply Rmult_le_reg_l with (r := powerRZ radix (Fexp q)); auto with real arith.
repeat rewrite (Rmult_comm (powerRZ radix (Fexp q))).
change (Float m' (Fexp q) <= Float (m + Zsucc (Fnum q)) (Fexp q))%R in |- *.
replace (m + Zsucc (Fnum q))%Z with (Zsucc (m + Fnum q));
 [ idtac | unfold Zsucc in |- *; ring ].
case (EvenBound (Float m (Fexp q)) (Zsucc (m + Fnum q))); auto.
apply OddSEven; auto.
apply Zle_trans with (m + Fnum q)%Z; auto with zarith.
case (Zle_next _ _ Rl4); auto.
intros H'3; Contradict E0.
rewrite <- H'3; auto.
apply EvenNOdd; auto.
apply EvenPlusInv1 with (n := radix); auto.
replace (radix + (Zpower_nat radix (S precision) - radix))%Z with
 (Zpower_nat radix (S precision)); [ idtac | ring ].
apply EvenExp; auto.
exists 1%Z; unfold radix in |- *; ring.
exists 1%Z; unfold radix in |- *; ring.
simpl in |- *; intros x H'3; elim H'3; intros H'6 H'7; rewrite <- H'7;
 clear H'3.
unfold FtoRradix in |- *; rewrite <- Hm'.
apply (ClosestMonotone b radix (Float m (Fexp q) + q)%R x); auto.
rewrite H'7.
unfold FtoRradix in |- *; rewrite <- Fplus_correct; auto.
unfold Fplus in |- *.
repeat rewrite Zmin_n_n; repeat rewrite <- Zminus_diag_reverse; simpl in |- *;
 auto.
rewrite Zpower_nat_O; repeat rewrite Zmult_1_r.
unfold FtoR in |- *; simpl in |- *.
repeat rewrite (fun x : Z => Rmult_comm x (powerRZ 2 (Fexp q))).
apply Rmult_lt_compat_l; auto with real arith.
apply Rlt_IZR; auto with zarith.
apply (RoundedModeProjectorIdem b radix (Closest b radix)); auto.
apply ClosestRoundedModeP with (precision := precision); auto.
apply (MDekkerAux2 b); auto.
cut (Closest b radix (Float m (Fexp q) + q) (Iplus (Float m (Fexp q)) q));
 auto.
unfold FtoRradix in |- *; rewrite <- Fplus_correct; auto with arith.
unfold Fplus in |- *.
repeat rewrite Zmin_n_n; repeat rewrite <- Zminus_diag_reverse; simpl in |- *;
 auto.
rewrite Zpower_nat_O; repeat rewrite Zmult_1_r.
case (EvenBound (Float m (Fexp q)) (m + Fnum q)); auto.
apply Zle_trans with (2 := Rl2); auto with zarith arith.
intros x (H'6, H'7); fold radix in |- *; simpl in H'7.
replace (FtoR radix (Float m (Fexp q)) + FtoR radix q)%R with
 (FtoR radix (Float (m + Fnum q) (Fexp q))).
fold FtoRradix in |- *; rewrite <- H'7.
intros H'3; apply sym_eq; apply (ClosestIdem b radix); auto with arith.
unfold FtoRradix, FtoR in |- *; simpl in |- *; rewrite plus_IZR; ring.
apply (MDekkerAux3 b); fold radix in |- *; auto.
replace (Fplus radix (Float m (Fexp q)) q) with
 (Fminus radix (Float m (Fexp q)) (Fopp q)).
apply BminusSameExp; auto.
apply oppBounded; auto.
unfold FtoRradix in H'4; rewrite H'4; auto.
replace 0%R with (-0)%R; auto with real.
rewrite Fopp_correct; apply Rlt_le; auto with real.
unfold Fminus in |- *; rewrite Fopp_Fopp; auto.
Qed.

Theorem ExtMDekker :
 forall p q : float,
 Fbounded b p ->
 Fbounded b q ->
 (Fexp q <= Fexp p)%Z -> Iminus (Iplus p q) p = (Iplus p q - p)%R :>R.
intros p q H' H'0 H'1.
case (Rle_or_lt 0 p); intros H1.
apply ExtMDekker1; auto.
apply (MDekkerAux5 b); auto.
apply ExtMDekker1; auto.
apply oppBounded; auto.
apply oppBounded; auto.
apply Rlt_le; unfold FtoRradix in |- *; rewrite Fopp_correct;
 replace 0%R with (-0)%R; auto with real.
Qed.

Theorem ExtDekker :
 forall p q : float,
 Fbounded b p ->
 Fbounded b q ->
 (Fexp q <= Fexp p)%Z ->
 Iminus q (Iminus (Iplus p q) p) = (p + q - Iplus p q)%R :>R.
intros p q H' H'0 H'1.
apply (MDekkerAux1 b precision); auto.
apply ExtMDekker; auto.
Qed.
End EFast.