src/HOL/Analysis/ex/Approximations.thy
author fleury <Mathias.Fleury@mpi-inf.mpg.de>
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 \<open>Numeric approximations to Constants\<close>
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 \<open>
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
\<close>
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
  "(\<Sum>k<(numeral n::nat). f k * x^k) = f 0 + (\<Sum>k<pred_numeral n. f (k+1) * x^k) * x"
eberlm@62085
    24
  "(\<Sum>k<Suc 0. f k * x^k) = (f 0 :: 'a :: semiring_1)"
eberlm@62085
    25
  "(\<Sum>k<(0::nat). f k * x^k) = 0"
eberlm@62085
    26
proof -
eberlm@62085
    27
  {
eberlm@62085
    28
    fix m :: nat
eberlm@62085
    29
    have "(\<Sum>k<Suc m. f k * x^k) = f 0 + (\<Sum>k=Suc 0..<Suc m. f k * x^k)"
eberlm@62085
    30
      by (subst atLeast0LessThan [symmetric], subst setsum_head_upt_Suc) simp_all
eberlm@62085
    31
    also have "(\<Sum>k=Suc 0..<Suc m. f k * x^k) = (\<Sum>k<m. f (k+1) * x^k) * x"
eberlm@62085
    32
      by (subst setsum_shift_bounds_Suc_ivl)
Mathias@63918
    33
         (simp add: setsum_distrib_right algebra_simps atLeast0LessThan power_commutes)
eberlm@62085
    34
    finally have "(\<Sum>k<Suc m. f k * x ^ k) = f 0 + (\<Sum>k<m. f (k + 1) * x ^ k) * x" .
eberlm@62085
    35
  }
eberlm@62089
    36
  from this[of "pred_numeral n"]
eberlm@62089
    37
    show "(\<Sum>k<numeral n. f k * x^k) = f 0 + (\<Sum>k<pred_numeral n. f (k+1) * x^k) * x"
eberlm@62085
    38
    by (simp add: numeral_eq_Suc)
eberlm@62085
    39
qed simp_all
eberlm@62085
    40
eberlm@62089
    41
lemma power_less_one:
eberlm@62085
    42
  assumes "n > 0" "x \<ge> 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 \<in> {a1..b1} \<Longrightarrow> y \<in> {a2..b2} \<Longrightarrow> a3 = a1 + a2 \<Longrightarrow> b3 = b1 + b2 \<Longrightarrow> x + y \<in> {a3..(b3::real)}"
eberlm@62085
    56
  "x \<in> {a1..b1} \<Longrightarrow> y \<in> {a2..b2} \<Longrightarrow> a3 = a1 - b2 \<Longrightarrow> b3 = b1 - a2 \<Longrightarrow> x - y \<in> {a3..(b3::real)}"
eberlm@62085
    57
  "c \<ge> 0 \<Longrightarrow> x \<in> {a..b} \<Longrightarrow> c * x \<in> {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
  "\<bar>x - a1\<bar> \<le> eps1 \<Longrightarrow> \<bar>a1 - a2\<bar> \<le> eps2 - eps1 \<Longrightarrow> \<bar>x - a2\<bar> \<le> (eps2 :: real)"
eberlm@62085
    62
  by simp
eberlm@62085
    63
eberlm@62085
    64
eberlm@62089
    65
eberlm@62089
    66
subsection \<open>Approximations of the exponential function\<close>
eberlm@62089
    67
eberlm@62089
    68
lemma two_power_fact_le_fact:
eberlm@62089
    69
  assumes "n \<ge> 1"
eberlm@62089
    70
  shows   "2^k * fact n \<le> (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 \<le> 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) \<le> 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 "\<dots> = 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 \<open>
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
\<close>
eberlm@62089
    88
lemma exp_approx:
eberlm@62089
    89
  assumes "n > 0" "0 \<le> x" "x < 2"
eberlm@62089
    90
  shows   "exp (x::real) - (\<Sum>k<n. x^k / fact k) \<in> {0..(2 * x^n / (2 - x)) / fact n}"
