src/HOL/Analysis/ex/Approximations.thy
 author fleury Mon Sep 19 20:06:21 2016 +0200 (2016-09-19) changeset 63918 6bf55e6e0b75 parent 63627 6ddb43c6b711 child 64267 b9a1486e79be permissions -rw-r--r--
left_distrib ~> distrib_right, right_distrib ~> distrib_left
 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)  wenzelm@63040  92  define approx where "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}"  haftmann@63417  188  have "fact n = (\i=1..n. real i)" by (simp add: fact_setprod)  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  haftmann@63417  192  also have "\ = fact k" by (simp add: fact_setprod)  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  Mathias@63918  202  by (subst setsum_distrib_left) (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  Mathias@63918  212  by (subst setsum_distrib_left) (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 -  wenzelm@63040  326  define l :: real  wenzelm@63040  327  where "l = 47388813395531028639296492901910937/82101866951584879688289000000000000"  wenzelm@63040  328  define u :: real  wenzelm@63040  329  where "u = 142196984054132045946501548559032969 / 246305600854754639064867000000000000"  eberlm@62085  330  have impI: "P \ Q" if Q for P Q using that by blast  wenzelm@63040  331  have hsum_63: "harm 63 = (310559566510213034489743057 / 65681493561267903750631200 :: real)"  eberlm@62085  332  by (simp add: harm_expand)  eberlm@62089  333  from harm_Suc[of 63] have hsum_64: "harm 64 =  eberlm@62089  334  623171679694215690971693339 / (131362987122535807501262400::real)"  eberlm@62085  335  by (subst (asm) hsum_63) simp  eberlm@62085  336  have "ln (64::real) = real (6::nat) * ln 2" by (subst ln_realpow[symmetric]) simp_all  eberlm@62089  337  hence "ln (real_of_nat (Suc 63)) \ {4.158883083293<..<4.158883083367}" using ln_2_64  nipkow@62390  338  by (simp add: abs_real_def split: if_split_asm)  eberlm@62085  339  from euler_mascheroni_bounds'[OF _ this]  eberlm@62089  340  have "(euler_mascheroni :: real) \ {l<.. \ {approx - e<..Approximation of pi\  eberlm@62085  351 eberlm@62085  352 eberlm@62085  353 subsubsection \Approximating the arctangent\  eberlm@62085  354 eberlm@62089  355 text\  eberlm@62089  356  The arctangent can be used to approximate pi. Fortunately, its Taylor series expansion  eberlm@62089  357  converges exponentially for small values, so we can get $\Theta(n)$ digits of precision  eberlm@62089  358  with $n$ summands of the expansion.  eberlm@62089  359 \  eberlm@62089  360 eberlm@62085  361 definition arctan_approx where  eberlm@62085  362  "arctan_approx n x = x * (\kx\ \ 1"  eberlm@62085  366  shows "(\k. (-1)^k * (1 / real (k*2+1) * x ^ (k*2+1))) sums arctan x"  eberlm@62085  367  using summable_arctan_series[OF assms] arctan_series[OF assms] by (simp add: sums_iff)  eberlm@62085  368 eberlm@62085  369 lemma arctan_approx:  eberlm@62085  370  assumes x: "0 \ x" "x < 1" and n: "even n"  eberlm@62085  371  shows "arctan x - arctan_approx n x \ {0..x^(2*n+1) / (1-x^4)}"  eberlm@62085  372 proof -  wenzelm@63040  373  define c where "c k = 1 / (1+(4*real k + 2*real n)) - x\<^sup>2 / (3+(4*real k + 2*real n))" for k  eberlm@62089  374  from assms have "(\k. (-1) ^ k * (1 / real (k * 2 + 1) * x^(k*2+1))) sums arctan x"  eberlm@62085  375  using arctan_series' by simp  eberlm@62089  376  also have "(\k. (-1) ^ k * (1 / real (k * 2 + 1) * x^(k*2+1))) =  eberlm@62085  377  (\k. x * ((- (x^2))^k / real (2*k+1)))"  eberlm@62085  378  by (simp add: power2_eq_square power_mult power_mult_distrib mult_ac power_minus')  eberlm@62085  379  finally have "(\k. x * ((- x\<^sup>2) ^ k / real (2 * k + 1))) sums arctan x" .  eberlm@62085  380  from sums_split_initial_segment[OF this, of n]  eberlm@62085  381  have "(\i. x * ((- x\<^sup>2) ^ (i + n) / real (2 * (i + n) + 1))) sums  eberlm@62085  382  (arctan x - arctan_approx n x)"  Mathias@63918  383  by (simp add: arctan_approx_def setsum_distrib_left)  eberlm@62085  384  from sums_group[OF this, of 2] assms  eberlm@62085  385  have sums: "(\k. x * (x\<^sup>2)^n * (x^4)^k * c k) sums (arctan x - arctan_approx n x)"  eberlm@62085  386  by (simp add: algebra_simps power_add power_mult [symmetric] c_def)  eberlm@62089  387 eberlm@62085  388  from assms have "0 \ arctan x - arctan_approx n x"  eberlm@62085  389  by (intro sums_le[OF _ sums_zero sums] allI mult_nonneg_nonneg)  eberlm@62085  390  (auto intro!: frac_le power_le_one simp: c_def)  eberlm@62085  391  moreover {  eberlm@62085  392  from assms have "c k \ 1 - 0" for k unfolding c_def  eberlm@62085  393  by (intro diff_mono divide_nonneg_nonneg add_nonneg_nonneg) auto  eberlm@62085  394  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  395  by (intro mult_left_mono mult_right_mono mult_nonneg_nonneg) simp_all  eberlm@62085  396  with assms have "arctan x - arctan_approx n x \ x * (x\<^sup>2)^n * (1 / (1 - x^4))"  eberlm@62085  397  by (intro sums_le[OF _ sums sums_mult[OF geometric_sums]] allI mult_left_mono)  eberlm@62085  398  (auto simp: power_less_one)  eberlm@62085  399  also have "x * (x^2)^n = x^(2*n+1)" by (simp add: power_mult power_add)  eberlm@62085  400  finally have "arctan x - arctan_approx n x \ x^(2*n+1) / (1 - x^4)" by simp  eberlm@62085  401  }  eberlm@62085  402  ultimately show ?thesis by simp  eberlm@62085  403 qed  eberlm@62085  404 eberlm@62089  405 lemma arctan_approx_def': "arctan_approx n (1/x) =  eberlm@62089  406  (\k2) ^ k)) / x"  eberlm@62085  407 proof -  eberlm@62085  408  have "(-1)^k / b = 1 / ((-1)^k * b)" for k :: nat and b :: real  eberlm@62085  409  by (cases "even k") auto  eberlm@62085  410  thus ?thesis by (simp add: arctan_approx_def field_simps power_minus')  eberlm@62085  411 qed  eberlm@62085  412 eberlm@62085  413 lemma expand_arctan_approx:  eberlm@62089  414  "(\k<(numeral n::nat). inverse (f k) * inverse (x ^ k)) =  eberlm@62085  415  inverse (f 0) + (\kkk<(0::nat). inverse (f k) * inverse (x^k)) = 0"  eberlm@62085  418 proof -  eberlm@62085  419  {  eberlm@62085  420  fix m :: nat  eberlm@62085  421  have "(\kk=Suc 0..k=Suc 0..kkkkkx*y::real\ < 1"  eberlm@62085  439  shows "arctan x - arctan y = arctan ((x - y) / (1 + x * y))"  eberlm@62085  440 proof -  eberlm@62085  441  have "arctan x - arctan y = arctan x + arctan (-y)" by (simp add: arctan_minus)  eberlm@62085  442  also from assms have "\ = arctan ((x - y) / (1 + x * y))" by (subst arctan_add_small) simp_all  eberlm@62085  443  finally show ?thesis .  lp15@59871  444 qed  lp15@59871  445 eberlm@62085  446 eberlm@62085  447 subsubsection \Machin-like formulae for pi\  eberlm@62085  448 eberlm@62085  449 text \  eberlm@62085  450  We first define a small proof method that can prove Machin-like formulae for @{term "pi"}  eberlm@62089  451  automatically. Unfortunately, this takes far too much time for larger formulae because  eberlm@62085  452  the numbers involved become too large.  eberlm@62085  453 \  eberlm@62085  454 eberlm@62085  455 definition "MACHIN_TAG a b \ a * (b::real)"  eberlm@62085  456 eberlm@62085  457 lemma numeral_horner_MACHIN_TAG:  eberlm@62085  458  "MACHIN_TAG Numeral1 x = x"  eberlm@62089  459  "MACHIN_TAG (numeral (Num.Bit0 (Num.Bit0 n))) x =  eberlm@62085  460  MACHIN_TAG 2 (MACHIN_TAG (numeral (Num.Bit0 n)) x)"  eberlm@62089  461  "MACHIN_TAG (numeral (Num.Bit0 (Num.Bit1 n))) x =  eberlm@62085  462  MACHIN_TAG 2 (MACHIN_TAG (numeral (Num.Bit1 n)) x)"  eberlm@62089  463  "MACHIN_TAG (numeral (Num.Bit1 n)) x =  eberlm@62085  464  MACHIN_TAG 2 (MACHIN_TAG (numeral n) x) + x"  eberlm@62085  465  unfolding numeral_Bit0 numeral_Bit1 ring_distribs one_add_one[symmetric] MACHIN_TAG_def  eberlm@62085  466  by (simp_all add: algebra_simps)  eberlm@62085  467 eberlm@62085  468 lemma tag_machin: "a * arctan b = MACHIN_TAG a (arctan b)" by (simp add: MACHIN_TAG_def)  eberlm@62085  469 eberlm@62085  470 lemma arctan_double': "\a::real\ < 1 \ MACHIN_TAG 2 (arctan a) = arctan (2 * a / (1 - a*a))"  eberlm@62085  471  unfolding MACHIN_TAG_def by (simp add: arctan_double power2_eq_square)  eberlm@62085  472 eberlm@62085  473 ML \  eberlm@62085  474  fun machin_term_conv ctxt ct =  eberlm@62085  475  let  eberlm@62085  476  val ctxt' = ctxt addsimps @{thms arctan_double' arctan_add_small}  eberlm@62085  477  in  eberlm@62085  478  case Thm.term_of ct of  eberlm@62089  479  Const (@{const_name MACHIN_TAG}, _) $_$  eberlm@62089  480  (Const (@{const_name "Transcendental.arctan"}, _) $_) =>  eberlm@62085  481  Simplifier.rewrite ctxt' ct  eberlm@62085  482  |  eberlm@62089  483  Const (@{const_name MACHIN_TAG}, _)$ _ $ eberlm@62089  484  (Const (@{const_name "Groups.plus"}, _)$  eberlm@62085  485  (Const (@{const_name "Transcendental.arctan"}, _) $_)$  eberlm@62089  486  (Const (@{const_name "Transcendental.arctan"}, _) $_)) =>  eberlm@62085  487  Simplifier.rewrite ctxt' ct  eberlm@62085  488  | _ => raise CTERM ("machin_conv", [ct])  eberlm@62085  489  end  eberlm@62085  490 eberlm@62089  491  fun machin_tac ctxt =  eberlm@62085  492  let val conv = Conv.top_conv (Conv.try_conv o machin_term_conv) ctxt  eberlm@62085  493  in  eberlm@62085  494  SELECT_GOAL (  eberlm@62089  495  Local_Defs.unfold_tac ctxt  eberlm@62085  496  @{thms tag_machin[THEN eq_reflection] numeral_horner_MACHIN_TAG[THEN eq_reflection]}  eberlm@62085  497  THEN REPEAT (CHANGED (HEADGOAL (CONVERSION conv))))  eberlm@62085  498  THEN' Simplifier.simp_tac (ctxt addsimps @{thms arctan_add_small arctan_diff_small})  eberlm@62085  499  end  eberlm@62085  500 \  eberlm@62085  501 eberlm@62085  502 method_setup machin = \Scan.succeed (SIMPLE_METHOD' o machin_tac)\  eberlm@62085  503 eberlm@62085  504 text \  eberlm@62089  505  We can now prove the standard'' Machin formula, which was already proven manually  eberlm@62085  506  in Isabelle, automatically.  eberlm@62085  507 }\  eberlm@62085  508 lemma "pi / 4 = (4::real) * arctan (1 / 5) - arctan (1 / 239)"  eberlm@62085  509  by machin  eberlm@62085  510 eberlm@62085  511 text \  eberlm@62085  512  We can also prove the following more complicated formula:  eberlm@62085  513 \  eberlm@62085  514 lemma machin': "pi/4 = (12::real) * arctan (1/18) + 8 * arctan (1/57) - 5 * arctan (1/239)"  eberlm@62085  515  by machin  eberlm@62085  516 eberlm@62085  517 eberlm@62085  518 eberlm@62085  519 subsubsection \Simple approximation of pi\  eberlm@62085  520 eberlm@62085  521 text \  eberlm@62085  522  We can use the simple Machin formula and the Taylor series expansion of the arctangent  eberlm@62089  523  to approximate pi. For a given even natural number$n$, we expand @{term "arctan (1/5)"}  eberlm@62085  524  to$3n$summands and @{term "arctan (1/239)"} to$n$summands. This gives us at least  eberlm@62085  525 $13n-2$bits of precision.  eberlm@62085  526 \  eberlm@62085  527 eberlm@62085  528 definition "pi_approx n = 16 * arctan_approx (3*n) (1/5) - 4 * arctan_approx n (1/239)"  eberlm@62085  529 eberlm@62085  530 lemma pi_approx:  eberlm@62085  531  fixes n :: nat assumes n: "even n" and "n > 0"  eberlm@62085  532  shows "\pi - pi_approx n\ \ inverse (2^(13*n - 2))"  eberlm@62085  533 proof -  eberlm@62085  534  from n have n': "even (3*n)" by simp  wenzelm@62175  535  \ \We apply the Machin formula\  eberlm@62085  536  from machin have "pi = 16 * arctan (1/5) - 4 * arctan (1/239::real)" by simp  wenzelm@62175  537  \ \Taylor series expansion of the arctangent\  eberlm@62085  538  also from arctan_approx[OF _ _ n', of "1/5"] arctan_approx[OF _ _ n, of "1/239"]  eberlm@62085  539  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  540  by (simp add: pi_approx_def)  wenzelm@62175  541  \ \Coarsening the bounds to make them a bit nicer\  eberlm@62085  542  also have "-4*((1/239::real)^(2*n+1) / (1-(1/239)^4)) = -((13651919 / 815702160) / 57121^n)"  eberlm@62085  543  by (simp add: power_mult power2_eq_square) (simp add: field_simps)  eberlm@62085  544  also have "16*(1/5)^(6*n+1) / (1-(1/5::real)^4) = (125/39) / 15625^n"  eberlm@62085  545  by (simp add: power_mult power2_eq_square) (simp add: field_simps)  eberlm@62089  546  also have "{-((13651919 / 815702160) / 57121^n) .. (125 / 39) / 15625^n} \  eberlm@62085  547  {- (4 / 2^(13*n)) .. 4 / (2^(13*n)::real)}"  eberlm@62085  548  by (subst atLeastatMost_subset_iff, intro disjI2 conjI le_imp_neg_le)  eberlm@62085  549  (rule frac_le; simp add: power_mult power_mono)+  eberlm@62085  550  finally have "abs (pi - pi_approx n) \ 4 / 2^(13*n)" by auto  eberlm@62085  551  also from \n > 0\ have "4 / 2^(13*n) = 1 / (2^(13*n - 2) :: real)"  eberlm@62085  552  by (cases n) (simp_all add: power_add)  eberlm@62085  553  finally show ?thesis by (simp add: divide_inverse)  eberlm@62085  554 qed  eberlm@62085  555 eberlm@62085  556 lemma pi_approx':  eberlm@62085  557  fixes n :: nat assumes n: "even n" and "n > 0" and "k \ 13*n - 2"  eberlm@62085  558  shows "\pi - pi_approx n\ \ inverse (2^k)"  eberlm@62085  559  using assms(3) by (intro order.trans[OF pi_approx[OF assms(1,2)]]) (simp_all add: field_simps)  eberlm@62085  560 eberlm@62085  561 text \We can now approximate pi to 22 decimals within a fraction of a second.\  eberlm@62085  562 lemma pi_approx_75: "abs (pi - 3.1415926535897932384626 :: real) \ inverse (10^22)"  eberlm@62085  563 proof -  wenzelm@63040  564  define a :: real  wenzelm@63040  565  where "a = 8295936325956147794769600190539918304 / 2626685325478320010006427764892578125"  wenzelm@63040  566  define b :: real  wenzelm@63040  567  where "b = 8428294561696506782041394632 / 503593538783547230635598424135"  wenzelm@62175  568  \ \The introduction of this constant prevents the simplifier from applying solvers that  eberlm@62085  569  we don't want. We want it to simply evaluate the terms to rational constants.}\  wenzelm@63040  570  define eq :: "real \ real \ bool" where "eq = op ="  eberlm@62089  571 wenzelm@62175  572  \ \Splitting the computation into several steps has the advantage that simplification can  eberlm@62085  573  be done in parallel\  eberlm@62085  574  have "abs (pi - pi_approx 6) \ inverse (2^76)" by (rule pi_approx') simp_all  eberlm@62089  575  also have "pi_approx 6 = 16 * arctan_approx (3 * 6) (1 / 5) - 4 * arctan_approx 6 (1 / 239)"  eberlm@62085  576  unfolding pi_approx_def by simp  eberlm@62085  577  also have [unfolded eq_def]: "eq (16 * arctan_approx (3 * 6) (1 / 5)) a"  eberlm@62085  578  by (simp add: arctan_approx_def' power2_eq_square,  eberlm@62085  579  simp add: expand_arctan_approx, unfold a_def eq_def, rule refl)  eberlm@62085  580  also have [unfolded eq_def]: "eq (4 * arctan_approx 6 (1 / 239::real)) b"  eberlm@62085  581  by (simp add: arctan_approx_def' power2_eq_square,  eberlm@62085  582  simp add: expand_arctan_approx, unfold b_def eq_def, rule refl)  eberlm@62089  583  also have [unfolded eq_def]:  eberlm@62085  584  "eq (a - b) (171331331860120333586637094112743033554946184594977368554649608 /  eberlm@62085  585  54536456744112171868276045488779391002026386559009552001953125)"  eberlm@62085  586  by (unfold a_def b_def, simp, unfold eq_def, rule refl)  eberlm@62085  587  finally show ?thesis by (rule approx_coarsen) simp  eberlm@62085  588 qed  eberlm@62085  589 eberlm@62085  590 text \  eberlm@62089  591  The previous estimate of pi in this file was based on approximating the root of the  eberlm@62089  592 $\sin(\pi/6)$in the interval$[0;4]$using the Taylor series expansion of the sine to  eberlm@62085  593  verify that it is between two given bounds.  eberlm@62089  594  This was much slower and much less precise. We can easily recover this coarser estimate from  eberlm@62085  595  the newer, precise estimate:  eberlm@62085  596 \  eberlm@62085  597 lemma pi_approx_32: "\pi - 13493037705/4294967296 :: real\ \ inverse(2 ^ 32)"  eberlm@62085  598  by (rule approx_coarsen[OF pi_approx_75]) simp  eberlm@62085  599 eberlm@62085  600 eberlm@62085  601 subsection \A more complicated approximation of pi\  eberlm@62085  602 eberlm@62085  603 text \  eberlm@62089  604  There are more complicated Machin-like formulae that have more terms with larger  eberlm@62085  605  denominators. Although they have more terms, each term requires fewer summands of the  eberlm@62085  606  Taylor series for the same precision, since it is evaluated closer to$0\$.  eberlm@62089  607 eberlm@62085  608  Using a good formula, one can therefore obtain the same precision with fewer operations.  eberlm@62089  609  The big formulae used for computations of pi in practice are too complicated for us to  eberlm@62085  610  prove here, but we can use the three-term Machin-like formula @{thm machin'}.  eberlm@62085  611 \  eberlm@62085  612 eberlm@62089  613 definition "pi_approx2 n = 48 * arctan_approx (6*n) (1/18::real) +  eberlm@62085  614  32 * arctan_approx (4*n) (1/57) - 20 * arctan_approx (3*n) (1/239)"  eberlm@62085  615 eberlm@62085  616 lemma pi_approx2:  eberlm@62085  617  fixes n :: nat assumes n: "even n" and "n > 0"  eberlm@62085  618  shows "\pi - pi_approx2 n\ \ inverse (2^(46*n - 1))"  eberlm@62085  619 proof -  eberlm@62085  620  from n have n': "even (6*n)" "even (4*n)" "even (3*n)" by simp_all  eberlm@62089  621  from machin' have "pi = 48 * arctan (1/18) + 32 * arctan (1/57) - 20 * arctan (1/239::real)"  eberlm@62085  622  by simp  eberlm@62085  623  hence "pi - pi_approx2 n = 48 * (arctan (1/18) - arctan_approx (6*n) (1/18)) +  eberlm@62085  624  32 * (arctan (1/57) - arctan_approx (4*n) (1/57)) -  eberlm@62085  625  20 * (arctan (1/239) - arctan_approx (3*n) (1/239))"  eberlm@62085  626  by (simp add: pi_approx2_def)  eberlm@62085  627  also have "\ \ {-((20/239/(1-(1/239)^4)) * (1/239)^(6*n))..  eberlm@62089  628  (48/18 / (1-(1/18)^4))*(1/18)^(12*n) + (32/57/(1-(1/57)^4)) * (1/57)^(8*n)}"  eberlm@62085  629  (is "_ \ {-?l..?u1 + ?u2}")  eberlm@62085  630  apply ((rule combine_bounds(1,2))+; (rule combine_bounds(3); (rule arctan_approx)?)?)  eberlm@62085  631  apply (simp_all add: n)  eberlm@62085  632  apply (simp_all add: divide_simps)?  eberlm@62085  633  done  eberlm@62085  634  also {  eberlm@62085  635  have "?l \ (1/8) * (1/2)^(46*n)"  eberlm@62085  636  unfolding power_mult by (intro mult_mono power_mono) (simp_all add: divide_simps)  eberlm@62085  637  also have "\ \ (1/2) ^ (46 * n - 1)"  eberlm@62085  638  by (cases n; simp_all add: power_add divide_simps)  eberlm@62085  639  finally have "?l \ (1/2) ^ (46 * n - 1)" .  eberlm@62085  640  moreover {  eberlm@62085  641  have "?u1 + ?u2 \ 4 * (1/2)^(48*n) + 1 * (1/2)^(46*n)"  eberlm@62085  642  unfolding power_mult by (intro add_mono mult_mono power_mono) (simp_all add: divide_simps)  eberlm@62089  643  also from \n > 0\ have "4 * (1/2::real)^(48*n) \ (1/2)^(46*n)"  eberlm@62085  644  by (cases n) (simp_all add: field_simps power_add)  eberlm@62085  645  also from \n > 0\ have "(1/2::real) ^ (46 * n) + 1 * (1 / 2) ^ (46 * n) = (1/2) ^ (46 * n - 1)"  eberlm@62085  646  by (cases n; simp_all add: power_add power_divide)  eberlm@62085  647  finally have "?u1 + ?u2 \ (1/2) ^ (46 * n - 1)" by - simp  eberlm@62085  648  }  eberlm@62085  649  ultimately have "{-?l..?u1 + ?u2} \ {-((1/2)^(46*n-1))..(1/2)^(46*n-1)}"  eberlm@62085  650  by (subst atLeastatMost_subset_iff) simp_all  eberlm@62085  651  }  eberlm@62085  652  finally have "\pi - pi_approx2 n\ \ ((1/2) ^ (46 * n - 1))" by auto  eberlm@62085  653  thus ?thesis by (simp add: divide_simps)  eberlm@62085  654 qed  eberlm@62085  655 eberlm@62085  656 lemma pi_approx2':  eberlm@62085  657  fixes n :: nat assumes n: "even n" and "n > 0" and "k \ 46*n - 1"  eberlm@62085  658  shows "\pi - pi_approx2 n\ \ inverse (2^k)"  eberlm@62085  659  using assms(3) by (intro order.trans[OF pi_approx2[OF assms(1,2)]]) (simp_all add: field_simps)  eberlm@62085  660 eberlm@62085  661 text \  eberlm@62089  662  We can now approximate pi to 54 decimals using this formula. The computations are much  eberlm@62089  663  slower now; this is mostly because we use arbitrary-precision rational numbers, whose  eberlm@62089  664  numerators and demoninators get very large. Using dyadic floating point numbers would be  eberlm@62085  665  much more economical.  eberlm@62085  666 \  eberlm@62089  667 lemma pi_approx_54_decimals:  eberlm@62085  668  "abs (pi - 3.141592653589793238462643383279502884197169399375105821 :: real) \ inverse (10^54)"  eberlm@62085  669  (is "abs (pi - ?pi') \ _")  eberlm@62085  670 proof -  wenzelm@63040  671  define a :: real  wenzelm@63040  672  where "a = 2829469759662002867886529831139137601191652261996513014734415222704732791803 /  wenzelm@63040  673  1062141879292765061960538947347721564047051545995266466660439319087625011200"  wenzelm@63040  674  define b :: real  wenzelm@63040  675  where "b = 13355545553549848714922837267299490903143206628621657811747118592 /  wenzelm@63040  676  23792006023392488526789546722992491355941103837356113731091180925"  wenzelm@63040  677  define c :: real  wenzelm@63040  678  where "c = 28274063397213534906669125255762067746830085389618481175335056 /  wenzelm@63040  679  337877029279505250241149903214554249587517250716358486542628059"  eberlm@62085  680  let ?pi'' = "3882327391761098513316067116522233897127356523627918964967729040413954225768920394233198626889767468122598417405434625348404038165437924058179155035564590497837027530349 /  eberlm@62089  681  1235783190199688165469648572769847552336447197542738425378629633275352407743112409829873464564018488572820294102599160968781449606552922108667790799771278860366957772800"  wenzelm@63040  682  define eq :: "real \ real \ bool" where "eq = op ="  eberlm@62089  683 eberlm@62085  684  have "abs (pi - pi_approx2 4) \ inverse (2^183)" by (rule pi_approx2') simp_all  eberlm@62085  685  also have "pi_approx2 4 = 48 * arctan_approx 24 (1 / 18) +  eberlm@62085  686  32 * arctan_approx 16 (1 / 57) -  eberlm@62089  687  20 * arctan_approx 12 (1 / 239)"  eberlm@62085  688  unfolding pi_approx2_def by simp  eberlm@62085  689  also have [unfolded eq_def]: "eq (48 * arctan_approx 24 (1 / 18)) a"  eberlm@62085  690  by (simp add: arctan_approx_def' power2_eq_square,  eberlm@62085  691  simp add: expand_arctan_approx, unfold a_def eq_def, rule refl)  eberlm@62085  692  also have [unfolded eq_def]: "eq (32 * arctan_approx 16 (1 / 57::real)) b"  eberlm@62085  693  by (simp add: arctan_approx_def' power2_eq_square,  eberlm@62085  694  simp add: expand_arctan_approx, unfold b_def eq_def, rule refl)  eberlm@62085  695  also have [unfolded eq_def]: "eq (20 * arctan_approx 12 (1 / 239::real)) c"  eberlm@62085  696  by (simp add: arctan_approx_def' power2_eq_square,  eberlm@62085  697  simp add: expand_arctan_approx, unfold c_def eq_def, rule refl)  eberlm@62085  698  also have [unfolded eq_def]:  eberlm@62085  699  "eq (a + b) (34326487387865555303797183505809267914709125998469664969258315922216638779011304447624792548723974104030355722677 /  eberlm@62089  700  10642967245546718617684989689985787964158885991018703366677373121531695267093031090059801733340658960857196134400)"  eberlm@62085  701  by (unfold a_def b_def c_def, simp, unfold eq_def, rule refl)  eberlm@62085  702  also have [unfolded eq_def]: "eq (\ - c) ?pi''"  eberlm@62085  703  by (unfold a_def b_def c_def, simp, unfold eq_def, rule refl)  wenzelm@62175  704  \ \This is incredibly slow because the numerators and denominators are huge.\  eberlm@62085  705  finally show ?thesis by (rule approx_coarsen) simp  eberlm@62085  706 qed  eberlm@62085  707 eberlm@62085  708 text \A 128 bit approximation of pi:\  eberlm@62085  709 lemma pi_approx_128:  eberlm@62085  710  "abs (pi - 1069028584064966747859680373161870783301 / 2^128) \ inverse (2^128)"  eberlm@62085  711  by (rule approx_coarsen[OF pi_approx_54_decimals]) simp  eberlm@62085  712 eberlm@62085  713 text \A 64 bit approximation of pi:\  eberlm@62089  714 lemma pi_approx_64:  eberlm@62085  715  "abs (pi - 57952155664616982739 / 2^64 :: real) \ inverse (2^64)"  eberlm@62085  716  by (rule approx_coarsen[OF pi_approx_54_decimals]) simp  eberlm@62089  717   eberlm@62089  718 text \  eberlm@62089  719  Again, going much farther with the simplifier takes a long time, but the code generator  eberlm@62089  720  can handle even two thousand decimal digits in under 20 seconds.  eberlm@62089  721 \  lp15@59871  722 lp15@59871  723 end