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