Average.v 53.5 KB
Newer Older
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1
Require Import Reals Fourier Psatz.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
2
From Flocq Require Import Core Plus_error.
3 4 5 6 7

Open Scope R_scope.

Section av1.

8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

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.

36 37 38 39 40 41 42 43 44 45

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)).
Notation round_flt :=(round radix2 (FLT_exp emin prec) ZnearestE). 
Notation ulp_flt :=(ulp radix2 (FLT_exp emin prec)).
46
Notation cexp := (cexp radix2 (FLT_exp emin prec)).
47
Notation pred_flt := (pred radix2 (FLT_exp emin prec)).
48 49 50 51 52 53

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...
54
destruct Fu as [uf H1 H2 H3].
55 56 57 58
exists (Float radix2 (Fnum uf) (Fexp uf+1)).
rewrite H1; unfold F2R; simpl.
rewrite bpow_plus, bpow_1.
simpl;ring.
59 60 61
easy.
apply Zle_trans with (1:=H3).
apply Zle_succ.
62 63 64 65 66 67 68
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...
69
destruct Fu as [[n e] H1 H2 H3].
70 71 72 73 74
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
75 76
change (bpow (-(1))) with (/2).
unfold Rdiv; ring.
77 78 79
easy.
cut (prec + emin < prec +e)%Z.
  simpl ; omega.
80 81 82 83 84 85 86 87
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
88 89 90
rewrite <- abs_IZR.
rewrite <- IZR_Zpower.
now apply IZR_lt.
91
now apply Zlt_le_weak.
92 93
Qed.

94 95 96 97 98 99 100 101 102 103 104 105 106
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.
107
assert (prec+emin < mag radix2 z)%Z.
108
apply lt_bpow with radix2.
109
destruct mag as (e,He); simpl.
110 111
apply Rle_lt_trans with (1:=Hz).
now apply He.
112
unfold cexp, FLT_exp.
113
replace ((mag radix2 (z/2))-prec)%Z with ((mag radix2 z -1) -prec)%Z.
114 115 116
rewrite Z.max_l; try omega.
rewrite Z.max_l; try omega.
apply Zplus_eq_compat; try reflexivity.
117 118
apply sym_eq, mag_unique.
destruct (mag radix2 z) as (e,He); simpl.
119 120 121 122 123 124 125 126 127
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
128
change (bpow 1) with 2%R.
129 130 131 132 133 134
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
135 136
change (bpow 1) with 2.
right; field.
137
now apply He.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
138
lra.
139 140 141 142 143 144 145
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
146 147
change (bpow 1) with 2.
field.
148
unfold Zminus; rewrite bpow_plus.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
149 150
change (bpow (-(1))) with (/2).
field.
151
Qed.
152 153 154 155

Lemma FLT_ulp_le_id: forall u, bpow emin <= u -> ulp_flt u <= u.
Proof with auto with typeclass_instances.
intros u H.
156 157
rewrite ulp_neq_0.
2: apply Rgt_not_eq, Rlt_le_trans with (2:=H), bpow_gt_0.
158
case (Rle_or_lt (bpow (emin+prec-1)) u); intros Hu.
159 160
unfold ulp; rewrite cexp_FLT_FLX.
unfold cexp, FLX_exp.
161
destruct (mag radix2 u) as (e,He); simpl.
162 163 164 165 166 167 168 169 170 171 172 173
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.
174
unfold ulp; rewrite cexp_FLT_FIX.
175 176 177 178 179 180 181 182 183 184 185 186 187
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
188
Proof.
189
intros u.
190 191 192 193
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
194
apply ulp_ge_0.
195 196
left; apply Rlt_plus_1.
rewrite 2!ulp_neq_0; trivial.
197
2: lra.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
198
change 2 at 2 with (bpow 1).
199
rewrite <- bpow_plus.
200 201 202
apply bpow_le.
case (Rle_or_lt (bpow (emin+prec-1)) (Rabs u)); intros Hu.
(* *)
203 204 205
rewrite cexp_FLT_FLX.
rewrite cexp_FLT_FLX; trivial.
unfold cexp, FLX_exp.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
206
change 2 with (bpow 1).
207
rewrite Rmult_comm, mag_mult_bpow.
208 209 210 211 212 213 214 215 216
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
217
rewrite <- (abs_IZR 2).
218
now apply IZR_le.
219
(* *)
220 221 222
rewrite cexp_FLT_FIX.
rewrite cexp_FLT_FIX; trivial.
unfold FIX_exp, cexp; omega.
223 224
apply Rlt_le_trans with (1:=Hu).
apply bpow_le; omega.
225
lra.
226
rewrite Rabs_mult.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
227
rewrite Rabs_pos_eq.
228
2: now apply IZR_le.
229
apply Rlt_le_trans with (2*bpow (emin + prec - 1)).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
230
apply Rmult_lt_compat_l with (1 := Rlt_0_2).
231
assumption.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
232
change 2 with (bpow 1).
233 234 235 236 237
rewrite <- bpow_plus.
apply bpow_le; omega.
Qed.


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

