Merge of Flocq version 2.2.0.

More precise modeling of NaNs.


git-svn-id: https://yquem.inria.fr/compcert/svn/compcert/trunk@2303 fca1b0fc-160b-0410-b1d3-a4f43f01ea2e

1 files changed, 339 insertions, 105 deletions

diff --git a/lib/Floats.v b/lib/Floats.v
index 94c19d2..0151cf0 100644
--- a/lib/Floats.v
+++ b/lib/Floats.v
@@ -38,20 +38,69 @@ Definition zero: float := B754_zero _ _ false. (**r the float [+0.0] *)
 Definition eq_dec: forall (f1 f2: float), {f1 = f2} + {f1 <> f2}.
 Proof.
   Ltac try_not_eq := try solve [right; congruence].
-  destruct f1, f2;
+  destruct f1 as [| |? []|], f2 as [| |? []|];
   try destruct b; try destruct b0;
-  try solve [left; auto]; try_not_eq;
+  try solve [left; auto]; try_not_eq.
+  destruct (positive_eq_dec x x0); try_not_eq;
+    subst; left; f_equal; f_equal; apply proof_irr.
+  destruct (positive_eq_dec x x0); try_not_eq;
+    subst; left; f_equal; f_equal; apply proof_irr.
+  destruct (positive_eq_dec m m0); try_not_eq;
+  destruct (Z_eq_dec e e1); try solve [right; intro H; inv H; congruence];
+  subst; left; rewrite (proof_irr e0 e2); auto.
   destruct (positive_eq_dec m m0); try_not_eq;
   destruct (Z_eq_dec e e1); try solve [right; intro H; inv H; congruence];
   subst; left; rewrite (proof_irr e0 e2); auto.
 Defined.
 
+(* Transform a Nan payload to a quiet Nan payload.
+   This is not part of the IEEE754 standard, but shared between all
+   architectures of Compcert. *)
+Program Definition transform_quiet_pl (pl:nan_pl 53) : nan_pl 53 :=
+  Pos.lor pl (nat_iter 51 xO xH).
+Next Obligation.
+  destruct pl.
+  simpl. rewrite Z.ltb_lt in *.
+  assert (forall x, S (Fcore_digits.digits2_Pnat x) = Pos.to_nat (Pos.size x)).
+  { induction x0; simpl; auto; rewrite IHx0; zify; omega. }
+  fold (Z.of_nat (S (Fcore_digits.digits2_Pnat (Pos.lor x 2251799813685248)))).
+  rewrite H, positive_nat_Z, Psize_log_inf, <- Zlog2_log_inf in *. clear H.
+  change (Z.pos (Pos.lor x 2251799813685248)) with (Z.lor (Z.pos x) 2251799813685248%Z).
+  rewrite Z.log2_lor by (zify; omega).
+  apply Z.max_case. auto. simpl. omega.
+Qed.
+
+Lemma nan_payload_fequal:
+  forall prec p1 e1 p2 e2, p1 = p2 -> (exist _ p1 e1:nan_pl prec) = exist _ p2 e2.
+Proof.
+  simpl; intros; subst. f_equal. apply Fcore_Zaux.eqbool_irrelevance.
+Qed.
+
+Lemma lor_idempotent:
+  forall x y, Pos.lor (Pos.lor x y) y = Pos.lor x y.
+Proof.
+  induction x; destruct y; simpl; f_equal; auto;
+  induction y; simpl; f_equal; auto.
+Qed.
+
+Lemma transform_quiet_pl_idempotent:
+  forall pl, transform_quiet_pl (transform_quiet_pl pl) = transform_quiet_pl pl.
+Proof.
+  intros []; simpl; intros. apply nan_payload_fequal.
+  simpl. apply lor_idempotent.
+Qed.
+
 (** Arithmetic operations *)
 
-Definition neg: float -> float := b64_opp. (**r opposite (change sign) *)
+(* The Nan payload operations for neg and abs is not part of the IEEE754
+   standard, but shared between all architectures of Compcert. *)
+Definition neg_pl (s:bool) (pl:nan_pl 53) := (negb s, pl).
+Definition abs_pl (s:bool) (pl:nan_pl 53) := (false, pl).
+
+Definition neg: float -> float := b64_opp neg_pl. (**r opposite (change sign) *)
 Definition abs (x: float): float := (**r absolute value (set sign to [+]) *)
   match x with
-  | B754_nan => x
+  | B754_nan s pl => let '(s, pl) := abs_pl s pl in B754_nan _ _ s pl
   | B754_infinity _ => B754_infinity _ _ false
   | B754_finite _ m e H => B754_finite _ _ false m e H
   | B754_zero _ => B754_zero _ _ false
@@ -71,9 +120,30 @@ Definition binary_normalize32_correct (m e:Z) (s:bool) :=
   binary_normalize_correct 24 128 eq_refl eq_refl mode_NE m e s.
 Global Opaque binary_normalize32_correct.
 
+(* The Nan payload operations for single <-> double conversions are not part of
+   the IEEE754 standard, but shared between all architectures of Compcert. *)
+Definition floatofbinary32_pl (s:bool) (pl:nan_pl 24) : (bool * nan_pl 53).
+  refine (s, transform_quiet_pl (exist _ (Pos.shiftl_nat (proj1_sig pl) 29) _)).
+  abstract (
+    destruct pl; unfold proj1_sig, Pos.shiftl_nat, nat_iter, Fcore_digits.digits2_Pnat;
+    fold (Fcore_digits.digits2_Pnat x);
+    rewrite Z.ltb_lt in *;
+    zify; omega).
+Defined.
+
+Definition binary32offloat_pl (s:bool) (pl:nan_pl 53) : (bool * nan_pl 24).
+  refine (s, exist _ (Pos.shiftr_nat (proj1_sig (transform_quiet_pl pl)) 29) _).
+  abstract (
+    destruct (transform_quiet_pl pl); unfold proj1_sig, Pos.shiftr_nat, nat_iter;
+    rewrite Z.ltb_lt in *;
+    assert (forall x, Fcore_digits.digits2_Pnat (Pos.div2 x) =
+                      (Fcore_digits.digits2_Pnat x - 1)%nat) by (destruct x0; simpl; zify; omega);
+    rewrite !H, <- !NPeano.Nat.sub_add_distr; zify; omega).
+Defined.
+
 Definition floatofbinary32 (f: binary32) : float := (**r single precision embedding in double precision *)
   match f with
