Bits.v 18.4 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) 2011-2013 Sylvie Boldo
6
#<br />#
7
Copyright (C) 2011-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 encoding of binary floating-point data *)
Guillaume Melquiond's avatar
Guillaume Melquiond committed
21
Require Import Core Digits Binary.
22 23 24

Section Binary_Bits.

25 26 27 28 29
Arguments exist {A} {P}.
Arguments B754_zero {prec} {emax}.
Arguments B754_infinity {prec} {emax}.
Arguments B754_nan {prec} {emax}.
Arguments B754_finite {prec} {emax}.
30

BOLDO Sylvie's avatar
BOLDO Sylvie committed
31
(** Number of bits for the fraction and exponent *)
32 33 34 35 36 37 38 39 40 41
Variable mw ew : Z.
Hypothesis Hmw : (0 < mw)%Z.
Hypothesis Hew : (0 < ew)%Z.

Let emax := Zpower 2 (ew - 1).
Let prec := (mw + 1)%Z.
Let emin := (3 - emax - prec)%Z.
Let binary_float := binary_float prec emax.

Let Hprec : (0 < prec)%Z.
42
Proof.
43 44 45 46 47 48
unfold prec.
apply Zle_lt_succ.
now apply Zlt_le_weak.
Qed.

Let Hm_gt_0 : (0 < 2^mw)%Z.
49
Proof.
50 51 52 53 54
apply (Zpower_gt_0 radix2).
now apply Zlt_le_weak.
Qed.

Let He_gt_0 : (0 < 2^ew)%Z.
55
Proof.
56 57 58 59 60 61 62
apply (Zpower_gt_0 radix2).
now apply Zlt_le_weak.
Qed.

Hypothesis Hmax : (prec < emax)%Z.

Definition join_bits (s : bool) m e :=
63
  (Z.shiftl ((if s then Zpower 2 ew else 0) + e) mw + m)%Z.
64

65 66 67 68 69 70 71 72
Lemma join_bits_range :
  forall s m e,
  (0 <= m < 2^mw)%Z ->
  (0 <= e < 2^ew)%Z ->
  (0 <= join_bits s m e < 2 ^ (mw + ew + 1))%Z.
Proof.
intros s m e Hm He.
unfold join_bits.
73
rewrite Z.shiftl_mul_pow2 by now apply Zlt_le_weak.
74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97
split.
- apply (Zplus_le_compat 0 _ 0) with (2 := proj1 Hm).
  rewrite <- (Zmult_0_l (2^mw)).
  apply Zmult_le_compat_r.
  case s.
  clear -He ; omega.
  now rewrite Zmult_0_l.
  clear -Hm ; omega.
- apply Zlt_le_trans with (((if s then 2 ^ ew else 0) + e + 1) * 2 ^ mw)%Z.
  rewrite (Zmult_plus_distr_l _ 1).
  apply Zplus_lt_compat_l.
  now rewrite Zmult_1_l.
  rewrite <- (Zplus_assoc mw), (Zplus_comm mw), Zpower_plus.
  apply Zmult_le_compat_r.
  rewrite Zpower_plus.
  change (2^1)%Z with 2%Z.
  case s ; clear -He ; omega.
  now apply Zlt_le_weak.
  easy.
  clear -Hm ; omega.
  clear -Hew ; omega.
  now apply Zlt_le_weak.
Qed.

98 99 100 101 102 103 104 105 106 107 108 109 110
Definition split_bits x :=
  let mm := Zpower 2 mw in
  let em := Zpower 2 ew in
  (Zle_bool (mm * em) x, Zmod x mm, Zmod (Zdiv x mm) em)%Z.

Theorem split_join_bits :
  forall s m e,
  (0 <= m < Zpower 2 mw)%Z ->
  (0 <= e < Zpower 2 ew)%Z ->
  split_bits (join_bits s m e) = (s, m, e).
