Average.v 53.6 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 88 89 90
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.
rewrite <- Z2R_abs.
rewrite <- Z2R_Zpower.
now apply Z2R_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.
217
change 2 with (Z2R 2).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
218 219
rewrite <- (Z2R_abs 2).
now apply (Z2R_le 1 2).
220
(* *)
221 222 223
rewrite cexp_FLT_FIX.
rewrite cexp_FLT_FIX; trivial.
unfold FIX_exp, cexp; omega.
224 225
apply Rlt_le_trans with (1:=Hu).
apply bpow_le; omega.
226
lra.
227
rewrite Rabs_mult.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
228 229
rewrite Rabs_pos_eq.
2: now apply (Z2R_le 0 2).
230
apply Rlt_le_trans with (2*bpow (emin + prec - 1)).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
231
apply Rmult_lt_compat_l with (1 := Rlt_0_2).
232
assumption.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
233
change 2 with (bpow 1).
234 235 236 237 238
rewrite <- bpow_plus.
apply bpow_le; omega.
Qed.


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

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


593

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

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

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

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


627

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

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

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

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


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



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


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

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



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


729 730 731



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

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


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



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

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

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

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


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

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


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


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


BOLDO Sylvie's avatar
BOLDO Sylvie committed
1173
Lemma avg_half_sub_no_underflow_aux2: forall u v, format u -> format v -> 
1174 1175
    (0 <= u /\ 0 <= v) \/ (u <= 0 /\ v <= 0) ->
    u <= v ->
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1176
   (bpow emin) <= Rabs ((u+v)/2) -> avg_half_sub u v <> 0.
1177
Proof with auto with typeclass_instances.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1178
clear Fx Fy a av x y; intros x y Fx Fy same_sign xLey H; unfold avg_half_sub.
1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190
intros J.
apply round_plus_eq_zero in J...
2: apply generic_format_round...
assert (H1:x <= 0).
apply Rplus_le_reg_r with (round_flt (round_flt (y - x) / 2)).
rewrite J, Rplus_0_l.
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
1191
lra.
1192 1193 1194 1195
destruct H1 as [H1|H1].
(* *)
destruct same_sign as [(H2,H3)|(_,H2)].
contradict H2; now apply Rlt_not_le.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1196
apply avg_half_sub_no_underflow_aux1 with (-x).
1197 1198 1199 1200 1201 1202 1203
now apply generic_format_opp.
rewrite <- Ropp_0; now apply Ropp_lt_contravar.
apply Rle_trans with (round_flt (round_flt (y - x) / 2)).
apply Rplus_le_reg_l with x.
rewrite J; right; ring.
apply round_le...
unfold Rdiv; apply Rmult_le_compat_r.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1204
lra.
1205 1206 1207 1208 1209 1210 1211 1212