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

Copyright (C) 2010 Sylvie Boldo
BOLDO Sylvie's avatar
BOLDO Sylvie committed
6
#<br />#
7 8 9 10 11 12 13 14 15 16 17 18 19
Copyright (C) 2010 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.
*)

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 32 33 34 35
Variable prec : Z.
Variable Hp : Zlt 0 prec.

Notation format := (generic_format beta (FLX_exp prec)).
Notation cexp := (canonic_exponent beta (FLX_exp prec)).

36
Variable choice : Z -> bool.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
37

38 39 40 41 42 43 44 45
Theorem format_add: 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
  -> format (x+y)%R.
intros 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, <- plus_F2R.
apply generic_format_canonic_exponent.
46 47 48
case_eq (Fplus beta fx fy).
intros mz ez Hz.
rewrite <- Hz.
49 50 51 52 53 54 55 56 57 58 59 60
apply Zle_trans with (Zmin (Fexp fx) (Fexp fy)).
unfold canonic_exponent, FLX_exp.
rewrite plus_F2R, <- Hx, <- Hy.
destruct (ln_beta beta (x+y)); simpl.
specialize (a H).
apply Zmin_case.
apply Zplus_le_reg_l with prec; ring_simplify.
apply (bpow_lt_bpow beta).
now apply Rle_lt_trans with (1:=proj1 a).
apply Zplus_le_reg_l with prec; ring_simplify.
apply (bpow_lt_bpow beta).
now apply Rle_lt_trans with (1:=proj1 a).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
61 62
rewrite <- Fexp_Fplus, Hz.
apply Zle_refl.
63
Qed.
64

65 66 67 68 69
Theorem format_nx: forall x, format x -> exists fx:float beta, (x=F2R fx)%R /\ Fexp fx = cexp x.
intros x; unfold generic_format.
exists (Float beta (Ztrunc (scaled_mantissa beta (FLX_exp prec) x)) (cexp x)).
split; auto.
Qed.
70

BOLDO Sylvie's avatar
BOLDO Sylvie committed
71
(** Remainder of the division in FLX *)
72 73
Theorem div_error_FLX :
  forall Zrnd x y,
BOLDO Sylvie's avatar
BOLDO Sylvie committed
74
  format x -> format y ->
75
  format (x - round beta (FLX_exp prec) Zrnd (x/y) * y)%R.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
76
Proof.
77 78 79
intros Zrnd x y Hx Hy.
destruct (Req_dec y 0) as [Zy|Zy].
now rewrite Zy, Rmult_0_r, Rminus_0_r.
80
case (Req_dec (round beta (FLX_exp prec) Zrnd (x/y)) 0); intros Hr.
81
rewrite Hr; ring_simplify (x-0*y)%R; assumption.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
82 83 84 85
assert (Zx: x <> R0).
contradict Hr.
rewrite Hr.
unfold Rdiv.
86
now rewrite Rmult_0_l, round_0.
87 88
destruct (format_nx x Hx) as (fx,(Hx1,Hx2)).
destruct (format_nx y Hy) as (fy,(Hy1,Hy2)).
89 90
destruct (format_nx (round beta (FLX_exp prec) Zrnd (x / y))) as (fr,(Hr1,Hr2)).
apply generic_format_round.
91 92 93 94
now apply FLX_exp_correct.
unfold Rminus; apply format_add with fx (Fopp beta (Fmult beta fr fy)); trivial.
now rewrite Fopp_F2R,mult_F2R, <- Hr1, <- Hy1.
(* *)
Guillaume Melquiond's avatar
Guillaume Melquiond committed
95 96 97 98 99 100 101 102 103 104
destruct (relative_error_FLX_ex beta prec Hp Zrnd (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.
105
now apply Rabs_pos_lt.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
106 107
apply Rlt_le_trans with (1 := Heps1).
change R1 with (bpow 0).
108
apply bpow_le.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
109 110 111 112 113 114 115 116 117
clear -Hp. omega.
rewrite Rmult_1_r.
rewrite Hx2.
unfold canonic_exponent.
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) Zrnd (x / y) * y))%R with
  (y * (-(round beta (FLX_exp prec) Zrnd (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.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
135
apply ulp_error_f.
136 137
now apply FLX_exp_correct.
clear; intros; unfold FLX_exp; omega.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
138 139 140
exact Hr.
unfold ulp; apply f_equal.
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 145
apply Rmult_le_compat_r.
apply bpow_ge_0.
rewrite Hy2; unfold canonic_exponent, 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
  format (x - Rsqr (round beta (FLX_exp prec) (rndN choice) (sqrt x)))%R.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
157
Proof.
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) (rndN choice) (sqrt x)) 0); intros Hr.
172 173
rewrite Hr; unfold Rsqr; ring_simplify (x-0*0)%R; assumption.
destruct (format_nx x Hx) as (fx,(Hx1,Hx2)).
174 175
destruct (format_nx (round beta (FLX_exp prec) (rndN choice) (sqrt x))) as (fr,(Hr1,Hr2)).
apply generic_format_round.
176 177 178 179 180 181
now apply FLX_exp_correct.
unfold Rminus; apply format_add with fx (Fopp beta (Fmult beta fr fr)); trivial.
unfold Rsqr; now rewrite Fopp_F2R,mult_F2R, <- Hr1.
(* *) 
apply Rle_lt_trans with x.
apply Rabs_Rminus_pos.
182
apply Rle_0_sqr.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
183 184 185 186 187 188 189
destruct (relative_error_N_FLX_ex beta prec Hp 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.
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).
229
rewrite Hx2; unfold canonic_exponent, 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 249 250 251 252 253 254 255
(* *)
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; trivial.
now apply FLX_exp_correct.
clear; intros; unfold FLX_exp; omega.
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.
256
rewrite bpow_plus, Rmult_assoc.
257 258 259 260 261 262 263 264 265 266 267
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_exponent; rewrite Hr1.
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 round_monotone_abs_l.
292 293 294
now apply FLX_exp_correct.
apply generic_format_bpow.
unfold FLX_exp; omega.
295
apply Es.
296
apply Rlt_le_trans with (1:=H).
297 298
apply bpow_le.
omega.
299
now apply Rlt_le.
300
Qed.
301

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