Tighten and prove correct the underflow/overflow bounds for parsing of FP literals.

This is a follow-up to commit 350354cb.
- Move Float.build_from_parsed to Fappli_IEEE_extra.Bparse
- Add early checks for overflow and underflow and prove them correct.
- Improve speed of Bparse by using a fast exponentiation (square-and-multiply).

cfrontend/C2C.ml

@@ -584,49 +584,37 @@ let checkFloatOverflow f =
       warning "Floating-point literal converts to Not-a-Number"
 
 let convertFloat f kind =
-  let iexp =(int_of_string f.C.exp)
-  and sign = match f.C.exp.[0] with '-' -> true | _ -> false in
-  (* Exp cannot be larger than 1023 or smaller than -1022 *)
-  if iexp > 1024  || iexp < -1034 then
-    let f =  if iexp > 1024 then
-      Fappli_IEEE.B754_infinity sign
-    else
-      Fappli_IEEE.B754_zero sign in
-    match kind with
-    | FFloat ->
-        checkFloatOverflow f;
-        Ctyping.econst_single f
-    | FDouble | FLongDouble ->
-        checkFloatOverflow f;
-        Ctyping.econst_float f
-  else
-    let mant = z_of_str f.C.hex (f.C.intPart ^ f.C.fracPart) 0 in
-    match mant with
+  let mant = z_of_str f.C.hex (f.C.intPart ^ f.C.fracPart) 0 in
+  match mant with
     | Z.Z0 ->
-        begin match kind with
-        | FFloat ->
-	    Ctyping.econst_single (Float.to_single Float.zero)
-        | FDouble | FLongDouble ->
-	    Ctyping.econst_float Float.zero
-        end
+      begin match kind with
+      | FFloat ->
+	  Ctyping.econst_single (Float.to_single Float.zero)
+      | FDouble | FLongDouble ->
+	  Ctyping.econst_float Float.zero
+      end
     | Z.Zpos mant ->
-        let sgExp = match f.C.exp.[0] with '+' | '-' -> true | _ -> false in
-        let exp = z_of_str false f.C.exp (if sgExp then 1 else 0) in
-        let exp = if f.C.exp.[0] = '-' then Z.neg exp else exp in
-        let shift_exp =
-	  (if f.C.hex then 4 else 1) * String.length f.C.fracPart in
-        let exp = Z.sub exp (Z.of_uint shift_exp) in
-        let base = P.of_int (if f.C.hex then 2 else 10) in
-        begin match kind with
-        | FFloat ->
-	    let f =  Float32.from_parsed base mant exp in
-            checkFloatOverflow f;
-            Ctyping.econst_single f
-        | FDouble | FLongDouble ->
-	    let f =  Float.from_parsed base mant exp in
-            checkFloatOverflow f;
-            Ctyping.econst_float f
-        end
+
+      let sgExp = match f.C.exp.[0] with '+' | '-' -> true | _ -> false in
+      let exp = z_of_str false f.C.exp (if sgExp then 1 else 0) in
+      let exp = if f.C.exp.[0] = '-' then Z.neg exp else exp in
+      let shift_exp =
+	(if f.C.hex then 4 else 1) * String.length f.C.fracPart in
+      let exp = Z.sub exp (Z.of_uint shift_exp) in
+
+      let base = P.of_int (if f.C.hex then 2 else 10) in
+
+      begin match kind with
+      | FFloat ->
+	  let f = Float32.from_parsed base mant exp in
+          checkFloatOverflow f;
+          Ctyping.econst_single f
+      | FDouble | FLongDouble ->
+	  let f = Float.from_parsed base mant exp in
+          checkFloatOverflow f;
+          Ctyping.econst_float f
+      end
+
     | Z.Zneg _ -> assert false
 
 (** Expressions *)

lib/Fappli_IEEE_extra.v

@@ -1110,6 +1110,224 @@ Proof.
   + erewrite NAN; eauto. 
 Qed.
 
