Commit 63602c45 authored by Guillaume Melquiond's avatar Guillaume Melquiond

Made proofs about Fsqrt_core a bit more reusable.

parent cffdf80d
......@@ -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 fexp m e) in
let '(m', e', l) := truncate beta fexp (Fsqrt beta fexp m e) in
Float beta (choice false m' l) e'
else Float beta 0 0.
......@@ -93,8 +93,8 @@ Proof.
intros [m e].
unfold sqrt.
case Zlt_bool_spec ; intros Hm.
generalize (Fsqrt_core_correct beta fexp m e Hm).
destruct Fsqrt_core as [[m' e'] l].
generalize (Fsqrt_correct beta fexp m e Hm).
destruct Fsqrt as [[m' e'] l].
intros [Hs1 Hs2].
rewrite (round_trunc_sign_any_correct' beta fexp rnd choice rnd_choice _ m' e' l).
destruct truncate as [[m'' e''] l'].
......
......@@ -45,92 +45,54 @@ The algorithm performs the following steps:
Complexity is fine as long as p1 <= 2p-1.
*)
Definition Fsqrt_core m e :=
let d := Zdigits beta m 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
Lemma mag_sqrt_F2R :
forall m1 e1,
(0 < m1)%Z ->
mag beta (sqrt (F2R (Float beta m1 e1))) = Zdiv2 (Zdigits beta m1 + e1 + 1) :> Z.
Proof.
intros m1 e1 Hm1.
rewrite <- (mag_F2R_Zdigits beta m1 e1) by now apply Zgt_not_eq.
apply mag_sqrt.
now apply F2R_gt_0.
Qed.
Definition Fsqrt_core m1 e1 e :=
let d1 := Zdigits beta m1 in
let m1' := (m1 * Zpower beta (e1 - 2 * e))%Z in
let (q, r) := Z.sqrtrem m1' 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).
(q, l).
Theorem Fsqrt_core_correct :
forall m e,
(0 < 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.
forall m1 e1 e,
(0 < m1)%Z ->
(2 * e <= e1)%Z ->
let '(m, l) := Fsqrt_core m1 e1 e in
inbetween_float beta m e (sqrt (F2R (Float beta m1 e1))) l.
Proof.
intros m e Hm.
intros m1 e1 e Hm1 He.
unfold Fsqrt_core.
set (d := Zdigits beta m).
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.
set (m' := Zmult _ _).
assert (0 <= m')%Z as Hm'.
apply Z.mul_nonneg_nonneg.
{ apply Z.mul_nonneg_nonneg.
now apply Zlt_le_weak.
apply Zpower_ge_0.
apply Zpower_ge_0. }
assert (sqrt (F2R (Float beta m1 e1)) = sqrt (IZR m') * bpow e)%R as Hf.
{ rewrite <- (sqrt_Rsqr (bpow e)) by apply bpow_ge_0.
rewrite <- sqrt_mult.
unfold Rsqr, m'.
rewrite mult_IZR, IZR_Zpower by omega.
rewrite Rmult_assoc, <- 2!bpow_plus.
now replace (_ + _)%Z with e1 by ring.
now apply IZR_le.
apply Rle_0_sqr. }
generalize (Z.sqrtrem_spec m' Hm').
destruct Z.sqrtrem as [q r].
intros [Hq Hr].
split.
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.
2: apply Rlt_le ; apply bpow_gt_0.
replace (2 * e')%Z with (e' + e')%Z by ring.
simpl Fexp.
rewrite bpow_plus.
fold (Rsqr (bpow e')).
rewrite sqrt_Rsqr by apply bpow_ge_0.
apply inbetween_mult_compat.
apply bpow_gt_0.
rewrite Hq.
......@@ -202,4 +164,33 @@ now apply IZR_le.
apply sqrt_ge_0.
Qed.
Definition Fsqrt m1 e1 :=
let e' := (Zdigits beta m1 + e1 + 1)%Z in
let e := Zmin (fexp (Z.div2 e')) (Z.div2 e1) in
let (m, l) := Fsqrt_core m1 e1 e in
(m, e, l).
Theorem Fsqrt_correct :
forall m1 e1,
(0 < m1)%Z ->
let '(m, e, l) := Fsqrt m1 e1 in
(e <= cexp beta fexp (sqrt (F2R (Float beta m1 e1))))%Z /\
inbetween_float beta m e (sqrt (F2R (Float beta m1 e1))) l.
Proof.
intros m1 e1 Hm1.
unfold Fsqrt.
set (e := Zmin _ _).
assert (2 * e <= e1)%Z as He.
{ assert (e <= Zdiv2 e1)%Z by apply Zle_min_r.
rewrite (Zdiv2_odd_eqn e1).
destruct Z.odd ; omega. }
generalize (Fsqrt_core_correct m1 e1 e Hm1 He).
destruct Fsqrt_core as [m l].
apply conj.
apply Zle_trans with (1 := Zle_min_l _ _).
unfold cexp.
rewrite (mag_sqrt_F2R m1 e1 Hm1).
apply Zle_refl.
Qed.
End Fcalc_sqrt.
......@@ -1818,7 +1818,7 @@ Proof.
intros m sx mx ex sy my ey.
unfold Fdiv_core_binary.
rewrite 2!Zdigits2_Zdigits.
match goal with |- context [Zmin ?m1 ?m2] => set (e' := Zmin m1 m2) end.
set (e' := Zmin _ _).
generalize (Fdiv_core_correct radix2 (Zpos mx) ex (Zpos my) ey e' eq_refl eq_refl).
unfold Fdiv_core.
rewrite Zle_bool_true by apply Zle_min_r.
......@@ -1969,23 +1969,36 @@ Lemma Bsqrt_correct_aux :
is_finite_FF z = true /\ sign_FF z = false.
Proof with auto with typeclass_instances.
intros m mx ex Hx.
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].
unfold Fsqrt_core_binary.
rewrite Zdigits2_Zdigits.
set (e' := Zmin _ _).
assert (2 * e' <= ex)%Z as He.
{ assert (e' <= Zdiv2 ex)%Z by apply Zle_min_r.
rewrite (Zdiv2_odd_eqn ex).
destruct Z.odd ; omega. }
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).
assert (mx' = Zpos mx * Zpower radix2 (ex - 2 * e'))%Z as <-.
{ unfold mx'.
destruct (ex - 2 * e')%Z as [|p|p].
now rewrite Zmult_1_r.
now rewrite Z.shiftl_mul_pow2.
easy. }
simpl.
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).
refine (_ (binary_round_aux_correct' m (sqrt (F2R (Float radix2 (Zpos mx) ex))) mz ez lz _ _ Pz)) ; cycle 1.
clearbody mx'.
destruct Z.sqrtrem as [mz r].
set (lz := if Zeq_bool r 0 then _ else _).
clearbody lz.
intros Bz.
refine (_ (binary_round_aux_correct' m (sqrt (F2R (Float radix2 (Zpos mx) ex))) mz e' lz _ _ _)) ; cycle 1.
now apply Rgt_not_eq, sqrt_lt_R0, F2R_gt_0.
rewrite Rabs_pos_eq.
exact Bz.
apply sqrt_ge_0.
apply Zle_trans with (1 := Zle_min_l _ _).
apply FLT_exp_monotone.
rewrite mag_sqrt_F2R by easy.
apply Zle_refl.
rewrite Rlt_bool_false by apply sqrt_ge_0.
rewrite Rlt_bool_true.
easy.
......
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