Fprop_div_sqrt_error.v 9.21 KB
Newer Older
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1
(**
2 3 4
This file is part of the Flocq formalization of floating-point
arithmetic in Coq: http://flocq.gforge.inria.fr/

5
Copyright (C) 2010-2013 Sylvie Boldo
BOLDO Sylvie's avatar
BOLDO Sylvie committed
6
#<br />#
7
Copyright (C) 2010-2013 Guillaume Melquiond
8 9 10 11 12 13 14 15 16 17 18 19

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.
*)

BOLDO Sylvie's avatar
BOLDO Sylvie committed
20
(** * Remainder of the division and square root are in the FLX format *)
BOLDO Sylvie's avatar
BOLDO Sylvie committed
21 22
Require Import Fcore.
Require Import Fcalc_ops.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
23
Require Import Fprop_relative.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
24 25 26 27 28 29

Section Fprop_divsqrt_error.

Variable beta : radix.
Notation bpow e := (bpow beta e).

30 31
Variable prec : Z.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
32 33 34
Theorem generic_format_plus_prec:
  forall fexp, (forall e, (fexp e  <= e - prec)%Z) ->
  forall x y (fx fy: float beta),
35
  (x = F2R fx)%R -> (y = F2R fy)%R -> (Rabs (x+y) < bpow (prec+Fexp fx))%R -> (Rabs (x+y) < bpow (prec+Fexp fy))%R
BOLDO Sylvie's avatar
BOLDO Sylvie committed
36 37
  -> generic_format beta fexp (x+y)%R.
intros fexp Hfexp x y fx fy Hx Hy H1 H2.
38 39
case (Req_dec (x+y) 0); intros H.
rewrite H; apply generic_format_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
40
rewrite Hx, Hy, <- F2R_plus.
41 42
apply generic_format_F2R.
intros _.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
43 44 45
case_eq (Fplus beta fx fy).
intros mz ez Hz.
rewrite <- Hz.
46
apply Zle_trans with (Zmin (Fexp fx) (Fexp fy)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
47 48 49
rewrite F2R_plus, <- Hx, <- Hy.
unfold canonic_exp.
apply Zle_trans with (1:=Hfexp _).
50
apply Zplus_le_reg_l with prec; ring_simplify.
51
apply ln_beta_le_bpow with (1 := H).
52
now apply Zmin_case.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
53 54
rewrite <- Fexp_Fplus, Hz.
apply Zle_refl.
55
Qed.
56

57
Theorem ex_Fexp_canonic: forall fexp, forall x, generic_format beta fexp x
BOLDO Sylvie's avatar
BOLDO Sylvie committed
58 59 60
  -> 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)).
61 62
split; auto.
Qed.
63

BOLDO Sylvie's avatar
BOLDO Sylvie committed
64 65 66 67 68 69 70 71 72

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.


BOLDO Sylvie's avatar
BOLDO Sylvie committed
73
(** Remainder of the division in FLX *)
74
Theorem div_error_FLX :
75
  forall rnd { Zrnd : Valid_rnd rnd } x y,
BOLDO Sylvie's avatar
BOLDO Sylvie committed
76
  format x -> format y ->
77 78 79
  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.
80 81
destruct (Req_dec y 0) as [Zy|Zy].
now rewrite Zy, Rmult_0_r, Rminus_0_r.
82
destruct (Req_dec (round beta (FLX_exp prec) rnd (x/y)) 0) as [Hr|Hr].
83
rewrite Hr; ring_simplify (x-0*y)%R; assumption.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
84 85 86 87
assert (Zx: x <> R0).
contradict Hr.
rewrite Hr.
unfold Rdiv.
88
now rewrite Rmult_0_l, round_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
89 90 91
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)).
92
apply generic_format_round...
BOLDO Sylvie's avatar
BOLDO Sylvie committed
93 94 95
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.
96
(* *)
97
destruct (relative_error_FLX_ex beta prec (prec_gt_0 prec) rnd (x / y)%R) as (eps,(Heps1,Heps2)).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
98 99 100 101 102 103
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.
104
now apply Rabs_pos_lt.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
105 106
apply Rlt_le_trans with (1 := Heps1).
change R1 with (bpow 0).
107
apply bpow_le.
108 109
generalize (prec_gt_0 prec).
clear ; omega.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
110 111
rewrite Rmult_1_r.
rewrite Hx2.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
112
unfold canonic_exp.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
113 114 115 116 117
destruct (ln_beta beta x) as (ex, Hex).
simpl.
specialize (Hex Zx).
apply Rlt_le.
apply Rlt_le_trans with (1 := proj2 Hex).
118
apply bpow_le.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
119 120 121
unfold FLX_exp.
ring_simplify.
apply Zle_refl.
122 123 124
(* *)
replace (Fexp (Fopp beta (Fmult beta fr fy))) with (Fexp fr + Fexp fy)%Z.
2: unfold Fopp, Fmult; destruct fr; destruct fy; now simpl.
125 126
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.
127 128 129 130 131 132
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.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
133
replace (bpow (Fexp fr)) with (ulp beta (FLX_exp prec) (F2R fr)).
134
rewrite <- Hr1.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
135
apply error_lt_ulp_round...
136 137 138 139
apply Rmult_integral_contrapositive_currified; try apply Rinv_neq_0_compat; assumption.
rewrite ulp_neq_0.
2: now rewrite <- Hr1.
apply f_equal.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
140
now rewrite Hr2, <- Hr1.
141
replace (prec+(Fexp fr+Fexp fy))%Z with ((prec+Fexp fy)+Fexp fr)%Z by ring.
142
rewrite bpow_plus.
143 144
apply Rmult_le_compat_r.
apply bpow_ge_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
145
rewrite Hy2; unfold canonic_exp, FLX_exp.
146
ring_simplify (prec + (ln_beta beta y - prec))%Z.
147 148 149 150
destruct (ln_beta beta y); simpl.
left; now apply a.
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
151
(** Remainder of the square in FLX (with p>1) and rounding to nearest *)
Guillaume Melquiond's avatar
Guillaume Melquiond committed
152 153
Variable Hp1 : Zlt 1 prec.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
154
Theorem sqrt_error_FLX_N :
155
  forall x, format x ->