Proof.
intros s m e Hm He.
unfold split_bits, join_bits.
111
rewrite Z.shiftl_mul_pow2 by now apply Zlt_le_weak.
112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156
apply f_equal2.
apply f_equal2.
(* *)
case s.
apply Zle_bool_true.
apply Zle_0_minus_le.
ring_simplify.
apply Zplus_le_0_compat.
apply Zmult_le_0_compat.
apply He.
now apply Zlt_le_weak.
apply Hm.
apply Zle_bool_false.
apply Zplus_lt_reg_l with (2^mw * (-e))%Z.
replace (2 ^ mw * - e + ((0 + e) * 2 ^ mw + m))%Z with (m * 1)%Z by ring.
rewrite <- Zmult_plus_distr_r.
apply Zlt_le_trans with (2^mw * 1)%Z.
now apply Zmult_lt_compat_r.
apply Zmult_le_compat_l.
clear -He. omega.
now apply Zlt_le_weak.
(* *)
rewrite Zplus_comm.
rewrite Z_mod_plus_full.
now apply Zmod_small.
(* *)
rewrite Z_div_plus_full_l.
rewrite Zdiv_small with (1 := Hm).
rewrite Zplus_0_r.
case s.
replace (2^ew + e)%Z with (e + 1 * 2^ew)%Z by ring.
rewrite Z_mod_plus_full.
now apply Zmod_small.
now apply Zmod_small.
now apply Zgt_not_eq.
Qed.

Theorem join_split_bits :
  forall x,
  (0 <= x < Zpower 2 (mw + ew + 1))%Z ->
  let '(s, m, e) := split_bits x in
  join_bits s m e = x.
Proof.
intros x Hx.
unfold split_bits, join_bits.
157
rewrite Z.shiftl_mul_pow2 by now apply Zlt_le_weak.
158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215
pattern x at 4 ; rewrite Z_div_mod_eq_full with x (2^mw)%Z.
apply (f_equal (fun v => (v + _)%Z)).
rewrite Zmult_comm.
apply f_equal.
pattern (x / (2^mw))%Z at 2 ; rewrite Z_div_mod_eq_full with (x / (2^mw))%Z (2^ew)%Z.
apply (f_equal (fun v => (v + _)%Z)).
replace (x / 2 ^ mw / 2 ^ ew)%Z with (if Zle_bool (2 ^ mw * 2 ^ ew) x then 1 else 0)%Z.
case Zle_bool.
now rewrite Zmult_1_r.
now rewrite Zmult_0_r.
rewrite Zdiv_Zdiv.
apply sym_eq.
case Zle_bool_spec ; intros Hs.
apply Zle_antisym.
cut (x / (2^mw * 2^ew) < 2)%Z. clear ; omega.
apply Zdiv_lt_upper_bound.
try apply Hx. (* 8.2/8.3 compatibility *)
now apply Zmult_lt_0_compat.
rewrite <- Zpower_exp ; try ( apply Zle_ge ; apply Zlt_le_weak ; assumption ).
change 2%Z at 1 with (Zpower 2 1).
rewrite <- Zpower_exp.
now rewrite Zplus_comm.
discriminate.
apply Zle_ge.
now apply Zplus_le_0_compat ; apply Zlt_le_weak.
apply Zdiv_le_lower_bound.
try apply Hx. (* 8.2/8.3 compatibility *)
now apply Zmult_lt_0_compat.
now rewrite Zmult_1_l.
apply Zdiv_small.
now split.
now apply Zlt_le_weak.
now apply Zlt_le_weak.
now apply Zgt_not_eq.
now apply Zgt_not_eq.
Qed.

Theorem split_bits_inj :
  forall x y,
  (0 <= x < Zpower 2 (mw + ew + 1))%Z ->
  (0 <= y < Zpower 2 (mw + ew + 1))%Z ->
  split_bits x = split_bits y ->
  x = y.
Proof.
intros x y Hx Hy.
generalize (join_split_bits x Hx) (join_split_bits y Hy).
destruct (split_bits x) as ((sx, mx), ex).
destruct (split_bits y) as ((sy, my), ey).
intros Jx Jy H. revert Jx Jy.
inversion_clear H.
intros Jx Jy.
now rewrite <- Jx.
Qed.

Definition bits_of_binary_float (x : binary_float) :=
  match x with
  | B754_zero sx => join_bits sx 0 0
  | B754_infinity sx => join_bits sx 0 (Zpower 2 ew - 1)
216
  | B754_nan sx plx _ => join_bits sx (Zpos plx) (Zpower 2 ew - 1)
217
  | B754_finite sx mx ex _ =>
