Library Float.Expansions.Fast2Sum

Require Export AllFloat.
Section Fast.
Variable b : Fbound.
Variable precision : nat.

Let radix := 2%Z.

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 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).

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.

Theorem IminusCorrect :
 forall p q : float,
 Fbounded b p -> Fbounded b q -> Closest b radix (p - q) (Iminus p q).
intros p q H' H'0.
rewrite IminusPlus.
unfold Rminus in |- *.
unfold FtoRradix in |- *; rewrite <- Fopp_correct.
apply IplusCorrect; auto.
apply oppBounded; auto.
Qed.
Hint Resolve IminusCorrect.

Theorem ErrorBoundedIplus :
 forall p q : float,
 Fbounded b p ->
 Fbounded b q ->
 exists error : float, error = (p + q - Iplus p q)%R :>R /\ Fbounded b error.
intros p q H' H'0.
case
 (errorBoundedPlus b radix precision)
  with (p := p) (q := q) (pq := Iplus p q); auto.
intros x H'1; elim H'1; intros H'2 H'3; elim H'3; intros H'4 H'5; exists x;
 auto.
Qed.

Theorem IplusOr :
 forall p q : float,
 Fbounded b p -> Fbounded b q -> q = 0%R :>R -> Iplus p q = p :>R.
intros p q H' H'0 H'1.
cut (Closest b radix (p + q) (Iplus p q)); auto.
rewrite H'1.
rewrite Rplus_0_r.
intros H'2; apply sym_eq.
unfold FtoRradix in |- *; apply (ClosestIdem b radix); auto.
Qed.

Theorem IminusId :
 forall p q : float,
 Fbounded b p -> Fbounded b q -> p = q :>R -> Iminus p q = 0%R :>R.