156 157
  format (x - Rsqr (round beta (FLX_exp prec) (Znearest choice) (sqrt x)))%R.
Proof with auto with typeclass_instances.
158 159 160 161
intros x Hx.
destruct (total_order_T x 0) as [[Hxz|Hxz]|Hxz].
unfold sqrt.
destruct (Rcase_abs x).
162
rewrite round_0...
163 164 165 166
unfold Rsqr.
now rewrite Rmult_0_l, Rminus_0_r.
elim (Rlt_irrefl 0).
now apply Rgt_ge_trans with x.
167
rewrite Hxz, sqrt_0, round_0...
168 169 170
unfold Rsqr.
rewrite Rmult_0_l, Rminus_0_r.
apply generic_format_0.
171
case (Req_dec (round beta (FLX_exp prec) (Znearest choice) (sqrt x)) 0); intros Hr.
172
rewrite Hr; unfold Rsqr; ring_simplify (x-0*0)%R; assumption.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
173 174
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)).
175
apply generic_format_round...
BOLDO Sylvie's avatar
BOLDO Sylvie committed
176 177 178
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.
179
(* *)
180
apply Rle_lt_trans with x.
181
apply Rabs_minus_le.
182
apply Rle_0_sqr.
183
destruct (relative_error_N_FLX_ex beta prec (prec_gt_0 prec) choice (sqrt x)) as (eps,(Heps1,Heps2)).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
184 185 186 187 188 189
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.
190
apply Rsqr_le_abs_1.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
191 192 193 194 195 196 197 198 199 200
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)).
201
apply bpow_le.
202
omega.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
203 204
replace (2 * (-1 + 5 / 4))%R with (/2)%R by field.
apply Rinv_le.
205 206
now apply (Z2R_lt 0 2).
apply (Z2R_le 2).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
207 208
unfold Zpower_pos. simpl.
rewrite Zmult_1_r.
209 210
apply Zle_bool_imp_le.
apply beta.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
211 212 213
apply Rgt_not_eq.
now apply (Z2R_lt 0 2).
unfold Rdiv.
214
apply Rmult_le_pos.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
215 216 217 218 219 220 221 222 223 224 225 226 227 228
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).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
229
rewrite Hx2; unfold canonic_exp, FLX_exp.
230
ring_simplify (prec + (ln_beta beta x - prec))%Z.
231 232
destruct (ln_beta beta x); simpl.
rewrite <- (Rabs_right x).
233 234 235
apply a.
now apply Rgt_not_eq.
now apply Rgt_ge.
236 237 238 239 240 241 242 243 244 245 246 247 248
(* *)
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.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
249
apply error_le_half_ulp_round...
250 251 252
right; rewrite ulp_neq_0.
2: now rewrite <- Hr1.
apply f_equal.
253
rewrite Hr2, <- Hr1; trivial.
254 255
rewrite Rmult_assoc, Rmult_comm.
replace (prec+(Fexp fr+Fexp fr))%Z with (Fexp fr + (prec+Fexp fr))%Z by ring.
256
rewrite bpow_plus, Rmult_assoc.
257 258 259 260 261 262 263 264 265
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.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
266
rewrite Hr2; unfold canonic_exp; rewrite Hr1.
267
unfold FLX_exp.
268
ring_simplify (prec + (ln_beta beta (F2R fr) - prec))%Z.
269 270 271 272 273 274 275 276 277 278
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.
279 280
apply Rgt_not_eq.
now apply sqrt_lt_R0.
281 282 283 284
destruct (ln_beta beta (sqrt x)) as (es,Es).
specialize (Es H0).
apply Rle_trans with (bpow es).
now apply Rlt_le.
285
apply bpow_le.
286
case (Zle_or_lt es (prec + Fexp fr)) ; trivial.
287
intros H1.
288
absurd (Rabs (F2R fr) < bpow (es - 1))%R.
289 290
apply Rle_not_lt.
rewrite <- Hr1.
291
apply abs_round_ge_generic...
292 293
apply generic_format_bpow.
unfold FLX_exp; omega.
294
apply Es.
295
apply Rlt_le_trans with (1:=H).
296 297
apply bpow_le.
omega.
298
now apply Rlt_le.
299
Qed.
300

BOLDO Sylvie's avatar
BOLDO Sylvie committed
301
End Fprop_divsqrt_error.