1 (* Author: Johannes Hoelzl <hoelzl@in.tum.de> 2008 / 2009 *)
3 header {* Prove unequations about real numbers by computation *}
6 imports Complex_Main Float Reflection Dense_Linear_Order Efficient_Nat
9 section "Horner Scheme"
11 subsection {* Define auxiliary helper @{text horner} function *}
13 fun horner :: "(nat \<Rightarrow> nat) \<Rightarrow> (nat \<Rightarrow> nat \<Rightarrow> nat) \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> real \<Rightarrow> real" where
14 "horner F G 0 i k x = 0" |
15 "horner F G (Suc n) i k x = 1 / real k - x * horner F G n (F i) (G i k) x"
17 lemma horner_schema': fixes x :: real and a :: "nat \<Rightarrow> real"
18 shows "a 0 - x * (\<Sum> i=0..<n. (-1)^i * a (Suc i) * x^i) = (\<Sum> i=0..<Suc n. (-1)^i * a i * x^i)"
20 have shift_pow: "\<And>i. - (x * ((-1)^i * a (Suc i) * x ^ i)) = (-1)^(Suc i) * a (Suc i) * x ^ (Suc i)" by auto
21 show ?thesis unfolding setsum_right_distrib shift_pow real_diff_def setsum_negf[symmetric] setsum_head_upt_Suc[OF zero_less_Suc]
22 setsum_reindex[OF inj_Suc, unfolded comp_def, symmetric, of "\<lambda> n. (-1)^n *a n * x^n"] by auto
25 lemma horner_schema: fixes f :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat" and F :: "nat \<Rightarrow> nat"
26 assumes f_Suc: "\<And>n. f (Suc n) = G ((F^n) s) (f n)"
27 shows "horner F G n ((F^j') s) (f j') x = (\<Sum> j = 0..< n. -1^j * (1 / real (f (j' + j))) * x^j)"
28 proof (induct n arbitrary: i k j')
31 show ?case unfolding horner.simps Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc]
32 using horner_schema'[of "\<lambda> j. 1 / real (f (j' + j))"] by auto
36 assumes "0 \<le> Ifloat x" and f_Suc: "\<And>n. f (Suc n) = G ((F^n) s) (f n)"
37 and lb_0: "\<And> i k x. lb 0 i k x = 0"
38 and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 (int k) - x * (ub n (F i) (G i k) x)"
39 and ub_0: "\<And> i k x. ub 0 i k x = 0"
40 and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 (int k) - x * (lb n (F i) (G i k) x)"
41 shows "Ifloat (lb n ((F^j') s) (f j') x) \<le> horner F G n ((F^j') s) (f j') (Ifloat x) \<and>
42 horner F G n ((F^j') s) (f j') (Ifloat x) \<le> Ifloat (ub n ((F^j') s) (f j') x)"
43 (is "?lb n j' \<le> ?horner n j' \<and> ?horner n j' \<le> ?ub n j'")
44 proof (induct n arbitrary: j')
45 case 0 thus ?case unfolding lb_0 ub_0 horner.simps by auto
48 have "?lb (Suc n) j' \<le> ?horner (Suc n) j'" unfolding lb_Suc ub_Suc horner.simps Ifloat_sub diff_def
50 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
51 from Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc, THEN conjunct2] `0 \<le> Ifloat x`
52 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))"
53 unfolding Ifloat_mult neg_le_iff_le by (rule mult_left_mono)
55 moreover have "?horner (Suc n) j' \<le> ?ub (Suc n) j'" unfolding ub_Suc ub_Suc horner.simps Ifloat_sub diff_def
57 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
58 from Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc, THEN conjunct1] `0 \<le> Ifloat x`
59 show "- (Ifloat x * horner F G n (F ((F ^ j') s)) (G ((F ^ j') s) (f j')) (Ifloat x)) \<le>
60 - Ifloat (x * lb n (F ((F ^ j') s)) (G ((F ^ j') s) (f j')) x)"
61 unfolding Ifloat_mult neg_le_iff_le by (rule mult_left_mono)
63 ultimately show ?case by blast
66 subsection "Theorems for floating point functions implementing the horner scheme"
70 Here @{term_type "f :: nat \<Rightarrow> nat"} is the sequence defining the Taylor series, the coefficients are
71 all alternating and reciprocs. We use @{term G} and @{term F} to describe the computation of @{term f}.
75 lemma horner_bounds: fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat"
76 assumes "0 \<le> Ifloat x" and f_Suc: "\<And>n. f (Suc n) = G ((F^n) s) (f n)"
77 and lb_0: "\<And> i k x. lb 0 i k x = 0"
78 and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 (int k) - x * (ub n (F i) (G i k) x)"
79 and ub_0: "\<And> i k x. ub 0 i k x = 0"
80 and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 (int k) - x * (lb n (F i) (G i k) x)"
81 shows "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
82 "(\<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")
85 using horner_bounds'[where lb=lb, OF `0 \<le> Ifloat x` f_Suc lb_0 lb_Suc ub_0 ub_Suc]
86 unfolding horner_schema[where f=f, OF f_Suc] .
87 thus "?lb" and "?ub" by auto
90 lemma horner_bounds_nonpos: fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat"
91 assumes "Ifloat x \<le> 0" and f_Suc: "\<And>n. f (Suc n) = G ((F^n) s) (f n)"
92 and lb_0: "\<And> i k x. lb 0 i k x = 0"
93 and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 (int k) + x * (ub n (F i) (G i k) x)"
94 and ub_0: "\<And> i k x. ub 0 i k x = 0"
95 and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 (int k) + x * (lb n (F i) (G i k) x)"
96 shows "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
97 "(\<Sum>j=0..<n. (1 / real (f (j' + j))) * (Ifloat x)^j) \<le> Ifloat (ub n ((F^j') s) (f j') x)" (is "?ub")
99 { fix x y z :: float have "x - y * z = x + - y * z"
100 by (cases x, cases y, cases z, simp add: plus_float.simps minus_float.simps uminus_float.simps times_float.simps algebra_simps)
101 } note diff_mult_minus = this
103 { fix x :: float have "- (- x) = x" by (cases x, auto simp add: uminus_float.simps) } note minus_minus = this
105 have move_minus: "Ifloat (-x) = -1 * Ifloat x" by auto
107 have sum_eq: "(\<Sum>j=0..<n. (1 / real (f (j' + j))) * (Ifloat x)^j) =
108 (\<Sum>j = 0..<n. -1 ^ j * (1 / real (f (j' + j))) * Ifloat (- x) ^ j)"
109 proof (rule setsum_cong, simp)
110 fix j assume "j \<in> {0 ..< n}"
111 show "1 / real (f (j' + j)) * Ifloat x ^ j = -1 ^ j * (1 / real (f (j' + j))) * Ifloat (- x) ^ j"
112 unfolding move_minus power_mult_distrib real_mult_assoc[symmetric]
113 unfolding real_mult_commute unfolding real_mult_assoc[of "-1^j", symmetric] power_mult_distrib[symmetric]
117 have "0 \<le> Ifloat (-x)" using assms by auto
118 from horner_bounds[where G=G and F=F and f=f and s=s and prec=prec
119 and lb="\<lambda> n i k x. lb n i k (-x)" and ub="\<lambda> n i k x. ub n i k (-x)", unfolded lb_Suc ub_Suc diff_mult_minus,
120 OF this f_Suc lb_0 refl ub_0 refl]
121 show "?lb" and "?ub" unfolding minus_minus sum_eq
125 subsection {* Selectors for next even or odd number *}
129 The horner scheme computes alternating series. To get the upper and lower bounds we need to
130 guarantee to access a even or odd member. To do this we use @{term get_odd} and @{term get_even}.
134 definition get_odd :: "nat \<Rightarrow> nat" where
135 "get_odd n = (if odd n then n else (Suc n))"
137 definition get_even :: "nat \<Rightarrow> nat" where
138 "get_even n = (if even n then n else (Suc n))"
140 lemma get_odd[simp]: "odd (get_odd n)" unfolding get_odd_def by (cases "odd n", auto)
141 lemma get_even[simp]: "even (get_even n)" unfolding get_even_def by (cases "even n", auto)
142 lemma get_odd_ex: "\<exists> k. Suc k = get_odd n \<and> odd (Suc k)"
143 proof (cases "odd n")
144 case True hence "0 < n" by (rule odd_pos)
145 from gr0_implies_Suc[OF this] obtain k where "Suc k = n" by auto
146 thus ?thesis unfolding get_odd_def if_P[OF True] using True[unfolded `Suc k = n`[symmetric]] by blast
148 case False hence "odd (Suc n)" by auto
149 thus ?thesis unfolding get_odd_def if_not_P[OF False] by blast
152 lemma get_even_double: "\<exists>i. get_even n = 2 * i" using get_even[unfolded even_mult_two_ex] .
153 lemma get_odd_double: "\<exists>i. get_odd n = 2 * i + 1" using get_odd[unfolded odd_Suc_mult_two_ex] by auto
155 section "Power function"
157 definition float_power_bnds :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where
158 "float_power_bnds n l u = (if odd n \<or> 0 < l then (l ^ n, u ^ n)
159 else if u < 0 then (u ^ n, l ^ n)
160 else (0, (max (-l) u) ^ n))"
162 lemma float_power_bnds: assumes "(l1, u1) = float_power_bnds n l u" and "x \<in> {Ifloat l .. Ifloat u}"
163 shows "x^n \<in> {Ifloat l1..Ifloat u1}"
164 proof (cases "even n")
167 proof (cases "0 < l")
168 case True hence "odd n \<or> 0 < l" and "0 \<le> Ifloat l" unfolding less_float_def by auto
169 have u1: "u1 = u ^ n" and l1: "l1 = l ^ n" using assms unfolding float_power_bnds_def if_P[OF `odd n \<or> 0 < l`] by auto
170 have "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
171 thus ?thesis using assms `0 < l` unfolding atLeastAtMost_iff l1 u1 float_power less_float_def by auto
173 case False hence P: "\<not> (odd n \<or> 0 < l)" using `even n` by auto
175 proof (cases "u < 0")
176 case True hence "0 \<le> - Ifloat u" and "- Ifloat u \<le> - x" and "0 \<le> - x" and "-x \<le> - Ifloat l" using assms unfolding less_float_def by auto
177 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]
178 unfolding power_minus_even[OF `even n`] by auto
179 moreover have u1: "u1 = l ^ n" and l1: "l1 = u ^ n" using assms unfolding float_power_bnds_def if_not_P[OF P] if_P[OF True] by auto
180 ultimately show ?thesis using float_power by auto
183 have "\<bar>x\<bar> \<le> Ifloat (max (-l) u)"
184 proof (cases "-l \<le> u")
185 case True thus ?thesis unfolding max_def if_P[OF True] using assms unfolding le_float_def by auto
187 case False thus ?thesis unfolding max_def if_not_P[OF False] using assms unfolding le_float_def by auto
189 hence x_abs: "\<bar>x\<bar> \<le> \<bar>Ifloat (max (-l) u)\<bar>" by auto
190 have u1: "u1 = (max (-l) u) ^ n" and l1: "l1 = 0" using assms unfolding float_power_bnds_def if_not_P[OF P] if_not_P[OF False] by auto
191 show ?thesis unfolding atLeastAtMost_iff l1 u1 float_power using zero_le_even_power[OF `even n`] power_mono_even[OF `even n` x_abs] by auto
195 case False hence "odd n \<or> 0 < l" by auto
196 have u1: "u1 = u ^ n" and l1: "l1 = l ^ n" using assms unfolding float_power_bnds_def if_P[OF `odd n \<or> 0 < l`] by auto
197 have "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
198 thus ?thesis unfolding atLeastAtMost_iff l1 u1 float_power less_float_def by auto
201 lemma bnds_power: "\<forall> x l u. (l1, u1) = float_power_bnds n l u \<and> x \<in> {Ifloat l .. Ifloat u} \<longrightarrow> Ifloat l1 \<le> x^n \<and> x^n \<le> Ifloat u1"
202 using float_power_bnds by auto
204 section "Square root"
208 The square root computation is implemented as newton iteration. As first first step we use the
209 nearest power of two greater than the square root.
213 fun sqrt_iteration :: "nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
214 "sqrt_iteration prec 0 (Float m e) = Float 1 ((e + bitlen m) div 2 + 1)" |
215 "sqrt_iteration prec (Suc m) x = (let y = sqrt_iteration prec m x
216 in Float 1 -1 * (y + float_divr prec x y))"
218 definition ub_sqrt :: "nat \<Rightarrow> float \<Rightarrow> float option" where
219 "ub_sqrt prec x = (if 0 < x then Some (sqrt_iteration prec prec x) else if x < 0 then None else Some 0)"
221 definition lb_sqrt :: "nat \<Rightarrow> float \<Rightarrow> float option" where
222 "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)"
224 lemma sqrt_ub_pos_pos_1:
225 assumes "sqrt x < b" and "0 < b" and "0 < x"
226 shows "sqrt x < (b + x / b)/2"
228 from assms have "0 < (b - sqrt x) ^ 2 " by simp
229 also have "\<dots> = b ^ 2 - 2 * b * sqrt x + (sqrt x) ^ 2" by algebra
230 also have "\<dots> = b ^ 2 - 2 * b * sqrt x + x" using assms by (simp add: real_sqrt_pow2)
231 finally have "0 < b ^ 2 - 2 * b * sqrt x + x" by assumption
232 hence "0 < b / 2 - sqrt x + x / (2 * b)" using assms
233 by (simp add: field_simps power2_eq_square)
234 thus ?thesis by (simp add: field_simps)
237 lemma sqrt_iteration_bound: assumes "0 < Ifloat x"
238 shows "sqrt (Ifloat x) < Ifloat (sqrt_iteration prec n x)"
244 hence "0 < m" using float_pos_m_pos[unfolded less_float_def] assms by auto
245 hence "0 < sqrt (real m)" by auto
247 have int_nat_bl: "int (nat (bitlen m)) = bitlen m" using bitlen_ge0 by auto
249 have "Ifloat x = (real m / 2^nat (bitlen m)) * pow2 (e + int (nat (bitlen m)))"
250 unfolding pow2_add pow2_int Float Ifloat.simps by auto
251 also have "\<dots> < 1 * pow2 (e + int (nat (bitlen m)))"
252 proof (rule mult_strict_right_mono, auto)
253 show "real m < 2^nat (bitlen m)" using bitlen_bounds[OF `0 < m`, THEN conjunct2]
254 unfolding real_of_int_less_iff[of m, symmetric] by auto
256 finally have "sqrt (Ifloat x) < sqrt (pow2 (e + bitlen m))" unfolding int_nat_bl by auto
257 also have "\<dots> \<le> pow2 ((e + bitlen m) div 2 + 1)"
259 let ?E = "e + bitlen m"
260 have E_mod_pow: "pow2 (?E mod 2) < 4"
261 proof (cases "?E mod 2 = 1")
262 case True thus ?thesis by auto
265 have "0 \<le> ?E mod 2" by auto
266 have "?E mod 2 < 2" by auto
267 from this[THEN zless_imp_add1_zle]
268 have "?E mod 2 \<le> 0" using False by auto
269 from xt1(5)[OF `0 \<le> ?E mod 2` this]
272 hence "sqrt (pow2 (?E mod 2)) < sqrt (2 * 2)" by auto
273 hence E_mod_pow: "sqrt (pow2 (?E mod 2)) < 2" unfolding real_sqrt_abs2 by auto
275 have E_eq: "pow2 ?E = pow2 (?E div 2 + ?E div 2 + ?E mod 2)" by auto
276 have "sqrt (pow2 ?E) = sqrt (pow2 (?E div 2) * pow2 (?E div 2) * pow2 (?E mod 2))"
277 unfolding E_eq unfolding pow2_add ..
278 also have "\<dots> = pow2 (?E div 2) * sqrt (pow2 (?E mod 2))"
279 unfolding real_sqrt_mult[of _ "pow2 (?E mod 2)"] real_sqrt_abs2 by auto
280 also have "\<dots> < pow2 (?E div 2) * 2"
281 by (rule mult_strict_left_mono, auto intro: E_mod_pow)
282 also have "\<dots> = pow2 (?E div 2 + 1)" unfolding zadd_commute[of _ 1] pow2_add1 by auto
283 finally show ?thesis by auto
286 unfolding Float sqrt_iteration.simps Ifloat.simps by auto
290 let ?b = "sqrt_iteration prec n x"
291 have "0 < sqrt (Ifloat x)" using `0 < Ifloat x` by auto
292 also have "\<dots> < Ifloat ?b" using Suc .
293 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
294 also have "\<dots> \<le> (Ifloat ?b + Ifloat (float_divr prec x ?b))/2" by (rule divide_right_mono, auto simp add: float_divr)
295 also have "\<dots> = Ifloat (Float 1 -1) * (Ifloat ?b + Ifloat (float_divr prec x ?b))" by auto
296 finally show ?case unfolding sqrt_iteration.simps Let_def Ifloat_mult Ifloat_add right_distrib .
299 lemma sqrt_iteration_lower_bound: assumes "0 < Ifloat x"
300 shows "0 < Ifloat (sqrt_iteration prec n x)" (is "0 < ?sqrt")
302 have "0 < sqrt (Ifloat x)" using assms by auto
303 also have "\<dots> < ?sqrt" using sqrt_iteration_bound[OF assms] .
304 finally show ?thesis .
307 lemma lb_sqrt_lower_bound: assumes "0 \<le> Ifloat x"
308 shows "0 \<le> Ifloat (the (lb_sqrt prec x))"
309 proof (cases "0 < x")
310 case True hence "0 < Ifloat x" and "0 \<le> x" using `0 \<le> Ifloat x` unfolding less_float_def le_float_def by auto
311 hence "0 < sqrt_iteration prec prec x" unfolding less_float_def using sqrt_iteration_lower_bound by auto
312 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
313 thus ?thesis unfolding lb_sqrt_def using True by auto
315 case False with `0 \<le> Ifloat x` have "Ifloat x = 0" unfolding less_float_def by auto
316 thus ?thesis unfolding lb_sqrt_def less_float_def by auto
319 lemma lb_sqrt_upper_bound: assumes "0 \<le> Ifloat x"
320 shows "Ifloat (the (lb_sqrt prec x)) \<le> sqrt (Ifloat x)"
321 proof (cases "0 < x")
322 case True hence "0 < Ifloat x" and "0 \<le> Ifloat x" unfolding less_float_def by auto
323 hence sqrt_gt0: "0 < sqrt (Ifloat x)" by auto
324 hence sqrt_ub: "sqrt (Ifloat x) < Ifloat (sqrt_iteration prec prec x)" using sqrt_iteration_bound by auto
326 have "Ifloat (float_divl prec x (sqrt_iteration prec prec x)) \<le> Ifloat x / Ifloat (sqrt_iteration prec prec x)" by (rule float_divl)
327 also have "\<dots> < Ifloat x / sqrt (Ifloat x)"
328 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]])
329 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
330 finally show ?thesis unfolding lb_sqrt_def if_P[OF `0 < x`] by auto
332 case False with `0 \<le> Ifloat x`
333 have "\<not> x < 0" unfolding less_float_def le_float_def by auto
334 show ?thesis unfolding lb_sqrt_def if_not_P[OF False] if_not_P[OF `\<not> x < 0`] using assms by auto
337 lemma lb_sqrt: assumes "Some y = lb_sqrt prec x"
338 shows "Ifloat y \<le> sqrt (Ifloat x)" and "0 \<le> Ifloat x"
340 show "0 \<le> Ifloat x"
342 assume "\<not> 0 \<le> Ifloat x"
343 hence "lb_sqrt prec x = None" unfolding lb_sqrt_def less_float_def by auto
344 thus False using assms by auto
346 from lb_sqrt_upper_bound[OF this, of prec]
347 show "Ifloat y \<le> sqrt (Ifloat x)" unfolding assms[symmetric] by auto
350 lemma ub_sqrt_lower_bound: assumes "0 \<le> Ifloat x"
351 shows "sqrt (Ifloat x) \<le> Ifloat (the (ub_sqrt prec x))"
352 proof (cases "0 < x")
353 case True hence "0 < Ifloat x" unfolding less_float_def by auto
354 hence "0 < sqrt (Ifloat x)" by auto
355 hence "sqrt (Ifloat x) < Ifloat (sqrt_iteration prec prec x)" using sqrt_iteration_bound by auto
356 thus ?thesis unfolding ub_sqrt_def if_P[OF `0 < x`] by auto
358 case False with `0 \<le> Ifloat x`
359 have "Ifloat x = 0" unfolding less_float_def le_float_def by auto
360 thus ?thesis unfolding ub_sqrt_def less_float_def le_float_def by auto
363 lemma ub_sqrt: assumes "Some y = ub_sqrt prec x"
364 shows "sqrt (Ifloat x) \<le> Ifloat y" and "0 \<le> Ifloat x"
366 show "0 \<le> Ifloat x"
368 assume "\<not> 0 \<le> Ifloat x"
369 hence "ub_sqrt prec x = None" unfolding ub_sqrt_def less_float_def by auto
370 thus False using assms by auto
372 from ub_sqrt_lower_bound[OF this, of prec]
373 show "sqrt (Ifloat x) \<le> Ifloat y" unfolding assms[symmetric] by auto
376 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"
377 proof (rule allI, rule allI, rule allI, rule impI)
379 assume "(Some l, Some u) = (lb_sqrt prec lx, ub_sqrt prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux}"
380 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
382 have "Ifloat lx \<le> x" and "x \<le> Ifloat ux" using x by auto
384 from lb_sqrt(1)[OF l] real_sqrt_le_mono[OF `Ifloat lx \<le> x`]
385 have "Ifloat l \<le> sqrt x" by (rule order_trans)
387 from real_sqrt_le_mono[OF `x \<le> Ifloat ux`] ub_sqrt(1)[OF u]
388 have "sqrt x \<le> Ifloat u" by (rule order_trans)
389 ultimately show "Ifloat l \<le> sqrt x \<and> sqrt x \<le> Ifloat u" ..
392 section "Arcus tangens and \<pi>"
394 subsection "Compute arcus tangens series"
398 As first step we implement the computation of the arcus tangens series. This is only valid in the range
399 @{term "{-1 :: real .. 1}"}. This is used to compute \<pi> and then the entire arcus tangens.
403 fun ub_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
404 and lb_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
405 "ub_arctan_horner prec 0 k x = 0"
406 | "ub_arctan_horner prec (Suc n) k x =
407 (rapprox_rat prec 1 (int k)) - x * (lb_arctan_horner prec n (k + 2) x)"
408 | "lb_arctan_horner prec 0 k x = 0"
409 | "lb_arctan_horner prec (Suc n) k x =
410 (lapprox_rat prec 1 (int k)) - x * (ub_arctan_horner prec n (k + 2) x)"
412 lemma arctan_0_1_bounds': assumes "0 \<le> Ifloat x" "Ifloat x \<le> 1" and "even n"
413 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))}"
415 let "?c i" = "-1^i * (1 / real (i * 2 + 1) * Ifloat x ^ (i * 2 + 1))"
416 let "?S n" = "\<Sum> i=0..<n. ?c i"
418 have "0 \<le> Ifloat (x * x)" by auto
419 from `even n` obtain m where "2 * m = n" unfolding even_mult_two_ex by auto
421 have "arctan (Ifloat x) \<in> { ?S n .. ?S (Suc n) }"
422 proof (cases "Ifloat x = 0")
424 hence "0 < Ifloat x" using `0 \<le> Ifloat x` by auto
425 hence prem: "0 < 1 / real (0 * 2 + (1::nat)) * Ifloat x ^ (0 * 2 + 1)" by auto
427 have "\<bar> Ifloat x \<bar> \<le> 1" using `0 \<le> Ifloat x` `Ifloat x \<le> 1` by auto
428 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`]
429 show ?thesis unfolding arctan_series[OF `\<bar> Ifloat x \<bar> \<le> 1`] Suc_plus1 .
431 note arctan_bounds = this[unfolded atLeastAtMost_iff]
433 have F: "\<And>n. 2 * Suc n + 1 = 2 * n + 1 + 2" by auto
435 note bounds = horner_bounds[where s=1 and f="\<lambda>i. 2 * i + 1" and j'=0
436 and lb="\<lambda>n i k x. lb_arctan_horner prec n k x"
437 and ub="\<lambda>n i k x. ub_arctan_horner prec n k x",
438 OF `0 \<le> Ifloat (x*x)` F lb_arctan_horner.simps ub_arctan_horner.simps]
440 { have "Ifloat (x * lb_arctan_horner prec n 1 (x*x)) \<le> ?S n"
441 using bounds(1) `0 \<le> Ifloat x`
442 unfolding Ifloat_mult power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]
443 unfolding real_mult_commute mult_commute[of _ "2::nat"] power_mult power2_eq_square[of "Ifloat x"]
444 by (auto intro!: mult_left_mono)
445 also have "\<dots> \<le> arctan (Ifloat x)" using arctan_bounds ..
446 finally have "Ifloat (x * lb_arctan_horner prec n 1 (x*x)) \<le> arctan (Ifloat x)" . }
448 { have "arctan (Ifloat x) \<le> ?S (Suc n)" using arctan_bounds ..
449 also have "\<dots> \<le> Ifloat (x * ub_arctan_horner prec (Suc n) 1 (x*x))"
450 using bounds(2)[of "Suc n"] `0 \<le> Ifloat x`
451 unfolding Ifloat_mult power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]
452 unfolding real_mult_commute mult_commute[of _ "2::nat"] power_mult power2_eq_square[of "Ifloat x"]
453 by (auto intro!: mult_left_mono)
454 finally have "arctan (Ifloat x) \<le> Ifloat (x * ub_arctan_horner prec (Suc n) 1 (x*x))" . }
455 ultimately show ?thesis by auto
458 lemma arctan_0_1_bounds: assumes "0 \<le> Ifloat x" "Ifloat x \<le> 1"
459 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))}"
460 proof (cases "even n")
462 obtain n' where "Suc n' = get_odd n" and "odd (Suc n')" using get_odd_ex by auto
463 hence "even n'" unfolding even_nat_Suc by auto
464 have "arctan (Ifloat x) \<le> Ifloat (x * ub_arctan_horner prec (get_odd n) 1 (x * x))"
465 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
467 have "Ifloat (x * lb_arctan_horner prec (get_even n) 1 (x * x)) \<le> arctan (Ifloat x)"
468 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
469 ultimately show ?thesis by auto
471 case False hence "0 < n" by (rule odd_pos)
472 from gr0_implies_Suc[OF this] obtain n' where "n = Suc n'" ..
473 from False[unfolded this even_nat_Suc]
474 have "even n'" and "even (Suc (Suc n'))" by auto
475 have "get_odd n = Suc n'" unfolding get_odd_def if_P[OF False] using `n = Suc n'` .
477 have "arctan (Ifloat x) \<le> Ifloat (x * ub_arctan_horner prec (get_odd n) 1 (x * x))"
478 unfolding `get_odd n = Suc n'` using arctan_0_1_bounds'[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1` `even n'`] by auto
480 have "Ifloat (x * lb_arctan_horner prec (get_even n) 1 (x * x)) \<le> arctan (Ifloat x)"
481 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
482 ultimately show ?thesis by auto
485 subsection "Compute \<pi>"
487 definition ub_pi :: "nat \<Rightarrow> float" where
488 "ub_pi prec = (let A = rapprox_rat prec 1 5 ;
489 B = lapprox_rat prec 1 239
490 in ((Float 1 2) * ((Float 1 2) * A * (ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (A * A)) -
491 B * (lb_arctan_horner prec (get_even (prec div 14 + 1)) 1 (B * B)))))"
493 definition lb_pi :: "nat \<Rightarrow> float" where
494 "lb_pi prec = (let A = lapprox_rat prec 1 5 ;
495 B = rapprox_rat prec 1 239
496 in ((Float 1 2) * ((Float 1 2) * A * (lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (A * A)) -
497 B * (ub_arctan_horner prec (get_odd (prec div 14 + 1)) 1 (B * B)))))"
499 lemma pi_boundaries: "pi \<in> {Ifloat (lb_pi n) .. Ifloat (ub_pi n)}"
501 have machin_pi: "pi = 4 * (4 * arctan (1 / 5) - arctan (1 / 239))" unfolding machin[symmetric] by auto
503 { fix prec n :: nat fix k :: int assume "1 < k" hence "0 \<le> k" and "0 < k" and "1 \<le> k" by auto
504 let ?k = "rapprox_rat prec 1 k"
505 have "1 div k = 0" using div_pos_pos_trivial[OF _ `1 < k`] by auto
507 have "0 \<le> Ifloat ?k" by (rule order_trans[OF _ rapprox_rat], auto simp add: `0 \<le> k`)
508 have "Ifloat ?k \<le> 1" unfolding rapprox_rat.simps(2)[OF zero_le_one `0 < k`]
509 by (rule rapprox_posrat_le1, auto simp add: `0 < k` `1 \<le> k`)
511 have "1 / real k \<le> Ifloat ?k" using rapprox_rat[where x=1 and y=k] by auto
512 hence "arctan (1 / real k) \<le> arctan (Ifloat ?k)" by (rule arctan_monotone')
513 also have "\<dots> \<le> Ifloat (?k * ub_arctan_horner prec (get_odd n) 1 (?k * ?k))"
514 using arctan_0_1_bounds[OF `0 \<le> Ifloat ?k` `Ifloat ?k \<le> 1`] by auto
515 finally have "arctan (1 / (real k)) \<le> Ifloat (?k * ub_arctan_horner prec (get_odd n) 1 (?k * ?k))" .
516 } note ub_arctan = this
518 { fix prec n :: nat fix k :: int assume "1 < k" hence "0 \<le> k" and "0 < k" by auto
519 let ?k = "lapprox_rat prec 1 k"
520 have "1 div k = 0" using div_pos_pos_trivial[OF _ `1 < k`] by auto
521 have "1 / real k \<le> 1" using `1 < k` by auto
523 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`)
524 have "\<And>n. Ifloat ?k \<le> 1" using lapprox_rat by (rule order_trans, auto simp add: `1 / real k \<le> 1`)
526 have "Ifloat ?k \<le> 1 / real k" using lapprox_rat[where x=1 and y=k] by auto
528 have "Ifloat (?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k)) \<le> arctan (Ifloat ?k)"
529 using arctan_0_1_bounds[OF `0 \<le> Ifloat ?k` `Ifloat ?k \<le> 1`] by auto
530 also have "\<dots> \<le> arctan (1 / real k)" using `Ifloat ?k \<le> 1 / real k` by (rule arctan_monotone')
531 finally have "Ifloat (?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k)) \<le> arctan (1 / (real k))" .
532 } note lb_arctan = this
534 have "pi \<le> Ifloat (ub_pi n)"
535 unfolding ub_pi_def machin_pi Let_def Ifloat_mult Ifloat_sub unfolding Float_num
536 using lb_arctan[of 239] ub_arctan[of 5]
537 by (auto intro!: mult_left_mono add_mono simp add: diff_minus simp del: lapprox_rat.simps rapprox_rat.simps)
539 have "Ifloat (lb_pi n) \<le> pi"
540 unfolding lb_pi_def machin_pi Let_def Ifloat_mult Ifloat_sub Float_num
541 using lb_arctan[of 5] ub_arctan[of 239]
542 by (auto intro!: mult_left_mono add_mono simp add: diff_minus simp del: lapprox_rat.simps rapprox_rat.simps)
543 ultimately show ?thesis by auto
546 subsection "Compute arcus tangens in the entire domain"
548 function lb_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" and ub_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" where
549 "lb_arctan prec x = (let ub_horner = \<lambda> x. x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x) ;
550 lb_horner = \<lambda> x. x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)
551 in (if x < 0 then - ub_arctan prec (-x) else
552 if x \<le> Float 1 -1 then lb_horner x else
553 if x \<le> Float 1 1 then Float 1 1 * lb_horner (float_divl prec x (1 + the (ub_sqrt prec (1 + x * x))))
554 else (let inv = float_divr prec 1 x
556 else lb_pi prec * Float 1 -1 - ub_horner inv)))"
558 | "ub_arctan prec x = (let lb_horner = \<lambda> x. x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x) ;
559 ub_horner = \<lambda> x. x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)
560 in (if x < 0 then - lb_arctan prec (-x) else
561 if x \<le> Float 1 -1 then ub_horner x else
562 if x \<le> Float 1 1 then let y = float_divr prec x (1 + the (lb_sqrt prec (1 + x * x)))
563 in if y > 1 then ub_pi prec * Float 1 -1
564 else Float 1 1 * ub_horner y
565 else ub_pi prec * Float 1 -1 - lb_horner (float_divl prec 1 x)))"
566 by pat_completeness auto
567 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)
569 declare ub_arctan_horner.simps[simp del]
570 declare lb_arctan_horner.simps[simp del]
572 lemma lb_arctan_bound': assumes "0 \<le> Ifloat x"
573 shows "Ifloat (lb_arctan prec x) \<le> arctan (Ifloat x)"
575 have "\<not> x < 0" and "0 \<le> x" unfolding less_float_def le_float_def using `0 \<le> Ifloat x` by auto
576 let "?ub_horner x" = "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)"
577 and "?lb_horner x" = "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)"
580 proof (cases "x \<le> Float 1 -1")
581 case True hence "Ifloat x \<le> 1" unfolding le_float_def Float_num by auto
582 show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True]
583 using arctan_0_1_bounds[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1`] by auto
585 case False hence "0 < Ifloat x" unfolding le_float_def Float_num by auto
586 let ?R = "1 + sqrt (1 + Ifloat x * Ifloat x)"
587 let ?fR = "1 + the (ub_sqrt prec (1 + x * x))"
588 let ?DIV = "float_divl prec x ?fR"
590 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
591 hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg)
593 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)
594 hence "?R \<le> Ifloat ?fR" by auto
595 hence "0 < ?fR" and "0 < Ifloat ?fR" unfolding less_float_def using `0 < ?R` by auto
597 have monotone: "Ifloat (float_divl prec x ?fR) \<le> Ifloat x / ?R"
599 have "Ifloat ?DIV \<le> Ifloat x / Ifloat ?fR" by (rule float_divl)
600 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]])
601 finally show ?thesis .
605 proof (cases "x \<le> Float 1 1")
608 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
609 also have "\<dots> \<le> Ifloat (the (ub_sqrt prec (1 + x * x)))" by (rule ub_sqrt_lower_bound, auto simp add: sqr_ge0)
610 finally have "Ifloat x \<le> Ifloat ?fR" by auto
611 moreover have "Ifloat ?DIV \<le> Ifloat x / Ifloat ?fR" by (rule float_divl)
612 ultimately have "Ifloat ?DIV \<le> 1" unfolding divide_le_eq_1_pos[OF `0 < Ifloat ?fR`, symmetric] by auto
614 have "0 \<le> Ifloat ?DIV" using float_divl_lower_bound[OF `0 \<le> x` `0 < ?fR`] unfolding le_float_def by auto
616 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
617 using arctan_0_1_bounds[OF `0 \<le> Ifloat ?DIV` `Ifloat ?DIV \<le> 1`] by auto
618 also have "\<dots> \<le> 2 * arctan (Ifloat x / ?R)"
619 using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono)
620 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 .
621 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] .
624 hence "2 < Ifloat x" unfolding le_float_def Float_num by auto
625 hence "1 \<le> Ifloat x" by auto
627 let "?invx" = "float_divr prec 1 x"
628 have "0 \<le> arctan (Ifloat x)" using arctan_monotone'[OF `0 \<le> Ifloat x`] using arctan_tan[of 0, unfolded tan_zero] by auto
631 proof (cases "1 < ?invx")
633 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]
634 using `0 \<le> arctan (Ifloat x)` by auto
637 hence "Ifloat ?invx \<le> 1" unfolding less_float_def by auto
638 have "0 \<le> Ifloat ?invx" by (rule order_trans[OF _ float_divr], auto simp add: `0 \<le> Ifloat x`)
640 have "1 / Ifloat x \<noteq> 0" and "0 < 1 / Ifloat x" using `0 < Ifloat x` by auto
642 have "arctan (1 / Ifloat x) \<le> arctan (Ifloat ?invx)" unfolding Ifloat_1[symmetric] by (rule arctan_monotone', rule float_divr)
643 also have "\<dots> \<le> Ifloat (?ub_horner ?invx)" using arctan_0_1_bounds[OF `0 \<le> Ifloat ?invx` `Ifloat ?invx \<le> 1`] by auto
644 finally have "pi / 2 - Ifloat (?ub_horner ?invx) \<le> arctan (Ifloat x)"
645 using `0 \<le> arctan (Ifloat x)` arctan_inverse[OF `1 / Ifloat x \<noteq> 0`]
646 unfolding real_sgn_pos[OF `0 < 1 / Ifloat x`] le_diff_eq by auto
648 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
650 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]
657 lemma ub_arctan_bound': assumes "0 \<le> Ifloat x"
658 shows "arctan (Ifloat x) \<le> Ifloat (ub_arctan prec x)"
660 have "\<not> x < 0" and "0 \<le> x" unfolding less_float_def le_float_def using `0 \<le> Ifloat x` by auto
662 let "?ub_horner x" = "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)"
663 and "?lb_horner x" = "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)"
666 proof (cases "x \<le> Float 1 -1")
667 case True hence "Ifloat x \<le> 1" unfolding le_float_def Float_num by auto
668 show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True]
669 using arctan_0_1_bounds[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1`] by auto
671 case False hence "0 < Ifloat x" unfolding le_float_def Float_num by auto
672 let ?R = "1 + sqrt (1 + Ifloat x * Ifloat x)"
673 let ?fR = "1 + the (lb_sqrt prec (1 + x * x))"
674 let ?DIV = "float_divr prec x ?fR"
676 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
677 hence "0 \<le> Ifloat (1 + x*x)" by auto
679 hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg)
681 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)
682 hence "Ifloat ?fR \<le> ?R" by auto
683 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)`])
685 have monotone: "Ifloat x / ?R \<le> Ifloat (float_divr prec x ?fR)"
687 from divide_left_mono[OF `Ifloat ?fR \<le> ?R` `0 \<le> Ifloat x` mult_pos_pos[OF divisor_gt0 `0 < Ifloat ?fR`]]
688 have "Ifloat x / ?R \<le> Ifloat x / Ifloat ?fR" .
689 also have "\<dots> \<le> Ifloat ?DIV" by (rule float_divr)
690 finally show ?thesis .
694 proof (cases "x \<le> Float 1 1")
697 proof (cases "?DIV > 1")
699 have "pi / 2 \<le> Ifloat (ub_pi prec * Float 1 -1)" unfolding Ifloat_mult Float_num times_divide_eq_right real_mult_1 using pi_boundaries by auto
700 from order_less_le_trans[OF arctan_ubound this, THEN less_imp_le]
701 show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_P[OF `x \<le> Float 1 1`] if_P[OF True] .
704 hence "Ifloat ?DIV \<le> 1" unfolding less_float_def by auto
706 have "0 \<le> Ifloat x / ?R" using `0 \<le> Ifloat x` `0 < ?R` unfolding real_0_le_divide_iff by auto
707 hence "0 \<le> Ifloat ?DIV" using monotone by (rule order_trans)
709 have "arctan (Ifloat x) = 2 * arctan (Ifloat x / ?R)" using arctan_half unfolding numeral_2_eq_2 power_Suc2 power_0 real_mult_1 .
710 also have "\<dots> \<le> 2 * arctan (Ifloat ?DIV)"
711 using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono)
712 also have "\<dots> \<le> Ifloat (Float 1 1 * ?ub_horner ?DIV)" unfolding Ifloat_mult[of "Float 1 1"] Float_num
713 using arctan_0_1_bounds[OF `0 \<le> Ifloat ?DIV` `Ifloat ?DIV \<le> 1`] by auto
714 finally show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_P[OF `x \<le> Float 1 1`] if_not_P[OF False] .
718 hence "2 < Ifloat x" unfolding le_float_def Float_num by auto
719 hence "1 \<le> Ifloat x" by auto
720 hence "0 < Ifloat x" by auto
721 hence "0 < x" unfolding less_float_def by auto
723 let "?invx" = "float_divl prec 1 x"
724 have "0 \<le> arctan (Ifloat x)" using arctan_monotone'[OF `0 \<le> Ifloat x`] using arctan_tan[of 0, unfolded tan_zero] by auto
726 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`])
727 have "0 \<le> Ifloat ?invx" unfolding Ifloat_0[symmetric] by (rule float_divl_lower_bound[unfolded le_float_def], auto simp add: `0 < x`)
729 have "1 / Ifloat x \<noteq> 0" and "0 < 1 / Ifloat x" using `0 < Ifloat x` by auto
731 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
732 also have "\<dots> \<le> arctan (1 / Ifloat x)" unfolding Ifloat_1[symmetric] by (rule arctan_monotone', rule float_divl)
733 finally have "arctan (Ifloat x) \<le> pi / 2 - Ifloat (?lb_horner ?invx)"
734 using `0 \<le> arctan (Ifloat x)` arctan_inverse[OF `1 / Ifloat x \<noteq> 0`]
735 unfolding real_sgn_pos[OF `0 < 1 / Ifloat x`] le_diff_eq by auto
737 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
739 show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_not_P[OF `\<not> x \<le> Float 1 1`] if_not_P[OF False]
745 lemma arctan_boundaries:
746 "arctan (Ifloat x) \<in> {Ifloat (lb_arctan prec x) .. Ifloat (ub_arctan prec x)}"
747 proof (cases "0 \<le> x")
748 case True hence "0 \<le> Ifloat x" unfolding le_float_def by auto
749 show ?thesis using ub_arctan_bound'[OF `0 \<le> Ifloat x`] lb_arctan_bound'[OF `0 \<le> Ifloat x`] unfolding atLeastAtMost_iff by auto
752 case False hence "x < 0" and "0 \<le> Ifloat ?mx" unfolding le_float_def less_float_def by auto
753 hence bounds: "Ifloat (lb_arctan prec ?mx) \<le> arctan (Ifloat ?mx) \<and> arctan (Ifloat ?mx) \<le> Ifloat (ub_arctan prec ?mx)"
754 using ub_arctan_bound'[OF `0 \<le> Ifloat ?mx`] lb_arctan_bound'[OF `0 \<le> Ifloat ?mx`] by auto
755 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`]
756 unfolding atLeastAtMost_iff using bounds[unfolded Ifloat_minus arctan_minus] by auto
759 lemma bnds_arctan: "\<forall> x lx ux. (l, u) = (lb_arctan prec lx, ub_arctan prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux} \<longrightarrow> Ifloat l \<le> arctan x \<and> arctan x \<le> Ifloat u"
760 proof (rule allI, rule allI, rule allI, rule impI)
762 assume "(l, u) = (lb_arctan prec lx, ub_arctan prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux}"
763 hence l: "lb_arctan prec lx = l " and u: "ub_arctan prec ux = u" and x: "x \<in> {Ifloat lx .. Ifloat ux}" by auto
765 { from arctan_boundaries[of lx prec, unfolded l]
766 have "Ifloat l \<le> arctan (Ifloat lx)" by (auto simp del: lb_arctan.simps)
767 also have "\<dots> \<le> arctan x" using x by (auto intro: arctan_monotone')
768 finally have "Ifloat l \<le> arctan x" .
770 { have "arctan x \<le> arctan (Ifloat ux)" using x by (auto intro: arctan_monotone')
771 also have "\<dots> \<le> Ifloat u" using arctan_boundaries[of ux prec, unfolded u] by (auto simp del: ub_arctan.simps)
772 finally have "arctan x \<le> Ifloat u" .
773 } ultimately show "Ifloat l \<le> arctan x \<and> arctan x \<le> Ifloat u" ..
776 section "Sinus and Cosinus"
778 subsection "Compute the cosinus and sinus series"
780 fun ub_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
781 and lb_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
782 "ub_sin_cos_aux prec 0 i k x = 0"
783 | "ub_sin_cos_aux prec (Suc n) i k x =
784 (rapprox_rat prec 1 (int k)) - x * (lb_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)"
785 | "lb_sin_cos_aux prec 0 i k x = 0"
786 | "lb_sin_cos_aux prec (Suc n) i k x =
787 (lapprox_rat prec 1 (int k)) - x * (ub_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)"
790 shows "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")
791 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")
793 have "0 \<le> Ifloat (x * x)" unfolding Ifloat_mult by auto
794 let "?f n" = "fact (2 * n)"
797 have F: "\<And>m. ((\<lambda>i. i + 2) ^ n) m = m + 2 * n" by (induct n arbitrary: m, auto)
798 have "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) ^ n) 1 * (((\<lambda>i. i + 2) ^ n) 1 + 1)"
799 unfolding F by auto } note f_eq = this
801 from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0,
802 OF `0 \<le> Ifloat (x * x)` f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps]
803 show "?lb" and "?ub" by (auto simp add: power_mult power2_eq_square[of "Ifloat x"])
806 lemma cos_boundaries: assumes "0 \<le> Ifloat x" and "Ifloat x \<le> pi / 2"
807 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))}"
808 proof (cases "Ifloat x = 0")
809 case False hence "Ifloat x \<noteq> 0" by auto
810 hence "0 < x" and "0 < Ifloat x" using `0 \<le> Ifloat x` unfolding less_float_def by auto
811 have "0 < x * x" using `0 < x` unfolding less_float_def Ifloat_mult Ifloat_0
812 using mult_pos_pos[where a="Ifloat x" and b="Ifloat x"] by auto
814 { fix x n have "(\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i))) * x^(2 * i))
815 = (\<Sum> i = 0 ..< 2 * n. (if even(i) then (-1 ^ (i div 2))/(real (fact i)) else 0) * x ^ i)" (is "?sum = ?ifsum")
817 have "?sum = ?sum + (\<Sum> j = 0 ..< n. 0)" by auto
819 (\<Sum> j = 0 ..< n. -1 ^ ((2 * j) div 2) / (real (fact (2 * j))) * x ^(2 * j)) + (\<Sum> j = 0 ..< n. 0)" by auto
820 also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. if even i then -1 ^ (i div 2) / (real (fact i)) * x ^ i else 0)"
821 unfolding sum_split_even_odd ..
822 also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. (if even i then -1 ^ (i div 2) / (real (fact i)) else 0) * x ^ i)"
823 by (rule setsum_cong2) auto
824 finally show ?thesis by assumption
825 qed } note morph_to_if_power = this
828 { fix n :: nat assume "0 < n"
829 hence "0 < 2 * n" by auto
830 obtain t where "0 < t" and "t < Ifloat x" and
831 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)
832 + (cos (t + 1/2 * real (2 * n) * pi) / real (fact (2*n))) * (Ifloat x)^(2*n)"
833 (is "_ = ?SUM + ?rest / ?fact * ?pow")
834 using Maclaurin_cos_expansion2[OF `0 < Ifloat x` `0 < 2 * n`] by auto
836 have "cos t * -1^n = cos t * cos (real n * pi) + sin t * sin (real n * pi)" by auto
837 also have "\<dots> = cos (t + real n * pi)" using cos_add by auto
838 also have "\<dots> = ?rest" by auto
839 finally have "cos t * -1^n = ?rest" .
841 have "t \<le> pi / 2" using `t < Ifloat x` and `Ifloat x \<le> pi / 2` by auto
842 hence "0 \<le> cos t" using `0 < t` and cos_ge_zero by auto
843 ultimately have even: "even n \<Longrightarrow> 0 \<le> ?rest" and odd: "odd n \<Longrightarrow> 0 \<le> - ?rest " by auto
845 have "0 < ?fact" by auto
846 have "0 < ?pow" using `0 < Ifloat x` by auto
850 have "Ifloat (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> ?SUM"
851 unfolding morph_to_if_power[symmetric] using cos_aux by auto
852 also have "\<dots> \<le> cos (Ifloat x)"
854 from even[OF `even n`] `0 < ?fact` `0 < ?pow`
855 have "0 \<le> (?rest / ?fact) * ?pow" by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
856 thus ?thesis unfolding cos_eq by auto
858 finally have "Ifloat (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> cos (Ifloat x)" .
863 have "cos (Ifloat x) \<le> ?SUM"
865 from `0 < ?fact` and `0 < ?pow` and odd[OF `odd n`]
866 have "0 \<le> (- ?rest) / ?fact * ?pow"
867 by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
868 thus ?thesis unfolding cos_eq by auto
870 also have "\<dots> \<le> Ifloat (ub_sin_cos_aux prec n 1 1 (x * x))"
871 unfolding morph_to_if_power[symmetric] using cos_aux by auto
872 finally have "cos (Ifloat x) \<le> Ifloat (ub_sin_cos_aux prec n 1 1 (x * x))" .
873 } note ub = this and lb
874 } note ub = this(1) and lb = this(2)
876 have "cos (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] .
877 moreover have "Ifloat (lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) \<le> cos (Ifloat x)"
878 proof (cases "0 < get_even n")
879 case True show ?thesis using lb[OF True get_even] .
882 hence "get_even n = 0" by auto
883 have "- (pi / 2) \<le> Ifloat x" by (rule order_trans[OF _ `0 < Ifloat x`[THEN less_imp_le]], auto)
884 with `Ifloat x \<le> pi / 2`
885 show ?thesis unfolding `get_even n = 0` lb_sin_cos_aux.simps Ifloat_minus Ifloat_0 using cos_ge_zero by auto
887 ultimately show ?thesis by auto
891 proof (cases "n = 0")
893 thus ?thesis unfolding `n = 0` get_even_def get_odd_def using `Ifloat x = 0` lapprox_rat[where x="-1" and y=1] by auto
895 case False with not0_implies_Suc obtain m where "n = Suc m" by blast
896 thus ?thesis unfolding `n = Suc m` get_even_def get_odd_def using `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)
900 lemma sin_aux: assumes "0 \<le> Ifloat x"
901 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")
902 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")
904 have "0 \<le> Ifloat (x * x)" unfolding Ifloat_mult by auto
905 let "?f n" = "fact (2 * n + 1)"
908 have F: "\<And>m. ((\<lambda>i. i + 2) ^ n) m = m + 2 * n" by (induct n arbitrary: m, auto)
909 have "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) ^ n) 2 * (((\<lambda>i. i + 2) ^ n) 2 + 1)"
910 unfolding F by auto } note f_eq = this
912 from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0,
913 OF `0 \<le> Ifloat (x * x)` f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps]
914 show "?lb" and "?ub" using `0 \<le> Ifloat x` unfolding Ifloat_mult
915 unfolding power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]
916 unfolding real_mult_commute
917 by (auto intro!: mult_left_mono simp add: power_mult power2_eq_square[of "Ifloat x"])
920 lemma sin_boundaries: assumes "0 \<le> Ifloat x" and "Ifloat x \<le> pi / 2"
921 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))}"
922 proof (cases "Ifloat x = 0")
923 case False hence "Ifloat x \<noteq> 0" by auto
924 hence "0 < x" and "0 < Ifloat x" using `0 \<le> Ifloat x` unfolding less_float_def by auto
925 have "0 < x * x" using `0 < x` unfolding less_float_def Ifloat_mult Ifloat_0
926 using mult_pos_pos[where a="Ifloat x" and b="Ifloat x"] by auto
928 { fix x n have "(\<Sum> j = 0 ..< n. -1 ^ (((2 * j + 1) - Suc 0) div 2) / (real (fact (2 * j + 1))) * x ^(2 * j + 1))
929 = (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * x ^ i)" (is "?SUM = _")
931 have pow: "!!i. x ^ (2 * i + 1) = x * x ^ (2 * i)" by auto
932 have "?SUM = (\<Sum> j = 0 ..< n. 0) + ?SUM" by auto
933 also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. if even i then 0 else -1 ^ ((i - Suc 0) div 2) / (real (fact i)) * x ^ i)"
934 unfolding sum_split_even_odd ..
935 also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. (if even i then 0 else -1 ^ ((i - Suc 0) div 2) / (real (fact i))) * x ^ i)"
936 by (rule setsum_cong2) auto
937 finally show ?thesis by assumption
938 qed } note setsum_morph = this
940 { fix n :: nat assume "0 < n"
941 hence "0 < 2 * n + 1" by auto
942 obtain t where "0 < t" and "t < Ifloat x" and
943 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)
944 + (sin (t + 1/2 * real (2 * n + 1) * pi) / real (fact (2*n + 1))) * (Ifloat x)^(2*n + 1)"
945 (is "_ = ?SUM + ?rest / ?fact * ?pow")
946 using Maclaurin_sin_expansion3[OF `0 < 2 * n + 1` `0 < Ifloat x`] by auto
948 have "?rest = cos t * -1^n" unfolding sin_add cos_add real_of_nat_add left_distrib right_distrib by auto
950 have "t \<le> pi / 2" using `t < Ifloat x` and `Ifloat x \<le> pi / 2` by auto
951 hence "0 \<le> cos t" using `0 < t` and cos_ge_zero by auto
952 ultimately have even: "even n \<Longrightarrow> 0 \<le> ?rest" and odd: "odd n \<Longrightarrow> 0 \<le> - ?rest " by auto
954 have "0 < ?fact" by (rule real_of_nat_fact_gt_zero)
955 have "0 < ?pow" using `0 < Ifloat x` by (rule zero_less_power)
959 have "Ifloat (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le>
960 (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * (Ifloat x) ^ i)"
961 using sin_aux[OF `0 \<le> Ifloat x`] unfolding setsum_morph[symmetric] by auto
962 also have "\<dots> \<le> ?SUM" by auto
963 also have "\<dots> \<le> sin (Ifloat x)"
965 from even[OF `even n`] `0 < ?fact` `0 < ?pow`
966 have "0 \<le> (?rest / ?fact) * ?pow" by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
967 thus ?thesis unfolding sin_eq by auto
969 finally have "Ifloat (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> sin (Ifloat x)" .
974 have "sin (Ifloat x) \<le> ?SUM"
976 from `0 < ?fact` and `0 < ?pow` and odd[OF `odd n`]
977 have "0 \<le> (- ?rest) / ?fact * ?pow"
978 by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
979 thus ?thesis unfolding sin_eq by auto
981 also have "\<dots> \<le> (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * (Ifloat x) ^ i)"
983 also have "\<dots> \<le> Ifloat (x * ub_sin_cos_aux prec n 2 1 (x * x))"
984 using sin_aux[OF `0 \<le> Ifloat x`] unfolding setsum_morph[symmetric] by auto
985 finally have "sin (Ifloat x) \<le> Ifloat (x * ub_sin_cos_aux prec n 2 1 (x * x))" .
986 } note ub = this and lb
987 } note ub = this(1) and lb = this(2)
989 have "sin (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] .
990 moreover have "Ifloat (x * lb_sin_cos_aux prec (get_even n) 2 1 (x * x)) \<le> sin (Ifloat x)"
991 proof (cases "0 < get_even n")
992 case True show ?thesis using lb[OF True get_even] .
995 hence "get_even n = 0" by auto
996 with `Ifloat x \<le> pi / 2` `0 \<le> Ifloat x`
997 show ?thesis unfolding `get_even n = 0` ub_sin_cos_aux.simps Ifloat_minus Ifloat_0 using sin_ge_zero by auto
999 ultimately show ?thesis by auto
1003 proof (cases "n = 0")
1005 thus ?thesis unfolding `n = 0` get_even_def get_odd_def using `Ifloat x = 0` lapprox_rat[where x="-1" and y=1] by auto
1007 case False with not0_implies_Suc obtain m where "n = Suc m" by blast
1008 thus ?thesis unfolding `n = Suc m` get_even_def get_odd_def using `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)
1012 subsection "Compute the cosinus in the entire domain"
1014 definition lb_cos :: "nat \<Rightarrow> float \<Rightarrow> float" where
1015 "lb_cos prec x = (let
1016 horner = \<lambda> x. lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 1 1 (x * x) ;
1017 half = \<lambda> x. if x < 0 then - 1 else Float 1 1 * x * x - 1
1018 in if x < Float 1 -1 then horner x
1019 else if x < 1 then half (horner (x * Float 1 -1))
1020 else half (half (horner (x * Float 1 -2))))"
1022 definition ub_cos :: "nat \<Rightarrow> float \<Rightarrow> float" where
1023 "ub_cos prec x = (let
1024 horner = \<lambda> x. ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x) ;
1025 half = \<lambda> x. Float 1 1 * x * x - 1
1026 in if x < Float 1 -1 then horner x
1027 else if x < 1 then half (horner (x * Float 1 -1))
1028 else half (half (horner (x * Float 1 -2))))"
1030 definition bnds_cos :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where
1031 "bnds_cos prec lx ux = (let lpi = lb_pi prec
1032 in if lx < -lpi \<or> ux > lpi then (Float -1 0, Float 1 0)
1033 else if ux \<le> 0 then (lb_cos prec (-lx), ub_cos prec (-ux))
1034 else if 0 \<le> lx then (lb_cos prec ux, ub_cos prec lx)
1035 else (min (lb_cos prec (-lx)) (lb_cos prec ux), Float 1 0))"
1037 lemma lb_cos: assumes "0 \<le> Ifloat x" and "Ifloat x \<le> pi"
1038 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) }")
1041 have "cos x = cos (x / 2 + x / 2)" by auto
1042 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"
1043 unfolding cos_add by auto
1044 also have "\<dots> = 2 * cos (x / 2) * cos (x / 2) - 1" by algebra
1045 finally have "cos x = 2 * cos (x / 2) * cos (x / 2) - 1" .
1046 } note x_half = this[symmetric]
1048 have "\<not> x < 0" using `0 \<le> Ifloat x` unfolding less_float_def by auto
1049 let "?ub_horner x" = "ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x)"
1050 let "?lb_horner x" = "lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 1 1 (x * x)"
1051 let "?ub_half x" = "Float 1 1 * x * x - 1"
1052 let "?lb_half x" = "if x < 0 then - 1 else Float 1 1 * x * x - 1"
1055 proof (cases "x < Float 1 -1")
1056 case True hence "Ifloat x \<le> pi / 2" unfolding less_float_def using pi_ge_two by auto
1057 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
1058 using cos_boundaries[OF `0 \<le> Ifloat x` `Ifloat x \<le> pi / 2`] .
1062 { fix y x :: float let ?x2 = "Ifloat (x * Float 1 -1)"
1063 assume "Ifloat y \<le> cos ?x2" and "-pi \<le> Ifloat x" and "Ifloat x \<le> pi"
1064 hence "- (pi / 2) \<le> ?x2" and "?x2 \<le> pi / 2" using pi_ge_two unfolding Ifloat_mult Float_num by auto
1065 hence "0 \<le> cos ?x2" by (rule cos_ge_zero)
1067 have "Ifloat (?lb_half y) \<le> cos (Ifloat x)"
1068 proof (cases "y < 0")
1069 case True show ?thesis using cos_ge_minus_one unfolding if_P[OF True] by auto
1072 hence "0 \<le> Ifloat y" unfolding less_float_def by auto
1073 from mult_mono[OF `Ifloat y \<le> cos ?x2` `Ifloat y \<le> cos ?x2` `0 \<le> cos ?x2` this]
1074 have "Ifloat y * Ifloat y \<le> cos ?x2 * cos ?x2" .
1075 hence "2 * Ifloat y * Ifloat y \<le> 2 * cos ?x2 * cos ?x2" by auto
1076 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
1077 thus ?thesis unfolding if_not_P[OF False] x_half Float_num Ifloat_mult Ifloat_sub by auto
1079 } note lb_half = this
1081 { fix y x :: float let ?x2 = "Ifloat (x * Float 1 -1)"
1082 assume ub: "cos ?x2 \<le> Ifloat y" and "- pi \<le> Ifloat x" and "Ifloat x \<le> pi"
1083 hence "- (pi / 2) \<le> ?x2" and "?x2 \<le> pi / 2" using pi_ge_two unfolding Ifloat_mult Float_num by auto
1084 hence "0 \<le> cos ?x2" by (rule cos_ge_zero)
1086 have "cos (Ifloat x) \<le> Ifloat (?ub_half y)"
1088 have "0 \<le> Ifloat y" using `0 \<le> cos ?x2` ub by (rule order_trans)
1089 from mult_mono[OF ub ub this `0 \<le> cos ?x2`]
1090 have "cos ?x2 * cos ?x2 \<le> Ifloat y * Ifloat y" .
1091 hence "2 * cos ?x2 * cos ?x2 \<le> 2 * Ifloat y * Ifloat y" by auto
1092 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
1093 thus ?thesis unfolding x_half Ifloat_mult Float_num Ifloat_sub by auto
1095 } note ub_half = this
1097 let ?x2 = "x * Float 1 -1"
1098 let ?x4 = "x * Float 1 -1 * Float 1 -1"
1100 have "-pi \<le> Ifloat x" using pi_ge_zero[THEN le_imp_neg_le, unfolded minus_zero] `0 \<le> Ifloat x` by (rule order_trans)
1103 proof (cases "x < 1")
1104 case True hence "Ifloat x \<le> 1" unfolding less_float_def by auto
1105 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
1106 from cos_boundaries[OF this]
1107 have lb: "Ifloat (?lb_horner ?x2) \<le> ?cos ?x2" and ub: "?cos ?x2 \<le> Ifloat (?ub_horner ?x2)" by auto
1109 have "Ifloat (?lb x) \<le> ?cos x"
1111 from lb_half[OF lb `-pi \<le> Ifloat x` `Ifloat x \<le> pi`]
1112 show ?thesis unfolding lb_cos_def[where x=x] Let_def using `\<not> x < 0` `\<not> x < Float 1 -1` `x < 1` by auto
1114 moreover have "?cos x \<le> Ifloat (?ub x)"
1116 from ub_half[OF ub `-pi \<le> Ifloat x` `Ifloat x \<le> pi`]
1117 show ?thesis unfolding ub_cos_def[where x=x] Let_def using `\<not> x < 0` `\<not> x < Float 1 -1` `x < 1` by auto
1119 ultimately show ?thesis by auto
1122 have "0 \<le> 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
1123 from cos_boundaries[OF this]
1124 have lb: "Ifloat (?lb_horner ?x4) \<le> ?cos ?x4" and ub: "?cos ?x4 \<le> Ifloat (?ub_horner ?x4)" by auto
1126 have eq_4: "?x2 * Float 1 -1 = x * Float 1 -2" by (cases x, auto simp add: times_float.simps)
1128 have "Ifloat (?lb x) \<le> ?cos x"
1130 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
1131 from lb_half[OF lb_half[OF lb this] `-pi \<le> Ifloat x` `Ifloat x \<le> pi`, unfolded eq_4]
1132 show ?thesis unfolding lb_cos_def[where x=x] if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x < Float 1 -1`] if_not_P[OF `\<not> x < 1`] Let_def .
1134 moreover have "?cos x \<le> Ifloat (?ub x)"
1136 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
1137 from ub_half[OF ub_half[OF ub this] `-pi \<le> Ifloat x` `Ifloat x \<le> pi`, unfolded eq_4]
1138 show ?thesis unfolding ub_cos_def[where x=x] if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x < Float 1 -1`] if_not_P[OF `\<not> x < 1`] Let_def .
1140 ultimately show ?thesis by auto
1145 lemma lb_cos_minus: assumes "-pi \<le> Ifloat x" and "Ifloat x \<le> 0"
1146 shows "cos (Ifloat (-x)) \<in> {Ifloat (lb_cos prec (-x)) .. Ifloat (ub_cos prec (-x))}"
1148 have "0 \<le> Ifloat (-x)" and "Ifloat (-x) \<le> pi" using `-pi \<le> Ifloat x` `Ifloat x \<le> 0` by auto
1149 from lb_cos[OF this] show ?thesis .
1152 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"
1153 proof (rule allI, rule allI, rule allI, rule impI)
1155 assume "(l, u) = bnds_cos prec lx ux \<and> x \<in> {Ifloat lx .. Ifloat ux}"
1156 hence bnds: "(l, u) = bnds_cos prec lx ux" and x: "x \<in> {Ifloat lx .. Ifloat ux}" by auto
1158 let ?lpi = "lb_pi prec"
1159 have [intro!]: "Ifloat lx \<le> Ifloat ux" using x by auto
1160 hence "lx \<le> ux" unfolding le_float_def .
1162 show "Ifloat l \<le> cos x \<and> cos x \<le> Ifloat u"
1163 proof (cases "lx < -?lpi \<or> ux > ?lpi")
1165 show ?thesis using bnds unfolding bnds_cos_def if_P[OF True] Let_def using cos_le_one cos_ge_minus_one by auto
1167 case False note not_out = this
1168 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
1170 from pi_boundaries[unfolded atLeastAtMost_iff, THEN conjunct1, THEN le_imp_neg_le] lpi_lx
1171 have "- pi \<le> Ifloat lx" by (rule order_trans)
1172 hence "- pi \<le> x" and "- pi \<le> Ifloat ux" and "x \<le> Ifloat ux" using x by auto
1174 from lpi_ux pi_boundaries[unfolded atLeastAtMost_iff, THEN conjunct1]
1175 have "Ifloat ux \<le> pi" by (rule order_trans)
1176 hence "x \<le> pi" and "Ifloat lx \<le> pi" and "Ifloat lx \<le> x" using x by auto
1178 note lb_cos_minus_bottom = lb_cos_minus[unfolded atLeastAtMost_iff, THEN conjunct1]
1179 note lb_cos_minus_top = lb_cos_minus[unfolded atLeastAtMost_iff, THEN conjunct2]
1180 note lb_cos_bottom = lb_cos[unfolded atLeastAtMost_iff, THEN conjunct1]
1181 note lb_cos_top = lb_cos[unfolded atLeastAtMost_iff, THEN conjunct2]
1184 proof (cases "ux \<le> 0")
1185 case True hence "Ifloat ux \<le> 0" unfolding le_float_def by auto
1186 hence "x \<le> 0" and "Ifloat lx \<le> 0" using x by auto
1188 { have "Ifloat (lb_cos prec (-lx)) \<le> cos (Ifloat (-lx))" using lb_cos_minus_bottom[OF `-pi \<le> Ifloat lx` `Ifloat lx \<le> 0`] .
1189 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`] .
1190 finally have "Ifloat (lb_cos prec (-lx)) \<le> cos x" . }
1192 { 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`] .
1193 also have "\<dots> \<le> Ifloat (ub_cos prec (-ux))" using lb_cos_minus_top[OF `-pi \<le> Ifloat ux` `Ifloat ux \<le> 0`] .
1194 finally have "cos x \<le> Ifloat (ub_cos prec (-ux))" . }
1195 ultimately show ?thesis using bnds unfolding bnds_cos_def Let_def if_not_P[OF not_out] if_P[OF True] by auto
1197 case False note not_ux = this
1200 proof (cases "0 \<le> lx")
1201 case True hence "0 \<le> Ifloat lx" unfolding le_float_def by auto
1202 hence "0 \<le> x" and "0 \<le> Ifloat ux" using x by auto
1204 { have "Ifloat (lb_cos prec ux) \<le> cos (Ifloat ux)" using lb_cos_bottom[OF `0 \<le> Ifloat ux` `Ifloat ux \<le> pi`] .
1205 also have "\<dots> \<le> cos x" using cos_monotone_0_pi'[OF `0 \<le> x` `x \<le> Ifloat ux` `Ifloat ux \<le> pi`] .
1206 finally have "Ifloat (lb_cos prec ux) \<le> cos x" . }
1208 { have "cos x \<le> cos (Ifloat lx)" using cos_monotone_0_pi'[OF `0 \<le> Ifloat lx` `Ifloat lx \<le> x` `x \<le> pi`] .
1209 also have "\<dots> \<le> Ifloat (ub_cos prec lx)" using lb_cos_top[OF `0 \<le> Ifloat lx` `Ifloat lx \<le> pi`] .
1210 finally have "cos x \<le> Ifloat (ub_cos prec lx)" . }
1211 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
1213 case False with not_ux
1214 have "Ifloat lx \<le> 0" and "0 \<le> Ifloat ux" unfolding le_float_def by auto
1216 have "Ifloat (min (lb_cos prec (-lx)) (lb_cos prec ux)) \<le> cos x"
1217 proof (cases "x \<le> 0")
1219 have "Ifloat (lb_cos prec (-lx)) \<le> cos (Ifloat (-lx))" using lb_cos_minus_bottom[OF `-pi \<le> Ifloat lx` `Ifloat lx \<le> 0`] .
1220 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`] .
1221 finally show ?thesis unfolding Ifloat_min by auto
1223 case False hence "0 \<le> x" by auto
1224 have "Ifloat (lb_cos prec ux) \<le> cos (Ifloat ux)" using lb_cos_bottom[OF `0 \<le> Ifloat ux` `Ifloat ux \<le> pi`] .
1225 also have "\<dots> \<le> cos x" using cos_monotone_0_pi'[OF `0 \<le> x` `x \<le> Ifloat ux` `Ifloat ux \<le> pi`] .
1226 finally show ?thesis unfolding Ifloat_min by auto
1228 moreover have "cos x \<le> Ifloat (Float 1 0)" by auto
1229 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
1235 subsection "Compute the sinus in the entire domain"
1237 function lb_sin :: "nat \<Rightarrow> float \<Rightarrow> float" and ub_sin :: "nat \<Rightarrow> float \<Rightarrow> float" where
1238 "lb_sin prec x = (let sqr_diff = \<lambda> x. if x > 1 then 0 else 1 - x * x
1239 in if x < 0 then - ub_sin prec (- x)
1240 else if x \<le> Float 1 -1 then x * lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 2 1 (x * x)
1241 else the (lb_sqrt prec (sqr_diff (ub_cos prec x))))" |
1243 "ub_sin prec x = (let sqr_diff = \<lambda> x. if x < 0 then 1 else 1 - x * x
1244 in if x < 0 then - lb_sin prec (- x)
1245 else if x \<le> Float 1 -1 then x * ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 2 1 (x * x)
1246 else the (ub_sqrt prec (sqr_diff (lb_cos prec x))))"
1247 by pat_completeness auto
1248 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)
1250 definition bnds_sin :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where
1251 "bnds_sin prec lx ux = (let
1253 half_pi = lpi * Float 1 -1
1254 in if lx \<le> - half_pi \<or> half_pi \<le> ux then (Float -1 0, Float 1 0)
1255 else (lb_sin prec lx, ub_sin prec ux))"
1257 lemma lb_sin: assumes "- (pi / 2) \<le> Ifloat x" and "Ifloat x \<le> pi / 2"
1258 shows "sin (Ifloat x) \<in> { Ifloat (lb_sin prec x) .. Ifloat (ub_sin prec x) }" (is "?sin x \<in> { ?lb x .. ?ub x}")
1260 { fix x :: float assume "0 \<le> Ifloat x" and "Ifloat x \<le> pi / 2"
1261 hence "\<not> (x < 0)" and "- (pi / 2) \<le> Ifloat x" unfolding less_float_def using pi_ge_two by auto
1263 have "Ifloat x \<le> pi" using `Ifloat x \<le> pi / 2` using pi_ge_two by auto
1265 have "?sin x \<in> { ?lb x .. ?ub x}"
1266 proof (cases "x \<le> Float 1 -1")
1267 case True from sin_boundaries[OF `0 \<le> Ifloat x` `Ifloat x \<le> pi / 2`]
1268 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 .
1271 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
1272 have "0 \<le> sin (Ifloat x)" using `0 \<le> Ifloat x` and `Ifloat x \<le> pi / 2` using sin_ge_zero by auto
1274 have "?sin x \<le> ?ub x"
1275 proof (cases "lb_cos prec x < 0")
1277 have "?sin x \<le> 1" using sin_le_one .
1278 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
1279 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 .
1281 case False hence "0 \<le> Ifloat (lb_cos prec x)" unfolding less_float_def by auto
1283 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
1284 also have "\<dots> \<le> sqrt (Ifloat (1 - lb_cos prec x * lb_cos prec x))"
1285 proof (rule real_sqrt_le_mono)
1286 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
1287 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)
1288 thus "1 - cos (Ifloat x) ^ 2 \<le> Ifloat (1 - lb_cos prec x * lb_cos prec x)" unfolding Ifloat_sub Ifloat_1 by auto
1290 also have "\<dots> \<le> Ifloat (the (ub_sqrt prec (1 - lb_cos prec x * lb_cos prec x)))"
1291 proof (rule ub_sqrt_lower_bound)
1292 have "Ifloat (lb_cos prec x) \<le> cos (Ifloat x)" using lb_cos[OF `0 \<le> Ifloat x` `Ifloat x \<le> pi`] by auto
1293 from mult_mono[OF order_trans[OF this cos_le_one] order_trans[OF this cos_le_one]]
1294 have "Ifloat (lb_cos prec x) * Ifloat (lb_cos prec x) \<le> 1" using `0 \<le> Ifloat (lb_cos prec x)` by auto
1295 thus "0 \<le> Ifloat (1 - lb_cos prec x * lb_cos prec x)" by auto
1297 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 .
1300 have "?lb x \<le> ?sin x"
1301 proof (cases "1 < ub_cos prec x")
1303 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
1304 by (rule order_trans[OF _ sin_ge_zero[OF `0 \<le> Ifloat x` `Ifloat x \<le> pi`]])
1305 (auto simp add: lb_sqrt_upper_bound[where prec=prec and x=0, unfolded Ifloat_0 real_sqrt_zero])
1307 case False hence "Ifloat (ub_cos prec x) \<le> 1" unfolding less_float_def by auto
1308 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
1310 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))"
1311 proof (rule lb_sqrt_upper_bound)
1312 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)`
1313 have "Ifloat (ub_cos prec x) * Ifloat (ub_cos prec x) \<le> 1" by auto
1314 thus "0 \<le> Ifloat (1 - ub_cos prec x * ub_cos prec x)" by auto
1316 also have "\<dots> \<le> sqrt (1 - cos (Ifloat x) ^ 2)"
1317 proof (rule real_sqrt_le_mono)
1318 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
1319 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)
1320 thus "Ifloat (1 - ub_cos prec x * ub_cos prec x) \<le> 1 - cos (Ifloat x) ^ 2" unfolding Ifloat_sub Ifloat_1 by auto
1322 also have "\<dots> = sin (Ifloat x)" unfolding sin_squared_eq[symmetric] real_sqrt_abs using `0 \<le> sin (Ifloat x)` by auto
1323 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 .
1325 ultimately show ?thesis by auto
1327 } note for_pos = this
1330 proof (cases "x < 0")
1332 hence "0 \<le> Ifloat (-x)" and "Ifloat (- x) \<le> pi / 2" using `-(pi/2) \<le> Ifloat x` unfolding less_float_def by auto
1333 from for_pos[OF this]
1334 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
1336 case False hence "0 \<le> Ifloat x" unfolding less_float_def by auto
1337 from for_pos[OF this `Ifloat x \<le> pi /2`]
1342 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"
1343 proof (rule allI, rule allI, rule allI, rule impI)
1345 assume "(l, u) = bnds_sin prec lx ux \<and> x \<in> {Ifloat lx .. Ifloat ux}"
1346 hence bnds: "(l, u) = bnds_sin prec lx ux" and x: "x \<in> {Ifloat lx .. Ifloat ux}" by auto
1347 show "Ifloat l \<le> sin x \<and> sin x \<le> Ifloat u"
1348 proof (cases "lx \<le> - (lb_pi prec * Float 1 -1) \<or> lb_pi prec * Float 1 -1 \<le> ux")
1349 case True show ?thesis using bnds unfolding bnds_sin_def if_P[OF True] Let_def by auto
1352 hence "- lb_pi prec * Float 1 -1 \<le> lx" and "ux \<le> lb_pi prec * Float 1 -1" unfolding le_float_def by auto
1353 moreover have "Ifloat (lb_pi prec * Float 1 -1) \<le> pi / 2" unfolding Ifloat_mult using pi_boundaries by auto
1354 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
1355 hence "- (pi / 2) \<le> Ifloat ux" and "Ifloat lx \<le> pi / 2" by auto
1357 have "- (pi / 2) \<le> x""x \<le> pi / 2" using `Ifloat ux \<le> pi / 2` `- (pi /2) \<le> Ifloat lx` x by auto
1359 { 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
1360 also have "\<dots> \<le> sin x" using sin_monotone_2pi' `- (pi / 2) \<le> Ifloat lx` x `x \<le> pi / 2` by auto
1361 finally have "Ifloat (lb_sin prec lx) \<le> sin x" . }
1363 { have "sin x \<le> sin (Ifloat ux)" using sin_monotone_2pi' `- (pi / 2) \<le> x` x `Ifloat ux \<le> pi / 2` by auto
1364 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
1365 finally have "sin x \<le> Ifloat (ub_sin prec ux)" . }
1367 show ?thesis using bnds unfolding bnds_sin_def if_not_P[OF False] Let_def by auto
1371 section "Exponential function"
1373 subsection "Compute the series of the exponential function"
1375 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
1376 "ub_exp_horner prec 0 i k x = 0" |
1377 "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" |
1378 "lb_exp_horner prec 0 i k x = 0" |
1379 "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"
1381 lemma bnds_exp_horner: assumes "Ifloat x \<le> 0"
1382 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) }"
1385 have F: "\<And> m. ((\<lambda>i. i + 1) ^ n) m = n + m" by (induct n, auto)
1386 have "fact (Suc n) = fact n * ((\<lambda>i. i + 1) ^ n) 1" unfolding F by auto } note f_eq = this
1388 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,
1389 OF assms f_eq lb_exp_horner.simps ub_exp_horner.simps]
1391 { 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)"
1392 using bounds(1) by auto
1393 also have "\<dots> \<le> exp (Ifloat x)"
1395 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)"
1396 using Maclaurin_exp_le by blast
1397 moreover have "0 \<le> exp t / real (fact (get_even n)) * (Ifloat x) ^ (get_even n)"
1398 by (auto intro!: mult_nonneg_nonneg divide_nonneg_pos simp add: get_even zero_le_even_power exp_gt_zero)
1399 ultimately show ?thesis
1400 using get_odd exp_gt_zero by (auto intro!: pordered_cancel_semiring_class.mult_nonneg_nonneg)
1402 finally have "Ifloat (lb_exp_horner prec (get_even n) 1 1 x) \<le> exp (Ifloat x)" .
1405 have x_less_zero: "Ifloat x ^ get_odd n \<le> 0"
1406 proof (cases "Ifloat x = 0")
1408 have "(get_odd n) \<noteq> 0" using get_odd[THEN odd_pos] by auto
1409 thus ?thesis unfolding True power_0_left by auto
1411 case False hence "Ifloat x < 0" using `Ifloat x \<le> 0` by auto
1412 show ?thesis by (rule less_imp_le, auto simp add: power_less_zero_eq get_odd `Ifloat x < 0`)
1415 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)"
1416 using Maclaurin_exp_le by blast
1417 moreover have "exp t / real (fact (get_odd n)) * (Ifloat x) ^ (get_odd n) \<le> 0"
1418 by (auto intro!: mult_nonneg_nonpos divide_nonpos_pos simp add: x_less_zero exp_gt_zero)
1419 ultimately have "exp (Ifloat x) \<le> (\<Sum>j = 0..<get_odd n. 1 / real (fact j) * Ifloat x ^ j)"
1420 using get_odd exp_gt_zero by (auto intro!: pordered_cancel_semiring_class.mult_nonneg_nonneg)
1421 also have "\<dots> \<le> Ifloat (ub_exp_horner prec (get_odd n) 1 1 x)"
1422 using bounds(2) by auto
1423 finally have "exp (Ifloat x) \<le> Ifloat (ub_exp_horner prec (get_odd n) 1 1 x)" .
1424 } ultimately show ?thesis by auto
1427 subsection "Compute the exponential function on the entire domain"
1429 function ub_exp :: "nat \<Rightarrow> float \<Rightarrow> float" and lb_exp :: "nat \<Rightarrow> float \<Rightarrow> float" where
1430 "lb_exp prec x = (if 0 < x then float_divl prec 1 (ub_exp prec (-x))
1432 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)
1433 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))
1435 "ub_exp prec x = (if 0 < x then float_divr prec 1 (lb_exp prec (-x))
1436 else if x < - 1 then (case floor_fl x of (Float m e) \<Rightarrow>
1437 (ub_exp_horner prec (get_odd (prec + 2)) 1 1 (float_divr prec x (- Float m e))) ^ (nat (-m) * 2 ^ nat e))
1438 else ub_exp_horner prec (get_odd (prec + 2)) 1 1 x)"
1439 by pat_completeness auto
1440 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)
1442 lemma exp_m1_ge_quarter: "(1 / 4 :: real) \<le> exp (- 1)"
1444 have eq4: "4 = Suc (Suc (Suc (Suc 0)))" by auto
1446 have "1 / 4 = Ifloat (Float 1 -2)" unfolding Float_num by auto
1447 also have "\<dots> \<le> Ifloat (lb_exp_horner 1 (get_even 4) 1 1 (- 1))"
1448 unfolding get_even_def eq4
1449 by (auto simp add: lapprox_posrat_def rapprox_posrat_def normfloat.simps)
1450 also have "\<dots> \<le> exp (Ifloat (- 1))" using bnds_exp_horner[where x="- 1"] by auto
1451 finally show ?thesis unfolding Ifloat_minus Ifloat_1 .
1454 lemma lb_exp_pos: assumes "\<not> 0 < x" shows "0 < lb_exp prec x"
1456 let "?lb_horner x" = "lb_exp_horner prec (get_even (prec + 2)) 1 1 x"
1457 let "?horner x" = "let y = ?lb_horner x in if y \<le> 0 then Float 1 -2 else y"
1458 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)
1459 moreover { fix x :: float fix num :: nat
1460 have "0 < Ifloat (?horner x) ^ num" using `0 < ?horner x`[unfolded less_float_def Ifloat_0] by (rule zero_less_power)
1461 also have "\<dots> = Ifloat ((?horner x) ^ num)" using float_power by auto
1462 finally have "0 < Ifloat ((?horner x) ^ num)" .
1464 ultimately show ?thesis
1465 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)
1468 lemma exp_boundaries': assumes "x \<le> 0"
1469 shows "exp (Ifloat x) \<in> { Ifloat (lb_exp prec x) .. Ifloat (ub_exp prec x)}"
1471 let "?lb_exp_horner x" = "lb_exp_horner prec (get_even (prec + 2)) 1 1 x"
1472 let "?ub_exp_horner x" = "ub_exp_horner prec (get_odd (prec + 2)) 1 1 x"
1474 have "Ifloat x \<le> 0" and "\<not> x > 0" using `x \<le> 0` unfolding le_float_def less_float_def by auto
1476 proof (cases "x < - 1")
1477 case False hence "- 1 \<le> Ifloat x" unfolding less_float_def by auto
1479 proof (cases "?lb_exp_horner x \<le> 0")
1480 from `\<not> x < - 1` have "- 1 \<le> Ifloat x" unfolding less_float_def by auto
1481 hence "exp (- 1) \<le> exp (Ifloat x)" unfolding exp_le_cancel_iff .
1482 from order_trans[OF exp_m1_ge_quarter this]
1483 have "Ifloat (Float 1 -2) \<le> exp (Ifloat x)" unfolding Float_num .
1485 ultimately show ?thesis using bnds_exp_horner `Ifloat x \<le> 0` `\<not> x > 0` `\<not> x < - 1` by auto
1487 case False thus ?thesis using bnds_exp_horner `Ifloat x \<le> 0` `\<not> x > 0` `\<not> x < - 1` by (auto simp add: Let_def)
1492 obtain m e where Float_floor: "floor_fl x = Float m e" by (cases "floor_fl x", auto)
1493 let ?num = "nat (- m) * 2 ^ nat e"
1495 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)
1496 hence "Ifloat (floor_fl x) < 0" unfolding Float_floor Ifloat.simps using zero_less_pow2[of xe] by auto
1498 unfolding less_float_def Ifloat_0 Float_floor Ifloat.simps
1499 unfolding pos_prod_lt[OF zero_less_pow2[of e], unfolded real_mult_commute] by auto
1500 hence "1 \<le> - m" by auto
1501 hence "0 < nat (- m)" by auto
1503 have "0 \<le> e" using floor_pos_exp Float_floor[symmetric] by auto
1504 hence "(0::nat) < 2 ^ nat e" by auto
1505 ultimately have "0 < ?num" by auto
1506 hence "real ?num \<noteq> 0" by auto
1507 have e_nat: "int (nat e) = e" using `0 \<le> e` by auto
1508 have num_eq: "real ?num = Ifloat (- floor_fl x)" using `0 < nat (- m)`
1509 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
1510 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 .
1511 hence "Ifloat (floor_fl x) < 0" unfolding less_float_def by auto
1513 have "exp (Ifloat x) \<le> Ifloat (ub_exp prec x)"
1515 have div_less_zero: "Ifloat (float_divr prec x (- floor_fl x)) \<le> 0"
1516 using float_divr_nonpos_pos_upper_bound[OF `x \<le> 0` `0 < - floor_fl x`] unfolding le_float_def Ifloat_0 .
1518 have "exp (Ifloat x) = exp (real ?num * (Ifloat x / real ?num))" using `real ?num \<noteq> 0` by auto
1519 also have "\<dots> = exp (Ifloat x / real ?num) ^ ?num" unfolding exp_real_of_nat_mult ..
1520 also have "\<dots> \<le> exp (Ifloat (float_divr prec x (- floor_fl x))) ^ ?num" unfolding num_eq
1521 by (rule power_mono, rule exp_le_cancel_iff[THEN iffD2], rule float_divr) auto
1522 also have "\<dots> \<le> Ifloat ((?ub_exp_horner (float_divr prec x (- floor_fl x))) ^ ?num)" unfolding float_power
1523 by (rule power_mono, rule bnds_exp_horner[OF div_less_zero, unfolded atLeastAtMost_iff, THEN conjunct2], auto)
1524 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 .
1527 have "Ifloat (lb_exp prec x) \<le> exp (Ifloat x)"
1529 let ?divl = "float_divl prec x (- Float m e)"
1530 let ?horner = "?lb_exp_horner ?divl"
1533 proof (cases "?horner \<le> 0")
1534 case False hence "0 \<le> Ifloat ?horner" unfolding le_float_def by auto
1536 have div_less_zero: "Ifloat (float_divl prec x (- floor_fl x)) \<le> 0"
1537 using `Ifloat (floor_fl x) < 0` `Ifloat x \<le> 0` by (auto intro!: order_trans[OF float_divl] divide_nonpos_neg)
1539 have "Ifloat ((?lb_exp_horner (float_divl prec x (- floor_fl x))) ^ ?num) \<le>
1540 exp (Ifloat (float_divl prec x (- floor_fl x))) ^ ?num" unfolding float_power
1541 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)
1542 also have "\<dots> \<le> exp (Ifloat x / real ?num) ^ ?num" unfolding num_eq
1543 using float_divl by (auto intro!: power_mono simp del: Ifloat_minus)
1544 also have "\<dots> = exp (real ?num * (Ifloat x / real ?num))" unfolding exp_real_of_nat_mult ..
1545 also have "\<dots> = exp (Ifloat x)" using `real ?num \<noteq> 0` by auto
1546 finally show ?thesis
1547 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
1550 have "Ifloat (floor_fl x) \<noteq> 0" and "Ifloat (floor_fl x) \<le> 0" using `Ifloat (floor_fl x) < 0` by auto
1551 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`]]
1552 have "- 1 \<le> Ifloat x / Ifloat (- floor_fl x)" unfolding Ifloat_minus by auto
1553 from order_trans[OF exp_m1_ge_quarter this[unfolded exp_le_cancel_iff[where x="- 1", symmetric]]]
1554 have "Ifloat (Float 1 -2) \<le> exp (Ifloat x / Ifloat (- floor_fl x))" unfolding Float_num .
1555 hence "Ifloat (Float 1 -2) ^ ?num \<le> exp (Ifloat x / Ifloat (- floor_fl x)) ^ ?num"
1556 by (auto intro!: power_mono simp add: Float_num)
1557 also have "\<dots> = exp (Ifloat x)" unfolding num_eq exp_real_of_nat_mult[symmetric] using `Ifloat (floor_fl x) \<noteq> 0` by auto
1558 finally show ?thesis
1559 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 .
1562 ultimately show ?thesis by auto
1566 lemma exp_boundaries: "exp (Ifloat x) \<in> { Ifloat (lb_exp prec x) .. Ifloat (ub_exp prec x)}"
1569 proof (cases "0 < x")
1570 case False hence "x \<le> 0" unfolding less_float_def le_float_def by auto
1571 from exp_boundaries'[OF this] show ?thesis .
1573 case True hence "-x \<le> 0" unfolding less_float_def le_float_def by auto
1575 have "Ifloat (lb_exp prec x) \<le> exp (Ifloat x)"
1577 from exp_boundaries'[OF `-x \<le> 0`]
1578 have ub_exp: "exp (- Ifloat x) \<le> Ifloat (ub_exp prec (-x))" unfolding atLeastAtMost_iff Ifloat_minus by auto
1580 have "Ifloat (float_divl prec 1 (ub_exp prec (-x))) \<le> Ifloat 1 / Ifloat (ub_exp prec (-x))" using float_divl .
1581 also have "Ifloat 1 / Ifloat (ub_exp prec (-x)) \<le> exp (Ifloat x)"
1582 using ub_exp[unfolded inverse_le_iff_le[OF order_less_le_trans[OF exp_gt_zero ub_exp] exp_gt_zero, symmetric]]
1583 unfolding exp_minus nonzero_inverse_inverse_eq[OF exp_not_eq_zero] inverse_eq_divide by auto
1584 finally show ?thesis unfolding lb_exp.simps if_P[OF True] .
1587 have "exp (Ifloat x) \<le> Ifloat (ub_exp prec x)"
1589 have "\<not> 0 < -x" using `0 < x` unfolding less_float_def by auto
1591 from exp_boundaries'[OF `-x \<le> 0`]
1592 have lb_exp: "Ifloat (lb_exp prec (-x)) \<le> exp (- Ifloat x)" unfolding atLeastAtMost_iff Ifloat_minus by auto
1594 have "exp (Ifloat x) \<le> Ifloat 1 / Ifloat (lb_exp prec (-x))"
1595 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]]
1596 unfolding exp_minus nonzero_inverse_inverse_eq[OF exp_not_eq_zero] inverse_eq_divide Ifloat_1 by auto
1597 also have "\<dots> \<le> Ifloat (float_divr prec 1 (lb_exp prec (-x)))" using float_divr .
1598 finally show ?thesis unfolding ub_exp.simps if_P[OF True] .
1600 ultimately show ?thesis by auto
1604 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"
1605 proof (rule allI, rule allI, rule allI, rule impI)
1607 assume "(l, u) = (lb_exp prec lx, ub_exp prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux}"
1608 hence l: "lb_exp prec lx = l " and u: "ub_exp prec ux = u" and x: "x \<in> {Ifloat lx .. Ifloat ux}" by auto
1610 { from exp_boundaries[of lx prec, unfolded l]
1611 have "Ifloat l \<le> exp (Ifloat lx)" by (auto simp del: lb_exp.simps)
1612 also have "\<dots> \<le> exp x" using x by auto
1613 finally have "Ifloat l \<le> exp x" .
1615 { have "exp x \<le> exp (Ifloat ux)" using x by auto
1616 also have "\<dots> \<le> Ifloat u" using exp_boundaries[of ux prec, unfolded u] by (auto simp del: ub_exp.simps)
1617 finally have "exp x \<le> Ifloat u" .
1618 } ultimately show "Ifloat l \<le> exp x \<and> exp x \<le> Ifloat u" ..
1623 subsection "Compute the logarithm series"
1625 fun ub_ln_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
1626 and lb_ln_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
1627 "ub_ln_horner prec 0 i x = 0" |
1628 "ub_ln_horner prec (Suc n) i x = rapprox_rat prec 1 (int i) - x * lb_ln_horner prec n (Suc i) x" |
1629 "lb_ln_horner prec 0 i x = 0" |
1630 "lb_ln_horner prec (Suc n) i x = lapprox_rat prec 1 (int i) - x * ub_ln_horner prec n (Suc i) x"
1633 assumes "0 \<le> x" and "x < 1"
1634 shows "(\<Sum>i=0..<2*n. -1^i * (1 / real (i + 1)) * x^(Suc i)) \<le> ln (x + 1)" (is "?lb")
1635 and "ln (x + 1) \<le> (\<Sum>i=0..<2*n + 1. -1^i * (1 / real (i + 1)) * x^(Suc i))" (is "?ub")
1637 let "?a n" = "(1/real (n +1)) * x^(Suc n)"
1639 have ln_eq: "(\<Sum> i. -1^i * ?a i) = ln (x + 1)"
1640 using ln_series[of "x + 1"] `0 \<le> x` `x < 1` by auto
1642 have "norm x < 1" using assms by auto
1643 have "?a ----> 0" unfolding Suc_plus1[symmetric] inverse_eq_divide[symmetric]
1644 using LIMSEQ_mult[OF LIMSEQ_inverse_real_of_nat LIMSEQ_Suc[OF LIMSEQ_power_zero[OF `norm x < 1`]]] by auto
1645 { fix n have "0 \<le> ?a n" by (rule mult_nonneg_nonneg, auto intro!: mult_nonneg_nonneg simp add: `0 \<le> x`) }
1646 { fix n have "?a (Suc n) \<le> ?a n" unfolding inverse_eq_divide[symmetric]
1647 proof (rule mult_mono)
1648 show "0 \<le> x ^ Suc (Suc n)" by (auto intro!: mult_nonneg_nonneg simp add: `0 \<le> x`)
1649 have "x ^ Suc (Suc n) \<le> x ^ Suc n * 1" unfolding power_Suc2 real_mult_assoc[symmetric]
1650 by (rule mult_left_mono, fact less_imp_le[OF `x < 1`], auto intro!: mult_nonneg_nonneg simp add: `0 \<le> x`)
1651 thus "x ^ Suc (Suc n) \<le> x ^ Suc n" by auto
1653 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]
1654 show "?lb" and "?ub" by auto
1657 lemma ln_float_bounds:
1658 assumes "0 \<le> Ifloat x" and "Ifloat x < 1"
1659 shows "Ifloat (x * lb_ln_horner prec (get_even n) 1 x) \<le> ln (Ifloat x + 1)" (is "?lb \<le> ?ln")
1660 and "ln (Ifloat x + 1) \<le> Ifloat (x * ub_ln_horner prec (get_odd n) 1 x)" (is "?ln \<le> ?ub")
1662 obtain ev where ev: "get_even n = 2 * ev" using get_even_double ..
1663 obtain od where od: "get_odd n = 2 * od + 1" using get_odd_double ..
1665 let "?s n" = "-1^n * (1 / real (1 + n)) * (Ifloat x)^(Suc n)"
1667 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
1668 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",
1669 OF `0 \<le> Ifloat x` refl lb_ln_horner.simps ub_ln_horner.simps] `0 \<le> Ifloat x`
1670 by (rule mult_right_mono)
1671 also have "\<dots> \<le> ?ln" using ln_bounds(1)[OF `0 \<le> Ifloat x` `Ifloat x < 1`] by auto
1672 finally show "?lb \<le> ?ln" .
1674 have "?ln \<le> setsum ?s {0 ..< 2 * od + 1}" using ln_bounds(2)[OF `0 \<le> Ifloat x` `Ifloat x < 1`] by auto
1675 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
1676 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",
1677 OF `0 \<le> Ifloat x` refl lb_ln_horner.simps ub_ln_horner.simps] `0 \<le> Ifloat x`
1678 by (rule mult_right_mono)
1679 finally show "?ln \<le> ?ub" .
1682 lemma ln_add: assumes "0 < x" and "0 < y" shows "ln (x + y) = ln x + ln (1 + y / x)"
1684 have "x \<noteq> 0" using assms by auto
1685 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
1687 have "0 < y / x" using assms divide_pos_pos by auto
1688 hence "0 < 1 + y / x" by auto
1689 ultimately show ?thesis using ln_mult assms by auto
1692 subsection "Compute the logarithm of 2"
1694 definition ub_ln2 where "ub_ln2 prec = (let third = rapprox_rat (max prec 1) 1 3
1695 in (Float 1 -1 * ub_ln_horner prec (get_odd prec) 1 (Float 1 -1)) +
1696 (third * ub_ln_horner prec (get_odd prec) 1 third))"
1697 definition lb_ln2 where "lb_ln2 prec = (let third = lapprox_rat prec 1 3
1698 in (Float 1 -1 * lb_ln_horner prec (get_even prec) 1 (Float 1 -1)) +
1699 (third * lb_ln_horner prec (get_even prec) 1 third))"
1701 lemma ub_ln2: "ln 2 \<le> Ifloat (ub_ln2 prec)" (is "?ub_ln2")
1702 and lb_ln2: "Ifloat (lb_ln2 prec) \<le> ln 2" (is "?lb_ln2")
1704 let ?uthird = "rapprox_rat (max prec 1) 1 3"
1705 let ?lthird = "lapprox_rat prec 1 3"
1707 have ln2_sum: "ln 2 = ln (1/2 + 1) + ln (1 / 3 + 1)"
1708 using ln_add[of "3 / 2" "1 / 2"] by auto
1709 have lb3: "Ifloat ?lthird \<le> 1 / 3" using lapprox_rat[of prec 1 3] by auto
1710 hence lb3_ub: "Ifloat ?lthird < 1" by auto
1711 have lb3_lb: "0 \<le> Ifloat ?lthird" using lapprox_rat_bottom[of 1 3] by auto
1712 have ub3: "1 / 3 \<le> Ifloat ?uthird" using rapprox_rat[of 1 3] by auto
1713 hence ub3_lb: "0 \<le> Ifloat ?uthird" by auto
1715 have lb2: "0 \<le> Ifloat (Float 1 -1)" and ub2: "Ifloat (Float 1 -1) < 1" unfolding Float_num by auto
1717 have "0 \<le> (1::int)" and "0 < (3::int)" by auto
1718 have ub3_ub: "Ifloat ?uthird < 1" unfolding rapprox_rat.simps(2)[OF `0 \<le> 1` `0 < 3`]
1719 by (rule rapprox_posrat_less1, auto)
1721 have third_gt0: "(0 :: real) < 1 / 3 + 1" by auto
1722 have uthird_gt0: "0 < Ifloat ?uthird + 1" using ub3_lb by auto
1723 have lthird_gt0: "0 < Ifloat ?lthird + 1" using lb3_lb by auto
1725 show ?ub_ln2 unfolding ub_ln2_def Let_def Ifloat_add ln2_sum Float_num(4)[symmetric]
1726 proof (rule add_mono, fact ln_float_bounds(2)[OF lb2 ub2])
1727 have "ln (1 / 3 + 1) \<le> ln (Ifloat ?uthird + 1)" unfolding ln_le_cancel_iff[OF third_gt0 uthird_gt0] using ub3 by auto
1728 also have "\<dots> \<le> Ifloat (?uthird * ub_ln_horner prec (get_odd prec) 1 ?uthird)"
1729 using ln_float_bounds(2)[OF ub3_lb ub3_ub] .
1730 finally show "ln (1 / 3 + 1) \<le> Ifloat (?uthird * ub_ln_horner prec (get_odd prec) 1 ?uthird)" .
1732 show ?lb_ln2 unfolding lb_ln2_def Let_def Ifloat_add ln2_sum Float_num(4)[symmetric]
1733 proof (rule add_mono, fact ln_float_bounds(1)[OF lb2 ub2])
1734 have "Ifloat (?lthird * lb_ln_horner prec (get_even prec) 1 ?lthird) \<le> ln (Ifloat ?lthird + 1)"
1735 using ln_float_bounds(1)[OF lb3_lb lb3_ub] .
1736 also have "\<dots> \<le> ln (1 / 3 + 1)" unfolding ln_le_cancel_iff[OF lthird_gt0 third_gt0] using lb3 by auto
1737 finally show "Ifloat (?lthird * lb_ln_horner prec (get_even prec) 1 ?lthird) \<le> ln (1 / 3 + 1)" .
1741 subsection "Compute the logarithm in the entire domain"
1743 function ub_ln :: "nat \<Rightarrow> float \<Rightarrow> float option" and lb_ln :: "nat \<Rightarrow> float \<Rightarrow> float option" where
1744 "ub_ln prec x = (if x \<le> 0 then None
1745 else if x < 1 then Some (- the (lb_ln prec (float_divl (max prec 1) 1 x)))
1746 else let horner = \<lambda>x. (x - 1) * ub_ln_horner prec (get_odd prec) 1 (x - 1) in
1747 if x < Float 1 1 then Some (horner x)
1748 else let l = bitlen (mantissa x) - 1 in
1749 Some (ub_ln2 prec * (Float (scale x + l) 0) + horner (Float (mantissa x) (- l))))" |
1750 "lb_ln prec x = (if x \<le> 0 then None
1751 else if x < 1 then Some (- the (ub_ln prec (float_divr prec 1 x)))
1752 else let horner = \<lambda>x. (x - 1) * lb_ln_horner prec (get_even prec) 1 (x - 1) in
1753 if x < Float 1 1 then Some (horner x)
1754 else let l = bitlen (mantissa x) - 1 in
1755 Some (lb_ln2 prec * (Float (scale x + l) 0) + horner (Float (mantissa x) (- l))))"
1756 by pat_completeness auto
1758 termination proof (relation "measure (\<lambda> v. let (prec, x) = sum_case id id v in (if x < 1 then 1 else 0))", auto)
1759 fix prec x assume "\<not> x \<le> 0" and "x < 1" and "float_divl (max prec (Suc 0)) 1 x < 1"
1760 hence "0 < x" and "0 < max prec (Suc 0)" unfolding less_float_def le_float_def by auto
1761 from float_divl_pos_less1_bound[OF `0 < x` `x < 1` `0 < max prec (Suc 0)`]
1762 show False using `float_divl (max prec (Suc 0)) 1 x < 1` unfolding less_float_def le_float_def by auto
1764 fix prec x assume "\<not> x \<le> 0" and "x < 1" and "float_divr prec 1 x < 1"
1765 hence "0 < x" unfolding less_float_def le_float_def by auto
1766 from float_divr_pos_less1_lower_bound[OF `0 < x` `x < 1`, of prec]
1767 show False using `float_divr prec 1 x < 1` unfolding less_float_def le_float_def by auto
1770 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))))"
1772 let ?B = "2^nat (bitlen m - 1)"
1773 have "0 < real m" and "\<And>X. (0 :: real) < 2^X" and "0 < (2 :: real)" and "m \<noteq> 0" using assms by auto
1774 hence "0 \<le> bitlen m - 1" using bitlen_ge1[OF `m \<noteq> 0`] by auto
1776 proof (cases "0 \<le> e")
1778 show ?thesis unfolding normalized_float[OF `m \<noteq> 0`]
1779 unfolding ln_div[OF `0 < real m` `0 < ?B`] real_of_int_add ln_realpow[OF `0 < 2`]
1780 unfolding Ifloat_ge0_exp[OF True] ln_mult[OF `0 < real m` `0 < 2^nat e`]
1781 ln_realpow[OF `0 < 2`] algebra_simps using `0 \<le> bitlen m - 1` True by auto
1783 case False hence "0 < -e" by auto
1784 hence pow_gt0: "(0::real) < 2^nat (-e)" by auto
1785 hence inv_gt0: "(0::real) < inverse (2^nat (-e))" by auto
1786 show ?thesis unfolding normalized_float[OF `m \<noteq> 0`]
1787 unfolding ln_div[OF `0 < real m` `0 < ?B`] real_of_int_add ln_realpow[OF `0 < 2`]
1788 unfolding Ifloat_nge0_exp[OF False] ln_mult[OF `0 < real m` inv_gt0] ln_inverse[OF pow_gt0]
1789 ln_realpow[OF `0 < 2`] algebra_simps using `0 \<le> bitlen m - 1` False by auto
1793 lemma ub_ln_lb_ln_bounds': assumes "1 \<le> x"
1794 shows "Ifloat (the (lb_ln prec x)) \<le> ln (Ifloat x) \<and> ln (Ifloat x) \<le> Ifloat (the (ub_ln prec x))"
1795 (is "?lb \<le> ?ln \<and> ?ln \<le> ?ub")
1796 proof (cases "x < Float 1 1")
1797 case True hence "Ifloat (x - 1) < 1" unfolding less_float_def Float_num by auto
1798 have "\<not> x \<le> 0" and "\<not> x < 1" using `1 \<le> x` unfolding less_float_def le_float_def by auto
1799 hence "0 \<le> Ifloat (x - 1)" using `1 \<le> x` unfolding less_float_def Float_num by auto
1800 show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps Let_def
1801 using ln_float_bounds[OF `0 \<le> Ifloat (x - 1)` `Ifloat (x - 1) < 1`] `\<not> x \<le> 0` `\<not> x < 1` True by auto
1804 have "\<not> x \<le> 0" and "\<not> x < 1" "0 < x" using `1 \<le> x` unfolding less_float_def le_float_def by auto
1808 let ?s = "Float (e + (bitlen m - 1)) 0"
1809 let ?x = "Float m (- (bitlen m - 1))"
1811 have "0 < m" and "m \<noteq> 0" using float_pos_m_pos `0 < x` Float by auto
1814 have "Ifloat (lb_ln2 prec * ?s) \<le> ln 2 * real (e + (bitlen m - 1))" (is "?lb2 \<le> _")
1815 unfolding Ifloat_mult Ifloat_ge0_exp[OF order_refl] nat_0 power_0 mult_1_right
1816 using lb_ln2[of prec]
1817 proof (rule mult_right_mono)
1818 have "1 \<le> Float m e" using `1 \<le> x` Float unfolding le_float_def by auto
1819 from float_gt1_scale[OF this]
1820 show "0 \<le> real (e + (bitlen m - 1))" by auto
1823 from bitlen_div[OF `0 < m`, unfolded normalized_float[OF `m \<noteq> 0`, symmetric]]
1824 have "0 \<le> Ifloat (?x - 1)" and "Ifloat (?x - 1) < 1" by auto
1825 from ln_float_bounds(1)[OF this]
1826 have "Ifloat ((?x - 1) * lb_ln_horner prec (get_even prec) 1 (?x - 1)) \<le> ln (Ifloat ?x)" (is "?lb_horner \<le> _") by auto
1827 ultimately have "?lb2 + ?lb_horner \<le> ln (Ifloat x)"
1828 unfolding Float ln_shifted_float[OF `0 < m`, of e] by auto
1832 from bitlen_div[OF `0 < m`, unfolded normalized_float[OF `m \<noteq> 0`, symmetric]]
1833 have "0 \<le> Ifloat (?x - 1)" and "Ifloat (?x - 1) < 1" by auto
1834 from ln_float_bounds(2)[OF this]
1835 have "ln (Ifloat ?x) \<le> Ifloat ((?x - 1) * ub_ln_horner prec (get_odd prec) 1 (?x - 1))" (is "_ \<le> ?ub_horner") by auto
1837 have "ln 2 * real (e + (bitlen m - 1)) \<le> Ifloat (ub_ln2 prec * ?s)" (is "_ \<le> ?ub2")
1838 unfolding Ifloat_mult Ifloat_ge0_exp[OF order_refl] nat_0 power_0 mult_1_right
1839 using ub_ln2[of prec]
1840 proof (rule mult_right_mono)
1841 have "1 \<le> Float m e" using `1 \<le> x` Float unfolding le_float_def by auto
1842 from float_gt1_scale[OF this]
1843 show "0 \<le> real (e + (bitlen m - 1))" by auto
1845 ultimately have "ln (Ifloat x) \<le> ?ub2 + ?ub_horner"
1846 unfolding Float ln_shifted_float[OF `0 < m`, of e] by auto
1848 ultimately show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps
1849 unfolding if_not_P[OF `\<not> x \<le> 0`] if_not_P[OF `\<not> x < 1`] if_not_P[OF False] Let_def
1850 unfolding scale.simps[of m e, unfolded Float[symmetric]] mantissa.simps[of m e, unfolded Float[symmetric]] Ifloat_add by auto
1854 lemma ub_ln_lb_ln_bounds: assumes "0 < x"
1855 shows "Ifloat (the (lb_ln prec x)) \<le> ln (Ifloat x) \<and> ln (Ifloat x) \<le> Ifloat (the (ub_ln prec x))"
1856 (is "?lb \<le> ?ln \<and> ?ln \<le> ?ub")
1857 proof (cases "x < 1")
1858 case False hence "1 \<le> x" unfolding less_float_def le_float_def by auto
1859 show ?thesis using ub_ln_lb_ln_bounds'[OF `1 \<le> x`] .
1861 case True have "\<not> x \<le> 0" using `0 < x` unfolding less_float_def le_float_def by auto
1863 have "0 < Ifloat x" and "Ifloat x \<noteq> 0" using `0 < x` unfolding less_float_def by auto
1864 hence A: "0 < 1 / Ifloat x" by auto
1867 let ?divl = "float_divl (max prec 1) 1 x"
1868 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
1869 hence B: "0 < Ifloat ?divl" unfolding le_float_def by auto
1871 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
1872 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
1873 from this ub_ln_lb_ln_bounds'[OF A', THEN conjunct1, THEN le_imp_neg_le]
1874 have "?ln \<le> Ifloat (- the (lb_ln prec ?divl))" unfolding Ifloat_minus by (rule order_trans)
1877 let ?divr = "float_divr prec 1 x"
1878 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
1879 hence B: "0 < Ifloat ?divr" unfolding le_float_def by auto
1881 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
1882 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
1883 from ub_ln_lb_ln_bounds'[OF A', THEN conjunct2, THEN le_imp_neg_le] this
1884 have "Ifloat (- the (ub_ln prec ?divr)) \<le> ?ln" unfolding Ifloat_minus by (rule order_trans)
1886 ultimately show ?thesis unfolding lb_ln.simps[where x=x] ub_ln.simps[where x=x]
1887 unfolding if_not_P[OF `\<not> x \<le> 0`] if_P[OF True] by auto
1890 lemma lb_ln: assumes "Some y = lb_ln prec x"
1891 shows "Ifloat y \<le> ln (Ifloat x)" and "0 < Ifloat x"
1895 assume "\<not> 0 < x" hence "x \<le> 0" unfolding le_float_def less_float_def by auto
1896 thus False using assms by auto
1898 thus "0 < Ifloat x" unfolding less_float_def by auto
1899 have "Ifloat (the (lb_ln prec x)) \<le> ln (Ifloat x)" using ub_ln_lb_ln_bounds[OF `0 < x`] ..
1900 thus "Ifloat y \<le> ln (Ifloat x)" unfolding assms[symmetric] by auto
1903 lemma ub_ln: assumes "Some y = ub_ln prec x"
1904 shows "ln (Ifloat x) \<le> Ifloat y" and "0 < Ifloat x"
1908 assume "\<not> 0 < x" hence "x \<le> 0" unfolding le_float_def less_float_def by auto
1909 thus False using assms by auto
1911 thus "0 < Ifloat x" unfolding less_float_def by auto
1912 have "ln (Ifloat x) \<le> Ifloat (the (ub_ln prec x))" using ub_ln_lb_ln_bounds[OF `0 < x`] ..
1913 thus "ln (Ifloat x) \<le> Ifloat y" unfolding assms[symmetric] by auto
1916 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"
1917 proof (rule allI, rule allI, rule allI, rule impI)
1919 assume "(Some l, Some u) = (lb_ln prec lx, ub_ln prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux}"
1920 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
1922 have "ln (Ifloat ux) \<le> Ifloat u" and "0 < Ifloat ux" using ub_ln u by auto
1923 have "Ifloat l \<le> ln (Ifloat lx)" and "0 < Ifloat lx" and "0 < x" using lb_ln[OF l] x by auto
1925 from ln_le_cancel_iff[OF `0 < Ifloat lx` `0 < x`] `Ifloat l \<le> ln (Ifloat lx)`
1926 have "Ifloat l \<le> ln x" using x unfolding atLeastAtMost_iff by auto
1928 from ln_le_cancel_iff[OF `0 < x` `0 < Ifloat ux`] `ln (Ifloat ux) \<le> Ifloat u`
1929 have "ln x \<le> Ifloat u" using x unfolding atLeastAtMost_iff by auto
1930 ultimately show "Ifloat l \<le> ln x \<and> ln x \<le> Ifloat u" ..
1934 section "Implement floatarith"
1936 subsection "Define syntax and semantics"
1939 = Add floatarith floatarith
1941 | Mult floatarith floatarith
1942 | Inverse floatarith
1947 | Max floatarith floatarith
1948 | Min floatarith floatarith
1953 | Power floatarith nat
1957 fun Ifloatarith :: "floatarith \<Rightarrow> real list \<Rightarrow> real"
1959 "Ifloatarith (Add a b) vs = (Ifloatarith a vs) + (Ifloatarith b vs)" |
1960 "Ifloatarith (Minus a) vs = - (Ifloatarith a vs)" |
1961 "Ifloatarith (Mult a b) vs = (Ifloatarith a vs) * (Ifloatarith b vs)" |
1962 "Ifloatarith (Inverse a) vs = inverse (Ifloatarith a vs)" |
1963 "Ifloatarith (Sin a) vs = sin (Ifloatarith a vs)" |
1964 "Ifloatarith (Cos a) vs = cos (Ifloatarith a vs)" |
1965 "Ifloatarith (Arctan a) vs = arctan (Ifloatarith a vs)" |
1966 "Ifloatarith (Min a b) vs = min (Ifloatarith a vs) (Ifloatarith b vs)" |
1967 "Ifloatarith (Max a b) vs = max (Ifloatarith a vs) (Ifloatarith b vs)" |
1968 "Ifloatarith (Abs a) vs = abs (Ifloatarith a vs)" |
1969 "Ifloatarith Pi vs = pi" |
1970 "Ifloatarith (Sqrt a) vs = sqrt (Ifloatarith a vs)" |
1971 "Ifloatarith (Exp a) vs = exp (Ifloatarith a vs)" |
1972 "Ifloatarith (Ln a) vs = ln (Ifloatarith a vs)" |
1973 "Ifloatarith (Power a n) vs = (Ifloatarith a vs)^n" |
1974 "Ifloatarith (Num f) vs = Ifloat f" |
1975 "Ifloatarith (Atom n) vs = vs ! n"
1977 subsection "Implement approximation function"
1979 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
1980 "lift_bin (Some (l1, u1)) (Some (l2, u2)) f = (case (f l1 u1 l2 u2) of (Some l, Some u) \<Rightarrow> Some (l, u)
1981 | t \<Rightarrow> None)" |
1982 "lift_bin a b f = None"
1984 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
1985 "lift_bin' (Some (l1, u1)) (Some (l2, u2)) f = Some (f l1 u1 l2 u2)" |
1986 "lift_bin' a b f = None"
1988 fun lift_un :: "(float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> ((float option) * (float option))) \<Rightarrow> (float * float) option" where
1989 "lift_un (Some (l1, u1)) f = (case (f l1 u1) of (Some l, Some u) \<Rightarrow> Some (l, u)
1990 | t \<Rightarrow> None)" |
1991 "lift_un b f = None"
1993 fun lift_un' :: "(float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> (float * float)) \<Rightarrow> (float * float) option" where
1994 "lift_un' (Some (l1, u1)) f = Some (f l1 u1)" |
1995 "lift_un' b f = None"
1997 fun bounded_by :: "real list \<Rightarrow> (float * float) list \<Rightarrow> bool " where
1998 bounded_by_Cons: "bounded_by (v#vs) ((l, u)#bs) = ((Ifloat l \<le> v \<and> v \<le> Ifloat u) \<and> bounded_by vs bs)" |
1999 bounded_by_Nil: "bounded_by [] [] = True" |
2000 "bounded_by _ _ = False"
2002 lemma bounded_by: assumes "bounded_by vs bs" and "i < length bs"
2003 shows "Ifloat (fst (bs ! i)) \<le> vs ! i \<and> vs ! i \<le> Ifloat (snd (bs ! i))"
2004 using `bounded_by vs bs` and `i < length bs`
2005 proof (induct arbitrary: i rule: bounded_by.induct)
2006 fix v :: real and vs :: "real list" and l u :: float and bs :: "(float * float) list" and i :: nat
2007 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))"
2008 assume bounded: "bounded_by (v # vs) ((l, u) # bs)" and length: "i < length ((l, u) # bs)"
2009 show "Ifloat (fst (((l, u) # bs) ! i)) \<le> (v # vs) ! i \<and> (v # vs) ! i \<le> Ifloat (snd (((l, u) # bs) ! i))"
2012 show ?thesis using bounded unfolding 0 nth_Cons_0 fst_conv snd_conv bounded_by.simps ..
2014 case (Suc i) with length have "i < length bs" by auto
2015 show ?thesis unfolding Suc nth_Cons_Suc bounded_by.simps
2016 using hyp[OF bounded[unfolded bounded_by.simps, THEN conjunct2] `i < length bs`] .
2020 fun approx approx' :: "nat \<Rightarrow> floatarith \<Rightarrow> (float * float) list \<Rightarrow> (float * float) option" where
2021 "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)" |
2022 "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))" |
2023 "approx prec (Minus a) bs = lift_un' (approx' prec a bs) (\<lambda> l u. (-u, -l))" |
2024 "approx prec (Mult a b) bs = lift_bin' (approx' prec a bs) (approx' prec b bs)
2025 (\<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,
2026 float_pprt a2 * float_pprt b2 + float_pprt a1 * float_nprt b2 + float_nprt a2 * float_pprt b1 + float_nprt a1 * float_nprt b1))" |
2027 "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))" |
2028 "approx prec (Sin a) bs = lift_un' (approx' prec a bs) (bnds_sin prec)" |
2029 "approx prec (Cos a) bs = lift_un' (approx' prec a bs) (bnds_cos prec)" |
2030 "approx prec Pi bs = Some (lb_pi prec, ub_pi prec)" |
2031 "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))" |
2032 "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))" |
2033 "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>))" |
2034 "approx prec (Arctan a) bs = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_arctan prec l, ub_arctan prec u))" |
2035 "approx prec (Sqrt a) bs = lift_un (approx' prec a bs) (\<lambda> l u. (lb_sqrt prec l, ub_sqrt prec u))" |
2036 "approx prec (Exp a) bs = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_exp prec l, ub_exp prec u))" |
2037 "approx prec (Ln a) bs = lift_un (approx' prec a bs) (\<lambda> l u. (lb_ln prec l, ub_ln prec u))" |
2038 "approx prec (Power a n) bs = lift_un' (approx' prec a bs) (float_power_bnds n)" |
2039 "approx prec (Num f) bs = Some (f, f)" |
2040 "approx prec (Atom i) bs = (if i < length bs then Some (bs ! i) else None)"
2043 assumes lift_bin'_Some: "Some (l, u) = lift_bin' a b f"
2044 shows "\<exists> l1 u1 l2 u2. Some (l1, u1) = a \<and> Some (l2, u2) = b"
2046 case None hence "None = lift_bin' a b f" unfolding None lift_bin'.simps ..
2047 thus ?thesis using lift_bin'_Some by auto
2052 case None hence "None = lift_bin' a b f" unfolding None lift_bin'.simps ..
2053 thus ?thesis using lift_bin'_Some by auto
2056 obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
2057 obtain lb ub where b': "b' = (lb, ub)" by (cases b', auto)
2058 thus ?thesis unfolding `a = Some a'` `b = Some b'` a' b' by auto
2063 assumes lift_bin'_Some: "Some (l, u) = lift_bin' (g a) (g b) f"
2064 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"
2065 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)"
2068 where Sa: "Some (l1, u1) = g a" and Sb: "Some (l2, u2) = g b" using lift_bin'_ex[OF assms(1)] by auto
2069 have lu: "(l, u) = f l1 u1 l2 u2" using lift_bin'_Some[unfolded Sa[symmetric] Sb[symmetric] lift_bin'.simps] by auto
2070 have "l = fst (f l1 u1 l2 u2)" and "u = snd (f l1 u1 l2 u2)" unfolding lu[symmetric] by auto
2071 thus ?thesis using Pa[OF Sa] Pb[OF Sb] by auto
2074 lemma approx_approx':
2075 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"
2076 and approx': "Some (l, u) = approx' prec a vs"
2077 shows "Ifloat l \<le> Ifloatarith a xs \<and> Ifloatarith a xs \<le> Ifloat u"
2079 obtain l' u' where S: "Some (l', u') = approx prec a vs"
2080 using approx' unfolding approx'.simps by (cases "approx prec a vs", auto)
2081 have l': "l = round_down prec l'" and u': "u = round_up prec u'"
2082 using approx' unfolding approx'.simps S[symmetric] by auto
2083 show ?thesis unfolding l' u'
2084 using order_trans[OF Pa[OF S, THEN conjunct2] round_up[of u']]
2085 using order_trans[OF round_down[of _ l'] Pa[OF S, THEN conjunct1]] by auto
2089 assumes lift_bin'_Some: "Some (l, u) = lift_bin' (approx' prec a bs) (approx' prec b bs) f"
2090 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")
2091 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"
2092 shows "\<exists> l1 u1 l2 u2. (Ifloat l1 \<le> Ifloatarith a xs \<and> Ifloatarith a xs \<le> Ifloat u1) \<and>
2093 (Ifloat l2 \<le> Ifloatarith b xs \<and> Ifloatarith b xs \<le> Ifloat u2) \<and>
2094 l = fst (f l1 u1 l2 u2) \<and> u = snd (f l1 u1 l2 u2)"
2096 { fix l u assume "Some (l, u) = approx' prec a bs"
2097 with approx_approx'[of prec a bs, OF _ this] Pa
2098 have "Ifloat l \<le> Ifloatarith a xs \<and> Ifloatarith a xs \<le> Ifloat u" by auto } note Pa = this
2099 { fix l u assume "Some (l, u) = approx' prec b bs"
2100 with approx_approx'[of prec b bs, OF _ this] Pb
2101 have "Ifloat l \<le> Ifloatarith b xs \<and> Ifloatarith b xs \<le> Ifloat u" by auto } note Pb = this
2103 from lift_bin'_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_bin'_Some, OF Pa Pb]
2104 show ?thesis by auto
2108 assumes lift_un'_Some: "Some (l, u) = lift_un' a f"
2109 shows "\<exists> l u. Some (l, u) = a"
2111 case None hence "None = lift_un' a f" unfolding None lift_un'.simps ..
2112 thus ?thesis using lift_un'_Some by auto
2115 obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
2116 thus ?thesis unfolding `a = Some a'` a' by auto
2120 assumes lift_un'_Some: "Some (l, u) = lift_un' (g a) f"
2121 and Pa: "\<And>l u. Some (l, u) = g a \<Longrightarrow> P l u a"
2122 shows "\<exists> l1 u1. P l1 u1 a \<and> l = fst (f l1 u1) \<and> u = snd (f l1 u1)"
2124 obtain l1 u1 where Sa: "Some (l1, u1) = g a" using lift_un'_ex[OF assms(1)] by auto
2125 have lu: "(l, u) = f l1 u1" using lift_un'_Some[unfolded Sa[symmetric] lift_un'.simps] by auto
2126 have "l = fst (f l1 u1)" and "u = snd (f l1 u1)" unfolding lu[symmetric] by auto
2127 thus ?thesis using Pa[OF Sa] by auto
2131 assumes lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f"
2132 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")
2133 shows "\<exists> l1 u1. (Ifloat l1 \<le> Ifloatarith a xs \<and> Ifloatarith a xs \<le> Ifloat u1) \<and>
2134 l = fst (f l1 u1) \<and> u = snd (f l1 u1)"
2136 { fix l u assume "Some (l, u) = approx' prec a bs"
2137 with approx_approx'[of prec a bs, OF _ this] Pa
2138 have "Ifloat l \<le> Ifloatarith a xs \<and> Ifloatarith a xs \<le> Ifloat u" by auto } note Pa = this
2139 from lift_un'_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_un'_Some, OF Pa]
2140 show ?thesis by auto
2143 lemma lift_un'_bnds:
2144 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"
2145 and lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f"
2146 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"
2147 shows "Ifloat l \<le> f' (Ifloatarith a xs) \<and> f' (Ifloatarith a xs) \<le> Ifloat u"
2149 from lift_un'[OF lift_un'_Some Pa]
2150 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
2151 hence "(l, u) = f l1 u1" and "Ifloatarith a xs \<in> {Ifloat l1 .. Ifloat u1}" by auto
2152 thus ?thesis using bnds by auto
2156 assumes lift_un_Some: "Some (l, u) = lift_un a f"
2157 shows "\<exists> l u. Some (l, u) = a"
2159 case None hence "None = lift_un a f" unfolding None lift_un.simps ..
2160 thus ?thesis using lift_un_Some by auto
2163 obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
2164 thus ?thesis unfolding `a = Some a'` a' by auto
2168 assumes lift_un_Some: "Some (l, u) = lift_un (g a) f"
2169 and Pa: "\<And>l u. Some (l, u) = g a \<Longrightarrow> P l u a"
2170 shows "\<exists> l1 u1. P l1 u1 a \<and> Some l = fst (f l1 u1) \<and> Some u = snd (f l1 u1)"
2172 obtain l1 u1 where Sa: "Some (l1, u1) = g a" using lift_un_ex[OF assms(1)] by auto
2173 have "fst (f l1 u1) \<noteq> None \<and> snd (f l1 u1) \<noteq> None"
2175 assume "\<not> (fst (f l1 u1) \<noteq> None \<and> snd (f l1 u1) \<noteq> None)"
2176 hence or: "fst (f l1 u1) = None \<or> snd (f l1 u1) = None" by auto
2177 hence "lift_un (g a) f = None"
2178 proof (cases "fst (f l1 u1) = None")
2180 then obtain b where b: "f l1 u1 = (None, b)" by (cases "f l1 u1", auto)
2181 thus ?thesis unfolding Sa[symmetric] lift_un.simps b by auto
2183 case False hence "snd (f l1 u1) = None" using or by auto
2184 with False obtain b where b: "f l1 u1 = (Some b, None)" by (cases "f l1 u1", auto)
2185 thus ?thesis unfolding Sa[symmetric] lift_un.simps b by auto
2187 thus False using lift_un_Some by auto
2189 then obtain a' b' where f: "f l1 u1 = (Some a', Some b')" by (cases "f l1 u1", auto)
2190 from lift_un_Some[unfolded Sa[symmetric] lift_un.simps f]
2191 have "Some l = fst (f l1 u1)" and "Some u = snd (f l1 u1)" unfolding f by auto
2192 thus ?thesis unfolding Sa[symmetric] lift_un.simps using Pa[OF Sa] by auto
2196 assumes lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f"
2197 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")
2198 shows "\<exists> l1 u1. (Ifloat l1 \<le> Ifloatarith a xs \<and> Ifloatarith a xs \<le> Ifloat u1) \<and>
2199 Some l = fst (f l1 u1) \<and> Some u = snd (f l1 u1)"
2201 { fix l u assume "Some (l, u) = approx' prec a bs"
2202 with approx_approx'[of prec a bs, OF _ this] Pa
2203 have "Ifloat l \<le> Ifloatarith a xs \<and> Ifloatarith a xs \<le> Ifloat u" by auto } note Pa = this
2204 from lift_un_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_un_Some, OF Pa]
2205 show ?thesis by auto
2209 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"
2210 and lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f"
2211 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"
2212 shows "Ifloat l \<le> f' (Ifloatarith a xs) \<and> f' (Ifloatarith a xs) \<le> Ifloat u"
2214 from lift_un[OF lift_un_Some Pa]
2215 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
2216 hence "(Some l, Some u) = f l1 u1" and "Ifloatarith a xs \<in> {Ifloat l1 .. Ifloat u1}" by auto
2217 thus ?thesis using bnds by auto
2221 assumes "bounded_by xs vs"
2222 and "Some (l, u) = approx prec arith vs" (is "_ = ?g arith")
2223 shows "Ifloat l \<le> Ifloatarith arith xs \<and> Ifloatarith arith xs \<le> Ifloat u" (is "?P l u arith")
2224 using `Some (l, u) = approx prec arith vs`
2225 proof (induct arith arbitrary: l u x)
2227 from lift_bin'[OF Add.prems[unfolded approx.simps]] Add.hyps
2228 obtain l1 u1 l2 u2 where "l = l1 + l2" and "u = u1 + u2"
2229 "Ifloat l1 \<le> Ifloatarith a xs" and "Ifloatarith a xs \<le> Ifloat u1"
2230 "Ifloat l2 \<le> Ifloatarith b xs" and "Ifloatarith b xs \<le> Ifloat u2" unfolding fst_conv snd_conv by blast
2231 thus ?case unfolding Ifloatarith.simps by auto
2234 from lift_un'[OF Minus.prems[unfolded approx.simps]] Minus.hyps
2235 obtain l1 u1 where "l = -u1" and "u = -l1"
2236 "Ifloat l1 \<le> Ifloatarith a xs" and "Ifloatarith a xs \<le> Ifloat u1" unfolding fst_conv snd_conv by blast
2237 thus ?case unfolding Ifloatarith.simps using Ifloat_minus by auto
2240 from lift_bin'[OF Mult.prems[unfolded approx.simps]] Mult.hyps
2242 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"
2243 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"
2244 and "Ifloat l1 \<le> Ifloatarith a xs" and "Ifloatarith a xs \<le> Ifloat u1"
2245 and "Ifloat l2 \<le> Ifloatarith b xs" and "Ifloatarith b xs \<le> Ifloat u2" unfolding fst_conv snd_conv by blast
2246 thus ?case unfolding Ifloatarith.simps l u Ifloat_add Ifloat_mult Ifloat_nprt Ifloat_pprt
2247 using mult_le_prts mult_ge_prts by auto
2250 from lift_un[OF Inverse.prems[unfolded approx.simps], unfolded if_distrib[of fst] if_distrib[of snd] fst_conv snd_conv] Inverse.hyps
2251 obtain l1 u1 where l': "Some l = (if 0 < l1 \<or> u1 < 0 then Some (float_divl prec 1 u1) else None)"
2252 and u': "Some u = (if 0 < l1 \<or> u1 < 0 then Some (float_divr prec 1 l1) else None)"
2253 and l1: "Ifloat l1 \<le> Ifloatarith a xs" and u1: "Ifloatarith a xs \<le> Ifloat u1" by blast
2254 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
2255 moreover have l1_le_u1: "Ifloat l1 \<le> Ifloat u1" using l1 u1 by auto
2256 ultimately have "Ifloat l1 \<noteq> 0" and "Ifloat u1 \<noteq> 0" unfolding less_float_def by auto
2258 have inv: "inverse (Ifloat u1) \<le> inverse (Ifloatarith a xs)
2259 \<and> inverse (Ifloatarith a xs) \<le> inverse (Ifloat l1)"
2260 proof (cases "0 < l1")
2261 case True hence "0 < Ifloat u1" and "0 < Ifloat l1" "0 < Ifloatarith a xs"
2262 unfolding less_float_def using l1_le_u1 l1 by auto
2264 unfolding inverse_le_iff_le[OF `0 < Ifloat u1` `0 < Ifloatarith a xs`]
2265 inverse_le_iff_le[OF `0 < Ifloatarith a xs` `0 < Ifloat l1`]
2268 case False hence "u1 < 0" using either by blast
2269 hence "Ifloat u1 < 0" and "Ifloat l1 < 0" "Ifloatarith a xs < 0"
2270 unfolding less_float_def using l1_le_u1 u1 by auto
2272 unfolding inverse_le_iff_le_neg[OF `Ifloat u1 < 0` `Ifloatarith a xs < 0`]
2273 inverse_le_iff_le_neg[OF `Ifloatarith a xs < 0` `Ifloat l1 < 0`]
2277 from l' have "l = float_divl prec 1 u1" by (cases "0 < l1 \<or> u1 < 0", auto)
2278 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
2279 also have "\<dots> \<le> inverse (Ifloatarith a xs)" using inv by auto
2280 finally have "Ifloat l \<le> inverse (Ifloatarith a xs)" .
2282 from u' have "u = float_divr prec 1 l1" by (cases "0 < l1 \<or> u1 < 0", auto)
2283 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
2284 hence "inverse (Ifloatarith a xs) \<le> Ifloat u" by (rule order_trans[OF inv[THEN conjunct2]])
2285 ultimately show ?case unfolding Ifloatarith.simps using l1 u1 by auto
2288 from lift_un'[OF Abs.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Abs.hyps
2289 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>"
2290 and l1: "Ifloat l1 \<le> Ifloatarith x xs" and u1: "Ifloatarith x xs \<le> Ifloat u1" by blast
2291 thus ?case unfolding l' u' by (cases "l1 < 0 \<and> 0 < u1", auto simp add: Ifloat_min Ifloat_max Ifloat_abs less_float_def)
2294 from lift_bin'[OF Min.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Min.hyps
2295 obtain l1 u1 l2 u2 where l': "l = min l1 l2" and u': "u = min u1 u2"
2296 and l1: "Ifloat l1 \<le> Ifloatarith a xs" and u1: "Ifloatarith a xs \<le> Ifloat u1"
2297 and l1: "Ifloat l2 \<le> Ifloatarith b xs" and u1: "Ifloatarith b xs \<le> Ifloat u2" by blast
2298 thus ?case unfolding l' u' by (auto simp add: Ifloat_min)
2301 from lift_bin'[OF Max.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Max.hyps
2302 obtain l1 u1 l2 u2 where l': "l = max l1 l2" and u': "u = max u1 u2"
2303 and l1: "Ifloat l1 \<le> Ifloatarith a xs" and u1: "Ifloatarith a xs \<le> Ifloat u1"
2304 and l1: "Ifloat l2 \<le> Ifloatarith b xs" and u1: "Ifloatarith b xs \<le> Ifloat u2" by blast
2305 thus ?case unfolding l' u' by (auto simp add: Ifloat_max)
2306 next case (Sin a) with lift_un'_bnds[OF bnds_sin] show ?case by auto
2307 next case (Cos a) with lift_un'_bnds[OF bnds_cos] show ?case by auto
2308 next case (Arctan a) with lift_un'_bnds[OF bnds_arctan] show ?case by auto
2309 next case Pi with pi_boundaries show ?case by auto
2310 next case (Sqrt a) with lift_un_bnds[OF bnds_sqrt] show ?case by auto
2311 next case (Exp a) with lift_un'_bnds[OF bnds_exp] show ?case by auto
2312 next case (Ln a) with lift_un_bnds[OF bnds_ln] show ?case by auto
2313 next case (Power a n) with lift_un'_bnds[OF bnds_power] show ?case by auto
2314 next case (Num f) thus ?case by auto
2318 proof (cases "n < length vs")
2320 with Atom have "vs ! n = (l, u)" by auto
2321 thus ?thesis using bounded_by[OF assms(1) True] by auto
2323 case False thus ?thesis using Atom by auto
2327 datatype ApproxEq = Less floatarith floatarith
2328 | LessEqual floatarith floatarith
2330 fun uneq :: "ApproxEq \<Rightarrow> real list \<Rightarrow> bool" where
2331 "uneq (Less a b) vs = (Ifloatarith a vs < Ifloatarith b vs)" |
2332 "uneq (LessEqual a b) vs = (Ifloatarith a vs \<le> Ifloatarith b vs)"
2334 fun uneq' :: "nat \<Rightarrow> ApproxEq \<Rightarrow> (float * float) list \<Rightarrow> bool" where
2335 "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)" |
2336 "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)"
2338 lemma uneq_approx: fixes m :: nat assumes "bounded_by vs bs" and "uneq' prec eq bs"
2343 proof (cases "\<exists> u l u' l'. approx prec a bs = Some (l, u) \<and>
2344 approx prec b bs = Some (l', u')")
2346 then obtain l u l' u' where a_approx: "approx prec a bs = Some (l, u)"
2347 and b_approx: "approx prec b bs = Some (l', u') " by auto
2348 with `uneq' prec eq bs` have "Ifloat u < Ifloat l'"
2349 unfolding Less uneq'.simps less_float_def by auto
2350 moreover from a_approx[symmetric] and b_approx[symmetric] and `bounded_by vs bs`
2351 have "Ifloatarith a vs \<le> Ifloat u" and "Ifloat l' \<le> Ifloatarith b vs"
2352 using approx by auto
2353 ultimately show ?thesis unfolding uneq.simps Less by auto
2356 hence "approx prec a bs = None \<or> approx prec b bs = None"
2357 unfolding not_Some_eq[symmetric] by auto
2358 hence "\<not> uneq' prec eq bs" unfolding Less uneq'.simps
2359 by (cases "approx prec a bs = None", auto)
2360 thus ?thesis using assms by auto
2363 case (LessEqual a b)
2365 proof (cases "\<exists> u l u' l'. approx prec a bs = Some (l, u) \<and>
2366 approx prec b bs = Some (l', u')")
2368 then obtain l u l' u' where a_approx: "approx prec a bs = Some (l, u)"
2369 and b_approx: "approx prec b bs = Some (l', u') " by auto
2370 with `uneq' prec eq bs` have "Ifloat u \<le> Ifloat l'"
2371 unfolding LessEqual uneq'.simps le_float_def by auto
2372 moreover from a_approx[symmetric] and b_approx[symmetric] and `bounded_by vs bs`
2373 have "Ifloatarith a vs \<le> Ifloat u" and "Ifloat l' \<le> Ifloatarith b vs"
2374 using approx by auto
2375 ultimately show ?thesis unfolding uneq.simps LessEqual by auto
2378 hence "approx prec a bs = None \<or> approx prec b bs = None"
2379 unfolding not_Some_eq[symmetric] by auto
2380 hence "\<not> uneq' prec eq bs" unfolding LessEqual uneq'.simps
2381 by (cases "approx prec a bs = None", auto)
2382 thus ?thesis using assms by auto
2386 lemma Ifloatarith_divide: "Ifloatarith (Mult a (Inverse b)) vs = (Ifloatarith a vs) / (Ifloatarith b vs)"
2387 unfolding real_divide_def Ifloatarith.simps ..
2389 lemma Ifloatarith_diff: "Ifloatarith (Add a (Minus b)) vs = (Ifloatarith a vs) - (Ifloatarith b vs)"
2390 unfolding real_diff_def Ifloatarith.simps ..
2392 lemma Ifloatarith_tan: "Ifloatarith (Mult (Sin a) (Inverse (Cos a))) vs = tan (Ifloatarith a vs)"
2393 unfolding tan_def Ifloatarith.simps real_divide_def ..
2395 lemma Ifloatarith_powr: "Ifloatarith (Exp (Mult b (Ln a))) vs = (Ifloatarith a vs) powr (Ifloatarith b vs)"
2396 unfolding powr_def Ifloatarith.simps ..
2398 lemma Ifloatarith_log: "Ifloatarith ((Mult (Ln x) (Inverse (Ln b)))) vs = log (Ifloatarith b vs) (Ifloatarith x vs)"
2399 unfolding log_def Ifloatarith.simps real_divide_def ..
2401 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
2403 subsection {* Implement proof method \texttt{approximation} *}
2405 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)
2406 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)
2407 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)"
2408 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)"
2409 and "Ifloat (Float (number_of A) (int B)) = (number_of A) * 2^B"
2410 and "Ifloat (Float 1 (int B)) = 2^B"
2411 and "Ifloat (Float (number_of A) (- int B)) = (number_of A) / 2^B"
2412 and "Ifloat (Float 1 (- int B)) = 1 / 2^B"
2413 by (auto simp add: Ifloat.simps pow2_def real_divide_def)
2415 lemmas bounded_by_equations = bounded_by_Cons bounded_by_Nil float_power bounded_divl bounded_divr bounded_num HOL.simp_thms
2416 lemmas uneq_equations = uneq.simps Ifloatarith.simps Ifloatarith_num Ifloatarith_divide Ifloatarith_diff Ifloatarith_tan Ifloatarith_powr Ifloatarith_log
2419 val uneq_equations = PureThy.get_thms @{theory} "uneq_equations";
2420 val bounded_by_equations = PureThy.get_thms @{theory} "bounded_by_equations";
2421 val bounded_by_simpset = (HOL_basic_ss addsimps bounded_by_equations)
2423 fun reify_uneq ctxt i = (fn st =>
2425 val to = HOLogic.dest_Trueprop (Logic.strip_imp_concl (List.nth (prems_of st, i - 1)))
2426 in (Reflection.genreify_tac ctxt uneq_equations (SOME to) i) st
2429 fun rule_uneq ctxt prec i thm = let
2430 fun conv_num typ = HOLogic.dest_number #> snd #> HOLogic.mk_number typ
2431 val to_natc = conv_num @{typ "nat"} #> Thm.cterm_of (ProofContext.theory_of ctxt)
2432 val to_nat = conv_num @{typ "nat"}
2433 val to_int = conv_num @{typ "int"}
2434 fun int_to_float A = @{term "Float"} $ to_int A $ @{term "0 :: int"}
2436 val prec' = to_nat prec
2438 fun bot_float (Const (@{const_name "times"}, _) $ mantisse $ (Const (@{const_name "pow2"}, _) $ exp))
2439 = @{term "Float"} $ to_int mantisse $ to_int exp
2440 | bot_float (Const (@{const_name "divide"}, _) $ mantisse $ (@{term "power 2 :: nat \<Rightarrow> real"} $ exp))
2441 = @{term "Float"} $ to_int mantisse $ (@{term "uminus :: int \<Rightarrow> int"} $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp))
2442 | bot_float (Const (@{const_name "times"}, _) $ mantisse $ (@{term "power 2 :: nat \<Rightarrow> real"} $ exp))
2443 = @{term "Float"} $ to_int mantisse $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp)
2444 | bot_float (Const (@{const_name "divide"}, _) $ A $ (@{term "power 10 :: nat \<Rightarrow> real"} $ exp))
2445 = @{term "float_divl"} $ prec' $ (int_to_float A) $ (@{term "power (Float 5 1)"} $ to_nat exp)
2446 | bot_float (Const (@{const_name "divide"}, _) $ A $ B)
2447 = @{term "float_divl"} $ prec' $ int_to_float A $ int_to_float B
2448 | bot_float (@{term "power 2 :: nat \<Rightarrow> real"} $ exp)
2449 = @{term "Float 1"} $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp)
2450 | bot_float A = int_to_float A
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 $ (@{term "power 2 :: nat \<Rightarrow> real"} $ exp))
2455 = @{term "Float"} $ to_int mantisse $ (@{term "uminus :: int \<Rightarrow> int"} $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp))
2456 | top_float (Const (@{const_name "times"}, _) $ mantisse $ (@{term "power 2 :: nat \<Rightarrow> real"} $ exp))
2457 = @{term "Float"} $ to_int mantisse $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp)
2458 | top_float (Const (@{const_name "divide"}, _) $ A $ (@{term "power 10 :: nat \<Rightarrow> real"} $ exp))
2459 = @{term "float_divr"} $ prec' $ (int_to_float A) $ (@{term "power (Float 5 1)"} $ to_nat exp)
2460 | top_float (Const (@{const_name "divide"}, _) $ A $ B)
2461 = @{term "float_divr"} $ prec' $ int_to_float A $ int_to_float B
2462 | top_float (@{term "power 2 :: nat \<Rightarrow> real"} $ exp)
2463 = @{term "Float 1"} $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp)
2464 | top_float A = int_to_float A
2466 val goal' : term = List.nth (prems_of thm, i - 1)
2468 fun lift_bnd (t as (Const (@{const_name "op &"}, _) $
2469 (Const (@{const_name "less_eq"}, _) $
2470 bottom $ (Free (name, _))) $
2471 (Const (@{const_name "less_eq"}, _) $ _ $ top)))
2472 = ((name, HOLogic.mk_prod (bot_float bottom, top_float top))
2473 handle TERM (txt, ts) => raise TERM ("Invalid bound number format: " ^
2474 (Syntax.string_of_term ctxt t), [t]))
2475 | lift_bnd t = raise TERM ("Premisse needs format '<num> <= <var> & <var> <= <num>', but found " ^
2476 (Syntax.string_of_term ctxt t), [t])
2477 val bound_eqs = map (HOLogic.dest_Trueprop #> lift_bnd) (Logic.strip_imp_prems goal')
2479 fun lift_var (Free (varname, _)) = (case AList.lookup (op =) bound_eqs varname of
2481 | NONE => raise TERM ("No bound equations found for " ^ varname, []))
2482 | lift_var t = raise TERM ("Can not convert expression " ^
2483 (Syntax.string_of_term ctxt t), [t])
2485 val _ $ vs = HOLogic.dest_Trueprop (Logic.strip_imp_concl goal')
2487 val bs = (HOLogic.dest_list #> map lift_var #> HOLogic.mk_list @{typ "float * float"}) vs
2488 val map = [(@{cpat "?prec::nat"}, to_natc prec),
2489 (@{cpat "?bs::(float * float) list"}, Thm.cterm_of (ProofContext.theory_of ctxt) bs)]
2490 in rtac (Thm.instantiate ([], map) @{thm "uneq_approx"}) i thm end
2492 val eval_tac = CSUBGOAL (fn (ct, i) => rtac (eval_oracle ct) i)
2494 fun gen_eval_tac conv ctxt = CONVERSION (Conv.params_conv (~1) (K (Conv.concl_conv (~1) conv)) ctxt)
2499 method_setup approximation = {*
2501 (fn prec => fn ctxt =>
2502 SIMPLE_METHOD' (fn i =>
2503 (DETERM (reify_uneq ctxt i)
2504 THEN rule_uneq ctxt prec i
2505 THEN Simplifier.asm_full_simp_tac bounded_by_simpset i
2506 THEN (TRY (filter_prems_tac (fn t => false) i))
2507 THEN (gen_eval_tac eval_oracle ctxt) i)))
2508 *} "real number approximation"