intros p q H' H'0 H'1.
cut (Closest b radix (p - q) (Iminus p q)); auto.
replace (p - q)%R with 0%R; [ idtac | rewrite <- H'1; ring ].
replace 0%R with (FtoRradix (Fzero (- dExp b))).
intros H'2; apply sym_eq.
unfold FtoRradix in |- *; apply (ClosestIdem b radix); auto.
apply FboundedFzero; auto.
unfold Fzero, FtoRradix, FtoR in |- *; simpl in |- *; ring.
Qed.

Theorem IminusOl :
 forall p q : float,
 Fbounded b p -> Fbounded b q -> p = 0%R :>R -> Iminus p q = (- q)%R :>R.
intros p q H' H'0 H'1.
cut (Closest b radix (p - q) (Iminus p q)); auto.
rewrite H'1.
replace (0 - q)%R with (FtoRradix (Fopp q)).
intros H'2; apply sym_eq.
unfold FtoRradix in |- *; rewrite <- Fopp_correct;
 apply (ClosestIdem b radix); auto.
apply oppBounded; auto.
unfold FtoRradix in |- *; rewrite Fopp_correct; auto; ring.
Qed.

Theorem IplusBounded :
 forall p q : float, Fbounded b p -> Fbounded b q -> Fbounded b (Iplus p q).
intros p q H' H'0; cut (Closest b radix (p + q) (Iplus p q)); auto.
intros H1; case H1; auto.
Qed.

Theorem IminusBounded :
 forall p q : float, Fbounded b p -> Fbounded b q -> Fbounded b (Iminus p q).
intros p q H' H'0; cut (Closest b radix (p - q) (Iminus p q)); auto.
intros H1; case H1; auto.
Qed.

Theorem IminusInv :
 forall p q r s : float,
 Fbounded b p ->
 Fbounded b q ->
 Fbounded b r ->
 Fbounded b s -> p = s :>R -> r = (s - q)%R :>R -> Iminus p q = r :>R.
intros p q r s H' H'0 H'1 H'2 H'3 H'4.
cut (Closest b radix (p - q) (Iminus p q)); auto.
rewrite H'3.
rewrite <- H'4.
intros H'5; apply sym_eq.
unfold FtoRradix in |- *; apply (ClosestIdem b radix); auto.
Qed.

Theorem IminusFminus :
 forall p q : float,
 Fbounded b p ->
 Fbounded b q ->
 Fbounded b (Fminus radix p q) -> Iminus p q = Fminus radix p q :>R.
intros p q H' H'0 H'1.
cut (Closest b radix (p - q) (Iminus p q)); auto.
intros H'2.
apply sym_eq; apply (ClosestIdem b radix); auto.
rewrite Fminus_correct; auto.
Qed.

Theorem MDekkerAux1 :
 forall p q : float,
 Iminus (Iplus p q) p = (Iplus p q - p)%R :>R ->
 Fbounded b p ->
 Fbounded b q -> Iminus q (Iminus (Iplus p q) p) = (p + q - Iplus p q)%R :>R.
intros p q H' H'0 H'1.
elim (ErrorBoundedIplus p q);
 [ intros error E; elim E; intros H'7 H'8; clear E | idtac | idtac ];
 auto.
cut
 (Closest b radix (q - Iminus (Iplus p q) p)
    (Iminus q (Iminus (Iplus p q) p))); auto.
rewrite H'.
replace (q - (Iplus p q - p))%R with (p + q - Iplus p q)%R; [ idtac | ring ].
rewrite <- H'7.
intros H'2.
apply sym_eq; apply (ClosestIdem b radix); auto.
apply IminusCorrect; auto.
repeat (try apply IplusBounded; try apply IminusBounded; auto); auto.
Qed.

Theorem MDekkerAux2 :
 forall p q : float,
 Iplus p q = (p + q)%R :>R ->
 Fbounded b p -> Fbounded b q -> Iminus (Iplus p q) p = (Iplus p q - p)%R :>R.
intros p q H' H'0 H'1.
cut (Closest b radix (Iplus p q - p) (Iminus (Iplus p q) p)); auto.
repeat rewrite H'.
replace (p + q - p)%R with (FtoRradix q); [ idtac | ring ].
intros H'2.
apply sym_eq; apply (ClosestIdem b radix); auto.
apply IminusCorrect; auto.
apply IplusBounded; auto.
Qed.

Theorem MDekkerAux3 :
 forall p q : float,
 Fbounded b (Fplus radix p q) ->
 Fbounded b p -> Fbounded b q -> Iminus (Iplus p q) p = (Iplus p q - p)%R :>R.
intros p q H' H'0 H'1.
apply MDekkerAux2; auto.
unfold FtoRradix in |- *; rewrite <- Fplus_correct; auto.
apply sym_eq; apply (ClosestIdem b radix); auto.
rewrite Fplus_correct; auto.
Qed.

Theorem MDekkerAux4 :
 forall p q : float,
 Fbounded b (Fminus radix (Iplus p q) p) ->
 Fbounded b p -> Fbounded b q -> Iminus (Iplus p q) p = (Iplus p q - p)%R :>R.
intros p q H' H'0 H'1.
unfold FtoRradix in |- *; rewrite <- Fminus_correct; auto.
apply sym_eq; apply (ClosestIdem b radix); auto.
rewrite Fminus_correct; auto.
apply IminusCorrect; auto.
repeat (try apply IplusBounded; try apply IminusBounded; auto); auto.
Qed.

Theorem Dekker1 :
 forall p q : float,
 (0 <= q)%R ->
 (q <= p)%R ->
 Fbounded b p -> Fbounded b q -> Iminus (Iplus p q) p = (Iplus p q - p)%R :>R.
intros p q H' H'0 H'1 H'2.
apply MDekkerAux4; auto.
apply oppBoundedInv; auto.
rewrite Fopp_Fminus; auto.
apply Sterbenz; auto.
apply IplusBounded; auto.
apply Rmult_le_reg_l with (r := INR 2); auto with real arith.
rewrite <- Rmult_assoc; rewrite Rinv_r; auto with real arith;
 rewrite Rmult_1_l.
apply (RoundedModeMult b radix) with (P := Closest b radix) (r := (p + q)%R);
 auto.
apply ClosestRoundedModeP with (precision := precision); auto.
replace (radix * FtoR radix p)%R with (p + p)%R;
 [ idtac | simpl in |- *; fold FtoRradix; ring ]; auto with real.
case H'; clear H'; intros H'.
apply Rmult_le_reg_l with (r := (/ radix)%R); auto with real.
rewrite <- Rmult_assoc; rewrite Rinv_l; auto with real; rewrite Rmult_1_l.
apply (FmultRadixInv b radix precision) with (y := (p + q)%R); auto.
apply Rmult_lt_reg_l with (r := INR 2); auto with real.
rewrite <- Rmult_assoc; rewrite Rinv_r; auto with real; rewrite Rmult_1_l.
apply Rle_lt_trans with (radix * p)%R.
apply Rledouble; auto.
apply Rlt_le; apply Rlt_le_trans with (FtoRradix q); auto.
apply Rmult_lt_reg_l with (r := (/ radix)%R); auto with real.
repeat rewrite <- Rmult_assoc; repeat rewrite Rinv_l; auto with real zarith;
 repeat rewrite Rmult_1_l.
pattern (FtoRradix p) at 1 in |- *; replace (FtoRradix p) with (p + 0)%R;
 [ idtac | ring ]; auto with real.
rewrite IplusOr; auto.
apply Rledouble; auto.
rewrite H'; auto.
Qed.

Theorem Dekker2 :
 forall p q : float,
 (0 <= p)%R ->
 (- q <= p)%R ->
 (p <= radix * - q)%R ->
 Fbounded b p -> Fbounded b q -> Iminus (Iplus p q) p = (Iplus p q - p)%R :>R.
intros p q H' H'0 H'1 H'2 H'3.
apply MDekkerAux3; auto.
rewrite <- (Fopp_Fopp q).
change (Fbounded b (Fminus radix p (Fopp q))) in |- *.
apply Sterbenz; auto.
apply oppBounded; auto.
apply Rmult_le_reg_l with (r := INR 2); auto with real arith.
rewrite <- Rmult_assoc; rewrite Rinv_r; auto with real arith;
 rewrite Rmult_1_l; auto.
rewrite Fopp_correct; auto.
apply Rle_trans with (FtoRradix p); auto.
apply Rledouble; auto.
rewrite Fopp_correct; auto.
Qed.

Theorem Dekker3 :
 forall p q : float,
 (q <= 0)%R ->
 (radix * - q < p)%R ->
 Fbounded b p -> Fbounded b q -> Iminus (Iplus p q) p = (Iplus p q - p)%R :>R.
intros p q H' H'0 H'1 H'2.
apply MDekkerAux4; auto.
apply Sterbenz; auto.
apply IplusBounded; auto.
apply (FmultRadixInv b radix precision) with (y := (p + q)%R); auto.
apply Rmult_lt_reg_l with (r := INR 2); auto with real arith.
rewrite <- Rmult_assoc; rewrite Rinv_r; auto with real arith;
 rewrite Rmult_1_l.
replace (2%nat * (p + q))%R with (2%nat * p + 2%nat * q)%R;
 [ idtac | simpl in |- *; ring ].
apply Rplus_lt_reg_r with (r := (- FtoR radix p)%R).
replace (- FtoR radix p + FtoR radix p)%R with 0%R;
 [ idtac | simpl in |- *; ring ].
replace (- FtoR radix p + (2%nat * p + 2%nat * q))%R with (p + 2%nat * q)%R;
 [ idtac | simpl in |- *; unfold FtoRradix in |- *; ring ].
apply Rplus_lt_reg_r with (r := (2%nat * - q)%R).
replace (2%nat * - q + 0)%R with (2%nat * - q)%R;
 [ idtac | simpl in |- *; ring ].
replace (2%nat * - q + (p + 2%nat * q))%R with (FtoRradix p);
 [ idtac | simpl in |- *; ring ]; auto.
apply (RoundedModeMult b radix) with (P := Closest b radix) (r := (p + q)%R);
 auto.
apply ClosestRoundedModeP with (precision := precision); auto.
apply Rle_trans with (FtoR radix p); auto.
replace (FtoR radix p) with (p + 0)%R; [ idtac | fold FtoRradix; ring ].
apply Rplus_le_compat_l; auto.
apply Rledouble; auto.
apply Rlt_le; auto.
apply Rle_lt_trans with (2 := H'0).
apply Rle_trans with (- q)%R; auto.
replace 0%R with (-0)%R; auto with real.
apply Rledouble; auto.
replace 0%R with (-0)%R; auto with real.
Qed.

Theorem MDekkerAux5 :
 forall p q : float,
 Fbounded b p ->
 Fbounded b q ->
 Iminus (Iplus (Fopp p) (Fopp q)) (Fopp p) =
 (Iplus (Fopp p) (Fopp q) - Fopp p)%R :>R ->
 Iminus (Iplus p q) p = (Iplus p q - p)%R :>R.
intros p q H' H'0 H'1.
rewrite <- (Fopp_Fopp (Iminus (Iplus p q) p)); auto.
repeat rewrite IminusPlus; auto.
rewrite IplusOp; auto.
rewrite <- IminusPlus.
rewrite IplusOp; auto.
unfold FtoRradix in |- *; rewrite Fopp_correct.
unfold FtoRradix in H'1; rewrite H'1.
rewrite <- IplusOp; auto.
repeat rewrite Fopp_correct.
ring.
Qed.

Theorem MDekker :
 forall p q : float,
 Fbounded b p ->
 Fbounded b q ->
 (Rabs q <= Rabs p)%R -> 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 Rl1.
rewrite (Rabs_right p) in H'1; auto with real.
case (Rle_or_lt 0 q); intros Rl2.
rewrite (Rabs_right q) in H'1; auto with real.
apply Dekker1; auto.
rewrite (Faux.Rabsolu_left1 q) in H'1; auto with real.
case (Rle_or_lt p (radix * - q)); intros Rl3.
apply Dekker2; auto.
apply Dekker3; auto.
apply Rlt_le; auto.
rewrite (Faux.Rabsolu_left1 p) in H'1; auto with real.
apply MDekkerAux5; auto.
case (Rle_or_lt 0 q); intros Rl2.
rewrite (Rabs_right q) in H'1; auto with real.
case (Rle_or_lt (- p) (radix * q)); intros Rl3.
apply Dekker2; auto with float.
unfold FtoRradix in |- *; rewrite Fopp_correct; auto.
replace 0%R with (-0)%R; auto with real.
unfold FtoRradix in |- *; repeat rewrite Fopp_correct; auto.
rewrite Ropp_involutive; auto.
unfold FtoRradix in |- *; repeat rewrite Fopp_correct; auto with real.
rewrite Ropp_involutive; auto.
apply Dekker3; auto with float.
unfold FtoRradix in |- *; rewrite Fopp_correct; auto.
replace 0%R with (-0)%R; auto with real.
unfold FtoRradix in |- *; repeat rewrite Fopp_correct; auto.
rewrite Ropp_involutive; auto.
rewrite (Faux.Rabsolu_left1 q) in H'1; auto with real.
apply Dekker1; auto with float.
unfold FtoRradix in |- *; rewrite Fopp_correct; auto.
replace 0%R with (-0)%R; auto with real.
unfold FtoRradix in |- *; repeat rewrite Fopp_correct; auto.
Qed.

Theorem Dekker :
 forall p q : float,
 Fbounded b p ->
 Fbounded b q ->
 (Rabs q <= Rabs p)%R ->
 Iminus q (Iminus (Iplus p q) p) = (p + q - Iplus p q)%R :>R.
intros p q H' H'0 H'1.
apply MDekkerAux1; auto.
apply MDekker; auto.
Qed.
End Fast.