From 0f919eb26c68d3882e612a1b3a9df45bee6d3624 Mon Sep 17 00:00:00 2001 From: Guillaume Melquiond Date: Wed, 13 Feb 2019 18:53:17 +0100 Subject: Upgrade embedded version of Flocq to 3.1. Main changes to CompCert outside of Flocq are as follows: - Minimal supported version of Coq is now 8.7, due to Flocq requirements. - Most modifications are due to Z2R being dropped in favor of IZR and to the way Flocq now handles NaNs. - CompCert now correctly handles NaNs for the Risc-V architecture (hopefully). --- flocq/Calc/Bracket.v | 672 +++++++++++++++++++++++++ flocq/Calc/Div.v | 159 ++++++ flocq/Calc/Fcalc_bracket.v | 688 -------------------------- flocq/Calc/Fcalc_digits.v | 63 --- flocq/Calc/Fcalc_div.v | 165 ------- flocq/Calc/Fcalc_ops.v | 163 ------ flocq/Calc/Fcalc_round.v | 1094 ----------------------------------------- flocq/Calc/Fcalc_sqrt.v | 244 --------- flocq/Calc/Operations.v | 164 +++++++ flocq/Calc/Round.v | 1171 ++++++++++++++++++++++++++++++++++++++++++++ flocq/Calc/Sqrt.v | 201 ++++++++ 11 files changed, 2367 insertions(+), 2417 deletions(-) create mode 100644 flocq/Calc/Bracket.v create mode 100644 flocq/Calc/Div.v delete mode 100644 flocq/Calc/Fcalc_bracket.v delete mode 100644 flocq/Calc/Fcalc_digits.v delete mode 100644 flocq/Calc/Fcalc_div.v delete mode 100644 flocq/Calc/Fcalc_ops.v delete mode 100644 flocq/Calc/Fcalc_round.v delete mode 100644 flocq/Calc/Fcalc_sqrt.v create mode 100644 flocq/Calc/Operations.v create mode 100644 flocq/Calc/Round.v create mode 100644 flocq/Calc/Sqrt.v (limited to 'flocq/Calc') diff --git a/flocq/Calc/Bracket.v b/flocq/Calc/Bracket.v new file mode 100644 index 00000000..83714e87 --- /dev/null +++ b/flocq/Calc/Bracket.v @@ -0,0 +1,672 @@ +(** +This file is part of the Flocq formalization of floating-point +arithmetic in Coq: http://flocq.gforge.inria.fr/ + +Copyright (C) 2010-2018 Sylvie Boldo +#
# +Copyright (C) 2010-2018 Guillaume Melquiond + +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. +*) + +(** * Locations: where a real number is positioned with respect to its rounded-down value in an arbitrary format. *) + +Require Import Raux Defs Float_prop. + +Section Fcalc_bracket. + +Variable d u : R. +Hypothesis Hdu : (d < u)%R. + +Inductive location := loc_Exact | loc_Inexact : comparison -> location. + +Variable x : R. + +Definition inbetween_loc := + match Rcompare x d with + | Gt => loc_Inexact (Rcompare x ((d + u) / 2)) + | _ => loc_Exact + end. + +(** Locates a real number with respect to the middle of two other numbers. *) + +Inductive inbetween : location -> Prop := + | inbetween_Exact : x = d -> inbetween loc_Exact + | inbetween_Inexact l : (d < x < u)%R -> Rcompare x ((d + u) / 2)%R = l -> inbetween (loc_Inexact l). + +Theorem inbetween_spec : + (d <= x < u)%R -> inbetween inbetween_loc. +Proof. +intros Hx. +unfold inbetween_loc. +destruct (Rcompare_spec x d) as [H|H|H]. +now elim Rle_not_lt with (1 := proj1 Hx). +now constructor. +constructor. +now split. +easy. +Qed. + +Theorem inbetween_unique : + forall l l', + inbetween l -> inbetween l' -> l = l'. +Proof. +intros l l' Hl Hl'. +inversion_clear Hl ; inversion_clear Hl'. +apply refl_equal. +rewrite H in H0. +elim Rlt_irrefl with (1 := proj1 H0). +rewrite H1 in H. +elim Rlt_irrefl with (1 := proj1 H). +apply f_equal. +now rewrite <- H0. +Qed. + +Section Fcalc_bracket_any. + +Variable l : location. + +Theorem inbetween_bounds : + inbetween l -> + (d <= x < u)%R. +Proof. +intros [Hx|l' Hx Hl] ; clear l. +rewrite Hx. +split. +apply Rle_refl. +exact Hdu. +now split ; try apply Rlt_le. +Qed. + +Theorem inbetween_bounds_not_Eq : + inbetween l -> + l <> loc_Exact -> + (d < x < u)%R. +Proof. +intros [Hx|l' Hx Hl] H. +now elim H. +exact Hx. +Qed. + +End Fcalc_bracket_any. + +Theorem inbetween_distance_inexact : + forall l, + inbetween (loc_Inexact l) -> + Rcompare (x - d) (u - x) = l. +Proof. +intros l Hl. +inversion_clear Hl as [|l' Hl' Hx]. +now rewrite Rcompare_middle. +Qed. + +Theorem inbetween_distance_inexact_abs : + forall l, + inbetween (loc_Inexact l) -> + Rcompare (Rabs (d - x)) (Rabs (u - x)) = l. +Proof. +intros l Hl. +rewrite Rabs_left1. +rewrite Rabs_pos_eq. +rewrite Ropp_minus_distr. +now apply inbetween_distance_inexact. +apply Rle_0_minus. +apply Rlt_le. +apply (inbetween_bounds _ Hl). +apply Rle_minus. +apply (inbetween_bounds _ Hl). +Qed. + +End Fcalc_bracket. + +Theorem inbetween_ex : + forall d u l, + (d < u)%R -> + exists x, + inbetween d u x l. +Proof. +intros d u [|l] Hdu. +exists d. +now constructor. +exists (d + match l with Lt => 1 | Eq => 2 | Gt => 3 end / 4 * (u - d))%R. +constructor. +(* *) +set (v := (match l with Lt => 1 | Eq => 2 | Gt => 3 end / 4)%R). +assert (0 < v < 1)%R. +split. +unfold v, Rdiv. +apply Rmult_lt_0_compat. +case l ; now apply IZR_lt. +apply Rinv_0_lt_compat. +now apply IZR_lt. +unfold v, Rdiv. +apply Rmult_lt_reg_r with 4%R. +now apply IZR_lt. +rewrite Rmult_assoc, Rinv_l. +rewrite Rmult_1_r, Rmult_1_l. +case l ; now apply IZR_lt. +apply Rgt_not_eq. +now apply IZR_lt. +split. +apply Rplus_lt_reg_r with (d * (v - 1))%R. +ring_simplify. +rewrite Rmult_comm. +now apply Rmult_lt_compat_l. +apply Rplus_lt_reg_l with (-u * v)%R. +replace (- u * v + (d + v * (u - d)))%R with (d * (1 - v))%R by ring. +replace (- u * v + u)%R with (u * (1 - v))%R by ring. +apply Rmult_lt_compat_r. +apply Rplus_lt_reg_r with v. +now ring_simplify. +exact Hdu. +(* *) +set (v := (match l with Lt => 1 | Eq => 2 | Gt => 3 end)%R). +rewrite <- (Rcompare_plus_r (- (d + u) / 2)). +rewrite <- (Rcompare_mult_r 4). +2: now apply IZR_lt. +replace (((d + u) / 2 + - (d + u) / 2) * 4)%R with ((u - d) * 0)%R by field. +replace ((d + v / 4 * (u - d) + - (d + u) / 2) * 4)%R with ((u - d) * (v - 2))%R by field. +rewrite Rcompare_mult_l. +2: now apply Rlt_Rminus. +rewrite <- (Rcompare_plus_r 2). +ring_simplify (v - 2 + 2)%R (0 + 2)%R. +unfold v. +case l ; apply Rcompare_IZR. +Qed. + +Section Fcalc_bracket_step. + +Variable start step : R. +Variable nb_steps : Z. +Variable Hstep : (0 < step)%R. + +Lemma ordered_steps : + forall k, + (start + IZR k * step < start + IZR (k + 1) * step)%R. +Proof. +intros k. +apply Rplus_lt_compat_l. +apply Rmult_lt_compat_r. +exact Hstep. +apply IZR_lt. +apply Zlt_succ. +Qed. + +Lemma middle_range : + forall k, + ((start + (start + IZR k * step)) / 2 = start + (IZR k / 2 * step))%R. +Proof. +intros k. +field. +Qed. + +Hypothesis (Hnb_steps : (1 < nb_steps)%Z). + +Lemma inbetween_step_not_Eq : + forall x k l l', + inbetween (start + IZR k * step) (start + IZR (k + 1) * step) x l -> + (0 < k < nb_steps)%Z -> + Rcompare x (start + (IZR nb_steps / 2 * step))%R = l' -> + inbetween start (start + IZR nb_steps * step) x (loc_Inexact l'). +Proof. +intros x k l l' Hx Hk Hl'. +constructor. +(* . *) +assert (Hx' := inbetween_bounds _ _ (ordered_steps _) _ _ Hx). +split. +apply Rlt_le_trans with (2 := proj1 Hx'). +rewrite <- (Rplus_0_r start) at 1. +apply Rplus_lt_compat_l. +apply Rmult_lt_0_compat. +now apply IZR_lt. +exact Hstep. +apply Rlt_le_trans with (1 := proj2 Hx'). +apply Rplus_le_compat_l. +apply Rmult_le_compat_r. +now apply Rlt_le. +apply IZR_le. +omega. +(* . *) +now rewrite middle_range. +Qed. + +Theorem inbetween_step_Lo : + forall x k l, + inbetween (start + IZR k * step) (start + IZR (k + 1) * step) x l -> + (0 < k)%Z -> (2 * k + 1 < nb_steps)%Z -> + inbetween start (start + IZR nb_steps * step) x (loc_Inexact Lt). +Proof. +intros x k l Hx Hk1 Hk2. +apply inbetween_step_not_Eq with (1 := Hx). +omega. +apply Rcompare_Lt. +assert (Hx' := inbetween_bounds _ _ (ordered_steps _) _ _ Hx). +apply Rlt_le_trans with (1 := proj2 Hx'). +apply Rcompare_not_Lt_inv. +rewrite Rcompare_plus_l, Rcompare_mult_r, Rcompare_half_l. +apply Rcompare_not_Lt. +rewrite <- mult_IZR. +apply IZR_le. +omega. +exact Hstep. +Qed. + +Theorem inbetween_step_Hi : + forall x k l, + inbetween (start + IZR k * step) (start + IZR (k + 1) * step) x l -> + (nb_steps < 2 * k)%Z -> (k < nb_steps)%Z -> + inbetween start (start + IZR nb_steps * step) x (loc_Inexact Gt). +Proof. +intros x k l Hx Hk1 Hk2. +apply inbetween_step_not_Eq with (1 := Hx). +omega. +apply Rcompare_Gt. +assert (Hx' := inbetween_bounds _ _ (ordered_steps _) _ _ Hx). +apply Rlt_le_trans with (2 := proj1 Hx'). +apply Rcompare_Lt_inv. +rewrite Rcompare_plus_l, Rcompare_mult_r, Rcompare_half_l. +apply Rcompare_Lt. +rewrite <- mult_IZR. +apply IZR_lt. +omega. +exact Hstep. +Qed. + +Theorem inbetween_step_Lo_not_Eq : + forall x l, + inbetween start (start + step) x l -> + l <> loc_Exact -> + inbetween start (start + IZR nb_steps * step) x (loc_Inexact Lt). +Proof. +intros x l Hx Hl. +assert (Hx' := inbetween_bounds_not_Eq _ _ _ _ Hx Hl). +constructor. +(* . *) +split. +apply Hx'. +apply Rlt_trans with (1 := proj2 Hx'). +apply Rplus_lt_compat_l. +rewrite <- (Rmult_1_l step) at 1. +apply Rmult_lt_compat_r. +exact Hstep. +now apply IZR_lt. +(* . *) +apply Rcompare_Lt. +apply Rlt_le_trans with (1 := proj2 Hx'). +rewrite middle_range. +apply Rcompare_not_Lt_inv. +rewrite <- (Rmult_1_l step) at 2. +rewrite Rcompare_plus_l, Rcompare_mult_r, Rcompare_half_l. +rewrite Rmult_1_r. +apply Rcompare_not_Lt. +apply IZR_le. +now apply (Zlt_le_succ 1). +exact Hstep. +Qed. + +Lemma middle_odd : + forall k, + (2 * k + 1 = nb_steps)%Z -> + (((start + IZR k * step) + (start + IZR (k + 1) * step))/2 = start + IZR nb_steps /2 * step)%R. +Proof. +intros k Hk. +rewrite <- Hk. +rewrite 2!plus_IZR, mult_IZR. +simpl. field. +Qed. + +Theorem inbetween_step_any_Mi_odd : + forall x k l, + inbetween (start + IZR k * step) (start + IZR (k + 1) * step) x (loc_Inexact l) -> + (2 * k + 1 = nb_steps)%Z -> + inbetween start (start + IZR nb_steps * step) x (loc_Inexact l). +Proof. +intros x k l Hx Hk. +apply inbetween_step_not_Eq with (1 := Hx). +omega. +inversion_clear Hx as [|l' _ Hl]. +now rewrite (middle_odd _ Hk) in Hl. +Qed. + +Theorem inbetween_step_Lo_Mi_Eq_odd : + forall x k, + inbetween (start + IZR k * step) (start + IZR (k + 1) * step) x loc_Exact -> + (2 * k + 1 = nb_steps)%Z -> + inbetween start (start + IZR nb_steps * step) x (loc_Inexact Lt). +Proof. +intros x k Hx Hk. +apply inbetween_step_not_Eq with (1 := Hx). +omega. +inversion_clear Hx as [Hl|]. +rewrite Hl. +rewrite Rcompare_plus_l, Rcompare_mult_r, Rcompare_half_r. +apply Rcompare_Lt. +rewrite <- mult_IZR. +apply IZR_lt. +rewrite <- Hk. +apply Zlt_succ. +exact Hstep. +Qed. + +Theorem inbetween_step_Hi_Mi_even : + forall x k l, + inbetween (start + IZR k * step) (start + IZR (k + 1) * step) x l -> + l <> loc_Exact -> + (2 * k = nb_steps)%Z -> + inbetween start (start + IZR nb_steps * step) x (loc_Inexact Gt). +Proof. +intros x k l Hx Hl Hk. +apply inbetween_step_not_Eq with (1 := Hx). +omega. +apply Rcompare_Gt. +assert (Hx' := inbetween_bounds_not_Eq _ _ _ _ Hx Hl). +apply Rle_lt_trans with (2 := proj1 Hx'). +apply Rcompare_not_Lt_inv. +rewrite Rcompare_plus_l, Rcompare_mult_r, Rcompare_half_r. +rewrite <- mult_IZR. +apply Rcompare_not_Lt. +apply IZR_le. +rewrite Hk. +apply Z.le_refl. +exact Hstep. +Qed. + +Theorem inbetween_step_Mi_Mi_even : + forall x k, + inbetween (start + IZR k * step) (start + IZR (k + 1) * step) x loc_Exact -> + (2 * k = nb_steps)%Z -> + inbetween start (start + IZR nb_steps * step) x (loc_Inexact Eq). +Proof. +intros x k Hx Hk. +apply inbetween_step_not_Eq with (1 := Hx). +omega. +apply Rcompare_Eq. +inversion_clear Hx as [Hx'|]. +rewrite Hx', <- Hk, mult_IZR. +field. +Qed. + +(** Computes a new location when the interval containing a real + number is split into nb_steps subintervals and the real is + in the k-th one. (Even radix.) *) + +Definition new_location_even k l := + if Zeq_bool k 0 then + match l with loc_Exact => l | _ => loc_Inexact Lt end + else + loc_Inexact + match Z.compare (2 * k) nb_steps with + | Lt => Lt + | Eq => match l with loc_Exact => Eq | _ => Gt end + | Gt => Gt + end. + +Theorem new_location_even_correct : + Z.even nb_steps = true -> + forall x k l, (0 <= k < nb_steps)%Z -> + inbetween (start + IZR k * step) (start + IZR (k + 1) * step) x l -> + inbetween start (start + IZR nb_steps * step) x (new_location_even k l). +Proof. +intros He x k l Hk Hx. +unfold new_location_even. +destruct (Zeq_bool_spec k 0) as [Hk0|Hk0]. +(* k = 0 *) +rewrite Hk0 in Hx. +rewrite Rmult_0_l, Rplus_0_r, Rmult_1_l in Hx. +set (l' := match l with loc_Exact => l | _ => loc_Inexact Lt end). +assert ((l = loc_Exact /\ l' = loc_Exact) \/ (l <> loc_Exact /\ l' = loc_Inexact Lt)). +unfold l' ; case l ; try (now left) ; right ; now split. +destruct H as [(H1,H2)|(H1,H2)] ; rewrite H2. +constructor. +rewrite H1 in Hx. +now inversion_clear Hx. +now apply inbetween_step_Lo_not_Eq with (2 := H1). +(* k <> 0 *) +destruct (Zcompare_spec (2 * k) nb_steps) as [Hk1|Hk1|Hk1]. +(* . 2 * k < nb_steps *) +apply inbetween_step_Lo with (1 := Hx). +omega. +destruct (Zeven_ex nb_steps). +rewrite He in H. +omega. +(* . 2 * k = nb_steps *) +set (l' := match l with loc_Exact => Eq | _ => Gt end). +assert ((l = loc_Exact /\ l' = Eq) \/ (l <> loc_Exact /\ l' = Gt)). +unfold l' ; case l ; try (now left) ; right ; now split. +destruct H as [(H1,H2)|(H1,H2)] ; rewrite H2. +rewrite H1 in Hx. +now apply inbetween_step_Mi_Mi_even with (1 := Hx). +now apply inbetween_step_Hi_Mi_even with (1 := Hx). +(* . 2 * k > nb_steps *) +apply inbetween_step_Hi with (1 := Hx). +exact Hk1. +apply Hk. +Qed. + +(** Computes a new location when the interval containing a real + number is split into nb_steps subintervals and the real is + in the k-th one. (Odd radix.) *) + +Definition new_location_odd k l := + if Zeq_bool k 0 then + match l with loc_Exact => l | _ => loc_Inexact Lt end + else + loc_Inexact + match Z.compare (2 * k + 1) nb_steps with + | Lt => Lt + | Eq => match l with loc_Inexact l => l | loc_Exact => Lt end + | Gt => Gt + end. + +Theorem new_location_odd_correct : + Z.even nb_steps = false -> + forall x k l, (0 <= k < nb_steps)%Z -> + inbetween (start + IZR k * step) (start + IZR (k + 1) * step) x l -> + inbetween start (start + IZR nb_steps * step) x (new_location_odd k l). +Proof. +intros Ho x k l Hk Hx. +unfold new_location_odd. +destruct (Zeq_bool_spec k 0) as [Hk0|Hk0]. +(* k = 0 *) +rewrite Hk0 in Hx. +rewrite Rmult_0_l, Rplus_0_r, Rmult_1_l in Hx. +set (l' := match l with loc_Exact => l | _ => loc_Inexact Lt end). +assert ((l = loc_Exact /\ l' = loc_Exact) \/ (l <> loc_Exact /\ l' = loc_Inexact Lt)). +unfold l' ; case l ; try (now left) ; right ; now split. +destruct H as [(H1,H2)|(H1,H2)] ; rewrite H2. +constructor. +rewrite H1 in Hx. +now inversion_clear Hx. +now apply inbetween_step_Lo_not_Eq with (2 := H1). +(* k <> 0 *) +destruct (Zcompare_spec (2 * k + 1) nb_steps) as [Hk1|Hk1|Hk1]. +(* . 2 * k + 1 < nb_steps *) +apply inbetween_step_Lo with (1 := Hx) (3 := Hk1). +omega. +(* . 2 * k + 1 = nb_steps *) +destruct l. +apply inbetween_step_Lo_Mi_Eq_odd with (1 := Hx) (2 := Hk1). +apply inbetween_step_any_Mi_odd with (1 := Hx) (2 := Hk1). +(* . 2 * k + 1 > nb_steps *) +apply inbetween_step_Hi with (1 := Hx). +destruct (Zeven_ex nb_steps). +rewrite Ho in H. +omega. +apply Hk. +Qed. + +Definition new_location := + if Z.even nb_steps then new_location_even else new_location_odd. + +Theorem new_location_correct : + forall x k l, (0 <= k < nb_steps)%Z -> + inbetween (start + IZR k * step) (start + IZR (k + 1) * step) x l -> + inbetween start (start + IZR nb_steps * step) x (new_location k l). +Proof. +intros x k l Hk Hx. +unfold new_location. +generalize (refl_equal nb_steps) (Z.le_lt_trans _ _ _ (proj1 Hk) (proj2 Hk)). +pattern nb_steps at 2 3 5 ; case nb_steps. +now intros _. +(* . *) +intros [p|p|] Hp _. +apply new_location_odd_correct with (2 := Hk) (3 := Hx). +now rewrite Hp. +apply new_location_even_correct with (2 := Hk) (3 := Hx). +now rewrite Hp. +now rewrite Hp in Hnb_steps. +(* . *) +now intros p _. +Qed. + +End Fcalc_bracket_step. + +Section Fcalc_bracket_scale. + +Lemma inbetween_mult_aux : + forall x d s, + ((x * s + d * s) / 2 = (x + d) / 2 * s)%R. +Proof. +intros x d s. +field. +Qed. + +Theorem inbetween_mult_compat : + forall x d u l s, + (0 < s)%R -> + inbetween x d u l -> + inbetween (x * s) (d * s) (u * s) l. +Proof. +intros x d u l s Hs [Hx|l' Hx Hl] ; constructor. +now rewrite Hx. +now split ; apply Rmult_lt_compat_r. +rewrite inbetween_mult_aux. +now rewrite Rcompare_mult_r. +Qed. + +Theorem inbetween_mult_reg : + forall x d u l s, + (0 < s)%R -> + inbetween (x * s) (d * s) (u * s) l -> + inbetween x d u l. +Proof. +intros x d u l s Hs [Hx|l' Hx Hl] ; constructor. +apply Rmult_eq_reg_r with (1 := Hx). +now apply Rgt_not_eq. +now split ; apply Rmult_lt_reg_r with s. +rewrite <- Rcompare_mult_r with (1 := Hs). +now rewrite inbetween_mult_aux in Hl. +Qed. + +End Fcalc_bracket_scale. + +Section Fcalc_bracket_generic. + +Variable beta : radix. +Notation bpow e := (bpow beta e). + +(** Specialization of inbetween for two consecutive floating-point numbers. *) + +Definition inbetween_float m e x l := + inbetween (F2R (Float beta m e)) (F2R (Float beta (m + 1) e)) x l. + +Theorem inbetween_float_bounds : + forall x m e l, + inbetween_float m e x l -> + (F2R (Float beta m e) <= x < F2R (Float beta (m + 1) e))%R. +Proof. +intros x m e l [Hx|l' Hx Hl]. +rewrite Hx. +split. +apply Rle_refl. +apply F2R_lt. +apply Zlt_succ. +split. +now apply Rlt_le. +apply Hx. +Qed. + +(** Specialization of inbetween for two consecutive integers. *) + +Definition inbetween_int m x l := + inbetween (IZR m) (IZR (m + 1)) x l. + +Theorem inbetween_float_new_location : + forall x m e l k, + (0 < k)%Z -> + inbetween_float m e x l -> + inbetween_float (Z.div m (Zpower beta k)) (e + k) x (new_location (Zpower beta k) (Zmod m (Zpower beta k)) l). +Proof. +intros x m e l k Hk Hx. +unfold inbetween_float in *. +assert (Hr: forall m, F2R (Float beta m (e + k)) = F2R (Float beta (m * Zpower beta k) e)). +clear -Hk. intros m. +rewrite (F2R_change_exp beta e). +apply (f_equal (fun r => F2R (Float beta (m * Zpower _ r) e))). +ring. +omega. +assert (Hp: (Zpower beta k > 0)%Z). +apply Z.lt_gt. +apply Zpower_gt_0. +now apply Zlt_le_weak. +(* . *) +rewrite 2!Hr. +rewrite Zmult_plus_distr_l, Zmult_1_l. +unfold F2R at 2. simpl. +rewrite plus_IZR, Rmult_plus_distr_r. +apply new_location_correct. +apply bpow_gt_0. +now apply Zpower_gt_1. +now apply Z_mod_lt. +rewrite <- 2!Rmult_plus_distr_r, <- 2!plus_IZR. +rewrite Zmult_comm, Zplus_assoc. +now rewrite <- Z_div_mod_eq. +Qed. + +Theorem inbetween_float_new_location_single : + forall x m e l, + inbetween_float m e x l -> + inbetween_float (Z.div m beta) (e + 1) x (new_location beta (Zmod m beta) l). +Proof. +intros x m e l Hx. +replace (radix_val beta) with (Zpower beta 1). +now apply inbetween_float_new_location. +apply Zmult_1_r. +Qed. + +Theorem inbetween_float_ex : + forall m e l, + exists x, + inbetween_float m e x l. +Proof. +intros m e l. +apply inbetween_ex. +apply F2R_lt. +apply Zlt_succ. +Qed. + +Theorem inbetween_float_unique : + forall x e m l m' l', + inbetween_float m e x l -> + inbetween_float m' e x l' -> + m = m' /\ l = l'. +Proof. +intros x e m l m' l' H H'. +refine ((fun Hm => conj Hm _) _). +rewrite <- Hm in H'. clear -H H'. +apply inbetween_unique with (1 := H) (2 := H'). +destruct (inbetween_float_bounds x m e l H) as (H1,H2). +destruct (inbetween_float_bounds x m' e l' H') as (H3,H4). +cut (m < m' + 1 /\ m' < m + 1)%Z. clear ; omega. +now split ; apply lt_F2R with beta e ; apply Rle_lt_trans with x. +Qed. + +End Fcalc_bracket_generic. diff --git a/flocq/Calc/Div.v b/flocq/Calc/Div.v new file mode 100644 index 00000000..65195562 --- /dev/null +++ b/flocq/Calc/Div.v @@ -0,0 +1,159 @@ +(** +This file is part of the Flocq formalization of floating-point +arithmetic in Coq: http://flocq.gforge.inria.fr/ + +Copyright (C) 2010-2018 Sylvie Boldo +#
# +Copyright (C) 2010-2018 Guillaume Melquiond + +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. +*) + +(** * Helper function and theorem for computing the rounded quotient of two floating-point numbers. *) + +Require Import Raux Defs Generic_fmt Float_prop Digits Bracket. + +Set Implicit Arguments. +Set Strongly Strict Implicit. + +Section Fcalc_div. + +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 quotient of the + input floating-point numbers. + +The algorithm performs the following steps: +- Shift dividend mantissa so that it has at least p2 + p digits. +- Perform the Euclidean division. +- Compute the position according to the division remainder. + +Complexity is fine as long as p1 <= 2p and p2 <= p. +*) + +Lemma mag_div_F2R : + forall m1 e1 m2 e2, + (0 < m1)%Z -> (0 < m2)%Z -> + let e := ((Zdigits beta m1 + e1) - (Zdigits beta m2 + e2))%Z in + (e <= mag beta (F2R (Float beta m1 e1) / F2R (Float beta m2 e2)) <= e + 1)%Z. +Proof. +intros m1 e1 m2 e2 Hm1 Hm2. +rewrite <- (mag_F2R_Zdigits beta m1 e1) by now apply Zgt_not_eq. +rewrite <- (mag_F2R_Zdigits beta m2 e2) by now apply Zgt_not_eq. +apply mag_div. +now apply F2R_neq_0, Zgt_not_eq. +now apply F2R_neq_0, Zgt_not_eq. +Qed. + +Definition Fdiv_core m1 e1 m2 e2 e := + let (m1', m2') := + if Zle_bool e (e1 - e2)%Z + then (m1 * Zpower beta (e1 - e2 - e), m2)%Z + else (m1, m2 * Zpower beta (e - (e1 - e2)))%Z in + let '(q, r) := Z.div_eucl m1' m2' in + (q, new_location m2' r loc_Exact). + +Theorem Fdiv_core_correct : + forall m1 e1 m2 e2 e, + (0 < m1)%Z -> (0 < m2)%Z -> + let '(m, l) := Fdiv_core m1 e1 m2 e2 e in + inbetween_float beta m e (F2R (Float beta m1 e1) / F2R (Float beta m2 e2)) l. +Proof. +intros m1 e1 m2 e2 e Hm1 Hm2. +unfold Fdiv_core. +match goal with |- context [if ?b then ?b1 else ?b2] => set (m12 := if b then b1 else b2) end. +case_eq m12 ; intros m1' m2' Hm. +assert ((F2R (Float beta m1 e1) / F2R (Float beta m2 e2) = IZR m1' / IZR m2' * bpow e)%R /\ (0 < m2')%Z) as [Hf Hm2']. +{ unfold F2R, Zminus ; simpl. + destruct (Zle_bool e (e1 - e2)) eqn:He' ; injection Hm ; intros ; subst. + - split ; try easy. + apply Zle_bool_imp_le in He'. + rewrite mult_IZR, IZR_Zpower by omega. + unfold Zminus ; rewrite 2!bpow_plus, 2!bpow_opp. + field. + repeat split ; try apply Rgt_not_eq, bpow_gt_0. + now apply IZR_neq, Zgt_not_eq. + - apply Z.leb_gt in He'. + split ; cycle 1. + { apply Z.mul_pos_pos with (1 := Hm2). + apply Zpower_gt_0 ; omega. } + rewrite mult_IZR, IZR_Zpower by omega. + unfold Zminus ; rewrite bpow_plus, bpow_opp, bpow_plus, bpow_opp. + field. + repeat split ; try apply Rgt_not_eq, bpow_gt_0. + now apply IZR_neq, Zgt_not_eq. } +clearbody m12 ; clear Hm Hm1 Hm2. +generalize (Z_div_mod m1' m2' (Z.lt_gt _ _ Hm2')). +destruct (Z.div_eucl m1' m2') as (q, r). +intros (Hq, Hr). +rewrite Hf. +unfold inbetween_float, F2R. simpl. +rewrite Hq, 2!plus_IZR, mult_IZR. +apply inbetween_mult_compat. + apply bpow_gt_0. +destruct (Z_lt_le_dec 1 m2') as [Hm2''|Hm2'']. +- replace 1%R with (IZR m2' * /IZR m2')%R. + apply new_location_correct ; try easy. + apply Rinv_0_lt_compat. + now apply IZR_lt. + constructor. + field. + now apply IZR_neq, Zgt_not_eq. + field. + now apply IZR_neq, Zgt_not_eq. +- assert (r = 0 /\ m2' = 1)%Z as [-> ->] by (clear -Hr Hm2'' ; omega). + unfold Rdiv. + rewrite Rmult_1_l, Rplus_0_r, Rinv_1, Rmult_1_r. + now constructor. +Qed. + +Definition Fdiv (x y : float beta) := + let (m1, e1) := x in + let (m2, e2) := y in + let e' := ((Zdigits beta m1 + e1) - (Zdigits beta m2 + e2))%Z in + let e := Z.min (Z.min (fexp e') (fexp (e' + 1))) (e1 - e2) in + let '(m, l) := Fdiv_core m1 e1 m2 e2 e in + (m, e, l). + +Theorem Fdiv_correct : + forall x y, + (0 < F2R x)%R -> (0 < F2R y)%R -> + let '(m, e, l) := Fdiv x y in + (e <= cexp beta fexp (F2R x / F2R y))%Z /\ + inbetween_float beta m e (F2R x / F2R y) l. +Proof. +intros [m1 e1] [m2 e2] Hm1 Hm2. +apply gt_0_F2R in Hm1. +apply gt_0_F2R in Hm2. +unfold Fdiv. +generalize (mag_div_F2R m1 e1 m2 e2 Hm1 Hm2). +set (e := Zminus _ _). +set (e' := Z.min (Z.min (fexp e) (fexp (e + 1))) (e1 - e2)). +intros [H1 H2]. +generalize (Fdiv_core_correct m1 e1 m2 e2 e' Hm1 Hm2). +destruct Fdiv_core as [m' l]. +apply conj. +apply Z.le_trans with (1 := Z.le_min_l _ _). +unfold cexp. +destruct (Zle_lt_or_eq _ _ H1) as [H|H]. +- replace (fexp (mag _ _)) with (fexp (e + 1)). + apply Z.le_min_r. + clear -H1 H2 H ; apply f_equal ; omega. +- replace (fexp (mag _ _)) with (fexp e). + apply Z.le_min_l. + clear -H1 H2 H ; apply f_equal ; omega. +Qed. + +End Fcalc_div. diff --git a/flocq/Calc/Fcalc_bracket.v b/flocq/Calc/Fcalc_bracket.v deleted file mode 100644 index 03ef7bd3..00000000 --- a/flocq/Calc/Fcalc_bracket.v +++ /dev/null @@ -1,688 +0,0 @@ -(** -This file is part of the Flocq formalization of floating-point -arithmetic in Coq: http://flocq.gforge.inria.fr/ - -Copyright (C) 2010-2013 Sylvie Boldo -#
# -Copyright (C) 2010-2013 Guillaume Melquiond - -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. -*) - -(** * Locations: where a real number is positioned with respect to its rounded-down value in an arbitrary format. *) - -Require Import Fcore_Raux. -Require Import Fcore_defs. -Require Import Fcore_float_prop. - -Section Fcalc_bracket. - -Variable d u : R. -Hypothesis Hdu : (d < u)%R. - -Inductive location := loc_Exact | loc_Inexact : comparison -> location. - -Variable x : R. - -Definition inbetween_loc := - match Rcompare x d with - | Gt => loc_Inexact (Rcompare x ((d + u) / 2)) - | _ => loc_Exact - end. - -(** Locates a real number with respect to the middle of two other numbers. *) - -Inductive inbetween : location -> Prop := - | inbetween_Exact : x = d -> inbetween loc_Exact - | inbetween_Inexact l : (d < x < u)%R -> Rcompare x ((d + u) / 2)%R = l -> inbetween (loc_Inexact l). - -Theorem inbetween_spec : - (d <= x < u)%R -> inbetween inbetween_loc. -Proof. -intros Hx. -unfold inbetween_loc. -destruct (Rcompare_spec x d) as [H|H|H]. -now elim Rle_not_lt with (1 := proj1 Hx). -now constructor. -constructor. -now split. -easy. -Qed. - -Theorem inbetween_unique : - forall l l', - inbetween l -> inbetween l' -> l = l'. -Proof. -intros l l' Hl Hl'. -inversion_clear Hl ; inversion_clear Hl'. -apply refl_equal. -rewrite H in H0. -elim Rlt_irrefl with (1 := proj1 H0). -rewrite H1 in H. -elim Rlt_irrefl with (1 := proj1 H). -apply f_equal. -now rewrite <- H0. -Qed. - -Section Fcalc_bracket_any. - -Variable l : location. - -Theorem inbetween_bounds : - inbetween l -> - (d <= x < u)%R. -Proof. -intros [Hx|l' Hx Hl] ; clear l. -rewrite Hx. -split. -apply Rle_refl. -exact Hdu. -now split ; try apply Rlt_le. -Qed. - -Theorem inbetween_bounds_not_Eq : - inbetween l -> - l <> loc_Exact -> - (d < x < u)%R. -Proof. -intros [Hx|l' Hx Hl] H. -now elim H. -exact Hx. -Qed. - -End Fcalc_bracket_any. - -Theorem inbetween_distance_inexact : - forall l, - inbetween (loc_Inexact l) -> - Rcompare (x - d) (u - x) = l. -Proof. -intros l Hl. -inversion_clear Hl as [|l' Hl' Hx]. -now rewrite Rcompare_middle. -Qed. - -Theorem inbetween_distance_inexact_abs : - forall l, - inbetween (loc_Inexact l) -> - Rcompare (Rabs (d - x)) (Rabs (u - x)) = l. -Proof. -intros l Hl. -rewrite Rabs_left1. -rewrite Rabs_pos_eq. -rewrite Ropp_minus_distr. -now apply inbetween_distance_inexact. -apply Rle_0_minus. -apply Rlt_le. -apply (inbetween_bounds _ Hl). -apply Rle_minus. -apply (inbetween_bounds _ Hl). -Qed. - -End Fcalc_bracket. - -Theorem inbetween_ex : - forall d u l, - (d < u)%R -> - exists x, - inbetween d u x l. -Proof. -intros d u [|l] Hdu. -exists d. -now constructor. -exists (d + match l with Lt => 1 | Eq => 2 | Gt => 3 end / 4 * (u - d))%R. -constructor. -(* *) -set (v := (match l with Lt => 1 | Eq => 2 | Gt => 3 end / 4)%R). -assert (0 < v < 1)%R. -split. -unfold v, Rdiv. -apply Rmult_lt_0_compat. -case l. -now apply (Z2R_lt 0 2). -now apply (Z2R_lt 0 1). -now apply (Z2R_lt 0 3). -apply Rinv_0_lt_compat. -now apply (Z2R_lt 0 4). -unfold v, Rdiv. -apply Rmult_lt_reg_r with 4%R. -now apply (Z2R_lt 0 4). -rewrite Rmult_assoc, Rinv_l. -rewrite Rmult_1_r, Rmult_1_l. -case l. -now apply (Z2R_lt 2 4). -now apply (Z2R_lt 1 4). -now apply (Z2R_lt 3 4). -apply Rgt_not_eq. -now apply (Z2R_lt 0 4). -split. -apply Rplus_lt_reg_r with (d * (v - 1))%R. -ring_simplify. -rewrite Rmult_comm. -now apply Rmult_lt_compat_l. -apply Rplus_lt_reg_l with (-u * v)%R. -replace (- u * v + (d + v * (u - d)))%R with (d * (1 - v))%R by ring. -replace (- u * v + u)%R with (u * (1 - v))%R by ring. -apply Rmult_lt_compat_r. -apply Rplus_lt_reg_r with v. -now ring_simplify. -exact Hdu. -(* *) -set (v := (match l with Lt => 1 | Eq => 2 | Gt => 3 end)%R). -rewrite <- (Rcompare_plus_r (- (d + u) / 2)). -rewrite <- (Rcompare_mult_r 4). -2: now apply (Z2R_lt 0 4). -replace (((d + u) / 2 + - (d + u) / 2) * 4)%R with ((u - d) * 0)%R by field. -replace ((d + v / 4 * (u - d) + - (d + u) / 2) * 4)%R with ((u - d) * (v - 2))%R by field. -rewrite Rcompare_mult_l. -2: now apply Rlt_Rminus. -rewrite <- (Rcompare_plus_r 2). -ring_simplify (v - 2 + 2)%R (0 + 2)%R. -unfold v. -case l. -exact (Rcompare_Z2R 2 2). -exact (Rcompare_Z2R 1 2). -exact (Rcompare_Z2R 3 2). -Qed. - -Section Fcalc_bracket_step. - -Variable start step : R. -Variable nb_steps : Z. -Variable Hstep : (0 < step)%R. - -Lemma ordered_steps : - forall k, - (start + Z2R k * step < start + Z2R (k + 1) * step)%R. -Proof. -intros k. -apply Rplus_lt_compat_l. -apply Rmult_lt_compat_r. -exact Hstep. -apply Z2R_lt. -apply Zlt_succ. -Qed. - -Lemma middle_range : - forall k, - ((start + (start + Z2R k * step)) / 2 = start + (Z2R k / 2 * step))%R. -Proof. -intros k. -field. -Qed. - -Hypothesis (Hnb_steps : (1 < nb_steps)%Z). - -Lemma inbetween_step_not_Eq : - forall x k l l', - inbetween (start + Z2R k * step) (start + Z2R (k + 1) * step) x l -> - (0 < k < nb_steps)%Z -> - Rcompare x (start + (Z2R nb_steps / 2 * step))%R = l' -> - inbetween start (start + Z2R nb_steps * step) x (loc_Inexact l'). -Proof. -intros x k l l' Hx Hk Hl'. -constructor. -(* . *) -assert (Hx' := inbetween_bounds _ _ (ordered_steps _) _ _ Hx). -split. -apply Rlt_le_trans with (2 := proj1 Hx'). -rewrite <- (Rplus_0_r start) at 1. -apply Rplus_lt_compat_l. -apply Rmult_lt_0_compat. -now apply (Z2R_lt 0). -exact Hstep. -apply Rlt_le_trans with (1 := proj2 Hx'). -apply Rplus_le_compat_l. -apply Rmult_le_compat_r. -now apply Rlt_le. -apply Z2R_le. -omega. -(* . *) -now rewrite middle_range. -Qed. - -Theorem inbetween_step_Lo : - forall x k l, - inbetween (start + Z2R k * step) (start + Z2R (k + 1) * step) x l -> - (0 < k)%Z -> (2 * k + 1 < nb_steps)%Z -> - inbetween start (start + Z2R nb_steps * step) x (loc_Inexact Lt). -Proof. -intros x k l Hx Hk1 Hk2. -apply inbetween_step_not_Eq with (1 := Hx). -omega. -apply Rcompare_Lt. -assert (Hx' := inbetween_bounds _ _ (ordered_steps _) _ _ Hx). -apply Rlt_le_trans with (1 := proj2 Hx'). -apply Rcompare_not_Lt_inv. -rewrite Rcompare_plus_l, Rcompare_mult_r, Rcompare_half_l. -apply Rcompare_not_Lt. -change 2%R with (Z2R 2). -rewrite <- Z2R_mult. -apply Z2R_le. -omega. -exact Hstep. -Qed. - -Theorem inbetween_step_Hi : - forall x k l, - inbetween (start + Z2R k * step) (start + Z2R (k + 1) * step) x l -> - (nb_steps < 2 * k)%Z -> (k < nb_steps)%Z -> - inbetween start (start + Z2R nb_steps * step) x (loc_Inexact Gt). -Proof. -intros x k l Hx Hk1 Hk2. -apply inbetween_step_not_Eq with (1 := Hx). -omega. -apply Rcompare_Gt. -assert (Hx' := inbetween_bounds _ _ (ordered_steps _) _ _ Hx). -apply Rlt_le_trans with (2 := proj1 Hx'). -apply Rcompare_Lt_inv. -rewrite Rcompare_plus_l, Rcompare_mult_r, Rcompare_half_l. -apply Rcompare_Lt. -change 2%R with (Z2R 2). -rewrite <- Z2R_mult. -apply Z2R_lt. -omega. -exact Hstep. -Qed. - -Theorem inbetween_step_Lo_not_Eq : - forall x l, - inbetween start (start + step) x l -> - l <> loc_Exact -> - inbetween start (start + Z2R nb_steps * step) x (loc_Inexact Lt). -Proof. -intros x l Hx Hl. -assert (Hx' := inbetween_bounds_not_Eq _ _ _ _ Hx Hl). -constructor. -(* . *) -split. -apply Hx'. -apply Rlt_trans with (1 := proj2 Hx'). -apply Rplus_lt_compat_l. -rewrite <- (Rmult_1_l step) at 1. -apply Rmult_lt_compat_r. -exact Hstep. -now apply (Z2R_lt 1). -(* . *) -apply Rcompare_Lt. -apply Rlt_le_trans with (1 := proj2 Hx'). -rewrite middle_range. -apply Rcompare_not_Lt_inv. -rewrite <- (Rmult_1_l step) at 2. -rewrite Rcompare_plus_l, Rcompare_mult_r, Rcompare_half_l. -rewrite Rmult_1_r. -apply Rcompare_not_Lt. -apply (Z2R_le 2). -now apply (Zlt_le_succ 1). -exact Hstep. -Qed. - -Lemma middle_odd : - forall k, - (2 * k + 1 = nb_steps)%Z -> - (((start + Z2R k * step) + (start + Z2R (k + 1) * step))/2 = start + Z2R nb_steps /2 * step)%R. -Proof. -intros k Hk. -rewrite <- Hk. -rewrite 2!Z2R_plus, Z2R_mult. -simpl. field. -Qed. - -Theorem inbetween_step_any_Mi_odd : - forall x k l, - inbetween (start + Z2R k * step) (start + Z2R (k + 1) * step) x (loc_Inexact l) -> - (2 * k + 1 = nb_steps)%Z -> - inbetween start (start + Z2R nb_steps * step) x (loc_Inexact l). -Proof. -intros x k l Hx Hk. -apply inbetween_step_not_Eq with (1 := Hx). -omega. -inversion_clear Hx as [|l' _ Hl]. -now rewrite (middle_odd _ Hk) in Hl. -Qed. - -Theorem inbetween_step_Lo_Mi_Eq_odd : - forall x k, - inbetween (start + Z2R k * step) (start + Z2R (k + 1) * step) x loc_Exact -> - (2 * k + 1 = nb_steps)%Z -> - inbetween start (start + Z2R nb_steps * step) x (loc_Inexact Lt). -Proof. -intros x k Hx Hk. -apply inbetween_step_not_Eq with (1 := Hx). -omega. -inversion_clear Hx as [Hl|]. -rewrite Hl. -rewrite Rcompare_plus_l, Rcompare_mult_r, Rcompare_half_r. -apply Rcompare_Lt. -change 2%R with (Z2R 2). -rewrite <- Z2R_mult. -apply Z2R_lt. -rewrite <- Hk. -apply Zlt_succ. -exact Hstep. -Qed. - -Theorem inbetween_step_Hi_Mi_even : - forall x k l, - inbetween (start + Z2R k * step) (start + Z2R (k + 1) * step) x l -> - l <> loc_Exact -> - (2 * k = nb_steps)%Z -> - inbetween start (start + Z2R nb_steps * step) x (loc_Inexact Gt). -Proof. -intros x k l Hx Hl Hk. -apply inbetween_step_not_Eq with (1 := Hx). -omega. -apply Rcompare_Gt. -assert (Hx' := inbetween_bounds_not_Eq _ _ _ _ Hx Hl). -apply Rle_lt_trans with (2 := proj1 Hx'). -apply Rcompare_not_Lt_inv. -rewrite Rcompare_plus_l, Rcompare_mult_r, Rcompare_half_r. -change 2%R with (Z2R 2). -rewrite <- Z2R_mult. -apply Rcompare_not_Lt. -apply Z2R_le. -rewrite Hk. -apply Zle_refl. -exact Hstep. -Qed. - -Theorem inbetween_step_Mi_Mi_even : - forall x k, - inbetween (start + Z2R k * step) (start + Z2R (k + 1) * step) x loc_Exact -> - (2 * k = nb_steps)%Z -> - inbetween start (start + Z2R nb_steps * step) x (loc_Inexact Eq). -Proof. -intros x k Hx Hk. -apply inbetween_step_not_Eq with (1 := Hx). -omega. -apply Rcompare_Eq. -inversion_clear Hx as [Hx'|]. -rewrite Hx', <- Hk, Z2R_mult. -simpl (Z2R 2). -field. -Qed. - -(** Computes a new location when the interval containing a real - number is split into nb_steps subintervals and the real is - in the k-th one. (Even radix.) *) - -Definition new_location_even k l := - if Zeq_bool k 0 then - match l with loc_Exact => l | _ => loc_Inexact Lt end - else - loc_Inexact - match Zcompare (2 * k) nb_steps with - | Lt => Lt - | Eq => match l with loc_Exact => Eq | _ => Gt end - | Gt => Gt - end. - -Theorem new_location_even_correct : - Zeven nb_steps = true -> - forall x k l, (0 <= k < nb_steps)%Z -> - inbetween (start + Z2R k * step) (start + Z2R (k + 1) * step) x l -> - inbetween start (start + Z2R nb_steps * step) x (new_location_even k l). -Proof. -intros He x k l Hk Hx. -unfold new_location_even. -destruct (Zeq_bool_spec k 0) as [Hk0|Hk0]. -(* k = 0 *) -rewrite Hk0 in Hx. -rewrite Rmult_0_l, Rplus_0_r, Rmult_1_l in Hx. -set (l' := match l with loc_Exact => l | _ => loc_Inexact Lt end). -assert ((l = loc_Exact /\ l' = loc_Exact) \/ (l <> loc_Exact /\ l' = loc_Inexact Lt)). -unfold l' ; case l ; try (now left) ; right ; now split. -destruct H as [(H1,H2)|(H1,H2)] ; rewrite H2. -constructor. -rewrite H1 in Hx. -now inversion_clear Hx. -now apply inbetween_step_Lo_not_Eq with (2 := H1). -(* k <> 0 *) -destruct (Zcompare_spec (2 * k) nb_steps) as [Hk1|Hk1|Hk1]. -(* . 2 * k < nb_steps *) -apply inbetween_step_Lo with (1 := Hx). -omega. -destruct (Zeven_ex nb_steps). -rewrite He in H. -omega. -(* . 2 * k = nb_steps *) -set (l' := match l with loc_Exact => Eq | _ => Gt end). -assert ((l = loc_Exact /\ l' = Eq) \/ (l <> loc_Exact /\ l' = Gt)). -unfold l' ; case l ; try (now left) ; right ; now split. -destruct H as [(H1,H2)|(H1,H2)] ; rewrite H2. -rewrite H1 in Hx. -now apply inbetween_step_Mi_Mi_even with (1 := Hx). -now apply inbetween_step_Hi_Mi_even with (1 := Hx). -(* . 2 * k > nb_steps *) -apply inbetween_step_Hi with (1 := Hx). -exact Hk1. -apply Hk. -Qed. - -(** Computes a new location when the interval containing a real - number is split into nb_steps subintervals and the real is - in the k-th one. (Odd radix.) *) - -Definition new_location_odd k l := - if Zeq_bool k 0 then - match l with loc_Exact => l | _ => loc_Inexact Lt end - else - loc_Inexact - match Zcompare (2 * k + 1) nb_steps with - | Lt => Lt - | Eq => match l with loc_Inexact l => l | loc_Exact => Lt end - | Gt => Gt - end. - -Theorem new_location_odd_correct : - Zeven nb_steps = false -> - forall x k l, (0 <= k < nb_steps)%Z -> - inbetween (start + Z2R k * step) (start + Z2R (k + 1) * step) x l -> - inbetween start (start + Z2R nb_steps * step) x (new_location_odd k l). -Proof. -intros Ho x k l Hk Hx. -unfold new_location_odd. -destruct (Zeq_bool_spec k 0) as [Hk0|Hk0]. -(* k = 0 *) -rewrite Hk0 in Hx. -rewrite Rmult_0_l, Rplus_0_r, Rmult_1_l in Hx. -set (l' := match l with loc_Exact => l | _ => loc_Inexact Lt end). -assert ((l = loc_Exact /\ l' = loc_Exact) \/ (l <> loc_Exact /\ l' = loc_Inexact Lt)). -unfold l' ; case l ; try (now left) ; right ; now split. -destruct H as [(H1,H2)|(H1,H2)] ; rewrite H2. -constructor. -rewrite H1 in Hx. -now inversion_clear Hx. -now apply inbetween_step_Lo_not_Eq with (2 := H1). -(* k <> 0 *) -destruct (Zcompare_spec (2 * k + 1) nb_steps) as [Hk1|Hk1|Hk1]. -(* . 2 * k + 1 < nb_steps *) -apply inbetween_step_Lo with (1 := Hx) (3 := Hk1). -omega. -(* . 2 * k + 1 = nb_steps *) -destruct l. -apply inbetween_step_Lo_Mi_Eq_odd with (1 := Hx) (2 := Hk1). -apply inbetween_step_any_Mi_odd with (1 := Hx) (2 := Hk1). -(* . 2 * k + 1 > nb_steps *) -apply inbetween_step_Hi with (1 := Hx). -destruct (Zeven_ex nb_steps). -rewrite Ho in H. -omega. -apply Hk. -Qed. - -Definition new_location := - if Zeven nb_steps then new_location_even else new_location_odd. - -Theorem new_location_correct : - forall x k l, (0 <= k < nb_steps)%Z -> - inbetween (start + Z2R k * step) (start + Z2R (k + 1) * step) x l -> - inbetween start (start + Z2R nb_steps * step) x (new_location k l). -Proof. -intros x k l Hk Hx. -unfold new_location. -generalize (refl_equal nb_steps) (Zle_lt_trans _ _ _ (proj1 Hk) (proj2 Hk)). -pattern nb_steps at 2 3 5 ; case nb_steps. -now intros _. -(* . *) -intros [p|p|] Hp _. -apply new_location_odd_correct with (2 := Hk) (3 := Hx). -now rewrite Hp. -apply new_location_even_correct with (2 := Hk) (3 := Hx). -now rewrite Hp. -now rewrite Hp in Hnb_steps. -(* . *) -now intros p _. -Qed. - -End Fcalc_bracket_step. - -Section Fcalc_bracket_scale. - -Lemma inbetween_mult_aux : - forall x d s, - ((x * s + d * s) / 2 = (x + d) / 2 * s)%R. -Proof. -intros x d s. -field. -Qed. - -Theorem inbetween_mult_compat : - forall x d u l s, - (0 < s)%R -> - inbetween x d u l -> - inbetween (x * s) (d * s) (u * s) l. -Proof. -intros x d u l s Hs [Hx|l' Hx Hl] ; constructor. -now rewrite Hx. -now split ; apply Rmult_lt_compat_r. -rewrite inbetween_mult_aux. -now rewrite Rcompare_mult_r. -Qed. - -Theorem inbetween_mult_reg : - forall x d u l s, - (0 < s)%R -> - inbetween (x * s) (d * s) (u * s) l -> - inbetween x d u l. -Proof. -intros x d u l s Hs [Hx|l' Hx Hl] ; constructor. -apply Rmult_eq_reg_r with (1 := Hx). -now apply Rgt_not_eq. -now split ; apply Rmult_lt_reg_r with s. -rewrite <- Rcompare_mult_r with (1 := Hs). -now rewrite inbetween_mult_aux in Hl. -Qed. - -End Fcalc_bracket_scale. - -Section Fcalc_bracket_generic. - -Variable beta : radix. -Notation bpow e := (bpow beta e). - -(** Specialization of inbetween for two consecutive floating-point numbers. *) - -Definition inbetween_float m e x l := - inbetween (F2R (Float beta m e)) (F2R (Float beta (m + 1) e)) x l. - -Theorem inbetween_float_bounds : - forall x m e l, - inbetween_float m e x l -> - (F2R (Float beta m e) <= x < F2R (Float beta (m + 1) e))%R. -Proof. -intros x m e l [Hx|l' Hx Hl]. -rewrite Hx. -split. -apply Rle_refl. -apply F2R_lt_compat. -apply Zlt_succ. -split. -now apply Rlt_le. -apply Hx. -Qed. - -(** Specialization of inbetween for two consecutive integers. *) - -Definition inbetween_int m x l := - inbetween (Z2R m) (Z2R (m + 1)) x l. - -Theorem inbetween_float_new_location : - forall x m e l k, - (0 < k)%Z -> - inbetween_float m e x l -> - inbetween_float (Zdiv m (Zpower beta k)) (e + k) x (new_location (Zpower beta k) (Zmod m (Zpower beta k)) l). -Proof. -intros x m e l k Hk Hx. -unfold inbetween_float in *. -assert (Hr: forall m, F2R (Float beta m (e + k)) = F2R (Float beta (m * Zpower beta k) e)). -clear -Hk. intros m. -rewrite (F2R_change_exp beta e). -apply (f_equal (fun r => F2R (Float beta (m * Zpower _ r) e))). -ring. -omega. -assert (Hp: (Zpower beta k > 0)%Z). -apply Zlt_gt. -apply Zpower_gt_0. -now apply Zlt_le_weak. -(* . *) -rewrite 2!Hr. -rewrite Zmult_plus_distr_l, Zmult_1_l. -unfold F2R at 2. simpl. -rewrite Z2R_plus, Rmult_plus_distr_r. -apply new_location_correct. -apply bpow_gt_0. -now apply Zpower_gt_1. -now apply Z_mod_lt. -rewrite <- 2!Rmult_plus_distr_r, <- 2!Z2R_plus. -rewrite Zmult_comm, Zplus_assoc. -now rewrite <- Z_div_mod_eq. -Qed. - -Theorem inbetween_float_new_location_single : - forall x m e l, - inbetween_float m e x l -> - inbetween_float (Zdiv m beta) (e + 1) x (new_location beta (Zmod m beta) l). -Proof. -intros x m e l Hx. -replace (radix_val beta) with (Zpower beta 1). -now apply inbetween_float_new_location. -apply Zmult_1_r. -Qed. - -Theorem inbetween_float_ex : - forall m e l, - exists x, - inbetween_float m e x l. -Proof. -intros m e l. -apply inbetween_ex. -apply F2R_lt_compat. -apply Zlt_succ. -Qed. - -Theorem inbetween_float_unique : - forall x e m l m' l', - inbetween_float m e x l -> - inbetween_float m' e x l' -> - m = m' /\ l = l'. -Proof. -intros x e m l m' l' H H'. -refine ((fun Hm => conj Hm _) _). -rewrite <- Hm in H'. clear -H H'. -apply inbetween_unique with (1 := H) (2 := H'). -destruct (inbetween_float_bounds x m e l H) as (H1,H2). -destruct (inbetween_float_bounds x m' e l' H') as (H3,H4). -cut (m < m' + 1 /\ m' < m + 1)%Z. clear ; omega. -now split ; apply F2R_lt_reg with beta e ; apply Rle_lt_trans with x. -Qed. - -End Fcalc_bracket_generic. diff --git a/flocq/Calc/Fcalc_digits.v b/flocq/Calc/Fcalc_digits.v deleted file mode 100644 index 45133e81..00000000 --- a/flocq/Calc/Fcalc_digits.v +++ /dev/null @@ -1,63 +0,0 @@ -(** -This file is part of the Flocq formalization of floating-point -arithmetic in Coq: http://flocq.gforge.inria.fr/ - -Copyright (C) 2010-2013 Sylvie Boldo -#
# -Copyright (C) 2010-2013 Guillaume Melquiond - -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. -*) - -(** * Functions for computing the number of digits of integers and related theorems. *) - -Require Import Fcore_Raux. -Require Import Fcore_defs. -Require Import Fcore_float_prop. -Require Import Fcore_digits. - -Section Fcalc_digits. - -Variable beta : radix. -Notation bpow e := (bpow beta e). - -Theorem Zdigits_ln_beta : - forall n, - n <> Z0 -> - Zdigits beta n = ln_beta beta (Z2R n). -Proof. -intros n Hn. -destruct (ln_beta beta (Z2R n)) as (e, He) ; simpl. -specialize (He (Z2R_neq _ _ Hn)). -rewrite <- Z2R_abs in He. -assert (Hd := Zdigits_correct beta n). -assert (Hd' := Zdigits_gt_0 beta n). -apply Zle_antisym ; apply (bpow_lt_bpow beta). -apply Rle_lt_trans with (2 := proj2 He). -rewrite <- Z2R_Zpower by omega. -now apply Z2R_le. -apply Rle_lt_trans with (1 := proj1 He). -rewrite <- Z2R_Zpower by omega. -now apply Z2R_lt. -Qed. - -Theorem ln_beta_F2R_Zdigits : - forall m e, m <> Z0 -> - (ln_beta beta (F2R (Float beta m e)) = Zdigits beta m + e :> Z)%Z. -Proof. -intros m e Hm. -rewrite ln_beta_F2R with (1 := Hm). -apply (f_equal (fun v => Zplus v e)). -apply sym_eq. -now apply Zdigits_ln_beta. -Qed. - -End Fcalc_digits. diff --git a/flocq/Calc/Fcalc_div.v b/flocq/Calc/Fcalc_div.v deleted file mode 100644 index c8f1f9fc..00000000 --- a/flocq/Calc/Fcalc_div.v +++ /dev/null @@ -1,165 +0,0 @@ -(** -This file is part of the Flocq formalization of floating-point -arithmetic in Coq: http://flocq.gforge.inria.fr/ - -Copyright (C) 2010-2013 Sylvie Boldo -#
# -Copyright (C) 2010-2013 Guillaume Melquiond - -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. -*) - -(** * Helper function and theorem for computing the rounded quotient of two floating-point numbers. *) - -Require Import Fcore_Raux. -Require Import Fcore_defs. -Require Import Fcore_float_prop. -Require Import Fcore_digits. -Require Import Fcalc_bracket. -Require Import Fcalc_digits. - -Section Fcalc_div. - -Variable beta : radix. -Notation bpow e := (bpow beta e). - -(** Computes a mantissa of precision p, the corresponding exponent, - and the position with respect to the real quotient of the - input floating-point numbers. - -The algorithm performs the following steps: -- Shift dividend mantissa so that it has at least p2 + p digits. -- Perform the Euclidean division. -- Compute the position according to the division remainder. - -Complexity is fine as long as p1 <= 2p and p2 <= p. -*) - -Definition Fdiv_core prec m1 e1 m2 e2 := - let d1 := Zdigits beta m1 in - let d2 := Zdigits beta m2 in - let e := (e1 - e2)%Z in - let (m, e') := - match (d2 + prec - d1)%Z with - | Zpos p => (m1 * Zpower_pos beta p, e + Zneg p)%Z - | _ => (m1, e) - end in - let '(q, r) := Zdiv_eucl m m2 in - (q, e', new_location m2 r loc_Exact). - -Theorem Fdiv_core_correct : - forall prec m1 e1 m2 e2, - (0 < prec)%Z -> - (0 < m1)%Z -> (0 < m2)%Z -> - let '(m, e, l) := Fdiv_core prec m1 e1 m2 e2 in - (prec <= Zdigits beta m)%Z /\ - inbetween_float beta m e (F2R (Float beta m1 e1) / F2R (Float beta m2 e2)) l. -Proof. -intros prec m1 e1 m2 e2 Hprec Hm1 Hm2. -unfold Fdiv_core. -set (d1 := Zdigits beta m1). -set (d2 := Zdigits beta m2). -case_eq - (match (d2 + prec - d1)%Z with - | Zpos p => ((m1 * Zpower_pos beta p)%Z, (e1 - e2 + Zneg p)%Z) - | _ => (m1, (e1 - e2)%Z) - end). -intros m' e' Hme. -(* . the shifted mantissa m' has enough digits *) -assert (Hs: F2R (Float beta m' (e' + e2)) = F2R (Float beta m1 e1) /\ (0 < m')%Z /\ (d2 + prec <= Zdigits beta m')%Z). -replace (d2 + prec)%Z with (d2 + prec - d1 + d1)%Z by ring. -destruct (d2 + prec - d1)%Z as [|p|p] ; - unfold d1 ; - inversion Hme. -ring_simplify (e1 - e2 + e2)%Z. -repeat split. -now rewrite <- H0. -apply Zle_refl. -replace (e1 - e2 + Zneg p + e2)%Z with (e1 - Zpos p)%Z by (unfold Zminus ; simpl ; ring). -fold (Zpower beta (Zpos p)). -split. -pattern (Zpos p) at 1 ; replace (Zpos p) with (e1 - (e1 - Zpos p))%Z by ring. -apply sym_eq. -apply F2R_change_exp. -assert (0 < Zpos p)%Z by easy. -omega. -split. -apply Zmult_lt_0_compat. -exact Hm1. -now apply Zpower_gt_0. -rewrite Zdigits_mult_Zpower. -rewrite Zplus_comm. -apply Zle_refl. -apply sym_not_eq. -now apply Zlt_not_eq. -easy. -split. -now ring_simplify (e1 - e2 + e2)%Z. -assert (Zneg p < 0)%Z by easy. -omega. -(* . *) -destruct Hs as (Hs1, (Hs2, Hs3)). -rewrite <- Hs1. -generalize (Z_div_mod m' m2 (Zlt_gt _ _ Hm2)). -destruct (Zdiv_eucl m' m2) as (q, r). -intros (Hq, Hr). -split. -(* . the result mantissa q has enough digits *) -cut (Zdigits beta m' <= d2 + Zdigits beta q)%Z. omega. -unfold d2. -rewrite Hq. -assert (Hq': (0 < q)%Z). -apply Zmult_lt_reg_r with (1 := Hm2). -assert (m2 < m')%Z. -apply lt_Zdigits with beta. -now apply Zlt_le_weak. -unfold d2 in Hs3. -clear -Hprec Hs3 ; omega. -cut (q * m2 = m' - r)%Z. clear -Hr H ; omega. -rewrite Hq. -ring. -apply Zle_trans with (Zdigits beta (m2 + q + m2 * q)). -apply Zdigits_le. -rewrite <- Hq. -now apply Zlt_le_weak. -clear -Hr Hq'. omega. -apply Zdigits_mult_strong ; apply Zlt_le_weak. -now apply Zle_lt_trans with r. -exact Hq'. -(* . the location is correctly computed *) -unfold inbetween_float, F2R. simpl. -rewrite bpow_plus, Z2R_plus. -rewrite Hq, Z2R_plus, Z2R_mult. -replace ((Z2R m2 * Z2R q + Z2R r) * (bpow e' * bpow e2) / (Z2R m2 * bpow e2))%R - with ((Z2R q + Z2R r / Z2R m2) * bpow e')%R. -apply inbetween_mult_compat. -apply bpow_gt_0. -destruct (Z_lt_le_dec 1 m2) as [Hm2''|Hm2'']. -replace (Z2R 1) with (Z2R m2 * /Z2R m2)%R. -apply new_location_correct ; try easy. -apply Rinv_0_lt_compat. -now apply (Z2R_lt 0). -now constructor. -apply Rinv_r. -apply Rgt_not_eq. -now apply (Z2R_lt 0). -assert (r = 0 /\ m2 = 1)%Z by (clear -Hr Hm2'' ; omega). -rewrite (proj1 H), (proj2 H). -unfold Rdiv. -rewrite Rmult_0_l, Rplus_0_r. -now constructor. -field. -split ; apply Rgt_not_eq. -apply bpow_gt_0. -now apply (Z2R_lt 0). -Qed. - -End Fcalc_div. diff --git a/flocq/Calc/Fcalc_ops.v b/flocq/Calc/Fcalc_ops.v deleted file mode 100644 index e834c044..00000000 --- a/flocq/Calc/Fcalc_ops.v +++ /dev/null @@ -1,163 +0,0 @@ -(** -This file is part of the Flocq formalization of floating-point -arithmetic in Coq: http://flocq.gforge.inria.fr/ - -Copyright (C) 2010-2013 Sylvie Boldo -#
# -Copyright (C) 2010-2013 Guillaume Melquiond - -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. -*) - -(** Basic operations on floats: alignment, addition, multiplication *) -Require Import Fcore_Raux. -Require Import Fcore_defs. -Require Import Fcore_float_prop. - -Section Float_ops. - -Variable beta : radix. - -Notation bpow e := (bpow beta e). - -Arguments Float {beta} Fnum Fexp. - -Definition Falign (f1 f2 : float beta) := - let '(Float m1 e1) := f1 in - let '(Float m2 e2) := f2 in - if Zle_bool e1 e2 - then (m1, (m2 * Zpower beta (e2 - e1))%Z, e1) - else ((m1 * Zpower beta (e1 - e2))%Z, m2, e2). - -Theorem Falign_spec : - forall f1 f2 : float beta, - let '(m1, m2, e) := Falign f1 f2 in - F2R f1 = @F2R beta (Float m1 e) /\ F2R f2 = @F2R beta (Float m2 e). -Proof. -unfold Falign. -intros (m1, e1) (m2, e2). -generalize (Zle_cases e1 e2). -case (Zle_bool e1 e2) ; intros He ; split ; trivial. -now rewrite <- F2R_change_exp. -rewrite <- F2R_change_exp. -apply refl_equal. -omega. -Qed. - -Theorem Falign_spec_exp: - forall f1 f2 : float beta, - snd (Falign f1 f2) = Zmin (Fexp f1) (Fexp f2). -Proof. -intros (m1,e1) (m2,e2). -unfold Falign; simpl. -generalize (Zle_cases e1 e2);case (Zle_bool e1 e2); intros He. -case (Zmin_spec e1 e2); intros (H1,H2); easy. -case (Zmin_spec e1 e2); intros (H1,H2); easy. -Qed. - -Definition Fopp (f1 : float beta) : float beta := - let '(Float m1 e1) := f1 in - Float (-m1)%Z e1. - -Theorem F2R_opp : - forall f1 : float beta, - (F2R (Fopp f1) = -F2R f1)%R. -intros (m1,e1). -apply F2R_Zopp. -Qed. - -Definition Fabs (f1 : float beta) : float beta := - let '(Float m1 e1) := f1 in - Float (Zabs m1)%Z e1. - -Theorem F2R_abs : - forall f1 : float beta, - (F2R (Fabs f1) = Rabs (F2R f1))%R. -intros (m1,e1). -apply F2R_Zabs. -Qed. - -Definition Fplus (f1 f2 : float beta) : float beta := - let '(m1, m2 ,e) := Falign f1 f2 in - Float (m1 + m2) e. - -Theorem F2R_plus : - forall f1 f2 : float beta, - F2R (Fplus f1 f2) = (F2R f1 + F2R f2)%R. -Proof. -intros f1 f2. -unfold Fplus. -generalize (Falign_spec f1 f2). -destruct (Falign f1 f2) as ((m1, m2), e). -intros (H1, H2). -rewrite H1, H2. -unfold F2R. simpl. -rewrite Z2R_plus. -apply Rmult_plus_distr_r. -Qed. - -Theorem Fplus_same_exp : - forall m1 m2 e, - Fplus (Float m1 e) (Float m2 e) = Float (m1 + m2) e. -Proof. -intros m1 m2 e. -unfold Fplus. -simpl. -now rewrite Zle_bool_refl, Zminus_diag, Zmult_1_r. -Qed. - -Theorem Fexp_Fplus : - forall f1 f2 : float beta, - Fexp (Fplus f1 f2) = Zmin (Fexp f1) (Fexp f2). -Proof. -intros f1 f2. -unfold Fplus. -rewrite <- Falign_spec_exp. -now destruct (Falign f1 f2) as ((p,q),e). -Qed. - -Definition Fminus (f1 f2 : float beta) := - Fplus f1 (Fopp f2). - -Theorem F2R_minus : - forall f1 f2 : float beta, - F2R (Fminus f1 f2) = (F2R f1 - F2R f2)%R. -Proof. -intros f1 f2; unfold Fminus. -rewrite F2R_plus, F2R_opp. -ring. -Qed. - -Theorem Fminus_same_exp : - forall m1 m2 e, - Fminus (Float m1 e) (Float m2 e) = Float (m1 - m2) e. -Proof. -intros m1 m2 e. -unfold Fminus. -apply Fplus_same_exp. -Qed. - -Definition Fmult (f1 f2 : float beta) : float beta := - let '(Float m1 e1) := f1 in - let '(Float m2 e2) := f2 in - Float (m1 * m2) (e1 + e2). - -Theorem F2R_mult : - forall f1 f2 : float beta, - F2R (Fmult f1 f2) = (F2R f1 * F2R f2)%R. -Proof. -intros (m1, e1) (m2, e2). -unfold Fmult, F2R. simpl. -rewrite Z2R_mult, bpow_plus. -ring. -Qed. - -End Float_ops. diff --git a/flocq/Calc/Fcalc_round.v b/flocq/Calc/Fcalc_round.v deleted file mode 100644 index 86422247..00000000 --- a/flocq/Calc/Fcalc_round.v +++ /dev/null @@ -1,1094 +0,0 @@ -(** -This file is part of the Flocq formalization of floating-point -arithmetic in Coq: http://flocq.gforge.inria.fr/ - -Copyright (C) 2010-2013 Sylvie Boldo -#
# -Copyright (C) 2010-2013 Guillaume Melquiond - -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. -*) - -(** * Helper function for computing the rounded value of a real number. *) - -Require Import Fcore. -Require Import Fcore_digits. -Require Import Fcalc_bracket. -Require Import Fcalc_digits. - -Section Fcalc_round. - -Variable beta : radix. -Notation bpow e := (bpow beta e). - -Section Fcalc_round_fexp. - -Variable fexp : Z -> Z. -Context { valid_exp : Valid_exp fexp }. -Notation format := (generic_format beta fexp). - -(** Relates location and rounding. *) - -Theorem inbetween_float_round : - forall rnd choice, - ( forall x m l, inbetween_int m x l -> rnd x = choice m l ) -> - forall x m l, - let e := canonic_exp beta fexp x in - inbetween_float beta m e x l -> - round beta fexp rnd x = F2R (Float beta (choice m l) e). -Proof. -intros rnd choice Hc x m l e Hl. -unfold round, F2R. simpl. -apply (f_equal (fun m => (Z2R m * bpow e)%R)). -apply Hc. -apply inbetween_mult_reg with (bpow e). -apply bpow_gt_0. -now rewrite scaled_mantissa_mult_bpow. -Qed. - -Definition cond_incr (b : bool) m := if b then (m + 1)%Z else m. - -Theorem inbetween_float_round_sign : - forall rnd choice, - ( forall x m l, inbetween_int m (Rabs x) l -> - rnd x = cond_Zopp (Rlt_bool x 0) (choice (Rlt_bool x 0) m l) ) -> - forall x m l, - let e := canonic_exp beta fexp x in - inbetween_float beta m e (Rabs x) l -> - round beta fexp rnd x = F2R (Float beta (cond_Zopp (Rlt_bool x 0) (choice (Rlt_bool x 0) m l)) e). -Proof. -intros rnd choice Hc x m l e Hx. -apply (f_equal (fun m => (Z2R m * bpow e)%R)). -simpl. -replace (Rlt_bool x 0) with (Rlt_bool (scaled_mantissa beta fexp x) 0). -(* *) -apply Hc. -apply inbetween_mult_reg with (bpow e). -apply bpow_gt_0. -rewrite <- (Rabs_right (bpow e)) at 3. -rewrite <- Rabs_mult. -now rewrite scaled_mantissa_mult_bpow. -apply Rle_ge. -apply bpow_ge_0. -(* *) -destruct (Rlt_bool_spec x 0) as [Zx|Zx] ; simpl. -apply Rlt_bool_true. -rewrite <- (Rmult_0_l (bpow (-e))). -apply Rmult_lt_compat_r with (2 := Zx). -apply bpow_gt_0. -apply Rlt_bool_false. -apply Rmult_le_pos with (1 := Zx). -apply bpow_ge_0. -Qed. - -(** Relates location and rounding down. *) - -Theorem inbetween_int_DN : - forall x m l, - inbetween_int m x l -> - Zfloor x = m. -Proof. -intros x m l Hl. -refine (Zfloor_imp m _ _). -apply inbetween_bounds with (2 := Hl). -apply Z2R_lt. -apply Zlt_succ. -Qed. - -Theorem inbetween_float_DN : - forall x m l, - let e := canonic_exp beta fexp x in - inbetween_float beta m e x l -> - round beta fexp Zfloor x = F2R (Float beta m e). -Proof. -apply inbetween_float_round with (choice := fun m l => m). -exact inbetween_int_DN. -Qed. - -Definition round_sign_DN s l := - match l with - | loc_Exact => false - | _ => s - end. - -Theorem inbetween_int_DN_sign : - forall x m l, - inbetween_int m (Rabs x) l -> - Zfloor x = cond_Zopp (Rlt_bool x 0) (cond_incr (round_sign_DN (Rlt_bool x 0) l) m). -Proof. -intros x m l Hl. -unfold Rabs in Hl. -destruct (Rcase_abs x) as [Zx|Zx] . -(* *) -rewrite Rlt_bool_true with (1 := Zx). -inversion_clear Hl ; simpl. -rewrite <- (Ropp_involutive x). -rewrite H, <- Z2R_opp. -apply Zfloor_Z2R. -apply Zfloor_imp. -split. -apply Rlt_le. -rewrite Z2R_opp. -apply Ropp_lt_cancel. -now rewrite Ropp_involutive. -ring_simplify (- (m + 1) + 1)%Z. -rewrite Z2R_opp. -apply Ropp_lt_cancel. -now rewrite Ropp_involutive. -(* *) -rewrite Rlt_bool_false. -inversion_clear Hl ; simpl. -rewrite H. -apply Zfloor_Z2R. -apply Zfloor_imp. -split. -now apply Rlt_le. -apply H. -now apply Rge_le. -Qed. - -Theorem inbetween_float_DN_sign : - forall x m l, - let e := canonic_exp beta fexp x in - inbetween_float beta m e (Rabs x) l -> - round beta fexp Zfloor x = F2R (Float beta (cond_Zopp (Rlt_bool x 0) (cond_incr (round_sign_DN (Rlt_bool x 0) l) m)) e). -Proof. -apply inbetween_float_round_sign with (choice := fun s m l => cond_incr (round_sign_DN s l) m). -exact inbetween_int_DN_sign. -Qed. - -(** Relates location and rounding up. *) - -Definition round_UP l := - match l with - | loc_Exact => false - | _ => true - end. - -Theorem inbetween_int_UP : - forall x m l, - inbetween_int m x l -> - Zceil x = cond_incr (round_UP l) m. -Proof. -intros x m l Hl. -assert (Hl': l = loc_Exact \/ (l <> loc_Exact /\ round_UP l = true)). -case l ; try (now left) ; now right ; split. -destruct Hl' as [Hl'|(Hl1, Hl2)]. -(* Exact *) -rewrite Hl'. -destruct Hl ; try easy. -rewrite H. -exact (Zceil_Z2R _). -(* not Exact *) -rewrite Hl2. -simpl. -apply Zceil_imp. -ring_simplify (m + 1 - 1)%Z. -refine (let H := _ in conj (proj1 H) (Rlt_le _ _ (proj2 H))). -apply inbetween_bounds_not_Eq with (2 := Hl1) (1 := Hl). -Qed. - -Theorem inbetween_float_UP : - forall x m l, - let e := canonic_exp beta fexp x in - inbetween_float beta m e x l -> - round beta fexp Zceil x = F2R (Float beta (cond_incr (round_UP l) m) e). -Proof. -apply inbetween_float_round with (choice := fun m l => cond_incr (round_UP l) m). -exact inbetween_int_UP. -Qed. - -Definition round_sign_UP s l := - match l with - | loc_Exact => false - | _ => negb s - end. - -Theorem inbetween_int_UP_sign : - forall x m l, - inbetween_int m (Rabs x) l -> - Zceil x = cond_Zopp (Rlt_bool x 0) (cond_incr (round_sign_UP (Rlt_bool x 0) l) m). -Proof. -intros x m l Hl. -unfold Rabs in Hl. -destruct (Rcase_abs x) as [Zx|Zx] . -(* *) -rewrite Rlt_bool_true with (1 := Zx). -simpl. -unfold Zceil. -apply f_equal. -inversion_clear Hl ; simpl. -rewrite H. -apply Zfloor_Z2R. -apply Zfloor_imp. -split. -now apply Rlt_le. -apply H. -(* *) -rewrite Rlt_bool_false. -simpl. -inversion_clear Hl ; simpl. -rewrite H. -apply Zceil_Z2R. -apply Zceil_imp. -split. -change (m + 1 - 1)%Z with (Zpred (Zsucc m)). -now rewrite <- Zpred_succ. -now apply Rlt_le. -now apply Rge_le. -Qed. - -Theorem inbetween_float_UP_sign : - forall x m l, - let e := canonic_exp beta fexp x in - inbetween_float beta m e (Rabs x) l -> - round beta fexp Zceil x = F2R (Float beta (cond_Zopp (Rlt_bool x 0) (cond_incr (round_sign_UP (Rlt_bool x 0) l) m)) e). -Proof. -apply inbetween_float_round_sign with (choice := fun s m l => cond_incr (round_sign_UP s l) m). -exact inbetween_int_UP_sign. -Qed. - -(** Relates location and rounding toward zero. *) - -Definition round_ZR (s : bool) l := - match l with - | loc_Exact => false - | _ => s - end. - -Theorem inbetween_int_ZR : - forall x m l, - inbetween_int m x l -> - Ztrunc x = cond_incr (round_ZR (Zlt_bool m 0) l) m. -Proof with auto with typeclass_instances. -intros x m l Hl. -inversion_clear Hl as [Hx|l' Hx Hl']. -(* Exact *) -rewrite Hx. -rewrite Zrnd_Z2R... -(* not Exact *) -unfold Ztrunc. -assert (Hm: Zfloor x = m). -apply Zfloor_imp. -exact (conj (Rlt_le _ _ (proj1 Hx)) (proj2 Hx)). -rewrite Zceil_floor_neq. -rewrite Hm. -unfold cond_incr. -simpl. -case Rlt_bool_spec ; intros Hx' ; - case Zlt_bool_spec ; intros Hm' ; try apply refl_equal. -elim Rlt_not_le with (1 := Hx'). -apply Rlt_le. -apply Rle_lt_trans with (2 := proj1 Hx). -now apply (Z2R_le 0). -elim Rle_not_lt with (1 := Hx'). -apply Rlt_le_trans with (1 := proj2 Hx). -apply (Z2R_le _ 0). -now apply Zlt_le_succ. -rewrite Hm. -now apply Rlt_not_eq. -Qed. - -Theorem inbetween_float_ZR : - forall x m l, - let e := canonic_exp beta fexp x in - inbetween_float beta m e x l -> - round beta fexp Ztrunc x = F2R (Float beta (cond_incr (round_ZR (Zlt_bool m 0) l) m) e). -Proof. -apply inbetween_float_round with (choice := fun m l => cond_incr (round_ZR (Zlt_bool m 0) l) m). -exact inbetween_int_ZR. -Qed. - -Theorem inbetween_int_ZR_sign : - forall x m l, - inbetween_int m (Rabs x) l -> - Ztrunc x = cond_Zopp (Rlt_bool x 0) m. -Proof. -intros x m l Hl. -simpl. -unfold Ztrunc. -destruct (Rlt_le_dec x 0) as [Zx|Zx]. -(* *) -rewrite Rlt_bool_true with (1 := Zx). -simpl. -unfold Zceil. -apply f_equal. -apply Zfloor_imp. -rewrite <- Rabs_left with (1 := Zx). -apply inbetween_bounds with (2 := Hl). -apply Z2R_lt. -apply Zlt_succ. -(* *) -rewrite Rlt_bool_false with (1 := Zx). -simpl. -apply Zfloor_imp. -rewrite <- Rabs_pos_eq with (1 := Zx). -apply inbetween_bounds with (2 := Hl). -apply Z2R_lt. -apply Zlt_succ. -Qed. - -Theorem inbetween_float_ZR_sign : - forall x m l, - let e := canonic_exp beta fexp x in - inbetween_float beta m e (Rabs x) l -> - round beta fexp Ztrunc x = F2R (Float beta (cond_Zopp (Rlt_bool x 0) m) e). -Proof. -apply inbetween_float_round_sign with (choice := fun s m l => m). -exact inbetween_int_ZR_sign. -Qed. - -(** Relates location and rounding to nearest. *) - -Definition round_N (p : bool) l := - match l with - | loc_Exact => false - | loc_Inexact Lt => false - | loc_Inexact Eq => p - | loc_Inexact Gt => true - end. - -Theorem inbetween_int_N : - forall choice x m l, - inbetween_int m x l -> - Znearest choice x = cond_incr (round_N (choice m) l) m. -Proof with auto with typeclass_instances. -intros choice x m l Hl. -inversion_clear Hl as [Hx|l' Hx Hl']. -(* Exact *) -rewrite Hx. -rewrite Zrnd_Z2R... -(* not Exact *) -unfold Znearest. -assert (Hm: Zfloor x = m). -apply Zfloor_imp. -exact (conj (Rlt_le _ _ (proj1 Hx)) (proj2 Hx)). -rewrite Zceil_floor_neq. -rewrite Hm. -replace (Rcompare (x - Z2R m) (/2)) with l'. -now case l'. -rewrite <- Hl'. -rewrite Z2R_plus. -rewrite <- (Rcompare_plus_r (- Z2R m) x). -apply f_equal. -simpl (Z2R 1). -field. -rewrite Hm. -now apply Rlt_not_eq. -Qed. - -Theorem inbetween_int_N_sign : - forall choice x m l, - inbetween_int m (Rabs x) l -> - Znearest choice x = cond_Zopp (Rlt_bool x 0) (cond_incr (round_N (if Rlt_bool x 0 then negb (choice (-(m + 1))%Z) else choice m) l) m). -Proof with auto with typeclass_instances. -intros choice x m l Hl. -simpl. -unfold Rabs in Hl. -destruct (Rcase_abs x) as [Zx|Zx]. -(* *) -rewrite Rlt_bool_true with (1 := Zx). -simpl. -rewrite <- (Ropp_involutive x). -rewrite Znearest_opp. -apply f_equal. -inversion_clear Hl as [Hx|l' Hx Hl']. -rewrite Hx. -apply Zrnd_Z2R... -assert (Hm: Zfloor (-x) = m). -apply Zfloor_imp. -exact (conj (Rlt_le _ _ (proj1 Hx)) (proj2 Hx)). -unfold Znearest. -rewrite Zceil_floor_neq. -rewrite Hm. -replace (Rcompare (- x - Z2R m) (/2)) with l'. -now case l'. -rewrite <- Hl'. -rewrite Z2R_plus. -rewrite <- (Rcompare_plus_r (- Z2R m) (-x)). -apply f_equal. -simpl (Z2R 1). -field. -rewrite Hm. -now apply Rlt_not_eq. -(* *) -generalize (Rge_le _ _ Zx). -clear Zx. intros Zx. -rewrite Rlt_bool_false with (1 := Zx). -simpl. -inversion_clear Hl as [Hx|l' Hx Hl']. -rewrite Hx. -apply Zrnd_Z2R... -assert (Hm: Zfloor x = m). -apply Zfloor_imp. -exact (conj (Rlt_le _ _ (proj1 Hx)) (proj2 Hx)). -unfold Znearest. -rewrite Zceil_floor_neq. -rewrite Hm. -replace (Rcompare (x - Z2R m) (/2)) with l'. -now case l'. -rewrite <- Hl'. -rewrite Z2R_plus. -rewrite <- (Rcompare_plus_r (- Z2R m) x). -apply f_equal. -simpl (Z2R 1). -field. -rewrite Hm. -now apply Rlt_not_eq. -Qed. - -(** Relates location and rounding to nearest even. *) - -Theorem inbetween_int_NE : - forall x m l, - inbetween_int m x l -> - ZnearestE x = cond_incr (round_N (negb (Zeven m)) l) m. -Proof. -intros x m l Hl. -now apply inbetween_int_N with (choice := fun x => negb (Zeven x)). -Qed. - -Theorem inbetween_float_NE : - forall x m l, - let e := canonic_exp beta fexp x in - inbetween_float beta m e x l -> - round beta fexp ZnearestE x = F2R (Float beta (cond_incr (round_N (negb (Zeven m)) l) m) e). -Proof. -apply inbetween_float_round with (choice := fun m l => cond_incr (round_N (negb (Zeven m)) l) m). -exact inbetween_int_NE. -Qed. - -Theorem inbetween_int_NE_sign : - forall x m l, - inbetween_int m (Rabs x) l -> - ZnearestE x = cond_Zopp (Rlt_bool x 0) (cond_incr (round_N (negb (Zeven m)) l) m). -Proof. -intros x m l Hl. -erewrite inbetween_int_N_sign with (choice := fun x => negb (Zeven x)). -2: eexact Hl. -apply f_equal. -case Rlt_bool. -rewrite Zeven_opp, Zeven_plus. -now case (Zeven m). -apply refl_equal. -Qed. - -Theorem inbetween_float_NE_sign : - forall x m l, - let e := canonic_exp beta fexp x in - inbetween_float beta m e (Rabs x) l -> - round beta fexp ZnearestE x = F2R (Float beta (cond_Zopp (Rlt_bool x 0) (cond_incr (round_N (negb (Zeven m)) l) m)) e). -Proof. -apply inbetween_float_round_sign with (choice := fun s m l => cond_incr (round_N (negb (Zeven m)) l) m). -exact inbetween_int_NE_sign. -Qed. - -(** Relates location and rounding to nearest away. *) - -Theorem inbetween_int_NA : - forall x m l, - inbetween_int m x l -> - ZnearestA x = cond_incr (round_N (Zle_bool 0 m) l) m. -Proof. -intros x m l Hl. -now apply inbetween_int_N with (choice := fun x => Zle_bool 0 x). -Qed. - -Theorem inbetween_float_NA : - forall x m l, - let e := canonic_exp beta fexp x in - inbetween_float beta m e x l -> - round beta fexp ZnearestA x = F2R (Float beta (cond_incr (round_N (Zle_bool 0 m) l) m) e). -Proof. -apply inbetween_float_round with (choice := fun m l => cond_incr (round_N (Zle_bool 0 m) l) m). -exact inbetween_int_NA. -Qed. - -Theorem inbetween_int_NA_sign : - forall x m l, - inbetween_int m (Rabs x) l -> - ZnearestA x = cond_Zopp (Rlt_bool x 0) (cond_incr (round_N true l) m). -Proof. -intros x m l Hl. -erewrite inbetween_int_N_sign with (choice := Zle_bool 0). -2: eexact Hl. -apply f_equal. -assert (Hm: (0 <= m)%Z). -apply Zlt_succ_le. -apply lt_Z2R. -apply Rle_lt_trans with (Rabs x). -apply Rabs_pos. -refine (proj2 (inbetween_bounds _ _ _ _ _ Hl)). -apply Z2R_lt. -apply Zlt_succ. -rewrite Zle_bool_true with (1 := Hm). -rewrite Zle_bool_false. -now case Rlt_bool. -omega. -Qed. - -Definition truncate_aux t k := - let '(m, e, l) := t in - let p := Zpower beta k in - (Zdiv m p, (e + k)%Z, new_location p (Zmod m p) l). - -Theorem truncate_aux_comp : - forall t k1 k2, - (0 < k1)%Z -> - (0 < k2)%Z -> - truncate_aux t (k1 + k2) = truncate_aux (truncate_aux t k1) k2. -Proof. -intros ((m,e),l) k1 k2 Hk1 Hk2. -unfold truncate_aux. -destruct (inbetween_float_ex beta m e l) as (x,Hx). -assert (B1 := inbetween_float_new_location _ _ _ _ _ _ Hk1 Hx). -assert (Hk3: (0 < k1 + k2)%Z). -change Z0 with (0 + 0)%Z. -now apply Zplus_lt_compat. -assert (B3 := inbetween_float_new_location _ _ _ _ _ _ Hk3 Hx). -assert (B2 := inbetween_float_new_location _ _ _ _ _ _ Hk2 B1). -rewrite Zplus_assoc in B3. -destruct (inbetween_float_unique _ _ _ _ _ _ _ B2 B3). -now rewrite H, H0, Zplus_assoc. -Qed. - -(** Given a triple (mantissa, exponent, position), this function - computes a triple with a canonic exponent, assuming the - original triple had enough precision. *) - -Definition truncate t := - let '(m, e, l) := t in - let k := (fexp (Zdigits beta m + e) - e)%Z in - if Zlt_bool 0 k then truncate_aux t k - else t. - -Theorem truncate_0 : - forall e l, - let '(m', e', l') := truncate (0, e, l)%Z in - m' = Z0. -Proof. -intros e l. -unfold truncate. -case Zlt_bool. -simpl. -apply Zdiv_0_l. -apply refl_equal. -Qed. - -Theorem generic_format_truncate : - forall m e l, - (0 <= m)%Z -> - let '(m', e', l') := truncate (m, e, l) in - format (F2R (Float beta m' e')). -Proof. -intros m e l Hm. -unfold truncate. -set (k := (fexp (Zdigits beta m + e) - e)%Z). -case Zlt_bool_spec ; intros Hk. -(* *) -unfold truncate_aux. -apply generic_format_F2R. -intros Hd. -unfold canonic_exp. -rewrite ln_beta_F2R_Zdigits with (1 := Hd). -rewrite Zdigits_div_Zpower with (1 := Hm). -replace (Zdigits beta m - k + (e + k))%Z with (Zdigits beta m + e)%Z by ring. -unfold k. -ring_simplify. -apply Zle_refl. -split. -now apply Zlt_le_weak. -apply Znot_gt_le. -contradict Hd. -apply Zdiv_small. -apply conj with (1 := Hm). -rewrite <- Zabs_eq with (1 := Hm). -apply Zpower_gt_Zdigits. -apply Zlt_le_weak. -now apply Zgt_lt. -(* *) -destruct (Zle_lt_or_eq _ _ Hm) as [Hm'|Hm']. -apply generic_format_F2R. -unfold canonic_exp. -rewrite ln_beta_F2R_Zdigits. -2: now apply Zgt_not_eq. -unfold k in Hk. clear -Hk. -omega. -rewrite <- Hm', F2R_0. -apply generic_format_0. -Qed. - -Theorem truncate_correct_format : - forall m e, - m <> Z0 -> - let x := F2R (Float beta m e) in - generic_format beta fexp x -> - (e <= fexp (Zdigits beta m + e))%Z -> - let '(m', e', l') := truncate (m, e, loc_Exact) in - x = F2R (Float beta m' e') /\ e' = canonic_exp beta fexp x. -Proof. -intros m e Hm x Fx He. -assert (Hc: canonic_exp beta fexp x = fexp (Zdigits beta m + e)). -unfold canonic_exp, x. -now rewrite ln_beta_F2R_Zdigits. -unfold truncate. -rewrite <- Hc. -set (k := (canonic_exp beta fexp x - e)%Z). -case Zlt_bool_spec ; intros Hk. -(* *) -unfold truncate_aux. -rewrite Fx at 1. -assert (H: (e + k)%Z = canonic_exp beta fexp x). -unfold k. ring. -refine (conj _ H). -rewrite <- H. -apply F2R_eq_compat. -replace (scaled_mantissa beta fexp x) with (Z2R (Zfloor (scaled_mantissa beta fexp x))). -rewrite Ztrunc_Z2R. -unfold scaled_mantissa. -rewrite <- H. -unfold x, F2R. simpl. -rewrite Rmult_assoc, <- bpow_plus. -ring_simplify (e + - (e + k))%Z. -clear -Hm Hk. -destruct k as [|k|k] ; try easy. -simpl. -apply Zfloor_div. -intros H. -generalize (Zpower_pos_gt_0 beta k) (Zle_bool_imp_le _ _ (radix_prop beta)). -omega. -rewrite scaled_mantissa_generic with (1 := Fx). -now rewrite Zfloor_Z2R. -(* *) -split. -apply refl_equal. -unfold k in Hk. -omega. -Qed. - -Theorem truncate_correct_partial : - forall x m e l, - (0 < x)%R -> - inbetween_float beta m e x l -> - (e <= fexp (Zdigits beta m + e))%Z -> - let '(m', e', l') := truncate (m, e, l) in - inbetween_float beta m' e' x l' /\ e' = canonic_exp beta fexp x. -Proof. -intros x m e l Hx H1 H2. -unfold truncate. -set (k := (fexp (Zdigits beta m + e) - e)%Z). -set (p := Zpower beta k). -(* *) -assert (Hx': (F2R (Float beta m e) <= x < F2R (Float beta (m + 1) e))%R). -apply inbetween_float_bounds with (1 := H1). -(* *) -assert (Hm: (0 <= m)%Z). -cut (0 < m + 1)%Z. omega. -apply F2R_lt_reg with beta e. -rewrite F2R_0. -apply Rlt_trans with (1 := Hx). -apply Hx'. -assert (He: (e + k = canonic_exp beta fexp x)%Z). -(* . *) -unfold canonic_exp. -destruct (Zle_lt_or_eq _ _ Hm) as [Hm'|Hm']. -(* .. 0 < m *) -rewrite ln_beta_F2R_bounds with (1 := Hm') (2 := Hx'). -assert (H: m <> Z0). -apply sym_not_eq. -now apply Zlt_not_eq. -rewrite ln_beta_F2R with (1 := H). -rewrite <- Zdigits_ln_beta with (1 := H). -unfold k. -ring. -(* .. m = 0 *) -rewrite <- Hm' in H2. -destruct (ln_beta beta x) as (ex, Hex). -simpl. -specialize (Hex (Rgt_not_eq _ _ Hx)). -unfold k. -ring_simplify. -rewrite <- Hm'. -simpl. -apply sym_eq. -apply valid_exp. -exact H2. -apply Zle_trans with e. -eapply bpow_lt_bpow. -apply Rle_lt_trans with (1 := proj1 Hex). -rewrite Rabs_pos_eq. -rewrite <- F2R_bpow. -rewrite <- Hm' in Hx'. -apply Hx'. -now apply Rlt_le. -exact H2. -(* . *) -generalize (Zlt_cases 0 k). -case (Zlt_bool 0 k) ; intros Hk ; unfold k in Hk. -split. -now apply inbetween_float_new_location. -exact He. -split. -exact H1. -rewrite <- He. -unfold k. -omega. -Qed. - -Theorem truncate_correct : - forall x m e l, - (0 <= x)%R -> - inbetween_float beta m e x l -> - (e <= fexp (Zdigits beta m + e))%Z \/ l = loc_Exact -> - let '(m', e', l') := truncate (m, e, l) in - inbetween_float beta m' e' x l' /\ - (e' = canonic_exp beta fexp x \/ (l' = loc_Exact /\ format x)). -Proof. -intros x m e l [Hx|Hx] H1 H2. -(* 0 < x *) -destruct (Zle_or_lt e (fexp (Zdigits beta m + e))) as [H3|H3]. -(* . enough digits *) -generalize (truncate_correct_partial x m e l Hx H1 H3). -destruct (truncate (m, e, l)) as ((m', e'), l'). -intros (H4, H5). -split. -exact H4. -now left. -(* . not enough digits but loc_Exact *) -destruct H2 as [H2|H2]. -elim (Zlt_irrefl e). -now apply Zle_lt_trans with (1 := H2). -rewrite H2 in H1 |- *. -unfold truncate. -generalize (Zlt_cases 0 (fexp (Zdigits beta m + e) - e)). -case Zlt_bool. -intros H. -apply False_ind. -omega. -intros _. -split. -exact H1. -right. -split. -apply refl_equal. -inversion_clear H1. -rewrite H. -apply generic_format_F2R. -intros Zm. -unfold canonic_exp. -rewrite ln_beta_F2R_Zdigits with (1 := Zm). -now apply Zlt_le_weak. -(* x = 0 *) -assert (Hm: m = Z0). -cut (m <= 0 < m + 1)%Z. omega. -assert (Hx': (F2R (Float beta m e) <= x < F2R (Float beta (m + 1) e))%R). -apply inbetween_float_bounds with (1 := H1). -rewrite <- Hx in Hx'. -split. -apply F2R_le_0_reg with (1 := proj1 Hx'). -apply F2R_gt_0_reg with (1 := proj2 Hx'). -rewrite Hm, <- Hx in H1 |- *. -clear -H1. -case H1. -(* . *) -intros _. -assert (exists e', truncate (Z0, e, loc_Exact) = (Z0, e', loc_Exact)). -unfold truncate, truncate_aux. -case Zlt_bool. -rewrite Zdiv_0_l, Zmod_0_l. -eexists. -apply f_equal. -unfold new_location. -now case Zeven. -now eexists. -destruct H as (e', H). -rewrite H. -split. -constructor. -apply sym_eq. -apply F2R_0. -right. -repeat split. -apply generic_format_0. -(* . *) -intros l' (H, _) _. -rewrite F2R_0 in H. -elim Rlt_irrefl with (1 := H). -Qed. - -Section round_dir. - -Variable rnd : R -> Z. -Context { valid_rnd : Valid_rnd rnd }. - -Variable choice : Z -> location -> Z. -Hypothesis inbetween_int_valid : - forall x m l, - inbetween_int m x l -> - rnd x = choice m l. - -Theorem round_any_correct : - forall x m e l, - inbetween_float beta m e x l -> - (e = canonic_exp beta fexp x \/ (l = loc_Exact /\ format x)) -> - round beta fexp rnd x = F2R (Float beta (choice m l) e). -Proof with auto with typeclass_instances. -intros x m e l Hin [He|(Hl,Hf)]. -rewrite He in Hin |- *. -apply inbetween_float_round with (2 := Hin). -exact inbetween_int_valid. -rewrite Hl in Hin. -inversion_clear Hin. -rewrite Hl. -replace (choice m loc_Exact) with m. -rewrite <- H. -apply round_generic... -rewrite <- (Zrnd_Z2R rnd m) at 1. -apply inbetween_int_valid. -now constructor. -Qed. - -(** Truncating a triple is sufficient to round a real number. *) - -Theorem round_trunc_any_correct : - forall x m e l, - (0 <= x)%R -> - inbetween_float beta m e x l -> - (e <= fexp (Zdigits beta m + e))%Z \/ l = loc_Exact -> - round beta fexp rnd x = let '(m', e', l') := truncate (m, e, l) in F2R (Float beta (choice m' l') e'). -Proof. -intros x m e l Hx Hl He. -generalize (truncate_correct x m e l Hx Hl He). -destruct (truncate (m, e, l)) as ((m', e'), l'). -intros (H1, H2). -now apply round_any_correct. -Qed. - -End round_dir. - -Section round_dir_sign. - -Variable rnd : R -> Z. -Context { valid_rnd : Valid_rnd rnd }. - -Variable choice : bool -> Z -> location -> Z. -Hypothesis inbetween_int_valid : - forall x m l, - inbetween_int m (Rabs x) l -> - rnd x = cond_Zopp (Rlt_bool x 0) (choice (Rlt_bool x 0) m l). - -Theorem round_sign_any_correct : - forall x m e l, - inbetween_float beta m e (Rabs x) l -> - (e = canonic_exp beta fexp x \/ (l = loc_Exact /\ format x)) -> - round beta fexp rnd x = F2R (Float beta (cond_Zopp (Rlt_bool x 0) (choice (Rlt_bool x 0) m l)) e). -Proof with auto with typeclass_instances. -intros x m e l Hin [He|(Hl,Hf)]. -rewrite He in Hin |- *. -apply inbetween_float_round_sign with (2 := Hin). -exact inbetween_int_valid. -rewrite Hl in Hin. -inversion_clear Hin. -rewrite Hl. -replace (choice (Rlt_bool x 0) m loc_Exact) with m. -(* *) -unfold Rabs in H. -destruct (Rcase_abs x) as [Zx|Zx]. -rewrite Rlt_bool_true with (1 := Zx). -simpl. -rewrite F2R_Zopp. -rewrite <- H, Ropp_involutive. -apply round_generic... -rewrite Rlt_bool_false. -simpl. -rewrite <- H. -now apply round_generic. -now apply Rge_le. -(* *) -destruct (Rlt_bool_spec x 0) as [Zx|Zx]. -(* . *) -apply Zopp_inj. -change (- m = cond_Zopp true (choice true m loc_Exact))%Z. -rewrite <- (Zrnd_Z2R rnd (-m)) at 1. -assert (Z2R (-m) < 0)%R. -rewrite Z2R_opp. -apply Ropp_lt_gt_0_contravar. -apply (Z2R_lt 0). -apply F2R_gt_0_reg with beta e. -rewrite <- H. -apply Rabs_pos_lt. -now apply Rlt_not_eq. -rewrite <- Rlt_bool_true with (1 := H0). -apply inbetween_int_valid. -constructor. -rewrite Rabs_left with (1 := H0). -rewrite Z2R_opp. -apply Ropp_involutive. -(* . *) -change (m = cond_Zopp false (choice false m loc_Exact))%Z. -rewrite <- (Zrnd_Z2R rnd m) at 1. -assert (0 <= Z2R m)%R. -apply (Z2R_le 0). -apply F2R_ge_0_reg with beta e. -rewrite <- H. -apply Rabs_pos. -rewrite <- Rlt_bool_false with (1 := H0). -apply inbetween_int_valid. -constructor. -now apply Rabs_pos_eq. -Qed. - -(** Truncating a triple is sufficient to round a real number. *) - -Theorem round_trunc_sign_any_correct : - forall x m e l, - inbetween_float beta m e (Rabs x) l -> - (e <= fexp (Zdigits beta m + e))%Z \/ l = loc_Exact -> - round beta fexp rnd x = let '(m', e', l') := truncate (m, e, l) in F2R (Float beta (cond_Zopp (Rlt_bool x 0) (choice (Rlt_bool x 0) m' l')) e'). -Proof. -intros x m e l Hl He. -generalize (truncate_correct (Rabs x) m e l (Rabs_pos _) Hl He). -destruct (truncate (m, e, l)) as ((m', e'), l'). -intros (H1, H2). -apply round_sign_any_correct. -exact H1. -destruct H2 as [H2|(H2,H3)]. -left. -now rewrite <- canonic_exp_abs. -right. -split. -exact H2. -unfold Rabs in H3. -destruct (Rcase_abs x) in H3. -rewrite <- Ropp_involutive. -now apply generic_format_opp. -exact H3. -Qed. - -End round_dir_sign. - -(** * Instances of the theorems above, for the usual rounding modes. *) - -Definition round_DN_correct := - round_any_correct _ (fun m _ => m) inbetween_int_DN. - -Definition round_trunc_DN_correct := - round_trunc_any_correct _ (fun m _ => m) inbetween_int_DN. - -Definition round_sign_DN_correct := - round_sign_any_correct _ (fun s m l => cond_incr (round_sign_DN s l) m) inbetween_int_DN_sign. - -Definition round_trunc_sign_DN_correct := - round_trunc_sign_any_correct _ (fun s m l => cond_incr (round_sign_DN s l) m) inbetween_int_DN_sign. - -Definition round_UP_correct := - round_any_correct _ (fun m l => cond_incr (round_UP l) m) inbetween_int_UP. - -Definition round_trunc_UP_correct := - round_trunc_any_correct _ (fun m l => cond_incr (round_UP l) m) inbetween_int_UP. - -Definition round_sign_UP_correct := - round_sign_any_correct _ (fun s m l => cond_incr (round_sign_UP s l) m) inbetween_int_UP_sign. - -Definition round_trunc_sign_UP_correct := - round_trunc_sign_any_correct _ (fun s m l => cond_incr (round_sign_UP s l) m) inbetween_int_UP_sign. - -Definition round_ZR_correct := - round_any_correct _ (fun m l => cond_incr (round_ZR (Zlt_bool m 0) l) m) inbetween_int_ZR. - -Definition round_trunc_ZR_correct := - round_trunc_any_correct _ (fun m l => cond_incr (round_ZR (Zlt_bool m 0) l) m) inbetween_int_ZR. - -Definition round_sign_ZR_correct := - round_sign_any_correct _ (fun s m l => m) inbetween_int_ZR_sign. - -Definition round_trunc_sign_ZR_correct := - round_trunc_sign_any_correct _ (fun s m l => m) inbetween_int_ZR_sign. - -Definition round_NE_correct := - round_any_correct _ (fun m l => cond_incr (round_N (negb (Zeven m)) l) m) inbetween_int_NE. - -Definition round_trunc_NE_correct := - round_trunc_any_correct _ (fun m l => cond_incr (round_N (negb (Zeven m)) l) m) inbetween_int_NE. - -Definition round_sign_NE_correct := - round_sign_any_correct _ (fun s m l => cond_incr (round_N (negb (Zeven m)) l) m) inbetween_int_NE_sign. - -Definition round_trunc_sign_NE_correct := - round_trunc_sign_any_correct _ (fun s m l => cond_incr (round_N (negb (Zeven m)) l) m) inbetween_int_NE_sign. - -Definition round_NA_correct := - round_any_correct _ (fun m l => cond_incr (round_N (Zle_bool 0 m) l) m) inbetween_int_NA. - -Definition round_trunc_NA_correct := - round_trunc_any_correct _ (fun m l => cond_incr (round_N (Zle_bool 0 m) l) m) inbetween_int_NA. - -Definition round_sign_NA_correct := - round_sign_any_correct _ (fun s m l => cond_incr (round_N true l) m) inbetween_int_NA_sign. - -Definition round_trunc_sign_NA_correct := - round_trunc_sign_any_correct _ (fun s m l => cond_incr (round_N true l) m) inbetween_int_NA_sign. - -End Fcalc_round_fexp. - -(** Specialization of truncate for FIX formats. *) - -Variable emin : Z. - -Definition truncate_FIX t := - let '(m, e, l) := t in - let k := (emin - e)%Z in - if Zlt_bool 0 k then - let p := Zpower beta k in - (Zdiv m p, (e + k)%Z, new_location p (Zmod m p) l) - else t. - -Theorem truncate_FIX_correct : - forall x m e l, - inbetween_float beta m e x l -> - (e <= emin)%Z \/ l = loc_Exact -> - let '(m', e', l') := truncate_FIX (m, e, l) in - inbetween_float beta m' e' x l' /\ - (e' = canonic_exp beta (FIX_exp emin) x \/ (l' = loc_Exact /\ generic_format beta (FIX_exp emin) x)). -Proof. -intros x m e l H1 H2. -unfold truncate_FIX. -set (k := (emin - e)%Z). -set (p := Zpower beta k). -unfold canonic_exp, FIX_exp. -generalize (Zlt_cases 0 k). -case (Zlt_bool 0 k) ; intros Hk. -(* shift *) -split. -now apply inbetween_float_new_location. -clear H2. -left. -unfold k. -ring. -(* no shift *) -split. -exact H1. -unfold k in Hk. -destruct H2 as [H2|H2]. -left. -omega. -right. -split. -exact H2. -rewrite H2 in H1. -inversion_clear H1. -rewrite H. -apply generic_format_F2R. -unfold canonic_exp. -omega. -Qed. - -End Fcalc_round. diff --git a/flocq/Calc/Fcalc_sqrt.v b/flocq/Calc/Fcalc_sqrt.v deleted file mode 100644 index 5f541d83..00000000 --- a/flocq/Calc/Fcalc_sqrt.v +++ /dev/null @@ -1,244 +0,0 @@ -(** -This file is part of the Flocq formalization of floating-point -arithmetic in Coq: http://flocq.gforge.inria.fr/ - -Copyright (C) 2010-2013 Sylvie Boldo -#
# -Copyright (C) 2010-2013 Guillaume Melquiond - -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. -*) - -(** * Helper functions and theorems for computing the rounded square root of a floating-point number. *) - -Require Import Fcore_Raux. -Require Import Fcore_defs. -Require Import Fcore_digits. -Require Import Fcore_float_prop. -Require Import Fcalc_bracket. -Require Import Fcalc_digits. - -Section Fcalc_sqrt. - -Variable beta : radix. -Notation bpow e := (bpow beta e). - -(** 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. - -The algorithm performs the following steps: -- Shift the mantissa so that it has at least 2p-1 digits; - shift it one digit more if the new exponent is not even. -- Compute the square root s (at least p digits) of the new - mantissa, and its remainder r. -- Compute the position according to the remainder: - -- r == 0 => Eq, - -- r <= s => Lo, - -- r >= s => Up. - -Complexity is fine as long as p1 <= 2p-1. -*) - -Definition Fsqrt_core prec 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 Zeven 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 (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). - -Theorem Fsqrt_core_correct : - forall prec m e, - (0 < m)%Z -> - let '(m', e', l) := Fsqrt_core prec m e in - (prec <= Zdigits beta m')%Z /\ - inbetween_float beta m' e' (sqrt (F2R (Float beta m e))) l. -Proof. -intros prec m e Hm. -unfold Fsqrt_core. -set (d := Zdigits beta m). -set (s := Zmax (2 * prec - d) 0). -(* . exponent *) -case_eq (if Zeven (e - s) then (s, (e - s)%Z) else ((s + 1)%Z, (e - s - 1)%Z)). -intros s' e' Hse. -assert (He: (Zeven e' = true /\ 0 <= s' /\ 2 * prec - d <= s' /\ s' + e' = e)%Z). -revert Hse. -case_eq (Zeven (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 Zeven_plus, 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. -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. -rewrite sqrt_mult. -2: now apply (Z2R_le 0) ; apply Zlt_le_weak. -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. -rewrite bpow_plus. -fold (Rsqr (bpow e2)). -rewrite sqrt_Rsqr. -2: apply Rlt_le ; apply bpow_gt_0. -apply inbetween_mult_compat. -apply bpow_gt_0. -rewrite Hq. -case Zeq_bool_spec ; intros Hr'. -(* .. r = 0 *) -rewrite Hr', Zplus_0_r, Z2R_mult. -fold (Rsqr (Z2R q)). -rewrite sqrt_Rsqr. -now constructor. -apply (Z2R_le 0). -omega. -(* .. r <> 0 *) -constructor. -split. -(* ... bounds *) -apply Rle_lt_trans with (sqrt (Z2R (q * q))). -rewrite Z2R_mult. -fold (Rsqr (Z2R q)). -rewrite sqrt_Rsqr. -apply Rle_refl. -apply (Z2R_le 0). -omega. -apply sqrt_lt_1. -rewrite Z2R_mult. -apply Rle_0_sqr. -rewrite <- Hq. -apply (Z2R_le 0). -now apply Zlt_le_weak. -apply Z2R_lt. -omega. -apply Rlt_le_trans with (sqrt (Z2R ((q + 1) * (q + 1)))). -apply sqrt_lt_1. -rewrite <- Hq. -apply (Z2R_le 0). -now apply Zlt_le_weak. -rewrite Z2R_mult. -apply Rle_0_sqr. -apply Z2R_lt. -ring_simplify. -omega. -rewrite Z2R_mult. -fold (Rsqr (Z2R (q + 1))). -rewrite sqrt_Rsqr. -apply Rle_refl. -apply (Z2R_le 0). -omega. -(* ... location *) -rewrite Rcompare_half_r. -rewrite <- Rcompare_sqr. -replace ((2 * sqrt (Z2R (q * q + r))) * (2 * sqrt (Z2R (q * q + r))))%R - with (4 * Rsqr (sqrt (Z2R (q * q + r))))%R by (unfold Rsqr ; ring). -rewrite Rsqr_sqrt. -change 4%R with (Z2R 4). -rewrite <- Z2R_plus, <- 2!Z2R_mult. -rewrite Rcompare_Z2R. -replace ((q + (q + 1)) * (q + (q + 1)))%Z with (4 * (q * q) + 4 * q + 1)%Z by ring. -generalize (Zle_cases r q). -case (Zle_bool r q) ; intros Hr''. -change (4 * (q * q + r) < 4 * (q * q) + 4 * q + 1)%Z. -omega. -change (4 * (q * q + r) > 4 * (q * q) + 4 * q + 1)%Z. -omega. -rewrite <- Hq. -apply (Z2R_le 0). -now apply Zlt_le_weak. -apply Rmult_le_pos. -now apply (Z2R_le 0 2). -apply sqrt_ge_0. -rewrite <- Z2R_plus. -apply (Z2R_le 0). -omega. -Qed. - -End Fcalc_sqrt. diff --git a/flocq/Calc/Operations.v b/flocq/Calc/Operations.v new file mode 100644 index 00000000..3416cb4e --- /dev/null +++ b/flocq/Calc/Operations.v @@ -0,0 +1,164 @@ +(** +This file is part of the Flocq formalization of floating-point +arithmetic in Coq: http://flocq.gforge.inria.fr/ + +Copyright (C) 2009-2018 Sylvie Boldo +#
# +Copyright (C) 2009-2018 Guillaume Melquiond + +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. +*) + +(** Basic operations on floats: alignment, addition, multiplication *) +Require Import Raux Defs Float_prop. + +Set Implicit Arguments. +Set Strongly Strict Implicit. + +Section Float_ops. + +Variable beta : radix. + +Notation bpow e := (bpow beta e). + +Arguments Float {beta}. + +Definition Falign (f1 f2 : float beta) := + let '(Float m1 e1) := f1 in + let '(Float m2 e2) := f2 in + if Zle_bool e1 e2 + then (m1, (m2 * Zpower beta (e2 - e1))%Z, e1) + else ((m1 * Zpower beta (e1 - e2))%Z, m2, e2). + +Theorem Falign_spec : + forall f1 f2 : float beta, + let '(m1, m2, e) := Falign f1 f2 in + F2R f1 = @F2R beta (Float m1 e) /\ F2R f2 = @F2R beta (Float m2 e). +Proof. +unfold Falign. +intros (m1, e1) (m2, e2). +generalize (Zle_cases e1 e2). +case (Zle_bool e1 e2) ; intros He ; split ; trivial. +now rewrite <- F2R_change_exp. +rewrite <- F2R_change_exp. +apply refl_equal. +omega. +Qed. + +Theorem Falign_spec_exp: + forall f1 f2 : float beta, + snd (Falign f1 f2) = Z.min (Fexp f1) (Fexp f2). +Proof. +intros (m1,e1) (m2,e2). +unfold Falign; simpl. +generalize (Zle_cases e1 e2);case (Zle_bool e1 e2); intros He. +case (Zmin_spec e1 e2); intros (H1,H2); easy. +case (Zmin_spec e1 e2); intros (H1,H2); easy. +Qed. + +Definition Fopp (f1 : float beta) : float beta := + let '(Float m1 e1) := f1 in + Float (-m1)%Z e1. + +Theorem F2R_opp : + forall f1 : float beta, + (F2R (Fopp f1) = -F2R f1)%R. +intros (m1,e1). +apply F2R_Zopp. +Qed. + +Definition Fabs (f1 : float beta) : float beta := + let '(Float m1 e1) := f1 in + Float (Z.abs m1)%Z e1. + +Theorem F2R_abs : + forall f1 : float beta, + (F2R (Fabs f1) = Rabs (F2R f1))%R. +intros (m1,e1). +apply F2R_Zabs. +Qed. + +Definition Fplus (f1 f2 : float beta) : float beta := + let '(m1, m2 ,e) := Falign f1 f2 in + Float (m1 + m2) e. + +Theorem F2R_plus : + forall f1 f2 : float beta, + F2R (Fplus f1 f2) = (F2R f1 + F2R f2)%R. +Proof. +intros f1 f2. +unfold Fplus. +generalize (Falign_spec f1 f2). +destruct (Falign f1 f2) as ((m1, m2), e). +intros (H1, H2). +rewrite H1, H2. +unfold F2R. simpl. +rewrite plus_IZR. +apply Rmult_plus_distr_r. +Qed. + +Theorem Fplus_same_exp : + forall m1 m2 e, + Fplus (Float m1 e) (Float m2 e) = Float (m1 + m2) e. +Proof. +intros m1 m2 e. +unfold Fplus. +simpl. +now rewrite Zle_bool_refl, Zminus_diag, Zmult_1_r. +Qed. + +Theorem Fexp_Fplus : + forall f1 f2 : float beta, + Fexp (Fplus f1 f2) = Z.min (Fexp f1) (Fexp f2). +Proof. +intros f1 f2. +unfold Fplus. +rewrite <- Falign_spec_exp. +now destruct (Falign f1 f2) as ((p,q),e). +Qed. + +Definition Fminus (f1 f2 : float beta) := + Fplus f1 (Fopp f2). + +Theorem F2R_minus : + forall f1 f2 : float beta, + F2R (Fminus f1 f2) = (F2R f1 - F2R f2)%R. +Proof. +intros f1 f2; unfold Fminus. +rewrite F2R_plus, F2R_opp. +ring. +Qed. + +Theorem Fminus_same_exp : + forall m1 m2 e, + Fminus (Float m1 e) (Float m2 e) = Float (m1 - m2) e. +Proof. +intros m1 m2 e. +unfold Fminus. +apply Fplus_same_exp. +Qed. + +Definition Fmult (f1 f2 : float beta) : float beta := + let '(Float m1 e1) := f1 in + let '(Float m2 e2) := f2 in + Float (m1 * m2) (e1 + e2). + +Theorem F2R_mult : + forall f1 f2 : float beta, + F2R (Fmult f1 f2) = (F2R f1 * F2R f2)%R. +Proof. +intros (m1, e1) (m2, e2). +unfold Fmult, F2R. simpl. +rewrite mult_IZR, bpow_plus. +ring. +Qed. + +End Float_ops. diff --git a/flocq/Calc/Round.v b/flocq/Calc/Round.v new file mode 100644 index 00000000..5bde6af4 --- /dev/null +++ b/flocq/Calc/Round.v @@ -0,0 +1,1171 @@ +(** +This file is part of the Flocq formalization of floating-point +arithmetic in Coq: http://flocq.gforge.inria.fr/ + +Copyright (C) 2010-2018 Sylvie Boldo +#
# +Copyright (C) 2010-2018 Guillaume Melquiond + +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. +*) + +(** * Helper function for computing the rounded value of a real number. *) + +Require Import Core Digits Float_prop Bracket. + +Section Fcalc_round. + +Variable beta : radix. +Notation bpow e := (bpow beta e). + +Section Fcalc_round_fexp. + +Variable fexp : Z -> Z. +Context { valid_exp : Valid_exp fexp }. +Notation format := (generic_format beta fexp). + +Theorem cexp_inbetween_float : + forall x m e l, + (0 < x)%R -> + inbetween_float beta m e x l -> + (e <= cexp beta fexp x \/ e <= fexp (Zdigits beta m + e))%Z -> + cexp beta fexp x = fexp (Zdigits beta m + e). +Proof. +intros x m e l Px Bx He. +unfold cexp. +apply inbetween_float_bounds in Bx. +assert (0 <= m)%Z as Hm. +{ apply Zlt_succ_le. + eapply gt_0_F2R. + apply Rlt_trans with (1 := Px). + apply Bx. } +destruct (Zle_lt_or_eq _ _ Hm) as [Hm'|<-]. + now erewrite <- mag_F2R_bounds_Zdigits with (1 := Hm'). +clear Hm. +assert (mag beta x <= e)%Z as Hx. +{ apply mag_le_bpow. + now apply Rgt_not_eq. + rewrite Rabs_pos_eq. + now rewrite <- F2R_bpow. + now apply Rlt_le. } +simpl in He |- *. +clear Bx. +destruct He as [He|He]. +- apply eq_sym, valid_exp with (2 := He). + now apply Z.le_trans with e. +- apply valid_exp with (1 := He). + now apply Z.le_trans with e. +Qed. + +Theorem cexp_inbetween_float_loc_Exact : + forall x m e l, + (0 <= x)%R -> + inbetween_float beta m e x l -> + (e <= cexp beta fexp x \/ l = loc_Exact <-> + e <= fexp (Zdigits beta m + e) \/ l = loc_Exact)%Z. +Proof. +intros x m e l Px Bx. +destruct Px as [Px|Px]. +- split ; (intros [H|H] ; [left|now right]). + rewrite <- cexp_inbetween_float with (1 := Px) (2 := Bx). + exact H. + now left. + rewrite cexp_inbetween_float with (1 := Px) (2 := Bx). + exact H. + now right. +- assert (H := Bx). + destruct Bx as [|c Bx _]. + now split ; right. + rewrite <- Px in Bx. + destruct Bx as [Bx1 Bx2]. + apply lt_0_F2R in Bx1. + apply gt_0_F2R in Bx2. + omega. +Qed. + +(** Relates location and rounding. *) + +Theorem inbetween_float_round : + forall rnd choice, + ( forall x m l, inbetween_int m x l -> rnd x = choice m l ) -> + forall x m l, + let e := cexp beta fexp x in + inbetween_float beta m e x l -> + round beta fexp rnd x = F2R (Float beta (choice m l) e). +Proof. +intros rnd choice Hc x m l e Hl. +unfold round, F2R. simpl. +apply (f_equal (fun m => (IZR m * bpow e)%R)). +apply Hc. +apply inbetween_mult_reg with (bpow e). +apply bpow_gt_0. +now rewrite scaled_mantissa_mult_bpow. +Qed. + +Definition cond_incr (b : bool) m := if b then (m + 1)%Z else m. + +Theorem inbetween_float_round_sign : + forall rnd choice, + ( forall x m l, inbetween_int m (Rabs x) l -> + rnd x = cond_Zopp (Rlt_bool x 0) (choice (Rlt_bool x 0) m l) ) -> + forall x m l, + let e := cexp beta fexp x in + inbetween_float beta m e (Rabs x) l -> + round beta fexp rnd x = F2R (Float beta (cond_Zopp (Rlt_bool x 0) (choice (Rlt_bool x 0) m l)) e). +Proof. +intros rnd choice Hc x m l e Hx. +apply (f_equal (fun m => (IZR m * bpow e)%R)). +simpl. +replace (Rlt_bool x 0) with (Rlt_bool (scaled_mantissa beta fexp x) 0). +(* *) +apply Hc. +apply inbetween_mult_reg with (bpow e). +apply bpow_gt_0. +rewrite <- (Rabs_right (bpow e)) at 3. +rewrite <- Rabs_mult. +now rewrite scaled_mantissa_mult_bpow. +apply Rle_ge. +apply bpow_ge_0. +(* *) +destruct (Rlt_bool_spec x 0) as [Zx|Zx] ; simpl. +apply Rlt_bool_true. +rewrite <- (Rmult_0_l (bpow (-e))). +apply Rmult_lt_compat_r with (2 := Zx). +apply bpow_gt_0. +apply Rlt_bool_false. +apply Rmult_le_pos with (1 := Zx). +apply bpow_ge_0. +Qed. + +(** Relates location and rounding down. *) + +Theorem inbetween_int_DN : + forall x m l, + inbetween_int m x l -> + Zfloor x = m. +Proof. +intros x m l Hl. +refine (Zfloor_imp m _ _). +apply inbetween_bounds with (2 := Hl). +apply IZR_lt. +apply Zlt_succ. +Qed. + +Theorem inbetween_float_DN : + forall x m l, + let e := cexp beta fexp x in + inbetween_float beta m e x l -> + round beta fexp Zfloor x = F2R (Float beta m e). +Proof. +apply inbetween_float_round with (choice := fun m l => m). +exact inbetween_int_DN. +Qed. + +Definition round_sign_DN s l := + match l with + | loc_Exact => false + | _ => s + end. + +Theorem inbetween_int_DN_sign : + forall x m l, + inbetween_int m (Rabs x) l -> + Zfloor x = cond_Zopp (Rlt_bool x 0) (cond_incr (round_sign_DN (Rlt_bool x 0) l) m). +Proof. +intros x m l Hl. +unfold Rabs in Hl. +destruct (Rcase_abs x) as [Zx|Zx] . +(* *) +rewrite Rlt_bool_true with (1 := Zx). +inversion_clear Hl ; simpl. +rewrite <- (Ropp_involutive x). +rewrite H, <- opp_IZR. +apply Zfloor_IZR. +apply Zfloor_imp. +split. +apply Rlt_le. +rewrite opp_IZR. +apply Ropp_lt_cancel. +now rewrite Ropp_involutive. +ring_simplify (- (m + 1) + 1)%Z. +rewrite opp_IZR. +apply Ropp_lt_cancel. +now rewrite Ropp_involutive. +(* *) +rewrite Rlt_bool_false. +inversion_clear Hl ; simpl. +rewrite H. +apply Zfloor_IZR. +apply Zfloor_imp. +split. +now apply Rlt_le. +apply H. +now apply Rge_le. +Qed. + +Theorem inbetween_float_DN_sign : + forall x m l, + let e := cexp beta fexp x in + inbetween_float beta m e (Rabs x) l -> + round beta fexp Zfloor x = F2R (Float beta (cond_Zopp (Rlt_bool x 0) (cond_incr (round_sign_DN (Rlt_bool x 0) l) m)) e). +Proof. +apply inbetween_float_round_sign with (choice := fun s m l => cond_incr (round_sign_DN s l) m). +exact inbetween_int_DN_sign. +Qed. + +(** Relates location and rounding up. *) + +Definition round_UP l := + match l with + | loc_Exact => false + | _ => true + end. + +Theorem inbetween_int_UP : + forall x m l, + inbetween_int m x l -> + Zceil x = cond_incr (round_UP l) m. +Proof. +intros x m l Hl. +assert (Hl': l = loc_Exact \/ (l <> loc_Exact /\ round_UP l = true)). +case l ; try (now left) ; now right ; split. +destruct Hl' as [Hl'|(Hl1, Hl2)]. +(* Exact *) +rewrite Hl'. +destruct Hl ; try easy. +rewrite H. +exact (Zceil_IZR _). +(* not Exact *) +rewrite Hl2. +simpl. +apply Zceil_imp. +ring_simplify (m + 1 - 1)%Z. +refine (let H := _ in conj (proj1 H) (Rlt_le _ _ (proj2 H))). +apply inbetween_bounds_not_Eq with (2 := Hl1) (1 := Hl). +Qed. + +Theorem inbetween_float_UP : + forall x m l, + let e := cexp beta fexp x in + inbetween_float beta m e x l -> + round beta fexp Zceil x = F2R (Float beta (cond_incr (round_UP l) m) e). +Proof. +apply inbetween_float_round with (choice := fun m l => cond_incr (round_UP l) m). +exact inbetween_int_UP. +Qed. + +Definition round_sign_UP s l := + match l with + | loc_Exact => false + | _ => negb s + end. + +Theorem inbetween_int_UP_sign : + forall x m l, + inbetween_int m (Rabs x) l -> + Zceil x = cond_Zopp (Rlt_bool x 0) (cond_incr (round_sign_UP (Rlt_bool x 0) l) m). +Proof. +intros x m l Hl. +unfold Rabs in Hl. +destruct (Rcase_abs x) as [Zx|Zx] . +(* *) +rewrite Rlt_bool_true with (1 := Zx). +simpl. +unfold Zceil. +apply f_equal. +inversion_clear Hl ; simpl. +rewrite H. +apply Zfloor_IZR. +apply Zfloor_imp. +split. +now apply Rlt_le. +apply H. +(* *) +rewrite Rlt_bool_false. +simpl. +inversion_clear Hl ; simpl. +rewrite H. +apply Zceil_IZR. +apply Zceil_imp. +split. +change (m + 1 - 1)%Z with (Z.pred (Z.succ m)). +now rewrite <- Zpred_succ. +now apply Rlt_le. +now apply Rge_le. +Qed. + +Theorem inbetween_float_UP_sign : + forall x m l, + let e := cexp beta fexp x in + inbetween_float beta m e (Rabs x) l -> + round beta fexp Zceil x = F2R (Float beta (cond_Zopp (Rlt_bool x 0) (cond_incr (round_sign_UP (Rlt_bool x 0) l) m)) e). +Proof. +apply inbetween_float_round_sign with (choice := fun s m l => cond_incr (round_sign_UP s l) m). +exact inbetween_int_UP_sign. +Qed. + +(** Relates location and rounding toward zero. *) + +Definition round_ZR (s : bool) l := + match l with + | loc_Exact => false + | _ => s + end. + +Theorem inbetween_int_ZR : + forall x m l, + inbetween_int m x l -> + Ztrunc x = cond_incr (round_ZR (Zlt_bool m 0) l) m. +Proof with auto with typeclass_instances. +intros x m l Hl. +inversion_clear Hl as [Hx|l' Hx Hl']. +(* Exact *) +rewrite Hx. +rewrite Zrnd_IZR... +(* not Exact *) +unfold Ztrunc. +assert (Hm: Zfloor x = m). +apply Zfloor_imp. +exact (conj (Rlt_le _ _ (proj1 Hx)) (proj2 Hx)). +rewrite Zceil_floor_neq. +rewrite Hm. +unfold cond_incr. +simpl. +case Rlt_bool_spec ; intros Hx' ; + case Zlt_bool_spec ; intros Hm' ; try apply refl_equal. +elim Rlt_not_le with (1 := Hx'). +apply Rlt_le. +apply Rle_lt_trans with (2 := proj1 Hx). +now apply IZR_le. +elim Rle_not_lt with (1 := Hx'). +apply Rlt_le_trans with (1 := proj2 Hx). +apply IZR_le. +now apply Zlt_le_succ. +rewrite Hm. +now apply Rlt_not_eq. +Qed. + +Theorem inbetween_float_ZR : + forall x m l, + let e := cexp beta fexp x in + inbetween_float beta m e x l -> + round beta fexp Ztrunc x = F2R (Float beta (cond_incr (round_ZR (Zlt_bool m 0) l) m) e). +Proof. +apply inbetween_float_round with (choice := fun m l => cond_incr (round_ZR (Zlt_bool m 0) l) m). +exact inbetween_int_ZR. +Qed. + +Theorem inbetween_int_ZR_sign : + forall x m l, + inbetween_int m (Rabs x) l -> + Ztrunc x = cond_Zopp (Rlt_bool x 0) m. +Proof. +intros x m l Hl. +simpl. +unfold Ztrunc. +destruct (Rlt_le_dec x 0) as [Zx|Zx]. +(* *) +rewrite Rlt_bool_true with (1 := Zx). +simpl. +unfold Zceil. +apply f_equal. +apply Zfloor_imp. +rewrite <- Rabs_left with (1 := Zx). +apply inbetween_bounds with (2 := Hl). +apply IZR_lt. +apply Zlt_succ. +(* *) +rewrite Rlt_bool_false with (1 := Zx). +simpl. +apply Zfloor_imp. +rewrite <- Rabs_pos_eq with (1 := Zx). +apply inbetween_bounds with (2 := Hl). +apply IZR_lt. +apply Zlt_succ. +Qed. + +Theorem inbetween_float_ZR_sign : + forall x m l, + let e := cexp beta fexp x in + inbetween_float beta m e (Rabs x) l -> + round beta fexp Ztrunc x = F2R (Float beta (cond_Zopp (Rlt_bool x 0) m) e). +Proof. +apply inbetween_float_round_sign with (choice := fun s m l => m). +exact inbetween_int_ZR_sign. +Qed. + +(** Relates location and rounding to nearest. *) + +Definition round_N (p : bool) l := + match l with + | loc_Exact => false + | loc_Inexact Lt => false + | loc_Inexact Eq => p + | loc_Inexact Gt => true + end. + +Theorem inbetween_int_N : + forall choice x m l, + inbetween_int m x l -> + Znearest choice x = cond_incr (round_N (choice m) l) m. +Proof with auto with typeclass_instances. +intros choice x m l Hl. +inversion_clear Hl as [Hx|l' Hx Hl']. +(* Exact *) +rewrite Hx. +rewrite Zrnd_IZR... +(* not Exact *) +unfold Znearest. +assert (Hm: Zfloor x = m). +apply Zfloor_imp. +exact (conj (Rlt_le _ _ (proj1 Hx)) (proj2 Hx)). +rewrite Zceil_floor_neq. +rewrite Hm. +replace (Rcompare (x - IZR m) (/2)) with l'. +now case l'. +rewrite <- Hl'. +rewrite plus_IZR. +rewrite <- (Rcompare_plus_r (- IZR m) x). +apply f_equal. +field. +rewrite Hm. +now apply Rlt_not_eq. +Qed. + +Theorem inbetween_int_N_sign : + forall choice x m l, + inbetween_int m (Rabs x) l -> + Znearest choice x = cond_Zopp (Rlt_bool x 0) (cond_incr (round_N (if Rlt_bool x 0 then negb (choice (-(m + 1))%Z) else choice m) l) m). +Proof with auto with typeclass_instances. +intros choice x m l Hl. +simpl. +unfold Rabs in Hl. +destruct (Rcase_abs x) as [Zx|Zx]. +(* *) +rewrite Rlt_bool_true with (1 := Zx). +simpl. +rewrite <- (Ropp_involutive x). +rewrite Znearest_opp. +apply f_equal. +inversion_clear Hl as [Hx|l' Hx Hl']. +rewrite Hx. +apply Zrnd_IZR... +assert (Hm: Zfloor (-x) = m). +apply Zfloor_imp. +exact (conj (Rlt_le _ _ (proj1 Hx)) (proj2 Hx)). +unfold Znearest. +rewrite Zceil_floor_neq. +rewrite Hm. +replace (Rcompare (- x - IZR m) (/2)) with l'. +now case l'. +rewrite <- Hl'. +rewrite plus_IZR. +rewrite <- (Rcompare_plus_r (- IZR m) (-x)). +apply f_equal. +field. +rewrite Hm. +now apply Rlt_not_eq. +(* *) +generalize (Rge_le _ _ Zx). +clear Zx. intros Zx. +rewrite Rlt_bool_false with (1 := Zx). +simpl. +inversion_clear Hl as [Hx|l' Hx Hl']. +rewrite Hx. +apply Zrnd_IZR... +assert (Hm: Zfloor x = m). +apply Zfloor_imp. +exact (conj (Rlt_le _ _ (proj1 Hx)) (proj2 Hx)). +unfold Znearest. +rewrite Zceil_floor_neq. +rewrite Hm. +replace (Rcompare (x - IZR m) (/2)) with l'. +now case l'. +rewrite <- Hl'. +rewrite plus_IZR. +rewrite <- (Rcompare_plus_r (- IZR m) x). +apply f_equal. +field. +rewrite Hm. +now apply Rlt_not_eq. +Qed. + +(** Relates location and rounding to nearest even. *) + +Theorem inbetween_int_NE : + forall x m l, + inbetween_int m x l -> + ZnearestE x = cond_incr (round_N (negb (Z.even m)) l) m. +Proof. +intros x m l Hl. +now apply inbetween_int_N with (choice := fun x => negb (Z.even x)). +Qed. + +Theorem inbetween_float_NE : + forall x m l, + let e := cexp beta fexp x in + inbetween_float beta m e x l -> + round beta fexp ZnearestE x = F2R (Float beta (cond_incr (round_N (negb (Z.even m)) l) m) e). +Proof. +apply inbetween_float_round with (choice := fun m l => cond_incr (round_N (negb (Z.even m)) l) m). +exact inbetween_int_NE. +Qed. + +Theorem inbetween_int_NE_sign : + forall x m l, + inbetween_int m (Rabs x) l -> + ZnearestE x = cond_Zopp (Rlt_bool x 0) (cond_incr (round_N (negb (Z.even m)) l) m). +Proof. +intros x m l Hl. +erewrite inbetween_int_N_sign with (choice := fun x => negb (Z.even x)). +2: eexact Hl. +apply f_equal. +case Rlt_bool. +rewrite Z.even_opp, Z.even_add. +now case (Z.even m). +apply refl_equal. +Qed. + +Theorem inbetween_float_NE_sign : + forall x m l, + let e := cexp beta fexp x in + inbetween_float beta m e (Rabs x) l -> + round beta fexp ZnearestE x = F2R (Float beta (cond_Zopp (Rlt_bool x 0) (cond_incr (round_N (negb (Z.even m)) l) m)) e). +Proof. +apply inbetween_float_round_sign with (choice := fun s m l => cond_incr (round_N (negb (Z.even m)) l) m). +exact inbetween_int_NE_sign. +Qed. + +(** Relates location and rounding to nearest away. *) + +Theorem inbetween_int_NA : + forall x m l, + inbetween_int m x l -> + ZnearestA x = cond_incr (round_N (Zle_bool 0 m) l) m. +Proof. +intros x m l Hl. +now apply inbetween_int_N with (choice := fun x => Zle_bool 0 x). +Qed. + +Theorem inbetween_float_NA : + forall x m l, + let e := cexp beta fexp x in + inbetween_float beta m e x l -> + round beta fexp ZnearestA x = F2R (Float beta (cond_incr (round_N (Zle_bool 0 m) l) m) e). +Proof. +apply inbetween_float_round with (choice := fun m l => cond_incr (round_N (Zle_bool 0 m) l) m). +exact inbetween_int_NA. +Qed. + +Theorem inbetween_int_NA_sign : + forall x m l, + inbetween_int m (Rabs x) l -> + ZnearestA x = cond_Zopp (Rlt_bool x 0) (cond_incr (round_N true l) m). +Proof. +intros x m l Hl. +erewrite inbetween_int_N_sign with (choice := Zle_bool 0). +2: eexact Hl. +apply f_equal. +assert (Hm: (0 <= m)%Z). +apply Zlt_succ_le. +apply lt_IZR. +apply Rle_lt_trans with (Rabs x). +apply Rabs_pos. +refine (proj2 (inbetween_bounds _ _ _ _ _ Hl)). +apply IZR_lt. +apply Zlt_succ. +rewrite Zle_bool_true with (1 := Hm). +rewrite Zle_bool_false. +now case Rlt_bool. +omega. +Qed. + +Definition truncate_aux t k := + let '(m, e, l) := t in + let p := Zpower beta k in + (Z.div m p, (e + k)%Z, new_location p (Zmod m p) l). + +Theorem truncate_aux_comp : + forall t k1 k2, + (0 < k1)%Z -> + (0 < k2)%Z -> + truncate_aux t (k1 + k2) = truncate_aux (truncate_aux t k1) k2. +Proof. +intros ((m,e),l) k1 k2 Hk1 Hk2. +unfold truncate_aux. +destruct (inbetween_float_ex beta m e l) as (x,Hx). +assert (B1 := inbetween_float_new_location _ _ _ _ _ _ Hk1 Hx). +assert (Hk3: (0 < k1 + k2)%Z). +change Z0 with (0 + 0)%Z. +now apply Zplus_lt_compat. +assert (B3 := inbetween_float_new_location _ _ _ _ _ _ Hk3 Hx). +assert (B2 := inbetween_float_new_location _ _ _ _ _ _ Hk2 B1). +rewrite Zplus_assoc in B3. +destruct (inbetween_float_unique _ _ _ _ _ _ _ B2 B3). +now rewrite H, H0, Zplus_assoc. +Qed. + +(** Given a triple (mantissa, exponent, position), this function + computes a triple with a canonic exponent, assuming the + original triple had enough precision. *) + +Definition truncate t := + let '(m, e, l) := t in + let k := (fexp (Zdigits beta m + e) - e)%Z in + if Zlt_bool 0 k then truncate_aux t k + else t. + +Theorem truncate_0 : + forall e l, + let '(m', e', l') := truncate (0, e, l)%Z in + m' = Z0. +Proof. +intros e l. +unfold truncate. +case Zlt_bool. +simpl. +apply Zdiv_0_l. +apply refl_equal. +Qed. + +Theorem generic_format_truncate : + forall m e l, + (0 <= m)%Z -> + let '(m', e', l') := truncate (m, e, l) in + format (F2R (Float beta m' e')). +Proof. +intros m e l Hm. +unfold truncate. +set (k := (fexp (Zdigits beta m + e) - e)%Z). +case Zlt_bool_spec ; intros Hk. +(* *) +unfold truncate_aux. +apply generic_format_F2R. +intros Hd. +unfold cexp. +rewrite mag_F2R_Zdigits with (1 := Hd). +rewrite Zdigits_div_Zpower with (1 := Hm). +replace (Zdigits beta m - k + (e + k))%Z with (Zdigits beta m + e)%Z by ring. +unfold k. +ring_simplify. +apply Z.le_refl. +split. +now apply Zlt_le_weak. +apply Znot_gt_le. +contradict Hd. +apply Zdiv_small. +apply conj with (1 := Hm). +rewrite <- Z.abs_eq with (1 := Hm). +apply Zpower_gt_Zdigits. +apply Zlt_le_weak. +now apply Z.gt_lt. +(* *) +destruct (Zle_lt_or_eq _ _ Hm) as [Hm'|Hm']. +apply generic_format_F2R. +unfold cexp. +rewrite mag_F2R_Zdigits. +2: now apply Zgt_not_eq. +unfold k in Hk. clear -Hk. +omega. +rewrite <- Hm', F2R_0. +apply generic_format_0. +Qed. + +Theorem truncate_correct_format : + forall m e, + m <> Z0 -> + let x := F2R (Float beta m e) in + generic_format beta fexp x -> + (e <= fexp (Zdigits beta m + e))%Z -> + let '(m', e', l') := truncate (m, e, loc_Exact) in + x = F2R (Float beta m' e') /\ e' = cexp beta fexp x. +Proof. +intros m e Hm x Fx He. +assert (Hc: cexp beta fexp x = fexp (Zdigits beta m + e)). +unfold cexp, x. +now rewrite mag_F2R_Zdigits. +unfold truncate. +rewrite <- Hc. +set (k := (cexp beta fexp x - e)%Z). +case Zlt_bool_spec ; intros Hk. +(* *) +unfold truncate_aux. +rewrite Fx at 1. +assert (H: (e + k)%Z = cexp beta fexp x). +unfold k. ring. +refine (conj _ H). +rewrite <- H. +apply F2R_eq. +replace (scaled_mantissa beta fexp x) with (IZR (Zfloor (scaled_mantissa beta fexp x))). +rewrite Ztrunc_IZR. +unfold scaled_mantissa. +rewrite <- H. +unfold x, F2R. simpl. +rewrite Rmult_assoc, <- bpow_plus. +ring_simplify (e + - (e + k))%Z. +clear -Hm Hk. +destruct k as [|k|k] ; try easy. +simpl. +apply Zfloor_div. +intros H. +generalize (Zpower_pos_gt_0 beta k) (Zle_bool_imp_le _ _ (radix_prop beta)). +omega. +rewrite scaled_mantissa_generic with (1 := Fx). +now rewrite Zfloor_IZR. +(* *) +split. +apply refl_equal. +unfold k in Hk. +omega. +Qed. + +Theorem truncate_correct_partial' : + forall x m e l, + (0 < x)%R -> + inbetween_float beta m e x l -> + (e <= cexp beta fexp x)%Z -> + let '(m', e', l') := truncate (m, e, l) in + inbetween_float beta m' e' x l' /\ e' = cexp beta fexp x. +Proof. +intros x m e l Hx H1 H2. +unfold truncate. +rewrite <- cexp_inbetween_float with (1 := Hx) (2 := H1) by now left. +generalize (Zlt_cases 0 (cexp beta fexp x - e)). +destruct Zlt_bool ; intros Hk. +- split. + now apply inbetween_float_new_location. + ring. +- apply (conj H1). + omega. +Qed. + +Theorem truncate_correct_partial : + forall x m e l, + (0 < x)%R -> + inbetween_float beta m e x l -> + (e <= fexp (Zdigits beta m + e))%Z -> + let '(m', e', l') := truncate (m, e, l) in + inbetween_float beta m' e' x l' /\ e' = cexp beta fexp x. +Proof. +intros x m e l Hx H1 H2. +apply truncate_correct_partial' with (1 := Hx) (2 := H1). +rewrite cexp_inbetween_float with (1 := Hx) (2 := H1). +exact H2. +now right. +Qed. + +Theorem truncate_correct' : + forall x m e l, + (0 <= x)%R -> + inbetween_float beta m e x l -> + (e <= cexp beta fexp x)%Z \/ l = loc_Exact -> + let '(m', e', l') := truncate (m, e, l) in + inbetween_float beta m' e' x l' /\ + (e' = cexp beta fexp x \/ (l' = loc_Exact /\ format x)). +Proof. +intros x m e l [Hx|Hx] H1 H2. +- destruct (Zle_or_lt e (fexp (Zdigits beta m + e))) as [H3|H3]. + + generalize (truncate_correct_partial x m e l Hx H1 H3). + destruct (truncate (m, e, l)) as [[m' e'] l']. + intros [H4 H5]. + apply (conj H4). + now left. + + destruct H2 as [H2|H2]. + generalize (truncate_correct_partial' x m e l Hx H1 H2). + destruct (truncate (m, e, l)) as [[m' e'] l']. + intros [H4 H5]. + apply (conj H4). + now left. + rewrite H2 in H1 |- *. + simpl. + generalize (Zlt_cases 0 (fexp (Zdigits beta m + e) - e)). + destruct Zlt_bool. + intros H. + apply False_ind. + omega. + intros _. + apply (conj H1). + right. + repeat split. + inversion_clear H1. + rewrite H. + apply generic_format_F2R. + intros Zm. + unfold cexp. + rewrite mag_F2R_Zdigits with (1 := Zm). + now apply Zlt_le_weak. +- assert (Hm: m = 0%Z). + cut (m <= 0 < m + 1)%Z. omega. + assert (F2R (Float beta m e) <= x < F2R (Float beta (m + 1) e))%R as Hx'. + apply inbetween_float_bounds with (1 := H1). + rewrite <- Hx in Hx'. + split. + apply le_0_F2R with (1 := proj1 Hx'). + apply gt_0_F2R with (1 := proj2 Hx'). + rewrite Hm, <- Hx in H1 |- *. + clear -H1. + destruct H1 as [_ | l' [H _] _]. + + assert (exists e', truncate (Z0, e, loc_Exact) = (Z0, e', loc_Exact)). + unfold truncate, truncate_aux. + case Zlt_bool. + rewrite Zdiv_0_l, Zmod_0_l. + eexists. + apply f_equal. + unfold new_location. + now case Z.even. + now eexists. + destruct H as [e' H]. + rewrite H. + split. + constructor. + apply eq_sym, F2R_0. + right. + repeat split. + apply generic_format_0. + + rewrite F2R_0 in H. + elim Rlt_irrefl with (1 := H). +Qed. + +Theorem truncate_correct : + forall x m e l, + (0 <= x)%R -> + inbetween_float beta m e x l -> + (e <= fexp (Zdigits beta m + e))%Z \/ l = loc_Exact -> + let '(m', e', l') := truncate (m, e, l) in + inbetween_float beta m' e' x l' /\ + (e' = cexp beta fexp x \/ (l' = loc_Exact /\ format x)). +Proof. +intros x m e l Hx H1 H2. +apply truncate_correct' with (1 := Hx) (2 := H1). +now apply cexp_inbetween_float_loc_Exact with (2 := H1). +Qed. + +Section round_dir. + +Variable rnd : R -> Z. +Context { valid_rnd : Valid_rnd rnd }. + +Variable choice : Z -> location -> Z. +Hypothesis inbetween_int_valid : + forall x m l, + inbetween_int m x l -> + rnd x = choice m l. + +Theorem round_any_correct : + forall x m e l, + inbetween_float beta m e x l -> + (e = cexp beta fexp x \/ (l = loc_Exact /\ format x)) -> + round beta fexp rnd x = F2R (Float beta (choice m l) e). +Proof with auto with typeclass_instances. +intros x m e l Hin [He|(Hl,Hf)]. +rewrite He in Hin |- *. +apply inbetween_float_round with (2 := Hin). +exact inbetween_int_valid. +rewrite Hl in Hin. +inversion_clear Hin. +rewrite Hl. +replace (choice m loc_Exact) with m. +rewrite <- H. +apply round_generic... +rewrite <- (Zrnd_IZR rnd m) at 1. +apply inbetween_int_valid. +now constructor. +Qed. + +(** Truncating a triple is sufficient to round a real number. *) + +Theorem round_trunc_any_correct : + forall x m e l, + (0 <= x)%R -> + inbetween_float beta m e x l -> + (e <= fexp (Zdigits beta m + e))%Z \/ l = loc_Exact -> + round beta fexp rnd x = let '(m', e', l') := truncate (m, e, l) in F2R (Float beta (choice m' l') e'). +Proof. +intros x m e l Hx Hl He. +generalize (truncate_correct x m e l Hx Hl He). +destruct (truncate (m, e, l)) as ((m', e'), l'). +intros (H1, H2). +now apply round_any_correct. +Qed. + +Theorem round_trunc_any_correct' : + forall x m e l, + (0 <= x)%R -> + inbetween_float beta m e x l -> + (e <= cexp beta fexp x)%Z \/ l = loc_Exact -> + round beta fexp rnd x = let '(m', e', l') := truncate (m, e, l) in F2R (Float beta (choice m' l') e'). +Proof. +intros x m e l Hx Hl He. +generalize (truncate_correct' x m e l Hx Hl He). +destruct (truncate (m, e, l)) as [[m' e'] l']. +intros [H1 H2]. +now apply round_any_correct. +Qed. + +End round_dir. + +Section round_dir_sign. + +Variable rnd : R -> Z. +Context { valid_rnd : Valid_rnd rnd }. + +Variable choice : bool -> Z -> location -> Z. +Hypothesis inbetween_int_valid : + forall x m l, + inbetween_int m (Rabs x) l -> + rnd x = cond_Zopp (Rlt_bool x 0) (choice (Rlt_bool x 0) m l). + +Theorem round_sign_any_correct : + forall x m e l, + inbetween_float beta m e (Rabs x) l -> + (e = cexp beta fexp x \/ (l = loc_Exact /\ format x)) -> + round beta fexp rnd x = F2R (Float beta (cond_Zopp (Rlt_bool x 0) (choice (Rlt_bool x 0) m l)) e). +Proof with auto with typeclass_instances. +intros x m e l Hin [He|(Hl,Hf)]. +rewrite He in Hin |- *. +apply inbetween_float_round_sign with (2 := Hin). +exact inbetween_int_valid. +rewrite Hl in Hin. +inversion_clear Hin. +rewrite Hl. +replace (choice (Rlt_bool x 0) m loc_Exact) with m. +(* *) +unfold Rabs in H. +destruct (Rcase_abs x) as [Zx|Zx]. +rewrite Rlt_bool_true with (1 := Zx). +simpl. +rewrite F2R_Zopp. +rewrite <- H, Ropp_involutive. +apply round_generic... +rewrite Rlt_bool_false. +simpl. +rewrite <- H. +now apply round_generic. +now apply Rge_le. +(* *) +destruct (Rlt_bool_spec x 0) as [Zx|Zx]. +(* . *) +apply Z.opp_inj. +change (- m = cond_Zopp true (choice true m loc_Exact))%Z. +rewrite <- (Zrnd_IZR rnd (-m)) at 1. +assert (IZR (-m) < 0)%R. +rewrite opp_IZR. +apply Ropp_lt_gt_0_contravar. +apply IZR_lt. +apply gt_0_F2R with beta e. +rewrite <- H. +apply Rabs_pos_lt. +now apply Rlt_not_eq. +rewrite <- Rlt_bool_true with (1 := H0). +apply inbetween_int_valid. +constructor. +rewrite Rabs_left with (1 := H0). +rewrite opp_IZR. +apply Ropp_involutive. +(* . *) +change (m = cond_Zopp false (choice false m loc_Exact))%Z. +rewrite <- (Zrnd_IZR rnd m) at 1. +assert (0 <= IZR m)%R. +apply IZR_le. +apply ge_0_F2R with beta e. +rewrite <- H. +apply Rabs_pos. +rewrite <- Rlt_bool_false with (1 := H0). +apply inbetween_int_valid. +constructor. +now apply Rabs_pos_eq. +Qed. + +(** Truncating a triple is sufficient to round a real number. *) + +Theorem round_trunc_sign_any_correct' : + forall x m e l, + inbetween_float beta m e (Rabs x) l -> + (e <= cexp beta fexp x)%Z \/ l = loc_Exact -> + round beta fexp rnd x = let '(m', e', l') := truncate (m, e, l) in F2R (Float beta (cond_Zopp (Rlt_bool x 0) (choice (Rlt_bool x 0) m' l')) e'). +Proof. +intros x m e l Hl He. +rewrite <- cexp_abs in He. +generalize (truncate_correct' (Rabs x) m e l (Rabs_pos _) Hl He). +destruct (truncate (m, e, l)) as [[m' e'] l']. +intros [H1 H2]. +apply round_sign_any_correct. +exact H1. +destruct H2 as [H2|[H2 H3]]. +left. +now rewrite <- cexp_abs. +right. +apply (conj H2). +now apply generic_format_abs_inv. +Qed. + +Theorem round_trunc_sign_any_correct : + forall x m e l, + inbetween_float beta m e (Rabs x) l -> + (e <= fexp (Zdigits beta m + e))%Z \/ l = loc_Exact -> + round beta fexp rnd x = let '(m', e', l') := truncate (m, e, l) in F2R (Float beta (cond_Zopp (Rlt_bool x 0) (choice (Rlt_bool x 0) m' l')) e'). +Proof. +intros x m e l Hl He. +apply round_trunc_sign_any_correct' with (1 := Hl). +rewrite <- cexp_abs. +apply cexp_inbetween_float_loc_Exact with (2 := Hl) (3 := He). +apply Rabs_pos. +Qed. + +End round_dir_sign. + +(** * Instances of the theorems above, for the usual rounding modes. *) + +Definition round_DN_correct := + round_any_correct _ (fun m _ => m) inbetween_int_DN. + +Definition round_trunc_DN_correct := + round_trunc_any_correct _ (fun m _ => m) inbetween_int_DN. + +Definition round_trunc_DN_correct' := + round_trunc_any_correct' _ (fun m _ => m) inbetween_int_DN. + +Definition round_sign_DN_correct := + round_sign_any_correct _ (fun s m l => cond_incr (round_sign_DN s l) m) inbetween_int_DN_sign. + +Definition round_trunc_sign_DN_correct := + round_trunc_sign_any_correct _ (fun s m l => cond_incr (round_sign_DN s l) m) inbetween_int_DN_sign. + +Definition round_trunc_sign_DN_correct' := + round_trunc_sign_any_correct' _ (fun s m l => cond_incr (round_sign_DN s l) m) inbetween_int_DN_sign. + +Definition round_UP_correct := + round_any_correct _ (fun m l => cond_incr (round_UP l) m) inbetween_int_UP. + +Definition round_trunc_UP_correct := + round_trunc_any_correct _ (fun m l => cond_incr (round_UP l) m) inbetween_int_UP. + +Definition round_trunc_UP_correct' := + round_trunc_any_correct' _ (fun m l => cond_incr (round_UP l) m) inbetween_int_UP. + +Definition round_sign_UP_correct := + round_sign_any_correct _ (fun s m l => cond_incr (round_sign_UP s l) m) inbetween_int_UP_sign. + +Definition round_trunc_sign_UP_correct := + round_trunc_sign_any_correct _ (fun s m l => cond_incr (round_sign_UP s l) m) inbetween_int_UP_sign. + +Definition round_trunc_sign_UP_correct' := + round_trunc_sign_any_correct' _ (fun s m l => cond_incr (round_sign_UP s l) m) inbetween_int_UP_sign. + +Definition round_ZR_correct := + round_any_correct _ (fun m l => cond_incr (round_ZR (Zlt_bool m 0) l) m) inbetween_int_ZR. + +Definition round_trunc_ZR_correct := + round_trunc_any_correct _ (fun m l => cond_incr (round_ZR (Zlt_bool m 0) l) m) inbetween_int_ZR. + +Definition round_trunc_ZR_correct' := + round_trunc_any_correct' _ (fun m l => cond_incr (round_ZR (Zlt_bool m 0) l) m) inbetween_int_ZR. + +Definition round_sign_ZR_correct := + round_sign_any_correct _ (fun s m l => m) inbetween_int_ZR_sign. + +Definition round_trunc_sign_ZR_correct := + round_trunc_sign_any_correct _ (fun s m l => m) inbetween_int_ZR_sign. + +Definition round_trunc_sign_ZR_correct' := + round_trunc_sign_any_correct' _ (fun s m l => m) inbetween_int_ZR_sign. + +Definition round_NE_correct := + round_any_correct _ (fun m l => cond_incr (round_N (negb (Z.even m)) l) m) inbetween_int_NE. + +Definition round_trunc_NE_correct := + round_trunc_any_correct _ (fun m l => cond_incr (round_N (negb (Z.even m)) l) m) inbetween_int_NE. + +Definition round_trunc_NE_correct' := + round_trunc_any_correct' _ (fun m l => cond_incr (round_N (negb (Z.even m)) l) m) inbetween_int_NE. + +Definition round_sign_NE_correct := + round_sign_any_correct _ (fun s m l => cond_incr (round_N (negb (Z.even m)) l) m) inbetween_int_NE_sign. + +Definition round_trunc_sign_NE_correct := + round_trunc_sign_any_correct _ (fun s m l => cond_incr (round_N (negb (Z.even m)) l) m) inbetween_int_NE_sign. + +Definition round_trunc_sign_NE_correct' := + round_trunc_sign_any_correct' _ (fun s m l => cond_incr (round_N (negb (Z.even m)) l) m) inbetween_int_NE_sign. + +Definition round_NA_correct := + round_any_correct _ (fun m l => cond_incr (round_N (Zle_bool 0 m) l) m) inbetween_int_NA. + +Definition round_trunc_NA_correct := + round_trunc_any_correct _ (fun m l => cond_incr (round_N (Zle_bool 0 m) l) m) inbetween_int_NA. + +Definition round_trunc_NA_correct' := + round_trunc_any_correct' _ (fun m l => cond_incr (round_N (Zle_bool 0 m) l) m) inbetween_int_NA. + +Definition round_sign_NA_correct := + round_sign_any_correct _ (fun s m l => cond_incr (round_N true l) m) inbetween_int_NA_sign. + +Definition round_trunc_sign_NA_correct := + round_trunc_sign_any_correct _ (fun s m l => cond_incr (round_N true l) m) inbetween_int_NA_sign. + +Definition round_trunc_sign_NA_correct' := + round_trunc_sign_any_correct' _ (fun s m l => cond_incr (round_N true l) m) inbetween_int_NA_sign. + +End Fcalc_round_fexp. + +(** Specialization of truncate for FIX formats. *) + +Variable emin : Z. + +Definition truncate_FIX t := + let '(m, e, l) := t in + let k := (emin - e)%Z in + if Zlt_bool 0 k then + let p := Zpower beta k in + (Z.div m p, (e + k)%Z, new_location p (Zmod m p) l) + else t. + +Theorem truncate_FIX_correct : + forall x m e l, + inbetween_float beta m e x l -> + (e <= emin)%Z \/ l = loc_Exact -> + let '(m', e', l') := truncate_FIX (m, e, l) in + inbetween_float beta m' e' x l' /\ + (e' = cexp beta (FIX_exp emin) x \/ (l' = loc_Exact /\ generic_format beta (FIX_exp emin) x)). +Proof. +intros x m e l H1 H2. +unfold truncate_FIX. +set (k := (emin - e)%Z). +set (p := Zpower beta k). +unfold cexp, FIX_exp. +generalize (Zlt_cases 0 k). +case (Zlt_bool 0 k) ; intros Hk. +(* shift *) +split. +now apply inbetween_float_new_location. +clear H2. +left. +unfold k. +ring. +(* no shift *) +split. +exact H1. +unfold k in Hk. +destruct H2 as [H2|H2]. +left. +omega. +right. +split. +exact H2. +rewrite H2 in H1. +inversion_clear H1. +rewrite H. +apply generic_format_F2R. +unfold cexp. +omega. +Qed. + +End Fcalc_round. diff --git a/flocq/Calc/Sqrt.v b/flocq/Calc/Sqrt.v new file mode 100644 index 00000000..8843d21e --- /dev/null +++ b/flocq/Calc/Sqrt.v @@ -0,0 +1,201 @@ +(** +This file is part of the Flocq formalization of floating-point +arithmetic in Coq: http://flocq.gforge.inria.fr/ + +Copyright (C) 2010-2018 Sylvie Boldo +#
# +Copyright (C) 2010-2018 Guillaume Melquiond + +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. +*) + +(** * Helper functions and theorems for computing the rounded square root of a floating-point number. *) + +Require Import Raux Defs Digits Generic_fmt Float_prop Bracket. + +Set Implicit Arguments. +Set Strongly Strict Implicit. + +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. + +The algorithm performs the following steps: +- Shift the mantissa so that it has at least 2p-1 digits; + shift it one digit more if the new exponent is not even. +- Compute the square root s (at least p digits) of the new + mantissa, and its remainder r. +- Compute the position according to the remainder: + -- r == 0 => Eq, + -- r <= s => Lo, + -- r >= s => Up. + +Complexity is fine as long as p1 <= 2p-1. +*) + +Lemma mag_sqrt_F2R : + forall m1 e1, + (0 < m1)%Z -> + mag beta (sqrt (F2R (Float beta m1 e1))) = Z.div2 (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, l). + +Theorem Fsqrt_core_correct : + 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 m1 e1 e Hm1 He. +unfold Fsqrt_core. +set (m' := Zmult _ _). +assert (0 <= m')%Z as Hm'. +{ apply Z.mul_nonneg_nonneg. + now apply Zlt_le_weak. + 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]. +rewrite Hf. +unfold inbetween_float, F2R. simpl Fnum. +apply inbetween_mult_compat. +apply bpow_gt_0. +rewrite Hq. +case Zeq_bool_spec ; intros Hr'. +(* .. r = 0 *) +rewrite Hr', Zplus_0_r, mult_IZR. +fold (Rsqr (IZR q)). +rewrite sqrt_Rsqr. +now constructor. +apply IZR_le. +clear -Hr ; omega. +(* .. r <> 0 *) +constructor. +split. +(* ... bounds *) +apply Rle_lt_trans with (sqrt (IZR (q * q))). +rewrite mult_IZR. +fold (Rsqr (IZR q)). +rewrite sqrt_Rsqr. +apply Rle_refl. +apply IZR_le. +clear -Hr ; omega. +apply sqrt_lt_1. +rewrite mult_IZR. +apply Rle_0_sqr. +rewrite <- Hq. +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. +now apply IZR_le. +rewrite mult_IZR. +apply Rle_0_sqr. +apply IZR_lt. +ring_simplify. +omega. +rewrite mult_IZR. +fold (Rsqr (IZR (q + 1))). +rewrite sqrt_Rsqr. +apply Rle_refl. +apply IZR_le. +clear -Hr ; omega. +(* ... location *) +rewrite Rcompare_half_r. +generalize (Rcompare_sqr (2 * sqrt (IZR (q * q + r))) (IZR q + IZR (q + 1))). +rewrite 2!Rabs_pos_eq. +intros <-. +replace ((2 * sqrt (IZR (q * q + r))) * (2 * sqrt (IZR (q * q + r))))%R + with (4 * Rsqr (sqrt (IZR (q * q + r))))%R by (unfold Rsqr ; ring). +rewrite Rsqr_sqrt. +rewrite <- plus_IZR, <- 2!mult_IZR. +rewrite Rcompare_IZR. +replace ((q + (q + 1)) * (q + (q + 1)))%Z with (4 * (q * q) + 4 * q + 1)%Z by ring. +generalize (Zle_cases r q). +case (Zle_bool r q) ; intros Hr''. +change (4 * (q * q + r) < 4 * (q * q) + 4 * q + 1)%Z. +omega. +change (4 * (q * q + r) > 4 * (q * q) + 4 * q + 1)%Z. +omega. +rewrite <- Hq. +now apply IZR_le. +rewrite <- plus_IZR. +apply IZR_le. +clear -Hr ; omega. +apply Rmult_le_pos. +now apply IZR_le. +apply sqrt_ge_0. +Qed. + +Definition Fsqrt (x : float beta) := + let (m1, e1) := x in + let e' := (Zdigits beta m1 + e1 + 1)%Z in + let e := Z.min (fexp (Z.div2 e')) (Z.div2 e1) in + let '(m, l) := Fsqrt_core m1 e1 e in + (m, e, l). + +Theorem Fsqrt_correct : + forall x, + (0 < F2R x)%R -> + let '(m, e, l) := Fsqrt x in + (e <= cexp beta fexp (sqrt (F2R x)))%Z /\ + inbetween_float beta m e (sqrt (F2R x)) l. +Proof. +intros [m1 e1] Hm1. +apply gt_0_F2R in Hm1. +unfold Fsqrt. +set (e := Z.min _ _). +assert (2 * e <= e1)%Z as He. +{ assert (e <= Z.div2 e1)%Z by apply Z.le_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 Z.le_trans with (1 := Z.le_min_l _ _). +unfold cexp. +rewrite (mag_sqrt_F2R m1 e1 Hm1). +apply Z.le_refl. +Qed. + +End Fcalc_sqrt. -- cgit