diff options
Diffstat (limited to 'flocq/Calc')
| -rw-r--r-- | flocq/Calc/Fcalc_bracket.v | 688 | ||||
| -rw-r--r-- | flocq/Calc/Fcalc_digits.v | 382 | ||||
| -rw-r--r-- | flocq/Calc/Fcalc_div.v | 165 | ||||
| -rw-r--r-- | flocq/Calc/Fcalc_ops.v | 161 | ||||
| -rw-r--r-- | flocq/Calc/Fcalc_round.v | 1093 | ||||
| -rw-r--r-- | flocq/Calc/Fcalc_sqrt.v | 346 |
6 files changed, 2835 insertions, 0 deletions
diff --git a/flocq/Calc/Fcalc_bracket.v b/flocq/Calc/Fcalc_bracket.v new file mode 100644 index 0000000..dd4bd97 --- /dev/null +++ b/flocq/Calc/Fcalc_bracket.v @@ -0,0 +1,688 @@ +(** +This file is part of the Flocq formalization of floating-point +arithmetic in Coq: http://flocq.gforge.inria.fr/ + +Copyright (C) 2010-2011 Sylvie Boldo +#<br /># +Copyright (C) 2010-2011 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_r 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 new file mode 100644 index 0000000..6210bac --- /dev/null +++ b/flocq/Calc/Fcalc_digits.v @@ -0,0 +1,382 @@ +(** +This file is part of the Flocq formalization of floating-point +arithmetic in Coq: http://flocq.gforge.inria.fr/ + +Copyright (C) 2010-2011 Sylvie Boldo +#<br /># +Copyright (C) 2010-2011 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. + +Theorem Zdigits_mult_Zpower : + forall m e, + m <> Z0 -> (0 <= e)%Z -> + Zdigits beta (m * Zpower beta e) = (Zdigits beta m + e)%Z. +Proof. +intros m e Hm He. +rewrite <- ln_beta_F2R_Zdigits with (1 := Hm). +rewrite Zdigits_ln_beta. +rewrite Z2R_mult. +now rewrite Z2R_Zpower with (1 := He). +contradict Hm. +apply Zmult_integral_l with (2 := Hm). +apply neq_Z2R. +rewrite Z2R_Zpower with (1 := He). +apply Rgt_not_eq. +apply bpow_gt_0. +Qed. + +Theorem Zdigits_Zpower : + forall e, + (0 <= e)%Z -> + Zdigits beta (Zpower beta e) = (e + 1)%Z. +Proof. +intros e He. +rewrite <- (Zmult_1_l (Zpower beta e)). +rewrite Zdigits_mult_Zpower ; try easy. +apply Zplus_comm. +Qed. + +Theorem Zdigits_le : + forall x y, + (0 <= x)%Z -> (x <= y)%Z -> + (Zdigits beta x <= Zdigits beta y)%Z. +Proof. +intros x y Hx Hxy. +case (Z_lt_le_dec 0 x). +clear Hx. intros Hx. +rewrite 2!Zdigits_ln_beta. +apply ln_beta_le. +now apply (Z2R_lt 0). +now apply Z2R_le. +apply Zgt_not_eq. +now apply Zlt_le_trans with x. +now apply Zgt_not_eq. +intros Hx'. +rewrite (Zle_antisym _ _ Hx' Hx). +apply Zdigits_ge_0. +Qed. + +Theorem lt_Zdigits : + forall x y, + (0 <= y)%Z -> + (Zdigits beta x < Zdigits beta y)%Z -> + (x < y)%Z. +Proof. +intros x y Hy. +cut (y <= x -> Zdigits beta y <= Zdigits beta x)%Z. omega. +now apply Zdigits_le. +Qed. + +Theorem Zpower_le_Zdigits : + forall e x, + (e < Zdigits beta x)%Z -> + (Zpower beta e <= Zabs x)%Z. +Proof. +intros e x Hex. +destruct (Zdigits_correct beta x) as (H1,H2). +apply Zle_trans with (2 := H1). +apply Zpower_le. +clear -Hex ; omega. +Qed. + +Theorem Zdigits_le_Zpower : + forall e x, + (Zabs x < Zpower beta e)%Z -> + (Zdigits beta x <= e)%Z. +Proof. +intros e x. +generalize (Zpower_le_Zdigits e x). +omega. +Qed. + +Theorem Zpower_gt_Zdigits : + forall e x, + (Zdigits beta x <= e)%Z -> + (Zabs x < Zpower beta e)%Z. +Proof. +intros e x Hex. +destruct (Zdigits_correct beta x) as (H1,H2). +apply Zlt_le_trans with (1 := H2). +now apply Zpower_le. +Qed. + +Theorem Zdigits_gt_Zpower : + forall e x, + (Zpower beta e <= Zabs x)%Z -> + (e < Zdigits beta x)%Z. +Proof. +intros e x Hex. +generalize (Zpower_gt_Zdigits e x). +omega. +Qed. + +(** Characterizes the number digits of a product. + +This strong version is needed for proofs of division and square root +algorithms, since they involve operation remainders. +*) + +Theorem Zdigits_mult_strong : + forall x y, + (0 <= x)%Z -> (0 <= y)%Z -> + (Zdigits beta (x + y + x * y) <= Zdigits beta x + Zdigits beta y)%Z. +Proof. +intros x y Hx Hy. +case (Z_lt_le_dec 0 x). +clear Hx. intros Hx. +case (Z_lt_le_dec 0 y). +clear Hy. intros Hy. +(* . *) +assert (Hxy: (0 < Z2R (x + y + x * y))%R). +apply (Z2R_lt 0). +change Z0 with (0 + 0 + 0)%Z. +apply Zplus_lt_compat. +now apply Zplus_lt_compat. +now apply Zmult_lt_0_compat. +rewrite 3!Zdigits_ln_beta ; try now (apply sym_not_eq ; apply Zlt_not_eq). +apply ln_beta_le_bpow with (1 := Rgt_not_eq _ _ Hxy). +rewrite Rabs_pos_eq with (1 := Rlt_le _ _ Hxy). +destruct (ln_beta beta (Z2R x)) as (ex, Hex). simpl. +specialize (Hex (Rgt_not_eq _ _ (Z2R_lt _ _ Hx))). +destruct (ln_beta beta (Z2R y)) as (ey, Hey). simpl. +specialize (Hey (Rgt_not_eq _ _ (Z2R_lt _ _ Hy))). +apply Rlt_le_trans with (Z2R (x + 1) * Z2R (y + 1))%R. +rewrite <- Z2R_mult. +apply Z2R_lt. +apply Zplus_lt_reg_r with (- (x + y + x * y + 1))%Z. +now ring_simplify. +rewrite bpow_plus. +apply Rmult_le_compat ; try (apply (Z2R_le 0) ; omega). +rewrite <- (Rmult_1_r (Z2R (x + 1))). +change (F2R (Float beta (x + 1) 0) <= bpow ex)%R. +apply F2R_p1_le_bpow. +exact Hx. +unfold F2R. rewrite Rmult_1_r. +apply Rle_lt_trans with (Rabs (Z2R x)). +apply RRle_abs. +apply Hex. +rewrite <- (Rmult_1_r (Z2R (y + 1))). +change (F2R (Float beta (y + 1) 0) <= bpow ey)%R. +apply F2R_p1_le_bpow. +exact Hy. +unfold F2R. rewrite Rmult_1_r. +apply Rle_lt_trans with (Rabs (Z2R y)). +apply RRle_abs. +apply Hey. +apply neq_Z2R. +now apply Rgt_not_eq. +(* . *) +intros Hy'. +rewrite (Zle_antisym _ _ Hy' Hy). +rewrite Zmult_0_r, 3!Zplus_0_r. +apply Zle_refl. +intros Hx'. +rewrite (Zle_antisym _ _ Hx' Hx). +rewrite Zmult_0_l, Zplus_0_r, 2!Zplus_0_l. +apply Zle_refl. +Qed. + +Theorem Zdigits_mult : + forall x y, + (Zdigits beta (x * y) <= Zdigits beta x + Zdigits beta y)%Z. +Proof. +intros x y. +rewrite <- Zdigits_abs. +rewrite <- (Zdigits_abs _ x). +rewrite <- (Zdigits_abs _ y). +apply Zle_trans with (Zdigits beta (Zabs x + Zabs y + Zabs x * Zabs y)). +apply Zdigits_le. +apply Zabs_pos. +rewrite Zabs_Zmult. +generalize (Zabs_pos x) (Zabs_pos y). +omega. +apply Zdigits_mult_strong ; apply Zabs_pos. +Qed. + +Theorem Zdigits_mult_ge : + forall x y, + (x <> 0)%Z -> (y <> 0)%Z -> + (Zdigits beta x + Zdigits beta y - 1 <= Zdigits beta (x * y))%Z. +Proof. +intros x y Zx Zy. +cut ((Zdigits beta x - 1) + (Zdigits beta y - 1) < Zdigits beta (x * y))%Z. omega. +apply Zdigits_gt_Zpower. +rewrite Zabs_Zmult. +rewrite Zpower_exp. +apply Zmult_le_compat. +apply Zpower_le_Zdigits. +apply Zlt_pred. +apply Zpower_le_Zdigits. +apply Zlt_pred. +apply Zpower_ge_0. +apply Zpower_ge_0. +generalize (Zdigits_gt_0 beta x). omega. +generalize (Zdigits_gt_0 beta y). omega. +Qed. + +Theorem Zdigits_div_Zpower : + forall m e, + (0 <= m)%Z -> + (0 <= e <= Zdigits beta m)%Z -> + Zdigits beta (m / Zpower beta e) = (Zdigits beta m - e)%Z. +Proof. +intros m e Hm He. +destruct (Zle_lt_or_eq 0 m Hm) as [Hm'|Hm']. +(* *) +destruct (Zle_lt_or_eq _ _ (proj2 He)) as [He'|He']. +(* . *) +unfold Zminus. +rewrite <- ln_beta_F2R_Zdigits. +2: now apply Zgt_not_eq. +assert (Hp: (0 < Zpower beta e)%Z). +apply lt_Z2R. +rewrite Z2R_Zpower with (1 := proj1 He). +apply bpow_gt_0. +rewrite Zdigits_ln_beta. +rewrite <- Zfloor_div with (1 := Zgt_not_eq _ _ Hp). +rewrite Z2R_Zpower with (1 := proj1 He). +unfold Rdiv. +rewrite <- bpow_opp. +change (Z2R m * bpow (-e))%R with (F2R (Float beta m (-e))). +destruct (ln_beta beta (F2R (Float beta m (-e)))) as (e', E'). +simpl. +specialize (E' (Rgt_not_eq _ _ (F2R_gt_0_compat beta (Float beta m (-e)) Hm'))). +apply ln_beta_unique. +rewrite Rabs_pos_eq in E'. +2: now apply F2R_ge_0_compat. +rewrite Rabs_pos_eq. +split. +assert (H': (0 <= e' - 1)%Z). +(* .. *) +cut (1 <= e')%Z. omega. +apply bpow_lt_bpow with beta. +apply Rle_lt_trans with (2 := proj2 E'). +simpl. +rewrite <- (Rinv_r (bpow e)). +rewrite <- bpow_opp. +unfold F2R. simpl. +apply Rmult_le_compat_r. +apply bpow_ge_0. +rewrite <- Z2R_Zpower with (1 := proj1 He). +apply Z2R_le. +rewrite <- Zabs_eq with (1 := Hm). +now apply Zpower_le_Zdigits. +apply Rgt_not_eq. +apply bpow_gt_0. +(* .. *) +rewrite <- Z2R_Zpower with (1 := H'). +apply Z2R_le. +apply Zfloor_lub. +rewrite Z2R_Zpower with (1 := H'). +apply E'. +apply Rle_lt_trans with (2 := proj2 E'). +apply Zfloor_lb. +apply (Z2R_le 0). +apply Zfloor_lub. +now apply F2R_ge_0_compat. +apply Zgt_not_eq. +apply Zlt_le_trans with (beta^e / beta^e)%Z. +rewrite Z_div_same_full. +apply refl_equal. +now apply Zgt_not_eq. +apply Z_div_le. +now apply Zlt_gt. +rewrite <- Zabs_eq with (1 := Hm). +now apply Zpower_le_Zdigits. +(* . *) +unfold Zminus. +rewrite He', Zplus_opp_r. +rewrite Zdiv_small. +apply refl_equal. +apply conj with (1 := Hm). +pattern m at 1 ; rewrite <- Zabs_eq with (1 := Hm). +apply Zpower_gt_Zdigits. +apply Zle_refl. +(* *) +revert He. +rewrite <- Hm'. +intros (H1, H2). +simpl. +now rewrite (Zle_antisym _ _ H2 H1). +Qed. + +End Fcalc_digits. + +Definition radix2 := Build_radix 2 (refl_equal _). + +Theorem Z_of_nat_S_digits2_Pnat : + forall m : positive, + Z_of_nat (S (digits2_Pnat m)) = Zdigits radix2 (Zpos m). +Proof. +intros m. +rewrite Zdigits_ln_beta. 2: easy. +apply sym_eq. +apply ln_beta_unique. +generalize (digits2_Pnat m) (digits2_Pnat_correct m). +intros d H. +simpl in H. +replace (Z_of_nat (S d) - 1)%Z with (Z_of_nat d). +rewrite <- Z2R_abs. +rewrite <- 2!Z2R_Zpower_nat. +split. +now apply Z2R_le. +now apply Z2R_lt. +rewrite inj_S. +apply Zpred_succ. +Qed. + diff --git a/flocq/Calc/Fcalc_div.v b/flocq/Calc/Fcalc_div.v new file mode 100644 index 0000000..6594e55 --- /dev/null +++ b/flocq/Calc/Fcalc_div.v @@ -0,0 +1,165 @@ +(** +This file is part of the Flocq formalization of floating-point +arithmetic in Coq: http://flocq.gforge.inria.fr/ + +Copyright (C) 2010-2011 Sylvie Boldo +#<br /># +Copyright (C) 2010-2011 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 new file mode 100644 index 0000000..15bb211 --- /dev/null +++ b/flocq/Calc/Fcalc_ops.v @@ -0,0 +1,161 @@ +(** +This file is part of the Flocq formalization of floating-point +arithmetic in Coq: http://flocq.gforge.inria.fr/ + +Copyright (C) 2010-2011 Sylvie Boldo +#<br /># +Copyright (C) 2010-2011 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). + +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 (Float beta m1 e) /\ F2R f2 = F2R (Float beta 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) := + let '(Float m1 e1) := f1 in + Float beta (-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) := + let '(Float m1 e1) := f1 in + Float beta (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) := + let '(m1, m2 ,e) := Falign f1 f2 in + Float beta (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 beta m1 e) (Float beta m2 e) = Float beta (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 beta m1 e) (Float beta m2 e) = Float beta (m1 - m2) e. +Proof. +intros m1 m2 e. +unfold Fminus. +apply Fplus_same_exp. +Qed. + +Definition Fmult (f1 f2 : float beta) := + let '(Float m1 e1) := f1 in + let '(Float m2 e2) := f2 in + Float beta (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 new file mode 100644 index 0000000..3d31aea --- /dev/null +++ b/flocq/Calc/Fcalc_round.v @@ -0,0 +1,1093 @@ +(** +This file is part of the Flocq formalization of floating-point +arithmetic in Coq: http://flocq.gforge.inria.fr/ + +Copyright (C) 2010-2011 Sylvie Boldo +#<br /># +Copyright (C) 2010-2011 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. +refine (let H := _ in conj _ H). +unfold k. ring. +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 new file mode 100644 index 0000000..e6fe74b --- /dev/null +++ b/flocq/Calc/Fcalc_sqrt.v @@ -0,0 +1,346 @@ +(** +This file is part of the Flocq formalization of floating-point +arithmetic in Coq: http://flocq.gforge.inria.fr/ + +Copyright (C) 2010-2011 Sylvie Boldo +#<br /># +Copyright (C) 2010-2011 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. + +Fixpoint Zsqrt_aux (p : positive) : Z * Z := + match p with + | xH => (1, 0)%Z + | xO xH => (1, 1)%Z + | xI xH => (1, 2)%Z + | xO (xO p) => + let (q, r) := Zsqrt_aux p in + let r' := (4 * r)%Z in + let d := (r' - (4 * q + 1))%Z in + if Zlt_bool d 0 then (2 * q, r')%Z else (2 * q + 1, d)%Z + | xO (xI p) => + let (q, r) := Zsqrt_aux p in + let r' := (4 * r + 2)%Z in + let d := (r' - (4 * q + 1))%Z in + if Zlt_bool d 0 then (2 * q, r')%Z else (2 * q + 1, d)%Z + | xI (xO p) => + let (q, r) := Zsqrt_aux p in + let r' := (4 * r + 1)%Z in + let d := (r' - (4 * q + 1))%Z in + if Zlt_bool d 0 then (2 * q, r')%Z else (2 * q + 1, d)%Z + | xI (xI p) => + let (q, r) := Zsqrt_aux p in + let r' := (4 * r + 3)%Z in + let d := (r' - (4 * q + 1))%Z in + if Zlt_bool d 0 then (2 * q, r')%Z else (2 * q + 1, d)%Z + end. + +Lemma Zsqrt_ind : + forall P : positive -> Prop, + P xH -> P (xO xH) -> P (xI xH) -> + ( forall p, P p -> P (xO (xO p)) /\ P (xO (xI p)) /\ P (xI (xO p)) /\ P (xI (xI p)) ) -> + forall p, P p. +Proof. +intros P H1 H2 H3 Hp. +fix 1. +intros [[p|p|]|[p|p|]|]. +refine (proj2 (proj2 (proj2 (Hp p _)))). +apply Zsqrt_ind. +refine (proj1 (proj2 (proj2 (Hp p _)))). +apply Zsqrt_ind. +exact H3. +refine (proj1 (proj2 (Hp p _))). +apply Zsqrt_ind. +refine (proj1 (Hp p _)). +apply Zsqrt_ind. +exact H2. +exact H1. +Qed. + +Lemma Zsqrt_aux_correct : + forall p, + let (q, r) := Zsqrt_aux p in + Zpos p = (q * q + r)%Z /\ (0 <= r <= 2 * q)%Z. +Proof. +intros p. +elim p using Zsqrt_ind ; clear p. +now repeat split. +now repeat split. +now repeat split. +intros p. +Opaque Zmult. simpl. Transparent Zmult. +destruct (Zsqrt_aux p) as (q, r). +intros (Hq, Hr). +change (Zpos p~0~0) with (4 * Zpos p)%Z. +change (Zpos p~0~1) with (4 * Zpos p + 1)%Z. +change (Zpos p~1~0) with (4 * Zpos p + 2)%Z. +change (Zpos p~1~1) with (4 * Zpos p + 3)%Z. +rewrite Hq. clear Hq. +repeat split. +generalize (Zlt_cases (4 * r - (4 * q + 1)) 0). +case Zlt_bool ; ( split ; [ ring | omega ] ). +generalize (Zlt_cases (4 * r + 2 - (4 * q + 1)) 0). +case Zlt_bool ; ( split ; [ ring | omega ] ). +generalize (Zlt_cases (4 * r + 1 - (4 * q + 1)) 0). +case Zlt_bool ; ( split ; [ ring | omega ] ). +generalize (Zlt_cases (4 * r + 3 - (4 * q + 1)) 0). +case Zlt_bool ; ( split ; [ ring | omega ] ). +Qed. + +(** Computes the integer square root and its remainder, but + without carrying a proof, contrarily to the operation of + the standard libary. *) + +Definition Zsqrt p := + match p with + | Zpos p => Zsqrt_aux p + | _ => (0, 0)%Z + end. + +Theorem Zsqrt_correct : + forall x, + (0 <= x)%Z -> + let (q, r) := Zsqrt x in + x = (q * q + r)%Z /\ (0 <= r <= 2 * q)%Z. +Proof. +unfold Zsqrt. +intros [|p|p] Hx. +now repeat split. +apply Zsqrt_aux_correct. +now elim Hx. +Qed. + +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) := Zsqrt 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 (Zsqrt_correct m' (Zlt_le_weak _ _ Hs3)). +destruct (Zsqrt 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. |