Fcore_ulp.v 25.8 KB
Newer Older
1
(**
2 3 4 5
This file is part of the Flocq formalization of floating-point
arithmetic in Coq: http://flocq.gforge.inria.fr/

Copyright (C) 2010 Sylvie Boldo
6
#<br />#
7 8 9 10 11 12 13 14 15 16 17 18 19
Copyright (C) 2010 Guillaume Melquiond

This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 3 of the License, or (at your option) any later version.

This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
COPYING file for more details.
*)

BOLDO Sylvie's avatar
BOLDO Sylvie committed
20
(** * Unit in the Last Place: our definition using fexp and its properties, successor and predecessor *)
21 22 23 24 25
Require Import Fcore_Raux.
Require Import Fcore_defs.
Require Import Fcore_rnd.
Require Import Fcore_generic_fmt.
Require Import Fcore_float_prop.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
26

27
Section Fcore_ulp.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
28 29 30

Variable beta : radix.

31
Notation bpow e := (bpow beta e).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
32 33 34

Variable fexp : Z -> Z.

35
Context { valid_exp : Valid_exp fexp }.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
36

BOLDO Sylvie's avatar
BOLDO Sylvie committed
37
Definition ulp x := bpow (canonic_exp beta fexp x).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
38

39
Notation F := (generic_format beta fexp).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
40

41 42 43 44 45
Theorem ulp_opp :
  forall x, ulp (- x) = ulp x.
Proof.
intros x.
unfold ulp.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
46
now rewrite canonic_exp_opp.
47 48 49 50 51 52 53
Qed.

Theorem ulp_abs :
  forall x, ulp (Rabs x) = ulp x.
Proof.
intros x.
unfold ulp.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
54
now rewrite canonic_exp_abs.
55 56
Qed.

57
Theorem ulp_le_id:
BOLDO Sylvie's avatar
BOLDO Sylvie committed
58 59 60 61 62 63 64 65 66 67 68 69 70 71
  forall x, 
    (0 < x)%R ->
    F x -> 
    (ulp x <= x)%R.
Proof.
intros x Zx Fx.
rewrite <- (Rmult_1_l (ulp x)).
pattern x at 2; rewrite Fx.
unfold F2R, ulp; simpl.
apply Rmult_le_compat_r.
apply bpow_ge_0.
replace 1%R with (Z2R (Zsucc 0)) by reflexivity.
apply Z2R_le.
apply Zlt_le_succ.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
72
apply F2R_gt_0_reg with beta (canonic_exp beta fexp x).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
73 74 75 76 77 78 79 80 81 82
now rewrite <- Fx.
Qed.

Theorem ulp_le_abs:
  forall x, 
    (x <> 0)%R ->
    F x -> 
    (ulp x <= Rabs x)%R.
Proof.
intros x Zx Fx.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
83
rewrite <- ulp_abs.
84
apply ulp_le_id.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
85 86
now apply Rabs_pos_lt.
now apply generic_format_abs.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
87 88
Qed.

89 90
Theorem ulp_DN_UP :
  forall x, ~ F x ->
91
  round beta fexp Zceil x = (round beta fexp Zfloor x + ulp x)%R.
92
Proof.
93
intros x Fx.
94
unfold round, ulp. simpl.
95 96
unfold F2R. simpl.
rewrite Zceil_floor_neq.
97
rewrite Z2R_plus. simpl.
98
ring.
99
intros H.
100
apply Fx.
101
unfold generic_format, F2R. simpl.
102 103 104
rewrite <- H.
rewrite Ztrunc_Z2R.
rewrite H.
105
now rewrite scaled_mantissa_mult_bpow.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
106 107
Qed.

108
(** The successor of x is x + ulp x *)
Guillaume Melquiond's avatar
Guillaume Melquiond committed
109 110 111 112 113 114 115 116
Theorem succ_le_bpow :
  forall x e, (0 < x)%R -> F x ->
  (x < bpow e)%R ->
  (x + ulp x <= bpow e)%R.
Proof.
intros x e Zx Fx Hx.
pattern x at 1 ; rewrite Fx.
unfold ulp, F2R. simpl.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
117
pattern (bpow (canonic_exp beta fexp x)) at 2 ; rewrite <- Rmult_1_l.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
118 119
rewrite <- Rmult_plus_distr_r.
change 1%R with (Z2R 1).
120
rewrite <- Z2R_plus.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
121
change (F2R (Float beta (Ztrunc (scaled_mantissa beta fexp x) + 1) (canonic_exp beta fexp x)) <= bpow e)%R.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
122
apply F2R_p1_le_bpow.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
123
apply F2R_gt_0_reg with beta (canonic_exp beta fexp x).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
124 125 126 127 128
now rewrite <- Fx.
now rewrite <- Fx.
Qed.

Theorem ln_beta_succ :
129 130
  forall x, (0 < x)%R -> F x ->
  forall eps, (0 <= eps < ulp x)%R ->
131
  ln_beta beta (x + eps) = ln_beta beta x :> Z.
