1 (* Author: Johannes Hoelzl <hoelzl@in.tum.de> 2008 / 2009 *)
3 header {* Prove Real Valued Inequalities by Computation *}
6 imports Complex_Main Float Reflection Dense_Linear_Order Efficient_Nat
9 section "Horner Scheme"
11 subsection {* Define auxiliary helper @{text horner} function *}
13 primrec horner :: "(nat \<Rightarrow> nat) \<Rightarrow> (nat \<Rightarrow> nat \<Rightarrow> nat) \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> real \<Rightarrow> real" where
14 "horner F G 0 i k x = 0" |
15 "horner F G (Suc n) i k x = 1 / real k - x * horner F G n (F i) (G i k) x"
17 lemma horner_schema': fixes x :: real and a :: "nat \<Rightarrow> real"
18 shows "a 0 - x * (\<Sum> i=0..<n. (-1)^i * a (Suc i) * x^i) = (\<Sum> i=0..<Suc n. (-1)^i * a i * x^i)"
20 have shift_pow: "\<And>i. - (x * ((-1)^i * a (Suc i) * x ^ i)) = (-1)^(Suc i) * a (Suc i) * x ^ (Suc i)" by auto
21 show ?thesis unfolding setsum_right_distrib shift_pow real_diff_def setsum_negf[symmetric] setsum_head_upt_Suc[OF zero_less_Suc]
22 setsum_reindex[OF inj_Suc, unfolded comp_def, symmetric, of "\<lambda> n. (-1)^n *a n * x^n"] by auto
25 lemma horner_schema: fixes f :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat" and F :: "nat \<Rightarrow> nat"
26 assumes f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"
27 shows "horner F G n ((F ^^ j') s) (f j') x = (\<Sum> j = 0..< n. -1 ^ j * (1 / real (f (j' + j))) * x ^ j)"
28 proof (induct n arbitrary: i k j')
31 show ?case unfolding horner.simps Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc]
32 using horner_schema'[of "\<lambda> j. 1 / real (f (j' + j))"] by auto
36 assumes "0 \<le> real x" and f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"
37 and lb_0: "\<And> i k x. lb 0 i k x = 0"
38 and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 (int k) - x * (ub n (F i) (G i k) x)"
39 and ub_0: "\<And> i k x. ub 0 i k x = 0"
40 and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 (int k) - x * (lb n (F i) (G i k) x)"
41 shows "real (lb n ((F ^^ j') s) (f j') x) \<le> horner F G n ((F ^^ j') s) (f j') (real x) \<and>
42 horner F G n ((F ^^ j') s) (f j') (real x) \<le> real (ub n ((F ^^ j') s) (f j') x)"
43 (is "?lb n j' \<le> ?horner n j' \<and> ?horner n j' \<le> ?ub n j'")
44 proof (induct n arbitrary: j')
45 case 0 thus ?case unfolding lb_0 ub_0 horner.simps by auto
48 have "?lb (Suc n) j' \<le> ?horner (Suc n) j'" unfolding lb_Suc ub_Suc horner.simps real_of_float_sub diff_def
50 show "real (lapprox_rat prec 1 (int (f j'))) \<le> 1 / real (f j')" using lapprox_rat[of prec 1 "int (f j')"] by auto
51 from Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc, THEN conjunct2] `0 \<le> real x`
52 show "- real (x * ub n (F ((F ^^ j') s)) (G ((F ^^ j') s) (f j')) x) \<le> - (real x * horner F G n (F ((F ^^ j') s)) (G ((F ^^ j') s) (f j')) (real x))"
53 unfolding real_of_float_mult neg_le_iff_le by (rule mult_left_mono)
55 moreover have "?horner (Suc n) j' \<le> ?ub (Suc n) j'" unfolding ub_Suc ub_Suc horner.simps real_of_float_sub diff_def
57 show "1 / real (f j') \<le> real (rapprox_rat prec 1 (int (f j')))" using rapprox_rat[of 1 "int (f j')" prec] by auto
58 from Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc, THEN conjunct1] `0 \<le> real x`
59 show "- (real x * horner F G n (F ((F ^^ j') s)) (G ((F ^^ j') s) (f j')) (real x)) \<le>
60 - real (x * lb n (F ((F ^^ j') s)) (G ((F ^^ j') s) (f j')) x)"
61 unfolding real_of_float_mult neg_le_iff_le by (rule mult_left_mono)
63 ultimately show ?case by blast
66 subsection "Theorems for floating point functions implementing the horner scheme"
70 Here @{term_type "f :: nat \<Rightarrow> nat"} is the sequence defining the Taylor series, the coefficients are
71 all alternating and reciprocs. We use @{term G} and @{term F} to describe the computation of @{term f}.
75 lemma horner_bounds: fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat"
76 assumes "0 \<le> real x" and f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"
77 and lb_0: "\<And> i k x. lb 0 i k x = 0"
78 and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 (int k) - x * (ub n (F i) (G i k) x)"
79 and ub_0: "\<And> i k x. ub 0 i k x = 0"
80 and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 (int k) - x * (lb n (F i) (G i k) x)"
81 shows "real (lb n ((F ^^ j') s) (f j') x) \<le> (\<Sum>j=0..<n. -1 ^ j * (1 / real (f (j' + j))) * real x ^ j)" (is "?lb") and
82 "(\<Sum>j=0..<n. -1 ^ j * (1 / real (f (j' + j))) * (real x ^ j)) \<le> real (ub n ((F ^^ j') s) (f j') x)" (is "?ub")
85 using horner_bounds'[where lb=lb, OF `0 \<le> real x` f_Suc lb_0 lb_Suc ub_0 ub_Suc]
86 unfolding horner_schema[where f=f, OF f_Suc] .
87 thus "?lb" and "?ub" by auto
90 lemma horner_bounds_nonpos: fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat"
91 assumes "real x \<le> 0" and f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"
92 and lb_0: "\<And> i k x. lb 0 i k x = 0"
93 and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 (int k) + x * (ub n (F i) (G i k) x)"
94 and ub_0: "\<And> i k x. ub 0 i k x = 0"
95 and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 (int k) + x * (lb n (F i) (G i k) x)"
96 shows "real (lb n ((F ^^ j') s) (f j') x) \<le> (\<Sum>j=0..<n. (1 / real (f (j' + j))) * real x ^ j)" (is "?lb") and
97 "(\<Sum>j=0..<n. (1 / real (f (j' + j))) * (real x ^ j)) \<le> real (ub n ((F ^^ j') s) (f j') x)" (is "?ub")
99 { fix x y z :: float have "x - y * z = x + - y * z"
100 by (cases x, cases y, cases z, simp add: plus_float.simps minus_float_def uminus_float.simps times_float.simps algebra_simps)
101 } note diff_mult_minus = this
103 { fix x :: float have "- (- x) = x" by (cases x, auto simp add: uminus_float.simps) } note minus_minus = this
105 have move_minus: "real (-x) = -1 * real x" by auto
107 have sum_eq: "(\<Sum>j=0..<n. (1 / real (f (j' + j))) * real x ^ j) =
108 (\<Sum>j = 0..<n. -1 ^ j * (1 / real (f (j' + j))) * real (- x) ^ j)"
109 proof (rule setsum_cong, simp)
110 fix j assume "j \<in> {0 ..< n}"
111 show "1 / real (f (j' + j)) * real x ^ j = -1 ^ j * (1 / real (f (j' + j))) * real (- x) ^ j"
112 unfolding move_minus power_mult_distrib real_mult_assoc[symmetric]
113 unfolding real_mult_commute unfolding real_mult_assoc[of "-1 ^ j", symmetric] power_mult_distrib[symmetric]
117 have "0 \<le> real (-x)" using assms by auto
118 from horner_bounds[where G=G and F=F and f=f and s=s and prec=prec
119 and lb="\<lambda> n i k x. lb n i k (-x)" and ub="\<lambda> n i k x. ub n i k (-x)", unfolded lb_Suc ub_Suc diff_mult_minus,
120 OF this f_Suc lb_0 refl ub_0 refl]
121 show "?lb" and "?ub" unfolding minus_minus sum_eq
125 subsection {* Selectors for next even or odd number *}
129 The horner scheme computes alternating series. To get the upper and lower bounds we need to
130 guarantee to access a even or odd member. To do this we use @{term get_odd} and @{term get_even}.
134 definition get_odd :: "nat \<Rightarrow> nat" where
135 "get_odd n = (if odd n then n else (Suc n))"
137 definition get_even :: "nat \<Rightarrow> nat" where
138 "get_even n = (if even n then n else (Suc n))"
140 lemma get_odd[simp]: "odd (get_odd n)" unfolding get_odd_def by (cases "odd n", auto)
141 lemma get_even[simp]: "even (get_even n)" unfolding get_even_def by (cases "even n", auto)
142 lemma get_odd_ex: "\<exists> k. Suc k = get_odd n \<and> odd (Suc k)"
143 proof (cases "odd n")
144 case True hence "0 < n" by (rule odd_pos)
145 from gr0_implies_Suc[OF this] obtain k where "Suc k = n" by auto
146 thus ?thesis unfolding get_odd_def if_P[OF True] using True[unfolded `Suc k = n`[symmetric]] by blast
148 case False hence "odd (Suc n)" by auto
149 thus ?thesis unfolding get_odd_def if_not_P[OF False] by blast
152 lemma get_even_double: "\<exists>i. get_even n = 2 * i" using get_even[unfolded even_mult_two_ex] .
153 lemma get_odd_double: "\<exists>i. get_odd n = 2 * i + 1" using get_odd[unfolded odd_Suc_mult_two_ex] by auto
155 section "Power function"
157 definition float_power_bnds :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where
158 "float_power_bnds n l u = (if odd n \<or> 0 < l then (l ^ n, u ^ n)
159 else if u < 0 then (u ^ n, l ^ n)
160 else (0, (max (-l) u) ^ n))"
162 lemma float_power_bnds: assumes "(l1, u1) = float_power_bnds n l u" and "x \<in> {real l .. real u}"
163 shows "x ^ n \<in> {real l1..real u1}"
164 proof (cases "even n")
167 proof (cases "0 < l")
168 case True hence "odd n \<or> 0 < l" and "0 \<le> real l" unfolding less_float_def by auto
169 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
170 have "real l ^ n \<le> x ^ n" and "x ^ n \<le> real u ^ n " using `0 \<le> real l` and assms unfolding atLeastAtMost_iff using power_mono[of "real l" x] power_mono[of x "real u"] by auto
171 thus ?thesis using assms `0 < l` unfolding atLeastAtMost_iff l1 u1 float_power less_float_def by auto
173 case False hence P: "\<not> (odd n \<or> 0 < l)" using `even n` by auto
175 proof (cases "u < 0")
176 case True hence "0 \<le> - real u" and "- real u \<le> - x" and "0 \<le> - x" and "-x \<le> - real l" using assms unfolding less_float_def by auto
177 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]
178 unfolding power_minus_even[OF `even n`] by auto
179 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
180 ultimately show ?thesis using float_power by auto
183 have "\<bar>x\<bar> \<le> real (max (-l) u)"
184 proof (cases "-l \<le> u")
185 case True thus ?thesis unfolding max_def if_P[OF True] using assms unfolding le_float_def by auto
187 case False thus ?thesis unfolding max_def if_not_P[OF False] using assms unfolding le_float_def by auto
189 hence x_abs: "\<bar>x\<bar> \<le> \<bar>real (max (-l) u)\<bar>" by auto
190 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
191 show ?thesis unfolding atLeastAtMost_iff l1 u1 float_power using zero_le_even_power[OF `even n`] power_mono_even[OF `even n` x_abs] by auto
195 case False hence "odd n \<or> 0 < l" by auto
196 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
197 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
198 thus ?thesis unfolding atLeastAtMost_iff l1 u1 float_power less_float_def by auto
201 lemma bnds_power: "\<forall> x l u. (l1, u1) = float_power_bnds n l u \<and> x \<in> {real l .. real u} \<longrightarrow> real l1 \<le> x ^ n \<and> x ^ n \<le> real u1"
202 using float_power_bnds by auto
204 section "Square root"
208 The square root computation is implemented as newton iteration. As first first step we use the
209 nearest power of two greater than the square root.
213 fun sqrt_iteration :: "nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
214 "sqrt_iteration prec 0 (Float m e) = Float 1 ((e + bitlen m) div 2 + 1)" |
215 "sqrt_iteration prec (Suc m) x = (let y = sqrt_iteration prec m x
216 in Float 1 -1 * (y + float_divr prec x y))"
218 function ub_sqrt lb_sqrt :: "nat \<Rightarrow> float \<Rightarrow> float" where
219 "ub_sqrt prec x = (if 0 < x then (sqrt_iteration prec prec x)
220 else if x < 0 then - lb_sqrt prec (- x)
222 "lb_sqrt prec x = (if 0 < x then (float_divl prec x (sqrt_iteration prec prec x))
223 else if x < 0 then - ub_sqrt prec (- x)
225 by pat_completeness auto
226 termination by (relation "measure (\<lambda> v. let (prec, x) = sum_case id id v in (if x < 0 then 1 else 0))", auto simp add: less_float_def)
228 declare lb_sqrt.simps[simp del]
229 declare ub_sqrt.simps[simp del]
231 lemma sqrt_ub_pos_pos_1:
232 assumes "sqrt x < b" and "0 < b" and "0 < x"
233 shows "sqrt x < (b + x / b)/2"
235 from assms have "0 < (b - sqrt x) ^ 2 " by simp
236 also have "\<dots> = b ^ 2 - 2 * b * sqrt x + (sqrt x) ^ 2" by algebra
237 also have "\<dots> = b ^ 2 - 2 * b * sqrt x + x" using assms by (simp add: real_sqrt_pow2)
238 finally have "0 < b ^ 2 - 2 * b * sqrt x + x" by assumption
239 hence "0 < b / 2 - sqrt x + x / (2 * b)" using assms
240 by (simp add: field_simps power2_eq_square)
241 thus ?thesis by (simp add: field_simps)
244 lemma sqrt_iteration_bound: assumes "0 < real x"
245 shows "sqrt (real x) < real (sqrt_iteration prec n x)"
251 hence "0 < m" using float_pos_m_pos[unfolded less_float_def] assms by auto
252 hence "0 < sqrt (real m)" by auto
254 have int_nat_bl: "int (nat (bitlen m)) = bitlen m" using bitlen_ge0 by auto
256 have "real x = (real m / 2^nat (bitlen m)) * pow2 (e + int (nat (bitlen m)))"
257 unfolding pow2_add pow2_int Float real_of_float_simp by auto
258 also have "\<dots> < 1 * pow2 (e + int (nat (bitlen m)))"
259 proof (rule mult_strict_right_mono, auto)
260 show "real m < 2^nat (bitlen m)" using bitlen_bounds[OF `0 < m`, THEN conjunct2]
261 unfolding real_of_int_less_iff[of m, symmetric] by auto
263 finally have "sqrt (real x) < sqrt (pow2 (e + bitlen m))" unfolding int_nat_bl by auto
264 also have "\<dots> \<le> pow2 ((e + bitlen m) div 2 + 1)"
266 let ?E = "e + bitlen m"
267 have E_mod_pow: "pow2 (?E mod 2) < 4"
268 proof (cases "?E mod 2 = 1")
269 case True thus ?thesis by auto
272 have "0 \<le> ?E mod 2" by auto
273 have "?E mod 2 < 2" by auto
274 from this[THEN zless_imp_add1_zle]
275 have "?E mod 2 \<le> 0" using False by auto
276 from xt1(5)[OF `0 \<le> ?E mod 2` this]
279 hence "sqrt (pow2 (?E mod 2)) < sqrt (2 * 2)" by auto
280 hence E_mod_pow: "sqrt (pow2 (?E mod 2)) < 2" unfolding real_sqrt_abs2 by auto
282 have E_eq: "pow2 ?E = pow2 (?E div 2 + ?E div 2 + ?E mod 2)" by auto
283 have "sqrt (pow2 ?E) = sqrt (pow2 (?E div 2) * pow2 (?E div 2) * pow2 (?E mod 2))"
284 unfolding E_eq unfolding pow2_add ..
285 also have "\<dots> = pow2 (?E div 2) * sqrt (pow2 (?E mod 2))"
286 unfolding real_sqrt_mult[of _ "pow2 (?E mod 2)"] real_sqrt_abs2 by auto
287 also have "\<dots> < pow2 (?E div 2) * 2"
288 by (rule mult_strict_left_mono, auto intro: E_mod_pow)
289 also have "\<dots> = pow2 (?E div 2 + 1)" unfolding zadd_commute[of _ 1] pow2_add1 by auto
290 finally show ?thesis by auto
293 unfolding Float sqrt_iteration.simps real_of_float_simp by auto
297 let ?b = "sqrt_iteration prec n x"
298 have "0 < sqrt (real x)" using `0 < real x` by auto
299 also have "\<dots> < real ?b" using Suc .
300 finally have "sqrt (real x) < (real ?b + real x / real ?b)/2" using sqrt_ub_pos_pos_1[OF Suc _ `0 < real x`] by auto
301 also have "\<dots> \<le> (real ?b + real (float_divr prec x ?b))/2" by (rule divide_right_mono, auto simp add: float_divr)
302 also have "\<dots> = real (Float 1 -1) * (real ?b + real (float_divr prec x ?b))" by auto
303 finally show ?case unfolding sqrt_iteration.simps Let_def real_of_float_mult real_of_float_add right_distrib .
306 lemma sqrt_iteration_lower_bound: assumes "0 < real x"
307 shows "0 < real (sqrt_iteration prec n x)" (is "0 < ?sqrt")
309 have "0 < sqrt (real x)" using assms by auto
310 also have "\<dots> < ?sqrt" using sqrt_iteration_bound[OF assms] .
311 finally show ?thesis .
314 lemma lb_sqrt_lower_bound: assumes "0 \<le> real x"
315 shows "0 \<le> real (lb_sqrt prec x)"
316 proof (cases "0 < x")
317 case True hence "0 < real x" and "0 \<le> x" using `0 \<le> real x` unfolding less_float_def le_float_def by auto
318 hence "0 < sqrt_iteration prec prec x" unfolding less_float_def using sqrt_iteration_lower_bound by auto
319 hence "0 \<le> real (float_divl prec x (sqrt_iteration prec prec x))" using float_divl_lower_bound[OF `0 \<le> x`] unfolding le_float_def by auto
320 thus ?thesis unfolding lb_sqrt.simps using True by auto
322 case False with `0 \<le> real x` have "real x = 0" unfolding less_float_def by auto
323 thus ?thesis unfolding lb_sqrt.simps less_float_def by auto
327 shows "sqrt (real x) \<in> { real (lb_sqrt prec x) .. real (ub_sqrt prec x) }"
329 { fix x :: float assume "0 < x"
330 hence "0 < real x" and "0 \<le> real x" unfolding less_float_def by auto
331 hence sqrt_gt0: "0 < sqrt (real x)" by auto
332 hence sqrt_ub: "sqrt (real x) < real (sqrt_iteration prec prec x)" using sqrt_iteration_bound by auto
334 have "real (float_divl prec x (sqrt_iteration prec prec x)) \<le>
335 real x / real (sqrt_iteration prec prec x)" by (rule float_divl)
336 also have "\<dots> < real x / sqrt (real x)"
337 by (rule divide_strict_left_mono[OF sqrt_ub `0 < real x`
338 mult_pos_pos[OF order_less_trans[OF sqrt_gt0 sqrt_ub] sqrt_gt0]])
339 also have "\<dots> = sqrt (real x)"
340 unfolding inverse_eq_iff_eq[of _ "sqrt (real x)", symmetric]
341 sqrt_divide_self_eq[OF `0 \<le> real x`, symmetric] by auto
342 finally have "real (lb_sqrt prec x) \<le> sqrt (real x)"
343 unfolding lb_sqrt.simps if_P[OF `0 < x`] by auto }
346 { fix x :: float assume "0 < x"
347 hence "0 < real x" unfolding less_float_def by auto
348 hence "0 < sqrt (real x)" by auto
349 hence "sqrt (real x) < real (sqrt_iteration prec prec x)"
350 using sqrt_iteration_bound by auto
351 hence "sqrt (real x) \<le> real (ub_sqrt prec x)"
352 unfolding ub_sqrt.simps if_P[OF `0 < x`] by auto }
356 proof (cases "0 < x")
357 case True with lb ub show ?thesis by auto
358 next case False show ?thesis
359 proof (cases "real x = 0")
360 case True thus ?thesis
361 by (auto simp add: less_float_def lb_sqrt.simps ub_sqrt.simps)
363 case False with `\<not> 0 < x` have "x < 0" and "0 < -x"
364 by (auto simp add: less_float_def)
367 show ?thesis using lb[OF `0 < -x`] ub[OF `0 < -x`]
368 by (auto simp add: real_sqrt_minus lb_sqrt.simps ub_sqrt.simps)
372 lemma bnds_sqrt: "\<forall> x lx ux. (l, u) = (lb_sqrt prec lx, ub_sqrt prec ux) \<and> x \<in> {real lx .. real ux} \<longrightarrow> real l \<le> sqrt x \<and> sqrt x \<le> real u"
373 proof ((rule allI) +, rule impI, erule conjE, rule conjI)
375 assume "(l, u) = (lb_sqrt prec lx, ub_sqrt prec ux)"
376 and x: "x \<in> {real lx .. real ux}"
377 hence l: "l = lb_sqrt prec lx " and u: "u = ub_sqrt prec ux" by auto
379 have "sqrt (real lx) \<le> sqrt x" using x by auto
380 from order_trans[OF _ this]
381 show "real l \<le> sqrt x" unfolding l using bnds_sqrt'[of lx prec] by auto
383 have "sqrt x \<le> sqrt (real ux)" using x by auto
384 from order_trans[OF this]
385 show "sqrt x \<le> real u" unfolding u using bnds_sqrt'[of ux prec] by auto
388 section "Arcus tangens and \<pi>"
390 subsection "Compute arcus tangens series"
394 As first step we implement the computation of the arcus tangens series. This is only valid in the range
395 @{term "{-1 :: real .. 1}"}. This is used to compute \<pi> and then the entire arcus tangens.
399 fun ub_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
400 and lb_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
401 "ub_arctan_horner prec 0 k x = 0"
402 | "ub_arctan_horner prec (Suc n) k x =
403 (rapprox_rat prec 1 (int k)) - x * (lb_arctan_horner prec n (k + 2) x)"
404 | "lb_arctan_horner prec 0 k x = 0"
405 | "lb_arctan_horner prec (Suc n) k x =
406 (lapprox_rat prec 1 (int k)) - x * (ub_arctan_horner prec n (k + 2) x)"
408 lemma arctan_0_1_bounds': assumes "0 \<le> real x" "real x \<le> 1" and "even n"
409 shows "arctan (real x) \<in> {real (x * lb_arctan_horner prec n 1 (x * x)) .. real (x * ub_arctan_horner prec (Suc n) 1 (x * x))}"
411 let "?c i" = "-1^i * (1 / real (i * 2 + 1) * real x ^ (i * 2 + 1))"
412 let "?S n" = "\<Sum> i=0..<n. ?c i"
414 have "0 \<le> real (x * x)" by auto
415 from `even n` obtain m where "2 * m = n" unfolding even_mult_two_ex by auto
417 have "arctan (real x) \<in> { ?S n .. ?S (Suc n) }"
418 proof (cases "real x = 0")
420 hence "0 < real x" using `0 \<le> real x` by auto
421 hence prem: "0 < 1 / real (0 * 2 + (1::nat)) * real x ^ (0 * 2 + 1)" by auto
423 have "\<bar> real x \<bar> \<le> 1" using `0 \<le> real x` `real x \<le> 1` by auto
424 from mp[OF summable_Leibniz(2)[OF zeroseq_arctan_series[OF this] monoseq_arctan_series[OF this]] prem, THEN spec, of m, unfolded `2 * m = n`]
425 show ?thesis unfolding arctan_series[OF `\<bar> real x \<bar> \<le> 1`] Suc_eq_plus1 .
427 note arctan_bounds = this[unfolded atLeastAtMost_iff]
429 have F: "\<And>n. 2 * Suc n + 1 = 2 * n + 1 + 2" by auto
431 note bounds = horner_bounds[where s=1 and f="\<lambda>i. 2 * i + 1" and j'=0
432 and lb="\<lambda>n i k x. lb_arctan_horner prec n k x"
433 and ub="\<lambda>n i k x. ub_arctan_horner prec n k x",
434 OF `0 \<le> real (x*x)` F lb_arctan_horner.simps ub_arctan_horner.simps]
436 { have "real (x * lb_arctan_horner prec n 1 (x*x)) \<le> ?S n"
437 using bounds(1) `0 \<le> real x`
438 unfolding real_of_float_mult power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]
439 unfolding real_mult_commute mult_commute[of _ "2::nat"] power_mult power2_eq_square[of "real x"]
440 by (auto intro!: mult_left_mono)
441 also have "\<dots> \<le> arctan (real x)" using arctan_bounds ..
442 finally have "real (x * lb_arctan_horner prec n 1 (x*x)) \<le> arctan (real x)" . }
444 { have "arctan (real x) \<le> ?S (Suc n)" using arctan_bounds ..
445 also have "\<dots> \<le> real (x * ub_arctan_horner prec (Suc n) 1 (x*x))"
446 using bounds(2)[of "Suc n"] `0 \<le> real x`
447 unfolding real_of_float_mult power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]
448 unfolding real_mult_commute mult_commute[of _ "2::nat"] power_mult power2_eq_square[of "real x"]
449 by (auto intro!: mult_left_mono)
450 finally have "arctan (real x) \<le> real (x * ub_arctan_horner prec (Suc n) 1 (x*x))" . }
451 ultimately show ?thesis by auto
454 lemma arctan_0_1_bounds: assumes "0 \<le> real x" "real x \<le> 1"
455 shows "arctan (real x) \<in> {real (x * lb_arctan_horner prec (get_even n) 1 (x * x)) .. real (x * ub_arctan_horner prec (get_odd n) 1 (x * x))}"
456 proof (cases "even n")
458 obtain n' where "Suc n' = get_odd n" and "odd (Suc n')" using get_odd_ex by auto
459 hence "even n'" unfolding even_Suc by auto
460 have "arctan (real x) \<le> real (x * ub_arctan_horner prec (get_odd n) 1 (x * x))"
461 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
463 have "real (x * lb_arctan_horner prec (get_even n) 1 (x * x)) \<le> arctan (real x)"
464 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
465 ultimately show ?thesis by auto
467 case False hence "0 < n" by (rule odd_pos)
468 from gr0_implies_Suc[OF this] obtain n' where "n = Suc n'" ..
469 from False[unfolded this even_Suc]
470 have "even n'" and "even (Suc (Suc n'))" by auto
471 have "get_odd n = Suc n'" unfolding get_odd_def if_P[OF False] using `n = Suc n'` .
473 have "arctan (real x) \<le> real (x * ub_arctan_horner prec (get_odd n) 1 (x * x))"
474 unfolding `get_odd n = Suc n'` using arctan_0_1_bounds'[OF `0 \<le> real x` `real x \<le> 1` `even n'`] by auto
476 have "real (x * lb_arctan_horner prec (get_even n) 1 (x * x)) \<le> arctan (real x)"
477 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
478 ultimately show ?thesis by auto
481 subsection "Compute \<pi>"
483 definition ub_pi :: "nat \<Rightarrow> float" where
484 "ub_pi prec = (let A = rapprox_rat prec 1 5 ;
485 B = lapprox_rat prec 1 239
486 in ((Float 1 2) * ((Float 1 2) * A * (ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (A * A)) -
487 B * (lb_arctan_horner prec (get_even (prec div 14 + 1)) 1 (B * B)))))"
489 definition lb_pi :: "nat \<Rightarrow> float" where
490 "lb_pi prec = (let A = lapprox_rat prec 1 5 ;
491 B = rapprox_rat prec 1 239
492 in ((Float 1 2) * ((Float 1 2) * A * (lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (A * A)) -
493 B * (ub_arctan_horner prec (get_odd (prec div 14 + 1)) 1 (B * B)))))"
495 lemma pi_boundaries: "pi \<in> {real (lb_pi n) .. real (ub_pi n)}"
497 have machin_pi: "pi = 4 * (4 * arctan (1 / 5) - arctan (1 / 239))" unfolding machin[symmetric] by auto
499 { fix prec n :: nat fix k :: int assume "1 < k" hence "0 \<le> k" and "0 < k" and "1 \<le> k" by auto
500 let ?k = "rapprox_rat prec 1 k"
501 have "1 div k = 0" using div_pos_pos_trivial[OF _ `1 < k`] by auto
503 have "0 \<le> real ?k" by (rule order_trans[OF _ rapprox_rat], auto simp add: `0 \<le> k`)
504 have "real ?k \<le> 1" unfolding rapprox_rat.simps(2)[OF zero_le_one `0 < k`]
505 by (rule rapprox_posrat_le1, auto simp add: `0 < k` `1 \<le> k`)
507 have "1 / real k \<le> real ?k" using rapprox_rat[where x=1 and y=k] by auto
508 hence "arctan (1 / real k) \<le> arctan (real ?k)" by (rule arctan_monotone')
509 also have "\<dots> \<le> real (?k * ub_arctan_horner prec (get_odd n) 1 (?k * ?k))"
510 using arctan_0_1_bounds[OF `0 \<le> real ?k` `real ?k \<le> 1`] by auto
511 finally have "arctan (1 / (real k)) \<le> real (?k * ub_arctan_horner prec (get_odd n) 1 (?k * ?k))" .
512 } note ub_arctan = this
514 { fix prec n :: nat fix k :: int assume "1 < k" hence "0 \<le> k" and "0 < k" by auto
515 let ?k = "lapprox_rat prec 1 k"
516 have "1 div k = 0" using div_pos_pos_trivial[OF _ `1 < k`] by auto
517 have "1 / real k \<le> 1" using `1 < k` by auto
519 have "\<And>n. 0 \<le> real ?k" using lapprox_rat_bottom[where x=1 and y=k, OF zero_le_one `0 < k`] by (auto simp add: `1 div k = 0`)
520 have "\<And>n. real ?k \<le> 1" using lapprox_rat by (rule order_trans, auto simp add: `1 / real k \<le> 1`)
522 have "real ?k \<le> 1 / real k" using lapprox_rat[where x=1 and y=k] by auto
524 have "real (?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k)) \<le> arctan (real ?k)"
525 using arctan_0_1_bounds[OF `0 \<le> real ?k` `real ?k \<le> 1`] by auto
526 also have "\<dots> \<le> arctan (1 / real k)" using `real ?k \<le> 1 / real k` by (rule arctan_monotone')
527 finally have "real (?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k)) \<le> arctan (1 / (real k))" .
528 } note lb_arctan = this
530 have "pi \<le> real (ub_pi n)"
531 unfolding ub_pi_def machin_pi Let_def real_of_float_mult real_of_float_sub unfolding Float_num
532 using lb_arctan[of 239] ub_arctan[of 5]
533 by (auto intro!: mult_left_mono add_mono simp add: diff_minus simp del: lapprox_rat.simps rapprox_rat.simps)
535 have "real (lb_pi n) \<le> pi"
536 unfolding lb_pi_def machin_pi Let_def real_of_float_mult real_of_float_sub Float_num
537 using lb_arctan[of 5] ub_arctan[of 239]
538 by (auto intro!: mult_left_mono add_mono simp add: diff_minus simp del: lapprox_rat.simps rapprox_rat.simps)
539 ultimately show ?thesis by auto
542 subsection "Compute arcus tangens in the entire domain"
544 function lb_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" and ub_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" where
545 "lb_arctan prec x = (let ub_horner = \<lambda> x. x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x) ;
546 lb_horner = \<lambda> x. x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)
547 in (if x < 0 then - ub_arctan prec (-x) else
548 if x \<le> Float 1 -1 then lb_horner x else
549 if x \<le> Float 1 1 then Float 1 1 * lb_horner (float_divl prec x (1 + ub_sqrt prec (1 + x * x)))
550 else (let inv = float_divr prec 1 x
552 else lb_pi prec * Float 1 -1 - ub_horner inv)))"
554 | "ub_arctan prec x = (let lb_horner = \<lambda> x. x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x) ;
555 ub_horner = \<lambda> x. x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)
556 in (if x < 0 then - lb_arctan prec (-x) else
557 if x \<le> Float 1 -1 then ub_horner x else
558 if x \<le> Float 1 1 then let y = float_divr prec x (1 + lb_sqrt prec (1 + x * x))
559 in if y > 1 then ub_pi prec * Float 1 -1
560 else Float 1 1 * ub_horner y
561 else ub_pi prec * Float 1 -1 - lb_horner (float_divl prec 1 x)))"
562 by pat_completeness auto
563 termination by (relation "measure (\<lambda> v. let (prec, x) = sum_case id id v in (if x < 0 then 1 else 0))", auto simp add: less_float_def)
565 declare ub_arctan_horner.simps[simp del]
566 declare lb_arctan_horner.simps[simp del]
568 lemma lb_arctan_bound': assumes "0 \<le> real x"
569 shows "real (lb_arctan prec x) \<le> arctan (real x)"
571 have "\<not> x < 0" and "0 \<le> x" unfolding less_float_def le_float_def using `0 \<le> real x` by auto
572 let "?ub_horner x" = "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)"
573 and "?lb_horner x" = "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)"
576 proof (cases "x \<le> Float 1 -1")
577 case True hence "real x \<le> 1" unfolding le_float_def Float_num by auto
578 show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True]
579 using arctan_0_1_bounds[OF `0 \<le> real x` `real x \<le> 1`] by auto
581 case False hence "0 < real x" unfolding le_float_def Float_num by auto
582 let ?R = "1 + sqrt (1 + real x * real x)"
583 let ?fR = "1 + ub_sqrt prec (1 + x * x)"
584 let ?DIV = "float_divl prec x ?fR"
586 have sqr_ge0: "0 \<le> 1 + real x * real x" using sum_power2_ge_zero[of 1 "real x", unfolded numeral_2_eq_2] by auto
587 hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg)
589 have "sqrt (real (1 + x * x)) \<le> real (ub_sqrt prec (1 + x * x))"
590 using bnds_sqrt'[of "1 + x * x"] by auto
592 hence "?R \<le> real ?fR" by auto
593 hence "0 < ?fR" and "0 < real ?fR" unfolding less_float_def using `0 < ?R` by auto
595 have monotone: "real (float_divl prec x ?fR) \<le> real x / ?R"
597 have "real ?DIV \<le> real x / real ?fR" by (rule float_divl)
598 also have "\<dots> \<le> real x / ?R" by (rule divide_left_mono[OF `?R \<le> real ?fR` `0 \<le> real x` mult_pos_pos[OF order_less_le_trans[OF divisor_gt0 `?R \<le> real ?fR`] divisor_gt0]])
599 finally show ?thesis .
603 proof (cases "x \<le> Float 1 1")
606 have "real x \<le> sqrt (real (1 + x * x))" using real_sqrt_sum_squares_ge2[where x=1, unfolded numeral_2_eq_2] by auto
607 also have "\<dots> \<le> real (ub_sqrt prec (1 + x * x))"
608 using bnds_sqrt'[of "1 + x * x"] by auto
609 finally have "real x \<le> real ?fR" by auto
610 moreover have "real ?DIV \<le> real x / real ?fR" by (rule float_divl)
611 ultimately have "real ?DIV \<le> 1" unfolding divide_le_eq_1_pos[OF `0 < real ?fR`, symmetric] by auto
613 have "0 \<le> real ?DIV" using float_divl_lower_bound[OF `0 \<le> x` `0 < ?fR`] unfolding le_float_def by auto
615 have "real (Float 1 1 * ?lb_horner ?DIV) \<le> 2 * arctan (real (float_divl prec x ?fR))" unfolding real_of_float_mult[of "Float 1 1"] Float_num
616 using arctan_0_1_bounds[OF `0 \<le> real ?DIV` `real ?DIV \<le> 1`] by auto
617 also have "\<dots> \<le> 2 * arctan (real x / ?R)"
618 using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono)
619 also have "2 * arctan (real x / ?R) = arctan (real x)" using arctan_half[symmetric] unfolding numeral_2_eq_2 power_Suc2 power_0 real_mult_1 .
620 finally show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_P[OF True] .
623 hence "2 < real x" unfolding le_float_def Float_num by auto
624 hence "1 \<le> real x" by auto
626 let "?invx" = "float_divr prec 1 x"
627 have "0 \<le> arctan (real x)" using arctan_monotone'[OF `0 \<le> real x`] using arctan_tan[of 0, unfolded tan_zero] by auto
630 proof (cases "1 < ?invx")
632 show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_not_P[OF False] if_P[OF True]
633 using `0 \<le> arctan (real x)` by auto
636 hence "real ?invx \<le> 1" unfolding less_float_def by auto
637 have "0 \<le> real ?invx" by (rule order_trans[OF _ float_divr], auto simp add: `0 \<le> real x`)
639 have "1 / real x \<noteq> 0" and "0 < 1 / real x" using `0 < real x` by auto
641 have "arctan (1 / real x) \<le> arctan (real ?invx)" unfolding real_of_float_1[symmetric] by (rule arctan_monotone', rule float_divr)
642 also have "\<dots> \<le> real (?ub_horner ?invx)" using arctan_0_1_bounds[OF `0 \<le> real ?invx` `real ?invx \<le> 1`] by auto
643 finally have "pi / 2 - real (?ub_horner ?invx) \<le> arctan (real x)"
644 using `0 \<le> arctan (real x)` arctan_inverse[OF `1 / real x \<noteq> 0`]
645 unfolding real_sgn_pos[OF `0 < 1 / real x`] le_diff_eq by auto
647 have "real (lb_pi prec * Float 1 -1) \<le> pi / 2" unfolding real_of_float_mult Float_num times_divide_eq_right real_mult_1 using pi_boundaries by auto
649 show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_not_P[OF `\<not> x \<le> Float 1 1`] if_not_P[OF False]
656 lemma ub_arctan_bound': assumes "0 \<le> real x"
657 shows "arctan (real x) \<le> real (ub_arctan prec x)"
659 have "\<not> x < 0" and "0 \<le> x" unfolding less_float_def le_float_def using `0 \<le> real x` by auto
661 let "?ub_horner x" = "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)"
662 and "?lb_horner x" = "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)"
665 proof (cases "x \<le> Float 1 -1")
666 case True hence "real x \<le> 1" unfolding le_float_def Float_num by auto
667 show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True]
668 using arctan_0_1_bounds[OF `0 \<le> real x` `real x \<le> 1`] by auto
670 case False hence "0 < real x" unfolding le_float_def Float_num by auto
671 let ?R = "1 + sqrt (1 + real x * real x)"
672 let ?fR = "1 + lb_sqrt prec (1 + x * x)"
673 let ?DIV = "float_divr prec x ?fR"
675 have sqr_ge0: "0 \<le> 1 + real x * real x" using sum_power2_ge_zero[of 1 "real x", unfolded numeral_2_eq_2] by auto
676 hence "0 \<le> real (1 + x*x)" by auto
678 hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg)
680 have "real (lb_sqrt prec (1 + x * x)) \<le> sqrt (real (1 + x * x))"
681 using bnds_sqrt'[of "1 + x * x"] by auto
682 hence "real ?fR \<le> ?R" by auto
683 have "0 < real ?fR" unfolding real_of_float_add real_of_float_1 by (rule order_less_le_trans[OF zero_less_one], auto simp add: lb_sqrt_lower_bound[OF `0 \<le> real (1 + x*x)`])
685 have monotone: "real x / ?R \<le> real (float_divr prec x ?fR)"
687 from divide_left_mono[OF `real ?fR \<le> ?R` `0 \<le> real x` mult_pos_pos[OF divisor_gt0 `0 < real ?fR`]]
688 have "real x / ?R \<le> real x / real ?fR" .
689 also have "\<dots> \<le> real ?DIV" by (rule float_divr)
690 finally show ?thesis .
694 proof (cases "x \<le> Float 1 1")
697 proof (cases "?DIV > 1")
699 have "pi / 2 \<le> real (ub_pi prec * Float 1 -1)" unfolding real_of_float_mult Float_num times_divide_eq_right real_mult_1 using pi_boundaries by auto
700 from order_less_le_trans[OF arctan_ubound this, THEN less_imp_le]
701 show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_P[OF `x \<le> Float 1 1`] if_P[OF True] .
704 hence "real ?DIV \<le> 1" unfolding less_float_def by auto
706 have "0 \<le> real x / ?R" using `0 \<le> real x` `0 < ?R` unfolding real_0_le_divide_iff by auto
707 hence "0 \<le> real ?DIV" using monotone by (rule order_trans)
709 have "arctan (real x) = 2 * arctan (real x / ?R)" using arctan_half unfolding numeral_2_eq_2 power_Suc2 power_0 real_mult_1 .
710 also have "\<dots> \<le> 2 * arctan (real ?DIV)"
711 using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono)
712 also have "\<dots> \<le> real (Float 1 1 * ?ub_horner ?DIV)" unfolding real_of_float_mult[of "Float 1 1"] Float_num
713 using arctan_0_1_bounds[OF `0 \<le> real ?DIV` `real ?DIV \<le> 1`] by auto
714 finally show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_P[OF `x \<le> Float 1 1`] if_not_P[OF False] .
718 hence "2 < real x" unfolding le_float_def Float_num by auto
719 hence "1 \<le> real x" by auto
720 hence "0 < real x" by auto
721 hence "0 < x" unfolding less_float_def by auto
723 let "?invx" = "float_divl prec 1 x"
724 have "0 \<le> arctan (real x)" using arctan_monotone'[OF `0 \<le> real x`] using arctan_tan[of 0, unfolded tan_zero] by auto
726 have "real ?invx \<le> 1" unfolding less_float_def by (rule order_trans[OF float_divl], auto simp add: `1 \<le> real x` divide_le_eq_1_pos[OF `0 < real x`])
727 have "0 \<le> real ?invx" unfolding real_of_float_0[symmetric] by (rule float_divl_lower_bound[unfolded le_float_def], auto simp add: `0 < x`)
729 have "1 / real x \<noteq> 0" and "0 < 1 / real x" using `0 < real x` by auto
731 have "real (?lb_horner ?invx) \<le> arctan (real ?invx)" using arctan_0_1_bounds[OF `0 \<le> real ?invx` `real ?invx \<le> 1`] by auto
732 also have "\<dots> \<le> arctan (1 / real x)" unfolding real_of_float_1[symmetric] by (rule arctan_monotone', rule float_divl)
733 finally have "arctan (real x) \<le> pi / 2 - real (?lb_horner ?invx)"
734 using `0 \<le> arctan (real x)` arctan_inverse[OF `1 / real x \<noteq> 0`]
735 unfolding real_sgn_pos[OF `0 < 1 / real x`] le_diff_eq by auto
737 have "pi / 2 \<le> real (ub_pi prec * Float 1 -1)" unfolding real_of_float_mult Float_num times_divide_eq_right mult_1_right using pi_boundaries by auto
739 show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_not_P[OF `\<not> x \<le> Float 1 1`] if_not_P[OF False]
745 lemma arctan_boundaries:
746 "arctan (real x) \<in> {real (lb_arctan prec x) .. real (ub_arctan prec x)}"
747 proof (cases "0 \<le> x")
748 case True hence "0 \<le> real x" unfolding le_float_def by auto
749 show ?thesis using ub_arctan_bound'[OF `0 \<le> real x`] lb_arctan_bound'[OF `0 \<le> real x`] unfolding atLeastAtMost_iff by auto
752 case False hence "x < 0" and "0 \<le> real ?mx" unfolding le_float_def less_float_def by auto
753 hence bounds: "real (lb_arctan prec ?mx) \<le> arctan (real ?mx) \<and> arctan (real ?mx) \<le> real (ub_arctan prec ?mx)"
754 using ub_arctan_bound'[OF `0 \<le> real ?mx`] lb_arctan_bound'[OF `0 \<le> real ?mx`] by auto
755 show ?thesis unfolding real_of_float_minus arctan_minus lb_arctan.simps[where x=x] ub_arctan.simps[where x=x] Let_def if_P[OF `x < 0`]
756 unfolding atLeastAtMost_iff using bounds[unfolded real_of_float_minus arctan_minus] by auto
759 lemma bnds_arctan: "\<forall> x lx ux. (l, u) = (lb_arctan prec lx, ub_arctan prec ux) \<and> x \<in> {real lx .. real ux} \<longrightarrow> real l \<le> arctan x \<and> arctan x \<le> real u"
760 proof (rule allI, rule allI, rule allI, rule impI)
762 assume "(l, u) = (lb_arctan prec lx, ub_arctan prec ux) \<and> x \<in> {real lx .. real ux}"
763 hence l: "lb_arctan prec lx = l " and u: "ub_arctan prec ux = u" and x: "x \<in> {real lx .. real ux}" by auto
765 { from arctan_boundaries[of lx prec, unfolded l]
766 have "real l \<le> arctan (real lx)" by (auto simp del: lb_arctan.simps)
767 also have "\<dots> \<le> arctan x" using x by (auto intro: arctan_monotone')
768 finally have "real l \<le> arctan x" .
770 { have "arctan x \<le> arctan (real ux)" using x by (auto intro: arctan_monotone')
771 also have "\<dots> \<le> real u" using arctan_boundaries[of ux prec, unfolded u] by (auto simp del: ub_arctan.simps)
772 finally have "arctan x \<le> real u" .
773 } ultimately show "real l \<le> arctan x \<and> arctan x \<le> real u" ..
776 section "Sinus and Cosinus"
778 subsection "Compute the cosinus and sinus series"
780 fun ub_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
781 and lb_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
782 "ub_sin_cos_aux prec 0 i k x = 0"
783 | "ub_sin_cos_aux prec (Suc n) i k x =
784 (rapprox_rat prec 1 (int k)) - x * (lb_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)"
785 | "lb_sin_cos_aux prec 0 i k x = 0"
786 | "lb_sin_cos_aux prec (Suc n) i k x =
787 (lapprox_rat prec 1 (int k)) - x * (ub_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)"
790 shows "real (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> (\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i))) * (real x)^(2 * i))" (is "?lb")
791 and "(\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i))) * (real x)^(2 * i)) \<le> real (ub_sin_cos_aux prec n 1 1 (x * x))" (is "?ub")
793 have "0 \<le> real (x * x)" unfolding real_of_float_mult by auto
794 let "?f n" = "fact (2 * n)"
797 have F: "\<And>m. ((\<lambda>i. i + 2) ^^ n) m = m + 2 * n" by (induct n arbitrary: m, auto)
798 have "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) ^^ n) 1 * (((\<lambda>i. i + 2) ^^ n) 1 + 1)"
799 unfolding F by auto } note f_eq = this
801 from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0,
802 OF `0 \<le> real (x * x)` f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps]
803 show "?lb" and "?ub" by (auto simp add: power_mult power2_eq_square[of "real x"])
806 lemma cos_boundaries: assumes "0 \<le> real x" and "real x \<le> pi / 2"
807 shows "cos (real x) \<in> {real (lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) .. real (ub_sin_cos_aux prec (get_odd n) 1 1 (x * x))}"
808 proof (cases "real x = 0")
809 case False hence "real x \<noteq> 0" by auto
810 hence "0 < x" and "0 < real x" using `0 \<le> real x` unfolding less_float_def by auto
811 have "0 < x * x" using `0 < x` unfolding less_float_def real_of_float_mult real_of_float_0
812 using mult_pos_pos[where a="real x" and b="real x"] by auto
814 { fix x n have "(\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i))) * x ^ (2 * i))
815 = (\<Sum> i = 0 ..< 2 * n. (if even(i) then (-1 ^ (i div 2))/(real (fact i)) else 0) * x ^ i)" (is "?sum = ?ifsum")
817 have "?sum = ?sum + (\<Sum> j = 0 ..< n. 0)" by auto
819 (\<Sum> j = 0 ..< n. -1 ^ ((2 * j) div 2) / (real (fact (2 * j))) * x ^(2 * j)) + (\<Sum> j = 0 ..< n. 0)" by auto
820 also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. if even i then -1 ^ (i div 2) / (real (fact i)) * x ^ i else 0)"
821 unfolding sum_split_even_odd ..
822 also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. (if even i then -1 ^ (i div 2) / (real (fact i)) else 0) * x ^ i)"
823 by (rule setsum_cong2) auto
824 finally show ?thesis by assumption
825 qed } note morph_to_if_power = this
828 { fix n :: nat assume "0 < n"
829 hence "0 < 2 * n" by auto
830 obtain t where "0 < t" and "t < real x" and
831 cos_eq: "cos (real x) = (\<Sum> i = 0 ..< 2 * n. (if even(i) then (-1 ^ (i div 2))/(real (fact i)) else 0) * (real x) ^ i)
832 + (cos (t + 1/2 * real (2 * n) * pi) / real (fact (2*n))) * (real x)^(2*n)"
833 (is "_ = ?SUM + ?rest / ?fact * ?pow")
834 using Maclaurin_cos_expansion2[OF `0 < real x` `0 < 2 * n`] by auto
836 have "cos t * -1^n = cos t * cos (real n * pi) + sin t * sin (real n * pi)" by auto
837 also have "\<dots> = cos (t + real n * pi)" using cos_add by auto
838 also have "\<dots> = ?rest" by auto
839 finally have "cos t * -1^n = ?rest" .
841 have "t \<le> pi / 2" using `t < real x` and `real x \<le> pi / 2` by auto
842 hence "0 \<le> cos t" using `0 < t` and cos_ge_zero by auto
843 ultimately have even: "even n \<Longrightarrow> 0 \<le> ?rest" and odd: "odd n \<Longrightarrow> 0 \<le> - ?rest " by auto
845 have "0 < ?fact" by auto
846 have "0 < ?pow" using `0 < real x` by auto
850 have "real (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> ?SUM"
851 unfolding morph_to_if_power[symmetric] using cos_aux by auto
852 also have "\<dots> \<le> cos (real x)"
854 from even[OF `even n`] `0 < ?fact` `0 < ?pow`
855 have "0 \<le> (?rest / ?fact) * ?pow" by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
856 thus ?thesis unfolding cos_eq by auto
858 finally have "real (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> cos (real x)" .
863 have "cos (real x) \<le> ?SUM"
865 from `0 < ?fact` and `0 < ?pow` and odd[OF `odd n`]
866 have "0 \<le> (- ?rest) / ?fact * ?pow"
867 by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
868 thus ?thesis unfolding cos_eq by auto
870 also have "\<dots> \<le> real (ub_sin_cos_aux prec n 1 1 (x * x))"
871 unfolding morph_to_if_power[symmetric] using cos_aux by auto
872 finally have "cos (real x) \<le> real (ub_sin_cos_aux prec n 1 1 (x * x))" .
873 } note ub = this and lb
874 } note ub = this(1) and lb = this(2)
876 have "cos (real x) \<le> real (ub_sin_cos_aux prec (get_odd n) 1 1 (x * x))" using ub[OF odd_pos[OF get_odd] get_odd] .
877 moreover have "real (lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) \<le> cos (real x)"
878 proof (cases "0 < get_even n")
879 case True show ?thesis using lb[OF True get_even] .
882 hence "get_even n = 0" by auto
883 have "- (pi / 2) \<le> real x" by (rule order_trans[OF _ `0 < real x`[THEN less_imp_le]], auto)
884 with `real x \<le> pi / 2`
885 show ?thesis unfolding `get_even n = 0` lb_sin_cos_aux.simps real_of_float_minus real_of_float_0 using cos_ge_zero by auto
887 ultimately show ?thesis by auto
891 proof (cases "n = 0")
893 thus ?thesis unfolding `n = 0` get_even_def get_odd_def using `real x = 0` lapprox_rat[where x="-1" and y=1] by auto
895 case False with not0_implies_Suc obtain m where "n = Suc m" by blast
896 thus ?thesis unfolding `n = Suc m` get_even_def get_odd_def using `real x = 0` rapprox_rat[where x=1 and y=1] lapprox_rat[where x=1 and y=1] by (cases "even (Suc m)", auto)
900 lemma sin_aux: assumes "0 \<le> real x"
901 shows "real (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> (\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i + 1))) * (real x)^(2 * i + 1))" (is "?lb")
902 and "(\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i + 1))) * (real x)^(2 * i + 1)) \<le> real (x * ub_sin_cos_aux prec n 2 1 (x * x))" (is "?ub")
904 have "0 \<le> real (x * x)" unfolding real_of_float_mult by auto
905 let "?f n" = "fact (2 * n + 1)"
908 have F: "\<And>m. ((\<lambda>i. i + 2) ^^ n) m = m + 2 * n" by (induct n arbitrary: m, auto)
909 have "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) ^^ n) 2 * (((\<lambda>i. i + 2) ^^ n) 2 + 1)"
910 unfolding F by auto } note f_eq = this
912 from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0,
913 OF `0 \<le> real (x * x)` f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps]
914 show "?lb" and "?ub" using `0 \<le> real x` unfolding real_of_float_mult
915 unfolding power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]
916 unfolding real_mult_commute
917 by (auto intro!: mult_left_mono simp add: power_mult power2_eq_square[of "real x"])
920 lemma sin_boundaries: assumes "0 \<le> real x" and "real x \<le> pi / 2"
921 shows "sin (real x) \<in> {real (x * lb_sin_cos_aux prec (get_even n) 2 1 (x * x)) .. real (x * ub_sin_cos_aux prec (get_odd n) 2 1 (x * x))}"
922 proof (cases "real x = 0")
923 case False hence "real x \<noteq> 0" by auto
924 hence "0 < x" and "0 < real x" using `0 \<le> real x` unfolding less_float_def by auto
925 have "0 < x * x" using `0 < x` unfolding less_float_def real_of_float_mult real_of_float_0
926 using mult_pos_pos[where a="real x" and b="real x"] by auto
928 { fix x n have "(\<Sum> j = 0 ..< n. -1 ^ (((2 * j + 1) - Suc 0) div 2) / (real (fact (2 * j + 1))) * x ^(2 * j + 1))
929 = (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * x ^ i)" (is "?SUM = _")
931 have pow: "!!i. x ^ (2 * i + 1) = x * x ^ (2 * i)" by auto
932 have "?SUM = (\<Sum> j = 0 ..< n. 0) + ?SUM" by auto
933 also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. if even i then 0 else -1 ^ ((i - Suc 0) div 2) / (real (fact i)) * x ^ i)"
934 unfolding sum_split_even_odd ..
935 also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. (if even i then 0 else -1 ^ ((i - Suc 0) div 2) / (real (fact i))) * x ^ i)"
936 by (rule setsum_cong2) auto
937 finally show ?thesis by assumption
938 qed } note setsum_morph = this
940 { fix n :: nat assume "0 < n"
941 hence "0 < 2 * n + 1" by auto
942 obtain t where "0 < t" and "t < real x" and
943 sin_eq: "sin (real x) = (\<Sum> i = 0 ..< 2 * n + 1. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * (real x) ^ i)
944 + (sin (t + 1/2 * real (2 * n + 1) * pi) / real (fact (2*n + 1))) * (real x)^(2*n + 1)"
945 (is "_ = ?SUM + ?rest / ?fact * ?pow")
946 using Maclaurin_sin_expansion3[OF `0 < 2 * n + 1` `0 < real x`] by auto
948 have "?rest = cos t * -1^n" unfolding sin_add cos_add real_of_nat_add left_distrib right_distrib by auto
950 have "t \<le> pi / 2" using `t < real x` and `real x \<le> pi / 2` by auto
951 hence "0 \<le> cos t" using `0 < t` and cos_ge_zero by auto
952 ultimately have even: "even n \<Longrightarrow> 0 \<le> ?rest" and odd: "odd n \<Longrightarrow> 0 \<le> - ?rest " by auto
954 have "0 < ?fact" by (rule real_of_nat_fact_gt_zero)
955 have "0 < ?pow" using `0 < real x` by (rule zero_less_power)
959 have "real (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le>
960 (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * (real x) ^ i)"
961 using sin_aux[OF `0 \<le> real x`] unfolding setsum_morph[symmetric] by auto
962 also have "\<dots> \<le> ?SUM" by auto
963 also have "\<dots> \<le> sin (real x)"
965 from even[OF `even n`] `0 < ?fact` `0 < ?pow`
966 have "0 \<le> (?rest / ?fact) * ?pow" by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
967 thus ?thesis unfolding sin_eq by auto
969 finally have "real (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> sin (real x)" .
974 have "sin (real x) \<le> ?SUM"
976 from `0 < ?fact` and `0 < ?pow` and odd[OF `odd n`]
977 have "0 \<le> (- ?rest) / ?fact * ?pow"
978 by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
979 thus ?thesis unfolding sin_eq by auto
981 also have "\<dots> \<le> (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * (real x) ^ i)"
983 also have "\<dots> \<le> real (x * ub_sin_cos_aux prec n 2 1 (x * x))"
984 using sin_aux[OF `0 \<le> real x`] unfolding setsum_morph[symmetric] by auto
985 finally have "sin (real x) \<le> real (x * ub_sin_cos_aux prec n 2 1 (x * x))" .
986 } note ub = this and lb
987 } note ub = this(1) and lb = this(2)
989 have "sin (real x) \<le> real (x * ub_sin_cos_aux prec (get_odd n) 2 1 (x * x))" using ub[OF odd_pos[OF get_odd] get_odd] .
990 moreover have "real (x * lb_sin_cos_aux prec (get_even n) 2 1 (x * x)) \<le> sin (real x)"
991 proof (cases "0 < get_even n")
992 case True show ?thesis using lb[OF True get_even] .
995 hence "get_even n = 0" by auto
996 with `real x \<le> pi / 2` `0 \<le> real x`
997 show ?thesis unfolding `get_even n = 0` ub_sin_cos_aux.simps real_of_float_minus real_of_float_0 using sin_ge_zero by auto
999 ultimately show ?thesis by auto
1003 proof (cases "n = 0")
1005 thus ?thesis unfolding `n = 0` get_even_def get_odd_def using `real x = 0` lapprox_rat[where x="-1" and y=1] by auto
1007 case False with not0_implies_Suc obtain m where "n = Suc m" by blast
1008 thus ?thesis unfolding `n = Suc m` get_even_def get_odd_def using `real x = 0` rapprox_rat[where x=1 and y=1] lapprox_rat[where x=1 and y=1] by (cases "even (Suc m)", auto)
1012 subsection "Compute the cosinus in the entire domain"
1014 definition lb_cos :: "nat \<Rightarrow> float \<Rightarrow> float" where
1015 "lb_cos prec x = (let
1016 horner = \<lambda> x. lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 1 1 (x * x) ;
1017 half = \<lambda> x. if x < 0 then - 1 else Float 1 1 * x * x - 1
1018 in if x < Float 1 -1 then horner x
1019 else if x < 1 then half (horner (x * Float 1 -1))
1020 else half (half (horner (x * Float 1 -2))))"
1022 definition ub_cos :: "nat \<Rightarrow> float \<Rightarrow> float" where
1023 "ub_cos prec x = (let
1024 horner = \<lambda> x. ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x) ;
1025 half = \<lambda> x. Float 1 1 * x * x - 1
1026 in if x < Float 1 -1 then horner x
1027 else if x < 1 then half (horner (x * Float 1 -1))
1028 else half (half (horner (x * Float 1 -2))))"
1030 lemma lb_cos: assumes "0 \<le> real x" and "real x \<le> pi"
1031 shows "cos (real x) \<in> {real (lb_cos prec x) .. real (ub_cos prec x)}" (is "?cos x \<in> { real (?lb x) .. real (?ub x) }")
1034 have "cos x = cos (x / 2 + x / 2)" by auto
1035 also have "\<dots> = cos (x / 2) * cos (x / 2) + sin (x / 2) * sin (x / 2) - sin (x / 2) * sin (x / 2) + cos (x / 2) * cos (x / 2) - 1"
1036 unfolding cos_add by auto
1037 also have "\<dots> = 2 * cos (x / 2) * cos (x / 2) - 1" by algebra
1038 finally have "cos x = 2 * cos (x / 2) * cos (x / 2) - 1" .
1039 } note x_half = this[symmetric]
1041 have "\<not> x < 0" using `0 \<le> real x` unfolding less_float_def by auto
1042 let "?ub_horner x" = "ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x)"
1043 let "?lb_horner x" = "lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 1 1 (x * x)"
1044 let "?ub_half x" = "Float 1 1 * x * x - 1"
1045 let "?lb_half x" = "if x < 0 then - 1 else Float 1 1 * x * x - 1"
1048 proof (cases "x < Float 1 -1")
1049 case True hence "real x \<le> pi / 2" unfolding less_float_def using pi_ge_two by auto
1050 show ?thesis unfolding lb_cos_def[where x=x] ub_cos_def[where x=x] if_not_P[OF `\<not> x < 0`] if_P[OF `x < Float 1 -1`] Let_def
1051 using cos_boundaries[OF `0 \<le> real x` `real x \<le> pi / 2`] .
1054 { fix y x :: float let ?x2 = "real (x * Float 1 -1)"
1055 assume "real y \<le> cos ?x2" and "-pi \<le> real x" and "real x \<le> pi"
1056 hence "- (pi / 2) \<le> ?x2" and "?x2 \<le> pi / 2" using pi_ge_two unfolding real_of_float_mult Float_num by auto
1057 hence "0 \<le> cos ?x2" by (rule cos_ge_zero)
1059 have "real (?lb_half y) \<le> cos (real x)"
1060 proof (cases "y < 0")
1061 case True show ?thesis using cos_ge_minus_one unfolding if_P[OF True] by auto
1064 hence "0 \<le> real y" unfolding less_float_def by auto
1065 from mult_mono[OF `real y \<le> cos ?x2` `real y \<le> cos ?x2` `0 \<le> cos ?x2` this]
1066 have "real y * real y \<le> cos ?x2 * cos ?x2" .
1067 hence "2 * real y * real y \<le> 2 * cos ?x2 * cos ?x2" by auto
1068 hence "2 * real y * real y - 1 \<le> 2 * cos (real x / 2) * cos (real x / 2) - 1" unfolding Float_num real_of_float_mult by auto
1069 thus ?thesis unfolding if_not_P[OF False] x_half Float_num real_of_float_mult real_of_float_sub by auto
1071 } note lb_half = this
1073 { fix y x :: float let ?x2 = "real (x * Float 1 -1)"
1074 assume ub: "cos ?x2 \<le> real y" and "- pi \<le> real x" and "real x \<le> pi"
1075 hence "- (pi / 2) \<le> ?x2" and "?x2 \<le> pi / 2" using pi_ge_two unfolding real_of_float_mult Float_num by auto
1076 hence "0 \<le> cos ?x2" by (rule cos_ge_zero)
1078 have "cos (real x) \<le> real (?ub_half y)"
1080 have "0 \<le> real y" using `0 \<le> cos ?x2` ub by (rule order_trans)
1081 from mult_mono[OF ub ub this `0 \<le> cos ?x2`]
1082 have "cos ?x2 * cos ?x2 \<le> real y * real y" .
1083 hence "2 * cos ?x2 * cos ?x2 \<le> 2 * real y * real y" by auto
1084 hence "2 * cos (real x / 2) * cos (real x / 2) - 1 \<le> 2 * real y * real y - 1" unfolding Float_num real_of_float_mult by auto
1085 thus ?thesis unfolding x_half real_of_float_mult Float_num real_of_float_sub by auto
1087 } note ub_half = this
1089 let ?x2 = "x * Float 1 -1"
1090 let ?x4 = "x * Float 1 -1 * Float 1 -1"
1092 have "-pi \<le> real x" using pi_ge_zero[THEN le_imp_neg_le, unfolded minus_zero] `0 \<le> real x` by (rule order_trans)
1095 proof (cases "x < 1")
1096 case True hence "real x \<le> 1" unfolding less_float_def by auto
1097 have "0 \<le> real ?x2" and "real ?x2 \<le> pi / 2" using pi_ge_two `0 \<le> real x` unfolding real_of_float_mult Float_num using assms by auto
1098 from cos_boundaries[OF this]
1099 have lb: "real (?lb_horner ?x2) \<le> ?cos ?x2" and ub: "?cos ?x2 \<le> real (?ub_horner ?x2)" by auto
1101 have "real (?lb x) \<le> ?cos x"
1103 from lb_half[OF lb `-pi \<le> real x` `real x \<le> pi`]
1104 show ?thesis unfolding lb_cos_def[where x=x] Let_def using `\<not> x < 0` `\<not> x < Float 1 -1` `x < 1` by auto
1106 moreover have "?cos x \<le> real (?ub x)"
1108 from ub_half[OF ub `-pi \<le> real x` `real x \<le> pi`]
1109 show ?thesis unfolding ub_cos_def[where x=x] Let_def using `\<not> x < 0` `\<not> x < Float 1 -1` `x < 1` by auto
1111 ultimately show ?thesis by auto
1114 have "0 \<le> real ?x4" and "real ?x4 \<le> pi / 2" using pi_ge_two `0 \<le> real x` `real x \<le> pi` unfolding real_of_float_mult Float_num by auto
1115 from cos_boundaries[OF this]
1116 have lb: "real (?lb_horner ?x4) \<le> ?cos ?x4" and ub: "?cos ?x4 \<le> real (?ub_horner ?x4)" by auto
1118 have eq_4: "?x2 * Float 1 -1 = x * Float 1 -2" by (cases x, auto simp add: times_float.simps)
1120 have "real (?lb x) \<le> ?cos x"
1122 have "-pi \<le> real ?x2" and "real ?x2 \<le> pi" unfolding real_of_float_mult Float_num using pi_ge_two `0 \<le> real x` `real x \<le> pi` by auto
1123 from lb_half[OF lb_half[OF lb this] `-pi \<le> real x` `real x \<le> pi`, unfolded eq_4]
1124 show ?thesis unfolding lb_cos_def[where x=x] if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x < Float 1 -1`] if_not_P[OF `\<not> x < 1`] Let_def .
1126 moreover have "?cos x \<le> real (?ub x)"
1128 have "-pi \<le> real ?x2" and "real ?x2 \<le> pi" unfolding real_of_float_mult Float_num using pi_ge_two `0 \<le> real x` `real x \<le> pi` by auto
1129 from ub_half[OF ub_half[OF ub this] `-pi \<le> real x` `real x \<le> pi`, unfolded eq_4]
1130 show ?thesis unfolding ub_cos_def[where x=x] if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x < Float 1 -1`] if_not_P[OF `\<not> x < 1`] Let_def .
1132 ultimately show ?thesis by auto
1137 lemma lb_cos_minus: assumes "-pi \<le> real x" and "real x \<le> 0"
1138 shows "cos (real (-x)) \<in> {real (lb_cos prec (-x)) .. real (ub_cos prec (-x))}"
1140 have "0 \<le> real (-x)" and "real (-x) \<le> pi" using `-pi \<le> real x` `real x \<le> 0` by auto
1141 from lb_cos[OF this] show ?thesis .
1144 definition bnds_cos :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where
1145 "bnds_cos prec lx ux = (let
1146 lpi = round_down prec (lb_pi prec) ;
1147 upi = round_up prec (ub_pi prec) ;
1148 k = floor_fl (float_divr prec (lx + lpi) (2 * lpi)) ;
1149 lx = lx - k * 2 * (if k < 0 then lpi else upi) ;
1150 ux = ux - k * 2 * (if k < 0 then upi else lpi)
1151 in if - lpi \<le> lx \<and> ux \<le> 0 then (lb_cos prec (-lx), ub_cos prec (-ux))
1152 else if 0 \<le> lx \<and> ux \<le> lpi then (lb_cos prec ux, ub_cos prec lx)
1153 else if - lpi \<le> lx \<and> ux \<le> lpi then (min (lb_cos prec (-lx)) (lb_cos prec ux), Float 1 0)
1154 else if 0 \<le> lx \<and> ux \<le> 2 * lpi then (Float -1 0, max (ub_cos prec lx) (ub_cos prec (- (ux - 2 * lpi))))
1155 else if -2 * lpi \<le> lx \<and> ux \<le> 0 then (Float -1 0, max (ub_cos prec (lx + 2 * lpi)) (ub_cos prec (-ux)))
1156 else (Float -1 0, Float 1 0))"
1159 obtains k :: int where "real k = real (floor_fl f)"
1161 assume *: "\<And> k :: int. real k = real (floor_fl f) \<Longrightarrow> thesis"
1162 obtain m e where fl: "Float m e = floor_fl f" by (cases "floor_fl f", auto)
1163 from floor_pos_exp[OF this]
1164 have "real (m* 2^(nat e)) = real (floor_fl f)"
1165 by (auto simp add: fl[symmetric] real_of_float_def pow2_def)
1166 from *[OF this] show thesis by blast
1169 lemma float_remove_real_numeral[simp]: "real (number_of k :: float) = number_of k"
1171 have "real (number_of k :: float) = real k"
1172 unfolding number_of_float_def real_of_float_def pow2_def by auto
1173 also have "\<dots> = real (number_of k :: int)"
1174 by (simp add: number_of_is_id)
1175 finally show ?thesis by auto
1178 lemma cos_periodic_nat[simp]: fixes n :: nat shows "cos (x + real n * 2 * pi) = cos x"
1179 proof (induct n arbitrary: x)
1181 have split_pi_off: "x + real (Suc n) * 2 * pi = (x + real n * 2 * pi) + 2 * pi"
1182 unfolding Suc_eq_plus1 real_of_nat_add real_of_one real_add_mult_distrib by auto
1183 show ?case unfolding split_pi_off using Suc by auto
1186 lemma cos_periodic_int[simp]: fixes i :: int shows "cos (x + real i * 2 * pi) = cos x"
1187 proof (cases "0 \<le> i")
1188 case True hence i_nat: "real i = real (nat i)" by auto
1189 show ?thesis unfolding i_nat by auto
1191 case False hence i_nat: "real i = - real (nat (-i))" by auto
1192 have "cos x = cos (x + real i * 2 * pi - real i * 2 * pi)" by auto
1193 also have "\<dots> = cos (x + real i * 2 * pi)"
1194 unfolding i_nat mult_minus_left diff_minus_eq_add by (rule cos_periodic_nat)
1195 finally show ?thesis by auto
1198 lemma bnds_cos: "\<forall> x lx ux. (l, u) = bnds_cos prec lx ux \<and> x \<in> {real lx .. real ux} \<longrightarrow> real l \<le> cos x \<and> cos x \<le> real u"
1199 proof ((rule allI | rule impI | erule conjE) +)
1201 assume bnds: "(l, u) = bnds_cos prec lx ux" and x: "x \<in> {real lx .. real ux}"
1203 let ?lpi = "round_down prec (lb_pi prec)"
1204 let ?upi = "round_up prec (ub_pi prec)"
1205 let ?k = "floor_fl (float_divr prec (lx + ?lpi) (2 * ?lpi))"
1206 let ?lx = "lx - ?k * 2 * (if ?k < 0 then ?lpi else ?upi)"
1207 let ?ux = "ux - ?k * 2 * (if ?k < 0 then ?upi else ?lpi)"
1209 obtain k :: int where k: "real k = real ?k" using floor_int .
1211 have upi: "pi \<le> real ?upi" and lpi: "real ?lpi \<le> pi"
1212 using round_up[of "ub_pi prec" prec] pi_boundaries[of prec]
1213 round_down[of prec "lb_pi prec"] by auto
1214 hence "real ?lx \<le> x - real k * 2 * pi \<and> x - real k * 2 * pi \<le> real ?ux"
1215 using x by (cases "k = 0") (auto intro!: add_mono
1216 simp add: real_diff_def k[symmetric] less_float_def)
1217 note lx = this[THEN conjunct1] and ux = this[THEN conjunct2]
1218 hence lx_less_ux: "real ?lx \<le> real ?ux" by (rule order_trans)
1220 { assume "- ?lpi \<le> ?lx" and x_le_0: "x - real k * 2 * pi \<le> 0"
1221 with lpi[THEN le_imp_neg_le] lx
1222 have pi_lx: "- pi \<le> real ?lx" and lx_0: "real ?lx \<le> 0"
1223 by (simp_all add: le_float_def)
1225 have "real (lb_cos prec (- ?lx)) \<le> cos (real (- ?lx))"
1226 using lb_cos_minus[OF pi_lx lx_0] by simp
1227 also have "\<dots> \<le> cos (x + real (-k) * 2 * pi)"
1228 using cos_monotone_minus_pi_0'[OF pi_lx lx x_le_0]
1229 by (simp only: real_of_float_minus real_of_int_minus
1230 cos_minus real_diff_def mult_minus_left)
1231 finally have "real (lb_cos prec (- ?lx)) \<le> cos x"
1232 unfolding cos_periodic_int . }
1233 note negative_lx = this
1235 { assume "0 \<le> ?lx" and pi_x: "x - real k * 2 * pi \<le> pi"
1237 have pi_lx: "real ?lx \<le> pi" and lx_0: "0 \<le> real ?lx"
1238 by (auto simp add: le_float_def)
1240 have "cos (x + real (-k) * 2 * pi) \<le> cos (real ?lx)"
1241 using cos_monotone_0_pi'[OF lx_0 lx pi_x]
1242 by (simp only: real_of_float_minus real_of_int_minus
1243 cos_minus real_diff_def mult_minus_left)
1244 also have "\<dots> \<le> real (ub_cos prec ?lx)"
1245 using lb_cos[OF lx_0 pi_lx] by simp
1246 finally have "cos x \<le> real (ub_cos prec ?lx)"
1247 unfolding cos_periodic_int . }
1248 note positive_lx = this
1250 { assume pi_x: "- pi \<le> x - real k * 2 * pi" and "?ux \<le> 0"
1252 have pi_ux: "- pi \<le> real ?ux" and ux_0: "real ?ux \<le> 0"
1253 by (simp_all add: le_float_def)
1255 have "cos (x + real (-k) * 2 * pi) \<le> cos (real (- ?ux))"
1256 using cos_monotone_minus_pi_0'[OF pi_x ux ux_0]
1257 by (simp only: real_of_float_minus real_of_int_minus
1258 cos_minus real_diff_def mult_minus_left)
1259 also have "\<dots> \<le> real (ub_cos prec (- ?ux))"
1260 using lb_cos_minus[OF pi_ux ux_0, of prec] by simp
1261 finally have "cos x \<le> real (ub_cos prec (- ?ux))"
1262 unfolding cos_periodic_int . }
1263 note negative_ux = this
1265 { assume "?ux \<le> ?lpi" and x_ge_0: "0 \<le> x - real k * 2 * pi"
1267 have pi_ux: "real ?ux \<le> pi" and ux_0: "0 \<le> real ?ux"
1268 by (simp_all add: le_float_def)
1270 have "real (lb_cos prec ?ux) \<le> cos (real ?ux)"
1271 using lb_cos[OF ux_0 pi_ux] by simp
1272 also have "\<dots> \<le> cos (x + real (-k) * 2 * pi)"
1273 using cos_monotone_0_pi'[OF x_ge_0 ux pi_ux]
1274 by (simp only: real_of_float_minus real_of_int_minus
1275 cos_minus real_diff_def mult_minus_left)
1276 finally have "real (lb_cos prec ?ux) \<le> cos x"
1277 unfolding cos_periodic_int . }
1278 note positive_ux = this
1280 show "real l \<le> cos x \<and> cos x \<le> real u"
1281 proof (cases "- ?lpi \<le> ?lx \<and> ?ux \<le> 0")
1283 have l: "l = lb_cos prec (-?lx)"
1284 and u: "u = ub_cos prec (-?ux)"
1285 by (auto simp add: bnds_cos_def Let_def)
1287 from True lpi[THEN le_imp_neg_le] lx ux
1288 have "- pi \<le> x - real k * 2 * pi"
1289 and "x - real k * 2 * pi \<le> 0"
1290 by (auto simp add: le_float_def)
1291 with True negative_ux negative_lx
1292 show ?thesis unfolding l u by simp
1293 next case False note 1 = this show ?thesis
1294 proof (cases "0 \<le> ?lx \<and> ?ux \<le> ?lpi")
1295 case True with bnds 1
1296 have l: "l = lb_cos prec ?ux"
1297 and u: "u = ub_cos prec ?lx"
1298 by (auto simp add: bnds_cos_def Let_def)
1301 have "0 \<le> x - real k * 2 * pi"
1302 and "x - real k * 2 * pi \<le> pi"
1303 by (auto simp add: le_float_def)
1304 with True positive_ux positive_lx
1305 show ?thesis unfolding l u by simp
1306 next case False note 2 = this show ?thesis
1307 proof (cases "- ?lpi \<le> ?lx \<and> ?ux \<le> ?lpi")
1308 case True note Cond = this with bnds 1 2
1309 have l: "l = min (lb_cos prec (-?lx)) (lb_cos prec ?ux)"
1310 and u: "u = Float 1 0"
1311 by (auto simp add: bnds_cos_def Let_def)
1313 show ?thesis unfolding u l using negative_lx positive_ux Cond
1314 by (cases "x - real k * 2 * pi < 0", simp_all add: real_of_float_min)
1315 next case False note 3 = this show ?thesis
1316 proof (cases "0 \<le> ?lx \<and> ?ux \<le> 2 * ?lpi")
1317 case True note Cond = this with bnds 1 2 3
1318 have l: "l = Float -1 0"
1319 and u: "u = max (ub_cos prec ?lx) (ub_cos prec (- (?ux - 2 * ?lpi)))"
1320 by (auto simp add: bnds_cos_def Let_def)
1322 have "cos x \<le> real u"
1323 proof (cases "x - real k * 2 * pi < pi")
1324 case True hence "x - real k * 2 * pi \<le> pi" by simp
1325 from positive_lx[OF Cond[THEN conjunct1] this]
1326 show ?thesis unfolding u by (simp add: real_of_float_max)
1328 case False hence "pi \<le> x - real k * 2 * pi" by simp
1329 hence pi_x: "- pi \<le> x - real k * 2 * pi - 2 * pi" by simp
1331 have "real ?ux \<le> 2 * pi" using Cond lpi by (auto simp add: le_float_def)
1332 hence "x - real k * 2 * pi - 2 * pi \<le> 0" using ux by simp
1334 have ux_0: "real (?ux - 2 * ?lpi) \<le> 0"
1335 using Cond by (auto simp add: le_float_def)
1337 from 2 and Cond have "\<not> ?ux \<le> ?lpi" by auto
1338 hence "- ?lpi \<le> ?ux - 2 * ?lpi" by (auto simp add: le_float_def)
1339 hence pi_ux: "- pi \<le> real (?ux - 2 * ?lpi)"
1340 using lpi[THEN le_imp_neg_le] by (auto simp add: le_float_def)
1342 have x_le_ux: "x - real k * 2 * pi - 2 * pi \<le> real (?ux - 2 * ?lpi)"
1343 using ux lpi by auto
1345 have "cos x = cos (x + real (-k) * 2 * pi + real (-1 :: int) * 2 * pi)"
1346 unfolding cos_periodic_int ..
1347 also have "\<dots> \<le> cos (real (?ux - 2 * ?lpi))"
1348 using cos_monotone_minus_pi_0'[OF pi_x x_le_ux ux_0]
1349 by (simp only: real_of_float_minus real_of_int_minus real_of_one
1350 number_of_Min real_diff_def mult_minus_left real_mult_1)
1351 also have "\<dots> = cos (real (- (?ux - 2 * ?lpi)))"
1352 unfolding real_of_float_minus cos_minus ..
1353 also have "\<dots> \<le> real (ub_cos prec (- (?ux - 2 * ?lpi)))"
1354 using lb_cos_minus[OF pi_ux ux_0] by simp
1355 finally show ?thesis unfolding u by (simp add: real_of_float_max)
1357 thus ?thesis unfolding l by auto
1358 next case False note 4 = this show ?thesis
1359 proof (cases "-2 * ?lpi \<le> ?lx \<and> ?ux \<le> 0")
1360 case True note Cond = this with bnds 1 2 3 4
1361 have l: "l = Float -1 0"
1362 and u: "u = max (ub_cos prec (?lx + 2 * ?lpi)) (ub_cos prec (-?ux))"
1363 by (auto simp add: bnds_cos_def Let_def)
1365 have "cos x \<le> real u"
1366 proof (cases "-pi < x - real k * 2 * pi")
1367 case True hence "-pi \<le> x - real k * 2 * pi" by simp
1368 from negative_ux[OF this Cond[THEN conjunct2]]
1369 show ?thesis unfolding u by (simp add: real_of_float_max)
1371 case False hence "x - real k * 2 * pi \<le> -pi" by simp
1372 hence pi_x: "x - real k * 2 * pi + 2 * pi \<le> pi" by simp
1374 have "-2 * pi \<le> real ?lx" using Cond lpi by (auto simp add: le_float_def)
1376 hence "0 \<le> x - real k * 2 * pi + 2 * pi" using lx by simp
1378 have lx_0: "0 \<le> real (?lx + 2 * ?lpi)"
1379 using Cond lpi by (auto simp add: le_float_def)
1381 from 1 and Cond have "\<not> -?lpi \<le> ?lx" by auto
1382 hence "?lx + 2 * ?lpi \<le> ?lpi" by (auto simp add: le_float_def)
1383 hence pi_lx: "real (?lx + 2 * ?lpi) \<le> pi"
1384 using lpi[THEN le_imp_neg_le] by (auto simp add: le_float_def)
1386 have lx_le_x: "real (?lx + 2 * ?lpi) \<le> x - real k * 2 * pi + 2 * pi"
1387 using lx lpi by auto
1389 have "cos x = cos (x + real (-k) * 2 * pi + real (1 :: int) * 2 * pi)"
1390 unfolding cos_periodic_int ..
1391 also have "\<dots> \<le> cos (real (?lx + 2 * ?lpi))"
1392 using cos_monotone_0_pi'[OF lx_0 lx_le_x pi_x]
1393 by (simp only: real_of_float_minus real_of_int_minus real_of_one
1394 number_of_Min real_diff_def mult_minus_left real_mult_1)
1395 also have "\<dots> \<le> real (ub_cos prec (?lx + 2 * ?lpi))"
1396 using lb_cos[OF lx_0 pi_lx] by simp
1397 finally show ?thesis unfolding u by (simp add: real_of_float_max)
1399 thus ?thesis unfolding l by auto
1401 case False with bnds 1 2 3 4 show ?thesis by (auto simp add: bnds_cos_def Let_def)
1405 section "Exponential function"
1407 subsection "Compute the series of the exponential function"
1409 fun ub_exp_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" and lb_exp_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
1410 "ub_exp_horner prec 0 i k x = 0" |
1411 "ub_exp_horner prec (Suc n) i k x = rapprox_rat prec 1 (int k) + x * lb_exp_horner prec n (i + 1) (k * i) x" |
1412 "lb_exp_horner prec 0 i k x = 0" |
1413 "lb_exp_horner prec (Suc n) i k x = lapprox_rat prec 1 (int k) + x * ub_exp_horner prec n (i + 1) (k * i) x"
1415 lemma bnds_exp_horner: assumes "real x \<le> 0"
1416 shows "exp (real x) \<in> { real (lb_exp_horner prec (get_even n) 1 1 x) .. real (ub_exp_horner prec (get_odd n) 1 1 x) }"
1419 have F: "\<And> m. ((\<lambda>i. i + 1) ^^ n) m = n + m" by (induct n, auto)
1420 have "fact (Suc n) = fact n * ((\<lambda>i. i + 1) ^^ n) 1" unfolding F by auto } note f_eq = this
1422 note bounds = horner_bounds_nonpos[where f="fact" and lb="lb_exp_horner prec" and ub="ub_exp_horner prec" and j'=0 and s=1,
1423 OF assms f_eq lb_exp_horner.simps ub_exp_horner.simps]
1425 { have "real (lb_exp_horner prec (get_even n) 1 1 x) \<le> (\<Sum>j = 0..<get_even n. 1 / real (fact j) * real x ^ j)"
1426 using bounds(1) by auto
1427 also have "\<dots> \<le> exp (real x)"
1429 obtain t where "\<bar>t\<bar> \<le> \<bar>real x\<bar>" and "exp (real x) = (\<Sum>m = 0..<get_even n. (real x) ^ m / real (fact m)) + exp t / real (fact (get_even n)) * (real x) ^ (get_even n)"
1430 using Maclaurin_exp_le by blast
1431 moreover have "0 \<le> exp t / real (fact (get_even n)) * (real x) ^ (get_even n)"
1432 by (auto intro!: mult_nonneg_nonneg divide_nonneg_pos simp add: get_even zero_le_even_power exp_gt_zero)
1433 ultimately show ?thesis
1434 using get_odd exp_gt_zero by (auto intro!: pordered_cancel_semiring_class.mult_nonneg_nonneg)
1436 finally have "real (lb_exp_horner prec (get_even n) 1 1 x) \<le> exp (real x)" .
1439 have x_less_zero: "real x ^ get_odd n \<le> 0"
1440 proof (cases "real x = 0")
1442 have "(get_odd n) \<noteq> 0" using get_odd[THEN odd_pos] by auto
1443 thus ?thesis unfolding True power_0_left by auto
1445 case False hence "real x < 0" using `real x \<le> 0` by auto
1446 show ?thesis by (rule less_imp_le, auto simp add: power_less_zero_eq get_odd `real x < 0`)
1449 obtain t where "\<bar>t\<bar> \<le> \<bar>real x\<bar>" and "exp (real x) = (\<Sum>m = 0..<get_odd n. (real x) ^ m / real (fact m)) + exp t / real (fact (get_odd n)) * (real x) ^ (get_odd n)"
1450 using Maclaurin_exp_le by blast
1451 moreover have "exp t / real (fact (get_odd n)) * (real x) ^ (get_odd n) \<le> 0"
1452 by (auto intro!: mult_nonneg_nonpos divide_nonpos_pos simp add: x_less_zero exp_gt_zero)
1453 ultimately have "exp (real x) \<le> (\<Sum>j = 0..<get_odd n. 1 / real (fact j) * real x ^ j)"
1454 using get_odd exp_gt_zero by (auto intro!: pordered_cancel_semiring_class.mult_nonneg_nonneg)
1455 also have "\<dots> \<le> real (ub_exp_horner prec (get_odd n) 1 1 x)"
1456 using bounds(2) by auto
1457 finally have "exp (real x) \<le> real (ub_exp_horner prec (get_odd n) 1 1 x)" .
1458 } ultimately show ?thesis by auto
1461 subsection "Compute the exponential function on the entire domain"
1463 function ub_exp :: "nat \<Rightarrow> float \<Rightarrow> float" and lb_exp :: "nat \<Rightarrow> float \<Rightarrow> float" where
1464 "lb_exp prec x = (if 0 < x then float_divl prec 1 (ub_exp prec (-x))
1466 horner = (\<lambda> x. let y = lb_exp_horner prec (get_even (prec + 2)) 1 1 x in if y \<le> 0 then Float 1 -2 else y)
1467 in if x < - 1 then (case floor_fl x of (Float m e) \<Rightarrow> (horner (float_divl prec x (- Float m e))) ^ (nat (-m) * 2 ^ nat e))
1469 "ub_exp prec x = (if 0 < x then float_divr prec 1 (lb_exp prec (-x))
1470 else if x < - 1 then (case floor_fl x of (Float m e) \<Rightarrow>
1471 (ub_exp_horner prec (get_odd (prec + 2)) 1 1 (float_divr prec x (- Float m e))) ^ (nat (-m) * 2 ^ nat e))
1472 else ub_exp_horner prec (get_odd (prec + 2)) 1 1 x)"
1473 by pat_completeness auto
1474 termination by (relation "measure (\<lambda> v. let (prec, x) = sum_case id id v in (if 0 < x then 1 else 0))", auto simp add: less_float_def)
1476 lemma exp_m1_ge_quarter: "(1 / 4 :: real) \<le> exp (- 1)"
1478 have eq4: "4 = Suc (Suc (Suc (Suc 0)))" by auto
1480 have "1 / 4 = real (Float 1 -2)" unfolding Float_num by auto
1481 also have "\<dots> \<le> real (lb_exp_horner 1 (get_even 4) 1 1 (- 1))"
1482 unfolding get_even_def eq4
1483 by (auto simp add: lapprox_posrat_def rapprox_posrat_def normfloat.simps)
1484 also have "\<dots> \<le> exp (real (- 1 :: float))" using bnds_exp_horner[where x="- 1"] by auto
1485 finally show ?thesis unfolding real_of_float_minus real_of_float_1 .
1488 lemma lb_exp_pos: assumes "\<not> 0 < x" shows "0 < lb_exp prec x"
1490 let "?lb_horner x" = "lb_exp_horner prec (get_even (prec + 2)) 1 1 x"
1491 let "?horner x" = "let y = ?lb_horner x in if y \<le> 0 then Float 1 -2 else y"
1492 have pos_horner: "\<And> x. 0 < ?horner x" unfolding Let_def by (cases "?lb_horner x \<le> 0", auto simp add: le_float_def less_float_def)
1493 moreover { fix x :: float fix num :: nat
1494 have "0 < real (?horner x) ^ num" using `0 < ?horner x`[unfolded less_float_def real_of_float_0] by (rule zero_less_power)
1495 also have "\<dots> = real ((?horner x) ^ num)" using float_power by auto
1496 finally have "0 < real ((?horner x) ^ num)" .
1498 ultimately show ?thesis
1499 unfolding lb_exp.simps if_not_P[OF `\<not> 0 < x`] Let_def
1500 by (cases "floor_fl x", cases "x < - 1", auto simp add: float_power le_float_def less_float_def)
1503 lemma exp_boundaries': assumes "x \<le> 0"
1504 shows "exp (real x) \<in> { real (lb_exp prec x) .. real (ub_exp prec x)}"
1506 let "?lb_exp_horner x" = "lb_exp_horner prec (get_even (prec + 2)) 1 1 x"
1507 let "?ub_exp_horner x" = "ub_exp_horner prec (get_odd (prec + 2)) 1 1 x"
1509 have "real x \<le> 0" and "\<not> x > 0" using `x \<le> 0` unfolding le_float_def less_float_def by auto
1511 proof (cases "x < - 1")
1512 case False hence "- 1 \<le> real x" unfolding less_float_def by auto
1514 proof (cases "?lb_exp_horner x \<le> 0")
1515 from `\<not> x < - 1` have "- 1 \<le> real x" unfolding less_float_def by auto
1516 hence "exp (- 1) \<le> exp (real x)" unfolding exp_le_cancel_iff .
1517 from order_trans[OF exp_m1_ge_quarter this]
1518 have "real (Float 1 -2) \<le> exp (real x)" unfolding Float_num .
1520 ultimately show ?thesis using bnds_exp_horner `real x \<le> 0` `\<not> x > 0` `\<not> x < - 1` by auto
1522 case False thus ?thesis using bnds_exp_horner `real x \<le> 0` `\<not> x > 0` `\<not> x < - 1` by (auto simp add: Let_def)
1527 obtain m e where Float_floor: "floor_fl x = Float m e" by (cases "floor_fl x", auto)
1528 let ?num = "nat (- m) * 2 ^ nat e"
1530 have "real (floor_fl x) < - 1" using floor_fl `x < - 1` unfolding le_float_def less_float_def real_of_float_minus real_of_float_1 by (rule order_le_less_trans)
1531 hence "real (floor_fl x) < 0" unfolding Float_floor real_of_float_simp using zero_less_pow2[of xe] by auto
1533 unfolding less_float_def real_of_float_0 Float_floor real_of_float_simp
1534 unfolding pos_prod_lt[OF zero_less_pow2[of e], unfolded real_mult_commute] by auto
1535 hence "1 \<le> - m" by auto
1536 hence "0 < nat (- m)" by auto
1538 have "0 \<le> e" using floor_pos_exp Float_floor[symmetric] by auto
1539 hence "(0::nat) < 2 ^ nat e" by auto
1540 ultimately have "0 < ?num" by auto
1541 hence "real ?num \<noteq> 0" by auto
1542 have e_nat: "int (nat e) = e" using `0 \<le> e` by auto
1543 have num_eq: "real ?num = real (- floor_fl x)" using `0 < nat (- m)`
1544 unfolding Float_floor real_of_float_minus real_of_float_simp real_of_nat_mult pow2_int[of "nat e", unfolded e_nat] realpow_real_of_nat[symmetric] by auto
1545 have "0 < - floor_fl x" using `0 < ?num`[unfolded real_of_nat_less_iff[symmetric]] unfolding less_float_def num_eq[symmetric] real_of_float_0 real_of_nat_zero .
1546 hence "real (floor_fl x) < 0" unfolding less_float_def by auto
1548 have "exp (real x) \<le> real (ub_exp prec x)"
1550 have div_less_zero: "real (float_divr prec x (- floor_fl x)) \<le> 0"
1551 using float_divr_nonpos_pos_upper_bound[OF `x \<le> 0` `0 < - floor_fl x`] unfolding le_float_def real_of_float_0 .
1553 have "exp (real x) = exp (real ?num * (real x / real ?num))" using `real ?num \<noteq> 0` by auto
1554 also have "\<dots> = exp (real x / real ?num) ^ ?num" unfolding exp_real_of_nat_mult ..
1555 also have "\<dots> \<le> exp (real (float_divr prec x (- floor_fl x))) ^ ?num" unfolding num_eq
1556 by (rule power_mono, rule exp_le_cancel_iff[THEN iffD2], rule float_divr) auto
1557 also have "\<dots> \<le> real ((?ub_exp_horner (float_divr prec x (- floor_fl x))) ^ ?num)" unfolding float_power
1558 by (rule power_mono, rule bnds_exp_horner[OF div_less_zero, unfolded atLeastAtMost_iff, THEN conjunct2], auto)
1559 finally show ?thesis unfolding ub_exp.simps if_not_P[OF `\<not> 0 < x`] if_P[OF `x < - 1`] float.cases Float_floor Let_def .
1562 have "real (lb_exp prec x) \<le> exp (real x)"
1564 let ?divl = "float_divl prec x (- Float m e)"
1565 let ?horner = "?lb_exp_horner ?divl"
1568 proof (cases "?horner \<le> 0")
1569 case False hence "0 \<le> real ?horner" unfolding le_float_def by auto
1571 have div_less_zero: "real (float_divl prec x (- floor_fl x)) \<le> 0"
1572 using `real (floor_fl x) < 0` `real x \<le> 0` by (auto intro!: order_trans[OF float_divl] divide_nonpos_neg)
1574 have "real ((?lb_exp_horner (float_divl prec x (- floor_fl x))) ^ ?num) \<le>
1575 exp (real (float_divl prec x (- floor_fl x))) ^ ?num" unfolding float_power
1576 using `0 \<le> real ?horner`[unfolded Float_floor[symmetric]] bnds_exp_horner[OF div_less_zero, unfolded atLeastAtMost_iff, THEN conjunct1] by (auto intro!: power_mono)
1577 also have "\<dots> \<le> exp (real x / real ?num) ^ ?num" unfolding num_eq
1578 using float_divl by (auto intro!: power_mono simp del: real_of_float_minus)
1579 also have "\<dots> = exp (real ?num * (real x / real ?num))" unfolding exp_real_of_nat_mult ..
1580 also have "\<dots> = exp (real x)" using `real ?num \<noteq> 0` by auto
1581 finally show ?thesis
1582 unfolding lb_exp.simps if_not_P[OF `\<not> 0 < x`] if_P[OF `x < - 1`] float.cases Float_floor Let_def if_not_P[OF False] by auto
1585 have "real (floor_fl x) \<noteq> 0" and "real (floor_fl x) \<le> 0" using `real (floor_fl x) < 0` by auto
1586 from divide_right_mono_neg[OF floor_fl[of x] `real (floor_fl x) \<le> 0`, unfolded divide_self[OF `real (floor_fl x) \<noteq> 0`]]
1587 have "- 1 \<le> real x / real (- floor_fl x)" unfolding real_of_float_minus by auto
1588 from order_trans[OF exp_m1_ge_quarter this[unfolded exp_le_cancel_iff[where x="- 1", symmetric]]]
1589 have "real (Float 1 -2) \<le> exp (real x / real (- floor_fl x))" unfolding Float_num .
1590 hence "real (Float 1 -2) ^ ?num \<le> exp (real x / real (- floor_fl x)) ^ ?num"
1591 by (auto intro!: power_mono simp add: Float_num)
1592 also have "\<dots> = exp (real x)" unfolding num_eq exp_real_of_nat_mult[symmetric] using `real (floor_fl x) \<noteq> 0` by auto
1593 finally show ?thesis
1594 unfolding lb_exp.simps if_not_P[OF `\<not> 0 < x`] if_P[OF `x < - 1`] float.cases Float_floor Let_def if_P[OF True] float_power .
1597 ultimately show ?thesis by auto
1601 lemma exp_boundaries: "exp (real x) \<in> { real (lb_exp prec x) .. real (ub_exp prec x)}"
1604 proof (cases "0 < x")
1605 case False hence "x \<le> 0" unfolding less_float_def le_float_def by auto
1606 from exp_boundaries'[OF this] show ?thesis .
1608 case True hence "-x \<le> 0" unfolding less_float_def le_float_def by auto
1610 have "real (lb_exp prec x) \<le> exp (real x)"
1612 from exp_boundaries'[OF `-x \<le> 0`]
1613 have ub_exp: "exp (- real x) \<le> real (ub_exp prec (-x))" unfolding atLeastAtMost_iff real_of_float_minus by auto
1615 have "real (float_divl prec 1 (ub_exp prec (-x))) \<le> 1 / real (ub_exp prec (-x))" using float_divl[where x=1] by auto
1616 also have "\<dots> \<le> exp (real x)"
1617 using ub_exp[unfolded inverse_le_iff_le[OF order_less_le_trans[OF exp_gt_zero ub_exp] exp_gt_zero, symmetric]]
1618 unfolding exp_minus nonzero_inverse_inverse_eq[OF exp_not_eq_zero] inverse_eq_divide by auto
1619 finally show ?thesis unfolding lb_exp.simps if_P[OF True] .
1622 have "exp (real x) \<le> real (ub_exp prec x)"
1624 have "\<not> 0 < -x" using `0 < x` unfolding less_float_def by auto
1626 from exp_boundaries'[OF `-x \<le> 0`]
1627 have lb_exp: "real (lb_exp prec (-x)) \<le> exp (- real x)" unfolding atLeastAtMost_iff real_of_float_minus by auto
1629 have "exp (real x) \<le> real (1 :: float) / real (lb_exp prec (-x))"
1630 using lb_exp[unfolded inverse_le_iff_le[OF exp_gt_zero lb_exp_pos[OF `\<not> 0 < -x`, unfolded less_float_def real_of_float_0],
1632 unfolding exp_minus nonzero_inverse_inverse_eq[OF exp_not_eq_zero] inverse_eq_divide real_of_float_1 by auto
1633 also have "\<dots> \<le> real (float_divr prec 1 (lb_exp prec (-x)))" using float_divr .
1634 finally show ?thesis unfolding ub_exp.simps if_P[OF True] .
1636 ultimately show ?thesis by auto
1640 lemma bnds_exp: "\<forall> x lx ux. (l, u) = (lb_exp prec lx, ub_exp prec ux) \<and> x \<in> {real lx .. real ux} \<longrightarrow> real l \<le> exp x \<and> exp x \<le> real u"
1641 proof (rule allI, rule allI, rule allI, rule impI)
1643 assume "(l, u) = (lb_exp prec lx, ub_exp prec ux) \<and> x \<in> {real lx .. real ux}"
1644 hence l: "lb_exp prec lx = l " and u: "ub_exp prec ux = u" and x: "x \<in> {real lx .. real ux}" by auto
1646 { from exp_boundaries[of lx prec, unfolded l]
1647 have "real l \<le> exp (real lx)" by (auto simp del: lb_exp.simps)
1648 also have "\<dots> \<le> exp x" using x by auto
1649 finally have "real l \<le> exp x" .
1651 { have "exp x \<le> exp (real ux)" using x by auto
1652 also have "\<dots> \<le> real u" using exp_boundaries[of ux prec, unfolded u] by (auto simp del: ub_exp.simps)
1653 finally have "exp x \<le> real u" .
1654 } ultimately show "real l \<le> exp x \<and> exp x \<le> real u" ..
1659 subsection "Compute the logarithm series"
1661 fun ub_ln_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
1662 and lb_ln_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
1663 "ub_ln_horner prec 0 i x = 0" |
1664 "ub_ln_horner prec (Suc n) i x = rapprox_rat prec 1 (int i) - x * lb_ln_horner prec n (Suc i) x" |
1665 "lb_ln_horner prec 0 i x = 0" |
1666 "lb_ln_horner prec (Suc n) i x = lapprox_rat prec 1 (int i) - x * ub_ln_horner prec n (Suc i) x"
1669 assumes "0 \<le> x" and "x < 1"
1670 shows "(\<Sum>i=0..<2*n. -1^i * (1 / real (i + 1)) * x ^ (Suc i)) \<le> ln (x + 1)" (is "?lb")
1671 and "ln (x + 1) \<le> (\<Sum>i=0..<2*n + 1. -1^i * (1 / real (i + 1)) * x ^ (Suc i))" (is "?ub")
1673 let "?a n" = "(1/real (n +1)) * x ^ (Suc n)"
1675 have ln_eq: "(\<Sum> i. -1^i * ?a i) = ln (x + 1)"
1676 using ln_series[of "x + 1"] `0 \<le> x` `x < 1` by auto
1678 have "norm x < 1" using assms by auto
1679 have "?a ----> 0" unfolding Suc_eq_plus1[symmetric] inverse_eq_divide[symmetric]
1680 using LIMSEQ_mult[OF LIMSEQ_inverse_real_of_nat LIMSEQ_Suc[OF LIMSEQ_power_zero[OF `norm x < 1`]]] by auto
1681 { fix n have "0 \<le> ?a n" by (rule mult_nonneg_nonneg, auto intro!: mult_nonneg_nonneg simp add: `0 \<le> x`) }
1682 { fix n have "?a (Suc n) \<le> ?a n" unfolding inverse_eq_divide[symmetric]
1683 proof (rule mult_mono)
1684 show "0 \<le> x ^ Suc (Suc n)" by (auto intro!: mult_nonneg_nonneg simp add: `0 \<le> x`)
1685 have "x ^ Suc (Suc n) \<le> x ^ Suc n * 1" unfolding power_Suc2 real_mult_assoc[symmetric]
1686 by (rule mult_left_mono, fact less_imp_le[OF `x < 1`], auto intro!: mult_nonneg_nonneg simp add: `0 \<le> x`)
1687 thus "x ^ Suc (Suc n) \<le> x ^ Suc n" by auto
1689 from summable_Leibniz'(2,4)[OF `?a ----> 0` `\<And>n. 0 \<le> ?a n`, OF `\<And>n. ?a (Suc n) \<le> ?a n`, unfolded ln_eq]
1690 show "?lb" and "?ub" by auto
1693 lemma ln_float_bounds:
1694 assumes "0 \<le> real x" and "real x < 1"
1695 shows "real (x * lb_ln_horner prec (get_even n) 1 x) \<le> ln (real x + 1)" (is "?lb \<le> ?ln")
1696 and "ln (real x + 1) \<le> real (x * ub_ln_horner prec (get_odd n) 1 x)" (is "?ln \<le> ?ub")
1698 obtain ev where ev: "get_even n = 2 * ev" using get_even_double ..
1699 obtain od where od: "get_odd n = 2 * od + 1" using get_odd_double ..
1701 let "?s n" = "-1^n * (1 / real (1 + n)) * (real x)^(Suc n)"
1703 have "?lb \<le> setsum ?s {0 ..< 2 * ev}" unfolding power_Suc2 real_mult_assoc[symmetric] real_of_float_mult setsum_left_distrib[symmetric] unfolding real_mult_commute[of "real x"] ev
1704 using horner_bounds(1)[where G="\<lambda> i k. Suc k" and F="\<lambda>x. x" and f="\<lambda>x. x" and lb="\<lambda>n i k x. lb_ln_horner prec n k x" and ub="\<lambda>n i k x. ub_ln_horner prec n k x" and j'=1 and n="2*ev",
1705 OF `0 \<le> real x` refl lb_ln_horner.simps ub_ln_horner.simps] `0 \<le> real x`
1706 by (rule mult_right_mono)
1707 also have "\<dots> \<le> ?ln" using ln_bounds(1)[OF `0 \<le> real x` `real x < 1`] by auto
1708 finally show "?lb \<le> ?ln" .
1710 have "?ln \<le> setsum ?s {0 ..< 2 * od + 1}" using ln_bounds(2)[OF `0 \<le> real x` `real x < 1`] by auto
1711 also have "\<dots> \<le> ?ub" unfolding power_Suc2 real_mult_assoc[symmetric] real_of_float_mult setsum_left_distrib[symmetric] unfolding real_mult_commute[of "real x"] od
1712 using horner_bounds(2)[where G="\<lambda> i k. Suc k" and F="\<lambda>x. x" and f="\<lambda>x. x" and lb="\<lambda>n i k x. lb_ln_horner prec n k x" and ub="\<lambda>n i k x. ub_ln_horner prec n k x" and j'=1 and n="2*od+1",
1713 OF `0 \<le> real x` refl lb_ln_horner.simps ub_ln_horner.simps] `0 \<le> real x`
1714 by (rule mult_right_mono)
1715 finally show "?ln \<le> ?ub" .
1718 lemma ln_add: assumes "0 < x" and "0 < y" shows "ln (x + y) = ln x + ln (1 + y / x)"
1720 have "x \<noteq> 0" using assms by auto
1721 have "x + y = x * (1 + y / x)" unfolding right_distrib times_divide_eq_right nonzero_mult_divide_cancel_left[OF `x \<noteq> 0`] by auto
1723 have "0 < y / x" using assms divide_pos_pos by auto
1724 hence "0 < 1 + y / x" by auto
1725 ultimately show ?thesis using ln_mult assms by auto
1728 subsection "Compute the logarithm of 2"
1730 definition ub_ln2 where "ub_ln2 prec = (let third = rapprox_rat (max prec 1) 1 3
1731 in (Float 1 -1 * ub_ln_horner prec (get_odd prec) 1 (Float 1 -1)) +
1732 (third * ub_ln_horner prec (get_odd prec) 1 third))"
1733 definition lb_ln2 where "lb_ln2 prec = (let third = lapprox_rat prec 1 3
1734 in (Float 1 -1 * lb_ln_horner prec (get_even prec) 1 (Float 1 -1)) +
1735 (third * lb_ln_horner prec (get_even prec) 1 third))"
1737 lemma ub_ln2: "ln 2 \<le> real (ub_ln2 prec)" (is "?ub_ln2")
1738 and lb_ln2: "real (lb_ln2 prec) \<le> ln 2" (is "?lb_ln2")
1740 let ?uthird = "rapprox_rat (max prec 1) 1 3"
1741 let ?lthird = "lapprox_rat prec 1 3"
1743 have ln2_sum: "ln 2 = ln (1/2 + 1) + ln (1 / 3 + 1)"
1744 using ln_add[of "3 / 2" "1 / 2"] by auto
1745 have lb3: "real ?lthird \<le> 1 / 3" using lapprox_rat[of prec 1 3] by auto
1746 hence lb3_ub: "real ?lthird < 1" by auto
1747 have lb3_lb: "0 \<le> real ?lthird" using lapprox_rat_bottom[of 1 3] by auto
1748 have ub3: "1 / 3 \<le> real ?uthird" using rapprox_rat[of 1 3] by auto
1749 hence ub3_lb: "0 \<le> real ?uthird" by auto
1751 have lb2: "0 \<le> real (Float 1 -1)" and ub2: "real (Float 1 -1) < 1" unfolding Float_num by auto
1753 have "0 \<le> (1::int)" and "0 < (3::int)" by auto
1754 have ub3_ub: "real ?uthird < 1" unfolding rapprox_rat.simps(2)[OF `0 \<le> 1` `0 < 3`]
1755 by (rule rapprox_posrat_less1, auto)
1757 have third_gt0: "(0 :: real) < 1 / 3 + 1" by auto
1758 have uthird_gt0: "0 < real ?uthird + 1" using ub3_lb by auto
1759 have lthird_gt0: "0 < real ?lthird + 1" using lb3_lb by auto
1761 show ?ub_ln2 unfolding ub_ln2_def Let_def real_of_float_add ln2_sum Float_num(4)[symmetric]
1762 proof (rule add_mono, fact ln_float_bounds(2)[OF lb2 ub2])
1763 have "ln (1 / 3 + 1) \<le> ln (real ?uthird + 1)" unfolding ln_le_cancel_iff[OF third_gt0 uthird_gt0] using ub3 by auto
1764 also have "\<dots> \<le> real (?uthird * ub_ln_horner prec (get_odd prec) 1 ?uthird)"
1765 using ln_float_bounds(2)[OF ub3_lb ub3_ub] .
1766 finally show "ln (1 / 3 + 1) \<le> real (?uthird * ub_ln_horner prec (get_odd prec) 1 ?uthird)" .
1768 show ?lb_ln2 unfolding lb_ln2_def Let_def real_of_float_add ln2_sum Float_num(4)[symmetric]
1769 proof (rule add_mono, fact ln_float_bounds(1)[OF lb2 ub2])
1770 have "real (?lthird * lb_ln_horner prec (get_even prec) 1 ?lthird) \<le> ln (real ?lthird + 1)"
1771 using ln_float_bounds(1)[OF lb3_lb lb3_ub] .
1772 also have "\<dots> \<le> ln (1 / 3 + 1)" unfolding ln_le_cancel_iff[OF lthird_gt0 third_gt0] using lb3 by auto
1773 finally show "real (?lthird * lb_ln_horner prec (get_even prec) 1 ?lthird) \<le> ln (1 / 3 + 1)" .
1777 subsection "Compute the logarithm in the entire domain"
1779 function ub_ln :: "nat \<Rightarrow> float \<Rightarrow> float option" and lb_ln :: "nat \<Rightarrow> float \<Rightarrow> float option" where
1780 "ub_ln prec x = (if x \<le> 0 then None
1781 else if x < 1 then Some (- the (lb_ln prec (float_divl (max prec 1) 1 x)))
1782 else let horner = \<lambda>x. x * ub_ln_horner prec (get_odd prec) 1 x in
1783 if x \<le> Float 3 -1 then Some (horner (x - 1))
1784 else if x < Float 1 1 then Some (horner (Float 1 -1) + horner (x * rapprox_rat prec 2 3 - 1))
1785 else let l = bitlen (mantissa x) - 1 in
1786 Some (ub_ln2 prec * (Float (scale x + l) 0) + horner (Float (mantissa x) (- l) - 1)))" |
1787 "lb_ln prec x = (if x \<le> 0 then None
1788 else if x < 1 then Some (- the (ub_ln prec (float_divr prec 1 x)))
1789 else let horner = \<lambda>x. x * lb_ln_horner prec (get_even prec) 1 x in
1790 if x \<le> Float 3 -1 then Some (horner (x - 1))
1791 else if x < Float 1 1 then Some (horner (Float 1 -1) +
1792 horner (max (x * lapprox_rat prec 2 3 - 1) 0))
1793 else let l = bitlen (mantissa x) - 1 in
1794 Some (lb_ln2 prec * (Float (scale x + l) 0) + horner (Float (mantissa x) (- l) - 1)))"
1795 by pat_completeness auto
1797 termination proof (relation "measure (\<lambda> v. let (prec, x) = sum_case id id v in (if x < 1 then 1 else 0))", auto)
1798 fix prec x assume "\<not> x \<le> 0" and "x < 1" and "float_divl (max prec (Suc 0)) 1 x < 1"
1799 hence "0 < x" and "0 < max prec (Suc 0)" unfolding less_float_def le_float_def by auto
1800 from float_divl_pos_less1_bound[OF `0 < x` `x < 1` `0 < max prec (Suc 0)`]
1801 show False using `float_divl (max prec (Suc 0)) 1 x < 1` unfolding less_float_def le_float_def by auto
1803 fix prec x assume "\<not> x \<le> 0" and "x < 1" and "float_divr prec 1 x < 1"
1804 hence "0 < x" unfolding less_float_def le_float_def by auto
1805 from float_divr_pos_less1_lower_bound[OF `0 < x` `x < 1`, of prec]
1806 show False using `float_divr prec 1 x < 1` unfolding less_float_def le_float_def by auto
1809 lemma ln_shifted_float: assumes "0 < m" shows "ln (real (Float m e)) = ln 2 * real (e + (bitlen m - 1)) + ln (real (Float m (- (bitlen m - 1))))"
1811 let ?B = "2^nat (bitlen m - 1)"
1812 have "0 < real m" and "\<And>X. (0 :: real) < 2^X" and "0 < (2 :: real)" and "m \<noteq> 0" using assms by auto
1813 hence "0 \<le> bitlen m - 1" using bitlen_ge1[OF `m \<noteq> 0`] by auto
1815 proof (cases "0 \<le> e")
1817 show ?thesis unfolding normalized_float[OF `m \<noteq> 0`]
1818 unfolding ln_div[OF `0 < real m` `0 < ?B`] real_of_int_add ln_realpow[OF `0 < 2`]
1819 unfolding real_of_float_ge0_exp[OF True] ln_mult[OF `0 < real m` `0 < 2^nat e`]
1820 ln_realpow[OF `0 < 2`] algebra_simps using `0 \<le> bitlen m - 1` True by auto
1822 case False hence "0 < -e" by auto
1823 hence pow_gt0: "(0::real) < 2^nat (-e)" by auto
1824 hence inv_gt0: "(0::real) < inverse (2^nat (-e))" by auto
1825 show ?thesis unfolding normalized_float[OF `m \<noteq> 0`]
1826 unfolding ln_div[OF `0 < real m` `0 < ?B`] real_of_int_add ln_realpow[OF `0 < 2`]
1827 unfolding real_of_float_nge0_exp[OF False] ln_mult[OF `0 < real m` inv_gt0] ln_inverse[OF pow_gt0]
1828 ln_realpow[OF `0 < 2`] algebra_simps using `0 \<le> bitlen m - 1` False by auto
1832 lemma ub_ln_lb_ln_bounds': assumes "1 \<le> x"
1833 shows "real (the (lb_ln prec x)) \<le> ln (real x) \<and> ln (real x) \<le> real (the (ub_ln prec x))"
1834 (is "?lb \<le> ?ln \<and> ?ln \<le> ?ub")
1835 proof (cases "x < Float 1 1")
1837 hence "real (x - 1) < 1" and "real x < 2" unfolding less_float_def Float_num by auto
1838 have "\<not> x \<le> 0" and "\<not> x < 1" using `1 \<le> x` unfolding less_float_def le_float_def by auto
1839 hence "0 \<le> real (x - 1)" using `1 \<le> x` unfolding less_float_def Float_num by auto
1841 have [simp]: "real (Float 3 -1) = 3 / 2" by (simp add: real_of_float_def pow2_def)
1844 proof (cases "x \<le> Float 3 -1")
1846 show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps Let_def
1847 using ln_float_bounds[OF `0 \<le> real (x - 1)` `real (x - 1) < 1`, of prec] `\<not> x \<le> 0` `\<not> x < 1` True
1850 case False hence *: "3 / 2 < real x" by (auto simp add: le_float_def)
1852 with ln_add[of "3 / 2" "real x - 3 / 2"]
1853 have add: "ln (real x) = ln (3 / 2) + ln (real x * 2 / 3)"
1854 by (auto simp add: algebra_simps diff_divide_distrib)
1856 let "?ub_horner x" = "x * ub_ln_horner prec (get_odd prec) 1 x"
1857 let "?lb_horner x" = "x * lb_ln_horner prec (get_even prec) 1 x"
1859 { have up: "real (rapprox_rat prec 2 3) \<le> 1"
1860 by (rule rapprox_rat_le1) simp_all
1861 have low: "2 / 3 \<le> real (rapprox_rat prec 2 3)"
1862 by (rule order_trans[OF _ rapprox_rat]) simp
1863 from mult_less_le_imp_less[OF * low] *
1864 have pos: "0 < real (x * rapprox_rat prec 2 3 - 1)" by auto
1866 have "ln (real x * 2/3)
1867 \<le> ln (real (x * rapprox_rat prec 2 3 - 1) + 1)"
1868 proof (rule ln_le_cancel_iff[symmetric, THEN iffD1])
1869 show "real x * 2 / 3 \<le> real (x * rapprox_rat prec 2 3 - 1) + 1"
1871 show "0 < real x * 2 / 3" using * by simp
1872 show "0 < real (x * rapprox_rat prec 2 3 - 1) + 1" using pos by auto
1874 also have "\<dots> \<le> real (?ub_horner (x * rapprox_rat prec 2 3 - 1))"
1875 proof (rule ln_float_bounds(2))
1876 from mult_less_le_imp_less[OF `real x < 2` up] low *
1877 show "real (x * rapprox_rat prec 2 3 - 1) < 1" by auto
1878 show "0 \<le> real (x * rapprox_rat prec 2 3 - 1)" using pos by auto
1880 finally have "ln (real x)
1881 \<le> real (?ub_horner (Float 1 -1))
1882 + real (?ub_horner (x * rapprox_rat prec 2 3 - 1))"
1883 using ln_float_bounds(2)[of "Float 1 -1" prec prec] add by auto }
1885 { let ?max = "max (x * lapprox_rat prec 2 3 - 1) 0"
1887 have up: "real (lapprox_rat prec 2 3) \<le> 2/3"
1888 by (rule order_trans[OF lapprox_rat], simp)
1890 have low: "0 \<le> real (lapprox_rat prec 2 3)"
1891 using lapprox_rat_bottom[of 2 3 prec] by simp
1893 have "real (?lb_horner ?max)
1894 \<le> ln (real ?max + 1)"
1895 proof (rule ln_float_bounds(1))
1896 from mult_less_le_imp_less[OF `real x < 2` up] * low
1897 show "real ?max < 1" by (cases "real (lapprox_rat prec 2 3) = 0",
1898 auto simp add: real_of_float_max)
1899 show "0 \<le> real ?max" by (auto simp add: real_of_float_max)
1901 also have "\<dots> \<le> ln (real x * 2/3)"
1902 proof (rule ln_le_cancel_iff[symmetric, THEN iffD1])
1903 show "0 < real ?max + 1" by (auto simp add: real_of_float_max)
1904 show "0 < real x * 2/3" using * by auto
1905 show "real ?max + 1 \<le> real x * 2/3" using * up
1906 by (cases "0 < real x * real (lapprox_posrat prec 2 3) - 1",
1907 auto simp add: real_of_float_max max_def)
1909 finally have "real (?lb_horner (Float 1 -1)) + real (?lb_horner ?max)
1911 using ln_float_bounds(1)[of "Float 1 -1" prec prec] add by auto }
1913 show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps Let_def
1914 using `\<not> x \<le> 0` `\<not> x < 1` True False by auto
1918 hence "\<not> x \<le> 0" and "\<not> x < 1" "0 < x" "\<not> x \<le> Float 3 -1"
1919 using `1 \<le> x` unfolding less_float_def le_float_def real_of_float_simp pow2_def
1924 let ?s = "Float (e + (bitlen m - 1)) 0"
1925 let ?x = "Float m (- (bitlen m - 1))"
1927 have "0 < m" and "m \<noteq> 0" using float_pos_m_pos `0 < x` Float by auto
1930 have "real (lb_ln2 prec * ?s) \<le> ln 2 * real (e + (bitlen m - 1))" (is "?lb2 \<le> _")
1931 unfolding real_of_float_mult real_of_float_ge0_exp[OF order_refl] nat_0 power_0 mult_1_right
1932 using lb_ln2[of prec]
1933 proof (rule mult_right_mono)
1934 have "1 \<le> Float m e" using `1 \<le> x` Float unfolding le_float_def by auto
1935 from float_gt1_scale[OF this]
1936 show "0 \<le> real (e + (bitlen m - 1))" by auto
1939 from bitlen_div[OF `0 < m`, unfolded normalized_float[OF `m \<noteq> 0`, symmetric]]
1940 have "0 \<le> real (?x - 1)" and "real (?x - 1) < 1" by auto
1941 from ln_float_bounds(1)[OF this]
1942 have "real ((?x - 1) * lb_ln_horner prec (get_even prec) 1 (?x - 1)) \<le> ln (real ?x)" (is "?lb_horner \<le> _") by auto
1943 ultimately have "?lb2 + ?lb_horner \<le> ln (real x)"
1944 unfolding Float ln_shifted_float[OF `0 < m`, of e] by auto
1948 from bitlen_div[OF `0 < m`, unfolded normalized_float[OF `m \<noteq> 0`, symmetric]]
1949 have "0 \<le> real (?x - 1)" and "real (?x - 1) < 1" by auto
1950 from ln_float_bounds(2)[OF this]
1951 have "ln (real ?x) \<le> real ((?x - 1) * ub_ln_horner prec (get_odd prec) 1 (?x - 1))" (is "_ \<le> ?ub_horner") by auto
1953 have "ln 2 * real (e + (bitlen m - 1)) \<le> real (ub_ln2 prec * ?s)" (is "_ \<le> ?ub2")
1954 unfolding real_of_float_mult real_of_float_ge0_exp[OF order_refl] nat_0 power_0 mult_1_right
1955 using ub_ln2[of prec]
1956 proof (rule mult_right_mono)
1957 have "1 \<le> Float m e" using `1 \<le> x` Float unfolding le_float_def by auto
1958 from float_gt1_scale[OF this]
1959 show "0 \<le> real (e + (bitlen m - 1))" by auto
1961 ultimately have "ln (real x) \<le> ?ub2 + ?ub_horner"
1962 unfolding Float ln_shifted_float[OF `0 < m`, of e] by auto
1964 ultimately show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps
1965 unfolding if_not_P[OF `\<not> x \<le> 0`] if_not_P[OF `\<not> x < 1`] if_not_P[OF False] if_not_P[OF `\<not> x \<le> Float 3 -1`] Let_def
1966 unfolding scale.simps[of m e, unfolded Float[symmetric]] mantissa.simps[of m e, unfolded Float[symmetric]] real_of_float_add
1971 lemma ub_ln_lb_ln_bounds: assumes "0 < x"
1972 shows "real (the (lb_ln prec x)) \<le> ln (real x) \<and> ln (real x) \<le> real (the (ub_ln prec x))"
1973 (is "?lb \<le> ?ln \<and> ?ln \<le> ?ub")
1974 proof (cases "x < 1")
1975 case False hence "1 \<le> x" unfolding less_float_def le_float_def by auto
1976 show ?thesis using ub_ln_lb_ln_bounds'[OF `1 \<le> x`] .
1978 case True have "\<not> x \<le> 0" using `0 < x` unfolding less_float_def le_float_def by auto
1980 have "0 < real x" and "real x \<noteq> 0" using `0 < x` unfolding less_float_def by auto
1981 hence A: "0 < 1 / real x" by auto
1984 let ?divl = "float_divl (max prec 1) 1 x"
1985 have A': "1 \<le> ?divl" using float_divl_pos_less1_bound[OF `0 < x` `x < 1`] unfolding le_float_def less_float_def by auto
1986 hence B: "0 < real ?divl" unfolding le_float_def by auto
1988 have "ln (real ?divl) \<le> ln (1 / real x)" unfolding ln_le_cancel_iff[OF B A] using float_divl[of _ 1 x] by auto
1989 hence "ln (real x) \<le> - ln (real ?divl)" unfolding nonzero_inverse_eq_divide[OF `real x \<noteq> 0`, symmetric] ln_inverse[OF `0 < real x`] by auto
1990 from this ub_ln_lb_ln_bounds'[OF A', THEN conjunct1, THEN le_imp_neg_le]
1991 have "?ln \<le> real (- the (lb_ln prec ?divl))" unfolding real_of_float_minus by (rule order_trans)
1994 let ?divr = "float_divr prec 1 x"
1995 have A': "1 \<le> ?divr" using float_divr_pos_less1_lower_bound[OF `0 < x` `x < 1`] unfolding le_float_def less_float_def by auto
1996 hence B: "0 < real ?divr" unfolding le_float_def by auto
1998 have "ln (1 / real x) \<le> ln (real ?divr)" unfolding ln_le_cancel_iff[OF A B] using float_divr[of 1 x] by auto
1999 hence "- ln (real ?divr) \<le> ln (real x)" unfolding nonzero_inverse_eq_divide[OF `real x \<noteq> 0`, symmetric] ln_inverse[OF `0 < real x`] by auto
2000 from ub_ln_lb_ln_bounds'[OF A', THEN conjunct2, THEN le_imp_neg_le] this
2001 have "real (- the (ub_ln prec ?divr)) \<le> ?ln" unfolding real_of_float_minus by (rule order_trans)
2003 ultimately show ?thesis unfolding lb_ln.simps[where x=x] ub_ln.simps[where x=x]
2004 unfolding if_not_P[OF `\<not> x \<le> 0`] if_P[OF True] by auto
2007 lemma lb_ln: assumes "Some y = lb_ln prec x"
2008 shows "real y \<le> ln (real x)" and "0 < real x"
2012 assume "\<not> 0 < x" hence "x \<le> 0" unfolding le_float_def less_float_def by auto
2013 thus False using assms by auto
2015 thus "0 < real x" unfolding less_float_def by auto
2016 have "real (the (lb_ln prec x)) \<le> ln (real x)" using ub_ln_lb_ln_bounds[OF `0 < x`] ..
2017 thus "real y \<le> ln (real x)" unfolding assms[symmetric] by auto
2020 lemma ub_ln: assumes "Some y = ub_ln prec x"
2021 shows "ln (real x) \<le> real y" and "0 < real x"
2025 assume "\<not> 0 < x" hence "x \<le> 0" unfolding le_float_def less_float_def by auto
2026 thus False using assms by auto
2028 thus "0 < real x" unfolding less_float_def by auto
2029 have "ln (real x) \<le> real (the (ub_ln prec x))" using ub_ln_lb_ln_bounds[OF `0 < x`] ..
2030 thus "ln (real x) \<le> real y" unfolding assms[symmetric] by auto
2033 lemma bnds_ln: "\<forall> x lx ux. (Some l, Some u) = (lb_ln prec lx, ub_ln prec ux) \<and> x \<in> {real lx .. real ux} \<longrightarrow> real l \<le> ln x \<and> ln x \<le> real u"
2034 proof (rule allI, rule allI, rule allI, rule impI)
2036 assume "(Some l, Some u) = (lb_ln prec lx, ub_ln prec ux) \<and> x \<in> {real lx .. real ux}"
2037 hence l: "Some l = lb_ln prec lx " and u: "Some u = ub_ln prec ux" and x: "x \<in> {real lx .. real ux}" by auto
2039 have "ln (real ux) \<le> real u" and "0 < real ux" using ub_ln u by auto
2040 have "real l \<le> ln (real lx)" and "0 < real lx" and "0 < x" using lb_ln[OF l] x by auto
2042 from ln_le_cancel_iff[OF `0 < real lx` `0 < x`] `real l \<le> ln (real lx)`
2043 have "real l \<le> ln x" using x unfolding atLeastAtMost_iff by auto
2045 from ln_le_cancel_iff[OF `0 < x` `0 < real ux`] `ln (real ux) \<le> real u`
2046 have "ln x \<le> real u" using x unfolding atLeastAtMost_iff by auto
2047 ultimately show "real l \<le> ln x \<and> ln x \<le> real u" ..
2050 section "Implement floatarith"
2052 subsection "Define syntax and semantics"
2055 = Add floatarith floatarith
2057 | Mult floatarith floatarith
2058 | Inverse floatarith
2062 | Max floatarith floatarith
2063 | Min floatarith floatarith
2068 | Power floatarith nat
2072 fun interpret_floatarith :: "floatarith \<Rightarrow> real list \<Rightarrow> real"
2074 "interpret_floatarith (Add a b) vs = (interpret_floatarith a vs) + (interpret_floatarith b vs)" |
2075 "interpret_floatarith (Minus a) vs = - (interpret_floatarith a vs)" |
2076 "interpret_floatarith (Mult a b) vs = (interpret_floatarith a vs) * (interpret_floatarith b vs)" |
2077 "interpret_floatarith (Inverse a) vs = inverse (interpret_floatarith a vs)" |
2078 "interpret_floatarith (Cos a) vs = cos (interpret_floatarith a vs)" |
2079 "interpret_floatarith (Arctan a) vs = arctan (interpret_floatarith a vs)" |
2080 "interpret_floatarith (Min a b) vs = min (interpret_floatarith a vs) (interpret_floatarith b vs)" |
2081 "interpret_floatarith (Max a b) vs = max (interpret_floatarith a vs) (interpret_floatarith b vs)" |
2082 "interpret_floatarith (Abs a) vs = abs (interpret_floatarith a vs)" |
2083 "interpret_floatarith Pi vs = pi" |
2084 "interpret_floatarith (Sqrt a) vs = sqrt (interpret_floatarith a vs)" |
2085 "interpret_floatarith (Exp a) vs = exp (interpret_floatarith a vs)" |
2086 "interpret_floatarith (Ln a) vs = ln (interpret_floatarith a vs)" |
2087 "interpret_floatarith (Power a n) vs = (interpret_floatarith a vs)^n" |
2088 "interpret_floatarith (Num f) vs = real f" |
2089 "interpret_floatarith (Atom n) vs = vs ! n"
2091 subsection "Implement approximation function"
2093 fun lift_bin' :: "(float * float) option \<Rightarrow> (float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> (float * float)) \<Rightarrow> (float * float) option" where
2094 "lift_bin' (Some (l1, u1)) (Some (l2, u2)) f = Some (f l1 u1 l2 u2)" |
2095 "lift_bin' a b f = None"
2097 fun lift_un :: "(float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> ((float option) * (float option))) \<Rightarrow> (float * float) option" where
2098 "lift_un (Some (l1, u1)) f = (case (f l1 u1) of (Some l, Some u) \<Rightarrow> Some (l, u)
2099 | t \<Rightarrow> None)" |
2100 "lift_un b f = None"
2102 fun lift_un' :: "(float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> (float * float)) \<Rightarrow> (float * float) option" where
2103 "lift_un' (Some (l1, u1)) f = Some (f l1 u1)" |
2104 "lift_un' b f = None"
2106 fun bounded_by :: "real list \<Rightarrow> (float * float) list \<Rightarrow> bool " where
2107 bounded_by_Cons: "bounded_by (v#vs) ((l, u)#bs) = ((real l \<le> v \<and> v \<le> real u) \<and> bounded_by vs bs)" |
2108 bounded_by_Nil: "bounded_by [] [] = True" |
2109 "bounded_by _ _ = False"
2111 lemma bounded_by: assumes "bounded_by vs bs" and "i < length bs"
2112 shows "real (fst (bs ! i)) \<le> vs ! i \<and> vs ! i \<le> real (snd (bs ! i))"
2113 using `bounded_by vs bs` and `i < length bs`
2114 proof (induct arbitrary: i rule: bounded_by.induct)
2115 fix v :: real and vs :: "real list" and l u :: float and bs :: "(float * float) list" and i :: nat
2116 assume hyp: "\<And>i. \<lbrakk>bounded_by vs bs; i < length bs\<rbrakk> \<Longrightarrow> real (fst (bs ! i)) \<le> vs ! i \<and> vs ! i \<le> real (snd (bs ! i))"
2117 assume bounded: "bounded_by (v # vs) ((l, u) # bs)" and length: "i < length ((l, u) # bs)"
2118 show "real (fst (((l, u) # bs) ! i)) \<le> (v # vs) ! i \<and> (v # vs) ! i \<le> real (snd (((l, u) # bs) ! i))"
2121 show ?thesis using bounded unfolding 0 nth_Cons_0 fst_conv snd_conv bounded_by.simps ..
2123 case (Suc i) with length have "i < length bs" by auto
2124 show ?thesis unfolding Suc nth_Cons_Suc bounded_by.simps
2125 using hyp[OF bounded[unfolded bounded_by.simps, THEN conjunct2] `i < length bs`] .
2129 fun approx approx' :: "nat \<Rightarrow> floatarith \<Rightarrow> (float * float) list \<Rightarrow> (float * float) option" where
2130 "approx' prec a bs = (case (approx prec a bs) of Some (l, u) \<Rightarrow> Some (round_down prec l, round_up prec u) | None \<Rightarrow> None)" |
2131 "approx prec (Add a b) bs = lift_bin' (approx' prec a bs) (approx' prec b bs) (\<lambda> l1 u1 l2 u2. (l1 + l2, u1 + u2))" |
2132 "approx prec (Minus a) bs = lift_un' (approx' prec a bs) (\<lambda> l u. (-u, -l))" |
2133 "approx prec (Mult a b) bs = lift_bin' (approx' prec a bs) (approx' prec b bs)
2134 (\<lambda> a1 a2 b1 b2. (float_nprt a1 * float_pprt b2 + float_nprt a2 * float_nprt b2 + float_pprt a1 * float_pprt b1 + float_pprt a2 * float_nprt b1,
2135 float_pprt a2 * float_pprt b2 + float_pprt a1 * float_nprt b2 + float_nprt a2 * float_pprt b1 + float_nprt a1 * float_nprt b1))" |
2136 "approx prec (Inverse a) bs = lift_un (approx' prec a bs) (\<lambda> l u. if (0 < l \<or> u < 0) then (Some (float_divl prec 1 u), Some (float_divr prec 1 l)) else (None, None))" |
2137 "approx prec (Cos a) bs = lift_un' (approx' prec a bs) (bnds_cos prec)" |
2138 "approx prec Pi bs = Some (lb_pi prec, ub_pi prec)" |
2139 "approx prec (Min a b) bs = lift_bin' (approx' prec a bs) (approx' prec b bs) (\<lambda> l1 u1 l2 u2. (min l1 l2, min u1 u2))" |
2140 "approx prec (Max a b) bs = lift_bin' (approx' prec a bs) (approx' prec b bs) (\<lambda> l1 u1 l2 u2. (max l1 l2, max u1 u2))" |
2141 "approx prec (Abs a) bs = lift_un' (approx' prec a bs) (\<lambda>l u. (if l < 0 \<and> 0 < u then 0 else min \<bar>l\<bar> \<bar>u\<bar>, max \<bar>l\<bar> \<bar>u\<bar>))" |
2142 "approx prec (Arctan a) bs = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_arctan prec l, ub_arctan prec u))" |
2143 "approx prec (Sqrt a) bs = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_sqrt prec l, ub_sqrt prec u))" |
2144 "approx prec (Exp a) bs = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_exp prec l, ub_exp prec u))" |
2145 "approx prec (Ln a) bs = lift_un (approx' prec a bs) (\<lambda> l u. (lb_ln prec l, ub_ln prec u))" |
2146 "approx prec (Power a n) bs = lift_un' (approx' prec a bs) (float_power_bnds n)" |
2147 "approx prec (Num f) bs = Some (f, f)" |
2148 "approx prec (Atom i) bs = (if i < length bs then Some (bs ! i) else None)"
2151 assumes lift_bin'_Some: "Some (l, u) = lift_bin' a b f"
2152 shows "\<exists> l1 u1 l2 u2. Some (l1, u1) = a \<and> Some (l2, u2) = b"
2154 case None hence "None = lift_bin' a b f" unfolding None lift_bin'.simps ..
2155 thus ?thesis using lift_bin'_Some by auto
2160 case None hence "None = lift_bin' a b f" unfolding None lift_bin'.simps ..
2161 thus ?thesis using lift_bin'_Some by auto
2164 obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
2165 obtain lb ub where b': "b' = (lb, ub)" by (cases b', auto)
2166 thus ?thesis unfolding `a = Some a'` `b = Some b'` a' b' by auto
2171 assumes lift_bin'_Some: "Some (l, u) = lift_bin' (g a) (g b) f"
2172 and Pa: "\<And>l u. Some (l, u) = g a \<Longrightarrow> P l u a" and Pb: "\<And>l u. Some (l, u) = g b \<Longrightarrow> P l u b"
2173 shows "\<exists> l1 u1 l2 u2. P l1 u1 a \<and> P l2 u2 b \<and> l = fst (f l1 u1 l2 u2) \<and> u = snd (f l1 u1 l2 u2)"
2176 where Sa: "Some (l1, u1) = g a" and Sb: "Some (l2, u2) = g b" using lift_bin'_ex[OF assms(1)] by auto
2177 have lu: "(l, u) = f l1 u1 l2 u2" using lift_bin'_Some[unfolded Sa[symmetric] Sb[symmetric] lift_bin'.simps] by auto
2178 have "l = fst (f l1 u1 l2 u2)" and "u = snd (f l1 u1 l2 u2)" unfolding lu[symmetric] by auto
2179 thus ?thesis using Pa[OF Sa] Pb[OF Sb] by auto
2182 lemma approx_approx':
2183 assumes Pa: "\<And>l u. Some (l, u) = approx prec a vs \<Longrightarrow> real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u"
2184 and approx': "Some (l, u) = approx' prec a vs"
2185 shows "real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u"
2187 obtain l' u' where S: "Some (l', u') = approx prec a vs"
2188 using approx' unfolding approx'.simps by (cases "approx prec a vs", auto)
2189 have l': "l = round_down prec l'" and u': "u = round_up prec u'"
2190 using approx' unfolding approx'.simps S[symmetric] by auto
2191 show ?thesis unfolding l' u'
2192 using order_trans[OF Pa[OF S, THEN conjunct2] round_up[of u']]
2193 using order_trans[OF round_down[of _ l'] Pa[OF S, THEN conjunct1]] by auto
2197 assumes lift_bin'_Some: "Some (l, u) = lift_bin' (approx' prec a bs) (approx' prec b bs) f"
2198 and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow> real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" (is "\<And>l u. _ = ?g a \<Longrightarrow> ?P l u a")
2199 and Pb: "\<And>l u. Some (l, u) = approx prec b bs \<Longrightarrow> real l \<le> interpret_floatarith b xs \<and> interpret_floatarith b xs \<le> real u"
2200 shows "\<exists> l1 u1 l2 u2. (real l1 \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u1) \<and>
2201 (real l2 \<le> interpret_floatarith b xs \<and> interpret_floatarith b xs \<le> real u2) \<and>
2202 l = fst (f l1 u1 l2 u2) \<and> u = snd (f l1 u1 l2 u2)"
2204 { fix l u assume "Some (l, u) = approx' prec a bs"
2205 with approx_approx'[of prec a bs, OF _ this] Pa
2206 have "real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" by auto } note Pa = this
2207 { fix l u assume "Some (l, u) = approx' prec b bs"
2208 with approx_approx'[of prec b bs, OF _ this] Pb
2209 have "real l \<le> interpret_floatarith b xs \<and> interpret_floatarith b xs \<le> real u" by auto } note Pb = this
2211 from lift_bin'_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_bin'_Some, OF Pa Pb]
2212 show ?thesis by auto
2216 assumes lift_un'_Some: "Some (l, u) = lift_un' a f"
2217 shows "\<exists> l u. Some (l, u) = a"
2219 case None hence "None = lift_un' a f" unfolding None lift_un'.simps ..
2220 thus ?thesis using lift_un'_Some by auto
2223 obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
2224 thus ?thesis unfolding `a = Some a'` a' by auto
2228 assumes lift_un'_Some: "Some (l, u) = lift_un' (g a) f"
2229 and Pa: "\<And>l u. Some (l, u) = g a \<Longrightarrow> P l u a"
2230 shows "\<exists> l1 u1. P l1 u1 a \<and> l = fst (f l1 u1) \<and> u = snd (f l1 u1)"
2232 obtain l1 u1 where Sa: "Some (l1, u1) = g a" using lift_un'_ex[OF assms(1)] by auto
2233 have lu: "(l, u) = f l1 u1" using lift_un'_Some[unfolded Sa[symmetric] lift_un'.simps] by auto
2234 have "l = fst (f l1 u1)" and "u = snd (f l1 u1)" unfolding lu[symmetric] by auto
2235 thus ?thesis using Pa[OF Sa] by auto
2239 assumes lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f"
2240 and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow> real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" (is "\<And>l u. _ = ?g a \<Longrightarrow> ?P l u a")
2241 shows "\<exists> l1 u1. (real l1 \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u1) \<and>
2242 l = fst (f l1 u1) \<and> u = snd (f l1 u1)"
2244 { fix l u assume "Some (l, u) = approx' prec a bs"
2245 with approx_approx'[of prec a bs, OF _ this] Pa
2246 have "real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" by auto } note Pa = this
2247 from lift_un'_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_un'_Some, OF Pa]
2248 show ?thesis by auto
2251 lemma lift_un'_bnds:
2252 assumes bnds: "\<forall> x lx ux. (l, u) = f lx ux \<and> x \<in> { real lx .. real ux } \<longrightarrow> real l \<le> f' x \<and> f' x \<le> real u"
2253 and lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f"
2254 and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow> real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u"
2255 shows "real l \<le> f' (interpret_floatarith a xs) \<and> f' (interpret_floatarith a xs) \<le> real u"
2257 from lift_un'[OF lift_un'_Some Pa]
2258 obtain l1 u1 where "real l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> real u1" and "l = fst (f l1 u1)" and "u = snd (f l1 u1)" by blast
2259 hence "(l, u) = f l1 u1" and "interpret_floatarith a xs \<in> {real l1 .. real u1}" by auto
2260 thus ?thesis using bnds by auto
2264 assumes lift_un_Some: "Some (l, u) = lift_un a f"
2265 shows "\<exists> l u. Some (l, u) = a"
2267 case None hence "None = lift_un a f" unfolding None lift_un.simps ..
2268 thus ?thesis using lift_un_Some by auto
2271 obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
2272 thus ?thesis unfolding `a = Some a'` a' by auto
2276 assumes lift_un_Some: "Some (l, u) = lift_un (g a) f"
2277 and Pa: "\<And>l u. Some (l, u) = g a \<Longrightarrow> P l u a"
2278 shows "\<exists> l1 u1. P l1 u1 a \<and> Some l = fst (f l1 u1) \<and> Some u = snd (f l1 u1)"
2280 obtain l1 u1 where Sa: "Some (l1, u1) = g a" using lift_un_ex[OF assms(1)] by auto
2281 have "fst (f l1 u1) \<noteq> None \<and> snd (f l1 u1) \<noteq> None"
2283 assume "\<not> (fst (f l1 u1) \<noteq> None \<and> snd (f l1 u1) \<noteq> None)"
2284 hence or: "fst (f l1 u1) = None \<or> snd (f l1 u1) = None" by auto
2285 hence "lift_un (g a) f = None"
2286 proof (cases "fst (f l1 u1) = None")
2288 then obtain b where b: "f l1 u1 = (None, b)" by (cases "f l1 u1", auto)
2289 thus ?thesis unfolding Sa[symmetric] lift_un.simps b by auto
2291 case False hence "snd (f l1 u1) = None" using or by auto
2292 with False obtain b where b: "f l1 u1 = (Some b, None)" by (cases "f l1 u1", auto)
2293 thus ?thesis unfolding Sa[symmetric] lift_un.simps b by auto
2295 thus False using lift_un_Some by auto
2297 then obtain a' b' where f: "f l1 u1 = (Some a', Some b')" by (cases "f l1 u1", auto)
2298 from lift_un_Some[unfolded Sa[symmetric] lift_un.simps f]
2299 have "Some l = fst (f l1 u1)" and "Some u = snd (f l1 u1)" unfolding f by auto
2300 thus ?thesis unfolding Sa[symmetric] lift_un.simps using Pa[OF Sa] by auto
2304 assumes lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f"
2305 and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow> real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" (is "\<And>l u. _ = ?g a \<Longrightarrow> ?P l u a")
2306 shows "\<exists> l1 u1. (real l1 \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u1) \<and>
2307 Some l = fst (f l1 u1) \<and> Some u = snd (f l1 u1)"
2309 { fix l u assume "Some (l, u) = approx' prec a bs"
2310 with approx_approx'[of prec a bs, OF _ this] Pa
2311 have "real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" by auto } note Pa = this
2312 from lift_un_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_un_Some, OF Pa]
2313 show ?thesis by auto
2317 assumes bnds: "\<forall> x lx ux. (Some l, Some u) = f lx ux \<and> x \<in> { real lx .. real ux } \<longrightarrow> real l \<le> f' x \<and> f' x \<le> real u"
2318 and lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f"
2319 and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow> real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u"
2320 shows "real l \<le> f' (interpret_floatarith a xs) \<and> f' (interpret_floatarith a xs) \<le> real u"
2322 from lift_un[OF lift_un_Some Pa]
2323 obtain l1 u1 where "real l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> real u1" and "Some l = fst (f l1 u1)" and "Some u = snd (f l1 u1)" by blast
2324 hence "(Some l, Some u) = f l1 u1" and "interpret_floatarith a xs \<in> {real l1 .. real u1}" by auto
2325 thus ?thesis using bnds by auto
2329 assumes "bounded_by xs vs"
2330 and "Some (l, u) = approx prec arith vs" (is "_ = ?g arith")
2331 shows "real l \<le> interpret_floatarith arith xs \<and> interpret_floatarith arith xs \<le> real u" (is "?P l u arith")
2332 using `Some (l, u) = approx prec arith vs`
2333 proof (induct arith arbitrary: l u x)
2335 from lift_bin'[OF Add.prems[unfolded approx.simps]] Add.hyps
2336 obtain l1 u1 l2 u2 where "l = l1 + l2" and "u = u1 + u2"
2337 "real l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> real u1"
2338 "real l2 \<le> interpret_floatarith b xs" and "interpret_floatarith b xs \<le> real u2" unfolding fst_conv snd_conv by blast
2339 thus ?case unfolding interpret_floatarith.simps by auto
2342 from lift_un'[OF Minus.prems[unfolded approx.simps]] Minus.hyps
2343 obtain l1 u1 where "l = -u1" and "u = -l1"
2344 "real l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> real u1" unfolding fst_conv snd_conv by blast
2345 thus ?case unfolding interpret_floatarith.simps using real_of_float_minus by auto
2348 from lift_bin'[OF Mult.prems[unfolded approx.simps]] Mult.hyps
2350 where l: "l = float_nprt l1 * float_pprt u2 + float_nprt u1 * float_nprt u2 + float_pprt l1 * float_pprt l2 + float_pprt u1 * float_nprt l2"
2351 and u: "u = float_pprt u1 * float_pprt u2 + float_pprt l1 * float_nprt u2 + float_nprt u1 * float_pprt l2 + float_nprt l1 * float_nprt l2"
2352 and "real l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> real u1"
2353 and "real l2 \<le> interpret_floatarith b xs" and "interpret_floatarith b xs \<le> real u2" unfolding fst_conv snd_conv by blast
2354 thus ?case unfolding interpret_floatarith.simps l u real_of_float_add real_of_float_mult real_of_float_nprt real_of_float_pprt
2355 using mult_le_prts mult_ge_prts by auto
2358 from lift_un[OF Inverse.prems[unfolded approx.simps], unfolded if_distrib[of fst] if_distrib[of snd] fst_conv snd_conv] Inverse.hyps
2359 obtain l1 u1 where l': "Some l = (if 0 < l1 \<or> u1 < 0 then Some (float_divl prec 1 u1) else None)"
2360 and u': "Some u = (if 0 < l1 \<or> u1 < 0 then Some (float_divr prec 1 l1) else None)"
2361 and l1: "real l1 \<le> interpret_floatarith a xs" and u1: "interpret_floatarith a xs \<le> real u1" by blast
2362 have either: "0 < l1 \<or> u1 < 0" proof (rule ccontr) assume P: "\<not> (0 < l1 \<or> u1 < 0)" show False using l' unfolding if_not_P[OF P] by auto qed
2363 moreover have l1_le_u1: "real l1 \<le> real u1" using l1 u1 by auto
2364 ultimately have "real l1 \<noteq> 0" and "real u1 \<noteq> 0" unfolding less_float_def by auto
2366 have inv: "inverse (real u1) \<le> inverse (interpret_floatarith a xs)
2367 \<and> inverse (interpret_floatarith a xs) \<le> inverse (real l1)"
2368 proof (cases "0 < l1")
2369 case True hence "0 < real u1" and "0 < real l1" "0 < interpret_floatarith a xs"
2370 unfolding less_float_def using l1_le_u1 l1 by auto
2372 unfolding inverse_le_iff_le[OF `0 < real u1` `0 < interpret_floatarith a xs`]
2373 inverse_le_iff_le[OF `0 < interpret_floatarith a xs` `0 < real l1`]
2376 case False hence "u1 < 0" using either by blast
2377 hence "real u1 < 0" and "real l1 < 0" "interpret_floatarith a xs < 0"
2378 unfolding less_float_def using l1_le_u1 u1 by auto
2380 unfolding inverse_le_iff_le_neg[OF `real u1 < 0` `interpret_floatarith a xs < 0`]
2381 inverse_le_iff_le_neg[OF `interpret_floatarith a xs < 0` `real l1 < 0`]
2385 from l' have "l = float_divl prec 1 u1" by (cases "0 < l1 \<or> u1 < 0", auto)
2386 hence "real l \<le> inverse (real u1)" unfolding nonzero_inverse_eq_divide[OF `real u1 \<noteq> 0`] using float_divl[of prec 1 u1] by auto
2387 also have "\<dots> \<le> inverse (interpret_floatarith a xs)" using inv by auto
2388 finally have "real l \<le> inverse (interpret_floatarith a xs)" .
2390 from u' have "u = float_divr prec 1 l1" by (cases "0 < l1 \<or> u1 < 0", auto)
2391 hence "inverse (real l1) \<le> real u" unfolding nonzero_inverse_eq_divide[OF `real l1 \<noteq> 0`] using float_divr[of 1 l1 prec] by auto
2392 hence "inverse (interpret_floatarith a xs) \<le> real u" by (rule order_trans[OF inv[THEN conjunct2]])
2393 ultimately show ?case unfolding interpret_floatarith.simps using l1 u1 by auto
2396 from lift_un'[OF Abs.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Abs.hyps
2397 obtain l1 u1 where l': "l = (if l1 < 0 \<and> 0 < u1 then 0 else min \<bar>l1\<bar> \<bar>u1\<bar>)" and u': "u = max \<bar>l1\<bar> \<bar>u1\<bar>"
2398 and l1: "real l1 \<le> interpret_floatarith x xs" and u1: "interpret_floatarith x xs \<le> real u1" by blast
2399 thus ?case unfolding l' u' by (cases "l1 < 0 \<and> 0 < u1", auto simp add: real_of_float_min real_of_float_max real_of_float_abs less_float_def)
2402 from lift_bin'[OF Min.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Min.hyps
2403 obtain l1 u1 l2 u2 where l': "l = min l1 l2" and u': "u = min u1 u2"
2404 and l1: "real l1 \<le> interpret_floatarith a xs" and u1: "interpret_floatarith a xs \<le> real u1"
2405 and l1: "real l2 \<le> interpret_floatarith b xs" and u1: "interpret_floatarith b xs \<le> real u2" by blast
2406 thus ?case unfolding l' u' by (auto simp add: real_of_float_min)
2409 from lift_bin'[OF Max.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Max.hyps
2410 obtain l1 u1 l2 u2 where l': "l = max l1 l2" and u': "u = max u1 u2"
2411 and l1: "real l1 \<le> interpret_floatarith a xs" and u1: "interpret_floatarith a xs \<le> real u1"
2412 and l1: "real l2 \<le> interpret_floatarith b xs" and u1: "interpret_floatarith b xs \<le> real u2" by blast
2413 thus ?case unfolding l' u' by (auto simp add: real_of_float_max)
2414 next case (Cos a) with lift_un'_bnds[OF bnds_cos] show ?case by auto
2415 next case (Arctan a) with lift_un'_bnds[OF bnds_arctan] show ?case by auto
2416 next case Pi with pi_boundaries show ?case by auto
2417 next case (Sqrt a) with lift_un'_bnds[OF bnds_sqrt] show ?case by auto
2418 next case (Exp a) with lift_un'_bnds[OF bnds_exp] show ?case by auto
2419 next case (Ln a) with lift_un_bnds[OF bnds_ln] show ?case by auto
2420 next case (Power a n) with lift_un'_bnds[OF bnds_power] show ?case by auto
2421 next case (Num f) thus ?case by auto
2425 proof (cases "n < length vs")
2427 with Atom have "vs ! n = (l, u)" by auto
2428 thus ?thesis using bounded_by[OF assms(1) True] by auto
2430 case False thus ?thesis using Atom by auto
2434 datatype inequality = Less floatarith floatarith
2435 | LessEqual floatarith floatarith
2437 fun interpret_inequality :: "inequality \<Rightarrow> real list \<Rightarrow> bool" where
2438 "interpret_inequality (Less a b) vs = (interpret_floatarith a vs < interpret_floatarith b vs)" |
2439 "interpret_inequality (LessEqual a b) vs = (interpret_floatarith a vs \<le> interpret_floatarith b vs)"
2441 fun approx_inequality :: "nat \<Rightarrow> inequality \<Rightarrow> (float * float) list \<Rightarrow> bool" where
2442 "approx_inequality prec (Less a b) bs = (case (approx prec a bs, approx prec b bs) of (Some (l, u), Some (l', u')) \<Rightarrow> u < l' | _ \<Rightarrow> False)" |
2443 "approx_inequality prec (LessEqual a b) bs = (case (approx prec a bs, approx prec b bs) of (Some (l, u), Some (l', u')) \<Rightarrow> u \<le> l' | _ \<Rightarrow> False)"
2445 lemma approx_inequality: fixes m :: nat assumes "bounded_by vs bs" and "approx_inequality prec eq bs"
2446 shows "interpret_inequality eq vs"
2450 proof (cases "\<exists> u l u' l'. approx prec a bs = Some (l, u) \<and>
2451 approx prec b bs = Some (l', u')")
2453 then obtain l u l' u' where a_approx: "approx prec a bs = Some (l, u)"
2454 and b_approx: "approx prec b bs = Some (l', u') " by auto
2455 with `approx_inequality prec eq bs` have "real u < real l'"
2456 unfolding Less approx_inequality.simps less_float_def by auto
2457 moreover from a_approx[symmetric] and b_approx[symmetric] and `bounded_by vs bs`
2458 have "interpret_floatarith a vs \<le> real u" and "real l' \<le> interpret_floatarith b vs"
2459 using approx by auto
2460 ultimately show ?thesis unfolding interpret_inequality.simps Less by auto
2463 hence "approx prec a bs = None \<or> approx prec b bs = None"
2464 unfolding not_Some_eq[symmetric] by auto
2465 hence "\<not> approx_inequality prec eq bs" unfolding Less approx_inequality.simps
2466 by (cases "approx prec a bs = None", auto)
2467 thus ?thesis using assms by auto
2470 case (LessEqual a b)
2472 proof (cases "\<exists> u l u' l'. approx prec a bs = Some (l, u) \<and>
2473 approx prec b bs = Some (l', u')")
2475 then obtain l u l' u' where a_approx: "approx prec a bs = Some (l, u)"
2476 and b_approx: "approx prec b bs = Some (l', u') " by auto
2477 with `approx_inequality prec eq bs` have "real u \<le> real l'"
2478 unfolding LessEqual approx_inequality.simps le_float_def by auto
2479 moreover from a_approx[symmetric] and b_approx[symmetric] and `bounded_by vs bs`
2480 have "interpret_floatarith a vs \<le> real u" and "real l' \<le> interpret_floatarith b vs"
2481 using approx by auto
2482 ultimately show ?thesis unfolding interpret_inequality.simps LessEqual by auto
2485 hence "approx prec a bs = None \<or> approx prec b bs = None"
2486 unfolding not_Some_eq[symmetric] by auto
2487 hence "\<not> approx_inequality prec eq bs" unfolding LessEqual approx_inequality.simps
2488 by (cases "approx prec a bs = None", auto)
2489 thus ?thesis using assms by auto
2493 lemma interpret_floatarith_divide: "interpret_floatarith (Mult a (Inverse b)) vs = (interpret_floatarith a vs) / (interpret_floatarith b vs)"
2494 unfolding real_divide_def interpret_floatarith.simps ..
2496 lemma interpret_floatarith_diff: "interpret_floatarith (Add a (Minus b)) vs = (interpret_floatarith a vs) - (interpret_floatarith b vs)"
2497 unfolding real_diff_def interpret_floatarith.simps ..
2499 lemma interpret_floatarith_sin: "interpret_floatarith (Cos (Add (Mult Pi (Num (Float 1 -1))) (Minus a))) vs =
2500 sin (interpret_floatarith a vs)"
2501 unfolding sin_cos_eq interpret_floatarith.simps
2502 interpret_floatarith_divide interpret_floatarith_diff real_diff_def
2505 lemma interpret_floatarith_tan:
2506 "interpret_floatarith (Mult (Cos (Add (Mult Pi (Num (Float 1 -1))) (Minus a))) (Inverse (Cos a))) vs =
2507 tan (interpret_floatarith a vs)"
2508 unfolding interpret_floatarith.simps(3,4) interpret_floatarith_sin tan_def real_divide_def
2511 lemma interpret_floatarith_powr: "interpret_floatarith (Exp (Mult b (Ln a))) vs = (interpret_floatarith a vs) powr (interpret_floatarith b vs)"
2512 unfolding powr_def interpret_floatarith.simps ..
2514 lemma interpret_floatarith_log: "interpret_floatarith ((Mult (Ln x) (Inverse (Ln b)))) vs = log (interpret_floatarith b vs) (interpret_floatarith x vs)"
2515 unfolding log_def interpret_floatarith.simps real_divide_def ..
2517 lemma interpret_floatarith_num:
2518 shows "interpret_floatarith (Num (Float 0 0)) vs = 0"
2519 and "interpret_floatarith (Num (Float 1 0)) vs = 1"
2520 and "interpret_floatarith (Num (Float (number_of a) 0)) vs = number_of a" by auto
2522 subsection {* Implement proof method \texttt{approximation} *}
2524 lemma bounded_divl: assumes "real a / real b \<le> x" shows "real (float_divl p a b) \<le> x" by (rule order_trans[OF _ assms], rule float_divl)
2525 lemma bounded_divr: assumes "x \<le> real a / real b" shows "x \<le> real (float_divr p a b)" by (rule order_trans[OF assms _], rule float_divr)
2526 lemma bounded_num: shows "real (Float 5 1) = 10" and "real (Float 0 0) = 0" and "real (Float 1 0) = 1" and "real (Float (number_of n) 0) = (number_of n)"
2527 and "0 * pow2 e = real (Float 0 e)" and "1 * pow2 e = real (Float 1 e)" and "number_of m * pow2 e = real (Float (number_of m) e)"
2528 and "real (Float (number_of A) (int B)) = (number_of A) * 2^B"
2529 and "real (Float 1 (int B)) = 2^B"
2530 and "real (Float (number_of A) (- int B)) = (number_of A) / 2^B"
2531 and "real (Float 1 (- int B)) = 1 / 2^B"
2532 by (auto simp add: real_of_float_simp pow2_def real_divide_def)
2534 lemmas bounded_by_equations = bounded_by_Cons bounded_by_Nil float_power bounded_divl bounded_divr bounded_num HOL.simp_thms
2535 lemmas interpret_inequality_equations = interpret_inequality.simps interpret_floatarith.simps interpret_floatarith_num
2536 interpret_floatarith_divide interpret_floatarith_diff interpret_floatarith_tan interpret_floatarith_powr interpret_floatarith_log
2537 interpret_floatarith_sin
2540 structure Float_Arith =
2543 @{code_datatype float = Float}
2544 @{code_datatype floatarith = Add | Minus | Mult | Inverse | Cos | Arctan
2545 | Abs | Max | Min | Pi | Sqrt | Exp | Ln | Power | Atom | Num }
2546 @{code_datatype inequality = Less | LessEqual }
2548 val approx_inequality = @{code approx_inequality}
2553 code_reserved Eval Float_Arith
2555 code_type float (Eval "Float'_Arith.float")
2556 code_const Float (Eval "Float'_Arith.Float/ (_,/ _)")
2558 code_type floatarith (Eval "Float'_Arith.floatarith")
2559 code_const Add and Minus and Mult and Inverse and Cos and Arctan and Abs and Max and Min and
2560 Pi and Sqrt and Exp and Ln and Power and Atom and Num
2561 (Eval "Float'_Arith.Add/ (_,/ _)" and "Float'_Arith.Minus" and "Float'_Arith.Mult/ (_,/ _)" and
2562 "Float'_Arith.Inverse" and "Float'_Arith.Cos" and
2563 "Float'_Arith.Arctan" and "Float'_Arith.Abs" and "Float'_Arith.Max/ (_,/ _)" and
2564 "Float'_Arith.Min/ (_,/ _)" and "Float'_Arith.Pi" and "Float'_Arith.Sqrt" and
2565 "Float'_Arith.Exp" and "Float'_Arith.Ln" and "Float'_Arith.Power/ (_,/ _)" and
2566 "Float'_Arith.Atom" and "Float'_Arith.Num")
2568 code_type inequality (Eval "Float'_Arith.inequality")
2569 code_const Less and LessEqual (Eval "Float'_Arith.Less/ (_,/ _)" and "Float'_Arith.LessEqual/ (_,/ _)")
2571 code_const approx_inequality (Eval "Float'_Arith.approx'_inequality")
2574 val ineq_equations = PureThy.get_thms @{theory} "interpret_inequality_equations";
2575 val bounded_by_equations = PureThy.get_thms @{theory} "bounded_by_equations";
2576 val bounded_by_simpset = (HOL_basic_ss addsimps bounded_by_equations)
2578 fun reify_ineq ctxt i = (fn st =>
2580 val to = HOLogic.dest_Trueprop (Logic.strip_imp_concl (List.nth (prems_of st, i - 1)))
2581 in (Reflection.genreify_tac ctxt ineq_equations (SOME to) i) st
2584 fun rule_ineq ctxt prec i thm = let
2585 fun conv_num typ = HOLogic.dest_number #> snd #> HOLogic.mk_number typ
2586 val to_natc = conv_num @{typ "nat"} #> Thm.cterm_of (ProofContext.theory_of ctxt)
2587 val to_nat = conv_num @{typ "nat"}
2588 val to_int = conv_num @{typ "int"}
2589 fun int_to_float A = @{term "Float"} $ to_int A $ @{term "0 :: int"}
2591 val prec' = to_nat prec
2593 fun bot_float (Const (@{const_name "times"}, _) $ mantisse $ (Const (@{const_name "pow2"}, _) $ exp))
2594 = @{term "Float"} $ to_int mantisse $ to_int exp
2595 | bot_float (Const (@{const_name "divide"}, _) $ mantisse $ (@{term "power 2 :: nat \<Rightarrow> real"} $ exp))
2596 = @{term "Float"} $ to_int mantisse $ (@{term "uminus :: int \<Rightarrow> int"} $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp))
2597 | bot_float (Const (@{const_name "times"}, _) $ mantisse $ (@{term "power 2 :: nat \<Rightarrow> real"} $ exp))
2598 = @{term "Float"} $ to_int mantisse $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp)
2599 | bot_float (Const (@{const_name "divide"}, _) $ A $ (@{term "power 10 :: nat \<Rightarrow> real"} $ exp))
2600 = @{term "float_divl"} $ prec' $ (int_to_float A) $ (@{term "power (Float 5 1)"} $ to_nat exp)
2601 | bot_float (Const (@{const_name "divide"}, _) $ A $ B)
2602 = @{term "float_divl"} $ prec' $ int_to_float A $ int_to_float B
2603 | bot_float (@{term "power 2 :: nat \<Rightarrow> real"} $ exp)
2604 = @{term "Float 1"} $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp)
2605 | bot_float A = int_to_float A
2607 fun top_float (Const (@{const_name "times"}, _) $ mantisse $ (Const (@{const_name "pow2"}, _) $ exp))
2608 = @{term "Float"} $ to_int mantisse $ to_int exp
2609 | top_float (Const (@{const_name "divide"}, _) $ mantisse $ (@{term "power 2 :: nat \<Rightarrow> real"} $ exp))
2610 = @{term "Float"} $ to_int mantisse $ (@{term "uminus :: int \<Rightarrow> int"} $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp))
2611 | top_float (Const (@{const_name "times"}, _) $ mantisse $ (@{term "power 2 :: nat \<Rightarrow> real"} $ exp))
2612 = @{term "Float"} $ to_int mantisse $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp)
2613 | top_float (Const (@{const_name "divide"}, _) $ A $ (@{term "power 10 :: nat \<Rightarrow> real"} $ exp))
2614 = @{term "float_divr"} $ prec' $ (int_to_float A) $ (@{term "power (Float 5 1)"} $ to_nat exp)
2615 | top_float (Const (@{const_name "divide"}, _) $ A $ B)
2616 = @{term "float_divr"} $ prec' $ int_to_float A $ int_to_float B
2617 | top_float (@{term "power 2 :: nat \<Rightarrow> real"} $ exp)
2618 = @{term "Float 1"} $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp)
2619 | top_float A = int_to_float A
2621 val goal' : term = List.nth (prems_of thm, i - 1)
2623 fun lift_bnd (t as (Const (@{const_name "op &"}, _) $
2624 (Const (@{const_name "less_eq"}, _) $
2625 bottom $ (Free (name, _))) $
2626 (Const (@{const_name "less_eq"}, _) $ _ $ top)))
2627 = ((name, HOLogic.mk_prod (bot_float bottom, top_float top))
2628 handle TERM (txt, ts) => raise TERM ("Invalid bound number format: " ^
2629 (Syntax.string_of_term ctxt t), [t]))
2630 | lift_bnd t = raise TERM ("Premisse needs format '<num> <= <var> & <var> <= <num>', but found " ^
2631 (Syntax.string_of_term ctxt t), [t])
2632 val bound_eqs = map (HOLogic.dest_Trueprop #> lift_bnd) (Logic.strip_imp_prems goal')
2634 fun lift_var (Free (varname, _)) = (case AList.lookup (op =) bound_eqs varname of
2636 | NONE => raise TERM ("No bound equations found for " ^ varname, []))
2637 | lift_var t = raise TERM ("Can not convert expression " ^
2638 (Syntax.string_of_term ctxt t), [t])
2640 val _ $ vs = HOLogic.dest_Trueprop (Logic.strip_imp_concl goal')
2642 val bs = (HOLogic.dest_list #> map lift_var #> HOLogic.mk_list @{typ "float * float"}) vs
2643 val map = [(@{cpat "?prec::nat"}, to_natc prec),
2644 (@{cpat "?bs::(float * float) list"}, Thm.cterm_of (ProofContext.theory_of ctxt) bs)]
2645 in rtac (Thm.instantiate ([], map) @{thm "approx_inequality"}) i thm end
2647 val eval_tac = CSUBGOAL (fn (ct, i) => rtac (eval_oracle ct) i)
2649 fun gen_eval_tac conv ctxt = CONVERSION (Conv.params_conv (~1) (K (Conv.concl_conv (~1) conv)) ctxt)
2654 method_setup approximation = {*
2656 (fn prec => fn ctxt =>
2657 SIMPLE_METHOD' (fn i =>
2658 (DETERM (reify_ineq ctxt i)
2659 THEN rule_ineq ctxt prec i
2660 THEN Simplifier.asm_full_simp_tac bounded_by_simpset i
2661 THEN (TRY (filter_prems_tac (fn t => false) i))
2662 THEN (gen_eval_tac eval_oracle ctxt) i)))
2663 *} "real number approximation"