diff options
Diffstat (limited to 'flocq/Prop/Div_sqrt_error.v')
| -rw-r--r-- | flocq/Prop/Div_sqrt_error.v | 872 |
1 files changed, 872 insertions, 0 deletions
diff --git a/flocq/Prop/Div_sqrt_error.v b/flocq/Prop/Div_sqrt_error.v new file mode 100644 index 00000000..76c7af95 --- /dev/null +++ b/flocq/Prop/Div_sqrt_error.v @@ -0,0 +1,872 @@ +(** +This file is part of the Flocq formalization of floating-point +arithmetic in Coq: http://flocq.gforge.inria.fr/ + +Copyright (C) 2010-2018 Sylvie Boldo +#<br /># +Copyright (C) 2010-2018 Guillaume Melquiond + +This library is free software; you can redistribute it and/or +modify it under the terms of the GNU Lesser General Public +License as published by the Free Software Foundation; either +version 3 of the License, or (at your option) any later version. + +This library is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +COPYING file for more details. +*) + +(** * Remainder of the division and square root are in the FLX format *) + +Require Import Psatz. +Require Import Core Operations Relative Sterbenz Mult_error. + +Section Fprop_divsqrt_error. + +Variable beta : radix. +Notation bpow e := (bpow beta e). + +Variable prec : Z. + +Lemma generic_format_plus_prec : + forall fexp, (forall e, (fexp e <= e - prec)%Z) -> + forall x y (fx fy: float beta), + (x = F2R fx)%R -> (y = F2R fy)%R -> (Rabs (x+y) < bpow (prec+Fexp fx))%R -> + (Rabs (x+y) < bpow (prec+Fexp fy))%R -> + generic_format beta fexp (x+y)%R. +Proof. +intros fexp Hfexp x y fx fy Hx Hy H1 H2. +case (Req_dec (x+y) 0); intros H. +rewrite H; apply generic_format_0. +rewrite Hx, Hy, <- F2R_plus. +apply generic_format_F2R. +intros _. +case_eq (Fplus fx fy). +intros mz ez Hz. +rewrite <- Hz. +apply Z.le_trans with (Z.min (Fexp fx) (Fexp fy)). +rewrite F2R_plus, <- Hx, <- Hy. +unfold cexp. +apply Z.le_trans with (1:=Hfexp _). +apply Zplus_le_reg_l with prec; ring_simplify. +apply mag_le_bpow with (1 := H). +now apply Z.min_case. +rewrite <- Fexp_Fplus, Hz. +apply Z.le_refl. +Qed. + +Context { prec_gt_0_ : Prec_gt_0 prec }. + +Notation format := (generic_format beta (FLX_exp prec)). +Notation cexp := (cexp beta (FLX_exp prec)). + +Variable choice : Z -> bool. + + +(** Remainder of the division in FLX *) +Theorem div_error_FLX : + forall rnd { Zrnd : Valid_rnd rnd } x y, + format x -> format y -> + format (x - round beta (FLX_exp prec) rnd (x/y) * y)%R. +Proof with auto with typeclass_instances. +intros rnd Zrnd x y Hx Hy. +destruct (Req_dec y 0) as [Zy|Zy]. +now rewrite Zy, Rmult_0_r, Rminus_0_r. +destruct (Req_dec (round beta (FLX_exp prec) rnd (x/y)) 0) as [Hr|Hr]. +rewrite Hr; ring_simplify (x-0*y)%R; assumption. +assert (Zx: x <> R0). +contradict Hr. +rewrite Hr. +unfold Rdiv. +now rewrite Rmult_0_l, round_0. +destruct (canonical_generic_format _ _ x Hx) as (fx,(Hx1,Hx2)). +destruct (canonical_generic_format _ _ y Hy) as (fy,(Hy1,Hy2)). +destruct (canonical_generic_format beta (FLX_exp prec) (round beta (FLX_exp prec) rnd (x / y))) as (fr,(Hr1,Hr2)). +apply generic_format_round... +unfold Rminus; apply generic_format_plus_prec with fx (Fopp (Fmult fr fy)); trivial. +intros e; apply Z.le_refl. +now rewrite F2R_opp, F2R_mult, <- Hr1, <- Hy1. +(* *) +destruct (relative_error_FLX_ex beta prec (prec_gt_0 prec) rnd (x / y)%R) as (eps,(Heps1,Heps2)). +rewrite Heps2. +rewrite <- Rabs_Ropp. +replace (-(x + - (x / y * (1 + eps) * y)))%R with (x * eps)%R by now field. +rewrite Rabs_mult. +apply Rlt_le_trans with (Rabs x * 1)%R. +apply Rmult_lt_compat_l. +now apply Rabs_pos_lt. +apply Rlt_le_trans with (1 := Heps1). +change 1%R with (bpow 0). +apply bpow_le. +generalize (prec_gt_0 prec). +clear ; omega. +rewrite Rmult_1_r. +rewrite Hx2, <- Hx1. +unfold cexp. +destruct (mag beta x) as (ex, Hex). +simpl. +specialize (Hex Zx). +apply Rlt_le. +apply Rlt_le_trans with (1 := proj2 Hex). +apply bpow_le. +unfold FLX_exp. +ring_simplify. +apply Z.le_refl. +(* *) +replace (Fexp (Fopp (Fmult fr fy))) with (Fexp fr + Fexp fy)%Z. +2: unfold Fopp, Fmult; destruct fr; destruct fy; now simpl. +replace (x + - (round beta (FLX_exp prec) rnd (x / y) * y))%R with + (y * (-(round beta (FLX_exp prec) rnd (x / y) - x/y)))%R. +2: field; assumption. +rewrite Rabs_mult. +apply Rlt_le_trans with (Rabs y * bpow (Fexp fr))%R. +apply Rmult_lt_compat_l. +now apply Rabs_pos_lt. +rewrite Rabs_Ropp. +replace (bpow (Fexp fr)) with (ulp beta (FLX_exp prec) (F2R fr)). +rewrite <- Hr1. +apply error_lt_ulp_round... +apply Rmult_integral_contrapositive_currified; try apply Rinv_neq_0_compat; assumption. +rewrite ulp_neq_0. +2: now rewrite <- Hr1. +apply f_equal. +now rewrite Hr2, <- Hr1. +replace (prec+(Fexp fr+Fexp fy))%Z with ((prec+Fexp fy)+Fexp fr)%Z by ring. +rewrite bpow_plus. +apply Rmult_le_compat_r. +apply bpow_ge_0. +rewrite Hy2, <- Hy1 ; unfold cexp, FLX_exp. +ring_simplify (prec + (mag beta y - prec))%Z. +destruct (mag beta y); simpl. +left; now apply a. +Qed. + +(** Remainder of the square in FLX (with p>1) and rounding to nearest *) +Variable Hp1 : Z.lt 1 prec. + +Theorem sqrt_error_FLX_N : + forall x, format x -> + format (x - Rsqr (round beta (FLX_exp prec) (Znearest choice) (sqrt x)))%R. +Proof with auto with typeclass_instances. +intros x Hx. +destruct (total_order_T x 0) as [[Hxz|Hxz]|Hxz]. +unfold sqrt. +destruct (Rcase_abs x). +rewrite round_0... +unfold Rsqr. +now rewrite Rmult_0_l, Rminus_0_r. +elim (Rlt_irrefl 0). +now apply Rgt_ge_trans with x. +rewrite Hxz, sqrt_0, round_0... +unfold Rsqr. +rewrite Rmult_0_l, Rminus_0_r. +apply generic_format_0. +case (Req_dec (round beta (FLX_exp prec) (Znearest choice) (sqrt x)) 0); intros Hr. +rewrite Hr; unfold Rsqr; ring_simplify (x-0*0)%R; assumption. +destruct (canonical_generic_format _ _ x Hx) as (fx,(Hx1,Hx2)). +destruct (canonical_generic_format beta (FLX_exp prec) (round beta (FLX_exp prec) (Znearest choice) (sqrt x))) as (fr,(Hr1,Hr2)). +apply generic_format_round... +unfold Rminus; apply generic_format_plus_prec with fx (Fopp (Fmult fr fr)); trivial. +intros e; apply Z.le_refl. +unfold Rsqr; now rewrite F2R_opp,F2R_mult, <- Hr1. +(* *) +apply Rle_lt_trans with x. +apply Rabs_minus_le. +apply Rle_0_sqr. +destruct (relative_error_N_FLX_ex beta prec (prec_gt_0 prec) choice (sqrt x)) as (eps,(Heps1,Heps2)). +rewrite Heps2. +rewrite Rsqr_mult, Rsqr_sqrt, Rmult_comm. 2: now apply Rlt_le. +apply Rmult_le_compat_r. +now apply Rlt_le. +apply Rle_trans with (5²/4²)%R. +rewrite <- Rsqr_div. +apply Rsqr_le_abs_1. +apply Rle_trans with (1 := Rabs_triang _ _). +rewrite Rabs_R1. +apply Rplus_le_reg_l with (-1)%R. +replace (-1 + (1 + Rabs eps))%R with (Rabs eps) by ring. +apply Rle_trans with (1 := Heps1). +rewrite Rabs_pos_eq. +apply Rmult_le_reg_l with 2%R. +now apply IZR_lt. +rewrite <- Rmult_assoc, Rinv_r, Rmult_1_l. +apply Rle_trans with (bpow (-1)). +apply bpow_le. +omega. +replace (2 * (-1 + 5 / 4))%R with (/2)%R by field. +apply Rinv_le. +now apply IZR_lt. +apply IZR_le. +unfold Zpower_pos. simpl. +rewrite Zmult_1_r. +apply Zle_bool_imp_le. +apply beta. +now apply IZR_neq. +unfold Rdiv. +apply Rmult_le_pos. +now apply IZR_le. +apply Rlt_le. +apply Rinv_0_lt_compat. +now apply IZR_lt. +now apply IZR_neq. +unfold Rsqr. +replace (5 * 5 / (4 * 4))%R with (25 * /16)%R by field. +apply Rmult_le_reg_r with 16%R. +now apply IZR_lt. +rewrite Rmult_assoc, Rinv_l, Rmult_1_r. +now apply (IZR_le _ 32). +now apply IZR_neq. +rewrite Hx2, <- Hx1; unfold cexp, FLX_exp. +ring_simplify (prec + (mag beta x - prec))%Z. +destruct (mag beta x); simpl. +rewrite <- (Rabs_right x). +apply a. +now apply Rgt_not_eq. +now apply Rgt_ge. +(* *) +replace (Fexp (Fopp (Fmult fr fr))) with (Fexp fr + Fexp fr)%Z. +2: unfold Fopp, Fmult; destruct fr; now simpl. +rewrite Hr1. +replace (x + - Rsqr (F2R fr))%R with (-((F2R fr - sqrt x)*(F2R fr + sqrt x)))%R. +2: rewrite <- (sqrt_sqrt x) at 3; auto. +2: unfold Rsqr; ring. +rewrite Rabs_Ropp, Rabs_mult. +apply Rle_lt_trans with ((/2*bpow (Fexp fr))* Rabs (F2R fr + sqrt x))%R. +apply Rmult_le_compat_r. +apply Rabs_pos. +apply Rle_trans with (/2*ulp beta (FLX_exp prec) (F2R fr))%R. +rewrite <- Hr1. +apply error_le_half_ulp_round... +right; rewrite ulp_neq_0. +2: now rewrite <- Hr1. +apply f_equal. +rewrite Hr2, <- Hr1; trivial. +rewrite Rmult_assoc, Rmult_comm. +replace (prec+(Fexp fr+Fexp fr))%Z with (Fexp fr + (prec+Fexp fr))%Z by ring. +rewrite bpow_plus, Rmult_assoc. +apply Rmult_lt_compat_l. +apply bpow_gt_0. +apply Rmult_lt_reg_l with (1 := Rlt_0_2). +apply Rle_lt_trans with (Rabs (F2R fr + sqrt x)). +right; field. +apply Rle_lt_trans with (1:=Rabs_triang _ _). +(* . *) +assert (Rabs (F2R fr) < bpow (prec + Fexp fr))%R. +rewrite Hr2. +unfold cexp, FLX_exp. +ring_simplify (prec + (mag beta (F2R fr) - prec))%Z. +destruct (mag beta (F2R fr)); simpl. +apply a. +rewrite <- Hr1; auto. +(* . *) +apply Rlt_le_trans with (bpow (prec + Fexp fr)+ Rabs (sqrt x))%R. +now apply Rplus_lt_compat_r. +(* . *) +replace (2 * bpow (prec + Fexp fr))%R with (bpow (prec + Fexp fr) + bpow (prec + Fexp fr))%R by ring. +apply Rplus_le_compat_l. +assert (sqrt x <> 0)%R. +apply Rgt_not_eq. +now apply sqrt_lt_R0. +destruct (mag beta (sqrt x)) as (es,Es). +specialize (Es H0). +apply Rle_trans with (bpow es). +now apply Rlt_le. +apply bpow_le. +case (Zle_or_lt es (prec + Fexp fr)) ; trivial. +intros H1. +absurd (Rabs (F2R fr) < bpow (es - 1))%R. +apply Rle_not_lt. +rewrite <- Hr1. +apply abs_round_ge_generic... +apply generic_format_bpow. +unfold FLX_exp; omega. +apply Es. +apply Rlt_le_trans with (1:=H). +apply bpow_le. +omega. +now apply Rlt_le. +Qed. + +Lemma sqrt_error_N_FLX_aux1 x (Fx : format x) (Px : (0 < x)%R) : + exists (mu : R) (e : Z), (format mu /\ x = mu * bpow (2 * e) :> R + /\ 1 <= mu < bpow 2)%R. +Proof. +set (e := ((mag beta x - 1) / 2)%Z). +set (mu := (x * bpow (-2 * e)%Z)%R). +assert (Hbe : (bpow (-2 * e) * bpow (2 * e) = 1)%R). +{ now rewrite <- bpow_plus; case e; simpl; [reflexivity| |]; intro p; + rewrite Z.pos_sub_diag. } +assert (Fmu : format mu); [now apply mult_bpow_exact_FLX|]. +exists mu, e; split; [exact Fmu|split; [|split]]. +{ set (e2 := (2 * e)%Z); simpl; unfold mu; rewrite Rmult_assoc. + now unfold e2; rewrite Hbe, Rmult_1_r. } +{ apply (Rmult_le_reg_r (bpow (2 * e))). + { apply bpow_gt_0. } + rewrite Rmult_1_l; set (e2 := (2 * e)%Z); simpl; unfold mu. + unfold e2; rewrite Rmult_assoc, Hbe, Rmult_1_r. + apply (Rle_trans _ (bpow (mag beta x - 1))). + { now apply bpow_le; unfold e; apply Z_mult_div_ge. } + set (l := mag _ _); rewrite <- (Rabs_pos_eq _ (Rlt_le _ _ Px)). + unfold l; apply bpow_mag_le. + intro Hx; revert Px; rewrite Hx; apply Rlt_irrefl. } +simpl; unfold mu; change (IZR _) with (bpow 2). +apply (Rmult_lt_reg_r (bpow (2 * e))); [now apply bpow_gt_0|]. +rewrite Rmult_assoc, Hbe, Rmult_1_r. +apply (Rlt_le_trans _ (bpow (mag beta x))). +{ rewrite <- (Rabs_pos_eq _ (Rlt_le _ _ Px)) at 1; apply bpow_mag_gt. } +rewrite <- bpow_plus; apply bpow_le; unfold e; set (mxm1 := (_ - 1)%Z). +replace (_ * _)%Z with (2 * (mxm1 / 2) + mxm1 mod 2 - mxm1 mod 2)%Z by ring. +rewrite <- Z.div_mod; [|now simpl]. +apply (Zplus_le_reg_r _ _ (mxm1 mod 2 - mag beta x)%Z). +unfold mxm1; destruct (Z.mod_bound_or (mag beta x - 1) 2); omega. +Qed. + +Notation u_ro := (u_ro beta prec). + +Lemma sqrt_error_N_FLX_aux2 x (Fx : format x) : + (1 <= x)%R -> + (x = 1 :> R \/ x = 1 + 2 * u_ro :> R \/ 1 + 4 * u_ro <= x)%R. +Proof. +intro HxGe1. +assert (Pu_ro : (0 <= u_ro)%R); [apply Rmult_le_pos; [lra|apply bpow_ge_0]|]. +destruct (Rle_or_lt x 1) as [HxLe1|HxGt1]; [now left; apply Rle_antisym|right]. +assert (F1 : format 1); [now apply generic_format_FLX_1|]. +assert (H2eps : (2 * u_ro = bpow (-prec + 1))%R). +{ unfold u_ro; rewrite bpow_plus; field. } +assert (HmuGe1p2eps : (1 + 2 * u_ro <= x)%R). +{ rewrite H2eps, <- succ_FLX_1. + now apply succ_le_lt; [now apply FLX_exp_valid| | |]. } +destruct (Rle_or_lt x (1 + 2 * u_ro)) as [HxLe1p2eps|HxGt1p2eps]; + [now left; apply Rle_antisym|right]. +assert (Hulp1p2eps : (ulp beta (FLX_exp prec) (1 + 2 * u_ro) = 2 * u_ro)%R). +{ destruct (ulp_succ_pos _ _ _ F1 Rlt_0_1) as [Hsucc|Hsucc]. + { now rewrite H2eps, <- succ_FLX_1, <- ulp_FLX_1. } + exfalso; revert Hsucc; apply Rlt_not_eq. + rewrite succ_FLX_1, mag_1, bpow_1, <- H2eps; simpl. + apply (Rlt_le_trans _ 2); [apply Rplus_lt_compat_l|]. + { unfold u_ro; rewrite <-Rmult_assoc, Rinv_r, Rmult_1_l; [|lra]. + change R1 with (bpow 0); apply bpow_lt; omega. } + apply IZR_le, Zle_bool_imp_le, radix_prop. } +assert (Hsucc1p2eps : + (succ beta (FLX_exp prec) (1 + 2 * u_ro) = 1 + 4 * u_ro)%R). +{ unfold succ; rewrite Rle_bool_true; [rewrite Hulp1p2eps; ring|]. + apply Rplus_le_le_0_compat; lra. } +rewrite <- Hsucc1p2eps. +apply succ_le_lt; [now apply FLX_exp_valid| |exact Fx|now simpl]. +rewrite H2eps, <- succ_FLX_1. +now apply generic_format_succ; [apply FLX_exp_valid|]. +Qed. + +Lemma sqrt_error_N_FLX_aux3 : + (u_ro / sqrt (1 + 4 * u_ro) <= 1 - 1 / sqrt (1 + 2 * u_ro))%R. +Proof. +assert (Pu_ro : (0 <= u_ro)%R); [apply Rmult_le_pos; [lra|apply bpow_ge_0]|]. +unfold Rdiv; apply (Rplus_le_reg_r (/ sqrt (1 + 2 * u_ro))); ring_simplify. +apply (Rmult_le_reg_r (sqrt (1 + 4 * u_ro) * sqrt (1 + 2 * u_ro))). +{ apply Rmult_lt_0_compat; apply sqrt_lt_R0; lra. } +field_simplify; [|split; apply Rgt_not_eq, Rlt_gt, sqrt_lt_R0; lra]. +unfold Rdiv; rewrite Rinv_1, !Rmult_1_r. +apply Rsqr_incr_0_var; [|now apply Rmult_le_pos; apply sqrt_pos]. +rewrite <-sqrt_mult; [|lra|lra]. +rewrite Rsqr_sqrt; [|apply Rmult_le_pos; lra]. +unfold Rsqr; ring_simplify; unfold pow; rewrite !Rmult_1_r. +rewrite !sqrt_def; [|lra|lra]. +apply (Rplus_le_reg_r (-u_ro * u_ro - 1 -4 * u_ro - 2 * u_ro ^ 3)). +ring_simplify; apply Rsqr_incr_0_var. +{ unfold Rsqr; ring_simplify. + unfold pow; rewrite !Rmult_1_r, !sqrt_def; [|lra|lra]. + apply (Rplus_le_reg_r (-32 * u_ro ^ 4 - 24 * u_ro ^ 3 - 4 * u_ro ^ 2)). + ring_simplify. + replace (_ + _)%R + with (((4 * u_ro ^ 2 - 28 * u_ro + 9) * u_ro + 4) * u_ro ^ 3)%R by ring. + apply Rmult_le_pos; [|now apply pow_le]. + assert (Heps_le_half : (u_ro <= 1 / 2)%R). + { unfold u_ro, Rdiv; rewrite Rmult_comm; apply Rmult_le_compat_r; [lra|]. + change 1%R with (bpow 0); apply bpow_le; omega. } + apply (Rle_trans _ (-8 * u_ro + 4)); [lra|]. + apply Rplus_le_compat_r, Rmult_le_compat_r; [apply Pu_ro|]. + now assert (H : (0 <= u_ro ^ 2)%R); [apply pow2_ge_0|lra]. } +assert (H : (u_ro ^ 3 <= u_ro ^ 2)%R). +{ unfold pow; rewrite <-!Rmult_assoc, Rmult_1_r. + apply Rmult_le_compat_l; [now apply Rmult_le_pos; apply Pu_ro|]. + now apply Rlt_le, u_ro_lt_1. } +now assert (H' : (0 <= u_ro ^ 2)%R); [apply pow2_ge_0|lra]. +Qed. + +Lemma om1ds1p2u_ro_pos : (0 <= 1 - 1 / sqrt (1 + 2 * u_ro))%R. +Proof. +unfold Rdiv; rewrite Rmult_1_l, <-Rinv_1 at 1. +apply Rle_0_minus, Rinv_le; [lra|]. +rewrite <- sqrt_1 at 1; apply sqrt_le_1_alt. +assert (H := u_ro_pos beta prec); lra. +Qed. + +Lemma om1ds1p2u_ro_le_u_rod1pu_ro : + (1 - 1 / sqrt (1 + 2 * u_ro) <= u_ro / (1 + u_ro))%R. +Proof. +assert (Pu_ro := u_ro_pos beta prec). +apply (Rmult_le_reg_r (sqrt (1 + 2 * u_ro) * (1 + u_ro))). +{ apply Rmult_lt_0_compat; [apply sqrt_lt_R0|]; lra. } +field_simplify; [|lra|intro H; apply sqrt_eq_0 in H; lra]. +unfold Rdiv, Rminus; rewrite Rinv_1, !Rmult_1_r, !Rplus_assoc. +rewrite <-(Rplus_0_r (sqrt _ * _)) at 2; apply Rplus_le_compat_l. +apply (Rplus_le_reg_r (1 + u_ro)); ring_simplify. +rewrite <-(sqrt_square (_ + 1)); [|lra]; apply sqrt_le_1_alt. +assert (H : (0 <= u_ro * u_ro)%R); [apply Rmult_le_pos|]; lra. +Qed. + +Lemma s1p2u_rom1_pos : (0 <= sqrt (1 + 2 * u_ro) - 1)%R. +apply (Rplus_le_reg_r 1); ring_simplify. +rewrite <-sqrt_1 at 1; apply sqrt_le_1_alt. +assert (H := u_ro_pos beta prec); lra. +Qed. + +Theorem sqrt_error_N_FLX x (Fx : format x) : + (Rabs (round beta (FLX_exp prec) (Znearest choice) (sqrt x) - sqrt x) + <= (1 - 1 / sqrt (1 + 2 * u_ro)) * Rabs (sqrt x))%R. +Proof. +assert (Peps := u_ro_pos beta prec). +assert (Peps' : (0 < u_ro)%R). +{ unfold u_ro; apply Rmult_lt_0_compat; [lra|apply bpow_gt_0]. } +assert (Pb := om1ds1p2u_ro_pos). +assert (Pb' := s1p2u_rom1_pos). +destruct (Rle_or_lt x 0) as [Nx|Px]. +{ rewrite (sqrt_neg _ Nx), round_0, Rabs_R0, Rmult_0_r; [|apply valid_rnd_N]. + now unfold Rminus; rewrite Rplus_0_l, Rabs_Ropp, Rabs_R0; right. } +destruct (sqrt_error_N_FLX_aux1 _ Fx Px) + as (mu, (e, (Fmu, (Hmu, (HmuGe1, HmuLtsqradix))))). +pose (t := sqrt x). +set (rt := round _ _ _ _). +assert (Ht : (t = sqrt mu * bpow e)%R). +{ unfold t; rewrite Hmu, sqrt_mult_alt; [|now apply (Rle_trans _ _ _ Rle_0_1)]. + now rewrite sqrt_bpow. } +destruct (sqrt_error_N_FLX_aux2 _ Fmu HmuGe1) as [Hmu'|[Hmu'|Hmu']]. +{ unfold rt; fold t; rewrite Ht, Hmu', sqrt_1, Rmult_1_l. + rewrite round_generic; [|now apply valid_rnd_N|]. + { rewrite Rminus_diag_eq, Rabs_R0; [|now simpl]. + now apply Rmult_le_pos; [|apply Rabs_pos]. } + apply generic_format_bpow'; [now apply FLX_exp_valid|]. + unfold FLX_exp; omega. } +{ assert (Hsqrtmu : (1 <= sqrt mu < 1 + u_ro)%R); [rewrite Hmu'; split|]. + { rewrite <- sqrt_1 at 1; apply sqrt_le_1_alt; lra. } + { rewrite <- sqrt_square; [|lra]; apply sqrt_lt_1_alt; split; [lra|]. + ring_simplify; assert (0 < u_ro ^ 2)%R; [apply pow_lt|]; lra. } + assert (Fbpowe : generic_format beta (FLX_exp prec) (bpow e)). + { apply generic_format_bpow; unfold FLX_exp; omega. } + assert (Hrt : rt = bpow e :> R). + { unfold rt; fold t; rewrite Ht; simpl; apply Rle_antisym. + { apply round_N_le_midp; [now apply FLX_exp_valid|exact Fbpowe|]. + apply (Rlt_le_trans _ ((1 + u_ro) * bpow e)). + { now apply Rmult_lt_compat_r; [apply bpow_gt_0|]. } + unfold succ; rewrite Rle_bool_true; [|now apply bpow_ge_0]. + rewrite ulp_bpow; unfold FLX_exp. + unfold Z.sub, u_ro; rewrite !bpow_plus; right; field. } + apply round_ge_generic; + [now apply FLX_exp_valid|now apply valid_rnd_N|exact Fbpowe|]. + rewrite <- (Rmult_1_l (bpow _)) at 1. + now apply Rmult_le_compat_r; [apply bpow_ge_0|]. } + fold t; rewrite Hrt, Ht, Hmu', <-(Rabs_pos_eq _ Pb), <-Rabs_mult. + rewrite Rabs_minus_sym; right; f_equal; field; lra. } +assert (Hsqrtmu : (1 + u_ro < sqrt mu)%R). +{ apply (Rlt_le_trans _ (sqrt (1 + 4 * u_ro))); [|now apply sqrt_le_1_alt]. + assert (P1peps : (0 <= 1 + u_ro)%R) + by now apply Rplus_le_le_0_compat; [lra|apply Peps]. + rewrite <- (sqrt_square (1 + u_ro)); [|lra]. + apply sqrt_lt_1_alt; split; [now apply Rmult_le_pos|]. + apply (Rplus_lt_reg_r (-1 - 2 * u_ro)); ring_simplify; simpl. + rewrite Rmult_1_r; apply Rmult_lt_compat_r; [apply Peps'|]. + now apply (Rlt_le_trans _ 1); [apply u_ro_lt_1|lra]. } +assert (Hulpt : (ulp beta (FLX_exp prec) t = 2 * u_ro * bpow e)%R). +{ unfold ulp; rewrite Req_bool_false; [|apply Rgt_not_eq, Rlt_gt]. + { unfold u_ro; rewrite <-Rmult_assoc, Rinv_r, Rmult_1_l, <-bpow_plus; [|lra]. + f_equal; unfold cexp, FLX_exp. + assert (Hmagt : (mag beta t = 1 + e :> Z)%Z). + { apply mag_unique. + unfold t; rewrite (Rabs_pos_eq _ (Rlt_le _ _ (sqrt_lt_R0 _ Px))). + fold t; split. + { rewrite Ht; replace (_ - _)%Z with e by ring. + rewrite <- (Rmult_1_l (bpow _)) at 1; apply Rmult_le_compat_r. + { apply bpow_ge_0. } + now rewrite <- sqrt_1; apply sqrt_le_1_alt. } + rewrite bpow_plus, bpow_1, Ht; simpl. + apply Rmult_lt_compat_r; [now apply bpow_gt_0|]. + rewrite <- sqrt_square. + { apply sqrt_lt_1_alt; split; [lra|]. + apply (Rlt_le_trans _ _ _ HmuLtsqradix); right. + now unfold bpow, Z.pow_pos; simpl; rewrite Zmult_1_r, mult_IZR. } + apply IZR_le, (Z.le_trans _ 2), Zle_bool_imp_le, radix_prop; omega. } + rewrite Hmagt; ring. } + rewrite Ht; apply Rmult_lt_0_compat; [|now apply bpow_gt_0]. + now apply (Rlt_le_trans _ 1); [lra|rewrite <- sqrt_1; apply sqrt_le_1_alt]. } +assert (Pt : (0 < t)%R). +{ rewrite Ht; apply Rmult_lt_0_compat; [lra|apply bpow_gt_0]. } +assert (H : (Rabs ((rt - sqrt x) / sqrt x) + <= 1 - 1 / sqrt (1 + 2 * u_ro))%R). +{ unfold Rdiv; rewrite Rabs_mult, (Rabs_pos_eq (/ _)); + [|now left; apply Rinv_0_lt_compat]. + apply (Rle_trans _ ((u_ro * bpow e) / t)). + { unfold Rdiv; apply Rmult_le_compat_r; [now left; apply Rinv_0_lt_compat|]. + apply (Rle_trans _ _ _ (error_le_half_ulp _ _ _ _)). + fold t; rewrite Hulpt; right; field. } + apply (Rle_trans _ (u_ro / sqrt (1 + 4 * u_ro))). + { apply (Rle_trans _ (u_ro * bpow e / (sqrt (1 + 4 * u_ro) * bpow e))). + { unfold Rdiv; apply Rmult_le_compat_l; + [now apply Rmult_le_pos; [apply Peps|apply bpow_ge_0]|]. + apply Rinv_le. + { apply Rmult_lt_0_compat; [apply sqrt_lt_R0; lra|apply bpow_gt_0]. } + now rewrite Ht; apply Rmult_le_compat_r; + [apply bpow_ge_0|apply sqrt_le_1_alt]. } + right; field; split; apply Rgt_not_eq, Rlt_gt; + [apply sqrt_lt_R0; lra|apply bpow_gt_0]. } + apply sqrt_error_N_FLX_aux3. } +revert H; unfold Rdiv; rewrite Rabs_mult, Rabs_Rinv; [|fold t; lra]; intro H. +apply (Rmult_le_reg_r (/ Rabs (sqrt x))); + [apply Rinv_0_lt_compat, Rabs_pos_lt; fold t; lra|]. +apply (Rle_trans _ _ _ H); right; field; split; [apply Rabs_no_R0;fold t|]; lra. +Qed. + +Theorem sqrt_error_N_FLX_ex x (Fx : format x) : + exists eps, + (Rabs eps <= 1 - 1 / sqrt (1 + 2 * u_ro))%R /\ + round beta (FLX_exp prec) (Znearest choice) (sqrt x) + = (sqrt x * (1 + eps))%R. +Proof. +now apply relative_error_le_conversion; + [apply valid_rnd_N|apply om1ds1p2u_ro_pos|apply sqrt_error_N_FLX]. +Qed. + +Lemma sqrt_error_N_round_ex_derive : + forall x rx, + (exists eps, + (Rabs eps <= 1 - 1 / sqrt (1 + 2 * u_ro))%R /\ rx = (x * (1 + eps))%R) -> + exists eps, + (Rabs eps <= sqrt (1 + 2 * u_ro) - 1)%R /\ x = (rx * (1 + eps))%R. +Proof. +intros x rx (d, (Bd, Hd)). +assert (H := Rabs_le_inv _ _ Bd). +assert (H' := om1ds1p2u_ro_le_u_rod1pu_ro). +assert (H'' := u_rod1pu_ro_le_u_ro beta prec). +assert (H''' := u_ro_lt_1 beta prec prec_gt_0_). +assert (Hpos := s1p2u_rom1_pos). +destruct (Req_dec rx 0) as [Zfx|Nzfx]. +{ exists 0%R; split; [now rewrite Rabs_R0|]. + rewrite Rplus_0_r, Rmult_1_r, Zfx; rewrite Zfx in Hd. + destruct (Rmult_integral _ _ (sym_eq Hd)); lra. } +destruct (Req_dec x 0) as [Zx|Nzx]. +{ now exfalso; revert Hd; rewrite Zx; rewrite Rmult_0_l. } +set (d' := ((x - rx) / rx)%R). +assert (Hd' : (Rabs d' <= sqrt (1 + 2 * u_ro) - 1)%R). +{ unfold d'; rewrite Hd. + replace (_ / _)%R with (- d / (1 + d))%R; [|now field; split; lra]. + unfold Rdiv; rewrite Rabs_mult, Rabs_Ropp. + rewrite (Rabs_pos_eq (/ _)); [|apply Rlt_le, Rinv_0_lt_compat; lra]. + apply (Rmult_le_reg_r (1 + d)); [lra|]. + rewrite Rmult_assoc, Rinv_l, Rmult_1_r; [|lra]. + apply (Rle_trans _ _ _ Bd). + apply (Rle_trans _ ((sqrt (1 + 2 * u_ro) - 1) * (1/sqrt (1 + 2 * u_ro)))); + [right; field|apply Rmult_le_compat_l]; lra. } +now exists d'; split; [exact Hd'|]; unfold d'; field. +Qed. + +(** sqrt(1 + 2 u_ro) - 1 <= u_ro *) +Theorem sqrt_error_N_FLX_round_ex : + forall x, + format x -> + exists eps, + (Rabs eps <= sqrt (1 + 2 * u_ro) - 1)%R /\ + sqrt x = (round beta (FLX_exp prec) (Znearest choice) (sqrt x) * (1 + eps))%R. +Proof. +now intros x Fx; apply sqrt_error_N_round_ex_derive, sqrt_error_N_FLX_ex. +Qed. + +Variable emin : Z. +Hypothesis Hemin : (emin <= 2 * (1 - prec))%Z. + +Theorem sqrt_error_N_FLT_ex : + forall x, + generic_format beta (FLT_exp emin prec) x -> + exists eps, + (Rabs eps <= 1 - 1 / sqrt (1 + 2 * u_ro))%R /\ + round beta (FLT_exp emin prec) (Znearest choice) (sqrt x) + = (sqrt x * (1 + eps))%R. +Proof. +intros x Fx. +assert (Heps := u_ro_pos). +assert (Pb := om1ds1p2u_ro_pos). +destruct (Rle_or_lt x 0) as [Nx|Px]. +{ exists 0%R; split; [now rewrite Rabs_R0|]. + now rewrite (sqrt_neg x Nx), round_0, Rmult_0_l; [|apply valid_rnd_N]. } +assert (Fx' := generic_format_FLX_FLT _ _ _ _ Fx). +destruct (sqrt_error_N_FLX_ex _ Fx') as (d, (Bd, Hd)). +exists d; split; [exact Bd|]; rewrite <-Hd; apply round_FLT_FLX. +apply (Rle_trans _ (bpow (emin / 2)%Z)). +{ apply bpow_le, Z.div_le_lower_bound; lia. } +apply (Rle_trans _ _ _ (sqrt_bpow_ge _ _)). +rewrite Rabs_pos_eq; [|now apply sqrt_pos]; apply sqrt_le_1_alt. +revert Fx; apply generic_format_ge_bpow; [|exact Px]. +intro e; unfold FLT_exp; apply Z.le_max_r. +Qed. + +(** sqrt(1 + 2 u_ro) - 1 <= u_ro *) +Theorem sqrt_error_N_FLT_round_ex : + forall x, + generic_format beta (FLT_exp emin prec) x -> + exists eps, + (Rabs eps <= sqrt (1 + 2 * u_ro) - 1)%R /\ + sqrt x + = (round beta (FLT_exp emin prec) (Znearest choice) (sqrt x) * (1 + eps))%R. +Proof. +now intros x Fx; apply sqrt_error_N_round_ex_derive, sqrt_error_N_FLT_ex. +Qed. + +End Fprop_divsqrt_error. + +Section format_REM_aux. + +Variable beta : radix. +Notation bpow e := (bpow beta e). + +Variable fexp : Z -> Z. +Context { valid_exp : Valid_exp fexp }. +Context { monotone_exp : Monotone_exp fexp }. + +Variable rnd : R -> Z. +Context { valid_rnd : Valid_rnd rnd }. + +Notation format := (generic_format beta fexp). + +Lemma format_REM_aux: + forall x y : R, + format x -> format y -> (0 <= x)%R -> (0 < y)%R -> + ((0 < x/y < /2)%R -> rnd (x/y) = 0%Z) -> + format (x - IZR (rnd (x/y))*y). +Proof with auto with typeclass_instances. +intros x y Fx Fy Hx Hy rnd_small. +pose (n:=rnd (x / y)). +assert (Hn:(IZR n = round beta (FIX_exp 0) rnd (x/y))%R). +unfold round, FIX_exp, cexp, scaled_mantissa, F2R; simpl. +now rewrite 2!Rmult_1_r. +assert (H:(0 <= n)%Z). +apply le_IZR; rewrite Hn; simpl. +apply Rle_trans with (round beta (FIX_exp 0) rnd 0). +right; apply sym_eq, round_0... +apply round_le... +apply Fourier_util.Rle_mult_inv_pos; assumption. +case (Zle_lt_or_eq 0 n); try exact H. +clear H; intros H. +case (Zle_lt_or_eq 1 n). +omega. +clear H; intros H. +set (ex := cexp beta fexp x). +set (ey := cexp beta fexp y). +set (mx := Ztrunc (scaled_mantissa beta fexp x)). +set (my := Ztrunc (scaled_mantissa beta fexp y)). +case (Zle_or_lt ey ex); intros Hexy. +(* ey <= ex *) +assert (H0:(x-IZR n *y = F2R (Float beta (mx*beta^(ex-ey) - n*my) ey))%R). +unfold Rminus; rewrite Rplus_comm. +replace (IZR n) with (F2R (Float beta n 0)). +rewrite Fx, Fy. +fold mx my ex ey. +rewrite <- F2R_mult. +rewrite <- F2R_opp. +rewrite <- F2R_plus. +unfold Fplus. simpl. +rewrite Zle_imp_le_bool with (1 := Hexy). +f_equal; f_equal; ring. +unfold F2R; simpl; ring. +fold n; rewrite H0. +apply generic_format_F2R. +rewrite <- H0; intros H3. +apply monotone_exp. +apply mag_le_abs. +rewrite H0; apply F2R_neq_0; easy. +apply Rmult_le_reg_l with (/Rabs y)%R. +apply Rinv_0_lt_compat. +apply Rabs_pos_lt. +now apply Rgt_not_eq. +rewrite Rinv_l. +2: apply Rgt_not_eq, Rabs_pos_lt. +2: now apply Rgt_not_eq. +rewrite <- Rabs_Rinv. +2: now apply Rgt_not_eq. +rewrite <- Rabs_mult. +replace (/y * (x - IZR n *y))%R with (-(IZR n - x/y))%R. +rewrite Rabs_Ropp. +rewrite Hn. +apply Rle_trans with (1:= error_le_ulp beta (FIX_exp 0) _ _). +rewrite ulp_FIX. +simpl; apply Rle_refl. +field. +now apply Rgt_not_eq. +(* ex < ey: impossible as 1 < n *) +absurd (1 < n)%Z; try easy. +apply Zle_not_lt. +apply le_IZR; simpl; rewrite Hn. +apply round_le_generic... +apply generic_format_FIX. +exists (Float beta 1 0); try easy. +unfold F2R; simpl; ring. +apply Rmult_le_reg_r with y; try easy. +unfold Rdiv; rewrite Rmult_assoc. +rewrite Rinv_l, Rmult_1_r, Rmult_1_l. +2: now apply Rgt_not_eq. +assert (mag beta x < mag beta y)%Z. +case (Zle_or_lt (mag beta y) (mag beta x)); try easy. +intros J; apply monotone_exp in J; clear -J Hexy. +unfold ex, ey, cexp in Hexy; omega. +left; apply lt_mag with beta; easy. +(* n = 1 -> Sterbenz + rnd_small *) +intros Hn'; fold n; rewrite <- Hn'. +rewrite Rmult_1_l. +case Hx; intros Hx'. +assert (J:(0 < x/y)%R). +apply Fourier_util.Rlt_mult_inv_pos; assumption. +apply sterbenz... +assert (H0:(Rabs (1 - x/y) < 1)%R). +rewrite Hn', Hn. +apply Rlt_le_trans with (ulp beta (FIX_exp 0) (round beta (FIX_exp 0) rnd (x / y)))%R. +apply error_lt_ulp_round... +now apply Rgt_not_eq. +rewrite ulp_FIX. +rewrite <- Hn, <- Hn'. +apply Rle_refl. +apply Rabs_lt_inv in H0. +split; apply Rmult_le_reg_l with (/y)%R; try now apply Rinv_0_lt_compat. +unfold Rdiv; rewrite <- Rmult_assoc. +rewrite Rinv_l. +2: now apply Rgt_not_eq. +rewrite Rmult_1_l, Rmult_comm; fold (x/y)%R. +case (Rle_or_lt (/2) (x/y)); try easy. +intros K. +elim Zlt_not_le with (1 := H). +apply Zeq_le. +apply rnd_small. +now split. +apply Ropp_le_cancel; apply Rplus_le_reg_l with 1%R. +apply Rle_trans with (1-x/y)%R. +2: right; unfold Rdiv; ring. +left; apply Rle_lt_trans with (2:=proj1 H0). +right; field. +now apply Rgt_not_eq. +rewrite <- Hx', Rminus_0_l. +now apply generic_format_opp. +(* n = 0 *) +clear H; intros H; fold n; rewrite <- H. +now rewrite Rmult_0_l, Rminus_0_r. +Qed. + +End format_REM_aux. + +Section format_REM. + +Variable beta : radix. +Notation bpow e := (bpow beta e). + +Variable fexp : Z -> Z. +Context { valid_exp : Valid_exp fexp }. +Context { monotone_exp : Monotone_exp fexp }. + +Notation format := (generic_format beta fexp). + +Theorem format_REM : + forall rnd : R -> Z, Valid_rnd rnd -> + forall x y : R, + ((Rabs (x/y) < /2)%R -> rnd (x/y)%R = 0%Z) -> + format x -> format y -> + format (x - IZR (rnd (x/y)%R) * y). +Proof with auto with typeclass_instances. +(* assume 0 < y *) +assert (H: forall rnd : R -> Z, Valid_rnd rnd -> + forall x y : R, + ((Rabs (x/y) < /2)%R -> rnd (x/y)%R = 0%Z) -> + format x -> format y -> (0 < y)%R -> + format (x - IZR (rnd (x/y)%R) * y)). +intros rnd valid_rnd x y Hrnd Fx Fy Hy. +case (Rle_or_lt 0 x); intros Hx. +apply format_REM_aux; try easy. +intros K. +apply Hrnd. +rewrite Rabs_pos_eq. +apply K. +apply Rlt_le, K. +replace (x - IZR (rnd (x/y)) * y)%R with + (- (-x - IZR (Zrnd_opp rnd (-x/y)) * y))%R. +apply generic_format_opp. +apply format_REM_aux; try easy... +now apply generic_format_opp. +apply Ropp_le_cancel; rewrite Ropp_0, Ropp_involutive; now left. +replace (- x / y)%R with (- (x/y))%R by (unfold Rdiv; ring). +intros K. +unfold Zrnd_opp. +rewrite Ropp_involutive, Hrnd. +easy. +rewrite Rabs_left. +apply K. +apply Ropp_lt_cancel. +now rewrite Ropp_0. +unfold Zrnd_opp. +replace (- (- x / y))%R with (x / y)%R by (unfold Rdiv; ring). +rewrite opp_IZR. +ring. +(* *) +intros rnd valid_rnd x y Hrnd Fx Fy. +case (Rle_or_lt 0 y); intros Hy. +destruct Hy as [Hy|Hy]. +now apply H. +now rewrite <- Hy, Rmult_0_r, Rminus_0_r. +replace (IZR (rnd (x/y)) * y)%R with + (IZR ((Zrnd_opp rnd) ((x / -y))) * -y)%R. +apply H; try easy... +replace (x / - y)%R with (- (x/y))%R. +intros K. +unfold Zrnd_opp. +rewrite Ropp_involutive, Hrnd. +easy. +now rewrite <- Rabs_Ropp. +field; now apply Rlt_not_eq. +now apply generic_format_opp. +apply Ropp_lt_cancel; now rewrite Ropp_0, Ropp_involutive. +unfold Zrnd_opp. +replace (- (x / - y))%R with (x/y)%R. +rewrite opp_IZR. +ring. +field; now apply Rlt_not_eq. +Qed. + +Theorem format_REM_ZR: + forall x y : R, + format x -> format y -> + format (x - IZR (Ztrunc (x/y)) * y). +Proof with auto with typeclass_instances. +intros x y Fx Fy. +apply format_REM; try easy... +intros K. +apply Z.abs_0_iff. +rewrite <- Ztrunc_abs. +rewrite Ztrunc_floor by apply Rabs_pos. +apply Zle_antisym. +replace 0%Z with (Zfloor (/2)). +apply Zfloor_le. +now apply Rlt_le. +apply Zfloor_imp. +simpl ; lra. +apply Zfloor_lub. +apply Rabs_pos. +Qed. + +Theorem format_REM_N : + forall choice, + forall x y : R, + format x -> format y -> + format (x - IZR (Znearest choice (x/y)) * y). +Proof with auto with typeclass_instances. +intros choice x y Fx Fy. +apply format_REM; try easy... +intros K. +apply Znearest_imp. +now rewrite Rminus_0_r. +Qed. + +End format_REM. |