eberlm@62089
    91
proof (unfold atLeastAtMost_iff, safe)
wenzelm@63040
    92
  define approx where "approx = (\<Sum>k<n. x^k / fact k)"
eberlm@62089
    93
  have "(\<lambda>k. 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: "(\<lambda>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) \<ge> 0"
eberlm@62089
   100
    by (intro sums_le[OF _ sums_zero sums]) auto
lp15@59871
   101
eberlm@62089
   102
  have "\<forall>k. x^n * (x^k / fact (n+k)) \<le> (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 "\<dots> \<le> 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 "\<dots> = (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)) \<le> (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 "(\<lambda>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 "(\<lambda>k. (x^n / fact n) * (x / 2)^k) sums \<dots>" .
eberlm@62089
   119
  }
eberlm@62089
   120
  ultimately show "(exp x - approx) \<le> (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 \<open>
eberlm@62089
   125
  The following variant gives a simpler error estimate for inputs between $0$ and $1$:  
eberlm@62089
   126
\<close>
eberlm@62089
   127
lemma exp_approx':
eberlm@62089
   128
  assumes "n > 0" "0 \<le> x" "x \<le> 1"
eberlm@62089
   129
  shows   "\<bar>exp (x::real) - (\<Sum>k\<le>n. x^k / fact k)\<bar> \<le> x ^ n / fact n"
eberlm@62085
   130
proof -
eberlm@62089
   131
  from assms have "x^n / (2 - x) \<le> x^n / 1" by (intro frac_le) simp_all 
eberlm@62089
   132
  hence "(2 * x^n / (2 - x)) / fact n \<le> 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) - (\<Sum>k<n. x^k / fact k) \<in> {0..2 * x^n / fact n}" by simp
eberlm@62089
   136
  moreover have "(\<Sum>k\<le>n. x^k / fact k) = (\<Sum>k<n. x^k / fact k) + x ^ n / fact n"
eberlm@62089
   137
    by (simp add: lessThan_Suc_atMost [symmetric])
eberlm@62089
   138
  ultimately show "\<bar>exp (x::real) - (\<Sum>k\<le>n. x^k / fact k)\<bar> \<le> x ^ n / fact n"
eberlm@62089
   139
    unfolding atLeastAtMost_iff by linarith
eberlm@62089
   140
qed
eberlm@62089
   141
eberlm@62089
   142
text \<open>
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
\<close>
eberlm@62089
   149
lemma exp_approx'':
eberlm@62089
   150
  assumes "n > 0" "0 \<le> x" "x \<le> 1"
eberlm@62089
   151
  shows   "\<bar>exp (x::real) - (\<Sum>k\<le>n. x^k / fact k)\<bar> \<le> 1 / fact n"
eberlm@62089
   152
proof -
eberlm@62089
   153
  from assms have "\<bar>exp x - (\<Sum>k\<le>n. x ^ k / fact k)\<bar> \<le> x ^ n / fact n"
