diff options
Diffstat (limited to 'flocq/Core/Round_NE.v')
| -rw-r--r-- | flocq/Core/Round_NE.v | 547 |
1 files changed, 547 insertions, 0 deletions
diff --git a/flocq/Core/Round_NE.v b/flocq/Core/Round_NE.v new file mode 100644 index 00000000..20b60ef5 --- /dev/null +++ b/flocq/Core/Round_NE.v @@ -0,0 +1,547 @@ +(** +This file is part of the Flocq formalization of floating-point +arithmetic in Coq: http://flocq.gforge.inria.fr/ + +Copyright (C) 2009-2018 Sylvie Boldo +#<br /># +Copyright (C) 2009-2018 Guillaume Melquiond + +This library is free software; you can redistribute it and/or +modify it under the terms of the GNU Lesser General Public +License as published by the Free Software Foundation; either +version 3 of the License, or (at your option) any later version. + +This library is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +COPYING file for more details. +*) + +(** * Rounding to nearest, ties to even: existence, unicity... *) +Require Import Raux Defs Round_pred Generic_fmt Float_prop Ulp. + +Notation ZnearestE := (Znearest (fun x => negb (Z.even x))). + +Section Fcore_rnd_NE. + +Variable beta : radix. + +Notation bpow e := (bpow beta e). + +Variable fexp : Z -> Z. + +Context { valid_exp : Valid_exp fexp }. + +Notation format := (generic_format beta fexp). +Notation canonical := (canonical beta fexp). + +Definition NE_prop (_ : R) f := + exists g : float beta, f = F2R g /\ canonical g /\ Z.even (Fnum g) = true. + +Definition Rnd_NE_pt := + Rnd_NG_pt format NE_prop. + +Definition DN_UP_parity_pos_prop := + forall x xd xu, + (0 < x)%R -> + ~ format x -> + canonical xd -> + canonical xu -> + F2R xd = round beta fexp Zfloor x -> + F2R xu = round beta fexp Zceil x -> + Z.even (Fnum xu) = negb (Z.even (Fnum xd)). + +Definition DN_UP_parity_prop := + forall x xd xu, + ~ format x -> + canonical xd -> + canonical xu -> + F2R xd = round beta fexp Zfloor x -> + F2R xu = round beta fexp Zceil x -> + Z.even (Fnum xu) = negb (Z.even (Fnum xd)). + +Lemma DN_UP_parity_aux : + DN_UP_parity_pos_prop -> + DN_UP_parity_prop. +Proof. +intros Hpos x xd xu Hfx Hd Hu Hxd Hxu. +destruct (total_order_T 0 x) as [[Hx|Hx]|Hx]. +(* . *) +exact (Hpos x xd xu Hx Hfx Hd Hu Hxd Hxu). +elim Hfx. +rewrite <- Hx. +apply generic_format_0. +(* . *) +assert (Hx': (0 < -x)%R). +apply Ropp_lt_cancel. +now rewrite Ropp_involutive, Ropp_0. +destruct xd as (md, ed). +destruct xu as (mu, eu). +simpl. +rewrite <- (Bool.negb_involutive (Z.even mu)). +apply f_equal. +apply sym_eq. +rewrite <- (Z.even_opp mu), <- (Z.even_opp md). +change (Z.even (Fnum (Float beta (-md) ed)) = negb (Z.even (Fnum (Float beta (-mu) eu)))). +apply (Hpos (-x)%R _ _ Hx'). +intros H. +apply Hfx. +rewrite <- Ropp_involutive. +now apply generic_format_opp. +now apply canonical_opp. +now apply canonical_opp. +rewrite round_DN_opp, F2R_Zopp. +now apply f_equal. +rewrite round_UP_opp, F2R_Zopp. +now apply f_equal. +Qed. + +Class Exists_NE := + exists_NE : Z.even beta = false \/ forall e, + ((fexp e < e)%Z -> (fexp (e + 1) < e)%Z) /\ ((e <= fexp e)%Z -> fexp (fexp e + 1) = fexp e). + +Context { exists_NE_ : Exists_NE }. + +Theorem DN_UP_parity_generic_pos : + DN_UP_parity_pos_prop. +Proof with auto with typeclass_instances. +intros x xd xu H0x Hfx Hd Hu Hxd Hxu. +destruct (mag beta x) as (ex, Hexa). +specialize (Hexa (Rgt_not_eq _ _ H0x)). +generalize Hexa. intros Hex. +rewrite (Rabs_pos_eq _ (Rlt_le _ _ H0x)) in Hex. +destruct (Zle_or_lt ex (fexp ex)) as [Hxe|Hxe]. +(* small x *) +assert (Hd3 : Fnum xd = Z0). +apply eq_0_F2R with beta (Fexp xd). +change (F2R xd = R0). +rewrite Hxd. +apply round_DN_small_pos with (1 := Hex) (2 := Hxe). +assert (Hu3 : xu = Float beta (1 * Zpower beta (fexp ex - fexp (fexp ex + 1))) (fexp (fexp ex + 1))). +apply canonical_unique with (1 := Hu). +apply (f_equal fexp). +rewrite <- F2R_change_exp. +now rewrite F2R_bpow, mag_bpow. +now apply valid_exp. +rewrite <- F2R_change_exp. +rewrite F2R_bpow. +apply sym_eq. +rewrite Hxu. +apply sym_eq. +apply round_UP_small_pos with (1 := Hex) (2 := Hxe). +now apply valid_exp. +rewrite Hd3, Hu3. +rewrite Zmult_1_l. +simpl. +destruct exists_NE_ as [H|H]. +apply Zeven_Zpower_odd with (2 := H). +apply Zle_minus_le_0. +now apply valid_exp. +rewrite (proj2 (H ex)). +now rewrite Zminus_diag. +exact Hxe. +(* large x *) +assert (Hd4: (bpow (ex - 1) <= Rabs (F2R xd) < bpow ex)%R). +rewrite Rabs_pos_eq. +rewrite Hxd. +split. +apply (round_DN_pt beta fexp x). +apply generic_format_bpow. +ring_simplify (ex - 1 + 1)%Z. +omega. +apply Hex. +apply Rle_lt_trans with (2 := proj2 Hex). +apply (round_DN_pt beta fexp x). +rewrite Hxd. +apply (round_DN_pt beta fexp x). +apply generic_format_0. +now apply Rlt_le. +assert (Hxe2 : (fexp (ex + 1) <= ex)%Z) by now apply valid_exp. +assert (Hud: (F2R xu = F2R xd + ulp beta fexp x)%R). +rewrite Hxu, Hxd. +now apply round_UP_DN_ulp. +destruct (total_order_T (bpow ex) (F2R xu)) as [[Hu2|Hu2]|Hu2]. +(* - xu > bpow ex *) +elim (Rlt_not_le _ _ Hu2). +rewrite Hxu. +apply round_bounded_large_pos... +(* - xu = bpow ex *) +assert (Hu3: xu = Float beta (1 * Zpower beta (ex - fexp (ex + 1))) (fexp (ex + 1))). +apply canonical_unique with (1 := Hu). +apply (f_equal fexp). +rewrite <- F2R_change_exp. +now rewrite F2R_bpow, mag_bpow. +now apply valid_exp. +rewrite <- Hu2. +apply sym_eq. +rewrite <- F2R_change_exp. +apply F2R_bpow. +exact Hxe2. +assert (Hd3: xd = Float beta (Zpower beta (ex - fexp ex) - 1) (fexp ex)). +assert (H: F2R xd = F2R (Float beta (Zpower beta (ex - fexp ex) - 1) (fexp ex))). +unfold F2R. simpl. +rewrite minus_IZR. +unfold Rminus. +rewrite Rmult_plus_distr_r. +rewrite IZR_Zpower, <- bpow_plus. +ring_simplify (ex - fexp ex + fexp ex)%Z. +rewrite Hu2, Hud. +rewrite ulp_neq_0;[idtac|now apply Rgt_not_eq]. +unfold cexp. +rewrite mag_unique with beta x ex. +unfold F2R. +simpl. ring. +rewrite Rabs_pos_eq. +exact Hex. +now apply Rlt_le. +apply Zle_minus_le_0. +now apply Zlt_le_weak. +apply canonical_unique with (1 := Hd) (3 := H). +apply (f_equal fexp). +rewrite <- H. +apply sym_eq. +now apply mag_unique. +rewrite Hd3, Hu3. +unfold Fnum. +rewrite Z.even_mul. simpl. +unfold Zminus at 2. +rewrite Z.even_add. +rewrite eqb_sym. simpl. +fold (negb (Z.even (beta ^ (ex - fexp ex)))). +rewrite Bool.negb_involutive. +rewrite (Z.even_pow beta (ex - fexp ex)). 2: omega. +destruct exists_NE_. +rewrite H. +apply Zeven_Zpower_odd with (2 := H). +now apply Zle_minus_le_0. +apply Z.even_pow. +specialize (H ex). +omega. +(* - xu < bpow ex *) +revert Hud. +rewrite ulp_neq_0;[idtac|now apply Rgt_not_eq]. +unfold F2R. +rewrite Hd, Hu. +unfold cexp. +rewrite mag_unique with beta (F2R xu) ex. +rewrite mag_unique with (1 := Hd4). +rewrite mag_unique with (1 := Hexa). +intros H. +replace (Fnum xu) with (Fnum xd + 1)%Z. +rewrite Z.even_add. +now apply eqb_sym. +apply sym_eq. +apply eq_IZR. +rewrite plus_IZR. +apply Rmult_eq_reg_r with (bpow (fexp ex)). +rewrite H. +simpl. ring. +apply Rgt_not_eq. +apply bpow_gt_0. +rewrite Rabs_pos_eq. +split. +apply Rle_trans with (1 := proj1 Hex). +rewrite Hxu. +apply (round_UP_pt beta fexp x). +exact Hu2. +apply Rlt_le. +apply Rlt_le_trans with (1 := H0x). +rewrite Hxu. +apply (round_UP_pt beta fexp x). +Qed. + +Theorem DN_UP_parity_generic : + DN_UP_parity_prop. +Proof. +apply DN_UP_parity_aux. +apply DN_UP_parity_generic_pos. +Qed. + +Theorem Rnd_NE_pt_total : + round_pred_total Rnd_NE_pt. +Proof. +apply satisfies_any_imp_NG. +now apply generic_format_satisfies_any. +intros x d u Hf Hd Hu. +generalize (proj1 Hd). +unfold generic_format. +set (ed := cexp beta fexp d). +set (md := Ztrunc (scaled_mantissa beta fexp d)). +intros Hd1. +case_eq (Z.even md) ; [ intros He | intros Ho ]. +right. +exists (Float beta md ed). +unfold Generic_fmt.canonical. +rewrite <- Hd1. +now repeat split. +left. +generalize (proj1 Hu). +unfold generic_format. +set (eu := cexp beta fexp u). +set (mu := Ztrunc (scaled_mantissa beta fexp u)). +intros Hu1. +rewrite Hu1. +eexists ; repeat split. +unfold Generic_fmt.canonical. +now rewrite <- Hu1. +rewrite (DN_UP_parity_generic x (Float beta md ed) (Float beta mu eu)). +simpl. +now rewrite Ho. +exact Hf. +unfold Generic_fmt.canonical. +now rewrite <- Hd1. +unfold Generic_fmt.canonical. +now rewrite <- Hu1. +rewrite <- Hd1. +apply Rnd_DN_pt_unique with (1 := Hd). +now apply round_DN_pt. +rewrite <- Hu1. +apply Rnd_UP_pt_unique with (1 := Hu). +now apply round_UP_pt. +Qed. + +Theorem Rnd_NE_pt_monotone : + round_pred_monotone Rnd_NE_pt. +Proof. +apply Rnd_NG_pt_monotone. +intros x d u Hd Hdn Hu Hun (cd, (Hd1, Hd2)) (cu, (Hu1, Hu2)). +destruct (Req_dec x d) as [Hx|Hx]. +rewrite <- Hx. +apply sym_eq. +apply Rnd_UP_pt_idempotent with (1 := Hu). +rewrite Hx. +apply Hd. +rewrite (DN_UP_parity_aux DN_UP_parity_generic_pos x cd cu) in Hu2 ; try easy. +now rewrite (proj2 Hd2) in Hu2. +intros Hf. +apply Hx. +apply sym_eq. +now apply Rnd_DN_pt_idempotent with (1 := Hd). +rewrite <- Hd1. +apply Rnd_DN_pt_unique with (1 := Hd). +now apply round_DN_pt. +rewrite <- Hu1. +apply Rnd_UP_pt_unique with (1 := Hu). +now apply round_UP_pt. +Qed. + +Theorem Rnd_NE_pt_round : + round_pred Rnd_NE_pt. +Proof. +split. +apply Rnd_NE_pt_total. +apply Rnd_NE_pt_monotone. +Qed. + +Lemma round_NE_pt_pos : + forall x, + (0 < x)%R -> + Rnd_NE_pt x (round beta fexp ZnearestE x). +Proof with auto with typeclass_instances. +intros x Hx. +split. +now apply round_N_pt. +unfold NE_prop. +set (mx := scaled_mantissa beta fexp x). +set (xr := round beta fexp ZnearestE x). +destruct (Req_dec (mx - IZR (Zfloor mx)) (/2)) as [Hm|Hm]. +(* midpoint *) +left. +exists (Float beta (Ztrunc (scaled_mantissa beta fexp xr)) (cexp beta fexp xr)). +split. +apply round_N_pt... +split. +unfold Generic_fmt.canonical. simpl. +apply f_equal. +apply round_N_pt... +simpl. +unfold xr, round, Znearest. +fold mx. +rewrite Hm. +rewrite Rcompare_Eq. 2: apply refl_equal. +case_eq (Z.even (Zfloor mx)) ; intros Hmx. +(* . even floor *) +change (Z.even (Ztrunc (scaled_mantissa beta fexp (round beta fexp Zfloor x))) = true). +destruct (Rle_or_lt (round beta fexp Zfloor x) 0) as [Hr|Hr]. +rewrite (Rle_antisym _ _ Hr). +unfold scaled_mantissa. +rewrite Rmult_0_l. +now rewrite Ztrunc_IZR. +rewrite <- (round_0 beta fexp Zfloor). +apply round_le... +now apply Rlt_le. +rewrite scaled_mantissa_DN... +now rewrite Ztrunc_IZR. +(* . odd floor *) +change (Z.even (Ztrunc (scaled_mantissa beta fexp (round beta fexp Zceil x))) = true). +destruct (mag beta x) as (ex, Hex). +specialize (Hex (Rgt_not_eq _ _ Hx)). +rewrite (Rabs_pos_eq _ (Rlt_le _ _ Hx)) in Hex. +destruct (Z_lt_le_dec (fexp ex) ex) as [He|He]. +(* .. large pos *) +assert (Hu := round_bounded_large_pos _ _ Zceil _ _ He Hex). +assert (Hfc: Zceil mx = (Zfloor mx + 1)%Z). +apply Zceil_floor_neq. +intros H. +rewrite H in Hm. +unfold Rminus in Hm. +rewrite Rplus_opp_r in Hm. +elim (Rlt_irrefl 0). +rewrite Hm at 2. +apply Rinv_0_lt_compat. +now apply IZR_lt. +destruct (proj2 Hu) as [Hu'|Hu']. +(* ... u <> bpow *) +unfold scaled_mantissa. +rewrite cexp_fexp_pos with (1 := conj (proj1 Hu) Hu'). +unfold round, F2R. simpl. +rewrite cexp_fexp_pos with (1 := Hex). +rewrite Rmult_assoc, <- bpow_plus, Zplus_opp_r, Rmult_1_r. +rewrite Ztrunc_IZR. +fold mx. +rewrite Hfc. +now rewrite Z.even_add, Hmx. +(* ... u = bpow *) +rewrite Hu'. +unfold scaled_mantissa, cexp. +rewrite mag_bpow. +rewrite <- bpow_plus, <- IZR_Zpower. +rewrite Ztrunc_IZR. +case_eq (Z.even beta) ; intros Hr. +destruct exists_NE_ as [Hs|Hs]. +now rewrite Hs in Hr. +destruct (Hs ex) as (H,_). +rewrite Z.even_pow. +exact Hr. +omega. +assert (Z.even (Zfloor mx) = true). 2: now rewrite H in Hmx. +replace (Zfloor mx) with (Zceil mx + -1)%Z by omega. +rewrite Z.even_add. +apply eqb_true. +unfold mx. +replace (Zceil (scaled_mantissa beta fexp x)) with (Zpower beta (ex - fexp ex)). +rewrite Zeven_Zpower_odd with (2 := Hr). +easy. +omega. +apply eq_IZR. +rewrite IZR_Zpower. 2: omega. +apply Rmult_eq_reg_r with (bpow (fexp ex)). +unfold Zminus. +rewrite bpow_plus. +rewrite Rmult_assoc, <- bpow_plus, Zplus_opp_l, Rmult_1_r. +pattern (fexp ex) ; rewrite <- cexp_fexp_pos with (1 := Hex). +now apply sym_eq. +apply Rgt_not_eq. +apply bpow_gt_0. +generalize (proj1 (valid_exp ex) He). +omega. +(* .. small pos *) +assert (Z.even (Zfloor mx) = true). 2: now rewrite H in Hmx. +unfold mx, scaled_mantissa. +rewrite cexp_fexp_pos with (1 := Hex). +now rewrite mantissa_DN_small_pos. +(* not midpoint *) +right. +intros g Hg. +destruct (Req_dec x g) as [Hxg|Hxg]. +rewrite <- Hxg. +apply sym_eq. +apply round_generic... +rewrite Hxg. +apply Hg. +set (d := round beta fexp Zfloor x). +set (u := round beta fexp Zceil x). +apply Rnd_N_pt_unique with (d := d) (u := u) (4 := Hg). +now apply round_DN_pt. +now apply round_UP_pt. +2: now apply round_N_pt. +rewrite <- (scaled_mantissa_mult_bpow beta fexp x). +unfold d, u, round, F2R. simpl. fold mx. +rewrite <- 2!Rmult_minus_distr_r. +intros H. +apply Rmult_eq_reg_r in H. +apply Hm. +apply Rcompare_Eq_inv. +rewrite Rcompare_floor_ceil_middle. +now apply Rcompare_Eq. +contradict Hxg. +apply sym_eq. +apply Rnd_N_pt_idempotent with (1 := Hg). +rewrite <- (scaled_mantissa_mult_bpow beta fexp x). +fold mx. +rewrite <- Hxg. +change (IZR (Zfloor mx) * bpow (cexp beta fexp x))%R with d. +now eapply round_DN_pt. +apply Rgt_not_eq. +apply bpow_gt_0. +Qed. + +Theorem round_NE_opp : + forall x, + round beta fexp ZnearestE (-x) = (- round beta fexp ZnearestE x)%R. +Proof. +intros x. +unfold round. simpl. +rewrite scaled_mantissa_opp, cexp_opp. +rewrite Znearest_opp. +rewrite <- F2R_Zopp. +apply (f_equal (fun v => F2R (Float beta (-v) _))). +set (m := scaled_mantissa beta fexp x). +unfold Znearest. +case Rcompare ; trivial. +apply (f_equal (fun (b : bool) => if b then Zceil m else Zfloor m)). +rewrite Bool.negb_involutive. +rewrite Z.even_opp. +rewrite Z.even_add. +now rewrite eqb_sym. +Qed. + +Lemma round_NE_abs: + forall x : R, + round beta fexp ZnearestE (Rabs x) = Rabs (round beta fexp ZnearestE x). +Proof with auto with typeclass_instances. +intros x. +apply sym_eq. +unfold Rabs at 2. +destruct (Rcase_abs x) as [Hx|Hx]. +rewrite round_NE_opp. +apply Rabs_left1. +rewrite <- (round_0 beta fexp ZnearestE). +apply round_le... +now apply Rlt_le. +apply Rabs_pos_eq. +rewrite <- (round_0 beta fexp ZnearestE). +apply round_le... +now apply Rge_le. +Qed. + +Theorem round_NE_pt : + forall x, + Rnd_NE_pt x (round beta fexp ZnearestE x). +Proof with auto with typeclass_instances. +intros x. +destruct (total_order_T x 0) as [[Hx|Hx]|Hx]. +apply Rnd_NG_pt_opp_inv. +apply generic_format_opp. +unfold NE_prop. +intros _ f ((mg,eg),(H1,(H2,H3))). +exists (Float beta (- mg) eg). +repeat split. +rewrite H1. +now rewrite F2R_Zopp. +now apply canonical_opp. +simpl. +now rewrite Z.even_opp. +rewrite <- round_NE_opp. +apply round_NE_pt_pos. +now apply Ropp_0_gt_lt_contravar. +rewrite Hx, round_0... +apply Rnd_NG_pt_refl. +apply generic_format_0. +now apply round_NE_pt_pos. +Qed. + +End Fcore_rnd_NE. + +(** Notations for backward-compatibility with Flocq 1.4. *) +Notation rndNE := ZnearestE (only parsing). |