Lemma round_plus_small_id: forall f h, format f -> (bpow (prec+emin) <= Rabs f)  
570
   -> Rabs h <= /4* ulp_flt f -> round_flt (f+h) = f.
571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591
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.


592

BOLDO Sylvie's avatar
BOLDO Sylvie committed
593
Definition avg_naive (x y : R) :=round_flt(round_flt(x+y)/2).
594 595 596 597 598 599

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

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

BOLDO Sylvie's avatar
BOLDO Sylvie committed
602
Lemma avg_naive_correct: av = round_flt a.
603 604 605 606
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
607
unfold av,a,avg_naive.
608 609 610 611 612 613 614 615 616
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
617
unfold av, avg_naive.
618 619 620 621 622 623 624 625
replace (round_flt (x + y)) with (x+y).
reflexivity.
apply sym_eq, round_generic...
apply FLT_format_plus_small...
left; assumption.
Qed.


626

BOLDO Sylvie's avatar
BOLDO Sylvie committed
627
Lemma avg_naive_symmetry: forall u v, avg_naive u v = avg_naive v u.
628
Proof.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
629
intros u v; unfold avg_naive.
630 631 632
rewrite Rplus_comm; reflexivity.
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
633
Lemma avg_naive_symmetry_Ropp: forall u v, avg_naive (-u) (-v) = - avg_naive u v.
634
Proof.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
635
intros u v; unfold avg_naive.
636 637 638 639 640 641
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
642
Lemma avg_naive_same_sign_1: 0 <= a -> 0 <= av.
643 644
Proof with auto with typeclass_instances.
intros H.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
645
rewrite avg_naive_correct.
646 647 648 649
apply round_ge_generic...
apply generic_format_0.
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
650
Lemma avg_naive_same_sign_2: a <= 0-> av <= 0.
651 652
Proof with auto with typeclass_instances.
intros H.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
653
rewrite avg_naive_correct.
654 655 656 657
apply round_le_generic...
apply generic_format_0.
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
658
Lemma avg_naive_between: Rmin x y <= av <= Rmax x y.
659
Proof with auto with typeclass_instances.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
660
rewrite avg_naive_correct.
661
split.
662
apply round_ge_generic...
663
now apply P_Rmin.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
664 665 666
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).
667 668 669 670
apply Rplus_le_compat.
apply Rmin_l.
apply Rmin_r.
(* *)
671
apply round_le_generic...
672
now apply Rmax_case.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
673 674 675
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.
676 677 678
apply Rplus_le_compat.
apply Rmax_l.
apply Rmax_r.
679 680 681
Qed.


BOLDO Sylvie's avatar
BOLDO Sylvie committed
682
Lemma avg_naive_zero: a = 0 -> av = 0.
683
Proof with auto with typeclass_instances.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
684
intros H1; rewrite avg_naive_correct, H1.
685 686 687 688 689
rewrite round_0...
Qed.



BOLDO Sylvie's avatar
BOLDO Sylvie committed
690
Lemma avg_naive_no_underflow: (bpow emin) <= Rabs a -> av <> 0.
691 692 693 694 695 696 697 698 699 700
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
701
rewrite avg_naive_correct.
702 703 704 705 706 707
apply abs_round_ge_generic...
apply FLT_format_bpow...
omega.
Qed.


BOLDO Sylvie's avatar
BOLDO Sylvie committed
708
Lemma avg_naive_correct_weak1: Rabs (av -a) <= /2*ulp_flt a.
709
Proof with auto with typeclass_instances.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
710
rewrite avg_naive_correct.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
711
apply error_le_half_ulp...
712 713
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
714
Lemma avg_naive_correct_weak2: Rabs (av -a) <= 3/2*ulp_flt a.
715
Proof with auto with typeclass_instances.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
716
apply Rle_trans with (1:=avg_naive_correct_weak1).
717
apply Rmult_le_compat_r.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
718
unfold ulp; apply ulp_ge_0.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
719
lra.
720 721 722 723 724 725 726 727
Qed.



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


728 729 730



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

BOLDO Sylvie's avatar
BOLDO Sylvie committed
733
Let av2:=avg_sum_half x y.
734 735


BOLDO Sylvie's avatar
BOLDO Sylvie committed
736
Lemma avg_sum_half_correct: bpow (emin +prec+prec+1) <= Rabs x -> av2 = round_flt a.
737 738 739 740
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
741
unfold av2, avg_sum_half.
742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772
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
773 774
change (bpow 1) with 2.
field.
775
apply Rle_ge, bpow_ge_0.
776
assert (K1: Rabs (y / 2) <= bpow (prec+emin-1)).
777 778
unfold Rdiv; rewrite Rabs_mult.
unfold Zminus; rewrite bpow_plus.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
779 780
simpl; rewrite (Rabs_pos_eq (/2)).
apply (Rmult_le_compat_r (/2)).
781 782 783
fourier.
now left.
fourier.
784
assert (K2:bpow (prec+emin-1) <= / 4 * ulp_flt (x / 2)).
785 786 787 788 789
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.
790 791
replace (/4) with (bpow (-2)) by reflexivity.
rewrite <- bpow_plus.
792
apply bpow_le.
793
unfold cexp, FLT_exp.
794 795
assert (emin+prec+prec+1 -1 < mag radix2 (x/2))%Z.
destruct (mag radix2 (x/2)) as (e,He).
796 797 798 799 800 801 802 803 804 805
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.
806
apply He; trivial.
807 808 809
rewrite Z.max_l.
omega.
omega.
810 811 812 813 814 815 816
(* . *)
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.
817 818 819
unfold a; apply sym_eq.
replace ((x+y)/2) with (x/2+y/2) by field.
apply round_plus_small_id; try assumption.
820
now apply Rle_trans with (2:=K2).
821 822 823 824
Qed.



825 826 827 828 829 830 831 832 833 834 835 836
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)).
837
Notation cexp := (cexp radix2 (FLT_exp emin prec)).
838

BOLDO Sylvie's avatar
BOLDO Sylvie committed
839
Definition avg_half_sub (x y : R) :=round_flt(x+round_flt(round_flt(y-x)/2)).
840 841 842 843 844 845

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

Let a:=(x+y)/2.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
846
Let av:=avg_half_sub x y.
847 848


BOLDO Sylvie's avatar
BOLDO Sylvie committed
849 850
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.
851 852 853 854 855 856 857 858 859 860
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
861
Lemma avg_half_sub_same_sign_1: 0 <= a -> 0 <= av.
862 863 864 865 866 867 868 869
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
870
apply Rmult_le_reg_l with (1 := Rlt_0_2).
871 872 873 874 875 876 877 878 879
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
880
lra.
881 882 883 884 885
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
886
Lemma avg_half_sub_same_sign_2: a <= 0-> av <= 0.
887 888 889 890 891 892 893 894
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
895
apply Rmult_le_reg_l with (1 := Rlt_0_2).
896 897 898 899 900 901 902 903 904
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
905
lra.
906 907 908 909 910 911 912 913
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
914 915
Lemma avg_half_sub_between_aux: forall u v, format u -> format v -> u <= v ->
    u <= avg_half_sub u v <= v.
