src/HOL/Multivariate_Analysis/ex/Approximations.thy
 author wenzelm Wed Jan 13 23:07:06 2016 +0100 (2016-01-13) changeset 62175 8ffc4d0e652d parent 62089 4d38c04957fc child 62390 842917225d56 permissions -rw-r--r--
isabelle update_cartouches -c -t;
 eberlm@62085 1 section \Numeric approximations to Constants\ lp15@59871 2 lp15@59871 3 theory Approximations eberlm@62085 4 imports "../Complex_Transcendental" "../Harmonic_Numbers" eberlm@62085 5 begin eberlm@62085 6 eberlm@62089 7 text \ eberlm@62089 8 In this theory, we will approximate some standard mathematical constants with high precision, eberlm@62089 9 using only Isabelle's simplifier. (no oracles, code generator, etc.) eberlm@62089 10 eberlm@62089 11 The constants we will look at are: $\pi$, $e$, $\ln 2$, and $\gamma$ (the Euler--Mascheroni eberlm@62089 12 constant). eberlm@62089 13 \ eberlm@62089 14 eberlm@62085 15 lemma eval_fact: eberlm@62085 16 "fact 0 = 1" eberlm@62085 17 "fact (Suc 0) = 1" eberlm@62085 18 "fact (numeral n) = numeral n * fact (pred_numeral n)" eberlm@62085 19 by (simp, simp, simp_all only: numeral_eq_Suc fact_Suc, eberlm@62085 20 simp only: numeral_eq_Suc [symmetric] of_nat_numeral) eberlm@62085 21 eberlm@62085 22 lemma setsum_poly_horner_expand: eberlm@62085 23 "(\k<(numeral n::nat). f k * x^k) = f 0 + (\kkk<(0::nat). f k * x^k) = 0" eberlm@62085 26 proof - eberlm@62085 27 { eberlm@62085 28 fix m :: nat eberlm@62085 29 have "(\kk=Suc 0..k=Suc 0..kkkkk 0" "x \ 0" "x < 1" eberlm@62085 43 shows "x ^ n < (1::'a::linordered_semidom)" eberlm@62085 44 proof - eberlm@62085 45 from assms consider "x > 0" | "x = 0" by force eberlm@62085 46 thus ?thesis eberlm@62085 47 proof cases eberlm@62085 48 assume "x > 0" eberlm@62085 49 with assms show ?thesis eberlm@62085 50 by (cases n) (simp, hypsubst, rule power_Suc_less_one) eberlm@62085 51 qed (insert assms, cases n, simp_all) eberlm@62085 52 qed eberlm@62085 53 eberlm@62085 54 lemma combine_bounds: eberlm@62085 55 "x \ {a1..b1} \ y \ {a2..b2} \ a3 = a1 + a2 \ b3 = b1 + b2 \ x + y \ {a3..(b3::real)}" eberlm@62085 56 "x \ {a1..b1} \ y \ {a2..b2} \ a3 = a1 - b2 \ b3 = b1 - a2 \ x - y \ {a3..(b3::real)}" eberlm@62085 57 "c \ 0 \ x \ {a..b} \ c * x \ {c*a..c*b}" eberlm@62085 58 by (auto simp: mult_left_mono) eberlm@62085 59 eberlm@62085 60 lemma approx_coarsen: eberlm@62085 61 "\x - a1\ \ eps1 \ \a1 - a2\ \ eps2 - eps1 \ \x - a2\ \ (eps2 :: real)" eberlm@62085 62 by simp eberlm@62085 63 eberlm@62085 64 eberlm@62089 65 eberlm@62089 66 subsection \Approximations of the exponential function\ eberlm@62089 67 eberlm@62089 68 lemma two_power_fact_le_fact: eberlm@62089 69 assumes "n \ 1" eberlm@62089 70 shows "2^k * fact n \ (fact (n + k) :: 'a :: {semiring_char_0,linordered_semidom})" eberlm@62089 71 proof (induction k) eberlm@62089 72 case (Suc k) eberlm@62089 73 have "2 ^ Suc k * fact n = 2 * (2 ^ k * fact n)" by (simp add: algebra_simps) eberlm@62089 74 also note Suc.IH eberlm@62089 75 also from assms have "of_nat 1 + of_nat 1 \ of_nat n + (of_nat (Suc k) :: 'a)" eberlm@62089 76 by (intro add_mono) (unfold of_nat_le_iff, simp_all) eberlm@62089 77 hence "2 * (fact (n + k) :: 'a) \ of_nat (n + Suc k) * fact (n + k)" eberlm@62089 78 by (intro mult_right_mono) (simp_all add: add_ac) eberlm@62089 79 also have "\ = fact (n + Suc k)" by simp eberlm@62089 80 finally show ?case by - (simp add: mult_left_mono) eberlm@62089 81 qed simp_all eberlm@62085 82 eberlm@62089 83 text \ eberlm@62089 84 We approximate the exponential function with inputs between $0$ and $2$ by its eberlm@62089 85 Taylor series expansion and bound the error term with $0$ from below and with a eberlm@62089 86 geometric series from above. eberlm@62089 87 \ eberlm@62089 88 lemma exp_approx: eberlm@62089 89 assumes "n > 0" "0 \ x" "x < 2" eberlm@62089 90 shows "exp (x::real) - (\k {0..(2 * x^n / (2 - x)) / fact n}" eberlm@62089 91 proof (unfold atLeastAtMost_iff, safe) eberlm@62089 92 def approx \ "(\kk. x^k / fact k) sums exp x" eberlm@62089 94 using exp_converges[of x] by (simp add: field_simps) eberlm@62089 95 from sums_split_initial_segment[OF this, of n] eberlm@62089 96 have sums: "(\k. x^n * (x^k / fact (n+k))) sums (exp x - approx)" eberlm@62089 97 by (simp add: approx_def algebra_simps power_add) eberlm@62089 98 eberlm@62089 99 from assms show "(exp x - approx) \ 0" eberlm@62089 100 by (intro sums_le[OF _ sums_zero sums]) auto lp15@59871 101 eberlm@62089 102 have "\k. x^n * (x^k / fact (n+k)) \ (x^n / fact n) * (x / 2)^k" eberlm@62089 103 proof eberlm@62089 104 fix k :: nat eberlm@62089 105 have "x^n * (x^k / fact (n + k)) = x^(n+k) / fact (n + k)" by (simp add: power_add) eberlm@62089 106 also from assms have "\ \ x^(n+k) / (2^k * fact n)" eberlm@62089 107 by (intro divide_left_mono two_power_fact_le_fact zero_le_power) simp_all eberlm@62089 108 also have "\ = (x^n / fact n) * (x / 2) ^ k" eberlm@62089 109 by (simp add: field_simps power_add) eberlm@62089 110 finally show "x^n * (x^k / fact (n+k)) \ (x^n / fact n) * (x / 2)^k" . eberlm@62089 111 qed eberlm@62089 112 moreover note sums eberlm@62089 113 moreover { eberlm@62089 114 from assms have "(\k. (x^n / fact n) * (x / 2)^k) sums ((x^n / fact n) * (1 / (1 - x / 2)))" eberlm@62089 115 by (intro sums_mult geometric_sums) simp_all eberlm@62089 116 also from assms have "((x^n / fact n) * (1 / (1 - x / 2))) = (2 * x^n / (2 - x)) / fact n" eberlm@62089 117 by (auto simp: divide_simps) eberlm@62089 118 finally have "(\k. (x^n / fact n) * (x / 2)^k) sums \" . eberlm@62089 119 } eberlm@62089 120 ultimately show "(exp x - approx) \ (2 * x^n / (2 - x)) / fact n" eberlm@62089 121 by (rule sums_le) eberlm@62089 122 qed eberlm@62089 123 eberlm@62089 124 text \ eberlm@62089 125 The following variant gives a simpler error estimate for inputs between $0$ and $1$: eberlm@62089 126 \ eberlm@62089 127 lemma exp_approx': eberlm@62089 128 assumes "n > 0" "0 \ x" "x \ 1" eberlm@62089 129 shows "\exp (x::real) - (\k\n. x^k / fact k)\ \ x ^ n / fact n" eberlm@62085 130 proof - eberlm@62089 131 from assms have "x^n / (2 - x) \ x^n / 1" by (intro frac_le) simp_all eberlm@62089 132 hence "(2 * x^n / (2 - x)) / fact n \ 2 * x^n / fact n" eberlm@62089 133 using assms by (simp add: divide_simps) eberlm@62089 134 with exp_approx[of n x] assms eberlm@62089 135 have "exp (x::real) - (\k {0..2 * x^n / fact n}" by simp eberlm@62089 136 moreover have "(\k\n. x^k / fact k) = (\kexp (x::real) - (\k\n. x^k / fact k)\ \ x ^ n / fact n" eberlm@62089 139 unfolding atLeastAtMost_iff by linarith eberlm@62089 140 qed eberlm@62089 141 eberlm@62089 142 text \ eberlm@62089 143 By adding $x^n / n!$ to the approximation (i.e. taking one more term from the eberlm@62089 144 Taylor series), one can get the error bound down to $x^n / n!$. eberlm@62089 145 eberlm@62089 146 This means that the number of accurate binary digits produced by the approximation is eberlm@62089 147 asymptotically equal to $(n \log n - n) / \log 2$ by Stirling's formula. eberlm@62089 148 \ eberlm@62089 149 lemma exp_approx'': eberlm@62089 150 assumes "n > 0" "0 \ x" "x \ 1" eberlm@62089 151 shows "\exp (x::real) - (\k\n. x^k / fact k)\ \ 1 / fact n" eberlm@62089 152 proof - eberlm@62089 153 from assms have "\exp x - (\k\n. x ^ k / fact k)\ \ x ^ n / fact n" eberlm@62089 154 by (rule exp_approx') eberlm@62089 155 also from assms have "\ \ 1 / fact n" by (simp add: divide_simps power_le_one) eberlm@62085 156 finally show ?thesis . eberlm@62085 157 qed lp15@59871 158 eberlm@62089 159 eberlm@62089 160 text \ eberlm@62089 161 We now define an approximation function for Euler's constant $e$. eberlm@62089 162 \ eberlm@62089 163 eberlm@62089 164 definition euler_approx :: "nat \ real" where eberlm@62089 165 "euler_approx n = (\k\n. inverse (fact k))" eberlm@62089 166 eberlm@62089 167 definition euler_approx_aux :: "nat \ nat" where eberlm@62089 168 "euler_approx_aux n = (\k\n. \{k + 1..n})" eberlm@62089 169 eberlm@62089 170 lemma exp_1_approx: eberlm@62089 171 "n > 0 \ \exp (1::real) - euler_approx n\ \ 1 / fact n" eberlm@62089 172 using exp_approx''[of n 1] by (simp add: euler_approx_def divide_simps) eberlm@62089 173 eberlm@62089 174 text \ eberlm@62089 175 The following allows us to compute the numerator and the denominator of the result eberlm@62089 176 separately, which greatly reduces the amount of rational number arithmetic that we eberlm@62089 177 have to do. eberlm@62089 178 \ eberlm@62089 179 lemma euler_approx_altdef [code]: eberlm@62089 180 "euler_approx n = real (euler_approx_aux n) / real (fact n)" eberlm@62085 181 proof - eberlm@62089 182 have "real (\k\n. \{k+1..n}) = (\k\n. \i=k+1..n. real i)" by simp eberlm@62089 183 also have "\ / fact n = (\k\n. 1 / (fact n / (\i=k+1..n. real i)))" eberlm@62089 184 by (simp add: setsum_divide_distrib) eberlm@62089 185 also have "\ = (\k\n. 1 / fact k)" eberlm@62089 186 proof (intro setsum.cong refl) eberlm@62089 187 fix k assume k: "k \ {..n}" eberlm@62089 188 have "fact n = (\i=1..n. real i)" by (rule fact_altdef) eberlm@62089 189 also from k have "{1..n} = {1..k} \ {k+1..n}" by auto eberlm@62089 190 also have "setprod real \ / (\i=k+1..n. real i) = (\i=1..k. real i)" eberlm@62089 191 by (subst nonzero_divide_eq_eq, simp, subst setprod.union_disjoint [symmetric]) auto eberlm@62089 192 also have "\ = fact k" by (simp only: fact_altdef) eberlm@62089 193 finally show "1 / (fact n / setprod real {k + 1..n}) = 1 / fact k" by simp eberlm@62089 194 qed eberlm@62089 195 also have "\ = euler_approx n" by (simp add: euler_approx_def field_simps) eberlm@62089 196 finally show ?thesis by (simp add: euler_approx_aux_def) eberlm@62085 197 qed eberlm@62085 198 eberlm@62089 199 lemma euler_approx_aux_Suc: eberlm@62089 200 "euler_approx_aux (Suc m) = 1 + Suc m * euler_approx_aux m" eberlm@62089 201 unfolding euler_approx_aux_def eberlm@62089 202 by (subst setsum_right_distrib) (simp add: atLeastAtMostSuc_conv) eberlm@62089 203 eberlm@62089 204 lemma eval_euler_approx_aux: eberlm@62089 205 "euler_approx_aux 0 = 1" eberlm@62089 206 "euler_approx_aux 1 = 2" eberlm@62089 207 "euler_approx_aux (Suc 0) = 2" eberlm@62089 208 "euler_approx_aux (numeral n) = 1 + numeral n * euler_approx_aux (pred_numeral n)" (is "?th") eberlm@62089 209 proof - eberlm@62089 210 have A: "euler_approx_aux (Suc m) = 1 + Suc m * euler_approx_aux m" for m :: nat eberlm@62089 211 unfolding euler_approx_aux_def eberlm@62089 212 by (subst setsum_right_distrib) (simp add: atLeastAtMostSuc_conv) eberlm@62089 213 show ?th by (subst numeral_eq_Suc, subst A, subst numeral_eq_Suc [symmetric]) simp eberlm@62089 214 qed (simp_all add: euler_approx_aux_def) eberlm@62089 215 eberlm@62089 216 lemma euler_approx_aux_code [code]: eberlm@62089 217 "euler_approx_aux n = (if n = 0 then 1 else 1 + n * euler_approx_aux (n - 1))" eberlm@62089 218 by (cases n) (simp_all add: eval_euler_approx_aux euler_approx_aux_Suc) eberlm@62089 219 eberlm@62089 220 lemmas eval_euler_approx = euler_approx_altdef eval_euler_approx_aux eberlm@62089 221 eberlm@62089 222 eberlm@62089 223 text \Approximations of $e$ to 60 decimals / 128 and 64 bits:\ eberlm@62089 224 eberlm@62089 225 lemma euler_60_decimals: eberlm@62089 226 "\exp 1 - 2.718281828459045235360287471352662497757247093699959574966968\ eberlm@62089 227 \ inverse (10^60::real)" eberlm@62089 228 by (rule approx_coarsen, rule exp_1_approx[of 48]) eberlm@62089 229 (simp_all add: eval_euler_approx eval_fact) eberlm@62089 230 eberlm@62089 231 lemma euler_128: eberlm@62089 232 "\exp 1 - 924983374546220337150911035843336795079 / 2 ^ 128\ \ inverse (2 ^ 128 :: real)" eberlm@62089 233 by (rule approx_coarsen[OF euler_60_decimals]) simp_all eberlm@62089 234 eberlm@62089 235 lemma euler_64: eberlm@62089 236 "\exp 1 - 50143449209799256683 / 2 ^ 64\ \ inverse (2 ^ 64 :: real)" eberlm@62089 237 by (rule approx_coarsen[OF euler_128]) simp_all eberlm@62089 238 eberlm@62089 239 text \ eberlm@62089 240 An approximation of $e$ to 60 decimals. This is about as far as we can go with the eberlm@62089 241 simplifier with this kind of setup; the exported code of the code generator, on the other eberlm@62089 242 hand, can easily approximate $e$ to 1000 decimals and verify that approximation within eberlm@62089 243 fractions of a second. eberlm@62089 244 \ eberlm@62089 245 eberlm@62089 246 (* (Uncommented because we don't want to use the code generator; eberlm@62089 247 don't forget to import Code\_Target\_Numeral)) *) eberlm@62089 248 (* eberlm@62089 249 lemma "\exp 1 - 2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952605956307381323286279434907632338298807531952510190115738341879307021540891499348841675092447614606680822648001684774118537423454424371075390777449920695517027618386062613313845830007520449338265602976067371132007093287091274437470472306969772093101416928368190255151086574637721112523897844250569536967707854499699679468644549059879316368892300987931277361782154249992295763514822082698951936680331825288693984964651058209392398294887933203625094431173012381970684161403970198376793206832823764648042953118023287825098194558153017567173613320698112509961818815930416903515988885193458072738667385894228792284998920868058257492796104841984443634632449684875602336248270419786232090021609902353043699418491463140934317381436405462531520961836908887070167683964243781405927145635490613031072085103837505101157477041718986106873969655212671546889570350354021\ eberlm@62089 250 \ inverse (10^1000::real)" eberlm@62089 251 by (rule approx_coarsen, rule exp_1_approx[of 450], simp) eval eberlm@62089 252 *) eberlm@62089 253 eberlm@62089 254 eberlm@62089 255 subsection \Approximation of $\ln 2$\ eberlm@62089 256 eberlm@62089 257 text \ eberlm@62089 258 The following three auxiliary constants allow us to force the simplifier to eberlm@62089 259 evaluate intermediate results, simulating call-by-value. eberlm@62089 260 \ eberlm@62089 261 eberlm@62089 262 definition "ln_approx_aux3 x' e n y d \ eberlm@62089 263 \(2 * y) * (\k \ e - d" eberlm@62089 264 definition "ln_approx_aux2 x' e n y \ eberlm@62089 265 ln_approx_aux3 x' e n y (y^(2*n+1) / (1 - y^2) / real (2*n+1))" eberlm@62089 266 definition "ln_approx_aux1 x' e n x \ eberlm@62089 267 ln_approx_aux2 x' e n ((x - 1) / (x + 1))" eberlm@62089 268 eberlm@62089 269 lemma ln_approx_abs'': eberlm@62089 270 fixes x :: real and n :: nat eberlm@62089 271 defines "y \ (x-1)/(x+1)" eberlm@62089 272 defines "approx \ (\k y^(2*n+1) / (1 - y^2) / of_nat (2*n+1)" eberlm@62089 274 assumes x: "x > 1" eberlm@62089 275 assumes A: "ln_approx_aux1 x' e n x" eberlm@62089 276 shows "\ln x - x'\ \ e" eberlm@62089 277 proof (rule approx_coarsen[OF ln_approx_abs[OF x, of n]], goal_cases) eberlm@62089 278 case 1 eberlm@62089 279 from A have "\2 * y * (\k2 ^ k) + d - x'\ \ e - d" eberlm@62089 280 by (simp only: ln_approx_aux3_def ln_approx_aux2_def ln_approx_aux1_def eberlm@62089 281 y_def [symmetric] d_def [symmetric]) eberlm@62089 282 also have "2 * y * (\k2 ^ k) = eberlm@62089 283 (\k eberlm@62089 290 We unfold the above three constants successively and then compute the eberlm@62089 291 sum using a Horner scheme. eberlm@62089 292 \ eberlm@62089 293 lemma ln_2_40_decimals: eberlm@62089 294 "\ln 2 - 0.6931471805599453094172321214581765680755\ eberlm@62089 295 \ inverse (10^40 :: real)" eberlm@62089 296 apply (rule ln_approx_abs''[where n = 40], simp) eberlm@62089 297 apply (simp, simp add: ln_approx_aux1_def) eberlm@62089 298 apply (simp add: ln_approx_aux2_def power2_eq_square power_divide) eberlm@62089 299 apply (simp add: ln_approx_aux3_def power2_eq_square) eberlm@62089 300 apply (simp add: setsum_poly_horner_expand) eberlm@62089 301 done eberlm@62089 302 eberlm@62089 303 lemma ln_2_128: eberlm@62089 304 "\ln 2 - 235865763225513294137944142764154484399 / 2 ^ 128\ \ inverse (2 ^ 128 :: real)" eberlm@62089 305 by (rule approx_coarsen[OF ln_2_40_decimals]) simp_all eberlm@62089 306 eberlm@62089 307 lemma ln_2_64: eberlm@62089 308 "\ln 2 - 12786308645202655660 / 2 ^ 64\ \ inverse (2 ^ 64 :: real)" eberlm@62089 309 by (rule approx_coarsen[OF ln_2_128]) simp_all eberlm@62089 310 eberlm@62085 311 eberlm@62085 312 eberlm@62085 313 subsection \Approximation of the Euler--Mascheroni constant\ eberlm@62085 314 eberlm@62089 315 text \ eberlm@62089 316 Unfortunatly, the best approximation we have formalised for the Euler--Mascheroni eberlm@62089 317 constant converges only quadratically. This is too slow to compute more than a eberlm@62089 318 few decimals, but we can get almost 4 decimals / 14 binary digits this way, eberlm@62089 319 which is not too bad. eberlm@62089 320 \ eberlm@62089 321 lemma euler_mascheroni_approx: eberlm@62085 322 defines "approx \ 0.577257 :: real" and "e \ 0.000063 :: real" eberlm@62085 323 shows "abs (euler_mascheroni - approx :: real) < e" eberlm@62085 324 (is "abs (_ - ?approx) < ?e") lp15@59871 325 proof - eberlm@62085 326 def l \ "47388813395531028639296492901910937/82101866951584879688289000000000000 :: real" eberlm@62085 327 def u \ "142196984054132045946501548559032969 / 246305600854754639064867000000000000 :: real" eberlm@62085 328 have impI: "P \ Q" if Q for P Q using that by blast eberlm@62085 329 have hsum_63: "harm 63 = (310559566510213034489743057 / 65681493561267903750631200 ::real)" eberlm@62085 330 by (simp add: harm_expand) eberlm@62089 331 from harm_Suc[of 63] have hsum_64: "harm 64 = eberlm@62089 332 623171679694215690971693339 / (131362987122535807501262400::real)" eberlm@62085 333 by (subst (asm) hsum_63) simp eberlm@62085 334 have "ln (64::real) = real (6::nat) * ln 2" by (subst ln_realpow[symmetric]) simp_all eberlm@62089 335 hence "ln (real_of_nat (Suc 63)) \ {4.158883083293<..<4.158883083367}" using ln_2_64 eberlm@62089 336 by (simp add: abs_real_def split: split_if_asm) eberlm@62085 337 from euler_mascheroni_bounds'[OF _ this] eberlm@62089 338 have "(euler_mascheroni :: real) \ {l<.. \ {approx - e<..Approximation of pi\ eberlm@62085 349 eberlm@62085 350 eberlm@62085 351 subsubsection \Approximating the arctangent\ eberlm@62085 352 eberlm@62089 353 text\ eberlm@62089 354 The arctangent can be used to approximate pi. Fortunately, its Taylor series expansion eberlm@62089 355 converges exponentially for small values, so we can get $\Theta(n)$ digits of precision eberlm@62089 356 with $n$ summands of the expansion. eberlm@62089 357 \ eberlm@62089 358 eberlm@62085 359 definition arctan_approx where eberlm@62085 360 "arctan_approx n x = x * (\kx\ \ 1" eberlm@62085 364 shows "(\k. (-1)^k * (1 / real (k*2+1) * x ^ (k*2+1))) sums arctan x" eberlm@62085 365 using summable_arctan_series[OF assms] arctan_series[OF assms] by (simp add: sums_iff) eberlm@62085 366 eberlm@62085 367 lemma arctan_approx: eberlm@62085 368 assumes x: "0 \ x" "x < 1" and n: "even n" eberlm@62085 369 shows "arctan x - arctan_approx n x \ {0..x^(2*n+1) / (1-x^4)}" eberlm@62085 370 proof - eberlm@62085 371 def c \ "\k. 1 / (1+(4*real k + 2*real n)) - x\<^sup>2 / (3+(4*real k + 2*real n))" eberlm@62089 372 from assms have "(\k. (-1) ^ k * (1 / real (k * 2 + 1) * x^(k*2+1))) sums arctan x" eberlm@62085 373 using arctan_series' by simp eberlm@62089 374 also have "(\k. (-1) ^ k * (1 / real (k * 2 + 1) * x^(k*2+1))) = eberlm@62085 375 (\k. x * ((- (x^2))^k / real (2*k+1)))" eberlm@62085 376 by (simp add: power2_eq_square power_mult power_mult_distrib mult_ac power_minus') eberlm@62085 377 finally have "(\k. x * ((- x\<^sup>2) ^ k / real (2 * k + 1))) sums arctan x" . eberlm@62085 378 from sums_split_initial_segment[OF this, of n] eberlm@62085 379 have "(\i. x * ((- x\<^sup>2) ^ (i + n) / real (2 * (i + n) + 1))) sums eberlm@62085 380 (arctan x - arctan_approx n x)" eberlm@62085 381 by (simp add: arctan_approx_def setsum_right_distrib) eberlm@62085 382 from sums_group[OF this, of 2] assms eberlm@62085 383 have sums: "(\k. x * (x\<^sup>2)^n * (x^4)^k * c k) sums (arctan x - arctan_approx n x)" eberlm@62085 384 by (simp add: algebra_simps power_add power_mult [symmetric] c_def) eberlm@62089 385 eberlm@62085 386 from assms have "0 \ arctan x - arctan_approx n x" eberlm@62085 387 by (intro sums_le[OF _ sums_zero sums] allI mult_nonneg_nonneg) eberlm@62085 388 (auto intro!: frac_le power_le_one simp: c_def) eberlm@62085 389 moreover { eberlm@62085 390 from assms have "c k \ 1 - 0" for k unfolding c_def eberlm@62085 391 by (intro diff_mono divide_nonneg_nonneg add_nonneg_nonneg) auto eberlm@62085 392 with assms have "x * x\<^sup>2 ^ n * (x ^ 4) ^ k * c k \ x * x\<^sup>2 ^ n * (x ^ 4) ^ k * 1" for k eberlm@62085 393 by (intro mult_left_mono mult_right_mono mult_nonneg_nonneg) simp_all eberlm@62085 394 with assms have "arctan x - arctan_approx n x \ x * (x\<^sup>2)^n * (1 / (1 - x^4))" eberlm@62085 395 by (intro sums_le[OF _ sums sums_mult[OF geometric_sums]] allI mult_left_mono) eberlm@62085 396 (auto simp: power_less_one) eberlm@62085 397 also have "x * (x^2)^n = x^(2*n+1)" by (simp add: power_mult power_add) eberlm@62085 398 finally have "arctan x - arctan_approx n x \ x^(2*n+1) / (1 - x^4)" by simp eberlm@62085 399 } eberlm@62085 400 ultimately show ?thesis by simp eberlm@62085 401 qed eberlm@62085 402 eberlm@62089 403 lemma arctan_approx_def': "arctan_approx n (1/x) = eberlm@62089 404 (\k2) ^ k)) / x" eberlm@62085 405 proof - eberlm@62085 406 have "(-1)^k / b = 1 / ((-1)^k * b)" for k :: nat and b :: real eberlm@62085 407 by (cases "even k") auto eberlm@62085 408 thus ?thesis by (simp add: arctan_approx_def field_simps power_minus') eberlm@62085 409 qed eberlm@62085 410 eberlm@62085 411 lemma expand_arctan_approx: eberlm@62089 412 "(\k<(numeral n::nat). inverse (f k) * inverse (x ^ k)) = eberlm@62085 413 inverse (f 0) + (\kkk<(0::nat). inverse (f k) * inverse (x^k)) = 0" eberlm@62085 416 proof - eberlm@62085 417 { eberlm@62085 418 fix m :: nat eberlm@62085 419 have "(\kk=Suc 0..k=Suc 0..kkkkkx*y::real\ < 1" eberlm@62085 437 shows "arctan x - arctan y = arctan ((x - y) / (1 + x * y))" eberlm@62085 438 proof - eberlm@62085 439 have "arctan x - arctan y = arctan x + arctan (-y)" by (simp add: arctan_minus) eberlm@62085 440 also from assms have "\ = arctan ((x - y) / (1 + x * y))" by (subst arctan_add_small) simp_all eberlm@62085 441 finally show ?thesis . lp15@59871 442 qed lp15@59871 443 eberlm@62085 444 eberlm@62085 445 subsubsection \Machin-like formulae for pi\ eberlm@62085 446 eberlm@62085 447 text \ eberlm@62085 448 We first define a small proof method that can prove Machin-like formulae for @{term "pi"} eberlm@62089 449 automatically. Unfortunately, this takes far too much time for larger formulae because eberlm@62085 450 the numbers involved become too large. eberlm@62085 451 \ eberlm@62085 452 eberlm@62085 453 definition "MACHIN_TAG a b \ a * (b::real)" eberlm@62085 454 eberlm@62085 455 lemma numeral_horner_MACHIN_TAG: eberlm@62085 456 "MACHIN_TAG Numeral1 x = x" eberlm@62089 457 "MACHIN_TAG (numeral (Num.Bit0 (Num.Bit0 n))) x = eberlm@62085 458 MACHIN_TAG 2 (MACHIN_TAG (numeral (Num.Bit0 n)) x)" eberlm@62089 459 "MACHIN_TAG (numeral (Num.Bit0 (Num.Bit1 n))) x = eberlm@62085 460 MACHIN_TAG 2 (MACHIN_TAG (numeral (Num.Bit1 n)) x)" eberlm@62089 461 "MACHIN_TAG (numeral (Num.Bit1 n)) x = eberlm@62085 462 MACHIN_TAG 2 (MACHIN_TAG (numeral n) x) + x" eberlm@62085 463 unfolding numeral_Bit0 numeral_Bit1 ring_distribs one_add_one[symmetric] MACHIN_TAG_def eberlm@62085 464 by (simp_all add: algebra_simps) eberlm@62085 465 eberlm@62085 466 lemma tag_machin: "a * arctan b = MACHIN_TAG a (arctan b)" by (simp add: MACHIN_TAG_def) eberlm@62085 467 eberlm@62085 468 lemma arctan_double': "\a::real\ < 1 \ MACHIN_TAG 2 (arctan a) = arctan (2 * a / (1 - a*a))" eberlm@62085 469 unfolding MACHIN_TAG_def by (simp add: arctan_double power2_eq_square) eberlm@62085 470 eberlm@62085 471 ML \ eberlm@62085 472 fun machin_term_conv ctxt ct = eberlm@62085 473 let eberlm@62085 474 val ctxt' = ctxt addsimps @{thms arctan_double' arctan_add_small} eberlm@62085 475 in eberlm@62085 476 case Thm.term_of ct of eberlm@62089 477 Const (@{const_name MACHIN_TAG}, _) $_$ eberlm@62089 478 (Const (@{const_name "Transcendental.arctan"}, _) $_) => eberlm@62085 479 Simplifier.rewrite ctxt' ct eberlm@62085 480 | eberlm@62089 481 Const (@{const_name MACHIN_TAG}, _)$ _ $eberlm@62089 482 (Const (@{const_name "Groups.plus"}, _)$ eberlm@62085 483 (Const (@{const_name "Transcendental.arctan"}, _) $_)$ eberlm@62089 484 (Const (@{const_name "Transcendental.arctan"}, _) $_)) => eberlm@62085 485 Simplifier.rewrite ctxt' ct eberlm@62085 486 | _ => raise CTERM ("machin_conv", [ct]) eberlm@62085 487 end eberlm@62085 488 eberlm@62089 489 fun machin_tac ctxt = eberlm@62085 490 let val conv = Conv.top_conv (Conv.try_conv o machin_term_conv) ctxt eberlm@62085 491 in eberlm@62085 492 SELECT_GOAL ( eberlm@62089 493 Local_Defs.unfold_tac ctxt eberlm@62085 494 @{thms tag_machin[THEN eq_reflection] numeral_horner_MACHIN_TAG[THEN eq_reflection]} eberlm@62085 495 THEN REPEAT (CHANGED (HEADGOAL (CONVERSION conv)))) eberlm@62085 496 THEN' Simplifier.simp_tac (ctxt addsimps @{thms arctan_add_small arctan_diff_small}) eberlm@62085 497 end eberlm@62085 498 \ eberlm@62085 499 eberlm@62085 500 method_setup machin = \Scan.succeed (SIMPLE_METHOD' o machin_tac)\ eberlm@62085 501 eberlm@62085 502 text \ eberlm@62089 503 We can now prove the standard'' Machin formula, which was already proven manually eberlm@62085 504 in Isabelle, automatically. eberlm@62085 505 }\ eberlm@62085 506 lemma "pi / 4 = (4::real) * arctan (1 / 5) - arctan (1 / 239)" eberlm@62085 507 by machin eberlm@62085 508 eberlm@62085 509 text \ eberlm@62085 510 We can also prove the following more complicated formula: eberlm@62085 511 \ eberlm@62085 512 lemma machin': "pi/4 = (12::real) * arctan (1/18) + 8 * arctan (1/57) - 5 * arctan (1/239)" eberlm@62085 513 by machin eberlm@62085 514 eberlm@62085 515 eberlm@62085 516 eberlm@62085 517 subsubsection \Simple approximation of pi\ eberlm@62085 518 eberlm@62085 519 text \ eberlm@62085 520 We can use the simple Machin formula and the Taylor series expansion of the arctangent eberlm@62089 521 to approximate pi. For a given even natural number$n$, we expand @{term "arctan (1/5)"} eberlm@62085 522 to$3n$summands and @{term "arctan (1/239)"} to$n$summands. This gives us at least eberlm@62085 523$13n-2$bits of precision. eberlm@62085 524 \ eberlm@62085 525 eberlm@62085 526 definition "pi_approx n = 16 * arctan_approx (3*n) (1/5) - 4 * arctan_approx n (1/239)" eberlm@62085 527 eberlm@62085 528 lemma pi_approx: eberlm@62085 529 fixes n :: nat assumes n: "even n" and "n > 0" eberlm@62085 530 shows "\pi - pi_approx n\ \ inverse (2^(13*n - 2))" eberlm@62085 531 proof - eberlm@62085 532 from n have n': "even (3*n)" by simp wenzelm@62175 533 \ \We apply the Machin formula\ eberlm@62085 534 from machin have "pi = 16 * arctan (1/5) - 4 * arctan (1/239::real)" by simp wenzelm@62175 535 \ \Taylor series expansion of the arctangent\ eberlm@62085 536 also from arctan_approx[OF _ _ n', of "1/5"] arctan_approx[OF _ _ n, of "1/239"] eberlm@62085 537 have "\ - pi_approx n \ {-4*((1/239)^(2*n+1) / (1-(1/239)^4))..16*(1/5)^(6*n+1) / (1-(1/5)^4)}" eberlm@62085 538 by (simp add: pi_approx_def) wenzelm@62175 539 \ \Coarsening the bounds to make them a bit nicer\ eberlm@62085 540 also have "-4*((1/239::real)^(2*n+1) / (1-(1/239)^4)) = -((13651919 / 815702160) / 57121^n)" eberlm@62085 541 by (simp add: power_mult power2_eq_square) (simp add: field_simps) eberlm@62085 542 also have "16*(1/5)^(6*n+1) / (1-(1/5::real)^4) = (125/39) / 15625^n" eberlm@62085 543 by (simp add: power_mult power2_eq_square) (simp add: field_simps) eberlm@62089 544 also have "{-((13651919 / 815702160) / 57121^n) .. (125 / 39) / 15625^n} \ eberlm@62085 545 {- (4 / 2^(13*n)) .. 4 / (2^(13*n)::real)}" eberlm@62085 546 by (subst atLeastatMost_subset_iff, intro disjI2 conjI le_imp_neg_le) eberlm@62085 547 (rule frac_le; simp add: power_mult power_mono)+ eberlm@62085 548 finally have "abs (pi - pi_approx n) \ 4 / 2^(13*n)" by auto eberlm@62085 549 also from \n > 0\ have "4 / 2^(13*n) = 1 / (2^(13*n - 2) :: real)" eberlm@62085 550 by (cases n) (simp_all add: power_add) eberlm@62085 551 finally show ?thesis by (simp add: divide_inverse) eberlm@62085 552 qed eberlm@62085 553 eberlm@62085 554 lemma pi_approx': eberlm@62085 555 fixes n :: nat assumes n: "even n" and "n > 0" and "k \ 13*n - 2" eberlm@62085 556 shows "\pi - pi_approx n\ \ inverse (2^k)" eberlm@62085 557 using assms(3) by (intro order.trans[OF pi_approx[OF assms(1,2)]]) (simp_all add: field_simps) eberlm@62085 558 eberlm@62085 559 text \We can now approximate pi to 22 decimals within a fraction of a second.\ eberlm@62085 560 lemma pi_approx_75: "abs (pi - 3.1415926535897932384626 :: real) \ inverse (10^22)" eberlm@62085 561 proof - eberlm@62085 562 def a \ "8295936325956147794769600190539918304 / 2626685325478320010006427764892578125 :: real" eberlm@62085 563 def b \ "8428294561696506782041394632 / 503593538783547230635598424135 :: real" wenzelm@62175 564 \ \The introduction of this constant prevents the simplifier from applying solvers that eberlm@62085 565 we don't want. We want it to simply evaluate the terms to rational constants.}\ eberlm@62085 566 def eq \ "op = :: real \ real \ bool" eberlm@62089 567 wenzelm@62175 568 \ \Splitting the computation into several steps has the advantage that simplification can eberlm@62085 569 be done in parallel\ eberlm@62085 570 have "abs (pi - pi_approx 6) \ inverse (2^76)" by (rule pi_approx') simp_all eberlm@62089 571 also have "pi_approx 6 = 16 * arctan_approx (3 * 6) (1 / 5) - 4 * arctan_approx 6 (1 / 239)" eberlm@62085 572 unfolding pi_approx_def by simp eberlm@62085 573 also have [unfolded eq_def]: "eq (16 * arctan_approx (3 * 6) (1 / 5)) a" eberlm@62085 574 by (simp add: arctan_approx_def' power2_eq_square, eberlm@62085 575 simp add: expand_arctan_approx, unfold a_def eq_def, rule refl) eberlm@62085 576 also have [unfolded eq_def]: "eq (4 * arctan_approx 6 (1 / 239::real)) b" eberlm@62085 577 by (simp add: arctan_approx_def' power2_eq_square, eberlm@62085 578 simp add: expand_arctan_approx, unfold b_def eq_def, rule refl) eberlm@62089 579 also have [unfolded eq_def]: eberlm@62085 580 "eq (a - b) (171331331860120333586637094112743033554946184594977368554649608 / eberlm@62085 581 54536456744112171868276045488779391002026386559009552001953125)" eberlm@62085 582 by (unfold a_def b_def, simp, unfold eq_def, rule refl) eberlm@62085 583 finally show ?thesis by (rule approx_coarsen) simp eberlm@62085 584 qed eberlm@62085 585 eberlm@62085 586 text \ eberlm@62089 587 The previous estimate of pi in this file was based on approximating the root of the eberlm@62089 588$\sin(\pi/6)$in the interval$[0;4]$using the Taylor series expansion of the sine to eberlm@62085 589 verify that it is between two given bounds. eberlm@62089 590 This was much slower and much less precise. We can easily recover this coarser estimate from eberlm@62085 591 the newer, precise estimate: eberlm@62085 592 \ eberlm@62085 593 lemma pi_approx_32: "\pi - 13493037705/4294967296 :: real\ \ inverse(2 ^ 32)" eberlm@62085 594 by (rule approx_coarsen[OF pi_approx_75]) simp eberlm@62085 595 eberlm@62085 596 eberlm@62085 597 subsection \A more complicated approximation of pi\ eberlm@62085 598 eberlm@62085 599 text \ eberlm@62089 600 There are more complicated Machin-like formulae that have more terms with larger eberlm@62085 601 denominators. Although they have more terms, each term requires fewer summands of the eberlm@62085 602 Taylor series for the same precision, since it is evaluated closer to$0\$. eberlm@62089 603 eberlm@62085 604 Using a good formula, one can therefore obtain the same precision with fewer operations. eberlm@62089 605 The big formulae used for computations of pi in practice are too complicated for us to eberlm@62085 606 prove here, but we can use the three-term Machin-like formula @{thm machin'}. eberlm@62085 607 \ eberlm@62085 608 eberlm@62089 609 definition "pi_approx2 n = 48 * arctan_approx (6*n) (1/18::real) + eberlm@62085 610 32 * arctan_approx (4*n) (1/57) - 20 * arctan_approx (3*n) (1/239)" eberlm@62085 611 eberlm@62085 612 lemma pi_approx2: eberlm@62085 613 fixes n :: nat assumes n: "even n" and "n > 0" eberlm@62085 614 shows "\pi - pi_approx2 n\ \ inverse (2^(46*n - 1))" eberlm@62085 615 proof - eberlm@62085 616 from n have n': "even (6*n)" "even (4*n)" "even (3*n)" by simp_all eberlm@62089 617 from machin' have "pi = 48 * arctan (1/18) + 32 * arctan (1/57) - 20 * arctan (1/239::real)" eberlm@62085 618 by simp eberlm@62085 619 hence "pi - pi_approx2 n = 48 * (arctan (1/18) - arctan_approx (6*n) (1/18)) + eberlm@62085 620 32 * (arctan (1/57) - arctan_approx (4*n) (1/57)) - eberlm@62085 621 20 * (arctan (1/239) - arctan_approx (3*n) (1/239))" eberlm@62085 622 by (simp add: pi_approx2_def) eberlm@62085 623 also have "\ \ {-((20/239/(1-(1/239)^4)) * (1/239)^(6*n)).. eberlm@62089 624 (48/18 / (1-(1/18)^4))*(1/18)^(12*n) + (32/57/(1-(1/57)^4)) * (1/57)^(8*n)}" eberlm@62085 625 (is "_ \ {-?l..?u1 + ?u2}") eberlm@62085 626 apply ((rule combine_bounds(1,2))+; (rule combine_bounds(3); (rule arctan_approx)?)?) eberlm@62085 627 apply (simp_all add: n) eberlm@62085 628 apply (simp_all add: divide_simps)? eberlm@62085 629 done eberlm@62085 630 also { eberlm@62085 631 have "?l \ (1/8) * (1/2)^(46*n)" eberlm@62085 632 unfolding power_mult by (intro mult_mono power_mono) (simp_all add: divide_simps) eberlm@62085 633 also have "\ \ (1/2) ^ (46 * n - 1)" eberlm@62085 634 by (cases n; simp_all add: power_add divide_simps) eberlm@62085 635 finally have "?l \ (1/2) ^ (46 * n - 1)" . eberlm@62085 636 moreover { eberlm@62085 637 have "?u1 + ?u2 \ 4 * (1/2)^(48*n) + 1 * (1/2)^(46*n)" eberlm@62085 638 unfolding power_mult by (intro add_mono mult_mono power_mono) (simp_all add: divide_simps) eberlm@62089 639 also from \n > 0\ have "4 * (1/2::real)^(48*n) \ (1/2)^(46*n)" eberlm@62085 640 by (cases n) (simp_all add: field_simps power_add) eberlm@62085 641 also from \n > 0\ have "(1/2::real) ^ (46 * n) + 1 * (1 / 2) ^ (46 * n) = (1/2) ^ (46 * n - 1)" eberlm@62085 642 by (cases n; simp_all add: power_add power_divide) eberlm@62085 643 finally have "?u1 + ?u2 \ (1/2) ^ (46 * n - 1)" by - simp eberlm@62085 644 } eberlm@62085 645 ultimately have "{-?l..?u1 + ?u2} \ {-((1/2)^(46*n-1))..(1/2)^(46*n-1)}" eberlm@62085 646 by (subst atLeastatMost_subset_iff) simp_all eberlm@62085 647 } eberlm@62085 648 finally have "\pi - pi_approx2 n\ \ ((1/2) ^ (46 * n - 1))" by auto eberlm@62085 649 thus ?thesis by (simp add: divide_simps) eberlm@62085 650 qed eberlm@62085 651 eberlm@62085 652 lemma pi_approx2': eberlm@62085 653 fixes n :: nat assumes n: "even n" and "n > 0" and "k \ 46*n - 1" eberlm@62085 654 shows "\pi - pi_approx2 n\ \ inverse (2^k)" eberlm@62085 655 using assms(3) by (intro order.trans[OF pi_approx2[OF assms(1,2)]]) (simp_all add: field_simps) eberlm@62085 656 eberlm@62085 657 text \ eberlm@62089 658 We can now approximate pi to 54 decimals using this formula. The computations are much eberlm@62089 659 slower now; this is mostly because we use arbitrary-precision rational numbers, whose eberlm@62089 660 numerators and demoninators get very large. Using dyadic floating point numbers would be eberlm@62085 661 much more economical. eberlm@62085 662 \ eberlm@62089 663 lemma pi_approx_54_decimals: eberlm@62085 664 "abs (pi - 3.141592653589793238462643383279502884197169399375105821 :: real) \ inverse (10^54)" eberlm@62085 665 (is "abs (pi - ?pi') \ _") eberlm@62085 666 proof - eberlm@62085 667 def a \ "2829469759662002867886529831139137601191652261996513014734415222704732791803 / eberlm@62085 668 1062141879292765061960538947347721564047051545995266466660439319087625011200 :: real" eberlm@62085 669 def b \ "13355545553549848714922837267299490903143206628621657811747118592 / eberlm@62085 670 23792006023392488526789546722992491355941103837356113731091180925 :: real" eberlm@62085 671 def c \ "28274063397213534906669125255762067746830085389618481175335056 / eberlm@62085 672 337877029279505250241149903214554249587517250716358486542628059 :: real" eberlm@62085 673 let ?pi'' = "3882327391761098513316067116522233897127356523627918964967729040413954225768920394233198626889767468122598417405434625348404038165437924058179155035564590497837027530349 / eberlm@62089 674 1235783190199688165469648572769847552336447197542738425378629633275352407743112409829873464564018488572820294102599160968781449606552922108667790799771278860366957772800" eberlm@62085 675 def eq \ "op = :: real \ real \ bool" eberlm@62089 676 eberlm@62085 677 have "abs (pi - pi_approx2 4) \ inverse (2^183)" by (rule pi_approx2') simp_all eberlm@62085 678 also have "pi_approx2 4 = 48 * arctan_approx 24 (1 / 18) + eberlm@62085 679 32 * arctan_approx 16 (1 / 57) - eberlm@62089 680 20 * arctan_approx 12 (1 / 239)" eberlm@62085 681 unfolding pi_approx2_def by simp eberlm@62085 682 also have [unfolded eq_def]: "eq (48 * arctan_approx 24 (1 / 18)) a" eberlm@62085 683 by (simp add: arctan_approx_def' power2_eq_square, eberlm@62085 684 simp add: expand_arctan_approx, unfold a_def eq_def, rule refl) eberlm@62085 685 also have [unfolded eq_def]: "eq (32 * arctan_approx 16 (1 / 57::real)) b" eberlm@62085 686 by (simp add: arctan_approx_def' power2_eq_square, eberlm@62085 687 simp add: expand_arctan_approx, unfold b_def eq_def, rule refl) eberlm@62085 688 also have [unfolded eq_def]: "eq (20 * arctan_approx 12 (1 / 239::real)) c" eberlm@62085 689 by (simp add: arctan_approx_def' power2_eq_square, eberlm@62085 690 simp add: expand_arctan_approx, unfold c_def eq_def, rule refl) eberlm@62085 691 also have [unfolded eq_def]: eberlm@62085 692 "eq (a + b) (34326487387865555303797183505809267914709125998469664969258315922216638779011304447624792548723974104030355722677 / eberlm@62089 693 10642967245546718617684989689985787964158885991018703366677373121531695267093031090059801733340658960857196134400)" eberlm@62085 694 by (unfold a_def b_def c_def, simp, unfold eq_def, rule refl) eberlm@62085 695 also have [unfolded eq_def]: "eq (\ - c) ?pi''" eberlm@62085 696 by (unfold a_def b_def c_def, simp, unfold eq_def, rule refl) wenzelm@62175 697 \ \This is incredibly slow because the numerators and denominators are huge.\ eberlm@62085 698 finally show ?thesis by (rule approx_coarsen) simp eberlm@62085 699 qed eberlm@62085 700 eberlm@62085 701 text \A 128 bit approximation of pi:\ eberlm@62085 702 lemma pi_approx_128: eberlm@62085 703 "abs (pi - 1069028584064966747859680373161870783301 / 2^128) \ inverse (2^128)" eberlm@62085 704 by (rule approx_coarsen[OF pi_approx_54_decimals]) simp eberlm@62085 705 eberlm@62085 706 text \A 64 bit approximation of pi:\ eberlm@62089 707 lemma pi_approx_64: eberlm@62085 708 "abs (pi - 57952155664616982739 / 2^64 :: real) \ inverse (2^64)" eberlm@62085 709 by (rule approx_coarsen[OF pi_approx_54_decimals]) simp eberlm@62089 710 eberlm@62089 711 text \ eberlm@62089 712 Again, going much farther with the simplifier takes a long time, but the code generator eberlm@62089 713 can handle even two thousand decimal digits in under 20 seconds. eberlm@62089 714 \ lp15@59871 715 lp15@59871 716 end