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 :
∀ p q : float,
Fbounded b p → Fbounded b q → Closest b radix (p + q) (Iplus p q).
Hypothesis IplusSym : ∀ p q : float, Iplus p q = Iplus q p.
Hypothesis
IplusOp : ∀ p q : float, Fopp (Iplus p q) = Iplus (Fopp p) (Fopp q).
Variable Iminus : float → float → float.
Hypothesis IminusPlus : ∀ 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 :
∀ 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 :
∀ p q : float,
Fbounded b p →
Fbounded b q →
∃ 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; ∃ x;
auto.
Qed.
Theorem IplusOr :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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.
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 :
∀ p q : float,
Fbounded b p → Fbounded b q → Closest b radix (p + q) (Iplus p q).
Hypothesis IplusSym : ∀ p q : float, Iplus p q = Iplus q p.
Hypothesis
IplusOp : ∀ p q : float, Fopp (Iplus p q) = Iplus (Fopp p) (Fopp q).
Variable Iminus : float → float → float.
Hypothesis IminusPlus : ∀ 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 :
∀ 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 :
∀ p q : float,
Fbounded b p →
Fbounded b q →
∃ 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; ∃ x;
auto.
Qed.
Theorem IplusOr :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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 :
∀ 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.