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

BOLDO Sylvie's avatar
BOLDO Sylvie committed
5
Copyright (C) 2010-2011 Sylvie Boldo
6
#<br />#
BOLDO Sylvie's avatar
BOLDO Sylvie committed
7
Copyright (C) 2010-2011 Guillaume Melquiond
8 9 10 11 12 13 14 15 16 17 18 19

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:
58
  forall x,
BOLDO Sylvie's avatar
BOLDO Sylvie committed
59
    (0 < x)%R ->
60
    F x ->
BOLDO Sylvie's avatar
BOLDO Sylvie committed
61 62 63 64 65 66 67 68 69 70 71
    (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
now rewrite <- Fx.
Qed.

Theorem ulp_le_abs:
77
  forall x,
BOLDO Sylvie's avatar
BOLDO Sylvie committed
78
    (x <> 0)%R ->
79
    F x ->
BOLDO Sylvie's avatar
BOLDO Sylvie committed
80 81 82
    (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
Qed.

496
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:
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 :
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.
660
assert (He:(fexp (e-1) <= e-1)%Z).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
661 662
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 :
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 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 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 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801
Theorem generic_format_plus_ulp :
  forall { monotone_exp : Monotone_exp fexp },
  forall x, (x <> 0)%R ->
  F x -> F (x + ulp x).
Proof with auto with typeclass_instances.
intros monotone_exp x Zx Fx.
destruct (Rtotal_order x 0) as [Hx|[Hx|Hx]].
rewrite <- Ropp_involutive.
apply generic_format_opp.
rewrite Ropp_plus_distr, <- ulp_opp.
apply generic_format_opp in Fx.
destruct (Req_dec (-x) (bpow (ln_beta beta (-x) - 1))) as [Hx'|Hx'].
rewrite Hx' in Fx |- *.
apply generic_format_bpow_inv' in Fx...
unfold ulp, canonic_exp.
rewrite ln_beta_bpow.
revert Fx.
generalize (ln_beta_val _ _ (ln_beta beta (-x)) - 1)%Z.
clear -monotone_exp valid_exp.
intros e He.
destruct (Zle_lt_or_eq _ _ He) as [He1|He1].
assert (He2 : e = (e - fexp (e + 1) + fexp (e + 1))%Z) by ring.
rewrite He2 at 1.
rewrite bpow_plus.
assert (Hb := Z2R_Zpower beta _ (Zle_minus_le_0 _ _ He)).
match goal with |- F (?a * ?b + - ?b) =>
  replace (a * b + -b)%R with ((a - 1) * b)%R by ring end.
rewrite <- Hb.
rewrite <- (Z2R_minus _ 1).
change (F (F2R (Float beta (Zpower beta (e - fexp (e + 1)) - 1) (fexp (e + 1))))).
apply generic_format_F2R.
intros Zb.
unfold canonic_exp.
rewrite ln_beta_F2R with (1 := Zb).
rewrite (ln_beta_unique beta _ (e - fexp (e + 1))).
apply monotone_exp.
rewrite <- He2.
apply Zle_succ.
rewrite Rabs_pos_eq.
rewrite Z2R_minus, Hb.
split.
apply Rplus_le_reg_r with (- bpow (e - fexp (e + 1) - 1) + Z2R 1)%R.
apply Rmult_le_reg_r with (bpow (-(e - fexp (e+1) - 1))).
apply bpow_gt_0.
ring_simplify.
apply Rle_trans with R1.
rewrite Rmult_1_l.
apply (bpow_le _ _ 0).
clear -He1 ; omega.
rewrite Ropp_mult_distr_l_reverse.
rewrite <- 2!bpow_plus.
ring_simplify (e - fexp (e + 1) - 1 + - (e - fexp (e + 1) - 1))%Z.
ring_simplify (- (e - fexp (e + 1) - 1) + (e - fexp (e + 1)))%Z.
simpl.
rewrite <- (Z2R_plus (-1) _).
apply (Z2R_le 1).
unfold Zpower_pos, iter_pos.
generalize (Zle_bool_imp_le _ _ (radix_prop beta)).
clear ; omega.
rewrite <- (Rplus_0_r (bpow (e - fexp (e + 1)))) at 2.
apply Rplus_lt_compat_l.
now apply (Z2R_lt (-1) 0).
rewrite Z2R_minus.
apply Rle_0_minus.
rewrite Hb.
apply (bpow_le _ 0).
now apply Zle_minus_le_0.
rewrite He1, Rplus_opp_r.
apply generic_format_0.
apply generic_format_pred_1 ; try easy.
rewrite <- Ropp_0.
now apply Ropp_lt_contravar.
now elim Zx.
now apply generic_format_succ.
Qed.

802
Lemma pred_plus_ulp_1 :
803
  forall x, (0 < x)%R -> F x ->
804
  x <> bpow (ln_beta beta x - 1) ->
BOLDO Sylvie's avatar
BOLDO Sylvie committed
805
  ((x - ulp x) + ulp (x-ulp x) = x)%R.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
806
Proof.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
807 808 809 810
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
811
unfold canonic_exp; apply f_equal.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
812 813 814 815 816 817 818 819 820 821 822 823 824
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
825
assert (ex-1 < canonic_exp beta fexp x  < ex)%Z.
826
split ; apply (lt_bpow beta) ; rewrite <- H3 ; easy.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
827 828 829 830 831 832 833 834 835 836 837 838
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
839
pattern (bpow (canonic_exp beta fexp x)) at 1; rewrite <- Rmult_1_l.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
840 841 842 843 844
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
845
apply F2R_gt_0_reg with beta (canonic_exp beta fexp x).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
846 847 848
now rewrite <- Fx.
Qed.

849
Lemma pred_plus_ulp_2 :
850
  forall x, (0 < x)%R -> F x ->
851
  let e := ln_beta_val beta x (ln_beta beta x) in
BOLDO Sylvie's avatar
BOLDO Sylvie committed
852 853 854
  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
855
Proof.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
856 857 858
intros x Zx Fx e Hxe Zp.
replace (ulp (x - bpow (fexp (e - 1)))) with (bpow (fexp (e - 1))).
ring.
859
assert (He:(fexp (e-1) <= e-1)%Z).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
860 861
apply generic_format_bpow_inv with beta; trivial.
rewrite <- Hxe; assumption.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
862 863 864
case (Zle_lt_or_eq _ _ He); clear He; intros He.
(* *)
unfold ulp; apply f_equal.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
865
unfold canonic_exp; apply f_equal.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
866 867 868 869 870 871 872
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.
873
apply Rplus_le_compat; apply bpow_le; omega.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
874 875 876 877 878
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.
879
replace (bpow 1) with (Z2R beta).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
880 881 882 883 884 885
apply Z2R_le.
apply <- Zle_is_le_bool.
now destruct beta.
simpl.
unfold Zpower_pos; simpl.
now rewrite Zmult_1_r.
886
rewrite <- bpow_plus.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
887 888 889 890 891 892 893 894 895
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.
896
apply bpow_le.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
897 898 899 900 901 902
omega.
(* *)
contradict Zp.
rewrite Hxe, He; ring.
Qed.

903
Theorem pred_plus_ulp :
BOLDO Sylvie's avatar
BOLDO Sylvie committed
904
  forall x, (0 < x)%R -> F x ->
BOLDO Sylvie's avatar
BOLDO Sylvie committed
905 906
  (pred x <> 0)%R ->
  (pred x + ulp (pred x) = x)%R.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
907
Proof.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
908
intros x Zx Fx.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
909 910
unfold pred.
case Req_bool_spec; intros H Zp.
911 912
now apply pred_plus_ulp_2.
now apply pred_plus_ulp_1.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
913
Qed.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
914

915
Theorem pred_lt_id :
916
  forall x,
917
  (pred x < x)%R.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
918
Proof.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
919
intros.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
920 921 922
unfold pred.
case Req_bool_spec; intros H.
(* *)