1 (* Title: HOL/Decision_Procs/Approximation.thy
2 Author: Johannes Hoelzl <hoelzl@in.tum.de> 2008 / 2009
5 header {* Prove unequations about real numbers by computation *}
8 imports Complex_Main Float Reflection Dense_Linear_Order Efficient_Nat
11 section "Horner Scheme"
13 subsection {* Define auxiliary helper @{text horner} function *}
15 fun horner :: "(nat \<Rightarrow> nat) \<Rightarrow> (nat \<Rightarrow> nat \<Rightarrow> nat) \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> real \<Rightarrow> real" where
16 "horner F G 0 i k x = 0" |
17 "horner F G (Suc n) i k x = 1 / real k - x * horner F G n (F i) (G i k) x"
19 lemma horner_schema': fixes x :: real and a :: "nat \<Rightarrow> real"
20 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)"
22 have shift_pow: "\<And>i. - (x * ((-1)^i * a (Suc i) * x ^ i)) = (-1)^(Suc i) * a (Suc i) * x ^ (Suc i)" by auto
23 show ?thesis unfolding setsum_right_distrib shift_pow real_diff_def setsum_negf[symmetric] setsum_head_upt_Suc[OF zero_less_Suc]
24 setsum_reindex[OF inj_Suc, unfolded comp_def, symmetric, of "\<lambda> n. (-1)^n *a n * x^n"] by auto
27 lemma horner_schema: fixes f :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat" and F :: "nat \<Rightarrow> nat"
28 assumes f_Suc: "\<And>n. f (Suc n) = G ((F^n) s) (f n)"
29 shows "horner F G n ((F^j') s) (f j') x = (\<Sum> j = 0..< n. -1^j * (1 / real (f (j' + j))) * x^j)"
30 proof (induct n arbitrary: i k j')
33 show ?case unfolding horner.simps Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc]
34 using horner_schema'[of "\<lambda> j. 1 / real (f (j' + j))"] by auto
38 assumes "0 \<le> Ifloat x" and f_Suc: "\<And>n. f (Suc n) = G ((F^n) s) (f n)"
39 and lb_0: "\<And> i k x. lb 0 i k x = 0"
40 and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 (int k) - x * (ub n (F i) (G i k) x)"
41 and ub_0: "\<And> i k x. ub 0 i k x = 0"
42 and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 (int k) - x * (lb n (F i) (G i k) x)"
43 shows "Ifloat (lb n ((F^j') s) (f j') x) \<le> horner F G n ((F^j') s) (f j') (Ifloat x) \<and>
44 horner F G n ((F^j') s) (f j') (Ifloat x) \<le> Ifloat (ub n ((F^j') s) (f j') x)"
45 (is "?lb n j' \<le> ?horner n j' \<and> ?horner n j' \<le> ?ub n j'")
46 proof (induct n arbitrary: j')
47 case 0 thus ?case unfolding lb_0 ub_0 horner.simps by auto
50 have "?lb (Suc n) j' \<le> ?horner (Suc n) j'" unfolding lb_Suc ub_Suc horner.simps Ifloat_sub diff_def
52 show "Ifloat (lapprox_rat prec 1 (int (f j'))) \<le> 1 / real (f j')" using lapprox_rat[of prec 1 "int (f j')"] by auto
53 from Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc, THEN conjunct2] `0 \<le> Ifloat x`
54 show "- Ifloat (x * ub n (F ((F ^ j') s)) (G ((F ^ j') s) (f j')) x) \<le> - (Ifloat x * horner F G n (F ((F ^ j') s)) (G ((F ^ j') s) (f j')) (Ifloat x))"
55 unfolding Ifloat_mult neg_le_iff_le by (rule mult_left_mono)
57 moreover have "?horner (Suc n) j' \<le> ?ub (Suc n) j'" unfolding ub_Suc ub_Suc horner.simps Ifloat_sub diff_def
59 show "1 / real (f j') \<le> Ifloat (rapprox_rat prec 1 (int (f j')))" using rapprox_rat[of 1 "int (f j')" prec] by auto
60 from Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc, THEN conjunct1] `0 \<le> Ifloat x`
61 show "- (Ifloat x * horner F G n (F ((F ^ j') s)) (G ((F ^ j') s) (f j')) (Ifloat x)) \<le>
62 - Ifloat (x * lb n (F ((F ^ j') s)) (G ((F ^ j') s) (f j')) x)"
63 unfolding Ifloat_mult neg_le_iff_le by (rule mult_left_mono)
65 ultimately show ?case by blast
68 subsection "Theorems for floating point functions implementing the horner scheme"
72 Here @{term_type "f :: nat \<Rightarrow> nat"} is the sequence defining the Taylor series, the coefficients are
73 all alternating and reciprocs. We use @{term G} and @{term F} to describe the computation of @{term f}.
77 lemma horner_bounds: fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat"
78 assumes "0 \<le> Ifloat x" and f_Suc: "\<And>n. f (Suc n) = G ((F^n) s) (f n)"
79 and lb_0: "\<And> i k x. lb 0 i k x = 0"
80 and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 (int k) - x * (ub n (F i) (G i k) x)"
81 and ub_0: "\<And> i k x. ub 0 i k x = 0"
82 and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 (int k) - x * (lb n (F i) (G i k) x)"
83 shows "Ifloat (lb n ((F^j') s) (f j') x) \<le> (\<Sum>j=0..<n. -1^j * (1 / real (f (j' + j))) * (Ifloat x)^j)" (is "?lb") and
84 "(\<Sum>j=0..<n. -1^j * (1 / real (f (j' + j))) * (Ifloat x)^j) \<le> Ifloat (ub n ((F^j') s) (f j') x)" (is "?ub")
87 using horner_bounds'[where lb=lb, OF `0 \<le> Ifloat x` f_Suc lb_0 lb_Suc ub_0 ub_Suc]
88 unfolding horner_schema[where f=f, OF f_Suc] .
89 thus "?lb" and "?ub" by auto
92 lemma horner_bounds_nonpos: fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat"
93 assumes "Ifloat x \<le> 0" and f_Suc: "\<And>n. f (Suc n) = G ((F^n) s) (f n)"
94 and lb_0: "\<And> i k x. lb 0 i k x = 0"
95 and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 (int k) + x * (ub n (F i) (G i k) x)"
96 and ub_0: "\<And> i k x. ub 0 i k x = 0"
97 and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 (int k) + x * (lb n (F i) (G i k) x)"
98 shows "Ifloat (lb n ((F^j') s) (f j') x) \<le> (\<Sum>j=0..<n. (1 / real (f (j' + j))) * (Ifloat x)^j)" (is "?lb") and
99 "(\<Sum>j=0..<n. (1 / real (f (j' + j))) * (Ifloat x)^j) \<le> Ifloat (ub n ((F^j') s) (f j') x)" (is "?ub")
101 { fix x y z :: float have "x - y * z = x + - y * z"
102 by (cases x, cases y, cases z, simp add: plus_float.simps minus_float.simps uminus_float.simps times_float.simps algebra_simps)
103 } note diff_mult_minus = this
105 { fix x :: float have "- (- x) = x" by (cases x, auto simp add: uminus_float.simps) } note minus_minus = this
107 have move_minus: "Ifloat (-x) = -1 * Ifloat x" by auto
109 have sum_eq: "(\<Sum>j=0..<n. (1 / real (f (j' + j))) * (Ifloat x)^j) =
110 (\<Sum>j = 0..<n. -1 ^ j * (1 / real (f (j' + j))) * Ifloat (- x) ^ j)"
111 proof (rule setsum_cong, simp)
112 fix j assume "j \<in> {0 ..< n}"
113 show "1 / real (f (j' + j)) * Ifloat x ^ j = -1 ^ j * (1 / real (f (j' + j))) * Ifloat (- x) ^ j"
114 unfolding move_minus power_mult_distrib real_mult_assoc[symmetric]
115 unfolding real_mult_commute unfolding real_mult_assoc[of "-1^j", symmetric] power_mult_distrib[symmetric]
119 have "0 \<le> Ifloat (-x)" using assms by auto
120 from horner_bounds[where G=G and F=F and f=f and s=s and prec=prec
121 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,
122 OF this f_Suc lb_0 refl ub_0 refl]
123 show "?lb" and "?ub" unfolding minus_minus sum_eq
127 subsection {* Selectors for next even or odd number *}
131 The horner scheme computes alternating series. To get the upper and lower bounds we need to
132 guarantee to access a even or odd member. To do this we use @{term get_odd} and @{term get_even}.
136 definition get_odd :: "nat \<Rightarrow> nat" where
137 "get_odd n = (if odd n then n else (Suc n))"
139 definition get_even :: "nat \<Rightarrow> nat" where
140 "get_even n = (if even n then n else (Suc n))"
142 lemma get_odd[simp]: "odd (get_odd n)" unfolding get_odd_def by (cases "odd n", auto)
143 lemma get_even[simp]: "even (get_even n)" unfolding get_even_def by (cases "even n", auto)
144 lemma get_odd_ex: "\<exists> k. Suc k = get_odd n \<and> odd (Suc k)"
145 proof (cases "odd n")
146 case True hence "0 < n" by (rule odd_pos)
147 from gr0_implies_Suc[OF this] obtain k where "Suc k = n" by auto
148 thus ?thesis unfolding get_odd_def if_P[OF True] using True[unfolded `Suc k = n`[symmetric]] by blast
150 case False hence "odd (Suc n)" by auto
151 thus ?thesis unfolding get_odd_def if_not_P[OF False] by blast
154 lemma get_even_double: "\<exists>i. get_even n = 2 * i" using get_even[unfolded even_mult_two_ex] .
155 lemma get_odd_double: "\<exists>i. get_odd n = 2 * i + 1" using get_odd[unfolded odd_Suc_mult_two_ex] by auto
157 section "Power function"
159 definition float_power_bnds :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where
160 "float_power_bnds n l u = (if odd n \<or> 0 < l then (l ^ n, u ^ n)
161 else if u < 0 then (u ^ n, l ^ n)
162 else (0, (max (-l) u) ^ n))"
164 lemma float_power_bnds: assumes "(l1, u1) = float_power_bnds n l u" and "x \<in> {Ifloat l .. Ifloat u}"
165 shows "x^n \<in> {Ifloat l1..Ifloat u1}"
166 proof (cases "even n")
169 proof (cases "0 < l")
170 case True hence "odd n \<or> 0 < l" and "0 \<le> Ifloat l" unfolding less_float_def by auto
171 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
172 have "Ifloat l^n \<le> x^n" and "x^n \<le> Ifloat u^n " using `0 \<le> Ifloat l` and assms unfolding atLeastAtMost_iff using power_mono[of "Ifloat l" x] power_mono[of x "Ifloat u"] by auto
173 thus ?thesis using assms `0 < l` unfolding atLeastAtMost_iff l1 u1 float_power less_float_def by auto
175 case False hence P: "\<not> (odd n \<or> 0 < l)" using `even n` by auto
177 proof (cases "u < 0")
178 case True hence "0 \<le> - Ifloat u" and "- Ifloat u \<le> - x" and "0 \<le> - x" and "-x \<le> - Ifloat l" using assms unfolding less_float_def by auto
179 hence "Ifloat u^n \<le> x^n" and "x^n \<le> Ifloat l^n" using power_mono[of "-x" "-Ifloat l" n] power_mono[of "-Ifloat u" "-x" n]
180 unfolding power_minus_even[OF `even n`] by auto
181 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
182 ultimately show ?thesis using float_power by auto
185 have "\<bar>x\<bar> \<le> Ifloat (max (-l) u)"
186 proof (cases "-l \<le> u")
187 case True thus ?thesis unfolding max_def if_P[OF True] using assms unfolding le_float_def by auto
189 case False thus ?thesis unfolding max_def if_not_P[OF False] using assms unfolding le_float_def by auto
191 hence x_abs: "\<bar>x\<bar> \<le> \<bar>Ifloat (max (-l) u)\<bar>" by auto
192 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
193 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
197 case False hence "odd n \<or> 0 < l" by auto
198 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
199 have "Ifloat l^n \<le> x^n" and "x^n \<le> Ifloat u^n " using assms unfolding atLeastAtMost_iff using power_mono_odd[OF False] by auto
200 thus ?thesis unfolding atLeastAtMost_iff l1 u1 float_power less_float_def by auto
203 lemma bnds_power: "\<forall> x l u. (l1, u1) = float_power_bnds n l u \<and> x \<in> {Ifloat l .. Ifloat u} \<longrightarrow> Ifloat l1 \<le> x^n \<and> x^n \<le> Ifloat u1"
204 using float_power_bnds by auto
206 section "Square root"
210 The square root computation is implemented as newton iteration. As first first step we use the
211 nearest power of two greater than the square root.
215 fun sqrt_iteration :: "nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
216 "sqrt_iteration prec 0 (Float m e) = Float 1 ((e + bitlen m) div 2 + 1)" |
217 "sqrt_iteration prec (Suc m) x = (let y = sqrt_iteration prec m x
218 in Float 1 -1 * (y + float_divr prec x y))"
220 definition ub_sqrt :: "nat \<Rightarrow> float \<Rightarrow> float option" where
221 "ub_sqrt prec x = (if 0 < x then Some (sqrt_iteration prec prec x) else if x < 0 then None else Some 0)"
223 definition lb_sqrt :: "nat \<Rightarrow> float \<Rightarrow> float option" where
224 "lb_sqrt prec x = (if 0 < x then Some (float_divl prec x (sqrt_iteration prec prec x)) else if x < 0 then None else Some 0)"
226 lemma sqrt_ub_pos_pos_1:
227 assumes "sqrt x < b" and "0 < b" and "0 < x"
228 shows "sqrt x < (b + x / b)/2"
230 from assms have "0 < (b - sqrt x) ^ 2 " by simp
231 also have "\<dots> = b ^ 2 - 2 * b * sqrt x + (sqrt x) ^ 2" by algebra
232 also have "\<dots> = b ^ 2 - 2 * b * sqrt x + x" using assms by (simp add: real_sqrt_pow2)
233 finally have "0 < b ^ 2 - 2 * b * sqrt x + x" by assumption
234 hence "0 < b / 2 - sqrt x + x / (2 * b)" using assms
235 by (simp add: field_simps power2_eq_square)
236 thus ?thesis by (simp add: field_simps)
239 lemma sqrt_iteration_bound: assumes "0 < Ifloat x"
240 shows "sqrt (Ifloat x) < Ifloat (sqrt_iteration prec n x)"
246 hence "0 < m" using float_pos_m_pos[unfolded less_float_def] assms by auto
247 hence "0 < sqrt (real m)" by auto
249 have int_nat_bl: "int (nat (bitlen m)) = bitlen m" using bitlen_ge0 by auto
251 have "Ifloat x = (real m / 2^nat (bitlen m)) * pow2 (e + int (nat (bitlen m)))"
252 unfolding pow2_add pow2_int Float Ifloat.simps by auto
253 also have "\<dots> < 1 * pow2 (e + int (nat (bitlen m)))"
254 proof (rule mult_strict_right_mono, auto)
255 show "real m < 2^nat (bitlen m)" using bitlen_bounds[OF `0 < m`, THEN conjunct2]
256 unfolding real_of_int_less_iff[of m, symmetric] by auto
258 finally have "sqrt (Ifloat x) < sqrt (pow2 (e + bitlen m))" unfolding int_nat_bl by auto
259 also have "\<dots> \<le> pow2 ((e + bitlen m) div 2 + 1)"
261 let ?E = "e + bitlen m"
262 have E_mod_pow: "pow2 (?E mod 2) < 4"
263 proof (cases "?E mod 2 = 1")
264 case True thus ?thesis by auto
267 have "0 \<le> ?E mod 2" by auto
268 have "?E mod 2 < 2" by auto
269 from this[THEN zless_imp_add1_zle]
270 have "?E mod 2 \<le> 0" using False by auto
271 from xt1(5)[OF `0 \<le> ?E mod 2` this]
274 hence "sqrt (pow2 (?E mod 2)) < sqrt (2 * 2)" by auto
275 hence E_mod_pow: "sqrt (pow2 (?E mod 2)) < 2" unfolding real_sqrt_abs2 by auto
277 have E_eq: "pow2 ?E = pow2 (?E div 2 + ?E div 2 + ?E mod 2)" by auto
278 have "sqrt (pow2 ?E) = sqrt (pow2 (?E div 2) * pow2 (?E div 2) * pow2 (?E mod 2))"
279 unfolding E_eq unfolding pow2_add ..
280 also have "\<dots> = pow2 (?E div 2) * sqrt (pow2 (?E mod 2))"
281 unfolding real_sqrt_mult[of _ "pow2 (?E mod 2)"] real_sqrt_abs2 by auto
282 also have "\<dots> < pow2 (?E div 2) * 2"
283 by (rule mult_strict_left_mono, auto intro: E_mod_pow)
284 also have "\<dots> = pow2 (?E div 2 + 1)" unfolding zadd_commute[of _ 1] pow2_add1 by auto
285 finally show ?thesis by auto
288 unfolding Float sqrt_iteration.simps Ifloat.simps by auto
292 let ?b = "sqrt_iteration prec n x"
293 have "0 < sqrt (Ifloat x)" using `0 < Ifloat x` by auto
294 also have "\<dots> < Ifloat ?b" using Suc .
295 finally have "sqrt (Ifloat x) < (Ifloat ?b + Ifloat x / Ifloat ?b)/2" using sqrt_ub_pos_pos_1[OF Suc _ `0 < Ifloat x`] by auto
296 also have "\<dots> \<le> (Ifloat ?b + Ifloat (float_divr prec x ?b))/2" by (rule divide_right_mono, auto simp add: float_divr)
297 also have "\<dots> = Ifloat (Float 1 -1) * (Ifloat ?b + Ifloat (float_divr prec x ?b))" by auto
298 finally show ?case unfolding sqrt_iteration.simps Let_def Ifloat_mult Ifloat_add right_distrib .
301 lemma sqrt_iteration_lower_bound: assumes "0 < Ifloat x"
302 shows "0 < Ifloat (sqrt_iteration prec n x)" (is "0 < ?sqrt")
304 have "0 < sqrt (Ifloat x)" using assms by auto
305 also have "\<dots> < ?sqrt" using sqrt_iteration_bound[OF assms] .
306 finally show ?thesis .
309 lemma lb_sqrt_lower_bound: assumes "0 \<le> Ifloat x"
310 shows "0 \<le> Ifloat (the (lb_sqrt prec x))"
311 proof (cases "0 < x")
312 case True hence "0 < Ifloat x" and "0 \<le> x" using `0 \<le> Ifloat x` unfolding less_float_def le_float_def by auto
313 hence "0 < sqrt_iteration prec prec x" unfolding less_float_def using sqrt_iteration_lower_bound by auto
314 hence "0 \<le> Ifloat (float_divl prec x (sqrt_iteration prec prec x))" using float_divl_lower_bound[OF `0 \<le> x`] unfolding le_float_def by auto
315 thus ?thesis unfolding lb_sqrt_def using True by auto
317 case False with `0 \<le> Ifloat x` have "Ifloat x = 0" unfolding less_float_def by auto
318 thus ?thesis unfolding lb_sqrt_def less_float_def by auto
321 lemma lb_sqrt_upper_bound: assumes "0 \<le> Ifloat x"
322 shows "Ifloat (the (lb_sqrt prec x)) \<le> sqrt (Ifloat x)"
323 proof (cases "0 < x")
324 case True hence "0 < Ifloat x" and "0 \<le> Ifloat x" unfolding less_float_def by auto
325 hence sqrt_gt0: "0 < sqrt (Ifloat x)" by auto
326 hence sqrt_ub: "sqrt (Ifloat x) < Ifloat (sqrt_iteration prec prec x)" using sqrt_iteration_bound by auto
328 have "Ifloat (float_divl prec x (sqrt_iteration prec prec x)) \<le> Ifloat x / Ifloat (sqrt_iteration prec prec x)" by (rule float_divl)
329 also have "\<dots> < Ifloat x / sqrt (Ifloat x)"
330 by (rule divide_strict_left_mono[OF sqrt_ub `0 < Ifloat x` mult_pos_pos[OF order_less_trans[OF sqrt_gt0 sqrt_ub] sqrt_gt0]])
331 also have "\<dots> = sqrt (Ifloat x)" unfolding inverse_eq_iff_eq[of _ "sqrt (Ifloat x)", symmetric] sqrt_divide_self_eq[OF `0 \<le> Ifloat x`, symmetric] by auto
332 finally show ?thesis unfolding lb_sqrt_def if_P[OF `0 < x`] by auto
334 case False with `0 \<le> Ifloat x`
335 have "\<not> x < 0" unfolding less_float_def le_float_def by auto
336 show ?thesis unfolding lb_sqrt_def if_not_P[OF False] if_not_P[OF `\<not> x < 0`] using assms by auto
339 lemma lb_sqrt: assumes "Some y = lb_sqrt prec x"
340 shows "Ifloat y \<le> sqrt (Ifloat x)" and "0 \<le> Ifloat x"
342 show "0 \<le> Ifloat x"
344 assume "\<not> 0 \<le> Ifloat x"
345 hence "lb_sqrt prec x = None" unfolding lb_sqrt_def less_float_def by auto
346 thus False using assms by auto
348 from lb_sqrt_upper_bound[OF this, of prec]
349 show "Ifloat y \<le> sqrt (Ifloat x)" unfolding assms[symmetric] by auto
352 lemma ub_sqrt_lower_bound: assumes "0 \<le> Ifloat x"
353 shows "sqrt (Ifloat x) \<le> Ifloat (the (ub_sqrt prec x))"
354 proof (cases "0 < x")
355 case True hence "0 < Ifloat x" unfolding less_float_def by auto
356 hence "0 < sqrt (Ifloat x)" by auto
357 hence "sqrt (Ifloat x) < Ifloat (sqrt_iteration prec prec x)" using sqrt_iteration_bound by auto
358 thus ?thesis unfolding ub_sqrt_def if_P[OF `0 < x`] by auto
360 case False with `0 \<le> Ifloat x`
361 have "Ifloat x = 0" unfolding less_float_def le_float_def by auto
362 thus ?thesis unfolding ub_sqrt_def less_float_def le_float_def by auto
365 lemma ub_sqrt: assumes "Some y = ub_sqrt prec x"
366 shows "sqrt (Ifloat x) \<le> Ifloat y" and "0 \<le> Ifloat x"
368 show "0 \<le> Ifloat x"
370 assume "\<not> 0 \<le> Ifloat x"
371 hence "ub_sqrt prec x = None" unfolding ub_sqrt_def less_float_def by auto
372 thus False using assms by auto
374 from ub_sqrt_lower_bound[OF this, of prec]
375 show "sqrt (Ifloat x) \<le> Ifloat y" unfolding assms[symmetric] by auto
378 lemma bnds_sqrt: "\<forall> x lx ux. (Some l, Some u) = (lb_sqrt prec lx, ub_sqrt prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux} \<longrightarrow> Ifloat l \<le> sqrt x \<and> sqrt x \<le> Ifloat u"
379 proof (rule allI, rule allI, rule allI, rule impI)
381 assume "(Some l, Some u) = (lb_sqrt prec lx, ub_sqrt prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux}"
382 hence l: "Some l = lb_sqrt prec lx " and u: "Some u = ub_sqrt prec ux" and x: "x \<in> {Ifloat lx .. Ifloat ux}" by auto
384 have "Ifloat lx \<le> x" and "x \<le> Ifloat ux" using x by auto
386 from lb_sqrt(1)[OF l] real_sqrt_le_mono[OF `Ifloat lx \<le> x`]
387 have "Ifloat l \<le> sqrt x" by (rule order_trans)
389 from real_sqrt_le_mono[OF `x \<le> Ifloat ux`] ub_sqrt(1)[OF u]
390 have "sqrt x \<le> Ifloat u" by (rule order_trans)
391 ultimately show "Ifloat l \<le> sqrt x \<and> sqrt x \<le> Ifloat u" ..
394 section "Arcus tangens and \<pi>"
396 subsection "Compute arcus tangens series"
400 As first step we implement the computation of the arcus tangens series. This is only valid in the range
401 @{term "{-1 :: real .. 1}"}. This is used to compute \<pi> and then the entire arcus tangens.
405 fun ub_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
406 and lb_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
407 "ub_arctan_horner prec 0 k x = 0"
408 | "ub_arctan_horner prec (Suc n) k x =
409 (rapprox_rat prec 1 (int k)) - x * (lb_arctan_horner prec n (k + 2) x)"
410 | "lb_arctan_horner prec 0 k x = 0"
411 | "lb_arctan_horner prec (Suc n) k x =
412 (lapprox_rat prec 1 (int k)) - x * (ub_arctan_horner prec n (k + 2) x)"
414 lemma arctan_0_1_bounds': assumes "0 \<le> Ifloat x" "Ifloat x \<le> 1" and "even n"
415 shows "arctan (Ifloat x) \<in> {Ifloat (x * lb_arctan_horner prec n 1 (x * x)) .. Ifloat (x * ub_arctan_horner prec (Suc n) 1 (x * x))}"
417 let "?c i" = "-1^i * (1 / real (i * 2 + 1) * Ifloat x ^ (i * 2 + 1))"
418 let "?S n" = "\<Sum> i=0..<n. ?c i"
420 have "0 \<le> Ifloat (x * x)" by auto
421 from `even n` obtain m where "2 * m = n" unfolding even_mult_two_ex by auto
423 have "arctan (Ifloat x) \<in> { ?S n .. ?S (Suc n) }"
424 proof (cases "Ifloat x = 0")
426 hence "0 < Ifloat x" using `0 \<le> Ifloat x` by auto
427 hence prem: "0 < 1 / real (0 * 2 + (1::nat)) * Ifloat x ^ (0 * 2 + 1)" by auto
429 have "\<bar> Ifloat x \<bar> \<le> 1" using `0 \<le> Ifloat x` `Ifloat x \<le> 1` by auto
430 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`]
431 show ?thesis unfolding arctan_series[OF `\<bar> Ifloat x \<bar> \<le> 1`] Suc_plus1 .
433 note arctan_bounds = this[unfolded atLeastAtMost_iff]
435 have F: "\<And>n. 2 * Suc n + 1 = 2 * n + 1 + 2" by auto
437 note bounds = horner_bounds[where s=1 and f="\<lambda>i. 2 * i + 1" and j'=0
438 and lb="\<lambda>n i k x. lb_arctan_horner prec n k x"
439 and ub="\<lambda>n i k x. ub_arctan_horner prec n k x",
440 OF `0 \<le> Ifloat (x*x)` F lb_arctan_horner.simps ub_arctan_horner.simps]
442 { have "Ifloat (x * lb_arctan_horner prec n 1 (x*x)) \<le> ?S n"
443 using bounds(1) `0 \<le> Ifloat x`
444 unfolding Ifloat_mult power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]
445 unfolding real_mult_commute mult_commute[of _ "2::nat"] power_mult power2_eq_square[of "Ifloat x"]
446 by (auto intro!: mult_left_mono)
447 also have "\<dots> \<le> arctan (Ifloat x)" using arctan_bounds ..
448 finally have "Ifloat (x * lb_arctan_horner prec n 1 (x*x)) \<le> arctan (Ifloat x)" . }
450 { have "arctan (Ifloat x) \<le> ?S (Suc n)" using arctan_bounds ..
451 also have "\<dots> \<le> Ifloat (x * ub_arctan_horner prec (Suc n) 1 (x*x))"
452 using bounds(2)[of "Suc n"] `0 \<le> Ifloat x`
453 unfolding Ifloat_mult power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]
454 unfolding real_mult_commute mult_commute[of _ "2::nat"] power_mult power2_eq_square[of "Ifloat x"]
455 by (auto intro!: mult_left_mono)
456 finally have "arctan (Ifloat x) \<le> Ifloat (x * ub_arctan_horner prec (Suc n) 1 (x*x))" . }
457 ultimately show ?thesis by auto
460 lemma arctan_0_1_bounds: assumes "0 \<le> Ifloat x" "Ifloat x \<le> 1"
461 shows "arctan (Ifloat x) \<in> {Ifloat (x * lb_arctan_horner prec (get_even n) 1 (x * x)) .. Ifloat (x * ub_arctan_horner prec (get_odd n) 1 (x * x))}"
462 proof (cases "even n")
464 obtain n' where "Suc n' = get_odd n" and "odd (Suc n')" using get_odd_ex by auto
465 hence "even n'" unfolding even_nat_Suc by auto
466 have "arctan (Ifloat x) \<le> Ifloat (x * ub_arctan_horner prec (get_odd n) 1 (x * x))"
467 unfolding `Suc n' = get_odd n`[symmetric] using arctan_0_1_bounds'[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1` `even n'`] by auto
469 have "Ifloat (x * lb_arctan_horner prec (get_even n) 1 (x * x)) \<le> arctan (Ifloat x)"
470 unfolding get_even_def if_P[OF True] using arctan_0_1_bounds'[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1` `even n`] by auto
471 ultimately show ?thesis by auto
473 case False hence "0 < n" by (rule odd_pos)
474 from gr0_implies_Suc[OF this] obtain n' where "n = Suc n'" ..
475 from False[unfolded this even_nat_Suc]
476 have "even n'" and "even (Suc (Suc n'))" by auto
477 have "get_odd n = Suc n'" unfolding get_odd_def if_P[OF False] using `n = Suc n'` .
479 have "arctan (Ifloat x) \<le> Ifloat (x * ub_arctan_horner prec (get_odd n) 1 (x * x))"
480 unfolding `get_odd n = Suc n'` using arctan_0_1_bounds'[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1` `even n'`] by auto
482 have "Ifloat (x * lb_arctan_horner prec (get_even n) 1 (x * x)) \<le> arctan (Ifloat x)"
483 unfolding get_even_def if_not_P[OF False] unfolding `n = Suc n'` using arctan_0_1_bounds'[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1` `even (Suc (Suc n'))`] by auto
484 ultimately show ?thesis by auto
487 subsection "Compute \<pi>"
489 definition ub_pi :: "nat \<Rightarrow> float" where
490 "ub_pi prec = (let A = rapprox_rat prec 1 5 ;
491 B = lapprox_rat prec 1 239
492 in ((Float 1 2) * ((Float 1 2) * A * (ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (A * A)) -
493 B * (lb_arctan_horner prec (get_even (prec div 14 + 1)) 1 (B * B)))))"
495 definition lb_pi :: "nat \<Rightarrow> float" where
496 "lb_pi prec = (let A = lapprox_rat prec 1 5 ;
497 B = rapprox_rat prec 1 239
498 in ((Float 1 2) * ((Float 1 2) * A * (lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (A * A)) -
499 B * (ub_arctan_horner prec (get_odd (prec div 14 + 1)) 1 (B * B)))))"
501 lemma pi_boundaries: "pi \<in> {Ifloat (lb_pi n) .. Ifloat (ub_pi n)}"
503 have machin_pi: "pi = 4 * (4 * arctan (1 / 5) - arctan (1 / 239))" unfolding machin[symmetric] by auto
505 { fix prec n :: nat fix k :: int assume "1 < k" hence "0 \<le> k" and "0 < k" and "1 \<le> k" by auto
506 let ?k = "rapprox_rat prec 1 k"
507 have "1 div k = 0" using div_pos_pos_trivial[OF _ `1 < k`] by auto
509 have "0 \<le> Ifloat ?k" by (rule order_trans[OF _ rapprox_rat], auto simp add: `0 \<le> k`)
510 have "Ifloat ?k \<le> 1" unfolding rapprox_rat.simps(2)[OF zero_le_one `0 < k`]
511 by (rule rapprox_posrat_le1, auto simp add: `0 < k` `1 \<le> k`)
513 have "1 / real k \<le> Ifloat ?k" using rapprox_rat[where x=1 and y=k] by auto
514 hence "arctan (1 / real k) \<le> arctan (Ifloat ?k)" by (rule arctan_monotone')
515 also have "\<dots> \<le> Ifloat (?k * ub_arctan_horner prec (get_odd n) 1 (?k * ?k))"
516 using arctan_0_1_bounds[OF `0 \<le> Ifloat ?k` `Ifloat ?k \<le> 1`] by auto
517 finally have "arctan (1 / (real k)) \<le> Ifloat (?k * ub_arctan_horner prec (get_odd n) 1 (?k * ?k))" .
518 } note ub_arctan = this
520 { fix prec n :: nat fix k :: int assume "1 < k" hence "0 \<le> k" and "0 < k" by auto
521 let ?k = "lapprox_rat prec 1 k"
522 have "1 div k = 0" using div_pos_pos_trivial[OF _ `1 < k`] by auto
523 have "1 / real k \<le> 1" using `1 < k` by auto
525 have "\<And>n. 0 \<le> Ifloat ?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`)
526 have "\<And>n. Ifloat ?k \<le> 1" using lapprox_rat by (rule order_trans, auto simp add: `1 / real k \<le> 1`)
528 have "Ifloat ?k \<le> 1 / real k" using lapprox_rat[where x=1 and y=k] by auto
530 have "Ifloat (?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k)) \<le> arctan (Ifloat ?k)"
531 using arctan_0_1_bounds[OF `0 \<le> Ifloat ?k` `Ifloat ?k \<le> 1`] by auto
532 also have "\<dots> \<le> arctan (1 / real k)" using `Ifloat ?k \<le> 1 / real k` by (rule arctan_monotone')
533 finally have "Ifloat (?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k)) \<le> arctan (1 / (real k))" .
534 } note lb_arctan = this
536 have "pi \<le> Ifloat (ub_pi n)"
537 unfolding ub_pi_def machin_pi Let_def Ifloat_mult Ifloat_sub unfolding Float_num
538 using lb_arctan[of 239] ub_arctan[of 5]
539 by (auto intro!: mult_left_mono add_mono simp add: diff_minus simp del: lapprox_rat.simps rapprox_rat.simps)
541 have "Ifloat (lb_pi n) \<le> pi"
542 unfolding lb_pi_def machin_pi Let_def Ifloat_mult Ifloat_sub Float_num
543 using lb_arctan[of 5] ub_arctan[of 239]
544 by (auto intro!: mult_left_mono add_mono simp add: diff_minus simp del: lapprox_rat.simps rapprox_rat.simps)
545 ultimately show ?thesis by auto
548 subsection "Compute arcus tangens in the entire domain"
550 function lb_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" and ub_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" where
551 "lb_arctan prec x = (let ub_horner = \<lambda> x. x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x) ;
552 lb_horner = \<lambda> x. x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)
553 in (if x < 0 then - ub_arctan prec (-x) else
554 if x \<le> Float 1 -1 then lb_horner x else
555 if x \<le> Float 1 1 then Float 1 1 * lb_horner (float_divl prec x (1 + the (ub_sqrt prec (1 + x * x))))
556 else (let inv = float_divr prec 1 x
558 else lb_pi prec * Float 1 -1 - ub_horner inv)))"
560 | "ub_arctan prec x = (let lb_horner = \<lambda> x. x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x) ;
561 ub_horner = \<lambda> x. x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)
562 in (if x < 0 then - lb_arctan prec (-x) else
563 if x \<le> Float 1 -1 then ub_horner x else
564 if x \<le> Float 1 1 then let y = float_divr prec x (1 + the (lb_sqrt prec (1 + x * x)))
565 in if y > 1 then ub_pi prec * Float 1 -1
566 else Float 1 1 * ub_horner y
567 else ub_pi prec * Float 1 -1 - lb_horner (float_divl prec 1 x)))"
568 by pat_completeness auto
569 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)
571 declare ub_arctan_horner.simps[simp del]
572 declare lb_arctan_horner.simps[simp del]
574 lemma lb_arctan_bound': assumes "0 \<le> Ifloat x"
575 shows "Ifloat (lb_arctan prec x) \<le> arctan (Ifloat x)"
577 have "\<not> x < 0" and "0 \<le> x" unfolding less_float_def le_float_def using `0 \<le> Ifloat x` by auto
578 let "?ub_horner x" = "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)"
579 and "?lb_horner x" = "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)"
582 proof (cases "x \<le> Float 1 -1")
583 case True hence "Ifloat x \<le> 1" unfolding le_float_def Float_num by auto
584 show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True]
585 using arctan_0_1_bounds[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1`] by auto
587 case False hence "0 < Ifloat x" unfolding le_float_def Float_num by auto
588 let ?R = "1 + sqrt (1 + Ifloat x * Ifloat x)"
589 let ?fR = "1 + the (ub_sqrt prec (1 + x * x))"
590 let ?DIV = "float_divl prec x ?fR"
592 have sqr_ge0: "0 \<le> 1 + Ifloat x * Ifloat x" using sum_power2_ge_zero[of 1 "Ifloat x", unfolded numeral_2_eq_2] by auto
593 hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg)
595 have "sqrt (Ifloat (1 + x * x)) \<le> Ifloat (the (ub_sqrt prec (1 + x * x)))" by (rule ub_sqrt_lower_bound, auto simp add: sqr_ge0)
596 hence "?R \<le> Ifloat ?fR" by auto
597 hence "0 < ?fR" and "0 < Ifloat ?fR" unfolding less_float_def using `0 < ?R` by auto
599 have monotone: "Ifloat (float_divl prec x ?fR) \<le> Ifloat x / ?R"
601 have "Ifloat ?DIV \<le> Ifloat x / Ifloat ?fR" by (rule float_divl)
602 also have "\<dots> \<le> Ifloat x / ?R" by (rule divide_left_mono[OF `?R \<le> Ifloat ?fR` `0 \<le> Ifloat x` mult_pos_pos[OF order_less_le_trans[OF divisor_gt0 `?R \<le> Ifloat ?fR`] divisor_gt0]])
603 finally show ?thesis .
607 proof (cases "x \<le> Float 1 1")
610 have "Ifloat x \<le> sqrt (Ifloat (1 + x * x))" using real_sqrt_sum_squares_ge2[where x=1, unfolded numeral_2_eq_2] by auto
611 also have "\<dots> \<le> Ifloat (the (ub_sqrt prec (1 + x * x)))" by (rule ub_sqrt_lower_bound, auto simp add: sqr_ge0)
612 finally have "Ifloat x \<le> Ifloat ?fR" by auto
613 moreover have "Ifloat ?DIV \<le> Ifloat x / Ifloat ?fR" by (rule float_divl)
614 ultimately have "Ifloat ?DIV \<le> 1" unfolding divide_le_eq_1_pos[OF `0 < Ifloat ?fR`, symmetric] by auto
616 have "0 \<le> Ifloat ?DIV" using float_divl_lower_bound[OF `0 \<le> x` `0 < ?fR`] unfolding le_float_def by auto
618 have "Ifloat (Float 1 1 * ?lb_horner ?DIV) \<le> 2 * arctan (Ifloat (float_divl prec x ?fR))" unfolding Ifloat_mult[of "Float 1 1"] Float_num
619 using arctan_0_1_bounds[OF `0 \<le> Ifloat ?DIV` `Ifloat ?DIV \<le> 1`] by auto
620 also have "\<dots> \<le> 2 * arctan (Ifloat x / ?R)"
621 using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono)
622 also have "2 * arctan (Ifloat x / ?R) = arctan (Ifloat x)" using arctan_half[symmetric] unfolding numeral_2_eq_2 power_Suc2 power_0 real_mult_1 .
623 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] .
626 hence "2 < Ifloat x" unfolding le_float_def Float_num by auto
627 hence "1 \<le> Ifloat x" by auto
629 let "?invx" = "float_divr prec 1 x"
630 have "0 \<le> arctan (Ifloat x)" using arctan_monotone'[OF `0 \<le> Ifloat x`] using arctan_tan[of 0, unfolded tan_zero] by auto
633 proof (cases "1 < ?invx")
635 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]
636 using `0 \<le> arctan (Ifloat x)` by auto
639 hence "Ifloat ?invx \<le> 1" unfolding less_float_def by auto
640 have "0 \<le> Ifloat ?invx" by (rule order_trans[OF _ float_divr], auto simp add: `0 \<le> Ifloat x`)
642 have "1 / Ifloat x \<noteq> 0" and "0 < 1 / Ifloat x" using `0 < Ifloat x` by auto
644 have "arctan (1 / Ifloat x) \<le> arctan (Ifloat ?invx)" unfolding Ifloat_1[symmetric] by (rule arctan_monotone', rule float_divr)
645 also have "\<dots> \<le> Ifloat (?ub_horner ?invx)" using arctan_0_1_bounds[OF `0 \<le> Ifloat ?invx` `Ifloat ?invx \<le> 1`] by auto
646 finally have "pi / 2 - Ifloat (?ub_horner ?invx) \<le> arctan (Ifloat x)"
647 using `0 \<le> arctan (Ifloat x)` arctan_inverse[OF `1 / Ifloat x \<noteq> 0`]
648 unfolding real_sgn_pos[OF `0 < 1 / Ifloat x`] le_diff_eq by auto
650 have "Ifloat (lb_pi prec * Float 1 -1) \<le> pi / 2" unfolding Ifloat_mult Float_num times_divide_eq_right real_mult_1 using pi_boundaries by auto
652 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]
659 lemma ub_arctan_bound': assumes "0 \<le> Ifloat x"
660 shows "arctan (Ifloat x) \<le> Ifloat (ub_arctan prec x)"
662 have "\<not> x < 0" and "0 \<le> x" unfolding less_float_def le_float_def using `0 \<le> Ifloat x` by auto
664 let "?ub_horner x" = "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)"
665 and "?lb_horner x" = "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)"
668 proof (cases "x \<le> Float 1 -1")
669 case True hence "Ifloat x \<le> 1" unfolding le_float_def Float_num by auto
670 show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True]
671 using arctan_0_1_bounds[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1`] by auto
673 case False hence "0 < Ifloat x" unfolding le_float_def Float_num by auto
674 let ?R = "1 + sqrt (1 + Ifloat x * Ifloat x)"
675 let ?fR = "1 + the (lb_sqrt prec (1 + x * x))"
676 let ?DIV = "float_divr prec x ?fR"
678 have sqr_ge0: "0 \<le> 1 + Ifloat x * Ifloat x" using sum_power2_ge_zero[of 1 "Ifloat x", unfolded numeral_2_eq_2] by auto
679 hence "0 \<le> Ifloat (1 + x*x)" by auto
681 hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg)
683 have "Ifloat (the (lb_sqrt prec (1 + x * x))) \<le> sqrt (Ifloat (1 + x * x))" by (rule lb_sqrt_upper_bound, auto simp add: sqr_ge0)
684 hence "Ifloat ?fR \<le> ?R" by auto
685 have "0 < Ifloat ?fR" unfolding Ifloat_add Ifloat_1 by (rule order_less_le_trans[OF zero_less_one], auto simp add: lb_sqrt_lower_bound[OF `0 \<le> Ifloat (1 + x*x)`])
687 have monotone: "Ifloat x / ?R \<le> Ifloat (float_divr prec x ?fR)"
689 from divide_left_mono[OF `Ifloat ?fR \<le> ?R` `0 \<le> Ifloat x` mult_pos_pos[OF divisor_gt0 `0 < Ifloat ?fR`]]
690 have "Ifloat x / ?R \<le> Ifloat x / Ifloat ?fR" .
691 also have "\<dots> \<le> Ifloat ?DIV" by (rule float_divr)
692 finally show ?thesis .
696 proof (cases "x \<le> Float 1 1")
699 proof (cases "?DIV > 1")
701 have "pi / 2 \<le> Ifloat (ub_pi prec * Float 1 -1)" unfolding Ifloat_mult Float_num times_divide_eq_right real_mult_1 using pi_boundaries by auto
702 from order_less_le_trans[OF arctan_ubound this, THEN less_imp_le]
703 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] .
706 hence "Ifloat ?DIV \<le> 1" unfolding less_float_def by auto
708 have "0 \<le> Ifloat x / ?R" using `0 \<le> Ifloat x` `0 < ?R` unfolding real_0_le_divide_iff by auto
709 hence "0 \<le> Ifloat ?DIV" using monotone by (rule order_trans)
711 have "arctan (Ifloat x) = 2 * arctan (Ifloat x / ?R)" using arctan_half unfolding numeral_2_eq_2 power_Suc2 power_0 real_mult_1 .
712 also have "\<dots> \<le> 2 * arctan (Ifloat ?DIV)"
713 using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono)
714 also have "\<dots> \<le> Ifloat (Float 1 1 * ?ub_horner ?DIV)" unfolding Ifloat_mult[of "Float 1 1"] Float_num
715 using arctan_0_1_bounds[OF `0 \<le> Ifloat ?DIV` `Ifloat ?DIV \<le> 1`] by auto
716 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] .
720 hence "2 < Ifloat x" unfolding le_float_def Float_num by auto
721 hence "1 \<le> Ifloat x" by auto
722 hence "0 < Ifloat x" by auto
723 hence "0 < x" unfolding less_float_def by auto
725 let "?invx" = "float_divl prec 1 x"
726 have "0 \<le> arctan (Ifloat x)" using arctan_monotone'[OF `0 \<le> Ifloat x`] using arctan_tan[of 0, unfolded tan_zero] by auto
728 have "Ifloat ?invx \<le> 1" unfolding less_float_def by (rule order_trans[OF float_divl], auto simp add: `1 \<le> Ifloat x` divide_le_eq_1_pos[OF `0 < Ifloat x`])
729 have "0 \<le> Ifloat ?invx" unfolding Ifloat_0[symmetric] by (rule float_divl_lower_bound[unfolded le_float_def], auto simp add: `0 < x`)
731 have "1 / Ifloat x \<noteq> 0" and "0 < 1 / Ifloat x" using `0 < Ifloat x` by auto
733 have "Ifloat (?lb_horner ?invx) \<le> arctan (Ifloat ?invx)" using arctan_0_1_bounds[OF `0 \<le> Ifloat ?invx` `Ifloat ?invx \<le> 1`] by auto
734 also have "\<dots> \<le> arctan (1 / Ifloat x)" unfolding Ifloat_1[symmetric] by (rule arctan_monotone', rule float_divl)
735 finally have "arctan (Ifloat x) \<le> pi / 2 - Ifloat (?lb_horner ?invx)"
736 using `0 \<le> arctan (Ifloat x)` arctan_inverse[OF `1 / Ifloat x \<noteq> 0`]
737 unfolding real_sgn_pos[OF `0 < 1 / Ifloat x`] le_diff_eq by auto
739 have "pi / 2 \<le> Ifloat (ub_pi prec * Float 1 -1)" unfolding Ifloat_mult Float_num times_divide_eq_right mult_1_right using pi_boundaries by auto
741 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]
747 lemma arctan_boundaries:
748 "arctan (Ifloat x) \<in> {Ifloat (lb_arctan prec x) .. Ifloat (ub_arctan prec x)}"
749 proof (cases "0 \<le> x")
750 case True hence "0 \<le> Ifloat x" unfolding le_float_def by auto
751 show ?thesis using ub_arctan_bound'[OF `0 \<le> Ifloat x`] lb_arctan_bound'[OF `0 \<le> Ifloat x`] unfolding atLeastAtMost_iff by auto
754 case False hence "x < 0" and "0 \<le> Ifloat ?mx" unfolding le_float_def less_float_def by auto
755 hence bounds: "Ifloat (lb_arctan prec ?mx) \<le> arctan (Ifloat ?mx) \<and> arctan (Ifloat ?mx) \<le> Ifloat (ub_arctan prec ?mx)"
756 using ub_arctan_bound'[OF `0 \<le> Ifloat ?mx`] lb_arctan_bound'[OF `0 \<le> Ifloat ?mx`] by auto
757 show ?thesis unfolding Ifloat_minus arctan_minus lb_arctan.simps[where x=x] ub_arctan.simps[where x=x] Let_def if_P[OF `x < 0`]
758 unfolding atLeastAtMost_iff using bounds[unfolded Ifloat_minus arctan_minus] by auto
761 lemma bnds_arctan: "\<forall> x lx ux. (l, u) = (lb_arctan prec lx, ub_arctan prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux} \<longrightarrow> Ifloat l \<le> arctan x \<and> arctan x \<le> Ifloat u"
762 proof (rule allI, rule allI, rule allI, rule impI)
764 assume "(l, u) = (lb_arctan prec lx, ub_arctan prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux}"
765 hence l: "lb_arctan prec lx = l " and u: "ub_arctan prec ux = u" and x: "x \<in> {Ifloat lx .. Ifloat ux}" by auto
767 { from arctan_boundaries[of lx prec, unfolded l]
768 have "Ifloat l \<le> arctan (Ifloat lx)" by (auto simp del: lb_arctan.simps)
769 also have "\<dots> \<le> arctan x" using x by (auto intro: arctan_monotone')
770 finally have "Ifloat l \<le> arctan x" .
772 { have "arctan x \<le> arctan (Ifloat ux)" using x by (auto intro: arctan_monotone')
773 also have "\<dots> \<le> Ifloat u" using arctan_boundaries[of ux prec, unfolded u] by (auto simp del: ub_arctan.simps)
774 finally have "arctan x \<le> Ifloat u" .
775 } ultimately show "Ifloat l \<le> arctan x \<and> arctan x \<le> Ifloat u" ..
778 section "Sinus and Cosinus"
780 subsection "Compute the cosinus and sinus series"
782 fun ub_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
783 and lb_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
784 "ub_sin_cos_aux prec 0 i k x = 0"
785 | "ub_sin_cos_aux prec (Suc n) i k x =
786 (rapprox_rat prec 1 (int k)) - x * (lb_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)"
787 | "lb_sin_cos_aux prec 0 i k x = 0"
788 | "lb_sin_cos_aux prec (Suc n) i k x =
789 (lapprox_rat prec 1 (int k)) - x * (ub_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)"
792 shows "Ifloat (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> (\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i))) * (Ifloat x)^(2 * i))" (is "?lb")
793 and "(\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i))) * (Ifloat x)^(2 * i)) \<le> Ifloat (ub_sin_cos_aux prec n 1 1 (x * x))" (is "?ub")
795 have "0 \<le> Ifloat (x * x)" unfolding Ifloat_mult by auto
796 let "?f n" = "fact (2 * n)"
799 have F: "\<And>m. ((\<lambda>i. i + 2) ^ n) m = m + 2 * n" by (induct n arbitrary: m, auto)
800 have "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) ^ n) 1 * (((\<lambda>i. i + 2) ^ n) 1 + 1)"
801 unfolding F by auto } note f_eq = this
803 from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0,
804 OF `0 \<le> Ifloat (x * x)` f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps]
805 show "?lb" and "?ub" by (auto simp add: power_mult power2_eq_square[of "Ifloat x"])
808 lemma cos_boundaries: assumes "0 \<le> Ifloat x" and "Ifloat x \<le> pi / 2"
809 shows "cos (Ifloat x) \<in> {Ifloat (lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) .. Ifloat (ub_sin_cos_aux prec (get_odd n) 1 1 (x * x))}"
810 proof (cases "Ifloat x = 0")
811 case False hence "Ifloat x \<noteq> 0" by auto
812 hence "0 < x" and "0 < Ifloat x" using `0 \<le> Ifloat x` unfolding less_float_def by auto
813 have "0 < x * x" using `0 < x` unfolding less_float_def Ifloat_mult Ifloat_0
814 using mult_pos_pos[where a="Ifloat x" and b="Ifloat x"] by auto
816 { fix x n have "(\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i))) * x^(2 * i))
817 = (\<Sum> i = 0 ..< 2 * n. (if even(i) then (-1 ^ (i div 2))/(real (fact i)) else 0) * x ^ i)" (is "?sum = ?ifsum")
819 have "?sum = ?sum + (\<Sum> j = 0 ..< n. 0)" by auto
821 (\<Sum> j = 0 ..< n. -1 ^ ((2 * j) div 2) / (real (fact (2 * j))) * x ^(2 * j)) + (\<Sum> j = 0 ..< n. 0)" by auto
822 also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. if even i then -1 ^ (i div 2) / (real (fact i)) * x ^ i else 0)"
823 unfolding sum_split_even_odd ..
824 also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. (if even i then -1 ^ (i div 2) / (real (fact i)) else 0) * x ^ i)"
825 by (rule setsum_cong2) auto
826 finally show ?thesis by assumption
827 qed } note morph_to_if_power = this
830 { fix n :: nat assume "0 < n"
831 hence "0 < 2 * n" by auto
832 obtain t where "0 < t" and "t < Ifloat x" and
833 cos_eq: "cos (Ifloat x) = (\<Sum> i = 0 ..< 2 * n. (if even(i) then (-1 ^ (i div 2))/(real (fact i)) else 0) * (Ifloat x) ^ i)
834 + (cos (t + 1/2 * real (2 * n) * pi) / real (fact (2*n))) * (Ifloat x)^(2*n)"
835 (is "_ = ?SUM + ?rest / ?fact * ?pow")
836 using Maclaurin_cos_expansion2[OF `0 < Ifloat x` `0 < 2 * n`] by auto
838 have "cos t * -1^n = cos t * cos (real n * pi) + sin t * sin (real n * pi)" by auto
839 also have "\<dots> = cos (t + real n * pi)" using cos_add by auto
840 also have "\<dots> = ?rest" by auto
841 finally have "cos t * -1^n = ?rest" .
843 have "t \<le> pi / 2" using `t < Ifloat x` and `Ifloat x \<le> pi / 2` by auto
844 hence "0 \<le> cos t" using `0 < t` and cos_ge_zero by auto
845 ultimately have even: "even n \<Longrightarrow> 0 \<le> ?rest" and odd: "odd n \<Longrightarrow> 0 \<le> - ?rest " by auto
847 have "0 < ?fact" by auto
848 have "0 < ?pow" using `0 < Ifloat x` by auto
852 have "Ifloat (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> ?SUM"
853 unfolding morph_to_if_power[symmetric] using cos_aux by auto
854 also have "\<dots> \<le> cos (Ifloat x)"
856 from even[OF `even n`] `0 < ?fact` `0 < ?pow`
857 have "0 \<le> (?rest / ?fact) * ?pow" by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
858 thus ?thesis unfolding cos_eq by auto
860 finally have "Ifloat (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> cos (Ifloat x)" .
865 have "cos (Ifloat x) \<le> ?SUM"
867 from `0 < ?fact` and `0 < ?pow` and odd[OF `odd n`]
868 have "0 \<le> (- ?rest) / ?fact * ?pow"
869 by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
870 thus ?thesis unfolding cos_eq by auto
872 also have "\<dots> \<le> Ifloat (ub_sin_cos_aux prec n 1 1 (x * x))"
873 unfolding morph_to_if_power[symmetric] using cos_aux by auto
874 finally have "cos (Ifloat x) \<le> Ifloat (ub_sin_cos_aux prec n 1 1 (x * x))" .
875 } note ub = this and lb
876 } note ub = this(1) and lb = this(2)
878 have "cos (Ifloat x) \<le> Ifloat (ub_sin_cos_aux prec (get_odd n) 1 1 (x * x))" using ub[OF odd_pos[OF get_odd] get_odd] .
879 moreover have "Ifloat (lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) \<le> cos (Ifloat x)"
880 proof (cases "0 < get_even n")
881 case True show ?thesis using lb[OF True get_even] .
884 hence "get_even n = 0" by auto
885 have "- (pi / 2) \<le> Ifloat x" by (rule order_trans[OF _ `0 < Ifloat x`[THEN less_imp_le]], auto)
886 with `Ifloat x \<le> pi / 2`
887 show ?thesis unfolding `get_even n = 0` lb_sin_cos_aux.simps Ifloat_minus Ifloat_0 using cos_ge_zero by auto
889 ultimately show ?thesis by auto
893 proof (cases "n = 0")
895 thus ?thesis unfolding `n = 0` get_even_def get_odd_def using `Ifloat x = 0` lapprox_rat[where x="-1" and y=1] by auto
897 case False with not0_implies_Suc obtain m where "n = Suc m" by blast
898 thus ?thesis unfolding `n = Suc m` get_even_def get_odd_def using `Ifloat x = 0` rapprox_rat[where x=1 and y=1] lapprox_rat[where x=1 and y=1] by (cases "even (Suc m)", auto)
902 lemma sin_aux: assumes "0 \<le> Ifloat x"
903 shows "Ifloat (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> (\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i + 1))) * (Ifloat x)^(2 * i + 1))" (is "?lb")
904 and "(\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i + 1))) * (Ifloat x)^(2 * i + 1)) \<le> Ifloat (x * ub_sin_cos_aux prec n 2 1 (x * x))" (is "?ub")
906 have "0 \<le> Ifloat (x * x)" unfolding Ifloat_mult by auto
907 let "?f n" = "fact (2 * n + 1)"
910 have F: "\<And>m. ((\<lambda>i. i + 2) ^ n) m = m + 2 * n" by (induct n arbitrary: m, auto)
911 have "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) ^ n) 2 * (((\<lambda>i. i + 2) ^ n) 2 + 1)"
912 unfolding F by auto } note f_eq = this
914 from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0,
915 OF `0 \<le> Ifloat (x * x)` f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps]
916 show "?lb" and "?ub" using `0 \<le> Ifloat x` unfolding Ifloat_mult
917 unfolding power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]
918 unfolding real_mult_commute
919 by (auto intro!: mult_left_mono simp add: power_mult power2_eq_square[of "Ifloat x"])
922 lemma sin_boundaries: assumes "0 \<le> Ifloat x" and "Ifloat x \<le> pi / 2"
923 shows "sin (Ifloat x) \<in> {Ifloat (x * lb_sin_cos_aux prec (get_even n) 2 1 (x * x)) .. Ifloat (x * ub_sin_cos_aux prec (get_odd n) 2 1 (x * x))}"
924 proof (cases "Ifloat x = 0")
925 case False hence "Ifloat x \<noteq> 0" by auto
926 hence "0 < x" and "0 < Ifloat x" using `0 \<le> Ifloat x` unfolding less_float_def by auto
927 have "0 < x * x" using `0 < x` unfolding less_float_def Ifloat_mult Ifloat_0
928 using mult_pos_pos[where a="Ifloat x" and b="Ifloat x"] by auto
930 { 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))
931 = (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * x ^ i)" (is "?SUM = _")
933 have pow: "!!i. x ^ (2 * i + 1) = x * x ^ (2 * i)" by auto
934 have "?SUM = (\<Sum> j = 0 ..< n. 0) + ?SUM" by auto
935 also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. if even i then 0 else -1 ^ ((i - Suc 0) div 2) / (real (fact i)) * x ^ i)"
936 unfolding sum_split_even_odd ..
937 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)"
938 by (rule setsum_cong2) auto
939 finally show ?thesis by assumption
940 qed } note setsum_morph = this
942 { fix n :: nat assume "0 < n"
943 hence "0 < 2 * n + 1" by auto
944 obtain t where "0 < t" and "t < Ifloat x" and
945 sin_eq: "sin (Ifloat x) = (\<Sum> i = 0 ..< 2 * n + 1. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * (Ifloat x) ^ i)
946 + (sin (t + 1/2 * real (2 * n + 1) * pi) / real (fact (2*n + 1))) * (Ifloat x)^(2*n + 1)"
947 (is "_ = ?SUM + ?rest / ?fact * ?pow")
948 using Maclaurin_sin_expansion3[OF `0 < 2 * n + 1` `0 < Ifloat x`] by auto
950 have "?rest = cos t * -1^n" unfolding sin_add cos_add real_of_nat_add left_distrib right_distrib by auto
952 have "t \<le> pi / 2" using `t < Ifloat x` and `Ifloat x \<le> pi / 2` by auto
953 hence "0 \<le> cos t" using `0 < t` and cos_ge_zero by auto
954 ultimately have even: "even n \<Longrightarrow> 0 \<le> ?rest" and odd: "odd n \<Longrightarrow> 0 \<le> - ?rest " by auto
956 have "0 < ?fact" by (rule real_of_nat_fact_gt_zero)
957 have "0 < ?pow" using `0 < Ifloat x` by (rule zero_less_power)
961 have "Ifloat (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le>
962 (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * (Ifloat x) ^ i)"
963 using sin_aux[OF `0 \<le> Ifloat x`] unfolding setsum_morph[symmetric] by auto
964 also have "\<dots> \<le> ?SUM" by auto
965 also have "\<dots> \<le> sin (Ifloat x)"
967 from even[OF `even n`] `0 < ?fact` `0 < ?pow`
968 have "0 \<le> (?rest / ?fact) * ?pow" by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
969 thus ?thesis unfolding sin_eq by auto
971 finally have "Ifloat (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> sin (Ifloat x)" .
976 have "sin (Ifloat x) \<le> ?SUM"
978 from `0 < ?fact` and `0 < ?pow` and odd[OF `odd n`]
979 have "0 \<le> (- ?rest) / ?fact * ?pow"
980 by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
981 thus ?thesis unfolding sin_eq by auto
983 also have "\<dots> \<le> (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * (Ifloat x) ^ i)"
985 also have "\<dots> \<le> Ifloat (x * ub_sin_cos_aux prec n 2 1 (x * x))"
986 using sin_aux[OF `0 \<le> Ifloat x`] unfolding setsum_morph[symmetric] by auto
987 finally have "sin (Ifloat x) \<le> Ifloat (x * ub_sin_cos_aux prec n 2 1 (x * x))" .
988 } note ub = this and lb
989 } note ub = this(1) and lb = this(2)
991 have "sin (Ifloat x) \<le> Ifloat (x * ub_sin_cos_aux prec (get_odd n) 2 1 (x * x))" using ub[OF odd_pos[OF get_odd] get_odd] .
992 moreover have "Ifloat (x * lb_sin_cos_aux prec (get_even n) 2 1 (x * x)) \<le> sin (Ifloat x)"
993 proof (cases "0 < get_even n")
994 case True show ?thesis using lb[OF True get_even] .
997 hence "get_even n = 0" by auto
998 with `Ifloat x \<le> pi / 2` `0 \<le> Ifloat x`
999 show ?thesis unfolding `get_even n = 0` ub_sin_cos_aux.simps Ifloat_minus Ifloat_0 using sin_ge_zero by auto
1001 ultimately show ?thesis by auto
1005 proof (cases "n = 0")
1007 thus ?thesis unfolding `n = 0` get_even_def get_odd_def using `Ifloat x = 0` lapprox_rat[where x="-1" and y=1] by auto
1009 case False with not0_implies_Suc obtain m where "n = Suc m" by blast
1010 thus ?thesis unfolding `n = Suc m` get_even_def get_odd_def using `Ifloat x = 0` rapprox_rat[where x=1 and y=1] lapprox_rat[where x=1 and y=1] by (cases "even (Suc m)", auto)
1014 subsection "Compute the cosinus in the entire domain"
1016 definition lb_cos :: "nat \<Rightarrow> float \<Rightarrow> float" where
1017 "lb_cos prec x = (let
1018 horner = \<lambda> x. lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 1 1 (x * x) ;
1019 half = \<lambda> x. if x < 0 then - 1 else Float 1 1 * x * x - 1
1020 in if x < Float 1 -1 then horner x
1021 else if x < 1 then half (horner (x * Float 1 -1))
1022 else half (half (horner (x * Float 1 -2))))"
1024 definition ub_cos :: "nat \<Rightarrow> float \<Rightarrow> float" where
1025 "ub_cos prec x = (let
1026 horner = \<lambda> x. ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x) ;
1027 half = \<lambda> x. Float 1 1 * x * x - 1
1028 in if x < Float 1 -1 then horner x
1029 else if x < 1 then half (horner (x * Float 1 -1))
1030 else half (half (horner (x * Float 1 -2))))"
1032 definition bnds_cos :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where
1033 "bnds_cos prec lx ux = (let lpi = lb_pi prec
1034 in if lx < -lpi \<or> ux > lpi then (Float -1 0, Float 1 0)
1035 else if ux \<le> 0 then (lb_cos prec (-lx), ub_cos prec (-ux))
1036 else if 0 \<le> lx then (lb_cos prec ux, ub_cos prec lx)
1037 else (min (lb_cos prec (-lx)) (lb_cos prec ux), Float 1 0))"
1039 lemma lb_cos: assumes "0 \<le> Ifloat x" and "Ifloat x \<le> pi"
1040 shows "cos (Ifloat x) \<in> {Ifloat (lb_cos prec x) .. Ifloat (ub_cos prec x)}" (is "?cos x \<in> { Ifloat (?lb x) .. Ifloat (?ub x) }")
1043 have "cos x = cos (x / 2 + x / 2)" by auto
1044 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"
1045 unfolding cos_add by auto
1046 also have "\<dots> = 2 * cos (x / 2) * cos (x / 2) - 1" by algebra
1047 finally have "cos x = 2 * cos (x / 2) * cos (x / 2) - 1" .
1048 } note x_half = this[symmetric]
1050 have "\<not> x < 0" using `0 \<le> Ifloat x` unfolding less_float_def by auto
1051 let "?ub_horner x" = "ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x)"
1052 let "?lb_horner x" = "lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 1 1 (x * x)"
1053 let "?ub_half x" = "Float 1 1 * x * x - 1"
1054 let "?lb_half x" = "if x < 0 then - 1 else Float 1 1 * x * x - 1"
1057 proof (cases "x < Float 1 -1")
1058 case True hence "Ifloat x \<le> pi / 2" unfolding less_float_def using pi_ge_two by auto
1059 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
1060 using cos_boundaries[OF `0 \<le> Ifloat x` `Ifloat x \<le> pi / 2`] .
1064 { fix y x :: float let ?x2 = "Ifloat (x * Float 1 -1)"
1065 assume "Ifloat y \<le> cos ?x2" and "-pi \<le> Ifloat x" and "Ifloat x \<le> pi"
1066 hence "- (pi / 2) \<le> ?x2" and "?x2 \<le> pi / 2" using pi_ge_two unfolding Ifloat_mult Float_num by auto
1067 hence "0 \<le> cos ?x2" by (rule cos_ge_zero)
1069 have "Ifloat (?lb_half y) \<le> cos (Ifloat x)"
1070 proof (cases "y < 0")
1071 case True show ?thesis using cos_ge_minus_one unfolding if_P[OF True] by auto
1074 hence "0 \<le> Ifloat y" unfolding less_float_def by auto
1075 from mult_mono[OF `Ifloat y \<le> cos ?x2` `Ifloat y \<le> cos ?x2` `0 \<le> cos ?x2` this]
1076 have "Ifloat y * Ifloat y \<le> cos ?x2 * cos ?x2" .
1077 hence "2 * Ifloat y * Ifloat y \<le> 2 * cos ?x2 * cos ?x2" by auto
1078 hence "2 * Ifloat y * Ifloat y - 1 \<le> 2 * cos (Ifloat x / 2) * cos (Ifloat x / 2) - 1" unfolding Float_num Ifloat_mult by auto
1079 thus ?thesis unfolding if_not_P[OF False] x_half Float_num Ifloat_mult Ifloat_sub by auto
1081 } note lb_half = this
1083 { fix y x :: float let ?x2 = "Ifloat (x * Float 1 -1)"
1084 assume ub: "cos ?x2 \<le> Ifloat y" and "- pi \<le> Ifloat x" and "Ifloat x \<le> pi"
1085 hence "- (pi / 2) \<le> ?x2" and "?x2 \<le> pi / 2" using pi_ge_two unfolding Ifloat_mult Float_num by auto
1086 hence "0 \<le> cos ?x2" by (rule cos_ge_zero)
1088 have "cos (Ifloat x) \<le> Ifloat (?ub_half y)"
1090 have "0 \<le> Ifloat y" using `0 \<le> cos ?x2` ub by (rule order_trans)
1091 from mult_mono[OF ub ub this `0 \<le> cos ?x2`]
1092 have "cos ?x2 * cos ?x2 \<le> Ifloat y * Ifloat y" .
1093 hence "2 * cos ?x2 * cos ?x2 \<le> 2 * Ifloat y * Ifloat y" by auto
1094 hence "2 * cos (Ifloat x / 2) * cos (Ifloat x / 2) - 1 \<le> 2 * Ifloat y * Ifloat y - 1" unfolding Float_num Ifloat_mult by auto
1095 thus ?thesis unfolding x_half Ifloat_mult Float_num Ifloat_sub by auto
1097 } note ub_half = this
1099 let ?x2 = "x * Float 1 -1"
1100 let ?x4 = "x * Float 1 -1 * Float 1 -1"
1102 have "-pi \<le> Ifloat x" using pi_ge_zero[THEN le_imp_neg_le, unfolded minus_zero] `0 \<le> Ifloat x` by (rule order_trans)
1105 proof (cases "x < 1")
1106 case True hence "Ifloat x \<le> 1" unfolding less_float_def by auto
1107 have "0 \<le> Ifloat ?x2" and "Ifloat ?x2 \<le> pi / 2" using pi_ge_two `0 \<le> Ifloat x` unfolding Ifloat_mult Float_num using assms by auto
1108 from cos_boundaries[OF this]
1109 have lb: "Ifloat (?lb_horner ?x2) \<le> ?cos ?x2" and ub: "?cos ?x2 \<le> Ifloat (?ub_horner ?x2)" by auto
1111 have "Ifloat (?lb x) \<le> ?cos x"
1113 from lb_half[OF lb `-pi \<le> Ifloat x` `Ifloat x \<le> pi`]
1114 show ?thesis unfolding lb_cos_def[where x=x] Let_def using `\<not> x < 0` `\<not> x < Float 1 -1` `x < 1` by auto
1116 moreover have "?cos x \<le> Ifloat (?ub x)"
1118 from ub_half[OF ub `-pi \<le> Ifloat x` `Ifloat x \<le> pi`]
1119 show ?thesis unfolding ub_cos_def[where x=x] Let_def using `\<not> x < 0` `\<not> x < Float 1 -1` `x < 1` by auto
1121 ultimately show ?thesis by auto
1124 have "0 \<le> Ifloat ?x4" and "Ifloat ?x4 \<le> pi / 2" using pi_ge_two `0 \<le> Ifloat x` `Ifloat x \<le> pi` unfolding Ifloat_mult Float_num by auto
1125 from cos_boundaries[OF this]
1126 have lb: "Ifloat (?lb_horner ?x4) \<le> ?cos ?x4" and ub: "?cos ?x4 \<le> Ifloat (?ub_horner ?x4)" by auto
1128 have eq_4: "?x2 * Float 1 -1 = x * Float 1 -2" by (cases x, auto simp add: times_float.simps)
1130 have "Ifloat (?lb x) \<le> ?cos x"
1132 have "-pi \<le> Ifloat ?x2" and "Ifloat ?x2 \<le> pi" unfolding Ifloat_mult Float_num using pi_ge_two `0 \<le> Ifloat x` `Ifloat x \<le> pi` by auto
1133 from lb_half[OF lb_half[OF lb this] `-pi \<le> Ifloat x` `Ifloat x \<le> pi`, unfolded eq_4]
1134 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 .
1136 moreover have "?cos x \<le> Ifloat (?ub x)"
1138 have "-pi \<le> Ifloat ?x2" and "Ifloat ?x2 \<le> pi" unfolding Ifloat_mult Float_num using pi_ge_two `0 \<le> Ifloat x` `Ifloat x \<le> pi` by auto
1139 from ub_half[OF ub_half[OF ub this] `-pi \<le> Ifloat x` `Ifloat x \<le> pi`, unfolded eq_4]
1140 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 .
1142 ultimately show ?thesis by auto
1147 lemma lb_cos_minus: assumes "-pi \<le> Ifloat x" and "Ifloat x \<le> 0"
1148 shows "cos (Ifloat (-x)) \<in> {Ifloat (lb_cos prec (-x)) .. Ifloat (ub_cos prec (-x))}"
1150 have "0 \<le> Ifloat (-x)" and "Ifloat (-x) \<le> pi" using `-pi \<le> Ifloat x` `Ifloat x \<le> 0` by auto
1151 from lb_cos[OF this] show ?thesis .
1154 lemma bnds_cos: "\<forall> x lx ux. (l, u) = bnds_cos prec lx ux \<and> x \<in> {Ifloat lx .. Ifloat ux} \<longrightarrow> Ifloat l \<le> cos x \<and> cos x \<le> Ifloat u"
1155 proof (rule allI, rule allI, rule allI, rule impI)
1157 assume "(l, u) = bnds_cos prec lx ux \<and> x \<in> {Ifloat lx .. Ifloat ux}"
1158 hence bnds: "(l, u) = bnds_cos prec lx ux" and x: "x \<in> {Ifloat lx .. Ifloat ux}" by auto
1160 let ?lpi = "lb_pi prec"
1161 have [intro!]: "Ifloat lx \<le> Ifloat ux" using x by auto
1162 hence "lx \<le> ux" unfolding le_float_def .
1164 show "Ifloat l \<le> cos x \<and> cos x \<le> Ifloat u"
1165 proof (cases "lx < -?lpi \<or> ux > ?lpi")
1167 show ?thesis using bnds unfolding bnds_cos_def if_P[OF True] Let_def using cos_le_one cos_ge_minus_one by auto
1169 case False note not_out = this
1170 hence lpi_lx: "- Ifloat ?lpi \<le> Ifloat lx" and lpi_ux: "Ifloat ux \<le> Ifloat ?lpi" unfolding le_float_def less_float_def by auto
1172 from pi_boundaries[unfolded atLeastAtMost_iff, THEN conjunct1, THEN le_imp_neg_le] lpi_lx
1173 have "- pi \<le> Ifloat lx" by (rule order_trans)
1174 hence "- pi \<le> x" and "- pi \<le> Ifloat ux" and "x \<le> Ifloat ux" using x by auto
1176 from lpi_ux pi_boundaries[unfolded atLeastAtMost_iff, THEN conjunct1]
1177 have "Ifloat ux \<le> pi" by (rule order_trans)
1178 hence "x \<le> pi" and "Ifloat lx \<le> pi" and "Ifloat lx \<le> x" using x by auto
1180 note lb_cos_minus_bottom = lb_cos_minus[unfolded atLeastAtMost_iff, THEN conjunct1]
1181 note lb_cos_minus_top = lb_cos_minus[unfolded atLeastAtMost_iff, THEN conjunct2]
1182 note lb_cos_bottom = lb_cos[unfolded atLeastAtMost_iff, THEN conjunct1]
1183 note lb_cos_top = lb_cos[unfolded atLeastAtMost_iff, THEN conjunct2]
1186 proof (cases "ux \<le> 0")
1187 case True hence "Ifloat ux \<le> 0" unfolding le_float_def by auto
1188 hence "x \<le> 0" and "Ifloat lx \<le> 0" using x by auto
1190 { have "Ifloat (lb_cos prec (-lx)) \<le> cos (Ifloat (-lx))" using lb_cos_minus_bottom[OF `-pi \<le> Ifloat lx` `Ifloat lx \<le> 0`] .
1191 also have "\<dots> \<le> cos x" unfolding Ifloat_minus cos_minus using cos_monotone_minus_pi_0'[OF `- pi \<le> Ifloat lx` `Ifloat lx \<le> x` `x \<le> 0`] .
1192 finally have "Ifloat (lb_cos prec (-lx)) \<le> cos x" . }
1194 { have "cos x \<le> cos (Ifloat (-ux))" unfolding Ifloat_minus cos_minus using cos_monotone_minus_pi_0'[OF `- pi \<le> x` `x \<le> Ifloat ux` `Ifloat ux \<le> 0`] .
1195 also have "\<dots> \<le> Ifloat (ub_cos prec (-ux))" using lb_cos_minus_top[OF `-pi \<le> Ifloat ux` `Ifloat ux \<le> 0`] .
1196 finally have "cos x \<le> Ifloat (ub_cos prec (-ux))" . }
1197 ultimately show ?thesis using bnds unfolding bnds_cos_def Let_def if_not_P[OF not_out] if_P[OF True] by auto
1199 case False note not_ux = this
1202 proof (cases "0 \<le> lx")
1203 case True hence "0 \<le> Ifloat lx" unfolding le_float_def by auto
1204 hence "0 \<le> x" and "0 \<le> Ifloat ux" using x by auto
1206 { have "Ifloat (lb_cos prec ux) \<le> cos (Ifloat ux)" using lb_cos_bottom[OF `0 \<le> Ifloat ux` `Ifloat ux \<le> pi`] .
1207 also have "\<dots> \<le> cos x" using cos_monotone_0_pi'[OF `0 \<le> x` `x \<le> Ifloat ux` `Ifloat ux \<le> pi`] .
1208 finally have "Ifloat (lb_cos prec ux) \<le> cos x" . }
1210 { have "cos x \<le> cos (Ifloat lx)" using cos_monotone_0_pi'[OF `0 \<le> Ifloat lx` `Ifloat lx \<le> x` `x \<le> pi`] .
1211 also have "\<dots> \<le> Ifloat (ub_cos prec lx)" using lb_cos_top[OF `0 \<le> Ifloat lx` `Ifloat lx \<le> pi`] .
1212 finally have "cos x \<le> Ifloat (ub_cos prec lx)" . }
1213 ultimately show ?thesis using bnds unfolding bnds_cos_def Let_def if_not_P[OF not_out] if_not_P[OF not_ux] if_P[OF True] by auto
1215 case False with not_ux
1216 have "Ifloat lx \<le> 0" and "0 \<le> Ifloat ux" unfolding le_float_def by auto
1218 have "Ifloat (min (lb_cos prec (-lx)) (lb_cos prec ux)) \<le> cos x"
1219 proof (cases "x \<le> 0")
1221 have "Ifloat (lb_cos prec (-lx)) \<le> cos (Ifloat (-lx))" using lb_cos_minus_bottom[OF `-pi \<le> Ifloat lx` `Ifloat lx \<le> 0`] .
1222 also have "\<dots> \<le> cos x" unfolding Ifloat_minus cos_minus using cos_monotone_minus_pi_0'[OF `- pi \<le> Ifloat lx` `Ifloat lx \<le> x` `x \<le> 0`] .
1223 finally show ?thesis unfolding Ifloat_min by auto
1225 case False hence "0 \<le> x" by auto
1226 have "Ifloat (lb_cos prec ux) \<le> cos (Ifloat ux)" using lb_cos_bottom[OF `0 \<le> Ifloat ux` `Ifloat ux \<le> pi`] .
1227 also have "\<dots> \<le> cos x" using cos_monotone_0_pi'[OF `0 \<le> x` `x \<le> Ifloat ux` `Ifloat ux \<le> pi`] .
1228 finally show ?thesis unfolding Ifloat_min by auto
1230 moreover have "cos x \<le> Ifloat (Float 1 0)" by auto
1231 ultimately show ?thesis using bnds unfolding bnds_cos_def Let_def if_not_P[OF not_out] if_not_P[OF not_ux] if_not_P[OF False] by auto
1237 subsection "Compute the sinus in the entire domain"
1239 function lb_sin :: "nat \<Rightarrow> float \<Rightarrow> float" and ub_sin :: "nat \<Rightarrow> float \<Rightarrow> float" where
1240 "lb_sin prec x = (let sqr_diff = \<lambda> x. if x > 1 then 0 else 1 - x * x
1241 in if x < 0 then - ub_sin prec (- x)
1242 else if x \<le> Float 1 -1 then x * lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 2 1 (x * x)
1243 else the (lb_sqrt prec (sqr_diff (ub_cos prec x))))" |
1245 "ub_sin prec x = (let sqr_diff = \<lambda> x. if x < 0 then 1 else 1 - x * x
1246 in if x < 0 then - lb_sin prec (- x)
1247 else if x \<le> Float 1 -1 then x * ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 2 1 (x * x)
1248 else the (ub_sqrt prec (sqr_diff (lb_cos prec x))))"
1249 by pat_completeness auto
1250 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)
1252 definition bnds_sin :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where
1253 "bnds_sin prec lx ux = (let
1255 half_pi = lpi * Float 1 -1
1256 in if lx \<le> - half_pi \<or> half_pi \<le> ux then (Float -1 0, Float 1 0)
1257 else (lb_sin prec lx, ub_sin prec ux))"
1259 lemma lb_sin: assumes "- (pi / 2) \<le> Ifloat x" and "Ifloat x \<le> pi / 2"
1260 shows "sin (Ifloat x) \<in> { Ifloat (lb_sin prec x) .. Ifloat (ub_sin prec x) }" (is "?sin x \<in> { ?lb x .. ?ub x}")
1262 { fix x :: float assume "0 \<le> Ifloat x" and "Ifloat x \<le> pi / 2"
1263 hence "\<not> (x < 0)" and "- (pi / 2) \<le> Ifloat x" unfolding less_float_def using pi_ge_two by auto
1265 have "Ifloat x \<le> pi" using `Ifloat x \<le> pi / 2` using pi_ge_two by auto
1267 have "?sin x \<in> { ?lb x .. ?ub x}"
1268 proof (cases "x \<le> Float 1 -1")
1269 case True from sin_boundaries[OF `0 \<le> Ifloat x` `Ifloat x \<le> pi / 2`]
1270 show ?thesis unfolding lb_sin.simps[of prec x] ub_sin.simps[of prec x] if_not_P[OF `\<not> (x < 0)`] if_P[OF True] Let_def .
1273 have "0 \<le> cos (Ifloat x)" using cos_ge_zero[OF _ `Ifloat x \<le> pi /2`] `0 \<le> Ifloat x` pi_ge_two by auto
1274 have "0 \<le> sin (Ifloat x)" using `0 \<le> Ifloat x` and `Ifloat x \<le> pi / 2` using sin_ge_zero by auto
1276 have "?sin x \<le> ?ub x"
1277 proof (cases "lb_cos prec x < 0")
1279 have "?sin x \<le> 1" using sin_le_one .
1280 also have "\<dots> \<le> Ifloat (the (ub_sqrt prec 1))" using ub_sqrt_lower_bound[where prec=prec and x=1] unfolding Ifloat_1 by auto
1281 finally show ?thesis unfolding ub_sin.simps if_not_P[OF `\<not> (x < 0)`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_P[OF True] Let_def .
1283 case False hence "0 \<le> Ifloat (lb_cos prec x)" unfolding less_float_def by auto
1285 have "sin (Ifloat x) = sqrt (1 - cos (Ifloat x) ^ 2)" unfolding sin_squared_eq[symmetric] real_sqrt_abs using `0 \<le> sin (Ifloat x)` by auto
1286 also have "\<dots> \<le> sqrt (Ifloat (1 - lb_cos prec x * lb_cos prec x))"
1287 proof (rule real_sqrt_le_mono)
1288 have "Ifloat (lb_cos prec x * lb_cos prec x) \<le> cos (Ifloat x) ^ 2" unfolding numeral_2_eq_2 power_Suc2 power_0 Ifloat_mult
1289 using `0 \<le> Ifloat (lb_cos prec x)` lb_cos[OF `0 \<le> Ifloat x` `Ifloat x \<le> pi`] `0 \<le> cos (Ifloat x)` by(auto intro!: mult_mono)
1290 thus "1 - cos (Ifloat x) ^ 2 \<le> Ifloat (1 - lb_cos prec x * lb_cos prec x)" unfolding Ifloat_sub Ifloat_1 by auto
1292 also have "\<dots> \<le> Ifloat (the (ub_sqrt prec (1 - lb_cos prec x * lb_cos prec x)))"
1293 proof (rule ub_sqrt_lower_bound)
1294 have "Ifloat (lb_cos prec x) \<le> cos (Ifloat x)" using lb_cos[OF `0 \<le> Ifloat x` `Ifloat x \<le> pi`] by auto
1295 from mult_mono[OF order_trans[OF this cos_le_one] order_trans[OF this cos_le_one]]
1296 have "Ifloat (lb_cos prec x) * Ifloat (lb_cos prec x) \<le> 1" using `0 \<le> Ifloat (lb_cos prec x)` by auto
1297 thus "0 \<le> Ifloat (1 - lb_cos prec x * lb_cos prec x)" by auto
1299 finally show ?thesis unfolding ub_sin.simps if_not_P[OF `\<not> (x < 0)`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_not_P[OF False] Let_def .
1302 have "?lb x \<le> ?sin x"
1303 proof (cases "1 < ub_cos prec x")
1305 show ?thesis unfolding lb_sin.simps if_not_P[OF `\<not> (x < 0)`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_P[OF True] Let_def
1306 by (rule order_trans[OF _ sin_ge_zero[OF `0 \<le> Ifloat x` `Ifloat x \<le> pi`]])
1307 (auto simp add: lb_sqrt_upper_bound[where prec=prec and x=0, unfolded Ifloat_0 real_sqrt_zero])
1309 case False hence "Ifloat (ub_cos prec x) \<le> 1" unfolding less_float_def by auto
1310 have "0 \<le> Ifloat (ub_cos prec x)" using order_trans[OF `0 \<le> cos (Ifloat x)`] lb_cos `0 \<le> Ifloat x` `Ifloat x \<le> pi` by auto
1312 have "Ifloat (the (lb_sqrt prec (1 - ub_cos prec x * ub_cos prec x))) \<le> sqrt (Ifloat (1 - ub_cos prec x * ub_cos prec x))"
1313 proof (rule lb_sqrt_upper_bound)
1314 from mult_mono[OF `Ifloat (ub_cos prec x) \<le> 1` `Ifloat (ub_cos prec x) \<le> 1`] `0 \<le> Ifloat (ub_cos prec x)`
1315 have "Ifloat (ub_cos prec x) * Ifloat (ub_cos prec x) \<le> 1" by auto
1316 thus "0 \<le> Ifloat (1 - ub_cos prec x * ub_cos prec x)" by auto
1318 also have "\<dots> \<le> sqrt (1 - cos (Ifloat x) ^ 2)"
1319 proof (rule real_sqrt_le_mono)
1320 have "cos (Ifloat x) ^ 2 \<le> Ifloat (ub_cos prec x * ub_cos prec x)" unfolding numeral_2_eq_2 power_Suc2 power_0 Ifloat_mult
1321 using `0 \<le> Ifloat (ub_cos prec x)` lb_cos[OF `0 \<le> Ifloat x` `Ifloat x \<le> pi`] `0 \<le> cos (Ifloat x)` by(auto intro!: mult_mono)
1322 thus "Ifloat (1 - ub_cos prec x * ub_cos prec x) \<le> 1 - cos (Ifloat x) ^ 2" unfolding Ifloat_sub Ifloat_1 by auto
1324 also have "\<dots> = sin (Ifloat x)" unfolding sin_squared_eq[symmetric] real_sqrt_abs using `0 \<le> sin (Ifloat x)` by auto
1325 finally show ?thesis unfolding lb_sin.simps if_not_P[OF `\<not> (x < 0)`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_not_P[OF False] Let_def .
1327 ultimately show ?thesis by auto
1329 } note for_pos = this
1332 proof (cases "x < 0")
1334 hence "0 \<le> Ifloat (-x)" and "Ifloat (- x) \<le> pi / 2" using `-(pi/2) \<le> Ifloat x` unfolding less_float_def by auto
1335 from for_pos[OF this]
1336 show ?thesis unfolding Ifloat_minus sin_minus lb_sin.simps[of prec x] ub_sin.simps[of prec x] if_P[OF True] Let_def atLeastAtMost_iff by auto
1338 case False hence "0 \<le> Ifloat x" unfolding less_float_def by auto
1339 from for_pos[OF this `Ifloat x \<le> pi /2`]
1344 lemma bnds_sin: "\<forall> x lx ux. (l, u) = bnds_sin prec lx ux \<and> x \<in> {Ifloat lx .. Ifloat ux} \<longrightarrow> Ifloat l \<le> sin x \<and> sin x \<le> Ifloat u"
1345 proof (rule allI, rule allI, rule allI, rule impI)
1347 assume "(l, u) = bnds_sin prec lx ux \<and> x \<in> {Ifloat lx .. Ifloat ux}"
1348 hence bnds: "(l, u) = bnds_sin prec lx ux" and x: "x \<in> {Ifloat lx .. Ifloat ux}" by auto
1349 show "Ifloat l \<le> sin x \<and> sin x \<le> Ifloat u"
1350 proof (cases "lx \<le> - (lb_pi prec * Float 1 -1) \<or> lb_pi prec * Float 1 -1 \<le> ux")
1351 case True show ?thesis using bnds unfolding bnds_sin_def if_P[OF True] Let_def by auto
1354 hence "- lb_pi prec * Float 1 -1 \<le> lx" and "ux \<le> lb_pi prec * Float 1 -1" unfolding le_float_def by auto
1355 moreover have "Ifloat (lb_pi prec * Float 1 -1) \<le> pi / 2" unfolding Ifloat_mult using pi_boundaries by auto
1356 ultimately have "- (pi / 2) \<le> Ifloat lx" and "Ifloat ux \<le> pi / 2" and "Ifloat lx \<le> Ifloat ux" unfolding le_float_def using x by auto
1357 hence "- (pi / 2) \<le> Ifloat ux" and "Ifloat lx \<le> pi / 2" by auto
1359 have "- (pi / 2) \<le> x""x \<le> pi / 2" using `Ifloat ux \<le> pi / 2` `- (pi /2) \<le> Ifloat lx` x by auto
1361 { have "Ifloat (lb_sin prec lx) \<le> sin (Ifloat lx)" using lb_sin[OF `- (pi / 2) \<le> Ifloat lx` `Ifloat lx \<le> pi / 2`] unfolding atLeastAtMost_iff by auto
1362 also have "\<dots> \<le> sin x" using sin_monotone_2pi' `- (pi / 2) \<le> Ifloat lx` x `x \<le> pi / 2` by auto
1363 finally have "Ifloat (lb_sin prec lx) \<le> sin x" . }
1365 { have "sin x \<le> sin (Ifloat ux)" using sin_monotone_2pi' `- (pi / 2) \<le> x` x `Ifloat ux \<le> pi / 2` by auto
1366 also have "\<dots> \<le> Ifloat (ub_sin prec ux)" using lb_sin[OF `- (pi / 2) \<le> Ifloat ux` `Ifloat ux \<le> pi / 2`] unfolding atLeastAtMost_iff by auto
1367 finally have "sin x \<le> Ifloat (ub_sin prec ux)" . }
1369 show ?thesis using bnds unfolding bnds_sin_def if_not_P[OF False] Let_def by auto
1373 section "Exponential function"
1375 subsection "Compute the series of the exponential function"
1377 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
1378 "ub_exp_horner prec 0 i k x = 0" |
1379 "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" |
1380 "lb_exp_horner prec 0 i k x = 0" |
1381 "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"
1383 lemma bnds_exp_horner: assumes "Ifloat x \<le> 0"
1384 shows "exp (Ifloat x) \<in> { Ifloat (lb_exp_horner prec (get_even n) 1 1 x) .. Ifloat (ub_exp_horner prec (get_odd n) 1 1 x) }"
1387 have F: "\<And> m. ((\<lambda>i. i + 1) ^ n) m = n + m" by (induct n, auto)
1388 have "fact (Suc n) = fact n * ((\<lambda>i. i + 1) ^ n) 1" unfolding F by auto } note f_eq = this
1390 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,
1391 OF assms f_eq lb_exp_horner.simps ub_exp_horner.simps]
1393 { have "Ifloat (lb_exp_horner prec (get_even n) 1 1 x) \<le> (\<Sum>j = 0..<get_even n. 1 / real (fact j) * Ifloat x ^ j)"
1394 using bounds(1) by auto
1395 also have "\<dots> \<le> exp (Ifloat x)"
1397 obtain t where "\<bar>t\<bar> \<le> \<bar>Ifloat x\<bar>" and "exp (Ifloat x) = (\<Sum>m = 0..<get_even n. (Ifloat x) ^ m / real (fact m)) + exp t / real (fact (get_even n)) * (Ifloat x) ^ (get_even n)"
1398 using Maclaurin_exp_le by blast
1399 moreover have "0 \<le> exp t / real (fact (get_even n)) * (Ifloat x) ^ (get_even n)"
1400 by (auto intro!: mult_nonneg_nonneg divide_nonneg_pos simp add: get_even zero_le_even_power exp_gt_zero)
1401 ultimately show ?thesis
1402 using get_odd exp_gt_zero by (auto intro!: pordered_cancel_semiring_class.mult_nonneg_nonneg)
1404 finally have "Ifloat (lb_exp_horner prec (get_even n) 1 1 x) \<le> exp (Ifloat x)" .
1407 have x_less_zero: "Ifloat x ^ get_odd n \<le> 0"
1408 proof (cases "Ifloat x = 0")
1410 have "(get_odd n) \<noteq> 0" using get_odd[THEN odd_pos] by auto
1411 thus ?thesis unfolding True power_0_left by auto
1413 case False hence "Ifloat x < 0" using `Ifloat x \<le> 0` by auto
1414 show ?thesis by (rule less_imp_le, auto simp add: power_less_zero_eq get_odd `Ifloat x < 0`)
1417 obtain t where "\<bar>t\<bar> \<le> \<bar>Ifloat x\<bar>" and "exp (Ifloat x) = (\<Sum>m = 0..<get_odd n. (Ifloat x) ^ m / real (fact m)) + exp t / real (fact (get_odd n)) * (Ifloat x) ^ (get_odd n)"
1418 using Maclaurin_exp_le by blast
1419 moreover have "exp t / real (fact (get_odd n)) * (Ifloat x) ^ (get_odd n) \<le> 0"
1420 by (auto intro!: mult_nonneg_nonpos divide_nonpos_pos simp add: x_less_zero exp_gt_zero)
1421 ultimately have "exp (Ifloat x) \<le> (\<Sum>j = 0..<get_odd n. 1 / real (fact j) * Ifloat x ^ j)"
1422 using get_odd exp_gt_zero by (auto intro!: pordered_cancel_semiring_class.mult_nonneg_nonneg)
1423 also have "\<dots> \<le> Ifloat (ub_exp_horner prec (get_odd n) 1 1 x)"
1424 using bounds(2) by auto
1425 finally have "exp (Ifloat x) \<le> Ifloat (ub_exp_horner prec (get_odd n) 1 1 x)" .
1426 } ultimately show ?thesis by auto
1429 subsection "Compute the exponential function on the entire domain"
1431 function ub_exp :: "nat \<Rightarrow> float \<Rightarrow> float" and lb_exp :: "nat \<Rightarrow> float \<Rightarrow> float" where
1432 "lb_exp prec x = (if 0 < x then float_divl prec 1 (ub_exp prec (-x))
1434 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)
1435 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))
1437 "ub_exp prec x = (if 0 < x then float_divr prec 1 (lb_exp prec (-x))
1438 else if x < - 1 then (case floor_fl x of (Float m e) \<Rightarrow>
1439 (ub_exp_horner prec (get_odd (prec + 2)) 1 1 (float_divr prec x (- Float m e))) ^ (nat (-m) * 2 ^ nat e))
1440 else ub_exp_horner prec (get_odd (prec + 2)) 1 1 x)"
1441 by pat_completeness auto
1442 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)
1444 lemma exp_m1_ge_quarter: "(1 / 4 :: real) \<le> exp (- 1)"
1446 have eq4: "4 = Suc (Suc (Suc (Suc 0)))" by auto
1448 have "1 / 4 = Ifloat (Float 1 -2)" unfolding Float_num by auto
1449 also have "\<dots> \<le> Ifloat (lb_exp_horner 1 (get_even 4) 1 1 (- 1))"
1450 unfolding get_even_def eq4
1451 by (auto simp add: lapprox_posrat_def rapprox_posrat_def normfloat.simps)
1452 also have "\<dots> \<le> exp (Ifloat (- 1))" using bnds_exp_horner[where x="- 1"] by auto
1453 finally show ?thesis unfolding Ifloat_minus Ifloat_1 .
1456 lemma lb_exp_pos: assumes "\<not> 0 < x" shows "0 < lb_exp prec x"
1458 let "?lb_horner x" = "lb_exp_horner prec (get_even (prec + 2)) 1 1 x"
1459 let "?horner x" = "let y = ?lb_horner x in if y \<le> 0 then Float 1 -2 else y"
1460 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)
1461 moreover { fix x :: float fix num :: nat
1462 have "0 < Ifloat (?horner x) ^ num" using `0 < ?horner x`[unfolded less_float_def Ifloat_0] by (rule zero_less_power)
1463 also have "\<dots> = Ifloat ((?horner x) ^ num)" using float_power by auto
1464 finally have "0 < Ifloat ((?horner x) ^ num)" .
1466 ultimately show ?thesis
1467 unfolding lb_exp.simps if_not_P[OF `\<not> 0 < x`] Let_def by (cases "floor_fl x", cases "x < - 1", auto simp add: le_float_def less_float_def normfloat)
1470 lemma exp_boundaries': assumes "x \<le> 0"
1471 shows "exp (Ifloat x) \<in> { Ifloat (lb_exp prec x) .. Ifloat (ub_exp prec x)}"
1473 let "?lb_exp_horner x" = "lb_exp_horner prec (get_even (prec + 2)) 1 1 x"
1474 let "?ub_exp_horner x" = "ub_exp_horner prec (get_odd (prec + 2)) 1 1 x"
1476 have "Ifloat x \<le> 0" and "\<not> x > 0" using `x \<le> 0` unfolding le_float_def less_float_def by auto
1478 proof (cases "x < - 1")
1479 case False hence "- 1 \<le> Ifloat x" unfolding less_float_def by auto
1481 proof (cases "?lb_exp_horner x \<le> 0")
1482 from `\<not> x < - 1` have "- 1 \<le> Ifloat x" unfolding less_float_def by auto
1483 hence "exp (- 1) \<le> exp (Ifloat x)" unfolding exp_le_cancel_iff .
1484 from order_trans[OF exp_m1_ge_quarter this]
1485 have "Ifloat (Float 1 -2) \<le> exp (Ifloat x)" unfolding Float_num .
1487 ultimately show ?thesis using bnds_exp_horner `Ifloat x \<le> 0` `\<not> x > 0` `\<not> x < - 1` by auto
1489 case False thus ?thesis using bnds_exp_horner `Ifloat x \<le> 0` `\<not> x > 0` `\<not> x < - 1` by (auto simp add: Let_def)
1494 obtain m e where Float_floor: "floor_fl x = Float m e" by (cases "floor_fl x", auto)
1495 let ?num = "nat (- m) * 2 ^ nat e"
1497 have "Ifloat (floor_fl x) < - 1" using floor_fl `x < - 1` unfolding le_float_def less_float_def Ifloat_minus Ifloat_1 by (rule order_le_less_trans)
1498 hence "Ifloat (floor_fl x) < 0" unfolding Float_floor Ifloat.simps using zero_less_pow2[of xe] by auto
1500 unfolding less_float_def Ifloat_0 Float_floor Ifloat.simps
1501 unfolding pos_prod_lt[OF zero_less_pow2[of e], unfolded real_mult_commute] by auto
1502 hence "1 \<le> - m" by auto
1503 hence "0 < nat (- m)" by auto
1505 have "0 \<le> e" using floor_pos_exp Float_floor[symmetric] by auto
1506 hence "(0::nat) < 2 ^ nat e" by auto
1507 ultimately have "0 < ?num" by auto
1508 hence "real ?num \<noteq> 0" by auto
1509 have e_nat: "int (nat e) = e" using `0 \<le> e` by auto
1510 have num_eq: "real ?num = Ifloat (- floor_fl x)" using `0 < nat (- m)`
1511 unfolding Float_floor Ifloat_minus Ifloat.simps real_of_nat_mult pow2_int[of "nat e", unfolded e_nat] realpow_real_of_nat[symmetric] by auto
1512 have "0 < - floor_fl x" using `0 < ?num`[unfolded real_of_nat_less_iff[symmetric]] unfolding less_float_def num_eq[symmetric] Ifloat_0 real_of_nat_zero .
1513 hence "Ifloat (floor_fl x) < 0" unfolding less_float_def by auto
1515 have "exp (Ifloat x) \<le> Ifloat (ub_exp prec x)"
1517 have div_less_zero: "Ifloat (float_divr prec x (- floor_fl x)) \<le> 0"
1518 using float_divr_nonpos_pos_upper_bound[OF `x \<le> 0` `0 < - floor_fl x`] unfolding le_float_def Ifloat_0 .
1520 have "exp (Ifloat x) = exp (real ?num * (Ifloat x / real ?num))" using `real ?num \<noteq> 0` by auto
1521 also have "\<dots> = exp (Ifloat x / real ?num) ^ ?num" unfolding exp_real_of_nat_mult ..
1522 also have "\<dots> \<le> exp (Ifloat (float_divr prec x (- floor_fl x))) ^ ?num" unfolding num_eq
1523 by (rule power_mono, rule exp_le_cancel_iff[THEN iffD2], rule float_divr) auto
1524 also have "\<dots> \<le> Ifloat ((?ub_exp_horner (float_divr prec x (- floor_fl x))) ^ ?num)" unfolding float_power
1525 by (rule power_mono, rule bnds_exp_horner[OF div_less_zero, unfolded atLeastAtMost_iff, THEN conjunct2], auto)
1526 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 .
1529 have "Ifloat (lb_exp prec x) \<le> exp (Ifloat x)"
1531 let ?divl = "float_divl prec x (- Float m e)"
1532 let ?horner = "?lb_exp_horner ?divl"
1535 proof (cases "?horner \<le> 0")
1536 case False hence "0 \<le> Ifloat ?horner" unfolding le_float_def by auto
1538 have div_less_zero: "Ifloat (float_divl prec x (- floor_fl x)) \<le> 0"
1539 using `Ifloat (floor_fl x) < 0` `Ifloat x \<le> 0` by (auto intro!: order_trans[OF float_divl] divide_nonpos_neg)
1541 have "Ifloat ((?lb_exp_horner (float_divl prec x (- floor_fl x))) ^ ?num) \<le>
1542 exp (Ifloat (float_divl prec x (- floor_fl x))) ^ ?num" unfolding float_power
1543 using `0 \<le> Ifloat ?horner`[unfolded Float_floor[symmetric]] bnds_exp_horner[OF div_less_zero, unfolded atLeastAtMost_iff, THEN conjunct1] by (auto intro!: power_mono)
1544 also have "\<dots> \<le> exp (Ifloat x / real ?num) ^ ?num" unfolding num_eq
1545 using float_divl by (auto intro!: power_mono simp del: Ifloat_minus)
1546 also have "\<dots> = exp (real ?num * (Ifloat x / real ?num))" unfolding exp_real_of_nat_mult ..
1547 also have "\<dots> = exp (Ifloat x)" using `real ?num \<noteq> 0` by auto
1548 finally show ?thesis
1549 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
1552 have "Ifloat (floor_fl x) \<noteq> 0" and "Ifloat (floor_fl x) \<le> 0" using `Ifloat (floor_fl x) < 0` by auto
1553 from divide_right_mono_neg[OF floor_fl[of x] `Ifloat (floor_fl x) \<le> 0`, unfolded divide_self[OF `Ifloat (floor_fl x) \<noteq> 0`]]
1554 have "- 1 \<le> Ifloat x / Ifloat (- floor_fl x)" unfolding Ifloat_minus by auto
1555 from order_trans[OF exp_m1_ge_quarter this[unfolded exp_le_cancel_iff[where x="- 1", symmetric]]]
1556 have "Ifloat (Float 1 -2) \<le> exp (Ifloat x / Ifloat (- floor_fl x))" unfolding Float_num .
1557 hence "Ifloat (Float 1 -2) ^ ?num \<le> exp (Ifloat x / Ifloat (- floor_fl x)) ^ ?num"
1558 by (auto intro!: power_mono simp add: Float_num)
1559 also have "\<dots> = exp (Ifloat x)" unfolding num_eq exp_real_of_nat_mult[symmetric] using `Ifloat (floor_fl x) \<noteq> 0` by auto
1560 finally show ?thesis
1561 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 .
1564 ultimately show ?thesis by auto
1568 lemma exp_boundaries: "exp (Ifloat x) \<in> { Ifloat (lb_exp prec x) .. Ifloat (ub_exp prec x)}"
1571 proof (cases "0 < x")
1572 case False hence "x \<le> 0" unfolding less_float_def le_float_def by auto
1573 from exp_boundaries'[OF this] show ?thesis .
1575 case True hence "-x \<le> 0" unfolding less_float_def le_float_def by auto
1577 have "Ifloat (lb_exp prec x) \<le> exp (Ifloat x)"
1579 from exp_boundaries'[OF `-x \<le> 0`]
1580 have ub_exp: "exp (- Ifloat x) \<le> Ifloat (ub_exp prec (-x))" unfolding atLeastAtMost_iff Ifloat_minus by auto
1582 have "Ifloat (float_divl prec 1 (ub_exp prec (-x))) \<le> Ifloat 1 / Ifloat (ub_exp prec (-x))" using float_divl .
1583 also have "Ifloat 1 / Ifloat (ub_exp prec (-x)) \<le> exp (Ifloat x)"
1584 using ub_exp[unfolded inverse_le_iff_le[OF order_less_le_trans[OF exp_gt_zero ub_exp] exp_gt_zero, symmetric]]
1585 unfolding exp_minus nonzero_inverse_inverse_eq[OF exp_not_eq_zero] inverse_eq_divide by auto
1586 finally show ?thesis unfolding lb_exp.simps if_P[OF True] .
1589 have "exp (Ifloat x) \<le> Ifloat (ub_exp prec x)"
1591 have "\<not> 0 < -x" using `0 < x` unfolding less_float_def by auto
1593 from exp_boundaries'[OF `-x \<le> 0`]
1594 have lb_exp: "Ifloat (lb_exp prec (-x)) \<le> exp (- Ifloat x)" unfolding atLeastAtMost_iff Ifloat_minus by auto
1596 have "exp (Ifloat x) \<le> Ifloat 1 / Ifloat (lb_exp prec (-x))"
1597 using lb_exp[unfolded inverse_le_iff_le[OF exp_gt_zero lb_exp_pos[OF `\<not> 0 < -x`, unfolded less_float_def Ifloat_0], symmetric]]
1598 unfolding exp_minus nonzero_inverse_inverse_eq[OF exp_not_eq_zero] inverse_eq_divide Ifloat_1 by auto
1599 also have "\<dots> \<le> Ifloat (float_divr prec 1 (lb_exp prec (-x)))" using float_divr .
1600 finally show ?thesis unfolding ub_exp.simps if_P[OF True] .
1602 ultimately show ?thesis by auto
1606 lemma bnds_exp: "\<forall> x lx ux. (l, u) = (lb_exp prec lx, ub_exp prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux} \<longrightarrow> Ifloat l \<le> exp x \<and> exp x \<le> Ifloat u"
1607 proof (rule allI, rule allI, rule allI, rule impI)
1609 assume "(l, u) = (lb_exp prec lx, ub_exp prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux}"
1610 hence l: "lb_exp prec lx = l " and u: "ub_exp prec ux = u" and x: "x \<in> {Ifloat lx .. Ifloat ux}" by auto
1612 { from exp_boundaries[of lx prec, unfolded l]
1613 have "Ifloat l \<le> exp (Ifloat lx)" by (auto simp del: lb_exp.simps)
1614 also have "\<dots> \<le> exp x" using x by auto
1615 finally have "Ifloat l \<le> exp x" .
1617 { have "exp x \<le> exp (Ifloat ux)" using x by auto
1618 also have "\<dots> \<le> Ifloat u" using exp_boundaries[of ux prec, unfolded u] by (auto simp del: ub_exp.simps)
1619 finally have "exp x \<le> Ifloat u" .
1620 } ultimately show "Ifloat l \<le> exp x \<and> exp x \<le> Ifloat u" ..
1625 subsection "Compute the logarithm series"
1627 fun ub_ln_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
1628 and lb_ln_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
1629 "ub_ln_horner prec 0 i x = 0" |
1630 "ub_ln_horner prec (Suc n) i x = rapprox_rat prec 1 (int i) - x * lb_ln_horner prec n (Suc i) x" |
1631 "lb_ln_horner prec 0 i x = 0" |
1632 "lb_ln_horner prec (Suc n) i x = lapprox_rat prec 1 (int i) - x * ub_ln_horner prec n (Suc i) x"
1635 assumes "0 \<le> x" and "x < 1"
1636 shows "(\<Sum>i=0..<2*n. -1^i * (1 / real (i + 1)) * x^(Suc i)) \<le> ln (x + 1)" (is "?lb")
1637 and "ln (x + 1) \<le> (\<Sum>i=0..<2*n + 1. -1^i * (1 / real (i + 1)) * x^(Suc i))" (is "?ub")
1639 let "?a n" = "(1/real (n +1)) * x^(Suc n)"
1641 have ln_eq: "(\<Sum> i. -1^i * ?a i) = ln (x + 1)"
1642 using ln_series[of "x + 1"] `0 \<le> x` `x < 1` by auto
1644 have "norm x < 1" using assms by auto
1645 have "?a ----> 0" unfolding Suc_plus1[symmetric] inverse_eq_divide[symmetric]
1646 using LIMSEQ_mult[OF LIMSEQ_inverse_real_of_nat LIMSEQ_Suc[OF LIMSEQ_power_zero[OF `norm x < 1`]]] by auto
1647 { fix n have "0 \<le> ?a n" by (rule mult_nonneg_nonneg, auto intro!: mult_nonneg_nonneg simp add: `0 \<le> x`) }
1648 { fix n have "?a (Suc n) \<le> ?a n" unfolding inverse_eq_divide[symmetric]
1649 proof (rule mult_mono)
1650 show "0 \<le> x ^ Suc (Suc n)" by (auto intro!: mult_nonneg_nonneg simp add: `0 \<le> x`)
1651 have "x ^ Suc (Suc n) \<le> x ^ Suc n * 1" unfolding power_Suc2 real_mult_assoc[symmetric]
1652 by (rule mult_left_mono, fact less_imp_le[OF `x < 1`], auto intro!: mult_nonneg_nonneg simp add: `0 \<le> x`)
1653 thus "x ^ Suc (Suc n) \<le> x ^ Suc n" by auto
1655 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]
1656 show "?lb" and "?ub" by auto
1659 lemma ln_float_bounds:
1660 assumes "0 \<le> Ifloat x" and "Ifloat x < 1"
1661 shows "Ifloat (x * lb_ln_horner prec (get_even n) 1 x) \<le> ln (Ifloat x + 1)" (is "?lb \<le> ?ln")
1662 and "ln (Ifloat x + 1) \<le> Ifloat (x * ub_ln_horner prec (get_odd n) 1 x)" (is "?ln \<le> ?ub")
1664 obtain ev where ev: "get_even n = 2 * ev" using get_even_double ..
1665 obtain od where od: "get_odd n = 2 * od + 1" using get_odd_double ..
1667 let "?s n" = "-1^n * (1 / real (1 + n)) * (Ifloat x)^(Suc n)"
1669 have "?lb \<le> setsum ?s {0 ..< 2 * ev}" unfolding power_Suc2 real_mult_assoc[symmetric] Ifloat_mult setsum_left_distrib[symmetric] unfolding real_mult_commute[of "Ifloat x"] ev
1670 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",
1671 OF `0 \<le> Ifloat x` refl lb_ln_horner.simps ub_ln_horner.simps] `0 \<le> Ifloat x`
1672 by (rule mult_right_mono)
1673 also have "\<dots> \<le> ?ln" using ln_bounds(1)[OF `0 \<le> Ifloat x` `Ifloat x < 1`] by auto
1674 finally show "?lb \<le> ?ln" .
1676 have "?ln \<le> setsum ?s {0 ..< 2 * od + 1}" using ln_bounds(2)[OF `0 \<le> Ifloat x` `Ifloat x < 1`] by auto
1677 also have "\<dots> \<le> ?ub" unfolding power_Suc2 real_mult_assoc[symmetric] Ifloat_mult setsum_left_distrib[symmetric] unfolding real_mult_commute[of "Ifloat x"] od
1678 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",
1679 OF `0 \<le> Ifloat x` refl lb_ln_horner.simps ub_ln_horner.simps] `0 \<le> Ifloat x`
1680 by (rule mult_right_mono)
1681 finally show "?ln \<le> ?ub" .
1684 lemma ln_add: assumes "0 < x" and "0 < y" shows "ln (x + y) = ln x + ln (1 + y / x)"
1686 have "x \<noteq> 0" using assms by auto
1687 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
1689 have "0 < y / x" using assms divide_pos_pos by auto
1690 hence "0 < 1 + y / x" by auto
1691 ultimately show ?thesis using ln_mult assms by auto
1694 subsection "Compute the logarithm of 2"
1696 definition ub_ln2 where "ub_ln2 prec = (let third = rapprox_rat (max prec 1) 1 3
1697 in (Float 1 -1 * ub_ln_horner prec (get_odd prec) 1 (Float 1 -1)) +
1698 (third * ub_ln_horner prec (get_odd prec) 1 third))"
1699 definition lb_ln2 where "lb_ln2 prec = (let third = lapprox_rat prec 1 3
1700 in (Float 1 -1 * lb_ln_horner prec (get_even prec) 1 (Float 1 -1)) +
1701 (third * lb_ln_horner prec (get_even prec) 1 third))"
1703 lemma ub_ln2: "ln 2 \<le> Ifloat (ub_ln2 prec)" (is "?ub_ln2")
1704 and lb_ln2: "Ifloat (lb_ln2 prec) \<le> ln 2" (is "?lb_ln2")
1706 let ?uthird = "rapprox_rat (max prec 1) 1 3"
1707 let ?lthird = "lapprox_rat prec 1 3"
1709 have ln2_sum: "ln 2 = ln (1/2 + 1) + ln (1 / 3 + 1)"
1710 using ln_add[of "3 / 2" "1 / 2"] by auto
1711 have lb3: "Ifloat ?lthird \<le> 1 / 3" using lapprox_rat[of prec 1 3] by auto
1712 hence lb3_ub: "Ifloat ?lthird < 1" by auto
1713 have lb3_lb: "0 \<le> Ifloat ?lthird" using lapprox_rat_bottom[of 1 3] by auto
1714 have ub3: "1 / 3 \<le> Ifloat ?uthird" using rapprox_rat[of 1 3] by auto
1715 hence ub3_lb: "0 \<le> Ifloat ?uthird" by auto
1717 have lb2: "0 \<le> Ifloat (Float 1 -1)" and ub2: "Ifloat (Float 1 -1) < 1" unfolding Float_num by auto
1719 have "0 \<le> (1::int)" and "0 < (3::int)" by auto
1720 have ub3_ub: "Ifloat ?uthird < 1" unfolding rapprox_rat.simps(2)[OF `0 \<le> 1` `0 < 3`]
1721 by (rule rapprox_posrat_less1, auto)
1723 have third_gt0: "(0 :: real) < 1 / 3 + 1" by auto
1724 have uthird_gt0: "0 < Ifloat ?uthird + 1" using ub3_lb by auto
1725 have lthird_gt0: "0 < Ifloat ?lthird + 1" using lb3_lb by auto
1727 show ?ub_ln2 unfolding ub_ln2_def Let_def Ifloat_add ln2_sum Float_num(4)[symmetric]
1728 proof (rule add_mono, fact ln_float_bounds(2)[OF lb2 ub2])
1729 have "ln (1 / 3 + 1) \<le> ln (Ifloat ?uthird + 1)" unfolding ln_le_cancel_iff[OF third_gt0 uthird_gt0] using ub3 by auto
1730 also have "\<dots> \<le> Ifloat (?uthird * ub_ln_horner prec (get_odd prec) 1 ?uthird)"
1731 using ln_float_bounds(2)[OF ub3_lb ub3_ub] .
1732 finally show "ln (1 / 3 + 1) \<le> Ifloat (?uthird * ub_ln_horner prec (get_odd prec) 1 ?uthird)" .
1734 show ?lb_ln2 unfolding lb_ln2_def Let_def Ifloat_add ln2_sum Float_num(4)[symmetric]
1735 proof (rule add_mono, fact ln_float_bounds(1)[OF lb2 ub2])
1736 have "Ifloat (?lthird * lb_ln_horner prec (get_even prec) 1 ?lthird) \<le> ln (Ifloat ?lthird + 1)"
1737 using ln_float_bounds(1)[OF lb3_lb lb3_ub] .
1738 also have "\<dots> \<le> ln (1 / 3 + 1)" unfolding ln_le_cancel_iff[OF lthird_gt0 third_gt0] using lb3 by auto
1739 finally show "Ifloat (?lthird * lb_ln_horner prec (get_even prec) 1 ?lthird) \<le> ln (1 / 3 + 1)" .
1743 subsection "Compute the logarithm in the entire domain"
1745 function ub_ln :: "nat \<Rightarrow> float \<Rightarrow> float option" and lb_ln :: "nat \<Rightarrow> float \<Rightarrow> float option" where
1746 "ub_ln prec x = (if x \<le> 0 then None
1747 else if x < 1 then Some (- the (lb_ln prec (float_divl (max prec 1) 1 x)))
1748 else let horner = \<lambda>x. (x - 1) * ub_ln_horner prec (get_odd prec) 1 (x - 1) in
1749 if x < Float 1 1 then Some (horner x)
1750 else let l = bitlen (mantissa x) - 1 in
1751 Some (ub_ln2 prec * (Float (scale x + l) 0) + horner (Float (mantissa x) (- l))))" |
1752 "lb_ln prec x = (if x \<le> 0 then None
1753 else if x < 1 then Some (- the (ub_ln prec (float_divr prec 1 x)))
1754 else let horner = \<lambda>x. (x - 1) * lb_ln_horner prec (get_even prec) 1 (x - 1) in
1755 if x < Float 1 1 then Some (horner x)
1756 else let l = bitlen (mantissa x) - 1 in
1757 Some (lb_ln2 prec * (Float (scale x + l) 0) + horner (Float (mantissa x) (- l))))"
1758 by pat_completeness auto
1760 termination proof (relation "measure (\<lambda> v. let (prec, x) = sum_case id id v in (if x < 1 then 1 else 0))", auto)
1761 fix prec x assume "\<not> x \<le> 0" and "x < 1" and "float_divl (max prec (Suc 0)) 1 x < 1"
1762 hence "0 < x" and "0 < max prec (Suc 0)" unfolding less_float_def le_float_def by auto
1763 from float_divl_pos_less1_bound[OF `0 < x` `x < 1` `0 < max prec (Suc 0)`]
1764 show False using `float_divl (max prec (Suc 0)) 1 x < 1` unfolding less_float_def le_float_def by auto
1766 fix prec x assume "\<not> x \<le> 0" and "x < 1" and "float_divr prec 1 x < 1"
1767 hence "0 < x" unfolding less_float_def le_float_def by auto
1768 from float_divr_pos_less1_lower_bound[OF `0 < x` `x < 1`, of prec]
1769 show False using `float_divr prec 1 x < 1` unfolding less_float_def le_float_def by auto
1772 lemma ln_shifted_float: assumes "0 < m" shows "ln (Ifloat (Float m e)) = ln 2 * real (e + (bitlen m - 1)) + ln (Ifloat (Float m (- (bitlen m - 1))))"
1774 let ?B = "2^nat (bitlen m - 1)"
1775 have "0 < real m" and "\<And>X. (0 :: real) < 2^X" and "0 < (2 :: real)" and "m \<noteq> 0" using assms by auto
1776 hence "0 \<le> bitlen m - 1" using bitlen_ge1[OF `m \<noteq> 0`] by auto
1778 proof (cases "0 \<le> e")
1780 show ?thesis unfolding normalized_float[OF `m \<noteq> 0`]
1781 unfolding ln_div[OF `0 < real m` `0 < ?B`] real_of_int_add ln_realpow[OF `0 < 2`]
1782 unfolding Ifloat_ge0_exp[OF True] ln_mult[OF `0 < real m` `0 < 2^nat e`]
1783 ln_realpow[OF `0 < 2`] algebra_simps using `0 \<le> bitlen m - 1` True by auto
1785 case False hence "0 < -e" by auto
1786 hence pow_gt0: "(0::real) < 2^nat (-e)" by auto
1787 hence inv_gt0: "(0::real) < inverse (2^nat (-e))" by auto
1788 show ?thesis unfolding normalized_float[OF `m \<noteq> 0`]
1789 unfolding ln_div[OF `0 < real m` `0 < ?B`] real_of_int_add ln_realpow[OF `0 < 2`]
1790 unfolding Ifloat_nge0_exp[OF False] ln_mult[OF `0 < real m` inv_gt0] ln_inverse[OF pow_gt0]
1791 ln_realpow[OF `0 < 2`] algebra_simps using `0 \<le> bitlen m - 1` False by auto
1795 lemma ub_ln_lb_ln_bounds': assumes "1 \<le> x"
1796 shows "Ifloat (the (lb_ln prec x)) \<le> ln (Ifloat x) \<and> ln (Ifloat x) \<le> Ifloat (the (ub_ln prec x))"
1797 (is "?lb \<le> ?ln \<and> ?ln \<le> ?ub")
1798 proof (cases "x < Float 1 1")
1799 case True hence "Ifloat (x - 1) < 1" unfolding less_float_def Float_num by auto
1800 have "\<not> x \<le> 0" and "\<not> x < 1" using `1 \<le> x` unfolding less_float_def le_float_def by auto
1801 hence "0 \<le> Ifloat (x - 1)" using `1 \<le> x` unfolding less_float_def Float_num by auto
1802 show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps Let_def
1803 using ln_float_bounds[OF `0 \<le> Ifloat (x - 1)` `Ifloat (x - 1) < 1`] `\<not> x \<le> 0` `\<not> x < 1` True by auto
1806 have "\<not> x \<le> 0" and "\<not> x < 1" "0 < x" using `1 \<le> x` unfolding less_float_def le_float_def by auto
1810 let ?s = "Float (e + (bitlen m - 1)) 0"
1811 let ?x = "Float m (- (bitlen m - 1))"
1813 have "0 < m" and "m \<noteq> 0" using float_pos_m_pos `0 < x` Float by auto
1816 have "Ifloat (lb_ln2 prec * ?s) \<le> ln 2 * real (e + (bitlen m - 1))" (is "?lb2 \<le> _")
1817 unfolding Ifloat_mult Ifloat_ge0_exp[OF order_refl] nat_0 power_0 mult_1_right
1818 using lb_ln2[of prec]
1819 proof (rule mult_right_mono)
1820 have "1 \<le> Float m e" using `1 \<le> x` Float unfolding le_float_def by auto
1821 from float_gt1_scale[OF this]
1822 show "0 \<le> real (e + (bitlen m - 1))" by auto
1825 from bitlen_div[OF `0 < m`, unfolded normalized_float[OF `m \<noteq> 0`, symmetric]]
1826 have "0 \<le> Ifloat (?x - 1)" and "Ifloat (?x - 1) < 1" by auto
1827 from ln_float_bounds(1)[OF this]
1828 have "Ifloat ((?x - 1) * lb_ln_horner prec (get_even prec) 1 (?x - 1)) \<le> ln (Ifloat ?x)" (is "?lb_horner \<le> _") by auto
1829 ultimately have "?lb2 + ?lb_horner \<le> ln (Ifloat x)"
1830 unfolding Float ln_shifted_float[OF `0 < m`, of e] by auto
1834 from bitlen_div[OF `0 < m`, unfolded normalized_float[OF `m \<noteq> 0`, symmetric]]
1835 have "0 \<le> Ifloat (?x - 1)" and "Ifloat (?x - 1) < 1" by auto
1836 from ln_float_bounds(2)[OF this]
1837 have "ln (Ifloat ?x) \<le> Ifloat ((?x - 1) * ub_ln_horner prec (get_odd prec) 1 (?x - 1))" (is "_ \<le> ?ub_horner") by auto
1839 have "ln 2 * real (e + (bitlen m - 1)) \<le> Ifloat (ub_ln2 prec * ?s)" (is "_ \<le> ?ub2")
1840 unfolding Ifloat_mult Ifloat_ge0_exp[OF order_refl] nat_0 power_0 mult_1_right
1841 using ub_ln2[of prec]
1842 proof (rule mult_right_mono)
1843 have "1 \<le> Float m e" using `1 \<le> x` Float unfolding le_float_def by auto
1844 from float_gt1_scale[OF this]
1845 show "0 \<le> real (e + (bitlen m - 1))" by auto
1847 ultimately have "ln (Ifloat x) \<le> ?ub2 + ?ub_horner"
1848 unfolding Float ln_shifted_float[OF `0 < m`, of e] by auto
1850 ultimately show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps
1851 unfolding if_not_P[OF `\<not> x \<le> 0`] if_not_P[OF `\<not> x < 1`] if_not_P[OF False] Let_def
1852 unfolding scale.simps[of m e, unfolded Float[symmetric]] mantissa.simps[of m e, unfolded Float[symmetric]] Ifloat_add by auto
1856 lemma ub_ln_lb_ln_bounds: assumes "0 < x"
1857 shows "Ifloat (the (lb_ln prec x)) \<le> ln (Ifloat x) \<and> ln (Ifloat x) \<le> Ifloat (the (ub_ln prec x))"
1858 (is "?lb \<le> ?ln \<and> ?ln \<le> ?ub")
1859 proof (cases "x < 1")
1860 case False hence "1 \<le> x" unfolding less_float_def le_float_def by auto
1861 show ?thesis using ub_ln_lb_ln_bounds'[OF `1 \<le> x`] .
1863 case True have "\<not> x \<le> 0" using `0 < x` unfolding less_float_def le_float_def by auto
1865 have "0 < Ifloat x" and "Ifloat x \<noteq> 0" using `0 < x` unfolding less_float_def by auto
1866 hence A: "0 < 1 / Ifloat x" by auto
1869 let ?divl = "float_divl (max prec 1) 1 x"
1870 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
1871 hence B: "0 < Ifloat ?divl" unfolding le_float_def by auto
1873 have "ln (Ifloat ?divl) \<le> ln (1 / Ifloat x)" unfolding ln_le_cancel_iff[OF B A] using float_divl[of _ 1 x] by auto
1874 hence "ln (Ifloat x) \<le> - ln (Ifloat ?divl)" unfolding nonzero_inverse_eq_divide[OF `Ifloat x \<noteq> 0`, symmetric] ln_inverse[OF `0 < Ifloat x`] by auto
1875 from this ub_ln_lb_ln_bounds'[OF A', THEN conjunct1, THEN le_imp_neg_le]
1876 have "?ln \<le> Ifloat (- the (lb_ln prec ?divl))" unfolding Ifloat_minus by (rule order_trans)
1879 let ?divr = "float_divr prec 1 x"
1880 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
1881 hence B: "0 < Ifloat ?divr" unfolding le_float_def by auto
1883 have "ln (1 / Ifloat x) \<le> ln (Ifloat ?divr)" unfolding ln_le_cancel_iff[OF A B] using float_divr[of 1 x] by auto
1884 hence "- ln (Ifloat ?divr) \<le> ln (Ifloat x)" unfolding nonzero_inverse_eq_divide[OF `Ifloat x \<noteq> 0`, symmetric] ln_inverse[OF `0 < Ifloat x`] by auto
1885 from ub_ln_lb_ln_bounds'[OF A', THEN conjunct2, THEN le_imp_neg_le] this
1886 have "Ifloat (- the (ub_ln prec ?divr)) \<le> ?ln" unfolding Ifloat_minus by (rule order_trans)
1888 ultimately show ?thesis unfolding lb_ln.simps[where x=x] ub_ln.simps[where x=x]
1889 unfolding if_not_P[OF `\<not> x \<le> 0`] if_P[OF True] by auto
1892 lemma lb_ln: assumes "Some y = lb_ln prec x"
1893 shows "Ifloat y \<le> ln (Ifloat x)" and "0 < Ifloat x"
1897 assume "\<not> 0 < x" hence "x \<le> 0" unfolding le_float_def less_float_def by auto
1898 thus False using assms by auto
1900 thus "0 < Ifloat x" unfolding less_float_def by auto
1901 have "Ifloat (the (lb_ln prec x)) \<le> ln (Ifloat x)" using ub_ln_lb_ln_bounds[OF `0 < x`] ..
1902 thus "Ifloat y \<le> ln (Ifloat x)" unfolding assms[symmetric] by auto
1905 lemma ub_ln: assumes "Some y = ub_ln prec x"
1906 shows "ln (Ifloat x) \<le> Ifloat y" and "0 < Ifloat x"
1910 assume "\<not> 0 < x" hence "x \<le> 0" unfolding le_float_def less_float_def by auto
1911 thus False using assms by auto
1913 thus "0 < Ifloat x" unfolding less_float_def by auto
1914 have "ln (Ifloat x) \<le> Ifloat (the (ub_ln prec x))" using ub_ln_lb_ln_bounds[OF `0 < x`] ..
1915 thus "ln (Ifloat x) \<le> Ifloat y" unfolding assms[symmetric] by auto
1918 lemma bnds_ln: "\<forall> x lx ux. (Some l, Some u) = (lb_ln prec lx, ub_ln prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux} \<longrightarrow> Ifloat l \<le> ln x \<and> ln x \<le> Ifloat u"
1919 proof (rule allI, rule allI, rule allI, rule impI)
1921 assume "(Some l, Some u) = (lb_ln prec lx, ub_ln prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux}"
1922 hence l: "Some l = lb_ln prec lx " and u: "Some u = ub_ln prec ux" and x: "x \<in> {Ifloat lx .. Ifloat ux}" by auto
1924 have "ln (Ifloat ux) \<le> Ifloat u" and "0 < Ifloat ux" using ub_ln u by auto
1925 have "Ifloat l \<le> ln (Ifloat lx)" and "0 < Ifloat lx" and "0 < x" using lb_ln[OF l] x by auto
1927 from ln_le_cancel_iff[OF `0 < Ifloat lx` `0 < x`] `Ifloat l \<le> ln (Ifloat lx)`
1928 have "Ifloat l \<le> ln x" using x unfolding atLeastAtMost_iff by auto
1930 from ln_le_cancel_iff[OF `0 < x` `0 < Ifloat ux`] `ln (Ifloat ux) \<le> Ifloat u`
1931 have "ln x \<le> Ifloat u" using x unfolding atLeastAtMost_iff by auto
1932 ultimately show "Ifloat l \<le> ln x \<and> ln x \<le> Ifloat u" ..
1936 section "Implement floatarith"
1938 subsection "Define syntax and semantics"
1941 = Add floatarith floatarith
1943 | Mult floatarith floatarith
1944 | Inverse floatarith
1949 | Max floatarith floatarith
1950 | Min floatarith floatarith
1955 | Power floatarith nat
1959 fun Ifloatarith :: "floatarith \<Rightarrow> real list \<Rightarrow> real"
1961 "Ifloatarith (Add a b) vs = (Ifloatarith a vs) + (Ifloatarith b vs)" |
1962 "Ifloatarith (Minus a) vs = - (Ifloatarith a vs)" |
1963 "Ifloatarith (Mult a b) vs = (Ifloatarith a vs) * (Ifloatarith b vs)" |
1964 "Ifloatarith (Inverse a) vs = inverse (Ifloatarith a vs)" |
1965 "Ifloatarith (Sin a) vs = sin (Ifloatarith a vs)" |
1966 "Ifloatarith (Cos a) vs = cos (Ifloatarith a vs)" |
1967 "Ifloatarith (Arctan a) vs = arctan (Ifloatarith a vs)" |
1968 "Ifloatarith (Min a b) vs = min (Ifloatarith a vs) (Ifloatarith b vs)" |
1969 "Ifloatarith (Max a b) vs = max (Ifloatarith a vs) (Ifloatarith b vs)" |
1970 "Ifloatarith (Abs a) vs = abs (Ifloatarith a vs)" |
1971 "Ifloatarith Pi vs = pi" |
1972 "Ifloatarith (Sqrt a) vs = sqrt (Ifloatarith a vs)" |
1973 "Ifloatarith (Exp a) vs = exp (Ifloatarith a vs)" |
1974 "Ifloatarith (Ln a) vs = ln (Ifloatarith a vs)" |
1975 "Ifloatarith (Power a n) vs = (Ifloatarith a vs)^n" |
1976 "Ifloatarith (Num f) vs = Ifloat f" |
1977 "Ifloatarith (Atom n) vs = vs ! n"
1979 subsection "Implement approximation function"
1981 fun lift_bin :: "(float * float) option \<Rightarrow> (float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> (float option * float option)) \<Rightarrow> (float * float) option" where
1982 "lift_bin (Some (l1, u1)) (Some (l2, u2)) f = (case (f l1 u1 l2 u2) of (Some l, Some u) \<Rightarrow> Some (l, u)
1983 | t \<Rightarrow> None)" |
1984 "lift_bin a b f = None"
1986 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
1987 "lift_bin' (Some (l1, u1)) (Some (l2, u2)) f = Some (f l1 u1 l2 u2)" |
1988 "lift_bin' a b f = None"
1990 fun lift_un :: "(float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> ((float option) * (float option))) \<Rightarrow> (float * float) option" where
1991 "lift_un (Some (l1, u1)) f = (case (f l1 u1) of (Some l, Some u) \<Rightarrow> Some (l, u)
1992 | t \<Rightarrow> None)" |
1993 "lift_un b f = None"
1995 fun lift_un' :: "(float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> (float * float)) \<Rightarrow> (float * float) option" where
1996 "lift_un' (Some (l1, u1)) f = Some (f l1 u1)" |
1997 "lift_un' b f = None"
1999 fun bounded_by :: "real list \<Rightarrow> (float * float) list \<Rightarrow> bool " where
2000 bounded_by_Cons: "bounded_by (v#vs) ((l, u)#bs) = ((Ifloat l \<le> v \<and> v \<le> Ifloat u) \<and> bounded_by vs bs)" |
2001 bounded_by_Nil: "bounded_by [] [] = True" |
2002 "bounded_by _ _ = False"
2004 lemma bounded_by: assumes "bounded_by vs bs" and "i < length bs"
2005 shows "Ifloat (fst (bs ! i)) \<le> vs ! i \<and> vs ! i \<le> Ifloat (snd (bs ! i))"
2006 using `bounded_by vs bs` and `i < length bs`
2007 proof (induct arbitrary: i rule: bounded_by.induct)
2008 fix v :: real and vs :: "real list" and l u :: float and bs :: "(float * float) list" and i :: nat
2009 assume hyp: "\<And>i. \<lbrakk>bounded_by vs bs; i < length bs\<rbrakk> \<Longrightarrow> Ifloat (fst (bs ! i)) \<le> vs ! i \<and> vs ! i \<le> Ifloat (snd (bs ! i))"
2010 assume bounded: "bounded_by (v # vs) ((l, u) # bs)" and length: "i < length ((l, u) # bs)"
2011 show "Ifloat (fst (((l, u) # bs) ! i)) \<le> (v # vs) ! i \<and> (v # vs) ! i \<le> Ifloat (snd (((l, u) # bs) ! i))"
2014 show ?thesis using bounded unfolding 0 nth_Cons_0 fst_conv snd_conv bounded_by.simps ..
2016 case (Suc i) with length have "i < length bs" by auto
2017 show ?thesis unfolding Suc nth_Cons_Suc bounded_by.simps
2018 using hyp[OF bounded[unfolded bounded_by.simps, THEN conjunct2] `i < length bs`] .
2022 fun approx approx' :: "nat \<Rightarrow> floatarith \<Rightarrow> (float * float) list \<Rightarrow> (float * float) option" where
2023 "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)" |
2024 "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))" |
2025 "approx prec (Minus a) bs = lift_un' (approx' prec a bs) (\<lambda> l u. (-u, -l))" |
2026 "approx prec (Mult a b) bs = lift_bin' (approx' prec a bs) (approx' prec b bs)
2027 (\<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,
2028 float_pprt a2 * float_pprt b2 + float_pprt a1 * float_nprt b2 + float_nprt a2 * float_pprt b1 + float_nprt a1 * float_nprt b1))" |
2029 "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))" |
2030 "approx prec (Sin a) bs = lift_un' (approx' prec a bs) (bnds_sin prec)" |
2031 "approx prec (Cos a) bs = lift_un' (approx' prec a bs) (bnds_cos prec)" |
2032 "approx prec Pi bs = Some (lb_pi prec, ub_pi prec)" |
2033 "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))" |
2034 "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))" |
2035 "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>))" |
2036 "approx prec (Arctan a) bs = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_arctan prec l, ub_arctan prec u))" |
2037 "approx prec (Sqrt a) bs = lift_un (approx' prec a bs) (\<lambda> l u. (lb_sqrt prec l, ub_sqrt prec u))" |
2038 "approx prec (Exp a) bs = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_exp prec l, ub_exp prec u))" |
2039 "approx prec (Ln a) bs = lift_un (approx' prec a bs) (\<lambda> l u. (lb_ln prec l, ub_ln prec u))" |
2040 "approx prec (Power a n) bs = lift_un' (approx' prec a bs) (float_power_bnds n)" |
2041 "approx prec (Num f) bs = Some (f, f)" |
2042 "approx prec (Atom i) bs = (if i < length bs then Some (bs ! i) else None)"
2045 assumes lift_bin'_Some: "Some (l, u) = lift_bin' a b f"
2046 shows "\<exists> l1 u1 l2 u2. Some (l1, u1) = a \<and> Some (l2, u2) = b"
2048 case None hence "None = lift_bin' a b f" unfolding None lift_bin'.simps ..
2049 thus ?thesis using lift_bin'_Some by auto
2054 case None hence "None = lift_bin' a b f" unfolding None lift_bin'.simps ..
2055 thus ?thesis using lift_bin'_Some by auto
2058 obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
2059 obtain lb ub where b': "b' = (lb, ub)" by (cases b', auto)
2060 thus ?thesis unfolding `a = Some a'` `b = Some b'` a' b' by auto
2065 assumes lift_bin'_Some: "Some (l, u) = lift_bin' (g a) (g b) f"
2066 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"
2067 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)"
2070 where Sa: "Some (l1, u1) = g a" and Sb: "Some (l2, u2) = g b" using lift_bin'_ex[OF assms(1)] by auto
2071 have lu: "(l, u) = f l1 u1 l2 u2" using lift_bin'_Some[unfolded Sa[symmetric] Sb[symmetric] lift_bin'.simps] by auto
2072 have "l = fst (f l1 u1 l2 u2)" and "u = snd (f l1 u1 l2 u2)" unfolding lu[symmetric] by auto
2073 thus ?thesis using Pa[OF Sa] Pb[OF Sb] by auto
2076 lemma approx_approx':
2077 assumes Pa: "\<And>l u. Some (l, u) = approx prec a vs \<Longrightarrow> Ifloat l \<le> Ifloatarith a xs \<and> Ifloatarith a xs \<le> Ifloat u"
2078 and approx': "Some (l, u) = approx' prec a vs"
2079 shows "Ifloat l \<le> Ifloatarith a xs \<and> Ifloatarith a xs \<le> Ifloat u"
2081 obtain l' u' where S: "Some (l', u') = approx prec a vs"
2082 using approx' unfolding approx'.simps by (cases "approx prec a vs", auto)
2083 have l': "l = round_down prec l'" and u': "u = round_up prec u'"
2084 using approx' unfolding approx'.simps S[symmetric] by auto
2085 show ?thesis unfolding l' u'
2086 using order_trans[OF Pa[OF S, THEN conjunct2] round_up[of u']]
2087 using order_trans[OF round_down[of _ l'] Pa[OF S, THEN conjunct1]] by auto
2091 assumes lift_bin'_Some: "Some (l, u) = lift_bin' (approx' prec a bs) (approx' prec b bs) f"
2092 and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow> Ifloat l \<le> Ifloatarith a xs \<and> Ifloatarith a xs \<le> Ifloat u" (is "\<And>l u. _ = ?g a \<Longrightarrow> ?P l u a")
2093 and Pb: "\<And>l u. Some (l, u) = approx prec b bs \<Longrightarrow> Ifloat l \<le> Ifloatarith b xs \<and> Ifloatarith b xs \<le> Ifloat u"
2094 shows "\<exists> l1 u1 l2 u2. (Ifloat l1 \<le> Ifloatarith a xs \<and> Ifloatarith a xs \<le> Ifloat u1) \<and>
2095 (Ifloat l2 \<le> Ifloatarith b xs \<and> Ifloatarith b xs \<le> Ifloat u2) \<and>
2096 l = fst (f l1 u1 l2 u2) \<and> u = snd (f l1 u1 l2 u2)"
2098 { fix l u assume "Some (l, u) = approx' prec a bs"
2099 with approx_approx'[of prec a bs, OF _ this] Pa
2100 have "Ifloat l \<le> Ifloatarith a xs \<and> Ifloatarith a xs \<le> Ifloat u" by auto } note Pa = this
2101 { fix l u assume "Some (l, u) = approx' prec b bs"
2102 with approx_approx'[of prec b bs, OF _ this] Pb
2103 have "Ifloat l \<le> Ifloatarith b xs \<and> Ifloatarith b xs \<le> Ifloat u" by auto } note Pb = this
2105 from lift_bin'_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_bin'_Some, OF Pa Pb]
2106 show ?thesis by auto
2110 assumes lift_un'_Some: "Some (l, u) = lift_un' a f"
2111 shows "\<exists> l u. Some (l, u) = a"
2113 case None hence "None = lift_un' a f" unfolding None lift_un'.simps ..
2114 thus ?thesis using lift_un'_Some by auto
2117 obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
2118 thus ?thesis unfolding `a = Some a'` a' by auto
2122 assumes lift_un'_Some: "Some (l, u) = lift_un' (g a) f"
2123 and Pa: "\<And>l u. Some (l, u) = g a \<Longrightarrow> P l u a"
2124 shows "\<exists> l1 u1. P l1 u1 a \<and> l = fst (f l1 u1) \<and> u = snd (f l1 u1)"
2126 obtain l1 u1 where Sa: "Some (l1, u1) = g a" using lift_un'_ex[OF assms(1)] by auto
2127 have lu: "(l, u) = f l1 u1" using lift_un'_Some[unfolded Sa[symmetric] lift_un'.simps] by auto
2128 have "l = fst (f l1 u1)" and "u = snd (f l1 u1)" unfolding lu[symmetric] by auto
2129 thus ?thesis using Pa[OF Sa] by auto
2133 assumes lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f"
2134 and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow> Ifloat l \<le> Ifloatarith a xs \<and> Ifloatarith a xs \<le> Ifloat u" (is "\<And>l u. _ = ?g a \<Longrightarrow> ?P l u a")
2135 shows "\<exists> l1 u1. (Ifloat l1 \<le> Ifloatarith a xs \<and> Ifloatarith a xs \<le> Ifloat u1) \<and>
2136 l = fst (f l1 u1) \<and> u = snd (f l1 u1)"
2138 { fix l u assume "Some (l, u) = approx' prec a bs"
2139 with approx_approx'[of prec a bs, OF _ this] Pa
2140 have "Ifloat l \<le> Ifloatarith a xs \<and> Ifloatarith a xs \<le> Ifloat u" by auto } note Pa = this
2141 from lift_un'_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_un'_Some, OF Pa]
2142 show ?thesis by auto
2145 lemma lift_un'_bnds:
2146 assumes bnds: "\<forall> x lx ux. (l, u) = f lx ux \<and> x \<in> { Ifloat lx .. Ifloat ux } \<longrightarrow> Ifloat l \<le> f' x \<and> f' x \<le> Ifloat u"
2147 and lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f"
2148 and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow> Ifloat l \<le> Ifloatarith a xs \<and> Ifloatarith a xs \<le> Ifloat u"
2149 shows "Ifloat l \<le> f' (Ifloatarith a xs) \<and> f' (Ifloatarith a xs) \<le> Ifloat u"
2151 from lift_un'[OF lift_un'_Some Pa]
2152 obtain l1 u1 where "Ifloat l1 \<le> Ifloatarith a xs" and "Ifloatarith a xs \<le> Ifloat u1" and "l = fst (f l1 u1)" and "u = snd (f l1 u1)" by blast
2153 hence "(l, u) = f l1 u1" and "Ifloatarith a xs \<in> {Ifloat l1 .. Ifloat u1}" by auto
2154 thus ?thesis using bnds by auto
2158 assumes lift_un_Some: "Some (l, u) = lift_un a f"
2159 shows "\<exists> l u. Some (l, u) = a"
2161 case None hence "None = lift_un a f" unfolding None lift_un.simps ..
2162 thus ?thesis using lift_un_Some by auto
2165 obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
2166 thus ?thesis unfolding `a = Some a'` a' by auto
2170 assumes lift_un_Some: "Some (l, u) = lift_un (g a) f"
2171 and Pa: "\<And>l u. Some (l, u) = g a \<Longrightarrow> P l u a"
2172 shows "\<exists> l1 u1. P l1 u1 a \<and> Some l = fst (f l1 u1) \<and> Some u = snd (f l1 u1)"
2174 obtain l1 u1 where Sa: "Some (l1, u1) = g a" using lift_un_ex[OF assms(1)] by auto
2175 have "fst (f l1 u1) \<noteq> None \<and> snd (f l1 u1) \<noteq> None"
2177 assume "\<not> (fst (f l1 u1) \<noteq> None \<and> snd (f l1 u1) \<noteq> None)"
2178 hence or: "fst (f l1 u1) = None \<or> snd (f l1 u1) = None" by auto
2179 hence "lift_un (g a) f = None"
2180 proof (cases "fst (f l1 u1) = None")
2182 then obtain b where b: "f l1 u1 = (None, b)" by (cases "f l1 u1", auto)
2183 thus ?thesis unfolding Sa[symmetric] lift_un.simps b by auto
2185 case False hence "snd (f l1 u1) = None" using or by auto
2186 with False obtain b where b: "f l1 u1 = (Some b, None)" by (cases "f l1 u1", auto)
2187 thus ?thesis unfolding Sa[symmetric] lift_un.simps b by auto
2189 thus False using lift_un_Some by auto
2191 then obtain a' b' where f: "f l1 u1 = (Some a', Some b')" by (cases "f l1 u1", auto)
2192 from lift_un_Some[unfolded Sa[symmetric] lift_un.simps f]
2193 have "Some l = fst (f l1 u1)" and "Some u = snd (f l1 u1)" unfolding f by auto
2194 thus ?thesis unfolding Sa[symmetric] lift_un.simps using Pa[OF Sa] by auto
2198 assumes lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f"
2199 and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow> Ifloat l \<le> Ifloatarith a xs \<and> Ifloatarith a xs \<le> Ifloat u" (is "\<And>l u. _ = ?g a \<Longrightarrow> ?P l u a")
2200 shows "\<exists> l1 u1. (Ifloat l1 \<le> Ifloatarith a xs \<and> Ifloatarith a xs \<le> Ifloat u1) \<and>
2201 Some l = fst (f l1 u1) \<and> Some u = snd (f l1 u1)"
2203 { fix l u assume "Some (l, u) = approx' prec a bs"
2204 with approx_approx'[of prec a bs, OF _ this] Pa
2205 have "Ifloat l \<le> Ifloatarith a xs \<and> Ifloatarith a xs \<le> Ifloat u" by auto } note Pa = this
2206 from lift_un_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_un_Some, OF Pa]
2207 show ?thesis by auto
2211 assumes bnds: "\<forall> x lx ux. (Some l, Some u) = f lx ux \<and> x \<in> { Ifloat lx .. Ifloat ux } \<longrightarrow> Ifloat l \<le> f' x \<and> f' x \<le> Ifloat u"
2212 and lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f"
2213 and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow> Ifloat l \<le> Ifloatarith a xs \<and> Ifloatarith a xs \<le> Ifloat u"
2214 shows "Ifloat l \<le> f' (Ifloatarith a xs) \<and> f' (Ifloatarith a xs) \<le> Ifloat u"
2216 from lift_un[OF lift_un_Some Pa]
2217 obtain l1 u1 where "Ifloat l1 \<le> Ifloatarith a xs" and "Ifloatarith a xs \<le> Ifloat u1" and "Some l = fst (f l1 u1)" and "Some u = snd (f l1 u1)" by blast
2218 hence "(Some l, Some u) = f l1 u1" and "Ifloatarith a xs \<in> {Ifloat l1 .. Ifloat u1}" by auto
2219 thus ?thesis using bnds by auto
2223 assumes "bounded_by xs vs"
2224 and "Some (l, u) = approx prec arith vs" (is "_ = ?g arith")
2225 shows "Ifloat l \<le> Ifloatarith arith xs \<and> Ifloatarith arith xs \<le> Ifloat u" (is "?P l u arith")
2226 using `Some (l, u) = approx prec arith vs`
2227 proof (induct arith arbitrary: l u x)
2229 from lift_bin'[OF Add.prems[unfolded approx.simps]] Add.hyps
2230 obtain l1 u1 l2 u2 where "l = l1 + l2" and "u = u1 + u2"
2231 "Ifloat l1 \<le> Ifloatarith a xs" and "Ifloatarith a xs \<le> Ifloat u1"
2232 "Ifloat l2 \<le> Ifloatarith b xs" and "Ifloatarith b xs \<le> Ifloat u2" unfolding fst_conv snd_conv by blast
2233 thus ?case unfolding Ifloatarith.simps by auto
2236 from lift_un'[OF Minus.prems[unfolded approx.simps]] Minus.hyps
2237 obtain l1 u1 where "l = -u1" and "u = -l1"
2238 "Ifloat l1 \<le> Ifloatarith a xs" and "Ifloatarith a xs \<le> Ifloat u1" unfolding fst_conv snd_conv by blast
2239 thus ?case unfolding Ifloatarith.simps using Ifloat_minus by auto
2242 from lift_bin'[OF Mult.prems[unfolded approx.simps]] Mult.hyps
2244 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"
2245 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"
2246 and "Ifloat l1 \<le> Ifloatarith a xs" and "Ifloatarith a xs \<le> Ifloat u1"
2247 and "Ifloat l2 \<le> Ifloatarith b xs" and "Ifloatarith b xs \<le> Ifloat u2" unfolding fst_conv snd_conv by blast
2248 thus ?case unfolding Ifloatarith.simps l u Ifloat_add Ifloat_mult Ifloat_nprt Ifloat_pprt
2249 using mult_le_prts mult_ge_prts by auto
2252 from lift_un[OF Inverse.prems[unfolded approx.simps], unfolded if_distrib[of fst] if_distrib[of snd] fst_conv snd_conv] Inverse.hyps
2253 obtain l1 u1 where l': "Some l = (if 0 < l1 \<or> u1 < 0 then Some (float_divl prec 1 u1) else None)"
2254 and u': "Some u = (if 0 < l1 \<or> u1 < 0 then Some (float_divr prec 1 l1) else None)"
2255 and l1: "Ifloat l1 \<le> Ifloatarith a xs" and u1: "Ifloatarith a xs \<le> Ifloat u1" by blast
2256 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
2257 moreover have l1_le_u1: "Ifloat l1 \<le> Ifloat u1" using l1 u1 by auto
2258 ultimately have "Ifloat l1 \<noteq> 0" and "Ifloat u1 \<noteq> 0" unfolding less_float_def by auto
2260 have inv: "inverse (Ifloat u1) \<le> inverse (Ifloatarith a xs)
2261 \<and> inverse (Ifloatarith a xs) \<le> inverse (Ifloat l1)"
2262 proof (cases "0 < l1")
2263 case True hence "0 < Ifloat u1" and "0 < Ifloat l1" "0 < Ifloatarith a xs"
2264 unfolding less_float_def using l1_le_u1 l1 by auto
2266 unfolding inverse_le_iff_le[OF `0 < Ifloat u1` `0 < Ifloatarith a xs`]
2267 inverse_le_iff_le[OF `0 < Ifloatarith a xs` `0 < Ifloat l1`]
2270 case False hence "u1 < 0" using either by blast
2271 hence "Ifloat u1 < 0" and "Ifloat l1 < 0" "Ifloatarith a xs < 0"
2272 unfolding less_float_def using l1_le_u1 u1 by auto
2274 unfolding inverse_le_iff_le_neg[OF `Ifloat u1 < 0` `Ifloatarith a xs < 0`]
2275 inverse_le_iff_le_neg[OF `Ifloatarith a xs < 0` `Ifloat l1 < 0`]
2279 from l' have "l = float_divl prec 1 u1" by (cases "0 < l1 \<or> u1 < 0", auto)
2280 hence "Ifloat l \<le> inverse (Ifloat u1)" unfolding nonzero_inverse_eq_divide[OF `Ifloat u1 \<noteq> 0`] using float_divl[of prec 1 u1] by auto
2281 also have "\<dots> \<le> inverse (Ifloatarith a xs)" using inv by auto
2282 finally have "Ifloat l \<le> inverse (Ifloatarith a xs)" .
2284 from u' have "u = float_divr prec 1 l1" by (cases "0 < l1 \<or> u1 < 0", auto)
2285 hence "inverse (Ifloat l1) \<le> Ifloat u" unfolding nonzero_inverse_eq_divide[OF `Ifloat l1 \<noteq> 0`] using float_divr[of 1 l1 prec] by auto
2286 hence "inverse (Ifloatarith a xs) \<le> Ifloat u" by (rule order_trans[OF inv[THEN conjunct2]])
2287 ultimately show ?case unfolding Ifloatarith.simps using l1 u1 by auto
2290 from lift_un'[OF Abs.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Abs.hyps
2291 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>"
2292 and l1: "Ifloat l1 \<le> Ifloatarith x xs" and u1: "Ifloatarith x xs \<le> Ifloat u1" by blast
2293 thus ?case unfolding l' u' by (cases "l1 < 0 \<and> 0 < u1", auto simp add: Ifloat_min Ifloat_max Ifloat_abs less_float_def)
2296 from lift_bin'[OF Min.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Min.hyps
2297 obtain l1 u1 l2 u2 where l': "l = min l1 l2" and u': "u = min u1 u2"
2298 and l1: "Ifloat l1 \<le> Ifloatarith a xs" and u1: "Ifloatarith a xs \<le> Ifloat u1"
2299 and l1: "Ifloat l2 \<le> Ifloatarith b xs" and u1: "Ifloatarith b xs \<le> Ifloat u2" by blast
2300 thus ?case unfolding l' u' by (auto simp add: Ifloat_min)
2303 from lift_bin'[OF Max.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Max.hyps
2304 obtain l1 u1 l2 u2 where l': "l = max l1 l2" and u': "u = max u1 u2"
2305 and l1: "Ifloat l1 \<le> Ifloatarith a xs" and u1: "Ifloatarith a xs \<le> Ifloat u1"
2306 and l1: "Ifloat l2 \<le> Ifloatarith b xs" and u1: "Ifloatarith b xs \<le> Ifloat u2" by blast
2307 thus ?case unfolding l' u' by (auto simp add: Ifloat_max)
2308 next case (Sin a) with lift_un'_bnds[OF bnds_sin] show ?case by auto
2309 next case (Cos a) with lift_un'_bnds[OF bnds_cos] show ?case by auto
2310 next case (Arctan a) with lift_un'_bnds[OF bnds_arctan] show ?case by auto
2311 next case Pi with pi_boundaries show ?case by auto
2312 next case (Sqrt a) with lift_un_bnds[OF bnds_sqrt] show ?case by auto
2313 next case (Exp a) with lift_un'_bnds[OF bnds_exp] show ?case by auto
2314 next case (Ln a) with lift_un_bnds[OF bnds_ln] show ?case by auto
2315 next case (Power a n) with lift_un'_bnds[OF bnds_power] show ?case by auto
2316 next case (Num f) thus ?case by auto
2320 proof (cases "n < length vs")
2322 with Atom have "vs ! n = (l, u)" by auto
2323 thus ?thesis using bounded_by[OF assms(1) True] by auto
2325 case False thus ?thesis using Atom by auto
2329 datatype ApproxEq = Less floatarith floatarith
2330 | LessEqual floatarith floatarith
2332 fun uneq :: "ApproxEq \<Rightarrow> real list \<Rightarrow> bool" where
2333 "uneq (Less a b) vs = (Ifloatarith a vs < Ifloatarith b vs)" |
2334 "uneq (LessEqual a b) vs = (Ifloatarith a vs \<le> Ifloatarith b vs)"
2336 fun uneq' :: "nat \<Rightarrow> ApproxEq \<Rightarrow> (float * float) list \<Rightarrow> bool" where
2337 "uneq' prec (Less a b) bs = (case (approx prec a bs, approx prec b bs) of (Some (l, u), Some (l', u')) \<Rightarrow> u < l' | _ \<Rightarrow> False)" |
2338 "uneq' prec (LessEqual a b) bs = (case (approx prec a bs, approx prec b bs) of (Some (l, u), Some (l', u')) \<Rightarrow> u \<le> l' | _ \<Rightarrow> False)"
2340 lemma uneq_approx: fixes m :: nat assumes "bounded_by vs bs" and "uneq' prec eq bs"
2345 proof (cases "\<exists> u l u' l'. approx prec a bs = Some (l, u) \<and>
2346 approx prec b bs = Some (l', u')")
2348 then obtain l u l' u' where a_approx: "approx prec a bs = Some (l, u)"
2349 and b_approx: "approx prec b bs = Some (l', u') " by auto
2350 with `uneq' prec eq bs` have "Ifloat u < Ifloat l'"
2351 unfolding Less uneq'.simps less_float_def by auto
2352 moreover from a_approx[symmetric] and b_approx[symmetric] and `bounded_by vs bs`
2353 have "Ifloatarith a vs \<le> Ifloat u" and "Ifloat l' \<le> Ifloatarith b vs"
2354 using approx by auto
2355 ultimately show ?thesis unfolding uneq.simps Less by auto
2358 hence "approx prec a bs = None \<or> approx prec b bs = None"
2359 unfolding not_Some_eq[symmetric] by auto
2360 hence "\<not> uneq' prec eq bs" unfolding Less uneq'.simps
2361 by (cases "approx prec a bs = None", auto)
2362 thus ?thesis using assms by auto
2365 case (LessEqual a b)
2367 proof (cases "\<exists> u l u' l'. approx prec a bs = Some (l, u) \<and>
2368 approx prec b bs = Some (l', u')")
2370 then obtain l u l' u' where a_approx: "approx prec a bs = Some (l, u)"
2371 and b_approx: "approx prec b bs = Some (l', u') " by auto
2372 with `uneq' prec eq bs` have "Ifloat u \<le> Ifloat l'"
2373 unfolding LessEqual uneq'.simps le_float_def by auto
2374 moreover from a_approx[symmetric] and b_approx[symmetric] and `bounded_by vs bs`
2375 have "Ifloatarith a vs \<le> Ifloat u" and "Ifloat l' \<le> Ifloatarith b vs"
2376 using approx by auto
2377 ultimately show ?thesis unfolding uneq.simps LessEqual by auto
2380 hence "approx prec a bs = None \<or> approx prec b bs = None"
2381 unfolding not_Some_eq[symmetric] by auto
2382 hence "\<not> uneq' prec eq bs" unfolding LessEqual uneq'.simps
2383 by (cases "approx prec a bs = None", auto)
2384 thus ?thesis using assms by auto
2388 lemma Ifloatarith_divide: "Ifloatarith (Mult a (Inverse b)) vs = (Ifloatarith a vs) / (Ifloatarith b vs)"
2389 unfolding real_divide_def Ifloatarith.simps ..
2391 lemma Ifloatarith_diff: "Ifloatarith (Add a (Minus b)) vs = (Ifloatarith a vs) - (Ifloatarith b vs)"
2392 unfolding real_diff_def Ifloatarith.simps ..
2394 lemma Ifloatarith_tan: "Ifloatarith (Mult (Sin a) (Inverse (Cos a))) vs = tan (Ifloatarith a vs)"
2395 unfolding tan_def Ifloatarith.simps real_divide_def ..
2397 lemma Ifloatarith_powr: "Ifloatarith (Exp (Mult b (Ln a))) vs = (Ifloatarith a vs) powr (Ifloatarith b vs)"
2398 unfolding powr_def Ifloatarith.simps ..
2400 lemma Ifloatarith_log: "Ifloatarith ((Mult (Ln x) (Inverse (Ln b)))) vs = log (Ifloatarith b vs) (Ifloatarith x vs)"
2401 unfolding log_def Ifloatarith.simps real_divide_def ..
2403 lemma Ifloatarith_num: shows "Ifloatarith (Num (Float 0 0)) vs = 0" and "Ifloatarith (Num (Float 1 0)) vs = 1" and "Ifloatarith (Num (Float (number_of a) 0)) vs = number_of a" by auto
2405 subsection {* Implement proof method \texttt{approximation} *}
2407 lemma bounded_divl: assumes "Ifloat a / Ifloat b \<le> x" shows "Ifloat (float_divl p a b) \<le> x" by (rule order_trans[OF _ assms], rule float_divl)
2408 lemma bounded_divr: assumes "x \<le> Ifloat a / Ifloat b" shows "x \<le> Ifloat (float_divr p a b)" by (rule order_trans[OF assms _], rule float_divr)
2409 lemma bounded_num: shows "Ifloat (Float 5 1) = 10" and "Ifloat (Float 0 0) = 0" and "Ifloat (Float 1 0) = 1" and "Ifloat (Float (number_of n) 0) = (number_of n)"
2410 and "0 * pow2 e = Ifloat (Float 0 e)" and "1 * pow2 e = Ifloat (Float 1 e)" and "number_of m * pow2 e = Ifloat (Float (number_of m) e)"
2411 by (auto simp add: Ifloat.simps pow2_def)
2413 lemmas bounded_by_equations = bounded_by_Cons bounded_by_Nil float_power bounded_divl bounded_divr bounded_num HOL.simp_thms
2414 lemmas uneq_equations = uneq.simps Ifloatarith.simps Ifloatarith_num Ifloatarith_divide Ifloatarith_diff Ifloatarith_tan Ifloatarith_powr Ifloatarith_log
2416 lemma "x div (0::int) = 0" by auto -- "What happens in the zero case for div"
2417 lemma "x mod (0::int) = x" by auto -- "What happens in the zero case for mod"
2419 text {* The following equations must hold for div & mod
2420 -- see "The Definition of Standard ML" by R. Milner, M. Tofte and R. Harper (pg. 79) *}
2421 lemma "d * (i div d) + i mod d = (i::int)" by auto
2422 lemma "0 < (d :: int) \<Longrightarrow> 0 \<le> i mod d \<and> i mod d < d" by auto
2423 lemma "(d :: int) < 0 \<Longrightarrow> d < i mod d \<and> i mod d \<le> 0" by auto
2426 val uneq_equations = PureThy.get_thms @{theory} "uneq_equations";
2427 val bounded_by_equations = PureThy.get_thms @{theory} "bounded_by_equations";
2428 val bounded_by_simpset = (HOL_basic_ss addsimps bounded_by_equations)
2430 fun reify_uneq ctxt i = (fn st =>
2432 val to = HOLogic.dest_Trueprop (Logic.strip_imp_concl (List.nth (prems_of st, i - 1)))
2433 in (Reflection.genreify_tac ctxt uneq_equations (SOME to) i) st
2436 fun rule_uneq ctxt prec i thm = let
2437 fun conv_num typ = HOLogic.dest_number #> snd #> HOLogic.mk_number typ
2438 val to_natc = conv_num @{typ "nat"} #> Thm.cterm_of (ProofContext.theory_of ctxt)
2439 val to_nat = conv_num @{typ "nat"}
2440 val to_int = conv_num @{typ "int"}
2442 val prec' = to_nat prec
2444 fun bot_float (Const (@{const_name "times"}, _) $ mantisse $ (Const (@{const_name "pow2"}, _) $ exp))
2445 = @{term "Float"} $ to_int mantisse $ to_int exp
2446 | bot_float (Const (@{const_name "divide"}, _) $ mantisse $ (Const (@{const_name "power"}, _) $ ten $ exp))
2447 = @{term "float_divl"} $ prec' $ (@{term "Float"} $ to_int mantisse $ @{term "0 :: int"}) $ (@{term "power (Float 5 1)"} $ to_nat exp)
2448 | bot_float (Const (@{const_name "divide"}, _) $ mantisse $ ten)
2449 = @{term "float_divl"} $ prec' $ (@{term "Float"} $ to_int mantisse $ @{term "0 :: int"}) $ @{term "Float 5 1"}
2450 | bot_float mantisse = @{term "Float"} $ to_int mantisse $ @{term "0 :: int"}
2452 fun top_float (Const (@{const_name "times"}, _) $ mantisse $ (Const (@{const_name "pow2"}, _) $ exp))
2453 = @{term "Float"} $ to_int mantisse $ to_int exp
2454 | top_float (Const (@{const_name "divide"}, _) $ mantisse $ (Const (@{const_name "power"}, _) $ ten $ exp))
2455 = @{term "float_divr"} $ prec' $ (@{term "Float"} $ to_int mantisse $ @{term "0 :: int"}) $ (@{term "power (Float 5 1)"} $ to_nat exp)
2456 | top_float (Const (@{const_name "divide"}, _) $ mantisse $ ten)
2457 = @{term "float_divr"} $ prec' $ (@{term "Float"} $ to_int mantisse $ @{term "0 :: int"}) $ @{term "Float 5 1"}
2458 | top_float mantisse = @{term "Float"} $ to_int mantisse $ @{term "0 :: int"}
2460 val goal' : term = List.nth (prems_of thm, i - 1)
2462 fun lift_bnd (t as (Const (@{const_name "op &"}, _) $
2463 (Const (@{const_name "less_eq"}, _) $
2464 bottom $ (Free (name, _))) $
2465 (Const (@{const_name "less_eq"}, _) $ _ $ top)))
2466 = ((name, HOLogic.mk_prod (bot_float bottom, top_float top))
2467 handle TERM (txt, ts) => raise TERM ("Premisse needs format '<num> <= <var> & <var> <= <num>', but found " ^
2468 (Syntax.string_of_term ctxt t), [t]))
2469 | lift_bnd t = raise TERM ("Premisse needs format '<num> <= <var> & <var> <= <num>', but found " ^
2470 (Syntax.string_of_term ctxt t), [t])
2471 val bound_eqs = map (HOLogic.dest_Trueprop #> lift_bnd) (Logic.strip_imp_prems goal')
2473 fun lift_var (Free (varname, _)) = (case AList.lookup (op =) bound_eqs varname of
2475 | NONE => raise TERM ("No bound equations found for " ^ varname, []))
2476 | lift_var t = raise TERM ("Can not convert expression " ^
2477 (Syntax.string_of_term ctxt t), [t])
2479 val _ $ vs = HOLogic.dest_Trueprop (Logic.strip_imp_concl goal')
2481 val bs = (HOLogic.dest_list #> map lift_var #> HOLogic.mk_list @{typ "float * float"}) vs
2482 val map = [(@{cpat "?prec::nat"}, to_natc prec),
2483 (@{cpat "?bs::(float * float) list"}, Thm.cterm_of (ProofContext.theory_of ctxt) bs)]
2484 in rtac (Thm.instantiate ([], map) @{thm "uneq_approx"}) i thm end
2486 val eval_tac = CSUBGOAL (fn (ct, i) => rtac (eval_oracle ct) i)
2488 fun gen_eval_tac conv ctxt = CONVERSION (Conv.params_conv (~1) (K (Conv.concl_conv (~1) conv)) ctxt)
2493 method_setup approximation = {* fn src =>
2494 Method.syntax Args.term src #>
2495 (fn (prec, ctxt) => let
2496 in Method.SIMPLE_METHOD' (fn i =>
2497 (DETERM (reify_uneq ctxt i)
2498 THEN rule_uneq ctxt prec i
2499 THEN Simplifier.asm_full_simp_tac bounded_by_simpset i
2500 THEN (TRY (filter_prems_tac (fn t => false) i))
2501 THEN (gen_eval_tac eval_oracle ctxt) i))
2503 *} "real number approximation"