eberlm@62085: section \Numeric approximations to Constants\ lp15@59871: lp15@59871: theory Approximations eberlm@62085: imports "../Complex_Transcendental" "../Harmonic_Numbers" eberlm@62085: begin eberlm@62085: eberlm@62089: text \ eberlm@62089: In this theory, we will approximate some standard mathematical constants with high precision, eberlm@62089: using only Isabelle's simplifier. (no oracles, code generator, etc.) eberlm@62089: eberlm@62089: The constants we will look at are: $\pi$, $e$, $\ln 2$, and $\gamma$ (the Euler--Mascheroni eberlm@62089: constant). eberlm@62089: \ eberlm@62089: eberlm@62085: lemma eval_fact: eberlm@62085: "fact 0 = 1" eberlm@62085: "fact (Suc 0) = 1" eberlm@62085: "fact (numeral n) = numeral n * fact (pred_numeral n)" eberlm@62085: by (simp, simp, simp_all only: numeral_eq_Suc fact_Suc, eberlm@62085: simp only: numeral_eq_Suc [symmetric] of_nat_numeral) eberlm@62085: eberlm@62085: lemma setsum_poly_horner_expand: eberlm@62085: "(\k<(numeral n::nat). f k * x^k) = f 0 + (\kkk<(0::nat). f k * x^k) = 0" eberlm@62085: proof - eberlm@62085: { eberlm@62085: fix m :: nat eberlm@62085: have "(\kk=Suc 0..k=Suc 0..kkkkk 0" "x \ 0" "x < 1" eberlm@62085: shows "x ^ n < (1::'a::linordered_semidom)" eberlm@62085: proof - eberlm@62085: from assms consider "x > 0" | "x = 0" by force eberlm@62085: thus ?thesis eberlm@62085: proof cases eberlm@62085: assume "x > 0" eberlm@62085: with assms show ?thesis eberlm@62085: by (cases n) (simp, hypsubst, rule power_Suc_less_one) eberlm@62085: qed (insert assms, cases n, simp_all) eberlm@62085: qed eberlm@62085: eberlm@62085: lemma combine_bounds: eberlm@62085: "x \ {a1..b1} \ y \ {a2..b2} \ a3 = a1 + a2 \ b3 = b1 + b2 \ x + y \ {a3..(b3::real)}" eberlm@62085: "x \ {a1..b1} \ y \ {a2..b2} \ a3 = a1 - b2 \ b3 = b1 - a2 \ x - y \ {a3..(b3::real)}" eberlm@62085: "c \ 0 \ x \ {a..b} \ c * x \ {c*a..c*b}" eberlm@62085: by (auto simp: mult_left_mono) eberlm@62085: eberlm@62085: lemma approx_coarsen: eberlm@62085: "\x - a1\ \ eps1 \ \a1 - a2\ \ eps2 - eps1 \ \x - a2\ \ (eps2 :: real)" eberlm@62085: by simp eberlm@62085: eberlm@62085: eberlm@62089: eberlm@62089: subsection \Approximations of the exponential function\ eberlm@62089: eberlm@62089: lemma two_power_fact_le_fact: eberlm@62089: assumes "n \ 1" eberlm@62089: shows "2^k * fact n \ (fact (n + k) :: 'a :: {semiring_char_0,linordered_semidom})" eberlm@62089: proof (induction k) eberlm@62089: case (Suc k) eberlm@62089: have "2 ^ Suc k * fact n = 2 * (2 ^ k * fact n)" by (simp add: algebra_simps) eberlm@62089: also note Suc.IH eberlm@62089: also from assms have "of_nat 1 + of_nat 1 \ of_nat n + (of_nat (Suc k) :: 'a)" eberlm@62089: by (intro add_mono) (unfold of_nat_le_iff, simp_all) eberlm@62089: hence "2 * (fact (n + k) :: 'a) \ of_nat (n + Suc k) * fact (n + k)" eberlm@62089: by (intro mult_right_mono) (simp_all add: add_ac) eberlm@62089: also have "\ = fact (n + Suc k)" by simp eberlm@62089: finally show ?case by - (simp add: mult_left_mono) eberlm@62089: qed simp_all eberlm@62085: eberlm@62089: text \ eberlm@62089: We approximate the exponential function with inputs between $0$ and $2$ by its eberlm@62089: Taylor series expansion and bound the error term with $0$ from below and with a eberlm@62089: geometric series from above. eberlm@62089: \ eberlm@62089: lemma exp_approx: eberlm@62089: assumes "n > 0" "0 \ x" "x < 2" eberlm@62089: shows "exp (x::real) - (\k {0..(2 * x^n / (2 - x)) / fact n}" eberlm@62089: proof (unfold atLeastAtMost_iff, safe) wenzelm@63040: define approx where "approx = (\kk. x^k / fact k) sums exp x" eberlm@62089: using exp_converges[of x] by (simp add: field_simps) eberlm@62089: from sums_split_initial_segment[OF this, of n] eberlm@62089: have sums: "(\k. x^n * (x^k / fact (n+k))) sums (exp x - approx)" eberlm@62089: by (simp add: approx_def algebra_simps power_add) eberlm@62089: eberlm@62089: from assms show "(exp x - approx) \ 0" eberlm@62089: by (intro sums_le[OF _ sums_zero sums]) auto lp15@59871: eberlm@62089: have "\k. x^n * (x^k / fact (n+k)) \ (x^n / fact n) * (x / 2)^k" eberlm@62089: proof eberlm@62089: fix k :: nat eberlm@62089: have "x^n * (x^k / fact (n + k)) = x^(n+k) / fact (n + k)" by (simp add: power_add) eberlm@62089: also from assms have "\ \ x^(n+k) / (2^k * fact n)" eberlm@62089: by (intro divide_left_mono two_power_fact_le_fact zero_le_power) simp_all eberlm@62089: also have "\ = (x^n / fact n) * (x / 2) ^ k" eberlm@62089: by (simp add: field_simps power_add) eberlm@62089: finally show "x^n * (x^k / fact (n+k)) \ (x^n / fact n) * (x / 2)^k" . eberlm@62089: qed eberlm@62089: moreover note sums eberlm@62089: moreover { eberlm@62089: from assms have "(\k. (x^n / fact n) * (x / 2)^k) sums ((x^n / fact n) * (1 / (1 - x / 2)))" eberlm@62089: by (intro sums_mult geometric_sums) simp_all eberlm@62089: also from assms have "((x^n / fact n) * (1 / (1 - x / 2))) = (2 * x^n / (2 - x)) / fact n" eberlm@62089: by (auto simp: divide_simps) eberlm@62089: finally have "(\k. (x^n / fact n) * (x / 2)^k) sums \" . eberlm@62089: } eberlm@62089: ultimately show "(exp x - approx) \ (2 * x^n / (2 - x)) / fact n" eberlm@62089: by (rule sums_le) eberlm@62089: qed eberlm@62089: eberlm@62089: text \ eberlm@62089: The following variant gives a simpler error estimate for inputs between $0$ and $1$: eberlm@62089: \ eberlm@62089: lemma exp_approx': eberlm@62089: assumes "n > 0" "0 \ x" "x \ 1" eberlm@62089: shows "\exp (x::real) - (\k\n. x^k / fact k)\ \ x ^ n / fact n" eberlm@62085: proof - eberlm@62089: from assms have "x^n / (2 - x) \ x^n / 1" by (intro frac_le) simp_all eberlm@62089: hence "(2 * x^n / (2 - x)) / fact n \ 2 * x^n / fact n" eberlm@62089: using assms by (simp add: divide_simps) eberlm@62089: with exp_approx[of n x] assms eberlm@62089: have "exp (x::real) - (\k {0..2 * x^n / fact n}" by simp eberlm@62089: moreover have "(\k\n. x^k / fact k) = (\kexp (x::real) - (\k\n. x^k / fact k)\ \ x ^ n / fact n" eberlm@62089: unfolding atLeastAtMost_iff by linarith eberlm@62089: qed eberlm@62089: eberlm@62089: text \ eberlm@62089: By adding $x^n / n!$ to the approximation (i.e. taking one more term from the eberlm@62089: Taylor series), one can get the error bound down to $x^n / n!$. eberlm@62089: eberlm@62089: This means that the number of accurate binary digits produced by the approximation is eberlm@62089: asymptotically equal to $(n \log n - n) / \log 2$ by Stirling's formula. eberlm@62089: \ eberlm@62089: lemma exp_approx'': eberlm@62089: assumes "n > 0" "0 \ x" "x \ 1" eberlm@62089: shows "\exp (x::real) - (\k\n. x^k / fact k)\ \ 1 / fact n" eberlm@62089: proof - eberlm@62089: from assms have "\exp x - (\k\n. x ^ k / fact k)\ \ x ^ n / fact n" eberlm@62089: by (rule exp_approx') eberlm@62089: also from assms have "\ \ 1 / fact n" by (simp add: divide_simps power_le_one) eberlm@62085: finally show ?thesis . eberlm@62085: qed lp15@59871: eberlm@62089: eberlm@62089: text \ eberlm@62089: We now define an approximation function for Euler's constant $e$. eberlm@62089: \ eberlm@62089: eberlm@62089: definition euler_approx :: "nat \ real" where eberlm@62089: "euler_approx n = (\k\n. inverse (fact k))" eberlm@62089: eberlm@62089: definition euler_approx_aux :: "nat \ nat" where eberlm@62089: "euler_approx_aux n = (\k\n. \{k + 1..n})" eberlm@62089: eberlm@62089: lemma exp_1_approx: eberlm@62089: "n > 0 \ \exp (1::real) - euler_approx n\ \ 1 / fact n" eberlm@62089: using exp_approx''[of n 1] by (simp add: euler_approx_def divide_simps) eberlm@62089: eberlm@62089: text \ eberlm@62089: The following allows us to compute the numerator and the denominator of the result eberlm@62089: separately, which greatly reduces the amount of rational number arithmetic that we eberlm@62089: have to do. eberlm@62089: \ eberlm@62089: lemma euler_approx_altdef [code]: eberlm@62089: "euler_approx n = real (euler_approx_aux n) / real (fact n)" eberlm@62085: proof - eberlm@62089: have "real (\k\n. \{k+1..n}) = (\k\n. \i=k+1..n. real i)" by simp eberlm@62089: also have "\ / fact n = (\k\n. 1 / (fact n / (\i=k+1..n. real i)))" eberlm@62089: by (simp add: setsum_divide_distrib) eberlm@62089: also have "\ = (\k\n. 1 / fact k)" eberlm@62089: proof (intro setsum.cong refl) eberlm@62089: fix k assume k: "k \ {..n}" haftmann@63417: have "fact n = (\i=1..n. real i)" by (simp add: fact_setprod) eberlm@62089: also from k have "{1..n} = {1..k} \ {k+1..n}" by auto eberlm@62089: also have "setprod real \ / (\i=k+1..n. real i) = (\i=1..k. real i)" eberlm@62089: by (subst nonzero_divide_eq_eq, simp, subst setprod.union_disjoint [symmetric]) auto haftmann@63417: also have "\ = fact k" by (simp add: fact_setprod) eberlm@62089: finally show "1 / (fact n / setprod real {k + 1..n}) = 1 / fact k" by simp eberlm@62089: qed eberlm@62089: also have "\ = euler_approx n" by (simp add: euler_approx_def field_simps) eberlm@62089: finally show ?thesis by (simp add: euler_approx_aux_def) eberlm@62085: qed eberlm@62085: eberlm@62089: lemma euler_approx_aux_Suc: eberlm@62089: "euler_approx_aux (Suc m) = 1 + Suc m * euler_approx_aux m" eberlm@62089: unfolding euler_approx_aux_def eberlm@62089: by (subst setsum_right_distrib) (simp add: atLeastAtMostSuc_conv) eberlm@62089: eberlm@62089: lemma eval_euler_approx_aux: eberlm@62089: "euler_approx_aux 0 = 1" eberlm@62089: "euler_approx_aux 1 = 2" eberlm@62089: "euler_approx_aux (Suc 0) = 2" eberlm@62089: "euler_approx_aux (numeral n) = 1 + numeral n * euler_approx_aux (pred_numeral n)" (is "?th") eberlm@62089: proof - eberlm@62089: have A: "euler_approx_aux (Suc m) = 1 + Suc m * euler_approx_aux m" for m :: nat eberlm@62089: unfolding euler_approx_aux_def eberlm@62089: by (subst setsum_right_distrib) (simp add: atLeastAtMostSuc_conv) eberlm@62089: show ?th by (subst numeral_eq_Suc, subst A, subst numeral_eq_Suc [symmetric]) simp eberlm@62089: qed (simp_all add: euler_approx_aux_def) eberlm@62089: eberlm@62089: lemma euler_approx_aux_code [code]: eberlm@62089: "euler_approx_aux n = (if n = 0 then 1 else 1 + n * euler_approx_aux (n - 1))" eberlm@62089: by (cases n) (simp_all add: eval_euler_approx_aux euler_approx_aux_Suc) eberlm@62089: eberlm@62089: lemmas eval_euler_approx = euler_approx_altdef eval_euler_approx_aux eberlm@62089: eberlm@62089: eberlm@62089: text \Approximations of $e$ to 60 decimals / 128 and 64 bits:\ eberlm@62089: eberlm@62089: lemma euler_60_decimals: eberlm@62089: "\exp 1 - 2.718281828459045235360287471352662497757247093699959574966968\ eberlm@62089: \ inverse (10^60::real)" eberlm@62089: by (rule approx_coarsen, rule exp_1_approx[of 48]) eberlm@62089: (simp_all add: eval_euler_approx eval_fact) eberlm@62089: eberlm@62089: lemma euler_128: eberlm@62089: "\exp 1 - 924983374546220337150911035843336795079 / 2 ^ 128\ \ inverse (2 ^ 128 :: real)" eberlm@62089: by (rule approx_coarsen[OF euler_60_decimals]) simp_all eberlm@62089: eberlm@62089: lemma euler_64: eberlm@62089: "\exp 1 - 50143449209799256683 / 2 ^ 64\ \ inverse (2 ^ 64 :: real)" eberlm@62089: by (rule approx_coarsen[OF euler_128]) simp_all eberlm@62089: eberlm@62089: text \ eberlm@62089: An approximation of $e$ to 60 decimals. This is about as far as we can go with the eberlm@62089: simplifier with this kind of setup; the exported code of the code generator, on the other eberlm@62089: hand, can easily approximate $e$ to 1000 decimals and verify that approximation within eberlm@62089: fractions of a second. eberlm@62089: \ eberlm@62089: eberlm@62089: (* (Uncommented because we don't want to use the code generator; eberlm@62089: don't forget to import Code\_Target\_Numeral)) *) eberlm@62089: (* eberlm@62089: lemma "\exp 1 - 2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952605956307381323286279434907632338298807531952510190115738341879307021540891499348841675092447614606680822648001684774118537423454424371075390777449920695517027618386062613313845830007520449338265602976067371132007093287091274437470472306969772093101416928368190255151086574637721112523897844250569536967707854499699679468644549059879316368892300987931277361782154249992295763514822082698951936680331825288693984964651058209392398294887933203625094431173012381970684161403970198376793206832823764648042953118023287825098194558153017567173613320698112509961818815930416903515988885193458072738667385894228792284998920868058257492796104841984443634632449684875602336248270419786232090021609902353043699418491463140934317381436405462531520961836908887070167683964243781405927145635490613031072085103837505101157477041718986106873969655212671546889570350354021\ eberlm@62089: \ inverse (10^1000::real)" eberlm@62089: by (rule approx_coarsen, rule exp_1_approx[of 450], simp) eval eberlm@62089: *) eberlm@62089: eberlm@62089: eberlm@62089: subsection \Approximation of $\ln 2$\ eberlm@62089: eberlm@62089: text \ eberlm@62089: The following three auxiliary constants allow us to force the simplifier to eberlm@62089: evaluate intermediate results, simulating call-by-value. eberlm@62089: \ eberlm@62089: eberlm@62089: definition "ln_approx_aux3 x' e n y d \ eberlm@62089: \(2 * y) * (\k \ e - d" eberlm@62089: definition "ln_approx_aux2 x' e n y \ eberlm@62089: ln_approx_aux3 x' e n y (y^(2*n+1) / (1 - y^2) / real (2*n+1))" eberlm@62089: definition "ln_approx_aux1 x' e n x \ eberlm@62089: ln_approx_aux2 x' e n ((x - 1) / (x + 1))" eberlm@62089: eberlm@62089: lemma ln_approx_abs'': eberlm@62089: fixes x :: real and n :: nat eberlm@62089: defines "y \ (x-1)/(x+1)" eberlm@62089: defines "approx \ (\k y^(2*n+1) / (1 - y^2) / of_nat (2*n+1)" eberlm@62089: assumes x: "x > 1" eberlm@62089: assumes A: "ln_approx_aux1 x' e n x" eberlm@62089: shows "\ln x - x'\ \ e" eberlm@62089: proof (rule approx_coarsen[OF ln_approx_abs[OF x, of n]], goal_cases) eberlm@62089: case 1 eberlm@62089: from A have "\2 * y * (\k2 ^ k) + d - x'\ \ e - d" eberlm@62089: by (simp only: ln_approx_aux3_def ln_approx_aux2_def ln_approx_aux1_def eberlm@62089: y_def [symmetric] d_def [symmetric]) eberlm@62089: also have "2 * y * (\k2 ^ k) = eberlm@62089: (\k eberlm@62089: We unfold the above three constants successively and then compute the eberlm@62089: sum using a Horner scheme. eberlm@62089: \ eberlm@62089: lemma ln_2_40_decimals: eberlm@62089: "\ln 2 - 0.6931471805599453094172321214581765680755\ eberlm@62089: \ inverse (10^40 :: real)" eberlm@62089: apply (rule ln_approx_abs''[where n = 40], simp) eberlm@62089: apply (simp, simp add: ln_approx_aux1_def) eberlm@62089: apply (simp add: ln_approx_aux2_def power2_eq_square power_divide) eberlm@62089: apply (simp add: ln_approx_aux3_def power2_eq_square) eberlm@62089: apply (simp add: setsum_poly_horner_expand) eberlm@62089: done eberlm@62089: eberlm@62089: lemma ln_2_128: eberlm@62089: "\ln 2 - 235865763225513294137944142764154484399 / 2 ^ 128\ \ inverse (2 ^ 128 :: real)" eberlm@62089: by (rule approx_coarsen[OF ln_2_40_decimals]) simp_all eberlm@62089: eberlm@62089: lemma ln_2_64: eberlm@62089: "\ln 2 - 12786308645202655660 / 2 ^ 64\ \ inverse (2 ^ 64 :: real)" eberlm@62089: by (rule approx_coarsen[OF ln_2_128]) simp_all eberlm@62089: eberlm@62085: eberlm@62085: eberlm@62085: subsection \Approximation of the Euler--Mascheroni constant\ eberlm@62085: eberlm@62089: text \ eberlm@62089: Unfortunatly, the best approximation we have formalised for the Euler--Mascheroni eberlm@62089: constant converges only quadratically. This is too slow to compute more than a eberlm@62089: few decimals, but we can get almost 4 decimals / 14 binary digits this way, eberlm@62089: which is not too bad. eberlm@62089: \ eberlm@62089: lemma euler_mascheroni_approx: eberlm@62085: defines "approx \ 0.577257 :: real" and "e \ 0.000063 :: real" eberlm@62085: shows "abs (euler_mascheroni - approx :: real) < e" eberlm@62085: (is "abs (_ - ?approx) < ?e") lp15@59871: proof - wenzelm@63040: define l :: real wenzelm@63040: where "l = 47388813395531028639296492901910937/82101866951584879688289000000000000" wenzelm@63040: define u :: real wenzelm@63040: where "u = 142196984054132045946501548559032969 / 246305600854754639064867000000000000" eberlm@62085: have impI: "P \ Q" if Q for P Q using that by blast wenzelm@63040: have hsum_63: "harm 63 = (310559566510213034489743057 / 65681493561267903750631200 :: real)" eberlm@62085: by (simp add: harm_expand) eberlm@62089: from harm_Suc[of 63] have hsum_64: "harm 64 = eberlm@62089: 623171679694215690971693339 / (131362987122535807501262400::real)" eberlm@62085: by (subst (asm) hsum_63) simp eberlm@62085: have "ln (64::real) = real (6::nat) * ln 2" by (subst ln_realpow[symmetric]) simp_all eberlm@62089: hence "ln (real_of_nat (Suc 63)) \ {4.158883083293<..<4.158883083367}" using ln_2_64 nipkow@62390: by (simp add: abs_real_def split: if_split_asm) eberlm@62085: from euler_mascheroni_bounds'[OF _ this] eberlm@62089: have "(euler_mascheroni :: real) \ {l<.. \ {approx - e<..Approximation of pi\ eberlm@62085: eberlm@62085: eberlm@62085: subsubsection \Approximating the arctangent\ eberlm@62085: eberlm@62089: text\ eberlm@62089: The arctangent can be used to approximate pi. Fortunately, its Taylor series expansion eberlm@62089: converges exponentially for small values, so we can get $\Theta(n)$ digits of precision eberlm@62089: with $n$ summands of the expansion. eberlm@62089: \ eberlm@62089: eberlm@62085: definition arctan_approx where eberlm@62085: "arctan_approx n x = x * (\kx\ \ 1" eberlm@62085: shows "(\k. (-1)^k * (1 / real (k*2+1) * x ^ (k*2+1))) sums arctan x" eberlm@62085: using summable_arctan_series[OF assms] arctan_series[OF assms] by (simp add: sums_iff) eberlm@62085: eberlm@62085: lemma arctan_approx: eberlm@62085: assumes x: "0 \ x" "x < 1" and n: "even n" eberlm@62085: shows "arctan x - arctan_approx n x \ {0..x^(2*n+1) / (1-x^4)}" eberlm@62085: proof - wenzelm@63040: 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: from assms have "(\k. (-1) ^ k * (1 / real (k * 2 + 1) * x^(k*2+1))) sums arctan x" eberlm@62085: using arctan_series' by simp eberlm@62089: also have "(\k. (-1) ^ k * (1 / real (k * 2 + 1) * x^(k*2+1))) = eberlm@62085: (\k. x * ((- (x^2))^k / real (2*k+1)))" eberlm@62085: by (simp add: power2_eq_square power_mult power_mult_distrib mult_ac power_minus') eberlm@62085: finally have "(\k. x * ((- x\<^sup>2) ^ k / real (2 * k + 1))) sums arctan x" . eberlm@62085: from sums_split_initial_segment[OF this, of n] eberlm@62085: have "(\i. x * ((- x\<^sup>2) ^ (i + n) / real (2 * (i + n) + 1))) sums eberlm@62085: (arctan x - arctan_approx n x)" eberlm@62085: by (simp add: arctan_approx_def setsum_right_distrib) eberlm@62085: from sums_group[OF this, of 2] assms eberlm@62085: have sums: "(\k. x * (x\<^sup>2)^n * (x^4)^k * c k) sums (arctan x - arctan_approx n x)" eberlm@62085: by (simp add: algebra_simps power_add power_mult [symmetric] c_def) eberlm@62089: eberlm@62085: from assms have "0 \ arctan x - arctan_approx n x" eberlm@62085: by (intro sums_le[OF _ sums_zero sums] allI mult_nonneg_nonneg) eberlm@62085: (auto intro!: frac_le power_le_one simp: c_def) eberlm@62085: moreover { eberlm@62085: from assms have "c k \ 1 - 0" for k unfolding c_def eberlm@62085: by (intro diff_mono divide_nonneg_nonneg add_nonneg_nonneg) auto eberlm@62085: with assms have "x * x\<^sup>2 ^ n * (x ^ 4) ^ k * c k \ x * x\<^sup>2 ^ n * (x ^ 4) ^ k * 1" for k eberlm@62085: by (intro mult_left_mono mult_right_mono mult_nonneg_nonneg) simp_all eberlm@62085: with assms have "arctan x - arctan_approx n x \ x * (x\<^sup>2)^n * (1 / (1 - x^4))" eberlm@62085: by (intro sums_le[OF _ sums sums_mult[OF geometric_sums]] allI mult_left_mono) eberlm@62085: (auto simp: power_less_one) eberlm@62085: also have "x * (x^2)^n = x^(2*n+1)" by (simp add: power_mult power_add) eberlm@62085: finally have "arctan x - arctan_approx n x \ x^(2*n+1) / (1 - x^4)" by simp eberlm@62085: } eberlm@62085: ultimately show ?thesis by simp eberlm@62085: qed eberlm@62085: eberlm@62089: lemma arctan_approx_def': "arctan_approx n (1/x) = eberlm@62089: (\k2) ^ k)) / x" eberlm@62085: proof - eberlm@62085: have "(-1)^k / b = 1 / ((-1)^k * b)" for k :: nat and b :: real eberlm@62085: by (cases "even k") auto eberlm@62085: thus ?thesis by (simp add: arctan_approx_def field_simps power_minus') eberlm@62085: qed eberlm@62085: eberlm@62085: lemma expand_arctan_approx: eberlm@62089: "(\k<(numeral n::nat). inverse (f k) * inverse (x ^ k)) = eberlm@62085: inverse (f 0) + (\kkk<(0::nat). inverse (f k) * inverse (x^k)) = 0" eberlm@62085: proof - eberlm@62085: { eberlm@62085: fix m :: nat eberlm@62085: have "(\kk=Suc 0..k=Suc 0..kkkkkx*y::real\ < 1" eberlm@62085: shows "arctan x - arctan y = arctan ((x - y) / (1 + x * y))" eberlm@62085: proof - eberlm@62085: have "arctan x - arctan y = arctan x + arctan (-y)" by (simp add: arctan_minus) eberlm@62085: also from assms have "\ = arctan ((x - y) / (1 + x * y))" by (subst arctan_add_small) simp_all eberlm@62085: finally show ?thesis . lp15@59871: qed lp15@59871: eberlm@62085: eberlm@62085: subsubsection \Machin-like formulae for pi\ eberlm@62085: eberlm@62085: text \ eberlm@62085: We first define a small proof method that can prove Machin-like formulae for @{term "pi"} eberlm@62089: automatically. Unfortunately, this takes far too much time for larger formulae because eberlm@62085: the numbers involved become too large. eberlm@62085: \ eberlm@62085: eberlm@62085: definition "MACHIN_TAG a b \ a * (b::real)" eberlm@62085: eberlm@62085: lemma numeral_horner_MACHIN_TAG: eberlm@62085: "MACHIN_TAG Numeral1 x = x" eberlm@62089: "MACHIN_TAG (numeral (Num.Bit0 (Num.Bit0 n))) x = eberlm@62085: MACHIN_TAG 2 (MACHIN_TAG (numeral (Num.Bit0 n)) x)" eberlm@62089: "MACHIN_TAG (numeral (Num.Bit0 (Num.Bit1 n))) x = eberlm@62085: MACHIN_TAG 2 (MACHIN_TAG (numeral (Num.Bit1 n)) x)" eberlm@62089: "MACHIN_TAG (numeral (Num.Bit1 n)) x = eberlm@62085: MACHIN_TAG 2 (MACHIN_TAG (numeral n) x) + x" eberlm@62085: unfolding numeral_Bit0 numeral_Bit1 ring_distribs one_add_one[symmetric] MACHIN_TAG_def eberlm@62085: by (simp_all add: algebra_simps) eberlm@62085: eberlm@62085: lemma tag_machin: "a * arctan b = MACHIN_TAG a (arctan b)" by (simp add: MACHIN_TAG_def) eberlm@62085: eberlm@62085: lemma arctan_double': "\a::real\ < 1 \ MACHIN_TAG 2 (arctan a) = arctan (2 * a / (1 - a*a))" eberlm@62085: unfolding MACHIN_TAG_def by (simp add: arctan_double power2_eq_square) eberlm@62085: eberlm@62085: ML \ eberlm@62085: fun machin_term_conv ctxt ct = eberlm@62085: let eberlm@62085: val ctxt' = ctxt addsimps @{thms arctan_double' arctan_add_small} eberlm@62085: in eberlm@62085: case Thm.term_of ct of eberlm@62089: Const (@{const_name MACHIN_TAG}, _) $_$ eberlm@62089: (Const (@{const_name "Transcendental.arctan"}, _) $_) => eberlm@62085: Simplifier.rewrite ctxt' ct eberlm@62085: | eberlm@62089: Const (@{const_name MACHIN_TAG}, _)$ _ $eberlm@62089: (Const (@{const_name "Groups.plus"}, _)$ eberlm@62085: (Const (@{const_name "Transcendental.arctan"}, _) $_)$ eberlm@62089: (Const (@{const_name "Transcendental.arctan"}, _) $_)) => eberlm@62085: Simplifier.rewrite ctxt' ct eberlm@62085: | _ => raise CTERM ("machin_conv", [ct]) eberlm@62085: end eberlm@62085: eberlm@62089: fun machin_tac ctxt = eberlm@62085: let val conv = Conv.top_conv (Conv.try_conv o machin_term_conv) ctxt eberlm@62085: in eberlm@62085: SELECT_GOAL ( eberlm@62089: Local_Defs.unfold_tac ctxt eberlm@62085: @{thms tag_machin[THEN eq_reflection] numeral_horner_MACHIN_TAG[THEN eq_reflection]} eberlm@62085: THEN REPEAT (CHANGED (HEADGOAL (CONVERSION conv)))) eberlm@62085: THEN' Simplifier.simp_tac (ctxt addsimps @{thms arctan_add_small arctan_diff_small}) eberlm@62085: end eberlm@62085: \ eberlm@62085: eberlm@62085: method_setup machin = \Scan.succeed (SIMPLE_METHOD' o machin_tac)\ eberlm@62085: eberlm@62085: text \ eberlm@62089: We can now prove the standard'' Machin formula, which was already proven manually eberlm@62085: in Isabelle, automatically. eberlm@62085: }\ eberlm@62085: lemma "pi / 4 = (4::real) * arctan (1 / 5) - arctan (1 / 239)" eberlm@62085: by machin eberlm@62085: eberlm@62085: text \ eberlm@62085: We can also prove the following more complicated formula: eberlm@62085: \ eberlm@62085: lemma machin': "pi/4 = (12::real) * arctan (1/18) + 8 * arctan (1/57) - 5 * arctan (1/239)" eberlm@62085: by machin eberlm@62085: eberlm@62085: eberlm@62085: eberlm@62085: subsubsection \Simple approximation of pi\ eberlm@62085: eberlm@62085: text \ eberlm@62085: We can use the simple Machin formula and the Taylor series expansion of the arctangent eberlm@62089: to approximate pi. For a given even natural number$n$, we expand @{term "arctan (1/5)"} eberlm@62085: to$3n$summands and @{term "arctan (1/239)"} to$n$summands. This gives us at least eberlm@62085:$13n-2$bits of precision. eberlm@62085: \ eberlm@62085: eberlm@62085: definition "pi_approx n = 16 * arctan_approx (3*n) (1/5) - 4 * arctan_approx n (1/239)" eberlm@62085: eberlm@62085: lemma pi_approx: eberlm@62085: fixes n :: nat assumes n: "even n" and "n > 0" eberlm@62085: shows "\pi - pi_approx n\ \ inverse (2^(13*n - 2))" eberlm@62085: proof - eberlm@62085: from n have n': "even (3*n)" by simp wenzelm@62175: \ \We apply the Machin formula\ eberlm@62085: from machin have "pi = 16 * arctan (1/5) - 4 * arctan (1/239::real)" by simp wenzelm@62175: \ \Taylor series expansion of the arctangent\ eberlm@62085: also from arctan_approx[OF _ _ n', of "1/5"] arctan_approx[OF _ _ n, of "1/239"] eberlm@62085: have "\ - pi_approx n \ {-4*((1/239)^(2*n+1) / (1-(1/239)^4))..16*(1/5)^(6*n+1) / (1-(1/5)^4)}" eberlm@62085: by (simp add: pi_approx_def) wenzelm@62175: \ \Coarsening the bounds to make them a bit nicer\ eberlm@62085: also have "-4*((1/239::real)^(2*n+1) / (1-(1/239)^4)) = -((13651919 / 815702160) / 57121^n)" eberlm@62085: by (simp add: power_mult power2_eq_square) (simp add: field_simps) eberlm@62085: also have "16*(1/5)^(6*n+1) / (1-(1/5::real)^4) = (125/39) / 15625^n" eberlm@62085: by (simp add: power_mult power2_eq_square) (simp add: field_simps) eberlm@62089: also have "{-((13651919 / 815702160) / 57121^n) .. (125 / 39) / 15625^n} \ eberlm@62085: {- (4 / 2^(13*n)) .. 4 / (2^(13*n)::real)}" eberlm@62085: by (subst atLeastatMost_subset_iff, intro disjI2 conjI le_imp_neg_le) eberlm@62085: (rule frac_le; simp add: power_mult power_mono)+ eberlm@62085: finally have "abs (pi - pi_approx n) \ 4 / 2^(13*n)" by auto eberlm@62085: also from \n > 0\ have "4 / 2^(13*n) = 1 / (2^(13*n - 2) :: real)" eberlm@62085: by (cases n) (simp_all add: power_add) eberlm@62085: finally show ?thesis by (simp add: divide_inverse) eberlm@62085: qed eberlm@62085: eberlm@62085: lemma pi_approx': eberlm@62085: fixes n :: nat assumes n: "even n" and "n > 0" and "k \ 13*n - 2" eberlm@62085: shows "\pi - pi_approx n\ \ inverse (2^k)" eberlm@62085: using assms(3) by (intro order.trans[OF pi_approx[OF assms(1,2)]]) (simp_all add: field_simps) eberlm@62085: eberlm@62085: text \We can now approximate pi to 22 decimals within a fraction of a second.\ eberlm@62085: lemma pi_approx_75: "abs (pi - 3.1415926535897932384626 :: real) \ inverse (10^22)" eberlm@62085: proof - wenzelm@63040: define a :: real wenzelm@63040: where "a = 8295936325956147794769600190539918304 / 2626685325478320010006427764892578125" wenzelm@63040: define b :: real wenzelm@63040: where "b = 8428294561696506782041394632 / 503593538783547230635598424135" wenzelm@62175: \ \The introduction of this constant prevents the simplifier from applying solvers that eberlm@62085: we don't want. We want it to simply evaluate the terms to rational constants.}\ wenzelm@63040: define eq :: "real \ real \ bool" where "eq = op =" eberlm@62089: wenzelm@62175: \ \Splitting the computation into several steps has the advantage that simplification can eberlm@62085: be done in parallel\ eberlm@62085: have "abs (pi - pi_approx 6) \ inverse (2^76)" by (rule pi_approx') simp_all eberlm@62089: also have "pi_approx 6 = 16 * arctan_approx (3 * 6) (1 / 5) - 4 * arctan_approx 6 (1 / 239)" eberlm@62085: unfolding pi_approx_def by simp eberlm@62085: also have [unfolded eq_def]: "eq (16 * arctan_approx (3 * 6) (1 / 5)) a" eberlm@62085: by (simp add: arctan_approx_def' power2_eq_square, eberlm@62085: simp add: expand_arctan_approx, unfold a_def eq_def, rule refl) eberlm@62085: also have [unfolded eq_def]: "eq (4 * arctan_approx 6 (1 / 239::real)) b" eberlm@62085: by (simp add: arctan_approx_def' power2_eq_square, eberlm@62085: simp add: expand_arctan_approx, unfold b_def eq_def, rule refl) eberlm@62089: also have [unfolded eq_def]: eberlm@62085: "eq (a - b) (171331331860120333586637094112743033554946184594977368554649608 / eberlm@62085: 54536456744112171868276045488779391002026386559009552001953125)" eberlm@62085: by (unfold a_def b_def, simp, unfold eq_def, rule refl) eberlm@62085: finally show ?thesis by (rule approx_coarsen) simp eberlm@62085: qed eberlm@62085: eberlm@62085: text \ eberlm@62089: The previous estimate of pi in this file was based on approximating the root of the eberlm@62089:$\sin(\pi/6)$in the interval$[0;4]$using the Taylor series expansion of the sine to eberlm@62085: verify that it is between two given bounds. eberlm@62089: This was much slower and much less precise. We can easily recover this coarser estimate from eberlm@62085: the newer, precise estimate: eberlm@62085: \ eberlm@62085: lemma pi_approx_32: "\pi - 13493037705/4294967296 :: real\ \ inverse(2 ^ 32)" eberlm@62085: by (rule approx_coarsen[OF pi_approx_75]) simp eberlm@62085: eberlm@62085: eberlm@62085: subsection \A more complicated approximation of pi\ eberlm@62085: eberlm@62085: text \ eberlm@62089: There are more complicated Machin-like formulae that have more terms with larger eberlm@62085: denominators. Although they have more terms, each term requires fewer summands of the eberlm@62085: Taylor series for the same precision, since it is evaluated closer to$0\$. eberlm@62089: eberlm@62085: Using a good formula, one can therefore obtain the same precision with fewer operations. eberlm@62089: The big formulae used for computations of pi in practice are too complicated for us to eberlm@62085: prove here, but we can use the three-term Machin-like formula @{thm machin'}. eberlm@62085: \ eberlm@62085: eberlm@62089: definition "pi_approx2 n = 48 * arctan_approx (6*n) (1/18::real) + eberlm@62085: 32 * arctan_approx (4*n) (1/57) - 20 * arctan_approx (3*n) (1/239)" eberlm@62085: eberlm@62085: lemma pi_approx2: eberlm@62085: fixes n :: nat assumes n: "even n" and "n > 0" eberlm@62085: shows "\pi - pi_approx2 n\ \ inverse (2^(46*n - 1))" eberlm@62085: proof - eberlm@62085: from n have n': "even (6*n)" "even (4*n)" "even (3*n)" by simp_all eberlm@62089: from machin' have "pi = 48 * arctan (1/18) + 32 * arctan (1/57) - 20 * arctan (1/239::real)" eberlm@62085: by simp eberlm@62085: hence "pi - pi_approx2 n = 48 * (arctan (1/18) - arctan_approx (6*n) (1/18)) + eberlm@62085: 32 * (arctan (1/57) - arctan_approx (4*n) (1/57)) - eberlm@62085: 20 * (arctan (1/239) - arctan_approx (3*n) (1/239))" eberlm@62085: by (simp add: pi_approx2_def) eberlm@62085: also have "\ \ {-((20/239/(1-(1/239)^4)) * (1/239)^(6*n)).. eberlm@62089: (48/18 / (1-(1/18)^4))*(1/18)^(12*n) + (32/57/(1-(1/57)^4)) * (1/57)^(8*n)}" eberlm@62085: (is "_ \ {-?l..?u1 + ?u2}") eberlm@62085: apply ((rule combine_bounds(1,2))+; (rule combine_bounds(3); (rule arctan_approx)?)?) eberlm@62085: apply (simp_all add: n) eberlm@62085: apply (simp_all add: divide_simps)? eberlm@62085: done eberlm@62085: also { eberlm@62085: have "?l \ (1/8) * (1/2)^(46*n)" eberlm@62085: unfolding power_mult by (intro mult_mono power_mono) (simp_all add: divide_simps) eberlm@62085: also have "\ \ (1/2) ^ (46 * n - 1)" eberlm@62085: by (cases n; simp_all add: power_add divide_simps) eberlm@62085: finally have "?l \ (1/2) ^ (46 * n - 1)" . eberlm@62085: moreover { eberlm@62085: have "?u1 + ?u2 \ 4 * (1/2)^(48*n) + 1 * (1/2)^(46*n)" eberlm@62085: unfolding power_mult by (intro add_mono mult_mono power_mono) (simp_all add: divide_simps) eberlm@62089: also from \n > 0\ have "4 * (1/2::real)^(48*n) \ (1/2)^(46*n)" eberlm@62085: by (cases n) (simp_all add: field_simps power_add) eberlm@62085: also from \n > 0\ have "(1/2::real) ^ (46 * n) + 1 * (1 / 2) ^ (46 * n) = (1/2) ^ (46 * n - 1)" eberlm@62085: by (cases n; simp_all add: power_add power_divide) eberlm@62085: finally have "?u1 + ?u2 \ (1/2) ^ (46 * n - 1)" by - simp eberlm@62085: } eberlm@62085: ultimately have "{-?l..?u1 + ?u2} \ {-((1/2)^(46*n-1))..(1/2)^(46*n-1)}" eberlm@62085: by (subst atLeastatMost_subset_iff) simp_all eberlm@62085: } eberlm@62085: finally have "\pi - pi_approx2 n\ \ ((1/2) ^ (46 * n - 1))" by auto eberlm@62085: thus ?thesis by (simp add: divide_simps) eberlm@62085: qed eberlm@62085: eberlm@62085: lemma pi_approx2': eberlm@62085: fixes n :: nat assumes n: "even n" and "n > 0" and "k \ 46*n - 1" eberlm@62085: shows "\pi - pi_approx2 n\ \ inverse (2^k)" eberlm@62085: using assms(3) by (intro order.trans[OF pi_approx2[OF assms(1,2)]]) (simp_all add: field_simps) eberlm@62085: eberlm@62085: text \ eberlm@62089: We can now approximate pi to 54 decimals using this formula. The computations are much eberlm@62089: slower now; this is mostly because we use arbitrary-precision rational numbers, whose eberlm@62089: numerators and demoninators get very large. Using dyadic floating point numbers would be eberlm@62085: much more economical. eberlm@62085: \ eberlm@62089: lemma pi_approx_54_decimals: eberlm@62085: "abs (pi - 3.141592653589793238462643383279502884197169399375105821 :: real) \ inverse (10^54)" eberlm@62085: (is "abs (pi - ?pi') \ _") eberlm@62085: proof - wenzelm@63040: define a :: real wenzelm@63040: where "a = 2829469759662002867886529831139137601191652261996513014734415222704732791803 / wenzelm@63040: 1062141879292765061960538947347721564047051545995266466660439319087625011200" wenzelm@63040: define b :: real wenzelm@63040: where "b = 13355545553549848714922837267299490903143206628621657811747118592 / wenzelm@63040: 23792006023392488526789546722992491355941103837356113731091180925" wenzelm@63040: define c :: real wenzelm@63040: where "c = 28274063397213534906669125255762067746830085389618481175335056 / wenzelm@63040: 337877029279505250241149903214554249587517250716358486542628059" eberlm@62085: let ?pi'' = "3882327391761098513316067116522233897127356523627918964967729040413954225768920394233198626889767468122598417405434625348404038165437924058179155035564590497837027530349 / eberlm@62089: 1235783190199688165469648572769847552336447197542738425378629633275352407743112409829873464564018488572820294102599160968781449606552922108667790799771278860366957772800" wenzelm@63040: define eq :: "real \ real \ bool" where "eq = op =" eberlm@62089: eberlm@62085: have "abs (pi - pi_approx2 4) \ inverse (2^183)" by (rule pi_approx2') simp_all eberlm@62085: also have "pi_approx2 4 = 48 * arctan_approx 24 (1 / 18) + eberlm@62085: 32 * arctan_approx 16 (1 / 57) - eberlm@62089: 20 * arctan_approx 12 (1 / 239)" eberlm@62085: unfolding pi_approx2_def by simp eberlm@62085: also have [unfolded eq_def]: "eq (48 * arctan_approx 24 (1 / 18)) a" eberlm@62085: by (simp add: arctan_approx_def' power2_eq_square, eberlm@62085: simp add: expand_arctan_approx, unfold a_def eq_def, rule refl) eberlm@62085: also have [unfolded eq_def]: "eq (32 * arctan_approx 16 (1 / 57::real)) b" eberlm@62085: by (simp add: arctan_approx_def' power2_eq_square, eberlm@62085: simp add: expand_arctan_approx, unfold b_def eq_def, rule refl) eberlm@62085: also have [unfolded eq_def]: "eq (20 * arctan_approx 12 (1 / 239::real)) c" eberlm@62085: by (simp add: arctan_approx_def' power2_eq_square, eberlm@62085: simp add: expand_arctan_approx, unfold c_def eq_def, rule refl) eberlm@62085: also have [unfolded eq_def]: eberlm@62085: "eq (a + b) (34326487387865555303797183505809267914709125998469664969258315922216638779011304447624792548723974104030355722677 / eberlm@62089: 10642967245546718617684989689985787964158885991018703366677373121531695267093031090059801733340658960857196134400)" eberlm@62085: by (unfold a_def b_def c_def, simp, unfold eq_def, rule refl) eberlm@62085: also have [unfolded eq_def]: "eq (\ - c) ?pi''" eberlm@62085: by (unfold a_def b_def c_def, simp, unfold eq_def, rule refl) wenzelm@62175: \ \This is incredibly slow because the numerators and denominators are huge.\ eberlm@62085: finally show ?thesis by (rule approx_coarsen) simp eberlm@62085: qed eberlm@62085: eberlm@62085: text \A 128 bit approximation of pi:\ eberlm@62085: lemma pi_approx_128: eberlm@62085: "abs (pi - 1069028584064966747859680373161870783301 / 2^128) \ inverse (2^128)" eberlm@62085: by (rule approx_coarsen[OF pi_approx_54_decimals]) simp eberlm@62085: eberlm@62085: text \A 64 bit approximation of pi:\ eberlm@62089: lemma pi_approx_64: eberlm@62085: "abs (pi - 57952155664616982739 / 2^64 :: real) \ inverse (2^64)" eberlm@62085: by (rule approx_coarsen[OF pi_approx_54_decimals]) simp eberlm@62089: eberlm@62089: text \ eberlm@62089: Again, going much farther with the simplifier takes a long time, but the code generator eberlm@62089: can handle even two thousand decimal digits in under 20 seconds. eberlm@62089: \ lp15@59871: lp15@59871: end