1 (* Author: Johannes Hoelzl, TU Muenchen
2 Coercions removed by Dmitriy Traytel *)
4 header {* Prove Real Valued Inequalities by Computation *}
9 "~~/src/HOL/Library/Float"
10 "~~/src/HOL/Library/Reflection"
11 "~~/src/HOL/Decision_Procs/Dense_Linear_Order"
12 "~~/src/HOL/Library/Efficient_Nat"
15 section "Horner Scheme"
17 subsection {* Define auxiliary helper @{text horner} function *}
19 primrec horner :: "(nat \<Rightarrow> nat) \<Rightarrow> (nat \<Rightarrow> nat \<Rightarrow> nat) \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> real \<Rightarrow> real" where
20 "horner F G 0 i k x = 0" |
21 "horner F G (Suc n) i k x = 1 / k - x * horner F G n (F i) (G i k) x"
23 lemma horner_schema': fixes x :: real and a :: "nat \<Rightarrow> real"
24 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)"
26 have shift_pow: "\<And>i. - (x * ((-1)^i * a (Suc i) * x ^ i)) = (-1)^(Suc i) * a (Suc i) * x ^ (Suc i)" by auto
27 show ?thesis unfolding setsum_right_distrib shift_pow diff_minus setsum_negf[symmetric] setsum_head_upt_Suc[OF zero_less_Suc]
28 setsum_reindex[OF inj_Suc, unfolded comp_def, symmetric, of "\<lambda> n. (-1)^n *a n * x^n"] by auto
31 lemma horner_schema: fixes f :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat" and F :: "nat \<Rightarrow> nat"
32 assumes f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"
33 shows "horner F G n ((F ^^ j') s) (f j') x = (\<Sum> j = 0..< n. -1 ^ j * (1 / (f (j' + j))) * x ^ j)"
34 proof (induct n arbitrary: i k j')
37 show ?case unfolding horner.simps Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc]
38 using horner_schema'[of "\<lambda> j. 1 / (f (j' + j))"] by auto
42 fixes lb :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" and ub :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
43 assumes "0 \<le> real x" and f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"
44 and lb_0: "\<And> i k x. lb 0 i k x = 0"
45 and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 k - x * (ub n (F i) (G i k) x)"
46 and ub_0: "\<And> i k x. ub 0 i k x = 0"
47 and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 k - x * (lb n (F i) (G i k) x)"
48 shows "(lb n ((F ^^ j') s) (f j') x) \<le> horner F G n ((F ^^ j') s) (f j') x \<and>
49 horner F G n ((F ^^ j') s) (f j') x \<le> (ub n ((F ^^ j') s) (f j') x)"
50 (is "?lb n j' \<le> ?horner n j' \<and> ?horner n j' \<le> ?ub n j'")
51 proof (induct n arbitrary: j')
52 case 0 thus ?case unfolding lb_0 ub_0 horner.simps by auto
55 have "?lb (Suc n) j' \<le> ?horner (Suc n) j'" unfolding lb_Suc ub_Suc horner.simps real_of_float_sub diff_minus
57 show "(lapprox_rat prec 1 (f j')) \<le> 1 / (f j')" using lapprox_rat[of prec 1 "f j'"] by auto
58 from Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc, THEN conjunct2] `0 \<le> real x`
59 show "- real (x * ub n (F ((F ^^ j') s)) (G ((F ^^ j') s) (f j')) x) \<le>
60 - (x * horner F G 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 moreover have "?horner (Suc n) j' \<le> ?ub (Suc n) j'" unfolding ub_Suc ub_Suc horner.simps real_of_float_sub diff_minus
65 show "1 / (f j') \<le> (rapprox_rat prec 1 (f j'))" using rapprox_rat[of 1 "f j'" prec] by auto
66 from Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc, THEN conjunct1] `0 \<le> real x`
67 show "- (x * horner F G n (F ((F ^^ j') s)) (G ((F ^^ j') s) (f j')) x) \<le>
68 - real (x * lb n (F ((F ^^ j') s)) (G ((F ^^ j') s) (f j')) x)"
69 unfolding real_of_float_mult neg_le_iff_le by (rule mult_left_mono)
71 ultimately show ?case by blast
74 subsection "Theorems for floating point functions implementing the horner scheme"
78 Here @{term_type "f :: nat \<Rightarrow> nat"} is the sequence defining the Taylor series, the coefficients are
79 all alternating and reciprocs. We use @{term G} and @{term F} to describe the computation of @{term f}.
83 lemma horner_bounds: fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat"
84 assumes "0 \<le> real x" and f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"
85 and lb_0: "\<And> i k x. lb 0 i k x = 0"
86 and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 k - x * (ub n (F i) (G i k) x)"
87 and ub_0: "\<And> i k x. ub 0 i k x = 0"
88 and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 k - x * (lb n (F i) (G i k) x)"
89 shows "(lb n ((F ^^ j') s) (f j') x) \<le> (\<Sum>j=0..<n. -1 ^ j * (1 / (f (j' + j))) * (x ^ j))" (is "?lb") and
90 "(\<Sum>j=0..<n. -1 ^ j * (1 / (f (j' + j))) * (x ^ j)) \<le> (ub n ((F ^^ j') s) (f j') x)" (is "?ub")
93 using horner_bounds'[where lb=lb, OF `0 \<le> real x` f_Suc lb_0 lb_Suc ub_0 ub_Suc]
94 unfolding horner_schema[where f=f, OF f_Suc] .
95 thus "?lb" and "?ub" by auto
98 lemma horner_bounds_nonpos: fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat"
99 assumes "real x \<le> 0" and f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"
100 and lb_0: "\<And> i k x. lb 0 i k x = 0"
101 and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 k + x * (ub n (F i) (G i k) x)"
102 and ub_0: "\<And> i k x. ub 0 i k x = 0"
103 and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 k + x * (lb n (F i) (G i k) x)"
104 shows "(lb n ((F ^^ j') s) (f j') x) \<le> (\<Sum>j=0..<n. (1 / (f (j' + j))) * real x ^ j)" (is "?lb") and
105 "(\<Sum>j=0..<n. (1 / (f (j' + j))) * real x ^ j) \<le> (ub n ((F ^^ j') s) (f j') x)" (is "?ub")
107 { fix x y z :: float have "x - y * z = x + - y * z"
108 by (cases x, cases y, cases z, simp add: plus_float.simps minus_float_def uminus_float.simps times_float.simps algebra_simps)
109 } note diff_mult_minus = this
111 { fix x :: float have "- (- x) = x" by (cases x, auto simp add: uminus_float.simps) } note minus_minus = this
113 have move_minus: "(-x) = -1 * real x" by auto (* coercion "inside" is necessary *)
115 have sum_eq: "(\<Sum>j=0..<n. (1 / (f (j' + j))) * real x ^ j) =
116 (\<Sum>j = 0..<n. -1 ^ j * (1 / (f (j' + j))) * real (- x) ^ j)"
117 proof (rule setsum_cong, simp)
118 fix j assume "j \<in> {0 ..< n}"
119 show "1 / (f (j' + j)) * real x ^ j = -1 ^ j * (1 / (f (j' + j))) * real (- x) ^ j"
120 unfolding move_minus power_mult_distrib mult_assoc[symmetric]
121 unfolding mult_commute unfolding mult_assoc[of "-1 ^ j", symmetric] power_mult_distrib[symmetric]
125 have "0 \<le> real (-x)" using assms by auto
126 from horner_bounds[where G=G and F=F and f=f and s=s and prec=prec
127 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,
128 OF this f_Suc lb_0 refl ub_0 refl]
129 show "?lb" and "?ub" unfolding minus_minus sum_eq
133 subsection {* Selectors for next even or odd number *}
137 The horner scheme computes alternating series. To get the upper and lower bounds we need to
138 guarantee to access a even or odd member. To do this we use @{term get_odd} and @{term get_even}.
142 definition get_odd :: "nat \<Rightarrow> nat" where
143 "get_odd n = (if odd n then n else (Suc n))"
145 definition get_even :: "nat \<Rightarrow> nat" where
146 "get_even n = (if even n then n else (Suc n))"
148 lemma get_odd[simp]: "odd (get_odd n)" unfolding get_odd_def by (cases "odd n", auto)
149 lemma get_even[simp]: "even (get_even n)" unfolding get_even_def by (cases "even n", auto)
150 lemma get_odd_ex: "\<exists> k. Suc k = get_odd n \<and> odd (Suc k)"
151 proof (cases "odd n")
152 case True hence "0 < n" by (rule odd_pos)
153 from gr0_implies_Suc[OF this] obtain k where "Suc k = n" by auto
154 thus ?thesis unfolding get_odd_def if_P[OF True] using True[unfolded `Suc k = n`[symmetric]] by blast
156 case False hence "odd (Suc n)" by auto
157 thus ?thesis unfolding get_odd_def if_not_P[OF False] by blast
160 lemma get_even_double: "\<exists>i. get_even n = 2 * i" using get_even[unfolded even_mult_two_ex] .
161 lemma get_odd_double: "\<exists>i. get_odd n = 2 * i + 1" using get_odd[unfolded odd_Suc_mult_two_ex] by auto
163 section "Power function"
165 definition float_power_bnds :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where
166 "float_power_bnds n l u = (if odd n \<or> 0 < l then (l ^ n, u ^ n)
167 else if u < 0 then (u ^ n, l ^ n)
168 else (0, (max (-l) u) ^ n))"
170 lemma float_power_bnds: fixes x :: real
171 assumes "(l1, u1) = float_power_bnds n l u" and "x \<in> {l .. u}"
172 shows "x ^ n \<in> {l1..u1}"
173 proof (cases "even n")
176 proof (cases "0 < l")
177 case True hence "odd n \<or> 0 < l" and "0 \<le> real l" unfolding less_float_def by auto
178 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
179 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 l x] power_mono[of x u] by auto
180 thus ?thesis using assms `0 < l` unfolding atLeastAtMost_iff l1 u1 float_power less_float_def by auto
182 case False hence P: "\<not> (odd n \<or> 0 < l)" using `even n` by auto
184 proof (cases "u < 0")
185 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
186 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]
187 unfolding power_minus_even[OF `even n`] by auto
188 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
189 ultimately show ?thesis using float_power by auto
192 have "\<bar>x\<bar> \<le> real (max (-l) u)"
193 proof (cases "-l \<le> u")
194 case True thus ?thesis unfolding max_def if_P[OF True] using assms unfolding le_float_def by auto
196 case False thus ?thesis unfolding max_def if_not_P[OF False] using assms unfolding le_float_def by auto
198 hence x_abs: "\<bar>x\<bar> \<le> \<bar>real (max (-l) u)\<bar>" by auto
199 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
200 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
204 case False hence "odd n \<or> 0 < l" by auto
205 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
206 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
207 thus ?thesis unfolding atLeastAtMost_iff l1 u1 float_power less_float_def by auto
210 lemma bnds_power: "\<forall> (x::real) l u. (l1, u1) = float_power_bnds n l u \<and> x \<in> {l .. u} \<longrightarrow> l1 \<le> x ^ n \<and> x ^ n \<le> u1"
211 using float_power_bnds by auto
213 section "Square root"
217 The square root computation is implemented as newton iteration. As first first step we use the
218 nearest power of two greater than the square root.
222 fun sqrt_iteration :: "nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
223 "sqrt_iteration prec 0 (Float m e) = Float 1 ((e + bitlen m) div 2 + 1)" |
224 "sqrt_iteration prec (Suc m) x = (let y = sqrt_iteration prec m x
225 in Float 1 -1 * (y + float_divr prec x y))"
227 function ub_sqrt lb_sqrt :: "nat \<Rightarrow> float \<Rightarrow> float" where
228 "ub_sqrt prec x = (if 0 < x then (sqrt_iteration prec prec x)
229 else if x < 0 then - lb_sqrt prec (- x)
231 "lb_sqrt prec x = (if 0 < x then (float_divl prec x (sqrt_iteration prec prec x))
232 else if x < 0 then - ub_sqrt prec (- x)
234 by pat_completeness auto
235 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)
237 declare lb_sqrt.simps[simp del]
238 declare ub_sqrt.simps[simp del]
240 lemma sqrt_ub_pos_pos_1:
241 assumes "sqrt x < b" and "0 < b" and "0 < x"
242 shows "sqrt x < (b + x / b)/2"
244 from assms have "0 < (b - sqrt x) ^ 2 " by simp
245 also have "\<dots> = b ^ 2 - 2 * b * sqrt x + (sqrt x) ^ 2" by algebra
246 also have "\<dots> = b ^ 2 - 2 * b * sqrt x + x" using assms by (simp add: real_sqrt_pow2)
247 finally have "0 < b ^ 2 - 2 * b * sqrt x + x" by assumption
248 hence "0 < b / 2 - sqrt x + x / (2 * b)" using assms
249 by (simp add: field_simps power2_eq_square)
250 thus ?thesis by (simp add: field_simps)
253 lemma sqrt_iteration_bound: assumes "0 < real x"
254 shows "sqrt x < (sqrt_iteration prec n x)"
260 hence "0 < m" using float_pos_m_pos[unfolded less_float_def] assms by auto
261 hence "0 < sqrt m" by auto
263 have int_nat_bl: "(nat (bitlen m)) = bitlen m" using bitlen_ge0 by auto
265 have "x = (m / 2^nat (bitlen m)) * pow2 (e + (nat (bitlen m)))"
266 unfolding pow2_add pow2_int Float real_of_float_simp by auto
267 also have "\<dots> < 1 * pow2 (e + nat (bitlen m))"
268 proof (rule mult_strict_right_mono, auto)
269 show "real m < 2^nat (bitlen m)" using bitlen_bounds[OF `0 < m`, THEN conjunct2]
270 unfolding real_of_int_less_iff[of m, symmetric] by auto
272 finally have "sqrt x < sqrt (pow2 (e + bitlen m))" unfolding int_nat_bl by auto
273 also have "\<dots> \<le> pow2 ((e + bitlen m) div 2 + 1)"
275 let ?E = "e + bitlen m"
276 have E_mod_pow: "pow2 (?E mod 2) < 4"
277 proof (cases "?E mod 2 = 1")
278 case True thus ?thesis by auto
281 have "0 \<le> ?E mod 2" by auto
282 have "?E mod 2 < 2" by auto
283 from this[THEN zless_imp_add1_zle]
284 have "?E mod 2 \<le> 0" using False by auto
285 from xt1(5)[OF `0 \<le> ?E mod 2` this]
288 hence "sqrt (pow2 (?E mod 2)) < sqrt (2 * 2)" by auto
289 hence E_mod_pow: "sqrt (pow2 (?E mod 2)) < 2" unfolding real_sqrt_abs2 by auto
291 have E_eq: "pow2 ?E = pow2 (?E div 2 + ?E div 2 + ?E mod 2)" by auto
292 have "sqrt (pow2 ?E) = sqrt (pow2 (?E div 2) * pow2 (?E div 2) * pow2 (?E mod 2))"
293 unfolding E_eq unfolding pow2_add ..
294 also have "\<dots> = pow2 (?E div 2) * sqrt (pow2 (?E mod 2))"
295 unfolding real_sqrt_mult[of _ "pow2 (?E mod 2)"] real_sqrt_abs2 by auto
296 also have "\<dots> < pow2 (?E div 2) * 2"
297 by (rule mult_strict_left_mono, auto intro: E_mod_pow)
298 also have "\<dots> = pow2 (?E div 2 + 1)" unfolding zadd_commute[of _ 1] pow2_add1 by auto
299 finally show ?thesis by auto
302 unfolding Float sqrt_iteration.simps real_of_float_simp by auto
306 let ?b = "sqrt_iteration prec n x"
307 have "0 < sqrt x" using `0 < real x` by auto
308 also have "\<dots> < real ?b" using Suc .
309 finally have "sqrt x < (?b + x / ?b)/2" using sqrt_ub_pos_pos_1[OF Suc _ `0 < real x`] by auto
310 also have "\<dots> \<le> (?b + (float_divr prec x ?b))/2" by (rule divide_right_mono, auto simp add: float_divr)
311 also have "\<dots> = (Float 1 -1) * (?b + (float_divr prec x ?b))" by auto
312 finally show ?case unfolding sqrt_iteration.simps Let_def real_of_float_mult real_of_float_add right_distrib .
315 lemma sqrt_iteration_lower_bound: assumes "0 < real x"
316 shows "0 < real (sqrt_iteration prec n x)" (is "0 < ?sqrt")
318 have "0 < sqrt x" using assms by auto
319 also have "\<dots> < ?sqrt" using sqrt_iteration_bound[OF assms] .
320 finally show ?thesis .
323 lemma lb_sqrt_lower_bound: assumes "0 \<le> real x"
324 shows "0 \<le> real (lb_sqrt prec x)"
325 proof (cases "0 < x")
326 case True hence "0 < real x" and "0 \<le> x" using `0 \<le> real x` unfolding less_float_def le_float_def by auto
327 hence "0 < sqrt_iteration prec prec x" unfolding less_float_def using sqrt_iteration_lower_bound by auto
328 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
329 thus ?thesis unfolding lb_sqrt.simps using True by auto
331 case False with `0 \<le> real x` have "real x = 0" unfolding less_float_def by auto
332 thus ?thesis unfolding lb_sqrt.simps less_float_def by auto
336 shows "sqrt x \<in> {(lb_sqrt prec x) .. (ub_sqrt prec x) }"
338 { fix x :: float assume "0 < x"
339 hence "0 < real x" and "0 \<le> real x" unfolding less_float_def by auto
340 hence sqrt_gt0: "0 < sqrt x" by auto
341 hence sqrt_ub: "sqrt x < sqrt_iteration prec prec x" using sqrt_iteration_bound by auto
343 have "(float_divl prec x (sqrt_iteration prec prec x)) \<le>
344 x / (sqrt_iteration prec prec x)" by (rule float_divl)
345 also have "\<dots> < x / sqrt x"
346 by (rule divide_strict_left_mono[OF sqrt_ub `0 < real x`
347 mult_pos_pos[OF order_less_trans[OF sqrt_gt0 sqrt_ub] sqrt_gt0]])
348 also have "\<dots> = sqrt x"
349 unfolding inverse_eq_iff_eq[of _ "sqrt x", symmetric]
350 sqrt_divide_self_eq[OF `0 \<le> real x`, symmetric] by auto
351 finally have "lb_sqrt prec x \<le> sqrt x"
352 unfolding lb_sqrt.simps if_P[OF `0 < x`] by auto }
355 { fix x :: float assume "0 < x"
356 hence "0 < real x" unfolding less_float_def by auto
357 hence "0 < sqrt x" by auto
358 hence "sqrt x < sqrt_iteration prec prec x"
359 using sqrt_iteration_bound by auto
360 hence "sqrt x \<le> ub_sqrt prec x"
361 unfolding ub_sqrt.simps if_P[OF `0 < x`] by auto }
365 proof (cases "0 < x")
366 case True with lb ub show ?thesis by auto
367 next case False show ?thesis
368 proof (cases "real x = 0")
369 case True thus ?thesis
370 by (auto simp add: less_float_def lb_sqrt.simps ub_sqrt.simps)
372 case False with `\<not> 0 < x` have "x < 0" and "0 < -x"
373 by (auto simp add: less_float_def)
376 show ?thesis using lb[OF `0 < -x`] ub[OF `0 < -x`]
377 by (auto simp add: real_sqrt_minus lb_sqrt.simps ub_sqrt.simps)
381 lemma bnds_sqrt: "\<forall> (x::real) lx ux. (l, u) = (lb_sqrt prec lx, ub_sqrt prec ux) \<and> x \<in> {lx .. ux} \<longrightarrow> l \<le> sqrt x \<and> sqrt x \<le> u"
382 proof ((rule allI) +, rule impI, erule conjE, rule conjI)
383 fix x :: real fix lx ux
384 assume "(l, u) = (lb_sqrt prec lx, ub_sqrt prec ux)"
385 and x: "x \<in> {lx .. ux}"
386 hence l: "l = lb_sqrt prec lx " and u: "u = ub_sqrt prec ux" by auto
388 have "sqrt lx \<le> sqrt x" using x by auto
389 from order_trans[OF _ this]
390 show "l \<le> sqrt x" unfolding l using bnds_sqrt'[of lx prec] by auto
392 have "sqrt x \<le> sqrt ux" using x by auto
393 from order_trans[OF this]
394 show "sqrt x \<le> u" unfolding u using bnds_sqrt'[of ux prec] by auto
397 section "Arcus tangens and \<pi>"
399 subsection "Compute arcus tangens series"
403 As first step we implement the computation of the arcus tangens series. This is only valid in the range
404 @{term "{-1 :: real .. 1}"}. This is used to compute \<pi> and then the entire arcus tangens.
408 fun ub_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
409 and lb_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
410 "ub_arctan_horner prec 0 k x = 0"
411 | "ub_arctan_horner prec (Suc n) k x =
412 (rapprox_rat prec 1 k) - x * (lb_arctan_horner prec n (k + 2) x)"
413 | "lb_arctan_horner prec 0 k x = 0"
414 | "lb_arctan_horner prec (Suc n) k x =
415 (lapprox_rat prec 1 k) - x * (ub_arctan_horner prec n (k + 2) x)"
417 lemma arctan_0_1_bounds': assumes "0 \<le> real x" "real x \<le> 1" and "even n"
418 shows "arctan x \<in> {(x * lb_arctan_horner prec n 1 (x * x)) .. (x * ub_arctan_horner prec (Suc n) 1 (x * x))}"
420 let "?c i" = "-1^i * (1 / (i * 2 + (1::nat)) * real x ^ (i * 2 + 1))"
421 let "?S n" = "\<Sum> i=0..<n. ?c i"
423 have "0 \<le> real (x * x)" by auto
424 from `even n` obtain m where "2 * m = n" unfolding even_mult_two_ex by auto
426 have "arctan x \<in> { ?S n .. ?S (Suc n) }"
427 proof (cases "real x = 0")
429 hence "0 < real x" using `0 \<le> real x` by auto
430 hence prem: "0 < 1 / (0 * 2 + (1::nat)) * real x ^ (0 * 2 + 1)" by auto
432 have "\<bar> real x \<bar> \<le> 1" using `0 \<le> real x` `real x \<le> 1` by auto
433 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`]
434 show ?thesis unfolding arctan_series[OF `\<bar> real x \<bar> \<le> 1`] Suc_eq_plus1 .
436 note arctan_bounds = this[unfolded atLeastAtMost_iff]
438 have F: "\<And>n. 2 * Suc n + 1 = 2 * n + 1 + 2" by auto
440 note bounds = horner_bounds[where s=1 and f="\<lambda>i. 2 * i + 1" and j'=0
441 and lb="\<lambda>n i k x. lb_arctan_horner prec n k x"
442 and ub="\<lambda>n i k x. ub_arctan_horner prec n k x",
443 OF `0 \<le> real (x*x)` F lb_arctan_horner.simps ub_arctan_horner.simps]
445 { have "(x * lb_arctan_horner prec n 1 (x*x)) \<le> ?S n"
446 using bounds(1) `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 also have "\<dots> \<le> arctan x" using arctan_bounds ..
451 finally have "(x * lb_arctan_horner prec n 1 (x*x)) \<le> arctan x" . }
453 { have "arctan x \<le> ?S (Suc n)" using arctan_bounds ..
454 also have "\<dots> \<le> (x * ub_arctan_horner prec (Suc n) 1 (x*x))"
455 using bounds(2)[of "Suc n"] `0 \<le> real x`
456 unfolding real_of_float_mult power_add power_one_right mult_assoc[symmetric] setsum_left_distrib[symmetric]
457 unfolding mult_commute[where 'a=real] mult_commute[of _ "2::nat"] power_mult power2_eq_square[of "real x"]
458 by (auto intro!: mult_left_mono)
459 finally have "arctan x \<le> (x * ub_arctan_horner prec (Suc n) 1 (x*x))" . }
460 ultimately show ?thesis by auto
463 lemma arctan_0_1_bounds: assumes "0 \<le> real x" "real x \<le> 1"
464 shows "arctan x \<in> {(x * lb_arctan_horner prec (get_even n) 1 (x * x)) .. (x * ub_arctan_horner prec (get_odd n) 1 (x * x))}"
465 proof (cases "even n")
467 obtain n' where "Suc n' = get_odd n" and "odd (Suc n')" using get_odd_ex by auto
468 hence "even n'" unfolding even_Suc by auto
469 have "arctan x \<le> x * ub_arctan_horner prec (get_odd n) 1 (x * x)"
470 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
472 have "x * lb_arctan_horner prec (get_even n) 1 (x * x) \<le> arctan x"
473 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
474 ultimately show ?thesis by auto
476 case False hence "0 < n" by (rule odd_pos)
477 from gr0_implies_Suc[OF this] obtain n' where "n = Suc n'" ..
478 from False[unfolded this even_Suc]
479 have "even n'" and "even (Suc (Suc n'))" by auto
480 have "get_odd n = Suc n'" unfolding get_odd_def if_P[OF False] using `n = Suc n'` .
482 have "arctan x \<le> x * ub_arctan_horner prec (get_odd n) 1 (x * x)"
483 unfolding `get_odd n = Suc n'` using arctan_0_1_bounds'[OF `0 \<le> real x` `real x \<le> 1` `even n'`] by auto
485 have "(x * lb_arctan_horner prec (get_even n) 1 (x * x)) \<le> arctan x"
486 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
487 ultimately show ?thesis by auto
490 subsection "Compute \<pi>"
492 definition ub_pi :: "nat \<Rightarrow> float" where
493 "ub_pi prec = (let A = rapprox_rat prec 1 5 ;
494 B = lapprox_rat prec 1 239
495 in ((Float 1 2) * ((Float 1 2) * A * (ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (A * A)) -
496 B * (lb_arctan_horner prec (get_even (prec div 14 + 1)) 1 (B * B)))))"
498 definition lb_pi :: "nat \<Rightarrow> float" where
499 "lb_pi prec = (let A = lapprox_rat prec 1 5 ;
500 B = rapprox_rat prec 1 239
501 in ((Float 1 2) * ((Float 1 2) * A * (lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (A * A)) -
502 B * (ub_arctan_horner prec (get_odd (prec div 14 + 1)) 1 (B * B)))))"
504 lemma pi_boundaries: "pi \<in> {(lb_pi n) .. (ub_pi n)}"
506 have machin_pi: "pi = 4 * (4 * arctan (1 / 5) - arctan (1 / 239))" unfolding machin[symmetric] by auto
508 { fix prec n :: nat fix k :: int assume "1 < k" hence "0 \<le> k" and "0 < k" and "1 \<le> k" by auto
509 let ?k = "rapprox_rat prec 1 k"
510 have "1 div k = 0" using div_pos_pos_trivial[OF _ `1 < k`] by auto
512 have "0 \<le> real ?k" by (rule order_trans[OF _ rapprox_rat], auto simp add: `0 \<le> k`)
513 have "real ?k \<le> 1" unfolding rapprox_rat.simps(2)[OF zero_le_one `0 < k`]
514 by (rule rapprox_posrat_le1, auto simp add: `0 < k` `1 \<le> k`)
516 have "1 / k \<le> ?k" using rapprox_rat[where x=1 and y=k] by auto
517 hence "arctan (1 / k) \<le> arctan ?k" by (rule arctan_monotone')
518 also have "\<dots> \<le> (?k * ub_arctan_horner prec (get_odd n) 1 (?k * ?k))"
519 using arctan_0_1_bounds[OF `0 \<le> real ?k` `real ?k \<le> 1`] by auto
520 finally have "arctan (1 / k) \<le> ?k * ub_arctan_horner prec (get_odd n) 1 (?k * ?k)" .
521 } note ub_arctan = this
523 { fix prec n :: nat fix k :: int assume "1 < k" hence "0 \<le> k" and "0 < k" by auto
524 let ?k = "lapprox_rat prec 1 k"
525 have "1 div k = 0" using div_pos_pos_trivial[OF _ `1 < k`] by auto
526 have "1 / k \<le> 1" using `1 < k` by auto
528 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`)
529 have "\<And>n. real ?k \<le> 1" using lapprox_rat by (rule order_trans, auto simp add: `1 / k \<le> 1`)
531 have "?k \<le> 1 / k" using lapprox_rat[where x=1 and y=k] by auto
533 have "?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k) \<le> arctan ?k"
534 using arctan_0_1_bounds[OF `0 \<le> real ?k` `real ?k \<le> 1`] by auto
535 also have "\<dots> \<le> arctan (1 / k)" using `?k \<le> 1 / k` by (rule arctan_monotone')
536 finally have "?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k) \<le> arctan (1 / k)" .
537 } note lb_arctan = this
539 have "pi \<le> ub_pi n"
540 unfolding ub_pi_def machin_pi Let_def real_of_float_mult real_of_float_sub unfolding Float_num
541 using lb_arctan[of 239] ub_arctan[of 5]
542 by (auto intro!: mult_left_mono add_mono simp add: diff_minus simp del: lapprox_rat.simps rapprox_rat.simps)
544 have "lb_pi n \<le> pi"
545 unfolding lb_pi_def machin_pi Let_def real_of_float_mult real_of_float_sub Float_num
546 using lb_arctan[of 5] ub_arctan[of 239]
547 by (auto intro!: mult_left_mono add_mono simp add: diff_minus simp del: lapprox_rat.simps rapprox_rat.simps)
548 ultimately show ?thesis by auto
551 subsection "Compute arcus tangens in the entire domain"
553 function lb_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" and ub_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" where
554 "lb_arctan prec x = (let ub_horner = \<lambda> x. x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x) ;
555 lb_horner = \<lambda> x. x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)
556 in (if x < 0 then - ub_arctan prec (-x) else
557 if x \<le> Float 1 -1 then lb_horner x else
558 if x \<le> Float 1 1 then Float 1 1 * lb_horner (float_divl prec x (1 + ub_sqrt prec (1 + x * x)))
559 else (let inv = float_divr prec 1 x
561 else lb_pi prec * Float 1 -1 - ub_horner inv)))"
563 | "ub_arctan prec x = (let lb_horner = \<lambda> x. x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x) ;
564 ub_horner = \<lambda> x. x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)
565 in (if x < 0 then - lb_arctan prec (-x) else
566 if x \<le> Float 1 -1 then ub_horner x else
567 if x \<le> Float 1 1 then let y = float_divr prec x (1 + lb_sqrt prec (1 + x * x))
568 in if y > 1 then ub_pi prec * Float 1 -1
569 else Float 1 1 * ub_horner y
570 else ub_pi prec * Float 1 -1 - lb_horner (float_divl prec 1 x)))"
571 by pat_completeness auto
572 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)
574 declare ub_arctan_horner.simps[simp del]
575 declare lb_arctan_horner.simps[simp del]
577 lemma lb_arctan_bound': assumes "0 \<le> real x"
578 shows "lb_arctan prec x \<le> arctan x"
580 have "\<not> x < 0" and "0 \<le> x" unfolding less_float_def le_float_def using `0 \<le> real x` by auto
581 let "?ub_horner x" = "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)"
582 and "?lb_horner x" = "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)"
585 proof (cases "x \<le> Float 1 -1")
586 case True hence "real x \<le> 1" unfolding le_float_def Float_num by auto
587 show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True]
588 using arctan_0_1_bounds[OF `0 \<le> real x` `real x \<le> 1`] by auto
590 case False hence "0 < real x" unfolding le_float_def Float_num by auto
591 let ?R = "1 + sqrt (1 + real x * real x)"
592 let ?fR = "1 + ub_sqrt prec (1 + x * x)"
593 let ?DIV = "float_divl prec x ?fR"
595 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
596 hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg)
598 have "sqrt (1 + x * x) \<le> ub_sqrt prec (1 + x * x)"
599 using bnds_sqrt'[of "1 + x * x"] by auto
601 hence "?R \<le> ?fR" by auto
602 hence "0 < ?fR" and "0 < real ?fR" unfolding less_float_def using `0 < ?R` by auto
604 have monotone: "(float_divl prec x ?fR) \<le> x / ?R"
606 have "?DIV \<le> real x / ?fR" by (rule float_divl)
607 also have "\<dots> \<le> x / ?R" by (rule divide_left_mono[OF `?R \<le> ?fR` `0 \<le> real x` mult_pos_pos[OF order_less_le_trans[OF divisor_gt0 `?R \<le> real ?fR`] divisor_gt0]])
608 finally show ?thesis .
612 proof (cases "x \<le> Float 1 1")
615 have "x \<le> sqrt (1 + x * x)" using real_sqrt_sum_squares_ge2[where x=1, unfolded numeral_2_eq_2] by auto
616 also have "\<dots> \<le> (ub_sqrt prec (1 + x * x))"
617 using bnds_sqrt'[of "1 + x * x"] by auto
618 finally have "real x \<le> ?fR" by auto
619 moreover have "?DIV \<le> real x / ?fR" by (rule float_divl)
620 ultimately have "real ?DIV \<le> 1" unfolding divide_le_eq_1_pos[OF `0 < real ?fR`, symmetric] by auto
622 have "0 \<le> real ?DIV" using float_divl_lower_bound[OF `0 \<le> x` `0 < ?fR`] unfolding le_float_def by auto
624 have "(Float 1 1 * ?lb_horner ?DIV) \<le> 2 * arctan (float_divl prec x ?fR)" unfolding real_of_float_mult[of "Float 1 1"] Float_num
625 using arctan_0_1_bounds[OF `0 \<le> real ?DIV` `real ?DIV \<le> 1`] by auto
626 also have "\<dots> \<le> 2 * arctan (x / ?R)"
627 using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono)
628 also have "2 * arctan (x / ?R) = arctan x" using arctan_half[symmetric] unfolding numeral_2_eq_2 power_Suc2 power_0 mult_1_left .
629 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] .
632 hence "2 < real x" unfolding le_float_def Float_num by auto
633 hence "1 \<le> real x" by auto
635 let "?invx" = "float_divr prec 1 x"
636 have "0 \<le> arctan x" using arctan_monotone'[OF `0 \<le> real x`] using arctan_tan[of 0, unfolded tan_zero] by auto
639 proof (cases "1 < ?invx")
641 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]
642 using `0 \<le> arctan x` by auto
645 hence "real ?invx \<le> 1" unfolding less_float_def by auto
646 have "0 \<le> real ?invx" by (rule order_trans[OF _ float_divr], auto simp add: `0 \<le> real x`)
648 have "1 / x \<noteq> 0" and "0 < 1 / x" using `0 < real x` by auto
650 have "arctan (1 / x) \<le> arctan ?invx" unfolding real_of_float_1[symmetric] by (rule arctan_monotone', rule float_divr)
651 also have "\<dots> \<le> (?ub_horner ?invx)" using arctan_0_1_bounds[OF `0 \<le> real ?invx` `real ?invx \<le> 1`] by auto
652 finally have "pi / 2 - (?ub_horner ?invx) \<le> arctan x"
653 using `0 \<le> arctan x` arctan_inverse[OF `1 / x \<noteq> 0`]
654 unfolding real_sgn_pos[OF `0 < 1 / real x`] le_diff_eq by auto
656 have "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
658 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]
665 lemma ub_arctan_bound': assumes "0 \<le> real x"
666 shows "arctan x \<le> ub_arctan prec x"
668 have "\<not> x < 0" and "0 \<le> x" unfolding less_float_def le_float_def using `0 \<le> real x` by auto
670 let "?ub_horner x" = "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)"
671 and "?lb_horner x" = "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)"
674 proof (cases "x \<le> Float 1 -1")
675 case True hence "real x \<le> 1" unfolding le_float_def Float_num by auto
676 show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True]
677 using arctan_0_1_bounds[OF `0 \<le> real x` `real x \<le> 1`] by auto
679 case False hence "0 < real x" unfolding le_float_def Float_num by auto
680 let ?R = "1 + sqrt (1 + real x * real x)"
681 let ?fR = "1 + lb_sqrt prec (1 + x * x)"
682 let ?DIV = "float_divr prec x ?fR"
684 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
685 hence "0 \<le> real (1 + x*x)" by auto
687 hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg)
689 have "lb_sqrt prec (1 + x * x) \<le> sqrt (1 + x * x)"
690 using bnds_sqrt'[of "1 + x * x"] by auto
691 hence "?fR \<le> ?R" by auto
692 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)`])
694 have monotone: "x / ?R \<le> (float_divr prec x ?fR)"
696 from divide_left_mono[OF `?fR \<le> ?R` `0 \<le> real x` mult_pos_pos[OF divisor_gt0 `0 < real ?fR`]]
697 have "x / ?R \<le> x / ?fR" .
698 also have "\<dots> \<le> ?DIV" by (rule float_divr)
699 finally show ?thesis .
703 proof (cases "x \<le> Float 1 1")
706 proof (cases "?DIV > 1")
708 have "pi / 2 \<le> 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
709 from order_less_le_trans[OF arctan_ubound this, THEN less_imp_le]
710 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] .
713 hence "real ?DIV \<le> 1" unfolding less_float_def by auto
715 have "0 \<le> x / ?R" using `0 \<le> real x` `0 < ?R` unfolding real_0_le_divide_iff by auto
716 hence "0 \<le> real ?DIV" using monotone by (rule order_trans)
718 have "arctan x = 2 * arctan (x / ?R)" using arctan_half unfolding numeral_2_eq_2 power_Suc2 power_0 mult_1_left .
719 also have "\<dots> \<le> 2 * arctan (?DIV)"
720 using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono)
721 also have "\<dots> \<le> (Float 1 1 * ?ub_horner ?DIV)" unfolding real_of_float_mult[of "Float 1 1"] Float_num
722 using arctan_0_1_bounds[OF `0 \<le> real ?DIV` `real ?DIV \<le> 1`] by auto
723 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] .
727 hence "2 < real x" unfolding le_float_def Float_num by auto
728 hence "1 \<le> real x" by auto
729 hence "0 < real x" by auto
730 hence "0 < x" unfolding less_float_def by auto
732 let "?invx" = "float_divl prec 1 x"
733 have "0 \<le> arctan x" using arctan_monotone'[OF `0 \<le> real x`] using arctan_tan[of 0, unfolded tan_zero] by auto
735 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`])
736 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`)
738 have "1 / x \<noteq> 0" and "0 < 1 / x" using `0 < real x` by auto
740 have "(?lb_horner ?invx) \<le> arctan (?invx)" using arctan_0_1_bounds[OF `0 \<le> real ?invx` `real ?invx \<le> 1`] by auto
741 also have "\<dots> \<le> arctan (1 / x)" unfolding real_of_float_1[symmetric] by (rule arctan_monotone', rule float_divl)
742 finally have "arctan x \<le> pi / 2 - (?lb_horner ?invx)"
743 using `0 \<le> arctan x` arctan_inverse[OF `1 / x \<noteq> 0`]
744 unfolding real_sgn_pos[OF `0 < 1 / x`] le_diff_eq by auto
746 have "pi / 2 \<le> 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
748 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]
754 lemma arctan_boundaries:
755 "arctan x \<in> {(lb_arctan prec x) .. (ub_arctan prec x)}"
756 proof (cases "0 \<le> x")
757 case True hence "0 \<le> real x" unfolding le_float_def by auto
758 show ?thesis using ub_arctan_bound'[OF `0 \<le> real x`] lb_arctan_bound'[OF `0 \<le> real x`] unfolding atLeastAtMost_iff by auto
761 case False hence "x < 0" and "0 \<le> real ?mx" unfolding le_float_def less_float_def by auto
762 hence bounds: "lb_arctan prec ?mx \<le> arctan ?mx \<and> arctan ?mx \<le> ub_arctan prec ?mx"
763 using ub_arctan_bound'[OF `0 \<le> real ?mx`] lb_arctan_bound'[OF `0 \<le> real ?mx`] by auto
764 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`]
765 unfolding atLeastAtMost_iff using bounds[unfolded real_of_float_minus arctan_minus] by auto
768 lemma bnds_arctan: "\<forall> (x::real) lx ux. (l, u) = (lb_arctan prec lx, ub_arctan prec ux) \<and> x \<in> {lx .. ux} \<longrightarrow> l \<le> arctan x \<and> arctan x \<le> u"
769 proof (rule allI, rule allI, rule allI, rule impI)
770 fix x :: real fix lx ux
771 assume "(l, u) = (lb_arctan prec lx, ub_arctan prec ux) \<and> x \<in> {lx .. ux}"
772 hence l: "lb_arctan prec lx = l " and u: "ub_arctan prec ux = u" and x: "x \<in> {lx .. ux}" by auto
774 { from arctan_boundaries[of lx prec, unfolded l]
775 have "l \<le> arctan lx" by (auto simp del: lb_arctan.simps)
776 also have "\<dots> \<le> arctan x" using x by (auto intro: arctan_monotone')
777 finally have "l \<le> arctan x" .
779 { have "arctan x \<le> arctan ux" using x by (auto intro: arctan_monotone')
780 also have "\<dots> \<le> u" using arctan_boundaries[of ux prec, unfolded u] by (auto simp del: ub_arctan.simps)
781 finally have "arctan x \<le> u" .
782 } ultimately show "l \<le> arctan x \<and> arctan x \<le> u" ..
785 section "Sinus and Cosinus"
787 subsection "Compute the cosinus and sinus series"
789 fun ub_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
790 and lb_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
791 "ub_sin_cos_aux prec 0 i k x = 0"
792 | "ub_sin_cos_aux prec (Suc n) i k x =
793 (rapprox_rat prec 1 k) - x * (lb_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)"
794 | "lb_sin_cos_aux prec 0 i k x = 0"
795 | "lb_sin_cos_aux prec (Suc n) i k x =
796 (lapprox_rat prec 1 k) - x * (ub_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)"
798 shows "(lb_sin_cos_aux prec n 1 1 (x * x)) \<le> (\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i))) * x ^(2 * i))" (is "?lb")
799 and "(\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i))) * x^(2 * i)) \<le> (ub_sin_cos_aux prec n 1 1 (x * x))" (is "?ub")
801 have "0 \<le> real (x * x)" unfolding real_of_float_mult by auto
802 let "?f n" = "fact (2 * n)"
805 have F: "\<And>m. ((\<lambda>i. i + 2) ^^ n) m = m + 2 * n" by (induct n arbitrary: m, auto)
806 have "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) ^^ n) 1 * (((\<lambda>i. i + 2) ^^ n) 1 + 1)"
807 unfolding F by auto } note f_eq = this
809 from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0,
810 OF `0 \<le> real (x * x)` f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps]
811 show "?lb" and "?ub" by (auto simp add: power_mult power2_eq_square[of "real x"])
814 lemma cos_boundaries: assumes "0 \<le> real x" and "x \<le> pi / 2"
815 shows "cos x \<in> {(lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) .. (ub_sin_cos_aux prec (get_odd n) 1 1 (x * x))}"
816 proof (cases "real x = 0")
817 case False hence "real x \<noteq> 0" by auto
818 hence "0 < x" and "0 < real x" using `0 \<le> real x` unfolding less_float_def by auto
819 have "0 < x * x" using `0 < x` unfolding less_float_def real_of_float_mult real_of_float_0
820 using mult_pos_pos[where a="real x" and b="real x"] by auto
822 { fix x n have "(\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i))) * x ^ (2 * i))
823 = (\<Sum> i = 0 ..< 2 * n. (if even(i) then (-1 ^ (i div 2))/(real (fact i)) else 0) * x ^ i)" (is "?sum = ?ifsum")
825 have "?sum = ?sum + (\<Sum> j = 0 ..< n. 0)" by auto
827 (\<Sum> j = 0 ..< n. -1 ^ ((2 * j) div 2) / (real (fact (2 * j))) * x ^(2 * j)) + (\<Sum> j = 0 ..< n. 0)" by auto
828 also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. if even i then -1 ^ (i div 2) / (real (fact i)) * x ^ i else 0)"
829 unfolding sum_split_even_odd ..
830 also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. (if even i then -1 ^ (i div 2) / (real (fact i)) else 0) * x ^ i)"
831 by (rule setsum_cong2) auto
832 finally show ?thesis by assumption
833 qed } note morph_to_if_power = this
836 { fix n :: nat assume "0 < n"
837 hence "0 < 2 * n" by auto
838 obtain t where "0 < t" and "t < real x" and
839 cos_eq: "cos x = (\<Sum> i = 0 ..< 2 * n. (if even(i) then (-1 ^ (i div 2))/(real (fact i)) else 0) * (real x) ^ i)
840 + (cos (t + 1/2 * (2 * n) * pi) / real (fact (2*n))) * (real x)^(2*n)"
841 (is "_ = ?SUM + ?rest / ?fact * ?pow")
842 using Maclaurin_cos_expansion2[OF `0 < real x` `0 < 2 * n`] by auto
844 have "cos t * -1^n = cos t * cos (n * pi) + sin t * sin (n * pi)" by auto
845 also have "\<dots> = cos (t + n * pi)" using cos_add by auto
846 also have "\<dots> = ?rest" by auto
847 finally have "cos t * -1^n = ?rest" .
849 have "t \<le> pi / 2" using `t < real x` and `x \<le> pi / 2` by auto
850 hence "0 \<le> cos t" using `0 < t` and cos_ge_zero by auto
851 ultimately have even: "even n \<Longrightarrow> 0 \<le> ?rest" and odd: "odd n \<Longrightarrow> 0 \<le> - ?rest " by auto
853 have "0 < ?fact" by auto
854 have "0 < ?pow" using `0 < real x` by auto
858 have "(lb_sin_cos_aux prec n 1 1 (x * x)) \<le> ?SUM"
859 unfolding morph_to_if_power[symmetric] using cos_aux by auto
860 also have "\<dots> \<le> cos x"
862 from even[OF `even n`] `0 < ?fact` `0 < ?pow`
863 have "0 \<le> (?rest / ?fact) * ?pow" by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
864 thus ?thesis unfolding cos_eq by auto
866 finally have "(lb_sin_cos_aux prec n 1 1 (x * x)) \<le> cos x" .
871 have "cos x \<le> ?SUM"
873 from `0 < ?fact` and `0 < ?pow` and odd[OF `odd n`]
874 have "0 \<le> (- ?rest) / ?fact * ?pow"
875 by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
876 thus ?thesis unfolding cos_eq by auto
878 also have "\<dots> \<le> (ub_sin_cos_aux prec n 1 1 (x * x))"
879 unfolding morph_to_if_power[symmetric] using cos_aux by auto
880 finally have "cos x \<le> (ub_sin_cos_aux prec n 1 1 (x * x))" .
881 } note ub = this and lb
882 } note ub = this(1) and lb = this(2)
884 have "cos x \<le> (ub_sin_cos_aux prec (get_odd n) 1 1 (x * x))" using ub[OF odd_pos[OF get_odd] get_odd] .
885 moreover have "(lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) \<le> cos x"
886 proof (cases "0 < get_even n")
887 case True show ?thesis using lb[OF True get_even] .
890 hence "get_even n = 0" by auto
891 have "- (pi / 2) \<le> x" by (rule order_trans[OF _ `0 < real x`[THEN less_imp_le]], auto)
892 with `x \<le> pi / 2`
893 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
895 ultimately show ?thesis by auto
899 proof (cases "n = 0")
901 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
903 case False with not0_implies_Suc obtain m where "n = Suc m" by blast
904 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)
908 lemma sin_aux: assumes "0 \<le> real x"
909 shows "(x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> (\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i + 1))) * x^(2 * i + 1))" (is "?lb")
910 and "(\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i + 1))) * x^(2 * i + 1)) \<le> (x * ub_sin_cos_aux prec n 2 1 (x * x))" (is "?ub")
912 have "0 \<le> real (x * x)" unfolding real_of_float_mult by auto
913 let "?f n" = "fact (2 * n + 1)"
916 have F: "\<And>m. ((\<lambda>i. i + 2) ^^ n) m = m + 2 * n" by (induct n arbitrary: m, auto)
917 have "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) ^^ n) 2 * (((\<lambda>i. i + 2) ^^ n) 2 + 1)"
918 unfolding F by auto } note f_eq = this
920 from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0,
921 OF `0 \<le> real (x * x)` f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps]
922 show "?lb" and "?ub" using `0 \<le> real x` unfolding real_of_float_mult
923 unfolding power_add power_one_right mult_assoc[symmetric] setsum_left_distrib[symmetric]
924 unfolding mult_commute[where 'a=real]
925 by (auto intro!: mult_left_mono simp add: power_mult power2_eq_square[of "real x"])
928 lemma sin_boundaries: assumes "0 \<le> real x" and "x \<le> pi / 2"
929 shows "sin x \<in> {(x * lb_sin_cos_aux prec (get_even n) 2 1 (x * x)) .. (x * ub_sin_cos_aux prec (get_odd n) 2 1 (x * x))}"
930 proof (cases "real x = 0")
931 case False hence "real x \<noteq> 0" by auto
932 hence "0 < x" and "0 < real x" using `0 \<le> real x` unfolding less_float_def by auto
933 have "0 < x * x" using `0 < x` unfolding less_float_def real_of_float_mult real_of_float_0
934 using mult_pos_pos[where a="real x" and b="real x"] by auto
936 { 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))
937 = (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * x ^ i)" (is "?SUM = _")
939 have pow: "!!i. x ^ (2 * i + 1) = x * x ^ (2 * i)" by auto
940 have "?SUM = (\<Sum> j = 0 ..< n. 0) + ?SUM" by auto
941 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)"
942 unfolding sum_split_even_odd ..
943 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)"
944 by (rule setsum_cong2) auto
945 finally show ?thesis by assumption
946 qed } note setsum_morph = this
948 { fix n :: nat assume "0 < n"
949 hence "0 < 2 * n + 1" by auto
950 obtain t where "0 < t" and "t < real x" and
951 sin_eq: "sin 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)
952 + (sin (t + 1/2 * (2 * n + 1) * pi) / real (fact (2*n + 1))) * (real x)^(2*n + 1)"
953 (is "_ = ?SUM + ?rest / ?fact * ?pow")
954 using Maclaurin_sin_expansion3[OF `0 < 2 * n + 1` `0 < real x`] by auto
956 have "?rest = cos t * -1^n" unfolding sin_add cos_add real_of_nat_add left_distrib right_distrib by auto
958 have "t \<le> pi / 2" using `t < real x` and `x \<le> pi / 2` by auto
959 hence "0 \<le> cos t" using `0 < t` and cos_ge_zero by auto
960 ultimately have even: "even n \<Longrightarrow> 0 \<le> ?rest" and odd: "odd n \<Longrightarrow> 0 \<le> - ?rest " by auto
962 have "0 < ?fact" by (rule real_of_nat_fact_gt_zero)
963 have "0 < ?pow" using `0 < real x` by (rule zero_less_power)
967 have "(x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le>
968 (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * (real x) ^ i)"
969 using sin_aux[OF `0 \<le> real x`] unfolding setsum_morph[symmetric] by auto
970 also have "\<dots> \<le> ?SUM" by auto
971 also have "\<dots> \<le> sin x"
973 from even[OF `even n`] `0 < ?fact` `0 < ?pow`
974 have "0 \<le> (?rest / ?fact) * ?pow" by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
975 thus ?thesis unfolding sin_eq by auto
977 finally have "(x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> sin x" .
982 have "sin x \<le> ?SUM"
984 from `0 < ?fact` and `0 < ?pow` and odd[OF `odd n`]
985 have "0 \<le> (- ?rest) / ?fact * ?pow"
986 by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
987 thus ?thesis unfolding sin_eq by auto
989 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)"
991 also have "\<dots> \<le> (x * ub_sin_cos_aux prec n 2 1 (x * x))"
992 using sin_aux[OF `0 \<le> real x`] unfolding setsum_morph[symmetric] by auto
993 finally have "sin x \<le> (x * ub_sin_cos_aux prec n 2 1 (x * x))" .
994 } note ub = this and lb
995 } note ub = this(1) and lb = this(2)
997 have "sin x \<le> (x * ub_sin_cos_aux prec (get_odd n) 2 1 (x * x))" using ub[OF odd_pos[OF get_odd] get_odd] .
998 moreover have "(x * lb_sin_cos_aux prec (get_even n) 2 1 (x * x)) \<le> sin x"
999 proof (cases "0 < get_even n")
1000 case True show ?thesis using lb[OF True get_even] .
1003 hence "get_even n = 0" by auto
1004 with `x \<le> pi / 2` `0 \<le> real x`
1005 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
1007 ultimately show ?thesis by auto
1011 proof (cases "n = 0")
1013 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
1015 case False with not0_implies_Suc obtain m where "n = Suc m" by blast
1016 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)
1020 subsection "Compute the cosinus in the entire domain"
1022 definition lb_cos :: "nat \<Rightarrow> float \<Rightarrow> float" where
1023 "lb_cos prec x = (let
1024 horner = \<lambda> x. lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 1 1 (x * x) ;
1025 half = \<lambda> x. if x < 0 then - 1 else 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 definition ub_cos :: "nat \<Rightarrow> float \<Rightarrow> float" where
1031 "ub_cos prec x = (let
1032 horner = \<lambda> x. ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x) ;
1033 half = \<lambda> x. Float 1 1 * x * x - 1
1034 in if x < Float 1 -1 then horner x
1035 else if x < 1 then half (horner (x * Float 1 -1))
1036 else half (half (horner (x * Float 1 -2))))"
1038 lemma lb_cos: assumes "0 \<le> real x" and "x \<le> pi"
1039 shows "cos x \<in> {(lb_cos prec x) .. (ub_cos prec x)}" (is "?cos x \<in> {(?lb x) .. (?ub x) }")
1042 have "cos x = cos (x / 2 + x / 2)" by auto
1043 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"
1044 unfolding cos_add by auto
1045 also have "\<dots> = 2 * cos (x / 2) * cos (x / 2) - 1" by algebra
1046 finally have "cos x = 2 * cos (x / 2) * cos (x / 2) - 1" .
1047 } note x_half = this[symmetric]
1049 have "\<not> x < 0" using `0 \<le> real x` unfolding less_float_def by auto
1050 let "?ub_horner x" = "ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x)"
1051 let "?lb_horner x" = "lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 1 1 (x * x)"
1052 let "?ub_half x" = "Float 1 1 * x * x - 1"
1053 let "?lb_half x" = "if x < 0 then - 1 else Float 1 1 * x * x - 1"
1056 proof (cases "x < Float 1 -1")
1057 case True hence "x \<le> pi / 2" unfolding less_float_def using pi_ge_two by auto
1058 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
1059 using cos_boundaries[OF `0 \<le> real x` `x \<le> pi / 2`] .
1062 { fix y x :: float let ?x2 = "(x * Float 1 -1)"
1063 assume "y \<le> cos ?x2" and "-pi \<le> x" and "x \<le> pi"
1064 hence "- (pi / 2) \<le> ?x2" and "?x2 \<le> pi / 2" using pi_ge_two unfolding real_of_float_mult Float_num by auto
1065 hence "0 \<le> cos ?x2" by (rule cos_ge_zero)
1067 have "(?lb_half y) \<le> cos x"
1068 proof (cases "y < 0")
1069 case True show ?thesis using cos_ge_minus_one unfolding if_P[OF True] by auto
1072 hence "0 \<le> real y" unfolding less_float_def by auto
1073 from mult_mono[OF `y \<le> cos ?x2` `y \<le> cos ?x2` `0 \<le> cos ?x2` this]
1074 have "real y * real y \<le> cos ?x2 * cos ?x2" .
1075 hence "2 * real y * real y \<le> 2 * cos ?x2 * cos ?x2" by auto
1076 hence "2 * real y * real y - 1 \<le> 2 * cos (x / 2) * cos (x / 2) - 1" unfolding Float_num real_of_float_mult by auto
1077 thus ?thesis unfolding if_not_P[OF False] x_half Float_num real_of_float_mult real_of_float_sub by auto
1079 } note lb_half = this
1081 { fix y x :: float let ?x2 = "(x * Float 1 -1)"
1082 assume ub: "cos ?x2 \<le> y" and "- pi \<le> x" and "x \<le> pi"
1083 hence "- (pi / 2) \<le> ?x2" and "?x2 \<le> pi / 2" using pi_ge_two unfolding real_of_float_mult Float_num by auto
1084 hence "0 \<le> cos ?x2" by (rule cos_ge_zero)
1086 have "cos x \<le> (?ub_half y)"
1088 have "0 \<le> real y" using `0 \<le> cos ?x2` ub by (rule order_trans)
1089 from mult_mono[OF ub ub this `0 \<le> cos ?x2`]
1090 have "cos ?x2 * cos ?x2 \<le> real y * real y" .
1091 hence "2 * cos ?x2 * cos ?x2 \<le> 2 * real y * real y" by auto
1092 hence "2 * cos (x / 2) * cos (x / 2) - 1 \<le> 2 * real y * real y - 1" unfolding Float_num real_of_float_mult by auto
1093 thus ?thesis unfolding x_half real_of_float_mult Float_num real_of_float_sub by auto
1095 } note ub_half = this
1097 let ?x2 = "x * Float 1 -1"
1098 let ?x4 = "x * Float 1 -1 * Float 1 -1"
1100 have "-pi \<le> x" using pi_ge_zero[THEN le_imp_neg_le, unfolded minus_zero] `0 \<le> real x` by (rule order_trans)
1103 proof (cases "x < 1")
1104 case True hence "real x \<le> 1" unfolding less_float_def by auto
1105 have "0 \<le> real ?x2" and "?x2 \<le> pi / 2" using pi_ge_two `0 \<le> real x` unfolding real_of_float_mult Float_num using assms by auto
1106 from cos_boundaries[OF this]
1107 have lb: "(?lb_horner ?x2) \<le> ?cos ?x2" and ub: "?cos ?x2 \<le> (?ub_horner ?x2)" by auto
1109 have "(?lb x) \<le> ?cos x"
1111 from lb_half[OF lb `-pi \<le> x` `x \<le> pi`]
1112 show ?thesis unfolding lb_cos_def[where x=x] Let_def using `\<not> x < 0` `\<not> x < Float 1 -1` `x < 1` by auto
1114 moreover have "?cos x \<le> (?ub x)"
1116 from ub_half[OF ub `-pi \<le> x` `x \<le> pi`]
1117 show ?thesis unfolding ub_cos_def[where x=x] Let_def using `\<not> x < 0` `\<not> x < Float 1 -1` `x < 1` by auto
1119 ultimately show ?thesis by auto
1122 have "0 \<le> real ?x4" and "?x4 \<le> pi / 2" using pi_ge_two `0 \<le> real x` `x \<le> pi` unfolding real_of_float_mult Float_num by auto
1123 from cos_boundaries[OF this]
1124 have lb: "(?lb_horner ?x4) \<le> ?cos ?x4" and ub: "?cos ?x4 \<le> (?ub_horner ?x4)" by auto
1126 have eq_4: "?x2 * Float 1 -1 = x * Float 1 -2" by (cases x, auto simp add: times_float.simps)
1128 have "(?lb x) \<le> ?cos x"
1130 have "-pi \<le> ?x2" and "?x2 \<le> pi" unfolding real_of_float_mult Float_num using pi_ge_two `0 \<le> real x` `x \<le> pi` by auto
1131 from lb_half[OF lb_half[OF lb this] `-pi \<le> x` `x \<le> pi`, unfolded eq_4]
1132 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 .
1134 moreover have "?cos x \<le> (?ub x)"
1136 have "-pi \<le> ?x2" and "?x2 \<le> pi" unfolding real_of_float_mult Float_num using pi_ge_two `0 \<le> real x` ` x \<le> pi` by auto
1137 from ub_half[OF ub_half[OF ub this] `-pi \<le> x` `x \<le> pi`, unfolded eq_4]
1138 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 .
1140 ultimately show ?thesis by auto
1145 lemma lb_cos_minus: assumes "-pi \<le> x" and "real x \<le> 0"
1146 shows "cos (real(-x)) \<in> {(lb_cos prec (-x)) .. (ub_cos prec (-x))}"
1148 have "0 \<le> real (-x)" and "(-x) \<le> pi" using `-pi \<le> x` `real x \<le> 0` by auto
1149 from lb_cos[OF this] show ?thesis .
1152 definition bnds_cos :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where
1153 "bnds_cos prec lx ux = (let
1154 lpi = round_down prec (lb_pi prec) ;
1155 upi = round_up prec (ub_pi prec) ;
1156 k = floor_fl (float_divr prec (lx + lpi) (2 * lpi)) ;
1157 lx = lx - k * 2 * (if k < 0 then lpi else upi) ;
1158 ux = ux - k * 2 * (if k < 0 then upi else lpi)
1159 in if - lpi \<le> lx \<and> ux \<le> 0 then (lb_cos prec (-lx), ub_cos prec (-ux))
1160 else if 0 \<le> lx \<and> ux \<le> lpi then (lb_cos prec ux, ub_cos prec lx)
1161 else if - lpi \<le> lx \<and> ux \<le> lpi then (min (lb_cos prec (-lx)) (lb_cos prec ux), Float 1 0)
1162 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))))
1163 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)))
1164 else (Float -1 0, Float 1 0))"
1167 obtains k :: int where "real k = (floor_fl f)"
1169 assume *: "\<And> k :: int. real k = (floor_fl f) \<Longrightarrow> thesis"
1170 obtain m e where fl: "Float m e = floor_fl f" by (cases "floor_fl f", auto)
1171 from floor_pos_exp[OF this]
1172 have "real (m* 2^(nat e)) = (floor_fl f)"
1173 by (auto simp add: fl[symmetric] real_of_float_def pow2_def)
1174 from *[OF this] show thesis by blast
1177 lemma float_remove_real_numeral[simp]: "(number_of k :: float) = (number_of k :: real)"
1179 have "(number_of k :: float) = real k"
1180 unfolding number_of_float_def real_of_float_def pow2_def by auto
1181 also have "\<dots> = (number_of k :: int)"
1182 by (simp add: number_of_is_id)
1183 finally show ?thesis by auto
1186 lemma cos_periodic_nat[simp]: fixes n :: nat shows "cos (x + n * (2 * pi)) = cos x"
1187 proof (induct n arbitrary: x)
1189 have split_pi_off: "x + (Suc n) * (2 * pi) = (x + n * (2 * pi)) + 2 * pi"
1190 unfolding Suc_eq_plus1 real_of_nat_add real_of_one left_distrib by auto
1191 show ?case unfolding split_pi_off using Suc by auto
1194 lemma cos_periodic_int[simp]: fixes i :: int shows "cos (x + i * (2 * pi)) = cos x"
1195 proof (cases "0 \<le> i")
1196 case True hence i_nat: "real i = nat i" by auto
1197 show ?thesis unfolding i_nat by auto
1199 case False hence i_nat: "i = - real (nat (-i))" by auto
1200 have "cos x = cos (x + i * (2 * pi) - i * (2 * pi))" by auto
1201 also have "\<dots> = cos (x + i * (2 * pi))"
1202 unfolding i_nat mult_minus_left diff_minus_eq_add by (rule cos_periodic_nat)
1203 finally show ?thesis by auto
1206 lemma bnds_cos: "\<forall> (x::real) lx ux. (l, u) = bnds_cos prec lx ux \<and> x \<in> {lx .. ux} \<longrightarrow> l \<le> cos x \<and> cos x \<le> u"
1207 proof ((rule allI | rule impI | erule conjE) +)
1208 fix x :: real fix lx ux
1209 assume bnds: "(l, u) = bnds_cos prec lx ux" and x: "x \<in> {lx .. ux}"
1211 let ?lpi = "round_down prec (lb_pi prec)"
1212 let ?upi = "round_up prec (ub_pi prec)"
1213 let ?k = "floor_fl (float_divr prec (lx + ?lpi) (2 * ?lpi))"
1214 let ?lx = "lx - ?k * 2 * (if ?k < 0 then ?lpi else ?upi)"
1215 let ?ux = "ux - ?k * 2 * (if ?k < 0 then ?upi else ?lpi)"
1217 obtain k :: int where k: "k = real ?k" using floor_int .
1219 have upi: "pi \<le> ?upi" and lpi: "?lpi \<le> pi"
1220 using round_up[of "ub_pi prec" prec] pi_boundaries[of prec]
1221 round_down[of prec "lb_pi prec"] by auto
1222 hence "?lx \<le> x - k * (2 * pi) \<and> x - k * (2 * pi) \<le> ?ux"
1223 using x by (cases "k = 0") (auto intro!: add_mono
1224 simp add: diff_minus k[symmetric] less_float_def)
1225 note lx = this[THEN conjunct1] and ux = this[THEN conjunct2]
1226 hence lx_less_ux: "?lx \<le> real ?ux" by (rule order_trans)
1228 { assume "- ?lpi \<le> ?lx" and x_le_0: "x - k * (2 * pi) \<le> 0"
1229 with lpi[THEN le_imp_neg_le] lx
1230 have pi_lx: "- pi \<le> ?lx" and lx_0: "real ?lx \<le> 0"
1231 by (simp_all add: le_float_def)
1233 have "(lb_cos prec (- ?lx)) \<le> cos (real (- ?lx))"
1234 using lb_cos_minus[OF pi_lx lx_0] by simp
1235 also have "\<dots> \<le> cos (x + (-k) * (2 * pi))"
1236 using cos_monotone_minus_pi_0'[OF pi_lx lx x_le_0]
1237 by (simp only: real_of_float_minus real_of_int_minus
1238 cos_minus diff_minus mult_minus_left)
1239 finally have "(lb_cos prec (- ?lx)) \<le> cos x"
1240 unfolding cos_periodic_int . }
1241 note negative_lx = this
1243 { assume "0 \<le> ?lx" and pi_x: "x - k * (2 * pi) \<le> pi"
1245 have pi_lx: "?lx \<le> pi" and lx_0: "0 \<le> real ?lx"
1246 by (auto simp add: le_float_def)
1248 have "cos (x + (-k) * (2 * pi)) \<le> cos ?lx"
1249 using cos_monotone_0_pi'[OF lx_0 lx pi_x]
1250 by (simp only: real_of_float_minus real_of_int_minus
1251 cos_minus diff_minus mult_minus_left)
1252 also have "\<dots> \<le> (ub_cos prec ?lx)"
1253 using lb_cos[OF lx_0 pi_lx] by simp
1254 finally have "cos x \<le> (ub_cos prec ?lx)"
1255 unfolding cos_periodic_int . }
1256 note positive_lx = this
1258 { assume pi_x: "- pi \<le> x - k * (2 * pi)" and "?ux \<le> 0"
1260 have pi_ux: "- pi \<le> ?ux" and ux_0: "real ?ux \<le> 0"
1261 by (simp_all add: le_float_def)
1263 have "cos (x + (-k) * (2 * pi)) \<le> cos (real (- ?ux))"
1264 using cos_monotone_minus_pi_0'[OF pi_x ux ux_0]
1265 by (simp only: real_of_float_minus real_of_int_minus
1266 cos_minus diff_minus mult_minus_left)
1267 also have "\<dots> \<le> (ub_cos prec (- ?ux))"
1268 using lb_cos_minus[OF pi_ux ux_0, of prec] by simp
1269 finally have "cos x \<le> (ub_cos prec (- ?ux))"
1270 unfolding cos_periodic_int . }
1271 note negative_ux = this
1273 { assume "?ux \<le> ?lpi" and x_ge_0: "0 \<le> x - k * (2 * pi)"
1275 have pi_ux: "?ux \<le> pi" and ux_0: "0 \<le> real ?ux"
1276 by (simp_all add: le_float_def)
1278 have "(lb_cos prec ?ux) \<le> cos ?ux"
1279 using lb_cos[OF ux_0 pi_ux] by simp
1280 also have "\<dots> \<le> cos (x + (-k) * (2 * pi))"
1281 using cos_monotone_0_pi'[OF x_ge_0 ux pi_ux]
1282 by (simp only: real_of_float_minus real_of_int_minus
1283 cos_minus diff_minus mult_minus_left)
1284 finally have "(lb_cos prec ?ux) \<le> cos x"
1285 unfolding cos_periodic_int . }
1286 note positive_ux = this
1288 show "l \<le> cos x \<and> cos x \<le> u"
1289 proof (cases "- ?lpi \<le> ?lx \<and> ?ux \<le> 0")
1291 have l: "l = lb_cos prec (-?lx)"
1292 and u: "u = ub_cos prec (-?ux)"
1293 by (auto simp add: bnds_cos_def Let_def)
1295 from True lpi[THEN le_imp_neg_le] lx ux
1296 have "- pi \<le> x - k * (2 * pi)"
1297 and "x - k * (2 * pi) \<le> 0"
1298 by (auto simp add: le_float_def)
1299 with True negative_ux negative_lx
1300 show ?thesis unfolding l u by simp
1301 next case False note 1 = this show ?thesis
1302 proof (cases "0 \<le> ?lx \<and> ?ux \<le> ?lpi")
1303 case True with bnds 1
1304 have l: "l = lb_cos prec ?ux"
1305 and u: "u = ub_cos prec ?lx"
1306 by (auto simp add: bnds_cos_def Let_def)
1309 have "0 \<le> x - k * (2 * pi)"
1310 and "x - k * (2 * pi) \<le> pi"
1311 by (auto simp add: le_float_def)
1312 with True positive_ux positive_lx
1313 show ?thesis unfolding l u by simp
1314 next case False note 2 = this show ?thesis
1315 proof (cases "- ?lpi \<le> ?lx \<and> ?ux \<le> ?lpi")
1316 case True note Cond = this with bnds 1 2
1317 have l: "l = min (lb_cos prec (-?lx)) (lb_cos prec ?ux)"
1318 and u: "u = Float 1 0"
1319 by (auto simp add: bnds_cos_def Let_def)
1321 show ?thesis unfolding u l using negative_lx positive_ux Cond
1322 by (cases "x - k * (2 * pi) < 0", simp_all add: real_of_float_min)
1323 next case False note 3 = this show ?thesis
1324 proof (cases "0 \<le> ?lx \<and> ?ux \<le> 2 * ?lpi")
1325 case True note Cond = this with bnds 1 2 3
1326 have l: "l = Float -1 0"
1327 and u: "u = max (ub_cos prec ?lx) (ub_cos prec (- (?ux - 2 * ?lpi)))"
1328 by (auto simp add: bnds_cos_def Let_def)
1330 have "cos x \<le> real u"
1331 proof (cases "x - k * (2 * pi) < pi")
1332 case True hence "x - k * (2 * pi) \<le> pi" by simp
1333 from positive_lx[OF Cond[THEN conjunct1] this]
1334 show ?thesis unfolding u by (simp add: real_of_float_max)
1336 case False hence "pi \<le> x - k * (2 * pi)" by simp
1337 hence pi_x: "- pi \<le> x - k * (2 * pi) - 2 * pi" by simp
1339 have "?ux \<le> 2 * pi" using Cond lpi by (auto simp add: le_float_def)
1340 hence "x - k * (2 * pi) - 2 * pi \<le> 0" using ux by simp
1342 have ux_0: "real (?ux - 2 * ?lpi) \<le> 0"
1343 using Cond by (auto simp add: le_float_def)
1345 from 2 and Cond have "\<not> ?ux \<le> ?lpi" by auto
1346 hence "- ?lpi \<le> ?ux - 2 * ?lpi" by (auto simp add: le_float_def)
1347 hence pi_ux: "- pi \<le> (?ux - 2 * ?lpi)"
1348 using lpi[THEN le_imp_neg_le] by (auto simp add: le_float_def)
1350 have x_le_ux: "x - k * (2 * pi) - 2 * pi \<le> (?ux - 2 * ?lpi)"
1351 using ux lpi by auto
1353 have "cos x = cos (x + (-k) * (2 * pi) + (-1::int) * (2 * pi))"
1354 unfolding cos_periodic_int ..
1355 also have "\<dots> \<le> cos ((?ux - 2 * ?lpi))"
1356 using cos_monotone_minus_pi_0'[OF pi_x x_le_ux ux_0]
1357 by (simp only: real_of_float_minus real_of_int_minus real_of_one
1358 number_of_Min diff_minus mult_minus_left mult_1_left)
1359 also have "\<dots> = cos ((- (?ux - 2 * ?lpi)))"
1360 unfolding real_of_float_minus cos_minus ..
1361 also have "\<dots> \<le> (ub_cos prec (- (?ux - 2 * ?lpi)))"
1362 using lb_cos_minus[OF pi_ux ux_0] by simp
1363 finally show ?thesis unfolding u by (simp add: real_of_float_max)
1365 thus ?thesis unfolding l by auto
1366 next case False note 4 = this show ?thesis
1367 proof (cases "-2 * ?lpi \<le> ?lx \<and> ?ux \<le> 0")
1368 case True note Cond = this with bnds 1 2 3 4
1369 have l: "l = Float -1 0"
1370 and u: "u = max (ub_cos prec (?lx + 2 * ?lpi)) (ub_cos prec (-?ux))"
1371 by (auto simp add: bnds_cos_def Let_def)
1373 have "cos x \<le> u"
1374 proof (cases "-pi < x - k * (2 * pi)")
1375 case True hence "-pi \<le> x - k * (2 * pi)" by simp
1376 from negative_ux[OF this Cond[THEN conjunct2]]
1377 show ?thesis unfolding u by (simp add: real_of_float_max)
1379 case False hence "x - k * (2 * pi) \<le> -pi" by simp
1380 hence pi_x: "x - k * (2 * pi) + 2 * pi \<le> pi" by simp
1382 have "-2 * pi \<le> ?lx" using Cond lpi by (auto simp add: le_float_def)
1384 hence "0 \<le> x - k * (2 * pi) + 2 * pi" using lx by simp
1386 have lx_0: "0 \<le> real (?lx + 2 * ?lpi)"
1387 using Cond lpi by (auto simp add: le_float_def)
1389 from 1 and Cond have "\<not> -?lpi \<le> ?lx" by auto
1390 hence "?lx + 2 * ?lpi \<le> ?lpi" by (auto simp add: le_float_def)
1391 hence pi_lx: "(?lx + 2 * ?lpi) \<le> pi"
1392 using lpi[THEN le_imp_neg_le] by (auto simp add: le_float_def)
1394 have lx_le_x: "(?lx + 2 * ?lpi) \<le> x - k * (2 * pi) + 2 * pi"
1395 using lx lpi by auto
1397 have "cos x = cos (x + (-k) * (2 * pi) + (1 :: int) * (2 * pi))"
1398 unfolding cos_periodic_int ..
1399 also have "\<dots> \<le> cos ((?lx + 2 * ?lpi))"
1400 using cos_monotone_0_pi'[OF lx_0 lx_le_x pi_x]
1401 by (simp only: real_of_float_minus real_of_int_minus real_of_one
1402 number_of_Min diff_minus mult_minus_left mult_1_left)
1403 also have "\<dots> \<le> (ub_cos prec (?lx + 2 * ?lpi))"
1404 using lb_cos[OF lx_0 pi_lx] by simp
1405 finally show ?thesis unfolding u by (simp add: real_of_float_max)
1407 thus ?thesis unfolding l by auto
1409 case False with bnds 1 2 3 4 show ?thesis by (auto simp add: bnds_cos_def Let_def)
1413 section "Exponential function"
1415 subsection "Compute the series of the exponential function"
1417 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
1418 "ub_exp_horner prec 0 i k x = 0" |
1419 "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" |
1420 "lb_exp_horner prec 0 i k x = 0" |
1421 "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"
1423 lemma bnds_exp_horner: assumes "real x \<le> 0"
1424 shows "exp x \<in> { lb_exp_horner prec (get_even n) 1 1 x .. ub_exp_horner prec (get_odd n) 1 1 x }"
1427 have F: "\<And> m. ((\<lambda>i. i + 1) ^^ n) m = n + m" by (induct n, auto)
1428 have "fact (Suc n) = fact n * ((\<lambda>i. i + 1) ^^ n) 1" unfolding F by auto } note f_eq = this
1430 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,
1431 OF assms f_eq lb_exp_horner.simps ub_exp_horner.simps]
1433 { have "lb_exp_horner prec (get_even n) 1 1 x \<le> (\<Sum>j = 0..<get_even n. 1 / real (fact j) * real x ^ j)"
1434 using bounds(1) by auto
1435 also have "\<dots> \<le> exp x"
1437 obtain t where "\<bar>t\<bar> \<le> \<bar>real x\<bar>" and "exp 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)"
1438 using Maclaurin_exp_le by blast
1439 moreover have "0 \<le> exp t / real (fact (get_even n)) * (real x) ^ (get_even n)"
1440 by (auto intro!: mult_nonneg_nonneg divide_nonneg_pos simp add: get_even zero_le_even_power exp_gt_zero)
1441 ultimately show ?thesis
1442 using get_odd exp_gt_zero by (auto intro!: mult_nonneg_nonneg)
1444 finally have "lb_exp_horner prec (get_even n) 1 1 x \<le> exp x" .
1447 have x_less_zero: "real x ^ get_odd n \<le> 0"
1448 proof (cases "real x = 0")
1450 have "(get_odd n) \<noteq> 0" using get_odd[THEN odd_pos] by auto
1451 thus ?thesis unfolding True power_0_left by auto
1453 case False hence "real x < 0" using `real x \<le> 0` by auto
1454 show ?thesis by (rule less_imp_le, auto simp add: power_less_zero_eq get_odd `real x < 0`)
1457 obtain t where "\<bar>t\<bar> \<le> \<bar>real x\<bar>" and "exp 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)"
1458 using Maclaurin_exp_le by blast
1459 moreover have "exp t / real (fact (get_odd n)) * (real x) ^ (get_odd n) \<le> 0"
1460 by (auto intro!: mult_nonneg_nonpos divide_nonpos_pos simp add: x_less_zero exp_gt_zero)
1461 ultimately have "exp x \<le> (\<Sum>j = 0..<get_odd n. 1 / real (fact j) * real x ^ j)"
1462 using get_odd exp_gt_zero by (auto intro!: mult_nonneg_nonneg)
1463 also have "\<dots> \<le> ub_exp_horner prec (get_odd n) 1 1 x"
1464 using bounds(2) by auto
1465 finally have "exp x \<le> ub_exp_horner prec (get_odd n) 1 1 x" .
1466 } ultimately show ?thesis by auto
1469 subsection "Compute the exponential function on the entire domain"
1471 function ub_exp :: "nat \<Rightarrow> float \<Rightarrow> float" and lb_exp :: "nat \<Rightarrow> float \<Rightarrow> float" where
1472 "lb_exp prec x = (if 0 < x then float_divl prec 1 (ub_exp prec (-x))
1474 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)
1475 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))
1477 "ub_exp prec x = (if 0 < x then float_divr prec 1 (lb_exp prec (-x))
1478 else if x < - 1 then (case floor_fl x of (Float m e) \<Rightarrow>
1479 (ub_exp_horner prec (get_odd (prec + 2)) 1 1 (float_divr prec x (- Float m e))) ^ (nat (-m) * 2 ^ nat e))
1480 else ub_exp_horner prec (get_odd (prec + 2)) 1 1 x)"
1481 by pat_completeness auto
1482 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)
1484 lemma exp_m1_ge_quarter: "(1 / 4 :: real) \<le> exp (- 1)"
1486 have eq4: "4 = Suc (Suc (Suc (Suc 0)))" by auto
1488 have "1 / 4 = (Float 1 -2)" unfolding Float_num by auto
1489 also have "\<dots> \<le> lb_exp_horner 1 (get_even 4) 1 1 (- 1)"
1490 unfolding get_even_def eq4
1491 by (auto simp add: lapprox_posrat_def rapprox_posrat_def normfloat.simps)
1492 also have "\<dots> \<le> exp (- 1 :: float)" using bnds_exp_horner[where x="- 1"] by auto
1493 finally show ?thesis unfolding real_of_float_minus real_of_float_1 .
1496 lemma lb_exp_pos: assumes "\<not> 0 < x" shows "0 < lb_exp prec x"
1498 let "?lb_horner x" = "lb_exp_horner prec (get_even (prec + 2)) 1 1 x"
1499 let "?horner x" = "let y = ?lb_horner x in if y \<le> 0 then Float 1 -2 else y"
1500 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)
1501 moreover { fix x :: float fix num :: nat
1502 have "0 < real (?horner x) ^ num" using `0 < ?horner x`[unfolded less_float_def real_of_float_0] by (rule zero_less_power)
1503 also have "\<dots> = (?horner x) ^ num" using float_power by auto
1504 finally have "0 < real ((?horner x) ^ num)" .
1506 ultimately show ?thesis
1507 unfolding lb_exp.simps if_not_P[OF `\<not> 0 < x`] Let_def
1508 by (cases "floor_fl x", cases "x < - 1", auto simp add: float_power le_float_def less_float_def)
1511 lemma exp_boundaries': assumes "x \<le> 0"
1512 shows "exp x \<in> { (lb_exp prec x) .. (ub_exp prec x)}"
1514 let "?lb_exp_horner x" = "lb_exp_horner prec (get_even (prec + 2)) 1 1 x"
1515 let "?ub_exp_horner x" = "ub_exp_horner prec (get_odd (prec + 2)) 1 1 x"
1517 have "real x \<le> 0" and "\<not> x > 0" using `x \<le> 0` unfolding le_float_def less_float_def by auto
1519 proof (cases "x < - 1")
1520 case False hence "- 1 \<le> real x" unfolding less_float_def by auto
1522 proof (cases "?lb_exp_horner x \<le> 0")
1523 from `\<not> x < - 1` have "- 1 \<le> real x" unfolding less_float_def by auto
1524 hence "exp (- 1) \<le> exp x" unfolding exp_le_cancel_iff .
1525 from order_trans[OF exp_m1_ge_quarter this]
1526 have "Float 1 -2 \<le> exp x" unfolding Float_num .
1528 ultimately show ?thesis using bnds_exp_horner `real x \<le> 0` `\<not> x > 0` `\<not> x < - 1` by auto
1530 case False thus ?thesis using bnds_exp_horner `real x \<le> 0` `\<not> x > 0` `\<not> x < - 1` by (auto simp add: Let_def)
1535 obtain m e where Float_floor: "floor_fl x = Float m e" by (cases "floor_fl x", auto)
1536 let ?num = "nat (- m) * 2 ^ nat e"
1538 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)
1539 hence "real (floor_fl x) < 0" unfolding Float_floor real_of_float_simp using zero_less_pow2[of xe] by auto
1541 unfolding less_float_def real_of_float_0 Float_floor real_of_float_simp
1542 unfolding pos_prod_lt[OF zero_less_pow2[of e], unfolded mult_commute] by auto
1543 hence "1 \<le> - m" by auto
1544 hence "0 < nat (- m)" by auto
1546 have "0 \<le> e" using floor_pos_exp Float_floor[symmetric] by auto
1547 hence "(0::nat) < 2 ^ nat e" by auto
1548 ultimately have "0 < ?num" by auto
1549 hence "real ?num \<noteq> 0" by auto
1550 have e_nat: "(nat e) = e" using `0 \<le> e` by auto
1551 have num_eq: "real ?num = - floor_fl x" using `0 < nat (- m)`
1552 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
1553 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 .
1554 hence "real (floor_fl x) < 0" unfolding less_float_def by auto
1556 have "exp x \<le> ub_exp prec x"
1558 have div_less_zero: "real (float_divr prec x (- floor_fl x)) \<le> 0"
1559 using float_divr_nonpos_pos_upper_bound[OF `x \<le> 0` `0 < - floor_fl x`] unfolding le_float_def real_of_float_0 .
1561 have "exp x = exp (?num * (x / ?num))" using `real ?num \<noteq> 0` by auto
1562 also have "\<dots> = exp (x / ?num) ^ ?num" unfolding exp_real_of_nat_mult ..
1563 also have "\<dots> \<le> exp (float_divr prec x (- floor_fl x)) ^ ?num" unfolding num_eq
1564 by (rule power_mono, rule exp_le_cancel_iff[THEN iffD2], rule float_divr) auto
1565 also have "\<dots> \<le> (?ub_exp_horner (float_divr prec x (- floor_fl x))) ^ ?num" unfolding float_power
1566 by (rule power_mono, rule bnds_exp_horner[OF div_less_zero, unfolded atLeastAtMost_iff, THEN conjunct2], auto)
1567 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 .
1570 have "lb_exp prec x \<le> exp x"
1572 let ?divl = "float_divl prec x (- Float m e)"
1573 let ?horner = "?lb_exp_horner ?divl"
1576 proof (cases "?horner \<le> 0")
1577 case False hence "0 \<le> real ?horner" unfolding le_float_def by auto
1579 have div_less_zero: "real (float_divl prec x (- floor_fl x)) \<le> 0"
1580 using `real (floor_fl x) < 0` `real x \<le> 0` by (auto intro!: order_trans[OF float_divl] divide_nonpos_neg)
1582 have "(?lb_exp_horner (float_divl prec x (- floor_fl x))) ^ ?num \<le>
1583 exp (float_divl prec x (- floor_fl x)) ^ ?num" unfolding float_power
1584 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)
1585 also have "\<dots> \<le> exp (x / ?num) ^ ?num" unfolding num_eq
1586 using float_divl by (auto intro!: power_mono simp del: real_of_float_minus)
1587 also have "\<dots> = exp (?num * (x / ?num))" unfolding exp_real_of_nat_mult ..
1588 also have "\<dots> = exp x" using `real ?num \<noteq> 0` by auto
1589 finally show ?thesis
1590 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
1593 have "real (floor_fl x) \<noteq> 0" and "real (floor_fl x) \<le> 0" using `real (floor_fl x) < 0` by auto
1594 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`]]
1595 have "- 1 \<le> x / (- floor_fl x)" unfolding real_of_float_minus by auto
1596 from order_trans[OF exp_m1_ge_quarter this[unfolded exp_le_cancel_iff[where x="- 1", symmetric]]]
1597 have "Float 1 -2 \<le> exp (x / (- floor_fl x))" unfolding Float_num .
1598 hence "real (Float 1 -2) ^ ?num \<le> exp (x / (- floor_fl x)) ^ ?num"
1599 by (auto intro!: power_mono simp add: Float_num)
1600 also have "\<dots> = exp x" unfolding num_eq exp_real_of_nat_mult[symmetric] using `real (floor_fl x) \<noteq> 0` by auto
1601 finally show ?thesis
1602 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 .
1605 ultimately show ?thesis by auto
1609 lemma exp_boundaries: "exp x \<in> { lb_exp prec x .. ub_exp prec x }"
1612 proof (cases "0 < x")
1613 case False hence "x \<le> 0" unfolding less_float_def le_float_def by auto
1614 from exp_boundaries'[OF this] show ?thesis .
1616 case True hence "-x \<le> 0" unfolding less_float_def le_float_def by auto
1618 have "lb_exp prec x \<le> exp x"
1620 from exp_boundaries'[OF `-x \<le> 0`]
1621 have ub_exp: "exp (- real x) \<le> ub_exp prec (-x)" unfolding atLeastAtMost_iff real_of_float_minus by auto
1623 have "float_divl prec 1 (ub_exp prec (-x)) \<le> 1 / ub_exp prec (-x)" using float_divl[where x=1] by auto
1624 also have "\<dots> \<le> exp x"
1625 using ub_exp[unfolded inverse_le_iff_le[OF order_less_le_trans[OF exp_gt_zero ub_exp] exp_gt_zero, symmetric]]
1626 unfolding exp_minus nonzero_inverse_inverse_eq[OF exp_not_eq_zero] inverse_eq_divide by auto
1627 finally show ?thesis unfolding lb_exp.simps if_P[OF True] .
1630 have "exp x \<le> ub_exp prec x"
1632 have "\<not> 0 < -x" using `0 < x` unfolding less_float_def by auto
1634 from exp_boundaries'[OF `-x \<le> 0`]
1635 have lb_exp: "lb_exp prec (-x) \<le> exp (- real x)" unfolding atLeastAtMost_iff real_of_float_minus by auto
1637 have "exp x \<le> (1 :: float) / lb_exp prec (-x)"
1638 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],
1640 unfolding exp_minus nonzero_inverse_inverse_eq[OF exp_not_eq_zero] inverse_eq_divide real_of_float_1 by auto
1641 also have "\<dots> \<le> float_divr prec 1 (lb_exp prec (-x))" using float_divr .
1642 finally show ?thesis unfolding ub_exp.simps if_P[OF True] .
1644 ultimately show ?thesis by auto
1648 lemma bnds_exp: "\<forall> (x::real) lx ux. (l, u) = (lb_exp prec lx, ub_exp prec ux) \<and> x \<in> {lx .. ux} \<longrightarrow> l \<le> exp x \<and> exp x \<le> u"
1649 proof (rule allI, rule allI, rule allI, rule impI)
1650 fix x::real and lx ux
1651 assume "(l, u) = (lb_exp prec lx, ub_exp prec ux) \<and> x \<in> {lx .. ux}"
1652 hence l: "lb_exp prec lx = l " and u: "ub_exp prec ux = u" and x: "x \<in> {lx .. ux}" by auto
1654 { from exp_boundaries[of lx prec, unfolded l]
1655 have "l \<le> exp lx" by (auto simp del: lb_exp.simps)
1656 also have "\<dots> \<le> exp x" using x by auto
1657 finally have "l \<le> exp x" .
1659 { have "exp x \<le> exp ux" using x by auto
1660 also have "\<dots> \<le> u" using exp_boundaries[of ux prec, unfolded u] by (auto simp del: ub_exp.simps)
1661 finally have "exp x \<le> u" .
1662 } ultimately show "l \<le> exp x \<and> exp x \<le> u" ..
1667 subsection "Compute the logarithm series"
1669 fun ub_ln_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
1670 and lb_ln_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
1671 "ub_ln_horner prec 0 i x = 0" |
1672 "ub_ln_horner prec (Suc n) i x = rapprox_rat prec 1 (int i) - x * lb_ln_horner prec n (Suc i) x" |
1673 "lb_ln_horner prec 0 i x = 0" |
1674 "lb_ln_horner prec (Suc n) i x = lapprox_rat prec 1 (int i) - x * ub_ln_horner prec n (Suc i) x"
1677 assumes "0 \<le> x" and "x < 1"
1678 shows "(\<Sum>i=0..<2*n. -1^i * (1 / real (i + 1)) * x ^ (Suc i)) \<le> ln (x + 1)" (is "?lb")
1679 and "ln (x + 1) \<le> (\<Sum>i=0..<2*n + 1. -1^i * (1 / real (i + 1)) * x ^ (Suc i))" (is "?ub")
1681 let "?a n" = "(1/real (n +1)) * x ^ (Suc n)"
1683 have ln_eq: "(\<Sum> i. -1^i * ?a i) = ln (x + 1)"
1684 using ln_series[of "x + 1"] `0 \<le> x` `x < 1` by auto
1686 have "norm x < 1" using assms by auto
1687 have "?a ----> 0" unfolding Suc_eq_plus1[symmetric] inverse_eq_divide[symmetric]
1688 using LIMSEQ_mult[OF LIMSEQ_inverse_real_of_nat LIMSEQ_Suc[OF LIMSEQ_power_zero[OF `norm x < 1`]]] by auto
1689 { fix n have "0 \<le> ?a n" by (rule mult_nonneg_nonneg, auto intro!: mult_nonneg_nonneg simp add: `0 \<le> x`) }
1690 { fix n have "?a (Suc n) \<le> ?a n" unfolding inverse_eq_divide[symmetric]
1691 proof (rule mult_mono)
1692 show "0 \<le> x ^ Suc (Suc n)" by (auto intro!: mult_nonneg_nonneg simp add: `0 \<le> x`)
1693 have "x ^ Suc (Suc n) \<le> x ^ Suc n * 1" unfolding power_Suc2 mult_assoc[symmetric]
1694 by (rule mult_left_mono, fact less_imp_le[OF `x < 1`], auto intro!: mult_nonneg_nonneg simp add: `0 \<le> x`)
1695 thus "x ^ Suc (Suc n) \<le> x ^ Suc n" by auto
1697 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]
1698 show "?lb" and "?ub" by auto
1701 lemma ln_float_bounds:
1702 assumes "0 \<le> real x" and "real x < 1"
1703 shows "x * lb_ln_horner prec (get_even n) 1 x \<le> ln (x + 1)" (is "?lb \<le> ?ln")
1704 and "ln (x + 1) \<le> x * ub_ln_horner prec (get_odd n) 1 x" (is "?ln \<le> ?ub")
1706 obtain ev where ev: "get_even n = 2 * ev" using get_even_double ..
1707 obtain od where od: "get_odd n = 2 * od + 1" using get_odd_double ..
1709 let "?s n" = "-1^n * (1 / real (1 + n)) * (real x)^(Suc n)"
1711 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
1712 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",
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 also have "\<dots> \<le> ?ln" using ln_bounds(1)[OF `0 \<le> real x` `real x < 1`] by auto
1716 finally show "?lb \<le> ?ln" .
1718 have "?ln \<le> setsum ?s {0 ..< 2 * od + 1}" using ln_bounds(2)[OF `0 \<le> real x` `real x < 1`] by auto
1719 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
1720 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",
1721 OF `0 \<le> real x` refl lb_ln_horner.simps ub_ln_horner.simps] `0 \<le> real x`
1722 by (rule mult_right_mono)
1723 finally show "?ln \<le> ?ub" .
1726 lemma ln_add: assumes "0 < x" and "0 < y" shows "ln (x + y) = ln x + ln (1 + y / x)"
1728 have "x \<noteq> 0" using assms by auto
1729 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
1731 have "0 < y / x" using assms divide_pos_pos by auto
1732 hence "0 < 1 + y / x" by auto
1733 ultimately show ?thesis using ln_mult assms by auto
1736 subsection "Compute the logarithm of 2"
1738 definition ub_ln2 where "ub_ln2 prec = (let third = rapprox_rat (max prec 1) 1 3
1739 in (Float 1 -1 * ub_ln_horner prec (get_odd prec) 1 (Float 1 -1)) +
1740 (third * ub_ln_horner prec (get_odd prec) 1 third))"
1741 definition lb_ln2 where "lb_ln2 prec = (let third = lapprox_rat prec 1 3
1742 in (Float 1 -1 * lb_ln_horner prec (get_even prec) 1 (Float 1 -1)) +
1743 (third * lb_ln_horner prec (get_even prec) 1 third))"
1745 lemma ub_ln2: "ln 2 \<le> ub_ln2 prec" (is "?ub_ln2")
1746 and lb_ln2: "lb_ln2 prec \<le> ln 2" (is "?lb_ln2")
1748 let ?uthird = "rapprox_rat (max prec 1) 1 3"
1749 let ?lthird = "lapprox_rat prec 1 3"
1751 have ln2_sum: "ln 2 = ln (1/2 + 1) + ln (1 / 3 + 1)"
1752 using ln_add[of "3 / 2" "1 / 2"] by auto
1753 have lb3: "?lthird \<le> 1 / 3" using lapprox_rat[of prec 1 3] by auto
1754 hence lb3_ub: "real ?lthird < 1" by auto
1755 have lb3_lb: "0 \<le> real ?lthird" using lapprox_rat_bottom[of 1 3] by auto
1756 have ub3: "1 / 3 \<le> ?uthird" using rapprox_rat[of 1 3] by auto
1757 hence ub3_lb: "0 \<le> real ?uthird" by auto
1759 have lb2: "0 \<le> real (Float 1 -1)" and ub2: "real (Float 1 -1) < 1" unfolding Float_num by auto
1761 have "0 \<le> (1::int)" and "0 < (3::int)" by auto
1762 have ub3_ub: "real ?uthird < 1" unfolding rapprox_rat.simps(2)[OF `0 \<le> 1` `0 < 3`]
1763 by (rule rapprox_posrat_less1, auto)
1765 have third_gt0: "(0 :: real) < 1 / 3 + 1" by auto
1766 have uthird_gt0: "0 < real ?uthird + 1" using ub3_lb by auto
1767 have lthird_gt0: "0 < real ?lthird + 1" using lb3_lb by auto
1769 show ?ub_ln2 unfolding ub_ln2_def Let_def real_of_float_add ln2_sum Float_num(4)[symmetric]
1770 proof (rule add_mono, fact ln_float_bounds(2)[OF lb2 ub2])
1771 have "ln (1 / 3 + 1) \<le> ln (real ?uthird + 1)" unfolding ln_le_cancel_iff[OF third_gt0 uthird_gt0] using ub3 by auto
1772 also have "\<dots> \<le> ?uthird * ub_ln_horner prec (get_odd prec) 1 ?uthird"
1773 using ln_float_bounds(2)[OF ub3_lb ub3_ub] .
1774 finally show "ln (1 / 3 + 1) \<le> ?uthird * ub_ln_horner prec (get_odd prec) 1 ?uthird" .
1776 show ?lb_ln2 unfolding lb_ln2_def Let_def real_of_float_add ln2_sum Float_num(4)[symmetric]
1777 proof (rule add_mono, fact ln_float_bounds(1)[OF lb2 ub2])
1778 have "?lthird * lb_ln_horner prec (get_even prec) 1 ?lthird \<le> ln (real ?lthird + 1)"
1779 using ln_float_bounds(1)[OF lb3_lb lb3_ub] .
1780 also have "\<dots> \<le> ln (1 / 3 + 1)" unfolding ln_le_cancel_iff[OF lthird_gt0 third_gt0] using lb3 by auto
1781 finally show "?lthird * lb_ln_horner prec (get_even prec) 1 ?lthird \<le> ln (1 / 3 + 1)" .
1785 subsection "Compute the logarithm in the entire domain"
1787 function ub_ln :: "nat \<Rightarrow> float \<Rightarrow> float option" and lb_ln :: "nat \<Rightarrow> float \<Rightarrow> float option" where
1788 "ub_ln prec x = (if x \<le> 0 then None
1789 else if x < 1 then Some (- the (lb_ln prec (float_divl (max prec 1) 1 x)))
1790 else let horner = \<lambda>x. x * ub_ln_horner prec (get_odd prec) 1 x in
1791 if x \<le> Float 3 -1 then Some (horner (x - 1))
1792 else if x < Float 1 1 then Some (horner (Float 1 -1) + horner (x * rapprox_rat prec 2 3 - 1))
1793 else let l = bitlen (mantissa x) - 1 in
1794 Some (ub_ln2 prec * (Float (scale x + l) 0) + horner (Float (mantissa x) (- l) - 1)))" |
1795 "lb_ln prec x = (if x \<le> 0 then None
1796 else if x < 1 then Some (- the (ub_ln prec (float_divr prec 1 x)))
1797 else let horner = \<lambda>x. x * lb_ln_horner prec (get_even prec) 1 x in
1798 if x \<le> Float 3 -1 then Some (horner (x - 1))
1799 else if x < Float 1 1 then Some (horner (Float 1 -1) +
1800 horner (max (x * lapprox_rat prec 2 3 - 1) 0))
1801 else let l = bitlen (mantissa x) - 1 in
1802 Some (lb_ln2 prec * (Float (scale x + l) 0) + horner (Float (mantissa x) (- l) - 1)))"
1803 by pat_completeness auto
1805 termination proof (relation "measure (\<lambda> v. let (prec, x) = sum_case id id v in (if x < 1 then 1 else 0))", auto)
1806 fix prec x assume "\<not> x \<le> 0" and "x < 1" and "float_divl (max prec (Suc 0)) 1 x < 1"
1807 hence "0 < x" and "0 < max prec (Suc 0)" unfolding less_float_def le_float_def by auto
1808 from float_divl_pos_less1_bound[OF `0 < x` `x < 1` `0 < max prec (Suc 0)`]
1809 show False using `float_divl (max prec (Suc 0)) 1 x < 1` unfolding less_float_def le_float_def by auto
1811 fix prec x assume "\<not> x \<le> 0" and "x < 1" and "float_divr prec 1 x < 1"
1812 hence "0 < x" unfolding less_float_def le_float_def by auto
1813 from float_divr_pos_less1_lower_bound[OF `0 < x` `x < 1`, of prec]
1814 show False using `float_divr prec 1 x < 1` unfolding less_float_def le_float_def by auto
1817 lemma ln_shifted_float: assumes "0 < m" shows "ln (Float m e) = ln 2 * (e + (bitlen m - 1)) + ln (Float m (- (bitlen m - 1)))"
1819 let ?B = "2^nat (bitlen m - 1)"
1820 have "0 < real m" and "\<And>X. (0 :: real) < 2^X" and "0 < (2 :: real)" and "m \<noteq> 0" using assms by auto
1821 hence "0 \<le> bitlen m - 1" using bitlen_ge1[OF `m \<noteq> 0`] by auto
1823 proof (cases "0 \<le> e")
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_ge0_exp[OF True] ln_mult[OF `0 < real m` `0 < 2^nat e`]
1828 ln_realpow[OF `0 < 2`] algebra_simps using `0 \<le> bitlen m - 1` True by auto
1830 case False hence "0 < -e" by auto
1831 hence pow_gt0: "(0::real) < 2^nat (-e)" by auto
1832 hence inv_gt0: "(0::real) < inverse (2^nat (-e))" by auto
1833 show ?thesis unfolding normalized_float[OF `m \<noteq> 0`]
1834 unfolding ln_div[OF `0 < real m` `0 < ?B`] real_of_int_add ln_realpow[OF `0 < 2`]
1835 unfolding real_of_float_nge0_exp[OF False] ln_mult[OF `0 < real m` inv_gt0] ln_inverse[OF pow_gt0]
1836 ln_realpow[OF `0 < 2`] algebra_simps using `0 \<le> bitlen m - 1` False by auto
1840 lemma ub_ln_lb_ln_bounds': assumes "1 \<le> x"
1841 shows "the (lb_ln prec x) \<le> ln x \<and> ln x \<le> the (ub_ln prec x)"
1842 (is "?lb \<le> ?ln \<and> ?ln \<le> ?ub")
1843 proof (cases "x < Float 1 1")
1845 hence "real (x - 1) < 1" and "real x < 2" unfolding less_float_def Float_num by auto
1846 have "\<not> x \<le> 0" and "\<not> x < 1" using `1 \<le> x` unfolding less_float_def le_float_def by auto
1847 hence "0 \<le> real (x - 1)" using `1 \<le> x` unfolding less_float_def Float_num by auto
1849 have [simp]: "(Float 3 -1) = 3 / 2" by (simp add: real_of_float_def pow2_def)
1852 proof (cases "x \<le> Float 3 -1")
1854 show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps Let_def
1855 using ln_float_bounds[OF `0 \<le> real (x - 1)` `real (x - 1) < 1`, of prec] `\<not> x \<le> 0` `\<not> x < 1` True
1858 case False hence *: "3 / 2 < x" by (auto simp add: le_float_def)
1860 with ln_add[of "3 / 2" "x - 3 / 2"]
1861 have add: "ln x = ln (3 / 2) + ln (real x * 2 / 3)"
1862 by (auto simp add: algebra_simps diff_divide_distrib)
1864 let "?ub_horner x" = "x * ub_ln_horner prec (get_odd prec) 1 x"
1865 let "?lb_horner x" = "x * lb_ln_horner prec (get_even prec) 1 x"
1867 { have up: "real (rapprox_rat prec 2 3) \<le> 1"
1868 by (rule rapprox_rat_le1) simp_all
1869 have low: "2 / 3 \<le> rapprox_rat prec 2 3"
1870 by (rule order_trans[OF _ rapprox_rat]) simp
1871 from mult_less_le_imp_less[OF * low] *
1872 have pos: "0 < real (x * rapprox_rat prec 2 3 - 1)" by auto
1874 have "ln (real x * 2/3)
1875 \<le> ln (real (x * rapprox_rat prec 2 3 - 1) + 1)"
1876 proof (rule ln_le_cancel_iff[symmetric, THEN iffD1])
1877 show "real x * 2 / 3 \<le> real (x * rapprox_rat prec 2 3 - 1) + 1"
1879 show "0 < real x * 2 / 3" using * by simp
1880 show "0 < real (x * rapprox_rat prec 2 3 - 1) + 1" using pos by auto
1882 also have "\<dots> \<le> ?ub_horner (x * rapprox_rat prec 2 3 - 1)"
1883 proof (rule ln_float_bounds(2))
1884 from mult_less_le_imp_less[OF `real x < 2` up] low *
1885 show "real (x * rapprox_rat prec 2 3 - 1) < 1" by auto
1886 show "0 \<le> real (x * rapprox_rat prec 2 3 - 1)" using pos by auto
1889 \<le> ?ub_horner (Float 1 -1)
1890 + ?ub_horner (x * rapprox_rat prec 2 3 - 1)"
1891 using ln_float_bounds(2)[of "Float 1 -1" prec prec] add by auto }
1893 { let ?max = "max (x * lapprox_rat prec 2 3 - 1) 0"
1895 have up: "lapprox_rat prec 2 3 \<le> 2/3"
1896 by (rule order_trans[OF lapprox_rat], simp)
1898 have low: "0 \<le> real (lapprox_rat prec 2 3)"
1899 using lapprox_rat_bottom[of 2 3 prec] by simp
1901 have "?lb_horner ?max
1902 \<le> ln (real ?max + 1)"
1903 proof (rule ln_float_bounds(1))
1904 from mult_less_le_imp_less[OF `real x < 2` up] * low
1905 show "real ?max < 1" by (cases "real (lapprox_rat prec 2 3) = 0",
1906 auto simp add: real_of_float_max)
1907 show "0 \<le> real ?max" by (auto simp add: real_of_float_max)
1909 also have "\<dots> \<le> ln (real x * 2/3)"
1910 proof (rule ln_le_cancel_iff[symmetric, THEN iffD1])
1911 show "0 < real ?max + 1" by (auto simp add: real_of_float_max)
1912 show "0 < real x * 2/3" using * by auto
1913 show "real ?max + 1 \<le> real x * 2/3" using * up
1914 by (cases "0 < real x * real (lapprox_posrat prec 2 3) - 1",
1915 auto simp add: real_of_float_max min_max.sup_absorb1)
1917 finally have "?lb_horner (Float 1 -1) + ?lb_horner ?max
1919 using ln_float_bounds(1)[of "Float 1 -1" prec prec] add by auto }
1921 show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps Let_def
1922 using `\<not> x \<le> 0` `\<not> x < 1` True False by auto
1926 hence "\<not> x \<le> 0" and "\<not> x < 1" "0 < x" "\<not> x \<le> Float 3 -1"
1927 using `1 \<le> x` unfolding less_float_def le_float_def real_of_float_simp pow2_def
1932 let ?s = "Float (e + (bitlen m - 1)) 0"
1933 let ?x = "Float m (- (bitlen m - 1))"
1935 have "0 < m" and "m \<noteq> 0" using float_pos_m_pos `0 < x` Float by auto
1938 have "lb_ln2 prec * ?s \<le> ln 2 * (e + (bitlen m - 1))" (is "?lb2 \<le> _")
1939 unfolding real_of_float_mult real_of_float_ge0_exp[OF order_refl] nat_0 power_0 mult_1_right
1940 using lb_ln2[of prec]
1941 proof (rule mult_right_mono)
1942 have "1 \<le> Float m e" using `1 \<le> x` Float unfolding le_float_def by auto
1943 from float_gt1_scale[OF this]
1944 show "0 \<le> real (e + (bitlen m - 1))" by auto
1947 from bitlen_div[OF `0 < m`, unfolded normalized_float[OF `m \<noteq> 0`, symmetric]]
1948 have "0 \<le> real (?x - 1)" and "real (?x - 1) < 1" by auto
1949 from ln_float_bounds(1)[OF this]
1950 have "(?x - 1) * lb_ln_horner prec (get_even prec) 1 (?x - 1) \<le> ln ?x" (is "?lb_horner \<le> _") by auto
1951 ultimately have "?lb2 + ?lb_horner \<le> ln x"
1952 unfolding Float ln_shifted_float[OF `0 < m`, of e] by auto
1956 from bitlen_div[OF `0 < m`, unfolded normalized_float[OF `m \<noteq> 0`, symmetric]]
1957 have "0 \<le> real (?x - 1)" and "real (?x - 1) < 1" by auto
1958 from ln_float_bounds(2)[OF this]
1959 have "ln ?x \<le> (?x - 1) * ub_ln_horner prec (get_odd prec) 1 (?x - 1)" (is "_ \<le> ?ub_horner") by auto
1961 have "ln 2 * (e + (bitlen m - 1)) \<le> ub_ln2 prec * ?s" (is "_ \<le> ?ub2")
1962 unfolding real_of_float_mult real_of_float_ge0_exp[OF order_refl] nat_0 power_0 mult_1_right
1963 using ub_ln2[of prec]
1964 proof (rule mult_right_mono)
1965 have "1 \<le> Float m e" using `1 \<le> x` Float unfolding le_float_def by auto
1966 from float_gt1_scale[OF this]
1967 show "0 \<le> real (e + (bitlen m - 1))" by auto
1969 ultimately have "ln x \<le> ?ub2 + ?ub_horner"
1970 unfolding Float ln_shifted_float[OF `0 < m`, of e] by auto
1972 ultimately show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps
1973 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
1974 unfolding scale.simps[of m e, unfolded Float[symmetric]] mantissa.simps[of m e, unfolded Float[symmetric]] real_of_float_add
1979 lemma ub_ln_lb_ln_bounds: assumes "0 < x"
1980 shows "the (lb_ln prec x) \<le> ln x \<and> ln x \<le> the (ub_ln prec x)"
1981 (is "?lb \<le> ?ln \<and> ?ln \<le> ?ub")
1982 proof (cases "x < 1")
1983 case False hence "1 \<le> x" unfolding less_float_def le_float_def by auto
1984 show ?thesis using ub_ln_lb_ln_bounds'[OF `1 \<le> x`] .
1986 case True have "\<not> x \<le> 0" using `0 < x` unfolding less_float_def le_float_def by auto
1988 have "0 < real x" and "real x \<noteq> 0" using `0 < x` unfolding less_float_def by auto
1989 hence A: "0 < 1 / real x" by auto
1992 let ?divl = "float_divl (max prec 1) 1 x"
1993 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
1994 hence B: "0 < real ?divl" unfolding le_float_def by auto
1996 have "ln ?divl \<le> ln (1 / x)" unfolding ln_le_cancel_iff[OF B A] using float_divl[of _ 1 x] by auto
1997 hence "ln x \<le> - ln ?divl" unfolding nonzero_inverse_eq_divide[OF `real x \<noteq> 0`, symmetric] ln_inverse[OF `0 < real x`] by auto
1998 from this ub_ln_lb_ln_bounds'[OF A', THEN conjunct1, THEN le_imp_neg_le]
1999 have "?ln \<le> - the (lb_ln prec ?divl)" unfolding real_of_float_minus by (rule order_trans)
2002 let ?divr = "float_divr prec 1 x"
2003 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
2004 hence B: "0 < real ?divr" unfolding le_float_def by auto
2006 have "ln (1 / x) \<le> ln ?divr" unfolding ln_le_cancel_iff[OF A B] using float_divr[of 1 x] by auto
2007 hence "- ln ?divr \<le> ln x" unfolding nonzero_inverse_eq_divide[OF `real x \<noteq> 0`, symmetric] ln_inverse[OF `0 < real x`] by auto
2008 from ub_ln_lb_ln_bounds'[OF A', THEN conjunct2, THEN le_imp_neg_le] this
2009 have "- the (ub_ln prec ?divr) \<le> ?ln" unfolding real_of_float_minus by (rule order_trans)
2011 ultimately show ?thesis unfolding lb_ln.simps[where x=x] ub_ln.simps[where x=x]
2012 unfolding if_not_P[OF `\<not> x \<le> 0`] if_P[OF True] by auto
2015 lemma lb_ln: assumes "Some y = lb_ln prec x"
2016 shows "y \<le> ln x" and "0 < real x"
2020 assume "\<not> 0 < x" hence "x \<le> 0" unfolding le_float_def less_float_def by auto
2021 thus False using assms by auto
2023 thus "0 < real x" unfolding less_float_def by auto
2024 have "the (lb_ln prec x) \<le> ln x" using ub_ln_lb_ln_bounds[OF `0 < x`] ..
2025 thus "y \<le> ln x" unfolding assms[symmetric] by auto
2028 lemma ub_ln: assumes "Some y = ub_ln prec x"
2029 shows "ln x \<le> y" and "0 < real x"
2033 assume "\<not> 0 < x" hence "x \<le> 0" unfolding le_float_def less_float_def by auto
2034 thus False using assms by auto
2036 thus "0 < real x" unfolding less_float_def by auto
2037 have "ln x \<le> the (ub_ln prec x)" using ub_ln_lb_ln_bounds[OF `0 < x`] ..
2038 thus "ln x \<le> y" unfolding assms[symmetric] by auto
2041 lemma bnds_ln: "\<forall> (x::real) lx ux. (Some l, Some u) = (lb_ln prec lx, ub_ln prec ux) \<and> x \<in> {lx .. ux} \<longrightarrow> l \<le> ln x \<and> ln x \<le> u"
2042 proof (rule allI, rule allI, rule allI, rule impI)
2043 fix x::real and lx ux
2044 assume "(Some l, Some u) = (lb_ln prec lx, ub_ln prec ux) \<and> x \<in> {lx .. ux}"
2045 hence l: "Some l = lb_ln prec lx " and u: "Some u = ub_ln prec ux" and x: "x \<in> {lx .. ux}" by auto
2047 have "ln ux \<le> u" and "0 < real ux" using ub_ln u by auto
2048 have "l \<le> ln lx" and "0 < real lx" and "0 < x" using lb_ln[OF l] x by auto
2050 from ln_le_cancel_iff[OF `0 < real lx` `0 < x`] `l \<le> ln lx`
2051 have "l \<le> ln x" using x unfolding atLeastAtMost_iff by auto
2053 from ln_le_cancel_iff[OF `0 < x` `0 < real ux`] `ln ux \<le> real u`
2054 have "ln x \<le> u" using x unfolding atLeastAtMost_iff by auto
2055 ultimately show "l \<le> ln x \<and> ln x \<le> u" ..
2058 section "Implement floatarith"
2060 subsection "Define syntax and semantics"
2063 = Add floatarith floatarith
2065 | Mult floatarith floatarith
2066 | Inverse floatarith
2070 | Max floatarith floatarith
2071 | Min floatarith floatarith
2076 | Power floatarith nat
2080 fun interpret_floatarith :: "floatarith \<Rightarrow> real list \<Rightarrow> real" where
2081 "interpret_floatarith (Add a b) vs = (interpret_floatarith a vs) + (interpret_floatarith b vs)" |
2082 "interpret_floatarith (Minus a) vs = - (interpret_floatarith a vs)" |
2083 "interpret_floatarith (Mult a b) vs = (interpret_floatarith a vs) * (interpret_floatarith b vs)" |
2084 "interpret_floatarith (Inverse a) vs = inverse (interpret_floatarith a vs)" |
2085 "interpret_floatarith (Cos a) vs = cos (interpret_floatarith a vs)" |
2086 "interpret_floatarith (Arctan a) vs = arctan (interpret_floatarith a vs)" |
2087 "interpret_floatarith (Min a b) vs = min (interpret_floatarith a vs) (interpret_floatarith b vs)" |
2088 "interpret_floatarith (Max a b) vs = max (interpret_floatarith a vs) (interpret_floatarith b vs)" |
2089 "interpret_floatarith (Abs a) vs = abs (interpret_floatarith a vs)" |
2090 "interpret_floatarith Pi vs = pi" |
2091 "interpret_floatarith (Sqrt a) vs = sqrt (interpret_floatarith a vs)" |
2092 "interpret_floatarith (Exp a) vs = exp (interpret_floatarith a vs)" |
2093 "interpret_floatarith (Ln a) vs = ln (interpret_floatarith a vs)" |
2094 "interpret_floatarith (Power a n) vs = (interpret_floatarith a vs)^n" |
2095 "interpret_floatarith (Num f) vs = f" |
2096 "interpret_floatarith (Var n) vs = vs ! n"
2098 lemma interpret_floatarith_divide: "interpret_floatarith (Mult a (Inverse b)) vs = (interpret_floatarith a vs) / (interpret_floatarith b vs)"
2099 unfolding divide_inverse interpret_floatarith.simps ..
2101 lemma interpret_floatarith_diff: "interpret_floatarith (Add a (Minus b)) vs = (interpret_floatarith a vs) - (interpret_floatarith b vs)"
2102 unfolding diff_minus interpret_floatarith.simps ..
2104 lemma interpret_floatarith_sin: "interpret_floatarith (Cos (Add (Mult Pi (Num (Float 1 -1))) (Minus a))) vs =
2105 sin (interpret_floatarith a vs)"
2106 unfolding sin_cos_eq interpret_floatarith.simps
2107 interpret_floatarith_divide interpret_floatarith_diff diff_minus
2110 lemma interpret_floatarith_tan:
2111 "interpret_floatarith (Mult (Cos (Add (Mult Pi (Num (Float 1 -1))) (Minus a))) (Inverse (Cos a))) vs =
2112 tan (interpret_floatarith a vs)"
2113 unfolding interpret_floatarith.simps(3,4) interpret_floatarith_sin tan_def divide_inverse
2116 lemma interpret_floatarith_powr: "interpret_floatarith (Exp (Mult b (Ln a))) vs = (interpret_floatarith a vs) powr (interpret_floatarith b vs)"
2117 unfolding powr_def interpret_floatarith.simps ..
2119 lemma interpret_floatarith_log: "interpret_floatarith ((Mult (Ln x) (Inverse (Ln b)))) vs = log (interpret_floatarith b vs) (interpret_floatarith x vs)"
2120 unfolding log_def interpret_floatarith.simps divide_inverse ..
2122 lemma interpret_floatarith_num:
2123 shows "interpret_floatarith (Num (Float 0 0)) vs = 0"
2124 and "interpret_floatarith (Num (Float 1 0)) vs = 1"
2125 and "interpret_floatarith (Num (Float (number_of a) 0)) vs = number_of a" by auto
2127 subsection "Implement approximation function"
2129 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
2130 "lift_bin' (Some (l1, u1)) (Some (l2, u2)) f = Some (f l1 u1 l2 u2)" |
2131 "lift_bin' a b f = None"
2133 fun lift_un :: "(float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> ((float option) * (float option))) \<Rightarrow> (float * float) option" where
2134 "lift_un (Some (l1, u1)) f = (case (f l1 u1) of (Some l, Some u) \<Rightarrow> Some (l, u)
2135 | t \<Rightarrow> None)" |
2136 "lift_un b f = None"
2138 fun lift_un' :: "(float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> (float * float)) \<Rightarrow> (float * float) option" where
2139 "lift_un' (Some (l1, u1)) f = Some (f l1 u1)" |
2140 "lift_un' b f = None"
2143 "bounded_by xs vs \<longleftrightarrow>
2144 (\<forall> i < length vs. case vs ! i of None \<Rightarrow> True
2145 | Some (l, u) \<Rightarrow> xs ! i \<in> { real l .. real u })"
2148 assumes "bounded_by xs vs"
2149 shows "\<And> i. i < length vs \<Longrightarrow> case vs ! i of None \<Rightarrow> True
2150 | Some (l, u) \<Rightarrow> xs ! i \<in> { real l .. real u }"
2151 using assms bounded_by_def by blast
2153 lemma bounded_by_update:
2154 assumes "bounded_by xs vs"
2155 and bnd: "xs ! i \<in> { real l .. real u }"
2156 shows "bounded_by xs (vs[i := Some (l,u)])"
2159 let ?vs = "vs[i := Some (l,u)]"
2160 assume "j < length ?vs" hence [simp]: "j < length vs" by simp
2161 have "case ?vs ! j of None \<Rightarrow> True | Some (l, u) \<Rightarrow> xs ! j \<in> { real l .. real u }"
2162 proof (cases "?vs ! j")
2165 proof (cases "i = j")
2167 thus ?thesis using `?vs ! j = Some b` and bnd by auto
2170 thus ?thesis using `bounded_by xs vs` unfolding bounded_by_def by auto
2173 thus ?thesis unfolding bounded_by_def by auto
2176 lemma bounded_by_None:
2177 shows "bounded_by xs (replicate (length xs) None)"
2178 unfolding bounded_by_def by auto
2180 fun approx approx' :: "nat \<Rightarrow> floatarith \<Rightarrow> (float * float) option list \<Rightarrow> (float * float) option" where
2181 "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)" |
2182 "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))" |
2183 "approx prec (Minus a) bs = lift_un' (approx' prec a bs) (\<lambda> l u. (-u, -l))" |
2184 "approx prec (Mult a b) bs = lift_bin' (approx' prec a bs) (approx' prec b bs)
2185 (\<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,
2186 float_pprt a2 * float_pprt b2 + float_pprt a1 * float_nprt b2 + float_nprt a2 * float_pprt b1 + float_nprt a1 * float_nprt b1))" |
2187 "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))" |
2188 "approx prec (Cos a) bs = lift_un' (approx' prec a bs) (bnds_cos prec)" |
2189 "approx prec Pi bs = Some (lb_pi prec, ub_pi prec)" |
2190 "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))" |
2191 "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))" |
2192 "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>))" |
2193 "approx prec (Arctan a) bs = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_arctan prec l, ub_arctan prec u))" |
2194 "approx prec (Sqrt a) bs = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_sqrt prec l, ub_sqrt prec u))" |
2195 "approx prec (Exp a) bs = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_exp prec l, ub_exp prec u))" |
2196 "approx prec (Ln a) bs = lift_un (approx' prec a bs) (\<lambda> l u. (lb_ln prec l, ub_ln prec u))" |
2197 "approx prec (Power a n) bs = lift_un' (approx' prec a bs) (float_power_bnds n)" |
2198 "approx prec (Num f) bs = Some (f, f)" |
2199 "approx prec (Var i) bs = (if i < length bs then bs ! i else None)"
2202 assumes lift_bin'_Some: "Some (l, u) = lift_bin' a b f"
2203 shows "\<exists> l1 u1 l2 u2. Some (l1, u1) = a \<and> Some (l2, u2) = b"
2205 case None hence "None = lift_bin' a b f" unfolding None lift_bin'.simps ..
2206 thus ?thesis using lift_bin'_Some by auto
2211 case None hence "None = lift_bin' a b f" unfolding None lift_bin'.simps ..
2212 thus ?thesis using lift_bin'_Some by auto
2215 obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
2216 obtain lb ub where b': "b' = (lb, ub)" by (cases b', auto)
2217 thus ?thesis unfolding `a = Some a'` `b = Some b'` a' b' by auto
2222 assumes lift_bin'_Some: "Some (l, u) = lift_bin' (g a) (g b) f"
2223 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"
2224 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)"
2227 where Sa: "Some (l1, u1) = g a" and Sb: "Some (l2, u2) = g b" using lift_bin'_ex[OF assms(1)] by auto
2228 have lu: "(l, u) = f l1 u1 l2 u2" using lift_bin'_Some[unfolded Sa[symmetric] Sb[symmetric] lift_bin'.simps] by auto
2229 have "l = fst (f l1 u1 l2 u2)" and "u = snd (f l1 u1 l2 u2)" unfolding lu[symmetric] by auto
2230 thus ?thesis using Pa[OF Sa] Pb[OF Sb] by auto
2233 lemma approx_approx':
2234 assumes Pa: "\<And>l u. Some (l, u) = approx prec a vs \<Longrightarrow> l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> u"
2235 and approx': "Some (l, u) = approx' prec a vs"
2236 shows "l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> u"
2238 obtain l' u' where S: "Some (l', u') = approx prec a vs"
2239 using approx' unfolding approx'.simps by (cases "approx prec a vs", auto)
2240 have l': "l = round_down prec l'" and u': "u = round_up prec u'"
2241 using approx' unfolding approx'.simps S[symmetric] by auto
2242 show ?thesis unfolding l' u'
2243 using order_trans[OF Pa[OF S, THEN conjunct2] round_up[of u']]
2244 using order_trans[OF round_down[of _ l'] Pa[OF S, THEN conjunct1]] by auto
2248 assumes lift_bin'_Some: "Some (l, u) = lift_bin' (approx' prec a bs) (approx' prec b bs) f"
2249 and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow> l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> u" (is "\<And>l u. _ = ?g a \<Longrightarrow> ?P l u a")
2250 and Pb: "\<And>l u. Some (l, u) = approx prec b bs \<Longrightarrow> l \<le> interpret_floatarith b xs \<and> interpret_floatarith b xs \<le> u"
2251 shows "\<exists> l1 u1 l2 u2. (l1 \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> u1) \<and>
2252 (l2 \<le> interpret_floatarith b xs \<and> interpret_floatarith b xs \<le> u2) \<and>
2253 l = fst (f l1 u1 l2 u2) \<and> u = snd (f l1 u1 l2 u2)"
2255 { fix l u assume "Some (l, u) = approx' prec a bs"
2256 with approx_approx'[of prec a bs, OF _ this] Pa
2257 have "l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> u" by auto } note Pa = this
2258 { fix l u assume "Some (l, u) = approx' prec b bs"
2259 with approx_approx'[of prec b bs, OF _ this] Pb
2260 have "l \<le> interpret_floatarith b xs \<and> interpret_floatarith b xs \<le> u" by auto } note Pb = this
2262 from lift_bin'_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_bin'_Some, OF Pa Pb]
2263 show ?thesis by auto
2267 assumes lift_un'_Some: "Some (l, u) = lift_un' a f"
2268 shows "\<exists> l u. Some (l, u) = a"
2270 case None hence "None = lift_un' a f" unfolding None lift_un'.simps ..
2271 thus ?thesis using lift_un'_Some by auto
2274 obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
2275 thus ?thesis unfolding `a = Some a'` a' by auto
2279 assumes lift_un'_Some: "Some (l, u) = lift_un' (g a) f"
2280 and Pa: "\<And>l u. Some (l, u) = g a \<Longrightarrow> P l u a"
2281 shows "\<exists> l1 u1. P l1 u1 a \<and> l = fst (f l1 u1) \<and> u = snd (f l1 u1)"
2283 obtain l1 u1 where Sa: "Some (l1, u1) = g a" using lift_un'_ex[OF assms(1)] by auto
2284 have lu: "(l, u) = f l1 u1" using lift_un'_Some[unfolded Sa[symmetric] lift_un'.simps] by auto
2285 have "l = fst (f l1 u1)" and "u = snd (f l1 u1)" unfolding lu[symmetric] by auto
2286 thus ?thesis using Pa[OF Sa] by auto
2290 assumes lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f"
2291 and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow> l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> u" (is "\<And>l u. _ = ?g a \<Longrightarrow> ?P l u a")
2292 shows "\<exists> l1 u1. (l1 \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> u1) \<and>
2293 l = fst (f l1 u1) \<and> u = snd (f l1 u1)"
2295 { fix l u assume "Some (l, u) = approx' prec a bs"
2296 with approx_approx'[of prec a bs, OF _ this] Pa
2297 have "l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> u" by auto } note Pa = this
2298 from lift_un'_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_un'_Some, OF Pa]
2299 show ?thesis by auto
2302 lemma lift_un'_bnds:
2303 assumes bnds: "\<forall> (x::real) lx ux. (l, u) = f lx ux \<and> x \<in> { lx .. ux } \<longrightarrow> l \<le> f' x \<and> f' x \<le> u"
2304 and lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f"
2305 and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow> l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> u"
2306 shows "real l \<le> f' (interpret_floatarith a xs) \<and> f' (interpret_floatarith a xs) \<le> real u"
2308 from lift_un'[OF lift_un'_Some Pa]
2309 obtain l1 u1 where "l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> u1" and "l = fst (f l1 u1)" and "u = snd (f l1 u1)" by blast
2310 hence "(l, u) = f l1 u1" and "interpret_floatarith a xs \<in> {l1 .. u1}" by auto
2311 thus ?thesis using bnds by auto
2315 assumes lift_un_Some: "Some (l, u) = lift_un a f"
2316 shows "\<exists> l u. Some (l, u) = a"
2318 case None hence "None = lift_un a f" unfolding None lift_un.simps ..
2319 thus ?thesis using lift_un_Some by auto
2322 obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
2323 thus ?thesis unfolding `a = Some a'` a' by auto
2327 assumes lift_un_Some: "Some (l, u) = lift_un (g a) f"
2328 and Pa: "\<And>l u. Some (l, u) = g a \<Longrightarrow> P l u a"
2329 shows "\<exists> l1 u1. P l1 u1 a \<and> Some l = fst (f l1 u1) \<and> Some u = snd (f l1 u1)"
2331 obtain l1 u1 where Sa: "Some (l1, u1) = g a" using lift_un_ex[OF assms(1)] by auto
2332 have "fst (f l1 u1) \<noteq> None \<and> snd (f l1 u1) \<noteq> None"
2334 assume "\<not> (fst (f l1 u1) \<noteq> None \<and> snd (f l1 u1) \<noteq> None)"
2335 hence or: "fst (f l1 u1) = None \<or> snd (f l1 u1) = None" by auto
2336 hence "lift_un (g a) f = None"
2337 proof (cases "fst (f l1 u1) = None")
2339 then obtain b where b: "f l1 u1 = (None, b)" by (cases "f l1 u1", auto)
2340 thus ?thesis unfolding Sa[symmetric] lift_un.simps b by auto
2342 case False hence "snd (f l1 u1) = None" using or by auto
2343 with False obtain b where b: "f l1 u1 = (Some b, None)" by (cases "f l1 u1", auto)
2344 thus ?thesis unfolding Sa[symmetric] lift_un.simps b by auto
2346 thus False using lift_un_Some by auto
2348 then obtain a' b' where f: "f l1 u1 = (Some a', Some b')" by (cases "f l1 u1", auto)
2349 from lift_un_Some[unfolded Sa[symmetric] lift_un.simps f]
2350 have "Some l = fst (f l1 u1)" and "Some u = snd (f l1 u1)" unfolding f by auto
2351 thus ?thesis unfolding Sa[symmetric] lift_un.simps using Pa[OF Sa] by auto
2355 assumes lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f"
2356 and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow> l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> u" (is "\<And>l u. _ = ?g a \<Longrightarrow> ?P l u a")
2357 shows "\<exists> l1 u1. (l1 \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> u1) \<and>
2358 Some l = fst (f l1 u1) \<and> Some u = snd (f l1 u1)"
2360 { fix l u assume "Some (l, u) = approx' prec a bs"
2361 with approx_approx'[of prec a bs, OF _ this] Pa
2362 have "l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> u" by auto } note Pa = this
2363 from lift_un_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_un_Some, OF Pa]
2364 show ?thesis by auto
2368 assumes bnds: "\<forall> (x::real) lx ux. (Some l, Some u) = f lx ux \<and> x \<in> { lx .. ux } \<longrightarrow> l \<le> f' x \<and> f' x \<le> u"
2369 and lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f"
2370 and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow> l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> u"
2371 shows "real l \<le> f' (interpret_floatarith a xs) \<and> f' (interpret_floatarith a xs) \<le> real u"
2373 from lift_un[OF lift_un_Some Pa]
2374 obtain l1 u1 where "l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> u1" and "Some l = fst (f l1 u1)" and "Some u = snd (f l1 u1)" by blast
2375 hence "(Some l, Some u) = f l1 u1" and "interpret_floatarith a xs \<in> {l1 .. u1}" by auto
2376 thus ?thesis using bnds by auto
2380 assumes "bounded_by xs vs"
2381 and "Some (l, u) = approx prec arith vs" (is "_ = ?g arith")
2382 shows "l \<le> interpret_floatarith arith xs \<and> interpret_floatarith arith xs \<le> u" (is "?P l u arith")
2383 using `Some (l, u) = approx prec arith vs`
2384 proof (induct arith arbitrary: l u x)
2386 from lift_bin'[OF Add.prems[unfolded approx.simps]] Add.hyps
2387 obtain l1 u1 l2 u2 where "l = l1 + l2" and "u = u1 + u2"
2388 "l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> u1"
2389 "l2 \<le> interpret_floatarith b xs" and "interpret_floatarith b xs \<le> u2" unfolding fst_conv snd_conv by blast
2390 thus ?case unfolding interpret_floatarith.simps by auto
2393 from lift_un'[OF Minus.prems[unfolded approx.simps]] Minus.hyps
2394 obtain l1 u1 where "l = -u1" and "u = -l1"
2395 "l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> u1" unfolding fst_conv snd_conv by blast
2396 thus ?case unfolding interpret_floatarith.simps using real_of_float_minus by auto
2399 from lift_bin'[OF Mult.prems[unfolded approx.simps]] Mult.hyps
2401 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"
2402 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"
2403 and "l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> u1"
2404 and "l2 \<le> interpret_floatarith b xs" and "interpret_floatarith b xs \<le> u2" unfolding fst_conv snd_conv by blast
2405 thus ?case unfolding interpret_floatarith.simps l u real_of_float_add real_of_float_mult real_of_float_nprt real_of_float_pprt
2406 using mult_le_prts mult_ge_prts by auto
2409 from lift_un[OF Inverse.prems[unfolded approx.simps], unfolded if_distrib[of fst] if_distrib[of snd] fst_conv snd_conv] Inverse.hyps
2410 obtain l1 u1 where l': "Some l = (if 0 < l1 \<or> u1 < 0 then Some (float_divl prec 1 u1) else None)"
2411 and u': "Some u = (if 0 < l1 \<or> u1 < 0 then Some (float_divr prec 1 l1) else None)"
2412 and l1: "l1 \<le> interpret_floatarith a xs" and u1: "interpret_floatarith a xs \<le> u1" by blast
2413 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
2414 moreover have l1_le_u1: "real l1 \<le> real u1" using l1 u1 by auto
2415 ultimately have "real l1 \<noteq> 0" and "real u1 \<noteq> 0" unfolding less_float_def by auto
2417 have inv: "inverse u1 \<le> inverse (interpret_floatarith a xs)
2418 \<and> inverse (interpret_floatarith a xs) \<le> inverse l1"
2419 proof (cases "0 < l1")
2420 case True hence "0 < real u1" and "0 < real l1" "0 < interpret_floatarith a xs"
2421 unfolding less_float_def using l1_le_u1 l1 by auto
2423 unfolding inverse_le_iff_le[OF `0 < real u1` `0 < interpret_floatarith a xs`]
2424 inverse_le_iff_le[OF `0 < interpret_floatarith a xs` `0 < real l1`]
2427 case False hence "u1 < 0" using either by blast
2428 hence "real u1 < 0" and "real l1 < 0" "interpret_floatarith a xs < 0"
2429 unfolding less_float_def using l1_le_u1 u1 by auto
2431 unfolding inverse_le_iff_le_neg[OF `real u1 < 0` `interpret_floatarith a xs < 0`]
2432 inverse_le_iff_le_neg[OF `interpret_floatarith a xs < 0` `real l1 < 0`]
2436 from l' have "l = float_divl prec 1 u1" by (cases "0 < l1 \<or> u1 < 0", auto)
2437 hence "l \<le> inverse u1" unfolding nonzero_inverse_eq_divide[OF `real u1 \<noteq> 0`] using float_divl[of prec 1 u1] by auto
2438 also have "\<dots> \<le> inverse (interpret_floatarith a xs)" using inv by auto
2439 finally have "l \<le> inverse (interpret_floatarith a xs)" .
2441 from u' have "u = float_divr prec 1 l1" by (cases "0 < l1 \<or> u1 < 0", auto)
2442 hence "inverse l1 \<le> u" unfolding nonzero_inverse_eq_divide[OF `real l1 \<noteq> 0`] using float_divr[of 1 l1 prec] by auto
2443 hence "inverse (interpret_floatarith a xs) \<le> u" by (rule order_trans[OF inv[THEN conjunct2]])
2444 ultimately show ?case unfolding interpret_floatarith.simps using l1 u1 by auto
2447 from lift_un'[OF Abs.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Abs.hyps
2448 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>"
2449 and l1: "l1 \<le> interpret_floatarith x xs" and u1: "interpret_floatarith x xs \<le> u1" by blast
2450 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)
2453 from lift_bin'[OF Min.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Min.hyps
2454 obtain l1 u1 l2 u2 where l': "l = min l1 l2" and u': "u = min u1 u2"
2455 and l1: "l1 \<le> interpret_floatarith a xs" and u1: "interpret_floatarith a xs \<le> u1"
2456 and l1: "l2 \<le> interpret_floatarith b xs" and u1: "interpret_floatarith b xs \<le> u2" by blast
2457 thus ?case unfolding l' u' by (auto simp add: real_of_float_min)
2460 from lift_bin'[OF Max.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Max.hyps
2461 obtain l1 u1 l2 u2 where l': "l = max l1 l2" and u': "u = max u1 u2"
2462 and l1: "l1 \<le> interpret_floatarith a xs" and u1: "interpret_floatarith a xs \<le> u1"
2463 and l1: "l2 \<le> interpret_floatarith b xs" and u1: "interpret_floatarith b xs \<le> u2" by blast
2464 thus ?case unfolding l' u' by (auto simp add: real_of_float_max)
2465 next case (Cos a) with lift_un'_bnds[OF bnds_cos] show ?case by auto
2466 next case (Arctan a) with lift_un'_bnds[OF bnds_arctan] show ?case by auto
2467 next case Pi with pi_boundaries show ?case by auto
2468 next case (Sqrt a) with lift_un'_bnds[OF bnds_sqrt] show ?case by auto
2469 next case (Exp a) with lift_un'_bnds[OF bnds_exp] show ?case by auto
2470 next case (Ln a) with lift_un_bnds[OF bnds_ln] show ?case by auto
2471 next case (Power a n) with lift_un'_bnds[OF bnds_power] show ?case by auto
2472 next case (Num f) thus ?case by auto
2475 from this[symmetric] `bounded_by xs vs`[THEN bounded_byE, of n]
2476 show ?case by (cases "n < length vs", auto)
2479 datatype form = Bound floatarith floatarith floatarith form
2480 | Assign floatarith floatarith form
2481 | Less floatarith floatarith
2482 | LessEqual floatarith floatarith
2483 | AtLeastAtMost floatarith floatarith floatarith
2485 fun interpret_form :: "form \<Rightarrow> real list \<Rightarrow> bool" where
2486 "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)" |
2487 "interpret_form (Assign x a f) vs = (interpret_floatarith x vs = interpret_floatarith a vs \<longrightarrow> interpret_form f vs)" |
2488 "interpret_form (Less a b) vs = (interpret_floatarith a vs < interpret_floatarith b vs)" |
2489 "interpret_form (LessEqual a b) vs = (interpret_floatarith a vs \<le> interpret_floatarith b vs)" |
2490 "interpret_form (AtLeastAtMost x a b) vs = (interpret_floatarith x vs \<in> { interpret_floatarith a vs .. interpret_floatarith b vs })"
2492 fun approx_form' and approx_form :: "nat \<Rightarrow> form \<Rightarrow> (float * float) option list \<Rightarrow> nat list \<Rightarrow> bool" where
2493 "approx_form' prec f 0 n l u bs ss = approx_form prec f (bs[n := Some (l, u)]) ss" |
2494 "approx_form' prec f (Suc s) n l u bs ss =
2495 (let m = (l + u) * Float 1 -1
2496 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))" |
2497 "approx_form prec (Bound (Var n) a b f) bs ss =
2498 (case (approx prec a bs, approx prec b bs)
2499 of (Some (l, _), Some (_, u)) \<Rightarrow> approx_form' prec f (ss ! n) n l u bs ss
2500 | _ \<Rightarrow> False)" |
2501 "approx_form prec (Assign (Var n) a f) bs ss =
2502 (case (approx prec a bs)
2503 of (Some (l, u)) \<Rightarrow> approx_form' prec f (ss ! n) n l u bs ss
2504 | _ \<Rightarrow> False)" |
2505 "approx_form prec (Less a b) bs ss =
2506 (case (approx prec a bs, approx prec b bs)
2507 of (Some (l, u), Some (l', u')) \<Rightarrow> u < l'
2508 | _ \<Rightarrow> False)" |
2509 "approx_form prec (LessEqual a b) bs ss =
2510 (case (approx prec a bs, approx prec b bs)
2511 of (Some (l, u), Some (l', u')) \<Rightarrow> u \<le> l'
2512 | _ \<Rightarrow> False)" |
2513 "approx_form prec (AtLeastAtMost x a b) bs ss =
2514 (case (approx prec x bs, approx prec a bs, approx prec b bs)
2515 of (Some (lx, ux), Some (l, u), Some (l', u')) \<Rightarrow> u \<le> lx \<and> ux \<le> l'
2516 | _ \<Rightarrow> False)" |
2517 "approx_form _ _ _ _ = False"
2519 lemma lazy_conj: "(if A then B else False) = (A \<and> B)" by simp
2521 lemma approx_form_approx_form':
2522 assumes "approx_form' prec f s n l u bs ss" and "(x::real) \<in> { l .. u }"
2523 obtains l' u' where "x \<in> { l' .. u' }"
2524 and "approx_form prec f (bs[n := Some (l', u')]) ss"
2525 using assms proof (induct s arbitrary: l u)
2527 from this(1)[of l u] this(2,3)
2532 let ?m = "(l + u) * Float 1 -1"
2533 have "real l \<le> ?m" and "?m \<le> real u"
2534 unfolding le_float_def using Suc.prems by auto
2536 with `x \<in> { l .. u }`
2537 have "x \<in> { l .. ?m} \<or> x \<in> { ?m .. u }" by auto
2540 assume *: "x \<in> { l .. ?m }"
2541 with Suc.hyps[OF _ _ *] Suc.prems
2542 show thesis by (simp add: Let_def lazy_conj)
2544 assume *: "x \<in> { ?m .. u }"
2545 with Suc.hyps[OF _ _ *] Suc.prems
2546 show thesis by (simp add: Let_def lazy_conj)
2550 lemma approx_form_aux:
2551 assumes "approx_form prec f vs ss"
2552 and "bounded_by xs vs"
2553 shows "interpret_form f xs"
2554 using assms proof (induct f arbitrary: vs)
2555 case (Bound x a b f)
2557 where x_eq: "x = Var n" by (cases x) auto
2559 with Bound.prems obtain l u' l' u
2560 where l_eq: "Some (l, u') = approx prec a vs"
2561 and u_eq: "Some (l', u) = approx prec b vs"
2562 and approx_form': "approx_form' prec f (ss ! n) n l u vs ss"
2563 by (cases "approx prec a vs", simp) (cases "approx prec b vs", auto)
2565 { assume "xs ! n \<in> { interpret_floatarith a xs .. interpret_floatarith b xs }"
2566 with approx[OF Bound.prems(2) l_eq] and approx[OF Bound.prems(2) u_eq]
2567 have "xs ! n \<in> { l .. u}" by auto
2569 from approx_form_approx_form'[OF approx_form' this]
2570 obtain lx ux where bnds: "xs ! n \<in> { lx .. ux }"
2571 and approx_form: "approx_form prec f (vs[n := Some (lx, ux)]) ss" .
2573 from `bounded_by xs vs` bnds
2574 have "bounded_by xs (vs[n := Some (lx, ux)])" by (rule bounded_by_update)
2575 with Bound.hyps[OF approx_form]
2576 have "interpret_form f xs" by blast }
2577 thus ?case using interpret_form.simps x_eq and interpret_floatarith.simps by simp
2581 where x_eq: "x = Var n" by (cases x) auto
2583 with Assign.prems obtain l u' l' u
2584 where bnd_eq: "Some (l, u) = approx prec a vs"
2585 and x_eq: "x = Var n"
2586 and approx_form': "approx_form' prec f (ss ! n) n l u vs ss"
2587 by (cases "approx prec a vs") auto
2589 { assume bnds: "xs ! n = interpret_floatarith a xs"
2590 with approx[OF Assign.prems(2) bnd_eq]
2591 have "xs ! n \<in> { l .. u}" by auto
2592 from approx_form_approx_form'[OF approx_form' this]
2593 obtain lx ux where bnds: "xs ! n \<in> { lx .. ux }"
2594 and approx_form: "approx_form prec f (vs[n := Some (lx, ux)]) ss" .
2596 from `bounded_by xs vs` bnds
2597 have "bounded_by xs (vs[n := Some (lx, ux)])" by (rule bounded_by_update)
2598 with Assign.hyps[OF approx_form]
2599 have "interpret_form f xs" by blast }
2600 thus ?case using interpret_form.simps x_eq and interpret_floatarith.simps by simp
2603 then obtain l u l' u'
2604 where l_eq: "Some (l, u) = approx prec a vs"
2605 and u_eq: "Some (l', u') = approx prec b vs"
2606 and inequality: "u < l'"
2607 by (cases "approx prec a vs", auto,
2608 cases "approx prec b vs", auto)
2609 from inequality[unfolded less_float_def] approx[OF Less.prems(2) l_eq] approx[OF Less.prems(2) u_eq]
2612 case (LessEqual a b)
2613 then obtain l u l' u'
2614 where l_eq: "Some (l, u) = approx prec a vs"
2615 and u_eq: "Some (l', u') = approx prec b vs"
2616 and inequality: "u \<le> l'"
2617 by (cases "approx prec a vs", auto,
2618 cases "approx prec b vs", auto)
2619 from inequality[unfolded le_float_def] approx[OF LessEqual.prems(2) l_eq] approx[OF LessEqual.prems(2) u_eq]
2622 case (AtLeastAtMost x a b)
2623 then obtain lx ux l u l' u'
2624 where x_eq: "Some (lx, ux) = approx prec x vs"
2625 and l_eq: "Some (l, u) = approx prec a vs"
2626 and u_eq: "Some (l', u') = approx prec b vs"
2627 and inequality: "u \<le> lx \<and> ux \<le> l'"
2628 by (cases "approx prec x vs", auto,
2629 cases "approx prec a vs", auto,
2630 cases "approx prec b vs", auto, blast)
2631 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]
2636 assumes "n = length xs"
2637 assumes "approx_form prec f (replicate n None) ss"
2638 shows "interpret_form f xs"
2639 using approx_form_aux[OF _ bounded_by_None] assms by auto
2641 subsection {* Implementing Taylor series expansion *}
2643 fun isDERIV :: "nat \<Rightarrow> floatarith \<Rightarrow> real list \<Rightarrow> bool" where
2644 "isDERIV x (Add a b) vs = (isDERIV x a vs \<and> isDERIV x b vs)" |
2645 "isDERIV x (Mult a b) vs = (isDERIV x a vs \<and> isDERIV x b vs)" |
2646 "isDERIV x (Minus a) vs = isDERIV x a vs" |
2647 "isDERIV x (Inverse a) vs = (isDERIV x a vs \<and> interpret_floatarith a vs \<noteq> 0)" |
2648 "isDERIV x (Cos a) vs = isDERIV x a vs" |
2649 "isDERIV x (Arctan a) vs = isDERIV x a vs" |
2650 "isDERIV x (Min a b) vs = False" |
2651 "isDERIV x (Max a b) vs = False" |
2652 "isDERIV x (Abs a) vs = False" |
2653 "isDERIV x Pi vs = True" |
2654 "isDERIV x (Sqrt a) vs = (isDERIV x a vs \<and> interpret_floatarith a vs > 0)" |
2655 "isDERIV x (Exp a) vs = isDERIV x a vs" |
2656 "isDERIV x (Ln a) vs = (isDERIV x a vs \<and> interpret_floatarith a vs > 0)" |
2657 "isDERIV x (Power a 0) vs = True" |
2658 "isDERIV x (Power a (Suc n)) vs = isDERIV x a vs" |
2659 "isDERIV x (Num f) vs = True" |
2660 "isDERIV x (Var n) vs = True"
2662 fun DERIV_floatarith :: "nat \<Rightarrow> floatarith \<Rightarrow> floatarith" where
2663 "DERIV_floatarith x (Add a b) = Add (DERIV_floatarith x a) (DERIV_floatarith x b)" |
2664 "DERIV_floatarith x (Mult a b) = Add (Mult a (DERIV_floatarith x b)) (Mult (DERIV_floatarith x a) b)" |
2665 "DERIV_floatarith x (Minus a) = Minus (DERIV_floatarith x a)" |
2666 "DERIV_floatarith x (Inverse a) = Minus (Mult (DERIV_floatarith x a) (Inverse (Power a 2)))" |
2667 "DERIV_floatarith x (Cos a) = Minus (Mult (Cos (Add (Mult Pi (Num (Float 1 -1))) (Minus a))) (DERIV_floatarith x a))" |
2668 "DERIV_floatarith x (Arctan a) = Mult (Inverse (Add (Num 1) (Power a 2))) (DERIV_floatarith x a)" |
2669 "DERIV_floatarith x (Min a b) = Num 0" |
2670 "DERIV_floatarith x (Max a b) = Num 0" |
2671 "DERIV_floatarith x (Abs a) = Num 0" |
2672 "DERIV_floatarith x Pi = Num 0" |
2673 "DERIV_floatarith x (Sqrt a) = (Mult (Inverse (Mult (Sqrt a) (Num 2))) (DERIV_floatarith x a))" |
2674 "DERIV_floatarith x (Exp a) = Mult (Exp a) (DERIV_floatarith x a)" |
2675 "DERIV_floatarith x (Ln a) = Mult (Inverse a) (DERIV_floatarith x a)" |
2676 "DERIV_floatarith x (Power a 0) = Num 0" |
2677 "DERIV_floatarith x (Power a (Suc n)) = Mult (Num (Float (int (Suc n)) 0)) (Mult (Power a n) (DERIV_floatarith x a))" |
2678 "DERIV_floatarith x (Num f) = Num 0" |
2679 "DERIV_floatarith x (Var n) = (if x = n then Num 1 else Num 0)"
2681 lemma DERIV_floatarith:
2682 assumes "n < length vs"
2683 assumes isDERIV: "isDERIV n f (vs[n := x])"
2684 shows "DERIV (\<lambda> x'. interpret_floatarith f (vs[n := x'])) x :>
2685 interpret_floatarith (DERIV_floatarith n f) (vs[n := x])"
2686 (is "DERIV (?i f) x :> _")
2687 using isDERIV proof (induct f arbitrary: x)
2688 case (Inverse a) thus ?case
2689 by (auto intro!: DERIV_intros
2690 simp add: algebra_simps power2_eq_square)
2691 next case (Cos a) thus ?case
2692 by (auto intro!: DERIV_intros
2693 simp del: interpret_floatarith.simps(5)
2694 simp add: interpret_floatarith_sin interpret_floatarith.simps(5)[of a])
2695 next case (Power a n) thus ?case
2696 by (cases n, auto intro!: DERIV_intros
2697 simp del: power_Suc simp add: real_eq_of_nat)
2698 next case (Ln a) thus ?case
2699 by (auto intro!: DERIV_intros simp add: divide_inverse)
2700 next case (Var i) thus ?case using `n < length vs` by auto
2701 qed (auto intro!: DERIV_intros)
2703 declare approx.simps[simp del]
2705 fun isDERIV_approx :: "nat \<Rightarrow> nat \<Rightarrow> floatarith \<Rightarrow> (float * float) option list \<Rightarrow> bool" where
2706 "isDERIV_approx prec x (Add a b) vs = (isDERIV_approx prec x a vs \<and> isDERIV_approx prec x b vs)" |
2707 "isDERIV_approx prec x (Mult a b) vs = (isDERIV_approx prec x a vs \<and> isDERIV_approx prec x b vs)" |
2708 "isDERIV_approx prec x (Minus a) vs = isDERIV_approx prec x a vs" |
2709 "isDERIV_approx prec x (Inverse a) vs =
2710 (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))" |
2711 "isDERIV_approx prec x (Cos a) vs = isDERIV_approx prec x a vs" |
2712 "isDERIV_approx prec x (Arctan a) vs = isDERIV_approx prec x a vs" |
2713 "isDERIV_approx prec x (Min a b) vs = False" |
2714 "isDERIV_approx prec x (Max a b) vs = False" |
2715 "isDERIV_approx prec x (Abs a) vs = False" |
2716 "isDERIV_approx prec x Pi vs = True" |
2717 "isDERIV_approx prec x (Sqrt a) vs =
2718 (isDERIV_approx prec x a vs \<and> (case approx prec a vs of Some (l, u) \<Rightarrow> 0 < l | None \<Rightarrow> False))" |
2719 "isDERIV_approx prec x (Exp a) vs = isDERIV_approx prec x a vs" |
2720 "isDERIV_approx prec x (Ln a) vs =
2721 (isDERIV_approx prec x a vs \<and> (case approx prec a vs of Some (l, u) \<Rightarrow> 0 < l | None \<Rightarrow> False))" |
2722 "isDERIV_approx prec x (Power a 0) vs = True" |
2723 "isDERIV_approx prec x (Power a (Suc n)) vs = isDERIV_approx prec x a vs" |
2724 "isDERIV_approx prec x (Num f) vs = True" |
2725 "isDERIV_approx prec x (Var n) vs = True"
2727 lemma isDERIV_approx:
2728 assumes "bounded_by xs vs"
2729 and isDERIV_approx: "isDERIV_approx prec x f vs"
2730 shows "isDERIV x f xs"
2731 using isDERIV_approx proof (induct f)
2733 then obtain l u where approx_Some: "Some (l, u) = approx prec a vs"
2734 and *: "0 < l \<or> u < 0"
2735 by (cases "approx prec a vs", auto)
2736 with approx[OF `bounded_by xs vs` approx_Some]
2737 have "interpret_floatarith a xs \<noteq> 0" unfolding less_float_def by auto
2738 thus ?case using Inverse by auto
2741 then obtain l u where approx_Some: "Some (l, u) = approx prec a vs"
2743 by (cases "approx prec a vs", auto)
2744 with approx[OF `bounded_by xs vs` approx_Some]
2745 have "0 < interpret_floatarith a xs" unfolding less_float_def by auto
2746 thus ?case using Ln by auto
2749 then obtain l u where approx_Some: "Some (l, u) = approx prec a vs"
2751 by (cases "approx prec a vs", auto)
2752 with approx[OF `bounded_by xs vs` approx_Some]
2753 have "0 < interpret_floatarith a xs" unfolding less_float_def by auto
2754 thus ?case using Sqrt by auto
2756 case (Power a n) thus ?case by (cases n, auto)
2759 lemma bounded_by_update_var:
2760 assumes "bounded_by xs vs" and "vs ! i = Some (l, u)"
2761 and bnd: "x \<in> { real l .. real u }"
2762 shows "bounded_by (xs[i := x]) vs"
2763 proof (cases "i < length xs")
2764 case False thus ?thesis using `bounded_by xs vs` by auto
2766 let ?xs = "xs[i := x]"
2767 case True hence "i < length ?xs" by auto
2769 assume "j < length vs"
2770 have "case vs ! j of None \<Rightarrow> True | Some (l, u) \<Rightarrow> ?xs ! j \<in> { real l .. real u }"
2771 proof (cases "vs ! j")
2774 proof (cases "i = j")
2776 thus ?thesis using `vs ! i = Some (l, u)` Some and bnd `i < length ?xs`
2780 thus ?thesis using `bounded_by xs vs`[THEN bounded_byE, OF `j < length vs`] Some
2784 thus ?thesis unfolding bounded_by_def by auto
2787 lemma isDERIV_approx':
2788 assumes "bounded_by xs vs"
2789 and vs_x: "vs ! x = Some (l, u)" and X_in: "X \<in> { real l .. real u }"
2790 and approx: "isDERIV_approx prec x f vs"
2791 shows "isDERIV x f (xs[x := X])"
2793 note bounded_by_update_var[OF `bounded_by xs vs` vs_x X_in] approx
2794 thus ?thesis by (rule isDERIV_approx)
2798 assumes "n < length xs" and bnd: "bounded_by xs vs"
2799 and isD: "isDERIV_approx prec n f vs"
2800 and app: "Some (l, u) = approx prec (DERIV_floatarith n f) vs" (is "_ = approx _ ?D _")
2801 shows "\<exists>(x::real). l \<le> x \<and> x \<le> u \<and>
2802 DERIV (\<lambda> x. interpret_floatarith f (xs[n := x])) (xs!n) :> x"
2803 (is "\<exists> x. _ \<and> _ \<and> DERIV (?i f) _ :> _")
2804 proof (rule exI[of _ "?i ?D (xs!n)"], rule conjI[OF _ conjI])
2805 let "?i f x" = "interpret_floatarith f (xs[n := x])"
2806 from approx[OF bnd app]
2807 show "l \<le> ?i ?D (xs!n)" and "?i ?D (xs!n) \<le> u"
2808 using `n < length xs` by auto
2809 from DERIV_floatarith[OF `n < length xs`, of f "xs!n"] isDERIV_approx[OF bnd isD]
2810 show "DERIV (?i f) (xs!n) :> (?i ?D (xs!n))" by simp
2813 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
2814 "lift_bin (Some (l1, u1)) (Some (l2, u2)) f = f l1 u1 l2 u2" |
2815 "lift_bin a b f = None"
2818 assumes lift_bin_Some: "Some (l, u) = lift_bin a b f"
2820 where "a = Some (l1, u1)"
2821 and "b = Some (l2, u2)"
2822 and "f l1 u1 l2 u2 = Some (l, u)"
2823 using assms by (cases a, simp, cases b, simp, auto)
2825 fun approx_tse where
2826 "approx_tse prec n 0 c k f bs = approx prec f bs" |
2827 "approx_tse prec n (Suc s) c k f bs =
2828 (if isDERIV_approx prec n f bs then
2829 lift_bin (approx prec f (bs[n := Some (c,c)]))
2830 (approx_tse prec n s c (Suc k) (DERIV_floatarith n f) bs)
2831 (\<lambda> l1 u1 l2 u2. approx prec
2833 (Mult (Inverse (Num (Float (int k) 0)))
2834 (Mult (Add (Var (Suc (Suc 0))) (Minus (Num c)))
2835 (Var (Suc 0))))) [Some (l1, u1), Some (l2, u2), bs!n])
2836 else approx prec f bs)"
2838 lemma bounded_by_Cons:
2839 assumes bnd: "bounded_by xs vs"
2840 and x: "x \<in> { real l .. real u }"
2841 shows "bounded_by (x#xs) ((Some (l, u))#vs)"
2843 { fix i assume *: "i < length ((Some (l, u))#vs)"
2844 have "case ((Some (l,u))#vs) ! i of Some (l, u) \<Rightarrow> (x#xs)!i \<in> { real l .. real u } | None \<Rightarrow> True"
2846 case 0 with x show ?thesis by auto
2848 case (Suc i) with * have "i < length vs" by auto
2849 from bnd[THEN bounded_byE, OF this]
2850 show ?thesis unfolding Suc nth_Cons_Suc .
2852 thus ?thesis by (auto simp add: bounded_by_def)
2855 lemma approx_tse_generic:
2856 assumes "bounded_by xs vs"
2857 and bnd_c: "bounded_by (xs[x := c]) vs" and "x < length vs" and "x < length xs"
2858 and bnd_x: "vs ! x = Some (lx, ux)"
2859 and ate: "Some (l, u) = approx_tse prec x s c k f vs"
2860 shows "\<exists> n. (\<forall> m < n. \<forall> (z::real) \<in> {lx .. ux}.
2861 DERIV (\<lambda> y. interpret_floatarith ((DERIV_floatarith x ^^ m) f) (xs[x := y])) z :>
2862 (interpret_floatarith ((DERIV_floatarith x ^^ (Suc m)) f) (xs[x := z])))
2863 \<and> (\<forall> (t::real) \<in> {lx .. ux}. (\<Sum> i = 0..<n. inverse (real (\<Prod> j \<in> {k..<k+i}. j)) *
2864 interpret_floatarith ((DERIV_floatarith x ^^ i) f) (xs[x := c]) *
2866 inverse (real (\<Prod> j \<in> {k..<k+n}. j)) *
2867 interpret_floatarith ((DERIV_floatarith x ^^ n) f) (xs[x := t]) *
2868 (xs!x - c)^n \<in> {l .. u})" (is "\<exists> n. ?taylor f k l u n")
2869 using ate proof (induct s arbitrary: k f l u)
2871 { fix t::real assume "t \<in> {lx .. ux}"
2872 note bounded_by_update_var[OF `bounded_by xs vs` bnd_x this]
2873 from approx[OF this 0[unfolded approx_tse.simps]]
2874 have "(interpret_floatarith f (xs[x := t])) \<in> {l .. u}"
2875 by (auto simp add: algebra_simps)
2876 } thus ?case by (auto intro!: exI[of _ 0])
2880 proof (cases "isDERIV_approx prec x f vs")
2882 note ap = Suc.prems[unfolded approx_tse.simps if_not_P[OF False]]
2884 { fix t::real assume "t \<in> {lx .. ux}"
2885 note bounded_by_update_var[OF `bounded_by xs vs` bnd_x this]
2886 from approx[OF this ap]
2887 have "(interpret_floatarith f (xs[x := t])) \<in> {l .. u}"
2888 by (auto simp add: algebra_simps)
2889 } thus ?thesis by (auto intro!: exI[of _ 0])
2894 where a: "Some (l1, u1) = approx prec f (vs[x := Some (c,c)])"
2895 and ate: "Some (l2, u2) = approx_tse prec x s c (Suc k) (DERIV_floatarith x f) vs"
2896 and final: "Some (l, u) = approx prec
2898 (Mult (Inverse (Num (Float (int k) 0)))
2899 (Mult (Add (Var (Suc (Suc 0))) (Minus (Num c)))
2900 (Var (Suc 0))))) [Some (l1, u1), Some (l2, u2), vs!x]"
2901 by (auto elim!: lift_bin) blast
2903 from bnd_c `x < length xs`
2904 have bnd: "bounded_by (xs[x:=c]) (vs[x:= Some (c,c)])"
2905 by (auto intro!: bounded_by_update)
2907 from approx[OF this a]
2908 have f_c: "interpret_floatarith ((DERIV_floatarith x ^^ 0) f) (xs[x := c]) \<in> { l1 .. u1 }"
2909 (is "?f 0 (real c) \<in> _")
2912 { fix f :: "'a \<Rightarrow> 'a" fix n :: nat fix x :: 'a
2913 have "(f ^^ Suc n) x = (f ^^ n) (f x)"
2914 by (induct n, auto) }
2915 note funpow_Suc = this[symmetric]
2916 from Suc.hyps[OF ate, unfolded this]
2918 where DERIV_hyp: "\<And> m z. \<lbrakk> m < n ; (z::real) \<in> { lx .. ux } \<rbrakk> \<Longrightarrow> DERIV (?f (Suc m)) z :> ?f (Suc (Suc m)) z"
2919 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) c * (xs!x - c)^i) +
2920 inverse (real (\<Prod> j \<in> {Suc k..<Suc k + n}. j)) * ?f (Suc n) t * (xs!x - c)^n \<in> {l2 .. u2}"
2921 (is "\<forall> t \<in> _. ?X (Suc k) f n t \<in> _")
2925 assume "m < Suc n" and bnd_z: "z \<in> { lx .. ux }"
2926 have "DERIV (?f m) z :> ?f (Suc m) z"
2929 with DERIV_floatarith[OF `x < length xs` isDERIV_approx'[OF `bounded_by xs vs` bnd_x bnd_z True]]
2930 show ?thesis by simp
2933 hence "m' < n" using `m < Suc n` by auto
2934 from DERIV_hyp[OF this bnd_z]
2935 show ?thesis using Suc by simp
2936 qed } note DERIV = this
2938 have "\<And> k i. k < i \<Longrightarrow> {k ..< i} = insert k {Suc k ..< i}" by auto
2939 hence setprod_head_Suc: "\<And> k i. \<Prod> {k ..< k + Suc i} = k * \<Prod> {Suc k ..< Suc k + i}" by auto
2940 have setsum_move0: "\<And> k F. setsum F {0..<Suc k} = F 0 + setsum (\<lambda> k. F (Suc k)) {0..<k}"
2941 unfolding setsum_shift_bounds_Suc_ivl[symmetric]
2942 unfolding setsum_head_upt_Suc[OF zero_less_Suc] ..
2943 def C \<equiv> "xs!x - c"
2945 { fix t::real assume t: "t \<in> {lx .. ux}"
2946 hence "bounded_by [xs!x] [vs!x]"
2947 using `bounded_by xs vs`[THEN bounded_byE, OF `x < length vs`]
2948 by (cases "vs!x", auto simp add: bounded_by_def)
2950 with hyp[THEN bspec, OF t] f_c
2951 have "bounded_by [?f 0 c, ?X (Suc k) f n t, xs!x] [Some (l1, u1), Some (l2, u2), vs!x]"
2952 by (auto intro!: bounded_by_Cons)
2953 from approx[OF this final, unfolded atLeastAtMost_iff[symmetric]]
2954 have "?X (Suc k) f n t * (xs!x - real c) * inverse k + ?f 0 c \<in> {l .. u}"
2955 by (auto simp add: algebra_simps)
2956 also have "?X (Suc k) f n t * (xs!x - real c) * inverse (real k) + ?f 0 c =
2957 (\<Sum> i = 0..<Suc n. inverse (real (\<Prod> j \<in> {k..<k+i}. j)) * ?f i c * (xs!x - c)^i) +
2958 inverse (real (\<Prod> j \<in> {k..<k+Suc n}. j)) * ?f (Suc n) t * (xs!x - c)^Suc n" (is "_ = ?T")
2959 unfolding funpow_Suc C_def[symmetric] setsum_move0 setprod_head_Suc
2960 by (auto simp add: algebra_simps)
2961 (simp only: mult_left_commute [of _ "inverse (real k)"] setsum_right_distrib [symmetric])
2962 finally have "?T \<in> {l .. u}" . }
2963 thus ?thesis using DERIV by blast
2967 lemma setprod_fact: "\<Prod> {1..<1 + k} = fact (k :: nat)"
2970 have "{ 1 ..< Suc (Suc k) } = insert (Suc k) { 1 ..< Suc k }" by auto
2971 hence "\<Prod> { 1 ..< Suc (Suc k) } = (Suc k) * \<Prod> { 1 ..< Suc k }" by auto
2972 thus ?case using Suc by auto
2976 assumes "bounded_by xs vs"
2977 and bnd_x: "vs ! x = Some (lx, ux)" and bnd_c: "real c \<in> {lx .. ux}"
2978 and "x < length vs" and "x < length xs"
2979 and ate: "Some (l, u) = approx_tse prec x s c 1 f vs"
2980 shows "interpret_floatarith f xs \<in> { l .. u }"
2982 def F \<equiv> "\<lambda> n z. interpret_floatarith ((DERIV_floatarith x ^^ n) f) (xs[x := z])"
2983 hence F0: "F 0 = (\<lambda> z. interpret_floatarith f (xs[x := z]))" by auto
2985 hence "bounded_by (xs[x := c]) vs" and "x < length vs" "x < length xs"
2986 using `bounded_by xs vs` bnd_x bnd_c `x < length vs` `x < length xs`
2987 by (auto intro!: bounded_by_update_var)
2989 from approx_tse_generic[OF `bounded_by xs vs` this bnd_x ate]
2991 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"
2992 and hyp: "\<And> (t::real). t \<in> {lx .. ux} \<Longrightarrow>
2993 (\<Sum> j = 0..<n. inverse (real (fact j)) * F j c * (xs!x - c)^j) +
2994 inverse (real (fact n)) * F n t * (xs!x - c)^n
2995 \<in> {l .. u}" (is "\<And> t. _ \<Longrightarrow> ?taylor t \<in> _")
2996 unfolding F_def atLeastAtMost_iff[symmetric] setprod_fact by blast
2998 have bnd_xs: "xs ! x \<in> { lx .. ux }"
2999 using `bounded_by xs vs`[THEN bounded_byE, OF `x < length vs`] bnd_x by auto
3003 case 0 thus ?thesis using hyp[OF bnd_xs] unfolding F_def by auto
3007 proof (cases "xs ! x = c")
3009 from True[symmetric] hyp[OF bnd_xs] Suc show ?thesis
3010 unfolding F_def Suc setsum_head_upt_Suc[OF zero_less_Suc] setsum_shift_bounds_Suc_ivl by auto
3014 have "lx \<le> real c" "real c \<le> ux" "lx \<le> xs!x" "xs!x \<le> ux"
3015 using Suc bnd_c `bounded_by xs vs`[THEN bounded_byE, OF `x < length vs`] bnd_x by auto
3016 from Taylor.taylor[OF zero_less_Suc, of F, OF F0 DERIV[unfolded Suc] this False]
3017 obtain t::real where t_bnd: "if xs ! x < c then xs ! x < t \<and> t < c else c < t \<and> t < xs ! x"
3018 and fl_eq: "interpret_floatarith f (xs[x := xs ! x]) =
3019 (\<Sum>m = 0..<Suc n'. F m c / real (fact m) * (xs ! x - c) ^ m) +
3020 F (Suc n') t / real (fact (Suc n')) * (xs ! x - c) ^ Suc n'"
3023 from t_bnd bnd_xs bnd_c have *: "t \<in> {lx .. ux}"
3024 by (cases "xs ! x < c", auto)
3026 have "interpret_floatarith f (xs[x := xs ! x]) = ?taylor t"
3027 unfolding fl_eq Suc by (auto simp add: algebra_simps divide_inverse)
3028 also have "\<dots> \<in> {l .. u}" using * by (rule hyp)
3029 finally show ?thesis by simp
3034 fun approx_tse_form' where
3035 "approx_tse_form' prec t f 0 l u cmp =
3036 (case approx_tse prec 0 t ((l + u) * Float 1 -1) 1 f [Some (l, u)]
3037 of Some (l, u) \<Rightarrow> cmp l u | None \<Rightarrow> False)" |
3038 "approx_tse_form' prec t f (Suc s) l u cmp =
3039 (let m = (l + u) * Float 1 -1
3040 in (if approx_tse_form' prec t f s l m cmp then
3041 approx_tse_form' prec t f s m u cmp else False))"
3043 lemma approx_tse_form':
3045 assumes "approx_tse_form' prec t f s l u cmp" and "x \<in> {l .. u}"
3046 shows "\<exists> l' u' ly uy. x \<in> { l' .. u' } \<and> real l \<le> l' \<and> u' \<le> real u \<and> cmp ly uy \<and>
3047 approx_tse prec 0 t ((l' + u') * Float 1 -1) 1 f [Some (l', u')] = Some (ly, uy)"
3048 using assms proof (induct s arbitrary: l u)
3051 where *: "approx_tse prec 0 t ((l + u) * Float 1 -1) 1 f [Some (l, u)] = Some (ly, uy)"
3052 and **: "cmp ly uy" by (auto elim!: option_caseE)
3053 with 0 show ?case by (auto intro!: exI)
3056 let ?m = "(l + u) * Float 1 -1"
3058 have l: "approx_tse_form' prec t f s l ?m cmp"
3059 and u: "approx_tse_form' prec t f s ?m u cmp"
3060 by (auto simp add: Let_def lazy_conj)
3062 have m_l: "real l \<le> ?m" and m_u: "?m \<le> real u"
3063 unfolding le_float_def using Suc.prems by auto
3065 with `x \<in> { l .. u }`
3066 have "x \<in> { l .. ?m} \<or> x \<in> { ?m .. u }" by auto
3069 assume "x \<in> { l .. ?m}"
3070 from Suc.hyps[OF l this]
3072 where "x \<in> { l' .. u' } \<and> real l \<le> l' \<and> real u' \<le> ?m \<and> cmp ly uy \<and>
3073 approx_tse prec 0 t ((l' + u') * Float 1 -1) 1 f [Some (l', u')] = Some (ly, uy)" by blast
3074 with m_u show ?thesis by (auto intro!: exI)
3076 assume "x \<in> { ?m .. u }"
3077 from Suc.hyps[OF u this]
3079 where "x \<in> { l' .. u' } \<and> ?m \<le> real l' \<and> u' \<le> real u \<and> cmp ly uy \<and>
3080 approx_tse prec 0 t ((l' + u') * Float 1 -1) 1 f [Some (l', u')] = Some (ly, uy)" by blast
3081 with m_u show ?thesis by (auto intro!: exI)
3085 lemma approx_tse_form'_less:
3087 assumes tse: "approx_tse_form' prec t (Add a (Minus b)) s l u (\<lambda> l u. 0 < l)"
3088 and x: "x \<in> {l .. u}"
3089 shows "interpret_floatarith b [x] < interpret_floatarith a [x]"
3091 from approx_tse_form'[OF tse x]
3093 where x': "x \<in> { l' .. u' }" and "l \<le> real l'"
3094 and "real u' \<le> u" and "0 < ly"
3095 and tse: "approx_tse prec 0 t ((l' + u') * Float 1 -1) 1 (Add a (Minus b)) [Some (l', u')] = Some (ly, uy)"
3098 hence "bounded_by [x] [Some (l', u')]" by (auto simp add: bounded_by_def)
3100 from approx_tse[OF this _ _ _ _ tse[symmetric], of l' u'] x'
3101 have "ly \<le> interpret_floatarith a [x] - interpret_floatarith b [x]"
3102 by (auto simp add: diff_minus)
3103 from order_less_le_trans[OF `0 < ly`[unfolded less_float_def] this]
3104 show ?thesis by auto
3107 lemma approx_tse_form'_le:
3109 assumes tse: "approx_tse_form' prec t (Add a (Minus b)) s l u (\<lambda> l u. 0 \<le> l)"
3110 and x: "x \<in> {l .. u}"
3111 shows "interpret_floatarith b [x] \<le> interpret_floatarith a [x]"
3113 from approx_tse_form'[OF tse x]
3115 where x': "x \<in> { l' .. u' }" and "l \<le> real l'"
3116 and "real u' \<le> u" and "0 \<le> ly"
3117 and tse: "approx_tse prec 0 t ((l' + u') * Float 1 -1) 1 (Add a (Minus b)) [Some (l', u')] = Some (ly, uy)"
3120 hence "bounded_by [x] [Some (l', u')]" by (auto simp add: bounded_by_def)
3122 from approx_tse[OF this _ _ _ _ tse[symmetric], of l' u'] x'
3123 have "ly \<le> interpret_floatarith a [x] - interpret_floatarith b [x]"
3124 by (auto simp add: diff_minus)
3125 from order_trans[OF `0 \<le> ly`[unfolded le_float_def] this]
3126 show ?thesis by auto
3130 "approx_tse_form prec t s f =
3132 of (Bound x a b f) \<Rightarrow> x = Var 0 \<and>
3133 (case (approx prec a [None], approx prec b [None])
3134 of (Some (l, u), Some (l', u')) \<Rightarrow>
3136 of Less lf rt \<Rightarrow> approx_tse_form' prec t (Add rt (Minus lf)) s l u' (\<lambda> l u. 0 < l)
3137 | LessEqual lf rt \<Rightarrow> approx_tse_form' prec t (Add rt (Minus lf)) s l u' (\<lambda> l u. 0 \<le> l)
3138 | AtLeastAtMost x lf rt \<Rightarrow>
3139 (if approx_tse_form' prec t (Add x (Minus lf)) s l u' (\<lambda> l u. 0 \<le> l) then
3140 approx_tse_form' prec t (Add rt (Minus x)) s l u' (\<lambda> l u. 0 \<le> l) else False)
3141 | _ \<Rightarrow> False)
3142 | _ \<Rightarrow> False)
3143 | _ \<Rightarrow> False)"
3145 lemma approx_tse_form:
3146 assumes "approx_tse_form prec t s f"
3147 shows "interpret_form f [x]"
3149 case (Bound i a b f') note f_def = this
3150 with assms obtain l u l' u'
3151 where a: "approx prec a [None] = Some (l, u)"
3152 and b: "approx prec b [None] = Some (l', u')"
3153 unfolding approx_tse_form_def by (auto elim!: option_caseE)
3155 from Bound assms have "i = Var 0" unfolding approx_tse_form_def by auto
3156 hence i: "interpret_floatarith i [x] = x" by auto
3158 { let "?f z" = "interpret_floatarith z [x]"
3159 assume "?f i \<in> { ?f a .. ?f b }"
3160 with approx[OF _ a[symmetric], of "[x]"] approx[OF _ b[symmetric], of "[x]"]
3161 have bnd: "x \<in> { l .. u'}" unfolding bounded_by_def i by auto
3163 have "interpret_form f' [x]"
3166 with Bound a b assms
3167 have "approx_tse_form' prec t (Add rt (Minus lf)) s l u' (\<lambda> l u. 0 < l)"
3168 unfolding approx_tse_form_def by auto
3169 from approx_tse_form'_less[OF this bnd]
3170 show ?thesis using Less by auto
3172 case (LessEqual lf rt)
3173 with Bound a b assms
3174 have "approx_tse_form' prec t (Add rt (Minus lf)) s l u' (\<lambda> l u. 0 \<le> l)"
3175 unfolding approx_tse_form_def by auto
3176 from approx_tse_form'_le[OF this bnd]
3177 show ?thesis using LessEqual by auto
3179 case (AtLeastAtMost x lf rt)
3180 with Bound a b assms
3181 have "approx_tse_form' prec t (Add rt (Minus x)) s l u' (\<lambda> l u. 0 \<le> l)"
3182 and "approx_tse_form' prec t (Add x (Minus lf)) s l u' (\<lambda> l u. 0 \<le> l)"
3183 unfolding approx_tse_form_def lazy_conj by auto
3184 from approx_tse_form'_le[OF this(1) bnd] approx_tse_form'_le[OF this(2) bnd]
3185 show ?thesis using AtLeastAtMost by auto
3187 case (Bound x a b f') with assms
3188 show ?thesis by (auto elim!: option_caseE simp add: f_def approx_tse_form_def)
3190 case (Assign x a f') with assms
3191 show ?thesis by (auto elim!: option_caseE simp add: f_def approx_tse_form_def)
3192 qed } thus ?thesis unfolding f_def by auto
3193 next case Assign with assms show ?thesis by (auto simp add: approx_tse_form_def)
3194 next case LessEqual with assms show ?thesis by (auto simp add: approx_tse_form_def)
3195 next case Less with assms show ?thesis by (auto simp add: approx_tse_form_def)
3196 next case AtLeastAtMost with assms show ?thesis by (auto simp add: approx_tse_form_def)
3199 text {* @{term approx_form_eval} is only used for the {\tt value}-command. *}
3201 fun approx_form_eval :: "nat \<Rightarrow> form \<Rightarrow> (float * float) option list \<Rightarrow> (float * float) option list" where
3202 "approx_form_eval prec (Bound (Var n) a b f) bs =
3203 (case (approx prec a bs, approx prec b bs)
3204 of (Some (l, _), Some (_, u)) \<Rightarrow> approx_form_eval prec f (bs[n := Some (l, u)])
3205 | _ \<Rightarrow> bs)" |
3206 "approx_form_eval prec (Assign (Var n) a f) bs =
3207 (case (approx prec a bs)
3208 of (Some (l, u)) \<Rightarrow> approx_form_eval prec f (bs[n := Some (l, u)])
3209 | _ \<Rightarrow> bs)" |
3210 "approx_form_eval prec (Less a b) bs = bs @ [approx prec a bs, approx prec b bs]" |
3211 "approx_form_eval prec (LessEqual a b) bs = bs @ [approx prec a bs, approx prec b bs]" |
3212 "approx_form_eval prec (AtLeastAtMost x a b) bs =
3213 bs @ [approx prec x bs, approx prec a bs, approx prec b bs]" |
3214 "approx_form_eval _ _ bs = bs"
3216 subsection {* Implement proof method \texttt{approximation} *}
3218 lemmas interpret_form_equations = interpret_form.simps interpret_floatarith.simps interpret_floatarith_num
3219 interpret_floatarith_divide interpret_floatarith_diff interpret_floatarith_tan interpret_floatarith_powr interpret_floatarith_log
3220 interpret_floatarith_sin
3222 oracle approximation_oracle = {* fn (thy, t) =>
3225 fun bad t = error ("Bad term: " ^ Syntax.string_of_term_global thy t);
3227 fun term_of_bool true = @{term True}
3228 | term_of_bool false = @{term False};
3230 fun term_of_float (@{code Float} (k, l)) =
3231 @{term Float} $ HOLogic.mk_number @{typ int} k $ HOLogic.mk_number @{typ int} l;
3233 fun term_of_float_float_option NONE = @{term "None :: (float \<times> float) option"}
3234 | term_of_float_float_option (SOME ff) = @{term "Some :: float \<times> float \<Rightarrow> _"}
3235 $ HOLogic.mk_prod (pairself term_of_float ff);
3237 val term_of_float_float_option_list =
3238 HOLogic.mk_list @{typ "(float \<times> float) option"} o map term_of_float_float_option;
3240 fun nat_of_term t = HOLogic.dest_nat t handle TERM _ => snd (HOLogic.dest_number t);
3242 fun float_of_term (@{term Float} $ k $ l) =
3243 @{code Float} (snd (HOLogic.dest_number k), snd (HOLogic.dest_number l))
3244 | float_of_term t = bad t;
3246 fun floatarith_of_term (@{term Add} $ a $ b) = @{code Add} (floatarith_of_term a, floatarith_of_term b)
3247 | floatarith_of_term (@{term Minus} $ a) = @{code Minus} (floatarith_of_term a)
3248 | floatarith_of_term (@{term Mult} $ a $ b) = @{code Mult} (floatarith_of_term a, floatarith_of_term b)
3249 | floatarith_of_term (@{term Inverse} $ a) = @{code Inverse} (floatarith_of_term a)
3250 | floatarith_of_term (@{term Cos} $ a) = @{code Cos} (floatarith_of_term a)
3251 | floatarith_of_term (@{term Arctan} $ a) = @{code Arctan} (floatarith_of_term a)
3252 | floatarith_of_term (@{term Abs} $ a) = @{code Abs} (floatarith_of_term a)
3253 | floatarith_of_term (@{term Max} $ a $ b) = @{code Max} (floatarith_of_term a, floatarith_of_term b)
3254 | floatarith_of_term (@{term Min} $ a $ b) = @{code Min} (floatarith_of_term a, floatarith_of_term b)
3255 | floatarith_of_term @{term Pi} = @{code Pi}
3256 | floatarith_of_term (@{term Sqrt} $ a) = @{code Sqrt} (floatarith_of_term a)
3257 | floatarith_of_term (@{term Exp} $ a) = @{code Exp} (floatarith_of_term a)
3258 | floatarith_of_term (@{term Ln} $ a) = @{code Ln} (floatarith_of_term a)
3259 | floatarith_of_term (@{term Power} $ a $ n) =
3260 @{code Power} (floatarith_of_term a, nat_of_term n)
3261 | floatarith_of_term (@{term Var} $ n) = @{code Var} (nat_of_term n)
3262 | floatarith_of_term (@{term Num} $ m) = @{code Num} (float_of_term m)
3263 | floatarith_of_term t = bad t;
3265 fun form_of_term (@{term Bound} $ a $ b $ c $ p) = @{code Bound}
3266 (floatarith_of_term a, floatarith_of_term b, floatarith_of_term c, form_of_term p)
3267 | form_of_term (@{term Assign} $ a $ b $ p) = @{code Assign}
3268 (floatarith_of_term a, floatarith_of_term b, form_of_term p)
3269 | form_of_term (@{term Less} $ a $ b) = @{code Less}
3270 (floatarith_of_term a, floatarith_of_term b)
3271 | form_of_term (@{term LessEqual} $ a $ b) = @{code LessEqual}
3272 (floatarith_of_term a, floatarith_of_term b)
3273 | form_of_term (@{term AtLeastAtMost} $ a $ b $ c) = @{code AtLeastAtMost}
3274 (floatarith_of_term a, floatarith_of_term b, floatarith_of_term c)
3275 | form_of_term t = bad t;
3277 fun float_float_option_of_term @{term "None :: (float \<times> float) option"} = NONE
3278 | float_float_option_of_term (@{term "Some :: float \<times> float \<Rightarrow> _"} $ ff) =
3279 SOME (pairself float_of_term (HOLogic.dest_prod ff))
3280 | float_float_option_of_term (@{term approx'} $ n $ a $ ffs) = @{code approx'}
3281 (nat_of_term n) (floatarith_of_term a) (float_float_option_list_of_term ffs)
3282 | float_float_option_of_term t = bad t
3283 and float_float_option_list_of_term
3284 (@{term "replicate :: _ \<Rightarrow> (float \<times> float) option \<Rightarrow> _"} $ n $ @{term "None :: (float \<times> float) option"}) =
3285 @{code replicate} (nat_of_term n) NONE
3286 | float_float_option_list_of_term (@{term approx_form_eval} $ n $ p $ ffs) =
3287 @{code approx_form_eval} (nat_of_term n) (form_of_term p) (float_float_option_list_of_term ffs)
3288 | float_float_option_list_of_term t = map float_float_option_of_term
3289 (HOLogic.dest_list t);
3291 val nat_list_of_term = map nat_of_term o HOLogic.dest_list ;
3293 fun bool_of_term (@{term approx_form} $ n $ p $ ffs $ ms) = @{code approx_form}
3294 (nat_of_term n) (form_of_term p) (float_float_option_list_of_term ffs) (nat_list_of_term ms)
3295 | bool_of_term (@{term approx_tse_form} $ m $ n $ q $ p) =
3296 @{code approx_tse_form} (nat_of_term m) (nat_of_term n) (nat_of_term q) (form_of_term p)
3297 | bool_of_term t = bad t;
3299 fun eval t = case fastype_of t
3301 (term_of_bool o bool_of_term) t
3302 | @{typ "(float \<times> float) option"} =>
3303 (term_of_float_float_option o float_float_option_of_term) t
3304 | @{typ "(float \<times> float) option list"} =>
3305 (term_of_float_float_option_list o float_float_option_list_of_term) t
3308 val normalize = eval o Envir.beta_norm o Pattern.eta_long [];
3310 in Thm.cterm_of thy (Logic.mk_equals (t, normalize t)) end
3314 fun reorder_bounds_tac prems i =
3316 fun variable_of_bound (Const (@{const_name Trueprop}, _) $
3317 (Const (@{const_name Set.member}, _) $
3318 Free (name, _) $ _)) = name
3319 | variable_of_bound (Const (@{const_name Trueprop}, _) $
3320 (Const (@{const_name HOL.eq}, _) $
3321 Free (name, _) $ _)) = name
3322 | variable_of_bound t = raise TERM ("variable_of_bound", [t])
3325 = map (` (variable_of_bound o prop_of)) prems
3327 fun add_deps (name, bnds)
3328 = Graph.add_deps_acyclic (name,
3329 remove (op =) name (Term.add_free_names (prop_of bnds) []))
3331 val order = Graph.empty
3332 |> fold Graph.new_node variable_bounds
3333 |> fold add_deps variable_bounds
3334 |> Graph.strong_conn |> map the_single |> rev
3335 |> map_filter (AList.lookup (op =) variable_bounds)
3337 fun prepend_prem th tac
3338 = tac THEN rtac (th RSN (2, @{thm mp})) i
3340 fold prepend_prem order all_tac
3343 fun approximation_conv ctxt ct =
3344 approximation_oracle (ProofContext.theory_of ctxt, Thm.term_of ct |> tap (tracing o Syntax.string_of_term ctxt));
3346 fun approximate ctxt t =
3347 approximation_oracle (ProofContext.theory_of ctxt, t)
3348 |> Thm.prop_of |> Logic.dest_equals |> snd;
3350 (* Should be in HOL.thy ? *)
3351 fun gen_eval_tac conv ctxt = CONVERSION
3352 (Object_Logic.judgment_conv (Conv.params_conv (~1) (K (Conv.concl_conv (~1) conv)) ctxt))
3355 val form_equations = @{thms interpret_form_equations};
3357 fun rewrite_interpret_form_tac ctxt prec splitting taylor i st = let
3358 fun lookup_splitting (Free (name, typ))
3359 = case AList.lookup (op =) splitting name
3360 of SOME s => HOLogic.mk_number @{typ nat} s
3361 | NONE => @{term "0 :: nat"}
3362 val vs = nth (prems_of st) (i - 1)
3363 |> Logic.strip_imp_concl
3364 |> HOLogic.dest_Trueprop
3365 |> Term.strip_comb |> snd |> List.last
3366 |> HOLogic.dest_list
3368 |> HOLogic.mk_number @{typ nat}
3369 |> Thm.cterm_of (ProofContext.theory_of ctxt)
3372 val n = vs |> length
3373 |> HOLogic.mk_number @{typ nat}
3374 |> Thm.cterm_of (ProofContext.theory_of ctxt)
3376 |> map lookup_splitting
3377 |> HOLogic.mk_list @{typ nat}
3378 |> Thm.cterm_of (ProofContext.theory_of ctxt)
3380 (rtac (Thm.instantiate ([], [(@{cpat "?n::nat"}, n),
3381 (@{cpat "?prec::nat"}, p),
3382 (@{cpat "?ss::nat list"}, s)])
3383 @{thm "approx_form"}) i
3384 THEN simp_tac @{simpset} i) st
3387 | SOME t => if length vs <> 1 then raise (TERM ("More than one variable used for taylor series expansion", [prop_of st]))
3390 |> HOLogic.mk_number @{typ nat}
3391 |> Thm.cterm_of (ProofContext.theory_of ctxt)
3392 val s = vs |> map lookup_splitting |> hd
3393 |> Thm.cterm_of (ProofContext.theory_of ctxt)
3395 rtac (Thm.instantiate ([], [(@{cpat "?s::nat"}, s),
3396 (@{cpat "?t::nat"}, t),
3397 (@{cpat "?prec::nat"}, p)])
3398 @{thm "approx_tse_form"}) i st
3402 (* copied from Tools/induct.ML should probably in args.ML *)
3403 val free = Args.context -- Args.term >> (fn (_, Free (n, t)) => n | (ctxt, t) =>
3404 error ("Bad free variable: " ^ Syntax.string_of_term ctxt t));
3408 lemma intervalE: "a \<le> x \<and> x \<le> b \<Longrightarrow> \<lbrakk> x \<in> { a .. b } \<Longrightarrow> P\<rbrakk> \<Longrightarrow> P"
3411 lemma meta_eqE: "x \<equiv> a \<Longrightarrow> \<lbrakk> x = a \<Longrightarrow> P\<rbrakk> \<Longrightarrow> P"
3414 method_setup approximation = {*
3417 Scan.optional (Scan.lift (Args.$$$ "splitting" |-- Args.colon)
3418 |-- Parse.and_list' (free --| Scan.lift (Args.$$$ "=") -- Scan.lift Parse.nat)) []
3420 Scan.option (Scan.lift (Args.$$$ "taylor" |-- Args.colon)
3421 |-- (free |-- Scan.lift (Args.$$$ "=") |-- Scan.lift Parse.nat))
3423 (fn ((prec, splitting), taylor) => fn ctxt =>
3424 SIMPLE_METHOD' (fn i =>
3425 REPEAT (FIRST' [etac @{thm intervalE},
3426 etac @{thm meta_eqE},
3427 rtac @{thm impI}] i)
3428 THEN Subgoal.FOCUS (fn {prems, ...} => reorder_bounds_tac prems i) @{context} i
3429 THEN DETERM (TRY (filter_prems_tac (K false) i))
3430 THEN DETERM (Reflection.genreify_tac ctxt form_equations NONE i)
3431 THEN rewrite_interpret_form_tac ctxt prec splitting taylor i
3432 THEN gen_eval_tac (approximation_conv ctxt) ctxt i))
3433 *} "real number approximation"
3436 fun calculated_subterms (@{const Trueprop} $ t) = calculated_subterms t
3437 | calculated_subterms (@{const HOL.implies} $ _ $ t) = calculated_subterms t
3438 | calculated_subterms (@{term "op <= :: real \<Rightarrow> real \<Rightarrow> bool"} $ t1 $ t2) = [t1, t2]
3439 | calculated_subterms (@{term "op < :: real \<Rightarrow> real \<Rightarrow> bool"} $ t1 $ t2) = [t1, t2]
3440 | calculated_subterms (@{term "op : :: real \<Rightarrow> real set \<Rightarrow> bool"} $ t1 $
3441 (@{term "atLeastAtMost :: real \<Rightarrow> real \<Rightarrow> real set"} $ t2 $ t3)) = [t1, t2, t3]
3442 | calculated_subterms t = raise TERM ("calculated_subterms", [t])
3444 fun dest_interpret_form (@{const "interpret_form"} $ b $ xs) = (b, xs)
3445 | dest_interpret_form t = raise TERM ("dest_interpret_form", [t])
3447 fun dest_interpret (@{const "interpret_floatarith"} $ b $ xs) = (b, xs)
3448 | dest_interpret t = raise TERM ("dest_interpret", [t])
3451 fun dest_float (@{const "Float"} $ m $ e) = (snd (HOLogic.dest_number m), snd (HOLogic.dest_number e))
3452 fun dest_ivl (Const (@{const_name "Some"}, _) $
3453 (Const (@{const_name Pair}, _) $ u $ l)) = SOME (dest_float u, dest_float l)
3454 | dest_ivl (Const (@{const_name "None"}, _)) = NONE
3455 | dest_ivl t = raise TERM ("dest_result", [t])
3457 fun mk_approx' prec t = (@{const "approx'"}
3458 $ HOLogic.mk_number @{typ nat} prec
3459 $ t $ @{term "[] :: (float * float) option list"})
3461 fun mk_approx_form_eval prec t xs = (@{const "approx_form_eval"}
3462 $ HOLogic.mk_number @{typ nat} prec
3465 fun float2_float10 prec round_down (m, e) = (
3467 val (m, e) = (if e < 0 then (m,e) else (m * Integer.pow e 2, 0))
3469 fun frac c p 0 digits cnt = (digits, cnt, 0)
3470 | frac c 0 r digits cnt = (digits, cnt, r)
3471 | frac c p r digits cnt = (let
3472 val (d, r) = Integer.div_mod (r * 10) (Integer.pow (~e) 2)
3473 in frac (c orelse d <> 0) (if d <> 0 orelse c then p - 1 else p) r
3474 (digits * 10 + d) (cnt + 1)
3477 val sgn = Int.sign m
3480 val round_down = (sgn = 1 andalso round_down) orelse
3481 (sgn = ~1 andalso not round_down)
3483 val (x, r) = Integer.div_mod m (Integer.pow (~e) 2)
3485 val p = ((if x = 0 then prec else prec - (IntInf.log2 x + 1)) * 3) div 10 + 1
3487 val (digits, e10, r) = if p > 0 then frac (x <> 0) p r 0 0 else (0,0,0)
3489 val digits = if round_down orelse r = 0 then digits else digits + 1
3491 in (sgn * (digits + x * (Integer.pow e10 10)), ~e10)
3494 fun mk_result prec (SOME (l, u)) = (let
3495 fun mk_float10 rnd x = (let val (m, e) = float2_float10 prec rnd x
3496 in if e = 0 then HOLogic.mk_number @{typ real} m
3497 else if e = 1 then @{term "divide :: real \<Rightarrow> real \<Rightarrow> real"} $
3498 HOLogic.mk_number @{typ real} m $
3500 else @{term "divide :: real \<Rightarrow> real \<Rightarrow> real"} $
3501 HOLogic.mk_number @{typ real} m $
3502 (@{term "power 10 :: nat \<Rightarrow> real"} $
3503 HOLogic.mk_number @{typ nat} (~e)) end)
3504 in @{term "atLeastAtMost :: real \<Rightarrow> real \<Rightarrow> real set"} $ mk_float10 true l $ mk_float10 false u end)
3505 | mk_result prec NONE = @{term "UNIV :: real set"}
3508 val t = Logic.varify_global t
3509 val m = map (fn (name, sort) => (name, @{typ real})) (Term.add_tvars t [])
3510 val t = Term.subst_TVars m t
3513 fun converted_result t =
3515 |> HOLogic.dest_Trueprop
3516 |> HOLogic.dest_eq |> snd
3518 fun apply_tactic context term tactic = cterm_of context term
3521 |> the |> prems_of |> hd
3523 fun prepare_form context term = apply_tactic context term (
3524 REPEAT (FIRST' [etac @{thm intervalE}, etac @{thm meta_eqE}, rtac @{thm impI}] 1)
3525 THEN Subgoal.FOCUS (fn {prems, ...} => reorder_bounds_tac prems 1) @{context} 1
3526 THEN DETERM (TRY (filter_prems_tac (K false) 1)))
3528 fun reify_form context term = apply_tactic context term
3529 (Reflection.genreify_tac @{context} form_equations NONE 1)
3531 fun approx_form prec ctxt t =
3533 |> prepare_form (ProofContext.theory_of ctxt)
3534 |> (fn arith_term =>
3535 reify_form (ProofContext.theory_of ctxt) arith_term
3536 |> HOLogic.dest_Trueprop |> dest_interpret_form
3537 |> (fn (data, xs) =>
3538 mk_approx_form_eval prec data (HOLogic.mk_list @{typ "(float * float) option"}
3539 (map (fn _ => @{term "None :: (float * float) option"}) (HOLogic.dest_list xs)))
3541 |> HOLogic.dest_list
3542 |> curry ListPair.zip (HOLogic.dest_list xs @ calculated_subterms arith_term)
3543 |> map (fn (elem, s) => @{term "op : :: real \<Rightarrow> real set \<Rightarrow> bool"} $ elem $ mk_result prec (dest_ivl s))
3544 |> foldr1 HOLogic.mk_conj))
3546 fun approx_arith prec ctxt t = realify t
3547 |> Reflection.genreif ctxt form_equations
3549 |> HOLogic.dest_Trueprop
3550 |> HOLogic.dest_eq |> snd
3551 |> dest_interpret |> fst
3557 fun approx prec ctxt t = if type_of t = @{typ prop} then approx_form prec ctxt t
3558 else if type_of t = @{typ bool} then approx_form prec ctxt (@{const Trueprop} $ t)
3559 else approx_arith prec ctxt t
3563 Value.add_evaluator ("approximate", approx 30)