-    | B754_nan => B754_nan _ _
+    | B754_nan s pl => let '(s, pl) := floatofbinary32_pl s pl in B754_nan _ _ s pl
     | B754_infinity s => B754_infinity _ _ s
     | B754_zero s => B754_zero _ _ s
     | B754_finite s m e _ =>
@@ -82,7 +152,7 @@ Definition floatofbinary32 (f: binary32) : float := (**r single precision embedd
 
 Definition binary32offloat (f: float) : binary32 := (**r conversion to single precision *)
   match f with
-    | B754_nan => B754_nan _ _
+    | B754_nan s pl => let '(s, pl) := binary32offloat_pl s pl in B754_nan _ _ s pl
     | B754_infinity s => B754_infinity _ _ s
     | B754_zero s => B754_zero _ _ s
     | B754_finite s m e _ =>
@@ -190,15 +260,33 @@ Definition singleoflong (n:int64): float := (**r conversion from signed 64-bit i
 Definition singleoflongu (n:int64): float:= (**r conversion from unsigned 64-bit int to single-precision float *)
   floatofbinary32 (binary_normalize32 (Int64.unsigned n) 0 false).
 
-Definition add: float -> float -> float := b64_plus mode_NE. (**r addition *)
-Definition sub: float -> float -> float := b64_minus mode_NE. (**r subtraction *)
-Definition mul: float -> float -> float := b64_mult mode_NE. (**r multiplication *)
-Definition div: float -> float -> float := b64_div mode_NE. (**r division *)
+Parameter binop_pl: binary64 -> binary64 -> bool*nan_pl 53.
+
+Definition binop_propagate1_prop binop_pl :=
+  forall f1 f2:float, is_nan _ _ f1 = true -> is_nan _ _ f2 = false ->
+    binop_pl f1 f1 = (binop_pl f1 f2 : bool*nan_pl 53).
+
+
+(* Proved in Nan.v for different architectures *)
+Hypothesis binop_propagate1: binop_propagate1_prop binop_pl.
+
+Definition binop_propagate2_prop binop_pl :=
+  forall f1 f2 f3:float, is_nan _ _ f1 = true ->
+                         is_nan _ _ f2 = false -> is_nan _ _ f3 = false ->
+    binop_pl f1 f2 = (binop_pl f1 f3 : bool*nan_pl 53).
+
+(* Proved in Nan.v for different architectures *)
+Hypothesis binop_propagate2: binop_propagate2_prop binop_pl.
+
+Definition add: float -> float -> float := b64_plus binop_pl mode_NE. (**r addition *)
+Definition sub: float -> float -> float := b64_minus binop_pl mode_NE. (**r subtraction *)
+Definition mul: float -> float -> float := b64_mult binop_pl mode_NE. (**r multiplication *)
+Definition div: float -> float -> float := b64_div binop_pl mode_NE. (**r division *)
 
 Definition order_float (f1 f2:float): option Datatypes.comparison :=
   match f1, f2 with
-    | B754_nan,_ | _,B754_nan => None
-    | B754_infinity true, B754_infinity true 
+    | B754_nan _ _,_ | _,B754_nan _ _ => None
+    | B754_infinity true, B754_infinity true
     | B754_infinity false, B754_infinity false => Some Eq
     | B754_infinity true, _ => Some Lt
     | B754_infinity false, _ => Some Gt
@@ -229,7 +317,7 @@ Definition order_float (f1 f2:float): option Datatypes.comparison :=
   end.
 
 Definition cmp (c:comparison) (f1 f2:float) : bool := (**r comparison *)
-  match c with                                                                        
+  match c with
   | Ceq =>
       match order_float f1 f2 with Some Eq => true | _ => false end
   | Cne =>
@@ -242,7 +330,7 @@ Definition cmp (c:comparison) (f1 f2:float) : bool := (**r comparison *)
       match order_float f1 f2 with Some Gt => true | _ => false end
   | Cge =>
       match order_float f1 f2 with Some(Gt|Eq) => true | _ => false end
-  end. 
+  end.
 
 (** Conversions between floats and their concrete in-memory representation
     as a sequence of 64 bits (double precision) or 32 bits (single precision). *)
@@ -264,29 +352,17 @@ Ltac compute_this val :=
   let x := fresh in set val as x in *; vm_compute in x; subst x.
 
 Ltac smart_omega :=
-  simpl radix_val in *; simpl Zpower in *;  
+  simpl radix_val in *; simpl Zpower in *;
   compute_this Int.modulus; compute_this Int.half_modulus;
   compute_this Int.max_unsigned;
   compute_this Int.min_signed; compute_this Int.max_signed;
   compute_this Int64.modulus; compute_this Int64.half_modulus;
   compute_this Int64.max_unsigned;
   compute_this (Zpower_pos 2 1024); compute_this (Zpower_pos 2 53); compute_this (Zpower_pos 2 52);
-  omega.
-
-Theorem addf_commut: forall f1 f2, add f1 f2 = add f2 f1.
-Proof.
-  intros.
-  destruct f1, f2; simpl; try reflexivity; try (destruct b, b0; reflexivity).
-  rewrite Zplus_comm; rewrite Zmin_comm; reflexivity.
-Qed.
-
-Theorem subf_addf_opp: forall f1 f2, sub f1 f2 = add f1 (neg f2).
-Proof.
-  destruct f1, f2; reflexivity.
-Qed.
+  zify; omega.
 
 Lemma floatofbinary32_exact :
-  forall f, is_finite_strict _ _ f = true -> 
+  forall f, is_finite_strict _ _ f = true ->
     is_finite_strict _ _ (floatofbinary32 f) = true /\ B2R _ _ f = B2R _ _ (floatofbinary32 f).
 Proof.
   destruct f as [ | | |s m e]; try discriminate; intro.
@@ -305,10 +381,12 @@ Proof.
   now apply bpow_lt.
 Qed.
 
-Lemma binary32offloatofbinary32 :
-  forall f, binary32offloat (floatofbinary32 f) = f.
+Lemma binary32offloatofbinary32_num :
+  forall f, is_nan _ _ f = false ->
+            binary32offloat (floatofbinary32 f) = f.
 Proof.
-  intro; pose proof (floatofbinary32_exact f); destruct f as [ | | |s m e]; try reflexivity.
+  intros f Hnan; pose proof (floatofbinary32_exact f); destruct f as [ | | |s m e]; try reflexivity.
+  discriminate.
   specialize (H eq_refl); destruct H.
   destruct (floatofbinary32 (B754_finite 24 128 s m e e0)) as [ | | |s1 m1 e1]; try discriminate.
   unfold binary32offloat.
@@ -324,22 +402,76 @@ Proof.
   now apply generic_format_B2R.
 Qed.
 
+Lemma floatofbinary32offloatofbinary32_pl:
+  forall s pl,
+    prod_rect (fun _ => _) floatofbinary32_pl (prod_rect (fun _ => _) binary32offloat_pl (floatofbinary32_pl s pl)) = floatofbinary32_pl s pl.
+Proof.
+  destruct pl. unfold binary32offloat_pl, floatofbinary32_pl.
+  unfold transform_quiet_pl, proj1_sig. simpl.
+  f_equal. apply nan_payload_fequal.
+  unfold Pos.shiftr_nat. simpl.
+  rewrite !lor_idempotent. reflexivity.
+Qed.
+
+Lemma floatofbinary32offloatofbinary32 :
+  forall f, floatofbinary32 (binary32offloat (floatofbinary32 f)) = floatofbinary32 f.
+Proof.
+  destruct f; try (rewrite binary32offloatofbinary32_num; tauto).
+  unfold floatofbinary32, binary32offloat.
+  rewrite <- floatofbinary32offloatofbinary32_pl at 2.
+  reflexivity.
+Qed.
+
+Lemma binary32offloatofbinary32offloat_pl:
+  forall s pl,
+    prod_rect (fun _ => _) binary32offloat_pl (prod_rect  (fun _ => _) floatofbinary32_pl (binary32offloat_pl s pl)) = binary32offloat_pl s pl.
+Proof.
+  destruct pl. unfold binary32offloat_pl, floatofbinary32_pl. unfold prod_rect.
+  f_equal. apply nan_payload_fequal.
+  rewrite transform_quiet_pl_idempotent.
+  unfold transform_quiet_pl, proj1_sig.
+  change 51 with (29+22).
+  clear - x. revert x. unfold Pos.shiftr_nat, Pos.shiftl_nat.
+  induction (29)%nat. intro. simpl. apply lor_idempotent.
+  intro.
+  rewrite !nat_iter_succ_r with (f:=Pos.div2).
+  destruct x; simpl; try apply IHn.
+  clear IHn. induction n. reflexivity.
+  rewrite !nat_iter_succ_r with (f:=Pos.div2). auto.
+Qed.
+
+Lemma binary32offloatofbinary32offloat :
+  forall f, binary32offloat (floatofbinary32 (binary32offloat f)) = binary32offloat f.
+Proof.
+  destruct f; try (rewrite binary32offloatofbinary32_num; simpl; tauto).
+  unfold floatofbinary32, binary32offloat.
+  rewrite <- binary32offloatofbinary32offloat_pl at 2.
+  reflexivity.
+  rewrite binary32offloatofbinary32_num; simpl. auto.
+  unfold binary_normalize32.
+  pose proof (binary_normalize32_correct (cond_Zopp b (Z.pos m)) e b).
+  destruct binary_normalize; auto. simpl in H.
+  destruct Rlt_bool in H. intuition.
+  unfold binary_overflow in H. destruct n.
+  destruct overflow_to_inf in H; discriminate.
+Qed.
+
 Theorem singleoffloat_idem:
   forall f, singleoffloat (singleoffloat f) = singleoffloat f.
 Proof.
-  intros; unfold singleoffloat; rewrite binary32offloatofbinary32; reflexivity.
+  intros; unfold singleoffloat; rewrite binary32offloatofbinary32offloat; reflexivity.
 Qed.
 
 Theorem singleoflong_idem:
   forall n, singleoffloat (singleoflong n) = singleoflong n.
 Proof.
-  intros; unfold singleoffloat, singleoflong; rewrite binary32offloatofbinary32; reflexivity.
+  intros; unfold singleoffloat, singleoflong. rewrite floatofbinary32offloatofbinary32; reflexivity.
 Qed.
 
 Theorem singleoflongu_idem:
   forall n, singleoffloat (singleoflongu n) = singleoflongu n.
 Proof.
-  intros; unfold singleoffloat, singleoflongu; rewrite binary32offloatofbinary32; reflexivity.
+  intros; unfold singleoffloat, singleoflongu. rewrite floatofbinary32offloatofbinary32; reflexivity.
 Qed.
 
 Definition is_single (f: float) : Prop := exists s, f = floatofbinary32 s.
@@ -353,8 +485,8 @@ Qed.
 Theorem singleoffloat_of_single:
   forall f, is_single f -> singleoffloat f = f.
 Proof.
-  intros. destruct H as [s EQ]. subst f. unfold singleoffloat. 
-  rewrite binary32offloatofbinary32; reflexivity.
+  intros. destruct H as [s EQ]. subst f. unfold singleoffloat.
+  apply floatofbinary32offloatofbinary32.
 Qed.
 
 Theorem is_single_dec: forall f, {is_single f} + {~is_single f}.
@@ -402,9 +534,9 @@ Proof.
   replace 0%R with (0*bpow radix2 e0)%R by ring; apply Rmult_lt_compat_r;
     [apply bpow_gt_0; reflexivity|now apply (Z2R_lt _ 0)].
   apply Rmult_gt_0_compat; [apply (Z2R_lt 0); reflexivity|now apply bpow_gt_0].
-  destruct b, b0; try (now apply_Rcompare; apply H5); inversion H3; 
+  destruct b, b0; try (now apply_Rcompare; apply H5); inversion H3;
     try (apply_Rcompare; apply H4; rewrite H, H1 in H7; assumption);
-    try (apply_Rcompare; do 2 rewrite Z2R_opp, Ropp_mult_distr_l_reverse; 
+    try (apply_Rcompare; do 2 rewrite Z2R_opp, Ropp_mult_distr_l_reverse;
       apply Ropp_lt_contravar; apply H4; rewrite H, H1 in H7; assumption);
     rewrite H7, Rcompare_mult_r, Rcompare_Z2R by (apply bpow_gt_0); reflexivity.
 Qed.
@@ -472,7 +604,16 @@ Proof.
   destruct f.
   simpl; try destruct b; vm_compute; split; congruence.
   simpl; try destruct b; vm_compute; split; congruence.
-  simpl; vm_compute; split; congruence.
+  destruct n as [p Hp].
+  simpl. rewrite Z.ltb_lt in Hp.
+  apply Zlt_succ_le with (m:=52) in Hp.
+  apply Zpower_le with (r:=radix2) in Hp.
+  edestruct Fcore_digits.digits2_Pnat_correct.
+  rewrite Zpower_nat_Z in H0.
+  eapply Z.lt_le_trans in Hp; eauto.
+  unfold join_bits; destruct b.
+  compute_this ((2 ^ 11 + 2047) * 2 ^ 52). smart_omega.
+  compute_this ((0 + 2047) * 2 ^ 52). smart_omega.
   unfold bits_of_binary_float, join_bits.
   destruct (andb_prop _ _ e0); apply Zle_bool_imp_le in H0; apply Zeq_bool_eq in H; unfold FLT_exp in H.
   match goal with [H:Zmax ?x ?y = e|-_] => pose proof (Zle_max_l x y); pose proof (Zle_max_r x y) end.
@@ -494,7 +635,17 @@ Proof.
   destruct (binary32offloat f).
   simpl; try destruct b; vm_compute; split; congruence.
   simpl; try destruct b; vm_compute; split; congruence.
-  simpl; vm_compute; split; congruence.
+  destruct n as [p Hp].
+  simpl. rewrite Z.ltb_lt in Hp.
+  apply Zlt_succ_le with (m:=23) in Hp.
+  apply Zpower_le with (r:=radix2) in Hp.
+  edestruct Fcore_digits.digits2_Pnat_correct.
+  rewrite Zpower_nat_Z in H0.
+  eapply Z.lt_le_trans in Hp; eauto.
+  compute_this (radix2^23).
+  unfold join_bits; destruct b.
+  compute_this ((2 ^ 8 + 255) * 2 ^ 23). smart_omega.
+  compute_this ((0 + 255) * 2 ^ 23). smart_omega.
   unfold bits_of_binary_float, join_bits.
   destruct (andb_prop _ _ e0); apply Zle_bool_imp_le in H0; apply Zeq_bool_eq in H.
   unfold FLT_exp in H.
@@ -512,13 +663,13 @@ Qed.
 Theorem bits_of_singleoffloat:
   forall f, bits_of_single (singleoffloat f) = bits_of_single f.
 Proof.
-  intro; unfold singleoffloat, bits_of_single; rewrite binary32offloatofbinary32; reflexivity.
+  intro; unfold singleoffloat, bits_of_single; rewrite binary32offloatofbinary32offloat; reflexivity.
 Qed.
 
 Theorem singleoffloat_of_bits:
   forall b, singleoffloat (single_of_bits b) = single_of_bits b.
 Proof.
-  intro; unfold singleoffloat, single_of_bits; rewrite binary32offloatofbinary32; reflexivity.
+  intro; unfold singleoffloat, single_of_bits; rewrite floatofbinary32offloatofbinary32; reflexivity.
 Qed.
 
 Theorem single_of_bits_is_single:
@@ -539,7 +690,7 @@ Lemma round_exact:
     round radix2 (FLT_exp (3 - 1024 - 53) 53)
       (round_mode mode_NE) (Z2R n) = Z2R n.
 Proof.
-  intros; rewrite round_generic; [reflexivity|now apply valid_rnd_round_mode|].  
+  intros; rewrite round_generic; [reflexivity|now apply valid_rnd_round_mode|].
   apply generic_format_FLT; exists (Float radix2 n 0).
   unfold F2R, Fnum, Fexp, bpow; rewrite Rmult_1_r; intuition.
   pose proof (Zabs_spec n); now smart_omega.
@@ -551,7 +702,7 @@ Lemma binary_normalize64_exact:
     is_finite _ _ (binary_normalize64 n 0 false) = true.
 Proof.
   intros; pose proof (binary_normalize64_correct n 0 false).
-  unfold F2R, Fnum, Fexp, bpow in H0; rewrite Rmult_1_r, round_exact, Rlt_bool_true in H0; try assumption.
+  unfold F2R, Fnum, Fexp, bpow in H0; rewrite Rmult_1_r, round_exact, Rlt_bool_true in H0; try now intuition.
   rewrite <- Z2R_abs; apply Z2R_lt; pose proof (Zabs_spec n); now smart_omega.
 Qed.
 
@@ -560,7 +711,7 @@ Theorem floatofintu_floatofint_1:
   Int.ltu x ox8000_0000 = true ->
   floatofintu x = floatofint x.
 Proof.
-  unfold floatofintu, floatofint, Int.signed, Int.ltu; intro. 
+  unfold floatofintu, floatofint, Int.signed, Int.ltu; intro.
   change (Int.unsigned ox8000_0000) with Int.half_modulus.
   destruct (zlt (Int.unsigned x) Int.half_modulus); now intuition.
 Qed.
@@ -571,19 +722,19 @@ Theorem floatofintu_floatofint_2:
   floatofintu x = add (floatofint (Int.sub x ox8000_0000))
                       (floatofintu ox8000_0000).
 Proof.
-  unfold floatofintu, floatofint, Int.signed, Int.ltu, Int.sub; intros. 
+  unfold floatofintu, floatofint, Int.signed, Int.ltu, Int.sub; intros.
   pose proof (Int.unsigned_range x).
   compute_this (Int.unsigned ox8000_0000).
   destruct (zlt (Int.unsigned x) 2147483648); try  discriminate.
   rewrite Int.unsigned_repr by smart_omega.
   destruct (zlt ((Int.unsigned x) - 2147483648) Int.half_modulus).
-  unfold add, b64_plus. 
-  match goal with [|- _ = Bplus _ _ _ _ _ ?x ?y] =>
-    pose proof (Bplus_correct 53 1024 eq_refl eq_refl mode_NE x y) end.
+  unfold add, b64_plus.
+  match goal with [|- _ = Bplus _ _ _ _ _ _ ?x ?y] =>
+    pose proof (Bplus_correct 53 1024 eq_refl eq_refl binop_pl mode_NE x y) end.
   do 2 rewrite (fun x H => proj1 (binary_normalize64_exact x H)) in H1 by smart_omega.
   do 2 rewrite (fun x H => proj2 (binary_normalize64_exact x H)) in H1 by smart_omega.
-  rewrite <- Z2R_plus, round_exact in H1 by smart_omega. 
-  rewrite Rlt_bool_true in H1; 
+  rewrite <- Z2R_plus, round_exact in H1 by smart_omega.
+  rewrite Rlt_bool_true in H1;
     replace (Int.unsigned x - 2147483648 + 2147483648) with (Int.unsigned x) in * by ring.
   apply B2R_inj.
   destruct (binary_normalize64_exact (Int.unsigned x)); [now smart_omega|].
@@ -616,7 +767,7 @@ Proof.
   now apply valid_rnd_round_mode.
   apply generic_format_FIX. exists (Float radix2 (cond_Zopp b (Zpos m)) 0). split; reflexivity.
   split; [reflexivity|].
-  rewrite round_generic, Z2R_mult, Z2R_Zpower_pos, <- bpow_powerRZ; 
+  rewrite round_generic, Z2R_mult, Z2R_Zpower_pos, <- bpow_powerRZ;
     [reflexivity|now apply valid_rnd_round_mode|apply generic_format_F2R; discriminate].
   rewrite (inbetween_float_ZR_sign _ _ _ ((Zpos m) / Zpower_pos radix2 p)
     (new_location (Zpower_pos radix2 p) (Zpos m mod Zpower_pos radix2 p) loc_Exact)).
@@ -685,7 +836,7 @@ Theorem intuoffloat_correct:
     end.
 Proof.
   intro; pose proof (Zoffloat_correct f); unfold intuoffloat; destruct (Zoffloat f).
-  pose proof (Zle_bool_spec 0 z); pose proof (Zle_bool_spec z Int.max_unsigned). 
+  pose proof (Zle_bool_spec 0 z); pose proof (Zle_bool_spec z Int.max_unsigned).
   compute_this Int.max_unsigned; destruct H.
   inversion H0. inversion H1.
   rewrite <- (Int.unsigned_repr z) in H2 by smart_omega; split; assumption.
@@ -713,14 +864,14 @@ Proof.
   simpl B2R; change 0%R with (Z2R 0); change (-1)%R with (Z2R (-1)); split; apply Z2R_lt; reflexivity.
   pose proof (Int.unsigned_range i).
   unfold round, scaled_mantissa, B2R, F2R, Fnum, Fexp in H0 |- *; simpl bpow in H0; do 2 rewrite Rmult_1_r in H0;
-    apply eq_Z2R in H0. 
+    apply eq_Z2R in H0.
   split; apply Rnot_le_lt; intro.
   rewrite Ztrunc_ceil in H0;
     [apply Zceil_le in H3; change (-1)%R with (Z2R (-1)) in H3; rewrite Zceil_Z2R in H3; omega|].
   eapply Rle_trans; [now apply H3|apply (Z2R_le (-1) 0); discriminate].
   rewrite Ztrunc_floor in H0; [apply Zfloor_le in H3; rewrite Zfloor_Z2R in H3; now smart_omega|].
   eapply Rle_trans; [|now apply H3]; apply (Z2R_le 0); discriminate.
-Qed.  
+Qed.
 
 Theorem intuoffloat_intoffloat_1:
   forall x n,
@@ -738,7 +889,7 @@ Proof.
   destruct H4; [rewrite H2 in H4; discriminate|].
   apply intuoffloat_interval in H0; exfalso; destruct H0, H4.
   eapply Rlt_le_trans in H0; [|now apply H4]; apply (lt_Z2R (-1)) in H0; discriminate.
-  apply Rcompare_Lt_inv in H1; eapply Rle_lt_trans in H1; [|now apply H4]. 
+  apply Rcompare_Lt_inv in H1; eapply Rle_lt_trans in H1; [|now apply H4].
   unfold floatofintu in H1; rewrite (fun x H => proj1 (binary_normalize64_exact x H)) in H1;
     [apply lt_Z2R in H1; discriminate|split; reflexivity].
 Qed.
@@ -769,16 +920,16 @@ Proof.
     [rewrite H in H2; simpl B2R in H2; apply (eq_Z2R 0) in H2; discriminate|reflexivity].
   reflexivity.
   rewrite H in H2; apply Rcompare_Gt_inv in H2; pose proof (intuoffloat_interval _ _ H1).
-  unfold sub, b64_minus. 
-  exploit (Bminus_correct 53 1024 eq_refl eq_refl mode_NE x (floatofintu ox8000_0000)); [assumption|reflexivity|]; intro.
+  unfold sub, b64_minus.
+  exploit (Bminus_correct 53 1024 eq_refl eq_refl binop_pl mode_NE x (floatofintu ox8000_0000)); [assumption|reflexivity|]; intro.
   rewrite H, round_generic in H6.
-  match goal with [H6:if Rlt_bool ?x ?y then _ else _|-_] => 
+  match goal with [H6:if Rlt_bool ?x ?y then _ else _|-_] =>
     pose proof (Rlt_bool_spec x y); destruct (Rlt_bool x y) end.
-  destruct H6.
+  destruct H6 as [? []].
   match goal with [|- _ ?y = _] => pose proof (intoffloat_correct y); destruct (intoffloat y) end.
-  destruct H9.
+  destruct H10.
   f_equal; rewrite <- (Int.repr_signed i); unfold Int.sub; f_equal; apply eq_Z2R. 
-  rewrite Z2R_minus, H10, H4.
+  rewrite Z2R_minus, H11, H4.
   unfold round, scaled_mantissa, F2R, Fexp, Fnum, round_mode; simpl bpow; repeat rewrite Rmult_1_r;
     rewrite <- Z2R_minus; f_equal.
   rewrite (Ztrunc_floor (B2R _ _ x)), <- Zfloor_minus, <- Ztrunc_floor;
@@ -786,11 +937,11 @@ Proof.
   left; eapply Rlt_trans; [|now apply H2]; apply (Z2R_lt 0); reflexivity.
   try (change (0 ?= 53) with Lt in H6,H8).  (* for Coq 8.4 *)
   try (change (53 ?= 1024) with Lt in H6,H8).  (* for Coq 8.4 *)
-  exfalso; simpl Zcompare in H6, H8; rewrite H6, H8 in H9.
-  destruct H9 as [|[]]; [discriminate|..].
-  eapply Rle_trans in H9; [|apply Rle_0_minus; left; assumption]; apply (le_Z2R 0) in H9; apply H9; reflexivity.
-  eapply Rle_lt_trans in H9; [|apply Rplus_lt_compat_r; now apply (proj2 H5)].
-  rewrite <- Z2R_opp, <- Z2R_plus in H9; apply lt_Z2R in H9; discriminate.
+  exfalso; simpl Zcompare in H6, H8; rewrite H6, H8 in H10.
+  destruct H10 as [|[]]; [discriminate|..].
+  eapply Rle_trans in H10; [|apply Rle_0_minus; left; assumption]; apply (le_Z2R 0) in H10; apply H10; reflexivity.
+  eapply Rle_lt_trans in H10; [|apply Rplus_lt_compat_r; now apply (proj2 H5)].
+  rewrite <- Z2R_opp, <- Z2R_plus in H10; apply lt_Z2R in H10; discriminate.
   exfalso; inversion H7; rewrite Rabs_right in H8.
   eapply Rle_lt_trans in H8. apply Rle_not_lt in H8; [assumption|apply (bpow_le _ 31); discriminate].
   change (bpow radix2 31) with (Z2R(Zsucc Int.max_unsigned - Int.unsigned ox8000_0000)); rewrite Z2R_minus.
@@ -816,22 +967,22 @@ Lemma split_bits_or:
 Proof.
   intros.
   transitivity (split_bits 52 11 (join_bits 52 11 false (Int.unsigned x) 1075)).
-  - f_equal. rewrite Int64.ofwords_add'. reflexivity. 
-  - apply split_join_bits. 
+  - f_equal. rewrite Int64.ofwords_add'. reflexivity.
+  - apply split_join_bits.
     compute; auto.
-    generalize (Int.unsigned_range x). 
-    compute_this Int.modulus; compute_this (2^52); omega. 
+    generalize (Int.unsigned_range x).
+    compute_this Int.modulus; compute_this (2^52); omega.
     compute_this (2^11); omega.
 Qed.
 
 Lemma from_words_value:
   forall x,
-    B2R _ _ (from_words ox4330_0000 x) = 
+    B2R _ _ (from_words ox4330_0000 x) =
     (bpow radix2 52 + Z2R (Int.unsigned x))%R /\
     is_finite _ _ (from_words ox4330_0000 x) = true.
 Proof.
-  intros; unfold from_words, double_of_bits, b64_of_bits, binary_float_of_bits. 
-  rewrite B2R_FF2B; unfold is_finite; rewrite match_FF2B;
+  intros; unfold from_words, double_of_bits, b64_of_bits, binary_float_of_bits.
+  rewrite B2R_FF2B. rewrite is_finite_FF2B.
   unfold binary_float_of_bits_aux; rewrite split_bits_or; simpl; pose proof (Int.unsigned_range x).
   destruct (Int.unsigned x + Zpower_pos 2 52) eqn:?.
   exfalso; now smart_omega.
@@ -848,37 +999,37 @@ Theorem floatofintu_from_words:
     sub (from_words ox4330_0000 x) (from_words ox4330_0000 Int.zero).
 Proof.
   intros; destruct (Int.eq_dec x Int.zero); [subst; vm_compute; reflexivity|].
-  assert (Int.unsigned x <> 0). 
+  assert (Int.unsigned x <> 0).
   intro; destruct n; rewrite <- (Int.repr_unsigned x), H; reflexivity.
   pose proof (Int.unsigned_range x).
   pose proof (binary_normalize64_exact (Int.unsigned x)). destruct H1; [smart_omega|].
   unfold floatofintu, sub, b64_minus.
-  match goal with [|- _ = Bminus _ _ _ _ _ ?x ?y] =>
-    pose proof (Bminus_correct 53 1024 eq_refl eq_refl mode_NE x y) end.
+  match goal with [|- _ = Bminus _ _ _ _ _ _ ?x ?y] =>
+    pose proof (Bminus_correct 53 1024 eq_refl eq_refl binop_pl mode_NE x y) end.
   apply (fun f x y => f x y) in H3; try apply (fun x => proj2 (from_words_value x)).
   do 2 rewrite (fun x => proj1 (from_words_value x)) in H3.
   rewrite Int.unsigned_zero in H3.
   replace (bpow radix2 52 + Z2R (Int.unsigned x) -
     (bpow radix2 52 + Z2R 0))%R with (Z2R (Int.unsigned x)) in H3 by (simpl; ring).
   rewrite round_exact in H3 by smart_omega.
-  match goal with [H3:if Rlt_bool ?x ?y then _ else _ |- _] => 
-    pose proof (Rlt_bool_spec x y); destruct (Rlt_bool x y) end; destruct H3. 
+  match goal with [H3:if Rlt_bool ?x ?y then _ else _ |- _] =>
+    pose proof (Rlt_bool_spec x y); destruct (Rlt_bool x y) end; destruct H3 as [? []].
   try (change (53 ?= 1024) with Lt in H3,H5).  (* for Coq 8.4 *)
   simpl Zcompare in *; apply B2R_inj;
     try match goal with [H':B2R _ _ ?f = _ , H'':is_finite _ _ ?f = true |- is_finite_strict _ _ ?f = true] => 
       destruct f; [
         simpl in H'; change 0%R with (Z2R 0) in H'; apply eq_Z2R in H'; now destruct (H (eq_sym H')) | 
         discriminate H'' | discriminate H'' | reflexivity
-      ] 
+      ]
     end.
   rewrite H3; assumption.
-  inversion H4; change (bpow radix2 1024) with (Z2R (radix2 ^ 1024)) in H6; rewrite <- Z2R_abs in H6.
-  apply le_Z2R in H6; pose proof (Zabs_spec (Int.unsigned x));
+  inversion H4; change (bpow radix2 1024) with (Z2R (radix2 ^ 1024)) in H5; rewrite <- Z2R_abs in H5.
+  apply le_Z2R in H5; pose proof (Zabs_spec (Int.unsigned x));
     exfalso; now smart_omega.
 Qed.
 
 Lemma ox8000_0000_signed_unsigned:
-  forall x, 
+  forall x,
     Int.unsigned (Int.add x ox8000_0000) = Int.signed x + Int.half_modulus.
 Proof.
   intro; unfold Int.signed, Int.add; pose proof (Int.unsigned_range x).
@@ -902,16 +1053,16 @@ Local Transparent Int.repr Int64.repr.
   pose proof (Int.signed_range x).
   pose proof (binary_normalize64_exact (Int.signed x)); destruct H1; [now smart_omega|].
   unfold floatofint, sub, b64_minus.
-  match goal with [|- _ = Bminus _ _ _ _ _ ?x ?y] =>
-    pose proof (Bminus_correct 53 1024 eq_refl eq_refl mode_NE x y) end.
+  match goal with [|- _ = Bminus _ _ _ _ _ _ ?x ?y] =>
+    pose proof (Bminus_correct 53 1024 eq_refl eq_refl binop_pl mode_NE x y) end.
   apply (fun f x y => f x y) in H3; try apply (fun x => proj2 (from_words_value x)).
   do 2 rewrite (fun x => proj1 (from_words_value x)) in H3.
   replace (bpow radix2 52 + Z2R (Int.unsigned (Int.add x ox8000_0000)) -
     (bpow radix2 52 + Z2R (Int.unsigned ox8000_0000)))%R with (Z2R (Int.signed x)) in H3
   by (rewrite ox8000_0000_signed_unsigned; rewrite Z2R_plus; simpl; ring).
   rewrite round_exact in H3 by smart_omega.
-  match goal with [H3:if Rlt_bool ?x ?y then _ else _ |- _] => 
-    pose proof (Rlt_bool_spec x y); destruct (Rlt_bool x y) end; destruct H3. 
+  match goal with [H3:if Rlt_bool ?x ?y then _ else _ |- _] =>
+    pose proof (Rlt_bool_spec x y); destruct (Rlt_bool x y) end; destruct H3 as [? []].
   try (change (0 ?= 53) with Lt in H3,H5).  (* for Coq 8.4 *)
   try (change (53 ?= 1024) with Lt in H3,H5).  (* for Coq 8.4 *)
   simpl Zcompare in *; apply B2R_inj;
@@ -919,11 +1070,11 @@ Local Transparent Int.repr Int64.repr.
       destruct f; [
         simpl in H'; change 0%R with (Z2R 0) in H'; apply eq_Z2R in H'; now destruct (H (eq_sym H')) | 
         discriminate H'' | discriminate H'' | reflexivity
