diff options
Diffstat (limited to 'flocq/Appli')
| -rw-r--r-- | flocq/Appli/Fappli_IEEE.v | 227 | ||||
| -rw-r--r-- | flocq/Appli/Fappli_IEEE_bits.v | 215 | ||||
| -rw-r--r-- | flocq/Appli/Fappli_double_round.v | 4554 | ||||
| -rw-r--r-- | flocq/Appli/Fappli_rnd_odd.v | 4 |
4 files changed, 4889 insertions, 111 deletions
diff --git a/flocq/Appli/Fappli_IEEE.v b/flocq/Appli/Fappli_IEEE.v index b44d711c..9b5826c1 100644 --- a/flocq/Appli/Fappli_IEEE.v +++ b/flocq/Appli/Fappli_IEEE.v @@ -46,6 +46,8 @@ End AnyRadix. 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 *) @@ -59,7 +61,7 @@ Instance fexp_correct : Valid_exp fexp := FLT_exp_valid emin prec. Instance fexp_monotone : Monotone_exp fexp := FLT_exp_monotone emin prec. Definition canonic_mantissa m e := - Zeq_bool (fexp (Z_of_nat (S (digits2_Pnat m)) + e)) e. + Zeq_bool (fexp (Zpos (digits2_pos m) + e)) e. Definition bounded m e := andb (canonic_mantissa m e) (Zle_bool e (emax - prec)). @@ -67,14 +69,15 @@ Definition bounded m e := Definition valid_binary x := match x with | F754_finite _ m e => bounded m e - | F754_nan _ pl => (Z_of_nat' (S (digits2_Pnat pl)) <? prec)%Z + | F754_nan _ pl => (Zpos (digits2_pos pl) <? prec)%Z | _ => true end. (** Basic type used for representing binary FP numbers. - Note that there is exactly one such object per FP datum. *) + Note that there is exactly one such object per FP datum. + NaNs do not have any payload. They cannot be distinguished. *) -Definition nan_pl := {pl | (Z_of_nat' (S (digits2_Pnat pl)) <? prec)%Z = true}. +Definition nan_pl := {pl | (Zpos (digits2_pos pl) <? prec)%Z = true}. Inductive binary_float := | B754_zero : bool -> binary_float @@ -88,7 +91,7 @@ Definition FF2B x := | F754_finite s m e => B754_finite s m e | F754_infinity s => fun _ => B754_infinity s | F754_zero s => fun _ => B754_zero s - | F754_nan b pl => fun H => B754_nan b (exist _ pl H) + | F754_nan b pl => fun H => B754_nan b (exist pl H) end. Definition B2FF x := @@ -99,8 +102,6 @@ Definition B2FF x := | B754_nan b (exist pl _) => F754_nan b pl end. -Definition radix2 := Build_radix 2 (refl_equal true). - Definition B2R f := match f with | B754_finite s m e _ => F2R (Float radix2 (cond_Zopp s (Zpos m)) e) @@ -183,7 +184,7 @@ pattern ex at 2 ; rewrite <- Hx. apply (f_equal fexp). rewrite ln_beta_F2R_Zdigits. rewrite <- Zdigits_abs. -rewrite Z_of_nat_S_digits2_Pnat. +rewrite Zpos_digits2_pos. now case sx. now case sx. Qed. @@ -419,6 +420,159 @@ simpl. now case opp_nan. Qed. +(** Absolute value *) + +Definition Babs abs_nan (x : binary_float) : binary_float := + match x with + | B754_nan sx plx => + let '(sres, plres) := abs_nan sx plx in B754_nan sres plres + | B754_infinity sx => B754_infinity false + | B754_finite sx mx ex Hx => B754_finite false mx ex Hx + | B754_zero sx => B754_zero false + end. + +Theorem B2R_Babs : + forall abs_nan x, + B2R (Babs abs_nan x) = Rabs (B2R x). +Proof. + intros abs_nan [sx|sx|sx plx|sx mx ex Hx]; apply sym_eq ; try apply Rabs_R0. + simpl. destruct abs_nan. simpl. apply Rabs_R0. + simpl. rewrite <- F2R_abs. now destruct sx. +Qed. + +Theorem is_finite_Babs : + forall abs_nan x, + is_finite (Babs abs_nan x) = is_finite x. +Proof. + intros abs_nan [| | |] ; try easy. + intros s pl. + simpl. + now case abs_nan. +Qed. + +Theorem Bsign_Babs : + forall abs_nan x, + is_nan x = false -> + Bsign (Babs abs_nan x) = false. +Proof. + now intros abs_nan [| | |]. +Qed. + +Theorem Babs_idempotent : + forall abs_nan (x: binary_float), + is_nan x = false -> + Babs abs_nan (Babs abs_nan x) = Babs abs_nan x. +Proof. + now intros abs_nan [sx|sx|sx plx|sx mx ex Hx]. +Qed. + +Theorem Babs_Bopp : + forall abs_nan opp_nan x, + is_nan x = false -> + Babs abs_nan (Bopp opp_nan x) = Babs abs_nan x. +Proof. + now intros abs_nan opp_nan [| | |]. +Qed. + +(** Comparison + +[Some c] means ordered as per [c]; [None] means unordered. *) + +Definition Bcompare (f1 f2 : binary_float) : option comparison := + match f1, f2 with + | B754_nan _ _,_ | _,B754_nan _ _ => None + | B754_infinity true, B754_infinity true + | B754_infinity false, B754_infinity false => Some Eq + | B754_infinity true, _ => Some Lt + | B754_infinity false, _ => Some Gt + | _, B754_infinity true => Some Gt + | _, B754_infinity false => Some Lt + | B754_finite true _ _ _, B754_zero _ => Some Lt + | B754_finite false _ _ _, B754_zero _ => Some Gt + | B754_zero _, B754_finite true _ _ _ => Some Gt + | B754_zero _, B754_finite false _ _ _ => Some Lt + | B754_zero _, B754_zero _ => Some Eq + | B754_finite s1 m1 e1 _, B754_finite s2 m2 e2 _ => + match s1, s2 with + | true, false => Some Lt + | false, true => Some Gt + | false, false => + match Zcompare e1 e2 with + | Lt => Some Lt + | Gt => Some Gt + | Eq => Some (Pcompare m1 m2 Eq) + end + | true, true => + match Zcompare e1 e2 with + | Lt => Some Gt + | Gt => Some Lt + | Eq => Some (CompOpp (Pcompare m1 m2 Eq)) + end + end + end. + +Theorem Bcompare_correct : + forall f1 f2, + is_finite f1 = true -> is_finite f2 = true -> + Bcompare f1 f2 = Some (Rcompare (B2R f1) (B2R f2)). +Proof. + Ltac apply_Rcompare := + match goal with + | [ |- Some Lt = Some (Rcompare _ _) ] => f_equal; symmetry; apply Rcompare_Lt + | [ |- Some Eq = Some (Rcompare _ _) ] => f_equal; symmetry; apply Rcompare_Eq + | [ |- Some Gt = Some (Rcompare _ _) ] => f_equal; symmetry; apply Rcompare_Gt + end. + unfold Bcompare; intros. + destruct f1, f2 ; try easy. + now rewrite Rcompare_Eq. + destruct b0 ; apply_Rcompare. + now apply F2R_lt_0_compat. + now apply F2R_gt_0_compat. + destruct b ; apply_Rcompare. + now apply F2R_lt_0_compat. + now apply F2R_gt_0_compat. + simpl. + clear H H0. + apply andb_prop in e0; destruct e0; apply (canonic_canonic_mantissa false) in H. + apply andb_prop in e2; destruct e2; apply (canonic_canonic_mantissa false) in H1. + pose proof (Zcompare_spec e e1); unfold canonic, Fexp in H1, H. + assert (forall m1 m2 e1 e2, + let x := (Z2R (Zpos m1) * bpow radix2 e1)%R in + let y := (Z2R (Zpos m2) * bpow radix2 e2)%R in + (canonic_exp radix2 fexp x < canonic_exp radix2 fexp y)%Z -> (x < y)%R). + { + intros; apply Rnot_le_lt; intro; apply (ln_beta_le radix2) in H5. + apply Zlt_not_le with (1 := H4). + now apply fexp_monotone. + now apply (F2R_gt_0_compat _ (Float radix2 (Zpos m2) e2)). + } + assert (forall m1 m2 e1 e2, (Z2R (- Zpos m1) * bpow radix2 e1 < Z2R (Zpos m2) * bpow radix2 e2)%R). + { + intros; apply (Rlt_trans _ 0%R). + now apply (F2R_lt_0_compat _ (Float radix2 (Zneg m1) e0)). + now apply (F2R_gt_0_compat _ (Float radix2 (Zpos m2) e2)). + } + unfold F2R, Fnum, Fexp. + destruct b, b0; try (now apply_Rcompare; apply H5); inversion H3; + try (apply_Rcompare; apply H4; rewrite H, H1 in H7; assumption); + try (apply_Rcompare; do 2 rewrite Z2R_opp, Ropp_mult_distr_l_reverse; + apply Ropp_lt_contravar; apply H4; rewrite H, H1 in H7; assumption); + rewrite H7, Rcompare_mult_r, Rcompare_Z2R by (apply bpow_gt_0); reflexivity. +Qed. + +Theorem Bcompare_swap : + forall x y, + Bcompare y x = match Bcompare x y with Some c => Some (CompOpp c) | None => None end. +Proof. + intros. + destruct x as [ ? | [] | ? ? | [] mx ex Bx ]; + destruct y as [ ? | [] | ? ? | [] my ey By ]; simpl; try easy. +- rewrite <- (Zcompare_antisym ex ey). destruct (ex ?= ey)%Z; try easy. + now rewrite (Pcompare_antisym mx my). +- rewrite <- (Zcompare_antisym ex ey). destruct (ex ?= ey)%Z; try easy. + now rewrite Pcompare_antisym. +Qed. + Theorem bounded_lt_emax : forall mx ex, bounded mx ex = true -> @@ -441,8 +595,7 @@ now apply F2R_gt_0_compat. apply bpow_le. rewrite H. 2: discriminate. revert H1. clear -H2. -rewrite Z_of_nat_S_digits2_Pnat. -change Fcalc_digits.radix2 with radix2. +rewrite Zpos_digits2_pos. unfold fexp, FLT_exp. generalize (Zdigits radix2 (Zpos mx)). intros ; zify ; subst. @@ -471,8 +624,7 @@ split. unfold canonic_mantissa. unfold canonic, Fexp in Cx. rewrite Cx at 2. -rewrite Z_of_nat_S_digits2_Pnat. -change Fcalc_digits.radix2 with radix2. +rewrite Zpos_digits2_pos. unfold canonic_exp. rewrite ln_beta_F2R_Zdigits. 2: discriminate. now apply -> Zeq_is_eq_bool. @@ -600,7 +752,7 @@ induction (nat_of_P n). simpl. rewrite Zplus_0_r. now destruct l as [|[| |]]. -simpl iter_nat. +simpl nat_rect. rewrite inj_S. unfold Zsucc. rewrite Zplus_assoc. @@ -617,18 +769,6 @@ intros (m, r, s) Hm. now destruct m as [|[m|m|]|m] ; try (now elim Hm) ; destruct r as [|] ; destruct s as [|]. Qed. -Definition Zdigits2 m := - match m with Z0 => m | Zpos p => Z_of_nat (S (digits2_Pnat p)) | Zneg p => Z_of_nat (S (digits2_Pnat p)) end. - -Theorem Zdigits2_Zdigits : - forall m, - Zdigits2 m = Zdigits radix2 m. -Proof. -unfold Zdigits2. -intros [|m|m] ; try apply Z_of_nat_S_digits2_Pnat. -easy. -Qed. - Definition shr_fexp m e l := shr (shr_record_of_loc m l) e (fexp (Zdigits2 m + e) - e). @@ -829,7 +969,7 @@ apply andb_true_intro. split. unfold canonic_mantissa. apply Zeq_bool_true. -rewrite Z_of_nat_S_digits2_Pnat. +rewrite Zpos_digits2_pos. rewrite <- ln_beta_F2R_Zdigits. apply sym_eq. now rewrite H3 in H4. @@ -849,8 +989,7 @@ unfold valid_binary, bounded. rewrite Zle_bool_refl. rewrite Bool.andb_true_r. apply Zeq_bool_true. -rewrite Z_of_nat_S_digits2_Pnat. -change Fcalc_digits.radix2 with radix2. +rewrite Zpos_digits2_pos. replace (Zdigits radix2 (Zpos (match (Zpower 2 prec - 1)%Z with Zpos p => p | _ => xH end))) with prec. unfold fexp, FLT_exp, emin. generalize (prec_gt_0 prec). @@ -942,7 +1081,7 @@ assert (forall m e, bounded m e = true -> fexp (Zdigits radix2 (Zpos m) + e) = e clear. intros m e Hb. destruct (andb_prop _ _ Hb) as (H,_). apply Zeq_bool_eq. -now rewrite <- Z_of_nat_S_digits2_Pnat. +now rewrite <- Zpos_digits2_pos. generalize (H _ _ Hx) (H _ _ Hy). clear x y sx sy Hx Hy H. unfold fexp, FLT_exp. @@ -1098,7 +1237,7 @@ apply refl_equal. Qed. Definition shl_align_fexp mx ex := - shl_align mx ex (fexp (Z_of_nat (S (digits2_Pnat mx)) + ex)). + shl_align mx ex (fexp (Zpos (digits2_pos mx) + ex)). Theorem shl_align_fexp_correct : forall mx ex, @@ -1108,8 +1247,8 @@ Theorem shl_align_fexp_correct : Proof. intros mx ex. unfold shl_align_fexp. -generalize (shl_align_correct mx ex (fexp (Z_of_nat (S (digits2_Pnat mx)) + ex))). -rewrite Z_of_nat_S_digits2_Pnat. +generalize (shl_align_correct mx ex (fexp (Zpos (digits2_pos mx) + ex))). +rewrite Zpos_digits2_pos. case shl_align. intros mx' ex' (H1, H2). split. @@ -1438,10 +1577,10 @@ Definition Fdiv_core_binary m1 e1 m2 e2 := let e := (e1 - e2)%Z in let (m, e') := match (d2 + prec - d1)%Z with - | Zpos p => (m1 * Zpower_pos 2 p, e + Zneg p)%Z + | Zpos p => (Z.shiftl m1 (Zpos p), e + Zneg p)%Z | _ => (m1, e) end in - let '(q, r) := Zdiv_eucl m m2 in + let '(q, r) := Zfast_div_eucl m m2 in (q, e', new_location m2 r loc_Exact). Lemma Bdiv_correct_aux : @@ -1463,7 +1602,6 @@ Lemma Bdiv_correct_aux : Proof. intros m sx mx ex sy my ey. replace (Fdiv_core_binary (Zpos mx) ex (Zpos my) ey) with (Fdiv_core radix2 prec (Zpos mx) ex (Zpos my) ey). -2: now unfold Fdiv_core_binary ; rewrite 2!Zdigits2_Zdigits. refine (_ (Fdiv_core_correct radix2 prec (Zpos mx) ex (Zpos my) ey _ _ _)) ; try easy. destruct (Fdiv_core radix2 prec (Zpos mx) ex (Zpos my) ey) as ((mz, ez), lz). intros (Pz, Bz). @@ -1540,6 +1678,15 @@ now apply F2R_ge_0_compat. apply Rlt_le. apply Rinv_0_lt_compat. now apply F2R_gt_0_compat. +(* *) +unfold Fdiv_core_binary, Fdiv_core. +rewrite 2!Zdigits2_Zdigits. +change 2%Z with (radix_val radix2). +destruct (Zdigits radix2 (Z.pos my) + prec - Zdigits radix2 (Z.pos mx))%Z as [|p|p]. +now rewrite Zfast_div_eucl_correct. +rewrite Z.shiftl_mul_pow2 by easy. +now rewrite Zfast_div_eucl_correct. +now rewrite Zfast_div_eucl_correct. Qed. Definition Bdiv div_nan m x y := @@ -1597,10 +1744,10 @@ Definition Fsqrt_core_binary m e := let (s', e'') := if Zeven e' then (s, e') else (s + 1, e' - 1)%Z in let m' := match s' with - | Zpos p => (m * Zpower_pos 2 p)%Z + | Zpos p => Z.shiftl m (Zpos p) | _ => m end in - let (q, r) := Zsqrt m' in + let (q, r) := Z.sqrtrem m' in let l := if Zeq_bool r 0 then loc_Exact else loc_Inexact (if Zle_bool r q then Lt else Gt) in @@ -1621,7 +1768,6 @@ Lemma Bsqrt_correct_aux : Proof with auto with typeclass_instances. intros m mx ex Hx. replace (Fsqrt_core_binary (Zpos mx) ex) with (Fsqrt_core radix2 prec (Zpos mx) ex). -2: now unfold Fsqrt_core_binary ; rewrite Zdigits2_Zdigits. simpl. refine (_ (Fsqrt_core_correct radix2 prec (Zpos mx) ex _)) ; try easy. destruct (Fsqrt_core radix2 prec (Zpos mx) ex) as ((mz, ez), lz). @@ -1716,6 +1862,11 @@ apply Rle_trans with R0. apply F2R_le_0_compat. now case mz. apply sqrt_ge_0. +(* *) +unfold Fsqrt_core, Fsqrt_core_binary. +rewrite Zdigits2_Zdigits. +destruct (if Zeven _ then _ else _) as [[|s'|s'] e''] ; try easy. +now rewrite Z.shiftl_mul_pow2. Qed. Definition Bsqrt sqrt_nan m x := diff --git a/flocq/Appli/Fappli_IEEE_bits.v b/flocq/Appli/Fappli_IEEE_bits.v index 937e8d43..5a77bf57 100644 --- a/flocq/Appli/Fappli_IEEE_bits.v +++ b/flocq/Appli/Fappli_IEEE_bits.v @@ -25,6 +25,12 @@ Require Import Fappli_IEEE. Section Binary_Bits. +Implicit Arguments exist [[A] [P]]. +Implicit Arguments B754_zero [[prec] [emax]]. +Implicit Arguments B754_infinity [[prec] [emax]]. +Implicit Arguments B754_nan [[prec] [emax]]. +Implicit Arguments B754_finite [[prec] [emax]]. + (** Number of bits for the fraction and exponent *) Variable mw ew : Z. Hypothesis Hmw : (0 < mw)%Z. @@ -54,7 +60,40 @@ Qed. Hypothesis Hmax : (prec < emax)%Z. Definition join_bits (s : bool) m e := - (((if s then Zpower 2 ew else 0) + e) * Zpower 2 mw + m)%Z. + (Z.shiftl ((if s then Zpower 2 ew else 0) + e) mw + m)%Z. + +Lemma join_bits_range : + forall s m e, + (0 <= m < 2^mw)%Z -> + (0 <= e < 2^ew)%Z -> + (0 <= join_bits s m e < 2 ^ (mw + ew + 1))%Z. +Proof. +intros s m e Hm He. +unfold join_bits. +rewrite Z.shiftl_mul_pow2 by now apply Zlt_le_weak. +split. +- apply (Zplus_le_compat 0 _ 0) with (2 := proj1 Hm). + rewrite <- (Zmult_0_l (2^mw)). + apply Zmult_le_compat_r. + case s. + clear -He ; omega. + now rewrite Zmult_0_l. + clear -Hm ; omega. +- apply Zlt_le_trans with (((if s then 2 ^ ew else 0) + e + 1) * 2 ^ mw)%Z. + rewrite (Zmult_plus_distr_l _ 1). + apply Zplus_lt_compat_l. + now rewrite Zmult_1_l. + rewrite <- (Zplus_assoc mw), (Zplus_comm mw), Zpower_plus. + apply Zmult_le_compat_r. + rewrite Zpower_plus. + change (2^1)%Z with 2%Z. + case s ; clear -He ; omega. + now apply Zlt_le_weak. + easy. + clear -Hm ; omega. + clear -Hew ; omega. + now apply Zlt_le_weak. +Qed. Definition split_bits x := let mm := Zpower 2 mw in @@ -69,6 +108,7 @@ Theorem split_join_bits : Proof. intros s m e Hm He. unfold split_bits, join_bits. +rewrite Z.shiftl_mul_pow2 by now apply Zlt_le_weak. apply f_equal2. apply f_equal2. (* *) @@ -114,6 +154,7 @@ Theorem join_split_bits : Proof. intros x Hx. unfold split_bits, join_bits. +rewrite Z.shiftl_mul_pow2 by now apply Zlt_le_weak. pattern x at 4 ; rewrite Z_div_mod_eq_full with x (2^mw)%Z. apply (f_equal (fun v => (v + _)%Z)). rewrite Zmult_comm. @@ -174,8 +215,9 @@ Definition bits_of_binary_float (x : binary_float) := | B754_infinity sx => join_bits sx 0 (Zpower 2 ew - 1) | B754_nan sx (exist plx _) => join_bits sx (Zpos plx) (Zpower 2 ew - 1) | B754_finite sx mx ex _ => - if Zle_bool (Zpower 2 mw) (Zpos mx) then - join_bits sx (Zpos mx - Zpower 2 mw) (ex - emin + 1) + let m := (Zpos mx - Zpower 2 mw)%Z in + if Zle_bool 0 m then + join_bits sx m (ex - emin + 1) else join_bits sx (Zpos mx) 0 end. @@ -186,8 +228,9 @@ Definition split_bits_of_binary_float (x : binary_float) := | B754_infinity sx => (sx, 0, Zpower 2 ew - 1)%Z | B754_nan sx (exist plx _) => (sx, Zpos plx, Zpower 2 ew - 1)%Z | B754_finite sx mx ex _ => - if Zle_bool (Zpower 2 mw) (Zpos mx) then - (sx, Zpos mx - Zpower 2 mw, ex - emin + 1)%Z + let m := (Zpos mx - Zpower 2 mw)%Z in + if Zle_bool 0 m then + (sx, m, ex - emin + 1)%Z else (sx, Zpos mx, 0)%Z end. @@ -200,6 +243,7 @@ intros [sx|sx|sx [plx Hplx]|sx mx ex Hx] ; try ( simpl ; apply split_join_bits ; split ; try apply Zle_refl ; try apply Zlt_pred ; trivial ; omega ). simpl. apply split_join_bits; split; try (zify; omega). destruct (digits2_Pnat_correct plx). +rewrite Zpos_digits2_pos, <- Z_of_nat_S_digits2_Pnat in Hplx. rewrite Zpower_nat_Z in H0. eapply Zlt_le_trans. apply H0. change 2%Z with (radix_val radix2). apply Zpower_le. @@ -210,13 +254,13 @@ unfold bits_of_binary_float, split_bits_of_binary_float. assert (Hf: (emin <= ex /\ Zdigits radix2 (Zpos mx) <= prec)%Z). destruct (andb_prop _ _ Hx) as (Hx', _). unfold canonic_mantissa in Hx'. -rewrite Z_of_nat_S_digits2_Pnat in Hx'. +rewrite Zpos_digits2_pos in Hx'. generalize (Zeq_bool_eq _ _ Hx'). unfold FLT_exp. -change (Fcalc_digits.radix2) with radix2. unfold emin. clear ; zify ; omega. -destruct (Zle_bool_spec (2^mw) (Zpos mx)) as [H|H] ; +case Zle_bool_spec ; intros H ; + [ apply -> Z.le_0_sub in H | apply -> Z.lt_sub_0 in H ] ; apply split_join_bits ; try now split. (* *) split. @@ -251,52 +295,45 @@ Qed. Theorem bits_of_binary_float_range: forall x, (0 <= bits_of_binary_float x < 2^(mw+ew+1))%Z. Proof. - intros. -Local Open Scope Z_scope. - assert (J: forall s m e, - 0 <= m < 2^mw -> 0 <= e < 2^ew -> - 0 <= join_bits s m e < 2^(mw+ew+1)). - { - intros. unfold join_bits. - set (se := (if s then 2 ^ ew else 0) + e). - assert (0 <= se < 2^(ew+1)). - { rewrite (Zpower_plus radix2) by omega. change (radix2^1) with 2. simpl. - unfold se. destruct s; omega. } - assert (0 <= se * 2^mw <= (2^(ew+1) - 1) * 2^mw). - { split. apply Zmult_gt_0_le_0_compat; omega. - apply Zmult_le_compat_r; omega. } - rewrite Z.mul_sub_distr_r in H2. - replace (mw + ew + 1) with ((ew + 1) + mw) by omega. - rewrite (Zpower_plus radix2) by omega. simpl. omega. - } - assert (D: forall p n, Z.of_nat (S (digits2_Pnat p)) <= n -> - 0 <= Z.pos p < 2^n). - { - intros. - generalize (digits2_Pnat_correct p). simpl. rewrite ! Zpower_nat_Z. intros [A B]. - split. zify; omega. eapply Zlt_le_trans. eassumption. - apply (Zpower_le radix2); auto. - } - destruct x; unfold bits_of_binary_float. -- apply J; omega. -- apply J; omega. -- destruct n as [pl pl_range]. apply Z.ltb_lt in pl_range. - apply J. apply D. unfold prec, Z_of_nat' in pl_range; omega. omega. -- unfold bounded in e0. apply Bool.andb_true_iff in e0; destruct e0 as [A B]. +unfold bits_of_binary_float. +intros [sx|sx|sx [pl pl_range]|sx mx ex H]. +- apply join_bits_range ; now split. +- apply join_bits_range. + now split. + clear -He_gt_0 ; omega. +- apply Z.ltb_lt in pl_range. + apply join_bits_range. + split. + easy. + apply (Zpower_gt_Zdigits radix2 _ (Zpos pl)). + apply Z.lt_succ_r. + now rewrite <- Zdigits2_Zdigits. + clear -He_gt_0 ; omega. +- unfold bounded in H. + apply Bool.andb_true_iff in H ; destruct H as [A B]. apply Z.leb_le in B. - unfold canonic_mantissa, FLT_exp in A. apply Zeq_bool_eq in A. - assert (G: Z.of_nat (S (digits2_Pnat m)) <= prec) by (zify; omega). - assert (M: emin <= e) by (unfold emin; zify; omega). - generalize (Zle_bool_spec (2^mw) (Z.pos m)); intro SPEC; inversion SPEC. - + apply J. - * split. omega. generalize (D _ _ G). unfold prec. - rewrite (Zpower_plus radix2) by omega. - change (radix2^1) with 2. simpl radix_val. omega. - * split. omega. unfold emin. replace (2^ew) with (2 * emax). omega. - symmetry. replace ew with (1 + (ew - 1)) by omega. - apply (Zpower_plus radix2); omega. - + apply J. zify; omega. omega. -Local Close Scope Z_scope. + unfold canonic_mantissa, FLT_exp in A. apply Zeq_bool_eq in A. + case Zle_bool_spec ; intros H. + + apply join_bits_range. + * split. + clear -H ; omega. + rewrite Zpos_digits2_pos in A. + cut (Zpos mx < 2 ^ prec)%Z. + unfold prec. + rewrite Zpower_plus by (clear -Hmw ; omega). + change (2^1)%Z with 2%Z. + clear ; omega. + apply (Zpower_gt_Zdigits radix2 _ (Zpos mx)). + clear -A ; zify ; omega. + * split. + unfold emin ; clear -A ; zify ; omega. + replace ew with ((ew - 1) + 1)%Z by ring. + rewrite Zpower_plus by (clear - Hew ; omega). + unfold emin, emax in *. + change (2^1)%Z with 2%Z. + clear -B ; omega. + + apply -> Z.lt_sub_0 in H. + apply join_bits_range ; now split. Qed. Definition binary_float_of_bits_aux x := @@ -360,7 +397,7 @@ case Zeq_bool_spec ; intros He2. case_eq (x mod 2 ^ mw)%Z; try easy. (* nan *) intros plx Eqplx. apply Z.ltb_lt. -rewrite Z_of_nat_S_digits2_Pnat. +rewrite Zpos_digits2_pos. assert (forall a b, a <= b -> a < b+1)%Z by (intros; omega). apply H. clear H. apply Zdigits_le_Zpower. simpl. rewrite <- Eqplx. edestruct Z_mod_lt; eauto. @@ -488,9 +525,8 @@ discriminate. clear -Hew ; omega. destruct (andb_prop _ _ Bx) as (H1, _). generalize (Zeq_bool_eq _ _ H1). -rewrite Z_of_nat_S_digits2_Pnat. +rewrite Zpos_digits2_pos. unfold FLT_exp, emin. -change Fcalc_digits.radix2 with radix2. generalize (Zdigits radix2 (Zpos mx)). clear. intros ; zify ; omega. @@ -500,9 +536,9 @@ simpl. apply f_equal. destruct (andb_prop _ _ Bx) as (H1, _). generalize (Zeq_bool_eq _ _ H1). -rewrite Z_of_nat_S_digits2_Pnat. +rewrite Zpos_digits2_pos. unfold FLT_exp, emin, prec. -change Fcalc_digits.radix2 with radix2. +apply -> Z.lt_sub_0 in Hm. generalize (Zdigits_le_Zpower radix2 _ (Zpos mx) Hm). generalize (Zdigits radix2 (Zpos mx)). clear. @@ -536,6 +572,7 @@ now rewrite He1 in Jx. intros px Hm Jx _. rewrite Zle_bool_false. now rewrite <- He1. +apply <- Z.lt_sub_0. now rewrite <- Hm. intros px Hm _ _. apply False_ind. @@ -556,7 +593,7 @@ intros p Hm Jx Cx. rewrite <- Hm. rewrite Zle_bool_true. now ring_simplify (mx + 2^mw - 2^mw)%Z (ex + emin - 1 - emin + 1)%Z. -now apply (Zplus_le_compat_r 0). +now ring_simplify. intros p Hm. apply False_ind. clear -Bm Hm ; zify ; omega. @@ -567,6 +604,8 @@ End Binary_Bits. (** Specialization for IEEE single precision operations *) Section B32_Bits. +Implicit Arguments B754_nan [[prec] [emax]]. + Definition binary32 := binary_float 24 128. Let Hprec : (0 < 24)%Z. @@ -577,12 +616,28 @@ Let Hprec_emax : (24 < 128)%Z. apply refl_equal. Qed. -Definition b32_opp := Bopp 24 128. -Definition b32_plus := Bplus _ _ Hprec Hprec_emax. -Definition b32_minus := Bminus _ _ Hprec Hprec_emax. -Definition b32_mult := Bmult _ _ Hprec Hprec_emax. -Definition b32_div := Bdiv _ _ Hprec Hprec_emax. -Definition b32_sqrt := Bsqrt _ _ Hprec Hprec_emax. +Definition default_nan_pl32 : bool * nan_pl 24 := + (false, exist _ (iter_nat 22 _ xO xH) (refl_equal true)). + +Definition unop_nan_pl32 (f : binary32) : bool * nan_pl 24 := + match f with + | B754_nan s pl => (s, pl) + | _ => default_nan_pl32 + end. + +Definition binop_nan_pl32 (f1 f2 : binary32) : bool * nan_pl 24 := + match f1, f2 with + | B754_nan s1 pl1, _ => (s1, pl1) + | _, B754_nan s2 pl2 => (s2, pl2) + | _, _ => default_nan_pl32 + end. + +Definition b32_opp := Bopp 24 128 pair. +Definition b32_plus := Bplus _ _ Hprec Hprec_emax binop_nan_pl32. +Definition b32_minus := Bminus _ _ Hprec Hprec_emax binop_nan_pl32. +Definition b32_mult := Bmult _ _ Hprec Hprec_emax binop_nan_pl32. +Definition b32_div := Bdiv _ _ Hprec Hprec_emax binop_nan_pl32. +Definition b32_sqrt := Bsqrt _ _ Hprec Hprec_emax unop_nan_pl32. Definition b32_of_bits : Z -> binary32 := binary_float_of_bits 23 8 (refl_equal _) (refl_equal _) (refl_equal _). Definition bits_of_b32 : binary32 -> Z := bits_of_binary_float 23 8. @@ -592,6 +647,8 @@ End B32_Bits. (** Specialization for IEEE double precision operations *) Section B64_Bits. +Implicit Arguments B754_nan [[prec] [emax]]. + Definition binary64 := binary_float 53 1024. Let Hprec : (0 < 53)%Z. @@ -602,12 +659,28 @@ Let Hprec_emax : (53 < 1024)%Z. apply refl_equal. Qed. -Definition b64_opp := Bopp 53 1024. -Definition b64_plus := Bplus _ _ Hprec Hprec_emax. -Definition b64_minus := Bminus _ _ Hprec Hprec_emax. -Definition b64_mult := Bmult _ _ Hprec Hprec_emax. -Definition b64_div := Bdiv _ _ Hprec Hprec_emax. -Definition b64_sqrt := Bsqrt _ _ Hprec Hprec_emax. +Definition default_nan_pl64 : bool * nan_pl 53 := + (false, exist _ (iter_nat 51 _ xO xH) (refl_equal true)). + +Definition unop_nan_pl64 (f : binary64) : bool * nan_pl 53 := + match f with + | B754_nan s pl => (s, pl) + | _ => default_nan_pl64 + end. + +Definition binop_nan_pl64 (pl1 pl2 : binary64) : bool * nan_pl 53 := + match pl1, pl2 with + | B754_nan s1 pl1, _ => (s1, pl1) + | _, B754_nan s2 pl2 => (s2, pl2) + | _, _ => default_nan_pl64 + end. + +Definition b64_opp := Bopp 53 1024 pair. +Definition b64_plus := Bplus _ _ Hprec Hprec_emax binop_nan_pl64. +Definition b64_minus := Bminus _ _ Hprec Hprec_emax binop_nan_pl64. +Definition b64_mult := Bmult _ _ Hprec Hprec_emax binop_nan_pl64. +Definition b64_div := Bdiv _ _ Hprec Hprec_emax binop_nan_pl64. +Definition b64_sqrt := Bsqrt _ _ Hprec Hprec_emax unop_nan_pl64. Definition b64_of_bits : Z -> binary64 := binary_float_of_bits 52 11 (refl_equal _) (refl_equal _) (refl_equal _). Definition bits_of_b64 : binary64 -> Z := bits_of_binary_float 52 11. diff --git a/flocq/Appli/Fappli_double_round.v b/flocq/Appli/Fappli_double_round.v new file mode 100644 index 00000000..f83abc47 --- /dev/null +++ b/flocq/Appli/Fappli_double_round.v @@ -0,0 +1,4554 @@ +(** * Conditions for innocuous double rounding. *) + +Require Import Fcore_Raux. +Require Import Fcore_defs. +Require Import Fcore_generic_fmt. +Require Import Fcalc_ops. +Require Import Fcore_ulp. + +Require Import Psatz. + +Open Scope R_scope. + +Section Double_round. + +Variable beta : radix. +Notation bpow e := (bpow beta e). +Notation ln_beta x := (ln_beta beta x). + +Definition double_round_eq fexp1 fexp2 choice1 choice2 x := + round beta fexp1 (Znearest choice1) (round beta fexp2 (Znearest choice2) x) + = round beta fexp1 (Znearest choice1) x. + +(** A little tactic to simplify terms of the form [bpow a * bpow b]. *) +Ltac bpow_simplify := + (* bpow ex * bpow ey ~~> bpow (ex + ey) *) + repeat + match goal with + | |- context [(Fcore_Raux.bpow _ _ * Fcore_Raux.bpow _ _)] => + rewrite <- bpow_plus + | |- context [(?X1 * Fcore_Raux.bpow _ _ * Fcore_Raux.bpow _ _)] => + rewrite (Rmult_assoc X1); rewrite <- bpow_plus + | |- context [(?X1 * (?X2 * Fcore_Raux.bpow _ _) * Fcore_Raux.bpow _ _)] => + rewrite <- (Rmult_assoc X1 X2); rewrite (Rmult_assoc (X1 * X2)); + rewrite <- bpow_plus + end; + (* ring_simplify arguments of bpow *) + repeat + match goal with + | |- context [(Fcore_Raux.bpow _ ?X)] => + progress ring_simplify X + end; + (* bpow 0 ~~> 1 *) + change (Fcore_Raux.bpow _ 0) with 1; + repeat + match goal with + | |- context [(_ * 1)] => + rewrite Rmult_1_r + end. + +Definition midp (fexp : Z -> Z) (x : R) := + round beta fexp Zfloor x + / 2 * ulp beta fexp x. + +Definition midp' (fexp : Z -> Z) (x : R) := + round beta fexp Zceil x - / 2 * ulp beta fexp x. + +Lemma double_round_lt_mid_further_place' : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + forall x, + 0 < x -> + (fexp2 (ln_beta x) <= fexp1 (ln_beta x) - 1)%Z -> + x < bpow (ln_beta x) - / 2 * ulp beta fexp2 x -> + x < midp fexp1 x - / 2 * ulp beta fexp2 x -> + double_round_eq fexp1 fexp2 choice1 choice2 x. +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 x Px Hf2f1 Hx1. +unfold double_round_eq. +set (x' := round beta fexp1 Zfloor x). +intro Hx2'. +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 (Pxx' : 0 <= x - x'). +{ apply Rle_0_minus. + apply round_DN_pt. + exact Vfexp1. } +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 _ _)). + replace (/ 2 * _) with (/ 2 * bpow (fexp2 (ln_beta x)) + + (/ 2 * (bpow (fexp1 (ln_beta x)) + - bpow (fexp2 (ln_beta x))))) by ring. + apply Rplus_le_lt_compat. + - exact Hr1. + - now rewrite Rabs_right; [|now apply Rle_ge]; apply Hx2. } +destruct (Req_dec x'' 0) as [Zx''|Nzx'']. +- (* x'' = 0 *) + rewrite Zx'' in Hr1 |- *. + rewrite round_0; [|now apply valid_rnd_N]. + unfold round, F2R, scaled_mantissa, canonic_exp; simpl. + rewrite (Znearest_imp _ _ 0); [now simpl; rewrite Rmult_0_l|]. + 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]. + rewrite <- Rabs_mult; rewrite Rmult_minus_distr_r. + rewrite Rmult_0_l. + bpow_simplify. + rewrite Rabs_minus_sym. + apply (Rle_lt_trans _ _ _ Hr1). + apply Rmult_lt_compat_l; [lra|]. + apply bpow_lt. + omega. +- (* x'' <> 0 *) + assert (Lx'' : ln_beta x'' = ln_beta x :> Z). + { apply Zle_antisym. + - apply ln_beta_le_bpow; [exact Nzx''|]. + replace x'' with (x'' - x + x) by ring. + apply (Rle_lt_trans _ _ _ (Rabs_triang _ _)). + replace (bpow _) with (/ 2 * bpow (fexp2 (ln_beta x)) + + (bpow (ln_beta x) + - / 2 * bpow (fexp2 (ln_beta x)))) by ring. + apply Rplus_le_lt_compat; [exact Hr1|]. + now rewrite Rabs_right; [|apply Rle_ge; apply Rlt_le]. + - unfold x'' in Nzx'' |- *. + now apply ln_beta_round_ge; [|apply valid_rnd_N|]. } + unfold round, F2R, scaled_mantissa, canonic_exp; simpl. + rewrite Lx''. + rewrite (Znearest_imp _ _ (Zfloor (scaled_mantissa beta fexp1 x))). + + rewrite (Znearest_imp _ _ (Zfloor (scaled_mantissa beta fexp1 x))); + [reflexivity|]. + 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]. + rewrite <- Rabs_mult. + rewrite Rmult_minus_distr_r. + fold x'. + bpow_simplify. + rewrite Rabs_right; [|now apply Rle_ge]. + 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. + + 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]. + rewrite <- Rabs_mult. + rewrite Rmult_minus_distr_r. + fold x'. + now bpow_simplify. +Qed. + +Lemma double_round_lt_mid_further_place : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + forall x, + 0 < x -> + (fexp2 (ln_beta x) <= fexp1 (ln_beta x) - 1)%Z -> + (fexp1 (ln_beta x) <= ln_beta x)%Z -> + x < midp fexp1 x - / 2 * ulp beta fexp2 x -> + double_round_eq fexp1 fexp2 choice1 choice2 x. +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 x Px Hf2f1 Hf1. +intro Hx2'. +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. } +revert Hx2. +unfold double_round_eq. +set (x' := round beta fexp1 Zfloor x). +intro Hx2. +assert (Pxx' : 0 <= x - x'). +{ apply Rle_0_minus. + 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']. +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. + 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. + rewrite <- (Rinv_r 2) at 3; [|lra]. + rewrite Rmult_comm; apply Rmult_le_compat_r; [lra|]. + apply Rle_trans with 1; [|lra]. + change 1 with (bpow 0); apply bpow_le. + omega. +- (* x' <> 0 *) + assert (Px' : 0 < x'). + { assert (0 <= x'); [|lra]. + unfold x'. + apply (Rmult_le_reg_r (bpow (- fexp1 (ln_beta x)))); + [now apply bpow_gt_0|]. + rewrite Rmult_0_l. + unfold round, F2R, canonic_exp; simpl; bpow_simplify. + change 0 with (Z2R 0); apply Z2R_le. + apply Zfloor_lub. + rewrite <- (Rabs_right x); [|now apply Rle_ge; apply Rlt_le]. + rewrite scaled_mantissa_abs. + apply Rabs_pos. } + assert (Hx' : x' <= bpow (ln_beta x) - ulp beta fexp1 x). + { apply (Rplus_le_reg_r (ulp beta fexp1 x)); ring_simplify. + rewrite <- ulp_DN. + - change (round _ _ _ _) with x'. + apply succ_le_bpow. + + exact Px'. + + change x' with (round beta fexp1 Zfloor x). + now apply generic_format_round; [|apply valid_rnd_DN]. + + apply Rle_lt_trans with x. + * now apply round_DN_pt. + * rewrite <- (Rabs_right x) at 1; [|now apply Rle_ge; apply Rlt_le]. + apply bpow_ln_beta_gt. + - 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]. + rewrite <- (Rmult_1_l (ulp _ _ _)) at 2. + apply Rmult_le_compat_r; [|lra]. + apply bpow_ge_0. +Qed. + +Lemma double_round_lt_mid_same_place : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> + forall (choice1 choice2 : Z -> bool), + forall x, + 0 < x -> + (fexp2 (ln_beta x) = fexp1 (ln_beta x))%Z -> + x < midp fexp1 x -> + double_round_eq fexp1 fexp2 choice1 choice2 x. +Proof. +intros fexp1 fexp2 Vfexp1 choice1 choice2 x Px Hf2f1. +intro Hx'. +assert (Hx : x - round beta fexp1 Zfloor x < / 2 * ulp beta fexp1 x). +{ now apply (Rplus_lt_reg_r (round beta fexp1 Zfloor x)); ring_simplify. } +revert Hx. +unfold double_round_eq. +set (x' := round beta fexp1 Zfloor x). +intro Hx. +assert (Pxx' : 0 <= x - x'). +{ apply Rle_0_minus. + apply round_DN_pt. + exact Vfexp1. } +assert (H : Rabs (x * bpow (- fexp1 (ln_beta x)) - + Z2R (Zfloor (x * bpow (- fexp1 (ln_beta x))))) < / 2). +{ apply (Rmult_lt_reg_r (bpow (fexp1 (ln_beta x)))); [now apply bpow_gt_0|]. + unfold scaled_mantissa, canonic_exp in Hx. + rewrite <- (Rabs_right (bpow (fexp1 _))) at 1; + [|now apply Rle_ge; apply bpow_ge_0]. + rewrite <- Rabs_mult. + rewrite Rmult_minus_distr_r. + bpow_simplify. + apply Rabs_lt. + change (Z2R _ * _) with x'. + split. + - apply Rlt_le_trans with 0; [|exact Pxx']. + rewrite <- Ropp_0. + apply Ropp_lt_contravar. + rewrite <- (Rmult_0_r (/ 2)). + apply Rmult_lt_compat_l; [lra|]. + apply bpow_gt_0. + - exact Hx. } +unfold round at 2. +unfold F2R, scaled_mantissa, canonic_exp; simpl. +rewrite Hf2f1. +rewrite (Znearest_imp _ _ (Zfloor (scaled_mantissa beta fexp1 x))). +- rewrite round_generic. + + unfold round, F2R, scaled_mantissa, canonic_exp; simpl. + now rewrite (Znearest_imp _ _ (Zfloor (x * bpow (- fexp1 (ln_beta x))))). + + now apply valid_rnd_N. + + fold (canonic_exp beta fexp1 x). + change (Z2R _ * bpow _) with (round beta fexp1 Zfloor x). + apply generic_format_round. + exact Vfexp1. + now apply valid_rnd_DN. +- now unfold scaled_mantissa, canonic_exp. +Qed. + +Lemma double_round_lt_mid : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + forall x, + 0 < x -> + (fexp2 (ln_beta x) <= fexp1 (ln_beta x))%Z -> + (fexp1 (ln_beta x) <= ln_beta x)%Z -> + x < midp fexp1 x -> + ((fexp2 (ln_beta x) <= fexp1 (ln_beta x) - 1)%Z -> + x < midp fexp1 x - / 2 * ulp beta fexp2 x) -> + double_round_eq fexp1 fexp2 choice1 choice2 x. +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 x Px Hf2f1 Hf1 Hx Hx'. +destruct (Zle_or_lt (fexp1 (ln_beta x)) (fexp2 (ln_beta x))) as [Hf2'|Hf2']. +- (* fexp1 (ln_beta x) <= fexp2 (ln_beta x) *) + assert (Hf2'' : (fexp2 (ln_beta x) = fexp1 (ln_beta x) :> Z)%Z); [omega|]. + now apply double_round_lt_mid_same_place. +- (* fexp2 (ln_beta x) < fexp1 (ln_beta x) *) + assert (Hf2'' : (fexp2 (ln_beta x) <= fexp1 (ln_beta x) - 1)%Z); [omega|]. + generalize (Hx' Hf2''); intro Hx''. + now apply double_round_lt_mid_further_place. +Qed. + +Lemma double_round_gt_mid_further_place' : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + forall x, + 0 < x -> + (fexp2 (ln_beta x) <= fexp1 (ln_beta x) - 1)%Z -> + round beta fexp2 (Znearest choice2) x < bpow (ln_beta x) -> + midp' fexp1 x + / 2 * ulp beta fexp2 x < x -> + double_round_eq fexp1 fexp2 choice1 choice2 x. +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 x Px Hf2f1. +intros Hx1 Hx2'. +assert (Hx2 : round beta fexp1 Zceil x - x + < / 2 * (ulp beta fexp1 x - ulp beta fexp2 x)). +{ apply (Rplus_lt_reg_r (- / 2 * ulp beta fexp1 x + x + + / 2 * ulp beta fexp2 x)); ring_simplify. + now unfold midp' in Hx2'. } +revert Hx1 Hx2. +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 (Px'x : 0 <= x' - x). +{ apply Rle_0_minus. + apply round_UP_pt. + exact Vfexp1. } +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 _ _)). + replace (/ 2 * _) with (/ 2 * bpow (fexp2 (ln_beta x)) + + (/ 2 * (bpow (fexp1 (ln_beta x)) + - bpow (fexp2 (ln_beta x))))) by ring. + apply Rplus_le_lt_compat. + - exact Hr1. + - rewrite Rabs_minus_sym. + now rewrite Rabs_right; [|now apply Rle_ge]; apply Hx2. } +destruct (Req_dec x'' 0) as [Zx''|Nzx'']. +- (* x'' = 0 *) + rewrite Zx'' in Hr1 |- *. + rewrite round_0; [|now apply valid_rnd_N]. + unfold round, F2R, scaled_mantissa, canonic_exp; simpl. + rewrite (Znearest_imp _ _ 0); [now simpl; rewrite Rmult_0_l|]. + 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]. + rewrite <- Rabs_mult; rewrite Rmult_minus_distr_r. + rewrite Rmult_0_l. + bpow_simplify. + rewrite Rabs_minus_sym. + apply (Rle_lt_trans _ _ _ Hr1). + apply Rmult_lt_compat_l; [lra|]. + apply bpow_lt. + omega. +- (* x'' <> 0 *) + assert (Lx'' : ln_beta x'' = ln_beta x :> Z). + { apply Zle_antisym. + - apply ln_beta_le_bpow; [exact Nzx''|]. + rewrite Rabs_right; [exact Hx1|apply Rle_ge]. + apply round_ge_generic. + + exact Vfexp2. + + now apply valid_rnd_N. + + apply generic_format_0. + + now apply Rlt_le. + - unfold x'' in Nzx'' |- *. + now apply ln_beta_round_ge; [|apply valid_rnd_N|]. } + unfold round, F2R, scaled_mantissa, canonic_exp; simpl. + rewrite Lx''. + rewrite (Znearest_imp _ _ (Zceil (scaled_mantissa beta fexp1 x))). + + rewrite (Znearest_imp _ _ (Zceil (scaled_mantissa beta fexp1 x))); + [reflexivity|]. + 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]. + rewrite <- Rabs_mult. + rewrite Rmult_minus_distr_r. + fold x'. + bpow_simplify. + rewrite Rabs_minus_sym. + rewrite Rabs_right; [|now apply Rle_ge]. + 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. + + 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]. + rewrite <- Rabs_mult. + rewrite Rmult_minus_distr_r. + fold x'. + now bpow_simplify. +Qed. + +Lemma double_round_gt_mid_further_place : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + forall x, + 0 < x -> + (fexp2 (ln_beta x) <= fexp1 (ln_beta x) - 1)%Z -> + (fexp1 (ln_beta x) <= ln_beta x)%Z -> + midp' fexp1 x + / 2 * ulp beta fexp2 x < x -> + double_round_eq fexp1 fexp2 choice1 choice2 x. +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 x Px Hf2f1 Hf1 Hx2'. +assert (Hx2 : round beta fexp1 Zceil x - x + < / 2 * (ulp beta fexp1 x - ulp beta fexp2 x)). +{ apply (Rplus_lt_reg_r (- / 2 * ulp beta fexp1 x + x + + / 2 * ulp beta fexp2 x)); ring_simplify. + now unfold midp' in Hx2'. } +revert Hx2. +unfold double_round_eq. +set (x' := round beta fexp1 Zfloor x). +intro Hx2. +set (x'' := round beta fexp2 (Znearest choice2) x). +destruct (Rlt_or_le x'' (bpow (ln_beta x))) as [Hx''|Hx'']; + [now apply double_round_gt_mid_further_place'|]. +(* bpow (ln_beta x) <= x'' *) +assert (Hx''pow : x'' = bpow (ln_beta x)). +{ assert (H'x'' : x'' < bpow (ln_beta x) + / 2 * ulp beta fexp2 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. + exact Vfexp2. + - apply Rplus_lt_compat_r. + rewrite <- Rabs_right at 1; [|now apply Rle_ge; apply Rlt_le]. + apply bpow_ln_beta_gt. } + apply Rle_antisym; [|exact Hx'']. + unfold x'', round, F2R, scaled_mantissa, canonic_exp; simpl. + apply (Rmult_le_reg_r (bpow (- fexp2 (ln_beta x)))); [now apply bpow_gt_0|]. + bpow_simplify. + rewrite <- (Z2R_Zpower _ (_ - _)); [|omega]. + apply Z2R_le. + apply Zlt_succ_le; unfold Z.succ. + apply lt_Z2R. + rewrite Z2R_plus; rewrite Z2R_Zpower; [|omega]. + apply (Rmult_lt_reg_r (bpow (fexp2 (ln_beta x)))); [now apply bpow_gt_0|]. + rewrite Rmult_plus_distr_r; rewrite Rmult_1_l. + bpow_simplify. + apply (Rlt_le_trans _ _ _ H'x''). + apply Rplus_le_compat_l. + rewrite <- (Rmult_1_l (Fcore_Raux.bpow _ _)). + 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. + exact Vfexp2. + - apply Rmult_lt_compat_l; [lra|]. + unfold ulp, 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). +{ rewrite Hx''pow. + rewrite ln_beta_bpow. + assert (fexp1 (ln_beta x + 1) <= ln_beta x)%Z; [|omega]. + destruct (Zle_or_lt (ln_beta x) (fexp1 (ln_beta x))) as [Hle|Hlt]; + [|now apply Vfexp1]. + assert (H : (ln_beta x = fexp1 (ln_beta x) :> Z)%Z); + [now apply Zle_antisym|]. + rewrite H. + now apply Vfexp1. } +rewrite (Znearest_imp _ _ (beta ^ (ln_beta x - fexp1 (ln_beta x'')))%Z). +- rewrite (Znearest_imp _ _ (beta ^ (ln_beta x - fexp1 (ln_beta x)))%Z). + + rewrite Z2R_Zpower; [|exact Hf]. + rewrite Z2R_Zpower; [|omega]. + now bpow_simplify. + + rewrite Z2R_Zpower; [|omega]. + 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]. + rewrite <- Rabs_mult. + rewrite Rmult_minus_distr_r. + bpow_simplify. + 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|]. + rewrite <- (Rabs_right (bpow (fexp1 _))) at 1; + [|now apply Rle_ge; apply bpow_ge_0]. + rewrite <- Rabs_mult. + rewrite Rmult_minus_distr_r. + bpow_simplify. + rewrite Rminus_diag_eq; [|exact Hx''pow]. + rewrite Rabs_R0. + rewrite <- (Rmult_0_r (/ 2)). + apply Rmult_lt_compat_l; [lra|apply bpow_gt_0]. +Qed. + +Lemma double_round_gt_mid_same_place : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> + forall (choice1 choice2 : Z -> bool), + forall x, + 0 < x -> + (fexp2 (ln_beta x) = fexp1 (ln_beta x))%Z -> + midp' fexp1 x < x -> + double_round_eq fexp1 fexp2 choice1 choice2 x. +Proof. +intros fexp1 fexp2 Vfexp1 choice1 choice2 x Px Hf2f1 Hx'. +assert (Hx : round beta fexp1 Zceil x - x < / 2 * ulp beta fexp1 x). +{ apply (Rplus_lt_reg_r (- / 2 * ulp beta fexp1 x + x)); ring_simplify. + now unfold midp' in Hx'. } +assert (H : Rabs (Z2R (Zceil (x * bpow (- fexp1 (ln_beta x)))) + - x * bpow (- fexp1 (ln_beta x))) < / 2). +{ apply (Rmult_lt_reg_r (bpow (fexp1 (ln_beta x)))); [now apply bpow_gt_0|]. + unfold scaled_mantissa, canonic_exp in Hx. + rewrite <- (Rabs_right (bpow (fexp1 _))) at 1; + [|now apply Rle_ge; apply bpow_ge_0]. + rewrite <- Rabs_mult. + rewrite Rmult_minus_distr_r. + bpow_simplify. + apply Rabs_lt. + split. + - apply Rlt_le_trans with 0. + + rewrite <- Ropp_0; apply Ropp_lt_contravar. + rewrite <- (Rmult_0_r (/ 2)). + apply Rmult_lt_compat_l; [lra|]. + apply bpow_gt_0. + + apply Rle_0_minus. + apply round_UP_pt. + exact Vfexp1. + - exact Hx. } +unfold double_round_eq, round at 2. +unfold F2R, scaled_mantissa, canonic_exp; simpl. +rewrite Hf2f1. +rewrite (Znearest_imp _ _ (Zceil (scaled_mantissa beta fexp1 x))). +- rewrite round_generic. + + unfold round, F2R, scaled_mantissa, canonic_exp; simpl. + now rewrite (Znearest_imp _ _ (Zceil (x * bpow (- fexp1 (ln_beta x))))); + [|rewrite Rabs_minus_sym]. + + now apply valid_rnd_N. + + fold (canonic_exp beta fexp1 x). + change (Z2R _ * bpow _) with (round beta fexp1 Zceil x). + apply generic_format_round. + exact Vfexp1. + now apply valid_rnd_UP. +- now rewrite Rabs_minus_sym. +Qed. + +Lemma double_round_gt_mid : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + forall x, + 0 < x -> + (fexp2 (ln_beta x) <= fexp1 (ln_beta x))%Z -> + (fexp1 (ln_beta x) <= ln_beta x)%Z -> + midp' fexp1 x < x -> + ((fexp2 (ln_beta x) <= fexp1 (ln_beta x) - 1)%Z -> + midp' fexp1 x + / 2 * ulp beta fexp2 x < x) -> + double_round_eq fexp1 fexp2 choice1 choice2 x. +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 x Px Hf2f1 Hf1 Hx Hx'. +destruct (Zle_or_lt (fexp1 (ln_beta x)) (fexp2 (ln_beta x))) as [Hf2'|Hf2']. +- (* fexp1 (ln_beta x) <= fexp2 (ln_beta x) *) + assert (Hf2'' : (fexp2 (ln_beta x) = fexp1 (ln_beta x) :> Z)%Z); [omega|]. + now apply double_round_gt_mid_same_place. +- (* fexp2 (ln_beta x) < fexp1 (ln_beta x) *) + assert (Hf2'' : (fexp2 (ln_beta x) <= fexp1 (ln_beta x) - 1)%Z); [omega|]. + generalize (Hx' Hf2''); intro Hx''. + now apply double_round_gt_mid_further_place. +Qed. + +Section Double_round_mult. + +Lemma ln_beta_mult_disj : + forall x y, + x <> 0 -> y <> 0 -> + ((ln_beta (x * y) = (ln_beta x + ln_beta y - 1)%Z :> Z) + \/ (ln_beta (x * y) = (ln_beta x + ln_beta y)%Z :> Z)). +Proof. +intros x y Zx Zy. +destruct (ln_beta_mult beta x y Zx Zy). +omega. +Qed. + +Definition double_round_mult_hyp fexp1 fexp2 := + (forall ex ey, (fexp2 (ex + ey) <= fexp1 ex + fexp1 ey)%Z) + /\ (forall ex ey, (fexp2 (ex + ey - 1) <= fexp1 ex + fexp1 ey)%Z). + +Lemma double_round_mult_aux : + forall (fexp1 fexp2 : Z -> Z), + double_round_mult_hyp fexp1 fexp2 -> + forall x y, + generic_format beta fexp1 x -> generic_format beta fexp1 y -> + generic_format beta fexp2 (x * y). +Proof. +intros fexp1 fexp2 Hfexp x y Fx Fy. +destruct (Req_dec x 0) as [Zx|Zx]. +- (* x = 0 *) + rewrite Zx. + rewrite Rmult_0_l. + now apply generic_format_0. +- (* x <> 0 *) + destruct (Req_dec y 0) as [Zy|Zy]. + + (* y = 0 *) + rewrite Zy. + rewrite Rmult_0_r. + now apply generic_format_0. + + (* y <> 0 *) + revert Fx Fy. + unfold generic_format. + unfold canonic_exp. + set (mx := Ztrunc (scaled_mantissa beta fexp1 x)). + set (my := Ztrunc (scaled_mantissa beta fexp1 y)). + unfold F2R; simpl. + intros Fx Fy. + set (fxy := Float beta (mx * my) (fexp1 (ln_beta x) + fexp1 (ln_beta y))). + assert (Hxy : x * y = F2R fxy). + { unfold fxy, F2R; simpl. + rewrite bpow_plus. + rewrite Z2R_mult. + rewrite Fx, Fy at 1. + ring. } + apply generic_format_F2R' with (f := fxy); [now rewrite Hxy|]. + intros _. + unfold canonic_exp, fxy; simpl. + destruct Hfexp as (Hfexp1, Hfexp2). + now destruct (ln_beta_mult_disj x y Zx Zy) as [Lxy|Lxy]; rewrite Lxy. +Qed. + +Variable rnd : R -> Z. +Context { valid_rnd : Valid_rnd rnd }. + +Theorem double_round_mult : + forall (fexp1 fexp2 : Z -> Z), + double_round_mult_hyp fexp1 fexp2 -> + forall x y, + generic_format beta fexp1 x -> generic_format beta fexp1 y -> + round beta fexp1 rnd (round beta fexp2 rnd (x * y)) + = round beta fexp1 rnd (x * y). +Proof. +intros fexp1 fexp2 Hfexp x y Fx Fy. +assert (Hxy : round beta fexp2 rnd (x * y) = x * y). +{ apply round_generic; [assumption|]. + now apply (double_round_mult_aux fexp1 fexp2). } +now rewrite Hxy at 1. +Qed. + +Section Double_round_mult_FLX. + +Require Import Fcore_FLX. + +Variable prec : Z. +Variable prec' : Z. + +Context { prec_gt_0_ : Prec_gt_0 prec }. +Context { prec_gt_0_' : Prec_gt_0 prec' }. + +Theorem double_round_mult_FLX : + (2 * prec <= prec')%Z -> + forall x y, + FLX_format beta prec x -> FLX_format beta prec y -> + round beta (FLX_exp prec) rnd (round beta (FLX_exp prec') rnd (x * y)) + = round beta (FLX_exp prec) rnd (x * y). +Proof. +intros Hprec x y Fx Fy. +apply double_round_mult; + [|now apply generic_format_FLX|now apply generic_format_FLX]. +unfold double_round_mult_hyp; split; intros ex ey; unfold FLX_exp; +omega. +Qed. + +End Double_round_mult_FLX. + +Section Double_round_mult_FLT. + +Require Import Fcore_FLX. +Require Import Fcore_FLT. + +Variable emin prec : Z. +Variable emin' prec' : Z. + +Context { prec_gt_0_ : Prec_gt_0 prec }. +Context { prec_gt_0_' : Prec_gt_0 prec' }. + +Theorem double_round_mult_FLT : + (emin' <= 2 * emin)%Z -> (2 * prec <= prec')%Z -> + forall x y, + FLT_format beta emin prec x -> FLT_format beta emin prec y -> + round beta (FLT_exp emin prec) rnd + (round beta (FLT_exp emin' prec') rnd (x * y)) + = round beta (FLT_exp emin prec) rnd (x * y). +Proof. +intros Hemin Hprec x y Fx Fy. +apply double_round_mult; + [|now apply generic_format_FLT|now apply generic_format_FLT]. +unfold double_round_mult_hyp; split; intros ex ey; +unfold FLT_exp; +generalize (Zmax_spec (ex + ey - prec') emin'); +generalize (Zmax_spec (ex + ey - 1 - prec') emin'); +generalize (Zmax_spec (ex - prec) emin); +generalize (Zmax_spec (ey - prec) emin); +omega. +Qed. + +End Double_round_mult_FLT. + +Section Double_round_mult_FTZ. + +Require Import Fcore_FTZ. + +Variable emin prec : Z. +Variable emin' prec' : Z. + +Context { prec_gt_0_ : Prec_gt_0 prec }. +Context { prec_gt_0_' : Prec_gt_0 prec' }. + +Theorem double_round_mult_FTZ : + (emin' + prec' <= 2 * emin + prec)%Z -> + (2 * prec <= prec')%Z -> + forall x y, + FTZ_format beta emin prec x -> FTZ_format beta emin prec y -> + round beta (FTZ_exp emin prec) rnd + (round beta (FTZ_exp emin' prec') rnd (x * y)) + = round beta (FTZ_exp emin prec) rnd (x * y). +Proof. +intros Hemin Hprec x y Fx Fy. +apply double_round_mult; + [|now apply generic_format_FTZ|now apply generic_format_FTZ]. +unfold double_round_mult_hyp; split; intros ex ey; +unfold FTZ_exp; +unfold Prec_gt_0 in *; +destruct (Z.ltb_spec (ex + ey - prec') emin'); +destruct (Z.ltb_spec (ex - prec) emin); +destruct (Z.ltb_spec (ey - prec) emin); +destruct (Z.ltb_spec (ex + ey - 1 - prec') emin'); +omega. +Qed. + +End Double_round_mult_FTZ. + +End Double_round_mult. + +Section Double_round_plus. + +Lemma ln_beta_plus_disj : + forall x y, + 0 < y -> y <= x -> + ((ln_beta (x + y) = ln_beta x :> Z) + \/ (ln_beta (x + y) = (ln_beta x + 1)%Z :> Z)). +Proof. +intros x y Py Hxy. +destruct (ln_beta_plus beta x y Py Hxy). +omega. +Qed. + +Lemma ln_beta_plus_separated : + forall fexp : Z -> Z, + forall x y, + 0 < x -> 0 <= y -> + generic_format beta fexp x -> + (ln_beta y <= fexp (ln_beta x))%Z -> + (ln_beta (x + y) = ln_beta x :> Z). +Proof. +intros fexp x y Px Nny Fx Hsep. +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|]. + split; [assumption|]. + unfold ulp, 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]. + + now apply bpow_le. +Qed. + +Lemma ln_beta_minus_disj : + forall x y, + 0 < x -> 0 < y -> + (ln_beta y <= ln_beta x - 2)%Z -> + ((ln_beta (x - y) = ln_beta x :> Z) + \/ (ln_beta (x - y) = (ln_beta x - 1)%Z :> Z)). +Proof. +intros x y Px Py Hln. +assert (Hxy : y < x); [now apply (ln_beta_lt_pos beta); [| |omega]|]. +generalize (ln_beta_minus beta x y Py Hxy); intro Hln2. +generalize (ln_beta_minus_lb beta x y Px Py Hln); intro Hln3. +omega. +Qed. + +Lemma ln_beta_minus_separated : + forall fexp : Z -> Z, Valid_exp fexp -> + forall x y, + 0 < x -> 0 < y -> y < x -> + bpow (ln_beta x - 1) < x -> + generic_format beta fexp x -> (ln_beta y <= fexp (ln_beta x))%Z -> + (ln_beta (x - y) = ln_beta x :> Z). +Proof. +intros fexp Vfexp x y Px Py Yltx Xgtpow Fx Ly. +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. + 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)). + - apply bpow_ln_beta_gt. + - now apply bpow_le. } + apply (Rplus_le_reg_r y); ring_simplify. + 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|]]. + 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]. + now apply ln_beta_generic_gt; [|now apply Rgt_not_eq|]. +- rewrite Rabs_right. + + apply Rlt_trans with x. + * rewrite <- (Rplus_0_r x) at 2. + apply Rplus_lt_compat_l. + rewrite <- Ropp_0. + now apply Ropp_lt_contravar. + * apply Rabs_lt_inv. + apply bpow_ln_beta_gt. + + lra. +Qed. + +Definition double_round_plus_hyp fexp1 fexp2 := + (forall ex ey, (fexp1 (ex + 1) - 1 <= ey)%Z -> (fexp2 ex <= fexp1 ey)%Z) + /\ (forall ex ey, (fexp1 (ex - 1) + 1 <= ey)%Z -> (fexp2 ex <= fexp1 ey)%Z) + /\ (forall ex ey, (fexp1 ex - 1 <= ey)%Z -> (fexp2 ex <= fexp1 ey)%Z) + /\ (forall ex ey, (ex - 1 <= ey)%Z -> (fexp2 ex <= fexp1 ey)%Z). + +Lemma double_round_plus_aux0_aux_aux : + forall (fexp1 fexp2 : Z -> Z), + forall x y, + (fexp1 (ln_beta x) <= fexp1 (ln_beta y))%Z -> + (fexp2 (ln_beta (x + y))%Z <= fexp1 (ln_beta x))%Z -> + (fexp2 (ln_beta (x + y))%Z <= fexp1 (ln_beta y))%Z -> + generic_format beta fexp1 x -> generic_format beta fexp1 y -> + generic_format beta fexp2 (x + y). +Proof. +intros fexp1 fexp2 x y Oxy Hlnx Hlny Fx Fy. +destruct (Req_dec x 0) as [Zx|Nzx]. +- (* x = 0 *) + rewrite Zx, Rplus_0_l in Hlny |- *. + now apply (generic_inclusion_ln_beta beta fexp1). +- (* x <> 0 *) + destruct (Req_dec y 0) as [Zy|Nzy]. + + (* y = 0 *) + rewrite Zy, Rplus_0_r in Hlnx |- *. + now apply (generic_inclusion_ln_beta beta fexp1). + + (* y <> 0 *) + revert Fx Fy. + unfold generic_format at -3, canonic_exp, F2R; simpl. + set (mx := Ztrunc (scaled_mantissa beta fexp1 x)). + set (my := Ztrunc (scaled_mantissa beta fexp1 y)). + intros Fx Fy. + set (fxy := Float beta (mx + my * (beta ^ (fexp1 (ln_beta y) + - fexp1 (ln_beta x)))) + (fexp1 (ln_beta x))). + assert (Hxy : x + y = F2R fxy). + { unfold fxy, F2R; simpl. + rewrite Z2R_plus. + rewrite Rmult_plus_distr_r. + rewrite <- Fx. + rewrite Z2R_mult. + rewrite Z2R_Zpower; [|omega]. + bpow_simplify. + now rewrite <- Fy. } + apply generic_format_F2R' with (f := fxy); [now rewrite Hxy|]. + intros _. + now unfold canonic_exp, fxy; simpl. +Qed. + +Lemma double_round_plus_aux0_aux : + forall (fexp1 fexp2 : Z -> Z), + forall x y, + (fexp2 (ln_beta (x + y))%Z <= fexp1 (ln_beta x))%Z -> + (fexp2 (ln_beta (x + y))%Z <= fexp1 (ln_beta y))%Z -> + generic_format beta fexp1 x -> generic_format beta fexp1 y -> + generic_format beta fexp2 (x + y). +Proof. +intros fexp1 fexp2 x y Hlnx Hlny Fx Fy. +destruct (Z.le_gt_cases (fexp1 (ln_beta x)) (fexp1 (ln_beta y))) as [Hle|Hgt]. +- now apply (double_round_plus_aux0_aux_aux fexp1). +- rewrite Rplus_comm in Hlnx, Hlny |- *. + now apply (double_round_plus_aux0_aux_aux fexp1); [omega| | | |]. +Qed. + +(* fexp1 (ln_beta x) - 1 <= ln_beta y : + * addition is exact in the largest precision (fexp2). *) +Lemma double_round_plus_aux0 : + forall (fexp1 fexp2 : Z -> Z), Valid_exp fexp1 -> + double_round_plus_hyp fexp1 fexp2 -> + forall x y, + (0 < x)%R -> (0 < y)%R -> (y <= x)%R -> + (fexp1 (ln_beta x) - 1 <= ln_beta y)%Z -> + generic_format beta fexp1 x -> generic_format beta fexp1 y -> + generic_format beta fexp2 (x + y). +Proof. +intros fexp1 fexp2 Vfexp1 Hexp x y Px Py Hyx Hln Fx Fy. +assert (Nny : (0 <= y)%R); [now apply Rlt_le|]. +destruct Hexp as (_,(Hexp2,(Hexp3,Hexp4))). +destruct (Z.le_gt_cases (ln_beta y) (fexp1 (ln_beta x))) as [Hle|Hgt]. +- (* ln_beta y <= fexp1 (ln_beta x) *) + assert (Lxy : ln_beta (x + y) = ln_beta x :> Z); + [now apply (ln_beta_plus_separated fexp1)|]. + apply (double_round_plus_aux0_aux fexp1); + [| |assumption|assumption]; rewrite Lxy. + + now apply Hexp4; omega. + + now apply Hexp3; omega. +- (* fexp1 (ln_beta x) < ln_beta y *) + apply (double_round_plus_aux0_aux fexp1); [| |assumption|assumption]. + destruct (ln_beta_plus_disj x y Py Hyx) as [Lxy|Lxy]; rewrite Lxy. + + now apply Hexp4; omega. + + apply Hexp2; apply (ln_beta_le beta y x Py) in Hyx. + replace (_ - _)%Z with (ln_beta x : Z) by ring. + omega. + + destruct (ln_beta_plus_disj x y Py Hyx) as [Lxy|Lxy]; rewrite Lxy. + * now apply Hexp3; omega. + * apply Hexp2. + replace (_ - _)%Z with (ln_beta x : Z) by ring. + omega. +Qed. + +Lemma double_round_plus_aux1_aux : + forall k, (0 < k)%Z -> + forall (fexp : Z -> Z), + forall x y, + 0 < x -> 0 < y -> + (ln_beta y <= fexp (ln_beta x) - k)%Z -> + (ln_beta (x + y) = ln_beta x :> Z) -> + generic_format beta fexp x -> + 0 < (x + y) - round beta fexp Zfloor (x + y) < bpow (fexp (ln_beta x) - k). +Proof. +assert (Hbeta : (2 <= beta)%Z). +{ destruct beta as (beta_val,beta_prop). + now apply Zle_bool_imp_le. } +intros k Hk fexp x y Px Py Hln Hlxy Fx. +revert Fx. +unfold round, generic_format, F2R, scaled_mantissa, canonic_exp; simpl. +rewrite Hlxy. +set (mx := Ztrunc (x * bpow (- fexp (ln_beta x)))). +intros Fx. +assert (R : (x + y) * bpow (- fexp (ln_beta x)) + = Z2R mx + y * bpow (- fexp (ln_beta x))). +{ rewrite Fx at 1. + rewrite Rmult_plus_distr_r. + now bpow_simplify. } +rewrite R. +assert (LB : 0 < y * bpow (- fexp (ln_beta x))). +{ rewrite <- (Rmult_0_r y). + now apply Rmult_lt_compat_l; [|apply bpow_gt_0]. } +assert (UB : y * bpow (- fexp (ln_beta x)) < / Z2R (beta ^ k)). +{ apply Rlt_le_trans with (bpow (ln_beta y) * bpow (- fexp (ln_beta x))). + - apply Rmult_lt_compat_r; [now apply bpow_gt_0|]. + rewrite <- (Rabs_right y) at 1; [|now apply Rle_ge; apply Rlt_le]. + apply bpow_ln_beta_gt. + - apply Rle_trans with (bpow (fexp (ln_beta x) - k) + * bpow (- fexp (ln_beta x)))%R. + + apply Rmult_le_compat_r; [now apply bpow_ge_0|]. + now apply bpow_le. + + bpow_simplify. + rewrite bpow_opp. + destruct k. + * omega. + * simpl; unfold Fcore_Raux.bpow, Z.pow_pos. + now apply Rle_refl. + * casetype False; apply (Zlt_irrefl 0). + apply (Zlt_trans _ _ _ Hk). + apply Zlt_neg_0. } +rewrite (Zfloor_imp mx). +{ split; ring_simplify. + - apply (Rmult_lt_reg_r (bpow (- fexp (ln_beta x)))); [now apply bpow_gt_0|]. + rewrite Rmult_minus_distr_r, Rmult_0_l. + bpow_simplify. + rewrite R; ring_simplify. + now apply Rmult_lt_0_compat; [|apply bpow_gt_0]. + - apply (Rmult_lt_reg_r (bpow (- fexp (ln_beta x)))); [now apply bpow_gt_0|]. + rewrite Rmult_minus_distr_r. + bpow_simplify. + rewrite R; ring_simplify. + apply (Rlt_le_trans _ _ _ UB). + rewrite bpow_opp. + apply Rinv_le; [now apply bpow_gt_0|]. + now rewrite Z2R_Zpower; [right|omega]. } +split. +- rewrite <- Rplus_0_r at 1; apply Rplus_le_compat_l. + now apply Rlt_le. +- rewrite Z2R_plus; apply Rplus_lt_compat_l. + apply (Rmult_lt_reg_r (bpow (fexp (ln_beta x)))); [now apply bpow_gt_0|]. + rewrite Rmult_1_l. + bpow_simplify. + apply Rlt_trans with (bpow (ln_beta y)). + + rewrite <- Rabs_right at 1; [|now apply Rle_ge; apply Rlt_le]. + apply bpow_ln_beta_gt. + + apply bpow_lt; omega. +Qed. + +(* ln_beta y <= fexp1 (ln_beta x) - 2 : double_round_lt_mid applies. *) +Lemma double_round_plus_aux1 : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_plus_hyp fexp1 fexp2 -> + forall x y, + 0 < x -> 0 < y -> + (ln_beta y <= fexp1 (ln_beta x) - 2)%Z -> + generic_format beta fexp1 x -> + double_round_eq fexp1 fexp2 choice1 choice2 (x + y). +Proof. +assert (Hbeta : (2 <= beta)%Z). +{ destruct beta as (beta_val,beta_prop). + now apply Zle_bool_imp_le. } +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x y Px Py Hly Fx. +assert (Lxy : ln_beta (x + y) = ln_beta x :> Z); + [now apply (ln_beta_plus_separated fexp1); [|apply Rlt_le| |omega]|]. +destruct Hexp as (_,(_,(_,Hexp4))). +assert (Hf2 : (fexp2 (ln_beta x) <= fexp1 (ln_beta x))%Z); + [now apply Hexp4; omega|]. +assert (Bpow2 : bpow (- 2) <= / 2 * / 2). +{ unfold Fcore_Raux.bpow, Z.pow_pos; simpl. + rewrite <- Rinv_mult_distr; [|lra|lra]. + apply Rinv_le; [lra|]. + change 4 with (Z2R (2 * 2)); apply Z2R_le; rewrite Zmult_1_r. + now apply Zmult_le_compat; omega. } +assert (P2 : (0 < 2)%Z) by omega. +unfold double_round_eq. +apply double_round_lt_mid. +- exact Vfexp1. +- exact Vfexp2. +- lra. +- now rewrite Lxy. +- rewrite Lxy. + assert (fexp1 (ln_beta x) < ln_beta x)%Z; [|omega]. + now apply ln_beta_generic_gt; [|apply Rgt_not_eq|]. +- unfold midp. + apply (Rplus_lt_reg_r (- round beta fexp1 Zfloor (x + y))). + 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. + 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. + intro Hf2'. + apply (Rmult_lt_reg_r (bpow (- fexp1 (ln_beta x)))); + [now apply bpow_gt_0|]. + apply (Rmult_lt_reg_r (bpow (fexp1 (ln_beta x)))); [now apply bpow_gt_0|]. + bpow_simplify. + apply (Rplus_lt_reg_r (- round beta fexp1 Zfloor (x + y))). + unfold midp; ring_simplify. + 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. + apply (Rle_trans _ _ _ Bpow2). + rewrite <- (Rmult_1_r (/ 2)) at 3; rewrite <- Rmult_minus_distr_l. + apply Rmult_le_compat_l; [lra|]. + apply (Rplus_le_reg_r (- 1)); ring_simplify. + replace (_ - _) with (- (/ 2)) by lra. + apply Ropp_le_contravar. + { apply Rle_trans with (bpow (- 1)). + - apply bpow_le; omega. + - unfold Fcore_Raux.bpow, Z.pow_pos; simpl. + apply Rinv_le; [lra|]. + change 2 with (Z2R 2); apply Z2R_le; omega. } +Qed. + +(* double_round_plus_aux{0,1} together *) +Lemma double_round_plus_aux2 : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_plus_hyp fexp1 fexp2 -> + forall x y, + 0 < x -> 0 < y -> y <= x -> + generic_format beta fexp1 x -> + generic_format beta fexp1 y -> + double_round_eq fexp1 fexp2 choice1 choice2 (x + y). +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x y Px Py Hyx Fx Fy. +unfold double_round_eq. +destruct (Zle_or_lt (ln_beta y) (fexp1 (ln_beta x) - 2)) as [Hly|Hly]. +- (* ln_beta y <= fexp1 (ln_beta x) - 2 *) + now apply double_round_plus_aux1. +- (* fexp1 (ln_beta x) - 2 < ln_beta y *) + rewrite (round_generic beta fexp2). + + reflexivity. + + now apply valid_rnd_N. + + assert (Hf1 : (fexp1 (ln_beta x) - 1 <= ln_beta y)%Z); [omega|]. + now apply (double_round_plus_aux0 fexp1). +Qed. + +Lemma double_round_plus_aux : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_plus_hyp fexp1 fexp2 -> + forall x y, + 0 <= x -> 0 <= y -> + generic_format beta fexp1 x -> + generic_format beta fexp1 y -> + double_round_eq fexp1 fexp2 choice1 choice2 (x + y). +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x y Nnx Nny Fx Fy. +unfold double_round_eq. +destruct (Req_dec x 0) as [Zx|Nzx]. +- (* x = 0 *) + destruct Hexp as (_,(_,(_,Hexp4))). + rewrite Zx; rewrite Rplus_0_l. + rewrite (round_generic beta fexp2). + + reflexivity. + + now apply valid_rnd_N. + + apply (generic_inclusion_ln_beta beta fexp1). + now intros _; apply Hexp4; omega. + exact Fy. +- (* x <> 0 *) + destruct (Req_dec y 0) as [Zy|Nzy]. + + (* y = 0 *) + destruct Hexp as (_,(_,(_,Hexp4))). + rewrite Zy; rewrite Rplus_0_r. + rewrite (round_generic beta fexp2). + * reflexivity. + * now apply valid_rnd_N. + * apply (generic_inclusion_ln_beta beta fexp1). + now intros _; apply Hexp4; omega. + exact Fx. + + (* y <> 0 *) + assert (Px : 0 < x); [lra|]. + assert (Py : 0 < y); [lra|]. + destruct (Rlt_or_le x y) as [H|H]. + * (* x < y *) + apply Rlt_le in H. + rewrite Rplus_comm. + now apply double_round_plus_aux2. + * now apply double_round_plus_aux2. +Qed. + +Lemma double_round_minus_aux0_aux : + forall (fexp1 fexp2 : Z -> Z), + forall x y, + (fexp2 (ln_beta (x - y))%Z <= fexp1 (ln_beta x))%Z -> + (fexp2 (ln_beta (x - y))%Z <= fexp1 (ln_beta y))%Z -> + generic_format beta fexp1 x -> generic_format beta fexp1 y -> + generic_format beta fexp2 (x - y). +Proof. +intros fexp1 fexp2 x y. +replace (x - y)%R with (x + (- y))%R; [|ring]. +intros Hlnx Hlny Fx Fy. +rewrite <- (ln_beta_opp beta y) in Hlny. +apply generic_format_opp in Fy. +now apply (double_round_plus_aux0_aux fexp1). +Qed. + +(* fexp1 (ln_beta x) - 1 <= ln_beta y : + * substraction is exact in the largest precision (fexp2). *) +Lemma double_round_minus_aux0 : + forall (fexp1 fexp2 : Z -> Z), + double_round_plus_hyp fexp1 fexp2 -> + forall x y, + 0 < y -> y < x -> + (fexp1 (ln_beta x) - 1 <= ln_beta y)%Z -> + generic_format beta fexp1 x -> generic_format beta fexp1 y -> + generic_format beta fexp2 (x - y). +Proof. +intros fexp1 fexp2 Hexp x y Py Hyx Hln Fx Fy. +assert (Px := Rlt_trans 0 y x Py Hyx). +destruct Hexp as (Hexp1,(_,(Hexp3,Hexp4))). +assert (Lyx : (ln_beta y <= ln_beta x)%Z); + [now apply ln_beta_le; [|apply Rlt_le]|]. +destruct (Z.lt_ge_cases (ln_beta x - 2) (ln_beta y)) as [Hlt|Hge]. +- (* ln_beta x - 2 < ln_beta y *) + assert (Hor : (ln_beta y = ln_beta x :> Z) + \/ (ln_beta y = ln_beta x - 1 :> Z)%Z); [omega|]. + destruct Hor as [Heq|Heqm1]. + + (* ln_beta y = ln_beta x *) + apply (double_round_minus_aux0_aux fexp1); [| |exact Fx|exact Fy]. + * apply Hexp4. + apply Zle_trans with (ln_beta (x - y)); [omega|]. + now apply ln_beta_minus. + * rewrite Heq. + apply Hexp4. + apply Zle_trans with (ln_beta (x - y)); [omega|]. + now apply ln_beta_minus. + + (* ln_beta y = ln_beta x - 1 *) + apply (double_round_minus_aux0_aux fexp1); [| |exact Fx|exact Fy]. + * apply Hexp4. + apply Zle_trans with (ln_beta (x - y)); [omega|]. + now apply ln_beta_minus. + * rewrite Heqm1. + apply Hexp4. + apply Zplus_le_compat_r. + now apply ln_beta_minus. +- (* ln_beta y <= ln_beta x - 2 *) + destruct (ln_beta_minus_disj x y Px Py Hge) as [Lxmy|Lxmy]. + + (* ln_beta (x - y) = ln_beta x *) + apply (double_round_minus_aux0_aux fexp1); [| |exact Fx|exact Fy]. + * apply Hexp4. + omega. + * now rewrite Lxmy; apply Hexp3. + + (* ln_beta (x - y) = ln_beta x - 1 *) + apply (double_round_minus_aux0_aux fexp1); [| |exact Fx|exact Fy]; + rewrite Lxmy. + * apply Hexp1. + replace (_ + _)%Z with (ln_beta x : Z); [|ring]. + now apply Zle_trans with (ln_beta y). + * apply Hexp1. + now replace (_ + _)%Z with (ln_beta x : Z); [|ring]. +Qed. + +(* ln_beta y <= fexp1 (ln_beta x) - 2, + * fexp1 (ln_beta (x - y)) - 1 <= ln_beta y : + * substraction is exact in the largest precision (fexp2). *) +Lemma double_round_minus_aux1 : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + double_round_plus_hyp fexp1 fexp2 -> + forall x y, + 0 < y -> y < x -> + (ln_beta y <= fexp1 (ln_beta x) - 2)%Z -> + (fexp1 (ln_beta (x - y)) - 1 <= ln_beta y)%Z -> + generic_format beta fexp1 x -> generic_format beta fexp1 y -> + generic_format beta fexp2 (x - y). +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 Hexp x y Py Hyx Hln Hln' Fx Fy. +assert (Px := Rlt_trans 0 y x Py Hyx). +destruct Hexp as (Hexp1,(Hexp2,(Hexp3,Hexp4))). +assert (Lyx : (ln_beta y <= ln_beta x)%Z); + [now apply ln_beta_le; [|apply Rlt_le]|]. +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|]|]. +apply (double_round_minus_aux0_aux fexp1); [| |exact Fx|exact Fy]. +- apply Zle_trans with (fexp1 (ln_beta (x - y))). + + apply Hexp4; omega. + + omega. +- now apply Hexp3. +Qed. + +Lemma double_round_minus_aux2_aux : + forall (fexp : Z -> Z), + Valid_exp fexp -> + forall x y, + 0 < y -> y < x -> + (ln_beta y <= fexp (ln_beta x) - 1)%Z -> + generic_format beta fexp x -> + generic_format beta fexp y -> + round beta fexp Zceil (x - y) - (x - y) <= y. +Proof. +intros fexp Vfexp x y Py Hxy Hly Fx Fy. +assert (Px := Rlt_trans 0 y x Py Hxy). +revert Fx. +unfold generic_format, F2R, scaled_mantissa, canonic_exp; simpl. +set (mx := Ztrunc (x * bpow (- fexp (ln_beta x)))). +intro Fx. +assert (Hfx : (fexp (ln_beta x) < ln_beta x)%Z); + [now apply ln_beta_generic_gt; [|apply Rgt_not_eq|]|]. +assert (Hfy : (fexp (ln_beta y) < ln_beta y)%Z); + [now apply ln_beta_generic_gt; [|apply Rgt_not_eq|]|]. +destruct (Rlt_or_le (bpow (ln_beta x - 1)) x) as [Hx|Hx]. +- (* bpow (ln_beta x - 1) < x *) + assert (Lxy : ln_beta (x - y) = ln_beta x :> Z); + [now apply (ln_beta_minus_separated fexp); [| | | | | |omega]|]. + assert (Rxy : round beta fexp Zceil (x - y) = x). + { unfold round, F2R, scaled_mantissa, canonic_exp; simpl. + rewrite Lxy. + apply eq_sym; rewrite Fx at 1; apply eq_sym. + apply Rmult_eq_compat_r. + apply f_equal. + rewrite Fx at 1. + rewrite Rmult_minus_distr_r. + bpow_simplify. + apply Zceil_imp. + split. + - unfold Zminus; rewrite Z2R_plus. + apply Rplus_lt_compat_l. + apply Ropp_lt_contravar; simpl. + apply (Rmult_lt_reg_r (bpow (fexp (ln_beta x)))); + [now apply bpow_gt_0|]. + rewrite Rmult_1_l; bpow_simplify. + apply Rlt_le_trans with (bpow (ln_beta y)). + + rewrite <- Rabs_right at 1; [|now apply Rle_ge; apply Rlt_le]. + apply bpow_ln_beta_gt. + + apply bpow_le. + omega. + - rewrite <- (Rplus_0_r (Z2R _)) at 2. + apply Rplus_le_compat_l. + rewrite <- Ropp_0; apply Ropp_le_contravar. + rewrite <- (Rmult_0_r y). + apply Rmult_le_compat_l; [now apply Rlt_le|]. + now apply bpow_ge_0. } + rewrite Rxy; ring_simplify. + apply Rle_refl. +- (* x <= bpow (ln_beta x - 1) *) + assert (Xpow : x = bpow (ln_beta x - 1)). + { apply Rle_antisym; [exact Hx|]. + destruct (ln_beta x) as (ex, Hex); simpl. + rewrite <- (Rabs_right x); [|now apply Rle_ge; apply Rlt_le]. + apply Hex. + now apply Rgt_not_eq. } + assert (Lxy : (ln_beta (x - y) = ln_beta x - 1 :> Z)%Z). + { apply Zle_antisym. + - apply ln_beta_le_bpow. + + apply Rminus_eq_contra. + now intro Hx'; rewrite Hx' in Hxy; apply (Rlt_irrefl y). + + rewrite Rabs_right; lra. + - apply (ln_beta_minus_lb beta x y Px Py). + omega. } + assert (Hfx1 : (fexp (ln_beta x - 1) < ln_beta x - 1)%Z); + [now apply (valid_exp_large fexp (ln_beta y)); [|omega]|]. + assert (Rxy : round beta fexp Zceil (x - y) <= x). + { rewrite Xpow at 2. + unfold round, F2R, scaled_mantissa, canonic_exp; simpl. + rewrite Lxy. + apply (Rmult_le_reg_r (bpow (- fexp (ln_beta x - 1)%Z))); + [now apply bpow_gt_0|]. + bpow_simplify. + rewrite <- (Z2R_Zpower beta (_ - _ - _)); [|omega]. + apply Z2R_le. + apply Zceil_glb. + rewrite Z2R_Zpower; [|omega]. + rewrite Xpow at 1. + rewrite Rmult_minus_distr_r. + bpow_simplify. + rewrite <- (Rplus_0_r (bpow _)) at 2. + apply Rplus_le_compat_l. + rewrite <- Ropp_0; apply Ropp_le_contravar. + rewrite <- (Rmult_0_r y). + apply Rmult_le_compat_l; [now apply Rlt_le|]. + now apply bpow_ge_0. } + lra. +Qed. + +(* ln_beta y <= fexp1 (ln_beta x) - 2 : + * ln_beta y <= fexp1 (ln_beta (x - y)) - 2 : + * double_round_gt_mid applies. *) +Lemma double_round_minus_aux2 : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_plus_hyp fexp1 fexp2 -> + forall x y, + 0 < y -> y < x -> + (ln_beta y <= fexp1 (ln_beta x) - 2)%Z -> + (ln_beta y <= fexp1 (ln_beta (x - y)) - 2)%Z -> + generic_format beta fexp1 x -> + generic_format beta fexp1 y -> + double_round_eq fexp1 fexp2 choice1 choice2 (x - y). +Proof. +assert (Hbeta : (2 <= beta)%Z). +{ destruct beta as (beta_val,beta_prop). + now apply Zle_bool_imp_le. } +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x y Py Hxy Hly Hly' Fx Fy. +assert (Px := Rlt_trans 0 y x Py Hxy). +destruct Hexp as (_,(_,(_,Hexp4))). +assert (Hf2 : (fexp2 (ln_beta x) <= fexp1 (ln_beta x))%Z); + [now apply Hexp4; omega|]. +assert (Hfx : (fexp1 (ln_beta x) < ln_beta x)%Z); + [now apply ln_beta_generic_gt; [|apply Rgt_not_eq|]|]. +assert (Bpow2 : bpow (- 2) <= / 2 * / 2). +{ unfold Fcore_Raux.bpow, Z.pow_pos; simpl. + rewrite <- Rinv_mult_distr; [|lra|lra]. + apply Rinv_le; [lra|]. + change 4 with (Z2R (2 * 2)); apply Z2R_le; rewrite Zmult_1_r. + now apply Zmult_le_compat; omega. } +assert (Ly : y < bpow (ln_beta y)). +{ apply Rabs_lt_inv. + apply bpow_ln_beta_gt. } +unfold double_round_eq. +apply double_round_gt_mid. +- exact Vfexp1. +- exact Vfexp2. +- lra. +- apply Hexp4; omega. +- assert (fexp1 (ln_beta (x - y)) < ln_beta (x - y))%Z; [|omega]. + apply (valid_exp_large fexp1 (ln_beta x - 1)). + + apply (valid_exp_large fexp1 (ln_beta y)); [|omega]. + now apply ln_beta_generic_gt; [|apply Rgt_not_eq|]. + + now apply ln_beta_minus_lb; [| |omega]. +- unfold midp'. + apply (Rplus_lt_reg_r (/ 2 * ulp beta fexp1 (x - y) - (x - y))). + ring_simplify. + replace (_ + _) with (round beta fexp1 Zceil (x - y) - (x - y)) by ring. + apply Rlt_le_trans with (bpow (fexp1 (ln_beta (x - y)) - 2)). + + apply Rle_lt_trans with y; + [now apply double_round_minus_aux2_aux; try assumption; omega|]. + apply (Rlt_le_trans _ _ _ Ly). + now apply bpow_le. + + unfold ulp, 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. + apply Rmult_le_compat. + * now apply bpow_ge_0. + * now apply bpow_ge_0. + * 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. + * apply bpow_le; omega. +- intro Hf2'. + unfold midp'. + apply (Rplus_lt_reg_r (/ 2 * ulp beta fexp1 (x - y) - (x - y) + - / 2 * ulp beta fexp2 (x - y))). + ring_simplify. + replace (_ + _) with (round beta fexp1 Zceil (x - y) - (x - y)) by ring. + apply Rle_lt_trans with y; + [now apply double_round_minus_aux2_aux; try assumption; omega|]. + apply (Rlt_le_trans _ _ _ Ly). + apply Rle_trans with (bpow (fexp1 (ln_beta (x - y)) - 2)); + [now apply bpow_le|]. + replace (_ - 2)%Z with (fexp1 (ln_beta (x - y)) - 1 - 1)%Z by ring. + unfold Zminus at 1; rewrite bpow_plus. + rewrite <- Rmult_minus_distr_l. + rewrite Rmult_comm; apply Rmult_le_compat. + + apply bpow_ge_0. + + apply bpow_ge_0. + + 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. + 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. + apply Rplus_le_compat_l. + now apply bpow_le. + * unfold Zminus; rewrite bpow_plus. + rewrite Rmult_comm; rewrite Rmult_assoc. + rewrite <- Rmult_1_r. + apply Rmult_le_compat_l; [now apply bpow_ge_0|]. + unfold Fcore_Raux.bpow, Z.pow_pos; simpl. + rewrite Zmult_1_r. + rewrite <- (Rinv_r 2) at 3; [|lra]. + rewrite Rmult_comm; apply Rmult_le_compat_l; [lra|]. + apply Rinv_le; [lra|]. + now change 2 with (Z2R 2); apply Z2R_le. +Qed. + +(* double_round_minus_aux{0,1,2} together *) +Lemma double_round_minus_aux3 : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_plus_hyp fexp1 fexp2 -> + forall x y, + 0 < y -> y <= x -> + generic_format beta fexp1 x -> + generic_format beta fexp1 y -> + double_round_eq fexp1 fexp2 choice1 choice2 (x - y). +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x y Py Hyx Fx Fy. +assert (Px := Rlt_le_trans 0 y x Py Hyx). +unfold double_round_eq. +destruct (Req_dec y x) as [Hy|Hy]. +- (* y = x *) + rewrite Hy; replace (x - x) with 0 by ring. + rewrite round_0. + + reflexivity. + + now apply valid_rnd_N. +- (* y < x *) + assert (Hyx' : y < x); [lra|]. + destruct (Zle_or_lt (ln_beta y) (fexp1 (ln_beta x) - 2)) as [Hly|Hly]. + + (* ln_beta y <= fexp1 (ln_beta x) - 2 *) + destruct (Zle_or_lt (ln_beta y) (fexp1 (ln_beta (x - y)) - 2)) + as [Hly'|Hly']. + * (* ln_beta y <= fexp1 (ln_beta (x - y)) - 2 *) + now apply double_round_minus_aux2. + * (* fexp1 (ln_beta (x - y)) - 2 < ln_beta y *) + { rewrite (round_generic beta fexp2). + - reflexivity. + - now apply valid_rnd_N. + - assert (Hf1 : (fexp1 (ln_beta (x - y)) - 1 <= ln_beta y)%Z); [omega|]. + now apply (double_round_minus_aux1 fexp1). } + + rewrite (round_generic beta fexp2). + * reflexivity. + * now apply valid_rnd_N. + * assert (Hf1 : (fexp1 (ln_beta x) - 1 <= ln_beta y)%Z); [omega|]. + now apply (double_round_minus_aux0 fexp1). +Qed. + +Lemma double_round_minus_aux : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_plus_hyp fexp1 fexp2 -> + forall x y, + 0 <= x -> 0 <= y -> + generic_format beta fexp1 x -> + generic_format beta fexp1 y -> + double_round_eq fexp1 fexp2 choice1 choice2 (x - y). +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x y Nnx Nny Fx Fy. +unfold double_round_eq. +destruct (Req_dec x 0) as [Zx|Nzx]. +- (* x = 0 *) + rewrite Zx; unfold Rminus; rewrite Rplus_0_l. + do 3 rewrite round_N_opp. + rewrite (round_generic beta fexp2). + * reflexivity. + * now apply valid_rnd_N. + * apply (generic_inclusion_ln_beta beta fexp1). + destruct Hexp as (_,(_,(_,Hexp4))). + now intros _; apply Hexp4; omega. + exact Fy. +- (* x <> 0 *) + destruct (Req_dec y 0) as [Zy|Nzy]. + + (* y = 0 *) + rewrite Zy; unfold Rminus; rewrite Ropp_0; rewrite Rplus_0_r. + rewrite (round_generic beta fexp2). + * reflexivity. + * now apply valid_rnd_N. + * apply (generic_inclusion_ln_beta beta fexp1). + destruct Hexp as (_,(_,(_,Hexp4))). + now intros _; apply Hexp4; omega. + exact Fx. + + (* y <> 0 *) + assert (Px : 0 < x); [lra|]. + assert (Py : 0 < y); [lra|]. + destruct (Rlt_or_le x y) as [H|H]. + * (* x < y *) + apply Rlt_le in H. + replace (x - y) with (- (y - x)) by ring. + do 3 rewrite round_N_opp. + apply Ropp_eq_compat. + now apply double_round_minus_aux3. + * (* y <= x *) + now apply double_round_minus_aux3. +Qed. + +Lemma double_round_plus : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_plus_hyp fexp1 fexp2 -> + forall x y, + generic_format beta fexp1 x -> + generic_format beta fexp1 y -> + double_round_eq fexp1 fexp2 choice1 choice2 (x + y). +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x y Fx Fy. +unfold double_round_eq. +destruct (Rlt_or_le x 0) as [Sx|Sx]; destruct (Rlt_or_le y 0) as [Sy|Sy]. +- (* x < 0, y < 0 *) + replace (x + y) with (- (- x - y)); [|ring]. + do 3 rewrite round_N_opp. + apply Ropp_eq_compat. + assert (Px : 0 <= - x); [lra|]. + assert (Py : 0 <= - y); [lra|]. + apply generic_format_opp in Fx. + apply generic_format_opp in Fy. + now apply double_round_plus_aux. +- (* x < 0, 0 <= y *) + replace (x + y) with (y - (- x)); [|ring]. + assert (Px : 0 <= - x); [lra|]. + apply generic_format_opp in Fx. + now apply double_round_minus_aux. +- (* 0 <= x, y < 0 *) + replace (x + y) with (x - (- y)); [|ring]. + assert (Py : 0 <= - y); [lra|]. + apply generic_format_opp in Fy. + now apply double_round_minus_aux. +- (* 0 <= x, 0 <= y *) + now apply double_round_plus_aux. +Qed. + +Lemma double_round_minus : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_plus_hyp fexp1 fexp2 -> + forall x y, + generic_format beta fexp1 x -> + generic_format beta fexp1 y -> + double_round_eq fexp1 fexp2 choice1 choice2 (x - y). +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x y Fx Fy. +unfold Rminus. +apply generic_format_opp in Fy. +now apply double_round_plus. +Qed. + +Section Double_round_plus_FLX. + +Require Import Fcore_FLX. + +Variable prec : Z. +Variable prec' : Z. + +Context { prec_gt_0_ : Prec_gt_0 prec }. +Context { prec_gt_0_' : Prec_gt_0 prec' }. + +Lemma FLX_double_round_plus_hyp : + (2 * prec + 1 <= prec')%Z -> + double_round_plus_hyp (FLX_exp prec) (FLX_exp prec'). +Proof. +intros Hprec. +unfold FLX_exp. +unfold double_round_plus_hyp; split; [|split; [|split]]; +intros ex ey; try omega. +unfold Prec_gt_0 in prec_gt_0_. +omega. +Qed. + +Theorem double_round_plus_FLX : + forall choice1 choice2, + (2 * prec + 1 <= prec')%Z -> + forall x y, + FLX_format beta prec x -> FLX_format beta prec y -> + double_round_eq (FLX_exp prec) (FLX_exp prec') choice1 choice2 (x + y). +Proof. +intros choice1 choice2 Hprec x y Fx Fy. +apply double_round_plus. +- now apply FLX_exp_valid. +- now apply FLX_exp_valid. +- now apply FLX_double_round_plus_hyp. +- now apply generic_format_FLX. +- now apply generic_format_FLX. +Qed. + +Theorem double_round_minus_FLX : + forall choice1 choice2, + (2 * prec + 1 <= prec')%Z -> + forall x y, + FLX_format beta prec x -> FLX_format beta prec y -> + double_round_eq (FLX_exp prec) (FLX_exp prec') choice1 choice2 (x - y). +Proof. +intros choice1 choice2 Hprec x y Fx Fy. +apply double_round_minus. +- now apply FLX_exp_valid. +- now apply FLX_exp_valid. +- now apply FLX_double_round_plus_hyp. +- now apply generic_format_FLX. +- now apply generic_format_FLX. +Qed. + +End Double_round_plus_FLX. + +Section Double_round_plus_FLT. + +Require Import Fcore_FLX. +Require Import Fcore_FLT. + +Variable emin prec : Z. +Variable emin' prec' : Z. + +Context { prec_gt_0_ : Prec_gt_0 prec }. +Context { prec_gt_0_' : Prec_gt_0 prec' }. + +Lemma FLT_double_round_plus_hyp : + (emin' <= emin)%Z -> (2 * prec + 1 <= prec')%Z -> + double_round_plus_hyp (FLT_exp emin prec) (FLT_exp emin' prec'). +Proof. +intros Hemin Hprec. +unfold FLT_exp. +unfold double_round_plus_hyp; split; [|split; [|split]]; intros ex ey. +- generalize (Zmax_spec (ex + 1 - prec) emin). + generalize (Zmax_spec (ex - prec') emin'). + generalize (Zmax_spec (ey - prec) emin). + omega. +- generalize (Zmax_spec (ex - 1 - prec) emin). + generalize (Zmax_spec (ex - prec') emin'). + generalize (Zmax_spec (ey - prec) emin). + omega. +- generalize (Zmax_spec (ex - prec) emin). + generalize (Zmax_spec (ex - prec') emin'). + generalize (Zmax_spec (ey - prec) emin). + omega. +- unfold Prec_gt_0 in prec_gt_0_. + generalize (Zmax_spec (ex - prec') emin'). + generalize (Zmax_spec (ey - prec) emin). + omega. +Qed. + +Theorem double_round_plus_FLT : + forall choice1 choice2, + (emin' <= emin)%Z -> (2 * prec + 1 <= prec')%Z -> + forall x y, + FLT_format beta emin prec x -> FLT_format beta emin prec y -> + double_round_eq (FLT_exp emin prec) (FLT_exp emin' prec') + choice1 choice2 (x + y). +Proof. +intros choice1 choice2 Hemin Hprec x y Fx Fy. +apply double_round_plus. +- now apply FLT_exp_valid. +- now apply FLT_exp_valid. +- now apply FLT_double_round_plus_hyp. +- now apply generic_format_FLT. +- now apply generic_format_FLT. +Qed. + +Theorem double_round_minus_FLT : + forall choice1 choice2, + (emin' <= emin)%Z -> (2 * prec + 1 <= prec')%Z -> + forall x y, + FLT_format beta emin prec x -> FLT_format beta emin prec y -> + double_round_eq (FLT_exp emin prec) (FLT_exp emin' prec') + choice1 choice2 (x - y). +Proof. +intros choice1 choice2 Hemin Hprec x y Fx Fy. +apply double_round_minus. +- now apply FLT_exp_valid. +- now apply FLT_exp_valid. +- now apply FLT_double_round_plus_hyp. +- now apply generic_format_FLT. +- now apply generic_format_FLT. +Qed. + +End Double_round_plus_FLT. + +Section Double_round_plus_FTZ. + +Require Import Fcore_FLX. +Require Import Fcore_FTZ. + +Variable emin prec : Z. +Variable emin' prec' : Z. + +Context { prec_gt_0_ : Prec_gt_0 prec }. +Context { prec_gt_0_' : Prec_gt_0 prec' }. + +Lemma FTZ_double_round_plus_hyp : + (emin' + prec' <= emin + 1)%Z -> (2 * prec + 1 <= prec')%Z -> + double_round_plus_hyp (FTZ_exp emin prec) (FTZ_exp emin' prec'). +Proof. +intros Hemin Hprec. +unfold FTZ_exp. +unfold Prec_gt_0 in *. +unfold double_round_plus_hyp; split; [|split; [|split]]; intros ex ey. +- destruct (Z.ltb_spec (ex + 1 - prec) emin); + destruct (Z.ltb_spec (ex - prec') emin'); + destruct (Z.ltb_spec (ey - prec) emin); + omega. +- destruct (Z.ltb_spec (ex - 1 - prec) emin); + destruct (Z.ltb_spec (ex - prec') emin'); + destruct (Z.ltb_spec (ey - prec) emin); + omega. +- destruct (Z.ltb_spec (ex - prec) emin); + destruct (Z.ltb_spec (ex - prec') emin'); + destruct (Z.ltb_spec (ey - prec) emin); + omega. +- destruct (Z.ltb_spec (ex - prec') emin'); + destruct (Z.ltb_spec (ey - prec) emin); + omega. +Qed. + +Theorem double_round_plus_FTZ : + forall choice1 choice2, + (emin' + prec' <= emin + 1)%Z -> (2 * prec + 1 <= prec')%Z -> + forall x y, + FTZ_format beta emin prec x -> FTZ_format beta emin prec y -> + double_round_eq (FTZ_exp emin prec) (FTZ_exp emin' prec') + choice1 choice2 (x + y). +Proof. +intros choice1 choice2 Hemin Hprec x y Fx Fy. +apply double_round_plus. +- now apply FTZ_exp_valid. +- now apply FTZ_exp_valid. +- now apply FTZ_double_round_plus_hyp. +- now apply generic_format_FTZ. +- now apply generic_format_FTZ. +Qed. + +Theorem double_round_minus_FTZ : + forall choice1 choice2, + (emin' + prec' <= emin + 1)%Z -> (2 * prec + 1 <= prec')%Z -> + forall x y, + FTZ_format beta emin prec x -> FTZ_format beta emin prec y -> + double_round_eq (FTZ_exp emin prec) (FTZ_exp emin' prec') + choice1 choice2 (x - y). +Proof. +intros choice1 choice2 Hemin Hprec x y Fx Fy. +apply double_round_minus. +- now apply FTZ_exp_valid. +- now apply FTZ_exp_valid. +- now apply FTZ_double_round_plus_hyp. +- now apply generic_format_FTZ. +- now apply generic_format_FTZ. +Qed. + +End Double_round_plus_FTZ. + +Section Double_round_plus_beta_ge_3. + +Definition double_round_plus_beta_ge_3_hyp fexp1 fexp2 := + (forall ex ey, (fexp1 (ex + 1) <= ey)%Z -> (fexp2 ex <= fexp1 ey)%Z) + /\ (forall ex ey, (fexp1 (ex - 1) + 1 <= ey)%Z -> (fexp2 ex <= fexp1 ey)%Z) + /\ (forall ex ey, (fexp1 ex <= ey)%Z -> (fexp2 ex <= fexp1 ey)%Z) + /\ (forall ex ey, (ex - 1 <= ey)%Z -> (fexp2 ex <= fexp1 ey)%Z). + +(* fexp1 (ln_beta x) <= ln_beta y : + * addition is exact in the largest precision (fexp2). *) +Lemma double_round_plus_beta_ge_3_aux0 : + forall (fexp1 fexp2 : Z -> Z), Valid_exp fexp1 -> + double_round_plus_beta_ge_3_hyp fexp1 fexp2 -> + forall x y, + (0 < y)%R -> (y <= x)%R -> + (fexp1 (ln_beta x) <= ln_beta y)%Z -> + generic_format beta fexp1 x -> generic_format beta fexp1 y -> + generic_format beta fexp2 (x + y). +Proof. +intros fexp1 fexp2 Vfexp1 Hexp x y Py Hyx Hln Fx Fy. +assert (Px := Rlt_le_trans 0 y x Py Hyx). +assert (Nny : (0 <= y)%R); [now apply Rlt_le|]. +destruct Hexp as (_,(Hexp2,(Hexp3,Hexp4))). +destruct (Z.le_gt_cases (ln_beta y) (fexp1 (ln_beta x))) as [Hle|Hgt]. +- (* ln_beta y <= fexp1 (ln_beta x) *) + assert (Lxy : ln_beta (x + y) = ln_beta x :> Z); + [now apply (ln_beta_plus_separated fexp1)|]. + apply (double_round_plus_aux0_aux fexp1); + [| |assumption|assumption]; rewrite Lxy. + + now apply Hexp4; omega. + + now apply Hexp3; omega. +- (* fexp1 (ln_beta x) < ln_beta y *) + apply (double_round_plus_aux0_aux fexp1); [| |assumption|assumption]. + destruct (ln_beta_plus_disj x y Py Hyx) as [Lxy|Lxy]; rewrite Lxy. + + now apply Hexp4; omega. + + apply Hexp2; apply (ln_beta_le beta y x Py) in Hyx. + replace (_ - _)%Z with (ln_beta x : Z) by ring. + omega. + + destruct (ln_beta_plus_disj x y Py Hyx) as [Lxy|Lxy]; rewrite Lxy. + * now apply Hexp3; omega. + * apply Hexp2. + replace (_ - _)%Z with (ln_beta x : Z) by ring. + omega. +Qed. + +(* ln_beta y <= fexp1 (ln_beta x) - 1 : double_round_lt_mid applies. *) +Lemma double_round_plus_beta_ge_3_aux1 : + (3 <= beta)%Z -> + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_plus_beta_ge_3_hyp fexp1 fexp2 -> + forall x y, + 0 < x -> 0 < y -> + (ln_beta y <= fexp1 (ln_beta x) - 1)%Z -> + generic_format beta fexp1 x -> + double_round_eq fexp1 fexp2 choice1 choice2 (x + y). +Proof. +intros Hbeta fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x y Px Py Hly Fx. +assert (Lxy : ln_beta (x + y) = ln_beta x :> Z); + [now apply (ln_beta_plus_separated fexp1); [|apply Rlt_le| |omega]|]. +destruct Hexp as (_,(_,(_,Hexp4))). +assert (Hf2 : (fexp2 (ln_beta x) <= fexp1 (ln_beta x))%Z); + [now apply Hexp4; omega|]. +assert (Bpow3 : bpow (- 1) <= / 3). +{ unfold Fcore_Raux.bpow, Z.pow_pos; simpl. + rewrite Zmult_1_r. + apply Rinv_le; [lra|]. + now change 3 with (Z2R 3); apply Z2R_le. } +assert (P1 : (0 < 1)%Z) by omega. +unfold double_round_eq. +apply double_round_lt_mid. +- exact Vfexp1. +- exact Vfexp2. +- lra. +- now rewrite Lxy. +- rewrite Lxy. + assert (fexp1 (ln_beta x) < ln_beta x)%Z; [|omega]. + now apply ln_beta_generic_gt; [|apply Rgt_not_eq|]. +- unfold midp. + apply (Rplus_lt_reg_r (- round beta fexp1 Zfloor (x + y))). + 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. + 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. + 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. + apply (Rmult_le_reg_r (bpow (- fexp1 (ln_beta x)))); + [now apply bpow_gt_0|]. + rewrite (Rmult_assoc (/ 2)). + rewrite Rmult_minus_distr_r. + bpow_simplify. + apply (Rle_trans _ _ _ Bpow3). + apply Rle_trans with (/ 2 * (2 / 3)); [lra|]. + apply Rmult_le_compat_l; [lra|]. + apply (Rplus_le_reg_r (- 1)); ring_simplify. + replace (_ - _) with (- (/ 3)) by lra. + apply Ropp_le_contravar. + now apply Rle_trans with (bpow (- 1)); [apply bpow_le; omega|]. +Qed. + +(* double_round_plus_beta_ge_3_aux{0,1} together *) +Lemma double_round_plus_beta_ge_3_aux2 : + (3 <= beta)%Z -> + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_plus_beta_ge_3_hyp fexp1 fexp2 -> + forall x y, + 0 < y -> y <= x -> + generic_format beta fexp1 x -> + generic_format beta fexp1 y -> + double_round_eq fexp1 fexp2 choice1 choice2 (x + y). +Proof. +intros Hbeta fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x y Py Hyx Fx Fy. +assert (Px := Rlt_le_trans 0 y x Py Hyx). +unfold double_round_eq. +destruct (Zle_or_lt (ln_beta y) (fexp1 (ln_beta x) - 1)) as [Hly|Hly]. +- (* ln_beta y <= fexp1 (ln_beta x) - 1 *) + now apply double_round_plus_beta_ge_3_aux1. +- (* fexp1 (ln_beta x) - 1 < ln_beta y *) + rewrite (round_generic beta fexp2). + + reflexivity. + + now apply valid_rnd_N. + + assert (Hf1 : (fexp1 (ln_beta x) <= ln_beta y)%Z); [omega|]. + now apply (double_round_plus_beta_ge_3_aux0 fexp1). +Qed. + +Lemma double_round_plus_beta_ge_3_aux : + (3 <= beta)%Z -> + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_plus_beta_ge_3_hyp fexp1 fexp2 -> + forall x y, + 0 <= x -> 0 <= y -> + generic_format beta fexp1 x -> + generic_format beta fexp1 y -> + double_round_eq fexp1 fexp2 choice1 choice2 (x + y). +Proof. +intros Hbeta fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x y Nnx Nny Fx Fy. +unfold double_round_eq. +destruct (Req_dec x 0) as [Zx|Nzx]. +- (* x = 0 *) + destruct Hexp as (_,(_,(_,Hexp4))). + rewrite Zx; rewrite Rplus_0_l. + rewrite (round_generic beta fexp2). + + reflexivity. + + now apply valid_rnd_N. + + apply (generic_inclusion_ln_beta beta fexp1). + now intros _; apply Hexp4; omega. + exact Fy. +- (* x <> 0 *) + destruct (Req_dec y 0) as [Zy|Nzy]. + + (* y = 0 *) + destruct Hexp as (_,(_,(_,Hexp4))). + rewrite Zy; rewrite Rplus_0_r. + rewrite (round_generic beta fexp2). + * reflexivity. + * now apply valid_rnd_N. + * apply (generic_inclusion_ln_beta beta fexp1). + now intros _; apply Hexp4; omega. + exact Fx. + + (* y <> 0 *) + assert (Px : 0 < x); [lra|]. + assert (Py : 0 < y); [lra|]. + destruct (Rlt_or_le x y) as [H|H]. + * (* x < y *) + apply Rlt_le in H. + rewrite Rplus_comm. + now apply double_round_plus_beta_ge_3_aux2. + * now apply double_round_plus_beta_ge_3_aux2. +Qed. + +(* fexp1 (ln_beta x) <= ln_beta y : + * substraction is exact in the largest precision (fexp2). *) +Lemma double_round_minus_beta_ge_3_aux0 : + forall (fexp1 fexp2 : Z -> Z), + double_round_plus_beta_ge_3_hyp fexp1 fexp2 -> + forall x y, + 0 < y -> y < x -> + (fexp1 (ln_beta x) <= ln_beta y)%Z -> + generic_format beta fexp1 x -> generic_format beta fexp1 y -> + generic_format beta fexp2 (x - y). +Proof. +intros fexp1 fexp2 Hexp x y Py Hyx Hln Fx Fy. +assert (Px := Rlt_trans 0 y x Py Hyx). +destruct Hexp as (Hexp1,(_,(Hexp3,Hexp4))). +assert (Lyx : (ln_beta y <= ln_beta x)%Z); + [now apply ln_beta_le; [|apply Rlt_le]|]. +destruct (Z.lt_ge_cases (ln_beta x - 2) (ln_beta y)) as [Hlt|Hge]. +- (* ln_beta x - 2 < ln_beta y *) + assert (Hor : (ln_beta y = ln_beta x :> Z) + \/ (ln_beta y = ln_beta x - 1 :> Z)%Z); [omega|]. + destruct Hor as [Heq|Heqm1]. + + (* ln_beta y = ln_beta x *) + apply (double_round_minus_aux0_aux fexp1); [| |exact Fx|exact Fy]. + * apply Hexp4. + apply Zle_trans with (ln_beta (x - y)); [omega|]. + now apply ln_beta_minus. + * rewrite Heq. + apply Hexp4. + apply Zle_trans with (ln_beta (x - y)); [omega|]. + now apply ln_beta_minus. + + (* ln_beta y = ln_beta x - 1 *) + apply (double_round_minus_aux0_aux fexp1); [| |exact Fx|exact Fy]. + * apply Hexp4. + apply Zle_trans with (ln_beta (x - y)); [omega|]. + now apply ln_beta_minus. + * rewrite Heqm1. + apply Hexp4. + apply Zplus_le_compat_r. + now apply ln_beta_minus. +- (* ln_beta y <= ln_beta x - 2 *) + destruct (ln_beta_minus_disj x y Px Py Hge) as [Lxmy|Lxmy]. + + (* ln_beta (x - y) = ln_beta x *) + apply (double_round_minus_aux0_aux fexp1); [| |exact Fx|exact Fy]. + * apply Hexp4. + omega. + * now rewrite Lxmy; apply Hexp3. + + (* ln_beta (x - y) = ln_beta x - 1 *) + apply (double_round_minus_aux0_aux fexp1); [| |exact Fx|exact Fy]; + rewrite Lxmy. + * apply Hexp1. + replace (_ + _)%Z with (ln_beta x : Z); [|ring]. + now apply Zle_trans with (ln_beta y). + * apply Hexp1. + now replace (_ + _)%Z with (ln_beta x : Z); [|ring]. +Qed. + +(* ln_beta y <= fexp1 (ln_beta x) - 1, + * fexp1 (ln_beta (x - y)) <= ln_beta y : + * substraction is exact in the largest precision (fexp2). *) +Lemma double_round_minus_beta_ge_3_aux1 : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + double_round_plus_beta_ge_3_hyp fexp1 fexp2 -> + forall x y, + 0 < y -> y < x -> + (ln_beta y <= fexp1 (ln_beta x) - 1)%Z -> + (fexp1 (ln_beta (x - y)) <= ln_beta y)%Z -> + generic_format beta fexp1 x -> generic_format beta fexp1 y -> + generic_format beta fexp2 (x - y). +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 Hexp x y Py Hyx Hln Hln' Fx Fy. +assert (Px := Rlt_trans 0 y x Py Hyx). +destruct Hexp as (Hexp1,(Hexp2,(Hexp3,Hexp4))). +assert (Lyx : (ln_beta y <= ln_beta x)%Z); + [now apply ln_beta_le; [|apply Rlt_le]|]. +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|]|]. +apply (double_round_minus_aux0_aux fexp1); [| |exact Fx|exact Fy]. +- apply Zle_trans with (fexp1 (ln_beta (x - y))). + + apply Hexp4; omega. + + omega. +- now apply Hexp3. +Qed. + +(* ln_beta y <= fexp1 (ln_beta x) - 1 : + * ln_beta y <= fexp1 (ln_beta (x - y)) - 1 : + * double_round_gt_mid applies. *) +Lemma double_round_minus_beta_ge_3_aux2 : + (3 <= beta)%Z -> + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_plus_beta_ge_3_hyp fexp1 fexp2 -> + forall x y, + 0 < y -> y < x -> + (ln_beta y <= fexp1 (ln_beta x) - 1)%Z -> + (ln_beta y <= fexp1 (ln_beta (x - y)) - 1)%Z -> + generic_format beta fexp1 x -> + generic_format beta fexp1 y -> + double_round_eq fexp1 fexp2 choice1 choice2 (x - y). +Proof. +intros Hbeta fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x y Py Hxy Hly Hly' Fx Fy. +assert (Px := Rlt_trans 0 y x Py Hxy). +destruct Hexp as (_,(_,(_,Hexp4))). +assert (Hf2 : (fexp2 (ln_beta x) <= fexp1 (ln_beta x))%Z); + [now apply Hexp4; omega|]. +assert (Hfx : (fexp1 (ln_beta x) < ln_beta x)%Z); + [now apply ln_beta_generic_gt; [|apply Rgt_not_eq|]|]. +assert (Bpow3 : bpow (- 1) <= / 3). +{ unfold Fcore_Raux.bpow, Z.pow_pos; simpl. + rewrite Zmult_1_r. + apply Rinv_le; [lra|]. + now change 3 with (Z2R 3); apply Z2R_le. } +assert (Ly : y < bpow (ln_beta y)). +{ apply Rabs_lt_inv. + apply bpow_ln_beta_gt. } +unfold double_round_eq. +apply double_round_gt_mid. +- exact Vfexp1. +- exact Vfexp2. +- lra. +- apply Hexp4; omega. +- assert (fexp1 (ln_beta (x - y)) < ln_beta (x - y))%Z; [|omega]. + apply (valid_exp_large fexp1 (ln_beta x - 1)). + + apply (valid_exp_large fexp1 (ln_beta y)); [|omega]. + now apply ln_beta_generic_gt; [|apply Rgt_not_eq|]. + + now apply ln_beta_minus_lb; [| |omega]. +- unfold midp'. + apply (Rplus_lt_reg_r (/ 2 * ulp beta fexp1 (x - y) - (x - y))). + ring_simplify. + replace (_ + _) with (round beta fexp1 Zceil (x - y) - (x - y)) by ring. + apply Rlt_le_trans with (bpow (fexp1 (ln_beta (x - y)) - 1)). + + apply Rle_lt_trans with y; + [now apply double_round_minus_aux2_aux|]. + apply (Rlt_le_trans _ _ _ Ly). + now apply bpow_le. + + unfold ulp, canonic_exp. + unfold Zminus at 1; rewrite bpow_plus. + rewrite Rmult_comm. + apply Rmult_le_compat_r; [now apply bpow_ge_0|]. + 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; omega. +- intro Hf2'. + unfold midp'. + apply (Rplus_lt_reg_r (/ 2 * (ulp beta fexp1 (x - y) + - ulp beta fexp2 (x - y)) - (x - y))). + ring_simplify; rewrite <- Rmult_minus_distr_l. + replace (_ + _) with (round beta fexp1 Zceil (x - y) - (x - y)) by ring. + apply Rle_lt_trans with y; + [now apply double_round_minus_aux2_aux|]. + apply (Rlt_le_trans _ _ _ Ly). + apply Rle_trans with (bpow (fexp1 (ln_beta (x - y)) - 1)); + [now apply bpow_le|]. + unfold ulp, canonic_exp. + apply (Rmult_le_reg_r (bpow (- fexp1 (ln_beta (x - y))))); + [now apply bpow_gt_0|]. + rewrite Rmult_assoc. + rewrite Rmult_minus_distr_r. + bpow_simplify. + apply Rle_trans with (/ 3). + + unfold Fcore_Raux.bpow, Z.pow_pos; simpl. + rewrite Zmult_1_r; apply Rinv_le; [lra|]. + now change 3 with (Z2R 3); apply Z2R_le. + + replace (/ 3) with (/ 2 * (2 / 3)) by field. + apply Rmult_le_compat_l; [lra|]. + apply (Rplus_le_reg_r (- 1)); ring_simplify. + replace (_ - _) with (- / 3) by field. + apply Ropp_le_contravar. + apply Rle_trans with (bpow (- 1)). + * apply bpow_le; omega. + * unfold Fcore_Raux.bpow, Z.pow_pos; simpl. + rewrite Zmult_1_r; apply Rinv_le; [lra|]. + now change 3 with (Z2R 3); apply Z2R_le. +Qed. + +(* double_round_minus_aux{0,1,2} together *) +Lemma double_round_minus_beta_ge_3_aux3 : + (3 <= beta)%Z -> + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_plus_beta_ge_3_hyp fexp1 fexp2 -> + forall x y, + 0 < y -> y <= x -> + generic_format beta fexp1 x -> + generic_format beta fexp1 y -> + double_round_eq fexp1 fexp2 choice1 choice2 (x - y). +Proof. +intros Hbeta fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x y Py Hyx Fx Fy. +assert (Px := Rlt_le_trans 0 y x Py Hyx). +unfold double_round_eq. +destruct (Req_dec y x) as [Hy|Hy]. +- (* y = x *) + rewrite Hy; replace (x - x) with 0 by ring. + rewrite round_0. + + reflexivity. + + now apply valid_rnd_N. +- (* y < x *) + assert (Hyx' : y < x); [lra|]. + destruct (Zle_or_lt (ln_beta y) (fexp1 (ln_beta x) - 1)) as [Hly|Hly]. + + (* ln_beta y <= fexp1 (ln_beta x) - 1 *) + destruct (Zle_or_lt (ln_beta y) (fexp1 (ln_beta (x - y)) - 1)) + as [Hly'|Hly']. + * (* ln_beta y <= fexp1 (ln_beta (x - y)) - 1 *) + now apply double_round_minus_beta_ge_3_aux2. + * (* fexp1 (ln_beta (x - y)) - 1 < ln_beta y *) + { rewrite (round_generic beta fexp2). + - reflexivity. + - now apply valid_rnd_N. + - assert (Hf1 : (fexp1 (ln_beta (x - y)) <= ln_beta y)%Z); [omega|]. + now apply (double_round_minus_beta_ge_3_aux1 fexp1). } + + rewrite (round_generic beta fexp2). + * reflexivity. + * now apply valid_rnd_N. + * assert (Hf1 : (fexp1 (ln_beta x) <= ln_beta y)%Z); [omega|]. + now apply (double_round_minus_beta_ge_3_aux0 fexp1). +Qed. + +Lemma double_round_minus_beta_ge_3_aux : + (3 <= beta)%Z -> + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_plus_beta_ge_3_hyp fexp1 fexp2 -> + forall x y, + 0 <= x -> 0 <= y -> + generic_format beta fexp1 x -> + generic_format beta fexp1 y -> + double_round_eq fexp1 fexp2 choice1 choice2 (x - y). +Proof. +intros Hbeta fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x y Nnx Nny Fx Fy. +unfold double_round_eq. +destruct (Req_dec x 0) as [Zx|Nzx]. +- (* x = 0 *) + rewrite Zx; unfold Rminus; rewrite Rplus_0_l. + do 3 rewrite round_N_opp. + rewrite (round_generic beta fexp2). + * reflexivity. + * now apply valid_rnd_N. + * apply (generic_inclusion_ln_beta beta fexp1). + destruct Hexp as (_,(_,(_,Hexp4))). + now intros _; apply Hexp4; omega. + exact Fy. +- (* x <> 0 *) + destruct (Req_dec y 0) as [Zy|Nzy]. + + (* y = 0 *) + rewrite Zy; unfold Rminus; rewrite Ropp_0; rewrite Rplus_0_r. + rewrite (round_generic beta fexp2). + * reflexivity. + * now apply valid_rnd_N. + * apply (generic_inclusion_ln_beta beta fexp1). + destruct Hexp as (_,(_,(_,Hexp4))). + now intros _; apply Hexp4; omega. + exact Fx. + + (* y <> 0 *) + assert (Px : 0 < x); [lra|]. + assert (Py : 0 < y); [lra|]. + destruct (Rlt_or_le x y) as [H|H]. + * (* x < y *) + apply Rlt_le in H. + replace (x - y) with (- (y - x)) by ring. + do 3 rewrite round_N_opp. + apply Ropp_eq_compat. + now apply double_round_minus_beta_ge_3_aux3. + * (* y <= x *) + now apply double_round_minus_beta_ge_3_aux3. +Qed. + +Lemma double_round_plus_beta_ge_3 : + (3 <= beta)%Z -> + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_plus_beta_ge_3_hyp fexp1 fexp2 -> + forall x y, + generic_format beta fexp1 x -> + generic_format beta fexp1 y -> + double_round_eq fexp1 fexp2 choice1 choice2 (x + y). +Proof. +intros Hbeta fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x y Fx Fy. +unfold double_round_eq. +destruct (Rlt_or_le x 0) as [Sx|Sx]; destruct (Rlt_or_le y 0) as [Sy|Sy]. +- (* x < 0, y < 0 *) + replace (x + y) with (- (- x - y)); [|ring]. + do 3 rewrite round_N_opp. + apply Ropp_eq_compat. + assert (Px : 0 <= - x); [lra|]. + assert (Py : 0 <= - y); [lra|]. + apply generic_format_opp in Fx. + apply generic_format_opp in Fy. + now apply double_round_plus_beta_ge_3_aux. +- (* x < 0, 0 <= y *) + replace (x + y) with (y - (- x)); [|ring]. + assert (Px : 0 <= - x); [lra|]. + apply generic_format_opp in Fx. + now apply double_round_minus_beta_ge_3_aux. +- (* 0 <= x, y < 0 *) + replace (x + y) with (x - (- y)); [|ring]. + assert (Py : 0 <= - y); [lra|]. + apply generic_format_opp in Fy. + now apply double_round_minus_beta_ge_3_aux. +- (* 0 <= x, 0 <= y *) + now apply double_round_plus_beta_ge_3_aux. +Qed. + +Lemma double_round_minus_beta_ge_3 : + (3 <= beta)%Z -> + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_plus_beta_ge_3_hyp fexp1 fexp2 -> + forall x y, + generic_format beta fexp1 x -> + generic_format beta fexp1 y -> + double_round_eq fexp1 fexp2 choice1 choice2 (x - y). +Proof. +intros Hbeta fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x y Fx Fy. +unfold Rminus. +apply generic_format_opp in Fy. +now apply double_round_plus_beta_ge_3. +Qed. + +Section Double_round_plus_beta_ge_3_FLX. + +Require Import Fcore_FLX. + +Variable prec : Z. +Variable prec' : Z. + +Context { prec_gt_0_ : Prec_gt_0 prec }. +Context { prec_gt_0_' : Prec_gt_0 prec' }. + +Lemma FLX_double_round_plus_beta_ge_3_hyp : + (2 * prec <= prec')%Z -> + double_round_plus_beta_ge_3_hyp (FLX_exp prec) (FLX_exp prec'). +Proof. +intros Hprec. +unfold FLX_exp. +unfold double_round_plus_beta_ge_3_hyp; split; [|split; [|split]]; +intros ex ey; try omega. +unfold Prec_gt_0 in prec_gt_0_. +omega. +Qed. + +Theorem double_round_plus_beta_ge_3_FLX : + (3 <= beta)%Z -> + forall choice1 choice2, + (2 * prec <= prec')%Z -> + forall x y, + FLX_format beta prec x -> FLX_format beta prec y -> + double_round_eq (FLX_exp prec) (FLX_exp prec') choice1 choice2 (x + y). +Proof. +intros Hbeta choice1 choice2 Hprec x y Fx Fy. +apply double_round_plus_beta_ge_3. +- exact Hbeta. +- now apply FLX_exp_valid. +- now apply FLX_exp_valid. +- now apply FLX_double_round_plus_beta_ge_3_hyp. +- now apply generic_format_FLX. +- now apply generic_format_FLX. +Qed. + +Theorem double_round_minus_beta_ge_3_FLX : + (3 <= beta)%Z -> + forall choice1 choice2, + (2 * prec <= prec')%Z -> + forall x y, + FLX_format beta prec x -> FLX_format beta prec y -> + double_round_eq (FLX_exp prec) (FLX_exp prec') choice1 choice2 (x - y). +Proof. +intros Hbeta choice1 choice2 Hprec x y Fx Fy. +apply double_round_minus_beta_ge_3. +- exact Hbeta. +- now apply FLX_exp_valid. +- now apply FLX_exp_valid. +- now apply FLX_double_round_plus_beta_ge_3_hyp. +- now apply generic_format_FLX. +- now apply generic_format_FLX. +Qed. + +End Double_round_plus_beta_ge_3_FLX. + +Section Double_round_plus_beta_ge_3_FLT. + +Require Import Fcore_FLX. +Require Import Fcore_FLT. + +Variable emin prec : Z. +Variable emin' prec' : Z. + +Context { prec_gt_0_ : Prec_gt_0 prec }. +Context { prec_gt_0_' : Prec_gt_0 prec' }. + +Lemma FLT_double_round_plus_beta_ge_3_hyp : + (emin' <= emin)%Z -> (2 * prec <= prec')%Z -> + double_round_plus_beta_ge_3_hyp (FLT_exp emin prec) (FLT_exp emin' prec'). +Proof. +intros Hemin Hprec. +unfold FLT_exp. +unfold double_round_plus_beta_ge_3_hyp; split; [|split; [|split]]; intros ex ey. +- generalize (Zmax_spec (ex + 1 - prec) emin). + generalize (Zmax_spec (ex - prec') emin'). + generalize (Zmax_spec (ey - prec) emin). + omega. +- generalize (Zmax_spec (ex - 1 - prec) emin). + generalize (Zmax_spec (ex - prec') emin'). + generalize (Zmax_spec (ey - prec) emin). + omega. +- generalize (Zmax_spec (ex - prec) emin). + generalize (Zmax_spec (ex - prec') emin'). + generalize (Zmax_spec (ey - prec) emin). + omega. +- unfold Prec_gt_0 in prec_gt_0_. + generalize (Zmax_spec (ex - prec') emin'). + generalize (Zmax_spec (ey - prec) emin). + omega. +Qed. + +Theorem double_round_plus_beta_ge_3_FLT : + (3 <= beta)%Z -> + forall choice1 choice2, + (emin' <= emin)%Z -> (2 * prec <= prec')%Z -> + forall x y, + FLT_format beta emin prec x -> FLT_format beta emin prec y -> + double_round_eq (FLT_exp emin prec) (FLT_exp emin' prec') + choice1 choice2 (x + y). +Proof. +intros Hbeta choice1 choice2 Hemin Hprec x y Fx Fy. +apply double_round_plus_beta_ge_3. +- exact Hbeta. +- now apply FLT_exp_valid. +- now apply FLT_exp_valid. +- now apply FLT_double_round_plus_beta_ge_3_hyp. +- now apply generic_format_FLT. +- now apply generic_format_FLT. +Qed. + +Theorem double_round_minus_beta_ge_3_FLT : + (3 <= beta)%Z -> + forall choice1 choice2, + (emin' <= emin)%Z -> (2 * prec <= prec')%Z -> + forall x y, + FLT_format beta emin prec x -> FLT_format beta emin prec y -> + double_round_eq (FLT_exp emin prec) (FLT_exp emin' prec') + choice1 choice2 (x - y). +Proof. +intros Hbeta choice1 choice2 Hemin Hprec x y Fx Fy. +apply double_round_minus_beta_ge_3. +- exact Hbeta. +- now apply FLT_exp_valid. +- now apply FLT_exp_valid. +- now apply FLT_double_round_plus_beta_ge_3_hyp. +- now apply generic_format_FLT. +- now apply generic_format_FLT. +Qed. + +End Double_round_plus_beta_ge_3_FLT. + +Section Double_round_plus_beta_ge_3_FTZ. + +Require Import Fcore_FLX. +Require Import Fcore_FTZ. + +Variable emin prec : Z. +Variable emin' prec' : Z. + +Context { prec_gt_0_ : Prec_gt_0 prec }. +Context { prec_gt_0_' : Prec_gt_0 prec' }. + +Lemma FTZ_double_round_plus_beta_ge_3_hyp : + (emin' + prec' <= emin + 1)%Z -> (2 * prec <= prec')%Z -> + double_round_plus_beta_ge_3_hyp (FTZ_exp emin prec) (FTZ_exp emin' prec'). +Proof. +intros Hemin Hprec. +unfold FTZ_exp. +unfold Prec_gt_0 in *. +unfold double_round_plus_beta_ge_3_hyp; split; [|split; [|split]]; intros ex ey. +- destruct (Z.ltb_spec (ex + 1 - prec) emin); + destruct (Z.ltb_spec (ex - prec') emin'); + destruct (Z.ltb_spec (ey - prec) emin); + omega. +- destruct (Z.ltb_spec (ex - 1 - prec) emin); + destruct (Z.ltb_spec (ex - prec') emin'); + destruct (Z.ltb_spec (ey - prec) emin); + omega. +- destruct (Z.ltb_spec (ex - prec) emin); + destruct (Z.ltb_spec (ex - prec') emin'); + destruct (Z.ltb_spec (ey - prec) emin); + omega. +- destruct (Z.ltb_spec (ex - prec') emin'); + destruct (Z.ltb_spec (ey - prec) emin); + omega. +Qed. + +Theorem double_round_plus_beta_ge_3_FTZ : + (3 <= beta)%Z -> + forall choice1 choice2, + (emin' + prec' <= emin + 1)%Z -> (2 * prec <= prec')%Z -> + forall x y, + FTZ_format beta emin prec x -> FTZ_format beta emin prec y -> + double_round_eq (FTZ_exp emin prec) (FTZ_exp emin' prec') + choice1 choice2 (x + y). +Proof. +intros Hbeta choice1 choice2 Hemin Hprec x y Fx Fy. +apply double_round_plus_beta_ge_3. +- exact Hbeta. +- now apply FTZ_exp_valid. +- now apply FTZ_exp_valid. +- now apply FTZ_double_round_plus_beta_ge_3_hyp. +- now apply generic_format_FTZ. +- now apply generic_format_FTZ. +Qed. + +Theorem double_round_minus_beta_ge_3_FTZ : + (3 <= beta)%Z -> + forall choice1 choice2, + (emin' + prec' <= emin + 1)%Z -> (2 * prec <= prec')%Z -> + forall x y, + FTZ_format beta emin prec x -> FTZ_format beta emin prec y -> + double_round_eq (FTZ_exp emin prec) (FTZ_exp emin' prec') + choice1 choice2 (x - y). +Proof. +intros Hbeta choice1 choice2 Hemin Hprec x y Fx Fy. +apply double_round_minus_beta_ge_3. +- exact Hbeta. +- now apply FTZ_exp_valid. +- now apply FTZ_exp_valid. +- now apply FTZ_double_round_plus_beta_ge_3_hyp. +- now apply generic_format_FTZ. +- now apply generic_format_FTZ. +Qed. + +End Double_round_plus_beta_ge_3_FTZ. + +End Double_round_plus_beta_ge_3. + +End Double_round_plus. + +Lemma double_round_mid_cases : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + forall x, + 0 < x -> + (fexp2 (ln_beta x) <= fexp1 (ln_beta x) - 1)%Z -> + (fexp1 (ln_beta x) <= ln_beta x)%Z -> + (Rabs (x - midp fexp1 x) <= / 2 * (ulp beta fexp2 x) -> + double_round_eq fexp1 fexp2 choice1 choice2 x) -> + double_round_eq fexp1 fexp2 choice1 choice2 x. +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 x Px Hf2f1 Hf1. +unfold double_round_eq, midp. +set (rd := round beta fexp1 Zfloor x). +set (u1 := ulp beta fexp1 x). +set (u2 := ulp beta fexp2 x). +intros Cmid. +destruct (generic_format_EM beta fexp1 x) as [Fx|Nfx]. +- (* generic_format beta fexp1 x *) + rewrite (round_generic beta fexp2); [reflexivity|now apply valid_rnd_N|]. + 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|]. + 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) *) + apply double_round_lt_mid_further_place; try assumption. + unfold midp. fold rd; fold u1; fold u2. + apply (Rplus_lt_reg_r (- rd)); ring_simplify. + now rewrite <- Rmult_minus_distr_l. + + (* / 2 * (u1 - u2) <= x - rd *) + { destruct (Rlt_or_le (/ 2 * (u1 + u2)) (x - rd)). + - (* / 2 * (u1 + u2) < x - rd *) + assert (round beta fexp1 Zceil x - x + < / 2 * (ulp beta fexp1 x - ulp beta fexp2 x)). + { rewrite Hceil; fold u1; fold u2. + lra. } + apply double_round_gt_mid_further_place; try assumption. + unfold midp'; lra. + - (* x - rd <= / 2 * (u1 + u2) *) + apply Cmid, Rabs_le; split; lra. } +Qed. + +Section Double_round_sqrt. + +Definition double_round_sqrt_hyp fexp1 fexp2 := + (forall ex, (2 * fexp1 ex <= fexp1 (2 * ex))%Z) + /\ (forall ex, (2 * fexp1 ex <= fexp1 (2 * ex - 1))%Z) + /\ (forall ex, (fexp1 (2 * ex) < 2 * ex)%Z -> + (fexp2 ex + ex <= 2 * fexp1 ex - 2)%Z). + +Lemma ln_beta_sqrt_disj : + forall x, + 0 < x -> + (ln_beta x = 2 * ln_beta (sqrt x) - 1 :> Z)%Z + \/ (ln_beta x = 2 * ln_beta (sqrt x) :> Z)%Z. +Proof. +intros x Px. +generalize (ln_beta_sqrt beta x Px). +intro H. +omega. +Qed. + +Lemma double_round_sqrt_aux : + forall fexp1 fexp2 : Z -> Z, + Valid_exp fexp1 -> Valid_exp fexp2 -> + double_round_sqrt_hyp fexp1 fexp2 -> + forall x, + 0 < x -> + (fexp2 (ln_beta (sqrt x)) <= fexp1 (ln_beta (sqrt x)) - 1)%Z -> + generic_format beta fexp1 x -> + / 2 * ulp beta fexp2 (sqrt x) < Rabs (sqrt x - midp fexp1 (sqrt x)). +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 := ulp beta fexp1 (sqrt x)). +set (u2 := ulp beta fexp2 (sqrt x)). +set (b := / 2 * (u1 - u2)). +set (b' := / 2 * (u1 + u2)). +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, canonic_exp; simpl. +set (mx := Ztrunc (x * bpow (- fexp1 (ln_beta x)))). +set (ma := Ztrunc (a * bpow (- fexp1 (ln_beta 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, ulp, canonic_exp. + apply bpow_lt. + omega. } +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 (ln_beta (sqrt x)) <= fexp1 (ln_beta (x)))%Z); + [destruct (ln_beta_sqrt_disj x Px) as [H'|H']; rewrite H'; apply Hexp|]. +assert (Hlx : (fexp1 (2 * ln_beta (sqrt x)) < 2 * ln_beta (sqrt x))%Z). +{ destruct (ln_beta_sqrt_disj x Px) as [Hlx|Hlx]. + - apply (valid_exp_large fexp1 (ln_beta x)); [|omega]. + now apply ln_beta_generic_gt; [|apply Rgt_not_eq|]. + - rewrite <- Hlx. + now apply ln_beta_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]. +- (* a = 0 *) + 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 (ln_beta 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. + change 0 with (Z2R 0); apply f_equal. + apply Zfloor_imp. + split; [now apply Rmult_le_pos; [apply Rlt_le|apply bpow_ge_0]|simpl]. + apply (Rmult_lt_reg_r (bpow (fexp1 (ln_beta x)))); [now apply bpow_gt_0|]. + rewrite Rmult_1_l; bpow_simplify. + apply Rlt_le_trans with (bpow (2 * fexp1 (ln_beta (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 (ln_beta (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, canonic_exp; apply bpow_lt; omega. +- (* a <> 0 *) + assert (Pa : 0 < a); [lra|]. + assert (Hla : (ln_beta a = ln_beta (sqrt x) :> Z)). + { unfold a; apply ln_beta_round_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 (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). + + apply Rplus_lt_compat_l. + rewrite <- (Rmult_1_l (ulp _ _ _)). + apply Rmult_lt_compat_r; [apply bpow_gt_0|lra]. + + apply (succ_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. + unfold u1, u2, ulp, canonic_exp; bpow_simplify; apply bpow_le. + now apply Hexp. + + apply Rmult_le_compat. + * apply bpow_ge_0. + * apply pow2_ge_0. + * unfold Fcore_Raux.bpow, Z.pow_pos; simpl; rewrite Zmult_1_r. + apply Rinv_le; [lra|]. + change 4 with (Z2R (2 * 2)%Z); apply Z2R_le, Zmult_le_compat; omega. + * 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, canonic_exp. + rewrite <- (Rmult_1_l (bpow _)); rewrite Fa; rewrite <- Rmult_plus_distr_r. + rewrite <- Rmult_assoc; rewrite (Rmult_comm _ (Z2R ma)). + rewrite <- (Rmult_assoc (Z2R ma)); bpow_simplify. + apply (Rmult_le_reg_r (bpow (- 2 * fexp1 (ln_beta (sqrt x))))); + [now apply bpow_gt_0|bpow_simplify]. + rewrite Fx at 1; bpow_simplify. + rewrite <- Z2R_Zpower; [|omega]. + change 1 with (Z2R 1); rewrite <- Z2R_plus; do 2 rewrite <- Z2R_mult. + apply Z2R_le, Zlt_succ_le, lt_Z2R. + unfold Z.succ; rewrite Z2R_plus; do 2 rewrite Z2R_mult; rewrite Z2R_plus. + rewrite Z2R_Zpower; [|omega]. + apply (Rmult_lt_reg_r (bpow (2 * fexp1 (ln_beta (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, canonic_exp; 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 (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. + + 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_ln_beta_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, canonic_exp. + bpow_simplify. + apply bpow_le. + rewrite Zplus_comm. + now apply Hexp. + * apply Rmult_le_compat_r; [now apply pow2_ge_0|]. + unfold Fcore_Raux.bpow; simpl; unfold Z.pow_pos; simpl. + rewrite Zmult_1_r. + apply Rinv_le; [lra|]. + change 2 with (Z2R 2); apply Z2R_le. + rewrite <- (Zmult_1_l 2). + apply Zmult_le_compat; omega. + + 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, canonic_exp; apply bpow_lt; omega. } + 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. + +(* --> 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, + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_sqrt_hyp fexp1 fexp2 -> + forall x, + generic_format beta fexp1 x -> + double_round_eq fexp1 fexp2 choice1 choice2 (sqrt x). +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x Fx. +unfold double_round_eq. +destruct (Rle_or_lt x 0) as [Npx|Px]. +- (* x <= 0 *) + rewrite (sqrt_neg _ Npx). + now rewrite round_0; [|apply valid_rnd_N]. +- (* 0 < x *) + assert (Hfx : (fexp1 (ln_beta x) < ln_beta x)%Z); + [now apply ln_beta_generic_gt; try assumption; lra|]. + assert (Hfsx : (fexp1 (ln_beta (sqrt x)) < ln_beta (sqrt x))%Z). + { destruct (Rle_or_lt x 1) as [Hx|Hx]. + - (* x <= 1 *) + apply (valid_exp_large fexp1 (ln_beta x)); [exact Hfx|]. + apply ln_beta_le; [exact Px|]. + rewrite <- (sqrt_def x) at 1; [|lra]. + rewrite <- Rmult_1_r. + apply Rmult_le_compat_l. + + apply sqrt_pos. + + rewrite <- sqrt_1. + now apply sqrt_le_1_alt. + - (* 1 < x *) + generalize ((proj1 (proj2 Hexp)) 1%Z). + replace (_ - 1)%Z with 1%Z by ring. + intro Hexp10. + assert (Hf0 : (fexp1 1 < 1)%Z); [omega|clear Hexp10]. + apply (valid_exp_large fexp1 1); [exact Hf0|]. + apply ln_beta_ge_bpow. + rewrite Zeq_minus; [|reflexivity]. + unfold Fcore_Raux.bpow; simpl. + apply Rabs_ge; right. + rewrite <- sqrt_1. + apply sqrt_le_1_alt. + now apply Rlt_le. } + assert (Hf2 : (fexp2 (ln_beta (sqrt x)) <= fexp1 (ln_beta (sqrt x)) - 1)%Z). + { assert (H : (fexp1 (2 * ln_beta (sqrt x)) < 2 * ln_beta (sqrt x))%Z). + { destruct (ln_beta_sqrt_disj x Px) as [Hlx|Hlx]. + - apply (valid_exp_large fexp1 (ln_beta x)); [|omega]. + now apply ln_beta_generic_gt; [|apply Rgt_not_eq|]. + - rewrite <- Hlx. + now apply ln_beta_generic_gt; [|apply Rgt_not_eq|]. } + generalize ((proj2 (proj2 Hexp)) (ln_beta (sqrt x)) H). + omega. } + apply double_round_mid_cases. + + exact Vfexp1. + + exact Vfexp2. + + now apply sqrt_lt_R0. + + omega. + + omega. + + intros Hmid; casetype False; apply (Rle_not_lt _ _ Hmid). + apply (double_round_sqrt_aux fexp1 fexp2 Vfexp1 Vfexp2 Hexp x Px Hf2 Fx). +Qed. + +Section Double_round_sqrt_FLX. + +Require Import Fcore_FLX. + +Variable prec : Z. +Variable prec' : Z. + +Context { prec_gt_0_ : Prec_gt_0 prec }. +Context { prec_gt_0_' : Prec_gt_0 prec' }. + +Lemma FLX_double_round_sqrt_hyp : + (2 * prec + 2 <= prec')%Z -> + double_round_sqrt_hyp (FLX_exp prec) (FLX_exp prec'). +Proof. +intros Hprec. +unfold FLX_exp. +unfold Prec_gt_0 in prec_gt_0_. +unfold double_round_sqrt_hyp; split; [|split]; intro ex; omega. +Qed. + +Theorem double_round_sqrt_FLX : + forall choice1 choice2, + (2 * prec + 2 <= prec')%Z -> + forall x, + FLX_format beta prec x -> + double_round_eq (FLX_exp prec) (FLX_exp prec') choice1 choice2 (sqrt x). +Proof. +intros choice1 choice2 Hprec x Fx. +apply double_round_sqrt. +- now apply FLX_exp_valid. +- now apply FLX_exp_valid. +- now apply FLX_double_round_sqrt_hyp. +- now apply generic_format_FLX. +Qed. + +End Double_round_sqrt_FLX. + +Section Double_round_sqrt_FLT. + +Require Import Fcore_FLX. +Require Import Fcore_FLT. + +Variable emin prec : Z. +Variable emin' prec' : Z. + +Context { prec_gt_0_ : Prec_gt_0 prec }. +Context { prec_gt_0_' : Prec_gt_0 prec' }. + +Lemma FLT_double_round_sqrt_hyp : + (emin <= 0)%Z -> + ((emin' <= emin - prec - 2)%Z + \/ (2 * emin' <= emin - 4 * prec - 2)%Z) -> + (2 * prec + 2 <= prec')%Z -> + double_round_sqrt_hyp (FLT_exp emin prec) (FLT_exp emin' prec'). +Proof. +intros Hemin Heminprec Hprec. +unfold FLT_exp. +unfold Prec_gt_0 in prec_gt_0_. +unfold double_round_sqrt_hyp; split; [|split]; intros ex. +- generalize (Zmax_spec (ex - prec) emin). + generalize (Zmax_spec (2 * ex - prec) emin). + omega. +- generalize (Zmax_spec (ex - prec) emin). + generalize (Zmax_spec (2 * ex - 1 - prec) emin). + omega. +- generalize (Zmax_spec (2 * ex - prec) emin). + generalize (Zmax_spec (ex - prec') emin'). + generalize (Zmax_spec (ex - prec) emin). + omega. +Qed. + +Theorem double_round_sqrt_FLT : + forall choice1 choice2, + (emin <= 0)%Z -> + ((emin' <= emin - prec - 2)%Z + \/ (2 * emin' <= emin - 4 * prec - 2)%Z) -> + (2 * prec + 2 <= prec')%Z -> + forall x, + FLT_format beta emin prec x -> + double_round_eq (FLT_exp emin prec) (FLT_exp emin' prec') + choice1 choice2 (sqrt x). +Proof. +intros choice1 choice2 Hemin Heminprec Hprec x Fx. +apply double_round_sqrt. +- now apply FLT_exp_valid. +- now apply FLT_exp_valid. +- now apply FLT_double_round_sqrt_hyp. +- now apply generic_format_FLT. +Qed. + +End Double_round_sqrt_FLT. + +Section Double_round_sqrt_FTZ. + +Require Import Fcore_FLX. +Require Import Fcore_FTZ. + +Variable emin prec : Z. +Variable emin' prec' : Z. + +Context { prec_gt_0_ : Prec_gt_0 prec }. +Context { prec_gt_0_' : Prec_gt_0 prec' }. + +Lemma FTZ_double_round_sqrt_hyp : + (2 * (emin' + prec') <= emin + prec <= 1)%Z -> + (2 * prec + 2 <= prec')%Z -> + double_round_sqrt_hyp (FTZ_exp emin prec) (FTZ_exp emin' prec'). +Proof. +intros Hemin Hprec. +unfold FTZ_exp. +unfold Prec_gt_0 in *. +unfold double_round_sqrt_hyp; split; [|split]; intros ex. +- destruct (Z.ltb_spec (ex - prec) emin); + destruct (Z.ltb_spec (2 * ex - prec) emin); + omega. +- destruct (Z.ltb_spec (ex - prec) emin); + destruct (Z.ltb_spec (2 * ex - 1 - prec) emin); + omega. +- intro H. + destruct (Zle_or_lt emin (2 * ex - prec)) as [H'|H']. + + destruct (Z.ltb_spec (ex - prec') emin'); + destruct (Z.ltb_spec (ex - prec) emin); + omega. + + casetype False. + rewrite (Zlt_bool_true _ _ H') in H. + omega. +Qed. + +Theorem double_round_sqrt_FTZ : + (4 <= beta)%Z -> + forall choice1 choice2, + (2 * (emin' + prec') <= emin + prec <= 1)%Z -> + (2 * prec + 2 <= prec')%Z -> + forall x, + FTZ_format beta emin prec x -> + double_round_eq (FTZ_exp emin prec) (FTZ_exp emin' prec') + choice1 choice2 (sqrt x). +Proof. +intros Hbeta choice1 choice2 Hemin Hprec x Fx. +apply double_round_sqrt. +- now apply FTZ_exp_valid. +- now apply FTZ_exp_valid. +- now apply FTZ_double_round_sqrt_hyp. +- now apply generic_format_FTZ. +Qed. + +End Double_round_sqrt_FTZ. + +Section Double_round_sqrt_beta_ge_4. + +Definition double_round_sqrt_beta_ge_4_hyp fexp1 fexp2 := + (forall ex, (2 * fexp1 ex <= fexp1 (2 * ex))%Z) + /\ (forall ex, (2 * fexp1 ex <= fexp1 (2 * ex - 1))%Z) + /\ (forall ex, (fexp1 (2 * ex) < 2 * ex)%Z -> + (fexp2 ex + ex <= 2 * fexp1 ex - 1)%Z). + +Lemma double_round_sqrt_beta_ge_4_aux : + (4 <= beta)%Z -> + forall fexp1 fexp2 : Z -> Z, + Valid_exp fexp1 -> Valid_exp fexp2 -> + double_round_sqrt_beta_ge_4_hyp fexp1 fexp2 -> + forall x, + 0 < x -> + (fexp2 (ln_beta (sqrt x)) <= fexp1 (ln_beta (sqrt x)) - 1)%Z -> + generic_format beta fexp1 x -> + / 2 * ulp beta fexp2 (sqrt x) < Rabs (sqrt x - midp fexp1 (sqrt x)). +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 (b := / 2 * (u1 - u2)). +set (b' := / 2 * (u1 + u2)). +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, canonic_exp; simpl. +set (mx := Ztrunc (x * bpow (- fexp1 (ln_beta x)))). +set (ma := Ztrunc (a * bpow (- fexp1 (ln_beta 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, ulp, canonic_exp. + apply bpow_lt. + omega. } +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 (ln_beta (sqrt x)) <= fexp1 (ln_beta (x)))%Z); + [destruct (ln_beta_sqrt_disj x Px) as [H'|H']; rewrite H'; apply Hexp|]. +assert (Hlx : (fexp1 (2 * ln_beta (sqrt x)) < 2 * ln_beta (sqrt x))%Z). +{ destruct (ln_beta_sqrt_disj x Px) as [Hlx|Hlx]. + - apply (valid_exp_large fexp1 (ln_beta x)); [|omega]. + now apply ln_beta_generic_gt; [|apply Rgt_not_eq|]. + - rewrite <- Hlx. + now apply ln_beta_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]. +- (* a = 0 *) + 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 (ln_beta 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. + change 0 with (Z2R 0); apply f_equal. + apply Zfloor_imp. + split; [now apply Rmult_le_pos; [apply Rlt_le|apply bpow_ge_0]|simpl]. + apply (Rmult_lt_reg_r (bpow (fexp1 (ln_beta x)))); [now apply bpow_gt_0|]. + rewrite Rmult_1_l; bpow_simplify. + apply Rlt_le_trans with (bpow (2 * fexp1 (ln_beta (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 (ln_beta (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, canonic_exp; apply bpow_lt; omega. +- (* a <> 0 *) + assert (Pa : 0 < a); [lra|]. + assert (Hla : (ln_beta a = ln_beta (sqrt x) :> Z)). + { unfold a; apply ln_beta_round_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 (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). + + apply Rplus_lt_compat_l. + rewrite <- (Rmult_1_l (ulp _ _ _)). + apply Rmult_lt_compat_r; [apply bpow_gt_0|lra]. + + apply (succ_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. + unfold u1, u2, ulp, canonic_exp; bpow_simplify; apply bpow_le. + now apply Hexp. + + apply Rmult_le_compat. + * apply bpow_ge_0. + * apply pow2_ge_0. + * unfold Fcore_Raux.bpow, Z.pow_pos; simpl; rewrite Zmult_1_r. + apply Rinv_le; [lra|]. + now change 4 with (Z2R 4); apply Z2R_le. + * 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, canonic_exp. + rewrite <- (Rmult_1_l (bpow _)); rewrite Fa; rewrite <- Rmult_plus_distr_r. + rewrite <- Rmult_assoc; rewrite (Rmult_comm _ (Z2R ma)). + rewrite <- (Rmult_assoc (Z2R ma)); bpow_simplify. + apply (Rmult_le_reg_r (bpow (- 2 * fexp1 (ln_beta (sqrt x))))); + [now apply bpow_gt_0|bpow_simplify]. + rewrite Fx at 1; bpow_simplify. + rewrite <- Z2R_Zpower; [|omega]. + change 1 with (Z2R 1); rewrite <- Z2R_plus; do 2 rewrite <- Z2R_mult. + apply Z2R_le, Zlt_succ_le, lt_Z2R. + unfold Z.succ; rewrite Z2R_plus; do 2 rewrite Z2R_mult; rewrite Z2R_plus. + rewrite Z2R_Zpower; [|omega]. + apply (Rmult_lt_reg_r (bpow (2 * fexp1 (ln_beta (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, canonic_exp; 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 (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. + + 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_ln_beta_gt. + - apply Rle_trans with (/ 2 * u1 ^ 2). + + apply Rle_trans with (bpow (- 1) * u1 ^ 2). + * unfold pow; rewrite Rmult_1_r. + unfold u2, u1, ulp, canonic_exp. + bpow_simplify. + apply bpow_le. + rewrite Zplus_comm. + now apply Hexp. + * apply Rmult_le_compat_r; [now apply pow2_ge_0|]. + unfold Fcore_Raux.bpow; simpl; unfold Z.pow_pos; simpl. + rewrite Zmult_1_r. + apply Rinv_le; [lra|]. + change 2 with (Z2R 2); apply Z2R_le; omega. + + 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, canonic_exp; apply bpow_lt; omega. } + 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. + +Lemma double_round_sqrt_beta_ge_4 : + (4 <= beta)%Z -> + forall fexp1 fexp2 : Z -> Z, + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_sqrt_beta_ge_4_hyp fexp1 fexp2 -> + forall x, + generic_format beta fexp1 x -> + double_round_eq fexp1 fexp2 choice1 choice2 (sqrt x). +Proof. +intros Hbeta fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x Fx. +unfold double_round_eq. +destruct (Rle_or_lt x 0) as [Npx|Px]. +- (* x <= 0 *) + assert (Hs : sqrt x = 0). + { 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; lra. } + rewrite Hs. + rewrite round_0. + + reflexivity. + + now apply valid_rnd_N. +- (* 0 < x *) + assert (Hfx : (fexp1 (ln_beta x) < ln_beta x)%Z); + [now apply ln_beta_generic_gt; try assumption; lra|]. + assert (Hfsx : (fexp1 (ln_beta (sqrt x)) < ln_beta (sqrt x))%Z). + { destruct (Rle_or_lt x 1) as [Hx|Hx]. + - (* x <= 1 *) + apply (valid_exp_large fexp1 (ln_beta x)); [exact Hfx|]. + apply ln_beta_le; [exact Px|]. + rewrite <- (sqrt_def x) at 1; [|lra]. + rewrite <- Rmult_1_r. + apply Rmult_le_compat_l. + + apply sqrt_pos. + + rewrite <- sqrt_1. + now apply sqrt_le_1_alt. + - (* 1 < x *) + generalize ((proj1 (proj2 Hexp)) 1%Z). + replace (_ - 1)%Z with 1%Z by ring. + intro Hexp10. + assert (Hf0 : (fexp1 1 < 1)%Z); [omega|clear Hexp10]. + apply (valid_exp_large fexp1 1); [exact Hf0|]. + apply ln_beta_ge_bpow. + rewrite Zeq_minus; [|reflexivity]. + unfold Fcore_Raux.bpow; simpl. + apply Rabs_ge; right. + rewrite <- sqrt_1. + apply sqrt_le_1_alt. + now apply Rlt_le. } + assert (Hf2 : (fexp2 (ln_beta (sqrt x)) <= fexp1 (ln_beta (sqrt x)) - 1)%Z). + { assert (H : (fexp1 (2 * ln_beta (sqrt x)) < 2 * ln_beta (sqrt x))%Z). + { destruct (ln_beta_sqrt_disj x Px) as [Hlx|Hlx]. + - apply (valid_exp_large fexp1 (ln_beta x)); [|omega]. + now apply ln_beta_generic_gt; [|apply Rgt_not_eq|]. + - rewrite <- Hlx. + now apply ln_beta_generic_gt; [|apply Rgt_not_eq|]. } + generalize ((proj2 (proj2 Hexp)) (ln_beta (sqrt x)) H). + omega. } + apply double_round_mid_cases. + + exact Vfexp1. + + exact Vfexp2. + + now apply sqrt_lt_R0. + + omega. + + omega. + + intros Hmid; casetype False; apply (Rle_not_lt _ _ Hmid). + apply (double_round_sqrt_beta_ge_4_aux Hbeta fexp1 fexp2 Vfexp1 Vfexp2 + Hexp x Px Hf2 Fx). +Qed. + +Section Double_round_sqrt_beta_ge_4_FLX. + +Require Import Fcore_FLX. + +Variable prec : Z. +Variable prec' : Z. + +Context { prec_gt_0_ : Prec_gt_0 prec }. +Context { prec_gt_0_' : Prec_gt_0 prec' }. + +Lemma FLX_double_round_sqrt_beta_ge_4_hyp : + (2 * prec + 1 <= prec')%Z -> + double_round_sqrt_beta_ge_4_hyp (FLX_exp prec) (FLX_exp prec'). +Proof. +intros Hprec. +unfold FLX_exp. +unfold Prec_gt_0 in prec_gt_0_. +unfold double_round_sqrt_beta_ge_4_hyp; split; [|split]; intro ex; omega. +Qed. + +Theorem double_round_sqrt_beta_ge_4_FLX : + (4 <= beta)%Z -> + forall choice1 choice2, + (2 * prec + 1 <= prec')%Z -> + forall x, + FLX_format beta prec x -> + double_round_eq (FLX_exp prec) (FLX_exp prec') choice1 choice2 (sqrt x). +Proof. +intros Hbeta choice1 choice2 Hprec x Fx. +apply double_round_sqrt_beta_ge_4. +- exact Hbeta. +- now apply FLX_exp_valid. +- now apply FLX_exp_valid. +- now apply FLX_double_round_sqrt_beta_ge_4_hyp. +- now apply generic_format_FLX. +Qed. + +End Double_round_sqrt_beta_ge_4_FLX. + +Section Double_round_sqrt_beta_ge_4_FLT. + +Require Import Fcore_FLX. +Require Import Fcore_FLT. + +Variable emin prec : Z. +Variable emin' prec' : Z. + +Context { prec_gt_0_ : Prec_gt_0 prec }. +Context { prec_gt_0_' : Prec_gt_0 prec' }. + +Lemma FLT_double_round_sqrt_beta_ge_4_hyp : + (emin <= 0)%Z -> + ((emin' <= emin - prec - 1)%Z + \/ (2 * emin' <= emin - 4 * prec)%Z) -> + (2 * prec + 1 <= prec')%Z -> + double_round_sqrt_beta_ge_4_hyp (FLT_exp emin prec) (FLT_exp emin' prec'). +Proof. +intros Hemin Heminprec Hprec. +unfold FLT_exp. +unfold Prec_gt_0 in prec_gt_0_. +unfold double_round_sqrt_beta_ge_4_hyp; split; [|split]; intros ex. +- generalize (Zmax_spec (ex - prec) emin). + generalize (Zmax_spec (2 * ex - prec) emin). + omega. +- generalize (Zmax_spec (ex - prec) emin). + generalize (Zmax_spec (2 * ex - 1 - prec) emin). + omega. +- generalize (Zmax_spec (2 * ex - prec) emin). + generalize (Zmax_spec (ex - prec') emin'). + generalize (Zmax_spec (ex - prec) emin). + omega. +Qed. + +Theorem double_round_sqrt_beta_ge_4_FLT : + (4 <= beta)%Z -> + forall choice1 choice2, + (emin <= 0)%Z -> + ((emin' <= emin - prec - 1)%Z + \/ (2 * emin' <= emin - 4 * prec)%Z) -> + (2 * prec + 1 <= prec')%Z -> + forall x, + FLT_format beta emin prec x -> + double_round_eq (FLT_exp emin prec) (FLT_exp emin' prec') + choice1 choice2 (sqrt x). +Proof. +intros Hbeta choice1 choice2 Hemin Heminprec Hprec x Fx. +apply double_round_sqrt_beta_ge_4. +- exact Hbeta. +- now apply FLT_exp_valid. +- now apply FLT_exp_valid. +- now apply FLT_double_round_sqrt_beta_ge_4_hyp. +- now apply generic_format_FLT. +Qed. + +End Double_round_sqrt_beta_ge_4_FLT. + +Section Double_round_sqrt_beta_ge_4_FTZ. + +Require Import Fcore_FLX. +Require Import Fcore_FTZ. + +Variable emin prec : Z. +Variable emin' prec' : Z. + +Context { prec_gt_0_ : Prec_gt_0 prec }. +Context { prec_gt_0_' : Prec_gt_0 prec' }. + +Lemma FTZ_double_round_sqrt_beta_ge_4_hyp : + (2 * (emin' + prec') <= emin + prec <= 1)%Z -> + (2 * prec + 1 <= prec')%Z -> + double_round_sqrt_beta_ge_4_hyp (FTZ_exp emin prec) (FTZ_exp emin' prec'). +Proof. +intros Hemin Hprec. +unfold FTZ_exp. +unfold Prec_gt_0 in *. +unfold double_round_sqrt_beta_ge_4_hyp; split; [|split]; intros ex. +- destruct (Z.ltb_spec (ex - prec) emin); + destruct (Z.ltb_spec (2 * ex - prec) emin); + omega. +- destruct (Z.ltb_spec (ex - prec) emin); + destruct (Z.ltb_spec (2 * ex - 1 - prec) emin); + omega. +- intro H. + destruct (Zle_or_lt emin (2 * ex - prec)) as [H'|H']. + + destruct (Z.ltb_spec (ex - prec') emin'); + destruct (Z.ltb_spec (ex - prec) emin); + omega. + + casetype False. + rewrite (Zlt_bool_true _ _ H') in H. + omega. +Qed. + +Theorem double_round_sqrt_beta_ge_4_FTZ : + (4 <= beta)%Z -> + forall choice1 choice2, + (2 * (emin' + prec') <= emin + prec <= 1)%Z -> + (2 * prec + 1 <= prec')%Z -> + forall x, + FTZ_format beta emin prec x -> + double_round_eq (FTZ_exp emin prec) (FTZ_exp emin' prec') + choice1 choice2 (sqrt x). +Proof. +intros Hbeta choice1 choice2 Hemin Hprec x Fx. +apply double_round_sqrt_beta_ge_4. +- exact Hbeta. +- now apply FTZ_exp_valid. +- now apply FTZ_exp_valid. +- now apply FTZ_double_round_sqrt_beta_ge_4_hyp. +- now apply generic_format_FTZ. +Qed. + +End Double_round_sqrt_beta_ge_4_FTZ. + +End Double_round_sqrt_beta_ge_4. + +End Double_round_sqrt. + +Section Double_round_div. + +Lemma double_round_eq_mid_beta_even : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + (exists n, (beta = 2 * n :> Z)%Z) -> + forall x, + 0 < x -> + (fexp2 (ln_beta x) <= fexp1 (ln_beta x) - 1)%Z -> + (fexp1 (ln_beta x) <= ln_beta x)%Z -> + x = midp fexp1 x -> + double_round_eq fexp1 fexp2 choice1 choice2 x. +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Ebeta x Px Hf2 Hf1. +unfold double_round_eq. +unfold midp. +set (rd := round beta fexp1 Zfloor x). +set (u := ulp beta fexp1 x). +intro H; apply (Rplus_eq_compat_l (- rd)) in H. +ring_simplify in H; revert H. +rewrite Rplus_comm; fold (Rminus x rd). +intro Xmid. +destruct Ebeta as (n,Ebeta). +assert (Hbeta : (2 <= beta)%Z). +{ destruct beta as (beta_val,beta_prop). + now apply Zle_bool_imp_le. } +apply (Rplus_eq_compat_l rd) in Xmid; ring_simplify in Xmid. +rewrite (round_generic beta fexp2); [reflexivity|now apply valid_rnd_N|]. +set (f := Float beta (Zfloor (scaled_mantissa beta fexp2 rd) + + n * beta ^ (fexp1 (ln_beta x) - 1 + - fexp2 (ln_beta x))) + (canonic_exp beta fexp2 x)). +assert (Hf : F2R f = x). +{ unfold f, F2R; simpl. + rewrite Z2R_plus. + rewrite Rmult_plus_distr_r. + rewrite Z2R_mult. + rewrite Z2R_Zpower; [|omega]. + unfold canonic_exp at 2; bpow_simplify. + unfold Zminus; rewrite bpow_plus. + rewrite (Rmult_comm _ (bpow (- 1))). + rewrite <- (Rmult_assoc (Z2R n)). + change (bpow (- 1)) with (/ Z2R (beta * 1)). + rewrite Zmult_1_r. + rewrite Ebeta. + rewrite (Z2R_mult 2). + rewrite Rinv_mult_distr; + [|simpl; lra|change 0 with (Z2R 0); apply Z2R_neq; omega]. + rewrite <- Rmult_assoc; rewrite (Rmult_comm (Z2R n)); + 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. + apply f_equal2; [|reflexivity]. + destruct (Req_dec rd 0) as [Zrd|Nzrd]. + - (* rd = 0 *) + rewrite Zrd. + rewrite scaled_mantissa_0. + change 0 with (Z2R 0) at 1; rewrite Zfloor_Z2R. + now rewrite Rmult_0_l. + - (* rd <> 0 *) + assert (Nnrd : 0 <= rd). + { apply round_DN_pt. + - exact Vfexp1. + - apply generic_format_0. + - now apply Rlt_le. } + assert (Prd : 0 < rd); [lra|]. + assert (Lrd : (ln_beta rd = ln_beta x :> Z)). + { apply Zle_antisym. + - apply ln_beta_le; [exact Prd|]. + now apply round_DN_pt. + - apply ln_beta_round_ge. + + exact Vfexp1. + + now apply valid_rnd_DN. + + exact Nzrd. } + unfold scaled_mantissa. + unfold rd at 1. + unfold round, F2R, scaled_mantissa, canonic_exp; simpl. + bpow_simplify. + rewrite Lrd. + rewrite <- (Z2R_Zpower _ (_ - _)); [|omega]. + rewrite <- Z2R_mult. + rewrite (Zfloor_imp (Zfloor (x * bpow (- fexp1 (ln_beta x))) * + beta ^ (fexp1 (ln_beta x) - fexp2 (ln_beta x)))). + + rewrite Z2R_mult. + rewrite Z2R_Zpower; [|omega]. + bpow_simplify. + now unfold rd. + + split; [now apply Rle_refl|]. + rewrite Z2R_plus. + simpl; lra. } +apply (generic_format_F2R' _ _ x f Hf). +intros _. +apply Zle_refl. +Qed. + +Lemma double_round_really_zero : + forall (fexp1 fexp2 : Z -> Z), + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + forall x, + 0 < x -> + (ln_beta x <= fexp1 (ln_beta x) - 2)%Z -> + double_round_eq fexp1 fexp2 choice1 choice2 x. +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 x Px Hf1. +assert (Hlx : bpow (ln_beta x - 1) <= x < bpow (ln_beta x)). +{ destruct (ln_beta x) as (ex,Hex); simpl. + rewrite <- (Rabs_right x); [|now apply Rle_ge; apply Rlt_le]. + apply Hex. + now apply Rgt_not_eq. } +unfold double_round_eq. +rewrite (round_N_really_small_pos beta fexp1 _ x (ln_beta x)); [|exact Hlx|omega]. +set (x'' := round beta fexp2 (Znearest choice2) x). +destruct (Req_dec x'' 0) as [Zx''|Nzx'']; + [now rewrite Zx''; rewrite round_0; [|apply valid_rnd_N]|]. +destruct (Zle_or_lt (fexp2 (ln_beta x)) (ln_beta x)). +- (* fexp2 (ln_beta x) <= ln_beta x *) + destruct (Rlt_or_le x'' (bpow (ln_beta x))). + + (* x'' < bpow (ln_beta x) *) + rewrite (round_N_really_small_pos beta fexp1 _ _ (ln_beta x)); + [reflexivity|split; [|exact H0]|omega]. + apply round_large_pos_ge_pow; [now apply valid_rnd_N| |now apply Hlx]. + fold x''; assert (0 <= x''); [|lra]; unfold x''. + rewrite <- (round_0 beta fexp2 (Znearest choice2)). + now apply round_le; [|apply valid_rnd_N|apply Rlt_le]. + + (* bpow (ln_beta x) <= x'' *) + assert (Hx'' : x'' = bpow (ln_beta x)). + { apply Rle_antisym; [|exact H0]. + rewrite <- (round_generic beta fexp2 (Znearest choice2) (bpow _)). + - now apply round_le; [|apply valid_rnd_N|apply Rlt_le]. + - now apply generic_format_bpow'. } + rewrite Hx''. + unfold round, F2R, scaled_mantissa, canonic_exp; simpl. + rewrite ln_beta_bpow. + assert (Hf11 : (fexp1 (ln_beta x + 1) = fexp1 (ln_beta x) :> Z)%Z); + [apply Vfexp1; omega|]. + rewrite Hf11. + apply (Rmult_eq_reg_r (bpow (- fexp1 (ln_beta x)))); + [|now apply Rgt_not_eq; apply bpow_gt_0]. + rewrite Rmult_0_l; bpow_simplify. + change 0 with (Z2R 0); apply f_equal. + apply Znearest_imp. + simpl; unfold Rminus; rewrite Ropp_0; rewrite Rplus_0_r. + rewrite Rabs_right; [|now apply Rle_ge; apply bpow_ge_0]. + apply Rle_lt_trans with (bpow (- 2)); [now apply bpow_le; omega|]. + unfold Fcore_Raux.bpow, Z.pow_pos; simpl; rewrite Zmult_1_r. + assert (Hbeta : (2 <= beta)%Z). + { destruct beta as (beta_val,beta_prop); simpl. + now apply Zle_bool_imp_le. } + apply Rinv_lt_contravar. + * apply Rmult_lt_0_compat; [lra|]. + rewrite Z2R_mult; apply Rmult_lt_0_compat; change 0 with (Z2R 0); + apply Z2R_lt; omega. + * change 2 with (Z2R 2); apply Z2R_lt. + apply (Zle_lt_trans _ _ _ Hbeta). + rewrite <- (Zmult_1_r beta) at 1. + apply Zmult_lt_compat_l; omega. +- (* ln_beta x < fexp2 (ln_beta x) *) + casetype False; apply Nzx''. + now apply (round_N_really_small_pos beta _ _ _ (ln_beta x)). +Qed. + +Lemma double_round_zero : + forall fexp1 fexp2 : Z -> Z, + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + forall x, + 0 < x -> + (fexp1 (ln_beta x) = ln_beta x + 1 :> Z)%Z -> + x < bpow (ln_beta x) - / 2 * ulp beta fexp2 x -> + double_round_eq fexp1 fexp2 choice1 choice2 x. +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 x Px Hf1. +unfold double_round_eq. +set (x'' := round beta fexp2 (Znearest choice2) x). +set (u1 := ulp beta fexp1 x). +set (u2 := ulp beta fexp2 x). +intro Hx. +assert (Hlx : bpow (ln_beta x - 1) <= x < bpow (ln_beta x)). +{ destruct (ln_beta x) as (ex,Hex); simpl. + rewrite <- (Rabs_right x); [|now apply Rle_ge; apply Rlt_le]. + apply Hex. + now apply Rgt_not_eq. } +rewrite (round_N_really_small_pos beta fexp1 choice1 x (ln_beta x)); + [|exact Hlx|omega]. +destruct (Req_dec x'' 0) as [Zx''|Nzx'']; + [now rewrite Zx''; rewrite round_0; [reflexivity|apply valid_rnd_N]|]. +rewrite (round_N_really_small_pos beta _ _ x'' (ln_beta x)); + [reflexivity| |omega]. +split. +- apply round_large_pos_ge_pow. + + now apply valid_rnd_N. + + assert (0 <= x''); [|now fold x''; lra]. + rewrite <- (round_0 beta fexp2 (Znearest choice2)). + now apply round_le; [|apply valid_rnd_N|apply Rlt_le]. + + apply Rle_trans with (Rabs x); + [|now rewrite Rabs_right; [apply Rle_refl|apply Rle_ge; apply Rlt_le]]. + destruct (ln_beta x) as (ex,Hex); simpl; apply Hex. + now apply Rgt_not_eq. +- replace x'' with (x + (x'' - x)) by ring. + 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. +Qed. + +Lemma double_round_all_mid_cases : + forall fexp1 fexp2 : Z -> Z, + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + forall x, + 0 < x -> + (fexp2 (ln_beta x) <= fexp1 (ln_beta x) - 1)%Z -> + ((fexp1 (ln_beta x) = ln_beta x + 1 :> Z)%Z -> + bpow (ln_beta x) - / 2 * ulp beta fexp2 x <= x -> + double_round_eq fexp1 fexp2 choice1 choice2 x) -> + ((fexp1 (ln_beta x) <= ln_beta x)%Z -> + midp fexp1 x - / 2 * ulp beta fexp2 x <= x < midp fexp1 x -> + double_round_eq fexp1 fexp2 choice1 choice2 x) -> + ((fexp1 (ln_beta x) <= ln_beta x)%Z -> + x = midp fexp1 x -> + double_round_eq fexp1 fexp2 choice1 choice2 x) -> + ((fexp1 (ln_beta x) <= ln_beta x)%Z -> + midp fexp1 x < x <= midp fexp1 x + / 2 * ulp beta fexp2 x -> + double_round_eq fexp1 fexp2 choice1 choice2 x) -> + double_round_eq fexp1 fexp2 choice1 choice2 x. +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 x Px Hf2. +set (x' := round beta fexp1 Zfloor x). +set (u1 := ulp beta fexp1 x). +set (u2 := ulp beta fexp2 x). +intros Cz Clt Ceq Cgt. +destruct (Ztrichotomy (ln_beta x) (fexp1 (ln_beta x) - 1)) as [Hlt|[Heq|Hgt]]. +- (* ln_beta x < fexp1 (ln_beta x) - 1 *) + assert (H : (ln_beta x <= fexp1 (ln_beta x) - 2)%Z) by omega. + now apply double_round_really_zero. +- (* ln_beta x = fexp1 (ln_beta x) - 1 *) + assert (H : (fexp1 (ln_beta x) = (ln_beta x + 1))%Z) by omega. + destruct (Rlt_or_le x (bpow (ln_beta x) - / 2 * u2)) as [Hlt'|Hge']. + + now apply double_round_zero. + + now apply Cz. +- (* ln_beta x > fexp1 (ln_beta x) - 1 *) + assert (H : (fexp1 (ln_beta x) <= ln_beta x)%Z) by omega. + destruct (Rtotal_order x (midp fexp1 x)) as [Hlt'|[Heq'|Hgt']]. + + (* x < midp fexp1 x *) + destruct (Rlt_or_le x (midp fexp1 x - / 2 * u2)) as [Hlt''|Hle'']. + * now apply double_round_lt_mid_further_place; [| | |omega| |]. + * now apply Clt; [|split]. + + (* x = midp fexp1 x *) + now apply Ceq. + + (* x > midp fexp1 x *) + destruct (Rle_or_lt x (midp fexp1 x + / 2 * u2)) as [Hlt''|Hle'']. + * now apply Cgt; [|split]. + * { destruct (generic_format_EM beta fexp1 x) as [Fx|Nfx]. + - (* generic_format beta fexp1 x *) + unfold double_round_eq; rewrite (round_generic beta fexp2); + [reflexivity|now apply valid_rnd_N|]. + 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|]. + assert (Hf2' : (fexp2 (ln_beta x) <= fexp1 (ln_beta x) - 1)%Z); + [omega|]. + assert (midp' fexp1 x + / 2 * ulp beta fexp2 x < x); + [|now apply double_round_gt_mid_further_place]. + revert Hle''; unfold midp, midp'; fold x'. + rewrite Hceil; fold u1; fold u2. + lra. } +Qed. + +Lemma ln_beta_div_disj : + forall x y : R, + 0 < x -> 0 < y -> + ((ln_beta (x / y) = ln_beta x - ln_beta y :> Z)%Z + \/ (ln_beta (x / y) = ln_beta x - ln_beta y + 1 :> Z)%Z). +Proof. +intros x y Px Py. +generalize (ln_beta_div beta x y Px Py). +omega. +Qed. + +Definition double_round_div_hyp fexp1 fexp2 := + (forall ex, (fexp2 ex <= fexp1 ex - 1)%Z) + /\ (forall ex ey, (fexp1 ex < ex)%Z -> (fexp1 ey < ey)%Z -> + (fexp1 (ex - ey) <= ex - ey + 1)%Z -> + (fexp2 (ex - ey) <= fexp1 ex - ey)%Z) + /\ (forall ex ey, (fexp1 ex < ex)%Z -> (fexp1 ey < ey)%Z -> + (fexp1 (ex - ey + 1) <= ex - ey + 1 + 1)%Z -> + (fexp2 (ex - ey + 1) <= fexp1 ex - ey)%Z) + /\ (forall ex ey, (fexp1 ex < ex)%Z -> (fexp1 ey < ey)%Z -> + (fexp1 (ex - ey) <= ex - ey)%Z -> + (fexp2 (ex - ey) <= fexp1 (ex - ey) + + fexp1 ey - ey)%Z) + /\ (forall ex ey, (fexp1 ex < ex)%Z -> (fexp1 ey < ey)%Z -> + (fexp1 (ex - ey) = ex - ey + 1)%Z -> + (fexp2 (ex - ey) <= ex - ey - ey + fexp1 ey)%Z). + +Lemma double_round_div_aux0 : + forall fexp1 fexp2 : Z -> Z, + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_div_hyp fexp1 fexp2 -> + forall x y, + 0 < x -> 0 < y -> + generic_format beta fexp1 x -> + generic_format beta fexp1 y -> + fexp1 (ln_beta (x / y)) = (ln_beta (x / y) + 1)%Z -> + ~ (bpow (ln_beta (x / y)) - / 2 * ulp beta fexp2 (x / y) <= x / y). +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x y Px Py Fx Fy Hf1. +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|]|]. +set (p := bpow (ln_beta (x / y))). +set (u2 := ulp beta fexp2 (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. +intro Hl. +assert (Hr : x / y < p); + [now apply Rabs_lt_inv; apply bpow_ln_beta_gt|]. +apply (Rlt_irrefl (p - / 2 * u2)). +apply (Rle_lt_trans _ _ _ Hl). +apply (Rmult_lt_reg_r y _ _ Py). +unfold Rdiv; rewrite Rmult_assoc. +rewrite Rinv_l; [|now apply Rgt_not_eq]; rewrite Rmult_1_r. +destruct (Zle_or_lt Z0 (fexp1 (ln_beta x) - ln_beta (x / y) + - fexp1 (ln_beta y))%Z) as [He|He]. +- (* ln_beta (x / y) + fexp1 (ln_beta y) <= fexp1 (ln_beta x) *) + apply Rle_lt_trans with (p * y - p * bpow (fexp1 (ln_beta y))). + + rewrite Fx; rewrite Fy at 1. + rewrite <- Rmult_assoc. + rewrite (Rmult_comm p). + unfold p; bpow_simplify. + apply (Rmult_le_reg_r (bpow (- ln_beta (x / y) - fexp1 (ln_beta y)))); + [now apply bpow_gt_0|]. + rewrite Rmult_minus_distr_r. + bpow_simplify. + rewrite <- Z2R_Zpower; [|exact He]. + rewrite <- Z2R_mult. + change 1 with (Z2R 1); rewrite <- Z2R_minus. + apply Z2R_le. + apply (Zplus_le_reg_r _ _ 1); ring_simplify. + apply Zlt_le_succ. + apply lt_Z2R. + rewrite Z2R_mult. + rewrite Z2R_Zpower; [|exact He]. + apply (Rmult_lt_reg_r (bpow (fexp1 (ln_beta y) + ln_beta (x / y)))); + [now apply bpow_gt_0|]. + bpow_simplify. + rewrite <- Fx. + rewrite bpow_plus. + rewrite <- Rmult_assoc; rewrite <- Fy. + fold p. + apply (Rmult_lt_reg_r (/ y)); [now apply Rinv_0_lt_compat|]. + field_simplify; lra. + + rewrite Rmult_minus_distr_r. + unfold Rminus; apply Rplus_lt_compat_l. + apply Ropp_lt_contravar. + apply Rlt_le_trans with (u2 * bpow (ln_beta y)). + * rewrite <- (Rmult_1_l (u2 * _)). + rewrite Rmult_assoc. + { apply Rmult_lt_compat. + - lra. + - now apply Rmult_le_pos; [apply bpow_ge_0|apply Rlt_le]. + - lra. + - apply Rmult_lt_compat_l; [now apply bpow_gt_0|]. + apply Rabs_lt_inv. + apply bpow_ln_beta_gt. } + * unfold u2, p, ulp, canonic_exp; bpow_simplify; apply bpow_le. + apply (Zplus_le_reg_r _ _ (- ln_beta y)); ring_simplify. + rewrite (Zplus_comm (- _)); fold (Zminus (ln_beta (x / y)) (ln_beta y)). + destruct (ln_beta_div_disj x y Px Py) as [Hxy|Hxy]; rewrite Hxy; + [now apply Hexp; [| |rewrite <- Hxy]|]. + replace (_ - _ + 1)%Z with ((ln_beta x + 1) - ln_beta y)%Z by ring. + apply Hexp. + { now assert (fexp1 (ln_beta x + 1) <= ln_beta x)%Z; + [apply valid_exp|omega]. } + { assumption. } + replace (_ + 1 - _)%Z with (ln_beta x - ln_beta y + 1)%Z by ring. + now rewrite <- Hxy. +- (* fexp1 (ln_beta x) < ln_beta (x / y) + fexp1 (ln_beta y) *) + apply Rle_lt_trans with (p * y - bpow (fexp1 (ln_beta x))). + + rewrite Fx at 1; rewrite Fy at 1. + apply (Rmult_le_reg_r (bpow (- fexp1 (ln_beta x)))); + [now apply bpow_gt_0|]. + rewrite Rmult_minus_distr_r. + bpow_simplify. + rewrite (Rmult_comm p). + unfold p; bpow_simplify. + rewrite <- Z2R_Zpower; [|omega]. + rewrite <- Z2R_mult. + change 1 with (Z2R 1); rewrite <- Z2R_minus. + apply Z2R_le. + apply (Zplus_le_reg_r _ _ 1); ring_simplify. + apply Zlt_le_succ. + apply lt_Z2R. + rewrite Z2R_mult. + rewrite Z2R_Zpower; [|omega]. + apply (Rmult_lt_reg_r (bpow (fexp1 (ln_beta x)))); + [now apply bpow_gt_0|bpow_simplify]. + rewrite <- Fx. + rewrite Zplus_comm; rewrite bpow_plus. + rewrite <- Rmult_assoc; rewrite <- Fy. + fold p. + apply (Rmult_lt_reg_r (/ y)); [now apply Rinv_0_lt_compat|]. + field_simplify; lra. + + rewrite Rmult_minus_distr_r. + unfold Rminus; apply Rplus_lt_compat_l. + apply Ropp_lt_contravar. + apply Rlt_le_trans with (u2 * bpow (ln_beta y)). + * rewrite <- (Rmult_1_l (u2 * _)). + rewrite Rmult_assoc. + { apply Rmult_lt_compat. + - lra. + - now apply Rmult_le_pos; [apply bpow_ge_0|apply Rlt_le]. + - lra. + - apply Rmult_lt_compat_l; [now apply bpow_gt_0|]. + apply Rabs_lt_inv. + apply bpow_ln_beta_gt. } + * unfold u2, p, ulp, canonic_exp; bpow_simplify; apply bpow_le. + apply (Zplus_le_reg_r _ _ (- ln_beta y)); ring_simplify. + rewrite (Zplus_comm (- _)); fold (Zminus (ln_beta (x / y)) (ln_beta y)). + destruct (ln_beta_div_disj x y Px Py) as [Hxy|Hxy]; rewrite Hxy; + apply Hexp; try assumption; rewrite <- Hxy; rewrite Hf1; apply Zle_refl. +Qed. + +Lemma double_round_div_aux1 : + forall fexp1 fexp2 : Z -> Z, + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_div_hyp fexp1 fexp2 -> + forall x y, + 0 < x -> 0 < y -> + generic_format beta fexp1 x -> + generic_format beta fexp1 y -> + (fexp1 (ln_beta (x / y)) <= ln_beta (x / y))%Z -> + ~ (midp fexp1 (x / y) - / 2 * ulp beta fexp2 (x / y) + <= x / y + < midp fexp1 (x / y)). +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x y Px Py Fx Fy Hf1. +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|]|]. +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))). +{ intro H; intro H'; apply H; split. + - apply (Rplus_le_reg_l (round beta fexp1 Zfloor (x / y))). + ring_simplify. + apply H'. + - 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 (x' := round beta fexp1 Zfloor (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. +intro Hlr. +apply (Rlt_irrefl (/ 2 * (u1 - u2))). +apply (Rle_lt_trans _ _ _ (proj1 Hlr)). +apply (Rplus_lt_reg_r x'); ring_simplify. +apply (Rmult_lt_reg_r y _ _ Py). +unfold Rdiv; rewrite Rmult_assoc. +rewrite Rinv_l; [|now apply Rgt_not_eq]; rewrite Rmult_1_r. +rewrite Rmult_minus_distr_r; rewrite Rmult_plus_distr_r. +apply (Rmult_lt_reg_l 2); [lra|]. +rewrite Rmult_minus_distr_l; rewrite Rmult_plus_distr_l. +do 5 rewrite <- Rmult_assoc. +rewrite Rinv_r; [|lra]; do 2 rewrite Rmult_1_l. +destruct (Zle_or_lt Z0 (fexp1 (ln_beta x) - fexp1 (ln_beta (x / y)) + - fexp1 (ln_beta y))%Z) as [He|He]. +- (* fexp1 (ln_beta (x / y)) + fexp1 (ln_beta y)) <= fexp1 (ln_beta x) *) + apply Rle_lt_trans with (2 * x' * y + u1 * y + - bpow (fexp1 (ln_beta (x / y)) + + fexp1 (ln_beta y))). + + rewrite Fx at 1; rewrite Fy at 1 2. + apply (Rmult_le_reg_r (bpow (- fexp1 (ln_beta (x / y)) + - fexp1 (ln_beta y)))); + [now apply bpow_gt_0|]. + rewrite Rmult_minus_distr_r; rewrite (Rmult_plus_distr_r (_ * _ * _)). + bpow_simplify. + replace (2 * x' * _ * _) + with (2 * Z2R my * x' * bpow (- fexp1 (ln_beta (x / y)))) by ring. + rewrite (Rmult_comm u1). + unfold x', u1, round, F2R, ulp, scaled_mantissa, canonic_exp; simpl. + bpow_simplify. + rewrite <- Z2R_Zpower; [|exact He]. + change 2 with (Z2R 2). + do 4 rewrite <- Z2R_mult. + rewrite <- Z2R_plus. + change 1 with (Z2R 1); rewrite <- Z2R_minus. + apply Z2R_le. + apply (Zplus_le_reg_r _ _ 1); ring_simplify. + apply Zlt_le_succ. + apply lt_Z2R. + rewrite Z2R_plus. + do 4 rewrite Z2R_mult; simpl. + rewrite Z2R_Zpower; [|exact He]. + apply (Rmult_lt_reg_r (bpow (fexp1 (ln_beta (x / y)) + + fexp1 (ln_beta y)))); + [now apply bpow_gt_0|bpow_simplify]. + rewrite Rmult_assoc. + rewrite <- Fx. + rewrite (Rmult_plus_distr_r _ _ (Fcore_Raux.bpow _ _)). + rewrite Rmult_assoc. + rewrite bpow_plus. + rewrite <- (Rmult_assoc (Z2R (Zfloor _))). + change (Z2R (Zfloor _) * _) with x'. + do 2 rewrite (Rmult_comm _ (bpow (fexp1 (ln_beta y)))). + rewrite Rmult_assoc. + do 2 rewrite <- (Rmult_assoc (Z2R my)). + rewrite <- Fy. + change (bpow _) with u1. + apply (Rmult_lt_reg_l (/ 2)); [lra|]. + rewrite Rmult_plus_distr_l. + do 4 rewrite <- Rmult_assoc. + rewrite Rinv_l; [|lra]; do 2 rewrite Rmult_1_l. + apply (Rplus_lt_reg_r (- y * x')); ring_simplify. + apply (Rmult_lt_reg_l (/ y)); [now apply Rinv_0_lt_compat|]. + rewrite Rmult_minus_distr_l. + do 3 rewrite <- Rmult_assoc. + rewrite Rinv_l; [|now apply Rgt_not_eq]; do 2 rewrite Rmult_1_l. + now rewrite Rmult_comm. + + apply Rplus_lt_compat_l. + apply Ropp_lt_contravar. + apply Rlt_le_trans with (u2 * bpow (ln_beta y)). + * { apply Rmult_lt_compat_l. + - apply bpow_gt_0. + - apply Rabs_lt_inv. + apply bpow_ln_beta_gt. } + * unfold u2, ulp, canonic_exp; bpow_simplify; apply bpow_le. + apply (Zplus_le_reg_r _ _ (- ln_beta y)); ring_simplify. + rewrite <- Zplus_assoc; rewrite (Zplus_comm (- _)). + destruct (ln_beta_div_disj x y Px Py) as [Hxy|Hxy]; rewrite Hxy; + [now apply Hexp; [| |rewrite <- Hxy]|]. + replace (_ - _ + 1)%Z with ((ln_beta x + 1) - ln_beta y)%Z by ring. + apply Hexp. + { now assert (fexp1 (ln_beta x + 1) <= ln_beta x)%Z; + [apply valid_exp|omega]. } + { assumption. } + replace (_ + 1 - _)%Z with (ln_beta x - ln_beta y + 1)%Z by ring. + now rewrite <- Hxy. +- (* fexp1 (ln_beta x) < fexp1 (ln_beta (x / y)) + fexp1 (ln_beta y) *) + apply Rle_lt_trans with (2 * x' * y + u1 * y - bpow (fexp1 (ln_beta x))). + + rewrite Fx at 1; rewrite Fy at 1 2. + apply (Rmult_le_reg_r (bpow (- fexp1 (ln_beta x)))); + [now apply bpow_gt_0|]. + rewrite Rmult_minus_distr_r; rewrite (Rmult_plus_distr_r (_ * _ * _)). + bpow_simplify. + replace (2 * x' * _ * _) + with (2 * Z2R my * x' * bpow (fexp1 (ln_beta y) - fexp1 (ln_beta x))) by ring. + rewrite (Rmult_comm u1). + unfold x', u1, round, F2R, ulp, scaled_mantissa, canonic_exp; simpl. + bpow_simplify. + rewrite <- (Z2R_Zpower _ (_ - _)%Z); [|omega]. + change 2 with (Z2R 2). + do 5 rewrite <- Z2R_mult. + rewrite <- Z2R_plus. + change 1 with (Z2R 1); rewrite <- Z2R_minus. + apply Z2R_le. + apply (Zplus_le_reg_r _ _ 1); ring_simplify. + apply Zlt_le_succ. + apply lt_Z2R. + rewrite Z2R_plus. + do 5 rewrite Z2R_mult; simpl. + rewrite Z2R_Zpower; [|omega]. + apply (Rmult_lt_reg_r (bpow (fexp1 (ln_beta x)))); + [now apply bpow_gt_0|]. + rewrite Rmult_assoc. + rewrite <- Fx. + rewrite (Rmult_plus_distr_r _ _ (Fcore_Raux.bpow _ _)). + bpow_simplify. + rewrite Rmult_assoc. + rewrite bpow_plus. + rewrite <- (Rmult_assoc (Z2R (Zfloor _))). + change (Z2R (Zfloor _) * _) with x'. + do 2 rewrite (Rmult_comm _ (bpow (fexp1 (ln_beta y)))). + rewrite Rmult_assoc. + do 2 rewrite <- (Rmult_assoc (Z2R my)). + rewrite <- Fy. + change (bpow _) with u1. + apply (Rmult_lt_reg_l (/ 2)); [lra|]. + rewrite Rmult_plus_distr_l. + do 4 rewrite <- Rmult_assoc. + rewrite Rinv_l; [|lra]; do 2 rewrite Rmult_1_l. + apply (Rplus_lt_reg_r (- y * x')); ring_simplify. + apply (Rmult_lt_reg_l (/ y)); [now apply Rinv_0_lt_compat|]. + rewrite Rmult_minus_distr_l. + do 3 rewrite <- Rmult_assoc. + rewrite Rinv_l; [|now apply Rgt_not_eq]; do 2 rewrite Rmult_1_l. + now rewrite Rmult_comm. + + apply Rplus_lt_compat_l. + apply Ropp_lt_contravar. + apply Rlt_le_trans with (u2 * bpow (ln_beta y)). + * { apply Rmult_lt_compat_l. + - apply bpow_gt_0. + - apply Rabs_lt_inv. + apply bpow_ln_beta_gt. } + * unfold u2, ulp, canonic_exp; bpow_simplify; apply bpow_le. + apply (Zplus_le_reg_r _ _ (- ln_beta y)); ring_simplify. + rewrite (Zplus_comm (- _)). + destruct (ln_beta_div_disj x y Px Py) as [Hxy|Hxy]; rewrite Hxy; + apply Hexp; try assumption; rewrite <- Hxy; omega. +Qed. + +Lemma double_round_div_aux2 : + forall fexp1 fexp2 : Z -> Z, + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + double_round_div_hyp fexp1 fexp2 -> + forall x y, + 0 < x -> 0 < y -> + generic_format beta fexp1 x -> + generic_format beta fexp1 y -> + (fexp1 (ln_beta (x / y)) <= ln_beta (x / y))%Z -> + ~ (midp fexp1 (x / y) + < x / y + <= midp fexp1 (x / y) + / 2 * ulp beta fexp2 (x / y)). +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Hexp x y Px Py Fx Fy Hf1. +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|]|]. +cut (~ (/ 2 * ulp beta fexp1 (x / y) + < x / y - round beta fexp1 Zfloor (x / y) + <= / 2 * (ulp beta fexp1 (x / y) + ulp beta fexp2 (x / y)))). +{ intro H; intro H'; apply H; split. + - apply (Rplus_lt_reg_l (round beta fexp1 Zfloor (x / y))). + ring_simplify. + apply H'. + - 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 (x' := round beta fexp1 Zfloor (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. +intro Hlr. +apply (Rlt_irrefl (/ 2 * (u1 + u2))). +apply Rlt_le_trans with (x / y - x'); [|now apply Hlr]. +apply (Rplus_lt_reg_r x'); ring_simplify. +apply (Rmult_lt_reg_r y _ _ Py). +unfold Rdiv; rewrite Rmult_assoc. +rewrite Rinv_l; [|now apply Rgt_not_eq]; rewrite Rmult_1_r. +do 2 rewrite Rmult_plus_distr_r. +apply (Rmult_lt_reg_l 2); [lra|]. +do 2 rewrite Rmult_plus_distr_l. +do 5 rewrite <- Rmult_assoc. +rewrite Rinv_r; [|lra]; do 2 rewrite Rmult_1_l. +destruct (Zle_or_lt Z0 (fexp1 (ln_beta x) - fexp1 (ln_beta (x / y)) + - fexp1 (ln_beta y))%Z) as [He|He]. +- (* fexp1 (ln_beta (x / y)) + fexp1 (ln_beta y) <= fexp1 (ln_beta x) *) + apply Rlt_le_trans with (u1 * y + bpow (fexp1 (ln_beta (x / y)) + + fexp1 (ln_beta y)) + + 2 * x' * y). + + apply Rplus_lt_compat_r, Rplus_lt_compat_l. + apply Rlt_le_trans with (u2 * bpow (ln_beta y)). + * { apply Rmult_lt_compat_l. + - apply bpow_gt_0. + - apply Rabs_lt_inv. + apply bpow_ln_beta_gt. } + * unfold u2, ulp, canonic_exp; bpow_simplify; apply bpow_le. + apply (Zplus_le_reg_r _ _ (- ln_beta y)); ring_simplify. + rewrite <- Zplus_assoc; rewrite (Zplus_comm (- _)). + destruct (ln_beta_div_disj x y Px Py) as [Hxy|Hxy]; rewrite Hxy; + [now apply Hexp; [| |rewrite <- Hxy]|]. + replace (_ - _ + 1)%Z with ((ln_beta x + 1) - ln_beta y)%Z by ring. + apply Hexp. + { now assert (fexp1 (ln_beta x + 1) <= ln_beta x)%Z; + [apply valid_exp|omega]. } + { assumption. } + replace (_ + 1 - _)%Z with (ln_beta x - ln_beta y + 1)%Z by ring. + now rewrite <- Hxy. + + apply Rge_le; rewrite Fx at 1; apply Rle_ge. + replace (u1 * y) with (u1 * (Z2R my * bpow (fexp1 (ln_beta y)))); + [|now apply eq_sym; rewrite Fy at 1]. + replace (2 * x' * y) with (2 * x' * (Z2R my * bpow (fexp1 (ln_beta y)))); + [|now apply eq_sym; rewrite Fy at 1]. + apply (Rmult_le_reg_r (bpow (- fexp1 (ln_beta (x / y)) + - fexp1 (ln_beta y)))); + [now apply bpow_gt_0|]. + do 2 rewrite Rmult_plus_distr_r. + bpow_simplify. + rewrite (Rmult_comm u1). + unfold u1, ulp, canonic_exp; bpow_simplify. + rewrite (Rmult_assoc 2). + rewrite (Rmult_comm x'). + rewrite (Rmult_assoc 2). + unfold x', round, F2R, scaled_mantissa, canonic_exp; simpl. + bpow_simplify. + rewrite <- (Z2R_Zpower _ (_ - _)%Z); [|exact He]. + change 2 with (Z2R 2). + do 4 rewrite <- Z2R_mult. + change 1 with (Z2R 1); do 2 rewrite <- Z2R_plus. + apply Z2R_le. + rewrite Zplus_comm, Zplus_assoc. + apply Zlt_le_succ. + apply lt_Z2R. + rewrite Z2R_plus. + do 4 rewrite Z2R_mult; simpl. + rewrite Z2R_Zpower; [|exact He]. + apply (Rmult_lt_reg_r (bpow (fexp1 (ln_beta y)))); + [now apply bpow_gt_0|]. + rewrite Rmult_plus_distr_r. + rewrite (Rmult_comm _ (Z2R _)). + do 2 rewrite Rmult_assoc. + rewrite <- Fy. + bpow_simplify. + unfold Zminus; rewrite bpow_plus. + rewrite (Rmult_assoc _ (Z2R mx)). + rewrite <- (Rmult_assoc (Z2R mx)). + rewrite <- Fx. + apply (Rmult_lt_reg_r (bpow (fexp1 (ln_beta (x / y))))); + [now apply bpow_gt_0|]. + rewrite Rmult_plus_distr_r. + bpow_simplify. + rewrite (Rmult_comm _ y). + do 2 rewrite Rmult_assoc. + change (Z2R _ * _) with x'. + change (bpow _) with u1. + apply (Rmult_lt_reg_l (/ 2)); [lra|]. + rewrite Rmult_plus_distr_l. + do 4 rewrite <- Rmult_assoc. + rewrite Rinv_l; [|lra]; do 2 rewrite Rmult_1_l. + apply (Rplus_lt_reg_r (- y * x')); ring_simplify. + apply (Rmult_lt_reg_l (/ y)); [now apply Rinv_0_lt_compat|]. + rewrite Rmult_plus_distr_l. + do 3 rewrite <- Rmult_assoc. + rewrite Ropp_mult_distr_r_reverse. + rewrite Ropp_mult_distr_l_reverse. + rewrite Rinv_l; [|now apply Rgt_not_eq]; do 2 rewrite Rmult_1_l. + rewrite (Rmult_comm (/ y)). + now rewrite (Rplus_comm (- x')). +- (* fexp1 (ln_beta x) < fexp1 (ln_beta (x / y)) + fexp1 (ln_beta y) *) + apply Rlt_le_trans with (2 * x' * y + u1 * y + bpow (fexp1 (ln_beta x))). + + rewrite Rplus_comm, Rplus_assoc; do 2 apply Rplus_lt_compat_l. + apply Rlt_le_trans with (u2 * bpow (ln_beta y)). + * apply Rmult_lt_compat_l. + now apply bpow_gt_0. + now apply Rabs_lt_inv; apply bpow_ln_beta_gt. + * unfold u2, ulp, canonic_exp; bpow_simplify; apply bpow_le. + apply (Zplus_le_reg_r _ _ (- ln_beta y)); ring_simplify. + rewrite (Zplus_comm (- _)). + destruct (ln_beta_div_disj x y Px Py) as [Hxy|Hxy]; rewrite Hxy; + apply Hexp; try assumption; rewrite <- Hxy; omega. + + apply Rge_le; rewrite Fx at 1; apply Rle_ge. + rewrite Fy at 1 2. + apply (Rmult_le_reg_r (bpow (- fexp1 (ln_beta x)))); + [now apply bpow_gt_0|]. + do 2 rewrite Rmult_plus_distr_r. + bpow_simplify. + replace (2 * x' * _ * _) + with (2 * Z2R my * x' * bpow (fexp1 (ln_beta y) - fexp1 (ln_beta x))) by ring. + rewrite (Rmult_comm u1). + unfold x', u1, round, F2R, ulp, scaled_mantissa, canonic_exp; simpl. + bpow_simplify. + rewrite <- (Z2R_Zpower _ (_ - _)%Z); [|omega]. + change 2 with (Z2R 2). + do 5 rewrite <- Z2R_mult. + change 1 with (Z2R 1); do 2 rewrite <- Z2R_plus. + apply Z2R_le. + apply Zlt_le_succ. + apply lt_Z2R. + rewrite Z2R_plus. + do 5 rewrite Z2R_mult; simpl. + rewrite Z2R_Zpower; [|omega]. + apply (Rmult_lt_reg_r (bpow (fexp1 (ln_beta x)))); + [now apply bpow_gt_0|]. + rewrite (Rmult_assoc _ (Z2R mx)). + rewrite <- Fx. + rewrite Rmult_plus_distr_r. + bpow_simplify. + rewrite bpow_plus. + rewrite Rmult_assoc. + rewrite <- (Rmult_assoc (Z2R _)). + change (Z2R _ * bpow _) with x'. + do 2 rewrite (Rmult_comm _ (bpow (fexp1 (ln_beta y)))). + rewrite Rmult_assoc. + do 2 rewrite <- (Rmult_assoc (Z2R my)). + rewrite <- Fy. + change (bpow _) with u1. + apply (Rmult_lt_reg_l (/ 2)); [lra|]. + rewrite Rmult_plus_distr_l. + do 4 rewrite <- Rmult_assoc. + rewrite Rinv_l; [|lra]; do 2 rewrite Rmult_1_l. + apply (Rplus_lt_reg_r (- y * x')); ring_simplify. + apply (Rmult_lt_reg_l (/ y)); [now apply Rinv_0_lt_compat|]. + rewrite Rmult_plus_distr_l. + do 3 rewrite <- Rmult_assoc. + rewrite Ropp_mult_distr_r_reverse. + rewrite Ropp_mult_distr_l_reverse. + rewrite Rinv_l; [|now apply Rgt_not_eq]; do 2 rewrite Rmult_1_l. + rewrite (Rmult_comm (/ y)). + now rewrite (Rplus_comm (- x')). +Qed. + +Lemma double_round_div_aux : + forall fexp1 fexp2 : Z -> Z, + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + (exists n, (beta = 2 * n :> Z)%Z) -> + double_round_div_hyp fexp1 fexp2 -> + forall x y, + 0 < x -> 0 < y -> + generic_format beta fexp1 x -> + generic_format beta fexp1 y -> + double_round_eq fexp1 fexp2 choice1 choice2 (x / y). +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Ebeta Hexp x y Px Py Fx Fy. +assert (Pxy : 0 < x / y). +{ apply Rmult_lt_0_compat; [exact Px|]. + now apply Rinv_0_lt_compat. } +apply double_round_all_mid_cases. +- exact Vfexp1. +- exact Vfexp2. +- exact Pxy. +- apply Hexp. +- intros Hf1 Hlxy. + casetype False. + now apply (double_round_div_aux0 fexp1 fexp2 _ _ choice1 choice2 Hexp x y). +- intros Hf1 Hlxy. + casetype False. + now apply (double_round_div_aux1 fexp1 fexp2 _ _ choice1 choice2 Hexp x y). +- intro H. + apply double_round_eq_mid_beta_even; try assumption. + apply Hexp. +- intros Hf1 Hlxy. + casetype False. + now apply (double_round_div_aux2 fexp1 fexp2 _ _ choice1 choice2 Hexp x y). +Qed. + +Lemma double_round_div : + forall fexp1 fexp2 : Z -> Z, + Valid_exp fexp1 -> Valid_exp fexp2 -> + forall (choice1 choice2 : Z -> bool), + (exists n, (beta = 2 * n :> Z)%Z) -> + double_round_div_hyp fexp1 fexp2 -> + forall x y, + y <> 0 -> + generic_format beta fexp1 x -> + generic_format beta fexp1 y -> + double_round_eq fexp1 fexp2 choice1 choice2 (x / y). +Proof. +intros fexp1 fexp2 Vfexp1 Vfexp2 choice1 choice2 Ebeta Hexp x y Nzy Fx Fy. +unfold double_round_eq. +destruct (Rtotal_order x 0) as [Nx|[Zx|Px]]. +- (* x < 0 *) + destruct (Rtotal_order y 0) as [Ny|[Zy|Py]]. + + (* y < 0 *) + rewrite <- (Ropp_involutive x). + rewrite <- (Ropp_involutive y). + rewrite Ropp_div. + unfold Rdiv; rewrite <- Ropp_inv_permute; [|lra]. + rewrite Ropp_mult_distr_r_reverse. + rewrite Ropp_involutive. + fold ((- x) / (- y)). + apply Ropp_lt_contravar in Nx. + apply Ropp_lt_contravar in Ny. + rewrite Ropp_0 in Nx, Ny. + apply generic_format_opp in Fx. + apply generic_format_opp in Fy. + now apply double_round_div_aux. + + (* y = 0 *) + now casetype False; apply Nzy. + + (* y > 0 *) + rewrite <- (Ropp_involutive x). + rewrite Ropp_div. + do 3 rewrite round_N_opp. + apply Ropp_eq_compat. + apply Ropp_lt_contravar in Nx. + rewrite Ropp_0 in Nx. + apply generic_format_opp in Fx. + now apply double_round_div_aux. +- (* x = 0 *) + rewrite Zx. + unfold Rdiv; rewrite Rmult_0_l. + now rewrite round_0; [|apply valid_rnd_N]. +- (* x > 0 *) + destruct (Rtotal_order y 0) as [Ny|[Zy|Py]]. + + (* y < 0 *) + rewrite <- (Ropp_involutive y). + unfold Rdiv; rewrite <- Ropp_inv_permute; [|lra]. + rewrite Ropp_mult_distr_r_reverse. + do 3 rewrite round_N_opp. + apply Ropp_eq_compat. + apply Ropp_lt_contravar in Ny. + rewrite Ropp_0 in Ny. + apply generic_format_opp in Fy. + now apply double_round_div_aux. + + (* y = 0 *) + now casetype False; apply Nzy. + + (* y > 0 *) + now apply double_round_div_aux. +Qed. + +Section Double_round_div_FLX. + +Require Import Fcore_FLX. + +Variable prec : Z. +Variable prec' : Z. + +Context { prec_gt_0_ : Prec_gt_0 prec }. +Context { prec_gt_0_' : Prec_gt_0 prec' }. + +Lemma FLX_double_round_div_hyp : + (2 * prec <= prec')%Z -> + double_round_div_hyp (FLX_exp prec) (FLX_exp prec'). +Proof. +intros Hprec. +unfold Prec_gt_0 in prec_gt_0_. +unfold FLX_exp. +unfold double_round_div_hyp. +split; [now intro ex; omega|]. +split; [|split; [|split]]; intros ex ey; omega. +Qed. + +Theorem double_round_div_FLX : + forall choice1 choice2, + (exists n, (beta = 2 * n :> Z)%Z) -> + (2 * prec <= prec')%Z -> + forall x y, + y <> 0 -> + FLX_format beta prec x -> FLX_format beta prec y -> + double_round_eq (FLX_exp prec) (FLX_exp prec') choice1 choice2 (x / y). +Proof. +intros choice1 choice2 Ebeta Hprec x y Nzy Fx Fy. +apply double_round_div. +- now apply FLX_exp_valid. +- now apply FLX_exp_valid. +- exact Ebeta. +- now apply FLX_double_round_div_hyp. +- exact Nzy. +- now apply generic_format_FLX. +- now apply generic_format_FLX. +Qed. + +End Double_round_div_FLX. + +Section Double_round_div_FLT. + +Require Import Fcore_FLX. +Require Import Fcore_FLT. + +Variable emin prec : Z. +Variable emin' prec' : Z. + +Context { prec_gt_0_ : Prec_gt_0 prec }. +Context { prec_gt_0_' : Prec_gt_0 prec' }. + +Lemma FLT_double_round_div_hyp : + (emin' <= emin - prec - 2)%Z -> + (2 * prec <= prec')%Z -> + double_round_div_hyp (FLT_exp emin prec) (FLT_exp emin' prec'). +Proof. +intros Hemin Hprec. +unfold FLT_exp. +unfold Prec_gt_0 in prec_gt_0_. +unfold double_round_div_hyp. +split; [intro ex|split; [|split; [|split]]; intros ex ey]. +- generalize (Zmax_spec (ex - prec') emin'). + generalize (Zmax_spec (ex - prec) emin). + omega. +- generalize (Zmax_spec (ex - prec) emin). + generalize (Zmax_spec (ey - prec) emin). + generalize (Zmax_spec (ex - ey - prec) emin). + generalize (Zmax_spec (ex - ey - prec') emin'). + omega. +- generalize (Zmax_spec (ex - prec) emin). + generalize (Zmax_spec (ey - prec) emin). + generalize (Zmax_spec (ex - ey + 1 - prec) emin). + generalize (Zmax_spec (ex - ey + 1 - prec') emin'). + omega. +- generalize (Zmax_spec (ex - prec) emin). + generalize (Zmax_spec (ey - prec) emin). + generalize (Zmax_spec (ex - ey - prec) emin). + generalize (Zmax_spec (ex - ey - prec') emin'). + omega. +- generalize (Zmax_spec (ex - prec) emin). + generalize (Zmax_spec (ey - prec) emin). + generalize (Zmax_spec (ex - ey - prec) emin). + generalize (Zmax_spec (ex - ey - prec') emin'). + omega. +Qed. + +Theorem double_round_div_FLT : + forall choice1 choice2, + (exists n, (beta = 2 * n :> Z)%Z) -> + (emin' <= emin - prec - 2)%Z -> + (2 * prec <= prec')%Z -> + forall x y, + y <> 0 -> + FLT_format beta emin prec x -> FLT_format beta emin prec y -> + double_round_eq (FLT_exp emin prec) (FLT_exp emin' prec') + choice1 choice2 (x / y). +Proof. +intros choice1 choice2 Ebeta Hemin Hprec x y Nzy Fx Fy. +apply double_round_div. +- now apply FLT_exp_valid. +- now apply FLT_exp_valid. +- exact Ebeta. +- now apply FLT_double_round_div_hyp. +- exact Nzy. +- now apply generic_format_FLT. +- now apply generic_format_FLT. +Qed. + +End Double_round_div_FLT. + +Section Double_round_div_FTZ. + +Require Import Fcore_FLX. +Require Import Fcore_FTZ. + +Variable emin prec : Z. +Variable emin' prec' : Z. + +Context { prec_gt_0_ : Prec_gt_0 prec }. +Context { prec_gt_0_' : Prec_gt_0 prec' }. + +Lemma FTZ_double_round_div_hyp : + (emin' + prec' <= emin - 1)%Z -> + (2 * prec <= prec')%Z -> + double_round_div_hyp (FTZ_exp emin prec) (FTZ_exp emin' prec'). +Proof. +intros Hemin Hprec. +unfold FTZ_exp. +unfold Prec_gt_0 in prec_gt_0_. +unfold Prec_gt_0 in prec_gt_0_. +unfold double_round_div_hyp. +split; [intro ex|split; [|split; [|split]]; intros ex ey]. +- destruct (Z.ltb_spec (ex - prec') emin'); + destruct (Z.ltb_spec (ex - prec) emin); + omega. +- destruct (Z.ltb_spec (ex - prec) emin); + destruct (Z.ltb_spec (ey - prec) emin); + destruct (Z.ltb_spec (ex - ey - prec) emin); + destruct (Z.ltb_spec (ex - ey - prec') emin'); + omega. +- destruct (Z.ltb_spec (ex - prec) emin); + destruct (Z.ltb_spec (ey - prec) emin); + destruct (Z.ltb_spec (ex - ey + 1 - prec) emin); + destruct (Z.ltb_spec (ex - ey + 1 - prec') emin'); + omega. +- destruct (Z.ltb_spec (ex - prec) emin); + destruct (Z.ltb_spec (ey - prec) emin); + destruct (Z.ltb_spec (ex - ey - prec) emin); + destruct (Z.ltb_spec (ex - ey - prec') emin'); + omega. +- destruct (Z.ltb_spec (ex - prec) emin); + destruct (Z.ltb_spec (ey - prec) emin); + destruct (Z.ltb_spec (ex - ey - prec) emin); + destruct (Z.ltb_spec (ex - ey - prec') emin'); + omega. +Qed. + +Theorem double_round_div_FTZ : + forall choice1 choice2, + (exists n, (beta = 2 * n :> Z)%Z) -> + (emin' + prec' <= emin - 1)%Z -> + (2 * prec <= prec')%Z -> + forall x y, + y <> 0 -> + FTZ_format beta emin prec x -> FTZ_format beta emin prec y -> + double_round_eq (FTZ_exp emin prec) (FTZ_exp emin' prec') + choice1 choice2 (x / y). +Proof. +intros choice1 choice2 Ebeta Hemin Hprec x y Nzy Fx Fy. +apply double_round_div. +- now apply FTZ_exp_valid. +- now apply FTZ_exp_valid. +- exact Ebeta. +- now apply FTZ_double_round_div_hyp. +- exact Nzy. +- now apply generic_format_FTZ. +- now apply generic_format_FTZ. +Qed. + +End Double_round_div_FTZ. + +End Double_round_div. + +End Double_round. diff --git a/flocq/Appli/Fappli_rnd_odd.v b/flocq/Appli/Fappli_rnd_odd.v index b4a2c522..b3244589 100644 --- a/flocq/Appli/Fappli_rnd_odd.v +++ b/flocq/Appli/Fappli_rnd_odd.v @@ -802,7 +802,7 @@ apply Hz1. Qed. -Theorem round_odd_prop_pos: +Theorem round_odd_prop_pos: round beta fexp (Znearest choice) (round beta fexpe Zrnd_odd x) = round beta fexp (Znearest choice) x. Proof with auto with typeclass_instances. @@ -945,7 +945,7 @@ Qed. -Theorem round_odd_prop: forall x, +Theorem round_odd_prop: forall x, round beta fexp (Znearest choice) (round beta fexpe Zrnd_odd x) = round beta fexp (Znearest choice) x. Proof with auto with typeclass_instances. |