From d4513f41c54869c9b4ba96ab79d50933a1d5c25a Mon Sep 17 00:00:00 2001 From: Guillaume Melquiond Date: Mon, 28 Dec 2020 15:53:36 +0100 Subject: Update Flocq to 3.4.0 (#383) --- flocq/IEEE754/Binary.v | 103 +++++++--- flocq/IEEE754/Bits.v | 138 ++++++++----- flocq/IEEE754/SpecFloatCompat.v | 435 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 603 insertions(+), 73 deletions(-) create mode 100644 flocq/IEEE754/SpecFloatCompat.v (limited to 'flocq/IEEE754') diff --git a/flocq/IEEE754/Binary.v b/flocq/IEEE754/Binary.v index ac38c761..35d15cb3 100644 --- a/flocq/IEEE754/Binary.v +++ b/flocq/IEEE754/Binary.v @@ -627,6 +627,52 @@ Proof. now rewrite Pcompare_antisym. Qed. +Theorem bounded_le_emax_minus_prec : + forall mx ex, + bounded mx ex = true -> + (F2R (Float radix2 (Zpos mx) ex) + <= bpow radix2 emax - bpow radix2 (emax - prec))%R. +Proof. +intros mx ex Hx. +destruct (andb_prop _ _ Hx) as (H1,H2). +generalize (Zeq_bool_eq _ _ H1). clear H1. intro H1. +generalize (Zle_bool_imp_le _ _ H2). clear H2. intro H2. +generalize (mag_F2R_Zdigits radix2 (Zpos mx) ex). +destruct (mag radix2 (F2R (Float radix2 (Zpos mx) ex))) as (e',Ex). +unfold mag_val. +intros H. +elim Ex; [|now apply Rgt_not_eq, F2R_gt_0]; intros _. +rewrite <-F2R_Zabs; simpl; clear Ex; intros Ex. +generalize (Rmult_lt_compat_r (bpow radix2 (-ex)) _ _ (bpow_gt_0 _ _) Ex). +unfold F2R; simpl; rewrite Rmult_assoc, <-!bpow_plus. +rewrite H; [|intro H'; discriminate H']. +rewrite <-Z.add_assoc, Z.add_opp_diag_r, Z.add_0_r, Rmult_1_r. +rewrite <-(IZR_Zpower _ _ (Zdigits_ge_0 _ _)); clear Ex; intro Ex. +generalize (Zlt_le_succ _ _ (lt_IZR _ _ Ex)); clear Ex; intro Ex. +generalize (IZR_le _ _ Ex). +rewrite succ_IZR; clear Ex; intro Ex. +generalize (Rplus_le_compat_r (-1) _ _ Ex); clear Ex; intro Ex. +ring_simplify in Ex; revert Ex. +rewrite (IZR_Zpower _ _ (Zdigits_ge_0 _ _)); intro Ex. +generalize (Rmult_le_compat_r (bpow radix2 ex) _ _ (bpow_ge_0 _ _) Ex). +intro H'; apply (Rle_trans _ _ _ H'). +rewrite Rmult_minus_distr_r, Rmult_1_l, <-bpow_plus. +revert H1; unfold fexp, FLT_exp; intro H1. +generalize (Z.le_max_l (Z.pos (digits2_pos mx) + ex - prec) emin). +rewrite H1; intro H1'. +generalize (proj1 (Z.le_sub_le_add_r _ _ _) H1'). +rewrite Zpos_digits2_pos; clear H1'; intro H1'. +apply (Rle_trans _ _ _ (Rplus_le_compat_r _ _ _ (bpow_le _ _ _ H1'))). +replace emax with (emax - prec - ex + (ex + prec))%Z at 1 by ring. +replace (emax - prec)%Z with (emax - prec - ex + ex)%Z at 2 by ring. +do 2 rewrite (bpow_plus _ (emax - prec - ex)). +rewrite <-Rmult_minus_distr_l. +rewrite <-(Rmult_1_l (_ + _)). +apply Rmult_le_compat_r. +{ apply Rle_0_minus, bpow_le; unfold Prec_gt_0 in prec_gt_0_; lia. } +change 1%R with (bpow radix2 0); apply bpow_le; lia. +Qed. + Theorem bounded_lt_emax : forall mx ex, bounded mx ex = true -> @@ -651,7 +697,7 @@ rewrite H. 2: discriminate. revert H1. clear -H2. rewrite Zpos_digits2_pos. unfold fexp, FLT_exp. -intros ; zify ; omega. +intros ; zify ; lia. Qed. Theorem bounded_ge_emin : @@ -679,7 +725,18 @@ unfold fexp, FLT_exp. clear -prec_gt_0_. unfold Prec_gt_0 in prec_gt_0_. clearbody emin. -intros ; zify ; omega. +intros ; zify ; lia. +Qed. + +Theorem abs_B2R_le_emax_minus_prec : + forall x, + (Rabs (B2R x) <= bpow radix2 emax - bpow radix2 (emax - prec))%R. +Proof. +intros [sx|sx|sx plx Hx|sx mx ex Hx] ; simpl ; + [rewrite Rabs_R0 ; apply Rle_0_minus, bpow_le ; + revert prec_gt_0_; unfold Prec_gt_0; lia..|]. +rewrite <- F2R_Zabs, abs_cond_Zopp. +now apply bounded_le_emax_minus_prec. Qed. Theorem abs_B2R_lt_emax : @@ -728,7 +785,7 @@ rewrite Cx. unfold cexp, fexp, FLT_exp. destruct (mag radix2 (F2R (Float radix2 (Zpos mx) ex))) as (e',Ex). simpl. apply Z.max_lub. -cut (e' - 1 < emax)%Z. clear ; omega. +cut (e' - 1 < emax)%Z. clear ; lia. apply lt_bpow with radix2. apply Rle_lt_trans with (2 := Bx). change (Zpos mx) with (Z.abs (Zpos mx)). @@ -738,7 +795,7 @@ apply Rgt_not_eq. now apply F2R_gt_0. unfold emin. generalize (prec_gt_0 prec). -clear -Hmax ; omega. +clear -Hmax ; lia. Qed. (** Truncation *) @@ -889,7 +946,7 @@ now inversion H. (* *) intros p Hp. assert (He: (e <= fexp (Zdigits radix2 m + e))%Z). -clear -Hp ; zify ; omega. +clear -Hp ; zify ; lia. destruct (inbetween_float_ex radix2 m e l) as (x, Hx). generalize (inbetween_shr x m e l (fexp (Zdigits radix2 m + e) - e) Hm Hx). assert (Hx0 : (0 <= x)%R). @@ -1091,18 +1148,18 @@ rewrite Zpos_digits2_pos. replace (Zdigits radix2 (Zpos (match (Zpower 2 prec - 1)%Z with Zpos p => p | _ => xH end))) with prec. unfold fexp, FLT_exp, emin. generalize (prec_gt_0 prec). -clear -Hmax ; zify ; omega. +clear -Hmax ; zify ; lia. change 2%Z with (radix_val radix2). case_eq (Zpower radix2 prec - 1)%Z. simpl Zdigits. generalize (Zpower_gt_1 radix2 prec (prec_gt_0 prec)). -clear ; omega. +clear ; lia. intros p Hp. apply Zle_antisym. -cut (prec - 1 < Zdigits radix2 (Zpos p))%Z. clear ; omega. +cut (prec - 1 < Zdigits radix2 (Zpos p))%Z. clear ; lia. apply Zdigits_gt_Zpower. simpl Z.abs. rewrite <- Hp. -cut (Zpower radix2 (prec - 1) < Zpower radix2 prec)%Z. clear ; omega. +cut (Zpower radix2 (prec - 1) < Zpower radix2 prec)%Z. clear ; lia. apply lt_IZR. rewrite 2!IZR_Zpower. 2: now apply Zlt_le_weak. apply bpow_lt. @@ -1113,7 +1170,7 @@ simpl Z.abs. rewrite <- Hp. apply Zlt_pred. intros p Hp. generalize (Zpower_gt_1 radix2 _ (prec_gt_0 prec)). -clear -Hp ; zify ; omega. +clear -Hp ; zify ; lia. apply Rnot_lt_le. intros Hx. generalize (refl_equal (bounded m2 e2)). @@ -1271,18 +1328,18 @@ rewrite Zpos_digits2_pos. replace (Zdigits radix2 (Zpos (match (Zpower 2 prec - 1)%Z with Zpos p => p | _ => xH end))) with prec. unfold fexp, FLT_exp, emin. generalize (prec_gt_0 prec). -clear -Hmax ; zify ; omega. +clear -Hmax ; zify ; lia. change 2%Z with (radix_val radix2). case_eq (Zpower radix2 prec - 1)%Z. simpl Zdigits. generalize (Zpower_gt_1 radix2 prec (prec_gt_0 prec)). -clear ; omega. +clear ; lia. intros p Hp. apply Zle_antisym. -cut (prec - 1 < Zdigits radix2 (Zpos p))%Z. clear ; omega. +cut (prec - 1 < Zdigits radix2 (Zpos p))%Z. clear ; lia. apply Zdigits_gt_Zpower. simpl Z.abs. rewrite <- Hp. -cut (Zpower radix2 (prec - 1) < Zpower radix2 prec)%Z. clear ; omega. +cut (Zpower radix2 (prec - 1) < Zpower radix2 prec)%Z. clear ; lia. apply lt_IZR. rewrite 2!IZR_Zpower. 2: now apply Zlt_le_weak. apply bpow_lt. @@ -1293,7 +1350,7 @@ simpl Z.abs. rewrite <- Hp. apply Zlt_pred. intros p Hp. generalize (Zpower_gt_1 radix2 _ (prec_gt_0 prec)). -clear -Hp ; zify ; omega. +clear -Hp ; zify ; lia. apply Rnot_lt_le. intros Hx. generalize (refl_equal (bounded m2 e2)). @@ -1370,7 +1427,7 @@ clear -Hmax. unfold emin. intros dx dy dxy Hx Hy Hxy. zify ; intros ; subst. -omega. +lia. (* *) case sx ; case sy. apply Rlt_bool_false. @@ -1479,7 +1536,7 @@ case_eq (ex' - ex)%Z ; simpl. intros H. now rewrite Zminus_eq with (1 := H). intros p. -clear -He ; zify ; omega. +clear -He ; zify ; lia. intros. apply refl_equal. Qed. @@ -1580,7 +1637,7 @@ now rewrite is_finite_FF2B. rewrite Bsign_FF2B, Rz''. rewrite Rcompare_Gt... apply F2R_gt_0. -simpl. zify; omega. +simpl. zify; lia. intros Hz' (Vz, Rz). rewrite B2FF_FF2B, Rz. apply f_equal. @@ -1599,7 +1656,7 @@ now rewrite is_finite_FF2B. rewrite Bsign_FF2B, Rz''. rewrite Rcompare_Lt... apply F2R_lt_0. -simpl. zify; omega. +simpl. zify; lia. intros Hz' (Vz, Rz). rewrite B2FF_FF2B, Rz. apply f_equal. @@ -2150,7 +2207,7 @@ set (e' := Z.min _ _). assert (2 * e' <= ex)%Z as He. { assert (e' <= Z.div2 ex)%Z by apply Z.le_min_r. rewrite (Zdiv2_odd_eqn ex). - destruct Z.odd ; omega. } + destruct Z.odd ; lia. } generalize (Fsqrt_core_correct radix2 (Zpos mx) ex e' eq_refl He). unfold Fsqrt_core. set (mx' := match (ex - 2 * e')%Z with Z0 => _ | _ => _ end). @@ -2187,7 +2244,7 @@ apply Rlt_le_trans with (1 := Heps). fold (bpow radix2 0). apply bpow_le. generalize (prec_gt_0 prec). -clear ; omega. +clear ; lia. apply Rsqr_incrst_0. 3: apply bpow_ge_0. rewrite Rsqr_mult. @@ -2211,7 +2268,7 @@ now apply IZR_le. change 4%R with (bpow radix2 2). apply bpow_le. generalize (prec_gt_0 prec). -clear -Hmax ; omega. +clear -Hmax ; lia. apply Rmult_le_pos. apply sqrt_ge_0. rewrite <- (Rplus_opp_r 1). @@ -2230,7 +2287,7 @@ unfold Rsqr. rewrite <- bpow_plus. apply bpow_le. unfold emin. -clear -Hmax ; omega. +clear -Hmax ; lia. apply generic_format_ge_bpow with fexp. intros. apply Z.le_max_r. diff --git a/flocq/IEEE754/Bits.v b/flocq/IEEE754/Bits.v index 3a84edfe..68bc541a 100644 --- a/flocq/IEEE754/Bits.v +++ b/flocq/IEEE754/Bits.v @@ -18,6 +18,8 @@ COPYING file for more details. *) (** * IEEE-754 encoding of binary floating-point data *) + +From Coq Require Import Lia. Require Import Core Digits Binary. Section Binary_Bits. @@ -43,10 +45,10 @@ Proof. intros s m e Hm He. assert (0 <= mw)%Z as Hmw. destruct mw as [|mw'|mw'] ; try easy. - clear -Hm ; simpl in Hm ; omega. + clear -Hm ; simpl in Hm ; lia. assert (0 <= ew)%Z as Hew. destruct ew as [|ew'|ew'] ; try easy. - clear -He ; simpl in He ; omega. + clear -He ; simpl in He ; lia. unfold join_bits. rewrite Z.shiftl_mul_pow2 by easy. split. @@ -54,9 +56,9 @@ split. rewrite <- (Zmult_0_l (2^mw)). apply Zmult_le_compat_r. case s. - clear -He ; omega. + clear -He ; lia. now rewrite Zmult_0_l. - clear -Hm ; omega. + clear -Hm ; lia. - apply Z.lt_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. @@ -65,9 +67,9 @@ split. apply Zmult_le_compat_r. rewrite Zpower_plus by easy. change (2^1)%Z with 2%Z. - case s ; clear -He ; omega. - clear -Hm ; omega. - clear -Hew ; omega. + case s ; clear -He ; lia. + clear -Hm ; lia. + clear -Hew ; lia. easy. Qed. @@ -85,10 +87,10 @@ Proof. intros s m e Hm He. assert (0 <= mw)%Z as Hmw. destruct mw as [|mw'|mw'] ; try easy. - clear -Hm ; simpl in Hm ; omega. + clear -Hm ; simpl in Hm ; lia. assert (0 <= ew)%Z as Hew. destruct ew as [|ew'|ew'] ; try easy. - clear -He ; simpl in He ; omega. + clear -He ; simpl in He ; lia. unfold split_bits, join_bits. rewrite Z.shiftl_mul_pow2 by easy. apply f_equal2 ; [apply f_equal2|]. @@ -99,7 +101,7 @@ apply f_equal2 ; [apply f_equal2|]. apply Zplus_le_0_compat. apply Zmult_le_0_compat. apply He. - clear -Hm ; omega. + clear -Hm ; lia. apply Hm. + apply Zle_bool_false. apply Zplus_lt_reg_l with (2^mw * (-e))%Z. @@ -108,12 +110,12 @@ apply f_equal2 ; [apply f_equal2|]. apply Z.lt_le_trans with (2^mw * 1)%Z. now apply Zmult_lt_compat_r. apply Zmult_le_compat_l. - clear -He ; omega. - clear -Hm ; omega. + clear -He ; lia. + clear -Hm ; lia. - rewrite Zplus_comm. rewrite Z_mod_plus_full. now apply Zmod_small. -- rewrite Z_div_plus_full_l by (clear -Hm ; omega). +- rewrite Z_div_plus_full_l by (clear -Hm ; lia). rewrite Zdiv_small with (1 := Hm). rewrite Zplus_0_r. case s. @@ -175,7 +177,7 @@ 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. +cut (x / (2^mw * 2^ew) < 2)%Z. clear ; lia. apply Zdiv_lt_upper_bound. now apply Zmult_lt_0_compat. rewrite <- Zpower_exp ; try ( apply Z.le_ge ; apply Zlt_le_weak ; assumption ). @@ -244,8 +246,8 @@ Theorem split_bits_of_binary_float_correct : split_bits (bits_of_binary_float x) = split_bits_of_binary_float x. Proof. intros [sx|sx|sx plx Hplx|sx mx ex Hx] ; - try ( simpl ; apply split_join_bits ; split ; try apply Z.le_refl ; try apply Zlt_pred ; trivial ; omega ). -simpl. apply split_join_bits; split; try (zify; omega). + try ( simpl ; apply split_join_bits ; split ; try apply Z.le_refl ; try apply Zlt_pred ; trivial ; lia ). +simpl. apply split_join_bits; split; try (zify; lia). destruct (digits2_Pnat_correct plx). unfold nan_pl in Hplx. rewrite Zpos_digits2_pos, <- Z_of_nat_S_digits2_Pnat in Hplx. @@ -253,7 +255,7 @@ rewrite Zpower_nat_Z in H0. eapply Z.lt_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. +unfold prec in *. zify; lia. (* *) unfold bits_of_binary_float, split_bits_of_binary_float. assert (Hf: (emin <= ex /\ Zdigits radix2 (Zpos mx) <= prec)%Z). @@ -263,14 +265,14 @@ rewrite Zpos_digits2_pos in Hx'. generalize (Zeq_bool_eq _ _ Hx'). unfold FLT_exp. unfold emin. -clear ; zify ; omega. +clear ; zify ; lia. case Zle_bool_spec ; intros H ; [ apply -> Z.le_0_sub in H | apply -> Z.lt_sub_0 in H ] ; apply split_join_bits ; try now split. (* *) split. -clear -He_gt_0 H ; omega. -cut (Zpos mx < 2 * 2^mw)%Z. clear ; omega. +clear -He_gt_0 H ; lia. +cut (Zpos mx < 2 * 2^mw)%Z. clear ; lia. replace (2 * 2^mw)%Z with (2^prec)%Z. apply (Zpower_gt_Zdigits radix2 _ (Zpos mx)). apply Hf. @@ -282,12 +284,12 @@ now apply Zlt_le_weak. (* *) split. generalize (proj1 Hf). -clear ; omega. +clear ; lia. destruct (andb_prop _ _ Hx) as (_, Hx'). unfold emin. replace (2^ew)%Z with (2 * emax)%Z. generalize (Zle_bool_imp_le _ _ Hx'). -clear ; omega. +clear ; lia. apply sym_eq. rewrite (Zsucc_pred ew). unfold Z.succ. @@ -305,7 +307,7 @@ intros [sx|sx|sx pl pl_range|sx mx ex H]. - apply join_bits_range ; now split. - apply join_bits_range. now split. - clear -He_gt_0 ; omega. + clear -He_gt_0 ; lia. - apply Z.ltb_lt in pl_range. apply join_bits_range. split. @@ -313,7 +315,7 @@ intros [sx|sx|sx pl pl_range|sx mx ex H]. apply (Zpower_gt_Zdigits radix2 _ (Zpos pl)). apply Z.lt_succ_r. now rewrite <- Zdigits2_Zdigits. - clear -He_gt_0 ; omega. + clear -He_gt_0 ; lia. - unfold bounded in H. apply Bool.andb_true_iff in H ; destruct H as [A B]. apply Z.leb_le in B. @@ -321,22 +323,22 @@ intros [sx|sx|sx pl pl_range|sx mx ex H]. case Zle_bool_spec ; intros H. + apply join_bits_range. * split. - clear -H ; omega. + clear -H ; lia. rewrite Zpos_digits2_pos in A. cut (Zpos mx < 2 ^ prec)%Z. unfold prec. - rewrite Zpower_plus by (clear -Hmw ; omega). + rewrite Zpower_plus by (clear -Hmw ; lia). change (2^1)%Z with 2%Z. - clear ; omega. + clear ; lia. apply (Zpower_gt_Zdigits radix2 _ (Zpos mx)). - clear -A ; zify ; omega. + clear -A ; zify ; lia. * split. - unfold emin ; clear -A ; zify ; omega. + unfold emin ; clear -A ; zify ; lia. replace ew with ((ew - 1) + 1)%Z by ring. - rewrite Zpower_plus by (clear - Hew ; omega). + rewrite Zpower_plus by (clear - Hew ; lia). unfold emin, emax in *. change (2^1)%Z with 2%Z. - clear -B ; omega. + clear -B ; lia. + apply -> Z.lt_sub_0 in H. apply join_bits_range ; now split. Qed. @@ -370,7 +372,7 @@ unfold binary_float_of_bits_aux, split_bits. assert (Hnan: nan_pl prec 1 = true). apply Z.ltb_lt. simpl. unfold prec. - clear -Hmw ; omega. + clear -Hmw ; lia. case Zeq_bool_spec ; intros He1. case_eq (x mod 2^mw)%Z ; try easy. (* subnormal *) @@ -389,7 +391,7 @@ unfold Fexp, FLT_exp. apply sym_eq. apply Zmax_right. clear -H Hprec. -unfold prec ; omega. +unfold prec ; lia. apply Rnot_le_lt. intros H0. refine (_ (mag_le radix2 _ _ _ H0)). @@ -397,20 +399,20 @@ rewrite mag_bpow. rewrite mag_F2R_Zdigits. 2: discriminate. unfold emin, prec. apply Zlt_not_le. -cut (0 < emax)%Z. clear -H Hew ; omega. +cut (0 < emax)%Z. clear -H Hew ; lia. apply (Zpower_gt_0 radix2). -clear -Hew ; omega. +clear -Hew ; lia. apply bpow_gt_0. case Zeq_bool_spec ; intros He2. case_eq (x mod 2 ^ mw)%Z; try easy. (* nan *) intros plx Eqplx. apply Z.ltb_lt. rewrite Zpos_digits2_pos. -assert (forall a b, a <= b -> a < b+1)%Z by (intros; omega). apply H. clear H. +assert (forall a b, a <= b -> a < b+1)%Z by (intros; lia). 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. +apply Z.lt_gt, Zpower_gt_0. lia. case_eq (x mod 2^mw + 2^mw)%Z ; try easy. (* normal *) intros px Hm. @@ -452,7 +454,7 @@ revert He1. fold ex. cut (0 <= ex)%Z. unfold emin. -clear ; intros H1 H2 ; omega. +clear ; intros H1 H2 ; lia. eapply Z_mod_lt. apply Z.lt_gt. apply (Zpower_gt_0 radix2). @@ -471,12 +473,12 @@ 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. +clear ; intros H1 H2 ; lia. replace ew with (1 + (ew - 1))%Z by ring. rewrite Zpower_exp. apply refl_equal. discriminate. -clear -Hew ; omega. +clear -Hew ; lia. eapply Z_mod_lt. apply Z.lt_gt. apply (Zpower_gt_0 radix2). @@ -503,13 +505,13 @@ apply refl_equal. simpl. rewrite Zeq_bool_false. now rewrite Zeq_bool_true. -cut (1 < 2^ew)%Z. clear ; omega. +cut (1 < 2^ew)%Z. clear ; lia. now apply (Zpower_gt_1 radix2). (* *) simpl. rewrite Zeq_bool_false. rewrite Zeq_bool_true; auto. -cut (1 < 2^ew)%Z. clear ; omega. +cut (1 < 2^ew)%Z. clear ; lia. now apply (Zpower_gt_1 radix2). (* *) unfold split_bits_of_binary_float. @@ -522,19 +524,19 @@ destruct (andb_prop _ _ Bx) as (_, H1). generalize (Zle_bool_imp_le _ _ H1). unfold emin. replace (2^ew)%Z with (2 * emax)%Z. -clear ; omega. +clear ; lia. replace ew with (1 + (ew - 1))%Z by ring. rewrite Zpower_exp. apply refl_equal. discriminate. -clear -Hew ; omega. +clear -Hew ; lia. destruct (andb_prop _ _ Bx) as (H1, _). generalize (Zeq_bool_eq _ _ H1). rewrite Zpos_digits2_pos. unfold FLT_exp, emin. generalize (Zdigits radix2 (Zpos mx)). clear. -intros ; zify ; omega. +intros ; zify ; lia. (* . *) rewrite Zeq_bool_true. 2: apply refl_equal. simpl. @@ -547,7 +549,7 @@ apply -> Z.lt_sub_0 in Hm. generalize (Zdigits_le_Zpower radix2 _ (Zpos mx) Hm). generalize (Zdigits radix2 (Zpos mx)). clear. -intros ; zify ; omega. +intros ; zify ; lia. Qed. Theorem bits_of_binary_float_of_bits : @@ -588,12 +590,12 @@ case Zeq_bool_spec ; intros He2. case_eq mx; intros Hm. now rewrite He2. now rewrite He2. -intros. zify; omega. +intros. zify; lia. (* normal *) case_eq (mx + 2 ^ mw)%Z. intros Hm. apply False_ind. -clear -Bm Hm ; omega. +clear -Bm Hm ; lia. intros p Hm Jx Cx. rewrite <- Hm. rewrite Zle_bool_true. @@ -601,7 +603,7 @@ now ring_simplify (mx + 2^mw - 2^mw)%Z (ex + emin - 1 - emin + 1)%Z. now ring_simplify. intros p Hm. apply False_ind. -clear -Bm Hm ; zify ; omega. +clear -Bm Hm ; zify ; lia. Qed. End Binary_Bits. @@ -623,6 +625,12 @@ Proof. apply refl_equal. Qed. +Let Hemax : (3 <= 128)%Z. +Proof. +intros H. +discriminate H. +Qed. + Definition default_nan_pl32 : { nan : binary32 | is_nan 24 128 nan = true } := exist _ (@B754_nan 24 128 false (iter_nat xO 22 xH) (refl_equal true)) (refl_equal true). @@ -639,16 +647,28 @@ Definition binop_nan_pl32 (f1 f2 : binary32) : { nan : binary32 | is_nan 24 128 | _, _ => default_nan_pl32 end. +Definition ternop_nan_pl32 (f1 f2 f3 : binary32) : { nan : binary32 | is_nan 24 128 nan = true } := + match f1, f2, f3 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) + | _, _, B754_nan s3 pl3 Hpl3 => exist _ (B754_nan s3 pl3 Hpl3) (refl_equal true) + | _, _, _ => default_nan_pl32 + end. + Definition b32_erase : binary32 -> binary32 := erase 24 128. Definition b32_opp : binary32 -> binary32 := Bopp 24 128 unop_nan_pl32. Definition b32_abs : binary32 -> binary32 := Babs 24 128 unop_nan_pl32. -Definition b32_sqrt : mode -> binary32 -> binary32 := Bsqrt _ _ Hprec Hprec_emax unop_nan_pl32. +Definition b32_pred : binary32 -> binary32 := Bpred _ _ Hprec Hprec_emax Hemax unop_nan_pl32. +Definition b32_succ : binary32 -> binary32 := Bsucc _ _ Hprec Hprec_emax Hemax unop_nan_pl32. +Definition b32_sqrt : mode -> binary32 -> binary32 := Bsqrt _ _ Hprec Hprec_emax unop_nan_pl32. Definition b32_plus : mode -> binary32 -> binary32 -> binary32 := Bplus _ _ Hprec Hprec_emax binop_nan_pl32. Definition b32_minus : mode -> binary32 -> binary32 -> binary32 := Bminus _ _ Hprec Hprec_emax binop_nan_pl32. Definition b32_mult : mode -> binary32 -> binary32 -> binary32 := Bmult _ _ Hprec Hprec_emax binop_nan_pl32. Definition b32_div : mode -> binary32 -> binary32 -> binary32 := Bdiv _ _ Hprec Hprec_emax binop_nan_pl32. +Definition b32_fma : mode -> binary32 -> binary32 -> binary32 -> binary32 := Bfma _ _ Hprec Hprec_emax ternop_nan_pl32. + Definition b32_compare : binary32 -> binary32 -> option comparison := Bcompare 24 128. 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. @@ -672,6 +692,12 @@ Proof. apply refl_equal. Qed. +Let Hemax : (3 <= 1024)%Z. +Proof. +intros H. +discriminate H. +Qed. + Definition default_nan_pl64 : { nan : binary64 | is_nan 53 1024 nan = true } := exist _ (@B754_nan 53 1024 false (iter_nat xO 51 xH) (refl_equal true)) (refl_equal true). @@ -688,9 +714,19 @@ Definition binop_nan_pl64 (f1 f2 : binary64) : { nan : binary64 | is_nan 53 1024 | _, _ => default_nan_pl64 end. +Definition ternop_nan_pl64 (f1 f2 f3 : binary64) : { nan : binary64 | is_nan 53 1024 nan = true } := + match f1, f2, f3 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) + | _, _, B754_nan s3 pl3 Hpl3 => exist _ (B754_nan s3 pl3 Hpl3) (refl_equal true) + | _, _, _ => default_nan_pl64 + end. + Definition b64_erase : binary64 -> binary64 := erase 53 1024. Definition b64_opp : binary64 -> binary64 := Bopp 53 1024 unop_nan_pl64. Definition b64_abs : binary64 -> binary64 := Babs 53 1024 unop_nan_pl64. +Definition b64_pred : binary64 -> binary64 := Bpred _ _ Hprec Hprec_emax Hemax unop_nan_pl64. +Definition b64_succ : binary64 -> binary64 := Bsucc _ _ Hprec Hprec_emax Hemax unop_nan_pl64. Definition b64_sqrt : mode -> binary64 -> binary64 := Bsqrt _ _ Hprec Hprec_emax unop_nan_pl64. Definition b64_plus : mode -> binary64 -> binary64 -> binary64 := Bplus _ _ Hprec Hprec_emax binop_nan_pl64. @@ -698,6 +734,8 @@ Definition b64_minus : mode -> binary64 -> binary64 -> binary64 := Bminus _ _ Hp Definition b64_mult : mode -> binary64 -> binary64 -> binary64 := Bmult _ _ Hprec Hprec_emax binop_nan_pl64. Definition b64_div : mode -> binary64 -> binary64 -> binary64 := Bdiv _ _ Hprec Hprec_emax binop_nan_pl64. +Definition b64_fma : mode -> binary64 -> binary64 -> binary64 -> binary64 := Bfma _ _ Hprec Hprec_emax ternop_nan_pl64. + Definition b64_compare : binary64 -> binary64 -> option comparison := Bcompare 53 1024. 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. diff --git a/flocq/IEEE754/SpecFloatCompat.v b/flocq/IEEE754/SpecFloatCompat.v new file mode 100644 index 00000000..e2ace4d5 --- /dev/null +++ b/flocq/IEEE754/SpecFloatCompat.v @@ -0,0 +1,435 @@ +(** +This file is part of the Flocq formalization of floating-point +arithmetic in Coq: http://flocq.gforge.inria.fr/ + +Copyright (C) 2018-2019 Guillaume Bertholon +#
# +Copyright (C) 2018-2019 Érik Martin-Dorel +#
# +Copyright (C) 2018-2019 Pierre Roux + +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. +*) + +Require Import ZArith. + +(** ** Inductive specification of floating-point numbers + +Similar to [IEEE754.Binary.full_float], but with no NaN payload. *) +Variant spec_float := + | S754_zero (s : bool) + | S754_infinity (s : bool) + | S754_nan + | S754_finite (s : bool) (m : positive) (e : Z). + +(** ** Parameterized definitions + +[prec] is the number of bits of the mantissa including the implicit one; +[emax] is the exponent of the infinities. + +For instance, Binary64 is defined by [prec = 53] and [emax = 1024]. *) +Section FloatOps. + Variable prec emax : Z. + + Definition emin := (3-emax-prec)%Z. + Definition fexp e := Z.max (e - prec) emin. + + Section Zdigits2. + Fixpoint digits2_pos (n : positive) : positive := + match n with + | xH => xH + | xO p => Pos.succ (digits2_pos p) + | xI p => Pos.succ (digits2_pos p) + end. + + Definition Zdigits2 n := + match n with + | Z0 => n + | Zpos p => Zpos (digits2_pos p) + | Zneg p => Zpos (digits2_pos p) + end. + End Zdigits2. + + Section ValidBinary. + Definition canonical_mantissa m e := + Zeq_bool (fexp (Zpos (digits2_pos m) + e)) e. + + Definition bounded m e := + andb (canonical_mantissa m e) (Zle_bool e (emax - prec)). + + Definition valid_binary x := + match x with + | S754_finite _ m e => bounded m e + | _ => true + end. + End ValidBinary. + + Section Iter. + Context {A : Type}. + Variable (f : A -> A). + + Fixpoint iter_pos (n : positive) (x : A) {struct n} : A := + match n with + | xI n' => iter_pos n' (iter_pos n' (f x)) + | xO n' => iter_pos n' (iter_pos n' x) + | xH => f x + end. + End Iter. + + Section Rounding. + Inductive location := loc_Exact | loc_Inexact : comparison -> location. + + Record shr_record := { shr_m : Z ; shr_r : bool ; shr_s : bool }. + + Definition shr_1 mrs := + let '(Build_shr_record m r s) := mrs in + let s := orb r s in + match m with + | Z0 => Build_shr_record Z0 false s + | Zpos xH => Build_shr_record Z0 true s + | Zpos (xO p) => Build_shr_record (Zpos p) false s + | Zpos (xI p) => Build_shr_record (Zpos p) true s + | Zneg xH => Build_shr_record Z0 true s + | Zneg (xO p) => Build_shr_record (Zneg p) false s + | Zneg (xI p) => Build_shr_record (Zneg p) true s + end. + + Definition loc_of_shr_record mrs := + match mrs with + | Build_shr_record _ false false => loc_Exact + | Build_shr_record _ false true => loc_Inexact Lt + | Build_shr_record _ true false => loc_Inexact Eq + | Build_shr_record _ true true => loc_Inexact Gt + end. + + Definition shr_record_of_loc m l := + match l with + | loc_Exact => Build_shr_record m false false + | loc_Inexact Lt => Build_shr_record m false true + | loc_Inexact Eq => Build_shr_record m true false + | loc_Inexact Gt => Build_shr_record m true true + end. + + Definition shr mrs e n := + match n with + | Zpos p => (iter_pos shr_1 p mrs, (e + n)%Z) + | _ => (mrs, e) + end. + + Definition shr_fexp m e l := + shr (shr_record_of_loc m l) e (fexp (Zdigits2 m + e) - e). + + Definition round_nearest_even mx lx := + match lx with + | loc_Exact => mx + | loc_Inexact Lt => mx + | loc_Inexact Eq => if Z.even mx then mx else (mx + 1)%Z + | loc_Inexact Gt => (mx + 1)%Z + end. + + Definition binary_round_aux sx mx ex lx := + let '(mrs', e') := shr_fexp mx ex lx in + let '(mrs'', e'') := shr_fexp (round_nearest_even (shr_m mrs') (loc_of_shr_record mrs')) e' loc_Exact in + match shr_m mrs'' with + | Z0 => S754_zero sx + | Zpos m => if Zle_bool e'' (emax - prec) then S754_finite sx m e'' else S754_infinity sx + | _ => S754_nan + end. + + Definition shl_align mx ex ex' := + match (ex' - ex)%Z with + | Zneg d => (shift_pos d mx, ex') + | _ => (mx, ex) + end. + + Definition binary_round sx mx ex := + let '(mz, ez) := shl_align mx ex (fexp (Zpos (digits2_pos mx) + ex))in + binary_round_aux sx (Zpos mz) ez loc_Exact. + + Definition binary_normalize m e szero := + match m with + | Z0 => S754_zero szero + | Zpos m => binary_round false m e + | Zneg m => binary_round true m e + end. + End Rounding. + + (** ** Define operations *) + + Definition SFopp x := + match x with + | S754_nan => S754_nan + | S754_infinity sx => S754_infinity (negb sx) + | S754_finite sx mx ex => S754_finite (negb sx) mx ex + | S754_zero sx => S754_zero (negb sx) + end. + + Definition SFabs x := + match x with + | S754_nan => S754_nan + | S754_infinity sx => S754_infinity false + | S754_finite sx mx ex => S754_finite false mx ex + | S754_zero sx => S754_zero false + end. + + Definition SFcompare f1 f2 := + match f1, f2 with + | S754_nan , _ | _, S754_nan => None + | S754_infinity s1, S754_infinity s2 => + Some match s1, s2 with + | true, true => Eq + | false, false => Eq + | true, false => Lt + | false, true => Gt + end + | S754_infinity s, _ => Some (if s then Lt else Gt) + | _, S754_infinity s => Some (if s then Gt else Lt) + | S754_finite s _ _, S754_zero _ => Some (if s then Lt else Gt) + | S754_zero _, S754_finite s _ _ => Some (if s then Gt else Lt) + | S754_zero _, S754_zero _ => Some Eq + | S754_finite s1 m1 e1, S754_finite s2 m2 e2 => + Some match s1, s2 with + | true, false => Lt + | false, true => Gt + | false, false => + match Z.compare e1 e2 with + | Lt => Lt + | Gt => Gt + | Eq => Pcompare m1 m2 Eq + end + | true, true => + match Z.compare e1 e2 with + | Lt => Gt + | Gt => Lt + | Eq => CompOpp (Pcompare m1 m2 Eq) + end + end + end. + + Definition SFeqb f1 f2 := + match SFcompare f1 f2 with + | Some Eq => true + | _ => false + end. + + Definition SFltb f1 f2 := + match SFcompare f1 f2 with + | Some Lt => true + | _ => false + end. + + Definition SFleb f1 f2 := + match SFcompare f1 f2 with + | Some (Lt | Eq) => true + | _ => false + end. + + Variant float_class : Set := + | PNormal | NNormal | PSubn | NSubn | PZero | NZero | PInf | NInf | NaN. + + Definition SFclassify f := + match f with + | S754_nan => NaN + | S754_infinity false => PInf + | S754_infinity true => NInf + | S754_zero false => NZero + | S754_zero true => PZero + | S754_finite false m _ => + if (digits2_pos m =? Z.to_pos prec)%positive then PNormal + else PSubn + | S754_finite true m _ => + if (digits2_pos m =? Z.to_pos prec)%positive then NNormal + else NSubn + end. + + Definition SFmul x y := + match x, y with + | S754_nan, _ | _, S754_nan => S754_nan + | S754_infinity sx, S754_infinity sy => S754_infinity (xorb sx sy) + | S754_infinity sx, S754_finite sy _ _ => S754_infinity (xorb sx sy) + | S754_finite sx _ _, S754_infinity sy => S754_infinity (xorb sx sy) + | S754_infinity _, S754_zero _ => S754_nan + | S754_zero _, S754_infinity _ => S754_nan + | S754_finite sx _ _, S754_zero sy => S754_zero (xorb sx sy) + | S754_zero sx, S754_finite sy _ _ => S754_zero (xorb sx sy) + | S754_zero sx, S754_zero sy => S754_zero (xorb sx sy) + | S754_finite sx mx ex, S754_finite sy my ey => + binary_round_aux (xorb sx sy) (Zpos (mx * my)) (ex + ey) loc_Exact + end. + + Definition cond_Zopp (b : bool) m := if b then Z.opp m else m. + + Definition SFadd x y := + match x, y with + | S754_nan, _ | _, S754_nan => S754_nan + | S754_infinity sx, S754_infinity sy => + if Bool.eqb sx sy then x else S754_nan + | S754_infinity _, _ => x + | _, S754_infinity _ => y + | S754_zero sx, S754_zero sy => + if Bool.eqb sx sy then x else + S754_zero false + | S754_zero _, _ => y + | _, S754_zero _ => x + | S754_finite sx mx ex, S754_finite sy my ey => + let ez := Z.min ex ey in + binary_normalize (Zplus (cond_Zopp sx (Zpos (fst (shl_align mx ex ez)))) (cond_Zopp sy (Zpos (fst (shl_align my ey ez))))) + ez false + end. + + Definition SFsub x y := + match x, y with + | S754_nan, _ | _, S754_nan => S754_nan + | S754_infinity sx, S754_infinity sy => + if Bool.eqb sx (negb sy) then x else S754_nan + | S754_infinity _, _ => x + | _, S754_infinity sy => S754_infinity (negb sy) + | S754_zero sx, S754_zero sy => + if Bool.eqb sx (negb sy) then x else + S754_zero false + | S754_zero _, S754_finite sy my ey => S754_finite (negb sy) my ey + | _, S754_zero _ => x + | S754_finite sx mx ex, S754_finite sy my ey => + let ez := Z.min ex ey in + binary_normalize (Zminus (cond_Zopp sx (Zpos (fst (shl_align mx ex ez)))) (cond_Zopp sy (Zpos (fst (shl_align my ey ez))))) + ez false + end. + + Definition new_location_even nb_steps k := + if Zeq_bool k 0 then loc_Exact + else loc_Inexact (Z.compare (2 * k) nb_steps). + + Definition new_location_odd nb_steps k := + if Zeq_bool k 0 then loc_Exact + else + loc_Inexact + match Z.compare (2 * k + 1) nb_steps with + | Lt => Lt + | Eq => Lt + | Gt => Gt + end. + + Definition new_location nb_steps := + if Z.even nb_steps then new_location_even nb_steps else new_location_odd nb_steps. + + Definition SFdiv_core_binary m1 e1 m2 e2 := + let d1 := Zdigits2 m1 in + let d2 := Zdigits2 m2 in + let e' := Z.min (fexp (d1 + e1 - (d2 + e2))) (e1 - e2) in + let s := (e1 - e2 - e')%Z in + let m' := + match s with + | Zpos _ => Z.shiftl m1 s + | Z0 => m1 + | Zneg _ => Z0 + end in + let '(q, r) := Z.div_eucl m' m2 in + (q, e', new_location m2 r). + + Definition SFdiv x y := + match x, y with + | S754_nan, _ | _, S754_nan => S754_nan + | S754_infinity sx, S754_infinity sy => S754_nan + | S754_infinity sx, S754_finite sy _ _ => S754_infinity (xorb sx sy) + | S754_finite sx _ _, S754_infinity sy => S754_zero (xorb sx sy) + | S754_infinity sx, S754_zero sy => S754_infinity (xorb sx sy) + | S754_zero sx, S754_infinity sy => S754_zero (xorb sx sy) + | S754_finite sx _ _, S754_zero sy => S754_infinity (xorb sx sy) + | S754_zero sx, S754_finite sy _ _ => S754_zero (xorb sx sy) + | S754_zero sx, S754_zero sy => S754_nan + | S754_finite sx mx ex, S754_finite sy my ey => + let '(mz, ez, lz) := SFdiv_core_binary (Zpos mx) ex (Zpos my) ey in + binary_round_aux (xorb sx sy) mz ez lz + end. + + Definition SFsqrt_core_binary m e := + let d := Zdigits2 m in + let e' := Z.min (fexp (Z.div2 (d + e + 1))) (Z.div2 e) in + let s := (e - 2 * e')%Z in + let m' := + match s with + | Zpos p => Z.shiftl m s + | Z0 => m + | Zneg _ => Z0 + end in + let (q, r) := Z.sqrtrem m' in + let l := + if Zeq_bool r 0 then loc_Exact + else loc_Inexact (if Zle_bool r q then Lt else Gt) in + (q, e', l). + + Definition SFsqrt x := + match x with + | S754_nan => S754_nan + | S754_infinity false => x + | S754_infinity true => S754_nan + | S754_finite true _ _ => S754_nan + | S754_zero _ => x + | S754_finite sx mx ex => + let '(mz, ez, lz) := SFsqrt_core_binary (Zpos mx) ex in + binary_round_aux false mz ez lz + end. + + Definition SFnormfr_mantissa f := + match f with + | S754_finite _ mx ex => + if Z.eqb ex (-prec) then Npos mx else 0%N + | _ => 0%N + end. + + Definition SFldexp f e := + match f with + | S754_finite sx mx ex => binary_round sx mx (ex+e) + | _ => f + end. + + Definition SFfrexp f := + match f with + | S754_finite sx mx ex => + if (Z.to_pos prec <=? digits2_pos mx)%positive then + (S754_finite sx mx (-prec), (ex+prec)%Z) + else + let d := (prec - Z.pos (digits2_pos mx))%Z in + (S754_finite sx (shift_pos (Z.to_pos d) mx) (-prec), (ex+prec-d)%Z) + | _ => (f, (-2*emax-prec)%Z) + end. + + Definition SFone := binary_round false 1 0. + + Definition SFulp x := SFldexp SFone (fexp (snd (SFfrexp x))). + + Definition SFpred_pos x := + match x with + | S754_finite _ mx _ => + let d := + if (mx~0 =? shift_pos (Z.to_pos prec) 1)%positive then + SFldexp SFone (fexp (snd (SFfrexp x) - 1)) + else + SFulp x in + SFsub x d + | _ => x + end. + + Definition SFmax_float := + S754_finite false (shift_pos (Z.to_pos prec) 1 - 1) (emax - prec). + + Definition SFsucc x := + match x with + | S754_zero _ => SFldexp SFone emin + | S754_infinity false => x + | S754_infinity true => SFopp SFmax_float + | S754_nan => x + | S754_finite false _ _ => SFadd x (SFulp x) + | S754_finite true _ _ => SFopp (SFpred_pos (SFopp x)) + end. + + Definition SFpred f := SFopp (SFsucc (SFopp f)). +End FloatOps. -- cgit