916 917 918 919 920 921 922 923 924 925 926 927 928 929 930
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
931
lra.
932 933 934 935 936 937 938 939 940 941 942 943 944 945
(* . *)
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
946
apply Rmult_le_reg_l with (1 := Rlt_0_2).
947 948 949 950 951 952 953 954
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
955
apply Rmult_le_reg_l with (1 := Rlt_0_2).
956 957 958 959 960 961 962
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
963
rewrite round_UP_DN_ulp.
964 965
apply Rplus_le_reg_l with (-round radix2 (FLT_exp emin prec) Zfloor (y - x)); ring_simplify.
apply round_DN_pt...
966
apply generic_format_ulp...
967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993
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
994
lra.
995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008
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)).
1009
unfold round, scaled_mantissa, cexp, FLT_exp.
1010
rewrite mag_bpow.
1011 1012 1013 1014 1015 1016 1017
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
1018
change (bpow (-1)) with (/2).
1019 1020 1021 1022 1023 1024
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
1025
simpl ; lra.
1026 1027 1028 1029
unfold Zminus; rewrite bpow_plus.
reflexivity.
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
1030
Lemma avg_half_sub_between: Rmin x y <= av <= Rmax x y.
1031 1032 1033 1034 1035
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
1036
now apply avg_half_sub_between_aux.
1037 1038 1039 1040
(* 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
1041
rewrite avg_half_sub_symmetry_Ropp.
1042
split; apply Ropp_le_contravar.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1043
apply avg_half_sub_between_aux.
1044 1045 1046
now apply generic_format_opp.
now apply generic_format_opp.
apply Ropp_le_contravar; now left.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1047
apply avg_half_sub_between_aux.
1048 1049 1050 1051 1052 1053
now apply generic_format_opp.
now apply generic_format_opp.
apply Ropp_le_contravar; now left.
Qed.


BOLDO Sylvie's avatar
BOLDO Sylvie committed
1054
Lemma avg_half_sub_zero: a = 0 -> av = 0.
1055 1056 1057 1058 1059 1060 1061 1062
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
1063
lra.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1064
unfold av, avg_half_sub.
1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077
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
1078
Lemma avg_half_sub_no_underflow_aux_aux: forall z:Z, (0 < z)%Z -> 
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1079
    (ZnearestE (IZR z / 2) < z)%Z.
1080 1081 1082
Proof.
intros z H1.
case (Zle_lt_or_eq 1 z); [omega|intros H2|intros H2].
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1083 1084 1085
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)))).
1086 1087 1088 1089 1090
right; ring.
apply Rle_lt_trans with (1:= RRle_abs _).
rewrite Rabs_Ropp.
apply Rle_lt_trans with (1:=Znearest_N (fun x => negb (Zeven x)) _).
apply Rle_lt_trans with (1*/2);[right; ring|idtac].
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1091
apply Rlt_le_trans with ((IZR z)*/2);[idtac|right; field].
1092
apply Rmult_lt_compat_r.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1093
lra.
1094
now apply IZR_lt.
1095 1096 1097 1098 1099 1100 1101 1102
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
1103
simpl; lra.
1104 1105 1106
Qed.


BOLDO Sylvie's avatar
BOLDO Sylvie committed
1107
Lemma avg_half_sub_no_underflow_aux1: forall f, format f -> 0 < f ->
1108 1109 1110 1111
  f <= round_flt (f/2) -> False.
Proof with auto with typeclass_instances.
intros f Ff Hf1 Hf2.
apply FLT_format_generic in Ff...
1112
destruct Ff as [g H1 H2 H3].
1113 1114 1115 1116 1117 1118 1119 1120
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
1121
lra.
1122 1123 1124 1125
apply generic_format_FLT.
exists (Float radix2 (Fnum g) (Fexp g-1)).
rewrite H1; unfold F2R; simpl.
unfold Zminus; rewrite bpow_plus.
1126 1127
apply Rmult_assoc.
easy.
1128 1129 1130 1131 1132 1133 1134 1135
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
1136 1137
apply IZR_lt.
replace (IZR (Fnum g) * bpow emin / 2 * bpow (- emin)) with (IZR (Fnum g) /2).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1138
apply avg_half_sub_no_underflow_aux_aux.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1139
apply lt_IZR.
1140 1141 1142 1143 1144
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
1145
unfold Rdiv; apply trans_eq with (IZR (Fnum g) * / 2 * (bpow (- emin)*bpow emin)).
1146 1147 1148 1149
rewrite <- bpow_plus.
ring_simplify (-emin+emin)%Z.
simpl; ring.
ring.
1150
apply sym_eq, cexp_FLT_FIX.
1151
apply Rgt_not_eq, Rlt_gt.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1152
lra.
1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163
rewrite H1; unfold F2R, Rdiv; simpl.
replace (/2) with (bpow (-1)) by reflexivity.
rewrite Rmult_assoc, <- bpow_plus.
rewrite Rabs_mult.
rewrite (Rabs_right (bpow _)).
2: apply Rle_ge, bpow_ge_0.
rewrite (Zplus_comm emin _).
rewrite (bpow_plus _ prec _).
apply Rmult_lt_compat.
apply Rabs_pos.
apply bpow_ge_0.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1164 1165
rewrite <- IZR_Zpower, <- abs_IZR.
now apply IZR_lt.
1166 1167 1168 1169 1170 1171
unfold Prec_gt_0 in prec_gt_0_; omega.
rewrite <- H4; apply bpow_lt.
omega.
Qed.


BOLDO Sylvie's avatar
BOLDO Sylvie committed
1172
Lemma avg_half_sub_no_underflow_aux2: forall u v, format u -> format v -> 
1173 1174
    (0 <= u /\ 0 <= v) \/ (u <= 0 /\ v <= 0) ->
    u <= v ->