132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
Proof.
intros x Zx Fx eps Heps.
destruct (ln_beta beta x) as (ex, He).
simpl.
specialize (He (Rgt_not_eq _ _ Zx)).
apply ln_beta_unique.
rewrite Rabs_pos_eq.
rewrite Rabs_pos_eq in He.
split.
apply Rle_trans with (1 := proj1 He).
pattern x at 1 ; rewrite <- Rplus_0_r.
now apply Rplus_le_compat_l.
apply Rlt_le_trans with (x + ulp x)%R.
now apply Rplus_lt_compat_l.
pattern x at 1 ; rewrite Fx.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
147
unfold ulp, F2R. simpl.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
148
pattern (bpow (canonic_exp beta fexp x)) at 2 ; rewrite <- Rmult_1_l.
149 150
rewrite <- Rmult_plus_distr_r.
change 1%R with (Z2R 1).
151
rewrite <- Z2R_plus.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
152
change (F2R (Float beta (Ztrunc (scaled_mantissa beta fexp x) + 1) (canonic_exp beta fexp x)) <= bpow ex)%R.
153
apply F2R_p1_le_bpow.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
154
apply F2R_gt_0_reg with beta (canonic_exp beta fexp x).
155 156 157 158 159 160
now rewrite <- Fx.
now rewrite <- Fx.
now apply Rlt_le.
apply Rplus_le_le_0_compat.
now apply Rlt_le.
apply Heps.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
161 162
Qed.

163
Theorem round_DN_succ :
Guillaume Melquiond's avatar
Guillaume Melquiond committed
164 165
  forall x, (0 < x)%R -> F x ->
  forall eps, (0 <= eps < ulp x)%R ->
166
  round beta fexp Zfloor (x + eps) = x.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
167 168 169
Proof.
intros x Zx Fx eps Heps.
pattern x at 2 ; rewrite Fx.
170
unfold round.
171
unfold scaled_mantissa. simpl.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
172
unfold canonic_exp at 1 2.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
173 174
rewrite ln_beta_succ ; trivial.
apply (f_equal (fun m => F2R (Float beta m _))).
175 176 177 178 179 180 181 182
rewrite Ztrunc_floor.
apply Zfloor_imp.
split.
apply (Rle_trans _ _ _ (Zfloor_lb _)).
apply Rmult_le_compat_r.
apply bpow_ge_0.
pattern x at 1 ; rewrite <- Rplus_0_r.
now apply Rplus_le_compat_l.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
183
apply Rlt_le_trans with ((x + ulp x) * bpow (- canonic_exp beta fexp x))%R.
184 185 186 187
apply Rmult_lt_compat_r.
apply bpow_gt_0.
now apply Rplus_lt_compat_l.
rewrite Rmult_plus_distr_r.
188
rewrite Z2R_plus.
189 190 191 192
apply Rplus_le_compat.
pattern x at 1 3 ; rewrite Fx.
unfold F2R. simpl.
rewrite Rmult_assoc.
193
rewrite <- bpow_plus.
194 195 196 197 198
rewrite Zplus_opp_r.
rewrite Rmult_1_r.
rewrite Zfloor_Z2R.
apply Rle_refl.
unfold ulp.
199
rewrite <- bpow_plus.
200 201 202 203 204 205 206
rewrite Zplus_opp_r.
apply Rle_refl.
apply Rmult_le_pos.
now apply Rlt_le.
apply bpow_ge_0.
Qed.

207 208 209 210 211 212 213 214 215 216
Theorem generic_format_succ :
  forall x, (0 < x)%R -> F x ->
  F (x + ulp x).