218 219 220
    let m := (Zpos mx - Zpower 2 mw)%Z in
    if Zle_bool 0 m then
      join_bits sx m (ex - emin + 1)
221 222 223 224 225 226 227 228
    else
      join_bits sx (Zpos mx) 0
  end.

Definition split_bits_of_binary_float (x : binary_float) :=
  match x with
  | B754_zero sx => (sx, 0, 0)%Z
  | B754_infinity sx => (sx, 0, Zpower 2 ew - 1)%Z
229
  | B754_nan sx plx _ => (sx, Zpos plx, Zpower 2 ew - 1)%Z
230
  | B754_finite sx mx ex _ =>
231 232 233
    let m := (Zpos mx - Zpower 2 mw)%Z in
    if Zle_bool 0 m then
      (sx, m, ex - emin + 1)%Z
234 235 236 237 238 239 240 241
    else
      (sx, Zpos mx, 0)%Z
  end.

Theorem split_bits_of_binary_float_correct :
  forall x,
  split_bits (bits_of_binary_float x) = split_bits_of_binary_float x.
Proof.
242
intros [sx|sx|sx plx Hplx|sx mx ex Hx] ;
243
  try ( simpl ; apply split_join_bits ; split ; try apply Zle_refl ; try apply Zlt_pred ; trivial ; omega ).
