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 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 mult_assoc[symmetric]
113 unfolding mult_commute unfolding 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 mult_assoc[symmetric] setsum_left_distrib[symmetric]
439 unfolding mult_commute[where 'a=real] 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 mult_assoc[symmetric] setsum_left_distrib[symmetric]
448 unfolding mult_commute[where 'a=real] 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 mult_1_left .
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 mult_1_left 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 mult_1_left 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 mult_1_left .
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 mult_assoc[symmetric] setsum_left_distrib[symmetric]
916 unfolding mult_commute[where 'a=real]
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 left_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: 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 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 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 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 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 diff_def mult_minus_left mult_1_left)
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 diff_def mult_minus_left mult_1_left)
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!: 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!: 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 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] real_of_nat_power 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 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 mult_assoc[symmetric] real_of_float_mult setsum_left_distrib[symmetric] unfolding 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 mult_assoc[symmetric] real_of_float_mult setsum_left_distrib[symmetric] unfolding 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 min_max.sup_absorb1)
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" where
2073 "interpret_floatarith (Add a b) vs = (interpret_floatarith a vs) + (interpret_floatarith b vs)" |
2074 "interpret_floatarith (Minus a) vs = - (interpret_floatarith a vs)" |
2075 "interpret_floatarith (Mult a b) vs = (interpret_floatarith a vs) * (interpret_floatarith b vs)" |
2076 "interpret_floatarith (Inverse a) vs = inverse (interpret_floatarith a vs)" |
2077 "interpret_floatarith (Cos a) vs = cos (interpret_floatarith a vs)" |
2078 "interpret_floatarith (Arctan a) vs = arctan (interpret_floatarith a vs)" |
2079 "interpret_floatarith (Min a b) vs = min (interpret_floatarith a vs) (interpret_floatarith b vs)" |
2080 "interpret_floatarith (Max a b) vs = max (interpret_floatarith a vs) (interpret_floatarith b vs)" |
2081 "interpret_floatarith (Abs a) vs = abs (interpret_floatarith a vs)" |
2082 "interpret_floatarith Pi vs = pi" |
2083 "interpret_floatarith (Sqrt a) vs = sqrt (interpret_floatarith a vs)" |
2084 "interpret_floatarith (Exp a) vs = exp (interpret_floatarith a vs)" |
2085 "interpret_floatarith (Ln a) vs = ln (interpret_floatarith a vs)" |
2086 "interpret_floatarith (Power a n) vs = (interpret_floatarith a vs)^n" |
2087 "interpret_floatarith (Num f) vs = real f" |
2088 "interpret_floatarith (Var n) vs = vs ! n"
2090 lemma interpret_floatarith_divide: "interpret_floatarith (Mult a (Inverse b)) vs = (interpret_floatarith a vs) / (interpret_floatarith b vs)"
2091 unfolding divide_inverse interpret_floatarith.simps ..
2093 lemma interpret_floatarith_diff: "interpret_floatarith (Add a (Minus b)) vs = (interpret_floatarith a vs) - (interpret_floatarith b vs)"
2094 unfolding diff_def interpret_floatarith.simps ..
2096 lemma interpret_floatarith_sin: "interpret_floatarith (Cos (Add (Mult Pi (Num (Float 1 -1))) (Minus a))) vs =
2097 sin (interpret_floatarith a vs)"
2098 unfolding sin_cos_eq interpret_floatarith.simps
2099 interpret_floatarith_divide interpret_floatarith_diff diff_def
2102 lemma interpret_floatarith_tan:
2103 "interpret_floatarith (Mult (Cos (Add (Mult Pi (Num (Float 1 -1))) (Minus a))) (Inverse (Cos a))) vs =
2104 tan (interpret_floatarith a vs)"
2105 unfolding interpret_floatarith.simps(3,4) interpret_floatarith_sin tan_def divide_inverse
2108 lemma interpret_floatarith_powr: "interpret_floatarith (Exp (Mult b (Ln a))) vs = (interpret_floatarith a vs) powr (interpret_floatarith b vs)"
2109 unfolding powr_def interpret_floatarith.simps ..
2111 lemma interpret_floatarith_log: "interpret_floatarith ((Mult (Ln x) (Inverse (Ln b)))) vs = log (interpret_floatarith b vs) (interpret_floatarith x vs)"
2112 unfolding log_def interpret_floatarith.simps divide_inverse ..
2114 lemma interpret_floatarith_num:
2115 shows "interpret_floatarith (Num (Float 0 0)) vs = 0"
2116 and "interpret_floatarith (Num (Float 1 0)) vs = 1"
2117 and "interpret_floatarith (Num (Float (number_of a) 0)) vs = number_of a" by auto
2119 subsection "Implement approximation function"
2121 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
2122 "lift_bin' (Some (l1, u1)) (Some (l2, u2)) f = Some (f l1 u1 l2 u2)" |
2123 "lift_bin' a b f = None"
2125 fun lift_un :: "(float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> ((float option) * (float option))) \<Rightarrow> (float * float) option" where
2126 "lift_un (Some (l1, u1)) f = (case (f l1 u1) of (Some l, Some u) \<Rightarrow> Some (l, u)
2127 | t \<Rightarrow> None)" |
2128 "lift_un b f = None"
2130 fun lift_un' :: "(float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> (float * float)) \<Rightarrow> (float * float) option" where
2131 "lift_un' (Some (l1, u1)) f = Some (f l1 u1)" |
2132 "lift_un' b f = None"
2135 "bounded_by xs vs \<longleftrightarrow>
2136 (\<forall> i < length vs. case vs ! i of None \<Rightarrow> True
2137 | Some (l, u) \<Rightarrow> xs ! i \<in> { real l .. real u })"
2140 assumes "bounded_by xs vs"
2141 shows "\<And> i. i < length vs \<Longrightarrow> case vs ! i of None \<Rightarrow> True
2142 | Some (l, u) \<Rightarrow> xs ! i \<in> { real l .. real u }"
2143 using assms bounded_by_def by blast
2145 lemma bounded_by_update:
2146 assumes "bounded_by xs vs"
2147 and bnd: "xs ! i \<in> { real l .. real u }"
2148 shows "bounded_by xs (vs[i := Some (l,u)])"
2151 let ?vs = "vs[i := Some (l,u)]"
2152 assume "j < length ?vs" hence [simp]: "j < length vs" by simp
2153 have "case ?vs ! j of None \<Rightarrow> True | Some (l, u) \<Rightarrow> xs ! j \<in> { real l .. real u }"
2154 proof (cases "?vs ! j")
2157 proof (cases "i = j")
2159 thus ?thesis using `?vs ! j = Some b` and bnd by auto
2162 thus ?thesis using `bounded_by xs vs` unfolding bounded_by_def by auto
2165 thus ?thesis unfolding bounded_by_def by auto
2168 lemma bounded_by_None:
2169 shows "bounded_by xs (replicate (length xs) None)"
2170 unfolding bounded_by_def by auto
2172 fun approx approx' :: "nat \<Rightarrow> floatarith \<Rightarrow> (float * float) option list \<Rightarrow> (float * float) option" where
2173 "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)" |
2174 "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))" |
2175 "approx prec (Minus a) bs = lift_un' (approx' prec a bs) (\<lambda> l u. (-u, -l))" |
2176 "approx prec (Mult a b) bs = lift_bin' (approx' prec a bs) (approx' prec b bs)
2177 (\<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,
2178 float_pprt a2 * float_pprt b2 + float_pprt a1 * float_nprt b2 + float_nprt a2 * float_pprt b1 + float_nprt a1 * float_nprt b1))" |
2179 "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))" |
2180 "approx prec (Cos a) bs = lift_un' (approx' prec a bs) (bnds_cos prec)" |
2181 "approx prec Pi bs = Some (lb_pi prec, ub_pi prec)" |
2182 "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))" |
2183 "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))" |
2184 "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>))" |
2185 "approx prec (Arctan a) bs = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_arctan prec l, ub_arctan prec u))" |
2186 "approx prec (Sqrt a) bs = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_sqrt prec l, ub_sqrt prec u))" |
2187 "approx prec (Exp a) bs = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_exp prec l, ub_exp prec u))" |
2188 "approx prec (Ln a) bs = lift_un (approx' prec a bs) (\<lambda> l u. (lb_ln prec l, ub_ln prec u))" |
2189 "approx prec (Power a n) bs = lift_un' (approx' prec a bs) (float_power_bnds n)" |
2190 "approx prec (Num f) bs = Some (f, f)" |
2191 "approx prec (Var i) bs = (if i < length bs then bs ! i else None)"
2194 assumes lift_bin'_Some: "Some (l, u) = lift_bin' a b f"
2195 shows "\<exists> l1 u1 l2 u2. Some (l1, u1) = a \<and> Some (l2, u2) = b"
2197 case None hence "None = lift_bin' a b f" unfolding None lift_bin'.simps ..
2198 thus ?thesis using lift_bin'_Some by auto
2203 case None hence "None = lift_bin' a b f" unfolding None lift_bin'.simps ..
2204 thus ?thesis using lift_bin'_Some by auto
2207 obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
2208 obtain lb ub where b': "b' = (lb, ub)" by (cases b', auto)
2209 thus ?thesis unfolding `a = Some a'` `b = Some b'` a' b' by auto
2214 assumes lift_bin'_Some: "Some (l, u) = lift_bin' (g a) (g b) f"
2215 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"
2216 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)"
2219 where Sa: "Some (l1, u1) = g a" and Sb: "Some (l2, u2) = g b" using lift_bin'_ex[OF assms(1)] by auto
2220 have lu: "(l, u) = f l1 u1 l2 u2" using lift_bin'_Some[unfolded Sa[symmetric] Sb[symmetric] lift_bin'.simps] by auto
2221 have "l = fst (f l1 u1 l2 u2)" and "u = snd (f l1 u1 l2 u2)" unfolding lu[symmetric] by auto
2222 thus ?thesis using Pa[OF Sa] Pb[OF Sb] by auto
2225 lemma approx_approx':
2226 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"
2227 and approx': "Some (l, u) = approx' prec a vs"
2228 shows "real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u"
2230 obtain l' u' where S: "Some (l', u') = approx prec a vs"
2231 using approx' unfolding approx'.simps by (cases "approx prec a vs", auto)
2232 have l': "l = round_down prec l'" and u': "u = round_up prec u'"
2233 using approx' unfolding approx'.simps S[symmetric] by auto
2234 show ?thesis unfolding l' u'
2235 using order_trans[OF Pa[OF S, THEN conjunct2] round_up[of u']]
2236 using order_trans[OF round_down[of _ l'] Pa[OF S, THEN conjunct1]] by auto
2240 assumes lift_bin'_Some: "Some (l, u) = lift_bin' (approx' prec a bs) (approx' prec b bs) f"
2241 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")
2242 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"
2243 shows "\<exists> l1 u1 l2 u2. (real l1 \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u1) \<and>
2244 (real l2 \<le> interpret_floatarith b xs \<and> interpret_floatarith b xs \<le> real u2) \<and>
2245 l = fst (f l1 u1 l2 u2) \<and> u = snd (f l1 u1 l2 u2)"
2247 { fix l u assume "Some (l, u) = approx' prec a bs"
2248 with approx_approx'[of prec a bs, OF _ this] Pa
2249 have "real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" by auto } note Pa = this
2250 { fix l u assume "Some (l, u) = approx' prec b bs"
2251 with approx_approx'[of prec b bs, OF _ this] Pb
2252 have "real l \<le> interpret_floatarith b xs \<and> interpret_floatarith b xs \<le> real u" by auto } note Pb = this
2254 from lift_bin'_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_bin'_Some, OF Pa Pb]
2255 show ?thesis by auto
2259 assumes lift_un'_Some: "Some (l, u) = lift_un' a f"
2260 shows "\<exists> l u. Some (l, u) = a"
2262 case None hence "None = lift_un' a f" unfolding None lift_un'.simps ..
2263 thus ?thesis using lift_un'_Some by auto
2266 obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
2267 thus ?thesis unfolding `a = Some a'` a' by auto
2271 assumes lift_un'_Some: "Some (l, u) = lift_un' (g a) f"
2272 and Pa: "\<And>l u. Some (l, u) = g a \<Longrightarrow> P l u a"
2273 shows "\<exists> l1 u1. P l1 u1 a \<and> l = fst (f l1 u1) \<and> u = snd (f l1 u1)"
2275 obtain l1 u1 where Sa: "Some (l1, u1) = g a" using lift_un'_ex[OF assms(1)] by auto
2276 have lu: "(l, u) = f l1 u1" using lift_un'_Some[unfolded Sa[symmetric] lift_un'.simps] by auto
2277 have "l = fst (f l1 u1)" and "u = snd (f l1 u1)" unfolding lu[symmetric] by auto
2278 thus ?thesis using Pa[OF Sa] by auto
2282 assumes lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f"
2283 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")
2284 shows "\<exists> l1 u1. (real l1 \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u1) \<and>
2285 l = fst (f l1 u1) \<and> u = snd (f l1 u1)"
2287 { fix l u assume "Some (l, u) = approx' prec a bs"
2288 with approx_approx'[of prec a bs, OF _ this] Pa
2289 have "real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" by auto } note Pa = this
2290 from lift_un'_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_un'_Some, OF Pa]
2291 show ?thesis by auto
2294 lemma lift_un'_bnds:
2295 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"
2296 and lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f"
2297 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"
2298 shows "real l \<le> f' (interpret_floatarith a xs) \<and> f' (interpret_floatarith a xs) \<le> real u"
2300 from lift_un'[OF lift_un'_Some Pa]
2301 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
2302 hence "(l, u) = f l1 u1" and "interpret_floatarith a xs \<in> {real l1 .. real u1}" by auto
2303 thus ?thesis using bnds by auto
2307 assumes lift_un_Some: "Some (l, u) = lift_un a f"
2308 shows "\<exists> l u. Some (l, u) = a"
2310 case None hence "None = lift_un a f" unfolding None lift_un.simps ..
2311 thus ?thesis using lift_un_Some by auto
2314 obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
2315 thus ?thesis unfolding `a = Some a'` a' by auto
2319 assumes lift_un_Some: "Some (l, u) = lift_un (g a) f"
2320 and Pa: "\<And>l u. Some (l, u) = g a \<Longrightarrow> P l u a"
2321 shows "\<exists> l1 u1. P l1 u1 a \<and> Some l = fst (f l1 u1) \<and> Some u = snd (f l1 u1)"
2323 obtain l1 u1 where Sa: "Some (l1, u1) = g a" using lift_un_ex[OF assms(1)] by auto
2324 have "fst (f l1 u1) \<noteq> None \<and> snd (f l1 u1) \<noteq> None"
2326 assume "\<not> (fst (f l1 u1) \<noteq> None \<and> snd (f l1 u1) \<noteq> None)"
2327 hence or: "fst (f l1 u1) = None \<or> snd (f l1 u1) = None" by auto
2328 hence "lift_un (g a) f = None"
2329 proof (cases "fst (f l1 u1) = None")
2331 then obtain b where b: "f l1 u1 = (None, b)" by (cases "f l1 u1", auto)
2332 thus ?thesis unfolding Sa[symmetric] lift_un.simps b by auto
2334 case False hence "snd (f l1 u1) = None" using or by auto
2335 with False obtain b where b: "f l1 u1 = (Some b, None)" by (cases "f l1 u1", auto)
2336 thus ?thesis unfolding Sa[symmetric] lift_un.simps b by auto
2338 thus False using lift_un_Some by auto
2340 then obtain a' b' where f: "f l1 u1 = (Some a', Some b')" by (cases "f l1 u1", auto)
2341 from lift_un_Some[unfolded Sa[symmetric] lift_un.simps f]
2342 have "Some l = fst (f l1 u1)" and "Some u = snd (f l1 u1)" unfolding f by auto
2343 thus ?thesis unfolding Sa[symmetric] lift_un.simps using Pa[OF Sa] by auto
2347 assumes lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f"
2348 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")
2349 shows "\<exists> l1 u1. (real l1 \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u1) \<and>
2350 Some l = fst (f l1 u1) \<and> Some u = snd (f l1 u1)"
2352 { fix l u assume "Some (l, u) = approx' prec a bs"
2353 with approx_approx'[of prec a bs, OF _ this] Pa
2354 have "real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" by auto } note Pa = this
2355 from lift_un_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_un_Some, OF Pa]
2356 show ?thesis by auto
2360 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"
2361 and lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f"
2362 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"
2363 shows "real l \<le> f' (interpret_floatarith a xs) \<and> f' (interpret_floatarith a xs) \<le> real u"
2365 from lift_un[OF lift_un_Some Pa]
2366 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
2367 hence "(Some l, Some u) = f l1 u1" and "interpret_floatarith a xs \<in> {real l1 .. real u1}" by auto
2368 thus ?thesis using bnds by auto
2372 assumes "bounded_by xs vs"
2373 and "Some (l, u) = approx prec arith vs" (is "_ = ?g arith")
2374 shows "real l \<le> interpret_floatarith arith xs \<and> interpret_floatarith arith xs \<le> real u" (is "?P l u arith")
2375 using `Some (l, u) = approx prec arith vs`
2376 proof (induct arith arbitrary: l u x)
2378 from lift_bin'[OF Add.prems[unfolded approx.simps]] Add.hyps
2379 obtain l1 u1 l2 u2 where "l = l1 + l2" and "u = u1 + u2"
2380 "real l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> real u1"
2381 "real l2 \<le> interpret_floatarith b xs" and "interpret_floatarith b xs \<le> real u2" unfolding fst_conv snd_conv by blast
2382 thus ?case unfolding interpret_floatarith.simps by auto
2385 from lift_un'[OF Minus.prems[unfolded approx.simps]] Minus.hyps
2386 obtain l1 u1 where "l = -u1" and "u = -l1"
2387 "real l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> real u1" unfolding fst_conv snd_conv by blast
2388 thus ?case unfolding interpret_floatarith.simps using real_of_float_minus by auto
2391 from lift_bin'[OF Mult.prems[unfolded approx.simps]] Mult.hyps
2393 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"
2394 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"
2395 and "real l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> real u1"
2396 and "real l2 \<le> interpret_floatarith b xs" and "interpret_floatarith b xs \<le> real u2" unfolding fst_conv snd_conv by blast
2397 thus ?case unfolding interpret_floatarith.simps l u real_of_float_add real_of_float_mult real_of_float_nprt real_of_float_pprt
2398 using mult_le_prts mult_ge_prts by auto
2401 from lift_un[OF Inverse.prems[unfolded approx.simps], unfolded if_distrib[of fst] if_distrib[of snd] fst_conv snd_conv] Inverse.hyps
2402 obtain l1 u1 where l': "Some l = (if 0 < l1 \<or> u1 < 0 then Some (float_divl prec 1 u1) else None)"
2403 and u': "Some u = (if 0 < l1 \<or> u1 < 0 then Some (float_divr prec 1 l1) else None)"
2404 and l1: "real l1 \<le> interpret_floatarith a xs" and u1: "interpret_floatarith a xs \<le> real u1" by blast
2405 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
2406 moreover have l1_le_u1: "real l1 \<le> real u1" using l1 u1 by auto
2407 ultimately have "real l1 \<noteq> 0" and "real u1 \<noteq> 0" unfolding less_float_def by auto
2409 have inv: "inverse (real u1) \<le> inverse (interpret_floatarith a xs)
2410 \<and> inverse (interpret_floatarith a xs) \<le> inverse (real l1)"
2411 proof (cases "0 < l1")
2412 case True hence "0 < real u1" and "0 < real l1" "0 < interpret_floatarith a xs"
2413 unfolding less_float_def using l1_le_u1 l1 by auto
2415 unfolding inverse_le_iff_le[OF `0 < real u1` `0 < interpret_floatarith a xs`]
2416 inverse_le_iff_le[OF `0 < interpret_floatarith a xs` `0 < real l1`]
2419 case False hence "u1 < 0" using either by blast
2420 hence "real u1 < 0" and "real l1 < 0" "interpret_floatarith a xs < 0"
2421 unfolding less_float_def using l1_le_u1 u1 by auto
2423 unfolding inverse_le_iff_le_neg[OF `real u1 < 0` `interpret_floatarith a xs < 0`]
2424 inverse_le_iff_le_neg[OF `interpret_floatarith a xs < 0` `real l1 < 0`]
2428 from l' have "l = float_divl prec 1 u1" by (cases "0 < l1 \<or> u1 < 0", auto)
2429 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
2430 also have "\<dots> \<le> inverse (interpret_floatarith a xs)" using inv by auto
2431 finally have "real l \<le> inverse (interpret_floatarith a xs)" .
2433 from u' have "u = float_divr prec 1 l1" by (cases "0 < l1 \<or> u1 < 0", auto)
2434 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
2435 hence "inverse (interpret_floatarith a xs) \<le> real u" by (rule order_trans[OF inv[THEN conjunct2]])
2436 ultimately show ?case unfolding interpret_floatarith.simps using l1 u1 by auto
2439 from lift_un'[OF Abs.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Abs.hyps
2440 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>"
2441 and l1: "real l1 \<le> interpret_floatarith x xs" and u1: "interpret_floatarith x xs \<le> real u1" by blast
2442 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)
2445 from lift_bin'[OF Min.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Min.hyps
2446 obtain l1 u1 l2 u2 where l': "l = min l1 l2" and u': "u = min u1 u2"
2447 and l1: "real l1 \<le> interpret_floatarith a xs" and u1: "interpret_floatarith a xs \<le> real u1"
2448 and l1: "real l2 \<le> interpret_floatarith b xs" and u1: "interpret_floatarith b xs \<le> real u2" by blast
2449 thus ?case unfolding l' u' by (auto simp add: real_of_float_min)
2452 from lift_bin'[OF Max.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Max.hyps
2453 obtain l1 u1 l2 u2 where l': "l = max l1 l2" and u': "u = max u1 u2"
2454 and l1: "real l1 \<le> interpret_floatarith a xs" and u1: "interpret_floatarith a xs \<le> real u1"
2455 and l1: "real l2 \<le> interpret_floatarith b xs" and u1: "interpret_floatarith b xs \<le> real u2" by blast
2456 thus ?case unfolding l' u' by (auto simp add: real_of_float_max)
2457 next case (Cos a) with lift_un'_bnds[OF bnds_cos] show ?case by auto
2458 next case (Arctan a) with lift_un'_bnds[OF bnds_arctan] show ?case by auto
2459 next case Pi with pi_boundaries show ?case by auto
2460 next case (Sqrt a) with lift_un'_bnds[OF bnds_sqrt] show ?case by auto
2461 next case (Exp a) with lift_un'_bnds[OF bnds_exp] show ?case by auto
2462 next case (Ln a) with lift_un_bnds[OF bnds_ln] show ?case by auto
2463 next case (Power a n) with lift_un'_bnds[OF bnds_power] show ?case by auto
2464 next case (Num f) thus ?case by auto
2467 from this[symmetric] `bounded_by xs vs`[THEN bounded_byE, of n]
2468 show ?case by (cases "n < length vs", auto)
2471 datatype form = Bound floatarith floatarith floatarith form
2472 | Assign floatarith floatarith form
2473 | Less floatarith floatarith
2474 | LessEqual floatarith floatarith
2475 | AtLeastAtMost floatarith floatarith floatarith
2477 fun interpret_form :: "form \<Rightarrow> real list \<Rightarrow> bool" where
2478 "interpret_form (Bound x a b f) vs = (interpret_floatarith x vs \<in> { interpret_floatarith a vs .. interpret_floatarith b vs } \<longrightarrow> interpret_form f vs)" |
2479 "interpret_form (Assign x a f) vs = (interpret_floatarith x vs = interpret_floatarith a vs \<longrightarrow> interpret_form f vs)" |
2480 "interpret_form (Less a b) vs = (interpret_floatarith a vs < interpret_floatarith b vs)" |
2481 "interpret_form (LessEqual a b) vs = (interpret_floatarith a vs \<le> interpret_floatarith b vs)" |
2482 "interpret_form (AtLeastAtMost x a b) vs = (interpret_floatarith x vs \<in> { interpret_floatarith a vs .. interpret_floatarith b vs })"
2484 fun approx_form' and approx_form :: "nat \<Rightarrow> form \<Rightarrow> (float * float) option list \<Rightarrow> nat list \<Rightarrow> bool" where
2485 "approx_form' prec f 0 n l u bs ss = approx_form prec f (bs[n := Some (l, u)]) ss" |
2486 "approx_form' prec f (Suc s) n l u bs ss =
2487 (let m = (l + u) * Float 1 -1
2488 in (if approx_form' prec f s n l m bs ss then approx_form' prec f s n m u bs ss else False))" |
2489 "approx_form prec (Bound (Var n) a b f) bs ss =
2490 (case (approx prec a bs, approx prec b bs)
2491 of (Some (l, _), Some (_, u)) \<Rightarrow> approx_form' prec f (ss ! n) n l u bs ss
2492 | _ \<Rightarrow> False)" |
2493 "approx_form prec (Assign (Var n) a f) bs ss =
2494 (case (approx prec a bs)
2495 of (Some (l, u)) \<Rightarrow> approx_form' prec f (ss ! n) n l u bs ss
2496 | _ \<Rightarrow> False)" |
2497 "approx_form prec (Less a b) bs ss =
2498 (case (approx prec a bs, approx prec b bs)
2499 of (Some (l, u), Some (l', u')) \<Rightarrow> u < l'
2500 | _ \<Rightarrow> False)" |
2501 "approx_form prec (LessEqual a b) bs ss =
2502 (case (approx prec a bs, approx prec b bs)
2503 of (Some (l, u), Some (l', u')) \<Rightarrow> u \<le> l'
2504 | _ \<Rightarrow> False)" |
2505 "approx_form prec (AtLeastAtMost x a b) bs ss =
2506 (case (approx prec x bs, approx prec a bs, approx prec b bs)
2507 of (Some (lx, ux), Some (l, u), Some (l', u')) \<Rightarrow> u \<le> lx \<and> ux \<le> l'
2508 | _ \<Rightarrow> False)" |
2509 "approx_form _ _ _ _ = False"
2511 lemma lazy_conj: "(if A then B else False) = (A \<and> B)" by simp
2513 lemma approx_form_approx_form':
2514 assumes "approx_form' prec f s n l u bs ss" and "x \<in> { real l .. real u }"
2515 obtains l' u' where "x \<in> { real l' .. real u' }"
2516 and "approx_form prec f (bs[n := Some (l', u')]) ss"
2517 using assms proof (induct s arbitrary: l u)
2519 from this(1)[of l u] this(2,3)
2524 let ?m = "(l + u) * Float 1 -1"
2525 have "real l \<le> real ?m" and "real ?m \<le> real u"
2526 unfolding le_float_def using Suc.prems by auto
2528 with `x \<in> { real l .. real u }`
2529 have "x \<in> { real l .. real ?m} \<or> x \<in> { real ?m .. real u }" by auto
2532 assume *: "x \<in> { real l .. real ?m }"
2533 with Suc.hyps[OF _ _ *] Suc.prems
2534 show thesis by (simp add: Let_def lazy_conj)
2536 assume *: "x \<in> { real ?m .. real u }"
2537 with Suc.hyps[OF _ _ *] Suc.prems
2538 show thesis by (simp add: Let_def lazy_conj)
2542 lemma approx_form_aux:
2543 assumes "approx_form prec f vs ss"
2544 and "bounded_by xs vs"
2545 shows "interpret_form f xs"
2546 using assms proof (induct f arbitrary: vs)
2547 case (Bound x a b f)
2549 where x_eq: "x = Var n" by (cases x) auto
2551 with Bound.prems obtain l u' l' u
2552 where l_eq: "Some (l, u') = approx prec a vs"
2553 and u_eq: "Some (l', u) = approx prec b vs"
2554 and approx_form': "approx_form' prec f (ss ! n) n l u vs ss"
2555 by (cases "approx prec a vs", simp,
2556 cases "approx prec b vs", auto) blast
2558 { assume "xs ! n \<in> { interpret_floatarith a xs .. interpret_floatarith b xs }"
2559 with approx[OF Bound.prems(2) l_eq] and approx[OF Bound.prems(2) u_eq]
2560 have "xs ! n \<in> { real l .. real u}" by auto
2562 from approx_form_approx_form'[OF approx_form' this]
2563 obtain lx ux where bnds: "xs ! n \<in> { real lx .. real ux }"
2564 and approx_form: "approx_form prec f (vs[n := Some (lx, ux)]) ss" .
2566 from `bounded_by xs vs` bnds
2567 have "bounded_by xs (vs[n := Some (lx, ux)])" by (rule bounded_by_update)
2568 with Bound.hyps[OF approx_form]
2569 have "interpret_form f xs" by blast }
2570 thus ?case using interpret_form.simps x_eq and interpret_floatarith.simps by simp
2574 where x_eq: "x = Var n" by (cases x) auto
2576 with Assign.prems obtain l u' l' u
2577 where bnd_eq: "Some (l, u) = approx prec a vs"
2578 and x_eq: "x = Var n"
2579 and approx_form': "approx_form' prec f (ss ! n) n l u vs ss"
2580 by (cases "approx prec a vs") auto
2582 { assume bnds: "xs ! n = interpret_floatarith a xs"
2583 with approx[OF Assign.prems(2) bnd_eq]
2584 have "xs ! n \<in> { real l .. real u}" by auto
2585 from approx_form_approx_form'[OF approx_form' this]
2586 obtain lx ux where bnds: "xs ! n \<in> { real lx .. real ux }"
2587 and approx_form: "approx_form prec f (vs[n := Some (lx, ux)]) ss" .
2589 from `bounded_by xs vs` bnds
2590 have "bounded_by xs (vs[n := Some (lx, ux)])" by (rule bounded_by_update)
2591 with Assign.hyps[OF approx_form]
2592 have "interpret_form f xs" by blast }
2593 thus ?case using interpret_form.simps x_eq and interpret_floatarith.simps by simp
2596 then obtain l u l' u'
2597 where l_eq: "Some (l, u) = approx prec a vs"
2598 and u_eq: "Some (l', u') = approx prec b vs"
2599 and inequality: "u < l'"
2600 by (cases "approx prec a vs", auto,
2601 cases "approx prec b vs", auto)
2602 from inequality[unfolded less_float_def] approx[OF Less.prems(2) l_eq] approx[OF Less.prems(2) u_eq]
2605 case (LessEqual a b)
2606 then obtain l u l' u'
2607 where l_eq: "Some (l, u) = approx prec a vs"
2608 and u_eq: "Some (l', u') = approx prec b vs"
2609 and inequality: "u \<le> l'"
2610 by (cases "approx prec a vs", auto,
2611 cases "approx prec b vs", auto)
2612 from inequality[unfolded le_float_def] approx[OF LessEqual.prems(2) l_eq] approx[OF LessEqual.prems(2) u_eq]
2615 case (AtLeastAtMost x a b)
2616 then obtain lx ux l u l' u'
2617 where x_eq: "Some (lx, ux) = approx prec x vs"
2618 and l_eq: "Some (l, u) = approx prec a vs"
2619 and u_eq: "Some (l', u') = approx prec b vs"
2620 and inequality: "u \<le> lx \<and> ux \<le> l'"
2621 by (cases "approx prec x vs", auto,
2622 cases "approx prec a vs", auto,
2623 cases "approx prec b vs", auto, blast)
2624 from inequality[unfolded le_float_def] approx[OF AtLeastAtMost.prems(2) l_eq] approx[OF AtLeastAtMost.prems(2) u_eq] approx[OF AtLeastAtMost.prems(2) x_eq]
2629 assumes "n = length xs"
2630 assumes "approx_form prec f (replicate n None) ss"
2631 shows "interpret_form f xs"
2632 using approx_form_aux[OF _ bounded_by_None] assms by auto
2634 subsection {* Implementing Taylor series expansion *}
2636 fun isDERIV :: "nat \<Rightarrow> floatarith \<Rightarrow> real list \<Rightarrow> bool" where
2637 "isDERIV x (Add a b) vs = (isDERIV x a vs \<and> isDERIV x b vs)" |
2638 "isDERIV x (Mult a b) vs = (isDERIV x a vs \<and> isDERIV x b vs)" |
2639 "isDERIV x (Minus a) vs = isDERIV x a vs" |
2640 "isDERIV x (Inverse a) vs = (isDERIV x a vs \<and> interpret_floatarith a vs \<noteq> 0)" |
2641 "isDERIV x (Cos a) vs = isDERIV x a vs" |
2642 "isDERIV x (Arctan a) vs = isDERIV x a vs" |
2643 "isDERIV x (Min a b) vs = False" |
2644 "isDERIV x (Max a b) vs = False" |
2645 "isDERIV x (Abs a) vs = False" |
2646 "isDERIV x Pi vs = True" |
2647 "isDERIV x (Sqrt a) vs = (isDERIV x a vs \<and> interpret_floatarith a vs > 0)" |
2648 "isDERIV x (Exp a) vs = isDERIV x a vs" |
2649 "isDERIV x (Ln a) vs = (isDERIV x a vs \<and> interpret_floatarith a vs > 0)" |
2650 "isDERIV x (Power a 0) vs = True" |
2651 "isDERIV x (Power a (Suc n)) vs = isDERIV x a vs" |
2652 "isDERIV x (Num f) vs = True" |
2653 "isDERIV x (Var n) vs = True"
2655 fun DERIV_floatarith :: "nat \<Rightarrow> floatarith \<Rightarrow> floatarith" where
2656 "DERIV_floatarith x (Add a b) = Add (DERIV_floatarith x a) (DERIV_floatarith x b)" |
2657 "DERIV_floatarith x (Mult a b) = Add (Mult a (DERIV_floatarith x b)) (Mult (DERIV_floatarith x a) b)" |
2658 "DERIV_floatarith x (Minus a) = Minus (DERIV_floatarith x a)" |
2659 "DERIV_floatarith x (Inverse a) = Minus (Mult (DERIV_floatarith x a) (Inverse (Power a 2)))" |
2660 "DERIV_floatarith x (Cos a) = Minus (Mult (Cos (Add (Mult Pi (Num (Float 1 -1))) (Minus a))) (DERIV_floatarith x a))" |
2661 "DERIV_floatarith x (Arctan a) = Mult (Inverse (Add (Num 1) (Power a 2))) (DERIV_floatarith x a)" |
2662 "DERIV_floatarith x (Min a b) = Num 0" |
2663 "DERIV_floatarith x (Max a b) = Num 0" |
2664 "DERIV_floatarith x (Abs a) = Num 0" |
2665 "DERIV_floatarith x Pi = Num 0" |
2666 "DERIV_floatarith x (Sqrt a) = (Mult (Inverse (Mult (Sqrt a) (Num 2))) (DERIV_floatarith x a))" |
2667 "DERIV_floatarith x (Exp a) = Mult (Exp a) (DERIV_floatarith x a)" |
2668 "DERIV_floatarith x (Ln a) = Mult (Inverse a) (DERIV_floatarith x a)" |
2669 "DERIV_floatarith x (Power a 0) = Num 0" |
2670 "DERIV_floatarith x (Power a (Suc n)) = Mult (Num (Float (int (Suc n)) 0)) (Mult (Power a n) (DERIV_floatarith x a))" |
2671 "DERIV_floatarith x (Num f) = Num 0" |
2672 "DERIV_floatarith x (Var n) = (if x = n then Num 1 else Num 0)"
2674 lemma DERIV_floatarith:
2675 assumes "n < length vs"
2676 assumes isDERIV: "isDERIV n f (vs[n := x])"
2677 shows "DERIV (\<lambda> x'. interpret_floatarith f (vs[n := x'])) x :>
2678 interpret_floatarith (DERIV_floatarith n f) (vs[n := x])"
2679 (is "DERIV (?i f) x :> _")
2680 using isDERIV proof (induct f arbitrary: x)
2681 case (Inverse a) thus ?case
2682 by (auto intro!: DERIV_intros
2683 simp add: algebra_simps power2_eq_square)
2684 next case (Cos a) thus ?case
2685 by (auto intro!: DERIV_intros
2686 simp del: interpret_floatarith.simps(5)
2687 simp add: interpret_floatarith_sin interpret_floatarith.simps(5)[of a])
2688 next case (Power a n) thus ?case
2689 by (cases n, auto intro!: DERIV_intros
2690 simp del: power_Suc simp add: real_eq_of_nat)
2691 next case (Ln a) thus ?case
2692 by (auto intro!: DERIV_intros simp add: divide_inverse)
2693 next case (Var i) thus ?case using `n < length vs` by auto
2694 qed (auto intro!: DERIV_intros)
2696 declare approx.simps[simp del]
2698 fun isDERIV_approx :: "nat \<Rightarrow> nat \<Rightarrow> floatarith \<Rightarrow> (float * float) option list \<Rightarrow> bool" where
2699 "isDERIV_approx prec x (Add a b) vs = (isDERIV_approx prec x a vs \<and> isDERIV_approx prec x b vs)" |
2700 "isDERIV_approx prec x (Mult a b) vs = (isDERIV_approx prec x a vs \<and> isDERIV_approx prec x b vs)" |
2701 "isDERIV_approx prec x (Minus a) vs = isDERIV_approx prec x a vs" |
2702 "isDERIV_approx prec x (Inverse a) vs =
2703 (isDERIV_approx prec x a vs \<and> (case approx prec a vs of Some (l, u) \<Rightarrow> 0 < l \<or> u < 0 | None \<Rightarrow> False))" |
2704 "isDERIV_approx prec x (Cos a) vs = isDERIV_approx prec x a vs" |
2705 "isDERIV_approx prec x (Arctan a) vs = isDERIV_approx prec x a vs" |
2706 "isDERIV_approx prec x (Min a b) vs = False" |
2707 "isDERIV_approx prec x (Max a b) vs = False" |
2708 "isDERIV_approx prec x (Abs a) vs = False" |
2709 "isDERIV_approx prec x Pi vs = True" |
2710 "isDERIV_approx prec x (Sqrt a) vs =
2711 (isDERIV_approx prec x a vs \<and> (case approx prec a vs of Some (l, u) \<Rightarrow> 0 < l | None \<Rightarrow> False))" |
2712 "isDERIV_approx prec x (Exp a) vs = isDERIV_approx prec x a vs" |
2713 "isDERIV_approx prec x (Ln a) vs =
2714 (isDERIV_approx prec x a vs \<and> (case approx prec a vs of Some (l, u) \<Rightarrow> 0 < l | None \<Rightarrow> False))" |
2715 "isDERIV_approx prec x (Power a 0) vs = True" |
2716 "isDERIV_approx prec x (Power a (Suc n)) vs = isDERIV_approx prec x a vs" |
2717 "isDERIV_approx prec x (Num f) vs = True" |
2718 "isDERIV_approx prec x (Var n) vs = True"
2720 lemma isDERIV_approx:
2721 assumes "bounded_by xs vs"
2722 and isDERIV_approx: "isDERIV_approx prec x f vs"
2723 shows "isDERIV x f xs"
2724 using isDERIV_approx proof (induct f)
2726 then obtain l u where approx_Some: "Some (l, u) = approx prec a vs"
2727 and *: "0 < l \<or> u < 0"
2728 by (cases "approx prec a vs", auto)
2729 with approx[OF `bounded_by xs vs` approx_Some]
2730 have "interpret_floatarith a xs \<noteq> 0" unfolding less_float_def by auto
2731 thus ?case using Inverse by auto
2734 then obtain l u where approx_Some: "Some (l, u) = approx prec a vs"
2736 by (cases "approx prec a vs", auto)
2737 with approx[OF `bounded_by xs vs` approx_Some]
2738 have "0 < interpret_floatarith a xs" unfolding less_float_def by auto
2739 thus ?case using Ln by auto
2742 then obtain l u where approx_Some: "Some (l, u) = approx prec a vs"
2744 by (cases "approx prec a vs", auto)
2745 with approx[OF `bounded_by xs vs` approx_Some]
2746 have "0 < interpret_floatarith a xs" unfolding less_float_def by auto
2747 thus ?case using Sqrt by auto
2749 case (Power a n) thus ?case by (cases n, auto)
2752 lemma bounded_by_update_var:
2753 assumes "bounded_by xs vs" and "vs ! i = Some (l, u)"
2754 and bnd: "x \<in> { real l .. real u }"
2755 shows "bounded_by (xs[i := x]) vs"
2756 proof (cases "i < length xs")
2757 case False thus ?thesis using `bounded_by xs vs` by auto
2759 let ?xs = "xs[i := x]"
2760 case True hence "i < length ?xs" by auto
2762 assume "j < length vs"
2763 have "case vs ! j of None \<Rightarrow> True | Some (l, u) \<Rightarrow> ?xs ! j \<in> { real l .. real u }"
2764 proof (cases "vs ! j")
2767 proof (cases "i = j")
2769 thus ?thesis using `vs ! i = Some (l, u)` Some and bnd `i < length ?xs`
2773 thus ?thesis using `bounded_by xs vs`[THEN bounded_byE, OF `j < length vs`] Some
2777 thus ?thesis unfolding bounded_by_def by auto
2780 lemma isDERIV_approx':
2781 assumes "bounded_by xs vs"
2782 and vs_x: "vs ! x = Some (l, u)" and X_in: "X \<in> { real l .. real u }"
2783 and approx: "isDERIV_approx prec x f vs"
2784 shows "isDERIV x f (xs[x := X])"
2786 note bounded_by_update_var[OF `bounded_by xs vs` vs_x X_in] approx
2787 thus ?thesis by (rule isDERIV_approx)
2791 assumes "n < length xs" and bnd: "bounded_by xs vs"
2792 and isD: "isDERIV_approx prec n f vs"
2793 and app: "Some (l, u) = approx prec (DERIV_floatarith n f) vs" (is "_ = approx _ ?D _")
2794 shows "\<exists>x. real l \<le> x \<and> x \<le> real u \<and>
2795 DERIV (\<lambda> x. interpret_floatarith f (xs[n := x])) (xs!n) :> x"
2796 (is "\<exists> x. _ \<and> _ \<and> DERIV (?i f) _ :> _")
2797 proof (rule exI[of _ "?i ?D (xs!n)"], rule conjI[OF _ conjI])
2798 let "?i f x" = "interpret_floatarith f (xs[n := x])"
2799 from approx[OF bnd app]
2800 show "real l \<le> ?i ?D (xs!n)" and "?i ?D (xs!n) \<le> real u"
2801 using `n < length xs` by auto
2802 from DERIV_floatarith[OF `n < length xs`, of f "xs!n"] isDERIV_approx[OF bnd isD]
2803 show "DERIV (?i f) (xs!n) :> (?i ?D (xs!n))" by simp
2806 fun lift_bin :: "(float * float) option \<Rightarrow> (float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> (float * float) option) \<Rightarrow> (float * float) option" where
2807 "lift_bin (Some (l1, u1)) (Some (l2, u2)) f = f l1 u1 l2 u2" |
2808 "lift_bin a b f = None"
2811 assumes lift_bin_Some: "Some (l, u) = lift_bin a b f"
2813 where "a = Some (l1, u1)"
2814 and "b = Some (l2, u2)"
2815 and "f l1 u1 l2 u2 = Some (l, u)"
2816 using assms by (cases a, simp, cases b, simp, auto)
2818 fun approx_tse where
2819 "approx_tse prec n 0 c k f bs = approx prec f bs" |
2820 "approx_tse prec n (Suc s) c k f bs =
2821 (if isDERIV_approx prec n f bs then
2822 lift_bin (approx prec f (bs[n := Some (c,c)]))
2823 (approx_tse prec n s c (Suc k) (DERIV_floatarith n f) bs)
2824 (\<lambda> l1 u1 l2 u2. approx prec
2826 (Mult (Inverse (Num (Float (int k) 0)))
2827 (Mult (Add (Var (Suc (Suc 0))) (Minus (Num c)))
2828 (Var (Suc 0))))) [Some (l1, u1), Some (l2, u2), bs!n])
2829 else approx prec f bs)"
2831 lemma bounded_by_Cons:
2832 assumes bnd: "bounded_by xs vs"
2833 and x: "x \<in> { real l .. real u }"
2834 shows "bounded_by (x#xs) ((Some (l, u))#vs)"
2836 { fix i assume *: "i < length ((Some (l, u))#vs)"
2837 have "case ((Some (l,u))#vs) ! i of Some (l, u) \<Rightarrow> (x#xs)!i \<in> { real l .. real u } | None \<Rightarrow> True"
2839 case 0 with x show ?thesis by auto
2841 case (Suc i) with * have "i < length vs" by auto
2842 from bnd[THEN bounded_byE, OF this]
2843 show ?thesis unfolding Suc nth_Cons_Suc .
2845 thus ?thesis by (auto simp add: bounded_by_def)
2848 lemma approx_tse_generic:
2849 assumes "bounded_by xs vs"
2850 and bnd_c: "bounded_by (xs[x := real c]) vs" and "x < length vs" and "x < length xs"
2851 and bnd_x: "vs ! x = Some (lx, ux)"
2852 and ate: "Some (l, u) = approx_tse prec x s c k f vs"
2853 shows "\<exists> n. (\<forall> m < n. \<forall> z \<in> {real lx .. real ux}.
2854 DERIV (\<lambda> y. interpret_floatarith ((DERIV_floatarith x ^^ m) f) (xs[x := y])) z :>
2855 (interpret_floatarith ((DERIV_floatarith x ^^ (Suc m)) f) (xs[x := z])))
2856 \<and> (\<forall> t \<in> {real lx .. real ux}. (\<Sum> i = 0..<n. inverse (real (\<Prod> j \<in> {k..<k+i}. j)) *
2857 interpret_floatarith ((DERIV_floatarith x ^^ i) f) (xs[x := real c]) *
2858 (xs!x - real c)^i) +
2859 inverse (real (\<Prod> j \<in> {k..<k+n}. j)) *
2860 interpret_floatarith ((DERIV_floatarith x ^^ n) f) (xs[x := t]) *
2861 (xs!x - real c)^n \<in> {real l .. real u})" (is "\<exists> n. ?taylor f k l u n")
2862 using ate proof (induct s arbitrary: k f l u)
2864 { fix t assume "t \<in> {real lx .. real ux}"
2865 note bounded_by_update_var[OF `bounded_by xs vs` bnd_x this]
2866 from approx[OF this 0[unfolded approx_tse.simps]]
2867 have "(interpret_floatarith f (xs[x := t])) \<in> {real l .. real u}"
2868 by (auto simp add: algebra_simps)
2869 } thus ?case by (auto intro!: exI[of _ 0])
2873 proof (cases "isDERIV_approx prec x f vs")
2875 note ap = Suc.prems[unfolded approx_tse.simps if_not_P[OF False]]
2877 { fix t assume "t \<in> {real lx .. real ux}"
2878 note bounded_by_update_var[OF `bounded_by xs vs` bnd_x this]
2879 from approx[OF this ap]
2880 have "(interpret_floatarith f (xs[x := t])) \<in> {real l .. real u}"
2881 by (auto simp add: algebra_simps)
2882 } thus ?thesis by (auto intro!: exI[of _ 0])
2887 where a: "Some (l1, u1) = approx prec f (vs[x := Some (c,c)])"
2888 and ate: "Some (l2, u2) = approx_tse prec x s c (Suc k) (DERIV_floatarith x f) vs"
2889 and final: "Some (l, u) = approx prec
2891 (Mult (Inverse (Num (Float (int k) 0)))
2892 (Mult (Add (Var (Suc (Suc 0))) (Minus (Num c)))
2893 (Var (Suc 0))))) [Some (l1, u1), Some (l2, u2), vs!x]"
2894 by (auto elim!: lift_bin) blast
2896 from bnd_c `x < length xs`
2897 have bnd: "bounded_by (xs[x:=real c]) (vs[x:= Some (c,c)])"
2898 by (auto intro!: bounded_by_update)
2900 from approx[OF this a]
2901 have f_c: "interpret_floatarith ((DERIV_floatarith x ^^ 0) f) (xs[x := real c]) \<in> { real l1 .. real u1 }"
2902 (is "?f 0 (real c) \<in> _")
2905 { fix f :: "'a \<Rightarrow> 'a" fix n :: nat fix x :: 'a
2906 have "(f ^^ Suc n) x = (f ^^ n) (f x)"
2907 by (induct n, auto) }
2908 note funpow_Suc = this[symmetric]
2909 from Suc.hyps[OF ate, unfolded this]
2911 where DERIV_hyp: "\<And> m z. \<lbrakk> m < n ; z \<in> { real lx .. real ux } \<rbrakk> \<Longrightarrow> DERIV (?f (Suc m)) z :> ?f (Suc (Suc m)) z"
2912 and hyp: "\<forall> t \<in> {real lx .. real ux}. (\<Sum> i = 0..<n. inverse (real (\<Prod> j \<in> {Suc k..<Suc k + i}. j)) * ?f (Suc i) (real c) * (xs!x - real c)^i) +
2913 inverse (real (\<Prod> j \<in> {Suc k..<Suc k + n}. j)) * ?f (Suc n) t * (xs!x - real c)^n \<in> {real l2 .. real u2}"
2914 (is "\<forall> t \<in> _. ?X (Suc k) f n t \<in> _")
2918 assume "m < Suc n" and bnd_z: "z \<in> { real lx .. real ux }"
2919 have "DERIV (?f m) z :> ?f (Suc m) z"
2922 with DERIV_floatarith[OF `x < length xs` isDERIV_approx'[OF `bounded_by xs vs` bnd_x bnd_z True]]
2923 show ?thesis by simp
2926 hence "m' < n" using `m < Suc n` by auto
2927 from DERIV_hyp[OF this bnd_z]
2928 show ?thesis using Suc by simp
2929 qed } note DERIV = this
2931 have "\<And> k i. k < i \<Longrightarrow> {k ..< i} = insert k {Suc k ..< i}" by auto
2932 hence setprod_head_Suc: "\<And> k i. \<Prod> {k ..< k + Suc i} = k * \<Prod> {Suc k ..< Suc k + i}" by auto
2933 have setsum_move0: "\<And> k F. setsum F {0..<Suc k} = F 0 + setsum (\<lambda> k. F (Suc k)) {0..<k}"
2934 unfolding setsum_shift_bounds_Suc_ivl[symmetric]
2935 unfolding setsum_head_upt_Suc[OF zero_less_Suc] ..
2936 def C \<equiv> "xs!x - real c"
2938 { fix t assume t: "t \<in> {real lx .. real ux}"
2939 hence "bounded_by [xs!x] [vs!x]"
2940 using `bounded_by xs vs`[THEN bounded_byE, OF `x < length vs`]
2941 by (cases "vs!x", auto simp add: bounded_by_def)
2943 with hyp[THEN bspec, OF t] f_c
2944 have "bounded_by [?f 0 (real c), ?X (Suc k) f n t, xs!x] [Some (l1, u1), Some (l2, u2), vs!x]"
2945 by (auto intro!: bounded_by_Cons)
2946 from approx[OF this final, unfolded atLeastAtMost_iff[symmetric]]
2947 have "?X (Suc k) f n t * (xs!x - real c) * inverse (real k) + ?f 0 (real c) \<in> {real l .. real u}"
2948 by (auto simp add: algebra_simps)
2949 also have "?X (Suc k) f n t * (xs!x - real c) * inverse (real k) + ?f 0 (real c) =
2950 (\<Sum> i = 0..<Suc n. inverse (real (\<Prod> j \<in> {k..<k+i}. j)) * ?f i (real c) * (xs!x - real c)^i) +
2951 inverse (real (\<Prod> j \<in> {k..<k+Suc n}. j)) * ?f (Suc n) t * (xs!x - real c)^Suc n" (is "_ = ?T")
2952 unfolding funpow_Suc C_def[symmetric] setsum_move0 setprod_head_Suc
2953 by (auto simp add: algebra_simps)
2954 (simp only: mult_left_commute [of _ "inverse (real k)"] setsum_right_distrib [symmetric])
2955 finally have "?T \<in> {real l .. real u}" . }
2956 thus ?thesis using DERIV by blast
2960 lemma setprod_fact: "\<Prod> {1..<1 + k} = fact (k :: nat)"
2963 have "{ 1 ..< Suc (Suc k) } = insert (Suc k) { 1 ..< Suc k }" by auto
2964 hence "\<Prod> { 1 ..< Suc (Suc k) } = (Suc k) * \<Prod> { 1 ..< Suc k }" by auto
2965 thus ?case using Suc by auto
2969 assumes "bounded_by xs vs"
2970 and bnd_x: "vs ! x = Some (lx, ux)" and bnd_c: "real c \<in> {real lx .. real ux}"
2971 and "x < length vs" and "x < length xs"
2972 and ate: "Some (l, u) = approx_tse prec x s c 1 f vs"
2973 shows "interpret_floatarith f xs \<in> { real l .. real u }"
2975 def F \<equiv> "\<lambda> n z. interpret_floatarith ((DERIV_floatarith x ^^ n) f) (xs[x := z])"
2976 hence F0: "F 0 = (\<lambda> z. interpret_floatarith f (xs[x := z]))" by auto
2978 hence "bounded_by (xs[x := real c]) vs" and "x < length vs" "x < length xs"
2979 using `bounded_by xs vs` bnd_x bnd_c `x < length vs` `x < length xs`
2980 by (auto intro!: bounded_by_update_var)
2982 from approx_tse_generic[OF `bounded_by xs vs` this bnd_x ate]
2984 where DERIV: "\<forall> m z. m < n \<and> real lx \<le> z \<and> z \<le> real ux \<longrightarrow> DERIV (F m) z :> F (Suc m) z"
2985 and hyp: "\<And> t. t \<in> {real lx .. real ux} \<Longrightarrow>
2986 (\<Sum> j = 0..<n. inverse (real (fact j)) * F j (real c) * (xs!x - real c)^j) +
2987 inverse (real (fact n)) * F n t * (xs!x - real c)^n
2988 \<in> {real l .. real u}" (is "\<And> t. _ \<Longrightarrow> ?taylor t \<in> _")
2989 unfolding F_def atLeastAtMost_iff[symmetric] setprod_fact by blast
2991 have bnd_xs: "xs ! x \<in> { real lx .. real ux }"
2992 using `bounded_by xs vs`[THEN bounded_byE, OF `x < length vs`] bnd_x by auto
2996 case 0 thus ?thesis using hyp[OF bnd_xs] unfolding F_def by auto
3000 proof (cases "xs ! x = real c")
3002 from True[symmetric] hyp[OF bnd_xs] Suc show ?thesis
3003 unfolding F_def Suc setsum_head_upt_Suc[OF zero_less_Suc] setsum_shift_bounds_Suc_ivl by auto
3007 have "real lx \<le> real c" "real c \<le> real ux" "real lx \<le> xs!x" "xs!x \<le> real ux"
3008 using Suc bnd_c `bounded_by xs vs`[THEN bounded_byE, OF `x < length vs`] bnd_x by auto
3009 from Taylor.taylor[OF zero_less_Suc, of F, OF F0 DERIV[unfolded Suc] this False]
3010 obtain t where t_bnd: "if xs ! x < real c then xs ! x < t \<and> t < real c else real c < t \<and> t < xs ! x"
3011 and fl_eq: "interpret_floatarith f (xs[x := xs ! x]) =
3012 (\<Sum>m = 0..<Suc n'. F m (real c) / real (fact m) * (xs ! x - real c) ^ m) +
3013 F (Suc n') t / real (fact (Suc n')) * (xs ! x - real c) ^ Suc n'"
3016 from t_bnd bnd_xs bnd_c have *: "t \<in> {real lx .. real ux}"
3017 by (cases "xs ! x < real c", auto)
3019 have "interpret_floatarith f (xs[x := xs ! x]) = ?taylor t"
3020 unfolding fl_eq Suc by (auto simp add: algebra_simps divide_inverse)
3021 also have "\<dots> \<in> {real l .. real u}" using * by (rule hyp)
3022 finally show ?thesis by simp
3027 fun approx_tse_form' where
3028 "approx_tse_form' prec t f 0 l u cmp =
3029 (case approx_tse prec 0 t ((l + u) * Float 1 -1) 1 f [Some (l, u)]
3030 of Some (l, u) \<Rightarrow> cmp l u | None \<Rightarrow> False)" |
3031 "approx_tse_form' prec t f (Suc s) l u cmp =
3032 (let m = (l + u) * Float 1 -1
3033 in (if approx_tse_form' prec t f s l m cmp then
3034 approx_tse_form' prec t f s m u cmp else False))"
3036 lemma approx_tse_form':
3037 assumes "approx_tse_form' prec t f s l u cmp" and "x \<in> {real l .. real u}"
3038 shows "\<exists> l' u' ly uy. x \<in> { real l' .. real u' } \<and> real l \<le> real l' \<and> real u' \<le> real u \<and> cmp ly uy \<and>
3039 approx_tse prec 0 t ((l' + u') * Float 1 -1) 1 f [Some (l', u')] = Some (ly, uy)"
3040 using assms proof (induct s arbitrary: l u)
3043 where *: "approx_tse prec 0 t ((l + u) * Float 1 -1) 1 f [Some (l, u)] = Some (ly, uy)"
3044 and **: "cmp ly uy" by (auto elim!: option_caseE)
3045 with 0 show ?case by (auto intro!: exI)
3048 let ?m = "(l + u) * Float 1 -1"
3050 have l: "approx_tse_form' prec t f s l ?m cmp"
3051 and u: "approx_tse_form' prec t f s ?m u cmp"
3052 by (auto simp add: Let_def lazy_conj)
3054 have m_l: "real l \<le> real ?m" and m_u: "real ?m \<le> real u"
3055 unfolding le_float_def using Suc.prems by auto
3057 with `x \<in> { real l .. real u }`
3058 have "x \<in> { real l .. real ?m} \<or> x \<in> { real ?m .. real u }" by auto
3061 assume "x \<in> { real l .. real ?m}"
3062 from Suc.hyps[OF l this]
3064 where "x \<in> { real l' .. real u' } \<and> real l \<le> real l' \<and> real u' \<le> real ?m \<and> cmp ly uy \<and>
3065 approx_tse prec 0 t ((l' + u') * Float 1 -1) 1 f [Some (l', u')] = Some (ly, uy)" by blast
3066 with m_u show ?thesis by (auto intro!: exI)
3068 assume "x \<in> { real ?m .. real u }"
3069 from Suc.hyps[OF u this]
3071 where "x \<in> { real l' .. real u' } \<and> real ?m \<le> real l' \<and> real u' \<le> real u \<and> cmp ly uy \<and>
3072 approx_tse prec 0 t ((l' + u') * Float 1 -1) 1 f [Some (l', u')] = Some (ly, uy)" by blast
3073 with m_u show ?thesis by (auto intro!: exI)
3077 lemma approx_tse_form'_less:
3078 assumes tse: "approx_tse_form' prec t (Add a (Minus b)) s l u (\<lambda> l u. 0 < l)"
3079 and x: "x \<in> {real l .. real u}"
3080 shows "interpret_floatarith b [x] < interpret_floatarith a [x]"
3082 from approx_tse_form'[OF tse x]
3084 where x': "x \<in> { real l' .. real u' }" and "real l \<le> real l'"
3085 and "real u' \<le> real u" and "0 < ly"
3086 and tse: "approx_tse prec 0 t ((l' + u') * Float 1 -1) 1 (Add a (Minus b)) [Some (l', u')] = Some (ly, uy)"
3089 hence "bounded_by [x] [Some (l', u')]" by (auto simp add: bounded_by_def)
3091 from approx_tse[OF this _ _ _ _ tse[symmetric], of l' u'] x'
3092 have "real ly \<le> interpret_floatarith a [x] - interpret_floatarith b [x]"
3093 by (auto simp add: diff_minus)
3094 from order_less_le_trans[OF `0 < ly`[unfolded less_float_def] this]
3095 show ?thesis by auto
3098 lemma approx_tse_form'_le:
3099 assumes tse: "approx_tse_form' prec t (Add a (Minus b)) s l u (\<lambda> l u. 0 \<le> l)"
3100 and x: "x \<in> {real l .. real u}"
3101 shows "interpret_floatarith b [x] \<le> interpret_floatarith a [x]"
3103 from approx_tse_form'[OF tse x]
3105 where x': "x \<in> { real l' .. real u' }" and "real l \<le> real l'"
3106 and "real u' \<le> real u" and "0 \<le> ly"
3107 and tse: "approx_tse prec 0 t ((l' + u') * Float 1 -1) 1 (Add a (Minus b)) [Some (l', u')] = Some (ly, uy)"
3110 hence "bounded_by [x] [Some (l', u')]" by (auto simp add: bounded_by_def)
3112 from approx_tse[OF this _ _ _ _ tse[symmetric], of l' u'] x'
3113 have "real ly \<le> interpret_floatarith a [x] - interpret_floatarith b [x]"
3114 by (auto simp add: diff_minus)
3115 from order_trans[OF `0 \<le> ly`[unfolded le_float_def] this]
3116 show ?thesis by auto
3120 "approx_tse_form prec t s f =
3122 of (Bound x a b f) \<Rightarrow> x = Var 0 \<and>
3123 (case (approx prec a [None], approx prec b [None])
3124 of (Some (l, u), Some (l', u')) \<Rightarrow>
3126 of Less lf rt \<Rightarrow> approx_tse_form' prec t (Add rt (Minus lf)) s l u' (\<lambda> l u. 0 < l)
3127 | LessEqual lf rt \<Rightarrow> approx_tse_form' prec t (Add rt (Minus lf)) s l u' (\<lambda> l u. 0 \<le> l)
3128 | AtLeastAtMost x lf rt \<Rightarrow>
3129 (if approx_tse_form' prec t (Add x (Minus lf)) s l u' (\<lambda> l u. 0 \<le> l) then
3130 approx_tse_form' prec t (Add rt (Minus x)) s l u' (\<lambda> l u. 0 \<le> l) else False)
3131 | _ \<Rightarrow> False)
3132 | _ \<Rightarrow> False)
3133 | _ \<Rightarrow> False)"
3135 lemma approx_tse_form:
3136 assumes "approx_tse_form prec t s f"
3137 shows "interpret_form f [x]"
3139 case (Bound i a b f') note f_def = this
3140 with assms obtain l u l' u'
3141 where a: "approx prec a [None] = Some (l, u)"
3142 and b: "approx prec b [None] = Some (l', u')"
3143 unfolding approx_tse_form_def by (auto elim!: option_caseE)
3145 from Bound assms have "i = Var 0" unfolding approx_tse_form_def by auto
3146 hence i: "interpret_floatarith i [x] = x" by auto
3148 { let "?f z" = "interpret_floatarith z [x]"
3149 assume "?f i \<in> { ?f a .. ?f b }"
3150 with approx[OF _ a[symmetric], of "[x]"] approx[OF _ b[symmetric], of "[x]"]
3151 have bnd: "x \<in> { real l .. real u'}" unfolding bounded_by_def i by auto
3153 have "interpret_form f' [x]"
3156 with Bound a b assms
3157 have "approx_tse_form' prec t (Add rt (Minus lf)) s l u' (\<lambda> l u. 0 < l)"
3158 unfolding approx_tse_form_def by auto
3159 from approx_tse_form'_less[OF this bnd]
3160 show ?thesis using Less by auto
3162 case (LessEqual lf rt)
3163 with Bound a b assms
3164 have "approx_tse_form' prec t (Add rt (Minus lf)) s l u' (\<lambda> l u. 0 \<le> l)"
3165 unfolding approx_tse_form_def by auto
3166 from approx_tse_form'_le[OF this bnd]
3167 show ?thesis using LessEqual by auto
3169 case (AtLeastAtMost x lf rt)
3170 with Bound a b assms
3171 have "approx_tse_form' prec t (Add rt (Minus x)) s l u' (\<lambda> l u. 0 \<le> l)"
3172 and "approx_tse_form' prec t (Add x (Minus lf)) s l u' (\<lambda> l u. 0 \<le> l)"
3173 unfolding approx_tse_form_def lazy_conj by auto
3174 from approx_tse_form'_le[OF this(1) bnd] approx_tse_form'_le[OF this(2) bnd]
3175 show ?thesis using AtLeastAtMost by auto
3177 case (Bound x a b f') with assms
3178 show ?thesis by (auto elim!: option_caseE simp add: f_def approx_tse_form_def)
3180 case (Assign x a f') with assms
3181 show ?thesis by (auto elim!: option_caseE simp add: f_def approx_tse_form_def)
3182 qed } thus ?thesis unfolding f_def by auto
3183 next case Assign with assms show ?thesis by (auto simp add: approx_tse_form_def)
3184 next case LessEqual with assms show ?thesis by (auto simp add: approx_tse_form_def)
3185 next case Less with assms show ?thesis by (auto simp add: approx_tse_form_def)
3186 next case AtLeastAtMost with assms show ?thesis by (auto simp add: approx_tse_form_def)
3189 text {* @{term approx_form_eval} is only used for the {\tt value}-command. *}
3191 fun approx_form_eval :: "nat \<Rightarrow> form \<Rightarrow> (float * float) option list \<Rightarrow> (float * float) option list" where
3192 "approx_form_eval prec (Bound (Var n) a b f) bs =
3193 (case (approx prec a bs, approx prec b bs)
3194 of (Some (l, _), Some (_, u)) \<Rightarrow> approx_form_eval prec f (bs[n := Some (l, u)])
3195 | _ \<Rightarrow> bs)" |
3196 "approx_form_eval prec (Assign (Var n) a f) bs =
3197 (case (approx prec a bs)
3198 of (Some (l, u)) \<Rightarrow> approx_form_eval prec f (bs[n := Some (l, u)])
3199 | _ \<Rightarrow> bs)" |
3200 "approx_form_eval prec (Less a b) bs = bs @ [approx prec a bs, approx prec b bs]" |
3201 "approx_form_eval prec (LessEqual a b) bs = bs @ [approx prec a bs, approx prec b bs]" |
3202 "approx_form_eval prec (AtLeastAtMost x a b) bs =
3203 bs @ [approx prec x bs, approx prec a bs, approx prec b bs]" |
3204 "approx_form_eval _ _ bs = bs"
3206 subsection {* Implement proof method \texttt{approximation} *}
3208 lemmas interpret_form_equations = interpret_form.simps interpret_floatarith.simps interpret_floatarith_num
3209 interpret_floatarith_divide interpret_floatarith_diff interpret_floatarith_tan interpret_floatarith_powr interpret_floatarith_log
3210 interpret_floatarith_sin
3212 code_reflect Float_Arith
3213 datatypes float = Float
3214 and floatarith = Add | Minus | Mult | Inverse | Cos | Arctan
3215 | Abs | Max | Min | Pi | Sqrt | Exp | Ln | Power | Var | Num
3216 and form = Bound | Assign | Less | LessEqual | AtLeastAtMost
3217 functions approx_form approx_tse_form approx' approx_form_eval
3220 fun reorder_bounds_tac prems i =
3222 fun variable_of_bound (Const ("Trueprop", _) $
3223 (Const (@{const_name "op :"}, _) $
3224 Free (name, _) $ _)) = name
3225 | variable_of_bound (Const ("Trueprop", _) $
3226 (Const ("op =", _) $
3227 Free (name, _) $ _)) = name
3228 | variable_of_bound t = raise TERM ("variable_of_bound", [t])
3231 = map (` (variable_of_bound o prop_of)) prems
3233 fun add_deps (name, bnds)
3234 = Graph.add_deps_acyclic (name,
3235 remove (op =) name (Term.add_free_names (prop_of bnds) []))
3237 val order = Graph.empty
3238 |> fold Graph.new_node variable_bounds
3239 |> fold add_deps variable_bounds
3240 |> Graph.strong_conn |> map the_single |> rev
3241 |> map_filter (AList.lookup (op =) variable_bounds)
3243 fun prepend_prem th tac
3244 = tac THEN rtac (th RSN (2, @{thm mp})) i
3246 fold prepend_prem order all_tac
3249 (* Should be in HOL.thy ? *)
3250 fun gen_eval_tac conv ctxt = CONVERSION (Conv.params_conv (~1) (K (Conv.concl_conv (~1) conv)) ctxt)
3253 val form_equations = PureThy.get_thms @{theory} "interpret_form_equations";
3255 fun rewrite_interpret_form_tac ctxt prec splitting taylor i st = let
3256 fun lookup_splitting (Free (name, typ))
3257 = case AList.lookup (op =) splitting name
3258 of SOME s => HOLogic.mk_number @{typ nat} s
3259 | NONE => @{term "0 :: nat"}
3260 val vs = nth (prems_of st) (i - 1)
3261 |> Logic.strip_imp_concl
3262 |> HOLogic.dest_Trueprop
3263 |> Term.strip_comb |> snd |> List.last
3264 |> HOLogic.dest_list
3266 |> HOLogic.mk_number @{typ nat}
3267 |> Thm.cterm_of (ProofContext.theory_of ctxt)
3270 val n = vs |> length
3271 |> HOLogic.mk_number @{typ nat}
3272 |> Thm.cterm_of (ProofContext.theory_of ctxt)
3274 |> map lookup_splitting
3275 |> HOLogic.mk_list @{typ nat}
3276 |> Thm.cterm_of (ProofContext.theory_of ctxt)
3278 (rtac (Thm.instantiate ([], [(@{cpat "?n::nat"}, n),
3279 (@{cpat "?prec::nat"}, p),
3280 (@{cpat "?ss::nat list"}, s)])
3281 @{thm "approx_form"}) i
3282 THEN simp_tac @{simpset} i) st
3285 | SOME t => if length vs <> 1 then raise (TERM ("More than one variable used for taylor series expansion", [prop_of st]))
3288 |> HOLogic.mk_number @{typ nat}
3289 |> Thm.cterm_of (ProofContext.theory_of ctxt)
3290 val s = vs |> map lookup_splitting |> hd
3291 |> Thm.cterm_of (ProofContext.theory_of ctxt)
3293 rtac (Thm.instantiate ([], [(@{cpat "?s::nat"}, s),
3294 (@{cpat "?t::nat"}, t),
3295 (@{cpat "?prec::nat"}, p)])
3296 @{thm "approx_tse_form"}) i st
3300 (* copied from Tools/induct.ML should probably in args.ML *)
3301 val free = Args.context -- Args.term >> (fn (_, Free (n, t)) => n | (ctxt, t) =>
3302 error ("Bad free variable: " ^ Syntax.string_of_term ctxt t));
3306 lemma intervalE: "a \<le> x \<and> x \<le> b \<Longrightarrow> \<lbrakk> x \<in> { a .. b } \<Longrightarrow> P\<rbrakk> \<Longrightarrow> P"
3309 lemma meta_eqE: "x \<equiv> a \<Longrightarrow> \<lbrakk> x = a \<Longrightarrow> P\<rbrakk> \<Longrightarrow> P"
3312 method_setup approximation = {*
3313 Scan.lift (OuterParse.nat)
3315 Scan.optional (Scan.lift (Args.$$$ "splitting" |-- Args.colon)
3316 |-- OuterParse.and_list' (free --| Scan.lift (Args.$$$ "=") -- Scan.lift OuterParse.nat)) []
3318 Scan.option (Scan.lift (Args.$$$ "taylor" |-- Args.colon)
3319 |-- (free |-- Scan.lift (Args.$$$ "=") |-- Scan.lift OuterParse.nat))
3321 (fn ((prec, splitting), taylor) => fn ctxt =>
3322 SIMPLE_METHOD' (fn i =>
3323 REPEAT (FIRST' [etac @{thm intervalE},
3324 etac @{thm meta_eqE},
3325 rtac @{thm impI}] i)
3326 THEN Subgoal.FOCUS (fn {prems, ...} => reorder_bounds_tac prems i) @{context} i
3327 THEN DETERM (TRY (filter_prems_tac (K false) i))
3328 THEN DETERM (Reflection.genreify_tac ctxt form_equations NONE i)
3329 THEN rewrite_interpret_form_tac ctxt prec splitting taylor i
3330 THEN gen_eval_tac eval_oracle ctxt i))
3331 *} "real number approximation"
3334 fun calculated_subterms (@{const Trueprop} $ t) = calculated_subterms t
3335 | calculated_subterms (@{const "op -->"} $ _ $ t) = calculated_subterms t
3336 | calculated_subterms (@{term "op <= :: real \<Rightarrow> real \<Rightarrow> bool"} $ t1 $ t2) = [t1, t2]
3337 | calculated_subterms (@{term "op < :: real \<Rightarrow> real \<Rightarrow> bool"} $ t1 $ t2) = [t1, t2]
3338 | calculated_subterms (@{term "op : :: real \<Rightarrow> real set \<Rightarrow> bool"} $ t1 $
3339 (@{term "atLeastAtMost :: real \<Rightarrow> real \<Rightarrow> real set"} $ t2 $ t3)) = [t1, t2, t3]
3340 | calculated_subterms t = raise TERM ("calculated_subterms", [t])
3342 fun dest_interpret_form (@{const "interpret_form"} $ b $ xs) = (b, xs)
3343 | dest_interpret_form t = raise TERM ("dest_interpret_form", [t])
3345 fun dest_interpret (@{const "interpret_floatarith"} $ b $ xs) = (b, xs)
3346 | dest_interpret t = raise TERM ("dest_interpret", [t])
3349 fun dest_float (@{const "Float"} $ m $ e) = (snd (HOLogic.dest_number m), snd (HOLogic.dest_number e))
3350 fun dest_ivl (Const (@{const_name "Some"}, _) $
3351 (Const (@{const_name "Pair"}, _) $ u $ l)) = SOME (dest_float u, dest_float l)
3352 | dest_ivl (Const (@{const_name "None"}, _)) = NONE
3353 | dest_ivl t = raise TERM ("dest_result", [t])
3355 fun mk_approx' prec t = (@{const "approx'"}
3356 $ HOLogic.mk_number @{typ nat} prec
3357 $ t $ @{term "[] :: (float * float) option list"})
3359 fun mk_approx_form_eval prec t xs = (@{const "approx_form_eval"}
3360 $ HOLogic.mk_number @{typ nat} prec
3363 fun float2_float10 prec round_down (m, e) = (
3365 val (m, e) = (if e < 0 then (m,e) else (m * Integer.pow e 2, 0))
3367 fun frac c p 0 digits cnt = (digits, cnt, 0)
3368 | frac c 0 r digits cnt = (digits, cnt, r)
3369 | frac c p r digits cnt = (let
3370 val (d, r) = Integer.div_mod (r * 10) (Integer.pow (~e) 2)
3371 in frac (c orelse d <> 0) (if d <> 0 orelse c then p - 1 else p) r
3372 (digits * 10 + d) (cnt + 1)
3375 val sgn = Int.sign m
3378 val round_down = (sgn = 1 andalso round_down) orelse
3379 (sgn = ~1 andalso not round_down)
3381 val (x, r) = Integer.div_mod m (Integer.pow (~e) 2)
3383 val p = ((if x = 0 then prec else prec - (IntInf.log2 x + 1)) * 3) div 10 + 1
3385 val (digits, e10, r) = if p > 0 then frac (x <> 0) p r 0 0 else (0,0,0)
3387 val digits = if round_down orelse r = 0 then digits else digits + 1
3389 in (sgn * (digits + x * (Integer.pow e10 10)), ~e10)
3392 fun mk_result prec (SOME (l, u)) = (let
3393 fun mk_float10 rnd x = (let val (m, e) = float2_float10 prec rnd x
3394 in if e = 0 then HOLogic.mk_number @{typ real} m
3395 else if e = 1 then @{term "divide :: real \<Rightarrow> real \<Rightarrow> real"} $
3396 HOLogic.mk_number @{typ real} m $
3398 else @{term "divide :: real \<Rightarrow> real \<Rightarrow> real"} $
3399 HOLogic.mk_number @{typ real} m $
3400 (@{term "power 10 :: nat \<Rightarrow> real"} $
3401 HOLogic.mk_number @{typ nat} (~e)) end)
3402 in @{term "atLeastAtMost :: real \<Rightarrow> real \<Rightarrow> real set"} $ mk_float10 true l $ mk_float10 false u end)
3403 | mk_result prec NONE = @{term "UNIV :: real set"}
3406 val t = Logic.varify_global t
3407 val m = map (fn (name, sort) => (name, @{typ real})) (Term.add_tvars t [])
3408 val t = Term.subst_TVars m t
3411 fun converted_result t =
3413 |> HOLogic.dest_Trueprop
3414 |> HOLogic.dest_eq |> snd
3416 fun apply_tactic context term tactic = cterm_of context term
3419 |> the |> prems_of |> hd
3421 fun prepare_form context term = apply_tactic context term (
3422 REPEAT (FIRST' [etac @{thm intervalE}, etac @{thm meta_eqE}, rtac @{thm impI}] 1)
3423 THEN Subgoal.FOCUS (fn {prems, ...} => reorder_bounds_tac prems 1) @{context} 1
3424 THEN DETERM (TRY (filter_prems_tac (K false) 1)))
3426 fun reify_form context term = apply_tactic context term
3427 (Reflection.genreify_tac @{context} form_equations NONE 1)
3429 fun approx_form prec ctxt t =
3431 |> prepare_form (ProofContext.theory_of ctxt)
3432 |> (fn arith_term =>
3433 reify_form (ProofContext.theory_of ctxt) arith_term
3434 |> HOLogic.dest_Trueprop |> dest_interpret_form
3435 |> (fn (data, xs) =>
3436 mk_approx_form_eval prec data (HOLogic.mk_list @{typ "(float * float) option"}
3437 (map (fn _ => @{term "None :: (float * float) option"}) (HOLogic.dest_list xs)))
3438 |> Codegen.eval_term @{theory}
3439 |> HOLogic.dest_list
3440 |> curry ListPair.zip (HOLogic.dest_list xs @ calculated_subterms arith_term)
3441 |> map (fn (elem, s) => @{term "op : :: real \<Rightarrow> real set \<Rightarrow> bool"} $ elem $ mk_result prec (dest_ivl s))
3442 |> foldr1 HOLogic.mk_conj))
3444 fun approx_arith prec ctxt t = realify t
3445 |> Reflection.genreif ctxt form_equations
3447 |> HOLogic.dest_Trueprop
3448 |> HOLogic.dest_eq |> snd
3449 |> dest_interpret |> fst
3451 |> Codegen.eval_term @{theory}
3455 fun approx prec ctxt t = if type_of t = @{typ prop} then approx_form prec ctxt t
3456 else if type_of t = @{typ bool} then approx_form prec ctxt (@{const Trueprop} $ t)
3457 else approx_arith prec ctxt t
3461 Value.add_evaluator ("approximate", approx 30)