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