Mise à jour terminée. Pour connaître les apports de la version 13.8.4 par rapport à notre ancienne version vous pouvez lire les "Release Notes" suivantes :
https://about.gitlab.com/releases/2021/02/11/security-release-gitlab-13-8-4-released/
https://about.gitlab.com/releases/2021/02/05/gitlab-13-8-3-released/

Fappli_IEEE.v 49.5 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/

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

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

(** * IEEE-754 arithmetic *)
21
Require Import Fcore.
22
Require Import Fcore_digits.
23
Require Import Fcalc_digits.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
24 25 26
Require Import Fcalc_round.
Require Import Fcalc_bracket.
Require Import Fcalc_ops.
27
Require Import Fcalc_div.
28
Require Import Fcalc_sqrt.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
29
Require Import Fprop_relative.
30

31 32 33 34 35
Section AnyRadix.

Inductive full_float :=
  | F754_zero : bool -> full_float
  | F754_infinity : bool -> full_float
36
  | F754_nan : bool -> positive -> full_float
37 38
  | F754_finite : bool -> positive -> Z -> full_float.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
39
Definition FF2R beta x :=
40
  match x with
BOLDO Sylvie's avatar
BOLDO Sylvie committed
41
  | F754_finite s m e => F2R (Float beta (cond_Zopp s (Zpos m)) e)
42 43 44 45 46
  | _ => R0
  end.

End AnyRadix.

47 48
Section Binary.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
49 50 51
(** prec is the number of bits of the mantissa including the implicit one
    emax is the exponent of the infinities
    Typically p=24 and emax = 128 in single precision *)
52
Variable prec emax : Z.
53
Context (prec_gt_0_ : Prec_gt_0 prec).
54
Hypothesis Hmax : (prec < emax)%Z.
55

56
Let emin := (3 - emax - prec)%Z.
57
Let fexp := FLT_exp emin prec.
58
Instance fexp_correct : Valid_exp fexp := FLT_exp_valid emin prec.
59
Instance fexp_monotone : Monotone_exp fexp := FLT_exp_monotone emin prec.
60

BOLDO Sylvie's avatar
BOLDO Sylvie committed
61
Definition canonic_mantissa m e :=
62 63 64
  Zeq_bool (fexp (Z_of_nat (S (digits2_Pnat m)) + e)) e.

Definition bounded m e :=
BOLDO Sylvie's avatar
BOLDO Sylvie committed
65
  andb (canonic_mantissa m e) (Zle_bool e (emax - prec)).
66

67 68 69
Definition valid_binary x :=
  match x with
  | F754_finite _ m e => bounded m e
