Fappli_IEEE.v 46.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 288 289 290 291 292 293 294
exact (proj1 (andb_prop _ _ Hy)).
(* *)
revert Hx.
rewrite Hs, (proj1 H), (proj2 H).
intros Hx.
apply f_equal.
apply eqbool_irrelevance.
Qed.

Definition is_finite f :=
  match f with
  | B754_finite _ _ _ _ => true
  | B754_zero _ => true
  | _ => false
  end.

295 296 297 298 299 300 301 302 303 304 305 306 307 308
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.

309 310 311 312
Theorem is_finite_FF_B2FF :
  forall x,
  is_finite_FF (B2FF x) = is_finite x.
Proof.
313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331
now intros [| |? []|].
Qed.

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.
332 333 334
now intros [| | |].
Qed.

335 336 337 338 339 340 341 342
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
343
  match x with
344 345
  | B754_nan sx plx =>
    let '(sres, plres) := opp_nan sx plx in B754_nan sres plres
Guillaume Melquiond's avatar
Guillaume Melquiond committed
346 347 348 349 350 351
  | 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 :
352 353 354
  forall opp_nan x,
  is_nan x = false ->
  Bopp opp_nan (Bopp opp_nan x) = x.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
355
Proof.
356
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
357 358 359
Qed.

Theorem B2R_Bopp :
360 361
  forall opp_nan x,
  B2R (Bopp opp_nan x) = (- B2R x)%R.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
362
Proof.
363 364
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
365
simpl.
366
rewrite <- F2R_opp.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
367 368 369
now case sx.
Qed.

370 371 372
Theorem is_finite_Bopp :
  forall opp_nan x,
  is_finite (Bopp opp_nan x) = is_finite x.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
373
Proof.
374 375 376 377
intros opp_nan [| | |] ; try easy.
intros s pl.
simpl.
now case opp_nan.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
378 379
Qed.

380 381 382
Theorem bounded_lt_emax :
  forall mx ex,
  bounded mx ex = true ->
383
  (F2R (Float radix2 (Zpos mx) ex) < bpow radix2 emax)%R.
