proposition: Lemma round_round_sqrt_aux : forall fexp1 fexp2 : Z -> Z, Valid_exp fexp1 -> Valid_exp fexp2 -> round_round_sqrt_hyp fexp1 fexp2 -> forall x, 0 < x -> (fexp2 (mag (sqrt x)) <= fexp1 (mag (sqrt x)) - 1)%Z -> generic_format beta fexp1 x -> / 2 * ulp beta fexp2 (sqrt x) < Rabs (sqrt x - midp fexp1 (sqrt x)).
proof: Proof. intros fexp1 fexp2 Vfexp1 Vfexp2 Hexp x Px Hf2 Fx. assert (Hbeta : (2 <= beta)%Z). { destruct beta as (beta_val,beta_prop). now apply Zle_bool_imp_le. } set (a := round beta fexp1 Zfloor (sqrt x)). set (u1 := bpow (fexp1 (mag (sqrt x)))). set (u2 := bpow (fexp2 (mag (sqrt x)))). set (b := / 2 * (u1 - u2)). set (b' := / 2 * (u1 + u2)). unfold midp; rewrite 2!ulp_neq_0; try now apply Rgt_not_eq, sqrt_lt_R0. apply Rnot_ge_lt; intro H; apply Rge_le in H. assert (Fa : generic_format beta fexp1 a). { unfold a. apply generic_format_round. - exact Vfexp1. - now apply valid_rnd_DN. } revert Fa; revert Fx. unfold generic_format, F2R, scaled_mantissa, cexp; simpl. set (mx := Ztrunc (x * bpow (- fexp1 (mag x)))). set (ma := Ztrunc (a * bpow (- fexp1 (mag a)))). intros Fx Fa. assert (Nna : 0 <= a). { rewrite <- (round_0 beta fexp1 Zfloor). unfold a; apply round_le. - exact Vfexp1. - now apply valid_rnd_DN. - apply sqrt_pos. } assert (Phu1 : 0 < / 2 * u1). { apply Rmult_lt_0_compat; [lra|apply bpow_gt_0]. } assert (Phu2 : 0 < / 2 * u2). { apply Rmult_lt_0_compat; [lra|apply bpow_gt_0]. } assert (Pb : 0 < b). { unfold b. rewrite <- (Rmult_0_r (/ 2)). apply Rmult_lt_compat_l; [lra|]. apply Rlt_Rminus. unfold u2, u1. apply bpow_lt. lia. } assert (Pb' : 0 < b'). { now unfold b'; rewrite Rmult_plus_distr_l; apply Rplus_lt_0_compat. } assert (Hr : sqrt x <= a + b'). { unfold b'; apply (Rplus_le_reg_r (- / 2 * u1 - a)); ring_simplify. replace (_ - _) with (sqrt x - (a + / 2 * u1)) by ring. now apply Rabs_le_inv. } assert (Hl : a + b <= sqrt x). { unfold b; apply (Rplus_le_reg_r (- / 2 * u1 - a)); ring_simplify. replace (_ + sqrt _) with (sqrt x - (a + / 2 * u1)) by ring. rewrite Ropp_mult_distr_l_reverse. now apply Rabs_le_inv in H; destruct H. } assert (Hf1 : (2 * fexp1 (mag (sqrt x)) <= fexp1 (mag (x)))%Z); [destruct (mag_sqrt_disj x Px) as [H'|H']; rewrite H'; apply Hexp|]. assert (Hlx : (fexp1 (2 * mag (sqrt x)) < 2 * mag (sqrt x))%Z). { destruct (mag_sqrt_disj x Px) as [Hlx|Hlx]. - apply (valid_exp_large fexp1 (mag x)); [|lia]. now apply mag_generic_gt; [|apply Rgt_not_eq|]. - rewrite <- Hlx. now apply mag_generic_gt; [|apply Rgt_not_eq|]. } assert (Hsl : a * a + u1 * a - u2 * a + b * b <= x). { replace (_ + _) with ((a + b) * (a + b)); [|now unfold b; field]. rewrite <- sqrt_def; [|now apply Rlt_le]. assert (H' : 0 <= a + b); [now apply Rlt_le, Rplus_le_lt_0_compat|]. now apply Rmult_le_compat. } assert (Hsr : x <= a * a + u1 * a + u2 * a + b' * b'). { replace (_ + _) with ((a + b') * (a + b')); [|now unfold b'; field]. rewrite <- (sqrt_def x); [|now apply Rlt_le]. assert (H' : 0 <= sqrt x); [now apply sqrt_pos|]. now apply Rmult_le_compat. } destruct (Req_dec a 0) as [Za|Nza]. - apply (Rlt_irrefl 0). apply Rlt_le_trans with (b * b); [now apply Rmult_lt_0_compat|]. apply Rle_trans with x. + revert Hsl; unfold Rminus; rewrite Za; do 3 rewrite Rmult_0_r. now rewrite Ropp_0; do 3 rewrite Rplus_0_l. + rewrite Fx. apply (Rmult_le_reg_r (bpow (- fexp1 (mag x)))); [now apply bpow_gt_0|]. rewrite Rmult_0_l; bpow_simplify. unfold mx. rewrite Ztrunc_floor; [|now apply Rmult_le_pos; [apply Rlt_le|apply bpow_ge_0]]. apply Req_le, IZR_eq. apply Zfloor_imp. split; [now apply Rmult_le_pos; [apply Rlt_le|apply bpow_ge_0]|simpl]. apply (Rmult_lt_reg_r (bpow (fexp1 (mag x)))); [now apply bpow_gt_0|]. rewrite Rmult_1_l; bpow_simplify. apply Rlt_le_trans with (bpow (2 * fexp1 (mag (sqrt x)))); [|now apply bpow_le]. change 2%Z with (1 + 1)%Z; rewrite Zmult_plus_distr_l; rewrite Zmult_1_l. rewrite bpow_plus. rewrite <- (sqrt_def x) at 1; [|now apply Rlt_le]. assert (sqrt x < bpow (fexp1 (mag (sqrt x)))); [|now apply Rmult_lt_compat; [apply sqrt_pos|apply sqrt_pos| |]]. apply (Rle_lt_trans _ _ _ Hr); rewrite Za; rewrite Rplus_0_l. unfold b'; change (bpow _) with u1. apply Rlt_le_trans with (/ 2 * (u1 + u1)); [|lra]. apply Rmult_lt_compat_l; [lra|]; apply Rplus_lt_compat_l. unfold u2, u1, ulp, cexp; apply bpow_lt; lia. - assert (Pa : 0 < a); [lra|]. assert (Hla : (mag a = mag (sqrt x) :> Z)). { unfold a; apply mag_DN. - exact Vfexp1. - now fold a. } assert (Hl' : 0 < - (u2 * a) + b * b). { apply (Rplus_lt_reg_r (u2 * a)); ring_simplify. unfold b; ring_simplify. apply (Rplus_lt_reg_r (/ 2 * u2 * u1)); field_simplify. replace (_ / 2) with (u2 * (a + / 2 * u1)) by field. replace (_ / 8) with (/ 4 * (u2 ^ 2 + u1 ^ 2)) by field. apply Rlt_le_trans with (u2 * bpow (mag (sqrt x))). - apply Rmult_lt_compat_l; [now unfold u2, ulp; apply bpow_gt_0|]. unfold u1; rewrite <- Hla. apply Rlt_le_trans with (a + bpow (fexp1 (mag a))). + apply Rplus_lt_compat_l. rewrite <- (Rmult_1_l (bpow _)) at 2. apply Rmult_lt_compat_r; [apply bpow_gt_0|lra]. + apply Rle_trans with (a+ ulp beta fexp1 a). right; now rewrite ulp_neq_0. apply (id_p_ulp_le_bpow _ _ _ _ Pa Fa). apply Rabs_lt_inv, bpow_mag_gt. - apply Rle_trans with (bpow (- 2) * u1 ^ 2). + unfold pow; rewrite Rmult_1_r. unfold u1, u2, ulp, cexp; bpow_simplify; apply bpow_le. now apply Hexp. + apply Rmult_le_compat. * apply bpow_ge_0. * apply pow2_ge_0. * unfold Raux.bpow, Z.pow_pos; simpl; rewrite Zmult_1_r. apply Rinv_le; [lra|]. change 4%Z with (2 * 2)%Z; apply IZR_le, Zmult_le_compat; lia. * rewrite <- (Rplus_0_l (u1 ^ 2)) at 1; apply Rplus_le_compat_r. apply pow2_ge_0. } assert (Hr' : x <= a * a + u1 * a). { rewrite Hla in Fa. rewrite <- Rmult_plus_distr_r. unfold u1, ulp, cexp. rewrite <- (Rmult_1_l (bpow _)); rewrite Fa; rewrite <- Rmult_plus_distr_r. rewrite <- Rmult_assoc; rewrite (Rmult_comm _ (IZR ma)). rewrite <- (Rmult_assoc (IZR ma)); bpow_simplify. apply (Rmult_le_reg_r (bpow (- 2 * fexp1 (mag (sqrt x))))); [now apply bpow_gt_0|bpow_simplify]. rewrite Fx at 1; bpow_simplify. rewrite <- IZR_Zpower by lia. rewrite <- plus_IZR, <- 2!mult_IZR. apply IZR_le, Zlt_succ_le, lt_IZR. unfold Z.succ; rewrite plus_IZR; do 2 rewrite mult_IZR; rewrite plus_IZR. rewrite IZR_Zpower by lia. apply (Rmult_lt_reg_r (bpow (2 * fexp1 (mag (sqrt x))))); [now apply bpow_gt_0|bpow_simplify]. rewrite <- Fx. change 2%Z with (1 + 1)%Z; rewrite Zmult_plus_distr_l; rewrite Zmult_1_l. rewrite bpow_plus; simpl. replace (_ * _) with (a * a + u1 * a + u1 * u1); [|unfold u1, ulp, cexp; rewrite Fa; ring]. apply (Rle_lt_trans _ _ _ Hsr). rewrite Rplus_assoc; apply Rplus_lt_compat_l. apply (Rplus_lt_reg_r (- b' * b' + / 2 * u1 * u2)); ring_simplify. replace (_ + _) with ((a + / 2 * u1) * u2) by ring. apply Rlt_le_trans with (bpow (mag (sqrt x)) * u2). - apply Rmult_lt_compat_r; [now unfold u2, ulp; apply bpow_gt_0|]. apply Rlt_le_trans with (a + u1); [lra|]. unfold u1; fold (cexp beta fexp1 (sqrt x)). rewrite <- cexp_DN; [|exact Vfexp1|exact Pa]; fold a. rewrite <- ulp_neq_0; trivial. apply id_p_ulp_le_bpow. + exact Pa. + now apply round_DN_pt. + apply Rle_lt_trans with (sqrt x). * now apply round_DN_pt. * apply Rabs_lt_inv. apply bpow_mag_gt. - apply Rle_trans with (/ 2 * u1 ^ 2). + apply Rle_trans with (bpow (- 2) * u1 ^ 2). * unfold pow; rewrite Rmult_1_r. unfold u2, u1, ulp, cexp. bpow_simplify. apply bpow_le. rewrite Zplus_comm. now apply Hexp. * apply Rmult_le_compat_r; [now apply pow2_ge_0|]. unfold Raux.bpow; simpl; unfold Z.pow_pos; simpl. rewrite Zmult_1_r. apply Rinv_le; [lra|]. apply IZR_le. rewrite <- (Zmult_1_l 2). apply Zmult_le_compat; lia. + assert (u2 ^ 2 < u1 ^ 2); [|unfold b'; lra]. unfold pow; do 2 rewrite Rmult_1_r. assert (H' : 0 <= u2); [unfold u2, ulp; apply bpow_ge_0|]. assert (u2 < u1); [|now apply Rmult_lt_compat]. unfold u1, u2, ulp, cexp; apply bpow_lt; lia. } apply (Rlt_irrefl (a * a + u1 * a)). apply Rlt_le_trans with (a * a + u1 * a - u2 * a + b * b). + rewrite <- (Rplus_0_r (a * a + _)) at 1. unfold Rminus; rewrite (Rplus_assoc _ _ (b * b)). now apply Rplus_lt_compat_l. + now apply Rle_trans with x. Qed.
imports: Coq(ZArith Reals Psatz) Core FTZ