eberlm@62089
   154
    by (rule exp_approx')
eberlm@62089
   155
  also from assms have "\<dots> \<le> 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 \<open>
eberlm@62089
   161
  We now define an approximation function for Euler's constant $e$.  
eberlm@62089
   162
\<close>
eberlm@62089
   163
eberlm@62089
   164
definition euler_approx :: "nat \<Rightarrow> real" where
eberlm@62089
   165
  "euler_approx n = (\<Sum>k\<le>n. inverse (fact k))"
eberlm@62089
   166
eberlm@62089
   167
definition euler_approx_aux :: "nat \<Rightarrow> nat" where
eberlm@62089
   168
  "euler_approx_aux n = (\<Sum>k\<le>n. \<Prod>{k + 1..n})"
eberlm@62089
   169
eberlm@62089
   170
lemma exp_1_approx:
eberlm@62089
   171
  "n > 0 \<Longrightarrow> \<bar>exp (1::real) - euler_approx n\<bar> \<le> 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 \<open>
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
\<close>
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 (\<Sum>k\<le>n. \<Prod>{k+1..n}) = (\<Sum>k\<le>n. \<Prod>i=k+1..n. real i)" by simp
eberlm@62089
   183
  also have "\<dots> / fact n = (\<Sum>k\<le>n. 1 / (fact n / (\<Prod>i=k+1..n. real i)))"
eberlm@62089
   184
    by (simp add: setsum_divide_distrib)
eberlm@62089
   185
  also have "\<dots> = (\<Sum>k\<le>n. 1 / fact k)"
eberlm@62089
   186
  proof (intro setsum.cong refl)
eberlm@62089
   187
    fix k assume k: "k \<in> {..n}"
haftmann@63417
   188
    have "fact n = (\<Prod>i=1..n. real i)" by (simp add: fact_setprod)
eberlm@62089
   189
    also from k have "{1..n} = {1..k} \<union> {k+1..n}" by auto
eberlm@62089
   190
    also have "setprod real \<dots> / (\<Prod>i=k+1..n. real i) = (\<Prod>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 "\<dots> = 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 "\<dots> = 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 \<open>Approximations of $e$ to 60 decimals / 128 and 64 bits:\<close>
eberlm@62089
   224
eberlm@62089
   225
lemma euler_60_decimals:
eberlm@62089
   226
  "\<bar>exp 1 - 2.718281828459045235360287471352662497757247093699959574966968\<bar> 
eberlm@62089
   227
      \<le> 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
  "\<bar>exp 1 - 924983374546220337150911035843336795079 / 2 ^ 128\<bar> \<le> 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
  "\<bar>exp 1 - 50143449209799256683 / 2 ^ 64\<bar> \<le> inverse (2 ^ 64 :: real)"
eberlm@62089
   237
  by (rule approx_coarsen[OF euler_128]) simp_all
eberlm@62089
   238
eberlm@62089
   239
text \<open>
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
\<close>
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 "\<bar>exp 1 - 2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952605956307381323286279434907632338298807531952510190115738341879307021540891499348841675092447614606680822648001684774118537423454424371075390777449920695517027618386062613313845830007520449338265602976067371132007093287091274437470472306969772093101416928368190255151086574637721112523897844250569536967707854499699679468644549059879316368892300987931277361782154249992295763514822082698951936680331825288693984964651058209392398294887933203625094431173012381970684161403970198376793206832823764648042953118023287825098194558153017567173613320698112509961818815930416903515988885193458072738667385894228792284998920868058257492796104841984443634632449684875602336248270419786232090021609902353043699418491463140934317381436405462531520961836908887070167683964243781405927145635490613031072085103837505101157477041718986106873969655212671546889570350354021\<bar>
eberlm@62089
   250
  \<le> 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 \<open>Approximation of $\ln 2$\<close>
eberlm@62089
   256
eberlm@62089
   257
text \<open>
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
\<close>
eberlm@62089
   261
eberlm@62089
   262
definition "ln_approx_aux3 x' e n y d \<longleftrightarrow> 
eberlm@62089
   263
  \<bar>(2 * y) * (\<Sum>k<n. inverse (real (2*k+1)) * (y^2)^k) + d - x'\<bar> \<le> e - d"
eberlm@62089
   264
definition "ln_approx_aux2 x' e n y \<longleftrightarrow> 
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 \<longleftrightarrow> 
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 \<equiv> (x-1)/(x+1)"
eberlm@62089
   272
  defines "approx \<equiv> (\<Sum>k<n. 2*y^(2*k+1) / of_nat (2*k+1))"
eberlm@62089
   273
  defines "d \<equiv> 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   "\<bar>ln x - x'\<bar> \<le> 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 "\<bar>2 * y * (\<Sum>k<n. inverse (real (2 * k + 1)) * y\<^sup>2 ^ k) + d - x'\<bar> \<le> 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 * (\<Sum>k<n. inverse (real (2 * k + 1)) * y\<^sup>2 ^ k) = 
eberlm@62089
   283
               (\<Sum>k<n. 2 * y^(2*k+1) / (real (2 * k + 1)))"
Mathias@63918
   284
    by (subst setsum_distrib_left, simp, subst power_mult) 
eberlm@62089
   285
       (simp_all add: divide_simps mult_ac power_mult)
eberlm@62089
   286
  finally show ?case by (simp only: d_def y_def approx_def) 
eberlm@62089
   287
qed
eberlm@62089
   288
eberlm@62089
   289
text \<open>
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
\<close>
eberlm@62089
   293
lemma ln_2_40_decimals: 
eberlm@62089
   294
  "\<bar>ln 2 - 0.6931471805599453094172321214581765680755\<bar> 
eberlm@62089
   295
      \<le> 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
  "\<bar>ln 2 - 235865763225513294137944142764154484399 / 2 ^ 128\<bar> \<le> 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
  "\<bar>ln 2 - 12786308645202655660 / 2 ^ 64\<bar> \<le> 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 \<open>Approximation of the Euler--Mascheroni constant\<close>
eberlm@62085
   314
eberlm@62089
   315
text \<open>
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
\<close>
eberlm@62089
   321
lemma euler_mascheroni_approx:
eberlm@62085
   322
  defines "approx \<equiv> 0.577257 :: real" and "e \<equiv> 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 \<longrightarrow> 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)) \<in> {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) \<in> {l<..<u}"
eberlm@62085
   341
    by (simp add: hsum_63 del: greaterThanLessThan_iff) (simp only: l_def u_def)
eberlm@62085
   342
  also have "\<dots> \<subseteq> {approx - e<..<approx + e}"
eberlm@62089
   343
    by (subst greaterThanLessThan_subseteq_greaterThanLessThan, rule impI)
eberlm@62085
   344
       (simp add: approx_def e_def u_def l_def)
eberlm@62085
   345
  finally show ?thesis by (simp add: abs_real_def)
eberlm@62085
   346
qed
eberlm@62085
   347
eberlm@62085
   348
eberlm@62089
   349
eberlm@62089
   350
subsection \<open>Approximation of pi\<close>
eberlm@62085
   351
eberlm@62085
   352
eberlm@62085
   353
subsubsection \<open>Approximating the arctangent\<close>
eberlm@62085
   354
eberlm@62089
   355
text\<open>
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
\<close>
eberlm@62089
   360
eberlm@62085
   361
definition arctan_approx where
eberlm@62085
   362
  "arctan_approx n x = x * (\<Sum>k<n. (-(x^2))^k / real (2*k+1))"
eberlm@62085
   363
eberlm@62085
   364
lemma arctan_series':
eberlm@62085
   365
  assumes "\<bar>x\<bar> \<le> 1"
eberlm@62085
   366
  shows "(\<lambda>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 \<le> x" "x < 1" and n: "even n"
eberlm@62085
   371
  shows   "arctan x - arctan_approx n x \<in> {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 "(\<lambda>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 "(\<lambda>k. (-1) ^ k * (1 / real (k * 2 + 1) * x^(k*2+1))) =
eberlm@62085
   377
                 (\<lambda>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 "(\<lambda>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 "(\<lambda>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: "(\<lambda>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 \<le> 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 \<le> 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 \<le> 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 \<le> 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 \<le> 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
  (\<Sum>k<n. inverse (real (2 * k + 1) * (- x\<^sup>2) ^ 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
  "(\<Sum>k<(numeral n::nat). inverse (f k) * inverse (x ^ k)) =
eberlm@62085
   415
     inverse (f 0) + (\<Sum>k<pred_numeral n. inverse (f (k+1)) * inverse (x^k)) / x"
eberlm@62085
   416
  "(\<Sum>k<Suc 0. inverse (f k) * inverse (x^k)) = inverse (f 0 :: 'a :: field)"
eberlm@62085
   417
  "(\<Sum>k<(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 "(\<Sum>k<Suc m. inverse (f k * x^k)) =
eberlm@62085
   422
             inverse (f 0) + (\<Sum>k=Suc 0..<Suc m. inverse (f k * x^k))"
eberlm@62085
   423
      by (subst atLeast0LessThan [symmetric], subst setsum_head_upt_Suc) simp_all
eberlm@62085
   424
    also have "(\<Sum>k=Suc 0..<Suc m. inverse (f k * x^k)) = (\<Sum>k<m. inverse (f (k+1) * x^k)) / x"
eberlm@62085
   425
      by (subst setsum_shift_bounds_Suc_ivl)
Mathias@63918
   426
         (simp add: setsum_distrib_left divide_inverse algebra_simps
eberlm@62085
   427
                    atLeast0LessThan power_commutes)
eberlm@62089
   428
    finally have "(\<Sum>k<Suc m. inverse (f k) * inverse (x ^ k)) =
eberlm@62085
   429
                      inverse (f 0) + (\<Sum>k<m. inverse (f (k + 1)) * inverse (x ^ k)) / x" by simp
eberlm@62085
   430
  }
eberlm@62089
   431
  from this[of "pred_numeral n"]
eberlm@62089
   432
    show "(\<Sum>k<numeral n. inverse (f k) * inverse (x^k)) =
eberlm@62089
   433
            inverse (f 0) + (\<Sum>k<pred_numeral n. inverse (f (k+1)) * inverse (x^k)) / x"
eberlm@62085
   434
    by (simp add: numeral_eq_Suc)
eberlm@62085
   435
qed simp_all
eberlm@62085
   436
eberlm@62089
   437
lemma arctan_diff_small:
eberlm@62085
   438
  assumes "\<bar>x*y::real\<bar> < 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 "\<dots> = 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 \<open>Machin-like formulae for pi\<close>
eberlm@62085
   448
eberlm@62085
   449
text \<open>
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
\<close>
eberlm@62085
   454
eberlm@62085
   455
definition "MACHIN_TAG a b \<equiv> 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': "\<bar>a::real\<bar> < 1 \<Longrightarrow> 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 \<open>
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
\<close>
eberlm@62085
   501
eberlm@62085
   502
method_setup machin = \<open>Scan.succeed (SIMPLE_METHOD' o machin_tac)\<close>
eberlm@62085
   503
eberlm@62085
   504
text \<open>
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
}\<close>
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 \<open>
eberlm@62085
   512
  We can also prove the following more complicated formula:
eberlm@62085
   513
\<close>
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 \<open>Simple approximation of pi\<close>
eberlm@62085
   520
eberlm@62085
   521
text \<open>
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
\<close>
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   "\<bar>pi - pi_approx n\<bar> \<le> inverse (2^(13*n - 2))"
eberlm@62085
   533
proof -
eberlm@62085
   534
  from n have n': "even (3*n)" by simp
wenzelm@62175
   535
  \<comment> \<open>We apply the Machin formula\<close>
eberlm@62085
   536
  from machin have "pi = 16 * arctan (1/5) - 4 * arctan (1/239::real)" by simp
wenzelm@62175
   537
  \<comment> \<open>Taylor series expansion of the arctangent\<close>
eberlm@62085
   538
  also from arctan_approx[OF _ _ n', of "1/5"] arctan_approx[OF _ _ n, of "1/239"]
eberlm@62085
   539
    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)}"
eberlm@62085
   540
    by (simp add: pi_approx_def)
wenzelm@62175
   541
  \<comment> \<open>Coarsening the bounds to make them a bit nicer\<close>
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} \<subseteq>
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) \<le> 4 / 2^(13*n)" by auto
eberlm@62085
   551
  also from \<open>n > 0\<close> 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 \<le> 13*n - 2"
eberlm@62085
   558
  shows   "\<bar>pi - pi_approx n\<bar> \<le> 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 \<open>We can now approximate pi to 22 decimals within a fraction of a second.\<close>
eberlm@62085
   562
lemma pi_approx_75: "abs (pi - 3.1415926535897932384626 :: real) \<le> 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
  \<comment> \<open>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.}\<close>
wenzelm@63040
   570
  define eq :: "real \<Rightarrow> real \<Rightarrow> bool" where "eq = op ="
eberlm@62089
   571
wenzelm@62175
   572
  \<comment> \<open>Splitting the computation into several steps has the advantage that simplification can
eberlm@62085
   573
      be done in parallel\<close>
eberlm@62085
   574
  have "abs (pi - pi_approx 6) \<le> 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 \<open>
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
\<close>
eberlm@62085
   597
lemma pi_approx_32: "\<bar>pi - 13493037705/4294967296 :: real\<bar> \<le> 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 \<open>A more complicated approximation of pi\<close>
eberlm@62085
   602
eberlm@62085
   603
text \<open>
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
\<close>
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   "\<bar>pi - pi_approx2 n\<bar> \<le> 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 "\<dots> \<in> {-((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 "_ \<in> {-?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 \<le> (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 "\<dots> \<le> (1/2) ^ (46 * n - 1)"
eberlm@62085
   638
      by (cases n; simp_all add: power_add divide_simps)
eberlm@62085
   639
    finally have "?l \<le> (1/2) ^ (46 * n - 1)" .
eberlm@62085
   640
    moreover {
eberlm@62085
   641
      have "?u1 + ?u2 \<le> 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 \<open>n > 0\<close> have "4 * (1/2::real)^(48*n) \<le> (1/2)^(46*n)"
eberlm@62085
   644
        by (cases n) (simp_all add: field_simps power_add)
eberlm@62085
   645
      also from \<open>n > 0\<close> 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 \<le> (1/2) ^ (46 * n - 1)" by - simp
eberlm@62085
   648
    }
eberlm@62085
   649
    ultimately have "{-?l..?u1 + ?u2} \<subseteq> {-((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 "\<bar>pi - pi_approx2 n\<bar> \<le> ((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 \<le> 46*n - 1"
eberlm@62085
   658
  shows   "\<bar>pi - pi_approx2 n\<bar> \<le> 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 \<open>
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
\<close>
eberlm@62089
   667
lemma pi_approx_54_decimals:
eberlm@62085
   668
  "abs (pi - 3.141592653589793238462643383279502884197169399375105821 :: real) \<le> inverse (10^54)"
eberlm@62085
   669
  (is "abs (pi - ?pi') \<le> _")
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 \<Rightarrow> real \<Rightarrow> bool" where "eq = op ="
eberlm@62089
   683
eberlm@62085
   684
  have "abs (pi - pi_approx2 4) \<le> 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 (\<dots> - c) ?pi''"
eberlm@62085
   703
    by (unfold a_def b_def c_def, simp, unfold eq_def, rule refl)
wenzelm@62175
   704
  \<comment> \<open>This is incredibly slow because the numerators and denominators are huge.\<close>
eberlm@62085
   705
  finally show ?thesis by (rule approx_coarsen) simp
eberlm@62085
   706
qed
eberlm@62085
   707
eberlm@62085
   708
text \<open>A 128 bit approximation of pi:\<close>
eberlm@62085
   709
lemma pi_approx_128:
eberlm@62085
   710
  "abs (pi - 1069028584064966747859680373161870783301 / 2^128) \<le> inverse (2^128)"
eberlm@62085
   711
  by (rule approx_coarsen[OF pi_approx_54_decimals]) simp
eberlm@62085
   712
eberlm@62085
   713
text \<open>A 64 bit approximation of pi:\<close>
eberlm@62089
   714
lemma pi_approx_64:
eberlm@62085
   715
  "abs (pi - 57952155664616982739 / 2^64 :: real) \<le> inverse (2^64)"
eberlm@62085
   716
  by (rule approx_coarsen[OF pi_approx_54_decimals]) simp
eberlm@62089
   717
  
eberlm@62089
   718
text \<open>
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
\<close>
lp15@59871
   722
lp15@59871
   723
end