diff options
Diffstat (limited to 'flocq/Appli/Fappli_IEEE.v')
-rw-r--r-- | flocq/Appli/Fappli_IEEE.v | 227 |
1 files changed, 189 insertions, 38 deletions
diff --git a/flocq/Appli/Fappli_IEEE.v b/flocq/Appli/Fappli_IEEE.v index b44d711c..9b5826c1 100644 --- a/flocq/Appli/Fappli_IEEE.v +++ b/flocq/Appli/Fappli_IEEE.v @@ -46,6 +46,8 @@ End AnyRadix. Section Binary. +Implicit Arguments exist [[A] [P]]. + (** prec is the number of bits of the mantissa including the implicit one emax is the exponent of the infinities Typically p=24 and emax = 128 in single precision *) @@ -59,7 +61,7 @@ Instance fexp_correct : Valid_exp fexp := FLT_exp_valid emin prec. Instance fexp_monotone : Monotone_exp fexp := FLT_exp_monotone emin prec. Definition canonic_mantissa m e := - Zeq_bool (fexp (Z_of_nat (S (digits2_Pnat m)) + e)) e. + Zeq_bool (fexp (Zpos (digits2_pos m) + e)) e. Definition bounded m e := andb (canonic_mantissa m e) (Zle_bool e (emax - prec)). @@ -67,14 +69,15 @@ Definition bounded m e := Definition valid_binary x := match x with | F754_finite _ m e => bounded m e - | F754_nan _ pl => (Z_of_nat' (S (digits2_Pnat pl)) <? prec)%Z + | F754_nan _ pl => (Zpos (digits2_pos pl) <? prec)%Z | _ => true end. (** Basic type used for representing binary FP numbers. - Note that there is exactly one such object per FP datum. *) + Note that there is exactly one such object per FP datum. + NaNs do not have any payload. They cannot be distinguished. *) -Definition nan_pl := {pl | (Z_of_nat' (S (digits2_Pnat pl)) <? prec)%Z = true}. +Definition nan_pl := {pl | (Zpos (digits2_pos pl) <? prec)%Z = true}. Inductive binary_float := | B754_zero : bool -> binary_float @@ -88,7 +91,7 @@ Definition FF2B x := | F754_finite s m e => B754_finite s m e | F754_infinity s => fun _ => B754_infinity s | F754_zero s => fun _ => B754_zero s - | F754_nan b pl => fun H => B754_nan b (exist _ pl H) + | F754_nan b pl => fun H => B754_nan b (exist pl H) end. Definition B2FF x := @@ -99,8 +102,6 @@ Definition B2FF x := | B754_nan b (exist pl _) => F754_nan b pl end. -Definition radix2 := Build_radix 2 (refl_equal true). - Definition B2R f := match f with | B754_finite s m e _ => F2R (Float radix2 (cond_Zopp s (Zpos m)) e) @@ -183,7 +184,7 @@ pattern ex at 2 ; rewrite <- Hx. apply (f_equal fexp). rewrite ln_beta_F2R_Zdigits. rewrite <- Zdigits_abs. -rewrite Z_of_nat_S_digits2_Pnat. +rewrite Zpos_digits2_pos. now case sx. now case sx. Qed. @@ -419,6 +420,159 @@ simpl. now case opp_nan. Qed. +(** Absolute value *) + +Definition Babs abs_nan (x : binary_float) : binary_float := + match x with + | B754_nan sx plx => + let '(sres, plres) := abs_nan sx plx in B754_nan sres plres + | B754_infinity sx => B754_infinity false + | B754_finite sx mx ex Hx => B754_finite false mx ex Hx + | B754_zero sx => B754_zero false + end. + +Theorem B2R_Babs : + forall abs_nan x, + B2R (Babs abs_nan x) = Rabs (B2R x). +Proof. + intros abs_nan [sx|sx|sx plx|sx mx ex Hx]; apply sym_eq ; try apply Rabs_R0. + simpl. destruct abs_nan. simpl. apply Rabs_R0. + simpl. rewrite <- F2R_abs. now destruct sx. +Qed. + +Theorem is_finite_Babs : + forall abs_nan x, + is_finite (Babs abs_nan x) = is_finite x. +Proof. + intros abs_nan [| | |] ; try easy. + intros s pl. + simpl. + now case abs_nan. +Qed. + +Theorem Bsign_Babs : + forall abs_nan x, + is_nan x = false -> + Bsign (Babs abs_nan x) = false. +Proof. + now intros abs_nan [| | |]. +Qed. + +Theorem Babs_idempotent : + forall abs_nan (x: binary_float), + is_nan x = false -> + Babs abs_nan (Babs abs_nan x) = Babs abs_nan x. +Proof. + now intros abs_nan [sx|sx|sx plx|sx mx ex Hx]. +Qed. + +Theorem Babs_Bopp : + forall abs_nan opp_nan x, + is_nan x = false -> + Babs abs_nan (Bopp opp_nan x) = Babs abs_nan x. +Proof. + now intros abs_nan opp_nan [| | |]. +Qed. + +(** Comparison + +[Some c] means ordered as per [c]; [None] means unordered. *) + +Definition Bcompare (f1 f2 : binary_float) : option comparison := + match f1, f2 with + | 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 + | _, B754_infinity true => Some Gt + | _, B754_infinity false => Some Lt + | B754_finite true _ _ _, B754_zero _ => Some Lt + | B754_finite false _ _ _, B754_zero _ => Some Gt + | B754_zero _, B754_finite true _ _ _ => Some Gt + | B754_zero _, B754_finite false _ _ _ => Some Lt + | B754_zero _, B754_zero _ => Some Eq + | B754_finite s1 m1 e1 _, B754_finite s2 m2 e2 _ => + match s1, s2 with + | true, false => Some Lt + | false, true => Some Gt + | false, false => + match Zcompare e1 e2 with + | Lt => Some Lt + | Gt => Some Gt + | Eq => Some (Pcompare m1 m2 Eq) + end + | true, true => + match Zcompare e1 e2 with + | Lt => Some Gt + | Gt => Some Lt + | Eq => Some (CompOpp (Pcompare m1 m2 Eq)) + end + end + end. + +Theorem Bcompare_correct : + forall f1 f2, + is_finite f1 = true -> is_finite f2 = true -> + Bcompare f1 f2 = Some (Rcompare (B2R f1) (B2R f2)). +Proof. + Ltac apply_Rcompare := + match goal with + | [ |- Some Lt = Some (Rcompare _ _) ] => f_equal; symmetry; apply Rcompare_Lt + | [ |- Some Eq = Some (Rcompare _ _) ] => f_equal; symmetry; apply Rcompare_Eq + | [ |- Some Gt = Some (Rcompare _ _) ] => f_equal; symmetry; apply Rcompare_Gt + end. + unfold Bcompare; intros. + destruct f1, f2 ; try easy. + now rewrite Rcompare_Eq. + destruct b0 ; apply_Rcompare. + now apply F2R_lt_0_compat. + now apply F2R_gt_0_compat. + destruct b ; apply_Rcompare. + now apply F2R_lt_0_compat. + now apply F2R_gt_0_compat. + simpl. + clear H H0. + apply andb_prop in e0; destruct e0; apply (canonic_canonic_mantissa false) in H. + apply andb_prop in e2; destruct e2; apply (canonic_canonic_mantissa false) in H1. + pose proof (Zcompare_spec e e1); unfold canonic, Fexp in H1, H. + assert (forall m1 m2 e1 e2, + let x := (Z2R (Zpos m1) * bpow radix2 e1)%R in + let y := (Z2R (Zpos m2) * bpow radix2 e2)%R in + (canonic_exp radix2 fexp x < canonic_exp radix2 fexp y)%Z -> (x < y)%R). + { + intros; apply Rnot_le_lt; intro; apply (ln_beta_le radix2) in H5. + apply Zlt_not_le with (1 := H4). + now apply fexp_monotone. + now apply (F2R_gt_0_compat _ (Float radix2 (Zpos m2) e2)). + } + assert (forall m1 m2 e1 e2, (Z2R (- Zpos m1) * bpow radix2 e1 < Z2R (Zpos m2) * bpow radix2 e2)%R). + { + intros; apply (Rlt_trans _ 0%R). + now apply (F2R_lt_0_compat _ (Float radix2 (Zneg m1) e0)). + now apply (F2R_gt_0_compat _ (Float radix2 (Zpos m2) e2)). + } + unfold F2R, Fnum, Fexp. + 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; + 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. + +Theorem Bcompare_swap : + forall x y, + Bcompare y x = match Bcompare x y with Some c => Some (CompOpp c) | None => None end. +Proof. + intros. + destruct x as [ ? | [] | ? ? | [] mx ex Bx ]; + destruct y as [ ? | [] | ? ? | [] my ey By ]; simpl; try easy. +- rewrite <- (Zcompare_antisym ex ey). destruct (ex ?= ey)%Z; try easy. + now rewrite (Pcompare_antisym mx my). +- rewrite <- (Zcompare_antisym ex ey). destruct (ex ?= ey)%Z; try easy. + now rewrite Pcompare_antisym. +Qed. + Theorem bounded_lt_emax : forall mx ex, bounded mx ex = true -> @@ -441,8 +595,7 @@ now apply F2R_gt_0_compat. apply bpow_le. rewrite H. 2: discriminate. revert H1. clear -H2. -rewrite Z_of_nat_S_digits2_Pnat. -change Fcalc_digits.radix2 with radix2. +rewrite Zpos_digits2_pos. unfold fexp, FLT_exp. generalize (Zdigits radix2 (Zpos mx)). intros ; zify ; subst. @@ -471,8 +624,7 @@ split. unfold canonic_mantissa. unfold canonic, Fexp in Cx. rewrite Cx at 2. -rewrite Z_of_nat_S_digits2_Pnat. -change Fcalc_digits.radix2 with radix2. +rewrite Zpos_digits2_pos. unfold canonic_exp. rewrite ln_beta_F2R_Zdigits. 2: discriminate. now apply -> Zeq_is_eq_bool. @@ -600,7 +752,7 @@ induction (nat_of_P n). simpl. rewrite Zplus_0_r. now destruct l as [|[| |]]. -simpl iter_nat. +simpl nat_rect. rewrite inj_S. unfold Zsucc. rewrite Zplus_assoc. @@ -617,18 +769,6 @@ intros (m, r, s) Hm. now destruct m as [|[m|m|]|m] ; try (now elim Hm) ; destruct r as [|] ; destruct s as [|]. Qed. -Definition Zdigits2 m := - match m with Z0 => m | Zpos p => Z_of_nat (S (digits2_Pnat p)) | Zneg p => Z_of_nat (S (digits2_Pnat p)) end. - -Theorem Zdigits2_Zdigits : - forall m, - Zdigits2 m = Zdigits radix2 m. -Proof. -unfold Zdigits2. -intros [|m|m] ; try apply Z_of_nat_S_digits2_Pnat. -easy. -Qed. - Definition shr_fexp m e l := shr (shr_record_of_loc m l) e (fexp (Zdigits2 m + e) - e). @@ -829,7 +969,7 @@ apply andb_true_intro. split. unfold canonic_mantissa. apply Zeq_bool_true. -rewrite Z_of_nat_S_digits2_Pnat. +rewrite Zpos_digits2_pos. rewrite <- ln_beta_F2R_Zdigits. apply sym_eq. now rewrite H3 in H4. @@ -849,8 +989,7 @@ unfold valid_binary, bounded. rewrite Zle_bool_refl. rewrite Bool.andb_true_r. apply Zeq_bool_true. -rewrite Z_of_nat_S_digits2_Pnat. -change Fcalc_digits.radix2 with radix2. +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). @@ -942,7 +1081,7 @@ assert (forall m e, bounded m e = true -> fexp (Zdigits radix2 (Zpos m) + e) = e clear. intros m e Hb. destruct (andb_prop _ _ Hb) as (H,_). apply Zeq_bool_eq. -now rewrite <- Z_of_nat_S_digits2_Pnat. +now rewrite <- Zpos_digits2_pos. generalize (H _ _ Hx) (H _ _ Hy). clear x y sx sy Hx Hy H. unfold fexp, FLT_exp. @@ -1098,7 +1237,7 @@ apply refl_equal. Qed. Definition shl_align_fexp mx ex := - shl_align mx ex (fexp (Z_of_nat (S (digits2_Pnat mx)) + ex)). + shl_align mx ex (fexp (Zpos (digits2_pos mx) + ex)). Theorem shl_align_fexp_correct : forall mx ex, @@ -1108,8 +1247,8 @@ Theorem shl_align_fexp_correct : Proof. intros mx ex. unfold shl_align_fexp. -generalize (shl_align_correct mx ex (fexp (Z_of_nat (S (digits2_Pnat mx)) + ex))). -rewrite Z_of_nat_S_digits2_Pnat. +generalize (shl_align_correct mx ex (fexp (Zpos (digits2_pos mx) + ex))). +rewrite Zpos_digits2_pos. case shl_align. intros mx' ex' (H1, H2). split. @@ -1438,10 +1577,10 @@ Definition Fdiv_core_binary m1 e1 m2 e2 := let e := (e1 - e2)%Z in let (m, e') := match (d2 + prec - d1)%Z with - | Zpos p => (m1 * Zpower_pos 2 p, e + Zneg p)%Z + | Zpos p => (Z.shiftl m1 (Zpos p), e + Zneg p)%Z | _ => (m1, e) end in - let '(q, r) := Zdiv_eucl m m2 in + let '(q, r) := Zfast_div_eucl m m2 in (q, e', new_location m2 r loc_Exact). Lemma Bdiv_correct_aux : @@ -1463,7 +1602,6 @@ Lemma Bdiv_correct_aux : Proof. intros m sx mx ex sy my ey. replace (Fdiv_core_binary (Zpos mx) ex (Zpos my) ey) with (Fdiv_core radix2 prec (Zpos mx) ex (Zpos my) ey). -2: now unfold Fdiv_core_binary ; rewrite 2!Zdigits2_Zdigits. refine (_ (Fdiv_core_correct radix2 prec (Zpos mx) ex (Zpos my) ey _ _ _)) ; try easy. destruct (Fdiv_core radix2 prec (Zpos mx) ex (Zpos my) ey) as ((mz, ez), lz). intros (Pz, Bz). @@ -1540,6 +1678,15 @@ now apply F2R_ge_0_compat. apply Rlt_le. apply Rinv_0_lt_compat. now apply F2R_gt_0_compat. +(* *) +unfold Fdiv_core_binary, Fdiv_core. +rewrite 2!Zdigits2_Zdigits. +change 2%Z with (radix_val radix2). +destruct (Zdigits radix2 (Z.pos my) + prec - Zdigits radix2 (Z.pos mx))%Z as [|p|p]. +now rewrite Zfast_div_eucl_correct. +rewrite Z.shiftl_mul_pow2 by easy. +now rewrite Zfast_div_eucl_correct. +now rewrite Zfast_div_eucl_correct. Qed. Definition Bdiv div_nan m x y := @@ -1597,10 +1744,10 @@ Definition Fsqrt_core_binary m e := let (s', e'') := if Zeven e' then (s, e') else (s + 1, e' - 1)%Z in let m' := match s' with - | Zpos p => (m * Zpower_pos 2 p)%Z + | Zpos p => Z.shiftl m (Zpos p) | _ => m end in - let (q, r) := Zsqrt m' 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 @@ -1621,7 +1768,6 @@ Lemma Bsqrt_correct_aux : Proof with auto with typeclass_instances. intros m mx ex Hx. replace (Fsqrt_core_binary (Zpos mx) ex) with (Fsqrt_core radix2 prec (Zpos mx) ex). -2: now unfold Fsqrt_core_binary ; rewrite Zdigits2_Zdigits. simpl. refine (_ (Fsqrt_core_correct radix2 prec (Zpos mx) ex _)) ; try easy. destruct (Fsqrt_core radix2 prec (Zpos mx) ex) as ((mz, ez), lz). @@ -1716,6 +1862,11 @@ apply Rle_trans with R0. apply F2R_le_0_compat. now case mz. apply sqrt_ge_0. +(* *) +unfold Fsqrt_core, Fsqrt_core_binary. +rewrite Zdigits2_Zdigits. +destruct (if Zeven _ then _ else _) as [[|s'|s'] e''] ; try easy. +now rewrite Z.shiftl_mul_pow2. Qed. Definition Bsqrt sqrt_nan m x := |