244 245
simpl. apply split_join_bits; split; try (zify; omega).
destruct (digits2_Pnat_correct plx).
246
unfold nan_pl in Hplx.
247
rewrite Zpos_digits2_pos, <- Z_of_nat_S_digits2_Pnat in Hplx.
248 249 250 251 252 253
rewrite Zpower_nat_Z in H0.
eapply Zlt_le_trans. apply H0.
change 2%Z with (radix_val radix2). apply Zpower_le.
rewrite Z.ltb_lt in Hplx.
unfold prec in *. zify; omega.
(* *)
254
unfold bits_of_binary_float, split_bits_of_binary_float.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
255
assert (Hf: (emin <= ex /\ Zdigits radix2 (Zpos mx) <= prec)%Z).
256
destruct (andb_prop _ _ Hx) as (Hx', _).
257
unfold canonical_mantissa in Hx'.
258
rewrite Zpos_digits2_pos in Hx'.
259 260 261 262
generalize (Zeq_bool_eq _ _ Hx').
unfold FLT_exp.
unfold emin.
clear ; zify ; omega.
263 264
case Zle_bool_spec ; intros H ;
  [ apply -> Z.le_0_sub in H | apply -> Z.lt_sub_0 in H ] ;
265 266 267 268 269 270
  apply split_join_bits ; try now split.
(* *)
split.
clear -He_gt_0 H ; omega.
cut (Zpos mx < 2 * 2^mw)%Z. clear ; omega.
replace (2 * 2^mw)%Z with (2^prec)%Z.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
271
apply (Zpower_gt_Zdigits radix2 _ (Zpos mx)).
272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295
apply Hf.
unfold prec.
rewrite Zplus_comm.
apply Zpower_exp ; apply Zle_ge.
discriminate.
now apply Zlt_le_weak.
(* *)
split.
generalize (proj1 Hf).
clear ; omega.
destruct (andb_prop _ _ Hx) as (_, Hx').
unfold emin.
replace (2^ew)%Z with (2 * emax)%Z.
generalize (Zle_bool_imp_le _ _ Hx').
clear ; omega.
apply sym_eq.
rewrite (Zsucc_pred ew).
unfold Zsucc.
rewrite Zplus_comm.
apply Zpower_exp ; apply Zle_ge.
discriminate.
now apply Zlt_0_le_0_pred.
Qed.

296 297 298
Theorem bits_of_binary_float_range:
  forall x, (0 <= bits_of_binary_float x < 2^(mw+ew+1))%Z.
Proof.
299
unfold bits_of_binary_float.
300
intros [sx|sx|sx pl pl_range|sx mx ex H].
301 302 303 304 305 306 307 308 309 310 311 312 313 314
- apply join_bits_range ; now split.
- apply join_bits_range.
  now split.
  clear -He_gt_0 ; omega.
- apply Z.ltb_lt in pl_range.
  apply join_bits_range.
  split.
  easy.
  apply (Zpower_gt_Zdigits radix2 _ (Zpos pl)).
  apply Z.lt_succ_r.
  now rewrite <- Zdigits2_Zdigits.
  clear -He_gt_0 ; omega.
- unfold bounded in H.
  apply Bool.andb_true_iff in H ; destruct H as [A B].
315
  apply Z.leb_le in B.
316
  unfold canonical_mantissa, FLT_exp in A. apply Zeq_bool_eq in A.
317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335
  case Zle_bool_spec ; intros H.
  + apply join_bits_range.
    * split.
      clear -H ; omega.
      rewrite Zpos_digits2_pos in A.
      cut (Zpos mx < 2 ^ prec)%Z.
      unfold prec.
      rewrite Zpower_plus by (clear -Hmw ; omega).
      change (2^1)%Z with 2%Z.
      clear ; omega.
      apply (Zpower_gt_Zdigits radix2 _ (Zpos mx)).
      clear -A ; zify ; omega.
    * split.
      unfold emin ; clear -A ; zify ; omega.
      replace ew with ((ew - 1) + 1)%Z by ring.
      rewrite Zpower_plus by (clear - Hew ; omega).
      unfold emin, emax in *.
      change (2^1)%Z with 2%Z.
      clear -B ; omega.
336 337
  + apply -> Z.lt_sub_0 in H.
    apply join_bits_range ; now split.
338 339
Qed.

340 341 342 343 344 345
Definition binary_float_of_bits_aux x :=
  let '(sx, mx, ex) := split_bits x in
  if Zeq_bool ex 0 then
    match mx with
    | Z0 => F754_zero sx
    | Zpos px => F754_finite sx px emin
346
    | Zneg _ => F754_nan false xH (* dummy *)
347 348
    end
  else if Zeq_bool ex (Zpower 2 ew - 1) then
349 350 351 352 353
    match mx with
      | Z0 => F754_infinity sx
      | Zpos plx => F754_nan sx plx
      | Zneg _ => F754_nan false xH (* dummy *)
    end
354 355 356
  else
    match (mx + Zpower 2 mw)%Z with
    | Zpos px => F754_finite sx px (ex + emin - 1)
357
    | _ => F754_nan false xH (* dummy *)
358 359 360 361 362 363 364 365
    end.

Lemma binary_float_of_bits_aux_correct :
  forall x,
  valid_binary prec emax (binary_float_of_bits_aux x) = true.
Proof.
intros x.
unfold binary_float_of_bits_aux, split_bits.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
366 367 368 369
assert (Hnan: nan_pl prec 1 = true).
  apply Z.ltb_lt.
  simpl. unfold prec.
  clear -Hmw ; omega.
370 371 372 373
case Zeq_bool_spec ; intros He1.
case_eq (x mod 2^mw)%Z ; try easy.
(* subnormal *)
intros px Hm.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
374 375
assert (Zdigits radix2 (Zpos px) <= mw)%Z.
apply Zdigits_le_Zpower.
376 377 378 379
simpl.
rewrite <- Hm.
eapply Z_mod_lt.
now apply Zlt_gt.
380 381
apply bounded_canonical_lt_emax ; try assumption.
unfold canonical, cexp.
382
fold emin.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
383
rewrite mag_F2R_Zdigits. 2: discriminate.
384 385 386 387 388 389 390
unfold Fexp, FLT_exp.
apply sym_eq.
apply Zmax_right.
clear -H Hprec.
unfold prec ; omega.
apply Rnot_le_lt.
intros H0.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
391 392 393
refine (_ (mag_le radix2 _ _ _ H0)).
rewrite mag_bpow.
rewrite mag_F2R_Zdigits. 2: discriminate.
394 395 396 397 398 399 400
unfold emin, prec.
apply Zlt_not_le.
cut (0 < emax)%Z. clear -H Hew ; omega.
apply (Zpower_gt_0 radix2).
clear -Hew ; omega.
apply bpow_gt_0.
case Zeq_bool_spec ; intros He2.
401 402 403
case_eq (x mod 2 ^ mw)%Z; try easy.
(* nan *)
intros plx Eqplx. apply Z.ltb_lt.
404
rewrite Zpos_digits2_pos.
405 406 407 408 409
assert (forall a b, a <= b -> a < b+1)%Z by (intros; omega). apply H. clear H.
apply Zdigits_le_Zpower. simpl.
rewrite <- Eqplx. edestruct Z_mod_lt; eauto.
change 2%Z with (radix_val radix2).
apply Z.lt_gt, Zpower_gt_0. omega.
410 411 412
case_eq (x mod 2^mw + 2^mw)%Z ; try easy.
(* normal *)
intros px Hm.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
413
assert (prec = Zdigits radix2 (Zpos px)).
414
(* . *)
Guillaume Melquiond's avatar
Guillaume Melquiond committed
415
rewrite Zdigits_mag. 2: discriminate.
416
apply sym_eq.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
417
apply mag_unique.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
418
rewrite <- abs_IZR.
419 420
unfold Zabs.
replace (prec - 1)%Z with mw by ( unfold prec ; ring ).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
421 422
rewrite <- IZR_Zpower with (1 := Zlt_le_weak _ _ Hmw).
rewrite <- IZR_Zpower. 2: now apply Zlt_le_weak.
423 424
rewrite <- Hm.
split.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
425
apply IZR_le.
426 427 428 429
change (radix2^mw)%Z with (0 + 2^mw)%Z.
apply Zplus_le_compat_r.
eapply Z_mod_lt.
now apply Zlt_gt.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
430
apply IZR_lt.
431 432 433 434 435 436 437
unfold prec.
rewrite Zpower_exp. 2: now apply Zle_ge ; apply Zlt_le_weak. 2: discriminate.
rewrite <- Zplus_diag_eq_mult_2.
apply Zplus_lt_compat_r.
eapply Z_mod_lt.
now apply Zlt_gt.
(* . *)
438 439
apply bounded_canonical_lt_emax ; try assumption.
unfold canonical, cexp.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
440
rewrite mag_F2R_Zdigits. 2: discriminate.
441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457
unfold Fexp, FLT_exp.
rewrite <- H.
set (ex := ((x / 2^mw) mod 2^ew)%Z).
replace (prec + (ex + emin - 1) - prec)%Z with (ex + emin - 1)%Z by ring.
apply sym_eq.
apply Zmax_left.
revert He1.
fold ex.
cut (0 <= ex)%Z.
unfold emin.
clear ; intros H1 H2 ; omega.
eapply Z_mod_lt.
apply Zlt_gt.
apply (Zpower_gt_0 radix2).
now apply Zlt_le_weak.
apply Rnot_le_lt.
intros H0.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
458 459 460
refine (_ (mag_le radix2 _ _ _ H0)).
rewrite mag_bpow.
rewrite mag_F2R_Zdigits. 2: discriminate.
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
rewrite <- H.
apply Zlt_not_le.
unfold emin.
apply Zplus_lt_reg_r with (emax - 1)%Z.
ring_simplify.
revert He2.
set (ex := ((x / 2^mw) mod 2^ew)%Z).
cut (ex < 2^ew)%Z.
replace (2^ew)%Z with (2 * emax)%Z.
clear ; intros H1 H2 ; omega.
replace ew with (1 + (ew - 1))%Z by ring.
rewrite Zpower_exp.
apply refl_equal.
discriminate.
clear -Hew ; omega.
eapply Z_mod_lt.
apply Zlt_gt.
apply (Zpower_gt_0 radix2).
now apply Zlt_le_weak.
apply bpow_gt_0.
Qed.

Definition binary_float_of_bits x :=
  FF2B prec emax _ (binary_float_of_bits_aux_correct x).

Theorem binary_float_of_bits_of_binary_float :
  forall x,
  binary_float_of_bits (bits_of_binary_float x) = x.
Proof.
intros x.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
491
apply B2FF_inj.
492 493 494 495
unfold binary_float_of_bits.
rewrite B2FF_FF2B.
unfold binary_float_of_bits_aux.
rewrite split_bits_of_binary_float_correct.
496
destruct x as [sx|sx|sx plx Hplx|sx mx ex Bx].
497 498 499 500 501 502 503 504 505 506
apply refl_equal.
(* *)
simpl.
rewrite Zeq_bool_false.
now rewrite Zeq_bool_true.
cut (1 < 2^ew)%Z. clear ; omega.
now apply (Zpower_gt_1 radix2).
(* *)
simpl.
rewrite Zeq_bool_false.
507
rewrite Zeq_bool_true; auto.
508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528
cut (1 < 2^ew)%Z. clear ; omega.
now apply (Zpower_gt_1 radix2).
(* *)
unfold split_bits_of_binary_float.
case Zle_bool_spec ; intros Hm.
(* . *)
rewrite Zeq_bool_false.
rewrite Zeq_bool_false.
now ring_simplify (Zpos mx - 2 ^ mw + 2 ^ mw)%Z (ex - emin + 1 + emin - 1)%Z.
destruct (andb_prop _ _ Bx) as (_, H1).
generalize (Zle_bool_imp_le _ _ H1).
unfold emin.
replace (2^ew)%Z with (2 * emax)%Z.
clear ; omega.
replace ew with (1 + (ew - 1))%Z by ring.
rewrite Zpower_exp.
apply refl_equal.
discriminate.
clear -Hew ; omega.
destruct (andb_prop _ _ Bx) as (H1, _).
generalize (Zeq_bool_eq _ _ H1).
529
rewrite Zpos_digits2_pos.
530
unfold FLT_exp, emin.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
531
generalize (Zdigits radix2 (Zpos mx)).
532 533 534 535 536 537 538 539
clear.
intros ; zify ; omega.
(* . *)
rewrite Zeq_bool_true. 2: apply refl_equal.
simpl.
apply f_equal.
destruct (andb_prop _ _ Bx) as (H1, _).
generalize (Zeq_bool_eq _ _ H1).
540
rewrite Zpos_digits2_pos.
541
unfold FLT_exp, emin, prec.
542
apply -> Z.lt_sub_0 in Hm.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
543 544
generalize (Zdigits_le_Zpower radix2 _ (Zpos mx) Hm).
generalize (Zdigits radix2 (Zpos mx)).
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
clear.
intros ; zify ; omega.
Qed.

Theorem bits_of_binary_float_of_bits :
  forall x,
  (0 <= x < 2^(mw+ew+1))%Z ->
  bits_of_binary_float (binary_float_of_bits x) = x.
Proof.
intros x Hx.
unfold binary_float_of_bits, bits_of_binary_float.
set (Cx := binary_float_of_bits_aux_correct x).
clearbody Cx.
rewrite match_FF2B.
revert Cx.
generalize (join_split_bits x Hx).
unfold binary_float_of_bits_aux.
case_eq (split_bits x).
intros (sx, mx) ex Sx.
assert (Bm: (0 <= mx < 2^mw)%Z).
inversion_clear Sx.
apply Z_mod_lt.
now apply Zlt_gt.
case Zeq_bool_spec ; intros He1.
(* subnormal *)
case_eq mx.
571
intros Hm Jx _.
572
now rewrite He1 in Jx.
573
intros px Hm Jx _.
574 575
rewrite Zle_bool_false.
now rewrite <- He1.
576
apply <- Z.lt_sub_0.
577
now rewrite <- Hm.
578
intros px Hm _ _.
579 580 581 582 583
apply False_ind.
apply Zle_not_lt with (1 := proj1 Bm).
now rewrite Hm.
case Zeq_bool_spec ; intros He2.
(* infinity/nan *)
584 585 586 587
case_eq mx; intros Hm.
now rewrite He2.
now rewrite He2.
intros. zify; omega.
588 589 590 591 592
(* normal *)
case_eq (mx + 2 ^ mw)%Z.
intros Hm.
apply False_ind.
clear -Bm Hm ; omega.
593
intros p Hm Jx Cx.
594 595 596
rewrite <- Hm.
rewrite Zle_bool_true.
now ring_simplify (mx + 2^mw - 2^mw)%Z (ex + emin - 1 - emin + 1)%Z.
597
now ring_simplify.
598 599 600 601 602 603 604
intros p Hm.
apply False_ind.
clear -Bm Hm ; zify ; omega.
Qed.

End Binary_Bits.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
605
(** Specialization for IEEE single precision operations *)
606 607
Section B32_Bits.

608
Arguments B754_nan {prec} {emax}.
609

610 611 612 613 614 615 616 617 618 619
Definition binary32 := binary_float 24 128.

Let Hprec : (0 < 24)%Z.
apply refl_equal.
Qed.

Let Hprec_emax : (24 < 128)%Z.
apply refl_equal.
Qed.

620 621
Definition default_nan_pl32 : { nan | is_nan 24 128 nan = true } :=
  exist _ (@B754_nan 24 128 false (iter_nat xO 22 xH) (refl_equal true)) (refl_equal true).
622

623
Definition unop_nan_pl32 (f : binary32) :=
624 625
  match f as f with
  | B754_nan s pl Hpl => exist _ (B754_nan s pl Hpl) (refl_equal true)
626 627 628
  | _ => default_nan_pl32
  end.

629
Definition binop_nan_pl32 (f1 f2 : binary32) :=
630
  match f1, f2 with
631 632
  | B754_nan s1 pl1 Hpl1, _ => exist _ (B754_nan s1 pl1 Hpl1) (refl_equal true)
  | _, B754_nan s2 pl2 Hpl2 => exist _ (B754_nan s2 pl2 Hpl2) (refl_equal true)
633 634 635
  | _, _ => default_nan_pl32
  end.

636
Definition b32_opp := Bopp 24 128 unop_nan_pl32.
637 638 639 640 641
Definition b32_plus := Bplus _ _ Hprec Hprec_emax binop_nan_pl32.
Definition b32_minus := Bminus _ _ Hprec Hprec_emax binop_nan_pl32.
Definition b32_mult := Bmult _ _ Hprec Hprec_emax binop_nan_pl32.
Definition b32_div := Bdiv _ _ Hprec Hprec_emax binop_nan_pl32.
Definition b32_sqrt := Bsqrt _ _ Hprec Hprec_emax unop_nan_pl32.
642 643 644 645 646 647

Definition b32_of_bits : Z -> binary32 := binary_float_of_bits 23 8 (refl_equal _) (refl_equal _) (refl_equal _).
Definition bits_of_b32 : binary32 -> Z := bits_of_binary_float 23 8.

End B32_Bits.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
648
(** Specialization for IEEE double precision operations *)
649 650
Section B64_Bits.

651
Arguments B754_nan {prec} {emax}.
652

653 654 655 656 657 658 659 660 661 662
Definition binary64 := binary_float 53 1024.

Let Hprec : (0 < 53)%Z.
apply refl_equal.
Qed.

Let Hprec_emax : (53 < 1024)%Z.
apply refl_equal.
Qed.

663 664
Definition default_nan_pl64 : { nan | is_nan 53 1024 nan = true } :=
  exist _ (@B754_nan 53 1024 false (iter_nat xO 51 xH) (refl_equal true)) (refl_equal true).
665

666
Definition unop_nan_pl64 (f : binary64) :=
667 668
  match f as f with
  | B754_nan s pl Hpl => exist _ (B754_nan s pl Hpl) (refl_equal true)
669 670 671
  | _ => default_nan_pl64
  end.

672 673 674 675
Definition binop_nan_pl64 (f1 f2 : binary64) :=
  match f1, f2 with
  | B754_nan s1 pl1 Hpl1, _ => exist _ (B754_nan s1 pl1 Hpl1) (refl_equal true)
  | _, B754_nan s2 pl2 Hpl2 => exist _ (B754_nan s2 pl2 Hpl2) (refl_equal true)
676 677 678
  | _, _ => default_nan_pl64
  end.

679
Definition b64_opp := Bopp 53 1024 unop_nan_pl64.
680 681 682 683 684
Definition b64_plus := Bplus _ _ Hprec Hprec_emax binop_nan_pl64.
Definition b64_minus := Bminus _ _ Hprec Hprec_emax binop_nan_pl64.
Definition b64_mult := Bmult _ _ Hprec Hprec_emax binop_nan_pl64.
Definition b64_div := Bdiv _ _ Hprec Hprec_emax binop_nan_pl64.
Definition b64_sqrt := Bsqrt _ _ Hprec Hprec_emax unop_nan_pl64.
685 686 687 688 689

Definition b64_of_bits : Z -> binary64 := binary_float_of_bits 52 11 (refl_equal _) (refl_equal _) (refl_equal _).
Definition bits_of_b64 : binary64 -> Z := bits_of_binary_float 52 11.

End B64_Bits.