52 lemma approx_coarsen: |
60 lemma approx_coarsen: |
53 "\<bar>x - a1\<bar> \<le> eps1 \<Longrightarrow> \<bar>a1 - a2\<bar> \<le> eps2 - eps1 \<Longrightarrow> \<bar>x - a2\<bar> \<le> (eps2 :: real)" |
61 "\<bar>x - a1\<bar> \<le> eps1 \<Longrightarrow> \<bar>a1 - a2\<bar> \<le> eps2 - eps1 \<Longrightarrow> \<bar>x - a2\<bar> \<le> (eps2 :: real)" |
54 by simp |
62 by simp |
55 |
63 |
56 |
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 def approx \<equiv> "(\<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 (rule fact_altdef) |
|
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 only: fact_altdef) |
|
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 |
57 subsection \<open>Approximation of $\ln 2$\<close> |
255 subsection \<open>Approximation of $\ln 2$\<close> |
58 |
256 |
59 context |
257 text \<open> |
60 begin |
258 The following three auxiliary constants allow us to force the simplifier to |
61 |
259 evaluate intermediate results, simulating call-by-value. |
62 qualified lemma ln_approx_abs': |
260 \<close> |
63 assumes "x > (1::real)" |
261 |
64 assumes "(x-1)/(x+1) = y" |
262 definition "ln_approx_aux3 x' e n y d \<longleftrightarrow> |
65 assumes "y^2 = ysqr" |
263 \<bar>(2 * y) * (\<Sum>k<n. inverse (real (2*k+1)) * (y^2)^k) + d - x'\<bar> \<le> e - d" |
66 assumes "(\<Sum>k<n. inverse (of_nat (2*k+1)) * ysqr^k) = approx" |
264 definition "ln_approx_aux2 x' e n y \<longleftrightarrow> |
67 assumes "y*ysqr^n / (1 - ysqr) / of_nat (2*n+1) = d" |
265 ln_approx_aux3 x' e n y (y^(2*n+1) / (1 - y^2) / real (2*n+1))" |
68 assumes "d \<le> e" |
266 definition "ln_approx_aux1 x' e n x \<longleftrightarrow> |
69 shows "abs (ln x - (2*y*approx + d)) \<le> e" |
267 ln_approx_aux2 x' e n ((x - 1) / (x + 1))" |
70 proof - |
268 |
71 note ln_approx_abs[OF assms(1), of n] |
269 lemma ln_approx_abs'': |
72 also note assms(2) |
270 fixes x :: real and n :: nat |
73 also have "y^(2*n+1) = y*ysqr^n" by (simp add: assms(3)[symmetric] power_mult) |
271 defines "y \<equiv> (x-1)/(x+1)" |
74 also note assms(3) |
272 defines "approx \<equiv> (\<Sum>k<n. 2*y^(2*k+1) / of_nat (2*k+1))" |
75 also note assms(5) |
273 defines "d \<equiv> y^(2*n+1) / (1 - y^2) / of_nat (2*n+1)" |
76 also note assms(5) |
274 assumes x: "x > 1" |
77 also note assms(6) |
275 assumes A: "ln_approx_aux1 x' e n x" |
78 also have "(\<Sum>k<n. 2*y^(2*k+1) / real_of_nat (2 * k + 1)) = (2*y) * approx" |
276 shows "\<bar>ln x - x'\<bar> \<le> e" |
79 apply (subst assms(4)[symmetric], subst setsum_right_distrib) |
277 proof (rule approx_coarsen[OF ln_approx_abs[OF x, of n]], goal_cases) |
80 apply (simp add: assms(3)[symmetric] power_mult) |
278 case 1 |
81 apply (simp add: mult_ac divide_simps)? |
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" |
82 done |
280 by (simp only: ln_approx_aux3_def ln_approx_aux2_def ln_approx_aux1_def |
83 finally show ?thesis . |
281 y_def [symmetric] d_def [symmetric]) |
84 qed |
282 also have "2 * y * (\<Sum>k<n. inverse (real (2 * k + 1)) * y\<^sup>2 ^ k) = |
85 |
283 (\<Sum>k<n. 2 * y^(2*k+1) / (real (2 * k + 1)))" |
86 lemma ln_2_approx: "\<bar>ln 2 - 0.69314718055\<bar> < inverse (2 ^ 36 :: real)" (is ?thesis1) |
284 by (subst setsum_right_distrib, simp, subst power_mult) |
87 and ln_2_bounds: "ln (2::real) \<in> {0.693147180549..0.693147180561}" (is ?thesis2) |
285 (simp_all add: divide_simps mult_ac power_mult) |
88 proof - |
286 finally show ?case by (simp only: d_def y_def approx_def) |
89 def approx \<equiv> "0.69314718055 :: real" and approx' \<equiv> "4465284211343447 / 6442043387911560 :: real" |
287 qed |
90 def d \<equiv> "inverse (195259926456::real)" |
288 |
91 have "dist (ln 2) approx \<le> dist (ln 2) approx' + dist approx' approx" by (rule dist_triangle) |
289 text \<open> |
92 also have "\<bar>ln (2::real) - (2 * (1/3) * (651187280816108 / 626309773824735) + |
290 We unfold the above three constants successively and then compute the |
93 inverse 195259926456)\<bar> \<le> inverse 195259926456" |
291 sum using a Horner scheme. |
94 proof (rule ln_approx_abs'[where n = 10]) |
292 \<close> |
95 show "(1/3::real)^2 = 1/9" by (simp add: power2_eq_square) |
293 lemma ln_2_40_decimals: |
96 qed (simp_all add: eval_nat_numeral) |
294 "\<bar>ln 2 - 0.6931471805599453094172321214581765680755\<bar> |
97 hence A: "dist (ln 2) approx' \<le> d" by (simp add: dist_real_def approx'_def d_def) |
295 \<le> inverse (10^40 :: real)" |
98 hence "dist (ln 2) approx' + dist approx' approx \<le> \<dots> + dist approx' approx" |
296 apply (rule ln_approx_abs''[where n = 40], simp) |
99 by (rule add_right_mono) |
297 apply (simp, simp add: ln_approx_aux1_def) |
100 also have "\<dots> < inverse (2 ^ 36)" by (simp add: dist_real_def approx'_def approx_def d_def) |
298 apply (simp add: ln_approx_aux2_def power2_eq_square power_divide) |
101 finally show ?thesis1 unfolding dist_real_def approx_def . |
299 apply (simp add: ln_approx_aux3_def power2_eq_square) |
102 |
300 apply (simp add: setsum_poly_horner_expand) |
103 from A have "ln 2 \<in> {approx' - d..approx' + d}" |
301 done |
104 by (simp add: dist_real_def abs_real_def split: split_if_asm) |
302 |
105 also have "\<dots> \<subseteq> {0.693147180549..0.693147180561}" |
303 lemma ln_2_128: |
106 by (subst atLeastatMost_subset_iff, rule disjI2) (simp add: approx'_def d_def) |
304 "\<bar>ln 2 - 235865763225513294137944142764154484399 / 2 ^ 128\<bar> \<le> inverse (2 ^ 128 :: real)" |
107 finally show ?thesis2 . |
305 by (rule approx_coarsen[OF ln_2_40_decimals]) simp_all |
108 qed |
306 |
109 |
307 lemma ln_2_64: |
110 end |
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 |
111 |
311 |
112 |
312 |
113 subsection \<open>Approximation of the Euler--Mascheroni constant\<close> |
313 subsection \<open>Approximation of the Euler--Mascheroni constant\<close> |
114 |
314 |
115 lemma euler_mascheroni_approx: |
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: |
116 defines "approx \<equiv> 0.577257 :: real" and "e \<equiv> 0.000063 :: real" |
322 defines "approx \<equiv> 0.577257 :: real" and "e \<equiv> 0.000063 :: real" |
117 shows "abs (euler_mascheroni - approx :: real) < e" |
323 shows "abs (euler_mascheroni - approx :: real) < e" |
118 (is "abs (_ - ?approx) < ?e") |
324 (is "abs (_ - ?approx) < ?e") |
119 proof - |
325 proof - |
120 def l \<equiv> "47388813395531028639296492901910937/82101866951584879688289000000000000 :: real" |
326 def l \<equiv> "47388813395531028639296492901910937/82101866951584879688289000000000000 :: real" |
121 def u \<equiv> "142196984054132045946501548559032969 / 246305600854754639064867000000000000 :: real" |
327 def u \<equiv> "142196984054132045946501548559032969 / 246305600854754639064867000000000000 :: real" |
122 have impI: "P \<longrightarrow> Q" if Q for P Q using that by blast |
328 have impI: "P \<longrightarrow> Q" if Q for P Q using that by blast |
123 have hsum_63: "harm 63 = (310559566510213034489743057 / 65681493561267903750631200 ::real)" |
329 have hsum_63: "harm 63 = (310559566510213034489743057 / 65681493561267903750631200 ::real)" |
124 by (simp add: harm_expand) |
330 by (simp add: harm_expand) |
125 from harm_Suc[of 63] have hsum_64: "harm 64 = |
331 from harm_Suc[of 63] have hsum_64: "harm 64 = |
126 623171679694215690971693339 / (131362987122535807501262400::real)" |
332 623171679694215690971693339 / (131362987122535807501262400::real)" |
127 by (subst (asm) hsum_63) simp |
333 by (subst (asm) hsum_63) simp |
128 have "ln (64::real) = real (6::nat) * ln 2" by (subst ln_realpow[symmetric]) simp_all |
334 have "ln (64::real) = real (6::nat) * ln 2" by (subst ln_realpow[symmetric]) simp_all |
129 hence "ln (real_of_nat (Suc 63)) \<in> {4.158883083293<..<4.158883083367}" using ln_2_bounds by simp |
335 hence "ln (real_of_nat (Suc 63)) \<in> {4.158883083293<..<4.158883083367}" using ln_2_64 |
|
336 by (simp add: abs_real_def split: split_if_asm) |
130 from euler_mascheroni_bounds'[OF _ this] |
337 from euler_mascheroni_bounds'[OF _ this] |
131 have "(euler_mascheroni :: real) \<in> {l<..<u}" |
338 have "(euler_mascheroni :: real) \<in> {l<..<u}" |
132 by (simp add: hsum_63 del: greaterThanLessThan_iff) (simp only: l_def u_def) |
339 by (simp add: hsum_63 del: greaterThanLessThan_iff) (simp only: l_def u_def) |
133 also have "\<dots> \<subseteq> {approx - e<..<approx + e}" |
340 also have "\<dots> \<subseteq> {approx - e<..<approx + e}" |
134 by (subst greaterThanLessThan_subseteq_greaterThanLessThan, rule impI) |
341 by (subst greaterThanLessThan_subseteq_greaterThanLessThan, rule impI) |
135 (simp add: approx_def e_def u_def l_def) |
342 (simp add: approx_def e_def u_def l_def) |
136 finally show ?thesis by (simp add: abs_real_def) |
343 finally show ?thesis by (simp add: abs_real_def) |
137 qed |
344 qed |
138 |
345 |
139 |
346 |
140 subsection \<open>Approximation to pi\<close> |
347 |
|
348 subsection \<open>Approximation of pi\<close> |
141 |
349 |
142 |
350 |
143 subsubsection \<open>Approximating the arctangent\<close> |
351 subsubsection \<open>Approximating the arctangent\<close> |
|
352 |
|
353 text\<open> |
|
354 The arctangent can be used to approximate pi. Fortunately, its Taylor series expansion |
|
355 converges exponentially for small values, so we can get $\Theta(n)$ digits of precision |
|
356 with $n$ summands of the expansion. |
|
357 \<close> |
144 |
358 |
145 definition arctan_approx where |
359 definition arctan_approx where |
146 "arctan_approx n x = x * (\<Sum>k<n. (-(x^2))^k / real (2*k+1))" |
360 "arctan_approx n x = x * (\<Sum>k<n. (-(x^2))^k / real (2*k+1))" |
147 |
361 |
148 lemma arctan_series': |
362 lemma arctan_series': |
345 text \<open>We can now approximate pi to 22 decimals within a fraction of a second.\<close> |
559 text \<open>We can now approximate pi to 22 decimals within a fraction of a second.\<close> |
346 lemma pi_approx_75: "abs (pi - 3.1415926535897932384626 :: real) \<le> inverse (10^22)" |
560 lemma pi_approx_75: "abs (pi - 3.1415926535897932384626 :: real) \<le> inverse (10^22)" |
347 proof - |
561 proof - |
348 def a \<equiv> "8295936325956147794769600190539918304 / 2626685325478320010006427764892578125 :: real" |
562 def a \<equiv> "8295936325956147794769600190539918304 / 2626685325478320010006427764892578125 :: real" |
349 def b \<equiv> "8428294561696506782041394632 / 503593538783547230635598424135 :: real" |
563 def b \<equiv> "8428294561696506782041394632 / 503593538783547230635598424135 :: real" |
350 -- \<open>The introduction of this constant prevents the simplifier from applying solvers that |
564 -- \<open>The introduction of this constant prevents the simplifier from applying solvers that |
351 we don't want. We want it to simply evaluate the terms to rational constants.}\<close> |
565 we don't want. We want it to simply evaluate the terms to rational constants.}\<close> |
352 def eq \<equiv> "op = :: real \<Rightarrow> real \<Rightarrow> bool" |
566 def eq \<equiv> "op = :: real \<Rightarrow> real \<Rightarrow> bool" |
353 |
567 |
354 -- \<open>Splitting the computation into several steps has the advantage that simplification can |
568 -- \<open>Splitting the computation into several steps has the advantage that simplification can |
355 be done in parallel\<close> |
569 be done in parallel\<close> |
356 have "abs (pi - pi_approx 6) \<le> inverse (2^76)" by (rule pi_approx') simp_all |
570 have "abs (pi - pi_approx 6) \<le> inverse (2^76)" by (rule pi_approx') simp_all |
357 also have "pi_approx 6 = 16 * arctan_approx (3 * 6) (1 / 5) - 4 * arctan_approx 6 (1 / 239)" |
571 also have "pi_approx 6 = 16 * arctan_approx (3 * 6) (1 / 5) - 4 * arctan_approx 6 (1 / 239)" |
358 unfolding pi_approx_def by simp |
572 unfolding pi_approx_def by simp |
359 also have [unfolded eq_def]: "eq (16 * arctan_approx (3 * 6) (1 / 5)) a" |
573 also have [unfolded eq_def]: "eq (16 * arctan_approx (3 * 6) (1 / 5)) a" |
360 by (simp add: arctan_approx_def' power2_eq_square, |
574 by (simp add: arctan_approx_def' power2_eq_square, |
361 simp add: expand_arctan_approx, unfold a_def eq_def, rule refl) |
575 simp add: expand_arctan_approx, unfold a_def eq_def, rule refl) |
362 also have [unfolded eq_def]: "eq (4 * arctan_approx 6 (1 / 239::real)) b" |
576 also have [unfolded eq_def]: "eq (4 * arctan_approx 6 (1 / 239::real)) b" |
363 by (simp add: arctan_approx_def' power2_eq_square, |
577 by (simp add: arctan_approx_def' power2_eq_square, |
364 simp add: expand_arctan_approx, unfold b_def eq_def, rule refl) |
578 simp add: expand_arctan_approx, unfold b_def eq_def, rule refl) |
365 also have [unfolded eq_def]: |
579 also have [unfolded eq_def]: |
366 "eq (a - b) (171331331860120333586637094112743033554946184594977368554649608 / |
580 "eq (a - b) (171331331860120333586637094112743033554946184594977368554649608 / |
367 54536456744112171868276045488779391002026386559009552001953125)" |
581 54536456744112171868276045488779391002026386559009552001953125)" |
368 by (unfold a_def b_def, simp, unfold eq_def, rule refl) |
582 by (unfold a_def b_def, simp, unfold eq_def, rule refl) |
369 finally show ?thesis by (rule approx_coarsen) simp |
583 finally show ?thesis by (rule approx_coarsen) simp |
370 qed |
584 qed |
371 |
585 |
372 text \<open> |
586 text \<open> |
373 The previous estimate of pi in this file was based on approximating the root of the |
587 The previous estimate of pi in this file was based on approximating the root of the |
374 $\sin(\pi/6)$ in the interval $[0;4]$ using the Taylor series expansion of the sine to |
588 $\sin(\pi/6)$ in the interval $[0;4]$ using the Taylor series expansion of the sine to |
375 verify that it is between two given bounds. |
589 verify that it is between two given bounds. |
376 This was much slower and much less precise. We can easily recover this coarser estimate from |
590 This was much slower and much less precise. We can easily recover this coarser estimate from |
377 the newer, precise estimate: |
591 the newer, precise estimate: |
378 \<close> |
592 \<close> |
379 lemma pi_approx_32: "\<bar>pi - 13493037705/4294967296 :: real\<bar> \<le> inverse(2 ^ 32)" |
593 lemma pi_approx_32: "\<bar>pi - 13493037705/4294967296 :: real\<bar> \<le> inverse(2 ^ 32)" |
380 by (rule approx_coarsen[OF pi_approx_75]) simp |
594 by (rule approx_coarsen[OF pi_approx_75]) simp |
381 |
595 |
382 |
596 |
383 subsection \<open>A more complicated approximation of pi\<close> |
597 subsection \<open>A more complicated approximation of pi\<close> |
384 |
598 |
385 text \<open> |
599 text \<open> |
386 There are more complicated Machin-like formulae that have more terms with larger |
600 There are more complicated Machin-like formulae that have more terms with larger |
387 denominators. Although they have more terms, each term requires fewer summands of the |
601 denominators. Although they have more terms, each term requires fewer summands of the |
388 Taylor series for the same precision, since it is evaluated closer to $0$. |
602 Taylor series for the same precision, since it is evaluated closer to $0$. |
389 |
603 |
390 Using a good formula, one can therefore obtain the same precision with fewer operations. |
604 Using a good formula, one can therefore obtain the same precision with fewer operations. |
391 The big formulae used for computations of pi in practice are too complicated for us to |
605 The big formulae used for computations of pi in practice are too complicated for us to |
392 prove here, but we can use the three-term Machin-like formula @{thm machin'}. |
606 prove here, but we can use the three-term Machin-like formula @{thm machin'}. |
393 \<close> |
607 \<close> |
394 |
608 |
395 definition "pi_approx2 n = 48 * arctan_approx (6*n) (1/18::real) + |
609 definition "pi_approx2 n = 48 * arctan_approx (6*n) (1/18::real) + |
396 32 * arctan_approx (4*n) (1/57) - 20 * arctan_approx (3*n) (1/239)" |
610 32 * arctan_approx (4*n) (1/57) - 20 * arctan_approx (3*n) (1/239)" |
397 |
611 |
398 lemma pi_approx2: |
612 lemma pi_approx2: |
399 fixes n :: nat assumes n: "even n" and "n > 0" |
613 fixes n :: nat assumes n: "even n" and "n > 0" |
400 shows "\<bar>pi - pi_approx2 n\<bar> \<le> inverse (2^(46*n - 1))" |
614 shows "\<bar>pi - pi_approx2 n\<bar> \<le> inverse (2^(46*n - 1))" |
401 proof - |
615 proof - |
402 from n have n': "even (6*n)" "even (4*n)" "even (3*n)" by simp_all |
616 from n have n': "even (6*n)" "even (4*n)" "even (3*n)" by simp_all |
403 from machin' have "pi = 48 * arctan (1/18) + 32 * arctan (1/57) - 20 * arctan (1/239::real)" |
617 from machin' have "pi = 48 * arctan (1/18) + 32 * arctan (1/57) - 20 * arctan (1/239::real)" |
404 by simp |
618 by simp |
405 hence "pi - pi_approx2 n = 48 * (arctan (1/18) - arctan_approx (6*n) (1/18)) + |
619 hence "pi - pi_approx2 n = 48 * (arctan (1/18) - arctan_approx (6*n) (1/18)) + |
406 32 * (arctan (1/57) - arctan_approx (4*n) (1/57)) - |
620 32 * (arctan (1/57) - arctan_approx (4*n) (1/57)) - |
407 20 * (arctan (1/239) - arctan_approx (3*n) (1/239))" |
621 20 * (arctan (1/239) - arctan_approx (3*n) (1/239))" |
408 by (simp add: pi_approx2_def) |
622 by (simp add: pi_approx2_def) |
409 also have "\<dots> \<in> {-((20/239/(1-(1/239)^4)) * (1/239)^(6*n)).. |
623 also have "\<dots> \<in> {-((20/239/(1-(1/239)^4)) * (1/239)^(6*n)).. |
410 (48/18 / (1-(1/18)^4))*(1/18)^(12*n) + (32/57/(1-(1/57)^4)) * (1/57)^(8*n)}" |
624 (48/18 / (1-(1/18)^4))*(1/18)^(12*n) + (32/57/(1-(1/57)^4)) * (1/57)^(8*n)}" |
411 (is "_ \<in> {-?l..?u1 + ?u2}") |
625 (is "_ \<in> {-?l..?u1 + ?u2}") |
412 apply ((rule combine_bounds(1,2))+; (rule combine_bounds(3); (rule arctan_approx)?)?) |
626 apply ((rule combine_bounds(1,2))+; (rule combine_bounds(3); (rule arctan_approx)?)?) |
413 apply (simp_all add: n) |
627 apply (simp_all add: n) |
414 apply (simp_all add: divide_simps)? |
628 apply (simp_all add: divide_simps)? |
415 done |
629 done |
439 fixes n :: nat assumes n: "even n" and "n > 0" and "k \<le> 46*n - 1" |
653 fixes n :: nat assumes n: "even n" and "n > 0" and "k \<le> 46*n - 1" |
440 shows "\<bar>pi - pi_approx2 n\<bar> \<le> inverse (2^k)" |
654 shows "\<bar>pi - pi_approx2 n\<bar> \<le> inverse (2^k)" |
441 using assms(3) by (intro order.trans[OF pi_approx2[OF assms(1,2)]]) (simp_all add: field_simps) |
655 using assms(3) by (intro order.trans[OF pi_approx2[OF assms(1,2)]]) (simp_all add: field_simps) |
442 |
656 |
443 text \<open> |
657 text \<open> |
444 We can now approximate pi to 54 decimals using this formula. The computations are much |
658 We can now approximate pi to 54 decimals using this formula. The computations are much |
445 slower now; this is mostly because we use arbitrary-precision rational numbers, whose |
659 slower now; this is mostly because we use arbitrary-precision rational numbers, whose |
446 numerators and demoninators get very large. Using dyadic floating point numbers would be |
660 numerators and demoninators get very large. Using dyadic floating point numbers would be |
447 much more economical. |
661 much more economical. |
448 \<close> |
662 \<close> |
449 lemma pi_approx_54_decimals: |
663 lemma pi_approx_54_decimals: |
450 "abs (pi - 3.141592653589793238462643383279502884197169399375105821 :: real) \<le> inverse (10^54)" |
664 "abs (pi - 3.141592653589793238462643383279502884197169399375105821 :: real) \<le> inverse (10^54)" |
451 (is "abs (pi - ?pi') \<le> _") |
665 (is "abs (pi - ?pi') \<le> _") |
452 proof - |
666 proof - |
453 def a \<equiv> "2829469759662002867886529831139137601191652261996513014734415222704732791803 / |
667 def a \<equiv> "2829469759662002867886529831139137601191652261996513014734415222704732791803 / |
454 1062141879292765061960538947347721564047051545995266466660439319087625011200 :: real" |
668 1062141879292765061960538947347721564047051545995266466660439319087625011200 :: real" |
455 def b \<equiv> "13355545553549848714922837267299490903143206628621657811747118592 / |
669 def b \<equiv> "13355545553549848714922837267299490903143206628621657811747118592 / |
456 23792006023392488526789546722992491355941103837356113731091180925 :: real" |
670 23792006023392488526789546722992491355941103837356113731091180925 :: real" |
457 def c \<equiv> "28274063397213534906669125255762067746830085389618481175335056 / |
671 def c \<equiv> "28274063397213534906669125255762067746830085389618481175335056 / |
458 337877029279505250241149903214554249587517250716358486542628059 :: real" |
672 337877029279505250241149903214554249587517250716358486542628059 :: real" |
459 let ?pi'' = "3882327391761098513316067116522233897127356523627918964967729040413954225768920394233198626889767468122598417405434625348404038165437924058179155035564590497837027530349 / |
673 let ?pi'' = "3882327391761098513316067116522233897127356523627918964967729040413954225768920394233198626889767468122598417405434625348404038165437924058179155035564590497837027530349 / |
460 1235783190199688165469648572769847552336447197542738425378629633275352407743112409829873464564018488572820294102599160968781449606552922108667790799771278860366957772800" |
674 1235783190199688165469648572769847552336447197542738425378629633275352407743112409829873464564018488572820294102599160968781449606552922108667790799771278860366957772800" |
461 def eq \<equiv> "op = :: real \<Rightarrow> real \<Rightarrow> bool" |
675 def eq \<equiv> "op = :: real \<Rightarrow> real \<Rightarrow> bool" |
462 |
676 |
463 have "abs (pi - pi_approx2 4) \<le> inverse (2^183)" by (rule pi_approx2') simp_all |
677 have "abs (pi - pi_approx2 4) \<le> inverse (2^183)" by (rule pi_approx2') simp_all |
464 also have "pi_approx2 4 = 48 * arctan_approx 24 (1 / 18) + |
678 also have "pi_approx2 4 = 48 * arctan_approx 24 (1 / 18) + |
465 32 * arctan_approx 16 (1 / 57) - |
679 32 * arctan_approx 16 (1 / 57) - |
466 20 * arctan_approx 12 (1 / 239)" |
680 20 * arctan_approx 12 (1 / 239)" |
467 unfolding pi_approx2_def by simp |
681 unfolding pi_approx2_def by simp |
468 also have [unfolded eq_def]: "eq (48 * arctan_approx 24 (1 / 18)) a" |
682 also have [unfolded eq_def]: "eq (48 * arctan_approx 24 (1 / 18)) a" |
469 by (simp add: arctan_approx_def' power2_eq_square, |
683 by (simp add: arctan_approx_def' power2_eq_square, |
470 simp add: expand_arctan_approx, unfold a_def eq_def, rule refl) |
684 simp add: expand_arctan_approx, unfold a_def eq_def, rule refl) |
471 also have [unfolded eq_def]: "eq (32 * arctan_approx 16 (1 / 57::real)) b" |
685 also have [unfolded eq_def]: "eq (32 * arctan_approx 16 (1 / 57::real)) b" |