Fcalc_digits.v 7.38 KB
 Guillaume Melquiond committed May 17, 2010 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 ``````Require Import Fcore. Section Fcalc_digits. Variable beta : radix. Notation bpow e := (bpow beta e). Fixpoint digits2_Pnat (n : positive) : nat := match n with | xH => O | xO p => S (digits2_Pnat p) | xI p => S (digits2_Pnat p) end. Theorem digits2_Pnat_correct : forall n, let d := digits2_Pnat n in (Zpower_nat 2 d <= Zpos n < Zpower_nat 2 (S d))%Z. Proof. intros n d. unfold d. clear. assert (Hp: forall m, (Zpower_nat 2 (S m) = 2 * Zpower_nat 2 m)%Z) by easy. induction n ; simpl. rewrite Zpos_xI, 2!Hp. omega. rewrite (Zpos_xO n), 2!Hp. omega. now split. Qed. Section digits_aux. Variable p : Z. Hypothesis Hp : (0 <= p)%Z. Fixpoint digits_aux (nb pow : Z) (n : nat) { struct n } : Z := match n with | O => nb `````` Guillaume Melquiond committed Oct 01, 2010 38 `````` | S n => if Zlt_bool p pow then nb else digits_aux (nb + 1) (Zmult beta pow) n `````` Guillaume Melquiond committed May 17, 2010 39 40 41 42 43 `````` end. Lemma digits_aux_invariant : forall k n nb, (0 <= nb)%Z -> `````` Guillaume Melquiond committed Oct 01, 2010 44 45 46 `````` (Zpower beta (nb + Z_of_nat k - 1) <= p)%Z -> digits_aux (nb + Z_of_nat k) (Zpower beta (nb + Z_of_nat k)) n = digits_aux nb (Zpower beta nb) (n + k). `````` Guillaume Melquiond committed May 17, 2010 47 48 49 50 51 52 53 54 55 56 ``````Proof. induction k ; intros n nb Hnb. now rewrite Zplus_0_r, plus_0_r. rewrite inj_S. unfold Zsucc. intros H. rewrite (Zplus_comm _ 1), Zplus_assoc. rewrite IHk. rewrite <- plus_n_Sm. simpl. `````` Guillaume Melquiond committed Oct 01, 2010 57 ``````generalize (Zlt_cases p (Zpower beta nb)). `````` Guillaume Melquiond committed May 17, 2010 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 ``````case Zlt_bool ; intros Hpp. elim (Zlt_irrefl p). apply Zlt_le_trans with (1 := Hpp). apply Zle_trans with (2 := H). replace (nb + (Z_of_nat k + 1) - 1)%Z with (nb + Z_of_nat k)%Z by ring. apply le_Z2R. rewrite Z2R_Zpower with (1 := Hnb). rewrite Z2R_Zpower. apply -> bpow_le. omega. omega. rewrite Zpower_exp. unfold Zpower at 2, Zpower_pos, iter_pos. rewrite Zmult_1_r. now rewrite Zmult_comm. now apply Zle_ge. easy. omega. now rewrite <- Zplus_assoc, (Zplus_comm 1). Qed. End digits_aux. Definition digits n := match n with | Z0 => Z0 `````` Guillaume Melquiond committed Oct 01, 2010 84 85 `````` | Zneg p => digits_aux (Zpos p) 1 beta (digits2_Pnat p) | Zpos p => digits_aux n 1 beta (digits2_Pnat p) `````` Guillaume Melquiond committed May 17, 2010 86 87 88 89 90 91 92 93 94 95 96 `````` end. Theorem digits_abs : forall n, digits (Zabs n) = digits n. Proof. now intros [|n|n]. Qed. Theorem digits_ln_beta : forall n, n <> Z0 -> `````` Guillaume Melquiond committed Oct 01, 2010 97 `````` digits n = ln_beta beta (Z2R n). `````` Guillaume Melquiond committed May 17, 2010 98 99 100 101 102 ``````Proof. intros n Hn. destruct (ln_beta beta (Z2R n)) as (e, He). simpl. specialize (He (Z2R_neq _ _ Hn)). `````` Guillaume Melquiond committed Oct 01, 2010 103 ``````rewrite <- Z2R_abs in He. `````` Guillaume Melquiond committed May 17, 2010 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 ``````assert (Hn': (0 < Zabs n)%Z). destruct n ; try easy. now elim Hn. rewrite <- digits_abs. destruct (Zabs n) as [|p|p] ; try (now elim Hn'). clear n Hn Hn'. simpl. assert (He1: (0 <= e - 1)%Z). apply Zlt_0_le_0_pred. apply <- bpow_lt. apply Rle_lt_trans with (2 := proj2 He). apply (Z2R_le 1). now apply (Zlt_le_succ 0). assert (He2: (Z_of_nat (Zabs_nat (e - 1)) = e - 1)%Z). rewrite inj_Zabs_nat. now apply Zabs_eq. `````` Guillaume Melquiond committed Oct 01, 2010 120 ``````replace (radix_val beta) with (Zpower beta 1). `````` Guillaume Melquiond committed May 17, 2010 121 122 123 124 125 126 127 ``````replace (digits2_Pnat p) with (digits2_Pnat p - Zabs_nat (e - 1) + Zabs_nat (e - 1)). rewrite <- digits_aux_invariant. rewrite He2. ring_simplify (1 + (e - 1))%Z. destruct (digits2_Pnat p - Zabs_nat (e - 1)) as [|n]. easy. simpl. `````` Guillaume Melquiond committed Oct 01, 2010 128 ``````generalize (Zlt_cases (Zpos p) (Zpower beta e)). `````` Guillaume Melquiond committed May 17, 2010 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 ``````case Zlt_bool ; intros Hpp. easy. elim Rlt_not_le with (1 := proj2 He). rewrite <- Z2R_Zpower. apply Z2R_le. now apply Zge_le. omega. easy. rewrite He2. ring_simplify (1 + (e - 1) - 1)%Z. apply le_Z2R. now rewrite Z2R_Zpower. rewrite plus_comm. apply le_plus_minus_r. apply inj_le_rev. rewrite He2. cut (e - 1 < Z_of_nat (S (digits2_Pnat p)))%Z. rewrite inj_S. omega. apply <- bpow_lt. apply Rle_lt_trans with (1 := proj1 He). `````` Guillaume Melquiond committed Aug 20, 2010 150 ``````rewrite <- Z2R_Zpower_nat. `````` Guillaume Melquiond committed May 17, 2010 151 ``````apply Z2R_lt. `````` Guillaume Melquiond committed Aug 20, 2010 152 ``````apply Zlt_le_trans with (Zpower_nat 2 (S (digits2_Pnat p))). `````` Guillaume Melquiond committed May 17, 2010 153 154 155 156 ``````exact (proj2 (digits2_Pnat_correct p)). clear. induction (S (digits2_Pnat p)). easy. `````` Guillaume Melquiond committed Oct 01, 2010 157 ``````change (2 * Zpower_nat 2 n <= beta * Zpower_nat beta n)%Z. `````` Guillaume Melquiond committed May 17, 2010 158 ``````apply Zmult_le_compat ; try easy. `````` Guillaume Melquiond committed Aug 19, 2010 159 ``````apply Zle_bool_imp_le. `````` Guillaume Melquiond committed May 17, 2010 160 161 162 163 164 ``````apply beta. now apply Zpower_NR0. apply Zmult_1_r. Qed. `````` Guillaume Melquiond committed May 17, 2010 165 166 167 168 169 170 171 172 173 174 175 176 ``````Theorem digits_ge_0 : forall n, (0 <= digits n)%Z. Proof. intros n. destruct (Z_eq_dec n 0) as [H|H]. now rewrite H. rewrite digits_ln_beta with (1 := H). destruct ln_beta as (e, He). simpl. apply <- bpow_le. apply Rlt_le. apply Rle_lt_trans with (Rabs (Z2R n)). simpl. `````` Guillaume Melquiond committed Oct 01, 2010 177 ``````rewrite <- Z2R_abs. `````` Guillaume Melquiond committed May 17, 2010 178 179 180 181 182 183 184 185 186 187 ``````apply (Z2R_le 1). apply (Zlt_le_succ 0). revert H. case n ; try easy. intros H. now elim H. apply He. now apply (Z2R_neq _ 0). Qed. `````` Guillaume Melquiond committed May 17, 2010 188 189 190 ``````Theorem digits_shift : forall m e, m <> Z0 -> (0 <= e)%Z -> `````` Guillaume Melquiond committed Oct 01, 2010 191 `````` digits (m * Zpower beta e) = (digits m + e)%Z. `````` Guillaume Melquiond committed May 17, 2010 192 193 194 ``````Proof. intros m e Hm He. rewrite 2!digits_ln_beta. `````` Guillaume Melquiond committed Oct 01, 2010 195 ``````rewrite Z2R_mult. `````` Guillaume Melquiond committed May 17, 2010 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 ``````rewrite Z2R_Zpower with (1 := He). change (Z2R m * bpow e)%R with (F2R (Float beta m e)). apply ln_beta_F2R. exact Hm. exact Hm. 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 digits_le : forall x y, `````` Guillaume Melquiond committed May 17, 2010 211 `````` (0 <= x)%Z -> (x <= y)%Z -> `````` Guillaume Melquiond committed May 17, 2010 212 213 214 `````` (digits x <= digits y)%Z. Proof. intros x y Hx Hxy. `````` Guillaume Melquiond committed May 17, 2010 215 216 ``````case (Z_lt_le_dec 0 x). clear Hx. intros Hx. `````` Guillaume Melquiond committed May 17, 2010 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 ``````assert (Hy: (y <> 0)%Z). apply sym_not_eq. apply Zlt_not_eq. now apply Zlt_le_trans with x. rewrite 2!digits_ln_beta. 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 (Z2R_neq _ _ Hy)). eapply bpow_lt_bpow. apply Rle_lt_trans with (1 := proj1 Hex). apply Rle_lt_trans with (Rabs (Z2R y)). rewrite Rabs_pos_eq. apply Rle_trans with (Z2R y). now apply Z2R_le. apply RRle_abs. apply (Z2R_le 0). now apply Zlt_le_weak. apply Hey. exact Hy. apply sym_not_eq. now apply Zlt_not_eq. `````` Guillaume Melquiond committed May 17, 2010 239 240 241 ``````intros Hx'. rewrite (Zle_antisym _ _ Hx' Hx). apply digits_ge_0. `````` Guillaume Melquiond committed May 17, 2010 242 243 244 245 ``````Qed. Theorem digits_lt : forall x y, `````` Guillaume Melquiond committed May 17, 2010 246 `````` (0 <= y)%Z -> `````` Guillaume Melquiond committed May 17, 2010 247 248 249 250 251 252 253 254 255 256 `````` (digits x < digits y)%Z -> (x < y)%Z. Proof. intros x y Hy. cut (y <= x -> digits y <= digits x)%Z. omega. now apply digits_le. Qed. Theorem digits_mult_strong : forall x y, `````` Guillaume Melquiond committed May 17, 2010 257 `````` (0 <= x)%Z -> (0 <= y)%Z -> `````` Guillaume Melquiond committed May 17, 2010 258 259 260 `````` (digits (x + y + x * y) <= digits x + digits y)%Z. Proof. intros x y Hx Hy. `````` Guillaume Melquiond committed May 17, 2010 261 262 263 264 265 ``````case (Z_lt_le_dec 0 x). clear Hx. intros Hx. case (Z_lt_le_dec 0 y). clear Hy. intros Hy. (* . *) `````` Guillaume Melquiond committed May 17, 2010 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 ``````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!digits_ln_beta ; try now (apply sym_not_eq ; apply Zlt_not_eq). destruct (ln_beta beta (Z2R (x + y + x * y))) as (exy, Hexy). simpl. specialize (Hexy (Rgt_not_eq _ _ 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))). eapply bpow_lt_bpow. apply Rlt_le_trans with (Z2R (x + 1) * Z2R (y + 1))%R. apply Rle_lt_trans with (Z2R (x + y + x * y)). rewrite <- (Rabs_pos_eq _ (Rlt_le _ _ Hxy)). apply Hexy. `````` Guillaume Melquiond committed Oct 01, 2010 284 ``````rewrite <- Z2R_mult. `````` Guillaume Melquiond committed May 17, 2010 285 286 287 ``````apply Z2R_lt. apply Zplus_lt_reg_r with (- (x + y + x * y + 1))%Z. now ring_simplify. `````` Guillaume Melquiond committed Oct 01, 2010 288 ``````rewrite bpow_plus. `````` Guillaume Melquiond committed May 17, 2010 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 ``````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. `````` Guillaume Melquiond committed May 17, 2010 308 309 310 311 312 313 314 315 316 ``````(* . *) 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. `````` Guillaume Melquiond committed May 17, 2010 317 318 ``````Qed. `````` Guillaume Melquiond committed Oct 01, 2010 319 ``End Fcalc_digits.``