70
  | F754_nan _ pl => (Z_of_nat' (S (digits2_Pnat pl)) <? prec)%Z
71 72 73
  | _ => true
  end.

74
(** Basic type used for representing binary FP numbers.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
75 76
    Note that there is exactly one such object per FP datum.
    NaNs do not have any payload. They cannot be distinguished. *)
77 78 79

Definition nan_pl := {pl | (Z_of_nat' (S (digits2_Pnat pl)) <? prec)%Z  = true}.

80 81 82
Inductive binary_float :=
  | B754_zero : bool -> binary_float
  | B754_infinity : bool -> binary_float
83
  | B754_nan : bool -> nan_pl -> binary_float
84 85 86
  | B754_finite : bool ->
    forall (m : positive) (e : Z), bounded m e = true -> binary_float.

87 88 89 90 91
Definition FF2B x :=
  match x as x return valid_binary x = true -> binary_float with
  | F754_finite s m e => B754_finite s m e
  | F754_infinity s => fun _ => B754_infinity s
  | F754_zero s => fun _ => B754_zero s
92
  | F754_nan b pl => fun H => B754_nan b (exist _ pl H)
93 94 95 96 97 98 99
  end.

Definition B2FF x :=
  match x with
  | B754_finite s m e _ => F754_finite s m e
  | B754_infinity s => F754_infinity s
  | B754_zero s => F754_zero s
100
  | B754_nan b (exist pl _) => F754_nan b pl
101 102
  end.

103 104 105 106
Definition radix2 := Build_radix 2 (refl_equal true).

Definition B2R f :=
  match f with
Guillaume Melquiond's avatar
Guillaume Melquiond committed
107
  | B754_finite s m e _ => F2R (Float radix2 (cond_Zopp s (Zpos m)) e)
108 109 110
  | _ => R0
  end.

111 112 113 114
Theorem FF2R_B2FF :
  forall x,
  FF2R radix2 (B2FF x) = B2R x.
Proof.
115
now intros [sx|sx|sx [plx Hplx]|sx mx ex Hx].
116 117 118 119 120 121
Qed.

Theorem B2FF_FF2B :
  forall x Hx,
  B2FF (FF2B x Hx) = x.
Proof.
122
now intros [sx|sx|sx plx|sx mx ex] Hx.
123 124
Qed.

125 126 127 128
Theorem valid_binary_B2FF :
  forall x,
  valid_binary (B2FF x) = true.
Proof.
129
now intros [sx|sx|sx [plx Hplx]|sx mx ex Hx].
130 131 132 133 134 135
Qed.

Theorem FF2B_B2FF :
  forall x H,
  FF2B (B2FF x) H = x.
Proof.
136 137 138
intros [sx|sx|sx [plx Hplx]|sx mx ex Hx] H ; try easy.
simpl. apply f_equal, f_equal, eqbool_irrelevance.
apply f_equal, eqbool_irrelevance.
139 140 141 142 143 144 145 146 147 148
Qed.

Theorem FF2B_B2FF_valid :
  forall x,
  FF2B (B2FF x) (valid_binary_B2FF x) = x.
Proof.
intros x.
apply FF2B_B2FF.
Qed.

149 150 151 152
Theorem B2R_FF2B :
  forall x Hx,
  B2R (FF2B x Hx) = FF2R radix2 x.
Proof.
153
now intros [sx|sx|sx plx|sx mx ex] Hx.
154 155
Qed.

156 157 158 159 160
Theorem match_FF2B :
  forall {T} fz fi fn ff x Hx,
  match FF2B x Hx return T with
  | B754_zero sx => fz sx
  | B754_infinity sx => fi sx
161
  | B754_nan b (exist p _) => fn b p
162 163 164 165 166
  | B754_finite sx mx ex _ => ff sx mx ex
  end =
  match x with
  | F754_zero sx => fz sx
  | F754_infinity sx => fi sx
167
  | F754_nan b p => fn b p
168 169 170
  | F754_finite sx mx ex => ff sx mx ex
  end.
Proof.
171
now intros T fz fi fn ff [sx|sx|sx plx|sx mx ex] Hx.
172 173
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
174
Theorem canonic_canonic_mantissa :
175
  forall (sx : bool) mx ex,
BOLDO Sylvie's avatar
BOLDO Sylvie committed
176
  canonic_mantissa mx ex = true ->
Guillaume Melquiond's avatar
Guillaume Melquiond committed
177
  canonic radix2 fexp (Float radix2 (cond_Zopp sx (Zpos mx)) ex).
178 179 180 181 182 183 184
Proof.
intros sx mx ex H.
assert (Hx := Zeq_bool_eq _ _ H). clear H.
apply sym_eq.
simpl.
pattern ex at 2 ; rewrite <- Hx.
apply (f_equal fexp).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
185 186
rewrite ln_beta_F2R_Zdigits.
rewrite <- Zdigits_abs.
187 188 189 190 191 192 193 194 195
rewrite Z_of_nat_S_digits2_Pnat.
now case sx.
now case sx.
Qed.

Theorem generic_format_B2R :
  forall x,
  generic_format radix2 fexp (B2R x).
Proof.
196
intros [sx|sx|sx plx|sx mx ex Hx] ; try apply generic_format_0.
197 198
simpl.
apply generic_format_canonic.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
199
apply canonic_canonic_mantissa.
200 201 202
now destruct (andb_prop _ _ Hx) as (H, _).
Qed.

Guillaume Melquiond's avatar
Guillaume Melquiond committed
203 204 205
Theorem FLT_format_B2R :
  forall x,
  FLT_format radix2 emin prec (B2R x).
206
Proof with auto with typeclass_instances.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
207 208 209 210 211
intros x.
apply FLT_format_generic...
apply generic_format_B2R.
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
212
Theorem B2FF_inj :
213 214 215 216
  forall x y : binary_float,
  B2FF x = B2FF y ->
  x = y.
Proof.
217
intros [sx|sx|sx [plx Hplx]|sx mx ex Hx] [sy|sy|sy [ply Hply]|sy my ey Hy] ; try easy.
218 219 220 221 222 223 224 225 226 227
(* *)
intros H.
now inversion H.
(* *)
intros H.
now inversion H.
(* *)
intros H.
inversion H.
clear H.
228 229 230 231 232 233 234 235
revert Hplx.
rewrite H2.
intros Hx.
apply f_equal, f_equal, eqbool_irrelevance.
(* *)
intros H.
inversion H.
clear H.
236 237 238
revert Hx.
rewrite H2, H3.
intros Hx.
239
apply f_equal, eqbool_irrelevance.
240 241
Qed.

242 243 244 245 246 247
Definition is_finite_strict f :=
  match f with
  | B754_finite _ _ _ _ => true
  | _ => false
  end.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
248
Theorem B2R_inj:
249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275
  forall x y : binary_float,
  is_finite_strict x = true ->
  is_finite_strict y = true ->
  B2R x = B2R y ->
  x = y.
Proof.
intros [sx|sx| |sx mx ex Hx] [sy|sy| |sy my ey Hy] ; try easy.
simpl.
intros _ _ Heq.
assert (Hs: sx = sy).
(* *)
revert Heq. clear.
case sx ; case sy ; try easy ;
  intros Heq ; apply False_ind ; revert Heq.
apply Rlt_not_eq.
apply Rlt_trans with R0.
now apply F2R_lt_0_compat.
now apply F2R_gt_0_compat.
apply Rgt_not_eq.
apply Rgt_trans with R0.
now apply F2R_gt_0_compat.
now apply F2R_lt_0_compat.
assert (mx = my /\ ex = ey).
(* *)
refine (_ (canonic_unicity _ fexp _ _ _ _ Heq)).
rewrite Hs.
now case sy ; intro H ; injection H ; split.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
276
apply canonic_canonic_mantissa.
277
exact (proj1 (andb_prop _ _ Hx)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
278
apply canonic_canonic_mantissa.
279 280 281 282 283 284 285 286 287
exact (proj1 (andb_prop _ _ Hy)).
(* *)
revert Hx.
rewrite Hs, (proj1 H), (proj2 H).
intros Hx.
apply f_equal.
apply eqbool_irrelevance.
Qed.

288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310
Definition Bsign x :=
  match x with
  | B754_nan s _ => s
  | B754_zero s => s
  | B754_infinity s => s
  | B754_finite s _ _ _ => s
  end.

Definition sign_FF x :=
  match x with
  | F754_nan s _ => s
  | F754_zero s => s
  | F754_infinity s => s
  | F754_finite s _ _ => s
  end.

Theorem Bsign_FF2B :
  forall x H,
  Bsign (FF2B x H) = sign_FF x.
Proof.
now intros [sx|sx|sx plx|sx mx ex] H.
Qed.

311 312 313 314 315 316 317
Definition is_finite f :=
  match f with
  | B754_finite _ _ _ _ => true
  | B754_zero _ => true
  | _ => false
  end.

318 319 320 321 322 323 324 325 326 327 328 329 330 331
Definition is_finite_FF f :=
  match f with
  | F754_finite _ _ _ => true
  | F754_zero _ => true
  | _ => false
  end.

Theorem is_finite_FF2B :
  forall x Hx,
  is_finite (FF2B x Hx) = is_finite_FF x.
Proof.
now intros [| | |].
Qed.

332 333 334 335
Theorem is_finite_FF_B2FF :
  forall x,
  is_finite_FF (B2FF x) = is_finite x.
Proof.
336 337 338
now intros [| |? []|].
Qed.

339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358
Theorem B2R_Bsign_inj:
  forall x y : binary_float,
    is_finite x = true ->
    is_finite y = true ->
    B2R x = B2R y ->
    Bsign x = Bsign y ->
    x = y.
Proof.
intros. destruct x, y; try (apply B2R_inj; now eauto).
- simpl in H2. congruence.
- symmetry in H1. apply Rmult_integral in H1.
  destruct H1. apply eq_Z2R with (n:=0%Z) in H1. destruct b0; discriminate H1.
  simpl in H1. pose proof (bpow_gt_0 radix2 e).
  rewrite H1 in H3. apply Rlt_irrefl in H3. destruct H3.
- apply Rmult_integral in H1.
  destruct H1. apply eq_Z2R with (n:=0%Z) in H1. destruct b; discriminate H1.
  simpl in H1. pose proof (bpow_gt_0 radix2 e).
  rewrite H1 in H3. apply Rlt_irrefl in H3. destruct H3.
Qed.

359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374
Definition is_nan f :=
  match f with
  | B754_nan _ _ => true
  | _ => false
  end.

Definition is_nan_FF f :=
  match f with
  | F754_nan _ _ => true
  | _ => false
  end.

Theorem is_nan_FF2B :
  forall x Hx,
  is_nan (FF2B x Hx) = is_nan_FF x.
Proof.
375 376 377
now intros [| | |].
Qed.

378 379 380 381 382 383 384 385
Theorem is_nan_FF_B2FF :
  forall x,
  is_nan_FF (B2FF x) = is_nan x.
Proof.
now intros [| |? []|].
Qed.

Definition Bopp opp_nan x :=
Guillaume Melquiond's avatar
Guillaume Melquiond committed
386
  match x with
387 388
  | B754_nan sx plx =>
    let '(sres, plres) := opp_nan sx plx in B754_nan sres plres
Guillaume Melquiond's avatar
Guillaume Melquiond committed
389 390 391 392 393 394
  | B754_infinity sx => B754_infinity (negb sx)
  | B754_finite sx mx ex Hx => B754_finite (negb sx) mx ex Hx
  | B754_zero sx => B754_zero (negb sx)
  end.

Theorem Bopp_involutive :
395 396 397
  forall opp_nan x,
  is_nan x = false ->
  Bopp opp_nan (Bopp opp_nan x) = x.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
398
Proof.
399
now intros opp_nan [sx|sx|sx plx|sx mx ex Hx] ; simpl ; try rewrite Bool.negb_involutive.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
400 401 402
Qed.

Theorem B2R_Bopp :
403 404
  forall opp_nan x,
  B2R (Bopp opp_nan x) = (- B2R x)%R.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
405
Proof.
406 407
intros opp_nan [sx|sx|sx plx|sx mx ex Hx]; apply sym_eq ; try apply Ropp_0.
simpl. destruct opp_nan. apply Ropp_0.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
408
simpl.
409
rewrite <- F2R_opp.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
410 411 412
now case sx.
Qed.

413 414 415
Theorem is_finite_Bopp :
  forall opp_nan x,
  is_finite (Bopp opp_nan x) = is_finite x.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
416
Proof.
417 418 419 420
intros opp_nan [| | |] ; try easy.
intros s pl.
simpl.
now case opp_nan.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
421 422
Qed.

423 424 425
Theorem bounded_lt_emax :
  forall mx ex,
  bounded mx ex = true ->
426
  (F2R (Float radix2 (Zpos mx) ex) < bpow radix2 emax)%R.
427
Proof.
428
intros mx ex Hx.
429 430 431
destruct (andb_prop _ _ Hx) as (H1,H2).
generalize (Zeq_bool_eq _ _ H1). clear H1. intro H1.
generalize (Zle_bool_imp_le _ _ H2). clear H2. intro H2.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
432
generalize (ln_beta_F2R_Zdigits radix2 (Zpos mx) ex).
433 434 435 436
destruct (ln_beta radix2 (F2R (Float radix2 (Zpos mx) ex))) as (e',Ex).
unfold ln_beta_val.
intros H.
apply Rlt_le_trans with (bpow radix2 e').
437
change (Zpos mx) with (Zabs (Zpos mx)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
438
rewrite F2R_Zabs.
439 440 441 442 443 444 445 446 447
apply Ex.
apply Rgt_not_eq.
now apply F2R_gt_0_compat.
apply bpow_le.
rewrite H. 2: discriminate.
revert H1. clear -H2.
rewrite Z_of_nat_S_digits2_Pnat.
change Fcalc_digits.radix2 with radix2.
unfold fexp, FLT_exp.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
448
generalize (Zdigits radix2 (Zpos mx)).
449 450
intros ; zify ; subst.
clear -H H2. clearbody emin.
451 452 453
omega.
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
454
Theorem abs_B2R_lt_emax :
455
  forall x,
456
  (Rabs (B2R x) < bpow radix2 emax)%R.
457
Proof.
458
intros [sx|sx|sx plx|sx mx ex Hx] ; simpl ; try ( rewrite Rabs_R0 ; apply bpow_gt_0 ).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
459
rewrite <- F2R_Zabs, abs_cond_Zopp.
460 461 462
now apply bounded_lt_emax.
Qed.

463 464 465
Theorem bounded_canonic_lt_emax :
  forall mx ex,
  canonic radix2 fexp (Float radix2 (Zpos mx) ex) ->
466
  (F2R (Float radix2 (Zpos mx) ex) < bpow radix2 emax)%R ->
467 468
  bounded mx ex = true.
Proof.
469
intros mx ex Cx Bx.
470 471
apply andb_true_intro.
split.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
472
unfold canonic_mantissa.
473 474 475 476
unfold canonic, Fexp in Cx.
rewrite Cx at 2.
rewrite Z_of_nat_S_digits2_Pnat.
change Fcalc_digits.radix2 with radix2.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
477 478
unfold canonic_exp.
rewrite ln_beta_F2R_Zdigits. 2: discriminate.
479 480 481 482
now apply -> Zeq_is_eq_bool.
apply Zle_bool_true.
unfold canonic, Fexp in Cx.
rewrite Cx.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
483
unfold canonic_exp, fexp, FLT_exp.
484
destruct (ln_beta radix2 (F2R (Float radix2 (Zpos mx) ex))) as (e',Ex). simpl.
485 486
apply Zmax_lub.
cut (e' - 1 < emax)%Z. clear ; omega.
487 488 489
apply lt_bpow with radix2.
apply Rle_lt_trans with (2 := Bx).
change (Zpos mx) with (Zabs (Zpos mx)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
490
rewrite F2R_Zabs.
491 492 493
apply Ex.
apply Rgt_not_eq.
now apply F2R_gt_0_compat.
494 495 496
unfold emin.
generalize (prec_gt_0 prec).
clear -Hmax ; omega.
497 498
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
499
(** mantissa, round and sticky bits *)
500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550
Record shr_record := { shr_m : Z ; shr_r : bool ; shr_s : bool }.

Definition shr_1 mrs :=
  let '(Build_shr_record m r s) := mrs in
  let s := orb r s in
  match m with
  | Z0 => Build_shr_record Z0 false s
  | Zpos xH => Build_shr_record Z0 true s
  | Zpos (xO p) => Build_shr_record (Zpos p) false s
  | Zpos (xI p) => Build_shr_record (Zpos p) true s
  | Zneg xH => Build_shr_record Z0 true s
  | Zneg (xO p) => Build_shr_record (Zneg p) false s
  | Zneg (xI p) => Build_shr_record (Zneg p) true s
  end.

Definition loc_of_shr_record mrs :=
  match mrs with
  | Build_shr_record _ false false => loc_Exact
  | Build_shr_record _ false true => loc_Inexact Lt
  | Build_shr_record _ true false => loc_Inexact Eq
  | Build_shr_record _ true true => loc_Inexact Gt
  end.

Definition shr_record_of_loc m l :=
  match l with
  | loc_Exact => Build_shr_record m false false
  | loc_Inexact Lt => Build_shr_record m false true
  | loc_Inexact Eq => Build_shr_record m true false
  | loc_Inexact Gt => Build_shr_record m true true
  end.

Theorem shr_m_shr_record_of_loc :
  forall m l,
  shr_m (shr_record_of_loc m l) = m.
Proof.
now intros m [|[| |]].
Qed.

Theorem loc_of_shr_record_of_loc :
  forall m l,
  loc_of_shr_record (shr_record_of_loc m l) = l.
Proof.
now intros m [|[| |]].
Qed.

Definition shr mrs e n :=
  match n with
  | Zpos p => (iter_pos p _ shr_1 mrs, (e + n)%Z)
  | _ => (mrs, e)
  end.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
551
Lemma inbetween_shr_1 :
552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620
  forall x mrs e,
  (0 <= shr_m mrs)%Z ->
  inbetween_float radix2 (shr_m mrs) e x (loc_of_shr_record mrs) ->
  inbetween_float radix2 (shr_m (shr_1 mrs)) (e + 1) x (loc_of_shr_record (shr_1 mrs)).
Proof.
intros x mrs e Hm Hl.
refine (_ (new_location_even_correct (F2R (Float radix2 (shr_m (shr_1 mrs)) (e + 1))) (bpow radix2 e) 2 _ _ _ x (if shr_r (shr_1 mrs) then 1 else 0) (loc_of_shr_record mrs) _ _)) ; try easy.
2: apply bpow_gt_0.
2: now case (shr_r (shr_1 mrs)) ; split.
change (Z2R 2) with (bpow radix2 1).
rewrite <- bpow_plus.
rewrite (Zplus_comm 1), <- (F2R_bpow radix2 (e + 1)).
unfold inbetween_float, F2R. simpl.
rewrite Z2R_plus, Rmult_plus_distr_r.
replace (new_location_even 2 (if shr_r (shr_1 mrs) then 1%Z else 0%Z) (loc_of_shr_record mrs)) with (loc_of_shr_record (shr_1 mrs)).
easy.
clear -Hm.
destruct mrs as (m, r, s).
now destruct m as [|[m|m|]|m] ; try (now elim Hm) ; destruct r as [|] ; destruct s as [|].
rewrite (F2R_change_exp radix2 e).
2: apply Zle_succ.
unfold F2R. simpl.
rewrite <- 2!Rmult_plus_distr_r, <- 2!Z2R_plus.
rewrite Zplus_assoc.
replace (shr_m (shr_1 mrs) * 2 ^ (e + 1 - e) + (if shr_r (shr_1 mrs) then 1%Z else 0%Z))%Z with (shr_m mrs).
exact Hl.
ring_simplify (e + 1 - e)%Z.
change (2^1)%Z with 2%Z.
rewrite Zmult_comm.
clear -Hm.
destruct mrs as (m, r, s).
now destruct m as [|[m|m|]|m] ; try (now elim Hm) ; destruct r as [|] ; destruct s as [|].
Qed.

Theorem inbetween_shr :
  forall x m e l n,
  (0 <= m)%Z ->
  inbetween_float radix2 m e x l ->
  let '(mrs, e') := shr (shr_record_of_loc m l) e n in
  inbetween_float radix2 (shr_m mrs) e' x (loc_of_shr_record mrs).
Proof.
intros x m e l n Hm Hl.
destruct n as [|n|n].
now destruct l as [|[| |]].
2: now destruct l as [|[| |]].
unfold shr.
rewrite iter_nat_of_P.
rewrite Zpos_eq_Z_of_nat_o_nat_of_P.
induction (nat_of_P n).
simpl.
rewrite Zplus_0_r.
now destruct l as [|[| |]].
simpl iter_nat.
rewrite inj_S.
unfold Zsucc.
rewrite  Zplus_assoc.
revert IHn0.
apply inbetween_shr_1.
clear -Hm.
induction n0.
now destruct l as [|[| |]].
simpl.
revert IHn0.
generalize (iter_nat n0 shr_record shr_1 (shr_record_of_loc m l)).
clear.
intros (m, r, s) Hm.
now destruct m as [|[m|m|]|m] ; try (now elim Hm) ; destruct r as [|] ; destruct s as [|].
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
621
Definition Zdigits2 m :=
622 623
  match m with Z0 => m | Zpos p => Z_of_nat (S (digits2_Pnat p)) | Zneg p => Z_of_nat (S (digits2_Pnat p)) end.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
624
Theorem Zdigits2_Zdigits :
625
  forall m,
BOLDO Sylvie's avatar
BOLDO Sylvie committed
626
  Zdigits2 m = Zdigits radix2 m.
627
Proof.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
628
unfold Zdigits2.
629 630 631 632
intros [|m|m] ; try apply Z_of_nat_S_digits2_Pnat.
easy.
Qed.

633
Definition shr_fexp m e l :=
BOLDO Sylvie's avatar
BOLDO Sylvie committed
634
  shr (shr_record_of_loc m l) e (fexp (Zdigits2 m + e) - e).
635 636 637 638 639 640 641 642 643 644 645

Theorem shr_truncate :
  forall m e l,
  (0 <= m)%Z ->
  shr_fexp m e l =
  let '(m', e', l') := truncate radix2 fexp (m, e, l) in (shr_record_of_loc m' l', e').
Proof.
intros m e l Hm.
case_eq (truncate radix2 fexp (m, e, l)).
intros (m', e') l'.
unfold shr_fexp.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
646
rewrite Zdigits2_Zdigits.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
647
case_eq (fexp (Zdigits radix2 m + e) - e)%Z.
648 649 650 651 652 653 654 655 656
(* *)
intros He.
unfold truncate.
rewrite He.
simpl.
intros H.
now inversion H.
(* *)
intros p Hp.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
657
assert (He: (e <= fexp (Zdigits radix2 m + e))%Z).
658 659
clear -Hp ; zify ; omega.
destruct (inbetween_float_ex radix2 m e l) as (x, Hx).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
660
generalize (inbetween_shr x m e l (fexp (Zdigits radix2 m + e) - e) Hm Hx).
661 662 663 664
assert (Hx0 : (0 <= x)%R).
apply Rle_trans with (F2R (Float radix2 m e)).
now apply F2R_ge_0_compat.
exact (proj1 (inbetween_float_bounds _ _ _ _ _ Hx)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
665
case_eq (shr (shr_record_of_loc m l) e (fexp (Zdigits radix2 m + e) - e)).
666
intros mrs e'' H3 H4 H1.
667
generalize (truncate_correct radix2 _ x m e l Hx0 Hx (or_introl _ He)).
668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690
rewrite H1.
intros (H2,_).
rewrite <- Hp, H3.
assert (e'' = e').
change (snd (mrs, e'') = snd (fst (m',e',l'))).
rewrite <- H1, <- H3.
unfold truncate.
now rewrite Hp.
rewrite H in H4 |- *.
apply (f_equal (fun v => (v, _))).
destruct (inbetween_float_unique _ _ _ _ _ _ _ H2 H4) as (H5, H6).
rewrite H5, H6.
case mrs.
now intros m0 [|] [|].
(* *)
intros p Hp.
unfold truncate.
rewrite Hp.
simpl.
intros H.
now inversion H.
Qed.

691 692 693 694
Inductive mode := mode_NE | mode_ZR | mode_DN | mode_UP | mode_NA.

Definition round_mode m :=
  match m with
695 696 697 698 699
  | mode_NE => ZnearestE
  | mode_ZR => Ztrunc
  | mode_DN => Zfloor
  | mode_UP => Zceil
  | mode_NA => ZnearestA
700
  end.
701 702 703 704 705 706 707 708 709 710

Definition choice_mode m sx mx lx :=
  match m with
  | mode_NE => cond_incr (round_N (negb (Zeven mx)) lx) mx
  | mode_ZR => mx
  | mode_DN => cond_incr (round_sign_DN sx lx) mx
  | mode_UP => cond_incr (round_sign_UP sx lx) mx
  | mode_NA => cond_incr (round_N true lx) mx
  end.

711 712 713 714 715
Global Instance valid_rnd_round_mode : forall m, Valid_rnd (round_mode m).
Proof.
destruct m ; unfold round_mode ; auto with typeclass_instances.
Qed.

716 717 718 719 720 721 722 723 724 725 726 727 728
Definition overflow_to_inf m s :=
  match m with
  | mode_NE => true
  | mode_NA => true
  | mode_ZR => false
  | mode_UP => negb s
  | mode_DN => s
  end.

Definition binary_overflow m s :=
  if overflow_to_inf m s then F754_infinity s
  else F754_finite s (match (Zpower 2 prec - 1)%Z with Zpos p => p | _ => xH end) (emax - prec).

BOLDO Sylvie's avatar
BOLDO Sylvie committed
729
Definition binary_round_aux mode sx mx ex lx :=
730 731 732
  let '(mrs', e') := shr_fexp (Zpos mx) ex lx in
  let '(mrs'', e'') := shr_fexp (choice_mode mode sx (shr_m mrs') (loc_of_shr_record mrs')) e' loc_Exact in
  match shr_m mrs'' with
733
  | Z0 => F754_zero sx
734
  | Zpos m => if Zle_bool e'' (emax - prec) then F754_finite sx m e'' else binary_overflow mode sx
735
  | _ => F754_nan false xH (* dummy *)
736 737
  end.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
738
Theorem binary_round_aux_correct :
739 740
  forall mode x mx ex lx,
  inbetween_float radix2 (Zpos mx) ex (Rabs x) lx ->
BOLDO Sylvie's avatar
BOLDO Sylvie committed
741
  (ex <= fexp (Zdigits radix2 (Zpos mx) + ex))%Z ->
BOLDO Sylvie's avatar
BOLDO Sylvie committed
742
  let z := binary_round_aux mode (Rlt_bool x 0) mx ex lx in
743
  valid_binary z = true /\
744
  if Rlt_bool (Rabs (round radix2 fexp (round_mode mode) x)) (bpow radix2 emax) then
745
    FF2R radix2 z = round radix2 fexp (round_mode mode) x /\
746
    is_finite_FF z = true /\ sign_FF z = Rlt_bool x 0
747
  else
748
    z = binary_overflow mode (Rlt_bool x 0).
749
Proof with auto with typeclass_instances.
750
intros m x mx ex lx Bx Ex z.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
751
unfold binary_round_aux in z.
752
revert z.
753
rewrite shr_truncate. 2: easy.
754 755
refine (_ (round_trunc_sign_any_correct _ _ (round_mode m) (choice_mode m) _ x (Zpos mx) ex lx Bx (or_introl _ Ex))).
refine (_ (truncate_correct_partial _ _ _ _ _ _ _ Bx Ex)).
756
destruct (truncate radix2 fexp (Zpos mx, ex, lx)) as ((m1, e1), l1).
757
rewrite loc_of_shr_record_of_loc, shr_m_shr_record_of_loc.
758 759 760 761 762 763 764 765 766 767 768 769 770 771
set (m1' := choice_mode m (Rlt_bool x 0) m1 l1).
intros (H1a,H1b) H1c.
rewrite H1c.
assert (Hm: (m1 <= m1')%Z).
(* . *)
unfold m1', choice_mode, cond_incr.
case m ;
  try apply Zle_refl ;
  match goal with |- (m1 <= if ?b then _ else _)%Z =>
    case b ; [ apply Zle_succ | apply Zle_refl ] end.
assert (Hr: Rabs (round radix2 fexp (round_mode m) x) = F2R (Float radix2 m1' e1)).
(* . *)
rewrite <- (Zabs_eq m1').
replace (Zabs m1') with (Zabs (cond_Zopp (Rlt_bool x 0) m1')).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
772
rewrite F2R_Zabs.
773 774 775 776 777 778 779 780 781 782 783 784
now apply f_equal.
apply abs_cond_Zopp.
apply Zle_trans with (2 := Hm).
apply Zlt_succ_le.
apply F2R_gt_0_reg with radix2 e1.
apply Rle_lt_trans with (1 := Rabs_pos x).
exact (proj2 (inbetween_float_bounds _ _ _ _ _ H1a)).
(* . *)
assert (Br: inbetween_float radix2 m1' e1 (Rabs (round radix2 fexp (round_mode m) x)) loc_Exact).
now apply inbetween_Exact.
destruct m1' as [|m1'|m1'].
(* . m1' = 0 *)
785
rewrite shr_truncate. 2: apply Zle_refl.
786 787
generalize (truncate_0 radix2 fexp e1 loc_Exact).
destruct (truncate radix2 fexp (Z0, e1, loc_Exact)) as ((m2, e2), l2).
788
rewrite shr_m_shr_record_of_loc.
789
intros Hm2.
790
rewrite Hm2.
791
intros z.
792 793
repeat split.
rewrite Rlt_bool_true.
794
repeat split.
795 796
apply sym_eq.
case Rlt_bool ; apply F2R_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
797
rewrite <- F2R_Zabs, abs_cond_Zopp, F2R_0.
798
apply bpow_gt_0.
799
(* . 0 < m1' *)
BOLDO Sylvie's avatar
BOLDO Sylvie committed
800 801
assert (He: (e1 <= fexp (Zdigits radix2 (Zpos m1') + e1))%Z).
rewrite <- ln_beta_F2R_Zdigits, <- Hr, ln_beta_abs.
802 803
2: discriminate.
rewrite H1b.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
804 805 806
rewrite canonic_exp_abs.
fold (canonic_exp radix2 fexp (round radix2 fexp (round_mode m) x)).
apply canonic_exp_round_ge...
807 808 809 810 811 812
rewrite H1c.
case (Rlt_bool x 0).
apply Rlt_not_eq.
now apply F2R_lt_0_compat.
apply Rgt_not_eq.
now apply F2R_gt_0_compat.
813
refine (_ (truncate_correct_partial _ _ _ _ _ _ _ Br He)).
814 815 816
2: now rewrite Hr ; apply F2R_gt_0_compat.
refine (_ (truncate_correct_format radix2 fexp (Zpos m1') e1 _ _ He)).
2: discriminate.
817
rewrite shr_truncate. 2: easy.
818
destruct (truncate radix2 fexp (Zpos m1', e1, loc_Exact)) as ((m2, e2), l2).
819
rewrite shr_m_shr_record_of_loc.
820 821 822 823 824
intros (H3,H4) (H2,_).
destruct m2 as [|m2|m2].
elim Rgt_not_eq with (2 := H3).
rewrite F2R_0.
now apply F2R_gt_0_compat.
825
rewrite F2R_cond_Zopp, H3, abs_cond_Ropp, <- F2R_abs.
826
simpl Zabs.
827 828 829 830
case_eq (Zle_bool e2 (emax - prec)) ; intros He2.
assert (bounded m2 e2 = true).
apply andb_true_intro.
split.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
831
unfold canonic_mantissa.
832 833
apply Zeq_bool_true.
rewrite Z_of_nat_S_digits2_Pnat.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
834
rewrite <- ln_beta_F2R_Zdigits.
835 836 837 838 839
apply sym_eq.
now rewrite H3 in H4.
discriminate.
exact He2.
apply (conj H).
840
rewrite Rlt_bool_true.
841
repeat split.
842 843
apply F2R_cond_Zopp.
now apply bounded_lt_emax.
844 845 846 847
rewrite (Rlt_bool_false _ (bpow radix2 emax)).
refine (conj _ (refl_equal _)).
unfold binary_overflow.
case overflow_to_inf.
848
apply refl_equal.
849 850 851 852 853 854
unfold valid_binary, bounded.
rewrite Zle_bool_refl.
rewrite Bool.andb_true_r.
apply Zeq_bool_true.
rewrite Z_of_nat_S_digits2_Pnat.
change Fcalc_digits.radix2 with radix2.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
855
replace (Zdigits radix2 (Zpos (match (Zpower 2 prec - 1)%Z with Zpos p => p | _ => xH end))) with prec.
856
unfold fexp, FLT_exp, emin.
857 858
generalize (prec_gt_0 prec).
clear -Hmax ; zify ; omega.
859 860
change 2%Z with (radix_val radix2).
case_eq (Zpower radix2 prec - 1)%Z.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
861
simpl Zdigits.
862
generalize (Zpower_gt_1 radix2 prec (prec_gt_0 prec)).
863 864 865
clear ; omega.
intros p Hp.
apply Zle_antisym.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
866 867
cut (prec - 1 < Zdigits radix2 (Zpos p))%Z. clear ; omega.
apply Zdigits_gt_Zpower.
868 869 870 871 872 873 874
simpl Zabs. rewrite <- Hp.
cut (Zpower radix2 (prec - 1) < Zpower radix2 prec)%Z. clear ; omega.
apply lt_Z2R.
rewrite 2!Z2R_Zpower. 2: now apply Zlt_le_weak.
apply bpow_lt.
apply Zlt_pred.
now apply Zlt_0_le_0_pred.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
875
apply Zdigits_le_Zpower.
876 877 878
simpl Zabs. rewrite <- Hp.
apply Zlt_pred.
intros p Hp.
879
generalize (Zpower_gt_1 radix2 _ (prec_gt_0 prec)).
880
clear -Hp ; zify ; omega.
881 882
apply Rnot_lt_le.
intros Hx.
883 884 885 886
generalize (refl_equal (bounded m2 e2)).
unfold bounded at 2.
rewrite He2.
rewrite Bool.andb_false_r.
887 888
rewrite bounded_canonic_lt_emax with (2 := Hx).
discriminate.
889 890 891 892 893 894 895
unfold canonic.
now rewrite <- H3.
elim Rgt_not_eq with (2 := H3).
apply Rlt_trans with R0.
now apply F2R_lt_0_compat.
now apply F2R_gt_0_compat.
rewrite <- Hr.
896 897
apply generic_format_abs...
apply generic_format_round...
898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913
(* . not m1' < 0 *)
elim Rgt_not_eq with (2 := Hr).
apply Rlt_le_trans with R0.
now apply F2R_lt_0_compat.
apply Rabs_pos.
(* *)
apply Rlt_le_trans with (2 := proj1 (inbetween_float_bounds _ _ _ _ _ Bx)).
now apply F2R_gt_0_compat.
(* all the modes are valid *)
clear. case m.
exact inbetween_int_NE_sign.
exact inbetween_int_ZR_sign.
exact inbetween_int_DN_sign.
exact inbetween_int_UP_sign.
exact inbetween_int_NA_sign.
Qed.
914

BOLDO Sylvie's avatar
BOLDO Sylvie committed
915 916
(** Multiplication *)

917 918 919 920
Lemma Bmult_correct_aux :
  forall m sx mx ex (Hx : bounded mx ex = true) sy my ey (Hy : bounded my ey = true),
  let x := F2R (Float radix2 (cond_Zopp sx (Zpos mx)) ex) in
  let y := F2R (Float radix2 (cond_Zopp sy (Zpos my)) ey) in
BOLDO Sylvie's avatar
BOLDO Sylvie committed
921
  let z := binary_round_aux m (xorb sx sy) (mx * my) (ex + ey) loc_Exact in
922 923
  valid_binary z = true /\
  if Rlt_bool (Rabs (round radix2 fexp (round_mode m) (x * y))) (bpow radix2 emax) then
924
    FF2R radix2 z = round radix2 fexp (round_mode m) (x * y) /\
925
    is_finite_FF z = true /\ sign_FF z = xorb sx sy
926
  else
927
    z = binary_overflow m (xorb sx sy).
928
Proof.
929 930
intros m sx mx ex Hx sy my ey Hy x y.
unfold x, y.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
931
rewrite <- F2R_mult.
932
simpl.
933
replace (xorb sx sy) with (Rlt_bool (F2R (Float radix2 (cond_Zopp sx (Zpos mx) * cond_Zopp sy (Zpos my)) (ex + ey))) 0).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
934
apply binary_round_aux_correct.
935
constructor.
936
rewrite <- F2R_abs.
937 938 939 940 941
apply F2R_eq_compat.
rewrite Zabs_Zmult.
now rewrite 2!abs_cond_Zopp.
(* *)
change (Zpos (mx * my)) with (Zpos mx * Zpos my)%Z.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
942
assert (forall m e, bounded m e = true -> fexp (Zdigits radix2 (Zpos m) + e) = e)%Z.
943 944 945 946 947
clear. intros m e Hb.
destruct (andb_prop _ _ Hb) as (H,_).
apply Zeq_bool_eq.
now rewrite <- Z_of_nat_S_digits2_Pnat.
generalize (H _ _ Hx) (H _ _ Hy).
948
clear x y sx sy Hx Hy H.
949
unfold fexp, FLT_exp.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
950
refine (_ (Zdigits_mult_ge radix2 (Zpos mx) (Zpos my) _ _)) ; try discriminate.
951
refine (_ (Zdigits_gt_0 radix2 (Zpos mx) _) (Zdigits_gt_0 radix2 (Zpos my) _)) ; try discriminate.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
952
generalize (Zdigits radix2 (Zpos mx)) (Zdigits radix2 (Zpos my)) (Zdigits radix2 (Zpos mx * Zpos my)).
953
clear -Hmax.
954
unfold emin.
955 956 957 958 959 960 961 962 963 964 965 966 967 968
intros dx dy dxy Hx Hy Hxy.
zify ; intros ; subst.
omega.
(* *)
case sx ; case sy.
apply Rlt_bool_false.
now apply F2R_ge_0_compat.
apply Rlt_bool_true.
now apply F2R_lt_0_compat.
apply Rlt_bool_true.
now apply F2R_lt_0_compat.
apply Rlt_bool_false.
now apply F2R_ge_0_compat.
Qed.
969

970 971
Definition Bmult mult_nan m x y :=
  let f pl := B754_nan (fst pl) (snd pl) in
972
  match x, y with
973
  | B754_nan _ _, _ | _, B754_nan _ _ => f (mult_nan x y)
974 975 976
  | B754_infinity sx, B754_infinity sy => B754_infinity (xorb sx sy)
  | B754_infinity sx, B754_finite sy _ _ _ => B754_infinity (xorb sx sy)
  | B754_finite sx _ _ _, B754_infinity sy => B754_infinity (xorb sx sy)
977 978
  | B754_infinity _, B754_zero _ => f (mult_nan x y)
  | B754_zero _, B754_infinity _ => f (mult_nan x y)
979 980 981 982 983 984 985 986
  | B754_finite sx _ _ _, B754_zero sy => B754_zero (xorb sx sy)
  | B754_zero sx, B754_finite sy _ _ _ => B754_zero (xorb sx sy)
  | B754_zero sx, B754_zero sy => B754_zero (xorb sx sy)
  | B754_finite sx mx ex Hx, B754_finite sy my ey Hy =>
    FF2B _ (proj1 (Bmult_correct_aux m sx mx ex Hx sy my ey Hy))
  end.

Theorem Bmult_correct :
987
  forall mult_nan m x y,
988
  if Rlt_bool (Rabs (round radix2 fexp (round_mode m) (B2R x * B2R y))) (bpow radix2 emax) then
989
    B2R (Bmult mult_nan m x y) = round radix2 fexp (round_mode m) (B2R x * B2R y) /\
990 991 992
    is_finite (Bmult mult_nan m x y) = andb (is_finite x) (is_finite y) /\
    (is_nan (Bmult mult_nan m x y) = false ->
      Bsign (Bmult mult_nan m x y) = xorb (Bsign x) (Bsign y))
993
  else
994
    B2FF (Bmult mult_nan m x y) = binary_overflow m (xorb (Bsign x) (Bsign y)).
995
Proof.
996
intros mult_nan m [sx|sx|sx plx|sx mx ex Hx] [sy|sy|sy ply|sy my ey Hy] ;
997
  try ( rewrite ?Rmult_0_r, ?Rmult_0_l, round_0, Rabs_R0, Rlt_bool_true ; [ now repeat constructor | apply bpow_gt_0 | now auto with typeclass_instances ] ).
998 999
simpl.
case Bmult_correct_aux.
1000 1001
intros H1.
case Rlt_bool.
1002
intros (H2, (H3, H4)).
1003
split.
1004
now rewrite B2R_FF2B.
1005
split.
1006
now rewrite is_finite_FF2B.
1007
rewrite Bsign_FF2B. auto.
1008
intros H2.