+(** ** Conversion from scientific notation *)
+
+(** Russian peasant exponentiation *)
+
+Fixpoint pos_pow (x y: positive) : positive :=
+  match y with
+  | xH => x
+  | xO y => Pos.square (pos_pow x y)
+  | xI y => Pos.mul x (Pos.square (pos_pow x y))
+  end.
+
+Lemma pos_pow_spec:
+  forall x y, Z.pos (pos_pow x y) = Z.pos x ^ Z.pos y.
+Proof.
+  intros x.
+  assert (REC: forall y a, Pos.iter y (Pos.mul x) a = Pos.mul (pos_pow x y) a).
+  { induction y; simpl; intros. 
+  - rewrite ! IHy, Pos.square_spec, ! Pos.mul_assoc. auto.
+  - rewrite ! IHy, Pos.square_spec, ! Pos.mul_assoc. auto.
+  - auto.
+  }
+  intros. simpl. rewrite <- Pos2Z.inj_pow_pos. unfold Pos.pow. rewrite REC. rewrite Pos.mul_1_r. auto.
+Qed. 
+
+(** Given a base [base], a mantissa [m] and an exponent [e], the following function
+  computes the FP number closest to [m * base ^ e], using round to odd, ties break to even.
+  The algorithm is naive, computing [base ^ |e|] exactly before doing a multiplication or
+  division with [m].  However, we treat specially very large or very small values of [e],
+  when the result is known to be [+infinity] or [0.0] respectively. *)
+
+Definition Bparse (base: positive) (m: positive) (e: Z): binary_float :=
+  match e with
+  | Z0 => 
+     BofZ (Zpos m)
+  | Zpos p =>
+     if e * Z.log2 (Zpos base) <? emax
+     then BofZ (Zpos m * Zpos (pos_pow base p))
+     else B754_infinity _ _ false
+  | Zneg p =>
+     if e * Z.log2 (Zpos base) + Z.log2_up (Zpos m) <? emin
+     then B754_zero _ _ false
+     else FF2B prec emax _ (proj1 (Bdiv_correct_aux prec emax prec_gt_0_ Hmax mode_NE
+                                     false m Z0 false (pos_pow base p) Z0))
+  end.
+
+(** Properties of [Z.log2] and [Z.log2_up]. *)
+
+Lemma Zpower_log:
+  forall (base: radix) n,
+  0 < n ->
+  2 ^ (n * Z.log2 base) <= base ^ n <= 2 ^ (n * Z.log2_up base).
+Proof.
+  intros.
+  assert (A: 0 < base) by apply radix_gt_0.
+  assert (B: 0 <= Z.log2 base) by apply Z.log2_nonneg.
+  assert (C: 0 <= Z.log2_up base) by apply Z.log2_up_nonneg.
+  destruct (Z.log2_spec base) as [D E]; auto.
+  destruct (Z.log2_up_spec base) as [F G]. apply radix_gt_1. 
+  assert (K: 0 <= 2 ^ Z.log2 base) by (apply Z.pow_nonneg; omega).
+  rewrite ! (Zmult_comm n). rewrite ! Z.pow_mul_r by omega.
+  split; apply Z.pow_le_mono_l; omega.
+Qed.
+
+Lemma bpow_log_pos:
+  forall (base: radix) n,
+  0 < n ->
+  (bpow radix2 (n * Z.log2 base)%Z <= bpow base n)%R.
+Proof.
+  intros. rewrite <- ! Z2R_Zpower. apply Z2R_le; apply Zpower_log; auto.
+  omega.
+  rewrite Zmult_comm; apply Zmult_gt_0_le_0_compat. omega. apply Z.log2_nonneg.
+Qed.
+
+Lemma bpow_log_neg:
+  forall (base: radix) n,
+  n < 0 ->
+  (bpow base n <= bpow radix2 (n * Z.log2 base)%Z)%R.
+Proof.
+  intros. set (m := -n). replace n with (-m) by (unfold m; omega).
+  rewrite ! Z.mul_opp_l, ! bpow_opp. apply Rinv_le. 
+  apply bpow_gt_0.
+  apply bpow_log_pos. unfold m; omega.
+Qed.
+
+(** Overflow and underflow conditions. *)
+
+Lemma round_integer_overflow:
+  forall (base: radix) e m,
+  0 < e ->
+  emax <= e * Z.log2 base ->
+  (bpow radix2 emax <= round radix2 fexp (round_mode mode_NE) (Z2R (Zpos m) * bpow base e))%R.
+Proof.
+  intros. 
+  rewrite <- (round_generic radix2 fexp (round_mode mode_NE) (bpow radix2 emax)); auto.
+  apply round_le; auto. apply fexp_correct; auto. apply valid_rnd_round_mode. 
+  rewrite <- (Rmult_1_l (bpow radix2 emax)). apply Rmult_le_compat.
+  apply Rle_0_1.
+  apply bpow_ge_0.
+  apply (Z2R_le 1). zify; omega.
+  eapply Rle_trans. eapply bpow_le. eassumption. apply bpow_log_pos; auto.
+  apply generic_format_FLT. exists (Float radix2 1 emax).
+  split. unfold F2R; simpl. ring.
+  split. simpl. apply (Zpower_gt_1 radix2); auto.
+  simpl. unfold emin; red in prec_gt_0_; omega.
+Qed.
+
+Lemma round_NE_underflows:
+  forall x,
+  (0 <= x <= bpow radix2 (emin - 1))%R ->
+  round radix2 fexp (round_mode mode_NE) x = 0%R.
+Proof.
+  intros.
+  set (eps := bpow radix2 (emin - 1)) in *.
+  assert (A: round radix2 fexp (round_mode mode_NE) eps = 0%R).
+  { unfold round. simpl. 
+    assert (E: canonic_exp radix2 fexp eps = emin).
+    { unfold canonic_exp, eps. rewrite ln_beta_bpow. unfold fexp, FLT_exp. zify; red in prec_gt_0_; omega. }
+    unfold scaled_mantissa; rewrite E. 
+    assert (P: (eps * bpow radix2 (-emin) = / 2)%R).
+    { unfold eps. rewrite <- bpow_plus. replace (emin - 1 + -emin) with (-1) by omega. auto. }
+    rewrite P. unfold Znearest.  
+    assert (F: Zfloor (/ 2)%R = 0).
+    { apply Zfloor_imp. 
+      split. apply Rlt_le. apply Rinv_0_lt_compat. apply (Z2R_lt 0 2). omega.
+      change (Z2R (0 + 1)) with 1%R. rewrite <- Rinv_1 at 3. apply Rinv_1_lt_contravar. apply Rle_refl. apply (Z2R_lt 1 2). omega.
+    }
+    rewrite F. change (Z2R 0) with 0%R. rewrite Rminus_0_r. rewrite Rcompare_Eq by auto. 
+    simpl. unfold F2R; simpl. apply Rmult_0_l. 
+  }
+  apply Rle_antisym.
+- rewrite <- A. apply round_le. apply fexp_correct; auto. apply valid_rnd_round_mode. tauto.
+- rewrite <- (round_0 radix2 fexp (round_mode mode_NE)). 
+  apply round_le. apply fexp_correct; auto. apply valid_rnd_round_mode. tauto.
+Qed.
+
+Lemma round_integer_underflow:
+  forall (base: radix) e m,
+  e < 0 ->
+  e * Z.log2 base + Z.log2_up (Zpos m) < emin ->
+  round radix2 fexp (round_mode mode_NE) (Z2R (Zpos m) * bpow base e) = 0%R.
+Proof.
+  intros. apply round_NE_underflows. split.
+- apply Rmult_le_pos. apply (Z2R_le 0). zify; omega. apply bpow_ge_0.
+- apply Rle_trans with (bpow radix2 (Z.log2_up (Z.pos m) + e * Z.log2 base)).
++ rewrite bpow_plus. apply Rmult_le_compat. 
+  apply (Z2R_le 0); zify; omega.
+  apply bpow_ge_0.
+  rewrite <- Z2R_Zpower. apply Z2R_le.
+  destruct (Z.eq_dec (Z.pos m) 1). 
+  rewrite e0. simpl. omega.
+  apply Z.log2_up_spec. zify; omega.
+  apply Z.log2_up_nonneg. 
+  apply bpow_log_neg. auto.
++ apply bpow_le. omega. 
+Qed.
+
+(** Correctness of Bparse *)
+
+Theorem Bparse_correct:
+  forall b m e (BASE: 2 <= Zpos b),
+  let base := {| radix_val := Zpos b; radix_prop := Zle_imp_le_bool _ _ BASE |} in
+  let r := round radix2 fexp (round_mode mode_NE) (Z2R (Zpos m) * bpow base e) in
+  if Rlt_bool (Rabs r) (bpow radix2 emax) then
+     B2R _ _ (Bparse b m e) = r
+  /\ is_finite _ _ (Bparse b m e) = true
+  /\ Bsign _ _ (Bparse b m e) = false
+  else
+    B2FF _ _ (Bparse b m e) = F754_infinity false.
+Proof.
+  intros. 
+  assert (A: forall x, @F2R radix2 {| Fnum := x; Fexp := 0 |} = Z2R x).
+  { intros. unfold F2R, Fnum; simpl. ring. }
+  unfold Bparse, r. destruct e as [ | e | e].
+- (* e = Z0 *)
+  change (bpow base 0) with 1%R. rewrite Rmult_1_r. 
+  exact (BofZ_correct (Z.pos m)).
+- (* e = Zpos e *)
+  destruct (Z.ltb_spec (Z.pos e * Z.log2 (Z.pos b)) emax).
++ (* no overflow *)
+  rewrite pos_pow_spec. rewrite <- Z2R_Zpower by (zify; omega). rewrite <- Z2R_mult.
+  replace false with (Z.pos m * Z.pos b ^ Z.pos e <? 0). 
+  exact (BofZ_correct (Z.pos m * Z.pos b ^ Z.pos e)).
+  rewrite Z.ltb_ge. rewrite Zmult_comm. apply Zmult_gt_0_le_0_compat. zify; omega.  apply (Zpower_ge_0 base). 
++ (* overflow *)
+  rewrite Rlt_bool_false. auto. eapply Rle_trans; [idtac|apply Rle_abs]. 
+  apply (round_integer_overflow base). zify; omega. auto. 
+- (* e = Zneg e *)
+  destruct (Z.ltb_spec (Z.neg e * Z.log2 (Z.pos b) + Z.log2_up (Z.pos m)) emin).
++ (* undeflow *)
+  rewrite round_integer_underflow; auto.
+  rewrite Rlt_bool_true. auto.  
+  replace (Rabs 0)%R with 0%R. apply bpow_gt_0. apply (Z2R_abs 0). 
+  zify; omega.
++ (* no underflow *)
+  generalize (Bdiv_correct_aux prec emax prec_gt_0_ Hmax mode_NE false m 0 false (pos_pow b e) 0).
+  set (f := match Fdiv_core_binary prec (Z.pos m) 0 (Z.pos (pos_pow b e)) 0 with
+      | (0, _, _) => F754_nan false 1
+      | (Z.pos mz0, ez, lz) =>
+          binary_round_aux prec emax mode_NE (xorb false false) mz0 ez lz
+      | (Z.neg _, _, _) => F754_nan false 1
+      end).
+  fold emin; fold fexp. rewrite ! A. unfold cond_Zopp. rewrite pos_pow_spec. 
+  assert (B: (Z2R (Z.pos m) / Z2R (Z.pos b ^ Z.pos e) =
+              Z2R (Z.pos m) * bpow base (Z.neg e))%R).
+  { change (Z.neg e) with (- (Z.pos e)). rewrite bpow_opp. auto. }
+  rewrite B. intros [P Q]. 
+  destruct (Rlt_bool
+     (Rabs
+        (round radix2 fexp (round_mode mode_NE)
+           (Z2R (Z.pos m) * bpow base (Z.neg e))))
+    (bpow radix2 emax)).
+* destruct Q as (Q1 & Q2 & Q3). 
+  split. rewrite B2R_FF2B, Q1. auto.
+  split. rewrite is_finite_FF2B. auto. 
+  rewrite Bsign_FF2B. auto.
+* rewrite B2FF_FF2B. auto.
+Qed.
+
 End Extra_ops.
 
 (** ** Conversions between two FP formats *)