-      ] 
+      ]
     end.
   rewrite H3; assumption.
-  inversion H4; unfold bpow in H6; rewrite <- Z2R_abs in H6; 
-    apply le_Z2R in H6; pose proof (Zabs_spec (Int.signed x)); exfalso; now smart_omega.
+  inversion H4; unfold bpow in H5; rewrite <- Z2R_abs in H5;
+    apply le_Z2R in H5; pose proof (Zabs_spec (Int.signed x)); exfalso; now smart_omega.
 Qed.
 
 (** Conversions from 32-bit integers to single-precision floats can
@@ -936,10 +1087,10 @@ Lemma is_finite_strict_ge_1:
   (1 <= Rabs (B2R _ _ f))%R ->
   is_finite_strict _ _ f = true.
 Proof.
-  intros. destruct f; auto. simpl in H0. 
+  intros. destruct f; auto. simpl in H0.
   change 0%R with (Z2R 0) in H0.
   change 1%R with (Z2R 1) in H0.
-  rewrite <- Z2R_abs in H0. 
+  rewrite <- Z2R_abs in H0.
   exploit le_Z2R; eauto.
 Qed.
 
@@ -954,29 +1105,27 @@ Proof.
   subst n; reflexivity.
   exploit binary_normalize64_exact; eauto. intros [A B].
   destruct (binary_normalize64 n 0 false) as [ | | | s m e] eqn:B64; simpl in *.
-- assert (0 = n) by (apply eq_Z2R; auto). subst n. simpl in GT. omegaContradiction. 
+- assert (0 = n) by (apply eq_Z2R; auto). subst n. simpl in GT. omegaContradiction.
 - discriminate.
 - discriminate.
 - set (n1 := cond_Zopp s (Z.pos m)) in *.
-  generalize (binary_normalize32_correct n1 e s). 
+  generalize (binary_normalize32_correct n1 e s).
   fold (binary_normalize32 n1 e s). intros C.
-  generalize (binary_normalize32_correct n 0 false). 
+  generalize (binary_normalize32_correct n 0 false).
   fold (binary_normalize32 n 0 false). intros D.
   assert (A': @F2R radix2 {| Fnum := n; Fexp := 0 |} = Z2R n).
-  {
-    unfold F2R. apply Rmult_1_r.
-  }
-  rewrite A in C. rewrite A' in D. 
+  { unfold F2R. apply Rmult_1_r. }
+  rewrite A in C. rewrite A' in D.
   destruct (Rlt_bool
          (Rabs
             (round radix2 (FLT_exp (3 - 128 - 24) 24) (round_mode mode_NE)
                (Z2R n))) (bpow radix2 128)).
-+ destruct C as [C1 C2]; destruct D as [D1 D2].
++ destruct C as [C1 [C2 _]]; destruct D as [D1 [D2 _]].
   assert (1 <= Rabs (round radix2 (FLT_exp (3 - 128 - 24) 24) (round_mode mode_NE) (Z2R n)))%R.
-  { apply abs_round_ge_generic. 
+  { apply abs_round_ge_generic.
     apply fexp_correct. red. omega.
     apply valid_rnd_round_mode.
-    apply generic_format_bpow with (e := 0). compute. congruence. 
+    apply generic_format_bpow with (e := 0). compute. congruence.
     rewrite <- Z2R_abs. change 1%R with (Z2R 1). apply Z2R_le. omega. }
   apply B2R_inj.
   apply is_finite_strict_ge_1; auto. rewrite C1; auto.
@@ -988,17 +1137,102 @@ Qed.
 Theorem singleofint_floatofint:
   forall n, singleofint n = singleoffloat (floatofint n).
 Proof.
-  intros. symmetry. apply single_float_of_int. 
+  intros. symmetry. apply single_float_of_int.
   generalize (Int.signed_range n). smart_omega.
 Qed.
 
 Theorem singleofintu_floatofintu:
   forall n, singleofintu n = singleoffloat (floatofintu n).
 Proof.
-  intros. symmetry. apply single_float_of_int. 
+  intros. symmetry. apply single_float_of_int.
   generalize (Int.unsigned_range n). smart_omega.
 Qed.
 
+Theorem mul2_add:
+  forall f, add f f = mul f (floatofint (Int.repr 2%Z)).
+Proof.
+  intros. unfold add, b64_plus, mul, b64_mult.
+  destruct (is_finite_strict _ _ f) eqn:EQFINST.
+  - assert (EQFIN:is_finite _ _ f = true) by (destruct f; simpl in *; congruence).
+    pose proof (Bplus_correct 53 1024 eq_refl eq_refl binop_pl mode_NE f f EQFIN EQFIN).
+    pose proof (Bmult_correct 53 1024 eq_refl eq_refl binop_pl mode_NE f
+                              (floatofint (Int.repr 2%Z))).
+    rewrite <- double, Rmult_comm in H.
+    replace (B2R 53 1024 (floatofint (Int.repr 2))) with 2%R in H0 by (compute; field).
+    destruct Rlt_bool.
+    + destruct H0 as [? []], H as [? []].
+      rewrite EQFIN in H1.
+      apply B2R_Bsign_inj; auto.
+      etransitivity. apply H. symmetry. apply H0.
+      etransitivity. apply H4. symmetry. etransitivity. apply H2.
+      destruct Bmult; try reflexivity; discriminate.
+      simpl. rewrite xorb_false_r.
+      erewrite <- Rmult_0_l, Rcompare_mult_r.
+      destruct f; try discriminate EQFINST.
+      simpl. unfold F2R.
+      erewrite <- Rmult_0_l, Rcompare_mult_r.
+      rewrite Rcompare_Z2R with (y:=0).
+      destruct b; reflexivity.
+      apply bpow_gt_0.
+      apply (Z2R_lt 0 2). omega.
+    + destruct H.
+      apply B2FF_inj.
+      etransitivity. apply H.
+      symmetry. etransitivity. apply H0.
+      f_equal. destruct Bsign; reflexivity.
+  - destruct f as [[]|[]| |]; try discriminate; try reflexivity.
+    simpl. erewrite binop_propagate1; reflexivity.
+Qed.
+
+Program Definition pow2_float (b:bool) (e:Z) (H:-1023 < e < 1023) : float :=
+  B754_finite _ _ b (nat_iter 52 xO xH) (e-52) _.
+Next Obligation.
+  unfold Fappli_IEEE.bounded, canonic_mantissa.
+  rewrite andb_true_iff, Zle_bool_true by omega. split; auto.
+  apply Zeq_bool_true. unfold FLT_exp. simpl Z.of_nat.
+  apply Z.max_case_strong; omega.
+Qed.
+
+Theorem mul_div_pow2:
+  forall b e f H H',
+    mul f (pow2_float b e H) = div f (pow2_float b (-e) H').
+Proof.
+  intros. unfold mul, b64_mult, div, b64_div.
+  pose proof (Bmult_correct 53 1024 eq_refl eq_refl binop_pl mode_NE f (pow2_float b e H)).
+  pose proof (Bdiv_correct 53 1024 eq_refl eq_refl binop_pl mode_NE f (pow2_float b (-e) H')).
+  lapply H1. clear H1. intro.
+  change (is_finite 53 1024 (pow2_float b e H)) with true in H0.
+  unfold Rdiv in H1.
+  replace (/ B2R 53 1024 (pow2_float b (-e) H'))%R
+    with (B2R 53 1024 (pow2_float b e H)) in H1.
+  destruct (is_finite _ _ f) eqn:EQFIN.
+  - destruct Rlt_bool.
+    + destruct H0 as [? []], H1 as [? []].
+      apply B2R_Bsign_inj; auto.
+      etransitivity. apply H0. symmetry. apply H1.
+      etransitivity. apply H3. destruct Bmult; try discriminate H2; reflexivity.
+      symmetry. etransitivity. apply H5. destruct Bdiv; try discriminate H4; reflexivity.
+      reflexivity.
+    + apply B2FF_inj.
+      etransitivity. apply H0. symmetry. etransitivity. apply H1.
+      reflexivity.
+  - destruct f; try discriminate EQFIN. reflexivity.
+    simpl. erewrite binop_propagate2. reflexivity.
+    reflexivity. reflexivity. reflexivity.
+  - simpl.
+    assert ((4503599627370496 * bpow radix2 (e - 52))%R =
+            (/ (4503599627370496 * bpow radix2 (- e - 52)))%R).
+    { etransitivity. symmetry. apply (bpow_plus radix2 52).
+      symmetry. etransitivity. apply f_equal. symmetry. apply (bpow_plus radix2 52).
+      rewrite <- bpow_opp. f_equal. ring. }
+    destruct b. unfold cond_Zopp.
+    rewrite !F2R_Zopp, <- Ropp_inv_permute. f_equal. auto.
+    intro. apply F2R_eq_0_reg in H3. omega.
+    apply H2.
+  - simpl. intro. apply F2R_eq_0_reg in H2.
+    destruct b; simpl in H2; omega.
+Qed.
+
 Global Opaque
   zero eq_dec neg abs singleoffloat intoffloat intuoffloat floatofint floatofintu
   add sub mul div cmp bits_of_double double_of_bits bits_of_single single_of_bits from_words.