diff options
Diffstat (limited to 'flocq')
29 files changed, 1908 insertions, 223 deletions
diff --git a/flocq/Appli/Fappli_IEEE.v b/flocq/Appli/Fappli_IEEE.v index 63b150f..5e9897f 100644 --- a/flocq/Appli/Fappli_IEEE.v +++ b/flocq/Appli/Fappli_IEEE.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public @@ -33,7 +33,7 @@ Section AnyRadix. Inductive full_float := | F754_zero : bool -> full_float | F754_infinity : bool -> full_float - | F754_nan : full_float + | F754_nan : bool -> positive -> full_float | F754_finite : bool -> positive -> Z -> full_float. Definition FF2R beta x := @@ -67,16 +67,20 @@ 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 | _ => true end. (** Basic type used for representing binary FP numbers. Note that there is exactly one such object per FP datum. NaNs do not have any payload. They cannot be distinguished. *) + +Definition nan_pl := {pl | (Z_of_nat' (S (digits2_Pnat pl)) <? prec)%Z = true}. + Inductive binary_float := | B754_zero : bool -> binary_float | B754_infinity : bool -> binary_float - | B754_nan : binary_float + | B754_nan : bool -> nan_pl -> binary_float | B754_finite : bool -> forall (m : positive) (e : Z), bounded m e = true -> binary_float. @@ -85,7 +89,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 => fun _ => B754_nan + | F754_nan b pl => fun H => B754_nan b (exist _ pl H) end. Definition B2FF x := @@ -93,7 +97,7 @@ Definition B2FF x := | B754_finite s m e _ => F754_finite s m e | B754_infinity s => F754_infinity s | B754_zero s => F754_zero s - | B754_nan => F754_nan + | B754_nan b (exist pl _) => F754_nan b pl end. Definition radix2 := Build_radix 2 (refl_equal true). @@ -108,30 +112,30 @@ Theorem FF2R_B2FF : forall x, FF2R radix2 (B2FF x) = B2R x. Proof. -now intros [sx|sx| |sx mx ex Hx]. +now intros [sx|sx|sx [plx Hplx]|sx mx ex Hx]. Qed. Theorem B2FF_FF2B : forall x Hx, B2FF (FF2B x Hx) = x. Proof. -now intros [sx|sx| |sx mx ex] Hx. +now intros [sx|sx|sx plx|sx mx ex] Hx. Qed. Theorem valid_binary_B2FF : forall x, valid_binary (B2FF x) = true. Proof. -now intros [sx|sx| |sx mx ex Hx]. +now intros [sx|sx|sx [plx Hplx]|sx mx ex Hx]. Qed. Theorem FF2B_B2FF : forall x H, FF2B (B2FF x) H = x. Proof. -intros [sx|sx| |sx mx ex Hx] H ; try easy. -apply f_equal. -apply eqbool_irrelevance. +intros [sx|sx|sx [plx Hplx]|sx mx ex Hx] H ; try easy. +simpl. apply f_equal, f_equal, eqbool_irrelevance. +apply f_equal, eqbool_irrelevance. Qed. Theorem FF2B_B2FF_valid : @@ -146,7 +150,7 @@ Theorem B2R_FF2B : forall x Hx, B2R (FF2B x Hx) = FF2R radix2 x. Proof. -now intros [sx|sx| |sx mx ex] Hx. +now intros [sx|sx|sx plx|sx mx ex] Hx. Qed. Theorem match_FF2B : @@ -154,17 +158,17 @@ Theorem match_FF2B : match FF2B x Hx return T with | B754_zero sx => fz sx | B754_infinity sx => fi sx - | B754_nan => fn + | B754_nan b (exist p _) => fn b p | B754_finite sx mx ex _ => ff sx mx ex end = match x with | F754_zero sx => fz sx | F754_infinity sx => fi sx - | F754_nan => fn + | F754_nan b p => fn b p | F754_finite sx mx ex => ff sx mx ex end. Proof. -now intros T fz fi fn ff [sx|sx| |sx mx ex] Hx. +now intros T fz fi fn ff [sx|sx|sx plx|sx mx ex] Hx. Qed. Theorem canonic_canonic_mantissa : @@ -189,19 +193,28 @@ Theorem generic_format_B2R : forall x, generic_format radix2 fexp (B2R x). Proof. -intros [sx|sx| |sx mx ex Hx] ; try apply generic_format_0. +intros [sx|sx|sx plx|sx mx ex Hx] ; try apply generic_format_0. simpl. apply generic_format_canonic. apply canonic_canonic_mantissa. now destruct (andb_prop _ _ Hx) as (H, _). Qed. +Theorem FLT_format_B2R : + forall x, + FLT_format radix2 emin prec (B2R x). +Proof with auto with typeclass_instances. +intros x. +apply FLT_format_generic... +apply generic_format_B2R. +Qed. + Theorem B2FF_inj : forall x y : binary_float, B2FF x = B2FF y -> x = y. Proof. -intros [sx|sx| |sx mx ex Hx] [sy|sy| |sy my ey Hy] ; try easy. +intros [sx|sx|sx [plx Hplx]|sx mx ex Hx] [sy|sy|sy [ply Hply]|sy my ey Hy] ; try easy. (* *) intros H. now inversion H. @@ -212,11 +225,18 @@ now inversion H. intros H. inversion H. clear H. +revert Hplx. +rewrite H2. +intros Hx. +apply f_equal, f_equal, eqbool_irrelevance. +(* *) +intros H. +inversion H. +clear H. revert Hx. rewrite H2, H3. intros Hx. -apply f_equal. -apply eqbool_irrelevance. +apply f_equal, eqbool_irrelevance. Qed. Definition is_finite_strict f := @@ -265,6 +285,29 @@ apply f_equal. apply eqbool_irrelevance. Qed. +Definition Bsign x := + match x with + | B754_nan s _ => s + | B754_zero s => s + | B754_infinity s => s + | B754_finite s _ _ _ => s + end. + +Definition sign_FF x := + match x with + | F754_nan s _ => s + | F754_zero s => s + | F754_infinity s => s + | F754_finite s _ _ => s + end. + +Theorem Bsign_FF2B : + forall x H, + Bsign (FF2B x H) = sign_FF x. +Proof. +now intros [sx|sx|sx plx|sx mx ex] H. +Qed. + Definition is_finite f := match f with | B754_finite _ _ _ _ => true @@ -290,42 +333,93 @@ Theorem is_finite_FF_B2FF : forall x, is_finite_FF (B2FF x) = is_finite x. Proof. +now intros [| |? []|]. +Qed. + +Theorem B2R_Bsign_inj: + forall x y : binary_float, + is_finite x = true -> + is_finite y = true -> + B2R x = B2R y -> + Bsign x = Bsign y -> + x = y. +Proof. +intros. destruct x, y; try (apply B2R_inj; now eauto). +- simpl in H2. congruence. +- symmetry in H1. apply Rmult_integral in H1. + destruct H1. apply eq_Z2R with (n:=0%Z) in H1. destruct b0; discriminate H1. + simpl in H1. pose proof (bpow_gt_0 radix2 e). + rewrite H1 in H3. apply Rlt_irrefl in H3. destruct H3. +- apply Rmult_integral in H1. + destruct H1. apply eq_Z2R with (n:=0%Z) in H1. destruct b; discriminate H1. + simpl in H1. pose proof (bpow_gt_0 radix2 e). + rewrite H1 in H3. apply Rlt_irrefl in H3. destruct H3. +Qed. + +Definition is_nan f := + match f with + | B754_nan _ _ => true + | _ => false + end. + +Definition is_nan_FF f := + match f with + | F754_nan _ _ => true + | _ => false + end. + +Theorem is_nan_FF2B : + forall x Hx, + is_nan (FF2B x Hx) = is_nan_FF x. +Proof. now intros [| | |]. Qed. -Definition Bopp x := +Theorem is_nan_FF_B2FF : + forall x, + is_nan_FF (B2FF x) = is_nan x. +Proof. +now intros [| |? []|]. +Qed. + +Definition Bopp opp_nan x := match x with - | B754_nan => x + | B754_nan sx plx => + let '(sres, plres) := opp_nan sx plx in B754_nan sres plres | B754_infinity sx => B754_infinity (negb sx) | B754_finite sx mx ex Hx => B754_finite (negb sx) mx ex Hx | B754_zero sx => B754_zero (negb sx) end. Theorem Bopp_involutive : - forall x, Bopp (Bopp x) = x. + forall opp_nan x, + is_nan x = false -> + Bopp opp_nan (Bopp opp_nan x) = x. Proof. -now intros [sx|sx| |sx mx ex Hx] ; simpl ; try rewrite Bool.negb_involutive. +now intros opp_nan [sx|sx|sx plx|sx mx ex Hx] ; simpl ; try rewrite Bool.negb_involutive. Qed. Theorem B2R_Bopp : - forall x, - B2R (Bopp x) = (- B2R x)%R. + forall opp_nan x, + B2R (Bopp opp_nan x) = (- B2R x)%R. Proof. -intros [sx|sx| |sx mx ex Hx] ; apply sym_eq ; try apply Ropp_0. +intros opp_nan [sx|sx|sx plx|sx mx ex Hx]; apply sym_eq ; try apply Ropp_0. +simpl. destruct opp_nan. apply Ropp_0. simpl. rewrite <- F2R_opp. now case sx. Qed. - -Theorem is_finite_Bopp: forall x, - is_finite (Bopp x) = is_finite x. +Theorem is_finite_Bopp : + forall opp_nan x, + is_finite (Bopp opp_nan x) = is_finite x. Proof. -now intros [| | |]. +intros opp_nan [| | |] ; try easy. +intros s pl. +simpl. +now case opp_nan. Qed. - - Theorem bounded_lt_emax : forall mx ex, bounded mx ex = true -> @@ -361,7 +455,7 @@ Theorem abs_B2R_lt_emax : forall x, (Rabs (B2R x) < bpow radix2 emax)%R. Proof. -intros [sx|sx| |sx mx ex Hx] ; simpl ; try ( rewrite Rabs_R0 ; apply bpow_gt_0 ). +intros [sx|sx|sx plx|sx mx ex Hx] ; simpl ; try ( rewrite Rabs_R0 ; apply bpow_gt_0 ). rewrite <- F2R_Zabs, abs_cond_Zopp. now apply bounded_lt_emax. Qed. @@ -638,7 +732,7 @@ Definition binary_round_aux mode sx mx ex lx := match shr_m mrs'' with | Z0 => F754_zero sx | Zpos m => if Zle_bool e'' (emax - prec) then F754_finite sx m e'' else binary_overflow mode sx - | _ => F754_nan (* dummy *) + | _ => F754_nan false xH (* dummy *) end. Theorem binary_round_aux_correct : @@ -649,7 +743,7 @@ Theorem binary_round_aux_correct : valid_binary z = true /\ if Rlt_bool (Rabs (round radix2 fexp (round_mode mode) x)) (bpow radix2 emax) then FF2R radix2 z = round radix2 fexp (round_mode mode) x /\ - is_finite_FF z = true + is_finite_FF z = true /\ sign_FF z = Rlt_bool x 0 else z = binary_overflow mode (Rlt_bool x 0). Proof with auto with typeclass_instances. @@ -818,29 +912,6 @@ exact inbetween_int_UP_sign. exact inbetween_int_NA_sign. Qed. -Definition Bsign x := - match x with - | B754_nan => false - | B754_zero s => s - | B754_infinity s => s - | B754_finite s _ _ _ => s - end. - -Definition sign_FF x := - match x with - | F754_nan => false - | F754_zero s => s - | F754_infinity s => s - | F754_finite s _ _ => s - end. - -Theorem Bsign_FF2B : - forall x H, - Bsign (FF2B x H) = sign_FF x. -Proof. -now intros [sx|sx| |sx mx ex] H. -Qed. - (** Multiplication *) Lemma Bmult_correct_aux : @@ -851,7 +922,7 @@ Lemma Bmult_correct_aux : valid_binary z = true /\ if Rlt_bool (Rabs (round radix2 fexp (round_mode m) (x * y))) (bpow radix2 emax) then FF2R radix2 z = round radix2 fexp (round_mode m) (x * y) /\ - is_finite_FF z = true + is_finite_FF z = true /\ sign_FF z = xorb sx sy else z = binary_overflow m (xorb sx sy). Proof. @@ -896,15 +967,15 @@ apply Rlt_bool_false. now apply F2R_ge_0_compat. Qed. -Definition Bmult m x y := +Definition Bmult mult_nan m x y := + let f pl := B754_nan (fst pl) (snd pl) in match x, y with - | B754_nan, _ => x - | _, B754_nan => y + | B754_nan _ _, _ | _, B754_nan _ _ => f (mult_nan x y) | B754_infinity sx, B754_infinity sy => B754_infinity (xorb sx sy) | B754_infinity sx, B754_finite sy _ _ _ => B754_infinity (xorb sx sy) | B754_finite sx _ _ _, B754_infinity sy => B754_infinity (xorb sx sy) - | B754_infinity _, B754_zero _ => B754_nan - | B754_zero _, B754_infinity _ => B754_nan + | B754_infinity _, B754_zero _ => f (mult_nan x y) + | B754_zero _, B754_infinity _ => f (mult_nan x y) | B754_finite sx _ _ _, B754_zero sy => B754_zero (xorb sx sy) | B754_zero sx, B754_finite sy _ _ _ => B754_zero (xorb sx sy) | B754_zero sx, B754_zero sy => B754_zero (xorb sx sy) @@ -913,36 +984,40 @@ Definition Bmult m x y := end. Theorem Bmult_correct : - forall m x y, + forall mult_nan m x y, if Rlt_bool (Rabs (round radix2 fexp (round_mode m) (B2R x * B2R y))) (bpow radix2 emax) then - B2R (Bmult m x y) = round radix2 fexp (round_mode m) (B2R x * B2R y) /\ - is_finite (Bmult m x y) = andb (is_finite x) (is_finite y) + B2R (Bmult mult_nan m x y) = round radix2 fexp (round_mode m) (B2R x * B2R y) /\ + is_finite (Bmult mult_nan m x y) = andb (is_finite x) (is_finite y) /\ + (is_nan (Bmult mult_nan m x y) = false -> + Bsign (Bmult mult_nan m x y) = xorb (Bsign x) (Bsign y)) else - B2FF (Bmult m x y) = binary_overflow m (xorb (Bsign x) (Bsign y)). + B2FF (Bmult mult_nan m x y) = binary_overflow m (xorb (Bsign x) (Bsign y)). Proof. -intros m [sx|sx| |sx mx ex Hx] [sy|sy| |sy my ey Hy] ; - try ( rewrite ?Rmult_0_r, ?Rmult_0_l, round_0, Rabs_R0, Rlt_bool_true ; [ split ; apply refl_equal | apply bpow_gt_0 | auto with typeclass_instances ] ). +intros mult_nan m [sx|sx|sx plx|sx mx ex Hx] [sy|sy|sy ply|sy my ey Hy] ; + try ( rewrite ?Rmult_0_r, ?Rmult_0_l, round_0, Rabs_R0, Rlt_bool_true ; [ now repeat constructor | apply bpow_gt_0 | now auto with typeclass_instances ] ). simpl. case Bmult_correct_aux. intros H1. case Rlt_bool. -intros (H2, H3). +intros (H2, (H3, H4)). split. now rewrite B2R_FF2B. +split. now rewrite is_finite_FF2B. +rewrite Bsign_FF2B. auto. intros H2. now rewrite B2FF_FF2B. Qed. -Definition Bmult_FF m x y := +Definition Bmult_FF mult_nan m x y := + let f pl := F754_nan (fst pl) (snd pl) in match x, y with - | F754_nan, _ => x - | _, F754_nan => y + | F754_nan _ _, _ | _, F754_nan _ _ => f (mult_nan x y) | F754_infinity sx, F754_infinity sy => F754_infinity (xorb sx sy) | F754_infinity sx, F754_finite sy _ _ => F754_infinity (xorb sx sy) | F754_finite sx _ _, F754_infinity sy => F754_infinity (xorb sx sy) - | F754_infinity _, F754_zero _ => F754_nan - | F754_zero _, F754_infinity _ => F754_nan + | F754_infinity _, F754_zero _ => f (mult_nan x y) + | F754_zero _, F754_infinity _ => f (mult_nan x y) | F754_finite sx _ _, F754_zero sy => F754_zero (xorb sx sy) | F754_zero sx, F754_finite sy _ _ => F754_zero (xorb sx sy) | F754_zero sx, F754_zero sy => F754_zero (xorb sx sy) @@ -951,14 +1026,20 @@ Definition Bmult_FF m x y := end. Theorem B2FF_Bmult : + forall mult_nan mult_nan_ff, forall m x y, - B2FF (Bmult m x y) = Bmult_FF m (B2FF x) (B2FF y). + mult_nan_ff (B2FF x) (B2FF y) = (let '(sr, exist plr _) := mult_nan x y in (sr, plr)) -> + B2FF (Bmult mult_nan m x y) = Bmult_FF mult_nan_ff m (B2FF x) (B2FF y). Proof. -intros m [sx|sx| |sx mx ex Hx] [sy|sy| |sy my ey Hy] ; try easy. +intros mult_nan mult_nan_ff m x y Hmult_nan. +unfold Bmult_FF. rewrite Hmult_nan. +destruct x as [sx|sx|sx [plx Hplx]|sx mx ex Hx], y as [sy|sy|sy [ply Hply]|sy my ey Hy] ; + simpl; try match goal with |- context [mult_nan ?x ?y] => + destruct (mult_nan x y) as [? []] end; + try easy. apply B2FF_FF2B. Qed. - Definition shl_align mx ex ex' := match (ex' - ex)%Z with | Zneg d => (shift_pos d mx, ex') @@ -1049,7 +1130,8 @@ Theorem binary_round_correct : let x := F2R (Float radix2 (cond_Zopp sx (Zpos mx)) ex) in if Rlt_bool (Rabs (round radix2 fexp (round_mode m) x)) (bpow radix2 emax) then FF2R radix2 z = round radix2 fexp (round_mode m) x /\ - is_finite_FF z = true + is_finite_FF z = true /\ + sign_FF z = sx else z = binary_overflow m sx. Proof. @@ -1084,22 +1166,35 @@ Theorem binary_normalize_correct : forall m mx ex szero, if Rlt_bool (Rabs (round radix2 fexp (round_mode m) (F2R (Float radix2 mx ex)))) (bpow radix2 emax) then B2R (binary_normalize m mx ex szero) = round radix2 fexp (round_mode m) (F2R (Float radix2 mx ex)) /\ - is_finite (binary_normalize m mx ex szero) = true + is_finite (binary_normalize m mx ex szero) = true /\ + Bsign (binary_normalize m mx ex szero) = + match Rcompare (F2R (Float radix2 mx ex)) 0 with + | Eq => szero + | Lt => true + | Gt => false + end else B2FF (binary_normalize m mx ex szero) = binary_overflow m (Rlt_bool (F2R (Float radix2 mx ex)) 0). Proof with auto with typeclass_instances. intros m mx ez szero. destruct mx as [|mz|mz] ; simpl. rewrite F2R_0, round_0, Rabs_R0, Rlt_bool_true... +split... split... +rewrite Rcompare_Eq... apply bpow_gt_0. (* . mz > 0 *) generalize (binary_round_correct m false mz ez). simpl. case Rlt_bool_spec. -intros _ (Vz, (Rz, Rz')). +intros _ (Vz, (Rz, (Rz', Rz''))). split. now rewrite B2R_FF2B. +split. now rewrite is_finite_FF2B. +rewrite Bsign_FF2B, Rz''. +rewrite Rcompare_Gt... +apply F2R_gt_0_compat. +simpl. zify; omega. intros Hz' (Vz, Rz). rewrite B2FF_FF2B, Rz. apply f_equal. @@ -1110,10 +1205,15 @@ now apply F2R_ge_0_compat. generalize (binary_round_correct m true mz ez). simpl. case Rlt_bool_spec. -intros _ (Vz, (Rz, Rz')). +intros _ (Vz, (Rz, (Rz', Rz''))). split. now rewrite B2R_FF2B. +split. now rewrite is_finite_FF2B. +rewrite Bsign_FF2B, Rz''. +rewrite Rcompare_Lt... +apply F2R_lt_0_compat. +simpl. zify; omega. intros Hz' (Vz, Rz). rewrite B2FF_FF2B, Rz. apply f_equal. @@ -1123,12 +1223,12 @@ now apply F2R_lt_0_compat. Qed. (** Addition *) -Definition Bplus m x y := +Definition Bplus plus_nan m x y := + let f pl := B754_nan (fst pl) (snd pl) in match x, y with - | B754_nan, _ => x - | _, B754_nan => y + | B754_nan _ _, _ | _, B754_nan _ _ => f (plus_nan x y) | B754_infinity sx, B754_infinity sy => - if Bool.eqb sx sy then x else B754_nan + if Bool.eqb sx sy then x else f (plus_nan x y) | B754_infinity _, _ => x | _, B754_infinity _ => y | B754_zero sx, B754_zero sy => @@ -1143,28 +1243,47 @@ Definition Bplus m x y := end. Theorem Bplus_correct : - forall m x y, + forall plus_nan m x y, is_finite x = true -> is_finite y = true -> if Rlt_bool (Rabs (round radix2 fexp (round_mode m) (B2R x + B2R y))) (bpow radix2 emax) then - B2R (Bplus m x y) = round radix2 fexp (round_mode m) (B2R x + B2R y) /\ - is_finite (Bplus m x y) = true + B2R (Bplus plus_nan m x y) = round radix2 fexp (round_mode m) (B2R x + B2R y) /\ + is_finite (Bplus plus_nan m x y) = true /\ + Bsign (Bplus plus_nan m x y) = + match Rcompare (B2R x + B2R y) 0 with + | Eq => match m with mode_DN => orb (Bsign x) (Bsign y) + | _ => andb (Bsign x) (Bsign y) end + | Lt => true + | Gt => false + end else - (B2FF (Bplus m x y) = binary_overflow m (Bsign x) /\ Bsign x = Bsign y). + (B2FF (Bplus plus_nan m x y) = binary_overflow m (Bsign x) /\ Bsign x = Bsign y). Proof with auto with typeclass_instances. -intros m [sx|sx| |sx mx ex Hx] [sy|sy| |sy my ey Hy] Fx Fy ; try easy. +intros plus_nan m [sx|sx| |sx mx ex Hx] [sy|sy| |sy my ey Hy] Fx Fy ; try easy. (* *) rewrite Rplus_0_r, round_0, Rabs_R0, Rlt_bool_true... simpl. -case (Bool.eqb sx sy) ; try easy. -now case m. +rewrite Rcompare_Eq by auto. +destruct sx, sy; try easy; now case m. apply bpow_gt_0. (* *) rewrite Rplus_0_l, round_generic, Rlt_bool_true... +split... split... +simpl. unfold F2R. +erewrite <- Rmult_0_l, Rcompare_mult_r. +rewrite Rcompare_Z2R with (y:=0%Z). +destruct sy... +apply bpow_gt_0. apply abs_B2R_lt_emax. apply generic_format_B2R. (* *) rewrite Rplus_0_r, round_generic, Rlt_bool_true... +split... split... +simpl. unfold F2R. +erewrite <- Rmult_0_l, Rcompare_mult_r. +rewrite Rcompare_Z2R with (y:=0%Z). +destruct sx... +apply bpow_gt_0. apply abs_B2R_lt_emax. apply generic_format_B2R. (* *) @@ -1264,7 +1383,19 @@ now apply F2R_le_0_compat. (* . *) generalize (binary_normalize_correct m mz ez szero). case Rlt_bool_spec. -easy. +split; try easy. split; try easy. +destruct (Rcompare_spec (F2R (beta:=radix2) {| Fnum := mz; Fexp := ez |}) 0); try easy. +rewrite H1 in Hp. +apply Rplus_opp_r_uniq in Hp. +rewrite <- F2R_Zopp in Hp. +eapply canonic_unicity in Hp. +inversion Hp. destruct sy, sx, m; try discriminate H3; easy. +apply canonic_canonic_mantissa. +apply Bool.andb_true_iff in Hy. easy. +replace (-cond_Zopp sx (Z.pos mx))%Z with (cond_Zopp (negb sx) (Z.pos mx)) + by (destruct sx; auto). +apply canonic_canonic_mantissa. +apply Bool.andb_true_iff in Hx. easy. intros Hz' Vz. specialize (Sz Hz'). split. @@ -1273,26 +1404,32 @@ now apply f_equal. apply Sz. Qed. -Definition Bminus m x y := Bplus m x (Bopp y). +Definition Bminus minus_nan m x y := Bplus minus_nan m x (Bopp pair y). Theorem Bminus_correct : - forall m x y, + forall minus_nan m x y, is_finite x = true -> is_finite y = true -> if Rlt_bool (Rabs (round radix2 fexp (round_mode m) (B2R x - B2R y))) (bpow radix2 emax) then - B2R (Bminus m x y) = round radix2 fexp (round_mode m) (B2R x - B2R y) /\ - is_finite (Bminus m x y) = true + B2R (Bminus minus_nan m x y) = round radix2 fexp (round_mode m) (B2R x - B2R y) /\ + is_finite (Bminus minus_nan m x y) = true /\ + Bsign (Bminus minus_nan m x y) = + match Rcompare (B2R x - B2R y) 0 with + | Eq => match m with mode_DN => orb (Bsign x) (negb (Bsign y)) + | _ => andb (Bsign x) (negb (Bsign y)) end + | Lt => true + | Gt => false + end else - (B2FF (Bminus m x y) = binary_overflow m (Bsign x) /\ Bsign x = negb (Bsign y)). + (B2FF (Bminus minus_nan m x y) = binary_overflow m (Bsign x) /\ Bsign x = negb (Bsign y)). Proof with auto with typeclass_instances. -intros m x y Fx Fy. -replace (negb (Bsign y)) with (Bsign (Bopp y)). +intros m minus_nan x y Fx Fy. +replace (negb (Bsign y)) with (Bsign (Bopp pair y)). unfold Rminus. -rewrite <- B2R_Bopp. +erewrite <- B2R_Bopp. apply Bplus_correct. exact Fx. -now rewrite is_finite_Bopp. -now destruct y as [ | | | ]. +rewrite is_finite_Bopp. auto. now destruct y as [ | | | ]. Qed. (** Division *) @@ -1316,12 +1453,12 @@ Lemma Bdiv_correct_aux : let '(mz, ez, lz) := Fdiv_core_binary (Zpos mx) ex (Zpos my) ey in match mz with | Zpos mz => binary_round_aux m (xorb sx sy) mz ez lz - | _ => F754_nan (* dummy *) + | _ => F754_nan false xH (* dummy *) end in valid_binary z = true /\ if Rlt_bool (Rabs (round radix2 fexp (round_mode m) (x / y))) (bpow radix2 emax) then FF2R radix2 z = round radix2 fexp (round_mode m) (x / y) /\ - is_finite_FF z = true + is_finite_FF z = true /\ sign_FF z = xorb sx sy else z = binary_overflow m (xorb sx sy). Proof. @@ -1406,45 +1543,49 @@ apply Rinv_0_lt_compat. now apply F2R_gt_0_compat. Qed. -Definition Bdiv m x y := +Definition Bdiv div_nan m x y := + let f pl := B754_nan (fst pl) (snd pl) in match x, y with - | B754_nan, _ => x - | _, B754_nan => y - | B754_infinity sx, B754_infinity sy => B754_nan + | B754_nan _ _, _ | _, B754_nan _ _ => f (div_nan x y) + | B754_infinity sx, B754_infinity sy => f (div_nan x y) | B754_infinity sx, B754_finite sy _ _ _ => B754_infinity (xorb sx sy) | B754_finite sx _ _ _, B754_infinity sy => B754_zero (xorb sx sy) | B754_infinity sx, B754_zero sy => B754_infinity (xorb sx sy) | B754_zero sx, B754_infinity sy => B754_zero (xorb sx sy) | B754_finite sx _ _ _, B754_zero sy => B754_infinity (xorb sx sy) | B754_zero sx, B754_finite sy _ _ _ => B754_zero (xorb sx sy) - | B754_zero sx, B754_zero sy => B754_nan + | B754_zero sx, B754_zero sy => f (div_nan x y) | B754_finite sx mx ex _, B754_finite sy my ey _ => FF2B _ (proj1 (Bdiv_correct_aux m sx mx ex sy my ey)) end. Theorem Bdiv_correct : - forall m x y, + forall div_nan m x y, B2R y <> R0 -> if Rlt_bool (Rabs (round radix2 fexp (round_mode m) (B2R x / B2R y))) (bpow radix2 emax) then - B2R (Bdiv m x y) = round radix2 fexp (round_mode m) (B2R x / B2R y) /\ - is_finite (Bdiv m x y) = is_finite x + B2R (Bdiv div_nan m x y) = round radix2 fexp (round_mode m) (B2R x / B2R y) /\ + is_finite (Bdiv div_nan m x y) = is_finite x /\ + (is_nan (Bdiv div_nan m x y) = false -> + Bsign (Bdiv div_nan m x y) = xorb (Bsign x) (Bsign y)) else - B2FF (Bdiv m x y) = binary_overflow m (xorb (Bsign x) (Bsign y)). + B2FF (Bdiv div_nan m x y) = binary_overflow m (xorb (Bsign x) (Bsign y)). Proof. -intros m x [sy|sy| |sy my ey Hy] Zy ; try now elim Zy. +intros div_nan m x [sy|sy|sy ply|sy my ey Hy] Zy ; try now elim Zy. revert x. unfold Rdiv. -intros [sx|sx| |sx mx ex Hx] ; - try ( rewrite Rmult_0_l, round_0, Rabs_R0, Rlt_bool_true ; [ split ; apply refl_equal | apply bpow_gt_0 | auto with typeclass_instances ] ). +intros [sx|sx|sx plx|sx mx ex Hx] ; + try ( rewrite Rmult_0_l, round_0, Rabs_R0, Rlt_bool_true ; [ now repeat constructor | apply bpow_gt_0 | auto with typeclass_instances ] ). simpl. case Bdiv_correct_aux. intros H1. unfold Rdiv. case Rlt_bool. -intros (H2, H3). +intros (H2, (H3, H4)). split. now rewrite B2R_FF2B. +split. now rewrite is_finite_FF2B. +rewrite Bsign_FF2B. congruence. intros H2. now rewrite B2FF_FF2B. Qed. @@ -1473,11 +1614,11 @@ Lemma Bsqrt_correct_aux : let '(mz, ez, lz) := Fsqrt_core_binary (Zpos mx) ex in match mz with | Zpos mz => binary_round_aux m false mz ez lz - | _ => F754_nan (* dummy *) + | _ => F754_nan false xH (* dummy *) end in valid_binary z = true /\ FF2R radix2 z = round radix2 fexp (round_mode m) (sqrt x) /\ - is_finite_FF z = true. + is_finite_FF z = true /\ sign_FF z = false. 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). @@ -1578,28 +1719,30 @@ now case mz. apply sqrt_ge_0. Qed. -Definition Bsqrt m x := +Definition Bsqrt sqrt_nan m x := + let f pl := B754_nan (fst pl) (snd pl) in match x with - | B754_nan => x + | B754_nan sx plx => f (sqrt_nan x) | B754_infinity false => x - | B754_infinity true => B754_nan - | B754_finite true _ _ _ => B754_nan + | B754_infinity true => f (sqrt_nan x) + | B754_finite true _ _ _ => f (sqrt_nan x) | B754_zero _ => x | B754_finite sx mx ex Hx => FF2B _ (proj1 (Bsqrt_correct_aux m mx ex Hx)) end. Theorem Bsqrt_correct : - forall m x, - B2R (Bsqrt m x) = round radix2 fexp (round_mode m) (sqrt (B2R x)) /\ - is_finite (Bsqrt m x) = match x with B754_zero _ => true | B754_finite false _ _ _ => true | _ => false end. + forall sqrt_nan m x, + B2R (Bsqrt sqrt_nan m x) = round radix2 fexp (round_mode m) (sqrt (B2R x)) /\ + is_finite (Bsqrt sqrt_nan m x) = match x with B754_zero _ => true | B754_finite false _ _ _ => true | _ => false end /\ + (is_nan (Bsqrt sqrt_nan m x) = false -> Bsign (Bsqrt sqrt_nan m x) = Bsign x). Proof. -intros m [sx|[|]| |sx mx ex Hx] ; try ( now simpl ; rewrite sqrt_0, round_0 ; auto with typeclass_instances ). +intros sqrt_nan m [sx|[|]| |sx mx ex Hx] ; try ( now simpl ; rewrite sqrt_0, round_0 ; intuition auto with typeclass_instances ). simpl. case Bsqrt_correct_aux. -intros H1 (H2, H3). +intros H1 (H2, (H3, H4)). case sx. -refine (conj _ (refl_equal false)). +refine (conj _ (conj (refl_equal false) _)). apply sym_eq. unfold sqrt. case Rcase_abs. @@ -1609,9 +1752,12 @@ auto with typeclass_instances. intros H. elim Rge_not_lt with (1 := H). now apply F2R_lt_0_compat. +easy. split. now rewrite B2R_FF2B. +split. now rewrite is_finite_FF2B. +intro. rewrite Bsign_FF2B. auto. Qed. End Binary. diff --git a/flocq/Appli/Fappli_IEEE_bits.v b/flocq/Appli/Fappli_IEEE_bits.v index 06ed21e..a41fba9 100644 --- a/flocq/Appli/Fappli_IEEE_bits.v +++ b/flocq/Appli/Fappli_IEEE_bits.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2011 Sylvie Boldo +Copyright (C) 2011-2013 Sylvie Boldo #<br /># -Copyright (C) 2011 Guillaume Melquiond +Copyright (C) 2011-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public @@ -172,7 +172,7 @@ Definition bits_of_binary_float (x : binary_float) := match x with | B754_zero sx => join_bits sx 0 0 | B754_infinity sx => join_bits sx 0 (Zpower 2 ew - 1) - | B754_nan => join_bits false (Zpower 2 mw - 1) (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) @@ -184,7 +184,7 @@ Definition split_bits_of_binary_float (x : binary_float) := match x with | B754_zero sx => (sx, 0, 0)%Z | B754_infinity sx => (sx, 0, Zpower 2 ew - 1)%Z - | B754_nan => (false, Zpower 2 mw - 1, 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 @@ -196,8 +196,16 @@ Theorem split_bits_of_binary_float_correct : forall x, split_bits (bits_of_binary_float x) = split_bits_of_binary_float x. Proof. -intros [sx|sx| |sx mx ex Hx] ; +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 Zpower_nat_Z in H0. +eapply Zlt_le_trans. apply H0. +change 2%Z with (radix_val radix2). apply Zpower_le. +rewrite Z.ltb_lt in Hplx. +unfold prec in *. zify; omega. +(* *) 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', _). @@ -246,14 +254,18 @@ Definition binary_float_of_bits_aux x := match mx with | Z0 => F754_zero sx | Zpos px => F754_finite sx px emin - | Zneg _ => F754_nan (* dummy *) + | Zneg _ => F754_nan false xH (* dummy *) end else if Zeq_bool ex (Zpower 2 ew - 1) then - if Zeq_bool mx 0 then F754_infinity sx else F754_nan + match mx with + | Z0 => F754_infinity sx + | Zpos plx => F754_nan sx plx + | Zneg _ => F754_nan false xH (* dummy *) + end else match (mx + Zpower 2 mw)%Z with | Zpos px => F754_finite sx px (ex + emin - 1) - | _ => F754_nan (* dummy *) + | _ => F754_nan false xH (* dummy *) end. Lemma binary_float_of_bits_aux_correct : @@ -292,9 +304,20 @@ cut (0 < emax)%Z. clear -H Hew ; omega. apply (Zpower_gt_0 radix2). clear -Hew ; omega. apply bpow_gt_0. +simpl. intros. rewrite Z.ltb_lt. unfold prec. zify; omega. case Zeq_bool_spec ; intros He2. -now case Zeq_bool. +case_eq (x mod 2 ^ mw)%Z; try easy. +(* nan *) +intros plx Eqplx. apply Z.ltb_lt. +rewrite Z_of_nat_S_digits2_Pnat. +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. +change 2%Z with (radix_val radix2). +apply Z.lt_gt, Zpower_gt_0. omega. +simpl. intros. rewrite Z.ltb_lt. unfold prec. zify; omega. case_eq (x mod 2^mw + 2^mw)%Z ; try easy. +simpl. intros. rewrite Z.ltb_lt. unfold prec. zify; omega. (* normal *) intros px Hm. assert (prec = Zdigits radix2 (Zpos px)). @@ -365,6 +388,7 @@ apply Zlt_gt. apply (Zpower_gt_0 radix2). now apply Zlt_le_weak. apply bpow_gt_0. +simpl. intros. rewrite Z.ltb_lt. unfold prec. zify; omega. Qed. Definition binary_float_of_bits x := @@ -380,7 +404,7 @@ unfold binary_float_of_bits. rewrite B2FF_FF2B. unfold binary_float_of_bits_aux. rewrite split_bits_of_binary_float_correct. -destruct x as [sx|sx| |sx mx ex Bx]. +destruct x as [sx|sx|sx [plx Hplx]|sx mx ex Bx]. apply refl_equal. (* *) simpl. @@ -391,12 +415,7 @@ now apply (Zpower_gt_1 radix2). (* *) simpl. rewrite Zeq_bool_false. -rewrite Zeq_bool_true. -rewrite Zeq_bool_false. -apply refl_equal. -cut (1 < 2^mw)%Z. clear ; omega. -now apply (Zpower_gt_1 radix2). -apply refl_equal. +rewrite Zeq_bool_true; auto. cut (1 < 2^ew)%Z. clear ; omega. now apply (Zpower_gt_1 radix2). (* *) @@ -442,7 +461,6 @@ Qed. Theorem bits_of_binary_float_of_bits : forall x, (0 <= x < 2^(mw+ew+1))%Z -> - binary_float_of_bits x <> B754_nan prec emax -> bits_of_binary_float (binary_float_of_bits x) = x. Proof. intros x Hx. @@ -462,28 +480,28 @@ now apply Zlt_gt. case Zeq_bool_spec ; intros He1. (* subnormal *) case_eq mx. -intros Hm Jx _ _. +intros Hm Jx _. now rewrite He1 in Jx. -intros px Hm Jx _ _. +intros px Hm Jx _. rewrite Zle_bool_false. now rewrite <- He1. now rewrite <- Hm. -intros px Hm _ _ _. +intros px Hm _ _. apply False_ind. apply Zle_not_lt with (1 := proj1 Bm). now rewrite Hm. case Zeq_bool_spec ; intros He2. (* infinity/nan *) -case Zeq_bool_spec ; intros Hm. -now rewrite Hm, He2. -intros _ Cx Nx. -now elim Nx. +case_eq mx; intros Hm. +now rewrite He2. +now rewrite He2. +intros. zify; omega. (* normal *) case_eq (mx + 2 ^ mw)%Z. intros Hm. apply False_ind. clear -Bm Hm ; omega. -intros p Hm Jx Cx _. +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. diff --git a/flocq/Appli/Fappli_rnd_odd.v b/flocq/Appli/Fappli_rnd_odd.v new file mode 100644 index 0000000..b4a2c52 --- /dev/null +++ b/flocq/Appli/Fappli_rnd_odd.v @@ -0,0 +1,979 @@ +(** +This file is part of the Flocq formalization of floating-point +arithmetic in Coq: http://flocq.gforge.inria.fr/ + +Copyright (C) 2010-2013 Sylvie Boldo +#<br /># +Copyright (C) 2010-2013 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 odd and its properties, including the equivalence + between rnd_NE and double rounding with rnd_odd and then rnd_NE *) + + +Require Import Fcore. +Require Import Fcalc_ops. + +Definition Zrnd_odd x := match Req_EM_T x (Z2R (Zfloor x)) with + | left _ => Zfloor x + | right _ => match (Zeven (Zfloor x)) with + | true => Zceil x + | false => Zfloor x + end + end. + + + +Global Instance valid_rnd_odd : Valid_rnd Zrnd_odd. +Proof. +split. +(* . *) +intros x y Hxy. +assert (Zfloor x <= Zrnd_odd y)%Z. +(* .. *) +apply Zle_trans with (Zfloor y). +now apply Zfloor_le. +unfold Zrnd_odd; destruct (Req_EM_T y (Z2R (Zfloor y))). +now apply Zle_refl. +case (Zeven (Zfloor y)). +apply le_Z2R. +apply Rle_trans with y. +apply Zfloor_lb. +apply Zceil_ub. +now apply Zle_refl. +unfold Zrnd_odd at 1. +(* . *) +destruct (Req_EM_T x (Z2R (Zfloor x))) as [Hx|Hx]. +(* .. *) +apply H. +(* .. *) +case_eq (Zeven (Zfloor x)); intros Hx2. +2: apply H. +unfold Zrnd_odd; destruct (Req_EM_T y (Z2R (Zfloor y))) as [Hy|Hy]. +apply Zceil_glb. +now rewrite <- Hy. +case_eq (Zeven (Zfloor y)); intros Hy2. +now apply Zceil_le. +apply Zceil_glb. +assert (H0:(Zfloor x <= Zfloor y)%Z) by now apply Zfloor_le. +case (Zle_lt_or_eq _ _ H0); intros H1. +apply Rle_trans with (1:=Zceil_ub _). +rewrite Zceil_floor_neq. +apply Z2R_le; omega. +now apply sym_not_eq. +contradict Hy2. +rewrite <- H1, Hx2; discriminate. +(* . *) +intros n; unfold Zrnd_odd. +rewrite Zfloor_Z2R, Zceil_Z2R. +destruct (Req_EM_T (Z2R n) (Z2R n)); trivial. +case (Zeven n); trivial. +Qed. + + + +Lemma Zrnd_odd_Zodd: forall x, x <> (Z2R (Zfloor x)) -> + (Zeven (Zrnd_odd x)) = false. +Proof. +intros x Hx; unfold Zrnd_odd. +destruct (Req_EM_T x (Z2R (Zfloor x))) as [H|H]. +now contradict H. +case_eq (Zeven (Zfloor x)). +(* difficult case *) +intros H'. +rewrite Zceil_floor_neq. +rewrite Zeven_plus, H'. +reflexivity. +now apply sym_not_eq. +trivial. +Qed. + + + + +Section Fcore_rnd_odd. + +Variable beta : radix. + +Notation bpow e := (bpow beta e). + +Variable fexp : Z -> Z. + +Context { valid_exp : Valid_exp fexp }. +Context { exists_NE_ : Exists_NE beta fexp }. + +Notation format := (generic_format beta fexp). +Notation canonic := (canonic beta fexp). +Notation cexp := (canonic_exp beta fexp). + + +Definition Rnd_odd_pt (x f : R) := + format f /\ ((f = x)%R \/ + ((Rnd_DN_pt format x f \/ Rnd_UP_pt format x f) /\ + exists g : float beta, f = F2R g /\ canonic g /\ Zeven (Fnum g) = false)). + +Definition Rnd_odd (rnd : R -> R) := + forall x : R, Rnd_odd_pt x (rnd x). + + +Theorem Rnd_odd_pt_sym : forall x f : R, + Rnd_odd_pt (-x) (-f) -> Rnd_odd_pt x f. +Proof with auto with typeclass_instances. +intros x f (H1,H2). +split. +replace f with (-(-f))%R by ring. +now apply generic_format_opp. +destruct H2. +left. +replace f with (-(-f))%R by ring. +rewrite H; ring. +right. +destruct H as (H2,(g,(Hg1,(Hg2,Hg3)))). +split. +destruct H2. +right. +replace f with (-(-f))%R by ring. +replace x with (-(-x))%R by ring. +apply Rnd_DN_UP_pt_sym... +apply generic_format_opp. +left. +replace f with (-(-f))%R by ring. +replace x with (-(-x))%R by ring. +apply Rnd_UP_DN_pt_sym... +apply generic_format_opp. +exists (Float beta (-Fnum g) (Fexp g)). +split. +rewrite F2R_Zopp. +replace f with (-(-f))%R by ring. +rewrite Hg1; reflexivity. +split. +now apply canonic_opp. +simpl. +now rewrite Zeven_opp. +Qed. + + +Theorem round_odd_opp : + forall x, + (round beta fexp Zrnd_odd (-x) = (- round beta fexp Zrnd_odd x))%R. +Proof. +intros x; unfold round. +rewrite <- F2R_Zopp. +unfold F2R; simpl. +apply f_equal2; apply f_equal. +rewrite scaled_mantissa_opp. +generalize (scaled_mantissa beta fexp x); intros r. +unfold Zrnd_odd. +case (Req_EM_T (- r) (Z2R (Zfloor (- r)))). +case (Req_EM_T r (Z2R (Zfloor r))). +intros Y1 Y2. +apply eq_Z2R. +now rewrite Z2R_opp, <- Y1, <-Y2. +intros Y1 Y2. +absurd (r=Z2R (Zfloor r)); trivial. +pattern r at 2; replace r with (-(-r))%R by ring. +rewrite Y2, <- Z2R_opp. +rewrite Zfloor_Z2R. +rewrite Z2R_opp, <- Y2. +ring. +case (Req_EM_T r (Z2R (Zfloor r))). +intros Y1 Y2. +absurd (-r=Z2R (Zfloor (-r)))%R; trivial. +pattern r at 2; rewrite Y1. +rewrite <- Z2R_opp, Zfloor_Z2R. +now rewrite Z2R_opp, <- Y1. +intros Y1 Y2. +unfold Zceil; rewrite Ropp_involutive. +replace (Zeven (Zfloor (- r))) with (negb (Zeven (Zfloor r))). +case (Zeven (Zfloor r)); simpl; ring. +apply trans_eq with (Zeven (Zceil r)). +rewrite Zceil_floor_neq. +rewrite Zeven_plus. +simpl; reflexivity. +now apply sym_not_eq. +rewrite <- (Zeven_opp (Zfloor (- r))). +reflexivity. +apply canonic_exp_opp. +Qed. + + + +Theorem round_odd_pt : + forall x, + Rnd_odd_pt x (round beta fexp Zrnd_odd x). +Proof with auto with typeclass_instances. +cut (forall x : R, (0 < x)%R -> Rnd_odd_pt x (round beta fexp Zrnd_odd x)). +intros H x; case (Rle_or_lt 0 x). +intros H1; destruct H1. +now apply H. +rewrite <- H0. +rewrite round_0... +split. +apply generic_format_0. +now left. +intros Hx; apply Rnd_odd_pt_sym. +rewrite <- round_odd_opp. +apply H. +auto with real. +(* *) +intros x Hxp. +generalize (generic_format_round beta fexp Zrnd_odd x). +set (o:=round beta fexp Zrnd_odd x). +intros Ho. +split. +assumption. +(* *) +case (Req_dec o x); intros Hx. +now left. +right. +assert (o=round beta fexp Zfloor x \/ o=round beta fexp Zceil x). +unfold o, round, F2R;simpl. +case (Zrnd_DN_or_UP Zrnd_odd (scaled_mantissa beta fexp x))... +intros H; rewrite H; now left. +intros H; rewrite H; now right. +split. +destruct H; rewrite H. +left; apply round_DN_pt... +right; apply round_UP_pt... +(* *) +unfold o, Zrnd_odd, round. +case (Req_EM_T (scaled_mantissa beta fexp x) + (Z2R (Zfloor (scaled_mantissa beta fexp x)))). +intros T. +absurd (o=x); trivial. +apply round_generic... +unfold generic_format, F2R; simpl. +rewrite <- (scaled_mantissa_mult_bpow beta fexp) at 1. +apply f_equal2; trivial; rewrite T at 1. +apply f_equal, sym_eq, Ztrunc_floor. +apply Rmult_le_pos. +now left. +apply bpow_ge_0. +intros L. +case_eq (Zeven (Zfloor (scaled_mantissa beta fexp x))). +(* . *) +generalize (generic_format_round beta fexp Zceil x). +unfold generic_format. +set (f:=round beta fexp Zceil x). +set (ef := canonic_exp beta fexp f). +set (mf := Ztrunc (scaled_mantissa beta fexp f)). +exists (Float beta mf ef). +unfold Fcore_generic_fmt.canonic. +rewrite <- H0. +repeat split; try assumption. +apply trans_eq with (negb (Zeven (Zfloor (scaled_mantissa beta fexp x)))). +2: rewrite H1; reflexivity. +apply trans_eq with (negb (Zeven (Fnum + (Float beta (Zfloor (scaled_mantissa beta fexp x)) (cexp x))))). +2: reflexivity. +case (Rle_lt_or_eq_dec 0 (round beta fexp Zfloor x)). +rewrite <- round_0 with beta fexp Zfloor... +apply round_le... +now left. +intros Y. +generalize (DN_UP_parity_generic beta fexp)... +unfold DN_UP_parity_prop. +intros T; apply T with x; clear T. +unfold generic_format. +rewrite <- (scaled_mantissa_mult_bpow beta fexp x) at 1. +unfold F2R; simpl. +apply Rmult_neq_compat_r. +apply Rgt_not_eq, bpow_gt_0. +rewrite Ztrunc_floor. +assumption. +apply Rmult_le_pos. +now left. +apply bpow_ge_0. +unfold Fcore_generic_fmt.canonic. +simpl. +apply sym_eq, canonic_exp_DN... +unfold Fcore_generic_fmt.canonic. +rewrite <- H0; reflexivity. +reflexivity. +apply trans_eq with (round beta fexp Ztrunc (round beta fexp Zceil x)). +reflexivity. +apply round_generic... +intros Y. +replace (Fnum {| Fnum := Zfloor (scaled_mantissa beta fexp x); Fexp := cexp x |}) + with (Fnum (Float beta 0 (fexp (ln_beta beta 0)))). +generalize (DN_UP_parity_generic beta fexp)... +unfold DN_UP_parity_prop. +intros T; apply T with x; clear T. +unfold generic_format. +rewrite <- (scaled_mantissa_mult_bpow beta fexp x) at 1. +unfold F2R; simpl. +apply Rmult_neq_compat_r. +apply Rgt_not_eq, bpow_gt_0. +rewrite Ztrunc_floor. +assumption. +apply Rmult_le_pos. +now left. +apply bpow_ge_0. +apply canonic_0. +unfold Fcore_generic_fmt.canonic. +rewrite <- H0; reflexivity. +rewrite <- Y; unfold F2R; simpl; ring. +apply trans_eq with (round beta fexp Ztrunc (round beta fexp Zceil x)). +reflexivity. +apply round_generic... +simpl. +apply eq_Z2R, Rmult_eq_reg_r with (bpow (cexp x)). +unfold round, F2R in Y; simpl in Y; rewrite <- Y. +simpl; ring. +apply Rgt_not_eq, bpow_gt_0. +(* . *) +intros Y. +case (Rle_lt_or_eq_dec 0 (round beta fexp Zfloor x)). +rewrite <- round_0 with beta fexp Zfloor... +apply round_le... +now left. +intros Hrx. +set (ef := canonic_exp beta fexp x). +set (mf := Zfloor (scaled_mantissa beta fexp x)). +exists (Float beta mf ef). +unfold Fcore_generic_fmt.canonic. +repeat split; try assumption. +simpl. +apply trans_eq with (cexp (round beta fexp Zfloor x )). +apply sym_eq, canonic_exp_DN... +reflexivity. +intros Hrx; contradict Y. +replace (Zfloor (scaled_mantissa beta fexp x)) with 0%Z. +simpl; discriminate. +apply eq_Z2R, Rmult_eq_reg_r with (bpow (cexp x)). +unfold round, F2R in Hrx; simpl in Hrx; rewrite <- Hrx. +simpl; ring. +apply Rgt_not_eq, bpow_gt_0. +Qed. + +End Fcore_rnd_odd. + +Section Odd_prop_aux. + +Variable beta : radix. +Hypothesis Even_beta: Zeven (radix_val beta)=true. + +Notation bpow e := (bpow beta e). + +Variable fexp : Z -> Z. +Variable fexpe : Z -> Z. + +Context { valid_exp : Valid_exp fexp }. +Context { exists_NE_ : Exists_NE beta fexp }. (* for underflow reason *) +Context { valid_expe : Valid_exp fexpe }. +Context { exists_NE_e : Exists_NE beta fexpe }. (* for defining rounding to odd *) + +Hypothesis fexpe_fexp: forall e, (fexpe e <= fexp e -2)%Z. + + +Lemma generic_format_fexpe_fexp: forall x, + generic_format beta fexp x -> generic_format beta fexpe x. +Proof. +intros x Hx. +apply generic_inclusion_ln_beta with fexp; trivial; intros Hx2. +generalize (fexpe_fexp (ln_beta beta x)). +omega. +Qed. + + + +Lemma exists_even_fexp_lt: forall (c:Z->Z), forall (x:R), + (exists f:float beta, F2R f = x /\ (c (ln_beta beta x) < Fexp f)%Z) -> + exists f:float beta, F2R f =x /\ canonic beta c f /\ Zeven (Fnum f) = true. +Proof with auto with typeclass_instances. +intros c x (g,(Hg1,Hg2)). +exists (Float beta + (Fnum g*Z.pow (radix_val beta) (Fexp g - c (ln_beta beta x))) + (c (ln_beta beta x))). +assert (F2R (Float beta + (Fnum g*Z.pow (radix_val beta) (Fexp g - c (ln_beta beta x))) + (c (ln_beta beta x))) = x). +unfold F2R; simpl. +rewrite Z2R_mult, Z2R_Zpower. +rewrite Rmult_assoc, <- bpow_plus. +rewrite <- Hg1; unfold F2R. +apply f_equal, f_equal. +ring. +omega. +split; trivial. +split. +unfold canonic, canonic_exp. +now rewrite H. +simpl. +rewrite Zeven_mult. +rewrite Zeven_Zpower. +rewrite Even_beta. +apply Bool.orb_true_intro. +now right. +omega. +Qed. + + +Variable choice:Z->bool. +Variable x:R. + + +Variable d u: float beta. +Hypothesis Hd: Rnd_DN_pt (generic_format beta fexp) x (F2R d). +Hypothesis Cd: canonic beta fexp d. +Hypothesis Hu: Rnd_UP_pt (generic_format beta fexp) x (F2R u). +Hypothesis Cu: canonic beta fexp u. + +Hypothesis xPos: (0 < x)%R. + + +Let m:= ((F2R d+F2R u)/2)%R. + + +Lemma d_eq: F2R d= round beta fexp Zfloor x. +Proof with auto with typeclass_instances. +apply Rnd_DN_pt_unicity with (generic_format beta fexp) x... +apply round_DN_pt... +Qed. + + +Lemma u_eq: F2R u= round beta fexp Zceil x. +Proof with auto with typeclass_instances. +apply Rnd_UP_pt_unicity with (generic_format beta fexp) x... +apply round_UP_pt... +Qed. + + +Lemma d_ge_0: (0 <= F2R d)%R. +Proof with auto with typeclass_instances. +rewrite d_eq; apply round_ge_generic... +apply generic_format_0. +now left. +Qed. + + + +Lemma ln_beta_d: (0< F2R d)%R -> + (ln_beta beta (F2R d) = ln_beta beta x :>Z). +Proof with auto with typeclass_instances. +intros Y. +rewrite d_eq; apply ln_beta_round_DN... +now rewrite <- d_eq. +Qed. + + +Lemma Fexp_d: (0 < F2R d)%R -> Fexp d =fexp (ln_beta beta x). +Proof with auto with typeclass_instances. +intros Y. +now rewrite Cd, <- ln_beta_d. +Qed. + + + +Lemma format_bpow_x: (0 < F2R d)%R + -> generic_format beta fexp (bpow (ln_beta beta x)). +Proof with auto with typeclass_instances. +intros Y. +apply generic_format_bpow. +apply valid_exp. +rewrite <- Fexp_d; trivial. +apply Zlt_le_trans with (ln_beta beta (F2R d))%Z. +rewrite Cd; apply ln_beta_generic_gt... +now apply Rgt_not_eq. +apply Hd. +apply ln_beta_le; trivial. +apply Hd. +Qed. + + +Lemma format_bpow_d: (0 < F2R d)%R -> + generic_format beta fexp (bpow (ln_beta beta (F2R d))). +Proof with auto with typeclass_instances. +intros Y; apply generic_format_bpow. +apply valid_exp. +apply ln_beta_generic_gt... +now apply Rgt_not_eq. +now apply generic_format_canonic. +Qed. + + +Lemma d_le_m: (F2R d <= m)%R. +apply Rmult_le_reg_l with 2%R. +auto with real. +apply Rplus_le_reg_l with (-F2R d)%R. +apply Rle_trans with (F2R d). +right; ring. +apply Rle_trans with (F2R u). +apply Rle_trans with x. +apply Hd. +apply Hu. +right; unfold m; field. +Qed. + +Lemma m_le_u: (m <= F2R u)%R. +apply Rmult_le_reg_l with 2%R. +auto with real. +apply Rplus_le_reg_l with (-F2R u)%R. +apply Rle_trans with (F2R d). +right; unfold m; field. +apply Rle_trans with (F2R u). +apply Rle_trans with x. +apply Hd. +apply Hu. +right; ring. +Qed. + +Lemma ln_beta_m: (0 < F2R d)%R -> (ln_beta beta m =ln_beta beta (F2R d) :>Z). +Proof with auto with typeclass_instances. +intros dPos; apply ln_beta_unique_pos. +split. +apply Rle_trans with (F2R d). +destruct (ln_beta beta (F2R d)) as (e,He). +simpl. +rewrite Rabs_right in He. +apply He. +now apply Rgt_not_eq. +apply Rle_ge; now left. +apply d_le_m. +case m_le_u; intros H. +apply Rlt_le_trans with (1:=H). +rewrite u_eq. +apply round_le_generic... +apply generic_format_bpow. +apply valid_exp. +apply ln_beta_generic_gt... +now apply Rgt_not_eq. +now apply generic_format_canonic. +case (Rle_or_lt x (bpow (ln_beta beta (F2R d)))); trivial; intros Z. +absurd ((bpow (ln_beta beta (F2R d)) <= (F2R d)))%R. +apply Rlt_not_le. +destruct (ln_beta beta (F2R d)) as (e,He). +simpl in *; rewrite Rabs_right in He. +apply He. +now apply Rgt_not_eq. +apply Rle_ge; now left. +apply Rle_trans with (round beta fexp Zfloor x). +2: right; apply sym_eq, d_eq. +apply round_ge_generic... +apply generic_format_bpow. +apply valid_exp. +apply ln_beta_generic_gt... +now apply Rgt_not_eq. +now apply generic_format_canonic. +now left. +replace m with (F2R d). +destruct (ln_beta beta (F2R d)) as (e,He). +simpl in *; rewrite Rabs_right in He. +apply He. +now apply Rgt_not_eq. +apply Rle_ge; now left. +assert (F2R d = F2R u). +apply Rmult_eq_reg_l with (/2)%R. +apply Rplus_eq_reg_l with (/2*F2R u)%R. +apply trans_eq with m. +unfold m, Rdiv; ring. +rewrite H; field. +auto with real. +apply Rgt_not_eq, Rlt_gt; auto with real. +unfold m; rewrite <- H0; field. +Qed. + + +Lemma ln_beta_m_0: (0 = F2R d)%R + -> (ln_beta beta m =ln_beta beta (F2R u)-1:>Z)%Z. +Proof with auto with typeclass_instances. +intros Y. +apply ln_beta_unique_pos. +unfold m; rewrite <- Y, Rplus_0_l. +rewrite u_eq. +destruct (ln_beta beta x) as (e,He). +rewrite Rabs_right in He. +rewrite round_UP_small_pos with (ex:=e). +rewrite ln_beta_bpow. +ring_simplify (fexp e + 1 - 1)%Z. +split. +unfold Zminus; rewrite bpow_plus. +unfold Rdiv; apply Rmult_le_compat_l. +apply bpow_ge_0. +simpl; unfold Z.pow_pos; simpl. +rewrite Zmult_1_r; apply Rinv_le. +auto with real. +apply (Z2R_le 2). +specialize (radix_gt_1 beta). +omega. +apply Rlt_le_trans with (bpow (fexp e)*1)%R. +2: right; ring. +unfold Rdiv; apply Rmult_lt_compat_l. +apply bpow_gt_0. +rewrite <- Rinv_1 at 3. +apply Rinv_lt; auto with real. +now apply He, Rgt_not_eq. +apply exp_small_round_0_pos with beta (Zfloor) x... +now apply He, Rgt_not_eq. +now rewrite <- d_eq, Y. +now left. +Qed. + + + + + +Lemma u'_eq: (0 < F2R d)%R -> exists f:float beta, F2R f = F2R u /\ (Fexp f = Fexp d)%Z. +Proof with auto with typeclass_instances. +intros Y. +rewrite u_eq; unfold round. +eexists; repeat split. +simpl; now rewrite Fexp_d. +Qed. + + + + +Lemma m_eq: (0 < F2R d)%R -> exists f:float beta, + F2R f = m /\ (Fexp f = fexp (ln_beta beta x) -1)%Z. +Proof with auto with typeclass_instances. +intros Y. +specialize (Zeven_ex (radix_val beta)); rewrite Even_beta. +intros (b, Hb); rewrite Zplus_0_r in Hb. +destruct u'_eq as (u', (Hu'1,Hu'2)); trivial. +exists (Fmult beta (Float beta b (-1)) (Fplus beta d u'))%R. +split. +rewrite F2R_mult, F2R_plus, Hu'1. +unfold m; rewrite Rmult_comm. +unfold Rdiv; apply f_equal. +unfold F2R; simpl; unfold Z.pow_pos; simpl. +rewrite Zmult_1_r, Hb, Z2R_mult. +simpl; field. +apply Rgt_not_eq, Rmult_lt_reg_l with (Z2R 2). +simpl; auto with real. +rewrite Rmult_0_r, <-Z2R_mult, <-Hb. +apply radix_pos. +apply trans_eq with (-1+Fexp (Fplus beta d u'))%Z. +unfold Fmult. +destruct (Fplus beta d u'); reflexivity. +rewrite Zplus_comm; unfold Zminus; apply f_equal2. +2: reflexivity. +rewrite Fexp_Fplus. +rewrite Z.min_l. +now rewrite Fexp_d. +rewrite Hu'2; omega. +Qed. + +Lemma m_eq_0: (0 = F2R d)%R -> exists f:float beta, + F2R f = m /\ (Fexp f = fexp (ln_beta beta (F2R u)) -1)%Z. +Proof with auto with typeclass_instances. +intros Y. +specialize (Zeven_ex (radix_val beta)); rewrite Even_beta. +intros (b, Hb); rewrite Zplus_0_r in Hb. +exists (Fmult beta (Float beta b (-1)) u)%R. +split. +rewrite F2R_mult; unfold m; rewrite <- Y, Rplus_0_l. +rewrite Rmult_comm. +unfold Rdiv; apply f_equal. +unfold F2R; simpl; unfold Z.pow_pos; simpl. +rewrite Zmult_1_r, Hb, Z2R_mult. +simpl; field. +apply Rgt_not_eq, Rmult_lt_reg_l with (Z2R 2). +simpl; auto with real. +rewrite Rmult_0_r, <-Z2R_mult, <-Hb. +apply radix_pos. +apply trans_eq with (-1+Fexp u)%Z. +unfold Fmult. +destruct u; reflexivity. +rewrite Zplus_comm, Cu; unfold Zminus; now apply f_equal2. +Qed. + +Lemma fexp_m_eq_0: (0 = F2R d)%R -> + (fexp (ln_beta beta (F2R u)-1) < fexp (ln_beta beta (F2R u))+1)%Z. +Proof with auto with typeclass_instances. +intros Y. +assert ((fexp (ln_beta beta (F2R u) - 1) <= fexp (ln_beta beta (F2R u))))%Z. +2: omega. +destruct (ln_beta beta x) as (e,He). +rewrite Rabs_right in He. +2: now left. +assert (e <= fexp e)%Z. +apply exp_small_round_0_pos with beta (Zfloor) x... +now apply He, Rgt_not_eq. +now rewrite <- d_eq, Y. +rewrite u_eq, round_UP_small_pos with (ex:=e); trivial. +2: now apply He, Rgt_not_eq. +rewrite ln_beta_bpow. +ring_simplify (fexp e + 1 - 1)%Z. +replace (fexp (fexp e)) with (fexp e). +case exists_NE_; intros V. +contradict V; rewrite Even_beta; discriminate. +rewrite (proj2 (V e)); omega. +apply sym_eq, valid_exp; omega. +Qed. + +Lemma Fm: generic_format beta fexpe m. +case (d_ge_0); intros Y. +(* *) +destruct m_eq as (g,(Hg1,Hg2)); trivial. +apply generic_format_F2R' with g. +now apply sym_eq. +intros H; unfold canonic_exp; rewrite Hg2. +rewrite ln_beta_m; trivial. +rewrite <- Fexp_d; trivial. +rewrite Cd. +unfold canonic_exp. +generalize (fexpe_fexp (ln_beta beta (F2R d))). +omega. +(* *) +destruct m_eq_0 as (g,(Hg1,Hg2)); trivial. +apply generic_format_F2R' with g. +assumption. +intros H; unfold canonic_exp; rewrite Hg2. +rewrite ln_beta_m_0; try assumption. +apply Zle_trans with (1:=fexpe_fexp _). +assert (fexp (ln_beta beta (F2R u)-1) < fexp (ln_beta beta (F2R u))+1)%Z;[idtac|omega]. +now apply fexp_m_eq_0. +Qed. + + + +Lemma Zm: + exists g : float beta, F2R g = m /\ canonic beta fexpe g /\ Zeven (Fnum g) = true. +Proof with auto with typeclass_instances. +case (d_ge_0); intros Y. +(* *) +destruct m_eq as (g,(Hg1,Hg2)); trivial. +apply exists_even_fexp_lt. +exists g; split; trivial. +rewrite Hg2. +rewrite ln_beta_m; trivial. +rewrite <- Fexp_d; trivial. +rewrite Cd. +unfold canonic_exp. +generalize (fexpe_fexp (ln_beta beta (F2R d))). +omega. +(* *) +destruct m_eq_0 as (g,(Hg1,Hg2)); trivial. +apply exists_even_fexp_lt. +exists g; split; trivial. +rewrite Hg2. +rewrite ln_beta_m_0; trivial. +apply Zle_lt_trans with (1:=fexpe_fexp _). +assert (fexp (ln_beta beta (F2R u)-1) < fexp (ln_beta beta (F2R u))+1)%Z;[idtac|omega]. +now apply fexp_m_eq_0. +Qed. + + +Lemma DN_odd_d_aux: forall z, (F2R d<= z< F2R u)%R -> + Rnd_DN_pt (generic_format beta fexp) z (F2R d). +Proof with auto with typeclass_instances. +intros z Hz1. +replace (F2R d) with (round beta fexp Zfloor z). +apply round_DN_pt... +case (Rnd_DN_UP_pt_split _ _ _ _ Hd Hu (round beta fexp Zfloor z)). +apply generic_format_round... +intros Y; apply Rle_antisym; trivial. +apply round_DN_pt... +apply Hd. +apply Hz1. +intros Y; absurd (z < z)%R. +auto with real. +apply Rlt_le_trans with (1:=proj2 Hz1), Rle_trans with (1:=Y). +apply round_DN_pt... +Qed. + +Lemma UP_odd_d_aux: forall z, (F2R d< z <= F2R u)%R -> + Rnd_UP_pt (generic_format beta fexp) z (F2R u). +Proof with auto with typeclass_instances. +intros z Hz1. +replace (F2R u) with (round beta fexp Zceil z). +apply round_UP_pt... +case (Rnd_DN_UP_pt_split _ _ _ _ Hd Hu (round beta fexp Zceil z)). +apply generic_format_round... +intros Y; absurd (z < z)%R. +auto with real. +apply Rle_lt_trans with (2:=proj1 Hz1), Rle_trans with (2:=Y). +apply round_UP_pt... +intros Y; apply Rle_antisym; trivial. +apply round_UP_pt... +apply Hu. +apply Hz1. +Qed. + + +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. +set (o:=round beta fexpe Zrnd_odd x). +case (generic_format_EM beta fexp x); intros Hx. +replace o with x; trivial. +unfold o; apply sym_eq, round_generic... +now apply generic_format_fexpe_fexp. +assert (K1:(F2R d <= o)%R). +apply round_ge_generic... +apply generic_format_fexpe_fexp, Hd. +apply Hd. +assert (K2:(o <= F2R u)%R). +apply round_le_generic... +apply generic_format_fexpe_fexp, Hu. +apply Hu. +assert (P:(x <> m -> o=m -> (forall P:Prop, P))). +intros Y1 Y2. +assert (H:(Rnd_odd_pt beta fexpe x o)). +apply round_odd_pt... +destruct H as (_,H); destruct H. +absurd (x=m)%R; try trivial. +now rewrite <- Y2, H. +destruct H as (_,(k,(Hk1,(Hk2,Hk3)))). +destruct Zm as (k',(Hk'1,(Hk'2,Hk'3))). +absurd (true=false). +discriminate. +rewrite <- Hk3, <- Hk'3. +apply f_equal, f_equal. +apply canonic_unicity with fexpe... +now rewrite Hk'1, <- Y2. +assert (generic_format beta fexp o -> (forall P:Prop, P)). +intros Y. +assert (H:(Rnd_odd_pt beta fexpe x o)). +apply round_odd_pt... +destruct H as (_,H); destruct H. +absurd (generic_format beta fexp x); trivial. +now rewrite <- H. +destruct H as (_,(k,(Hk1,(Hk2,Hk3)))). +destruct (exists_even_fexp_lt fexpe o) as (k',(Hk'1,(Hk'2,Hk'3))). +eexists; split. +apply sym_eq, Y. +simpl; unfold canonic_exp. +apply Zle_lt_trans with (1:=fexpe_fexp _). +omega. +absurd (true=false). +discriminate. +rewrite <- Hk3, <- Hk'3. +apply f_equal, f_equal. +apply canonic_unicity with fexpe... +now rewrite Hk'1, <- Hk1. +case K1; clear K1; intros K1. +2: apply H; rewrite <- K1; apply Hd. +case K2; clear K2; intros K2. +2: apply H; rewrite K2; apply Hu. +case (Rle_or_lt x m); intros Y;[destruct Y|idtac]. +(* . *) +apply trans_eq with (F2R d). +apply round_N_DN_betw with (F2R u)... +apply DN_odd_d_aux; split; try left; assumption. +apply UP_odd_d_aux; split; try left; assumption. +split. +apply round_ge_generic... +apply generic_format_fexpe_fexp, Hd. +apply Hd. +assert (o <= (F2R d + F2R u) / 2)%R. +apply round_le_generic... +apply Fm. +now left. +destruct H1; trivial. +apply P. +now apply Rlt_not_eq. +trivial. +apply sym_eq, round_N_DN_betw with (F2R u)... +split. +apply Hd. +exact H0. +(* . *) +replace o with x. +reflexivity. +apply sym_eq, round_generic... +rewrite H0; apply Fm. +(* . *) +apply trans_eq with (F2R u). +apply round_N_UP_betw with (F2R d)... +apply DN_odd_d_aux; split; try left; assumption. +apply UP_odd_d_aux; split; try left; assumption. +split. +assert ((F2R d + F2R u) / 2 <= o)%R. +apply round_ge_generic... +apply Fm. +now left. +destruct H0; trivial. +apply P. +now apply Rgt_not_eq. +rewrite <- H0; trivial. +apply round_le_generic... +apply generic_format_fexpe_fexp, Hu. +apply Hu. +apply sym_eq, round_N_UP_betw with (F2R d)... +split. +exact Y. +apply Hu. +Qed. + + +End Odd_prop_aux. + +Section Odd_prop. + +Variable beta : radix. +Hypothesis Even_beta: Zeven (radix_val beta)=true. + +Variable fexp : Z -> Z. +Variable fexpe : Z -> Z. +Variable choice:Z->bool. + +Context { valid_exp : Valid_exp fexp }. +Context { exists_NE_ : Exists_NE beta fexp }. (* for underflow reason *) +Context { valid_expe : Valid_exp fexpe }. +Context { exists_NE_e : Exists_NE beta fexpe }. (* for defining rounding to odd *) + +Hypothesis fexpe_fexp: forall e, (fexpe e <= fexp e -2)%Z. + + +Theorem canonizer: forall f, generic_format beta fexp f + -> exists g : float beta, f = F2R g /\ canonic beta fexp g. +Proof with auto with typeclass_instances. +intros f Hf. +exists (Float beta (Ztrunc (scaled_mantissa beta fexp f)) (canonic_exp beta fexp f)). +assert (L:(f = F2R (Float beta (Ztrunc (scaled_mantissa beta fexp f)) (canonic_exp beta fexp f)))). +apply trans_eq with (round beta fexp Ztrunc f). +apply sym_eq, round_generic... +reflexivity. +split; trivial. +unfold canonic; rewrite <- L. +reflexivity. +Qed. + + + + +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. +intros x. +case (total_order_T x 0); intros H; [case H; clear H; intros H | idtac]. +rewrite <- (Ropp_involutive x). +rewrite round_odd_opp. +rewrite 2!round_N_opp. +apply f_equal. +destruct (canonizer (round beta fexp Zfloor (-x))) as (d,(Hd1,Hd2)). +apply generic_format_round... +destruct (canonizer (round beta fexp Zceil (-x))) as (u,(Hu1,Hu2)). +apply generic_format_round... +apply round_odd_prop_pos with d u... +rewrite <- Hd1; apply round_DN_pt... +rewrite <- Hu1; apply round_UP_pt... +auto with real. +(* . *) +rewrite H; repeat rewrite round_0... +(* . *) +destruct (canonizer (round beta fexp Zfloor x)) as (d,(Hd1,Hd2)). +apply generic_format_round... +destruct (canonizer (round beta fexp Zceil x)) as (u,(Hu1,Hu2)). +apply generic_format_round... +apply round_odd_prop_pos with d u... +rewrite <- Hd1; apply round_DN_pt... +rewrite <- Hu1; apply round_UP_pt... +Qed. + + +End Odd_prop. diff --git a/flocq/Calc/Fcalc_bracket.v b/flocq/Calc/Fcalc_bracket.v index dd4bd97..90a8588 100644 --- a/flocq/Calc/Fcalc_bracket.v +++ b/flocq/Calc/Fcalc_bracket.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/flocq/Calc/Fcalc_digits.v b/flocq/Calc/Fcalc_digits.v index 6210bac..4f76cc2 100644 --- a/flocq/Calc/Fcalc_digits.v +++ b/flocq/Calc/Fcalc_digits.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/flocq/Calc/Fcalc_div.v b/flocq/Calc/Fcalc_div.v index 6594e55..c8f1f9f 100644 --- a/flocq/Calc/Fcalc_div.v +++ b/flocq/Calc/Fcalc_div.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/flocq/Calc/Fcalc_ops.v b/flocq/Calc/Fcalc_ops.v index 15bb211..7ece683 100644 --- a/flocq/Calc/Fcalc_ops.v +++ b/flocq/Calc/Fcalc_ops.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/flocq/Calc/Fcalc_round.v b/flocq/Calc/Fcalc_round.v index 3d31aea..a1bcb84 100644 --- a/flocq/Calc/Fcalc_round.v +++ b/flocq/Calc/Fcalc_round.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/flocq/Calc/Fcalc_sqrt.v b/flocq/Calc/Fcalc_sqrt.v index e6fe74b..2ed3234 100644 --- a/flocq/Calc/Fcalc_sqrt.v +++ b/flocq/Calc/Fcalc_sqrt.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/flocq/Core/Fcore.v b/flocq/Core/Fcore.v index 23ebb39..2a5a5f0 100644 --- a/flocq/Core/Fcore.v +++ b/flocq/Core/Fcore.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/flocq/Core/Fcore_FIX.v b/flocq/Core/Fcore_FIX.v index f185ddf..a3b8d4d 100644 --- a/flocq/Core/Fcore_FIX.v +++ b/flocq/Core/Fcore_FIX.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/flocq/Core/Fcore_FLT.v b/flocq/Core/Fcore_FLT.v index 4ad4797..c057b6c 100644 --- a/flocq/Core/Fcore_FLT.v +++ b/flocq/Core/Fcore_FLT.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/flocq/Core/Fcore_FLX.v b/flocq/Core/Fcore_FLX.v index 62ecb6f..800540f 100644 --- a/flocq/Core/Fcore_FLX.v +++ b/flocq/Core/Fcore_FLX.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/flocq/Core/Fcore_FTZ.v b/flocq/Core/Fcore_FTZ.v index 5356c11..5f3e533 100644 --- a/flocq/Core/Fcore_FTZ.v +++ b/flocq/Core/Fcore_FTZ.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/flocq/Core/Fcore_Raux.v b/flocq/Core/Fcore_Raux.v index 748e36e..d811019 100644 --- a/flocq/Core/Fcore_Raux.v +++ b/flocq/Core/Fcore_Raux.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public @@ -121,6 +121,18 @@ rewrite <- 3!(Rmult_comm r). apply Rmult_minus_distr_l. Qed. +Lemma Rmult_neq_reg_r: forall r1 r2 r3:R, (r2 * r1 <> r3 * r1)%R -> r2 <> r3. +intros r1 r2 r3 H1 H2. +apply H1; rewrite H2; ring. +Qed. + +Lemma Rmult_neq_compat_r: forall r1 r2 r3:R, (r1 <> 0)%R -> (r2 <> r3)%R + -> (r2 *r1 <> r3*r1)%R. +intros r1 r2 r3 H H1 H2. +now apply H1, Rmult_eq_reg_r with r1. +Qed. + + Theorem Rmult_min_distr_r : forall r r1 r2 : R, (0 <= r)%R -> diff --git a/flocq/Core/Fcore_Zaux.v b/flocq/Core/Fcore_Zaux.v index af0d837..7ba28ca 100644 --- a/flocq/Core/Fcore_Zaux.v +++ b/flocq/Core/Fcore_Zaux.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2011 Sylvie Boldo +Copyright (C) 2011-2013 Sylvie Boldo #<br /># -Copyright (C) 2011 Guillaume Melquiond +Copyright (C) 2011-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/flocq/Core/Fcore_defs.v b/flocq/Core/Fcore_defs.v index fda3a85..ad8cf4f 100644 --- a/flocq/Core/Fcore_defs.v +++ b/flocq/Core/Fcore_defs.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/flocq/Core/Fcore_digits.v b/flocq/Core/Fcore_digits.v index 2ae076e..02c7a0e 100644 --- a/flocq/Core/Fcore_digits.v +++ b/flocq/Core/Fcore_digits.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2011 Sylvie Boldo +Copyright (C) 2011-2013 Sylvie Boldo #<br /># -Copyright (C) 2011 Guillaume Melquiond +Copyright (C) 2011-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/flocq/Core/Fcore_float_prop.v b/flocq/Core/Fcore_float_prop.v index 746f7a6..e1535bc 100644 --- a/flocq/Core/Fcore_float_prop.v +++ b/flocq/Core/Fcore_float_prop.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/flocq/Core/Fcore_generic_fmt.v b/flocq/Core/Fcore_generic_fmt.v index b1db47c..b04cf3d 100644 --- a/flocq/Core/Fcore_generic_fmt.v +++ b/flocq/Core/Fcore_generic_fmt.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public @@ -165,6 +165,22 @@ now rewrite Ztrunc_Z2R. now apply Zle_left. Qed. +Lemma generic_format_F2R': forall (x:R) (f:float beta), + F2R f = x -> ((x <> 0)%R -> + (canonic_exp x <= Fexp f)%Z) -> + generic_format x. +Proof. +intros x f H1 H2. +rewrite <- H1; destruct f as (m,e). +apply generic_format_F2R. +simpl in *; intros H3. +rewrite H1; apply H2. +intros Y; apply H3. +apply F2R_eq_0_reg with beta e. +now rewrite H1. +Qed. + + Theorem canonic_opp : forall m e, canonic (Float beta m e) -> @@ -175,6 +191,26 @@ unfold canonic. now rewrite F2R_Zopp, canonic_exp_opp. Qed. +Theorem canonic_abs : + forall m e, + canonic (Float beta m e) -> + canonic (Float beta (Zabs m) e). +Proof. +intros m e H. +unfold canonic. +now rewrite F2R_Zabs, canonic_exp_abs. +Qed. + +Theorem canonic_0: canonic (Float beta 0 (fexp (ln_beta beta 0%R))). +Proof. +unfold canonic; simpl; unfold canonic_exp. +replace (F2R {| Fnum := 0; Fexp := fexp (ln_beta beta 0) |}) with 0%R. +reflexivity. +unfold F2R; simpl; ring. +Qed. + + + Theorem canonic_unicity : forall f1 f2, canonic f1 -> @@ -756,6 +792,24 @@ refine (let H := _ in conj (proj1 H) (Rlt_le _ _ (proj2 H))). now apply mantissa_small_pos. Qed. + +Theorem exp_small_round_0_pos : + forall x ex, + (bpow (ex - 1) <= x < bpow ex)%R -> + round x =R0 -> (ex <= fexp ex)%Z . +Proof. +intros x ex H H1. +case (Zle_or_lt ex (fexp ex)); trivial; intros V. +contradict H1. +apply Rgt_not_eq. +apply Rlt_le_trans with (bpow (ex-1)). +apply bpow_gt_0. +apply (round_bounded_large_pos); assumption. +Qed. + + + + Theorem generic_format_round_pos : forall x, (0 < x)%R -> @@ -1014,6 +1068,24 @@ intros rnd' Hr x. apply round_bounded_large_pos... Qed. +Theorem exp_small_round_0 : + forall rnd {Hr : Valid_rnd rnd} x ex, + (bpow (ex - 1) <= Rabs x < bpow ex)%R -> + round rnd x =R0 -> (ex <= fexp ex)%Z . +Proof. +intros rnd Hr x ex H1 H2. +generalize Rabs_R0. +rewrite <- H2 at 1. +apply (round_abs_abs' (fun t rt => forall (ex : Z), +(bpow (ex - 1) <= t < bpow ex)%R -> +rt = 0%R -> (ex <= fexp ex)%Z)); trivial. +clear; intros rnd Hr x Hx. +now apply exp_small_round_0_pos. +Qed. + + + + Section monotone_abs. Variable rnd : R -> Z. @@ -1283,6 +1355,33 @@ rewrite <- mantissa_DN_small_pos with (1 := Hx) (2 := He). now rewrite <- canonic_exp_fexp_pos with (1 := Hx). Qed. + +Theorem round_DN_UP_lt : + forall x, ~ generic_format x -> + (round Zfloor x < x < round Zceil x)%R. +Proof with auto with typeclass_instances. +intros x Fx. +assert (Hx:(round Zfloor x <= x <= round Zceil x)%R). +split. +apply round_DN_pt. +apply round_UP_pt. +split. + destruct Hx as [Hx _]. + apply Rnot_le_lt; intro Hle. + assert (x = round Zfloor x) by now apply Rle_antisym. + rewrite H in Fx. + contradict Fx. + apply generic_format_round... +destruct Hx as [_ Hx]. +apply Rnot_le_lt; intro Hle. +assert (x = round Zceil x) by now apply Rle_antisym. +rewrite H in Fx. +contradict Fx. +apply generic_format_round... +Qed. + + + Theorem round_UP_small_pos : forall x ex, (bpow (ex - 1) <= x < bpow ex)%R -> diff --git a/flocq/Core/Fcore_rnd.v b/flocq/Core/Fcore_rnd.v index 6b4d807..94c9420 100644 --- a/flocq/Core/Fcore_rnd.v +++ b/flocq/Core/Fcore_rnd.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/flocq/Core/Fcore_rnd_ne.v b/flocq/Core/Fcore_rnd_ne.v index 0b0776e..6829c0c 100644 --- a/flocq/Core/Fcore_rnd_ne.v +++ b/flocq/Core/Fcore_rnd_ne.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public @@ -499,6 +499,25 @@ rewrite Zeven_plus. 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). diff --git a/flocq/Core/Fcore_ulp.v b/flocq/Core/Fcore_ulp.v index 492fac6..07ef3ec 100644 --- a/flocq/Core/Fcore_ulp.v +++ b/flocq/Core/Fcore_ulp.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public @@ -1139,4 +1139,278 @@ apply le_pred_lt_aux ; try easy. now split. Qed. + +Theorem pred_succ : forall { monotone_exp : Monotone_exp fexp }, + forall x, F x -> (0 < x)%R -> pred (x + ulp x)=x. +Proof. +intros L x Fx Hx. +assert (x <= pred (x + ulp x))%R. +apply le_pred_lt. +assumption. +now apply generic_format_succ. +replace 0%R with (0+0)%R by ring; apply Rplus_lt_compat; try apply Hx. +apply bpow_gt_0. +apply Rplus_lt_reg_r with (-x)%R; ring_simplify. +apply bpow_gt_0. +apply Rle_antisym; trivial. +apply Rplus_le_reg_r with (ulp (pred (x + ulp x))). +rewrite pred_plus_ulp. +apply Rplus_le_compat_l. +now apply ulp_le. +replace 0%R with (0+0)%R by ring; apply Rplus_lt_compat; try apply Hx. +apply bpow_gt_0. +now apply generic_format_succ. +apply Rgt_not_eq. +now apply Rlt_le_trans with x. +Qed. + + +Theorem lt_UP_le_DN : + forall x y, F y -> + (y < round beta fexp Zceil x -> y <= round beta fexp Zfloor x)%R. +Proof with auto with typeclass_instances. +intros x y Fy Hlt. +apply round_DN_pt... +apply Rnot_lt_le. +contradict Hlt. +apply RIneq.Rle_not_lt. +apply round_UP_pt... +now apply Rlt_le. +Qed. + +Theorem pred_UP_le_DN : + forall x, (0 < round beta fexp Zceil x)%R -> + (pred (round beta fexp Zceil x) <= round beta fexp Zfloor x)%R. +Proof with auto with typeclass_instances. +intros x Pxu. +destruct (generic_format_EM beta fexp x) as [Fx|Fx]. +rewrite !round_generic... +now apply Rlt_le; apply pred_lt_id. +assert (let u := round beta fexp Zceil x in pred u < u)%R as Hup. + now apply pred_lt_id. +apply lt_UP_le_DN... +apply generic_format_pred... +now apply round_UP_pt. +Qed. + +Theorem pred_UP_eq_DN : + forall x, (0 < round beta fexp Zceil x)%R -> ~ F x -> + (pred (round beta fexp Zceil x) = round beta fexp Zfloor x)%R. +Proof with auto with typeclass_instances. +intros x Px Fx. +apply Rle_antisym. +now apply pred_UP_le_DN. +apply le_pred_lt; try apply generic_format_round... +pose proof round_DN_UP_lt _ _ _ Fx as HE. +now apply Rlt_trans with (1 := proj1 HE) (2 := proj2 HE). +Qed. + + + + + +(** Properties of rounding to nearest and ulp *) + +Theorem rnd_N_le_half_an_ulp: forall choice u v, + F u -> (0 < u)%R -> (v < u + (ulp u)/2)%R + -> (round beta fexp (Znearest choice) v <= u)%R. +Proof with auto with typeclass_instances. +intros choice u v Fu Hu H. +(* . *) +assert (0 < ulp u / 2)%R. +unfold Rdiv; apply Rmult_lt_0_compat. +unfold ulp; apply bpow_gt_0. +auto with real. +(* . *) +assert (ulp u / 2 < ulp u)%R. +apply Rlt_le_trans with (ulp u *1)%R;[idtac|right; ring]. +unfold Rdiv; apply Rmult_lt_compat_l. +apply bpow_gt_0. +apply Rmult_lt_reg_l with 2%R. +auto with real. +apply Rle_lt_trans with 1%R. +right; field. +rewrite Rmult_1_r; auto with real. +(* *) +apply Rnd_N_pt_monotone with F v (u + ulp u / 2)%R... +apply round_N_pt... +apply Rnd_DN_pt_N with (u+ulp u)%R. +pattern u at 3; replace u with (round beta fexp Zfloor (u + ulp u / 2)). +apply round_DN_pt... +apply round_DN_succ; try assumption. +split; try left; assumption. +replace (u+ulp u)%R with (round beta fexp Zceil (u + ulp u / 2)). +apply round_UP_pt... +apply round_UP_succ; try assumption... +split; try left; assumption. +right; field. +Qed. + + +Theorem rnd_N_ge_half_an_ulp_pred: forall choice u v, + F u -> (0 < pred u)%R -> (u - (ulp (pred u))/2 < v)%R + -> (u <= round beta fexp (Znearest choice) v)%R. +Proof with auto with typeclass_instances. +intros choice u v Fu Hu H. +(* . *) +assert (0 < u)%R. +apply Rlt_trans with (1:= Hu). +apply pred_lt_id. +assert (0 < ulp (pred u) / 2)%R. +unfold Rdiv; apply Rmult_lt_0_compat. +unfold ulp; apply bpow_gt_0. +auto with real. +assert (ulp (pred u) / 2 < ulp (pred u))%R. +apply Rlt_le_trans with (ulp (pred u) *1)%R;[idtac|right; ring]. +unfold Rdiv; apply Rmult_lt_compat_l. +apply bpow_gt_0. +apply Rmult_lt_reg_l with 2%R. +auto with real. +apply Rle_lt_trans with 1%R. +right; field. +rewrite Rmult_1_r; auto with real. +(* *) +apply Rnd_N_pt_monotone with F (u - ulp (pred u) / 2)%R v... +2: apply round_N_pt... +apply Rnd_UP_pt_N with (pred u). +pattern (pred u) at 2; replace (pred u) with (round beta fexp Zfloor (u - ulp (pred u) / 2)). +apply round_DN_pt... +replace (u - ulp (pred u) / 2)%R with (pred u + ulp (pred u) / 2)%R. +apply round_DN_succ; try assumption. +apply generic_format_pred; assumption. +split; [left|idtac]; assumption. +pattern u at 3; rewrite <- (pred_plus_ulp u); try assumption. +field. +now apply Rgt_not_eq. +pattern u at 3; replace u with (round beta fexp Zceil (u - ulp (pred u) / 2)). +apply round_UP_pt... +replace (u - ulp (pred u) / 2)%R with (pred u + ulp (pred u) / 2)%R. +apply trans_eq with (pred u +ulp(pred u))%R. +apply round_UP_succ; try assumption... +apply generic_format_pred; assumption. +split; [idtac|left]; assumption. +apply pred_plus_ulp; try assumption. +now apply Rgt_not_eq. +pattern u at 3; rewrite <- (pred_plus_ulp u); try assumption. +field. +now apply Rgt_not_eq. +pattern u at 4; rewrite <- (pred_plus_ulp u); try assumption. +right; field. +now apply Rgt_not_eq. +Qed. + + +Theorem rnd_N_ge_half_an_ulp: forall choice u v, + F u -> (0 < u)%R -> (u <> bpow (ln_beta beta u - 1))%R + -> (u - (ulp u)/2 < v)%R + -> (u <= round beta fexp (Znearest choice) v)%R. +Proof with auto with typeclass_instances. +intros choice u v Fu Hupos Hu H. +(* *) +assert (bpow (ln_beta beta u-1) <= pred u)%R. +apply le_pred_lt; try assumption. +apply generic_format_bpow. +assert (canonic_exp beta fexp u < ln_beta beta u)%Z. +apply ln_beta_generic_gt; try assumption. +now apply Rgt_not_eq. +unfold canonic_exp in H0. +ring_simplify (ln_beta beta u - 1 + 1)%Z. +omega. +destruct ln_beta as (e,He); simpl in *. +assert (bpow (e - 1) <= Rabs u)%R. +apply He. +now apply Rgt_not_eq. +rewrite Rabs_right in H0. +case H0; auto. +intros T; contradict T. +now apply sym_not_eq. +apply Rle_ge; now left. +assert (Hu2:(ulp (pred u) = ulp u)). +unfold ulp, canonic_exp. +apply f_equal; apply f_equal. +apply ln_beta_unique. +rewrite Rabs_right. +split. +assumption. +apply Rlt_trans with (1:=pred_lt_id _). +destruct ln_beta as (e,He); simpl in *. +rewrite Rabs_right in He. +apply He. +now apply Rgt_not_eq. +apply Rle_ge; now left. +apply Rle_ge, pred_ge_0; assumption. +apply rnd_N_ge_half_an_ulp_pred; try assumption. +apply Rlt_le_trans with (2:=H0). +apply bpow_gt_0. +rewrite Hu2; assumption. +Qed. + + +Lemma round_N_DN_betw: forall choice x d u, + Rnd_DN_pt (generic_format beta fexp) x d -> + Rnd_UP_pt (generic_format beta fexp) x u -> + (d<=x<(d+u)/2)%R -> + round beta fexp (Znearest choice) x = d. +Proof with auto with typeclass_instances. +intros choice x d u Hd Hu H. +apply Rnd_N_pt_unicity with (generic_format beta fexp) x d u; try assumption. +intros Y. +absurd (x < (d+u)/2)%R; try apply H. +apply Rle_not_lt; right. +apply Rplus_eq_reg_r with (-x)%R. +apply trans_eq with (- (x-d)/2 + (u-x)/2)%R. +field. +rewrite Y; field. +apply round_N_pt... +apply Rnd_DN_UP_pt_N with d u... +apply Hd. +right; apply trans_eq with (-(d-x))%R;[idtac|ring]. +apply Rabs_left1. +apply Rplus_le_reg_l with x; ring_simplify. +apply H. +rewrite Rabs_left1. +apply Rplus_le_reg_l with (d+x)%R. +apply Rmult_le_reg_l with (/2)%R. +auto with real. +apply Rle_trans with x. +right; field. +apply Rle_trans with ((d+u)/2)%R. +now left. +right; field. +apply Rplus_le_reg_l with x; ring_simplify. +apply H. +Qed. + + +Lemma round_N_UP_betw: forall choice x d u, + Rnd_DN_pt (generic_format beta fexp) x d -> + Rnd_UP_pt (generic_format beta fexp) x u -> + ((d+u)/2 < x <= u)%R -> + round beta fexp (Znearest choice) x = u. +Proof with auto with typeclass_instances. +intros choice x d u Hd Hu H. +rewrite <- (Ropp_involutive (round beta fexp (Znearest choice) x )), + <- (Ropp_involutive u) . +apply f_equal. +rewrite <- (Ropp_involutive x) . +rewrite round_N_opp, Ropp_involutive. +apply round_N_DN_betw with (-d)%R. +replace u with (round beta fexp Zceil x). +rewrite <- round_DN_opp. +apply round_DN_pt... +apply Rnd_UP_pt_unicity with (generic_format beta fexp) x... +apply round_UP_pt... +replace d with (round beta fexp Zfloor x). +rewrite <- round_UP_opp. +apply round_UP_pt... +apply Rnd_DN_pt_unicity with (generic_format beta fexp) x... +apply round_DN_pt... +split. +apply Ropp_le_contravar, H. +apply Rlt_le_trans with (-((d + u) / 2))%R. +apply Ropp_lt_contravar, H. +unfold Rdiv; right; ring. +Qed. + + End Fcore_ulp. diff --git a/flocq/Flocq_version.v b/flocq/Flocq_version.v index 662d83a..b375857 100644 --- a/flocq/Flocq_version.v +++ b/flocq/Flocq_version.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2011 Sylvie Boldo +Copyright (C) 2011-2013 Sylvie Boldo #<br /># -Copyright (C) 2011 Guillaume Melquiond +Copyright (C) 2011-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public @@ -28,4 +28,4 @@ Definition Flocq_version := Eval vm_compute in | String h t => parse t major (minor + N_of_ascii h - N_of_ascii "0"%char)%N | Empty_string => (major * 100 + minor)%N end in - parse "2.1.0"%string N0 N0. + parse "2.2.0"%string N0 N0. diff --git a/flocq/Prop/Fprop_Sterbenz.v b/flocq/Prop/Fprop_Sterbenz.v index 7d2c2e7..7260d2e 100644 --- a/flocq/Prop/Fprop_Sterbenz.v +++ b/flocq/Prop/Fprop_Sterbenz.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/flocq/Prop/Fprop_div_sqrt_error.v b/flocq/Prop/Fprop_div_sqrt_error.v index 84a0694..ec00ca4 100644 --- a/flocq/Prop/Fprop_div_sqrt_error.v +++ b/flocq/Prop/Fprop_div_sqrt_error.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/flocq/Prop/Fprop_mult_error.v b/flocq/Prop/Fprop_mult_error.v index e3e094c..e84e80b 100644 --- a/flocq/Prop/Fprop_mult_error.v +++ b/flocq/Prop/Fprop_mult_error.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/flocq/Prop/Fprop_plus_error.v b/flocq/Prop/Fprop_plus_error.v index d9dee7c..ddae698 100644 --- a/flocq/Prop/Fprop_plus_error.v +++ b/flocq/Prop/Fprop_plus_error.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public diff --git a/flocq/Prop/Fprop_relative.v b/flocq/Prop/Fprop_relative.v index 8df7336..a8cc1ff 100644 --- a/flocq/Prop/Fprop_relative.v +++ b/flocq/Prop/Fprop_relative.v @@ -2,9 +2,9 @@ This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/ -Copyright (C) 2010-2011 Sylvie Boldo +Copyright (C) 2010-2013 Sylvie Boldo #<br /># -Copyright (C) 2010-2011 Guillaume Melquiond +Copyright (C) 2010-2013 Guillaume Melquiond This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public @@ -35,9 +35,9 @@ Section relative_error_conversion. Variable rnd : R -> Z. Context { valid_rnd : Valid_rnd rnd }. -Lemma relative_error_lt_conversion : +Lemma relative_error_lt_conversion' : forall x b, (0 < b)%R -> - (Rabs (round beta fexp rnd x - x) < b * Rabs x)%R -> + (x <> 0 -> Rabs (round beta fexp rnd x - x) < b * Rabs x)%R -> exists eps, (Rabs eps < b)%R /\ round beta fexp rnd x = (x * (1 + eps))%R. Proof with auto with typeclass_instances. @@ -50,6 +50,7 @@ now rewrite Rabs_R0. rewrite Hx0, Rmult_0_l. apply round_0... (* *) +specialize (Hxb Hx0). exists ((round beta fexp rnd x - x) / x)%R. split. 2: now field. unfold Rdiv. @@ -61,6 +62,19 @@ rewrite Rinv_l with (1 := Hx0). now rewrite Rabs_R1, Rmult_1_r. Qed. +(* TODO: remove *) +Lemma relative_error_lt_conversion : + forall x b, (0 < b)%R -> + (Rabs (round beta fexp rnd x - x) < b * Rabs x)%R -> + exists eps, + (Rabs eps < b)%R /\ round beta fexp rnd x = (x * (1 + eps))%R. +Proof. +intros x b Hb0 Hxb. +apply relative_error_lt_conversion'. +exact Hb0. +now intros _. +Qed. + Lemma relative_error_le_conversion : forall x b, (0 <= b)%R -> (Rabs (round beta fexp rnd x - x) <= b * Rabs x)%R -> @@ -154,18 +168,28 @@ rewrite F2R_0, F2R_Zabs. now apply Rabs_pos_lt. Qed. -Theorem relative_error_F2R_emin_ex : +Theorem relative_error_F2R_emin_ex' : forall m, let x := F2R (Float beta m emin) in - (x <> 0)%R -> exists eps, (Rabs eps < bpow (-p + 1))%R /\ round beta fexp rnd x = (x * (1 + eps))%R. Proof with auto with typeclass_instances. -intros m x Hx. -apply relative_error_lt_conversion... +intros m x. +apply relative_error_lt_conversion'... apply bpow_gt_0. now apply relative_error_F2R_emin. Qed. +(* TODO: remove *) +Theorem relative_error_F2R_emin_ex : + forall m, let x := F2R (Float beta m emin) in + (x <> 0)%R -> + exists eps, + (Rabs eps < bpow (-p + 1))%R /\ round beta fexp rnd x = (x * (1 + eps))%R. +Proof with auto with typeclass_instances. +intros m x _. +apply relative_error_F2R_emin_ex'. +Qed. + Theorem relative_error_round : (0 < p)%Z -> forall x, @@ -404,6 +428,7 @@ Qed. Variable rnd : R -> Z. Context { valid_rnd : Valid_rnd rnd }. +(* TODO: remove *) Theorem relative_error_FLT_F2R_emin : forall m, let x := F2R (Float beta m (emin + prec - 1)) in (x <> 0)%R -> @@ -424,14 +449,49 @@ apply relative_error with (emin + prec - 1)%Z... apply relative_error_FLT_aux. Qed. +Theorem relative_error_FLT_F2R_emin' : + forall m, let x := F2R (Float beta m emin) in + (x <> 0)%R -> + (Rabs (round beta (FLT_exp emin prec) rnd x - x) < bpow (-prec + 1) * Rabs x)%R. +Proof with auto with typeclass_instances. +intros m x Zx. +destruct (Rlt_or_le (Rabs x) (bpow (emin + prec - 1))) as [Hx|Hx]. +rewrite round_generic... +unfold Rminus. +rewrite Rplus_opp_r, Rabs_R0. +apply Rmult_lt_0_compat. +apply bpow_gt_0. +now apply Rabs_pos_lt. +apply generic_format_FLT_FIX... +apply Rlt_le. +apply Rlt_le_trans with (1 := Hx). +apply bpow_le. +apply Zle_pred. +apply generic_format_FIX. +now exists (Float beta m emin). +now apply relative_error_FLT. +Qed. + +Theorem relative_error_FLT_F2R_emin_ex' : + forall m, let x := F2R (Float beta m emin) in + exists eps, + (Rabs eps < bpow (-prec + 1))%R /\ round beta (FLT_exp emin prec) rnd x = (x * (1 + eps))%R. +Proof with auto with typeclass_instances. +intros m x. +apply relative_error_lt_conversion'... +apply bpow_gt_0. +now apply relative_error_FLT_F2R_emin'. +Qed. + +(* TODO: remove *) Theorem relative_error_FLT_F2R_emin_ex : forall m, let x := F2R (Float beta m (emin + prec - 1)) in (x <> 0)%R -> exists eps, (Rabs eps < bpow (-prec + 1))%R /\ round beta (FLT_exp emin prec) rnd x = (x * (1 + eps))%R. Proof with auto with typeclass_instances. -intros m x Hx. -apply relative_error_lt_conversion... +intros m x _. +apply relative_error_lt_conversion'... apply bpow_gt_0. now apply relative_error_FLT_F2R_emin. Qed. @@ -488,6 +548,32 @@ apply relative_error_N_round with (emin + prec - 1)%Z... apply relative_error_FLT_aux. Qed. +Theorem relative_error_N_FLT_F2R_emin' : + forall m, let x := F2R (Float beta m emin) in + (Rabs (round beta (FLT_exp emin prec) (Znearest choice) x - x) <= /2 * bpow (-prec + 1) * Rabs x)%R. +Proof with auto with typeclass_instances. +intros m x. +destruct (Rlt_or_le (Rabs x) (bpow (emin + prec - 1))) as [Hx|Hx]. +rewrite round_generic... +unfold Rminus. +rewrite Rplus_opp_r, Rabs_R0. +apply Rmult_le_pos. +apply Rmult_le_pos. +apply Rlt_le. +apply (RinvN_pos 1). +apply bpow_ge_0. +apply Rabs_pos. +apply generic_format_FLT_FIX... +apply Rlt_le. +apply Rlt_le_trans with (1 := Hx). +apply bpow_le. +apply Zle_pred. +apply generic_format_FIX. +now exists (Float beta m emin). +now apply relative_error_N_FLT. +Qed. + +(* TODO: remove *) Theorem relative_error_N_FLT_F2R_emin : forall m, let x := F2R (Float beta m (emin + prec - 1)) in (Rabs (round beta (FLT_exp emin prec) (Znearest choice) x - x) <= /2 * bpow (-prec + 1) * Rabs x)%R. @@ -497,6 +583,21 @@ apply relative_error_N_F2R_emin... apply relative_error_FLT_aux. Qed. +Theorem relative_error_N_FLT_F2R_emin_ex' : + forall m, let x := F2R (Float beta m emin) in + exists eps, + (Rabs eps <= /2 * bpow (-prec + 1))%R /\ round beta (FLT_exp emin prec) (Znearest choice) x = (x * (1 + eps))%R. +Proof with auto with typeclass_instances. +intros m x. +apply relative_error_le_conversion... +apply Rmult_le_pos. +apply Rlt_le. +apply (RinvN_pos 1). +apply bpow_ge_0. +now apply relative_error_N_FLT_F2R_emin'. +Qed. + +(* TODO: remove *) Theorem relative_error_N_FLT_F2R_emin_ex : forall m, let x := F2R (Float beta m (emin + prec - 1)) in exists eps, @@ -512,6 +613,33 @@ apply bpow_gt_0. now apply relative_error_N_FLT_F2R_emin. Qed. +Theorem relative_error_N_FLT_round_F2R_emin' : + forall m, let x := F2R (Float beta m emin) in + (Rabs (round beta (FLT_exp emin prec) (Znearest choice) x - x) <= /2 * bpow (-prec + 1) * Rabs (round beta (FLT_exp emin prec) (Znearest choice) x))%R. +Proof with auto with typeclass_instances. +intros m x. +destruct (Rlt_or_le (Rabs x) (bpow (emin + prec - 1))) as [Hx|Hx]. +rewrite round_generic... +unfold Rminus. +rewrite Rplus_opp_r, Rabs_R0. +apply Rmult_le_pos. +apply Rmult_le_pos. +apply Rlt_le. +apply (RinvN_pos 1). +apply bpow_ge_0. +apply Rabs_pos. +apply generic_format_FLT_FIX... +apply Rlt_le. +apply Rlt_le_trans with (1 := Hx). +apply bpow_le. +apply Zle_pred. +apply generic_format_FIX. +now exists (Float beta m emin). +apply relative_error_N_round with (emin := (emin + prec - 1)%Z)... +apply relative_error_FLT_aux. +Qed. + +(* TODO: remove *) Theorem relative_error_N_FLT_round_F2R_emin : forall m, let x := F2R (Float beta m (emin + prec - 1)) in (Rabs (round beta (FLT_exp emin prec) (Znearest choice) x - x) <= /2 * bpow (-prec + 1) * Rabs (round beta (FLT_exp emin prec) (Znearest choice) x))%R. @@ -606,18 +734,28 @@ apply He. Qed. (** 1+#ε# property in any rounding in FLX *) -Theorem relative_error_FLX_ex : +Theorem relative_error_FLX_ex' : forall x, - (x <> 0)%R -> exists eps, (Rabs eps < bpow (-prec + 1))%R /\ round beta (FLX_exp prec) rnd x = (x * (1 + eps))%R. Proof with auto with typeclass_instances. -intros x Hx. -apply relative_error_lt_conversion... +intros x. +apply relative_error_lt_conversion'... apply bpow_gt_0. now apply relative_error_FLX. Qed. +(* TODO: remove *) +Theorem relative_error_FLX_ex : + forall x, + (x <> 0)%R -> + exists eps, + (Rabs eps < bpow (-prec + 1))%R /\ round beta (FLX_exp prec) rnd x = (x * (1 + eps))%R. +Proof with auto with typeclass_instances. +intros x _. +apply relative_error_FLX_ex'. +Qed. + Theorem relative_error_FLX_round : forall x, (x <> 0)%R -> |