lib/Floats.v

@@ -92,100 +92,6 @@ Proof.
   destruct x as [[]|]; simpl; intros; discriminate.
 Qed.
 
-Section FP_PARSING.
-
-Variables prec emax: Z.
-Context (prec_gt_0 : Prec_gt_0 prec).
-Hypothesis Hmax : (prec < emax)%Z.
-
-(** Function used to convert literals into FP numbers during parsing.
-   [intPart] is the mantissa
-   [expPart] is the exponent
-   [base] is the base for the exponent (usually 10 or 16).
-   The result is [intPart * base^expPart] rounded to the nearest FP number,
-   ties to even. *)
-
-Definition build_from_parsed (base:positive) (intPart:positive) (expPart:Z) : binary_float prec emax :=
-  match expPart with
-    | Z0 =>
-      binary_normalize prec emax prec_gt_0 Hmax mode_NE (Zpos intPart) Z0 false
-    | Zpos p =>
-      binary_normalize prec emax prec_gt_0 Hmax mode_NE ((Zpos intPart) * Zpower_pos (Zpos base) p) Z0 false
-    | Zneg p =>
-      FF2B prec emax _ (proj1 (Bdiv_correct_aux prec emax prec_gt_0 Hmax mode_NE
-                                  false intPart Z0 false (Pos.pow base p) Z0))
-  end.
-
-Let emin := (3 - emax - prec)%Z.
-Let fexp := FLT_exp emin prec.
-
-Theorem build_from_parsed_correct:
-  forall base m e (BASE: 2 <= Zpos base),
-  let base_r := {| radix_val := Zpos base; radix_prop := Zle_imp_le_bool _ _ BASE |} in
-  let r := round radix2 fexp (round_mode mode_NE) (Z2R (Zpos m) * bpow base_r e) in
-  if Rlt_bool (Rabs r) (bpow radix2 emax) then
-     B2R _ _ (build_from_parsed base m e) = r
-  /\ is_finite _ _ (build_from_parsed base m e) = true
-  /\ Bsign _ _ (build_from_parsed base m e) = false
-  else
-    B2FF _ _ (build_from_parsed base m e) = F754_infinity false.
-Proof.
-  intros. 
-  assert (A: forall x, @F2R radix2 {| Fnum := x; Fexp := 0 |} = Z2R x).
-  { intros. unfold F2R, Fnum; simpl. ring. }
-  unfold build_from_parsed, r; destruct e. 
-- change (bpow base_r 0) with (1%R). rewrite Rmult_1_r.  
-  generalize (binary_normalize_correct _ _ _ Hmax mode_NE (Zpos m) 0 false).
-  fold emin; fold fexp. rewrite ! A.
-  destruct (Rlt_bool (Rabs (round radix2 fexp (round_mode mode_NE) (Z2R (Z.pos m))))
-                     (bpow radix2 emax)).
-  + intros (P & Q & R). split. auto. split. auto. rewrite R. rewrite Rcompare_Gt; auto. 
-    apply (Z2R_lt 0). compute; auto.
-  + intros P; rewrite P. unfold binary_overflow.
-    replace (Rlt_bool (Z2R (Z.pos m)) 0%R) with false. auto. 
-    rewrite Rlt_bool_false; auto. apply (Z2R_le 0). xomega.
-- generalize (binary_normalize_correct _ _ _ Hmax mode_NE (Zpos m * Z.pow_pos (Zpos base) p) 0 false).
-  fold emin; fold fexp. rewrite ! A.
-  assert (B: Z2R (Z.pos m * Z.pow_pos (Z.pos base) p) = (Z2R (Z.pos m) * bpow base_r (Z.pos p))%R).
-  { unfold bpow. apply Z2R_mult. }
-  rewrite B. 
-  destruct (Rlt_bool
-     (Rabs
-        (round radix2 fexp (round_mode mode_NE)
-           (Z2R (Z.pos m) * bpow base_r (Z.pos p)))) (bpow radix2 emax)).
-  + intros (P & Q & R). split. auto. split. auto. rewrite R. rewrite Rcompare_Gt; auto. 
-    apply Rmult_lt_0_compat. apply (Z2R_lt 0). xomega. apply bpow_gt_0. 
-  + intros P. rewrite P. unfold binary_overflow. 
-    replace (Rlt_bool (Z2R (Z.pos m) * bpow base_r (Z.pos p)) 0) with false.
-    auto.
-    rewrite Rlt_bool_false; auto. apply Rlt_le. apply Rmult_lt_0_compat. apply (Z2R_lt 0). xomega. apply bpow_gt_0. 
-- generalize (Bdiv_correct_aux prec emax prec_gt_0 Hmax mode_NE false m 0 false (base ^ p) 0).
-  set (f := 
-      match Fdiv_core_binary prec (Z.pos m) 0 (Z.pos (base ^ p)) 0 with
-      | (0, _, _) => F754_nan false 1
-      | (Z.pos mz0, ez, lz) =>
-          binary_round_aux prec emax mode_NE (xorb false false) mz0 ez lz
-      | (Z.neg _, _, _) => F754_nan false 1
-      end).
-  fold emin; fold fexp. rewrite ! A. unfold cond_Zopp. 
-  assert (B: (Z2R (Z.pos m) / Z2R (Z.pos (base ^ p)) =
-              Z2R (Z.pos m) * bpow base_r (Z.neg p))%R).
-  { change (Z.neg p) with (- (Z.pos p)). rewrite bpow_opp. unfold Rdiv. f_equal. f_equal. 
-    unfold bpow. f_equal. simpl. apply Pos2Z.inj_pow_pos. }
-  rewrite ! B. 
-  destruct (Rlt_bool
-     (Rabs
-        (round radix2 fexp (round_mode mode_NE)
-           (Z2R (Z.pos m) * bpow base_r (Z.neg p)))) (bpow radix2 emax)).
-  + intros (P & Q & R & S).
-    split. rewrite B2R_FF2B. exact Q. 
-    split. rewrite is_finite_FF2B. auto.
-    rewrite Bsign_FF2B. auto. 
-  + intros (P & Q). rewrite B2FF_FF2B. auto.
-Qed.
-
-End FP_PARSING.
-
 Local Notation __ := (refl_equal Datatypes.Lt).
 
 Local Hint Extern 1 (Prec_gt_0 _) => exact (refl_equal Datatypes.Lt).
@@ -339,7 +245,7 @@ Definition of_longu (n:int64): float:= (**r conversion from unsigned 64-bit int
   BofZ 53 1024 __ __ (Int64.unsigned n).
 
 Definition from_parsed (base:positive) (intPart:positive) (expPart:Z) : float :=
-  build_from_parsed 53 1024 __ __ base intPart expPart.
+  Bparse 53 1024 __ __ base intPart expPart.
 
 (** Conversions between floats and their concrete in-memory representation
     as a sequence of 64 bits. *)
@@ -1004,7 +910,7 @@ Definition of_longu (n:int64): float32 := (**r conversion from unsigned 64-bit i
   BofZ 24 128 __ __ (Int64.unsigned n).
 
 Definition from_parsed (base:positive) (intPart:positive) (expPart:Z) : float32 :=
-  build_from_parsed 24 128 __ __ base intPart expPart.
+  Bparse 24 128 __ __ base intPart expPart.
 
 (** Conversions between floats and their concrete in-memory representation
     as a sequence of 32 bits. *)