From 5312915c1b29929f82e1f8de80609a277584913f Mon Sep 17 00:00:00 2001 From: xleroy Date: Thu, 28 Jun 2012 07:59:03 +0000 Subject: Use Flocq for floats git-svn-id: https://yquem.inria.fr/compcert/svn/compcert/trunk@1939 fca1b0fc-160b-0410-b1d3-a4f43f01ea2e --- flocq/Prop/Fprop_Sterbenz.v | 169 +++++++++ flocq/Prop/Fprop_div_sqrt_error.v | 299 ++++++++++++++++ flocq/Prop/Fprop_mult_error.v | 237 +++++++++++++ flocq/Prop/Fprop_plus_error.v | 234 +++++++++++++ flocq/Prop/Fprop_relative.v | 699 ++++++++++++++++++++++++++++++++++++++ 5 files changed, 1638 insertions(+) create mode 100644 flocq/Prop/Fprop_Sterbenz.v create mode 100644 flocq/Prop/Fprop_div_sqrt_error.v create mode 100644 flocq/Prop/Fprop_mult_error.v create mode 100644 flocq/Prop/Fprop_plus_error.v create mode 100644 flocq/Prop/Fprop_relative.v (limited to 'flocq/Prop') diff --git a/flocq/Prop/Fprop_Sterbenz.v b/flocq/Prop/Fprop_Sterbenz.v new file mode 100644 index 00000000..7d2c2e74 --- /dev/null +++ b/flocq/Prop/Fprop_Sterbenz.v @@ -0,0 +1,169 @@ +(** +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-2011 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. +*) + +(** * Sterbenz conditions for exact subtraction *) + +Require Import Fcore_Raux. +Require Import Fcore_defs. +Require Import Fcore_generic_fmt. +Require Import Fcalc_ops. + +Section Fprop_Sterbenz. + +Variable beta : radix. +Notation bpow e := (bpow beta e). + +Variable fexp : Z -> Z. +Context { valid_exp : Valid_exp fexp }. +Context { monotone_exp : Monotone_exp fexp }. +Notation format := (generic_format beta fexp). + +Theorem generic_format_plus : + forall x y, + format x -> format y -> + (Rabs (x + y) < bpow (Zmin (ln_beta beta x) (ln_beta beta y)))%R -> + format (x + y)%R. +Proof. +intros x y Fx Fy Hxy. +destruct (Req_dec (x + y) 0) as [Zxy|Zxy]. +rewrite Zxy. +apply generic_format_0. +destruct (Req_dec x R0) as [Zx|Zx]. +now rewrite Zx, Rplus_0_l. +destruct (Req_dec y R0) as [Zy|Zy]. +now rewrite Zy, Rplus_0_r. +revert Hxy. +destruct (ln_beta beta x) as (ex, Ex). simpl. +specialize (Ex Zx). +destruct (ln_beta beta y) as (ey, Ey). simpl. +specialize (Ey Zy). +intros Hxy. +set (fx := Float beta (Ztrunc (scaled_mantissa beta fexp x)) (fexp ex)). +assert (Hx: x = F2R fx). +rewrite Fx at 1. +unfold canonic_exp. +now rewrite ln_beta_unique with (1 := Ex). +set (fy := Float beta (Ztrunc (scaled_mantissa beta fexp y)) (fexp ey)). +assert (Hy: y = F2R fy). +rewrite Fy at 1. +unfold canonic_exp. +now rewrite ln_beta_unique with (1 := Ey). +rewrite Hx, Hy. +rewrite <- F2R_plus. +apply generic_format_F2R. +intros _. +case_eq (Fplus beta fx fy). +intros mxy exy Pxy. +rewrite <- Pxy, F2R_plus, <- Hx, <- Hy. +unfold canonic_exp. +replace exy with (fexp (Zmin ex ey)). +apply monotone_exp. +now apply ln_beta_le_bpow. +replace exy with (Fexp (Fplus beta fx fy)) by exact (f_equal Fexp Pxy). +rewrite Fexp_Fplus. +simpl. clear -monotone_exp. +apply sym_eq. +destruct (Zmin_spec ex ey) as [(H1,H2)|(H1,H2)] ; rewrite H2. +apply Zmin_l. +now apply monotone_exp. +apply Zmin_r. +apply monotone_exp. +apply Zlt_le_weak. +now apply Zgt_lt. +Qed. + +Theorem generic_format_plus_weak : + forall x y, + format x -> format y -> + (Rabs (x + y) <= Rmin (Rabs x) (Rabs y))%R -> + format (x + y)%R. +Proof. +intros x y Fx Fy Hxy. +destruct (Req_dec x R0) as [Zx|Zx]. +now rewrite Zx, Rplus_0_l. +destruct (Req_dec y R0) as [Zy|Zy]. +now rewrite Zy, Rplus_0_r. +apply generic_format_plus ; try assumption. +apply Rle_lt_trans with (1 := Hxy). +unfold Rmin. +destruct (Rle_dec (Rabs x) (Rabs y)) as [Hxy'|Hxy']. +rewrite Zmin_l. +destruct (ln_beta beta x) as (ex, Hx). +now apply Hx. +now apply ln_beta_le_abs. +rewrite Zmin_r. +destruct (ln_beta beta y) as (ex, Hy). +now apply Hy. +apply ln_beta_le_abs. +exact Zy. +apply Rlt_le. +now apply Rnot_le_lt. +Qed. + +Lemma sterbenz_aux : + forall x y, format x -> format y -> + (y <= x <= 2 * y)%R -> + format (x - y)%R. +Proof. +intros x y Hx Hy (Hxy1, Hxy2). +unfold Rminus. +apply generic_format_plus_weak. +exact Hx. +now apply generic_format_opp. +rewrite Rabs_pos_eq. +rewrite Rabs_Ropp. +rewrite Rmin_comm. +assert (Hy0: (0 <= y)%R). +apply Rplus_le_reg_r with y. +apply Rle_trans with x. +now rewrite Rplus_0_l. +now rewrite Rmult_plus_distr_r, Rmult_1_l in Hxy2. +rewrite Rabs_pos_eq with (1 := Hy0). +rewrite Rabs_pos_eq. +unfold Rmin. +destruct (Rle_dec y x) as [Hyx|Hyx]. +apply Rplus_le_reg_r with y. +now ring_simplify. +now elim Hyx. +now apply Rle_trans with y. +now apply Rle_0_minus. +Qed. + +Theorem sterbenz : + forall x y, format x -> format y -> + (y / 2 <= x <= 2 * y)%R -> + format (x - y)%R. +Proof. +intros x y Hx Hy (Hxy1, Hxy2). +destruct (Rle_or_lt x y) as [Hxy|Hxy]. +rewrite <- Ropp_minus_distr. +apply generic_format_opp. +apply sterbenz_aux ; try easy. +split. +exact Hxy. +apply Rcompare_not_Lt_inv. +rewrite <- Rcompare_half_r. +now apply Rcompare_not_Lt. +apply sterbenz_aux ; try easy. +split. +now apply Rlt_le. +exact Hxy2. +Qed. + +End Fprop_Sterbenz. diff --git a/flocq/Prop/Fprop_div_sqrt_error.v b/flocq/Prop/Fprop_div_sqrt_error.v new file mode 100644 index 00000000..84a06942 --- /dev/null +++ b/flocq/Prop/Fprop_div_sqrt_error.v @@ -0,0 +1,299 @@ +(** +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-2011 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. +*) + +(** * Remainder of the division and square root are in the FLX format *) +Require Import Fcore. +Require Import Fcalc_ops. +Require Import Fprop_relative. + +Section Fprop_divsqrt_error. + +Variable beta : radix. +Notation bpow e := (bpow beta e). + +Variable prec : Z. + +Theorem generic_format_plus_prec: + forall fexp, (forall e, (fexp e <= e - prec)%Z) -> + forall x y (fx fy: float beta), + (x = F2R fx)%R -> (y = F2R fy)%R -> (Rabs (x+y) < bpow (prec+Fexp fx))%R -> (Rabs (x+y) < bpow (prec+Fexp fy))%R + -> generic_format beta fexp (x+y)%R. +intros fexp Hfexp x y fx fy Hx Hy H1 H2. +case (Req_dec (x+y) 0); intros H. +rewrite H; apply generic_format_0. +rewrite Hx, Hy, <- F2R_plus. +apply generic_format_F2R. +intros _. +case_eq (Fplus beta fx fy). +intros mz ez Hz. +rewrite <- Hz. +apply Zle_trans with (Zmin (Fexp fx) (Fexp fy)). +rewrite F2R_plus, <- Hx, <- Hy. +unfold canonic_exp. +apply Zle_trans with (1:=Hfexp _). +apply Zplus_le_reg_l with prec; ring_simplify. +apply ln_beta_le_bpow with (1 := H). +now apply Zmin_case. +rewrite <- Fexp_Fplus, Hz. +apply Zle_refl. +Qed. + +Theorem ex_Fexp_canonic: forall fexp, forall x, generic_format beta fexp x + -> exists fx:float beta, (x=F2R fx)%R /\ Fexp fx = canonic_exp beta fexp x. +intros fexp x; unfold generic_format. +exists (Float beta (Ztrunc (scaled_mantissa beta fexp x)) (canonic_exp beta fexp x)). +split; auto. +Qed. + + +Context { prec_gt_0_ : Prec_gt_0 prec }. + +Notation format := (generic_format beta (FLX_exp prec)). +Notation cexp := (canonic_exp beta (FLX_exp prec)). + +Variable choice : Z -> bool. + + +(** Remainder of the division in FLX *) +Theorem div_error_FLX : + forall rnd { Zrnd : Valid_rnd rnd } x y, + format x -> format y -> + format (x - round beta (FLX_exp prec) rnd (x/y) * y)%R. +Proof with auto with typeclass_instances. +intros rnd Zrnd x y Hx Hy. +destruct (Req_dec y 0) as [Zy|Zy]. +now rewrite Zy, Rmult_0_r, Rminus_0_r. +destruct (Req_dec (round beta (FLX_exp prec) rnd (x/y)) 0) as [Hr|Hr]. +rewrite Hr; ring_simplify (x-0*y)%R; assumption. +assert (Zx: x <> R0). +contradict Hr. +rewrite Hr. +unfold Rdiv. +now rewrite Rmult_0_l, round_0. +destruct (ex_Fexp_canonic _ x Hx) as (fx,(Hx1,Hx2)). +destruct (ex_Fexp_canonic _ y Hy) as (fy,(Hy1,Hy2)). +destruct (ex_Fexp_canonic (FLX_exp prec) (round beta (FLX_exp prec) rnd (x / y))) as (fr,(Hr1,Hr2)). +apply generic_format_round... +unfold Rminus; apply generic_format_plus_prec with fx (Fopp beta (Fmult beta fr fy)); trivial. +intros e; apply Zle_refl. +now rewrite F2R_opp, F2R_mult, <- Hr1, <- Hy1. +(* *) +destruct (relative_error_FLX_ex beta prec (prec_gt_0 prec) rnd (x / y)%R) as (eps,(Heps1,Heps2)). +apply Rmult_integral_contrapositive_currified. +exact Zx. +now apply Rinv_neq_0_compat. +rewrite Heps2. +rewrite <- Rabs_Ropp. +replace (-(x + - (x / y * (1 + eps) * y)))%R with (x * eps)%R by now field. +rewrite Rabs_mult. +apply Rlt_le_trans with (Rabs x * 1)%R. +apply Rmult_lt_compat_l. +now apply Rabs_pos_lt. +apply Rlt_le_trans with (1 := Heps1). +change R1 with (bpow 0). +apply bpow_le. +generalize (prec_gt_0 prec). +clear ; omega. +rewrite Rmult_1_r. +rewrite Hx2. +unfold canonic_exp. +destruct (ln_beta beta x) as (ex, Hex). +simpl. +specialize (Hex Zx). +apply Rlt_le. +apply Rlt_le_trans with (1 := proj2 Hex). +apply bpow_le. +unfold FLX_exp. +ring_simplify. +apply Zle_refl. +(* *) +replace (Fexp (Fopp beta (Fmult beta fr fy))) with (Fexp fr + Fexp fy)%Z. +2: unfold Fopp, Fmult; destruct fr; destruct fy; now simpl. +replace (x + - (round beta (FLX_exp prec) rnd (x / y) * y))%R with + (y * (-(round beta (FLX_exp prec) rnd (x / y) - x/y)))%R. +2: field; assumption. +rewrite Rabs_mult. +apply Rlt_le_trans with (Rabs y * bpow (Fexp fr))%R. +apply Rmult_lt_compat_l. +now apply Rabs_pos_lt. +rewrite Rabs_Ropp. +replace (bpow (Fexp fr)) with (ulp beta (FLX_exp prec) (F2R fr)). +rewrite <- Hr1. +apply ulp_error_f... +unfold ulp; apply f_equal. +now rewrite Hr2, <- Hr1. +replace (prec+(Fexp fr+Fexp fy))%Z with ((prec+Fexp fy)+Fexp fr)%Z by ring. +rewrite bpow_plus. +apply Rmult_le_compat_r. +apply bpow_ge_0. +rewrite Hy2; unfold canonic_exp, FLX_exp. +ring_simplify (prec + (ln_beta beta y - prec))%Z. +destruct (ln_beta beta y); simpl. +left; now apply a. +Qed. + +(** Remainder of the square in FLX (with p>1) and rounding to nearest *) +Variable Hp1 : Zlt 1 prec. + +Theorem sqrt_error_FLX_N : + forall x, format x -> + format (x - Rsqr (round beta (FLX_exp prec) (Znearest choice) (sqrt x)))%R. +Proof with auto with typeclass_instances. +intros x Hx. +destruct (total_order_T x 0) as [[Hxz|Hxz]|Hxz]. +unfold sqrt. +destruct (Rcase_abs x). +rewrite round_0... +unfold Rsqr. +now rewrite Rmult_0_l, Rminus_0_r. +elim (Rlt_irrefl 0). +now apply Rgt_ge_trans with x. +rewrite Hxz, sqrt_0, round_0... +unfold Rsqr. +rewrite Rmult_0_l, Rminus_0_r. +apply generic_format_0. +case (Req_dec (round beta (FLX_exp prec) (Znearest choice) (sqrt x)) 0); intros Hr. +rewrite Hr; unfold Rsqr; ring_simplify (x-0*0)%R; assumption. +destruct (ex_Fexp_canonic _ x Hx) as (fx,(Hx1,Hx2)). +destruct (ex_Fexp_canonic (FLX_exp prec) (round beta (FLX_exp prec) (Znearest choice) (sqrt x))) as (fr,(Hr1,Hr2)). +apply generic_format_round... +unfold Rminus; apply generic_format_plus_prec with fx (Fopp beta (Fmult beta fr fr)); trivial. +intros e; apply Zle_refl. +unfold Rsqr; now rewrite F2R_opp,F2R_mult, <- Hr1. +(* *) +apply Rle_lt_trans with x. +apply Rabs_minus_le. +apply Rle_0_sqr. +destruct (relative_error_N_FLX_ex beta prec (prec_gt_0 prec) choice (sqrt x)) as (eps,(Heps1,Heps2)). +rewrite Heps2. +rewrite Rsqr_mult, Rsqr_sqrt, Rmult_comm. 2: now apply Rlt_le. +apply Rmult_le_compat_r. +now apply Rlt_le. +apply Rle_trans with (5²/4²)%R. +rewrite <- Rsqr_div. +apply Rsqr_le_abs_1. +apply Rle_trans with (1 := Rabs_triang _ _). +rewrite Rabs_R1. +apply Rplus_le_reg_l with (-1)%R. +rewrite <- Rplus_assoc, Rplus_opp_l, Rplus_0_l. +apply Rle_trans with (1 := Heps1). +rewrite Rabs_pos_eq. +apply Rmult_le_reg_l with 2%R. +now apply (Z2R_lt 0 2). +rewrite <- Rmult_assoc, Rinv_r, Rmult_1_l. +apply Rle_trans with (bpow (-1)). +apply bpow_le. +omega. +replace (2 * (-1 + 5 / 4))%R with (/2)%R by field. +apply Rinv_le. +now apply (Z2R_lt 0 2). +apply (Z2R_le 2). +unfold Zpower_pos. simpl. +rewrite Zmult_1_r. +apply Zle_bool_imp_le. +apply beta. +apply Rgt_not_eq. +now apply (Z2R_lt 0 2). +unfold Rdiv. +apply Rmult_le_pos. +now apply (Z2R_le 0 5). +apply Rlt_le. +apply Rinv_0_lt_compat. +now apply (Z2R_lt 0 4). +apply Rgt_not_eq. +now apply (Z2R_lt 0 4). +unfold Rsqr. +replace (5 * 5 / (4 * 4))%R with (25 * /16)%R by field. +apply Rmult_le_reg_r with 16%R. +now apply (Z2R_lt 0 16). +rewrite Rmult_assoc, Rinv_l, Rmult_1_r. +now apply (Z2R_le 25 32). +apply Rgt_not_eq. +now apply (Z2R_lt 0 16). +rewrite Hx2; unfold canonic_exp, FLX_exp. +ring_simplify (prec + (ln_beta beta x - prec))%Z. +destruct (ln_beta beta x); simpl. +rewrite <- (Rabs_right x). +apply a. +now apply Rgt_not_eq. +now apply Rgt_ge. +(* *) +replace (Fexp (Fopp beta (Fmult beta fr fr))) with (Fexp fr + Fexp fr)%Z. +2: unfold Fopp, Fmult; destruct fr; now simpl. +rewrite Hr1. +replace (x + - Rsqr (F2R fr))%R with (-((F2R fr - sqrt x)*(F2R fr + sqrt x)))%R. +2: rewrite <- (sqrt_sqrt x) at 3; auto. +2: unfold Rsqr; ring. +rewrite Rabs_Ropp, Rabs_mult. +apply Rle_lt_trans with ((/2*bpow (Fexp fr))* Rabs (F2R fr + sqrt x))%R. +apply Rmult_le_compat_r. +apply Rabs_pos. +apply Rle_trans with (/2*ulp beta (FLX_exp prec) (F2R fr))%R. +rewrite <- Hr1. +apply ulp_half_error_f... +right; unfold ulp; apply f_equal. +rewrite Hr2, <- Hr1; trivial. +rewrite Rmult_assoc, Rmult_comm. +replace (prec+(Fexp fr+Fexp fr))%Z with (Fexp fr + (prec+Fexp fr))%Z by ring. +rewrite bpow_plus, Rmult_assoc. +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 (Rabs (F2R fr + sqrt x)). +right; field. +apply Rle_lt_trans with (1:=Rabs_triang _ _). +(* . *) +assert (Rabs (F2R fr) < bpow (prec + Fexp fr))%R. +rewrite Hr2; unfold canonic_exp; rewrite Hr1. +unfold FLX_exp. +ring_simplify (prec + (ln_beta beta (F2R fr) - prec))%Z. +destruct (ln_beta beta (F2R fr)); simpl. +apply a. +rewrite <- Hr1; auto. +(* . *) +apply Rlt_le_trans with (bpow (prec + Fexp fr)+ Rabs (sqrt x))%R. +now apply Rplus_lt_compat_r. +(* . *) +rewrite Rmult_plus_distr_r, Rmult_1_l. +apply Rplus_le_compat_l. +assert (sqrt x <> 0)%R. +apply Rgt_not_eq. +now apply sqrt_lt_R0. +destruct (ln_beta beta (sqrt x)) as (es,Es). +specialize (Es H0). +apply Rle_trans with (bpow es). +now apply Rlt_le. +apply bpow_le. +case (Zle_or_lt es (prec + Fexp fr)) ; trivial. +intros H1. +absurd (Rabs (F2R fr) < bpow (es - 1))%R. +apply Rle_not_lt. +rewrite <- Hr1. +apply abs_round_ge_generic... +apply generic_format_bpow. +unfold FLX_exp; omega. +apply Es. +apply Rlt_le_trans with (1:=H). +apply bpow_le. +omega. +now apply Rlt_le. +Qed. + +End Fprop_divsqrt_error. diff --git a/flocq/Prop/Fprop_mult_error.v b/flocq/Prop/Fprop_mult_error.v new file mode 100644 index 00000000..e3e094c6 --- /dev/null +++ b/flocq/Prop/Fprop_mult_error.v @@ -0,0 +1,237 @@ +(** +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-2011 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. +*) + +(** * Error of the multiplication is in the FLX/FLT format *) +Require Import Fcore. +Require Import Fcalc_ops. + +Section Fprop_mult_error. + +Variable beta : radix. +Notation bpow e := (bpow beta e). + +Variable prec : Z. +Context { prec_gt_0_ : Prec_gt_0 prec }. + +Notation format := (generic_format beta (FLX_exp prec)). +Notation cexp := (canonic_exp beta (FLX_exp prec)). + +Variable rnd : R -> Z. +Context { valid_rnd : Valid_rnd rnd }. + +(** Auxiliary result that provides the exponent *) +Lemma mult_error_FLX_aux: + forall x y, + format x -> format y -> + (round beta (FLX_exp prec) rnd (x * y) - (x * y) <> 0)%R -> + exists f:float beta, + (F2R f = round beta (FLX_exp prec) rnd (x * y) - (x * y))%R + /\ (canonic_exp beta (FLX_exp prec) (F2R f) <= Fexp f)%Z + /\ (Fexp f = cexp x + cexp y)%Z. +Proof with auto with typeclass_instances. +intros x y Hx Hy Hz. +set (f := (round beta (FLX_exp prec) rnd (x * y))). +destruct (Req_dec (x * y) 0) as [Hxy0|Hxy0]. +contradict Hz. +rewrite Hxy0. +rewrite round_0... +ring. +destruct (ln_beta beta (x * y)) as (exy, Hexy). +specialize (Hexy Hxy0). +destruct (ln_beta beta (f - x * y)) as (er, Her). +specialize (Her Hz). +destruct (ln_beta beta x) as (ex, Hex). +assert (Hx0: (x <> 0)%R). +contradict Hxy0. +now rewrite Hxy0, Rmult_0_l. +specialize (Hex Hx0). +destruct (ln_beta beta y) as (ey, Hey). +assert (Hy0: (y <> 0)%R). +contradict Hxy0. +now rewrite Hxy0, Rmult_0_r. +specialize (Hey Hy0). +(* *) +assert (Hc1: (cexp (x * y)%R - prec <= cexp x + cexp y)%Z). +unfold canonic_exp, FLX_exp. +rewrite ln_beta_unique with (1 := Hex). +rewrite ln_beta_unique with (1 := Hey). +rewrite ln_beta_unique with (1 := Hexy). +cut (exy - 1 < ex + ey)%Z. omega. +apply (lt_bpow beta). +apply Rle_lt_trans with (1 := proj1 Hexy). +rewrite Rabs_mult. +rewrite bpow_plus. +apply Rmult_le_0_lt_compat. +apply Rabs_pos. +apply Rabs_pos. +apply Hex. +apply Hey. +(* *) +assert (Hc2: (cexp x + cexp y <= cexp (x * y)%R)%Z). +unfold canonic_exp, FLX_exp. +rewrite ln_beta_unique with (1 := Hex). +rewrite ln_beta_unique with (1 := Hey). +rewrite ln_beta_unique with (1 := Hexy). +cut ((ex - 1) + (ey - 1) < exy)%Z. +generalize (prec_gt_0 prec). +clear ; omega. +apply (lt_bpow beta). +apply Rle_lt_trans with (2 := proj2 Hexy). +rewrite Rabs_mult. +rewrite bpow_plus. +apply Rmult_le_compat. +apply bpow_ge_0. +apply bpow_ge_0. +apply Hex. +apply Hey. +(* *) +assert (Hr: ((F2R (Float beta (- (Ztrunc (scaled_mantissa beta (FLX_exp prec) x) * + Ztrunc (scaled_mantissa beta (FLX_exp prec) y)) + rnd (scaled_mantissa beta (FLX_exp prec) (x * y)) * + beta ^ (cexp (x * y)%R - (cexp x + cexp y))) (cexp x + cexp y))) = f - x * y)%R). +rewrite Hx at 6. +rewrite Hy at 6. +rewrite <- F2R_mult. +simpl. +unfold f, round, Rminus. +rewrite <- F2R_opp, Rplus_comm, <- F2R_plus. +unfold Fplus. simpl. +now rewrite Zle_imp_le_bool with (1 := Hc2). +(* *) +exists (Float beta (- (Ztrunc (scaled_mantissa beta (FLX_exp prec) x) * + Ztrunc (scaled_mantissa beta (FLX_exp prec) y)) + rnd (scaled_mantissa beta (FLX_exp prec) (x * y)) * + beta ^ (cexp (x * y)%R - (cexp x + cexp y))) (cexp x + cexp y)). +split;[assumption|split]. +rewrite Hr. +simpl. +clear Hr. +apply Zle_trans with (cexp (x * y)%R - prec)%Z. +unfold canonic_exp, FLX_exp. +apply Zplus_le_compat_r. +rewrite ln_beta_unique with (1 := Hexy). +apply ln_beta_le_bpow with (1 := Hz). +replace (bpow (exy - prec)) with (ulp beta (FLX_exp prec) (x * y)). +apply ulp_error... +unfold ulp, canonic_exp. +now rewrite ln_beta_unique with (1 := Hexy). +apply Hc1. +reflexivity. +Qed. + +(** Error of the multiplication in FLX *) +Theorem mult_error_FLX : + forall x y, + format x -> format y -> + format (round beta (FLX_exp prec) rnd (x * y) - (x * y))%R. +Proof. +intros x y Hx Hy. +destruct (Req_dec (round beta (FLX_exp prec) rnd (x * y) - x * y) 0) as [Hr0|Hr0]. +rewrite Hr0. +apply generic_format_0. +destruct (mult_error_FLX_aux x y Hx Hy Hr0) as ((m,e),(H1,(H2,H3))). +rewrite <- H1. +now apply generic_format_F2R. +Qed. + +End Fprop_mult_error. + +Section Fprop_mult_error_FLT. + +Variable beta : radix. +Notation bpow e := (bpow beta e). + +Variable emin prec : Z. +Context { prec_gt_0_ : Prec_gt_0 prec }. +Variable Hpemin: (emin <= prec)%Z. + +Notation format := (generic_format beta (FLT_exp emin prec)). +Notation cexp := (canonic_exp beta (FLT_exp emin prec)). + +Variable rnd : R -> Z. +Context { valid_rnd : Valid_rnd rnd }. + +(** Error of the multiplication in FLT with underflow requirements *) +Theorem mult_error_FLT : + forall x y, + format x -> format y -> + (x*y = 0)%R \/ (bpow (emin + 2*prec - 1) <= Rabs (x * y))%R -> + format (round beta (FLT_exp emin prec) rnd (x * y) - (x * y))%R. +Proof with auto with typeclass_instances. +clear Hpemin. +intros x y Hx Hy Hxy. +set (f := (round beta (FLT_exp emin prec) rnd (x * y))). +destruct (Req_dec (f - x * y) 0) as [Hr0|Hr0]. +rewrite Hr0. +apply generic_format_0. +destruct Hxy as [Hxy|Hxy]. +unfold f. +rewrite Hxy. +rewrite round_0... +ring_simplify (0 - 0)%R. +apply generic_format_0. +destruct (mult_error_FLX_aux beta prec rnd x y) as ((m,e),(H1,(H2,H3))). +now apply generic_format_FLX_FLT with emin. +now apply generic_format_FLX_FLT with emin. +rewrite <- (round_FLT_FLX beta emin). +assumption. +apply Rle_trans with (2:=Hxy). +apply bpow_le. +generalize (prec_gt_0 prec). +clear ; omega. +rewrite <- (round_FLT_FLX beta emin) in H1. +2:apply Rle_trans with (2:=Hxy). +2:apply bpow_le ; generalize (prec_gt_0 prec) ; clear ; omega. +unfold f; rewrite <- H1. +apply generic_format_F2R. +intros _. +simpl in H2, H3. +unfold canonic_exp, FLT_exp. +case (Zmax_spec (ln_beta beta (F2R (Float beta m e)) - prec) emin); + intros (M1,M2); rewrite M2. +apply Zle_trans with (2:=H2). +unfold canonic_exp, FLX_exp. +apply Zle_refl. +rewrite H3. +unfold canonic_exp, FLX_exp. +assert (Hxy0:(x*y <> 0)%R). +contradict Hr0. +unfold f. +rewrite Hr0. +rewrite round_0... +ring. +assert (Hx0: (x <> 0)%R). +contradict Hxy0. +now rewrite Hxy0, Rmult_0_l. +assert (Hy0: (y <> 0)%R). +contradict Hxy0. +now rewrite Hxy0, Rmult_0_r. +destruct (ln_beta beta x) as (ex,Ex) ; simpl. +specialize (Ex Hx0). +destruct (ln_beta beta y) as (ey,Ey) ; simpl. +specialize (Ey Hy0). +assert (emin + 2 * prec -1 < ex + ey)%Z. +2: omega. +apply (lt_bpow beta). +apply Rle_lt_trans with (1:=Hxy). +rewrite Rabs_mult, bpow_plus. +apply Rmult_le_0_lt_compat; try apply Rabs_pos. +apply Ex. +apply Ey. +Qed. + +End Fprop_mult_error_FLT. diff --git a/flocq/Prop/Fprop_plus_error.v b/flocq/Prop/Fprop_plus_error.v new file mode 100644 index 00000000..d9dee7ce --- /dev/null +++ b/flocq/Prop/Fprop_plus_error.v @@ -0,0 +1,234 @@ +(** +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-2011 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. +*) + +(** * Error of the rounded-to-nearest addition is representable. *) + +Require Import Fcore_Raux. +Require Import Fcore_defs. +Require Import Fcore_float_prop. +Require Import Fcore_generic_fmt. +Require Import Fcalc_ops. + +Section Fprop_plus_error. + +Variable beta : radix. +Notation bpow e := (bpow beta e). + +Variable fexp : Z -> Z. +Context { valid_exp : Valid_exp fexp }. + +Section round_repr_same_exp. + +Variable rnd : R -> Z. +Context { valid_rnd : Valid_rnd rnd }. + +Theorem round_repr_same_exp : + forall m e, + exists m', + round beta fexp rnd (F2R (Float beta m e)) = F2R (Float beta m' e). +Proof with auto with typeclass_instances. +intros m e. +set (e' := canonic_exp beta fexp (F2R (Float beta m e))). +unfold round, scaled_mantissa. fold e'. +destruct (Zle_or_lt e' e) as [He|He]. +exists m. +unfold F2R at 2. simpl. +rewrite Rmult_assoc, <- bpow_plus. +rewrite <- Z2R_Zpower. 2: omega. +rewrite <- Z2R_mult, Zrnd_Z2R... +unfold F2R. simpl. +rewrite Z2R_mult. +rewrite Rmult_assoc. +rewrite Z2R_Zpower. 2: omega. +rewrite <- bpow_plus. +apply (f_equal (fun v => Z2R m * bpow v)%R). +ring. +exists ((rnd (Z2R m * bpow (e - e'))) * Zpower beta (e' - e))%Z. +unfold F2R. simpl. +rewrite Z2R_mult. +rewrite Z2R_Zpower. 2: omega. +rewrite 2!Rmult_assoc. +rewrite <- 2!bpow_plus. +apply (f_equal (fun v => _ * bpow v)%R). +ring. +Qed. + +End round_repr_same_exp. + +Context { monotone_exp : Monotone_exp fexp }. +Notation format := (generic_format beta fexp). + +Variable choice : Z -> bool. + +Lemma plus_error_aux : + forall x y, + (canonic_exp beta fexp x <= canonic_exp beta fexp y)%Z -> + format x -> format y -> + format (round beta fexp (Znearest choice) (x + y) - (x + y))%R. +Proof. +intros x y. +set (ex := canonic_exp beta fexp x). +set (ey := canonic_exp beta fexp y). +intros He Hx Hy. +destruct (Req_dec (round beta fexp (Znearest choice) (x + y) - (x + y)) R0) as [H0|H0]. +rewrite H0. +apply generic_format_0. +set (mx := Ztrunc (scaled_mantissa beta fexp x)). +set (my := Ztrunc (scaled_mantissa beta fexp y)). +(* *) +assert (Hxy: (x + y)%R = F2R (Float beta (mx + my * beta ^ (ey - ex)) ex)). +rewrite Hx, Hy. +fold mx my ex ey. +rewrite <- F2R_plus. +unfold Fplus. simpl. +now rewrite Zle_imp_le_bool with (1 := He). +(* *) +rewrite Hxy. +destruct (round_repr_same_exp (Znearest choice) (mx + my * beta ^ (ey - ex)) ex) as (mxy, Hxy'). +rewrite Hxy'. +assert (H: (F2R (Float beta mxy ex) - F2R (Float beta (mx + my * beta ^ (ey - ex)) ex))%R = + F2R (Float beta (mxy - (mx + my * beta ^ (ey - ex))) ex)). +now rewrite <- F2R_minus, Fminus_same_exp. +rewrite H. +apply generic_format_F2R. +intros _. +apply monotone_exp. +rewrite <- H, <- Hxy', <- Hxy. +apply ln_beta_le_abs. +exact H0. +pattern x at 3 ; replace x with (-(y - (x + y)))%R by ring. +rewrite Rabs_Ropp. +now apply (round_N_pt beta _ choice (x + y)). +Qed. + +(** Error of the addition *) +Theorem plus_error : + forall x y, + format x -> format y -> + format (round beta fexp (Znearest choice) (x + y) - (x + y))%R. +Proof. +intros x y Hx Hy. +destruct (Zle_or_lt (canonic_exp beta fexp x) (canonic_exp beta fexp y)). +now apply plus_error_aux. +rewrite Rplus_comm. +apply plus_error_aux ; try easy. +now apply Zlt_le_weak. +Qed. + +End Fprop_plus_error. + +Section Fprop_plus_zero. + +Variable beta : radix. +Notation bpow e := (bpow beta e). + +Variable fexp : Z -> Z. +Context { valid_exp : Valid_exp fexp }. +Context { exp_not_FTZ : Exp_not_FTZ fexp }. +Notation format := (generic_format beta fexp). + +Section round_plus_eq_zero_aux. + +Variable rnd : R -> Z. +Context { valid_rnd : Valid_rnd rnd }. + +Lemma round_plus_eq_zero_aux : + forall x y, + (canonic_exp beta fexp x <= canonic_exp beta fexp y)%Z -> + format x -> format y -> + (0 <= x + y)%R -> + round beta fexp rnd (x + y) = R0 -> + (x + y = 0)%R. +Proof with auto with typeclass_instances. +intros x y He Hx Hy Hp Hxy. +destruct (Req_dec (x + y) 0) as [H0|H0]. +exact H0. +destruct (ln_beta beta (x + y)) as (exy, Hexy). +simpl. +specialize (Hexy H0). +destruct (Zle_or_lt exy (fexp exy)) as [He'|He']. +(* . *) +assert (H: (x + y)%R = F2R (Float beta (Ztrunc (x * bpow (- fexp exy)) + + Ztrunc (y * bpow (- fexp exy))) (fexp exy))). +rewrite (subnormal_exponent beta fexp exy x He' Hx) at 1. +rewrite (subnormal_exponent beta fexp exy y He' Hy) at 1. +now rewrite <- F2R_plus, Fplus_same_exp. +rewrite H in Hxy. +rewrite round_generic in Hxy... +now rewrite <- H in Hxy. +apply generic_format_F2R. +intros _. +rewrite <- H. +unfold canonic_exp. +rewrite ln_beta_unique with (1 := Hexy). +apply Zle_refl. +(* . *) +elim Rle_not_lt with (1 := round_le beta _ rnd _ _ (proj1 Hexy)). +rewrite (Rabs_pos_eq _ Hp). +rewrite Hxy. +rewrite round_generic... +apply bpow_gt_0. +apply generic_format_bpow. +apply Zlt_succ_le. +now rewrite (Zsucc_pred exy) in He'. +Qed. + +End round_plus_eq_zero_aux. + +Variable rnd : R -> Z. +Context { valid_rnd : Valid_rnd rnd }. + +(** rnd(x+y)=0 -> x+y = 0 provided this is not a FTZ format *) +Theorem round_plus_eq_zero : + forall x y, + format x -> format y -> + round beta fexp rnd (x + y) = R0 -> + (x + y = 0)%R. +Proof with auto with typeclass_instances. +intros x y Hx Hy. +destruct (Rle_or_lt R0 (x + y)) as [H1|H1]. +(* . *) +revert H1. +destruct (Zle_or_lt (canonic_exp beta fexp x) (canonic_exp beta fexp y)) as [H2|H2]. +now apply round_plus_eq_zero_aux. +rewrite Rplus_comm. +apply round_plus_eq_zero_aux ; try easy. +now apply Zlt_le_weak. +(* . *) +revert H1. +rewrite <- (Ropp_involutive (x + y)), Ropp_plus_distr, <- Ropp_0. +intros H1. +rewrite round_opp. +intros Hxy. +apply f_equal. +cut (round beta fexp (Zrnd_opp rnd) (- x + - y) = 0)%R. +cut (0 <= -x + -y)%R. +destruct (Zle_or_lt (canonic_exp beta fexp (-x)) (canonic_exp beta fexp (-y))) as [H2|H2]. +apply round_plus_eq_zero_aux ; try apply generic_format_opp... +rewrite Rplus_comm. +apply round_plus_eq_zero_aux ; try apply generic_format_opp... +now apply Zlt_le_weak. +apply Rlt_le. +now apply Ropp_lt_cancel. +rewrite <- (Ropp_involutive (round _ _ _ _)). +rewrite Hxy. +apply Ropp_involutive. +Qed. + +End Fprop_plus_zero. diff --git a/flocq/Prop/Fprop_relative.v b/flocq/Prop/Fprop_relative.v new file mode 100644 index 00000000..8df73361 --- /dev/null +++ b/flocq/Prop/Fprop_relative.v @@ -0,0 +1,699 @@ +(** +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-2011 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. +*) + +(** * Relative error of the roundings *) +Require Import Fcore. + +Section Fprop_relative. + +Variable beta : radix. +Notation bpow e := (bpow beta e). + +Section Fprop_relative_generic. + +Variable fexp : Z -> Z. +Context { prop_exp : Valid_exp fexp }. + +Section relative_error_conversion. + +Variable rnd : R -> Z. +Context { valid_rnd : Valid_rnd rnd }. + +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 with auto with typeclass_instances. +intros x b Hb0 Hxb. +destruct (Req_dec x 0) as [Hx0|Hx0]. +(* *) +exists R0. +split. +now rewrite Rabs_R0. +rewrite Hx0, Rmult_0_l. +apply round_0... +(* *) +exists ((round beta fexp rnd x - x) / x)%R. +split. 2: now field. +unfold Rdiv. +rewrite Rabs_mult. +apply Rmult_lt_reg_r with (Rabs x). +now apply Rabs_pos_lt. +rewrite Rmult_assoc, <- Rabs_mult. +rewrite Rinv_l with (1 := Hx0). +now rewrite Rabs_R1, Rmult_1_r. +Qed. + +Lemma relative_error_le_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 with auto with typeclass_instances. +intros x b Hb0 Hxb. +destruct (Req_dec x 0) as [Hx0|Hx0]. +(* *) +exists R0. +split. +now rewrite Rabs_R0. +rewrite Hx0, Rmult_0_l. +apply round_0... +(* *) +exists ((round beta fexp rnd x - x) / x)%R. +split. 2: now field. +unfold Rdiv. +rewrite Rabs_mult. +apply Rmult_le_reg_r with (Rabs x). +now apply Rabs_pos_lt. +rewrite Rmult_assoc, <- Rabs_mult. +rewrite Rinv_l with (1 := Hx0). +now rewrite Rabs_R1, Rmult_1_r. +Qed. + +End relative_error_conversion. + +Variable emin p : Z. +Hypothesis Hmin : forall k, (emin < k)%Z -> (p <= k - fexp k)%Z. + +Variable rnd : R -> Z. +Context { valid_rnd : Valid_rnd rnd }. + +Theorem relative_error : + forall x, + (bpow emin <= Rabs x)%R -> + (Rabs (round beta fexp rnd x - x) < bpow (-p + 1) * Rabs x)%R. +Proof. +intros x Hx. +apply Rlt_le_trans with (ulp beta fexp x)%R. +now apply ulp_error. +unfold ulp, canonic_exp. +assert (Hx': (x <> 0)%R). +intros H. +apply Rlt_not_le with (2 := Hx). +rewrite H, Rabs_R0. +apply bpow_gt_0. +destruct (ln_beta beta x) as (ex, He). +simpl. +specialize (He Hx'). +apply Rle_trans with (bpow (-p + 1) * bpow (ex - 1))%R. +rewrite <- bpow_plus. +apply bpow_le. +assert (emin < ex)%Z. +apply (lt_bpow beta). +apply Rle_lt_trans with (2 := proj2 He). +exact Hx. +generalize (Hmin ex). +omega. +apply Rmult_le_compat_l. +apply bpow_ge_0. +apply He. +Qed. + +(** 1+#ε# property in any rounding *) +Theorem relative_error_ex : + forall x, + (bpow emin <= Rabs x)%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 x Hx. +apply relative_error_lt_conversion... +apply bpow_gt_0. +now apply relative_error. +Qed. + +Theorem relative_error_F2R_emin : + forall m, let x := F2R (Float beta m emin) in + (x <> 0)%R -> + (Rabs (round beta fexp rnd x - x) < bpow (-p + 1) * Rabs x)%R. +Proof. +intros m x Hx. +apply relative_error. +unfold x. +rewrite <- F2R_Zabs. +apply bpow_le_F2R. +apply F2R_lt_reg with beta emin. +rewrite F2R_0, F2R_Zabs. +now apply Rabs_pos_lt. +Qed. + +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... +apply bpow_gt_0. +now apply relative_error_F2R_emin. +Qed. + +Theorem relative_error_round : + (0 < p)%Z -> + forall x, + (bpow emin <= Rabs x)%R -> + (Rabs (round beta fexp rnd x - x) < bpow (-p + 1) * Rabs (round beta fexp rnd x))%R. +Proof with auto with typeclass_instances. +intros Hp x Hx. +apply Rlt_le_trans with (ulp beta fexp x)%R. +now apply ulp_error. +assert (Hx': (x <> 0)%R). +intros H. +apply Rlt_not_le with (2 := Hx). +rewrite H, Rabs_R0. +apply bpow_gt_0. +unfold ulp, canonic_exp. +destruct (ln_beta beta x) as (ex, He). +simpl. +specialize (He Hx'). +assert (He': (emin < ex)%Z). +apply (lt_bpow beta). +apply Rle_lt_trans with (2 := proj2 He). +exact Hx. +apply Rle_trans with (bpow (-p + 1) * bpow (ex - 1))%R. +rewrite <- bpow_plus. +apply bpow_le. +generalize (Hmin ex). +omega. +apply Rmult_le_compat_l. +apply bpow_ge_0. +generalize He. +apply round_abs_abs... +clear rnd valid_rnd x Hx Hx' He. +intros rnd valid_rnd x Hx. +rewrite <- (round_generic beta fexp rnd (bpow (ex - 1))). +now apply round_le. +apply generic_format_bpow. +ring_simplify (ex - 1 + 1)%Z. +generalize (Hmin ex). +omega. +Qed. + +Theorem relative_error_round_F2R_emin : + (0 < p)%Z -> + forall m, let x := F2R (Float beta m emin) in + (x <> 0)%R -> + (Rabs (round beta fexp rnd x - x) < bpow (-p + 1) * Rabs (round beta fexp rnd x))%R. +Proof. +intros Hp m x Hx. +apply relative_error_round. +exact Hp. +unfold x. +rewrite <- F2R_Zabs. +apply bpow_le_F2R. +apply F2R_lt_reg with beta emin. +rewrite F2R_0, F2R_Zabs. +now apply Rabs_pos_lt. +Qed. + +Variable choice : Z -> bool. + +Theorem relative_error_N : + forall x, + (bpow emin <= Rabs x)%R -> + (Rabs (round beta fexp (Znearest choice) x - x) <= /2 * bpow (-p + 1) * Rabs x)%R. +Proof. +intros x Hx. +apply Rle_trans with (/2 * ulp beta fexp x)%R. +now apply ulp_half_error. +rewrite Rmult_assoc. +apply Rmult_le_compat_l. +apply Rlt_le. +apply Rinv_0_lt_compat. +now apply (Z2R_lt 0 2). +assert (Hx': (x <> 0)%R). +intros H. +apply Rlt_not_le with (2 := Hx). +rewrite H, Rabs_R0. +apply bpow_gt_0. +unfold ulp, canonic_exp. +destruct (ln_beta beta x) as (ex, He). +simpl. +specialize (He Hx'). +apply Rle_trans with (bpow (-p + 1) * bpow (ex - 1))%R. +rewrite <- bpow_plus. +apply bpow_le. +assert (emin < ex)%Z. +apply (lt_bpow beta). +apply Rle_lt_trans with (2 := proj2 He). +exact Hx. +generalize (Hmin ex). +omega. +apply Rmult_le_compat_l. +apply bpow_ge_0. +apply He. +Qed. + +(** 1+#ε# property in rounding to nearest *) +Theorem relative_error_N_ex : + forall x, + (bpow emin <= Rabs x)%R -> + exists eps, + (Rabs eps <= /2 * bpow (-p + 1))%R /\ round beta fexp (Znearest choice) x = (x * (1 + eps))%R. +Proof with auto with typeclass_instances. +intros x Hx. +apply relative_error_le_conversion... +apply Rlt_le. +apply Rmult_lt_0_compat. +apply Rinv_0_lt_compat. +now apply (Z2R_lt 0 2). +apply bpow_gt_0. +now apply relative_error_N. +Qed. + +Theorem relative_error_N_F2R_emin : + forall m, let x := F2R (Float beta m emin) in + (Rabs (round beta fexp (Znearest choice) x - x) <= /2 * bpow (-p + 1) * Rabs x)%R. +Proof with auto with typeclass_instances. +intros m x. +destruct (Req_dec x 0) as [Hx|Hx]. +(* . *) +rewrite Hx, round_0... +unfold Rminus. +rewrite Rplus_0_l, Rabs_Ropp, Rabs_R0. +rewrite Rmult_0_r. +apply Rle_refl. +(* . *) +apply relative_error_N. +unfold x. +rewrite <- F2R_Zabs. +apply bpow_le_F2R. +apply F2R_lt_reg with beta emin. +rewrite F2R_0, F2R_Zabs. +now apply Rabs_pos_lt. +Qed. + +Theorem relative_error_N_F2R_emin_ex : + forall m, let x := F2R (Float beta m emin) in + exists eps, + (Rabs eps <= /2 * bpow (-p + 1))%R /\ round beta fexp (Znearest choice) x = (x * (1 + eps))%R. +Proof with auto with typeclass_instances. +intros m x. +apply relative_error_le_conversion... +apply Rlt_le. +apply Rmult_lt_0_compat. +apply Rinv_0_lt_compat. +now apply (Z2R_lt 0 2). +apply bpow_gt_0. +now apply relative_error_N_F2R_emin. +Qed. + +Theorem relative_error_N_round : + (0 < p)%Z -> + forall x, + (bpow emin <= Rabs x)%R -> + (Rabs (round beta fexp (Znearest choice) x - x) <= /2 * bpow (-p + 1) * Rabs (round beta fexp (Znearest choice) x))%R. +Proof with auto with typeclass_instances. +intros Hp x Hx. +apply Rle_trans with (/2 * ulp beta fexp x)%R. +now apply ulp_half_error. +rewrite Rmult_assoc. +apply Rmult_le_compat_l. +apply Rlt_le. +apply Rinv_0_lt_compat. +now apply (Z2R_lt 0 2). +assert (Hx': (x <> 0)%R). +intros H. +apply Rlt_not_le with (2 := Hx). +rewrite H, Rabs_R0. +apply bpow_gt_0. +unfold ulp, canonic_exp. +destruct (ln_beta beta x) as (ex, He). +simpl. +specialize (He Hx'). +assert (He': (emin < ex)%Z). +apply (lt_bpow beta). +apply Rle_lt_trans with (2 := proj2 He). +exact Hx. +apply Rle_trans with (bpow (-p + 1) * bpow (ex - 1))%R. +rewrite <- bpow_plus. +apply bpow_le. +generalize (Hmin ex). +omega. +apply Rmult_le_compat_l. +apply bpow_ge_0. +generalize He. +apply round_abs_abs... +clear rnd valid_rnd x Hx Hx' He. +intros rnd valid_rnd x Hx. +rewrite <- (round_generic beta fexp rnd (bpow (ex - 1))). +now apply round_le. +apply generic_format_bpow. +ring_simplify (ex - 1 + 1)%Z. +generalize (Hmin ex). +omega. +Qed. + +Theorem relative_error_N_round_F2R_emin : + (0 < p)%Z -> + forall m, let x := F2R (Float beta m emin) in + (Rabs (round beta fexp (Znearest choice) x - x) <= /2 * bpow (-p + 1) * Rabs (round beta fexp (Znearest choice) x))%R. +Proof with auto with typeclass_instances. +intros Hp m x. +destruct (Req_dec x 0) as [Hx|Hx]. +(* . *) +rewrite Hx, round_0... +unfold Rminus. +rewrite Rplus_0_l, Rabs_Ropp, Rabs_R0. +rewrite Rmult_0_r. +apply Rle_refl. +(* . *) +apply relative_error_N_round with (1 := Hp). +unfold x. +rewrite <- F2R_Zabs. +apply bpow_le_F2R. +apply F2R_lt_reg with beta emin. +rewrite F2R_0, F2R_Zabs. +now apply Rabs_pos_lt. +Qed. + +End Fprop_relative_generic. + +Section Fprop_relative_FLT. + +Variable emin prec : Z. +Variable Hp : Zlt 0 prec. + +Lemma relative_error_FLT_aux : + forall k, (emin + prec - 1 < k)%Z -> (prec <= k - FLT_exp emin prec k)%Z. +Proof. +intros k Hk. +unfold FLT_exp. +generalize (Zmax_spec (k - prec) emin). +omega. +Qed. + +Variable rnd : R -> Z. +Context { valid_rnd : Valid_rnd rnd }. + +Theorem relative_error_FLT_F2R_emin : + forall m, let x := F2R (Float beta m (emin + prec - 1)) 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 Hx. +apply relative_error_F2R_emin... +apply relative_error_FLT_aux. +Qed. + +Theorem relative_error_FLT : + forall x, + (bpow (emin + prec - 1) <= Rabs x)%R -> + (Rabs (round beta (FLT_exp emin prec) rnd x - x) < bpow (-prec + 1) * Rabs x)%R. +Proof with auto with typeclass_instances. +intros x Hx. +apply relative_error with (emin + prec - 1)%Z... +apply relative_error_FLT_aux. +Qed. + +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... +apply bpow_gt_0. +now apply relative_error_FLT_F2R_emin. +Qed. + +(** 1+#ε# property in any rounding in FLT *) +Theorem relative_error_FLT_ex : + forall x, + (bpow (emin + prec - 1) <= Rabs x)%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 x Hx. +apply relative_error_lt_conversion... +apply bpow_gt_0. +now apply relative_error_FLT. +Qed. + +Variable choice : Z -> bool. + +Theorem relative_error_N_FLT : + forall x, + (bpow (emin + prec - 1) <= Rabs x)%R -> + (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 x Hx. +apply relative_error_N with (emin + prec - 1)%Z... +apply relative_error_FLT_aux. +Qed. + +(** 1+#ε# property in rounding to nearest in FLT *) +Theorem relative_error_N_FLT_ex : + forall x, + (bpow (emin + prec - 1) <= Rabs x)%R -> + 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 x Hx. +apply relative_error_le_conversion... +apply Rlt_le. +apply Rmult_lt_0_compat. +apply Rinv_0_lt_compat. +now apply (Z2R_lt 0 2). +apply bpow_gt_0. +now apply relative_error_N_FLT. +Qed. + +Theorem relative_error_N_FLT_round : + forall x, + (bpow (emin + prec - 1) <= Rabs x)%R -> + (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 x Hx. +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 + prec - 1)) 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. +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 + prec - 1)) 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 Rlt_le. +apply Rmult_lt_0_compat. +apply Rinv_0_lt_compat. +now apply (Z2R_lt 0 2). +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 + 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. +Proof with auto with typeclass_instances. +intros m x. +apply relative_error_N_round_F2R_emin... +apply relative_error_FLT_aux. +Qed. + +Lemma error_N_FLT_aux : + forall x, + (0 < x)%R -> + exists eps, exists eta, + (Rabs eps <= /2 * bpow (-prec + 1))%R /\ + (Rabs eta <= /2 * bpow (emin))%R /\ + (eps*eta=0)%R /\ + round beta (FLT_exp emin prec) (Znearest choice) x = (x * (1 + eps) + eta)%R. +Proof. +intros x Hx2. +case (Rle_or_lt (bpow (emin+prec)) x); intros Hx. +(* *) +destruct relative_error_N_ex with (FLT_exp emin prec) (emin+prec)%Z prec choice x + as (eps,(Heps1,Heps2)). +now apply FLT_exp_valid. +intros; unfold FLT_exp. +rewrite Zmax_left; omega. +rewrite Rabs_right;[assumption|apply Rle_ge; now left]. +exists eps; exists 0%R. +split;[assumption|split]. +rewrite Rabs_R0; apply Rmult_le_pos. +auto with real. +apply bpow_ge_0. +split;[apply Rmult_0_r|idtac]. +now rewrite Rplus_0_r. +(* *) +exists 0%R; exists (round beta (FLT_exp emin prec) (Znearest choice) x - x)%R. +split. +rewrite Rabs_R0; apply Rmult_le_pos. +auto with real. +apply bpow_ge_0. +split. +apply Rle_trans with (/2*ulp beta (FLT_exp emin prec) x)%R. +apply ulp_half_error. +now apply FLT_exp_valid. +apply Rmult_le_compat_l; auto with real. +unfold ulp. +apply bpow_le. +unfold FLT_exp, canonic_exp. +rewrite Zmax_right. +omega. +destruct (ln_beta beta x) as (e,He); simpl. +assert (e-1 < emin+prec)%Z. +apply (lt_bpow beta). +apply Rle_lt_trans with (2:=Hx). +rewrite <- (Rabs_right x). +apply He; auto with real. +apply Rle_ge; now left. +omega. +split;ring. +Qed. + +End Fprop_relative_FLT. + +Section Fprop_relative_FLX. + +Variable prec : Z. +Variable Hp : Zlt 0 prec. + +Lemma relative_error_FLX_aux : + forall k, (prec <= k - FLX_exp prec k)%Z. +Proof. +intros k. +unfold FLX_exp. +omega. +Qed. + +Variable rnd : R -> Z. +Context { valid_rnd : Valid_rnd rnd }. + +Theorem relative_error_FLX : + forall x, + (x <> 0)%R -> + (Rabs (round beta (FLX_exp prec) rnd x - x) < bpow (-prec + 1) * Rabs x)%R. +Proof with auto with typeclass_instances. +intros x Hx. +destruct (ln_beta beta x) as (ex, He). +specialize (He Hx). +apply relative_error with (ex - 1)%Z... +intros k _. +apply relative_error_FLX_aux. +apply He. +Qed. + +(** 1+#ε# property in any rounding in FLX *) +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... +apply bpow_gt_0. +now apply relative_error_FLX. +Qed. + +Theorem relative_error_FLX_round : + forall x, + (x <> 0)%R -> + (Rabs (round beta (FLX_exp prec) rnd x - x) < bpow (-prec + 1) * Rabs (round beta (FLX_exp prec) rnd x))%R. +Proof with auto with typeclass_instances. +intros x Hx. +destruct (ln_beta beta x) as (ex, He). +specialize (He Hx). +apply relative_error_round with (ex - 1)%Z... +intros k _. +apply relative_error_FLX_aux. +apply He. +Qed. + +Variable choice : Z -> bool. + +Theorem relative_error_N_FLX : + forall x, + (Rabs (round beta (FLX_exp prec) (Znearest choice) x - x) <= /2 * bpow (-prec + 1) * Rabs x)%R. +Proof with auto with typeclass_instances. +intros x. +destruct (Req_dec x 0) as [Hx|Hx]. +(* . *) +rewrite Hx, round_0... +unfold Rminus. +rewrite Rplus_0_l, Rabs_Ropp, Rabs_R0. +rewrite Rmult_0_r. +apply Rle_refl. +(* . *) +destruct (ln_beta beta x) as (ex, He). +specialize (He Hx). +apply relative_error_N with (ex - 1)%Z... +intros k _. +apply relative_error_FLX_aux. +apply He. +Qed. + +(** 1+#ε# property in rounding to nearest in FLX *) +Theorem relative_error_N_FLX_ex : + forall x, + exists eps, + (Rabs eps <= /2 * bpow (-prec + 1))%R /\ round beta (FLX_exp prec) (Znearest choice) x = (x * (1 + eps))%R. +Proof with auto with typeclass_instances. +intros x. +apply relative_error_le_conversion... +apply Rlt_le. +apply Rmult_lt_0_compat. +apply Rinv_0_lt_compat. +now apply (Z2R_lt 0 2). +apply bpow_gt_0. +now apply relative_error_N_FLX. +Qed. + +Theorem relative_error_N_FLX_round : + forall x, + (Rabs (round beta (FLX_exp prec) (Znearest choice) x - x) <= /2 * bpow (-prec + 1) * Rabs (round beta (FLX_exp prec) (Znearest choice) x))%R. +Proof with auto with typeclass_instances. +intros x. +destruct (Req_dec x 0) as [Hx|Hx]. +(* . *) +rewrite Hx, round_0... +unfold Rminus. +rewrite Rplus_0_l, Rabs_Ropp, Rabs_R0. +rewrite Rmult_0_r. +apply Rle_refl. +(* . *) +destruct (ln_beta beta x) as (ex, He). +specialize (He Hx). +apply relative_error_N_round with (ex - 1)%Z. +now apply FLX_exp_valid. +intros k _. +apply relative_error_FLX_aux. +exact Hp. +apply He. +Qed. + +End Fprop_relative_FLX. + +End Fprop_relative. \ No newline at end of file -- cgit