384
Proof.
385
intros mx ex Hx.
386 387 388
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
389
generalize (ln_beta_F2R_Zdigits radix2 (Zpos mx) ex).
390 391 392 393
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').
394
change (Zpos mx) with (Zabs (Zpos mx)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
395
rewrite F2R_Zabs.
396 397 398 399 400 401 402 403 404
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
405
generalize (Zdigits radix2 (Zpos mx)).
406 407
intros ; zify ; subst.
clear -H H2. clearbody emin.
408 409 410
omega.
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
411
Theorem abs_B2R_lt_emax :
412
  forall x,
413
  (Rabs (B2R x) < bpow radix2 emax)%R.
414
Proof.
415
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
416
rewrite <- F2R_Zabs, abs_cond_Zopp.
417 418 419
now apply bounded_lt_emax.
Qed.

420 421 422
Theorem bounded_canonic_lt_emax :
  forall mx ex,
  canonic radix2 fexp (Float radix2 (Zpos mx) ex) ->
423
  (F2R (Float radix2 (Zpos mx) ex) < bpow radix2 emax)%R ->
424 425
  bounded mx ex = true.
Proof.
426
intros mx ex Cx Bx.
427 428
apply andb_true_intro.
split.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
429
unfold canonic_mantissa.
430 431 432 433
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
434 435
unfold canonic_exp.
rewrite ln_beta_F2R_Zdigits. 2: discriminate.
436 437 438 439
now apply -> Zeq_is_eq_bool.
apply Zle_bool_true.
unfold canonic, Fexp in Cx.
rewrite Cx.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
440
unfold canonic_exp, fexp, FLT_exp.
441
destruct (ln_beta radix2 (F2R (Float radix2 (Zpos mx) ex))) as (e',Ex). simpl.
442 443
apply Zmax_lub.
cut (e' - 1 < emax)%Z. clear ; omega.
444 445 446
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
447
rewrite F2R_Zabs.
448 449 450
apply Ex.
apply Rgt_not_eq.
now apply F2R_gt_0_compat.
451 452 453
unfold emin.
generalize (prec_gt_0 prec).
clear -Hmax ; omega.
454 455
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
456
(** mantissa, round and sticky bits *)
457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507
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
508
Lemma inbetween_shr_1 :
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 551 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
  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
578
Definition Zdigits2 m :=
579 580
  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
581
Theorem Zdigits2_Zdigits :
582
  forall m,
BOLDO Sylvie's avatar
BOLDO Sylvie committed
583
  Zdigits2 m = Zdigits radix2 m.
584
Proof.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
585
unfold Zdigits2.
586 587 588 589
intros [|m|m] ; try apply Z_of_nat_S_digits2_Pnat.
easy.
Qed.

590
Definition shr_fexp m e l :=
BOLDO Sylvie's avatar
BOLDO Sylvie committed
591
  shr (shr_record_of_loc m l) e (fexp (Zdigits2 m + e) - e).
592 593 594 595 596 597 598 599 600 601 602

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
603
rewrite Zdigits2_Zdigits.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
604
case_eq (fexp (Zdigits radix2 m + e) - e)%Z.
605 606 607 608 609 610 611 612 613
(* *)
intros He.
unfold truncate.
rewrite He.
simpl.
intros H.
now inversion H.
(* *)
intros p Hp.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
614
assert (He: (e <= fexp (Zdigits radix2 m + e))%Z).
615 616
clear -Hp ; zify ; omega.
destruct (inbetween_float_ex radix2 m e l) as (x, Hx).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
617
generalize (inbetween_shr x m e l (fexp (Zdigits radix2 m + e) - e) Hm Hx).
618 619 620 621
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
622
case_eq (shr (shr_record_of_loc m l) e (fexp (Zdigits radix2 m + e) - e)).
623
intros mrs e'' H3 H4 H1.
624
generalize (truncate_correct radix2 _ x m e l Hx0 Hx (or_introl _ He)).
625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647
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.

648 649 650 651
Inductive mode := mode_NE | mode_ZR | mode_DN | mode_UP | mode_NA.

Definition round_mode m :=
  match m with
652 653 654 655 656
  | mode_NE => ZnearestE
  | mode_ZR => Ztrunc
  | mode_DN => Zfloor
  | mode_UP => Zceil
  | mode_NA => ZnearestA
657
  end.
658 659 660 661 662 663 664 665 666 667

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.

668 669 670 671 672
Global Instance valid_rnd_round_mode : forall m, Valid_rnd (round_mode m).
Proof.
destruct m ; unfold round_mode ; auto with typeclass_instances.
Qed.

673 674 675 676 677 678 679 680 681 682 683 684 685
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
686
Definition binary_round_aux mode sx mx ex lx :=
687 688 689
  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
690
  | Z0 => F754_zero sx
691
  | Zpos m => if Zle_bool e'' (emax - prec) then F754_finite sx m e'' else binary_overflow mode sx
692
  | _ => F754_nan false xH (* dummy *)
693 694
  end.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
695
Theorem binary_round_aux_correct :
696 697
  forall mode x mx ex lx,
  inbetween_float radix2 (Zpos mx) ex (Rabs x) lx ->
BOLDO Sylvie's avatar
BOLDO Sylvie committed
698
  (ex <= fexp (Zdigits radix2 (Zpos mx) + ex))%Z ->
BOLDO Sylvie's avatar
BOLDO Sylvie committed
699
  let z := binary_round_aux mode (Rlt_bool x 0) mx ex lx in
700
  valid_binary z = true /\
701
  if Rlt_bool (Rabs (round radix2 fexp (round_mode mode) x)) (bpow radix2 emax) then
702 703
    FF2R radix2 z = round radix2 fexp (round_mode mode) x /\
    is_finite_FF z = true
704
  else
705
    z = binary_overflow mode (Rlt_bool x 0).
706
Proof with auto with typeclass_instances.
707
intros m x mx ex lx Bx Ex z.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
708
unfold binary_round_aux in z.
709
revert z.
710
rewrite shr_truncate. 2: easy.
711 712
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)).
713
destruct (truncate radix2 fexp (Zpos mx, ex, lx)) as ((m1, e1), l1).
714
rewrite loc_of_shr_record_of_loc, shr_m_shr_record_of_loc.
715 716 717 718 719 720 721 722 723 724 725 726 727 728
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
729
rewrite F2R_Zabs.
730 731 732 733 734 735 736 737 738 739 740 741
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 *)
742
rewrite shr_truncate. 2: apply Zle_refl.
743 744
generalize (truncate_0 radix2 fexp e1 loc_Exact).
destruct (truncate radix2 fexp (Z0, e1, loc_Exact)) as ((m2, e2), l2).
745
rewrite shr_m_shr_record_of_loc.
746
intros Hm2.
747
rewrite Hm2.
748
intros z.
749 750
repeat split.
rewrite Rlt_bool_true.
751
repeat split.
752 753
apply sym_eq.
case Rlt_bool ; apply F2R_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
754
rewrite <- F2R_Zabs, abs_cond_Zopp, F2R_0.
755
apply bpow_gt_0.
756
(* . 0 < m1' *)
BOLDO Sylvie's avatar
BOLDO Sylvie committed
757 758
assert (He: (e1 <= fexp (Zdigits radix2 (Zpos m1') + e1))%Z).
rewrite <- ln_beta_F2R_Zdigits, <- Hr, ln_beta_abs.
759 760
2: discriminate.
rewrite H1b.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
761 762 763
rewrite canonic_exp_abs.
fold (canonic_exp radix2 fexp (round radix2 fexp (round_mode m) x)).
apply canonic_exp_round_ge...
764 765 766 767 768 769
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.
770
refine (_ (truncate_correct_partial _ _ _ _ _ _ _ Br He)).
771 772 773
2: now rewrite Hr ; apply F2R_gt_0_compat.
refine (_ (truncate_correct_format radix2 fexp (Zpos m1') e1 _ _ He)).
2: discriminate.
774
rewrite shr_truncate. 2: easy.
775
destruct (truncate radix2 fexp (Zpos m1', e1, loc_Exact)) as ((m2, e2), l2).
776
rewrite shr_m_shr_record_of_loc.
777 778 779 780 781
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.
782
rewrite F2R_cond_Zopp, H3, abs_cond_Ropp, <- F2R_abs.
783
simpl Zabs.
784 785 786 787
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
788
unfold canonic_mantissa.
789 790
apply Zeq_bool_true.
rewrite Z_of_nat_S_digits2_Pnat.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
791
rewrite <- ln_beta_F2R_Zdigits.
792 793 794 795 796
apply sym_eq.
now rewrite H3 in H4.
discriminate.
exact He2.
apply (conj H).
797
rewrite Rlt_bool_true.
798
repeat split.
799 800
apply F2R_cond_Zopp.
now apply bounded_lt_emax.
801 802 803 804
rewrite (Rlt_bool_false _ (bpow radix2 emax)).
refine (conj _ (refl_equal _)).
unfold binary_overflow.
case overflow_to_inf.
805
apply refl_equal.
806 807 808 809 810 811
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
812
replace (Zdigits radix2 (Zpos (match (Zpower 2 prec - 1)%Z with Zpos p => p | _ => xH end))) with prec.
813
unfold fexp, FLT_exp, emin.
814 815
generalize (prec_gt_0 prec).
clear -Hmax ; zify ; omega.
816 817
change 2%Z with (radix_val radix2).
case_eq (Zpower radix2 prec - 1)%Z.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
818
simpl Zdigits.
819
generalize (Zpower_gt_1 radix2 prec (prec_gt_0 prec)).
820 821 822
clear ; omega.
intros p Hp.
apply Zle_antisym.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
823 824
cut (prec - 1 < Zdigits radix2 (Zpos p))%Z. clear ; omega.
apply Zdigits_gt_Zpower.
825 826 827 828 829 830 831
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
832
apply Zdigits_le_Zpower.
833 834 835
simpl Zabs. rewrite <- Hp.
apply Zlt_pred.
intros p Hp.
836
generalize (Zpower_gt_1 radix2 _ (prec_gt_0 prec)).
837
clear -Hp ; zify ; omega.
838 839
apply Rnot_lt_le.
intros Hx.
840 841 842 843
generalize (refl_equal (bounded m2 e2)).
unfold bounded at 2.
rewrite He2.
rewrite Bool.andb_false_r.
844 845
rewrite bounded_canonic_lt_emax with (2 := Hx).
discriminate.
846 847 848 849 850 851 852
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.
853 854
apply generic_format_abs...
apply generic_format_round...
855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870
(* . 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.
871

872 873
Definition Bsign x :=
  match x with
874
  | B754_nan s _ => s
875 876 877 878 879
  | B754_zero s => s
  | B754_infinity s => s
  | B754_finite s _ _ _ => s
  end.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
880
Definition sign_FF x :=
881
  match x with
882
  | F754_nan s _ => s
883 884 885 886 887 888 889
  | F754_zero s => s
  | F754_infinity s => s
  | F754_finite s _ _ => s
  end.

Theorem Bsign_FF2B :
  forall x H,
BOLDO Sylvie's avatar
BOLDO Sylvie committed
890
  Bsign (FF2B x H) = sign_FF x.
891
Proof.
892
now intros [sx|sx|sx plx|sx mx ex] H.
893 894
Qed.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
895 896
(** Multiplication *)

897 898 899 900
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
901
  let z := binary_round_aux m (xorb sx sy) (mx * my) (ex + ey) loc_Exact in
902 903
  valid_binary z = true /\
  if Rlt_bool (Rabs (round radix2 fexp (round_mode m) (x * y))) (bpow radix2 emax) then
904 905
    FF2R radix2 z = round radix2 fexp (round_mode m) (x * y) /\
    is_finite_FF z = true
906
  else
907
    z = binary_overflow m (xorb sx sy).
908
Proof.
909 910
intros m sx mx ex Hx sy my ey Hy x y.
unfold x, y.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
911
rewrite <- F2R_mult.
912
simpl.
913
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
914
apply binary_round_aux_correct.
915
constructor.
916
rewrite <- F2R_abs.
917 918 919 920 921
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
922
assert (forall m e, bounded m e = true -> fexp (Zdigits radix2 (Zpos m) + e) = e)%Z.
923 924 925 926 927
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).
928
clear x y sx sy Hx Hy H.
929
unfold fexp, FLT_exp.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
930
refine (_ (Zdigits_mult_ge radix2 (Zpos mx) (Zpos my) _ _)) ; try discriminate.
931
refine (_ (Zdigits_gt_0 radix2 (Zpos mx) _) (Zdigits_gt_0 radix2 (Zpos my) _)) ; try discriminate.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
932
generalize (Zdigits radix2 (Zpos mx)) (Zdigits radix2 (Zpos my)) (Zdigits radix2 (Zpos mx * Zpos my)).
933
clear -Hmax.
934
unfold emin.
935 936 937 938 939 940 941 942 943 944 945 946 947 948
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.
949

950 951
Definition Bmult mult_nan m x y :=
  let f pl := B754_nan (fst pl) (snd pl) in
952
  match x, y with
953
  | B754_nan _ _, _ | _, B754_nan _ _ => f (mult_nan x y)
954 955 956
  | 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)
957 958
  | B754_infinity _, B754_zero _ => f (mult_nan x y)
  | B754_zero _, B754_infinity _ => f (mult_nan x y)
959 960 961 962 963 964 965 966
  | 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 :
967
  forall mult_nan m x y,
968
  if Rlt_bool (Rabs (round radix2 fexp (round_mode m) (B2R x * B2R y))) (bpow radix2 emax) then
969 970
    B2R (Bmult mult_nan m x y) = round radix2 fexp (round_mode m) (B2R x * B2R y) /\
    is_finite (Bmult mult_nan m x y) = andb (is_finite x) (is_finite y)
971
  else
972
    B2FF (Bmult mult_nan m x y) = binary_overflow m (xorb (Bsign x) (Bsign y)).
973
Proof.
974
intros mult_nan m [sx|sx|sx plx|sx mx ex Hx] [sy|sy|sy ply|sy my ey Hy] ;
975
  try ( rewrite ?Rmult_0_r, ?Rmult_0_l, round_0, Rabs_R0, Rlt_bool_true ; [ split ; apply refl_equal | apply bpow_gt_0 | auto with typeclass_instances ] ).
976 977
simpl.
case Bmult_correct_aux.
978 979 980 981
intros H1.
case Rlt_bool.
intros (H2, H3).
split.
982
now rewrite B2R_FF2B.
983 984
now rewrite is_finite_FF2B.
intros H2.
985 986 987
now rewrite B2FF_FF2B.
Qed.

988 989
Definition Bmult_FF mult_nan m x y :=
  let f pl := F754_nan (fst pl) (snd pl) in
990
  match x, y with
991
  | F754_nan _ _, _ | _, F754_nan _ _ => f (mult_nan x y)
992 993 994
  | F754_infinity sx, F754_infinity sy => F754_infinity (xorb sx sy)
  | F754_infinity sx, F754_finite sy _ _ => F754_infinity (xorb sx sy)
  | F754_finite sx _ _, F754_infinity sy => F754_infinity (xorb sx sy)
995 996
  | F754_infinity _, F754_zero _ => f (mult_nan x y)
  | F754_zero _, F754_infinity _ => f (mult_nan x y)