Proof.
intros x Zx Fx.
destruct (ln_beta beta x) as (ex, Ex).
specialize (Ex (Rgt_not_eq _ _ Zx)).
assert (Ex' := Ex).
rewrite Rabs_pos_eq in Ex'.
destruct (succ_le_bpow x ex) ; try easy.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
217
unfold generic_format, scaled_mantissa, canonic_exp.
218 219
rewrite ln_beta_unique with beta (x + ulp x)%R ex.
pattern x at 1 3 ; rewrite Fx.
220
unfold ulp, scaled_mantissa.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
221
rewrite canonic_exp_fexp with (1 := Ex).
222
unfold F2R. simpl.
223 224
rewrite Rmult_plus_distr_r.
rewrite Rmult_assoc.
225
rewrite <- bpow_plus, Zplus_opp_r, Rmult_1_r.
226
change (bpow 0) with (Z2R 1).
227
rewrite <- Z2R_plus.
228
rewrite Ztrunc_Z2R.
229
rewrite Z2R_plus.
230 231 232 233 234 235 236 237 238 239 240 241 242 243
rewrite Rmult_plus_distr_r.
now rewrite Rmult_1_l.
rewrite Rabs_pos_eq.
split.
apply Rle_trans with (1 := proj1 Ex').
pattern x at 1 ; rewrite <- Rplus_0_r.
apply Rplus_le_compat_l.
apply bpow_ge_0.
exact H.
apply Rplus_le_le_0_compat.
now apply Rlt_le.
apply bpow_ge_0.
rewrite H.
apply generic_format_bpow.
244
apply valid_exp.
245 246 247 248 249 250 251
destruct (Zle_or_lt ex (fexp ex)) ; trivial.
elim Rlt_not_le with (1 := Zx).
rewrite Fx.
replace (Ztrunc (scaled_mantissa beta fexp x)) with Z0.
rewrite F2R_0.
apply Rle_refl.
unfold scaled_mantissa.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
252
rewrite canonic_exp_fexp with (1 := Ex).
253 254 255 256 257 258 259 260
destruct (mantissa_small_pos beta fexp x ex) ; trivial.
rewrite Ztrunc_floor.
apply sym_eq.
apply Zfloor_imp.
split.
now apply Rlt_le.
exact H2.
now apply Rlt_le.
261 262 263
now apply Rlt_le.
Qed.

264
Theorem round_UP_succ :
265 266
  forall x, (0 < x)%R -> F x ->
  forall eps, (0 < eps <= ulp x)%R ->
267 268
  round beta fexp Zceil (x + eps) = (x + ulp x)%R.
Proof with auto with typeclass_instances.
269 270 271 272 273
intros x Zx Fx eps (Heps1,[Heps2|Heps2]).
assert (Heps: (0 <= eps < ulp x)%R).
split.
now apply Rlt_le.
exact Heps2.
274
assert (Hd := round_DN_succ x Zx Fx eps Heps).
275 276
rewrite ulp_DN_UP.
rewrite Hd.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
277
unfold ulp, canonic_exp.
278 279
now rewrite ln_beta_succ.
intros Fs.
280
rewrite round_generic in Hd...
281 282 283 284
apply Rgt_not_eq with (2 := Hd).
pattern x at 2 ; rewrite <- Rplus_0_r.
now apply Rplus_lt_compat_l.
rewrite Heps2.
285
apply round_generic...
286 287 288
now apply generic_format_succ.
Qed.

289
Theorem succ_le_lt:
BOLDO Sylvie's avatar
BOLDO Sylvie committed
290 291 292 293
  forall x y,
  F x -> F y ->
  (0 < x < y)%R ->
  (x + ulp x <= y)%R.
294
Proof with auto with typeclass_instances.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
295 296 297 298 299 300 301 302
intros x y Hx Hy H.
case (Rle_or_lt (ulp x) (y-x)); intros H1.
apply Rplus_le_reg_r with (-x)%R.
now ring_simplify (x+ulp x + -x)%R.
replace y with (x+(y-x))%R by ring.
absurd (x < y)%R.
2: apply H.
apply Rle_not_lt; apply Req_le.
303
rewrite <- round_DN_succ with (eps:=(y-x)%R); try easy.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
304
ring_simplify (x+(y-x))%R.
305 306
apply sym_eq.
apply round_generic...
BOLDO Sylvie's avatar
BOLDO Sylvie committed
307 308 309 310
split; trivial.
apply Rlt_le; now apply Rlt_Rminus.
Qed.

311
(** Error of a rounding, expressed in number of ulps *)
Guillaume Melquiond's avatar
Guillaume Melquiond committed
312
Theorem ulp_error :
313 314 315 316
  forall rnd { Zrnd : Valid_rnd rnd } x,
  (Rabs (round beta fexp rnd x - x) < ulp x)%R.
Proof with auto with typeclass_instances.
intros rnd Zrnd x.
317
destruct (generic_format_EM beta fexp x) as [Hx|Hx].
318
(* x = rnd x *)
319
rewrite round_generic...
320 321 322
unfold Rminus.
rewrite Rplus_opp_r, Rabs_R0.
apply bpow_gt_0.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
323
(* x <> rnd x *)
324
destruct (round_DN_or_UP beta fexp rnd x) as [H|H] ; rewrite H ; clear H.
325 326 327
(* . *)
rewrite Rabs_left1.
rewrite Ropp_minus_distr.
328
apply Rplus_lt_reg_r with (round beta fexp Zfloor x).
329 330
rewrite <- ulp_DN_UP with (1 := Hx).
ring_simplify.
331 332
assert (Hu: (x <= round beta fexp Zceil x)%R).
apply round_UP_pt...
333 334 335 336
destruct Hu as [Hu|Hu].
exact Hu.
elim Hx.
rewrite Hu.
337
apply generic_format_round...
338
apply Rle_minus.
339
apply round_DN_pt...
340
(* . *)
Guillaume Melquiond's avatar
Guillaume Melquiond committed
341
rewrite Rabs_pos_eq.
342
rewrite ulp_DN_UP with (1 := Hx).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
343
apply Rplus_lt_reg_r with (x - ulp x)%R.
344
ring_simplify.
345 346
assert (Hd: (round beta fexp Zfloor x <= x)%R).
apply round_DN_pt...
347 348 349 350
destruct Hd as [Hd|Hd].
exact Hd.
elim Hx.
rewrite <- Hd.
351
apply generic_format_round...
Guillaume Melquiond's avatar
Guillaume Melquiond committed
352
apply Rle_0_minus.
353
apply round_UP_pt...
Guillaume Melquiond's avatar
Guillaume Melquiond committed
354 355
Qed.

Guillaume Melquiond's avatar
Guillaume Melquiond committed
356 357
Theorem ulp_half_error :
  forall choice x,
358 359
  (Rabs (round beta fexp (Znearest choice) x - x) <= /2 * ulp x)%R.
Proof with auto with typeclass_instances.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
360
intros choice x.
361
destruct (generic_format_EM beta fexp x) as [Hx|Hx].
Guillaume Melquiond's avatar
Guillaume Melquiond committed
362
(* x = rnd x *)
363
rewrite round_generic...
Guillaume Melquiond's avatar
Guillaume Melquiond committed
364
unfold Rminus.
365
rewrite Rplus_opp_r, Rabs_R0.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
366 367 368 369
apply Rmult_le_pos.
apply Rlt_le.
apply Rinv_0_lt_compat.
now apply (Z2R_lt 0 2).
370
apply bpow_ge_0.
371
(* x <> rnd x *)
372
set (d := round beta fexp Zfloor x).
373
destruct (round_N_pt beta fexp choice x) as (Hr1, Hr2).
374 375 376 377
destruct (Rle_or_lt (x - d) (d + ulp x - x)) as [H|H].
(* . rnd(x) = rndd(x) *)
apply Rle_trans with (Rabs (d - x)).
apply Hr2.
378
apply (round_DN_pt beta fexp x).
379 380 381 382 383 384 385 386 387
rewrite Rabs_left1.
rewrite Ropp_minus_distr.
apply Rmult_le_reg_r with 2%R.
now apply (Z2R_lt 0 2).
apply Rplus_le_reg_r with (d - x)%R.
ring_simplify.
apply Rle_trans with (1 := H).
right. field.
apply Rle_minus.
388
apply (round_DN_pt beta fexp x).
389
(* . rnd(x) = rndu(x) *)
390
assert (Hu: (d + ulp x)%R = round beta fexp Zceil x).
391 392 393 394 395
unfold d.
now rewrite <- ulp_DN_UP.
apply Rle_trans with (Rabs (d + ulp x - x)).
apply Hr2.
rewrite Hu.
396
apply (round_UP_pt beta fexp x).
397 398 399 400 401 402 403 404 405 406
rewrite Rabs_pos_eq.
apply Rmult_le_reg_r with 2%R.
now apply (Z2R_lt 0 2).
apply Rplus_le_reg_r with (- (d + ulp x - x))%R.
ring_simplify.
apply Rlt_le.
apply Rlt_le_trans with (1 := H).
right. field.
apply Rle_0_minus.
rewrite Hu.
407
apply (round_UP_pt beta fexp x).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
408 409
Qed.

410
Theorem ulp_le :
411
  forall { Hm : Monotone_exp fexp },
Guillaume Melquiond's avatar
Guillaume Melquiond committed
412 413 414 415 416
  forall x y: R,
  (0 < x)%R -> (x <= y)%R ->
  (ulp x <= ulp y)%R.
Proof.
intros Hm x y Hx Hxy.
417
apply bpow_le.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
418
apply Hm.
419
now apply ln_beta_le.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
420 421
Qed.

422
Theorem ulp_bpow :
423 424 425
  forall e, ulp (bpow e) = bpow (fexp (e + 1)).
intros e.
unfold ulp.
426
apply f_equal.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
427
apply canonic_exp_fexp.
428 429
rewrite Rabs_pos_eq.
split.
430 431
ring_simplify (e + 1 - 1)%Z.
apply Rle_refl.
432
apply bpow_lt.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
433
apply Zlt_succ.
434
apply bpow_ge_0.
435 436
Qed.

437
Theorem ulp_DN :
438
  forall x,
439 440
  (0 < round beta fexp Zfloor x)%R ->
  ulp (round beta fexp Zfloor x) = ulp x.
441
Proof.
442
intros x Hd.
443
unfold ulp.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
444
now rewrite canonic_exp_DN with (2 := Hd).
445 446
Qed.

447
Theorem ulp_error_f :
448 449 450 451 452 453
  forall { Hm : Monotone_exp fexp } rnd { Zrnd : Valid_rnd rnd } x,
  (round beta fexp rnd x <> 0)%R ->
  (Rabs (round beta fexp rnd x - x) < ulp (round beta fexp rnd x))%R.
Proof with auto with typeclass_instances.
intros Hm rnd Zrnd x Hfx.
case (round_DN_or_UP beta fexp rnd x); intros Hx.
454
(* *)
455
case (Rle_or_lt 0 (round beta fexp Zfloor x)).
456 457
intros H; destruct H.
rewrite Hx at 2.
458
rewrite ulp_DN; trivial.
459
apply ulp_error...
460 461 462
rewrite Hx in Hfx; contradict Hfx; auto with real.
intros H.
apply Rlt_le_trans with (1:=ulp_error _ _).
463
rewrite <- (ulp_opp x), <- (ulp_opp (round beta fexp rnd x)).
464
apply ulp_le; trivial.
465 466 467 468
apply Ropp_0_gt_lt_contravar; apply Rlt_gt.
case (Rle_or_lt 0 x); trivial.
intros H1; contradict H.
apply Rle_not_lt.
469
apply round_ge_generic...
470 471
apply generic_format_0.
apply Ropp_le_contravar; rewrite Hx.
472
apply (round_DN_pt beta fexp x).
473
(* *)
474
rewrite Hx; case (Rle_or_lt 0 (round beta fexp Zceil x)).
475 476
intros H; destruct H.
apply Rlt_le_trans with (1:=ulp_error _ _).
477
apply ulp_le; trivial.
478 479 480
case (Rle_or_lt x 0); trivial.
intros H1; contradict H.
apply Rle_not_lt.
481
apply round_le_generic...
482
apply generic_format_0.
483
apply round_UP_pt...
484 485
rewrite Hx in Hfx; contradict Hfx; auto with real.
intros H.
486
rewrite <- (ulp_opp (round beta fexp Zceil x)).
487
rewrite <- round_DN_opp.
488
rewrite ulp_DN; trivial.
489
replace (round beta fexp Zceil x - x)%R with (-((round beta fexp Zfloor (-x) - (-x))))%R.
490
rewrite Rabs_Ropp.
491
apply ulp_error...
492 493
rewrite round_DN_opp; ring.
rewrite round_DN_opp; apply Ropp_0_gt_lt_contravar; apply Rlt_gt; assumption.
494 495 496
Qed.

Theorem ulp_half_error_f : 
497
  forall { Hm : Monotone_exp fexp },
498
  forall choice x,
499 500 501
  (round beta fexp (Znearest choice) x <> 0)%R ->
  (Rabs (round beta fexp (Znearest choice) x - x) <= /2 * ulp (round beta fexp (Znearest choice) x))%R.
Proof with auto with typeclass_instances.
502
intros Hm choice x Hfx.
503
case (round_DN_or_UP beta fexp (Znearest choice) x); intros Hx.
504
(* *)
505
case (Rle_or_lt 0 (round beta fexp Zfloor x)).
506 507
intros H; destruct H.
rewrite Hx at 2.
508
rewrite ulp_DN; trivial.
509 510 511 512 513 514
apply ulp_half_error.
rewrite Hx in Hfx; contradict Hfx; auto with real.
intros H.
apply Rle_trans with (1:=ulp_half_error _ _).
apply Rmult_le_compat_l.
auto with real.
515
rewrite <- (ulp_opp x), <- (ulp_opp (round beta fexp (Znearest choice) x)).
516
apply ulp_le; trivial.
517 518 519 520
apply Ropp_0_gt_lt_contravar; apply Rlt_gt.
case (Rle_or_lt 0 x); trivial.
intros H1; contradict H.
apply Rle_not_lt.
521
apply round_ge_generic...
522 523
apply generic_format_0.
apply Ropp_le_contravar; rewrite Hx.
524
apply (round_DN_pt beta fexp x).
525
(* *)
526
case (Rle_or_lt 0 (round beta fexp Zceil x)).
527 528 529 530
intros H; destruct H.
apply Rle_trans with (1:=ulp_half_error _ _).
apply Rmult_le_compat_l.
auto with real.
531
apply ulp_le; trivial.
532 533 534
case (Rle_or_lt x 0); trivial.
intros H1; contradict H.
apply Rle_not_lt.
535
apply round_le_generic...
536
apply generic_format_0.
537
rewrite Hx; apply (round_UP_pt beta fexp x).
538 539
rewrite Hx in Hfx; contradict Hfx; auto with real.
intros H.
540
rewrite Hx at 2; rewrite <- (ulp_opp (round beta fexp Zceil x)).
541
rewrite <- round_DN_opp.
542
rewrite ulp_DN; trivial.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
543
pattern x at 1 2; rewrite <- Ropp_involutive.
544
rewrite round_N_opp.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
545 546
unfold Rminus.
rewrite <- Ropp_plus_distr, Rabs_Ropp.
547
apply ulp_half_error.
548
rewrite round_DN_opp; apply Ropp_0_gt_lt_contravar; apply Rlt_gt; assumption.
549 550
Qed.

551
(** Predecessor *)
BOLDO Sylvie's avatar
BOLDO Sylvie committed
552
Definition pred x :=
553 554 555 556
  if Req_bool x (bpow (ln_beta beta x - 1)) then
    (x - bpow (fexp (ln_beta beta x - 1)))%R
  else
    (x - ulp x)%R.
BOLDO Sylvie's avatar
Fpred  
BOLDO Sylvie committed
557

BOLDO Sylvie's avatar
BOLDO Sylvie committed
558 559 560 561 562 563 564 565 566 567
Theorem pred_ge_bpow :
  forall x e,  F x ->
  x <> ulp x ->
  (bpow e < x)%R ->
  (bpow e <= x - ulp x)%R.
Proof.
intros x e Fx Hx' Hx.
(* *)
assert (1 <= Ztrunc (scaled_mantissa beta fexp x))%Z.
assert (0 <  Ztrunc (scaled_mantissa beta fexp x))%Z.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
568
apply F2R_gt_0_reg with beta (canonic_exp beta fexp x).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
569 570 571 572 573 574 575 576
rewrite <- Fx.
apply Rle_lt_trans with (2:=Hx).
apply bpow_ge_0.
omega.
case (Zle_lt_or_eq _ _ H); intros Hm.
(* *)
pattern x at 1 ; rewrite Fx.
unfold ulp, F2R. simpl.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
577
pattern (bpow (canonic_exp beta fexp x)) at 2 ; rewrite <- Rmult_1_l.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
578 579
rewrite <- Rmult_minus_distr_r.
change 1%R with (Z2R 1).
580
rewrite <- Z2R_minus.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
581
change (bpow e <= F2R (Float beta (Ztrunc (scaled_mantissa beta fexp x) - 1) (canonic_exp beta fexp x)))%R.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
582 583 584 585 586 587 588 589 590 591
apply bpow_le_F2R_m1; trivial.
now rewrite <- Fx.
(* *)
contradict Hx'.
pattern x at 1; rewrite Fx.
rewrite  <- Hm.
unfold ulp, F2R; simpl.
now rewrite Rmult_1_l.
Qed.

592
Lemma generic_format_pred_1:
BOLDO Sylvie's avatar
BOLDO Sylvie committed
593
  forall x, (0 < x)%R -> F x -> 
594
  x <> bpow (ln_beta beta x - 1) ->
BOLDO Sylvie's avatar
BOLDO Sylvie committed
595
  F (x - ulp x).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
596
Proof.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
597 598 599 600 601 602 603 604 605
intros x Zx Fx Hx.
destruct (ln_beta beta x) as (ex, Ex).
simpl in Hx.
specialize (Ex (Rgt_not_eq _ _ Zx)).
assert (Ex' : (bpow (ex - 1) < x < bpow ex)%R).
rewrite Rabs_pos_eq in Ex.
destruct Ex as (H,H'); destruct H; split; trivial.
contradict Hx; easy.
now apply Rlt_le.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
606
unfold generic_format, scaled_mantissa, canonic_exp.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
607 608 609
rewrite ln_beta_unique with beta (x - ulp x)%R ex.
pattern x at 1 3 ; rewrite Fx.
unfold ulp, scaled_mantissa.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
610
rewrite canonic_exp_fexp with (1 := Ex).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
611 612 613
unfold F2R. simpl.
rewrite Rmult_minus_distr_r.
rewrite Rmult_assoc.
614
rewrite <- bpow_plus, Zplus_opp_r, Rmult_1_r.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
615
change (bpow 0) with (Z2R 1).
616
rewrite <- Z2R_minus.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
617
rewrite Ztrunc_Z2R.
618
rewrite Z2R_minus.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
619 620 621 622 623 624
rewrite Rmult_minus_distr_r.
now rewrite Rmult_1_l.
rewrite Rabs_pos_eq.
split.
apply pred_ge_bpow; trivial.
unfold ulp; intro H.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
625
assert (ex-1 < canonic_exp beta fexp x  < ex)%Z.
626
split ; apply (lt_bpow beta) ; rewrite <- H ; easy.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
627 628 629 630 631 632 633 634 635 636 637
clear -H0. omega.
apply Ex'.
apply Rle_lt_trans with (2 := proj2 Ex').
pattern x at 3 ; rewrite <- Rplus_0_r.
apply Rplus_le_compat_l.
rewrite <-Ropp_0.
apply Ropp_le_contravar.
apply bpow_ge_0.
apply Rle_0_minus.
pattern x at 2; rewrite Fx.
unfold ulp, F2R; simpl.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
638
pattern (bpow (canonic_exp beta fexp x)) at 1; rewrite <- Rmult_1_l.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
639 640 641 642 643
apply Rmult_le_compat_r.
apply bpow_ge_0.
replace 1%R with (Z2R 1) by reflexivity.
apply Z2R_le.
assert (0 <  Ztrunc (scaled_mantissa beta fexp x))%Z.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
644
apply F2R_gt_0_reg with beta (canonic_exp beta fexp x).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
645 646 647 648 649 650
rewrite <- Fx.
apply Rle_lt_trans with (2:=proj1 Ex').
apply bpow_ge_0.
omega.
Qed.

651
Lemma generic_format_pred_2 :
BOLDO Sylvie's avatar
BOLDO Sylvie committed
652
  forall x, (0 < x)%R -> F x -> 
653
  let e := ln_beta_val beta x (ln_beta beta x) in
BOLDO Sylvie's avatar
BOLDO Sylvie committed
654
  x =  bpow (e - 1) ->
655
  F (x - bpow (fexp (e - 1))).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
656
Proof.
BOLDO Sylvie's avatar
Fpred  
BOLDO Sylvie committed
657 658
intros x Zx Fx e Hx.
pose (f:=(x - bpow (fexp (e - 1)))%R).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
659
fold f.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
660 661 662
assert (He:(fexp (e-1) <= e-1)%Z). 
apply generic_format_bpow_inv with beta; trivial.
rewrite <- Hx; assumption.
BOLDO Sylvie's avatar
Fpred  
BOLDO Sylvie committed
663
case (Zle_lt_or_eq _ _ He); clear He; intros He.
664
assert (f = F2R (Float beta (Zpower beta (e-1-(fexp (e-1))) -1) (fexp (e-1))))%R.
BOLDO Sylvie's avatar
Fpred  
BOLDO Sylvie committed
665 666
unfold f; rewrite Hx.
unfold F2R; simpl.
667
rewrite Z2R_minus, Z2R_Zpower.
BOLDO Sylvie's avatar
Fpred  
BOLDO Sylvie committed
668
rewrite Rmult_minus_distr_r, Rmult_1_l.
669
rewrite <- bpow_plus.
BOLDO Sylvie's avatar
Fpred  
BOLDO Sylvie committed
670 671 672
now replace (e - 1 - fexp (e - 1) + fexp (e - 1))%Z with (e-1)%Z by ring.
omega.
rewrite H.
673 674
apply generic_format_F2R.
intros _.
BOLDO Sylvie's avatar
Fpred  
BOLDO Sylvie committed
675
apply Zeq_le.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
676
apply canonic_exp_fexp.
BOLDO Sylvie's avatar
Fpred  
BOLDO Sylvie committed
677 678 679 680 681 682 683
rewrite <- H.
unfold f; rewrite Hx.
rewrite Rabs_right.
split.
apply Rplus_le_reg_l with (bpow (fexp (e-1))).
ring_simplify.
apply Rle_trans with (bpow (e - 2) + bpow (e - 2))%R.
684
apply Rplus_le_compat ; apply bpow_le ; omega.
BOLDO Sylvie's avatar
Fpred  
BOLDO Sylvie committed
685 686 687 688 689
apply Rle_trans with (2*bpow (e - 2))%R;[right; ring|idtac].
apply Rle_trans with (bpow 1*bpow (e - 2))%R.
apply Rmult_le_compat_r.
apply bpow_ge_0.
replace 2%R with (Z2R 2) by reflexivity.
690
replace (bpow 1) with (Z2R beta).
BOLDO Sylvie's avatar
Fpred  
BOLDO Sylvie committed
691 692 693 694 695 696
apply Z2R_le.
apply <- Zle_is_le_bool.
now destruct beta.
simpl.
unfold Zpower_pos; simpl.
now rewrite Zmult_1_r.
697
rewrite <- bpow_plus.
BOLDO Sylvie's avatar
Fpred  
BOLDO Sylvie committed
698 699 700 701 702 703 704 705
replace (1+(e-2))%Z with (e-1)%Z by ring.
now right.
rewrite <- Rplus_0_r.
apply Rplus_lt_compat_l.
rewrite <- Ropp_0.
apply Ropp_lt_contravar.
apply bpow_gt_0.
apply Rle_ge; apply Rle_0_minus.
706
apply bpow_le.
BOLDO Sylvie's avatar
Fpred  
BOLDO Sylvie committed
707 708 709 710 711 712 713
omega.
replace f with 0%R.
apply generic_format_0.
unfold f.
rewrite Hx, He.
ring.
Qed.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
714

715
Theorem generic_format_pred :
BOLDO Sylvie's avatar
Fpred  
BOLDO Sylvie committed
716
  forall x, (0 < x)%R -> F x -> 
BOLDO Sylvie's avatar
BOLDO Sylvie committed
717
  F (pred x).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
718
Proof.
BOLDO Sylvie's avatar
Fpred  
BOLDO Sylvie committed
719
intros x Zx Fx.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
720 721
unfold pred.
case Req_bool_spec; intros H.
722 723
now apply generic_format_pred_2.
now apply generic_format_pred_1.
BOLDO Sylvie's avatar
Fpred  
BOLDO Sylvie committed
724
Qed.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
725

726
Lemma pred_plus_ulp_1 :
BOLDO Sylvie's avatar
BOLDO Sylvie committed
727
  forall x, (0 < x)%R -> F x -> 
728
  x <> bpow (ln_beta beta x - 1) ->
BOLDO Sylvie's avatar
BOLDO Sylvie committed
729
  ((x - ulp x) + ulp (x-ulp x) = x)%R.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
730
Proof.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
731 732 733 734
intros x Zx Fx Hx.
replace (ulp (x - ulp x)) with (ulp x).
ring.
unfold ulp at 1 2; apply f_equal.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
735
unfold canonic_exp; apply f_equal.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
736 737 738 739 740 741 742 743 744 745 746 747 748
destruct (ln_beta beta x) as (ex, Hex).
simpl in *.
assert (x <> 0)%R by auto with real.
specialize (Hex H).
apply sym_eq.
apply ln_beta_unique.
rewrite Rabs_right.
rewrite Rabs_right in Hex.
2: apply Rle_ge; apply Rlt_le; easy.
split.
destruct Hex as ([H1|H1],H2).
apply pred_ge_bpow; trivial.
unfold ulp; intros H3.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
749
assert (ex-1 < canonic_exp beta fexp x  < ex)%Z.
750
split ; apply (lt_bpow beta) ; rewrite <- H3 ; easy.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
751 752 753 754 755 756 757 758 759 760 761 762
omega.
contradict Hx; auto with real.
apply Rle_lt_trans with (2:=proj2 Hex).
rewrite <- Rplus_0_r.
apply Rplus_le_compat_l.
rewrite <- Ropp_0.
apply Ropp_le_contravar.
apply bpow_ge_0.
apply Rle_ge.
apply Rle_0_minus.
pattern x at 2; rewrite Fx.
unfold ulp, F2R; simpl.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
763
pattern (bpow (canonic_exp beta fexp x)) at 1; rewrite <- Rmult_1_l.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
764 765 766 767 768
apply Rmult_le_compat_r.
apply bpow_ge_0.
replace 1%R with (Z2R (Zsucc 0)) by reflexivity.
apply Z2R_le.
apply Zlt_le_succ.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
769
apply F2R_gt_0_reg with beta (canonic_exp beta fexp x).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
770 771 772
now rewrite <- Fx.
Qed.

773
Lemma pred_plus_ulp_2 :
BOLDO Sylvie's avatar
BOLDO Sylvie committed
774
  forall x, (0 < x)%R -> F x -> 
775
  let e := ln_beta_val beta x (ln_beta beta x) in
BOLDO Sylvie's avatar
BOLDO Sylvie committed
776 777 778
  x =  bpow (e - 1) ->
  (x - bpow (fexp (e-1)) <> 0)%R ->
  ((x - bpow (fexp (e-1))) + ulp (x - bpow (fexp (e-1))) = x)%R.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
779
Proof.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
780 781 782
intros x Zx Fx e Hxe Zp.
replace (ulp (x - bpow (fexp (e - 1)))) with (bpow (fexp (e - 1))).
ring.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
783 784 785
assert (He:(fexp (e-1) <= e-1)%Z). 
apply generic_format_bpow_inv with beta; trivial.
rewrite <- Hxe; assumption.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
786 787 788
case (Zle_lt_or_eq _ _ He); clear He; intros He.
(* *)
unfold ulp; apply f_equal.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
789
unfold canonic_exp; apply f_equal.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
790 791 792 793 794 795 796
apply sym_eq.
apply ln_beta_unique.
rewrite Rabs_right.
split.
apply Rplus_le_reg_l with (bpow (fexp (e-1))).
ring_simplify.
apply Rle_trans with (bpow (e - 2) + bpow (e - 2))%R.
797
apply Rplus_le_compat; apply bpow_le; omega.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
798 799 800 801 802
apply Rle_trans with (2*bpow (e - 2))%R;[right; ring|idtac].
apply Rle_trans with (bpow 1*bpow (e - 2))%R.
apply Rmult_le_compat_r.
apply bpow_ge_0.
replace 2%R with (Z2R 2) by reflexivity.
803
replace (bpow 1) with (Z2R beta).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
804 805 806 807 808 809
apply Z2R_le.
apply <- Zle_is_le_bool.
now destruct beta.
simpl.
unfold Zpower_pos; simpl.
now rewrite Zmult_1_r.
810
rewrite <- bpow_plus.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
811 812 813 814 815 816 817 818 819
replace (1+(e-2))%Z with (e-1)%Z by ring.
now right.
rewrite <- Rplus_0_r, Hxe.
apply Rplus_lt_compat_l.
rewrite <- Ropp_0.
apply Ropp_lt_contravar.
apply bpow_gt_0.
apply Rle_ge; apply Rle_0_minus.
rewrite Hxe.
820
apply bpow_le.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
821 822 823 824 825 826
omega.
(* *)
contradict Zp.
rewrite Hxe, He; ring.
Qed.

827
Theorem pred_plus_ulp :
BOLDO Sylvie's avatar
BOLDO Sylvie committed
828
  forall x, (0 < x)%R -> F x ->
BOLDO Sylvie's avatar
BOLDO Sylvie committed
829 830
  (pred x <> 0)%R ->
  (pred x + ulp (pred x) = x)%R.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
831
Proof.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
832
intros x Zx Fx.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
833 834
unfold pred.
case Req_bool_spec; intros H Zp.
835 836
now apply pred_plus_ulp_2.
now apply pred_plus_ulp_1.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
837
Qed.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
838

839
Theorem pred_lt_id :
BOLDO Sylvie's avatar
BOLDO Sylvie committed
840
  forall x, 
841
  (pred x < x)%R.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
842
Proof.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
843
intros.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
844 845 846
unfold pred.
case Req_bool_spec; intros H.
(* *)
BOLDO Sylvie's avatar
BOLDO Sylvie committed
847 848 849 850 851
rewrite <- Rplus_0_r.
apply Rplus_lt_compat_l.
rewrite <- Ropp_0.
apply Ropp_lt_contravar.
apply bpow_gt_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
852
(* *)
BOLDO Sylvie's avatar
BOLDO Sylvie committed
853 854 855 856 857 858
rewrite <- Rplus_0_r.
apply Rplus_lt_compat_l.
rewrite <- Ropp_0.
apply Ropp_lt_contravar.
apply bpow_gt_0.
Qed.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
859

860
Theorem pred_ge_0 :
BOLDO Sylvie's avatar
BOLDO Sylvie committed
861
  forall x, 
862
  (0 < x)%R -> F x -> (0 <= pred x)%R.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
863 864 865 866 867 868
intros x Zx Fx.
unfold pred.
case Req_bool_spec; intros H.
(* *)
apply Rle_0_minus.
rewrite H.
869
apply bpow_le.
870
destruct (ln_beta beta x) as (ex,Ex) ; simpl.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
871
rewrite ln_beta_bpow.
872
ring_simplify (ex - 1 + 1 - 1)%Z.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
873 874 875 876
apply generic_format_bpow_inv with beta; trivial.
simpl in H.
rewrite <- H; assumption.
apply Rle_0_minus.
877
now apply ulp_le_id.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
878 879
Qed.

880
Theorem round_UP_pred :
BOLDO Sylvie's avatar
BOLDO Sylvie committed
881 882
  forall x, (0 < pred x)%R -> F x ->
  forall eps, (0 < eps <= ulp (pred x) )%R ->
883
  round beta fexp Zceil (pred x + eps) = x.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
884 885
Proof.
intros x Hx Fx eps Heps.
886
rewrite round_UP_succ; trivial.
887
apply pred_plus_ulp; trivial.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
888
apply Rlt_trans with (1:=Hx).
889
apply pred_lt_id.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
890
now apply Rgt_not_eq.
891
apply generic_format_pred; trivial.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
892
apply Rlt_trans with (1:=Hx).