src/HOL/Analysis/ex/Approximations.thy
 changeset 63627 6ddb43c6b711 parent 63417 c184ec919c70 child 63918 6bf55e6e0b75
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/HOL/Analysis/ex/Approximations.thy	Mon Aug 08 14:13:14 2016 +0200
1.3 @@ -0,0 +1,723 @@
1.4 +section \<open>Numeric approximations to Constants\<close>
1.5 +
1.6 +theory Approximations
1.7 +imports "../Complex_Transcendental" "../Harmonic_Numbers"
1.8 +begin
1.9 +
1.10 +text \<open>
1.11 +  In this theory, we will approximate some standard mathematical constants with high precision,
1.12 +  using only Isabelle's simplifier. (no oracles, code generator, etc.)
1.13 +
1.14 +  The constants we will look at are: $\pi$, $e$, $\ln 2$, and $\gamma$ (the Euler--Mascheroni
1.15 +  constant).
1.16 +\<close>
1.17 +
1.18 +lemma eval_fact:
1.19 +  "fact 0 = 1"
1.20 +  "fact (Suc 0) = 1"
1.21 +  "fact (numeral n) = numeral n * fact (pred_numeral n)"
1.22 +  by (simp, simp, simp_all only: numeral_eq_Suc fact_Suc,
1.23 +      simp only: numeral_eq_Suc [symmetric] of_nat_numeral)
1.24 +
1.25 +lemma setsum_poly_horner_expand:
1.26 +  "(\<Sum>k<(numeral n::nat). f k * x^k) = f 0 + (\<Sum>k<pred_numeral n. f (k+1) * x^k) * x"
1.27 +  "(\<Sum>k<Suc 0. f k * x^k) = (f 0 :: 'a :: semiring_1)"
1.28 +  "(\<Sum>k<(0::nat). f k * x^k) = 0"
1.29 +proof -
1.30 +  {
1.31 +    fix m :: nat
1.32 +    have "(\<Sum>k<Suc m. f k * x^k) = f 0 + (\<Sum>k=Suc 0..<Suc m. f k * x^k)"
1.33 +      by (subst atLeast0LessThan [symmetric], subst setsum_head_upt_Suc) simp_all
1.34 +    also have "(\<Sum>k=Suc 0..<Suc m. f k * x^k) = (\<Sum>k<m. f (k+1) * x^k) * x"
1.35 +      by (subst setsum_shift_bounds_Suc_ivl)
1.36 +         (simp add: setsum_left_distrib algebra_simps atLeast0LessThan power_commutes)
1.37 +    finally have "(\<Sum>k<Suc m. f k * x ^ k) = f 0 + (\<Sum>k<m. f (k + 1) * x ^ k) * x" .
1.38 +  }
1.39 +  from this[of "pred_numeral n"]
1.40 +    show "(\<Sum>k<numeral n. f k * x^k) = f 0 + (\<Sum>k<pred_numeral n. f (k+1) * x^k) * x"
1.41 +    by (simp add: numeral_eq_Suc)
1.42 +qed simp_all
1.43 +
1.44 +lemma power_less_one:
1.45 +  assumes "n > 0" "x \<ge> 0" "x < 1"
1.46 +  shows   "x ^ n < (1::'a::linordered_semidom)"
1.47 +proof -
1.48 +  from assms consider "x > 0" | "x = 0" by force
1.49 +  thus ?thesis
1.50 +  proof cases
1.51 +    assume "x > 0"
1.52 +    with assms show ?thesis
1.53 +      by (cases n) (simp, hypsubst, rule power_Suc_less_one)
1.54 +  qed (insert assms, cases n, simp_all)
1.55 +qed
1.56 +
1.57 +lemma combine_bounds:
1.58 +  "x \<in> {a1..b1} \<Longrightarrow> y \<in> {a2..b2} \<Longrightarrow> a3 = a1 + a2 \<Longrightarrow> b3 = b1 + b2 \<Longrightarrow> x + y \<in> {a3..(b3::real)}"
1.59 +  "x \<in> {a1..b1} \<Longrightarrow> y \<in> {a2..b2} \<Longrightarrow> a3 = a1 - b2 \<Longrightarrow> b3 = b1 - a2 \<Longrightarrow> x - y \<in> {a3..(b3::real)}"
1.60 +  "c \<ge> 0 \<Longrightarrow> x \<in> {a..b} \<Longrightarrow> c * x \<in> {c*a..c*b}"
1.61 +  by (auto simp: mult_left_mono)
1.62 +
1.63 +lemma approx_coarsen:
1.64 +  "\<bar>x - a1\<bar> \<le> eps1 \<Longrightarrow> \<bar>a1 - a2\<bar> \<le> eps2 - eps1 \<Longrightarrow> \<bar>x - a2\<bar> \<le> (eps2 :: real)"
1.65 +  by simp
1.66 +
1.67 +
1.68 +
1.69 +subsection \<open>Approximations of the exponential function\<close>
1.70 +
1.71 +lemma two_power_fact_le_fact:
1.72 +  assumes "n \<ge> 1"
1.73 +  shows   "2^k * fact n \<le> (fact (n + k) :: 'a :: {semiring_char_0,linordered_semidom})"
1.74 +proof (induction k)
1.75 +  case (Suc k)
1.76 +  have "2 ^ Suc k * fact n = 2 * (2 ^ k * fact n)" by (simp add: algebra_simps)
1.77 +  also note Suc.IH
1.78 +  also from assms have "of_nat 1 + of_nat 1 \<le> of_nat n + (of_nat (Suc k) :: 'a)"
1.79 +    by (intro add_mono) (unfold of_nat_le_iff, simp_all)
1.80 +  hence "2 * (fact (n + k) :: 'a) \<le> of_nat (n + Suc k) * fact (n + k)"
1.82 +  also have "\<dots> = fact (n + Suc k)" by simp
1.83 +  finally show ?case by - (simp add: mult_left_mono)
1.84 +qed simp_all
1.85 +
1.86 +text \<open>
1.87 +  We approximate the exponential function with inputs between $0$ and $2$ by its
1.88 +  Taylor series expansion and bound the error term with $0$ from below and with a
1.89 +  geometric series from above.
1.90 +\<close>
1.91 +lemma exp_approx:
1.92 +  assumes "n > 0" "0 \<le> x" "x < 2"
1.93 +  shows   "exp (x::real) - (\<Sum>k<n. x^k / fact k) \<in> {0..(2 * x^n / (2 - x)) / fact n}"
1.94 +proof (unfold atLeastAtMost_iff, safe)
1.95 +  define approx where "approx = (\<Sum>k<n. x^k / fact k)"
1.96 +  have "(\<lambda>k. x^k / fact k) sums exp x"
1.97 +    using exp_converges[of x] by (simp add: field_simps)
1.98 +  from sums_split_initial_segment[OF this, of n]
1.99 +    have sums: "(\<lambda>k. x^n * (x^k / fact (n+k))) sums (exp x - approx)"
1.101 +
1.102 +  from assms show "(exp x - approx) \<ge> 0"
1.103 +    by (intro sums_le[OF _ sums_zero sums]) auto
1.104 +
1.105 +  have "\<forall>k. x^n * (x^k / fact (n+k)) \<le> (x^n / fact n) * (x / 2)^k"
1.106 +  proof
1.107 +    fix k :: nat
1.108 +    have "x^n * (x^k / fact (n + k)) = x^(n+k) / fact (n + k)" by (simp add: power_add)
1.109 +    also from assms have "\<dots> \<le> x^(n+k) / (2^k * fact n)"
1.110 +      by (intro divide_left_mono two_power_fact_le_fact zero_le_power) simp_all
1.111 +    also have "\<dots> = (x^n / fact n) * (x / 2) ^ k"
1.113 +    finally show "x^n * (x^k / fact (n+k)) \<le> (x^n / fact n) * (x / 2)^k" .
1.114 +  qed
1.115 +  moreover note sums
1.116 +  moreover {
1.117 +    from assms have "(\<lambda>k. (x^n / fact n) * (x / 2)^k) sums ((x^n / fact n) * (1 / (1 - x / 2)))"
1.118 +      by (intro sums_mult geometric_sums) simp_all
1.119 +    also from assms have "((x^n / fact n) * (1 / (1 - x / 2))) = (2 * x^n / (2 - x)) / fact n"
1.120 +      by (auto simp: divide_simps)
1.121 +    finally have "(\<lambda>k. (x^n / fact n) * (x / 2)^k) sums \<dots>" .
1.122 +  }
1.123 +  ultimately show "(exp x - approx) \<le> (2 * x^n / (2 - x)) / fact n"
1.124 +    by (rule sums_le)
1.125 +qed
1.126 +
1.127 +text \<open>
1.128 +  The following variant gives a simpler error estimate for inputs between $0$ and $1$:
1.129 +\<close>
1.130 +lemma exp_approx':
1.131 +  assumes "n > 0" "0 \<le> x" "x \<le> 1"
1.132 +  shows   "\<bar>exp (x::real) - (\<Sum>k\<le>n. x^k / fact k)\<bar> \<le> x ^ n / fact n"
1.133 +proof -
1.134 +  from assms have "x^n / (2 - x) \<le> x^n / 1" by (intro frac_le) simp_all
1.135 +  hence "(2 * x^n / (2 - x)) / fact n \<le> 2 * x^n / fact n"
1.136 +    using assms by (simp add: divide_simps)
1.137 +  with exp_approx[of n x] assms
1.138 +    have "exp (x::real) - (\<Sum>k<n. x^k / fact k) \<in> {0..2 * x^n / fact n}" by simp
1.139 +  moreover have "(\<Sum>k\<le>n. x^k / fact k) = (\<Sum>k<n. x^k / fact k) + x ^ n / fact n"
1.140 +    by (simp add: lessThan_Suc_atMost [symmetric])
1.141 +  ultimately show "\<bar>exp (x::real) - (\<Sum>k\<le>n. x^k / fact k)\<bar> \<le> x ^ n / fact n"
1.142 +    unfolding atLeastAtMost_iff by linarith
1.143 +qed
1.144 +
1.145 +text \<open>
1.146 +  By adding $x^n / n!$ to the approximation (i.e. taking one more term from the
1.147 +  Taylor series), one can get the error bound down to $x^n / n!$.
1.148 +
1.149 +  This means that the number of accurate binary digits produced by the approximation is
1.150 +  asymptotically equal to $(n \log n - n) / \log 2$ by Stirling's formula.
1.151 +\<close>
1.152 +lemma exp_approx'':
1.153 +  assumes "n > 0" "0 \<le> x" "x \<le> 1"
1.154 +  shows   "\<bar>exp (x::real) - (\<Sum>k\<le>n. x^k / fact k)\<bar> \<le> 1 / fact n"
1.155 +proof -
1.156 +  from assms have "\<bar>exp x - (\<Sum>k\<le>n. x ^ k / fact k)\<bar> \<le> x ^ n / fact n"
1.157 +    by (rule exp_approx')
1.158 +  also from assms have "\<dots> \<le> 1 / fact n" by (simp add: divide_simps power_le_one)
1.159 +  finally show ?thesis .
1.160 +qed
1.161 +
1.162 +
1.163 +text \<open>
1.164 +  We now define an approximation function for Euler's constant $e$.
1.165 +\<close>
1.166 +
1.167 +definition euler_approx :: "nat \<Rightarrow> real" where
1.168 +  "euler_approx n = (\<Sum>k\<le>n. inverse (fact k))"
1.169 +
1.170 +definition euler_approx_aux :: "nat \<Rightarrow> nat" where
1.171 +  "euler_approx_aux n = (\<Sum>k\<le>n. \<Prod>{k + 1..n})"
1.172 +
1.173 +lemma exp_1_approx:
1.174 +  "n > 0 \<Longrightarrow> \<bar>exp (1::real) - euler_approx n\<bar> \<le> 1 / fact n"
1.175 +  using exp_approx''[of n 1] by (simp add: euler_approx_def divide_simps)
1.176 +
1.177 +text \<open>
1.178 +  The following allows us to compute the numerator and the denominator of the result
1.179 +  separately, which greatly reduces the amount of rational number arithmetic that we
1.180 +  have to do.
1.181 +\<close>
1.182 +lemma euler_approx_altdef [code]:
1.183 +  "euler_approx n = real (euler_approx_aux n) / real (fact n)"
1.184 +proof -
1.185 +  have "real (\<Sum>k\<le>n. \<Prod>{k+1..n}) = (\<Sum>k\<le>n. \<Prod>i=k+1..n. real i)" by simp
1.186 +  also have "\<dots> / fact n = (\<Sum>k\<le>n. 1 / (fact n / (\<Prod>i=k+1..n. real i)))"
1.187 +    by (simp add: setsum_divide_distrib)
1.188 +  also have "\<dots> = (\<Sum>k\<le>n. 1 / fact k)"
1.189 +  proof (intro setsum.cong refl)
1.190 +    fix k assume k: "k \<in> {..n}"
1.191 +    have "fact n = (\<Prod>i=1..n. real i)" by (simp add: fact_setprod)
1.192 +    also from k have "{1..n} = {1..k} \<union> {k+1..n}" by auto
1.193 +    also have "setprod real \<dots> / (\<Prod>i=k+1..n. real i) = (\<Prod>i=1..k. real i)"
1.194 +      by (subst nonzero_divide_eq_eq, simp, subst setprod.union_disjoint [symmetric]) auto
1.195 +    also have "\<dots> = fact k" by (simp add: fact_setprod)
1.196 +    finally show "1 / (fact n / setprod real {k + 1..n}) = 1 / fact k" by simp
1.197 +  qed
1.198 +  also have "\<dots> = euler_approx n" by (simp add: euler_approx_def field_simps)
1.199 +  finally show ?thesis by (simp add: euler_approx_aux_def)
1.200 +qed
1.201 +
1.202 +lemma euler_approx_aux_Suc:
1.203 +  "euler_approx_aux (Suc m) = 1 + Suc m * euler_approx_aux m"
1.204 +  unfolding euler_approx_aux_def
1.205 +  by (subst setsum_right_distrib) (simp add: atLeastAtMostSuc_conv)
1.206 +
1.207 +lemma eval_euler_approx_aux:
1.208 +  "euler_approx_aux 0 = 1"
1.209 +  "euler_approx_aux 1 = 2"
1.210 +  "euler_approx_aux (Suc 0) = 2"
1.211 +  "euler_approx_aux (numeral n) = 1 + numeral n * euler_approx_aux (pred_numeral n)" (is "?th")
1.212 +proof -
1.213 +  have A: "euler_approx_aux (Suc m) = 1 + Suc m * euler_approx_aux m" for m :: nat
1.214 +    unfolding euler_approx_aux_def
1.215 +    by (subst setsum_right_distrib) (simp add: atLeastAtMostSuc_conv)
1.216 +  show ?th by (subst numeral_eq_Suc, subst A, subst numeral_eq_Suc [symmetric]) simp
1.218 +
1.219 +lemma euler_approx_aux_code [code]:
1.220 +  "euler_approx_aux n = (if n = 0 then 1 else 1 + n * euler_approx_aux (n - 1))"
1.221 +  by (cases n) (simp_all add: eval_euler_approx_aux euler_approx_aux_Suc)
1.222 +
1.223 +lemmas eval_euler_approx = euler_approx_altdef eval_euler_approx_aux
1.224 +
1.225 +
1.226 +text \<open>Approximations of $e$ to 60 decimals / 128 and 64 bits:\<close>
1.227 +
1.228 +lemma euler_60_decimals:
1.229 +  "\<bar>exp 1 - 2.718281828459045235360287471352662497757247093699959574966968\<bar>
1.230 +      \<le> inverse (10^60::real)"
1.231 +  by (rule approx_coarsen, rule exp_1_approx[of 48])
1.232 +     (simp_all add: eval_euler_approx eval_fact)
1.233 +
1.234 +lemma euler_128:
1.235 +  "\<bar>exp 1 - 924983374546220337150911035843336795079 / 2 ^ 128\<bar> \<le> inverse (2 ^ 128 :: real)"
1.236 +  by (rule approx_coarsen[OF euler_60_decimals]) simp_all
1.237 +
1.238 +lemma euler_64:
1.239 +  "\<bar>exp 1 - 50143449209799256683 / 2 ^ 64\<bar> \<le> inverse (2 ^ 64 :: real)"
1.240 +  by (rule approx_coarsen[OF euler_128]) simp_all
1.241 +
1.242 +text \<open>
1.243 +  An approximation of $e$ to 60 decimals. This is about as far as we can go with the
1.244 +  simplifier with this kind of setup; the exported code of the code generator, on the other
1.245 +  hand, can easily approximate $e$ to 1000 decimals and verify that approximation within
1.246 +  fractions of a second.
1.247 +\<close>
1.248 +
1.249 +(* (Uncommented because we don't want to use the code generator;
1.250 +   don't forget to import Code\_Target\_Numeral)) *)
1.251 +(*
1.252 +lemma "\<bar>exp 1 - 2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952605956307381323286279434907632338298807531952510190115738341879307021540891499348841675092447614606680822648001684774118537423454424371075390777449920695517027618386062613313845830007520449338265602976067371132007093287091274437470472306969772093101416928368190255151086574637721112523897844250569536967707854499699679468644549059879316368892300987931277361782154249992295763514822082698951936680331825288693984964651058209392398294887933203625094431173012381970684161403970198376793206832823764648042953118023287825098194558153017567173613320698112509961818815930416903515988885193458072738667385894228792284998920868058257492796104841984443634632449684875602336248270419786232090021609902353043699418491463140934317381436405462531520961836908887070167683964243781405927145635490613031072085103837505101157477041718986106873969655212671546889570350354021\<bar>
1.253 +  \<le> inverse (10^1000::real)"
1.254 +  by (rule approx_coarsen, rule exp_1_approx[of 450], simp) eval
1.255 +*)
1.256 +
1.257 +
1.258 +subsection \<open>Approximation of $\ln 2$\<close>
1.259 +
1.260 +text \<open>
1.261 +  The following three auxiliary constants allow us to force the simplifier to
1.262 +  evaluate intermediate results, simulating call-by-value.
1.263 +\<close>
1.264 +
1.265 +definition "ln_approx_aux3 x' e n y d \<longleftrightarrow>
1.266 +  \<bar>(2 * y) * (\<Sum>k<n. inverse (real (2*k+1)) * (y^2)^k) + d - x'\<bar> \<le> e - d"
1.267 +definition "ln_approx_aux2 x' e n y \<longleftrightarrow>
1.268 +  ln_approx_aux3 x' e n y (y^(2*n+1) / (1 - y^2) / real (2*n+1))"
1.269 +definition "ln_approx_aux1 x' e n x \<longleftrightarrow>
1.270 +  ln_approx_aux2 x' e n ((x - 1) / (x + 1))"
1.271 +
1.272 +lemma ln_approx_abs'':
1.273 +  fixes x :: real and n :: nat
1.274 +  defines "y \<equiv> (x-1)/(x+1)"
1.275 +  defines "approx \<equiv> (\<Sum>k<n. 2*y^(2*k+1) / of_nat (2*k+1))"
1.276 +  defines "d \<equiv> y^(2*n+1) / (1 - y^2) / of_nat (2*n+1)"
1.277 +  assumes x: "x > 1"
1.278 +  assumes A: "ln_approx_aux1 x' e n x"
1.279 +  shows   "\<bar>ln x - x'\<bar> \<le> e"
1.280 +proof (rule approx_coarsen[OF ln_approx_abs[OF x, of n]], goal_cases)
1.281 +  case 1
1.282 +  from A have "\<bar>2 * y * (\<Sum>k<n. inverse (real (2 * k + 1)) * y\<^sup>2 ^ k) + d - x'\<bar> \<le> e - d"
1.283 +    by (simp only: ln_approx_aux3_def ln_approx_aux2_def ln_approx_aux1_def
1.284 +                   y_def [symmetric] d_def [symmetric])
1.285 +  also have "2 * y * (\<Sum>k<n. inverse (real (2 * k + 1)) * y\<^sup>2 ^ k) =
1.286 +               (\<Sum>k<n. 2 * y^(2*k+1) / (real (2 * k + 1)))"
1.287 +    by (subst setsum_right_distrib, simp, subst power_mult)
1.288 +       (simp_all add: divide_simps mult_ac power_mult)
1.289 +  finally show ?case by (simp only: d_def y_def approx_def)
1.290 +qed
1.291 +
1.292 +text \<open>
1.293 +  We unfold the above three constants successively and then compute the
1.294 +  sum using a Horner scheme.
1.295 +\<close>
1.296 +lemma ln_2_40_decimals:
1.297 +  "\<bar>ln 2 - 0.6931471805599453094172321214581765680755\<bar>
1.298 +      \<le> inverse (10^40 :: real)"
1.299 +  apply (rule ln_approx_abs''[where n = 40], simp)
1.300 +  apply (simp, simp add: ln_approx_aux1_def)
1.301 +  apply (simp add: ln_approx_aux2_def power2_eq_square power_divide)
1.302 +  apply (simp add: ln_approx_aux3_def power2_eq_square)
1.303 +  apply (simp add: setsum_poly_horner_expand)
1.304 +  done
1.305 +
1.306 +lemma ln_2_128:
1.307 +  "\<bar>ln 2 - 235865763225513294137944142764154484399 / 2 ^ 128\<bar> \<le> inverse (2 ^ 128 :: real)"
1.308 +  by (rule approx_coarsen[OF ln_2_40_decimals]) simp_all
1.309 +
1.310 +lemma ln_2_64:
1.311 +  "\<bar>ln 2 - 12786308645202655660 / 2 ^ 64\<bar> \<le> inverse (2 ^ 64 :: real)"
1.312 +  by (rule approx_coarsen[OF ln_2_128]) simp_all
1.313 +
1.314 +
1.315 +
1.316 +subsection \<open>Approximation of the Euler--Mascheroni constant\<close>
1.317 +
1.318 +text \<open>
1.319 +  Unfortunatly, the best approximation we have formalised for the Euler--Mascheroni
1.320 +  constant converges only quadratically. This is too slow to compute more than a
1.321 +  few decimals, but we can get almost 4 decimals / 14 binary digits this way,
1.322 +  which is not too bad.
1.323 +\<close>
1.324 +lemma euler_mascheroni_approx:
1.325 +  defines "approx \<equiv> 0.577257 :: real" and "e \<equiv> 0.000063 :: real"
1.326 +  shows   "abs (euler_mascheroni - approx :: real) < e"
1.327 +  (is "abs (_ - ?approx) < ?e")
1.328 +proof -
1.329 +  define l :: real
1.330 +    where "l = 47388813395531028639296492901910937/82101866951584879688289000000000000"
1.331 +  define u :: real
1.332 +    where "u = 142196984054132045946501548559032969 / 246305600854754639064867000000000000"
1.333 +  have impI: "P \<longrightarrow> Q" if Q for P Q using that by blast
1.334 +  have hsum_63: "harm 63 = (310559566510213034489743057 / 65681493561267903750631200 :: real)"
1.335 +    by (simp add: harm_expand)
1.336 +  from harm_Suc[of 63] have hsum_64: "harm 64 =
1.337 +          623171679694215690971693339 / (131362987122535807501262400::real)"
1.338 +    by (subst (asm) hsum_63) simp
1.339 +  have "ln (64::real) = real (6::nat) * ln 2" by (subst ln_realpow[symmetric]) simp_all
1.340 +  hence "ln (real_of_nat (Suc 63)) \<in> {4.158883083293<..<4.158883083367}" using ln_2_64
1.341 +    by (simp add: abs_real_def split: if_split_asm)
1.342 +  from euler_mascheroni_bounds'[OF _ this]
1.343 +    have "(euler_mascheroni :: real) \<in> {l<..<u}"
1.344 +    by (simp add: hsum_63 del: greaterThanLessThan_iff) (simp only: l_def u_def)
1.345 +  also have "\<dots> \<subseteq> {approx - e<..<approx + e}"
1.346 +    by (subst greaterThanLessThan_subseteq_greaterThanLessThan, rule impI)
1.347 +       (simp add: approx_def e_def u_def l_def)
1.348 +  finally show ?thesis by (simp add: abs_real_def)
1.349 +qed
1.350 +
1.351 +
1.352 +
1.353 +subsection \<open>Approximation of pi\<close>
1.354 +
1.355 +
1.356 +subsubsection \<open>Approximating the arctangent\<close>
1.357 +
1.358 +text\<open>
1.359 +  The arctangent can be used to approximate pi. Fortunately, its Taylor series expansion
1.360 +  converges exponentially for small values, so we can get $\Theta(n)$ digits of precision
1.361 +  with $n$ summands of the expansion.
1.362 +\<close>
1.363 +
1.364 +definition arctan_approx where
1.365 +  "arctan_approx n x = x * (\<Sum>k<n. (-(x^2))^k / real (2*k+1))"
1.366 +
1.367 +lemma arctan_series':
1.368 +  assumes "\<bar>x\<bar> \<le> 1"
1.369 +  shows "(\<lambda>k. (-1)^k * (1 / real (k*2+1) * x ^ (k*2+1))) sums arctan x"
1.370 +  using summable_arctan_series[OF assms] arctan_series[OF assms] by (simp add: sums_iff)
1.371 +
1.372 +lemma arctan_approx:
1.373 +  assumes x: "0 \<le> x" "x < 1" and n: "even n"
1.374 +  shows   "arctan x - arctan_approx n x \<in> {0..x^(2*n+1) / (1-x^4)}"
1.375 +proof -
1.376 +  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
1.377 +  from assms have "(\<lambda>k. (-1) ^ k * (1 / real (k * 2 + 1) * x^(k*2+1))) sums arctan x"
1.378 +    using arctan_series' by simp
1.379 +  also have "(\<lambda>k. (-1) ^ k * (1 / real (k * 2 + 1) * x^(k*2+1))) =
1.380 +                 (\<lambda>k. x * ((- (x^2))^k / real (2*k+1)))"
1.381 +    by (simp add: power2_eq_square power_mult power_mult_distrib mult_ac power_minus')
1.382 +  finally have "(\<lambda>k. x * ((- x\<^sup>2) ^ k / real (2 * k + 1))) sums arctan x" .
1.383 +  from sums_split_initial_segment[OF this, of n]
1.384 +    have "(\<lambda>i. x * ((- x\<^sup>2) ^ (i + n) / real (2 * (i + n) + 1))) sums
1.385 +            (arctan x - arctan_approx n x)"
1.386 +    by (simp add: arctan_approx_def setsum_right_distrib)
1.387 +  from sums_group[OF this, of 2] assms
1.388 +    have sums: "(\<lambda>k. x * (x\<^sup>2)^n * (x^4)^k * c k) sums (arctan x - arctan_approx n x)"
1.390 +
1.391 +  from assms have "0 \<le> arctan x - arctan_approx n x"
1.392 +    by (intro sums_le[OF _ sums_zero sums] allI mult_nonneg_nonneg)
1.393 +       (auto intro!: frac_le power_le_one simp: c_def)
1.394 +  moreover {
1.395 +    from assms have "c k \<le> 1 - 0" for k unfolding c_def
1.396 +      by (intro diff_mono divide_nonneg_nonneg add_nonneg_nonneg) auto
1.397 +    with assms have "x * x\<^sup>2 ^ n * (x ^ 4) ^ k * c k \<le> x * x\<^sup>2 ^ n * (x ^ 4) ^ k * 1" for k
1.398 +      by (intro mult_left_mono mult_right_mono mult_nonneg_nonneg) simp_all
1.399 +    with assms have "arctan x - arctan_approx n x \<le> x * (x\<^sup>2)^n * (1 / (1 - x^4))"
1.400 +      by (intro sums_le[OF _ sums sums_mult[OF geometric_sums]] allI mult_left_mono)
1.401 +         (auto simp: power_less_one)
1.402 +    also have "x * (x^2)^n = x^(2*n+1)" by (simp add: power_mult power_add)
1.403 +    finally have "arctan x - arctan_approx n x \<le> x^(2*n+1) / (1 - x^4)" by simp
1.404 +  }
1.405 +  ultimately show ?thesis by simp
1.406 +qed
1.407 +
1.408 +lemma arctan_approx_def': "arctan_approx n (1/x) =
1.409 +  (\<Sum>k<n. inverse (real (2 * k + 1) * (- x\<^sup>2) ^ k)) / x"
1.410 +proof -
1.411 +  have "(-1)^k / b = 1 / ((-1)^k * b)" for k :: nat and b :: real
1.412 +    by (cases "even k") auto
1.413 +  thus ?thesis by (simp add: arctan_approx_def  field_simps power_minus')
1.414 +qed
1.415 +
1.416 +lemma expand_arctan_approx:
1.417 +  "(\<Sum>k<(numeral n::nat). inverse (f k) * inverse (x ^ k)) =
1.418 +     inverse (f 0) + (\<Sum>k<pred_numeral n. inverse (f (k+1)) * inverse (x^k)) / x"
1.419 +  "(\<Sum>k<Suc 0. inverse (f k) * inverse (x^k)) = inverse (f 0 :: 'a :: field)"
1.420 +  "(\<Sum>k<(0::nat). inverse (f k) * inverse (x^k)) = 0"
1.421 +proof -
1.422 +  {
1.423 +    fix m :: nat
1.424 +    have "(\<Sum>k<Suc m. inverse (f k * x^k)) =
1.425 +             inverse (f 0) + (\<Sum>k=Suc 0..<Suc m. inverse (f k * x^k))"
1.426 +      by (subst atLeast0LessThan [symmetric], subst setsum_head_upt_Suc) simp_all
1.427 +    also have "(\<Sum>k=Suc 0..<Suc m. inverse (f k * x^k)) = (\<Sum>k<m. inverse (f (k+1) * x^k)) / x"
1.428 +      by (subst setsum_shift_bounds_Suc_ivl)
1.429 +         (simp add: setsum_right_distrib divide_inverse algebra_simps
1.430 +                    atLeast0LessThan power_commutes)
1.431 +    finally have "(\<Sum>k<Suc m. inverse (f k) * inverse (x ^ k)) =
1.432 +                      inverse (f 0) + (\<Sum>k<m. inverse (f (k + 1)) * inverse (x ^ k)) / x" by simp
1.433 +  }
1.434 +  from this[of "pred_numeral n"]
1.435 +    show "(\<Sum>k<numeral n. inverse (f k) * inverse (x^k)) =
1.436 +            inverse (f 0) + (\<Sum>k<pred_numeral n. inverse (f (k+1)) * inverse (x^k)) / x"
1.437 +    by (simp add: numeral_eq_Suc)
1.438 +qed simp_all
1.439 +
1.440 +lemma arctan_diff_small:
1.441 +  assumes "\<bar>x*y::real\<bar> < 1"
1.442 +  shows   "arctan x - arctan y = arctan ((x - y) / (1 + x * y))"
1.443 +proof -
1.444 +  have "arctan x - arctan y = arctan x + arctan (-y)" by (simp add: arctan_minus)
1.445 +  also from assms have "\<dots> = arctan ((x - y) / (1 + x * y))" by (subst arctan_add_small) simp_all
1.446 +  finally show ?thesis .
1.447 +qed
1.448 +
1.449 +
1.450 +subsubsection \<open>Machin-like formulae for pi\<close>
1.451 +
1.452 +text \<open>
1.453 +  We first define a small proof method that can prove Machin-like formulae for @{term "pi"}
1.454 +  automatically. Unfortunately, this takes far too much time for larger formulae because
1.455 +  the numbers involved become too large.
1.456 +\<close>
1.457 +
1.458 +definition "MACHIN_TAG a b \<equiv> a * (b::real)"
1.459 +
1.460 +lemma numeral_horner_MACHIN_TAG:
1.461 +  "MACHIN_TAG Numeral1 x = x"
1.462 +  "MACHIN_TAG (numeral (Num.Bit0 (Num.Bit0 n))) x =
1.463 +     MACHIN_TAG 2 (MACHIN_TAG (numeral (Num.Bit0 n)) x)"
1.464 +  "MACHIN_TAG (numeral (Num.Bit0 (Num.Bit1 n))) x =
1.465 +     MACHIN_TAG 2 (MACHIN_TAG (numeral (Num.Bit1 n)) x)"
1.466 +  "MACHIN_TAG (numeral (Num.Bit1 n)) x =
1.467 +     MACHIN_TAG 2 (MACHIN_TAG (numeral n) x) + x"
1.468 +  unfolding numeral_Bit0 numeral_Bit1 ring_distribs one_add_one[symmetric] MACHIN_TAG_def
1.469 +     by (simp_all add: algebra_simps)
1.470 +
1.471 +lemma tag_machin: "a * arctan b = MACHIN_TAG a (arctan b)" by (simp add: MACHIN_TAG_def)
1.472 +
1.473 +lemma arctan_double': "\<bar>a::real\<bar> < 1 \<Longrightarrow> MACHIN_TAG 2 (arctan a) = arctan (2 * a / (1 - a*a))"
1.474 +  unfolding MACHIN_TAG_def by (simp add: arctan_double power2_eq_square)
1.475 +
1.476 +ML \<open>
1.477 +  fun machin_term_conv ctxt ct =
1.478 +    let
1.480 +    in
1.481 +      case Thm.term_of ct of
1.482 +        Const (@{const_name MACHIN_TAG}, _) $_$
1.483 +          (Const (@{const_name "Transcendental.arctan"}, _) $_) => 1.484 + Simplifier.rewrite ctxt' ct 1.485 + | 1.486 + Const (@{const_name MACHIN_TAG}, _)$ _ $1.487 + (Const (@{const_name "Groups.plus"}, _)$
1.488 +            (Const (@{const_name "Transcendental.arctan"}, _) $_)$
1.489 +            (Const (@{const_name "Transcendental.arctan"}, _) $_)) => 1.490 + Simplifier.rewrite ctxt' ct 1.491 + | _ => raise CTERM ("machin_conv", [ct]) 1.492 + end 1.493 + 1.494 + fun machin_tac ctxt = 1.495 + let val conv = Conv.top_conv (Conv.try_conv o machin_term_conv) ctxt 1.496 + in 1.497 + SELECT_GOAL ( 1.498 + Local_Defs.unfold_tac ctxt 1.499 + @{thms tag_machin[THEN eq_reflection] numeral_horner_MACHIN_TAG[THEN eq_reflection]} 1.500 + THEN REPEAT (CHANGED (HEADGOAL (CONVERSION conv)))) 1.501 + THEN' Simplifier.simp_tac (ctxt addsimps @{thms arctan_add_small arctan_diff_small}) 1.502 + end 1.503 +\<close> 1.504 + 1.505 +method_setup machin = \<open>Scan.succeed (SIMPLE_METHOD' o machin_tac)\<close> 1.506 + 1.507 +text \<open> 1.508 + We can now prove the standard'' Machin formula, which was already proven manually 1.509 + in Isabelle, automatically. 1.510 +}\<close> 1.511 +lemma "pi / 4 = (4::real) * arctan (1 / 5) - arctan (1 / 239)" 1.512 + by machin 1.513 + 1.514 +text \<open> 1.515 + We can also prove the following more complicated formula: 1.516 +\<close> 1.517 +lemma machin': "pi/4 = (12::real) * arctan (1/18) + 8 * arctan (1/57) - 5 * arctan (1/239)" 1.518 + by machin 1.519 + 1.520 + 1.521 + 1.522 +subsubsection \<open>Simple approximation of pi\<close> 1.523 + 1.524 +text \<open> 1.525 + We can use the simple Machin formula and the Taylor series expansion of the arctangent 1.526 + to approximate pi. For a given even natural number$n$, we expand @{term "arctan (1/5)"} 1.527 + to$3n$summands and @{term "arctan (1/239)"} to$n$summands. This gives us at least 1.528 +$13n-2$bits of precision. 1.529 +\<close> 1.530 + 1.531 +definition "pi_approx n = 16 * arctan_approx (3*n) (1/5) - 4 * arctan_approx n (1/239)" 1.532 + 1.533 +lemma pi_approx: 1.534 + fixes n :: nat assumes n: "even n" and "n > 0" 1.535 + shows "\<bar>pi - pi_approx n\<bar> \<le> inverse (2^(13*n - 2))" 1.536 +proof - 1.537 + from n have n': "even (3*n)" by simp 1.538 + \<comment> \<open>We apply the Machin formula\<close> 1.539 + from machin have "pi = 16 * arctan (1/5) - 4 * arctan (1/239::real)" by simp 1.540 + \<comment> \<open>Taylor series expansion of the arctangent\<close> 1.541 + also from arctan_approx[OF _ _ n', of "1/5"] arctan_approx[OF _ _ n, of "1/239"] 1.542 + have "\<dots> - pi_approx n \<in> {-4*((1/239)^(2*n+1) / (1-(1/239)^4))..16*(1/5)^(6*n+1) / (1-(1/5)^4)}" 1.543 + by (simp add: pi_approx_def) 1.544 + \<comment> \<open>Coarsening the bounds to make them a bit nicer\<close> 1.545 + also have "-4*((1/239::real)^(2*n+1) / (1-(1/239)^4)) = -((13651919 / 815702160) / 57121^n)" 1.546 + by (simp add: power_mult power2_eq_square) (simp add: field_simps) 1.547 + also have "16*(1/5)^(6*n+1) / (1-(1/5::real)^4) = (125/39) / 15625^n" 1.548 + by (simp add: power_mult power2_eq_square) (simp add: field_simps) 1.549 + also have "{-((13651919 / 815702160) / 57121^n) .. (125 / 39) / 15625^n} \<subseteq> 1.550 + {- (4 / 2^(13*n)) .. 4 / (2^(13*n)::real)}" 1.551 + by (subst atLeastatMost_subset_iff, intro disjI2 conjI le_imp_neg_le) 1.552 + (rule frac_le; simp add: power_mult power_mono)+ 1.553 + finally have "abs (pi - pi_approx n) \<le> 4 / 2^(13*n)" by auto 1.554 + also from \<open>n > 0\<close> have "4 / 2^(13*n) = 1 / (2^(13*n - 2) :: real)" 1.555 + by (cases n) (simp_all add: power_add) 1.556 + finally show ?thesis by (simp add: divide_inverse) 1.557 +qed 1.558 + 1.559 +lemma pi_approx': 1.560 + fixes n :: nat assumes n: "even n" and "n > 0" and "k \<le> 13*n - 2" 1.561 + shows "\<bar>pi - pi_approx n\<bar> \<le> inverse (2^k)" 1.562 + using assms(3) by (intro order.trans[OF pi_approx[OF assms(1,2)]]) (simp_all add: field_simps) 1.563 + 1.564 +text \<open>We can now approximate pi to 22 decimals within a fraction of a second.\<close> 1.565 +lemma pi_approx_75: "abs (pi - 3.1415926535897932384626 :: real) \<le> inverse (10^22)" 1.566 +proof - 1.567 + define a :: real 1.568 + where "a = 8295936325956147794769600190539918304 / 2626685325478320010006427764892578125" 1.569 + define b :: real 1.570 + where "b = 8428294561696506782041394632 / 503593538783547230635598424135" 1.571 + \<comment> \<open>The introduction of this constant prevents the simplifier from applying solvers that 1.572 + we don't want. We want it to simply evaluate the terms to rational constants.}\<close> 1.573 + define eq :: "real \<Rightarrow> real \<Rightarrow> bool" where "eq = op =" 1.574 + 1.575 + \<comment> \<open>Splitting the computation into several steps has the advantage that simplification can 1.576 + be done in parallel\<close> 1.577 + have "abs (pi - pi_approx 6) \<le> inverse (2^76)" by (rule pi_approx') simp_all 1.578 + also have "pi_approx 6 = 16 * arctan_approx (3 * 6) (1 / 5) - 4 * arctan_approx 6 (1 / 239)" 1.579 + unfolding pi_approx_def by simp 1.580 + also have [unfolded eq_def]: "eq (16 * arctan_approx (3 * 6) (1 / 5)) a" 1.581 + by (simp add: arctan_approx_def' power2_eq_square, 1.582 + simp add: expand_arctan_approx, unfold a_def eq_def, rule refl) 1.583 + also have [unfolded eq_def]: "eq (4 * arctan_approx 6 (1 / 239::real)) b" 1.584 + by (simp add: arctan_approx_def' power2_eq_square, 1.585 + simp add: expand_arctan_approx, unfold b_def eq_def, rule refl) 1.586 + also have [unfolded eq_def]: 1.587 + "eq (a - b) (171331331860120333586637094112743033554946184594977368554649608 / 1.588 + 54536456744112171868276045488779391002026386559009552001953125)" 1.589 + by (unfold a_def b_def, simp, unfold eq_def, rule refl) 1.590 + finally show ?thesis by (rule approx_coarsen) simp 1.591 +qed 1.592 + 1.593 +text \<open> 1.594 + The previous estimate of pi in this file was based on approximating the root of the 1.595 +$\sin(\pi/6)$in the interval$[0;4]$using the Taylor series expansion of the sine to 1.596 + verify that it is between two given bounds. 1.597 + This was much slower and much less precise. We can easily recover this coarser estimate from 1.598 + the newer, precise estimate: 1.599 +\<close> 1.600 +lemma pi_approx_32: "\<bar>pi - 13493037705/4294967296 :: real\<bar> \<le> inverse(2 ^ 32)" 1.601 + by (rule approx_coarsen[OF pi_approx_75]) simp 1.602 + 1.603 + 1.604 +subsection \<open>A more complicated approximation of pi\<close> 1.605 + 1.606 +text \<open> 1.607 + There are more complicated Machin-like formulae that have more terms with larger 1.608 + denominators. Although they have more terms, each term requires fewer summands of the 1.609 + Taylor series for the same precision, since it is evaluated closer to$0\$.
1.610 +
1.611 +  Using a good formula, one can therefore obtain the same precision with fewer operations.
1.612 +  The big formulae used for computations of pi in practice are too complicated for us to
1.613 +  prove here, but we can use the three-term Machin-like formula @{thm machin'}.
1.614 +\<close>
1.615 +
1.616 +definition "pi_approx2 n = 48 * arctan_approx (6*n) (1/18::real) +
1.617 +                             32 * arctan_approx (4*n) (1/57) - 20 * arctan_approx (3*n) (1/239)"
1.618 +
1.619 +lemma pi_approx2:
1.620 +  fixes n :: nat assumes n: "even n" and "n > 0"
1.621 +  shows   "\<bar>pi - pi_approx2 n\<bar> \<le> inverse (2^(46*n - 1))"
1.622 +proof -
1.623 +  from n have n': "even (6*n)" "even (4*n)" "even (3*n)" by simp_all
1.624 +  from machin' have "pi = 48 * arctan (1/18) + 32 * arctan (1/57) - 20 * arctan (1/239::real)"
1.625 +    by simp
1.626 +  hence "pi - pi_approx2 n = 48 * (arctan (1/18) - arctan_approx (6*n) (1/18)) +
1.627 +                                 32 * (arctan (1/57) - arctan_approx (4*n) (1/57)) -
1.628 +                                 20 * (arctan (1/239) - arctan_approx (3*n) (1/239))"
1.629 +    by (simp add: pi_approx2_def)
1.630 +  also have "\<dots> \<in> {-((20/239/(1-(1/239)^4)) * (1/239)^(6*n))..
1.631 +              (48/18 / (1-(1/18)^4))*(1/18)^(12*n) + (32/57/(1-(1/57)^4)) * (1/57)^(8*n)}"
1.632 +    (is "_ \<in> {-?l..?u1 + ?u2}")
1.633 +    apply ((rule combine_bounds(1,2))+; (rule combine_bounds(3); (rule arctan_approx)?)?)
1.634 +    apply (simp_all add: n)
1.635 +    apply (simp_all add: divide_simps)?
1.636 +    done
1.637 +  also {
1.638 +    have "?l \<le> (1/8) * (1/2)^(46*n)"
1.639 +      unfolding power_mult by (intro mult_mono power_mono) (simp_all add: divide_simps)
1.640 +    also have "\<dots> \<le> (1/2) ^ (46 * n - 1)"
1.642 +    finally have "?l \<le> (1/2) ^ (46 * n - 1)" .
1.643 +    moreover {
1.644 +      have "?u1 + ?u2 \<le> 4 * (1/2)^(48*n) + 1 * (1/2)^(46*n)"
1.645 +        unfolding power_mult by (intro add_mono mult_mono power_mono) (simp_all add: divide_simps)
1.646 +      also from \<open>n > 0\<close> have "4 * (1/2::real)^(48*n) \<le> (1/2)^(46*n)"
1.648 +      also from \<open>n > 0\<close> have "(1/2::real) ^ (46 * n) + 1 * (1 / 2) ^ (46 * n) = (1/2) ^ (46 * n - 1)"
1.650 +      finally have "?u1 + ?u2 \<le> (1/2) ^ (46 * n - 1)" by - simp
1.651 +    }
1.652 +    ultimately have "{-?l..?u1 + ?u2} \<subseteq> {-((1/2)^(46*n-1))..(1/2)^(46*n-1)}"
1.653 +      by (subst atLeastatMost_subset_iff) simp_all
1.654 +  }
1.655 +  finally have "\<bar>pi - pi_approx2 n\<bar> \<le> ((1/2) ^ (46 * n - 1))" by auto
1.656 +  thus ?thesis by (simp add: divide_simps)
1.657 +qed
1.658 +
1.659 +lemma pi_approx2':
1.660 +  fixes n :: nat assumes n: "even n" and "n > 0" and "k \<le> 46*n - 1"
1.661 +  shows   "\<bar>pi - pi_approx2 n\<bar> \<le> inverse (2^k)"
1.662 +  using assms(3) by (intro order.trans[OF pi_approx2[OF assms(1,2)]]) (simp_all add: field_simps)
1.663 +
1.664 +text \<open>
1.665 +  We can now approximate pi to 54 decimals using this formula. The computations are much
1.666 +  slower now; this is mostly because we use arbitrary-precision rational numbers, whose
1.667 +  numerators and demoninators get very large. Using dyadic floating point numbers would be
1.668 +  much more economical.
1.669 +\<close>
1.670 +lemma pi_approx_54_decimals:
1.671 +  "abs (pi - 3.141592653589793238462643383279502884197169399375105821 :: real) \<le> inverse (10^54)"
1.672 +  (is "abs (pi - ?pi') \<le> _")
1.673 +proof -
1.674 +  define a :: real
1.675 +    where "a = 2829469759662002867886529831139137601191652261996513014734415222704732791803 /
1.676 +           1062141879292765061960538947347721564047051545995266466660439319087625011200"
1.677 +  define b :: real
1.678 +    where "b = 13355545553549848714922837267299490903143206628621657811747118592 /
1.679 +           23792006023392488526789546722992491355941103837356113731091180925"
1.680 +  define c :: real
1.681 +    where "c = 28274063397213534906669125255762067746830085389618481175335056 /
1.682 +           337877029279505250241149903214554249587517250716358486542628059"
1.683 +  let ?pi'' = "3882327391761098513316067116522233897127356523627918964967729040413954225768920394233198626889767468122598417405434625348404038165437924058179155035564590497837027530349 /
1.684 +               1235783190199688165469648572769847552336447197542738425378629633275352407743112409829873464564018488572820294102599160968781449606552922108667790799771278860366957772800"
1.685 +  define eq :: "real \<Rightarrow> real \<Rightarrow> bool" where "eq = op ="
1.686 +
1.687 +  have "abs (pi - pi_approx2 4) \<le> inverse (2^183)" by (rule pi_approx2') simp_all
1.688 +  also have "pi_approx2 4 = 48 * arctan_approx 24 (1 / 18) +
1.689 +                            32 * arctan_approx 16 (1 / 57) -
1.690 +                            20 * arctan_approx 12 (1 / 239)"
1.691 +    unfolding pi_approx2_def by simp
1.692 +  also have [unfolded eq_def]: "eq (48 * arctan_approx 24 (1 / 18)) a"
1.693 +    by (simp add: arctan_approx_def' power2_eq_square,
1.694 +        simp add: expand_arctan_approx, unfold a_def eq_def, rule refl)
1.695 +  also have [unfolded eq_def]: "eq (32 * arctan_approx 16 (1 / 57::real)) b"
1.696 +    by (simp add: arctan_approx_def' power2_eq_square,
1.697 +        simp add: expand_arctan_approx, unfold b_def eq_def, rule refl)
1.698 +  also have [unfolded eq_def]: "eq (20 * arctan_approx 12 (1 / 239::real)) c"
1.699 +    by (simp add: arctan_approx_def' power2_eq_square,
1.700 +        simp add: expand_arctan_approx, unfold c_def eq_def, rule refl)
1.701 +  also have [unfolded eq_def]:
1.702 +    "eq (a + b) (34326487387865555303797183505809267914709125998469664969258315922216638779011304447624792548723974104030355722677 /
1.703 +                 10642967245546718617684989689985787964158885991018703366677373121531695267093031090059801733340658960857196134400)"
1.704 +    by (unfold a_def b_def c_def, simp, unfold eq_def, rule refl)
1.705 +  also have [unfolded eq_def]: "eq (\<dots> - c) ?pi''"
1.706 +    by (unfold a_def b_def c_def, simp, unfold eq_def, rule refl)
1.707 +  \<comment> \<open>This is incredibly slow because the numerators and denominators are huge.\<close>
1.708 +  finally show ?thesis by (rule approx_coarsen) simp
1.709 +qed
1.710 +
1.711 +text \<open>A 128 bit approximation of pi:\<close>
1.712 +lemma pi_approx_128:
1.713 +  "abs (pi - 1069028584064966747859680373161870783301 / 2^128) \<le> inverse (2^128)"
1.714 +  by (rule approx_coarsen[OF pi_approx_54_decimals]) simp
1.715 +
1.716 +text \<open>A 64 bit approximation of pi:\<close>
1.717 +lemma pi_approx_64:
1.718 +  "abs (pi - 57952155664616982739 / 2^64 :: real) \<le> inverse (2^64)"
1.719 +  by (rule approx_coarsen[OF pi_approx_54_decimals]) simp
1.720 +
1.721 +text \<open>
1.722 +  Again, going much farther with the simplifier takes a long time, but the code generator
1.723 +  can handle even two thousand decimal digits in under 20 seconds.
1.724 +\<close>
1.725 +
1.726 +end