From 0af966a42eb60e9af43f9a450d924758a83946c6 Mon Sep 17 00:00:00 2001 From: Guillaume Melquiond Date: Tue, 22 Sep 2015 15:41:50 +0200 Subject: Upgrade to Flocq 2.5.0. Note: this version of Flocq is compatible with both Coq 8.4 and 8.5. --- flocq/Appli/Fappli_IEEE.v | 38 +++++--- flocq/Appli/Fappli_IEEE_bits.v | 4 +- flocq/Appli/Fappli_double_round.v | 179 ++++++++++++++++++++++---------------- flocq/Appli/Fappli_rnd_odd.v | 98 +++++++++++++++++---- 4 files changed, 211 insertions(+), 108 deletions(-) (limited to 'flocq/Appli') diff --git a/flocq/Appli/Fappli_IEEE.v b/flocq/Appli/Fappli_IEEE.v index 9b5826c1..23999a50 100644 --- a/flocq/Appli/Fappli_IEEE.v +++ b/flocq/Appli/Fappli_IEEE.v @@ -48,9 +48,9 @@ Section Binary. Implicit Arguments exist [[A] [P]]. -(** prec is the number of bits of the mantissa including the implicit one - emax is the exponent of the infinities - Typically p=24 and emax = 128 in single precision *) +(** [prec] is the number of bits of the mantissa including the implicit one; + [emax] is the exponent of the infinities. + For instance, binary32 is defined by [prec = 24] and [emax = 128]. *) Variable prec emax : Z. Context (prec_gt_0_ : Prec_gt_0 prec). Hypothesis Hmax : (prec < emax)%Z. @@ -74,8 +74,7 @@ Definition valid_binary x := end. (** Basic type used for representing binary FP numbers. - Note that there is exactly one such object per FP datum. - NaNs do not have any payload. They cannot be distinguished. *) + Note that there is exactly one such object per FP datum. *) Definition nan_pl := {pl | (Zpos (digits2_pos pl) @@ -647,7 +648,8 @@ generalize (prec_gt_0 prec). clear -Hmax ; omega. Qed. -(** mantissa, round and sticky bits *) +(** Truncation *) + Record shr_record := { shr_m : Z ; shr_r : bool ; shr_s : bool }. Definition shr_1 mrs := @@ -695,7 +697,7 @@ Qed. Definition shr mrs e n := match n with - | Zpos p => (iter_pos p _ shr_1 mrs, (e + n)%Z) + | Zpos p => (iter_pos shr_1 p mrs, (e + n)%Z) | _ => (mrs, e) end. @@ -746,24 +748,24 @@ destruct n as [|n|n]. now destruct l as [|[| |]]. 2: now destruct l as [|[| |]]. unfold shr. -rewrite iter_nat_of_P. +rewrite iter_pos_nat. rewrite Zpos_eq_Z_of_nat_o_nat_of_P. induction (nat_of_P n). simpl. rewrite Zplus_0_r. now destruct l as [|[| |]]. -simpl nat_rect. +rewrite iter_nat_S. rewrite inj_S. unfold Zsucc. -rewrite Zplus_assoc. +rewrite Zplus_assoc. revert IHn0. apply inbetween_shr_1. clear -Hm. induction n0. now destruct l as [|[| |]]. -simpl. +rewrite iter_nat_S. revert IHn0. -generalize (iter_nat n0 shr_record shr_1 (shr_record_of_loc m l)). +generalize (iter_nat shr_1 n0 (shr_record_of_loc m l)). clear. intros (m, r, s) Hm. now destruct m as [|[m|m|]|m] ; try (now elim Hm) ; destruct r as [|] ; destruct s as [|]. @@ -827,6 +829,8 @@ intros H. now inversion H. Qed. +(** Rounding modes *) + Inductive mode := mode_NE | mode_ZR | mode_DN | mode_UP | mode_NA. Definition round_mode m := @@ -927,7 +931,6 @@ destruct (truncate radix2 fexp (Z0, e1, loc_Exact)) as ((m2, e2), l2). rewrite shr_m_shr_record_of_loc. intros Hm2. rewrite Hm2. -intros z. repeat split. rewrite Rlt_bool_true. repeat split. @@ -1178,6 +1181,8 @@ destruct x as [sx|sx|sx [plx Hplx]|sx mx ex Hx], y as [sy|sy|sy [ply Hply]|sy my apply B2FF_FF2B. Qed. +(** Normalization and rounding *) + Definition shl_align mx ex ex' := match (ex' - ex)%Z with | Zneg d => (shift_pos d mx, ex') @@ -1361,6 +1366,7 @@ now apply F2R_lt_0_compat. Qed. (** Addition *) + Definition Bplus plus_nan m x y := let f pl := B754_nan (fst pl) (snd pl) in match x, y with @@ -1475,7 +1481,7 @@ elim Rle_not_lt with (1 := Bz). generalize (bounded_lt_emax _ _ Hx) (bounded_lt_emax _ _ Hy) (andb_prop _ _ Hx) (andb_prop _ _ Hy). intros Bx By (Hx',_) (Hy',_). generalize (canonic_canonic_mantissa sx _ _ Hx') (canonic_canonic_mantissa sy _ _ Hy'). -clear -Bx By Hs. +clear -Bx By Hs prec_gt_0_. intros Cx Cy. destruct sx. (* ... *) @@ -1542,6 +1548,8 @@ now apply f_equal. apply Sz. Qed. +(** Subtraction *) + Definition Bminus minus_nan m x y := Bplus minus_nan m x (Bopp pair y). Theorem Bminus_correct : @@ -1571,6 +1579,7 @@ rewrite is_finite_Bopp. auto. now destruct y as [ | | | ]. Qed. (** Division *) + Definition Fdiv_core_binary m1 e1 m2 e2 := let d1 := Zdigits2 m1 in let d2 := Zdigits2 m2 in @@ -1737,6 +1746,7 @@ now rewrite B2FF_FF2B. Qed. (** Square root *) + Definition Fsqrt_core_binary m e := let d := Zdigits2 m in let s := Zmax (2 * prec - d) 0 in diff --git a/flocq/Appli/Fappli_IEEE_bits.v b/flocq/Appli/Fappli_IEEE_bits.v index 5a77bf57..87aa1046 100644 --- a/flocq/Appli/Fappli_IEEE_bits.v +++ b/flocq/Appli/Fappli_IEEE_bits.v @@ -617,7 +617,7 @@ apply refl_equal. Qed. Definition default_nan_pl32 : bool * nan_pl 24 := - (false, exist _ (iter_nat 22 _ xO xH) (refl_equal true)). + (false, exist _ (iter_nat xO 22 xH) (refl_equal true)). Definition unop_nan_pl32 (f : binary32) : bool * nan_pl 24 := match f with @@ -660,7 +660,7 @@ apply refl_equal. Qed. Definition default_nan_pl64 : bool * nan_pl 53 := - (false, exist _ (iter_nat 51 _ xO xH) (refl_equal true)). + (false, exist _ (iter_nat xO 51 xH) (refl_equal true)). Definition unop_nan_pl64 (f : binary64) : bool * nan_pl 53 := match f with diff --git a/flocq/Appli/Fappli_double_round.v b/flocq/Appli/Fappli_double_round.v index f83abc47..3ff6c31a 100644 --- a/flocq/Appli/Fappli_double_round.v +++ b/flocq/Appli/Fappli_double_round.v @@ -72,12 +72,15 @@ assert (Hx2 : x - round beta fexp1 Zfloor x < / 2 * (ulp beta fexp1 x - ulp beta fexp2 x)). { now apply (Rplus_lt_reg_r (round beta fexp1 Zfloor x)); ring_simplify. } set (x'' := round beta fexp2 (Znearest choice2) x). -assert (Hr1 : Rabs (x'' - x) <= / 2 * bpow (fexp2 (ln_beta x))); - [now unfold x''; apply ulp_half_error|]. +assert (Hr1 : Rabs (x'' - x) <= / 2 * bpow (fexp2 (ln_beta x))). +apply Rle_trans with (/ 2 * ulp beta fexp2 x). +now unfold x''; apply error_le_half_ulp... +rewrite ulp_neq_0;[now right|now apply Rgt_not_eq]. assert (Pxx' : 0 <= x - x'). { apply Rle_0_minus. apply round_DN_pt. exact Vfexp1. } +rewrite 2!ulp_neq_0 in Hx2; try (apply Rgt_not_eq; assumption). assert (Hr2 : Rabs (x'' - x') < / 2 * bpow (fexp1 (ln_beta x))). { replace (x'' - x') with (x'' - x + (x - x')) by ring. apply (Rle_lt_trans _ _ _ (Rabs_triang _ _)). @@ -114,6 +117,7 @@ destruct (Req_dec x'' 0) as [Zx''|Nzx'']. + (bpow (ln_beta x) - / 2 * bpow (fexp2 (ln_beta x)))) by ring. apply Rplus_le_lt_compat; [exact Hr1|]. + rewrite ulp_neq_0 in Hx1;[idtac| now apply Rgt_not_eq]. now rewrite Rabs_right; [|apply Rle_ge; apply Rlt_le]. - unfold x'' in Nzx'' |- *. now apply ln_beta_round_ge; [|apply valid_rnd_N|]. } @@ -168,12 +172,14 @@ assert (Pxx' : 0 <= x - x'). apply round_DN_pt. exact Vfexp1. } assert (x < bpow (ln_beta x) - / 2 * bpow (fexp2 (ln_beta x))); - [|now apply double_round_lt_mid_further_place']. + [|apply double_round_lt_mid_further_place'; try assumption]... +2: rewrite ulp_neq_0;[assumption|now apply Rgt_not_eq]. destruct (Req_dec x' 0) as [Zx'|Nzx']. - (* x' = 0 *) rewrite Zx' in Hx2; rewrite Rminus_0_r in Hx2. apply (Rlt_le_trans _ _ _ Hx2). rewrite Rmult_minus_distr_l. + rewrite 2!ulp_neq_0;[idtac|now apply Rgt_not_eq|now apply Rgt_not_eq]. apply Rplus_le_compat_r. apply (Rmult_le_reg_r (bpow (- ln_beta x))); [now apply bpow_gt_0|]. unfold ulp, canonic_exp; bpow_simplify. @@ -199,7 +205,7 @@ destruct (Req_dec x' 0) as [Zx'|Nzx']. { apply (Rplus_le_reg_r (ulp beta fexp1 x)); ring_simplify. rewrite <- ulp_DN. - change (round _ _ _ _) with x'. - apply succ_le_bpow. + apply id_p_ulp_le_bpow. + exact Px'. + change x' with (round beta fexp1 Zfloor x). now apply generic_format_round; [|apply valid_rnd_DN]. @@ -210,10 +216,14 @@ destruct (Req_dec x' 0) as [Zx'|Nzx']. - exact Vfexp1. - exact Px'. } fold (canonic_exp beta fexp2 x); fold (ulp beta fexp2 x). - assert (/ 2 * ulp beta fexp1 x <= ulp beta fexp1 x); [|lra]. + assert (/ 2 * ulp beta fexp1 x <= ulp beta fexp1 x). rewrite <- (Rmult_1_l (ulp _ _ _)) at 2. apply Rmult_le_compat_r; [|lra]. - apply bpow_ge_0. + apply ulp_ge_0. + rewrite 2!ulp_neq_0 in Hx2;[|now apply Rgt_not_eq|now apply Rgt_not_eq]. + rewrite ulp_neq_0 in Hx';[|now apply Rgt_not_eq]. + rewrite ulp_neq_0 in H;[|now apply Rgt_not_eq]. + lra. Qed. Lemma double_round_lt_mid_same_place : @@ -256,7 +266,7 @@ assert (H : Rabs (x * bpow (- fexp1 (ln_beta x)) - rewrite <- (Rmult_0_r (/ 2)). apply Rmult_lt_compat_l; [lra|]. apply bpow_gt_0. - - exact Hx. } + - rewrite ulp_neq_0 in Hx;try apply Rgt_not_eq; assumption. } unfold round at 2. unfold F2R, scaled_mantissa, canonic_exp; simpl. rewrite Hf2f1. @@ -320,8 +330,10 @@ unfold double_round_eq. set (x' := round beta fexp1 Zceil x). set (x'' := round beta fexp2 (Znearest choice2) x). intros Hx1 Hx2. -assert (Hr1 : Rabs (x'' - x) <= / 2 * bpow (fexp2 (ln_beta x))); - [now unfold x''; apply ulp_half_error|]. +assert (Hr1 : Rabs (x'' - x) <= / 2 * bpow (fexp2 (ln_beta x))). + apply Rle_trans with (/2* ulp beta fexp2 x). + now unfold x''; apply error_le_half_ulp... + rewrite ulp_neq_0;[now right|now apply Rgt_not_eq]. assert (Px'x : 0 <= x' - x). { apply Rle_0_minus. apply round_UP_pt. @@ -335,6 +347,7 @@ assert (Hr2 : Rabs (x'' - x') < / 2 * bpow (fexp1 (ln_beta x))). apply Rplus_le_lt_compat. - exact Hr1. - rewrite Rabs_minus_sym. + rewrite 2!ulp_neq_0 in Hx2; try (apply Rgt_not_eq; assumption). now rewrite Rabs_right; [|now apply Rle_ge]; apply Hx2. } destruct (Req_dec x'' 0) as [Zx''|Nzx'']. - (* x'' = 0 *) @@ -382,7 +395,8 @@ destruct (Req_dec x'' 0) as [Zx''|Nzx'']. apply (Rlt_le_trans _ _ _ Hx2). apply Rmult_le_compat_l; [lra|]. generalize (bpow_ge_0 beta (fexp2 (ln_beta x))). - unfold ulp, canonic_exp; lra. + rewrite 2!ulp_neq_0; try (apply Rgt_not_eq; assumption). + unfold canonic_exp; lra. + apply (Rmult_lt_reg_r (bpow (fexp1 (ln_beta x)))); [now apply bpow_gt_0|]. rewrite <- (Rabs_right (bpow (fexp1 _))) at 1; [|now apply Rle_ge; apply bpow_ge_0]. @@ -422,7 +436,7 @@ assert (Hx''pow : x'' = bpow (ln_beta x)). { apply Rle_lt_trans with (x + / 2 * ulp beta fexp2 x). - apply (Rplus_le_reg_r (- x)); ring_simplify. apply Rabs_le_inv. - apply ulp_half_error. + apply error_le_half_ulp. exact Vfexp2. - apply Rplus_lt_compat_r. rewrite <- Rabs_right at 1; [|now apply Rle_ge; apply Rlt_le]. @@ -442,15 +456,17 @@ assert (Hx''pow : x'' = bpow (ln_beta x)). apply (Rlt_le_trans _ _ _ H'x''). apply Rplus_le_compat_l. rewrite <- (Rmult_1_l (Fcore_Raux.bpow _ _)). + rewrite ulp_neq_0;[idtac|now apply Rgt_not_eq]. apply Rmult_le_compat_r; [now apply bpow_ge_0|]. lra. } assert (Hr : Rabs (x - x'') < / 2 * ulp beta fexp1 x). { apply Rle_lt_trans with (/ 2 * ulp beta fexp2 x). - rewrite Rabs_minus_sym. - apply ulp_half_error. + apply error_le_half_ulp. exact Vfexp2. - apply Rmult_lt_compat_l; [lra|]. - unfold ulp, canonic_exp; apply bpow_lt. + rewrite 2!ulp_neq_0; try now apply Rgt_not_eq. + unfold canonic_exp; apply bpow_lt. omega. } unfold round, F2R, scaled_mantissa, canonic_exp; simpl. assert (Hf : (0 <= ln_beta x - fexp1 (ln_beta x''))%Z). @@ -475,6 +491,7 @@ rewrite (Znearest_imp _ _ (beta ^ (ln_beta x - fexp1 (ln_beta x'')))%Z). rewrite <- Rabs_mult. rewrite Rmult_minus_distr_r. bpow_simplify. + rewrite ulp_neq_0 in Hr;[idtac|now apply Rgt_not_eq]. rewrite <- Hx''pow; exact Hr. - rewrite Z2R_Zpower; [|exact Hf]. apply (Rmult_lt_reg_r (bpow (fexp1 (ln_beta x'')))); [now apply bpow_gt_0|]. @@ -522,7 +539,7 @@ assert (H : Rabs (Z2R (Zceil (x * bpow (- fexp1 (ln_beta x)))) + apply Rle_0_minus. apply round_UP_pt. exact Vfexp1. - - exact Hx. } + - rewrite ulp_neq_0 in Hx;[exact Hx|now apply Rgt_not_eq]. } unfold double_round_eq, round at 2. unfold F2R, scaled_mantissa, canonic_exp; simpl. rewrite Hf2f1. @@ -761,9 +778,10 @@ destruct (Req_dec y 0) as [Zy|Nzy]. - (* y = 0 *) now rewrite Zy; rewrite Rplus_0_r. - (* y <> 0 *) - apply (ln_beta_succ beta fexp); [assumption|assumption|]. + apply (ln_beta_plus_eps beta fexp); [assumption|assumption|]. split; [assumption|]. - unfold ulp, canonic_exp. + rewrite ulp_neq_0;[idtac|now apply Rgt_not_eq]. + unfold canonic_exp. destruct (ln_beta y) as (ey, Hey); simpl in *. apply Rlt_le_trans with (bpow ey). + now rewrite <- (Rabs_right y); [apply Hey|apply Rle_ge]. @@ -797,8 +815,7 @@ apply ln_beta_unique. split. - apply Rabs_ge; right. assert (Hy : y < ulp beta fexp (bpow (ln_beta x - 1))). - { unfold ulp, canonic_exp. - rewrite ln_beta_bpow. + { rewrite ulp_bpow. replace (_ + _)%Z with (ln_beta x : Z) by ring. rewrite <- (Rabs_right y); [|now apply Rle_ge; apply Rlt_le]. apply Rlt_le_trans with (bpow (ln_beta y)). @@ -808,7 +825,8 @@ split. apply Rle_trans with (bpow (ln_beta x - 1) + ulp beta fexp (bpow (ln_beta x - 1))). + now apply Rplus_le_compat_l; apply Rlt_le. - + apply succ_le_lt; [|exact Fx|now split; [apply bpow_gt_0|]]. + + rewrite <- succ_eq_pos;[idtac|apply bpow_ge_0]. + apply succ_le_lt; [apply Vfexp|idtac|exact Fx|assumption]. apply (generic_format_bpow beta fexp (ln_beta x - 1)). replace (_ + _)%Z with (ln_beta x : Z) by ring. assert (fexp (ln_beta x) < ln_beta x)%Z; [|omega]. @@ -1039,13 +1057,15 @@ apply double_round_lt_mid. apply (Rlt_le_trans _ _ _ (proj2 (double_round_plus_aux1_aux 2 P2 fexp1 x y Px Py Hly Lxy Fx))). ring_simplify. - unfold ulp, canonic_exp; rewrite Lxy. + rewrite ulp_neq_0;[idtac|now apply Rgt_not_eq, Rplus_lt_0_compat]. + unfold canonic_exp; rewrite Lxy. apply (Rmult_le_reg_r (bpow (- fexp1 (ln_beta x)))); [now apply bpow_gt_0|]. bpow_simplify. apply (Rle_trans _ _ _ Bpow2). rewrite <- (Rmult_1_r (/ 2)) at 3. apply Rmult_le_compat_l; lra. -- unfold ulp, round, F2R, scaled_mantissa, canonic_exp; simpl; rewrite Lxy. +- rewrite ulp_neq_0;[idtac|now apply Rgt_not_eq, Rplus_lt_0_compat]. + unfold round, F2R, scaled_mantissa, canonic_exp; simpl; rewrite Lxy. intro Hf2'. apply (Rmult_lt_reg_r (bpow (- fexp1 (ln_beta x)))); [now apply bpow_gt_0|]. @@ -1056,7 +1076,8 @@ apply double_round_lt_mid. apply (Rlt_le_trans _ _ _ (proj2 (double_round_plus_aux1_aux 2 P2 fexp1 x y Px Py Hly Lxy Fx))). apply (Rmult_le_reg_r (bpow (- fexp1 (ln_beta x)))); [now apply bpow_gt_0|]. - unfold ulp, canonic_exp; rewrite Lxy, Rmult_minus_distr_r; bpow_simplify. + rewrite ulp_neq_0;[idtac|now apply Rgt_not_eq, Rplus_lt_0_compat]. + unfold canonic_exp; rewrite Lxy, Rmult_minus_distr_r; bpow_simplify. apply (Rle_trans _ _ _ Bpow2). rewrite <- (Rmult_1_r (/ 2)) at 3; rewrite <- Rmult_minus_distr_l. apply Rmult_le_compat_l; [lra|]. @@ -1391,7 +1412,8 @@ apply double_round_gt_mid. [now apply double_round_minus_aux2_aux; try assumption; omega|]. apply (Rlt_le_trans _ _ _ Ly). now apply bpow_le. - + unfold ulp, canonic_exp. + + rewrite ulp_neq_0;[idtac|now apply sym_not_eq, Rlt_not_eq, Rgt_minus]. + unfold canonic_exp. replace (_ - 2)%Z with (fexp1 (ln_beta (x - y)) - 1 - 1)%Z by ring. unfold Zminus at 1; rewrite bpow_plus. rewrite Rmult_comm. @@ -1423,7 +1445,8 @@ apply double_round_gt_mid. + unfold Fcore_Raux.bpow, Z.pow_pos; simpl. rewrite Zmult_1_r; apply Rinv_le; [lra|]. now change 2 with (Z2R 2); apply Z2R_le. - + unfold ulp, canonic_exp. + + rewrite 2!ulp_neq_0; try now apply Rgt_not_eq, Rgt_minus. + unfold canonic_exp. apply (Rplus_le_reg_r (bpow (fexp2 (ln_beta (x - y))))); ring_simplify. apply Rle_trans with (2 * bpow (fexp1 (ln_beta (x - y)) - 1)). * rewrite Rmult_plus_distr_r; rewrite Rmult_1_l. @@ -1868,19 +1891,22 @@ apply double_round_lt_mid. apply (Rlt_le_trans _ _ _ (proj2 (double_round_plus_aux1_aux 1 P1 fexp1 x y Px Py Hly Lxy Fx))). ring_simplify. - unfold ulp, canonic_exp; rewrite Lxy. + rewrite ulp_neq_0;[idtac|now apply Rgt_not_eq, Rplus_lt_0_compat]. + unfold canonic_exp; rewrite Lxy. apply (Rmult_le_reg_r (bpow (- fexp1 (ln_beta x)))); [now apply bpow_gt_0|]. bpow_simplify. apply (Rle_trans _ _ _ Bpow3); lra. -- unfold ulp, round, F2R, scaled_mantissa, canonic_exp; simpl; rewrite Lxy. +- rewrite ulp_neq_0;[idtac|now apply Rgt_not_eq, Rplus_lt_0_compat]. + unfold round, F2R, scaled_mantissa, canonic_exp; simpl; rewrite Lxy. intro Hf2'. unfold midp. apply (Rplus_lt_reg_r (- round beta fexp1 Zfloor (x + y))); ring_simplify. rewrite <- Rmult_minus_distr_l. apply (Rlt_le_trans _ _ _ (proj2 (double_round_plus_aux1_aux 1 P1 fexp1 x y Px Py Hly Lxy Fx))). - unfold ulp, canonic_exp; rewrite Lxy. + rewrite ulp_neq_0;[idtac|now apply Rgt_not_eq, Rplus_lt_0_compat]. + unfold canonic_exp; rewrite Lxy. apply (Rmult_le_reg_r (bpow (- fexp1 (ln_beta x)))); [now apply bpow_gt_0|]. rewrite (Rmult_assoc (/ 2)). @@ -2106,7 +2132,8 @@ apply double_round_gt_mid. [now apply double_round_minus_aux2_aux|]. apply (Rlt_le_trans _ _ _ Ly). now apply bpow_le. - + unfold ulp, canonic_exp. + + rewrite ulp_neq_0;[idtac|now apply Rgt_not_eq, Rgt_minus]. + unfold canonic_exp. unfold Zminus at 1; rewrite bpow_plus. rewrite Rmult_comm. apply Rmult_le_compat_r; [now apply bpow_ge_0|]. @@ -2124,7 +2151,8 @@ apply double_round_gt_mid. apply (Rlt_le_trans _ _ _ Ly). apply Rle_trans with (bpow (fexp1 (ln_beta (x - y)) - 1)); [now apply bpow_le|]. - unfold ulp, canonic_exp. + rewrite 2!ulp_neq_0; try now apply Rgt_not_eq, Rgt_minus. + unfold canonic_exp. apply (Rmult_le_reg_r (bpow (- fexp1 (ln_beta (x - y))))); [now apply bpow_gt_0|]. rewrite Rmult_assoc. @@ -2533,7 +2561,7 @@ destruct (generic_format_EM beta fexp1 x) as [Fx|Nfx]. now apply (generic_inclusion_ln_beta beta fexp1); [omega|]. - (* ~ generic_format beta fexp1 x *) assert (Hceil : round beta fexp1 Zceil x = rd + u1); - [now apply ulp_DN_UP|]. + [now apply round_UP_DN_ulp|]. assert (Hf2' : (fexp2 (ln_beta x) <= fexp1 (ln_beta x) - 1)%Z); [omega|]. destruct (Rlt_or_le (x - rd) (/ 2 * (u1 - u2))). + (* x - rd < / 2 * (u1 - u2) *) @@ -2589,10 +2617,11 @@ 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 := ulp beta fexp1 (sqrt x)). -set (u2 := ulp beta fexp2 (sqrt x)). +set (u1 := bpow (fexp1 (ln_beta (sqrt x)))). +set (u2 := bpow (fexp2 (ln_beta (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. @@ -2619,7 +2648,7 @@ assert (Pb : 0 < b). rewrite <- (Rmult_0_r (/ 2)). apply Rmult_lt_compat_l; [lra|]. apply Rlt_Rminus. - unfold u2, u1, ulp, canonic_exp. + unfold u2, u1. apply bpow_lt. omega. } assert (Pb' : 0 < b'). @@ -2686,7 +2715,7 @@ destruct (Req_dec a 0) as [Za|Nza]. - (* a <> 0 *) assert (Pa : 0 < a); [lra|]. assert (Hla : (ln_beta a = ln_beta (sqrt x) :> Z)). - { unfold a; apply ln_beta_round_DN. + { unfold a; apply ln_beta_DN. - exact Vfexp1. - now fold a. } assert (Hl' : 0 < - (u2 * a) + b * b). @@ -2697,12 +2726,14 @@ destruct (Req_dec a 0) as [Za|Nza]. replace (_ / 8) with (/ 4 * (u2 ^ 2 + u1 ^ 2)) by field. apply Rlt_le_trans with (u2 * bpow (ln_beta (sqrt x))). - apply Rmult_lt_compat_l; [now unfold u2, ulp; apply bpow_gt_0|]. - unfold u1, ulp, canonic_exp; rewrite <- Hla. - apply Rlt_le_trans with (a + ulp beta fexp1 a). + unfold u1; rewrite <- Hla. + apply Rlt_le_trans with (a + bpow (fexp1 (ln_beta a))). + apply Rplus_lt_compat_l. - rewrite <- (Rmult_1_l (ulp _ _ _)). + rewrite <- (Rmult_1_l (bpow _)) at 2. apply Rmult_lt_compat_r; [apply bpow_gt_0|lra]. - + apply (succ_le_bpow _ _ _ _ Pa Fa). + + 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_ln_beta_gt. - apply Rle_trans with (bpow (- 2) * u1 ^ 2). + unfold pow; rewrite Rmult_1_r. @@ -2745,9 +2776,10 @@ destruct (Req_dec a 0) as [Za|Nza]. apply Rlt_le_trans with (bpow (ln_beta (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. - rewrite <- ulp_DN; [|exact Vfexp1|exact Pa]; fold a. - apply succ_le_bpow. + unfold u1; fold (canonic_exp beta fexp1 (sqrt x)). + rewrite <- canonic_exp_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). @@ -2782,21 +2814,6 @@ destruct (Req_dec a 0) as [Za|Nza]. + now apply Rle_trans with x. Qed. -(* --> Fcore_Raux *) -Lemma sqrt_neg : forall x, x <= 0 -> sqrt x = 0. -Proof. -intros x Npx. -destruct (Req_dec x 0) as [Zx|Nzx]. -- (* x = 0 *) - rewrite Zx. - exact sqrt_0. -- (* x < 0 *) - unfold sqrt. - destruct Rcase_abs. - + reflexivity. - + casetype False. - now apply Nzx, Rle_antisym; [|apply Rge_le]. -Qed. Lemma double_round_sqrt : forall fexp1 fexp2 : Z -> Z, @@ -3028,10 +3045,11 @@ Lemma double_round_sqrt_beta_ge_4_aux : Proof. intros Hbeta fexp1 fexp2 Vfexp1 Vfexp2 Hexp x Px Hf2 Fx. set (a := round beta fexp1 Zfloor (sqrt x)). -set (u1 := ulp beta fexp1 (sqrt x)). -set (u2 := ulp beta fexp2 (sqrt x)). +set (u1 := bpow (fexp1 (ln_beta (sqrt x)))). +set (u2 := bpow (fexp2 (ln_beta (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. @@ -3125,7 +3143,7 @@ destruct (Req_dec a 0) as [Za|Nza]. - (* a <> 0 *) assert (Pa : 0 < a); [lra|]. assert (Hla : (ln_beta a = ln_beta (sqrt x) :> Z)). - { unfold a; apply ln_beta_round_DN. + { unfold a; apply ln_beta_DN. - exact Vfexp1. - now fold a. } assert (Hl' : 0 < - (u2 * a) + b * b). @@ -3136,12 +3154,13 @@ destruct (Req_dec a 0) as [Za|Nza]. replace (_ / 8) with (/ 4 * (u2 ^ 2 + u1 ^ 2)) by field. apply Rlt_le_trans with (u2 * bpow (ln_beta (sqrt x))). - apply Rmult_lt_compat_l; [now unfold u2, ulp; apply bpow_gt_0|]. - unfold u1, ulp, canonic_exp; rewrite <- Hla. + unfold u1; rewrite <- Hla. apply Rlt_le_trans with (a + ulp beta fexp1 a). + apply Rplus_lt_compat_l. rewrite <- (Rmult_1_l (ulp _ _ _)). + rewrite ulp_neq_0; trivial. apply Rmult_lt_compat_r; [apply bpow_gt_0|lra]. - + apply (succ_le_bpow _ _ _ _ Pa Fa). + + apply (id_p_ulp_le_bpow _ _ _ _ Pa Fa). apply Rabs_lt_inv, bpow_ln_beta_gt. - apply Rle_trans with (bpow (- 1) * u1 ^ 2). + unfold pow; rewrite Rmult_1_r. @@ -3184,9 +3203,10 @@ destruct (Req_dec a 0) as [Za|Nza]. apply Rlt_le_trans with (bpow (ln_beta (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. - rewrite <- ulp_DN; [|exact Vfexp1|exact Pa]; fold a. - apply succ_le_bpow. + unfold u1; fold (canonic_exp beta fexp1 (sqrt x)). + rewrite <- canonic_exp_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). @@ -3504,9 +3524,11 @@ assert (Hf : F2R f = x). rewrite (Rmult_assoc _ (Z2R n)). rewrite Rinv_r; [rewrite Rmult_1_r|change 0 with (Z2R 0); apply Z2R_neq; omega]. - simpl; fold (canonic_exp beta fexp1 x); fold (ulp beta fexp1 x); fold u. - rewrite Xmid at 2. + simpl; fold (canonic_exp beta fexp1 x). + rewrite <- 2!ulp_neq_0; try now apply Rgt_not_eq. + fold u; rewrite Xmid at 2. apply f_equal2; [|reflexivity]. + rewrite ulp_neq_0; try now apply Rgt_not_eq. destruct (Req_dec rd 0) as [Zrd|Nzrd]. - (* rd = 0 *) rewrite Zrd. @@ -3657,7 +3679,7 @@ split. replace (bpow _) with (bpow (ln_beta x) - / 2 * u2 + / 2 * u2) by ring. apply Rplus_lt_le_compat; [exact Hx|]. apply Rabs_le_inv. - now apply ulp_half_error. + now apply error_le_half_ulp. Qed. Lemma double_round_all_mid_cases : @@ -3714,7 +3736,7 @@ destruct (Ztrichotomy (ln_beta x) (fexp1 (ln_beta x) - 1)) as [Hlt|[Heq|Hgt]]. now apply (generic_inclusion_ln_beta beta fexp1); [omega|]. - (* ~ generic_format beta fexp1 x *) assert (Hceil : round beta fexp1 Zceil x = x' + u1); - [now apply ulp_DN_UP|]. + [now apply round_UP_DN_ulp|]. assert (Hf2' : (fexp2 (ln_beta x) <= fexp1 (ln_beta x) - 1)%Z); [omega|]. assert (midp' fexp1 x + / 2 * ulp beta fexp2 x < x); @@ -3769,12 +3791,15 @@ assert (Hfx : (fexp1 (ln_beta x) < ln_beta x)%Z); assert (Hfy : (fexp1 (ln_beta y) < ln_beta y)%Z); [now apply ln_beta_generic_gt; [|apply Rgt_not_eq|]|]. set (p := bpow (ln_beta (x / y))). -set (u2 := ulp beta fexp2 (x / y)). +set (u2 := bpow (fexp2 (ln_beta (x / y)))). revert Fx Fy. unfold generic_format, F2R, scaled_mantissa, canonic_exp; simpl. set (mx := Ztrunc (x * bpow (- fexp1 (ln_beta x)))). set (my := Ztrunc (y * bpow (- fexp1 (ln_beta y)))). intros Fx Fy. +rewrite ulp_neq_0. +2: apply Rmult_integral_contrapositive_currified; [now apply Rgt_not_eq|idtac]. +2: now apply Rinv_neq_0_compat, Rgt_not_eq. intro Hl. assert (Hr : x / y < p); [now apply Rabs_lt_inv; apply bpow_ln_beta_gt|]. @@ -3903,6 +3928,9 @@ assert (Hfx : (fexp1 (ln_beta x) < ln_beta x)%Z); [now apply ln_beta_generic_gt; [|apply Rgt_not_eq|]|]. assert (Hfy : (fexp1 (ln_beta y) < ln_beta y)%Z); [now apply ln_beta_generic_gt; [|apply Rgt_not_eq|]|]. +assert (S : (x / y <> 0)%R). +apply Rmult_integral_contrapositive_currified; [now apply Rgt_not_eq|idtac]. +now apply Rinv_neq_0_compat, Rgt_not_eq. cut (~ (/ 2 * (ulp beta fexp1 (x / y) - ulp beta fexp2 (x / y)) <= x / y - round beta fexp1 Zfloor (x / y) < / 2 * ulp beta fexp1 (x / y))). @@ -3913,9 +3941,10 @@ cut (~ (/ 2 * (ulp beta fexp1 (x / y) - ulp beta fexp2 (x / y)) - apply (Rplus_lt_reg_l (round beta fexp1 Zfloor (x / y))). ring_simplify. apply H'. } -set (u1 := ulp beta fexp1 (x / y)). -set (u2 := ulp beta fexp2 (x / y)). +set (u1 := bpow (fexp1 (ln_beta (x / y)))). +set (u2 := bpow (fexp2 (ln_beta (x / y)))). set (x' := round beta fexp1 Zfloor (x / y)). +rewrite 2!ulp_neq_0; trivial. revert Fx Fy. unfold generic_format, F2R, scaled_mantissa, canonic_exp; simpl. set (mx := Ztrunc (x * bpow (- fexp1 (ln_beta x)))). @@ -4098,9 +4127,13 @@ cut (~ (/ 2 * ulp beta fexp1 (x / y) - apply (Rplus_le_reg_l (round beta fexp1 Zfloor (x / y))). ring_simplify. apply H'. } -set (u1 := ulp beta fexp1 (x / y)). -set (u2 := ulp beta fexp2 (x / y)). +set (u1 := bpow (fexp1 (ln_beta (x / y)))). +set (u2 := bpow (fexp2 (ln_beta (x / y)))). set (x' := round beta fexp1 Zfloor (x / y)). +assert (S : (x / y <> 0)%R). +apply Rmult_integral_contrapositive_currified; [now apply Rgt_not_eq|idtac]. +now apply Rinv_neq_0_compat, Rgt_not_eq. +rewrite 2!ulp_neq_0; trivial. revert Fx Fy. unfold generic_format, F2R, scaled_mantissa, canonic_exp; simpl. set (mx := Ztrunc (x * bpow (- fexp1 (ln_beta x)))). diff --git a/flocq/Appli/Fappli_rnd_odd.v b/flocq/Appli/Fappli_rnd_odd.v index b3244589..4741047f 100644 --- a/flocq/Appli/Fappli_rnd_odd.v +++ b/flocq/Appli/Fappli_rnd_odd.v @@ -356,6 +356,80 @@ simpl; ring. apply Rgt_not_eq, bpow_gt_0. Qed. + + +Theorem Rnd_odd_pt_unicity : + forall x f1 f2 : R, + Rnd_odd_pt x f1 -> Rnd_odd_pt x f2 -> + f1 = f2. +Proof. +intros x f1 f2 (Ff1,H1) (Ff2,H2). +(* *) +case (generic_format_EM beta fexp x); intros L. +apply trans_eq with x. +case H1; try easy. +intros (J,_); case J; intros J'. +apply Rnd_DN_pt_idempotent with format; assumption. +apply Rnd_UP_pt_idempotent with format; assumption. +case H2; try easy. +intros (J,_); case J; intros J'; apply sym_eq. +apply Rnd_DN_pt_idempotent with format; assumption. +apply Rnd_UP_pt_idempotent with format; assumption. +(* *) +destruct H1 as [H1|(H1,H1')]. +contradict L; now rewrite <- H1. +destruct H2 as [H2|(H2,H2')]. +contradict L; now rewrite <- H2. +destruct H1 as [H1|H1]; destruct H2 as [H2|H2]. +apply Rnd_DN_pt_unicity with format x; assumption. +destruct H1' as (ff,(K1,(K2,K3))). +destruct H2' as (gg,(L1,(L2,L3))). +absurd (true = false); try discriminate. +rewrite <- L3. +apply trans_eq with (negb (Zeven (Fnum ff))). +rewrite K3; easy. +apply sym_eq. +generalize (DN_UP_parity_generic beta fexp). +unfold DN_UP_parity_prop; intros T; apply (T x); clear T; try assumption... +rewrite <- K1; apply Rnd_DN_pt_unicity with (generic_format beta fexp) x; try easy... +now apply round_DN_pt... +rewrite <- L1; apply Rnd_UP_pt_unicity with (generic_format beta fexp) x; try easy... +now apply round_UP_pt... +(* *) +destruct H1' as (ff,(K1,(K2,K3))). +destruct H2' as (gg,(L1,(L2,L3))). +absurd (true = false); try discriminate. +rewrite <- K3. +apply trans_eq with (negb (Zeven (Fnum gg))). +rewrite L3; easy. +apply sym_eq. +generalize (DN_UP_parity_generic beta fexp). +unfold DN_UP_parity_prop; intros T; apply (T x); clear T; try assumption... +rewrite <- L1; apply Rnd_DN_pt_unicity with (generic_format beta fexp) x; try easy... +now apply round_DN_pt... +rewrite <- K1; apply Rnd_UP_pt_unicity with (generic_format beta fexp) x; try easy... +now apply round_UP_pt... +apply Rnd_UP_pt_unicity with format x; assumption. +Qed. + + + +Theorem Rnd_odd_pt_monotone : + round_pred_monotone (Rnd_odd_pt). +Proof with auto with typeclass_instances. +intros x y f g H1 H2 Hxy. +apply Rle_trans with (round beta fexp Zrnd_odd x). +right; apply Rnd_odd_pt_unicity with x; try assumption. +apply round_odd_pt. +apply Rle_trans with (round beta fexp Zrnd_odd y). +apply round_le... +right; apply Rnd_odd_pt_unicity with y; try assumption. +apply round_odd_pt. +Qed. + + + + End Fcore_rnd_odd. Section Odd_prop_aux. @@ -462,7 +536,7 @@ Lemma ln_beta_d: (0< F2R d)%R -> (ln_beta beta (F2R d) = ln_beta beta x :>Z). Proof with auto with typeclass_instances. intros Y. -rewrite d_eq; apply ln_beta_round_DN... +rewrite d_eq; apply ln_beta_DN... now rewrite <- d_eq. Qed. @@ -861,13 +935,9 @@ case K2; clear K2; intros K2. case (Rle_or_lt x m); intros Y;[destruct Y|idtac]. (* . *) apply trans_eq with (F2R d). -apply round_N_DN_betw with (F2R u)... +apply round_N_eq_DN_pt with (F2R u)... apply DN_odd_d_aux; split; try left; assumption. apply UP_odd_d_aux; split; try left; assumption. -split. -apply round_ge_generic... -apply generic_format_fexpe_fexp, Hd. -apply Hd. assert (o <= (F2R d + F2R u) / 2)%R. apply round_le_generic... apply Fm. @@ -876,10 +946,7 @@ destruct H1; trivial. apply P. now apply Rlt_not_eq. trivial. -apply sym_eq, round_N_DN_betw with (F2R u)... -split. -apply Hd. -exact H0. +apply sym_eq, round_N_eq_DN_pt with (F2R u)... (* . *) replace o with x. reflexivity. @@ -887,10 +954,9 @@ apply sym_eq, round_generic... rewrite H0; apply Fm. (* . *) apply trans_eq with (F2R u). -apply round_N_UP_betw with (F2R d)... +apply round_N_eq_UP_pt with (F2R d)... apply DN_odd_d_aux; split; try left; assumption. apply UP_odd_d_aux; split; try left; assumption. -split. assert ((F2R d + F2R u) / 2 <= o)%R. apply round_ge_generic... apply Fm. @@ -899,13 +965,7 @@ destruct H0; trivial. apply P. now apply Rgt_not_eq. rewrite <- H0; trivial. -apply round_le_generic... -apply generic_format_fexpe_fexp, Hu. -apply Hu. -apply sym_eq, round_N_UP_betw with (F2R d)... -split. -exact Y. -apply Hu. +apply sym_eq, round_N_eq_UP_pt with (F2R d)... Qed. -- cgit