130 "get_even n = (if even n then n else (Suc n))" |
130 "get_even n = (if even n then n else (Suc n))" |
131 |
131 |
132 lemma get_odd[simp]: "odd (get_odd n)" unfolding get_odd_def by (cases "odd n", auto) |
132 lemma get_odd[simp]: "odd (get_odd n)" unfolding get_odd_def by (cases "odd n", auto) |
133 lemma get_even[simp]: "even (get_even n)" unfolding get_even_def by (cases "even n", auto) |
133 lemma get_even[simp]: "even (get_even n)" unfolding get_even_def by (cases "even n", auto) |
134 lemma get_odd_ex: "\<exists> k. Suc k = get_odd n \<and> odd (Suc k)" |
134 lemma get_odd_ex: "\<exists> k. Suc k = get_odd n \<and> odd (Suc k)" |
135 proof (cases "odd n") |
135 by (auto simp: get_odd_def odd_pos intro!: exI[of _ "n - 1"]) |
136 case True hence "0 < n" by (rule odd_pos) |
|
137 from gr0_implies_Suc[OF this] obtain k where "Suc k = n" by auto |
|
138 thus ?thesis unfolding get_odd_def if_P[OF True] using True[unfolded `Suc k = n`[symmetric]] by blast |
|
139 next |
|
140 case False hence "odd (Suc n)" by auto |
|
141 thus ?thesis unfolding get_odd_def if_not_P[OF False] by blast |
|
142 qed |
|
143 |
136 |
144 lemma get_even_double: "\<exists>i. get_even n = 2 * i" using get_even[unfolded even_mult_two_ex] . |
137 lemma get_even_double: "\<exists>i. get_even n = 2 * i" using get_even[unfolded even_mult_two_ex] . |
145 lemma get_odd_double: "\<exists>i. get_odd n = 2 * i + 1" using get_odd[unfolded odd_Suc_mult_two_ex] by auto |
138 lemma get_odd_double: "\<exists>i. get_odd n = 2 * i + 1" using get_odd[unfolded odd_Suc_mult_two_ex] by auto |
146 |
139 |
147 section "Power function" |
140 section "Power function" |
149 definition float_power_bnds :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where |
142 definition float_power_bnds :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where |
150 "float_power_bnds n l u = (if odd n \<or> 0 < l then (l ^ n, u ^ n) |
143 "float_power_bnds n l u = (if odd n \<or> 0 < l then (l ^ n, u ^ n) |
151 else if u < 0 then (u ^ n, l ^ n) |
144 else if u < 0 then (u ^ n, l ^ n) |
152 else (0, (max (-l) u) ^ n))" |
145 else (0, (max (-l) u) ^ n))" |
153 |
146 |
154 lemma float_power_bnds: fixes x :: real |
147 lemma float_power_bnds: "(l1, u1) = float_power_bnds n l u \<Longrightarrow> x \<in> {l .. u} \<Longrightarrow> (x::real) ^ n \<in> {l1..u1}" |
155 assumes "(l1, u1) = float_power_bnds n l u" and "x \<in> {l .. u}" |
148 by (auto simp: float_power_bnds_def max_def split: split_if_asm |
156 shows "x ^ n \<in> {l1..u1}" |
149 intro: power_mono_odd power_mono power_mono_even zero_le_even_power) |
157 proof (cases "even n") |
|
158 case True |
|
159 show ?thesis |
|
160 proof (cases "0 < l") |
|
161 case True hence "odd n \<or> 0 < l" and "0 \<le> real l" by auto |
|
162 have u1: "u1 = u ^ n" and l1: "l1 = l ^ n" using assms |
|
163 unfolding float_power_bnds_def if_P[OF `odd n \<or> 0 < l`] by auto |
|
164 have "real l ^ n \<le> x ^ n" and "x ^ n \<le> real u ^ n " using `0 \<le> real l` assms |
|
165 by (auto simp: power_mono) |
|
166 thus ?thesis using assms `0 < l` unfolding l1 u1 by auto |
|
167 next |
|
168 case False hence P: "\<not> (odd n \<or> 0 < l)" using `even n` by auto |
|
169 show ?thesis |
|
170 proof (cases "u < 0") |
|
171 case True hence "0 \<le> - real u" and "- real u \<le> - x" and "0 \<le> - x" and "-x \<le> - real l" using assms by auto |
|
172 hence "real u ^ n \<le> x ^ n" and "x ^ n \<le> real l ^ n" using power_mono[of "-x" "-real l" n] power_mono[of "-real u" "-x" n] |
|
173 unfolding power_minus_even[OF `even n`] by auto |
|
174 moreover have u1: "u1 = l ^ n" and l1: "l1 = u ^ n" using assms unfolding float_power_bnds_def if_not_P[OF P] if_P[OF True] by auto |
|
175 ultimately show ?thesis by auto |
|
176 next |
|
177 case False |
|
178 have "\<bar>x\<bar> \<le> real (max (-l) u)" |
|
179 proof (cases "-l \<le> u") |
|
180 case True thus ?thesis unfolding max_def if_P[OF True] using assms by auto |
|
181 next |
|
182 case False thus ?thesis unfolding max_def if_not_P[OF False] using assms by auto |
|
183 qed |
|
184 hence x_abs: "\<bar>x\<bar> \<le> \<bar>real (max (-l) u)\<bar>" by auto |
|
185 have u1: "u1 = (max (-l) u) ^ n" and l1: "l1 = 0" using assms unfolding float_power_bnds_def if_not_P[OF P] if_not_P[OF False] by auto |
|
186 show ?thesis unfolding atLeastAtMost_iff l1 u1 using zero_le_even_power[OF `even n`] power_mono_even[OF `even n` x_abs] by auto |
|
187 qed |
|
188 qed |
|
189 next |
|
190 case False hence "odd n \<or> 0 < l" by auto |
|
191 have u1: "u1 = u ^ n" and l1: "l1 = l ^ n" using assms unfolding float_power_bnds_def if_P[OF `odd n \<or> 0 < l`] by auto |
|
192 have "real l ^ n \<le> x ^ n" and "x ^ n \<le> real u ^ n " using assms unfolding atLeastAtMost_iff using power_mono_odd[OF False] by auto |
|
193 thus ?thesis unfolding atLeastAtMost_iff l1 u1 less_float_def by auto |
|
194 qed |
|
195 |
150 |
196 lemma bnds_power: "\<forall> (x::real) l u. (l1, u1) = float_power_bnds n l u \<and> x \<in> {l .. u} \<longrightarrow> l1 \<le> x ^ n \<and> x ^ n \<le> u1" |
151 lemma bnds_power: "\<forall> (x::real) l u. (l1, u1) = float_power_bnds n l u \<and> x \<in> {l .. u} \<longrightarrow> l1 \<le> x ^ n \<and> x ^ n \<le> u1" |
197 using float_power_bnds by auto |
152 using float_power_bnds by auto |
198 |
153 |
199 section "Square root" |
154 section "Square root" |
354 hence "sqrt x \<le> ub_sqrt prec x" |
309 hence "sqrt x \<le> ub_sqrt prec x" |
355 unfolding ub_sqrt.simps if_P[OF `0 < x`] by auto } |
310 unfolding ub_sqrt.simps if_P[OF `0 < x`] by auto } |
356 note ub = this |
311 note ub = this |
357 |
312 |
358 show ?thesis |
313 show ?thesis |
359 proof (cases "0 < x") |
314 using lb[of "-x"] ub[of "-x"] lb[of x] ub[of x] |
360 case True with lb ub show ?thesis by auto |
315 by (auto simp add: lb_sqrt.simps ub_sqrt.simps real_sqrt_minus) |
361 next case False show ?thesis |
|
362 proof (cases "real x = 0") |
|
363 case True thus ?thesis |
|
364 by (auto simp add: lb_sqrt.simps ub_sqrt.simps) |
|
365 next |
|
366 case False with `\<not> 0 < x` have "x < 0" and "0 < -x" |
|
367 by auto |
|
368 |
|
369 with `\<not> 0 < x` |
|
370 show ?thesis using lb[OF `0 < -x`] ub[OF `0 < -x`] |
|
371 by (auto simp add: real_sqrt_minus lb_sqrt.simps ub_sqrt.simps) |
|
372 qed qed |
|
373 qed |
316 qed |
374 |
317 |
375 lemma bnds_sqrt: "\<forall> (x::real) lx ux. (l, u) = (lb_sqrt prec lx, ub_sqrt prec ux) \<and> x \<in> {lx .. ux} \<longrightarrow> l \<le> sqrt x \<and> sqrt x \<le> u" |
318 lemma bnds_sqrt: "\<forall> (x::real) lx ux. (l, u) = (lb_sqrt prec lx, ub_sqrt prec ux) \<and> x \<in> {lx .. ux} \<longrightarrow> l \<le> sqrt x \<and> sqrt x \<le> u" |
376 proof ((rule allI) +, rule impI, erule conjE, rule conjI) |
319 proof ((rule allI) +, rule impI, erule conjE, rule conjI) |
377 fix x :: real fix lx ux |
320 fix x :: real fix lx ux |
410 |
353 |
411 lemma arctan_0_1_bounds': |
354 lemma arctan_0_1_bounds': |
412 assumes "0 \<le> real x" "real x \<le> 1" and "even n" |
355 assumes "0 \<le> real x" "real x \<le> 1" and "even n" |
413 shows "arctan x \<in> {(x * lb_arctan_horner prec n 1 (x * x)) .. (x * ub_arctan_horner prec (Suc n) 1 (x * x))}" |
356 shows "arctan x \<in> {(x * lb_arctan_horner prec n 1 (x * x)) .. (x * ub_arctan_horner prec (Suc n) 1 (x * x))}" |
414 proof - |
357 proof - |
415 let "?c i" = "-1^i * (1 / (i * 2 + (1::nat)) * real x ^ (i * 2 + 1))" |
358 let ?c = "\<lambda>i. -1^i * (1 / (i * 2 + (1::nat)) * real x ^ (i * 2 + 1))" |
416 let "?S n" = "\<Sum> i=0..<n. ?c i" |
359 let ?S = "\<lambda>n. \<Sum> i=0..<n. ?c i" |
417 |
360 |
418 have "0 \<le> real (x * x)" by auto |
361 have "0 \<le> real (x * x)" by auto |
419 from `even n` obtain m where "2 * m = n" unfolding even_mult_two_ex by auto |
362 from `even n` obtain m where "2 * m = n" unfolding even_mult_two_ex by auto |
420 |
363 |
421 have "arctan x \<in> { ?S n .. ?S (Suc n) }" |
364 have "arctan x \<in> { ?S n .. ?S (Suc n) }" |
455 ultimately show ?thesis by auto |
398 ultimately show ?thesis by auto |
456 qed |
399 qed |
457 |
400 |
458 lemma arctan_0_1_bounds: assumes "0 \<le> real x" "real x \<le> 1" |
401 lemma arctan_0_1_bounds: assumes "0 \<le> real x" "real x \<le> 1" |
459 shows "arctan x \<in> {(x * lb_arctan_horner prec (get_even n) 1 (x * x)) .. (x * ub_arctan_horner prec (get_odd n) 1 (x * x))}" |
402 shows "arctan x \<in> {(x * lb_arctan_horner prec (get_even n) 1 (x * x)) .. (x * ub_arctan_horner prec (get_odd n) 1 (x * x))}" |
460 proof (cases "even n") |
403 using |
461 case True |
404 arctan_0_1_bounds'[OF assms, of n prec] |
462 obtain n' where "Suc n' = get_odd n" and "odd (Suc n')" using get_odd_ex by auto |
405 arctan_0_1_bounds'[OF assms, of "n + 1" prec] |
463 hence "even n'" unfolding even_Suc by auto |
406 arctan_0_1_bounds'[OF assms, of "n - 1" prec] |
464 have "arctan x \<le> x * ub_arctan_horner prec (get_odd n) 1 (x * x)" |
407 by (auto simp: get_even_def get_odd_def odd_pos simp del: ub_arctan_horner.simps lb_arctan_horner.simps) |
465 unfolding `Suc n' = get_odd n`[symmetric] using arctan_0_1_bounds'[OF `0 \<le> real x` `real x \<le> 1` `even n'`] by auto |
|
466 moreover |
|
467 have "x * lb_arctan_horner prec (get_even n) 1 (x * x) \<le> arctan x" |
|
468 unfolding get_even_def if_P[OF True] using arctan_0_1_bounds'[OF `0 \<le> real x` `real x \<le> 1` `even n`] by auto |
|
469 ultimately show ?thesis by auto |
|
470 next |
|
471 case False hence "0 < n" by (rule odd_pos) |
|
472 from gr0_implies_Suc[OF this] obtain n' where "n = Suc n'" .. |
|
473 from False[unfolded this even_Suc] |
|
474 have "even n'" and "even (Suc (Suc n'))" by auto |
|
475 have "get_odd n = Suc n'" unfolding get_odd_def if_P[OF False] using `n = Suc n'` . |
|
476 |
|
477 have "arctan x \<le> x * ub_arctan_horner prec (get_odd n) 1 (x * x)" |
|
478 unfolding `get_odd n = Suc n'` using arctan_0_1_bounds'[OF `0 \<le> real x` `real x \<le> 1` `even n'`] by auto |
|
479 moreover |
|
480 have "(x * lb_arctan_horner prec (get_even n) 1 (x * x)) \<le> arctan x" |
|
481 unfolding get_even_def if_not_P[OF False] unfolding `n = Suc n'` using arctan_0_1_bounds'[OF `0 \<le> real x` `real x \<le> 1` `even (Suc (Suc n'))`] by auto |
|
482 ultimately show ?thesis by auto |
|
483 qed |
|
484 |
408 |
485 subsection "Compute \<pi>" |
409 subsection "Compute \<pi>" |
486 |
410 |
487 definition ub_pi :: "nat \<Rightarrow> float" where |
411 definition ub_pi :: "nat \<Rightarrow> float" where |
488 "ub_pi prec = (let A = rapprox_rat prec 1 5 ; |
412 "ub_pi prec = (let A = rapprox_rat prec 1 5 ; |
528 using arctan_0_1_bounds[OF `0 \<le> real ?k` `real ?k \<le> 1`] by auto |
452 using arctan_0_1_bounds[OF `0 \<le> real ?k` `real ?k \<le> 1`] by auto |
529 also have "\<dots> \<le> arctan (1 / k)" using `?k \<le> 1 / k` by (rule arctan_monotone') |
453 also have "\<dots> \<le> arctan (1 / k)" using `?k \<le> 1 / k` by (rule arctan_monotone') |
530 finally have "?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k) \<le> arctan (1 / k)" . |
454 finally have "?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k) \<le> arctan (1 / k)" . |
531 } note lb_arctan = this |
455 } note lb_arctan = this |
532 |
456 |
533 have "pi \<le> ub_pi n" |
457 have "pi \<le> ub_pi n \<and> lb_pi n \<le> pi" |
534 unfolding ub_pi_def machin_pi Let_def unfolding Float_num |
458 unfolding lb_pi_def ub_pi_def machin_pi Let_def unfolding Float_num |
535 using lb_arctan[of 239] ub_arctan[of 5] powr_realpow[of 2 2] |
459 using lb_arctan[of 5] ub_arctan[of 239] lb_arctan[of 239] ub_arctan[of 5] powr_realpow[of 2 2] |
536 by (auto intro!: mult_left_mono add_mono simp add: uminus_add_conv_diff [symmetric] simp del: uminus_add_conv_diff) |
460 by (auto intro!: mult_left_mono add_mono simp add: uminus_add_conv_diff [symmetric] simp del: uminus_add_conv_diff) |
537 moreover |
461 then show ?thesis by auto |
538 have "lb_pi n \<le> pi" |
|
539 unfolding lb_pi_def machin_pi Let_def Float_num |
|
540 using lb_arctan[of 5] ub_arctan[of 239] powr_realpow[of 2 2] |
|
541 by (auto intro!: mult_left_mono add_mono simp add: uminus_add_conv_diff [symmetric] simp del: uminus_add_conv_diff) |
|
542 ultimately show ?thesis by auto |
|
543 qed |
462 qed |
544 |
463 |
545 subsection "Compute arcus tangens in the entire domain" |
464 subsection "Compute arcus tangens in the entire domain" |
546 |
465 |
547 function lb_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" and ub_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" where |
466 function lb_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" and ub_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" where |
1806 using powr_gt_zero[of 2 "exponent x"] |
1725 using powr_gt_zero[of 2 "exponent x"] |
1807 apply simp |
1726 apply simp |
1808 done |
1727 done |
1809 |
1728 |
1810 lemma Float_pos_eq_mantissa_pos: "Float m e > 0 \<longleftrightarrow> m > 0" |
1729 lemma Float_pos_eq_mantissa_pos: "Float m e > 0 \<longleftrightarrow> m > 0" |
1811 apply (auto simp add: zero_less_mult_iff zero_float_def powr_gt_zero[of 2 "exponent x"] dest: less_zeroE) |
|
1812 using powr_gt_zero[of 2 "e"] |
1730 using powr_gt_zero[of 2 "e"] |
1813 apply simp |
1731 by (auto simp add: zero_less_mult_iff zero_float_def simp del: powr_gt_zero dest: less_zeroE) |
1814 done |
|
1815 |
1732 |
1816 lemma Float_representation_aux: |
1733 lemma Float_representation_aux: |
1817 fixes m e |
1734 fixes m e |
1818 defines "x \<equiv> Float m e" |
1735 defines "x \<equiv> Float m e" |
1819 assumes "x > 0" |
1736 assumes "x > 0" |