diff options
| -rw-r--r-- | cfrontend/C2C.ml | 70 | ||||
| -rw-r--r-- | lib/Fappli_IEEE_extra.v | 218 | ||||
| -rw-r--r-- | lib/Floats.v | 98 |
3 files changed, 249 insertions, 137 deletions
diff --git a/cfrontend/C2C.ml b/cfrontend/C2C.ml index 389fc62d..b919c1d4 100644 --- a/cfrontend/C2C.ml +++ b/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 *) diff --git a/lib/Fappli_IEEE_extra.v b/lib/Fappli_IEEE_extra.v index 92ed11ae..3de7b103 100644 --- a/lib/Fappli_IEEE_extra.v +++ b/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 *) diff --git a/lib/Floats.v b/lib/Floats.v index f86632b9..e893e3e7 100644 --- a/lib/Floats.v +++ b/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. *) |