Average.v 54.1 KB
Newer Older
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
(**
This example is part of the Flocq formalization of floating-point
arithmetic in Coq: http://flocq.gforge.inria.fr/

Copyright (C) 2014-2018 Sylvie Boldo

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

Guillaume Melquiond's avatar
Guillaume Melquiond committed
18
Require Import Reals Fourier Psatz.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
19
From Flocq Require Import Core Plus_error.
20 21 22 23 24

Open Scope R_scope.

Section av1.

25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52

Lemma Rmin_Rmax_overflow: forall x y z M, Rabs x <= M -> Rabs y <= M ->
        Rmin x y <= z <= Rmax x y -> Rabs z <= M.
Proof.
intros x y z M Hx Hy H.
case (Rle_or_lt 0 z); intros Hz.
rewrite Rabs_right.
apply Rle_trans with (1:=proj2 H).
generalize (proj2 H).
apply Rmax_case_strong.
intros; apply Rle_trans with (2:=Hx).
apply RRle_abs.
intros; apply Rle_trans with (2:=Hy).
apply RRle_abs.
now apply Rle_ge.
rewrite Rabs_left; try assumption.
apply Rle_trans with (Rmax (-x) (-y)).
rewrite Rmax_opp.
apply Ropp_le_contravar, H.
apply Rmax_case_strong.
intros; apply Rle_trans with (2:=Hx).
rewrite <- Rabs_Ropp.
apply RRle_abs.
intros; apply Rle_trans with (2:=Hy).
rewrite <- Rabs_Ropp.
apply RRle_abs.
Qed.

53 54 55 56 57 58 59 60

Definition radix2 := Build_radix 2 (refl_equal true).
Notation bpow e := (bpow radix2 e).

Variable emin prec : Z.
Context { prec_gt_0_ : Prec_gt_0 prec }.

Notation format := (generic_format radix2 (FLT_exp emin prec)).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
61
Notation round_flt :=(round radix2 (FLT_exp emin prec) ZnearestE).
62
Notation ulp_flt :=(ulp radix2 (FLT_exp emin prec)).
63
Notation cexp := (cexp radix2 (FLT_exp emin prec)).
64
Notation pred_flt := (pred radix2 (FLT_exp emin prec)).
65 66 67 68 69 70

Lemma FLT_format_double: forall u, format u -> format (2*u).
Proof with auto with typeclass_instances.
intros u Fu.
apply generic_format_FLT.
apply FLT_format_generic in Fu...
71
destruct Fu as [uf H1 H2 H3].
72 73 74 75
exists (Float radix2 (Fnum uf) (Fexp uf+1)).
rewrite H1; unfold F2R; simpl.
rewrite bpow_plus, bpow_1.
simpl;ring.
76 77 78
easy.
apply Zle_trans with (1:=H3).
apply Zle_succ.
79 80 81 82 83 84 85
Qed.

Lemma FLT_format_half: forall u, 
   format u -> bpow (prec+emin) <= Rabs u -> format (u/2).
Proof with auto with typeclass_instances.
intros u Fu H.
apply FLT_format_generic in Fu...
86
destruct Fu as [[n e] H1 H2 H3].
87 88 89 90 91
simpl in H1, H2, H3.
apply generic_format_FLT.
exists (Float radix2 n (e-1)).
rewrite H1; unfold F2R; simpl.
unfold Zminus; rewrite bpow_plus.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
92 93
change (bpow (-(1))) with (/2).
unfold Rdiv; ring.
94 95 96
easy.
cut (prec + emin < prec +e)%Z.
  simpl ; omega.
97 98 99 100 101 102 103 104
apply lt_bpow with radix2.
apply Rle_lt_trans with (1:=H).
rewrite H1; unfold F2R; simpl.
rewrite Rabs_mult; rewrite (Rabs_right (bpow e)).
2: apply Rle_ge, bpow_ge_0.
rewrite bpow_plus.
apply Rmult_lt_compat_r.
apply bpow_gt_0.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
105 106 107
rewrite <- abs_IZR.
rewrite <- IZR_Zpower.
now apply IZR_lt.
108
now apply Zlt_le_weak.
109 110
Qed.

111 112 113 114 115 116 117 118 119 120 121 122 123
Lemma FLT_round_half: forall z, bpow (prec+emin) <= Rabs z -> 
   round_flt (z/2)= round_flt z /2.
Proof with auto with typeclass_instances.
intros z Hz.
apply Rmult_eq_reg_l with 2.
2: apply sym_not_eq; auto with real.
apply trans_eq with (round_flt z).
2: field.
assert (z <> 0)%R.
intros K; contradict Hz.
rewrite K, Rabs_R0; apply Rlt_not_le.
apply bpow_gt_0.
assert (cexp (z/2) = cexp z -1)%Z.
124
assert (prec+emin < mag radix2 z)%Z.
125
apply lt_bpow with radix2.
126
destruct mag as (e,He); simpl.
127 128
apply Rle_lt_trans with (1:=Hz).
now apply He.
129
unfold cexp, FLT_exp.
130
replace ((mag radix2 (z/2))-prec)%Z with ((mag radix2 z -1) -prec)%Z.
131 132 133
rewrite Z.max_l; try omega.
rewrite Z.max_l; try omega.
apply Zplus_eq_compat; try reflexivity.
134 135
apply sym_eq, mag_unique.
destruct (mag radix2 z) as (e,He); simpl.
136 137 138 139 140 141 142 143 144
unfold Rdiv; rewrite Rabs_mult.
rewrite (Rabs_right (/2)).
split.
apply Rmult_le_reg_l with (bpow 1).
apply bpow_gt_0.
rewrite <- bpow_plus.
replace (1+(e-1-1))%Z with (e-1)%Z by ring.
apply Rle_trans with (Rabs z).
now apply He.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
145
change (bpow 1) with 2%R.
146 147 148 149 150 151
right; simpl; field.
apply Rmult_lt_reg_l with (bpow 1).
apply bpow_gt_0.
rewrite <- bpow_plus.
replace (1+(e-1))%Z with e by ring.
apply Rle_lt_trans with (Rabs z).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
152 153
change (bpow 1) with 2.
right; field.
154
now apply He.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
155
lra.
156 157 158 159 160 161 162
unfold round, scaled_mantissa, F2R.
rewrite H0; simpl.
rewrite Rmult_comm, Rmult_assoc.
apply f_equal2.
apply f_equal, f_equal.
replace (-(cexp z -1))%Z with (-cexp z +1)%Z by ring.
rewrite bpow_plus.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
163 164
change (bpow 1) with 2.
field.
165
unfold Zminus; rewrite bpow_plus.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
166 167
change (bpow (-(1))) with (/2).
field.
168
Qed.
169 170 171 172

Lemma FLT_ulp_le_id: forall u, bpow emin <= u -> ulp_flt u <= u.
Proof with auto with typeclass_instances.
intros u H.
173 174
rewrite ulp_neq_0.
2: apply Rgt_not_eq, Rlt_le_trans with (2:=H), bpow_gt_0.
175
case (Rle_or_lt (bpow (emin+prec-1)) u); intros Hu.
176 177
unfold ulp; rewrite cexp_FLT_FLX.
unfold cexp, FLX_exp.
178
destruct (mag radix2 u) as (e,He); simpl.
179 180 181 182 183 184 185 186 187 188 189 190
apply Rle_trans with (bpow (e-1)).
apply bpow_le.
unfold Prec_gt_0 in prec_gt_0_; omega.
rewrite <- (Rabs_right u).
apply He.
apply Rgt_not_eq, Rlt_gt.
apply Rlt_le_trans with (2:=Hu).
apply bpow_gt_0.
apply Rle_ge, Rle_trans with (2:=Hu), bpow_ge_0.
rewrite Rabs_right.
assumption.
apply Rle_ge, Rle_trans with (2:=Hu), bpow_ge_0.
191
unfold ulp; rewrite cexp_FLT_FIX.
192 193 194 195 196 197 198 199 200 201 202 203 204
apply H.
apply Rgt_not_eq, Rlt_gt.
apply Rlt_le_trans with (2:=H).
apply bpow_gt_0.
rewrite Rabs_right.
apply Rlt_le_trans with (1:=Hu).
apply bpow_le; omega.
apply Rle_ge, Rle_trans with (2:=H), bpow_ge_0.
Qed.



Lemma FLT_ulp_double: forall u, ulp_flt (2*u) <= 2*ulp_flt(u).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
205
Proof.
206
intros u.
207 208 209 210
case (Req_bool_spec u 0); intros Hu'.
rewrite Hu', Rmult_0_r.
rewrite <- (Rmult_1_l (ulp_flt 0)) at 1.
apply Rmult_le_compat_r.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
211
apply ulp_ge_0.
212 213
left; apply Rlt_plus_1.
rewrite 2!ulp_neq_0; trivial.
214
2: lra.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
215
change 2 at 2 with (bpow 1).
216
rewrite <- bpow_plus.
217 218 219
apply bpow_le.
case (Rle_or_lt (bpow (emin+prec-1)) (Rabs u)); intros Hu.
(* *)
220 221 222
rewrite cexp_FLT_FLX.
rewrite cexp_FLT_FLX; trivial.
unfold cexp, FLX_exp.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
223
change 2 with (bpow 1).
224
rewrite Rmult_comm, mag_mult_bpow.
225 226 227 228 229 230 231 232 233
omega.
intros H; contradict Hu.
apply Rlt_not_le; rewrite H, Rabs_R0.
apply bpow_gt_0.
apply Rle_trans with (1:=Hu).
rewrite Rabs_mult.
pattern (Rabs u) at 1; rewrite <- (Rmult_1_l (Rabs u)).
apply Rmult_le_compat_r.
apply Rabs_pos.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
234
rewrite <- (abs_IZR 2).
235
now apply IZR_le.
236
(* *)
237 238 239
rewrite cexp_FLT_FIX.
rewrite cexp_FLT_FIX; trivial.
unfold FIX_exp, cexp; omega.
240 241
apply Rlt_le_trans with (1:=Hu).
apply bpow_le; omega.
242
lra.
243
rewrite Rabs_mult.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
244
rewrite Rabs_pos_eq.
245
2: now apply IZR_le.
246
apply Rlt_le_trans with (2*bpow (emin + prec - 1)).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
247
apply Rmult_lt_compat_l with (1 := Rlt_0_2).
248
assumption.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
249
change 2 with (bpow 1).
250 251 252 253 254
rewrite <- bpow_plus.
apply bpow_le; omega.
Qed.


255
Lemma round_plus_small_id_aux: forall f h, format f -> (bpow (prec+emin) <= f) -> 0 < f 
256
   -> Rabs h <= /4* ulp_flt f -> round_flt (f+h) = f.
257 258 259 260 261 262
Proof with auto with typeclass_instances.
intros f h Ff H1 H2 Hh.
case (Rle_or_lt 0 h); intros H3;[destruct H3|idtac].
(* 0 < h *)
rewrite Rabs_right in Hh.
2: now apply Rle_ge, Rlt_le.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
263
apply round_N_eq_DN_pt with (f+ ulp_flt f)...
264
pattern f at 2; rewrite <- (round_DN_plus_eps_pos radix2 (FLT_exp emin prec) f) with (eps:=h); try assumption.
265
apply round_DN_pt...
266
now left.
267 268
split.
now left.
269
apply Rle_lt_trans with (1:=Hh).
270 271
rewrite <- (Rmult_1_l (ulp_flt f)) at 2.
apply Rmult_lt_compat_r.
272
rewrite ulp_neq_0; try now apply Rgt_not_eq.
273 274
apply bpow_gt_0.
fourier.
275
rewrite <- (round_UP_plus_eps_pos radix2 (FLT_exp emin prec) f) with (eps:=h); try assumption.
276
apply round_UP_pt...
277
now left.
278
split; trivial.
279
apply Rle_trans with (1:=Hh).
280 281
rewrite <- (Rmult_1_l (ulp_flt f)) at 2.
apply Rmult_le_compat_r.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
282
apply ulp_ge_0.
283 284 285 286
fourier.
apply Rplus_lt_reg_l with (-f); ring_simplify.
apply Rlt_le_trans with (/2*ulp_flt f).
2: right; field.
287
apply Rle_lt_trans with (1:=Hh).
288
apply Rmult_lt_compat_r.
289
rewrite ulp_neq_0; try now apply Rgt_not_eq.
290 291 292 293 294 295 296 297 298 299 300
apply bpow_gt_0.
fourier.
(* h = 0 *)
rewrite <- H, Rplus_0_r.
apply round_generic...
(* h < 0 *)
(*  - assertions *)
rewrite Rabs_left in Hh; try assumption.
assert (0 < pred_flt f).
apply Rlt_le_trans with (bpow emin).
apply bpow_gt_0.
301
apply pred_ge_gt...
302 303 304 305 306
apply FLT_format_bpow...
omega.
apply Rlt_le_trans with (2:=H1).
apply bpow_lt.
unfold Prec_gt_0 in prec_gt_0_; omega.
307 308
assert (M:(prec + emin +1 <= mag radix2 f)%Z).
apply mag_ge_bpow.
309 310 311 312 313
replace (prec+emin+1-1)%Z with (prec+emin)%Z by ring.
rewrite Rabs_right; try assumption.
apply Rle_ge; now left.
assert (T1:(ulp_flt (pred_flt f) = ulp_flt f) 
     \/ ( ulp_flt (pred_flt f) = /2* ulp_flt f 
314
               /\ f = bpow (mag radix2 f -1))).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
315
generalize H; rewrite pred_eq_pos; [idtac|now left].
316
unfold pred_pos; case Req_bool_spec; intros K HH.
317 318
(**)
right; split; try assumption.
319
rewrite ulp_neq_0;[idtac|now apply Rgt_not_eq].
320
apply trans_eq with (bpow (mag radix2 f- prec -1)).
321
apply f_equal.
322
unfold cexp.
323
apply trans_eq with (FLT_exp emin prec (mag radix2 f -1)%Z).
324 325 326 327
apply f_equal.
unfold FLT_exp.
rewrite Z.max_l.
2: omega.
328
apply mag_unique.
329 330
rewrite Rabs_right.
split.
331
apply Rplus_le_reg_l with (bpow (mag radix2 f -1-prec)).
332
ring_simplify.
333
apply Rle_trans with (bpow (mag radix2 f - 1 - 1) + bpow (mag radix2 f - 1 - 1)).
334 335 336
apply Rplus_le_compat_r.
apply bpow_le.
unfold Prec_gt_0 in prec_gt_0_; omega.
337
apply Rle_trans with (bpow 1*bpow (mag radix2 f - 1 - 1)).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
338 339
change (bpow 1) with 2.
right; ring.
340
rewrite <- bpow_plus.
341
apply Rle_trans with (bpow (mag radix2 f -1)).
342 343 344
apply bpow_le; omega.
rewrite <- K; now right.
rewrite <- K.
345
apply Rplus_lt_reg_l with (-f+bpow (mag radix2 f-1-prec)); ring_simplify.
346 347 348
apply bpow_gt_0.
apply Rle_ge.
rewrite K at 1.
349
apply Rplus_le_reg_l with (bpow (mag radix2 f - 1 - prec)).
350 351 352 353 354 355
ring_simplify.
apply bpow_le.
unfold Prec_gt_0 in prec_gt_0_; omega.
unfold FLT_exp.
rewrite Z.max_l;[ring|omega].
replace (/2) with (bpow (-1)) by reflexivity.
356 357
rewrite ulp_neq_0; try now apply Rgt_not_eq.
rewrite <- bpow_plus.
358
apply f_equal.
359
unfold cexp, FLT_exp.
360 361 362
rewrite Z.max_l;[ring|omega].
(**)
left.
363 364
assert (bpow (mag radix2 f -1) < f).
destruct (mag radix2 f); simpl in *.
365
destruct a.
366
now apply Rgt_not_eq.
367 368 369 370 371
rewrite Rabs_right in H0.
destruct H0; try assumption.
contradict H0.
now apply sym_not_eq.
apply Rle_ge; now left.
372
assert (bpow (mag radix2 f -1) + ulp_flt (bpow (mag radix2 f-1)) <= f).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
373
rewrite <- succ_eq_pos;[idtac|apply bpow_ge_0].
374
apply succ_le_lt...
375 376 377 378 379 380
apply FLT_format_bpow...
unfold Prec_gt_0 in prec_gt_0_;omega.
rewrite ulp_bpow in H4.
unfold FLT_exp in H4.
rewrite Z.max_l in H4.
2: omega.
381
replace (mag radix2 f - 1 + 1 - prec)%Z with  (mag radix2 f - prec)%Z in H4 by ring.
382 383
rewrite ulp_neq_0; try now apply Rgt_not_eq.
rewrite ulp_neq_0 at 2; try now apply Rgt_not_eq.
384
unfold cexp.
385
apply f_equal; apply f_equal.
386 387
replace (ulp_flt f) with (bpow (mag radix2 f -prec)).
apply mag_unique.
388 389
rewrite Rabs_right.
split.
390
apply Rplus_le_reg_l with (bpow (mag radix2 f -prec)).
391 392 393
ring_simplify.
apply Rle_trans with (2:=H4); right; ring.
apply Rlt_trans with f.
394
apply Rplus_lt_reg_l with (-f+bpow (mag radix2 f - prec)).
395 396 397
ring_simplify.
apply bpow_gt_0.
apply Rle_lt_trans with (1:=RRle_abs _).
398
apply bpow_mag_gt.
399
apply Rle_ge.
400
apply Rplus_le_reg_l with (bpow (mag radix2 f - prec)).
401 402 403 404
ring_simplify.
left; apply Rle_lt_trans with (2:=H0).
apply bpow_le.
unfold Prec_gt_0 in prec_gt_0_;omega.
405
rewrite ulp_neq_0; try now apply Rgt_not_eq.
406
unfold cexp, FLT_exp.
407 408 409 410 411 412
rewrite Z.max_l.
reflexivity.
omega.
assert (T: (ulp_flt (pred_flt f) = ulp_flt f \/ 
              (ulp_flt (pred_flt f) = / 2 * ulp_flt f /\ - h < / 4 * ulp_flt f))
         \/ (ulp_flt (pred_flt f) = / 2 * ulp_flt f /\
413
              f = bpow (mag radix2 f - 1) /\
414 415 416 417 418 419 420 421 422
              - h = / 4 * ulp_flt f) ).
destruct T1.
left; now left.
case Hh; intros P.
left; right.
split; try apply H0; assumption.
right.
split; try split; try apply H0; assumption.
clear T1.
423
(*  - end of assertions *)
424 425
destruct T.
(* normal case *)
BOLDO Sylvie's avatar
BOLDO Sylvie committed
426
apply round_N_eq_UP_pt with (pred_flt f)...
427
rewrite <- (round_DN_minus_eps_pos radix2 (FLT_exp emin prec) f) with (eps:=-h); try assumption.
428 429 430 431
replace (f--h) with (f+h) by ring.
apply round_DN_pt...
split.
auto with real.
432 433 434
apply Rle_trans with (1:=Hh).
apply Rle_trans with (/2*ulp_flt f).
apply Rmult_le_compat_r.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
435
apply ulp_ge_0.
436 437 438 439
fourier.
case H0.
intros Y; rewrite Y.
rewrite <- (Rmult_1_l (ulp_flt f)) at 2.
440
apply Rmult_le_compat_r.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
441
apply ulp_ge_0.
442
fourier.
443
intros Y; rewrite (proj1 Y); now right.
444
replace (f+h) with (pred_flt f + (f-pred_flt f+h)) by ring.
445
pattern f at 4; rewrite <- (round_UP_pred_plus_eps_pos radix2 (FLT_exp emin prec) f) with (eps:=(f - pred_flt f + h)); try assumption.
446 447 448 449
apply round_UP_pt...
replace (f-pred_flt f) with (ulp_flt (pred_flt f)).
split.
apply Rplus_lt_reg_l with (-h); ring_simplify.
450 451 452 453
case H0; [intros Y|intros (Y1,Y2)].
apply Rle_lt_trans with (1:=Hh).
rewrite Y.
rewrite <- (Rmult_1_l (ulp_flt f)) at 2.
454
apply Rmult_lt_compat_r.
455
rewrite ulp_neq_0;[apply bpow_gt_0|now apply Rgt_not_eq].
456
fourier.
457 458 459
apply Rlt_le_trans with (1:=Y2).
rewrite Y1.
apply Rmult_le_compat_r.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
460
apply ulp_ge_0.
461
fourier.
462 463
apply Rplus_le_reg_l with (-ulp_flt (pred_flt f)); ring_simplify.
now left.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
464
rewrite pred_eq_pos; try now left.
465
pattern f at 2; rewrite <- (pred_pos_plus_ulp radix2 (FLT_exp emin prec) f)...
466 467 468 469 470 471
ring.
apply Rplus_lt_reg_l with (-f); ring_simplify.
apply Rle_lt_trans with (-(/2 * ulp_flt (pred_flt f))).
right.
apply trans_eq with ((pred_flt f - f) / 2).
field.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
472
rewrite pred_eq_pos; try now left.
473
pattern f at 2; rewrite <- (pred_pos_plus_ulp radix2 (FLT_exp emin prec) f)...
474 475 476
field.
replace h with (--h) by ring.
apply Ropp_lt_contravar.
477 478 479 480
case H0;[intros Y|intros (Y1,Y2)].
apply Rle_lt_trans with (1:=Hh).
rewrite Y.
apply Rmult_lt_compat_r.
481
rewrite ulp_neq_0; try apply bpow_gt_0; now apply Rgt_not_eq.
482
fourier.
483 484 485 486 487
apply Rlt_le_trans with (1:=Y2).
rewrite Y1.
right; field.
(* complex case: even choosing *)
elim H0; intros  T1 (T2,T3); clear H0.
488
assert (pred_flt f = bpow (mag radix2 f - 1) - bpow (mag radix2 f - 1 -prec)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
489
rewrite pred_eq_pos; try now left.
490
unfold pred_pos; case Req_bool_spec.
491 492 493 494 495 496 497 498 499
intros _; rewrite <- T2.
apply f_equal, f_equal.
unfold FLT_exp.
rewrite Z.max_l.
ring.
omega.
intros Y; now contradict T2.
assert (round radix2 (FLT_exp emin prec) Zfloor (f+h) = pred_flt f).
replace (f+h) with (f-(-h)) by ring.
500
apply round_DN_minus_eps_pos...
501 502 503 504
split.
auto with real.
rewrite T3, T1.
apply Rmult_le_compat_r.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
505
apply ulp_ge_0.
506 507 508
fourier.
assert (round radix2 (FLT_exp emin prec) Zceil (f+h) = f).
replace (f+h) with (pred_flt f + /2*ulp_flt (pred_flt f)).
509
apply round_UP_pred_plus_eps_pos...
510 511 512
split.
apply Rmult_lt_0_compat.
fourier.
513
rewrite ulp_neq_0; try now apply Rgt_not_eq.
514 515 516
apply bpow_gt_0.
rewrite <- (Rmult_1_l (ulp_flt (pred_flt f))) at 2.
apply Rmult_le_compat_r.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
517
apply ulp_ge_0.
518 519 520
fourier.
rewrite T1, H0, <- T2.
replace h with (--h) by ring; rewrite T3.
521
replace (bpow (mag radix2 f - 1 - prec)) with (/2*ulp_flt f).
522 523 524 525 526 527 528 529
field.
replace (/2) with (bpow (-1)) by reflexivity.
rewrite T2 at 1.
rewrite ulp_bpow, <- bpow_plus.
apply f_equal; unfold FLT_exp.
rewrite Z.max_l.
ring.
omega.
530
assert ((Z.even (Zfloor (scaled_mantissa radix2 (FLT_exp emin prec) (f + h)))) = false).
531 532
replace (Zfloor (scaled_mantissa radix2 (FLT_exp emin prec) (f + h)))
   with (Zpower radix2 prec -1)%Z.
533 534 535
unfold Zminus; rewrite Z.even_add.
rewrite Z.even_opp.
rewrite Z.even_pow.
536 537
reflexivity.
unfold Prec_gt_0 in prec_gt_0_; omega.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
538
apply eq_IZR.
539 540 541 542 543
rewrite <- scaled_mantissa_DN...
2: rewrite H4; assumption.
rewrite H4.
unfold scaled_mantissa.
rewrite bpow_opp.
544
rewrite <- ulp_neq_0; try now apply Rgt_not_eq.
545 546 547
rewrite T1.
rewrite Rinv_mult_distr.
2: apply Rgt_not_eq; fourier.
548 549
2: apply Rgt_not_eq; rewrite ulp_neq_0; try apply bpow_gt_0.
2: now apply Rgt_not_eq.
550 551 552 553 554 555 556 557 558 559 560 561
rewrite Rinv_involutive.
2: apply Rgt_not_eq; fourier.
rewrite T2 at 2.
rewrite ulp_bpow.
rewrite <- bpow_opp.
unfold FLT_exp at 2.
rewrite Z.max_l.
2: omega.
replace 2 with (bpow 1) by reflexivity.
rewrite <- bpow_plus.
rewrite H0.
rewrite Rmult_minus_distr_r, <- 2!bpow_plus.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
562
rewrite minus_IZR.
563
apply f_equal2.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
564
rewrite IZR_Zpower.
565 566 567 568 569 570 571 572 573 574 575 576
apply f_equal.
ring.
unfold Prec_gt_0 in prec_gt_0_; omega.
apply trans_eq with (bpow 0).
reflexivity.
apply f_equal.
ring.
rewrite round_N_middle.
rewrite H5.
rewrite H6.
reflexivity.
rewrite H5, H4.
577
pattern f at 1; rewrite <- (pred_pos_plus_ulp radix2 (FLT_exp emin prec) f); try assumption.
578
ring_simplify.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
579
rewrite <- pred_eq_pos;[idtac|now left].
580 581 582 583
rewrite T1.
replace h with (--h) by ring.
rewrite T3.
field.
584 585 586
Qed.

Lemma round_plus_small_id: forall f h, format f -> (bpow (prec+emin) <= Rabs f)  
587
   -> Rabs h <= /4* ulp_flt f -> round_flt (f+h) = f.
588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608
intros f h Ff H1 H2.
case (Rle_or_lt 0 f); intros V.
case V; clear V; intros V.
apply round_plus_small_id_aux; try assumption.
rewrite Rabs_right in H1; try assumption.
apply Rle_ge; now left.
contradict H1.
rewrite <- V, Rabs_R0.
apply Rlt_not_le, bpow_gt_0.
rewrite <- (Ropp_involutive f), <- (Ropp_involutive h).
replace (--f + --h) with (-(-f+-h)) by ring.
rewrite round_NE_opp.
apply f_equal.
apply round_plus_small_id_aux.
now apply generic_format_opp.
rewrite Rabs_left in H1; try assumption.
auto with real.
now rewrite Rabs_Ropp, ulp_opp.
Qed.


609

BOLDO Sylvie's avatar
BOLDO Sylvie committed
610
Definition avg_naive (x y : R) :=round_flt(round_flt(x+y)/2).
611 612 613 614 615 616

Variables x y:R.
Hypothesis Fx: format x.
Hypothesis Fy: format y.

Let a:=(x+y)/2.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
617
Let av:=avg_naive x y.
618

BOLDO Sylvie's avatar
BOLDO Sylvie committed
619
Lemma avg_naive_correct: av = round_flt a.
620 621 622 623
Proof with auto with typeclass_instances.
case (Rle_or_lt (bpow (prec + emin)) (Rabs (x+y))).
(* normal case: division by 2 is exact *)
intros H.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
624
unfold av,a,avg_naive.
625 626 627 628 629 630 631 632 633
rewrite round_generic...
now apply sym_eq, FLT_round_half.
apply FLT_format_half.
apply generic_format_round...
apply abs_round_ge_generic...
apply FLT_format_bpow...
unfold Prec_gt_0 in prec_gt_0_; omega.
(* subnormal case: addition is exact, but division by 2 is not *)
intros H.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
634
unfold av, avg_naive.
635 636 637 638 639 640 641 642
replace (round_flt (x + y)) with (x+y).
reflexivity.
apply sym_eq, round_generic...
apply FLT_format_plus_small...
left; assumption.
Qed.


643

BOLDO Sylvie's avatar
BOLDO Sylvie committed
644
Lemma avg_naive_symmetry: forall u v, avg_naive u v = avg_naive v u.
645
Proof.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
646
intros u v; unfold avg_naive.
647 648 649
rewrite Rplus_comm; reflexivity.
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
650
Lemma avg_naive_symmetry_Ropp: forall u v, avg_naive (-u) (-v) = - avg_naive u v.
651
Proof.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
652
intros u v; unfold avg_naive.
653 654 655 656 657 658
replace (-u+-v) with (-(u+v)) by ring.
rewrite round_NE_opp.
replace (- round_flt (u + v) / 2) with (- (round_flt (u + v) / 2)) by (unfold Rdiv; ring).
now rewrite round_NE_opp.
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
659
Lemma avg_naive_same_sign_1: 0 <= a -> 0 <= av.
660 661
Proof with auto with typeclass_instances.
intros H.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
662
rewrite avg_naive_correct.
663 664 665 666
apply round_ge_generic...
apply generic_format_0.
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
667
Lemma avg_naive_same_sign_2: a <= 0-> av <= 0.
668 669
Proof with auto with typeclass_instances.
intros H.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
670
rewrite avg_naive_correct.
671 672 673 674
apply round_le_generic...
apply generic_format_0.
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
675
Lemma avg_naive_between: Rmin x y <= av <= Rmax x y.
676
Proof with auto with typeclass_instances.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
677
rewrite avg_naive_correct.
678
split.
679
apply round_ge_generic...
680
now apply P_Rmin.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
681 682 683
apply Rmult_le_reg_l with (1 := Rlt_0_2).
replace (2 * Rmin x y) with (Rmin x y + Rmin x y) by ring.
replace (2 * a) with (x + y) by (unfold a; field).
684 685 686 687
apply Rplus_le_compat.
apply Rmin_l.
apply Rmin_r.
(* *)
688
apply round_le_generic...
689
now apply Rmax_case.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
690 691 692
apply Rmult_le_reg_l with (1 := Rlt_0_2).
replace (2 * a) with (x + y) by (unfold a; field).
replace (2 * Rmax x y) with (Rmax x y + Rmax x y) by ring.
693 694 695
apply Rplus_le_compat.
apply Rmax_l.
apply Rmax_r.
696 697 698
Qed.


BOLDO Sylvie's avatar
BOLDO Sylvie committed
699
Lemma avg_naive_zero: a = 0 -> av = 0.
700
Proof with auto with typeclass_instances.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
701
intros H1; rewrite avg_naive_correct, H1.
702 703 704 705 706
rewrite round_0...
Qed.



BOLDO Sylvie's avatar
BOLDO Sylvie committed
707
Lemma avg_naive_no_underflow: (bpow emin) <= Rabs a -> av <> 0.
708 709 710 711 712 713 714 715 716 717
Proof with auto with typeclass_instances.
intros H.
(* *)
cut (bpow emin <= Rabs av).
intros H1 H2.
rewrite H2 in H1; rewrite Rabs_R0 in H1.
contradict H1.
apply Rlt_not_le.
apply bpow_gt_0.
(* *)
BOLDO Sylvie's avatar
BOLDO Sylvie committed
718
rewrite avg_naive_correct.
719 720 721 722 723 724
apply abs_round_ge_generic...
apply FLT_format_bpow...
omega.
Qed.


BOLDO Sylvie's avatar
BOLDO Sylvie committed
725
Lemma avg_naive_correct_weak1: Rabs (av -a) <= /2*ulp_flt a.
726
Proof with auto with typeclass_instances.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
727
rewrite avg_naive_correct.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
728
apply error_le_half_ulp...
729 730
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
731
Lemma avg_naive_correct_weak2: Rabs (av -a) <= 3/2*ulp_flt a.
732
Proof with auto with typeclass_instances.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
733
apply Rle_trans with (1:=avg_naive_correct_weak1).
734
apply Rmult_le_compat_r.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
735
unfold ulp; apply ulp_ge_0.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
736
lra.
737 738 739 740 741 742 743 744
Qed.



(* Hypothesis diff_sign: (0 <= x /\ y <= 0) \/ (x <= 0 /\ 0 <= y).
  is useless for properties: only useful for preventing overflow *)


745 746 747



BOLDO Sylvie's avatar
BOLDO Sylvie committed
748
Definition avg_sum_half (x y : R) :=round_flt(round_flt(x/2) + round_flt(y/2)).
749

BOLDO Sylvie's avatar
BOLDO Sylvie committed
750
Let av2:=avg_sum_half x y.
751 752


BOLDO Sylvie's avatar
BOLDO Sylvie committed
753
Lemma avg_sum_half_correct: bpow (emin +prec+prec+1) <= Rabs x -> av2 = round_flt a.
754 755 756 757
Proof with auto with typeclass_instances.
intros Hx.
assert (G:(0 < prec)%Z).
unfold Prec_gt_0 in prec_gt_0_; assumption.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
758
unfold av2, avg_sum_half.
759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789
replace (round_flt (x/2)) with (x/2).
2: apply sym_eq, round_generic...
2: apply FLT_format_half; try assumption.
2: apply Rle_trans with (2:=Hx).
2: apply bpow_le; omega.
case (Rle_or_lt (bpow (prec + emin)) (Rabs y)).
(* y is big enough so that y/2 is correct *)
intros Hy.
replace (round_flt (y/2)) with (y/2).
apply f_equal; unfold a; field.
apply sym_eq, round_generic...
apply FLT_format_half; assumption.
(* y is a subnormal, then it is too small to impact the result *)
intros Hy.
assert (format (x/2)).
apply FLT_format_half.
assumption.
apply Rle_trans with (2:=Hx).
apply bpow_le.
omega.
assert (bpow (prec+emin) <= Rabs (x/2)).
apply Rmult_le_reg_l with (bpow 1).
apply bpow_gt_0.
rewrite <- bpow_plus.
apply Rle_trans with (Rabs x).
apply Rle_trans with (2:=Hx).
apply bpow_le.
omega.
rewrite <- (Rabs_right (bpow 1)).
rewrite <- Rabs_mult.
right; apply f_equal.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
790 791
change (bpow 1) with 2.
field.
792
apply Rle_ge, bpow_ge_0.
793
assert (K1: Rabs (y / 2) <= bpow (prec+emin-1)).
794 795
unfold Rdiv; rewrite Rabs_mult.
unfold Zminus; rewrite bpow_plus.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
796 797
simpl; rewrite (Rabs_pos_eq (/2)).
apply (Rmult_le_compat_r (/2)).
798 799 800
fourier.
now left.
fourier.
801
assert (K2:bpow (prec+emin-1) <= / 4 * ulp_flt (x / 2)).
802 803 804 805 806
assert (Z: x/2 <> 0).
intros K; contradict H0.
rewrite K, Rabs_R0.
apply Rlt_not_le, bpow_gt_0.
rewrite ulp_neq_0; trivial.
807 808
replace (/4) with (bpow (-2)) by reflexivity.
rewrite <- bpow_plus.
809
apply bpow_le.
810
unfold cexp, FLT_exp.
811 812
assert (emin+prec+prec+1 -1 < mag radix2 (x/2))%Z.
destruct (mag radix2 (x/2)) as (e,He).
813 814 815 816 817 818 819 820 821 822
simpl.
apply lt_bpow with radix2.
apply Rle_lt_trans with (Rabs (x/2)).
unfold Rdiv; rewrite Rabs_mult.
unfold Zminus; rewrite bpow_plus.
simpl; rewrite (Rabs_right (/2)).
apply Rmult_le_compat_r.
fourier.
exact Hx.
fourier.
823
apply He; trivial.
824 825 826
rewrite Z.max_l.
omega.
omega.
827 828 829 830 831 832 833
(* . *)
apply trans_eq with (x/2).
apply round_plus_small_id; try assumption.
apply Rle_trans with (2:=K2).
apply abs_round_le_generic...
apply FLT_format_bpow...
omega.
834 835 836
unfold a; apply sym_eq.
replace ((x+y)/2) with (x/2+y/2) by field.
apply round_plus_small_id; try assumption.
837
now apply Rle_trans with (2:=K2).
838 839 840 841
Qed.



842 843 844 845 846 847 848 849 850 851 852 853
End av1.

Section av3.

Notation bpow e := (bpow radix2 e).

Variable emin prec : Z.
Context { prec_gt_0_ : Prec_gt_0 prec }.

Notation format := (generic_format radix2 (FLT_exp emin prec)).
Notation round_flt :=(round radix2 (FLT_exp emin prec) ZnearestE). 
Notation ulp_flt :=(ulp radix2 (FLT_exp emin prec)).
854
Notation cexp := (cexp radix2 (FLT_exp emin prec)).
855

BOLDO Sylvie's avatar
BOLDO Sylvie committed
856
Definition avg_half_sub (x y : R) :=round_flt(x+round_flt(round_flt(y-x)/2)).
857 858 859 860 861 862

Variables x y:R.
Hypothesis Fx: format x.
Hypothesis Fy: format y.

Let a:=(x+y)/2.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
863
Let av:=avg_half_sub x y.
864 865


BOLDO Sylvie's avatar
BOLDO Sylvie committed
866 867
Lemma avg_half_sub_symmetry_Ropp: forall u v, avg_half_sub (-u) (-v) = - avg_half_sub u v.
intros u v; unfold avg_half_sub.
868 869 870 871 872 873 874 875 876 877
replace (-v--u) with (-(v-u)) by ring.
rewrite round_NE_opp.
replace (- round_flt (v-u) / 2) with (- (round_flt (v-u) / 2)) by (unfold Rdiv; ring).
rewrite round_NE_opp.
replace (- u + - round_flt (round_flt (v - u) / 2)) with
   (-(u+round_flt (round_flt (v - u) / 2))) by ring.
apply round_NE_opp.
Qed.


BOLDO Sylvie's avatar
BOLDO Sylvie committed
878
Lemma avg_half_sub_same_sign_1: 0 <= a -> 0 <= av.
879 880 881 882 883 884 885 886
Proof with auto with typeclass_instances.
intros H.
apply round_ge_generic...
apply generic_format_0.
apply Rplus_le_reg_l with (-x).
ring_simplify.
apply round_ge_generic...
now apply generic_format_opp.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
887
apply Rmult_le_reg_l with (1 := Rlt_0_2).
888 889 890 891 892 893 894 895 896
apply Rle_trans with (-(2*x)).
right; ring.
apply Rle_trans with (round_flt (y - x)).
2: right; field.
apply round_ge_generic...
apply generic_format_opp.
now apply FLT_format_double...
apply Rplus_le_reg_l with (2*x).
apply Rmult_le_reg_r with (/2).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
897
lra.
898 899 900 901 902
apply Rle_trans with 0;[right; ring|idtac].
apply Rle_trans with (1:=H).
right; unfold a, Rdiv; ring.
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
903
Lemma avg_half_sub_same_sign_2: a <= 0-> av <= 0.
904 905 906 907 908 909 910 911
Proof with auto with typeclass_instances.
intros H.
apply round_le_generic...
apply generic_format_0.
apply Rplus_le_reg_l with (-x).
ring_simplify.
apply round_le_generic...
now apply generic_format_opp.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
912
apply Rmult_le_reg_l with (1 := Rlt_0_2).
913 914 915 916 917 918 919 920 921
apply Rle_trans with (-(2*x)).
2: right; ring.
apply Rle_trans with (round_flt (y - x)).
right; field.
apply round_le_generic...
apply generic_format_opp.
now apply FLT_format_double...
apply Rplus_le_reg_l with (2*x).
apply Rmult_le_reg_r with (/2).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
922
lra.
923 924 925 926 927 928 929 930
apply Rle_trans with 0;[idtac|right; ring].
apply Rle_trans with (2:=H).
right; unfold a, Rdiv; ring.
Qed.




BOLDO Sylvie's avatar
BOLDO Sylvie committed
931 932
Lemma avg_half_sub_between_aux: forall u v, format u -> format v -> u <= v ->
    u <= avg_half_sub u v <= v.
933 934 935 936 937 938 939 940 941 942 943 944 945 946 947
Proof with auto with typeclass_instances.
clear Fx Fy a av x y.
intros x y Fx Fy M.
split.
(* . *)
apply round_ge_generic...
apply Rplus_le_reg_l with (-x).
ring_simplify.
apply round_ge_generic...
apply generic_format_0.
unfold Rdiv; apply Rmult_le_pos.
apply round_ge_generic...
apply generic_format_0.
apply Rplus_le_reg_l with x.
now ring_simplify.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
948
lra.
949 950 951 952 953 954 955 956 957 958 959 960 961 962
(* . *)
apply round_le_generic...
assert (H:(0 <= round radix2 (FLT_exp emin prec) Zfloor (y-x))).
apply round_ge_generic...
apply generic_format_0.
apply Rplus_le_reg_l with x.
now ring_simplify.
destruct H as [H|H].
(* .. *)
pattern y at 2; replace y with (x + (y-x)) by ring.
apply Rplus_le_compat_l.
case (generic_format_EM radix2 (FLT_exp emin prec) (y-x)); intros K.
apply round_le_generic...
rewrite round_generic...
Guillaume Melquiond's avatar
Guillaume Melquiond committed
963
apply Rmult_le_reg_l with (1 := Rlt_0_2).
964 965 966 967 968 969 970 971
apply Rplus_le_reg_l with (2*x-y).
apply Rle_trans with x.
right; field.
apply Rle_trans with (1:=M).
right; field.
apply Rle_trans with (round radix2 (FLT_exp emin prec) Zfloor (y - x)).
apply round_le_generic...
apply generic_format_round...
Guillaume Melquiond's avatar
Guillaume Melquiond committed
972
apply Rmult_le_reg_l with (1 := Rlt_0_2).
973 974 975 976 977 978 979
apply Rle_trans with (round_flt (y - x)).
right; field.
case (round_DN_or_UP radix2 (FLT_exp emin prec) ZnearestE (y-x));
   intros H1; rewrite H1.
apply Rplus_le_reg_l with (-round radix2 (FLT_exp emin prec) Zfloor (y - x)).
ring_simplify.
now left.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
980
rewrite round_UP_DN_ulp.
981 982
apply Rplus_le_reg_l with (-round radix2 (FLT_exp emin prec) Zfloor (y - x)); ring_simplify.
apply round_DN_pt...
983
apply generic_format_ulp...
984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010
case (Rle_or_lt (bpow (emin + prec - 1))  (y-x)); intros P.
apply FLT_ulp_le_id...
apply Rle_trans with (2:=P).
apply bpow_le; unfold Prec_gt_0 in prec_gt_0_; omega.
contradict K.
apply FLT_format_plus_small...
now apply generic_format_opp.
rewrite Rabs_right.
apply Rle_trans with (bpow (emin+prec-1)).
left; exact P.
apply bpow_le; omega.
apply Rle_ge; apply Rplus_le_reg_l with x; now ring_simplify.
assumption.
apply round_DN_pt...
(* .. *)
case M; intros H1.
2: rewrite H1; replace (y-y) with 0 by ring.
2: rewrite round_0...
2: unfold Rdiv; rewrite Rmult_0_l.
2: rewrite round_0...
2: right; ring.
apply Rle_trans with (x+0).
2: rewrite Rplus_0_r; assumption.
apply Rplus_le_compat_l.
replace 0 with (round_flt (bpow emin/2)).
apply round_le...
unfold Rdiv; apply Rmult_le_compat_r.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1011
lra.
1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025
apply round_le_generic...
apply FLT_format_bpow...
omega.
case (Rle_or_lt (y-x) (bpow emin)); trivial.
intros H2.
contradict H.
apply Rlt_not_eq.
apply Rlt_le_trans with (bpow emin).
apply bpow_gt_0.
apply round_DN_pt...
apply FLT_format_bpow...
omega.
now left.
replace (bpow emin /2) with (bpow (emin-1)).
1026
unfold round, scaled_mantissa, cexp, FLT_exp.
1027
rewrite mag_bpow.
1028 1029 1030 1031 1032 1033 1034
replace (emin - 1 + 1 - prec)%Z with (emin-prec)%Z by ring.
rewrite Z.max_r.
2: unfold Prec_gt_0 in prec_gt_0_; omega.
rewrite <- bpow_plus.
replace (emin-1+-emin)%Z with (-1)%Z by ring.
replace (ZnearestE (bpow (-1))) with 0%Z.
unfold F2R; simpl; ring.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1035
change (bpow (-1)) with (/2).
1036 1037 1038 1039 1040 1041
simpl; unfold Znearest.
replace (Zfloor (/2)) with 0%Z.
rewrite Rcompare_Eq.
reflexivity.
simpl; ring.
apply sym_eq, Zfloor_imp.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1042
simpl ; lra.
1043 1044 1045 1046
unfold Zminus; rewrite bpow_plus.
reflexivity.
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
1047
Lemma avg_half_sub_between: Rmin x y <= av <= Rmax x y.
1048 1049 1050 1051 1052
Proof with auto with typeclass_instances.
case (Rle_or_lt x y); intros M.
(* x <= y *)
rewrite Rmin_left; try exact M.
rewrite Rmax_right; try exact M.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1053
now apply avg_half_sub_between_aux.
1054 1055 1056 1057
(* y < x *)
rewrite Rmin_right; try now left.
rewrite Rmax_left; try now left.
unfold av; rewrite <- (Ropp_involutive x); rewrite <- (Ropp_involutive y).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1058
rewrite avg_half_sub_symmetry_Ropp.
1059
split; apply Ropp_le_contravar.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1060
apply avg_half_sub_between_aux.
1061 1062 1063
now apply generic_format_opp.
now apply generic_format_opp.
apply Ropp_le_contravar; now left.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1064
apply avg_half_sub_between_aux.
1065 1066 1067 1068 1069 1070
now apply generic_format_opp.
now apply generic_format_opp.
apply Ropp_le_contravar; now left.
Qed.


BOLDO Sylvie's avatar
BOLDO Sylvie committed
1071
Lemma avg_half_sub_zero: a = 0 -> av = 0.
1072 1073 1074 1075 1076 1077 1078 1079
Proof with auto with typeclass_instances.
intros H.
assert (y=-x).
apply Rplus_eq_reg_l with x.
apply Rmult_eq_reg_r with (/2).
apply trans_eq with a.
reflexivity.
rewrite H; ring.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1080
lra.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1081
unfold av, avg_half_sub.
1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094
rewrite H0.
replace (-x-x) with (-(2*x)) by ring.
rewrite round_generic with (x:=(-(2*x)))...
replace (-(2*x)/2) with (-x) by field.
rewrite round_generic with (x:=-x)...
replace (x+-x) with 0 by ring.
apply round_0...
now apply generic_format_opp.
apply generic_format_opp.
now apply FLT_format_double.
Qed.


BOLDO Sylvie's avatar
BOLDO Sylvie committed
1095
Lemma avg_half_sub_no_underflow_aux_aux: forall z:Z, (0 < z)%Z -> 
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1096
    (ZnearestE (IZR z / 2) < z)%Z.
1097 1098 1099
Proof.
intros z H1.
case (Zle_lt_or_eq 1 z); [omega|intros H2|intros H2].
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1100 1101 1102
apply lt_IZR.
apply Rplus_lt_reg_r with (- ((IZR z)/2)).
apply Rle_lt_trans with (-(((IZR z) /2) - IZR (ZnearestE (IZR z / 2)))).
1103 1104 1105
right; ring.
apply Rle_lt_trans with (1:= RRle_abs _).
rewrite Rabs_Ropp.
1106
apply Rle_lt_trans with (1:=Znearest_half (fun x => negb (Z.even x)) _).
1107
apply Rle_lt_trans with (1*/2);[right; ring|idtac].
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1108
apply Rlt_le_trans with ((IZR z)*/2);[idtac|right; field].
1109
apply Rmult_lt_compat_r.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1110
lra.
1111
now apply IZR_lt.
1112 1113 1114 1115 1116 1117 1118 1119
rewrite <- H2.
unfold Znearest; simpl.
replace (Zfloor (1 / 2)) with 0%Z.
rewrite Rcompare_Eq.
simpl; omega.
simpl; field.
unfold Rdiv; rewrite Rmult_1_l.
apply sym_eq, Zfloor_imp.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1120
simpl; lra.
1121 1122 1123
Qed.


BOLDO Sylvie's avatar
BOLDO Sylvie committed
1124
Lemma avg_half_sub_no_underflow_aux1: forall f, format f -> 0 < f ->
1125 1126 1127 1128
  f <= round_flt (f/2) -> False.
Proof with auto with typeclass_instances.
intros f Ff Hf1 Hf2.
apply FLT_format_generic in Ff...
1129
destruct Ff as [g H1 H2 H3].
1130 1131 1132 1133 1134 1135 1136 1137
case (Zle_lt_or_eq emin (Fexp g)); try exact H3; intros H4.
contradict Hf2.
apply Rlt_not_le.
rewrite round_generic...
apply Rplus_lt_reg_l with (-(f/2)).
apply Rle_lt_trans with 0;[right; ring|idtac].
apply Rlt_le_trans with (f*/2);[idtac|right;field].
apply Rmult_lt_0_compat; try assumption.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1138
lra.
1139 1140 1141 1142
apply generic_format_FLT.
exists (Float radix2 (Fnum g) (Fexp g-1)).
rewrite H1; unfold F2R; simpl.
unfold Zminus; rewrite bpow_plus.
1143 1144
apply Rmult_assoc.
easy.
1145 1146 1147 1148 1149 1150 1151 1152
simpl; omega.
contradict Hf2; apply Rlt_not_le.
unfold round, scaled_mantissa.
replace (cexp (f/2)) with emin.
rewrite H1; unfold F2R; simpl.
rewrite <- H4.
apply Rmult_lt_compat_r.
apply bpow_gt_0.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1153 1154
apply IZR_lt.
replace (IZR (Fnum g) * bpow emin / 2 * bpow (- emin)) with (IZR (Fnum g) /2).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1155
apply avg_half_sub_no_underflow_aux_aux.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1156
apply lt_IZR.
1157 1158 1159 1160 1161
apply Rmult_lt_reg_r with (bpow (Fexp g)).
apply bpow_gt_0.
rewrite Rmult_0_l.
apply Rlt_le_trans with (1:=Hf1).
right; rewrite H1; reflexivity.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1162
unfold Rdiv; apply trans_eq with (