Commit 9e2ea1b4 authored by Guillaume Melquiond's avatar Guillaume Melquiond

Make Fsqrt_core generic with respect to format.

parent 3424cf9b
......@@ -82,7 +82,7 @@ Qed.
Definition sqrt (x : float beta) :=
let (m, e) := x in
if Zlt_bool 0 m then
let '(m', e', l) := truncate beta fexp (Fsqrt_core beta prec m e) in
let '(m', e', l) := truncate beta fexp (Fsqrt_core beta fexp m e) in
Float beta (choice false m' l) e'
else Float beta 0 0.
......@@ -93,16 +93,14 @@ Proof.
intros [m e].
unfold sqrt.
case Zlt_bool_spec ; intros Hm.
generalize (Fsqrt_core_correct beta prec m e Hm).
generalize (Fsqrt_core_correct beta fexp m e Hm).
destruct Fsqrt_core as [[m' e'] l].
intros [Hs1 Hs2].
rewrite (round_trunc_sign_any_correct beta fexp rnd choice rnd_choice _ m' e' l).
rewrite (round_trunc_sign_any_correct' beta fexp rnd choice rnd_choice _ m' e' l).
destruct truncate as [[m'' e''] l'].
now rewrite Rlt_bool_false by apply sqrt_ge_0.
now rewrite Rabs_pos_eq by apply sqrt_ge_0.
left.
apply Zle_trans with (2 := fexp_prec _).
clear -Hs1 ; omega.
now left.
destruct (Req_dec (F2R (Float beta m e)) 0) as [Hz|Hz].
rewrite Hz, sqrt_0, F2R_0.
now apply round_0.
......
......@@ -19,13 +19,15 @@ COPYING file for more details.
(** * Helper functions and theorems for computing the rounded square root of a floating-point number. *)
Require Import Raux Defs Digits Float_prop Bracket.
Require Import Raux Defs Digits Generic_fmt Float_prop Bracket.
Section Fcalc_sqrt.
Variable beta : radix.
Notation bpow e := (bpow beta e).
Variable fexp : Z -> Z.
(** Computes a mantissa of precision p, the corresponding exponent,
and the position with respect to the real square root of the
input floating-point number.
......@@ -43,126 +45,92 @@ The algorithm performs the following steps:
Complexity is fine as long as p1 <= 2p-1.
*)
Definition Fsqrt_core prec m e :=
Definition Fsqrt_core m e :=
let d := Zdigits beta m in
let s := Zmax (2 * prec - d) 0 in
let e' := (e - s)%Z in
let (s', e'') := if Z.even e' then (s, e') else (s + 1, e' - 1)%Z in
let m' :=
match s' with
| Zpos p => (m * Zpower_pos beta p)%Z
| _ => m
end in
let e' := Zmin (fexp (Z.div2 (d + e + 1))) (Z.div2 e) in
let m' := (m * Zpower beta (e - 2 * e'))%Z 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, Zdiv2 e'', l).
(q, e', l).
Theorem Fsqrt_core_correct :
forall prec m e,
forall m e,
(0 < m)%Z ->
let '(m', e', l) := Fsqrt_core prec m e in
(prec <= Zdigits beta m')%Z /\
let '(m', e', l) := Fsqrt_core m e in
(e' <= cexp beta fexp (sqrt (F2R (Float beta m e))))%Z /\
inbetween_float beta m' e' (sqrt (F2R (Float beta m e))) l.
Proof.
intros prec m e Hm.
intros m e Hm.
unfold Fsqrt_core.
set (d := Zdigits beta m).
set (s := Zmax (2 * prec - d) 0).
(* . exponent *)
case_eq (if Z.even (e - s) then (s, (e - s)%Z) else ((s + 1)%Z, (e - s - 1)%Z)).
intros s' e' Hse.
assert (He: (Z.even e' = true /\ 0 <= s' /\ 2 * prec - d <= s' /\ s' + e' = e)%Z).
revert Hse.
case_eq (Z.even (e - s)) ; intros He Hse ; inversion Hse.
repeat split.
exact He.
unfold s.
apply Zle_max_r.
apply Zle_max_l.
ring.
assert (H: (Zmax (2 * prec - d) 0 <= s + 1)%Z).
fold s.
apply Zle_succ.
repeat split.
unfold Zminus at 1.
now rewrite Z.even_add, He.
apply Zle_trans with (2 := H).
apply Zle_max_r.
apply Zle_trans with (2 := H).
apply Zle_max_l.
ring.
clear -Hm He.
destruct He as (He1, (He2, (He3, He4))).
(* . shift *)
set (m' := match s' with
| Z0 => m
| Zpos p => (m * Zpower_pos beta p)%Z
| Zneg _ => m
end).
assert (Hs: F2R (Float beta m' e') = F2R (Float beta m e) /\ (2 * prec <= Zdigits beta m')%Z /\ (0 < m')%Z).
rewrite <- He4.
unfold m'.
destruct s' as [|p|p].
repeat split ; try easy.
fold d.
omega.
fold (Zpower beta (Zpos p)).
split.
replace (Zpos p) with (Zpos p + e' - e')%Z by ring.
rewrite <- F2R_change_exp.
apply (f_equal (fun v => F2R (Float beta m v))).
ring.
assert (0 < Zpos p)%Z by easy.
omega.
split.
rewrite Zdigits_mult_Zpower.
fold d.
omega.
apply sym_not_eq.
now apply Zlt_not_eq.
easy.
apply Zmult_lt_0_compat.
exact Hm.
now apply Zpower_gt_0.
now elim He2.
clearbody m'.
destruct Hs as (Hs1, (Hs2, Hs3)).
generalize (Z.sqrtrem_spec m' (Zlt_le_weak _ _ Hs3)).
destruct (Z.sqrtrem m') as (q, r).
intros (Hq, Hr).
rewrite <- Hs1. clear Hs1.
set (e'' := Z.div2 (d + e + 1)).
set (e' := Zmin (fexp e'') (Z.div2 e)).
set (s := (e - 2 * e')%Z).
set (m' := (m * Zpower beta s)%Z).
assert (bpow (e'' - 1) <= sqrt (F2R (Float beta m e)) < bpow e'')%R as Hs.
{ unfold e'', d. clear -Hm.
rewrite <- (mag_F2R_Zdigits beta m e (Zgt_not_eq _ _ Hm)).
destruct mag as [e' Hs'].
simpl.
assert (bpow (e' - 1) <= F2R (Float beta m e) < bpow e')%R as Hs.
rewrite Rabs_pos_eq in Hs'.
now apply Hs', F2R_neq_0, Zgt_not_eq.
now apply F2R_ge_0, Zlt_le_weak.
clear Hs'.
split.
- rewrite <- (sqrt_Rsqr (bpow _)) by now apply bpow_ge_0.
apply sqrt_le_1_alt.
unfold Rsqr.
rewrite <- bpow_plus.
apply Rle_trans with (2 := proj1 Hs).
apply bpow_le.
replace (e' - 1)%Z with (e' + 1 - 2)%Z by ring.
rewrite (Zdiv2_odd_eqn (e' + 1)) at 3.
case Z.odd ; clear ; omega.
- rewrite <- (sqrt_Rsqr (bpow _)) by now apply bpow_ge_0.
apply sqrt_lt_1_alt.
split.
now apply F2R_ge_0, Zlt_le_weak.
unfold Rsqr.
rewrite <- bpow_plus.
apply Rlt_le_trans with (1 := proj2 Hs).
apply bpow_le.
replace e' with (e' + 1 - 1)%Z at 1 by ring.
rewrite (Zdiv2_odd_eqn (e' + 1)) at 1.
case Z.odd ; clear ; omega. }
assert (F2R (Float beta m e) = F2R (Float beta m' (2 * e'))) as Hf.
apply F2R_change_exp.
unfold e'.
rewrite <- Z.mul_min_distr_nonneg_l by easy.
apply Zle_trans with (1 := Z.le_min_r _ _).
rewrite (Zdiv2_odd_eqn e) at 2.
case Z.odd ; clear ; omega.
assert (0 <= m')%Z as Hm'.
apply Z.mul_nonneg_nonneg.
now apply Zlt_le_weak.
apply Zpower_ge_0.
generalize (Z.sqrtrem_spec m' Hm').
destruct Z.sqrtrem as [q r].
intros [Hq Hr].
split.
(* . mantissa width *)
apply Zmult_le_reg_r with 2%Z.
easy.
rewrite Zmult_comm.
apply Zle_trans with (1 := Hs2).
rewrite Hq.
apply Zle_trans with (Zdigits beta (q + q + q * q)).
apply Zdigits_le.
rewrite <- Hq.
now apply Zlt_le_weak.
omega.
replace (Zdigits beta q * 2)%Z with (Zdigits beta q + Zdigits beta q)%Z by ring.
apply Zdigits_mult_strong.
omega.
omega.
(* . round *)
unfold inbetween_float, F2R. simpl.
unfold cexp.
rewrite (mag_unique _ _ e'').
apply Z.le_min_l.
rewrite Rabs_pos_eq.
exact Hs.
apply sqrt_ge_0.
rewrite Hf.
unfold inbetween_float, F2R. simpl Fnum.
rewrite sqrt_mult.
2: now apply IZR_le, Zlt_le_weak.
2: now apply IZR_le.
2: apply Rlt_le ; apply bpow_gt_0.
destruct (Zeven_ex e') as (e2, Hev).
rewrite He1, Zplus_0_r in Hev. clear He1.
rewrite Hev.
replace (Zdiv2 (2 * e2)) with e2 by now case e2.
replace (2 * e2)%Z with (e2 + e2)%Z by ring.
replace (2 * e')%Z with (e' + e')%Z by ring.
simpl Fexp.
rewrite bpow_plus.
fold (Rsqr (bpow e2)).
rewrite sqrt_Rsqr.
2: apply Rlt_le ; apply bpow_gt_0.
fold (Rsqr (bpow e')).
rewrite sqrt_Rsqr by apply bpow_ge_0.
apply inbetween_mult_compat.
apply bpow_gt_0.
rewrite Hq.
......@@ -189,15 +157,13 @@ apply sqrt_lt_1.
rewrite mult_IZR.
apply Rle_0_sqr.
rewrite <- Hq.
apply IZR_le.
now apply Zlt_le_weak.
now apply IZR_le.
apply IZR_lt.
omega.
apply Rlt_le_trans with (sqrt (IZR ((q + 1) * (q + 1)))).
apply sqrt_lt_1.
rewrite <- Hq.
apply IZR_le.
now apply Zlt_le_weak.
now apply IZR_le.
rewrite mult_IZR.
apply Rle_0_sqr.
apply IZR_lt.
......@@ -227,8 +193,7 @@ omega.
change (4 * (q * q + r) > 4 * (q * q) + 4 * q + 1)%Z.
omega.
rewrite <- Hq.
apply IZR_le.
now apply Zlt_le_weak.
now apply IZR_le.
rewrite <- plus_IZR.
apply IZR_le.
clear -Hr ; omega.
......
......@@ -615,6 +615,34 @@ clear -H H2. clearbody emin.
omega.
Qed.
Theorem bounded_ge_emin :
forall mx ex,
bounded mx ex = true ->
(bpow radix2 emin <= F2R (Float radix2 (Zpos mx) ex))%R.
Proof.
intros mx ex Hx.
destruct (andb_prop _ _ Hx) as [H1 _].
apply Zeq_bool_eq in H1.
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.
assert (H0 : Zpos mx <> 0%Z) by easy.
rewrite Rabs_pos_eq in Ex by now apply F2R_ge_0.
refine (Rle_trans _ _ _ _ (proj1 (Ex _))).
2: now apply F2R_neq_0.
apply bpow_le.
rewrite H by easy.
revert H1.
rewrite Zpos_digits2_pos.
generalize (Zdigits radix2 (Zpos mx)) (Zdigits_gt_0 radix2 (Zpos mx) H0).
unfold fexp, FLT_exp.
clear -prec_gt_0_.
unfold Prec_gt_0 in prec_gt_0_.
clearbody emin.
intros ; zify ; omega.
Qed.
Theorem abs_B2R_lt_emax :
forall x,
(Rabs (B2R x) < bpow radix2 emax)%R.
......@@ -1921,19 +1949,19 @@ Qed.
Definition Fsqrt_core_binary m e :=
let d := Zdigits2 m in
let s := Zmax (2 * prec - d) 0 in
let e' := (e - s)%Z in
let (s', e'') := if Z.even e' then (s, e') else (s + 1, e' - 1)%Z in
let e' := Zmin (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 (Zpos p)
| _ => 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, Zdiv2 e'', l).
(q, e', l).
Lemma Bsqrt_correct_aux :
forall m mx ex (Hx : bounded mx ex = true),
......@@ -1949,21 +1977,54 @@ Lemma Bsqrt_correct_aux :
is_finite_FF z = true /\ sign_FF z = false.
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).
assert (Fsqrt_core_binary (Zpos mx) ex = Fsqrt_core radix2 fexp (Zpos mx) ex) as ->.
unfold Fsqrt_core, Fsqrt_core_binary.
rewrite Zdigits2_Zdigits.
set (e' := Zmin (fexp (Z.div2 (Zdigits radix2 (Zpos mx) + ex + 1))) (Z.div2 ex)).
destruct (ex - 2 * e')%Z as [|s|s].
now rewrite Zmult_1_r.
now rewrite Z.shiftl_mul_pow2.
easy.
simpl.
refine (_ (Fsqrt_core_correct radix2 prec (Zpos mx) ex _)) ; try easy.
destruct (Fsqrt_core radix2 prec (Zpos mx) ex) as ((mz, ez), lz).
refine (_ (Fsqrt_core_correct radix2 fexp (Zpos mx) ex _)) ; try easy.
destruct (Fsqrt_core radix2 fexp (Zpos mx) ex) as ((mz, ez), lz).
intros (Pz, Bz).
destruct mz as [|mz|mz].
(* . mz = 0 *)
elim (Zlt_irrefl prec).
now apply Zle_lt_trans with Z0.
(* . mz > 0 *)
refine (_ (binary_round_aux_correct m (sqrt (F2R (Float radix2 (Zpos mx) ex))) mz ez lz _ _)).
rewrite Rlt_bool_false. 2: apply sqrt_ge_0.
- apply inbetween_float_bounds in Bz.
elim (Zlt_irrefl ez).
apply Zle_lt_trans with (1 := Pz).
apply lt_bpow with radix2.
rewrite <- (F2R_bpow radix2 ez).
apply Rle_lt_trans with (2 := proj2 Bz).
clear -Hx prec_gt_0_ Hmax.
apply bounded_ge_emin in Hx.
unfold cexp.
destruct mag as [e He].
simpl.
refine (_ (He _)).
2: now apply Rgt_not_eq, sqrt_lt_R0, F2R_gt_0.
clear He.
rewrite Rabs_pos_eq by apply sqrt_ge_0.
intros He.
apply Rle_trans with (2 := proj1 He).
apply bpow_le.
apply sqrt_le_1_alt in Hx.
apply (fun H => Rle_lt_trans _ _ _ H (proj2 He)) in Hx.
rewrite <- (sqrt_Rsqr (bpow radix2 e)) in Hx by apply bpow_ge_0.
apply sqrt_lt_0_alt in Hx.
unfold Rsqr in Hx.
rewrite <- bpow_plus in Hx.
apply lt_bpow in Hx.
unfold fexp, FLT_exp.
unfold Prec_gt_0 in prec_gt_0_.
revert Hx.
unfold emin.
intros ; zify ; omega.
- refine (_ (binary_round_aux_correct' m (sqrt (F2R (Float radix2 (Zpos mx) ex))) mz ez lz _ _)).
rewrite Rlt_bool_false.
2: apply sqrt_ge_0.
rewrite Rlt_bool_true.
easy.
(* .. *)
rewrite Rabs_pos_eq.
refine (_ (relative_error_FLT_ex radix2 emin prec (prec_gt_0 prec) (round_mode m) (sqrt (F2R (Float radix2 (Zpos mx) ex))) _)).
fold fexp.
......@@ -2038,17 +2099,11 @@ unfold fexp, FLT_exp.
clear.
intros ; zify ; subst.
omega.
(* . mz < 0 *)
elim Rlt_not_le with (1 := proj2 (inbetween_float_bounds _ _ _ _ _ Bz)).
- elim Rlt_not_le with (1 := proj2 (inbetween_float_bounds _ _ _ _ _ Bz)).
apply Rle_trans with R0.
apply F2R_le_0.
now case mz.
apply sqrt_ge_0.
(* *)
unfold Fsqrt_core, Fsqrt_core_binary.
rewrite Zdigits2_Zdigits.
destruct (if Z.even _ then _ else _) as [[|s'|s'] e''] ; try easy.
now rewrite Z.shiftl_mul_pow2.
Qed.
Definition Bsqrt sqrt_nan m x :=
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment