|
1 (* |
|
2 Author: Johannes Hoelzl, TU Muenchen |
|
3 Coercions removed by Dmitriy Traytel |
|
4 |
|
5 This file contains only general material about computing lower/upper bounds |
|
6 on real functions. Approximation.thy contains the actual approximation algorithm |
|
7 and the approximation oracle. This is in order to make a clear separation between |
|
8 "morally immaculate" material about upper/lower bounds and the trusted oracle/reflection. |
|
9 *) |
|
10 |
|
11 theory Approximation_Bounds |
|
12 imports |
|
13 Complex_Main |
|
14 "~~/src/HOL/Library/Float" |
|
15 Dense_Linear_Order |
|
16 begin |
|
17 |
|
18 declare powr_neg_one [simp] |
|
19 declare powr_neg_numeral [simp] |
|
20 |
|
21 section "Horner Scheme" |
|
22 |
|
23 subsection \<open>Define auxiliary helper \<open>horner\<close> function\<close> |
|
24 |
|
25 primrec horner :: "(nat \<Rightarrow> nat) \<Rightarrow> (nat \<Rightarrow> nat \<Rightarrow> nat) \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> real \<Rightarrow> real" where |
|
26 "horner F G 0 i k x = 0" | |
|
27 "horner F G (Suc n) i k x = 1 / k - x * horner F G n (F i) (G i k) x" |
|
28 |
|
29 lemma horner_schema': |
|
30 fixes x :: real and a :: "nat \<Rightarrow> real" |
|
31 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)" |
|
32 proof - |
|
33 have shift_pow: "\<And>i. - (x * ((-1)^i * a (Suc i) * x ^ i)) = (-1)^(Suc i) * a (Suc i) * x ^ (Suc i)" |
|
34 by auto |
|
35 show ?thesis |
|
36 unfolding sum_distrib_left shift_pow uminus_add_conv_diff [symmetric] sum_negf[symmetric] |
|
37 sum_head_upt_Suc[OF zero_less_Suc] |
|
38 sum.reindex[OF inj_Suc, unfolded comp_def, symmetric, of "\<lambda> n. (-1)^n *a n * x^n"] by auto |
|
39 qed |
|
40 |
|
41 lemma horner_schema: |
|
42 fixes f :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat" and F :: "nat \<Rightarrow> nat" |
|
43 assumes f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)" |
|
44 shows "horner F G n ((F ^^ j') s) (f j') x = (\<Sum> j = 0..< n. (- 1) ^ j * (1 / (f (j' + j))) * x ^ j)" |
|
45 proof (induct n arbitrary: j') |
|
46 case 0 |
|
47 then show ?case by auto |
|
48 next |
|
49 case (Suc n) |
|
50 show ?case unfolding horner.simps Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc] |
|
51 using horner_schema'[of "\<lambda> j. 1 / (f (j' + j))"] by auto |
|
52 qed |
|
53 |
|
54 lemma horner_bounds': |
|
55 fixes lb :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" and ub :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" |
|
56 assumes "0 \<le> real_of_float x" and f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)" |
|
57 and lb_0: "\<And> i k x. lb 0 i k x = 0" |
|
58 and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = float_plus_down prec |
|
59 (lapprox_rat prec 1 k) |
|
60 (- float_round_up prec (x * (ub n (F i) (G i k) x)))" |
|
61 and ub_0: "\<And> i k x. ub 0 i k x = 0" |
|
62 and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = float_plus_up prec |
|
63 (rapprox_rat prec 1 k) |
|
64 (- float_round_down prec (x * (lb n (F i) (G i k) x)))" |
|
65 shows "(lb n ((F ^^ j') s) (f j') x) \<le> horner F G n ((F ^^ j') s) (f j') x \<and> |
|
66 horner F G n ((F ^^ j') s) (f j') x \<le> (ub n ((F ^^ j') s) (f j') x)" |
|
67 (is "?lb n j' \<le> ?horner n j' \<and> ?horner n j' \<le> ?ub n j'") |
|
68 proof (induct n arbitrary: j') |
|
69 case 0 |
|
70 thus ?case unfolding lb_0 ub_0 horner.simps by auto |
|
71 next |
|
72 case (Suc n) |
|
73 thus ?case using lapprox_rat[of prec 1 "f j'"] using rapprox_rat[of 1 "f j'" prec] |
|
74 Suc[where j'="Suc j'"] \<open>0 \<le> real_of_float x\<close> |
|
75 by (auto intro!: add_mono mult_left_mono float_round_down_le float_round_up_le |
|
76 order_trans[OF add_mono[OF _ float_plus_down_le]] |
|
77 order_trans[OF _ add_mono[OF _ float_plus_up_le]] |
|
78 simp add: lb_Suc ub_Suc field_simps f_Suc) |
|
79 qed |
|
80 |
|
81 subsection "Theorems for floating point functions implementing the horner scheme" |
|
82 |
|
83 text \<open> |
|
84 |
|
85 Here @{term_type "f :: nat \<Rightarrow> nat"} is the sequence defining the Taylor series, the coefficients are |
|
86 all alternating and reciprocs. We use @{term G} and @{term F} to describe the computation of @{term f}. |
|
87 |
|
88 \<close> |
|
89 |
|
90 lemma horner_bounds: |
|
91 fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat" |
|
92 assumes "0 \<le> real_of_float x" and f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)" |
|
93 and lb_0: "\<And> i k x. lb 0 i k x = 0" |
|
94 and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = float_plus_down prec |
|
95 (lapprox_rat prec 1 k) |
|
96 (- float_round_up prec (x * (ub n (F i) (G i k) x)))" |
|
97 and ub_0: "\<And> i k x. ub 0 i k x = 0" |
|
98 and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = float_plus_up prec |
|
99 (rapprox_rat prec 1 k) |
|
100 (- float_round_down prec (x * (lb n (F i) (G i k) x)))" |
|
101 shows "(lb n ((F ^^ j') s) (f j') x) \<le> (\<Sum>j=0..<n. (- 1) ^ j * (1 / (f (j' + j))) * (x ^ j))" |
|
102 (is "?lb") |
|
103 and "(\<Sum>j=0..<n. (- 1) ^ j * (1 / (f (j' + j))) * (x ^ j)) \<le> (ub n ((F ^^ j') s) (f j') x)" |
|
104 (is "?ub") |
|
105 proof - |
|
106 have "?lb \<and> ?ub" |
|
107 using horner_bounds'[where lb=lb, OF \<open>0 \<le> real_of_float x\<close> f_Suc lb_0 lb_Suc ub_0 ub_Suc] |
|
108 unfolding horner_schema[where f=f, OF f_Suc] by simp |
|
109 thus "?lb" and "?ub" by auto |
|
110 qed |
|
111 |
|
112 lemma horner_bounds_nonpos: |
|
113 fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat" |
|
114 assumes "real_of_float x \<le> 0" and f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)" |
|
115 and lb_0: "\<And> i k x. lb 0 i k x = 0" |
|
116 and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = float_plus_down prec |
|
117 (lapprox_rat prec 1 k) |
|
118 (float_round_down prec (x * (ub n (F i) (G i k) x)))" |
|
119 and ub_0: "\<And> i k x. ub 0 i k x = 0" |
|
120 and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = float_plus_up prec |
|
121 (rapprox_rat prec 1 k) |
|
122 (float_round_up prec (x * (lb n (F i) (G i k) x)))" |
|
123 shows "(lb n ((F ^^ j') s) (f j') x) \<le> (\<Sum>j=0..<n. (1 / (f (j' + j))) * real_of_float x ^ j)" (is "?lb") |
|
124 and "(\<Sum>j=0..<n. (1 / (f (j' + j))) * real_of_float x ^ j) \<le> (ub n ((F ^^ j') s) (f j') x)" (is "?ub") |
|
125 proof - |
|
126 have diff_mult_minus: "x - y * z = x + - y * z" for x y z :: float by simp |
|
127 have sum_eq: "(\<Sum>j=0..<n. (1 / (f (j' + j))) * real_of_float x ^ j) = |
|
128 (\<Sum>j = 0..<n. (- 1) ^ j * (1 / (f (j' + j))) * real_of_float (- x) ^ j)" |
|
129 by (auto simp add: field_simps power_mult_distrib[symmetric]) |
|
130 have "0 \<le> real_of_float (-x)" using assms by auto |
|
131 from horner_bounds[where G=G and F=F and f=f and s=s and prec=prec |
|
132 and lb="\<lambda> n i k x. lb n i k (-x)" and ub="\<lambda> n i k x. ub n i k (-x)", |
|
133 unfolded lb_Suc ub_Suc diff_mult_minus, |
|
134 OF this f_Suc lb_0 _ ub_0 _] |
|
135 show "?lb" and "?ub" unfolding minus_minus sum_eq |
|
136 by (auto simp: minus_float_round_up_eq minus_float_round_down_eq) |
|
137 qed |
|
138 |
|
139 |
|
140 subsection \<open>Selectors for next even or odd number\<close> |
|
141 |
|
142 text \<open> |
|
143 The horner scheme computes alternating series. To get the upper and lower bounds we need to |
|
144 guarantee to access a even or odd member. To do this we use @{term get_odd} and @{term get_even}. |
|
145 \<close> |
|
146 |
|
147 definition get_odd :: "nat \<Rightarrow> nat" where |
|
148 "get_odd n = (if odd n then n else (Suc n))" |
|
149 |
|
150 definition get_even :: "nat \<Rightarrow> nat" where |
|
151 "get_even n = (if even n then n else (Suc n))" |
|
152 |
|
153 lemma get_odd[simp]: "odd (get_odd n)" |
|
154 unfolding get_odd_def by (cases "odd n") auto |
|
155 |
|
156 lemma get_even[simp]: "even (get_even n)" |
|
157 unfolding get_even_def by (cases "even n") auto |
|
158 |
|
159 lemma get_odd_ex: "\<exists> k. Suc k = get_odd n \<and> odd (Suc k)" |
|
160 by (auto simp: get_odd_def odd_pos intro!: exI[of _ "n - 1"]) |
|
161 |
|
162 lemma get_even_double: "\<exists>i. get_even n = 2 * i" |
|
163 using get_even by (blast elim: evenE) |
|
164 |
|
165 lemma get_odd_double: "\<exists>i. get_odd n = 2 * i + 1" |
|
166 using get_odd by (blast elim: oddE) |
|
167 |
|
168 |
|
169 section "Power function" |
|
170 |
|
171 definition float_power_bnds :: "nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where |
|
172 "float_power_bnds prec n l u = |
|
173 (if 0 < l then (power_down_fl prec l n, power_up_fl prec u n) |
|
174 else if odd n then |
|
175 (- power_up_fl prec \<bar>l\<bar> n, |
|
176 if u < 0 then - power_down_fl prec \<bar>u\<bar> n else power_up_fl prec u n) |
|
177 else if u < 0 then (power_down_fl prec \<bar>u\<bar> n, power_up_fl prec \<bar>l\<bar> n) |
|
178 else (0, power_up_fl prec (max \<bar>l\<bar> \<bar>u\<bar>) n))" |
|
179 |
|
180 lemma le_minus_power_downI: "0 \<le> x \<Longrightarrow> x ^ n \<le> - a \<Longrightarrow> a \<le> - power_down prec x n" |
|
181 by (subst le_minus_iff) (auto intro: power_down_le power_mono_odd) |
|
182 |
|
183 lemma float_power_bnds: |
|
184 "(l1, u1) = float_power_bnds prec n l u \<Longrightarrow> x \<in> {l .. u} \<Longrightarrow> (x::real) ^ n \<in> {l1..u1}" |
|
185 by (auto |
|
186 simp: float_power_bnds_def max_def real_power_up_fl real_power_down_fl minus_le_iff |
|
187 split: if_split_asm |
|
188 intro!: power_up_le power_down_le le_minus_power_downI |
|
189 intro: power_mono_odd power_mono power_mono_even zero_le_even_power) |
|
190 |
|
191 lemma bnds_power: |
|
192 "\<forall>(x::real) l u. (l1, u1) = float_power_bnds prec n l u \<and> x \<in> {l .. u} \<longrightarrow> |
|
193 l1 \<le> x ^ n \<and> x ^ n \<le> u1" |
|
194 using float_power_bnds by auto |
|
195 |
|
196 section \<open>Approximation utility functions\<close> |
|
197 |
|
198 definition bnds_mult :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float \<times> float" where |
|
199 "bnds_mult prec a1 a2 b1 b2 = |
|
200 (float_plus_down prec (nprt a1 * pprt b2) |
|
201 (float_plus_down prec (nprt a2 * nprt b2) |
|
202 (float_plus_down prec (pprt a1 * pprt b1) (pprt a2 * nprt b1))), |
|
203 float_plus_up prec (pprt a2 * pprt b2) |
|
204 (float_plus_up prec (pprt a1 * nprt b2) |
|
205 (float_plus_up prec (nprt a2 * pprt b1) (nprt a1 * nprt b1))))" |
|
206 |
|
207 lemma bnds_mult: |
|
208 fixes prec :: nat and a1 aa2 b1 b2 :: float |
|
209 assumes "(l, u) = bnds_mult prec a1 a2 b1 b2" |
|
210 assumes "a \<in> {real_of_float a1..real_of_float a2}" |
|
211 assumes "b \<in> {real_of_float b1..real_of_float b2}" |
|
212 shows "a * b \<in> {real_of_float l..real_of_float u}" |
|
213 proof - |
|
214 from assms have "real_of_float l \<le> a * b" |
|
215 by (intro order.trans[OF _ mult_ge_prts[of a1 a a2 b1 b b2]]) |
|
216 (auto simp: bnds_mult_def intro!: float_plus_down_le) |
|
217 moreover from assms have "real_of_float u \<ge> a * b" |
|
218 by (intro order.trans[OF mult_le_prts[of a1 a a2 b1 b b2]]) |
|
219 (auto simp: bnds_mult_def intro!: float_plus_up_le) |
|
220 ultimately show ?thesis by simp |
|
221 qed |
|
222 |
|
223 definition map_bnds :: "(nat \<Rightarrow> float \<Rightarrow> float) \<Rightarrow> (nat \<Rightarrow> float \<Rightarrow> float) \<Rightarrow> |
|
224 nat \<Rightarrow> (float \<times> float) \<Rightarrow> (float \<times> float)" where |
|
225 "map_bnds lb ub prec = (\<lambda>(l,u). (lb prec l, ub prec u))" |
|
226 |
|
227 lemma map_bnds: |
|
228 assumes "(lf, uf) = map_bnds lb ub prec (l, u)" |
|
229 assumes "mono f" |
|
230 assumes "x \<in> {real_of_float l..real_of_float u}" |
|
231 assumes "real_of_float (lb prec l) \<le> f (real_of_float l)" |
|
232 assumes "real_of_float (ub prec u) \<ge> f (real_of_float u)" |
|
233 shows "f x \<in> {real_of_float lf..real_of_float uf}" |
|
234 proof - |
|
235 from assms have "real_of_float lf = real_of_float (lb prec l)" |
|
236 by (simp add: map_bnds_def) |
|
237 also have "real_of_float (lb prec l) \<le> f (real_of_float l)" by fact |
|
238 also from assms have "\<dots> \<le> f x" |
|
239 by (intro monoD[OF \<open>mono f\<close>]) auto |
|
240 finally have lf: "real_of_float lf \<le> f x" . |
|
241 |
|
242 from assms have "f x \<le> f (real_of_float u)" |
|
243 by (intro monoD[OF \<open>mono f\<close>]) auto |
|
244 also have "\<dots> \<le> real_of_float (ub prec u)" by fact |
|
245 also from assms have "\<dots> = real_of_float uf" |
|
246 by (simp add: map_bnds_def) |
|
247 finally have uf: "f x \<le> real_of_float uf" . |
|
248 |
|
249 from lf uf show ?thesis by simp |
|
250 qed |
|
251 |
|
252 |
|
253 section "Square root" |
|
254 |
|
255 text \<open> |
|
256 The square root computation is implemented as newton iteration. As first first step we use the |
|
257 nearest power of two greater than the square root. |
|
258 \<close> |
|
259 |
|
260 fun sqrt_iteration :: "nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where |
|
261 "sqrt_iteration prec 0 x = Float 1 ((bitlen \<bar>mantissa x\<bar> + exponent x) div 2 + 1)" | |
|
262 "sqrt_iteration prec (Suc m) x = (let y = sqrt_iteration prec m x |
|
263 in Float 1 (- 1) * float_plus_up prec y (float_divr prec x y))" |
|
264 |
|
265 lemma compute_sqrt_iteration_base[code]: |
|
266 shows "sqrt_iteration prec n (Float m e) = |
|
267 (if n = 0 then Float 1 ((if m = 0 then 0 else bitlen \<bar>m\<bar> + e) div 2 + 1) |
|
268 else (let y = sqrt_iteration prec (n - 1) (Float m e) in |
|
269 Float 1 (- 1) * float_plus_up prec y (float_divr prec (Float m e) y)))" |
|
270 using bitlen_Float by (cases n) simp_all |
|
271 |
|
272 function ub_sqrt lb_sqrt :: "nat \<Rightarrow> float \<Rightarrow> float" where |
|
273 "ub_sqrt prec x = (if 0 < x then (sqrt_iteration prec prec x) |
|
274 else if x < 0 then - lb_sqrt prec (- x) |
|
275 else 0)" | |
|
276 "lb_sqrt prec x = (if 0 < x then (float_divl prec x (sqrt_iteration prec prec x)) |
|
277 else if x < 0 then - ub_sqrt prec (- x) |
|
278 else 0)" |
|
279 by pat_completeness auto |
|
280 termination by (relation "measure (\<lambda> v. let (prec, x) = case_sum id id v in (if x < 0 then 1 else 0))", auto) |
|
281 |
|
282 declare lb_sqrt.simps[simp del] |
|
283 declare ub_sqrt.simps[simp del] |
|
284 |
|
285 lemma sqrt_ub_pos_pos_1: |
|
286 assumes "sqrt x < b" and "0 < b" and "0 < x" |
|
287 shows "sqrt x < (b + x / b)/2" |
|
288 proof - |
|
289 from assms have "0 < (b - sqrt x)\<^sup>2 " by simp |
|
290 also have "\<dots> = b\<^sup>2 - 2 * b * sqrt x + (sqrt x)\<^sup>2" by algebra |
|
291 also have "\<dots> = b\<^sup>2 - 2 * b * sqrt x + x" using assms by simp |
|
292 finally have "0 < b\<^sup>2 - 2 * b * sqrt x + x" . |
|
293 hence "0 < b / 2 - sqrt x + x / (2 * b)" using assms |
|
294 by (simp add: field_simps power2_eq_square) |
|
295 thus ?thesis by (simp add: field_simps) |
|
296 qed |
|
297 |
|
298 lemma sqrt_iteration_bound: |
|
299 assumes "0 < real_of_float x" |
|
300 shows "sqrt x < sqrt_iteration prec n x" |
|
301 proof (induct n) |
|
302 case 0 |
|
303 show ?case |
|
304 proof (cases x) |
|
305 case (Float m e) |
|
306 hence "0 < m" |
|
307 using assms |
|
308 apply (auto simp: sign_simps) |
|
309 by (meson not_less powr_ge_pzero) |
|
310 hence "0 < sqrt m" by auto |
|
311 |
|
312 have int_nat_bl: "(nat (bitlen m)) = bitlen m" |
|
313 using bitlen_nonneg by auto |
|
314 |
|
315 have "x = (m / 2^nat (bitlen m)) * 2 powr (e + (nat (bitlen m)))" |
|
316 unfolding Float by (auto simp: powr_realpow[symmetric] field_simps powr_add) |
|
317 also have "\<dots> < 1 * 2 powr (e + nat (bitlen m))" |
|
318 proof (rule mult_strict_right_mono, auto) |
|
319 show "m < 2^nat (bitlen m)" |
|
320 using bitlen_bounds[OF \<open>0 < m\<close>, THEN conjunct2] |
|
321 unfolding of_int_less_iff[of m, symmetric] by auto |
|
322 qed |
|
323 finally have "sqrt x < sqrt (2 powr (e + bitlen m))" |
|
324 unfolding int_nat_bl by auto |
|
325 also have "\<dots> \<le> 2 powr ((e + bitlen m) div 2 + 1)" |
|
326 proof - |
|
327 let ?E = "e + bitlen m" |
|
328 have E_mod_pow: "2 powr (?E mod 2) < 4" |
|
329 proof (cases "?E mod 2 = 1") |
|
330 case True |
|
331 thus ?thesis by auto |
|
332 next |
|
333 case False |
|
334 have "0 \<le> ?E mod 2" by auto |
|
335 have "?E mod 2 < 2" by auto |
|
336 from this[THEN zless_imp_add1_zle] |
|
337 have "?E mod 2 \<le> 0" using False by auto |
|
338 from xt1(5)[OF \<open>0 \<le> ?E mod 2\<close> this] |
|
339 show ?thesis by auto |
|
340 qed |
|
341 hence "sqrt (2 powr (?E mod 2)) < sqrt (2 * 2)" |
|
342 by (auto simp del: real_sqrt_four) |
|
343 hence E_mod_pow: "sqrt (2 powr (?E mod 2)) < 2" by auto |
|
344 |
|
345 have E_eq: "2 powr ?E = 2 powr (?E div 2 + ?E div 2 + ?E mod 2)" |
|
346 by auto |
|
347 have "sqrt (2 powr ?E) = sqrt (2 powr (?E div 2) * 2 powr (?E div 2) * 2 powr (?E mod 2))" |
|
348 unfolding E_eq unfolding powr_add[symmetric] by (metis of_int_add) |
|
349 also have "\<dots> = 2 powr (?E div 2) * sqrt (2 powr (?E mod 2))" |
|
350 unfolding real_sqrt_mult[of _ "2 powr (?E mod 2)"] real_sqrt_abs2 by auto |
|
351 also have "\<dots> < 2 powr (?E div 2) * 2 powr 1" |
|
352 by (rule mult_strict_left_mono) (auto intro: E_mod_pow) |
|
353 also have "\<dots> = 2 powr (?E div 2 + 1)" |
|
354 unfolding add.commute[of _ 1] powr_add[symmetric] by simp |
|
355 finally show ?thesis by auto |
|
356 qed |
|
357 finally show ?thesis using \<open>0 < m\<close> |
|
358 unfolding Float |
|
359 by (subst compute_sqrt_iteration_base) (simp add: ac_simps) |
|
360 qed |
|
361 next |
|
362 case (Suc n) |
|
363 let ?b = "sqrt_iteration prec n x" |
|
364 have "0 < sqrt x" |
|
365 using \<open>0 < real_of_float x\<close> by auto |
|
366 also have "\<dots> < real_of_float ?b" |
|
367 using Suc . |
|
368 finally have "sqrt x < (?b + x / ?b)/2" |
|
369 using sqrt_ub_pos_pos_1[OF Suc _ \<open>0 < real_of_float x\<close>] by auto |
|
370 also have "\<dots> \<le> (?b + (float_divr prec x ?b))/2" |
|
371 by (rule divide_right_mono, auto simp add: float_divr) |
|
372 also have "\<dots> = (Float 1 (- 1)) * (?b + (float_divr prec x ?b))" |
|
373 by simp |
|
374 also have "\<dots> \<le> (Float 1 (- 1)) * (float_plus_up prec ?b (float_divr prec x ?b))" |
|
375 by (auto simp add: algebra_simps float_plus_up_le) |
|
376 finally show ?case |
|
377 unfolding sqrt_iteration.simps Let_def distrib_left . |
|
378 qed |
|
379 |
|
380 lemma sqrt_iteration_lower_bound: |
|
381 assumes "0 < real_of_float x" |
|
382 shows "0 < real_of_float (sqrt_iteration prec n x)" (is "0 < ?sqrt") |
|
383 proof - |
|
384 have "0 < sqrt x" using assms by auto |
|
385 also have "\<dots> < ?sqrt" using sqrt_iteration_bound[OF assms] . |
|
386 finally show ?thesis . |
|
387 qed |
|
388 |
|
389 lemma lb_sqrt_lower_bound: |
|
390 assumes "0 \<le> real_of_float x" |
|
391 shows "0 \<le> real_of_float (lb_sqrt prec x)" |
|
392 proof (cases "0 < x") |
|
393 case True |
|
394 hence "0 < real_of_float x" and "0 \<le> x" |
|
395 using \<open>0 \<le> real_of_float x\<close> by auto |
|
396 hence "0 < sqrt_iteration prec prec x" |
|
397 using sqrt_iteration_lower_bound by auto |
|
398 hence "0 \<le> real_of_float (float_divl prec x (sqrt_iteration prec prec x))" |
|
399 using float_divl_lower_bound[OF \<open>0 \<le> x\<close>] unfolding less_eq_float_def by auto |
|
400 thus ?thesis |
|
401 unfolding lb_sqrt.simps using True by auto |
|
402 next |
|
403 case False |
|
404 with \<open>0 \<le> real_of_float x\<close> have "real_of_float x = 0" by auto |
|
405 thus ?thesis |
|
406 unfolding lb_sqrt.simps by auto |
|
407 qed |
|
408 |
|
409 lemma bnds_sqrt': "sqrt x \<in> {(lb_sqrt prec x) .. (ub_sqrt prec x)}" |
|
410 proof - |
|
411 have lb: "lb_sqrt prec x \<le> sqrt x" if "0 < x" for x :: float |
|
412 proof - |
|
413 from that have "0 < real_of_float x" and "0 \<le> real_of_float x" by auto |
|
414 hence sqrt_gt0: "0 < sqrt x" by auto |
|
415 hence sqrt_ub: "sqrt x < sqrt_iteration prec prec x" |
|
416 using sqrt_iteration_bound by auto |
|
417 have "(float_divl prec x (sqrt_iteration prec prec x)) \<le> |
|
418 x / (sqrt_iteration prec prec x)" by (rule float_divl) |
|
419 also have "\<dots> < x / sqrt x" |
|
420 by (rule divide_strict_left_mono[OF sqrt_ub \<open>0 < real_of_float x\<close> |
|
421 mult_pos_pos[OF order_less_trans[OF sqrt_gt0 sqrt_ub] sqrt_gt0]]) |
|
422 also have "\<dots> = sqrt x" |
|
423 unfolding inverse_eq_iff_eq[of _ "sqrt x", symmetric] |
|
424 sqrt_divide_self_eq[OF \<open>0 \<le> real_of_float x\<close>, symmetric] by auto |
|
425 finally show ?thesis |
|
426 unfolding lb_sqrt.simps if_P[OF \<open>0 < x\<close>] by auto |
|
427 qed |
|
428 have ub: "sqrt x \<le> ub_sqrt prec x" if "0 < x" for x :: float |
|
429 proof - |
|
430 from that have "0 < real_of_float x" by auto |
|
431 hence "0 < sqrt x" by auto |
|
432 hence "sqrt x < sqrt_iteration prec prec x" |
|
433 using sqrt_iteration_bound by auto |
|
434 then show ?thesis |
|
435 unfolding ub_sqrt.simps if_P[OF \<open>0 < x\<close>] by auto |
|
436 qed |
|
437 show ?thesis |
|
438 using lb[of "-x"] ub[of "-x"] lb[of x] ub[of x] |
|
439 by (auto simp add: lb_sqrt.simps ub_sqrt.simps real_sqrt_minus) |
|
440 qed |
|
441 |
|
442 lemma bnds_sqrt: "\<forall>(x::real) lx ux. |
|
443 (l, u) = (lb_sqrt prec lx, ub_sqrt prec ux) \<and> x \<in> {lx .. ux} \<longrightarrow> l \<le> sqrt x \<and> sqrt x \<le> u" |
|
444 proof ((rule allI) +, rule impI, erule conjE, rule conjI) |
|
445 fix x :: real |
|
446 fix lx ux |
|
447 assume "(l, u) = (lb_sqrt prec lx, ub_sqrt prec ux)" |
|
448 and x: "x \<in> {lx .. ux}" |
|
449 hence l: "l = lb_sqrt prec lx " and u: "u = ub_sqrt prec ux" by auto |
|
450 |
|
451 have "sqrt lx \<le> sqrt x" using x by auto |
|
452 from order_trans[OF _ this] |
|
453 show "l \<le> sqrt x" unfolding l using bnds_sqrt'[of lx prec] by auto |
|
454 |
|
455 have "sqrt x \<le> sqrt ux" using x by auto |
|
456 from order_trans[OF this] |
|
457 show "sqrt x \<le> u" unfolding u using bnds_sqrt'[of ux prec] by auto |
|
458 qed |
|
459 |
|
460 |
|
461 section "Arcus tangens and \<pi>" |
|
462 |
|
463 subsection "Compute arcus tangens series" |
|
464 |
|
465 text \<open> |
|
466 As first step we implement the computation of the arcus tangens series. This is only valid in the range |
|
467 @{term "{-1 :: real .. 1}"}. This is used to compute \<pi> and then the entire arcus tangens. |
|
468 \<close> |
|
469 |
|
470 fun ub_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" |
|
471 and lb_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where |
|
472 "ub_arctan_horner prec 0 k x = 0" |
|
473 | "ub_arctan_horner prec (Suc n) k x = float_plus_up prec |
|
474 (rapprox_rat prec 1 k) (- float_round_down prec (x * (lb_arctan_horner prec n (k + 2) x)))" |
|
475 | "lb_arctan_horner prec 0 k x = 0" |
|
476 | "lb_arctan_horner prec (Suc n) k x = float_plus_down prec |
|
477 (lapprox_rat prec 1 k) (- float_round_up prec (x * (ub_arctan_horner prec n (k + 2) x)))" |
|
478 |
|
479 lemma arctan_0_1_bounds': |
|
480 assumes "0 \<le> real_of_float y" "real_of_float y \<le> 1" |
|
481 and "even n" |
|
482 shows "arctan (sqrt y) \<in> |
|
483 {(sqrt y * lb_arctan_horner prec n 1 y) .. (sqrt y * ub_arctan_horner prec (Suc n) 1 y)}" |
|
484 proof - |
|
485 let ?c = "\<lambda>i. (- 1) ^ i * (1 / (i * 2 + (1::nat)) * sqrt y ^ (i * 2 + 1))" |
|
486 let ?S = "\<lambda>n. \<Sum> i=0..<n. ?c i" |
|
487 |
|
488 have "0 \<le> sqrt y" using assms by auto |
|
489 have "sqrt y \<le> 1" using assms by auto |
|
490 from \<open>even n\<close> obtain m where "2 * m = n" by (blast elim: evenE) |
|
491 |
|
492 have "arctan (sqrt y) \<in> { ?S n .. ?S (Suc n) }" |
|
493 proof (cases "sqrt y = 0") |
|
494 case True |
|
495 then show ?thesis by simp |
|
496 next |
|
497 case False |
|
498 hence "0 < sqrt y" using \<open>0 \<le> sqrt y\<close> by auto |
|
499 hence prem: "0 < 1 / (0 * 2 + (1::nat)) * sqrt y ^ (0 * 2 + 1)" by auto |
|
500 |
|
501 have "\<bar> sqrt y \<bar> \<le> 1" using \<open>0 \<le> sqrt y\<close> \<open>sqrt y \<le> 1\<close> by auto |
|
502 from mp[OF summable_Leibniz(2)[OF zeroseq_arctan_series[OF this] |
|
503 monoseq_arctan_series[OF this]] prem, THEN spec, of m, unfolded \<open>2 * m = n\<close>] |
|
504 show ?thesis unfolding arctan_series[OF \<open>\<bar> sqrt y \<bar> \<le> 1\<close>] Suc_eq_plus1 atLeast0LessThan . |
|
505 qed |
|
506 note arctan_bounds = this[unfolded atLeastAtMost_iff] |
|
507 |
|
508 have F: "\<And>n. 2 * Suc n + 1 = 2 * n + 1 + 2" by auto |
|
509 |
|
510 note bounds = horner_bounds[where s=1 and f="\<lambda>i. 2 * i + 1" and j'=0 |
|
511 and lb="\<lambda>n i k x. lb_arctan_horner prec n k x" |
|
512 and ub="\<lambda>n i k x. ub_arctan_horner prec n k x", |
|
513 OF \<open>0 \<le> real_of_float y\<close> F lb_arctan_horner.simps ub_arctan_horner.simps] |
|
514 |
|
515 have "(sqrt y * lb_arctan_horner prec n 1 y) \<le> arctan (sqrt y)" |
|
516 proof - |
|
517 have "(sqrt y * lb_arctan_horner prec n 1 y) \<le> ?S n" |
|
518 using bounds(1) \<open>0 \<le> sqrt y\<close> |
|
519 apply (simp only: power_add power_one_right mult.assoc[symmetric] sum_distrib_right[symmetric]) |
|
520 apply (simp only: mult.commute[where 'a=real] mult.commute[of _ "2::nat"] power_mult) |
|
521 apply (auto intro!: mult_left_mono) |
|
522 done |
|
523 also have "\<dots> \<le> arctan (sqrt y)" using arctan_bounds .. |
|
524 finally show ?thesis . |
|
525 qed |
|
526 moreover |
|
527 have "arctan (sqrt y) \<le> (sqrt y * ub_arctan_horner prec (Suc n) 1 y)" |
|
528 proof - |
|
529 have "arctan (sqrt y) \<le> ?S (Suc n)" using arctan_bounds .. |
|
530 also have "\<dots> \<le> (sqrt y * ub_arctan_horner prec (Suc n) 1 y)" |
|
531 using bounds(2)[of "Suc n"] \<open>0 \<le> sqrt y\<close> |
|
532 apply (simp only: power_add power_one_right mult.assoc[symmetric] sum_distrib_right[symmetric]) |
|
533 apply (simp only: mult.commute[where 'a=real] mult.commute[of _ "2::nat"] power_mult) |
|
534 apply (auto intro!: mult_left_mono) |
|
535 done |
|
536 finally show ?thesis . |
|
537 qed |
|
538 ultimately show ?thesis by auto |
|
539 qed |
|
540 |
|
541 lemma arctan_0_1_bounds: |
|
542 assumes "0 \<le> real_of_float y" "real_of_float y \<le> 1" |
|
543 shows "arctan (sqrt y) \<in> |
|
544 {(sqrt y * lb_arctan_horner prec (get_even n) 1 y) .. |
|
545 (sqrt y * ub_arctan_horner prec (get_odd n) 1 y)}" |
|
546 using |
|
547 arctan_0_1_bounds'[OF assms, of n prec] |
|
548 arctan_0_1_bounds'[OF assms, of "n + 1" prec] |
|
549 arctan_0_1_bounds'[OF assms, of "n - 1" prec] |
|
550 by (auto simp: get_even_def get_odd_def odd_pos |
|
551 simp del: ub_arctan_horner.simps lb_arctan_horner.simps) |
|
552 |
|
553 lemma arctan_lower_bound: |
|
554 assumes "0 \<le> x" |
|
555 shows "x / (1 + x\<^sup>2) \<le> arctan x" (is "?l x \<le> _") |
|
556 proof - |
|
557 have "?l x - arctan x \<le> ?l 0 - arctan 0" |
|
558 using assms |
|
559 by (intro DERIV_nonpos_imp_nonincreasing[where f="\<lambda>x. ?l x - arctan x"]) |
|
560 (auto intro!: derivative_eq_intros simp: add_nonneg_eq_0_iff field_simps) |
|
561 thus ?thesis by simp |
|
562 qed |
|
563 |
|
564 lemma arctan_divide_mono: "0 < x \<Longrightarrow> x \<le> y \<Longrightarrow> arctan y / y \<le> arctan x / x" |
|
565 by (rule DERIV_nonpos_imp_nonincreasing[where f="\<lambda>x. arctan x / x"]) |
|
566 (auto intro!: derivative_eq_intros divide_nonpos_nonneg |
|
567 simp: inverse_eq_divide arctan_lower_bound) |
|
568 |
|
569 lemma arctan_mult_mono: "0 \<le> x \<Longrightarrow> x \<le> y \<Longrightarrow> x * arctan y \<le> y * arctan x" |
|
570 using arctan_divide_mono[of x y] by (cases "x = 0") (simp_all add: field_simps) |
|
571 |
|
572 lemma arctan_mult_le: |
|
573 assumes "0 \<le> x" "x \<le> y" "y * z \<le> arctan y" |
|
574 shows "x * z \<le> arctan x" |
|
575 proof (cases "x = 0") |
|
576 case True |
|
577 then show ?thesis by simp |
|
578 next |
|
579 case False |
|
580 with assms have "z \<le> arctan y / y" by (simp add: field_simps) |
|
581 also have "\<dots> \<le> arctan x / x" using assms \<open>x \<noteq> 0\<close> by (auto intro!: arctan_divide_mono) |
|
582 finally show ?thesis using assms \<open>x \<noteq> 0\<close> by (simp add: field_simps) |
|
583 qed |
|
584 |
|
585 lemma arctan_le_mult: |
|
586 assumes "0 < x" "x \<le> y" "arctan x \<le> x * z" |
|
587 shows "arctan y \<le> y * z" |
|
588 proof - |
|
589 from assms have "arctan y / y \<le> arctan x / x" by (auto intro!: arctan_divide_mono) |
|
590 also have "\<dots> \<le> z" using assms by (auto simp: field_simps) |
|
591 finally show ?thesis using assms by (simp add: field_simps) |
|
592 qed |
|
593 |
|
594 lemma arctan_0_1_bounds_le: |
|
595 assumes "0 \<le> x" "x \<le> 1" "0 < real_of_float xl" "real_of_float xl \<le> x * x" "x * x \<le> real_of_float xu" "real_of_float xu \<le> 1" |
|
596 shows "arctan x \<in> |
|
597 {x * lb_arctan_horner p1 (get_even n) 1 xu .. x * ub_arctan_horner p2 (get_odd n) 1 xl}" |
|
598 proof - |
|
599 from assms have "real_of_float xl \<le> 1" "sqrt (real_of_float xl) \<le> x" "x \<le> sqrt (real_of_float xu)" "0 \<le> real_of_float xu" |
|
600 "0 \<le> real_of_float xl" "0 < sqrt (real_of_float xl)" |
|
601 by (auto intro!: real_le_rsqrt real_le_lsqrt simp: power2_eq_square) |
|
602 from arctan_0_1_bounds[OF \<open>0 \<le> real_of_float xu\<close> \<open>real_of_float xu \<le> 1\<close>] |
|
603 have "sqrt (real_of_float xu) * real_of_float (lb_arctan_horner p1 (get_even n) 1 xu) \<le> arctan (sqrt (real_of_float xu))" |
|
604 by simp |
|
605 from arctan_mult_le[OF \<open>0 \<le> x\<close> \<open>x \<le> sqrt _\<close> this] |
|
606 have "x * real_of_float (lb_arctan_horner p1 (get_even n) 1 xu) \<le> arctan x" . |
|
607 moreover |
|
608 from arctan_0_1_bounds[OF \<open>0 \<le> real_of_float xl\<close> \<open>real_of_float xl \<le> 1\<close>] |
|
609 have "arctan (sqrt (real_of_float xl)) \<le> sqrt (real_of_float xl) * real_of_float (ub_arctan_horner p2 (get_odd n) 1 xl)" |
|
610 by simp |
|
611 from arctan_le_mult[OF \<open>0 < sqrt xl\<close> \<open>sqrt xl \<le> x\<close> this] |
|
612 have "arctan x \<le> x * real_of_float (ub_arctan_horner p2 (get_odd n) 1 xl)" . |
|
613 ultimately show ?thesis by simp |
|
614 qed |
|
615 |
|
616 lemma arctan_0_1_bounds_round: |
|
617 assumes "0 \<le> real_of_float x" "real_of_float x \<le> 1" |
|
618 shows "arctan x \<in> |
|
619 {real_of_float x * lb_arctan_horner p1 (get_even n) 1 (float_round_up (Suc p2) (x * x)) .. |
|
620 real_of_float x * ub_arctan_horner p3 (get_odd n) 1 (float_round_down (Suc p4) (x * x))}" |
|
621 using assms |
|
622 apply (cases "x > 0") |
|
623 apply (intro arctan_0_1_bounds_le) |
|
624 apply (auto simp: float_round_down.rep_eq float_round_up.rep_eq |
|
625 intro!: truncate_up_le1 mult_le_one truncate_down_le truncate_up_le truncate_down_pos |
|
626 mult_pos_pos) |
|
627 done |
|
628 |
|
629 |
|
630 subsection "Compute \<pi>" |
|
631 |
|
632 definition ub_pi :: "nat \<Rightarrow> float" where |
|
633 "ub_pi prec = |
|
634 (let |
|
635 A = rapprox_rat prec 1 5 ; |
|
636 B = lapprox_rat prec 1 239 |
|
637 in ((Float 1 2) * float_plus_up prec |
|
638 ((Float 1 2) * float_round_up prec (A * (ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 |
|
639 (float_round_down (Suc prec) (A * A))))) |
|
640 (- float_round_down prec (B * (lb_arctan_horner prec (get_even (prec div 14 + 1)) 1 |
|
641 (float_round_up (Suc prec) (B * B)))))))" |
|
642 |
|
643 definition lb_pi :: "nat \<Rightarrow> float" where |
|
644 "lb_pi prec = |
|
645 (let |
|
646 A = lapprox_rat prec 1 5 ; |
|
647 B = rapprox_rat prec 1 239 |
|
648 in ((Float 1 2) * float_plus_down prec |
|
649 ((Float 1 2) * float_round_down prec (A * (lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 |
|
650 (float_round_up (Suc prec) (A * A))))) |
|
651 (- float_round_up prec (B * (ub_arctan_horner prec (get_odd (prec div 14 + 1)) 1 |
|
652 (float_round_down (Suc prec) (B * B)))))))" |
|
653 |
|
654 lemma pi_boundaries: "pi \<in> {(lb_pi n) .. (ub_pi n)}" |
|
655 proof - |
|
656 have machin_pi: "pi = 4 * (4 * arctan (1 / 5) - arctan (1 / 239))" |
|
657 unfolding machin[symmetric] by auto |
|
658 |
|
659 { |
|
660 fix prec n :: nat |
|
661 fix k :: int |
|
662 assume "1 < k" hence "0 \<le> k" and "0 < k" and "1 \<le> k" by auto |
|
663 let ?k = "rapprox_rat prec 1 k" |
|
664 let ?kl = "float_round_down (Suc prec) (?k * ?k)" |
|
665 have "1 div k = 0" using div_pos_pos_trivial[OF _ \<open>1 < k\<close>] by auto |
|
666 |
|
667 have "0 \<le> real_of_float ?k" by (rule order_trans[OF _ rapprox_rat]) (auto simp add: \<open>0 \<le> k\<close>) |
|
668 have "real_of_float ?k \<le> 1" |
|
669 by (auto simp add: \<open>0 < k\<close> \<open>1 \<le> k\<close> less_imp_le |
|
670 intro!: mult_le_one order_trans[OF _ rapprox_rat] rapprox_rat_le1) |
|
671 have "1 / k \<le> ?k" using rapprox_rat[where x=1 and y=k] by auto |
|
672 hence "arctan (1 / k) \<le> arctan ?k" by (rule arctan_monotone') |
|
673 also have "\<dots> \<le> (?k * ub_arctan_horner prec (get_odd n) 1 ?kl)" |
|
674 using arctan_0_1_bounds_round[OF \<open>0 \<le> real_of_float ?k\<close> \<open>real_of_float ?k \<le> 1\<close>] |
|
675 by auto |
|
676 finally have "arctan (1 / k) \<le> ?k * ub_arctan_horner prec (get_odd n) 1 ?kl" . |
|
677 } note ub_arctan = this |
|
678 |
|
679 { |
|
680 fix prec n :: nat |
|
681 fix k :: int |
|
682 assume "1 < k" hence "0 \<le> k" and "0 < k" by auto |
|
683 let ?k = "lapprox_rat prec 1 k" |
|
684 let ?ku = "float_round_up (Suc prec) (?k * ?k)" |
|
685 have "1 div k = 0" using div_pos_pos_trivial[OF _ \<open>1 < k\<close>] by auto |
|
686 have "1 / k \<le> 1" using \<open>1 < k\<close> by auto |
|
687 have "0 \<le> real_of_float ?k" using lapprox_rat_nonneg[where x=1 and y=k, OF zero_le_one \<open>0 \<le> k\<close>] |
|
688 by (auto simp add: \<open>1 div k = 0\<close>) |
|
689 have "0 \<le> real_of_float (?k * ?k)" by simp |
|
690 have "real_of_float ?k \<le> 1" using lapprox_rat by (rule order_trans, auto simp add: \<open>1 / k \<le> 1\<close>) |
|
691 hence "real_of_float (?k * ?k) \<le> 1" using \<open>0 \<le> real_of_float ?k\<close> by (auto intro!: mult_le_one) |
|
692 |
|
693 have "?k \<le> 1 / k" using lapprox_rat[where x=1 and y=k] by auto |
|
694 |
|
695 have "?k * lb_arctan_horner prec (get_even n) 1 ?ku \<le> arctan ?k" |
|
696 using arctan_0_1_bounds_round[OF \<open>0 \<le> real_of_float ?k\<close> \<open>real_of_float ?k \<le> 1\<close>] |
|
697 by auto |
|
698 also have "\<dots> \<le> arctan (1 / k)" using \<open>?k \<le> 1 / k\<close> by (rule arctan_monotone') |
|
699 finally have "?k * lb_arctan_horner prec (get_even n) 1 ?ku \<le> arctan (1 / k)" . |
|
700 } note lb_arctan = this |
|
701 |
|
702 have "pi \<le> ub_pi n " |
|
703 unfolding ub_pi_def machin_pi Let_def times_float.rep_eq Float_num |
|
704 using lb_arctan[of 239] ub_arctan[of 5] powr_realpow[of 2 2] |
|
705 by (intro mult_left_mono float_plus_up_le float_plus_down_le) |
|
706 (auto intro!: mult_left_mono float_round_down_le float_round_up_le diff_mono) |
|
707 moreover have "lb_pi n \<le> pi" |
|
708 unfolding lb_pi_def machin_pi Let_def times_float.rep_eq Float_num |
|
709 using lb_arctan[of 5] ub_arctan[of 239] |
|
710 by (intro mult_left_mono float_plus_up_le float_plus_down_le) |
|
711 (auto intro!: mult_left_mono float_round_down_le float_round_up_le diff_mono) |
|
712 ultimately show ?thesis by auto |
|
713 qed |
|
714 |
|
715 |
|
716 subsection "Compute arcus tangens in the entire domain" |
|
717 |
|
718 function lb_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" and ub_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" where |
|
719 "lb_arctan prec x = |
|
720 (let |
|
721 ub_horner = \<lambda> x. float_round_up prec |
|
722 (x * |
|
723 ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (float_round_down (Suc prec) (x * x))); |
|
724 lb_horner = \<lambda> x. float_round_down prec |
|
725 (x * |
|
726 lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (float_round_up (Suc prec) (x * x))) |
|
727 in |
|
728 if x < 0 then - ub_arctan prec (-x) |
|
729 else if x \<le> Float 1 (- 1) then lb_horner x |
|
730 else if x \<le> Float 1 1 then |
|
731 Float 1 1 * |
|
732 lb_horner |
|
733 (float_divl prec x |
|
734 (float_plus_up prec 1 |
|
735 (ub_sqrt prec (float_plus_up prec 1 (float_round_up prec (x * x)))))) |
|
736 else let inv = float_divr prec 1 x in |
|
737 if inv > 1 then 0 |
|
738 else float_plus_down prec (lb_pi prec * Float 1 (- 1)) ( - ub_horner inv))" |
|
739 |
|
740 | "ub_arctan prec x = |
|
741 (let |
|
742 lb_horner = \<lambda> x. float_round_down prec |
|
743 (x * |
|
744 lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (float_round_up (Suc prec) (x * x))) ; |
|
745 ub_horner = \<lambda> x. float_round_up prec |
|
746 (x * |
|
747 ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (float_round_down (Suc prec) (x * x))) |
|
748 in if x < 0 then - lb_arctan prec (-x) |
|
749 else if x \<le> Float 1 (- 1) then ub_horner x |
|
750 else if x \<le> Float 1 1 then |
|
751 let y = float_divr prec x |
|
752 (float_plus_down |
|
753 (Suc prec) 1 (lb_sqrt prec (float_plus_down prec 1 (float_round_down prec (x * x))))) |
|
754 in if y > 1 then ub_pi prec * Float 1 (- 1) else Float 1 1 * ub_horner y |
|
755 else float_plus_up prec (ub_pi prec * Float 1 (- 1)) ( - lb_horner (float_divl prec 1 x)))" |
|
756 by pat_completeness auto |
|
757 termination |
|
758 by (relation "measure (\<lambda> v. let (prec, x) = case_sum id id v in (if x < 0 then 1 else 0))", auto) |
|
759 |
|
760 declare ub_arctan_horner.simps[simp del] |
|
761 declare lb_arctan_horner.simps[simp del] |
|
762 |
|
763 lemma lb_arctan_bound': |
|
764 assumes "0 \<le> real_of_float x" |
|
765 shows "lb_arctan prec x \<le> arctan x" |
|
766 proof - |
|
767 have "\<not> x < 0" and "0 \<le> x" |
|
768 using \<open>0 \<le> real_of_float x\<close> by (auto intro!: truncate_up_le ) |
|
769 |
|
770 let "?ub_horner x" = |
|
771 "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (float_round_down (Suc prec) (x * x))" |
|
772 and "?lb_horner x" = |
|
773 "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (float_round_up (Suc prec) (x * x))" |
|
774 |
|
775 show ?thesis |
|
776 proof (cases "x \<le> Float 1 (- 1)") |
|
777 case True |
|
778 hence "real_of_float x \<le> 1" by simp |
|
779 from arctan_0_1_bounds_round[OF \<open>0 \<le> real_of_float x\<close> \<open>real_of_float x \<le> 1\<close>] |
|
780 show ?thesis |
|
781 unfolding lb_arctan.simps Let_def if_not_P[OF \<open>\<not> x < 0\<close>] if_P[OF True] using \<open>0 \<le> x\<close> |
|
782 by (auto intro!: float_round_down_le) |
|
783 next |
|
784 case False |
|
785 hence "0 < real_of_float x" by auto |
|
786 let ?R = "1 + sqrt (1 + real_of_float x * real_of_float x)" |
|
787 let ?sxx = "float_plus_up prec 1 (float_round_up prec (x * x))" |
|
788 let ?fR = "float_plus_up prec 1 (ub_sqrt prec ?sxx)" |
|
789 let ?DIV = "float_divl prec x ?fR" |
|
790 |
|
791 have divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg) |
|
792 |
|
793 have "sqrt (1 + x*x) \<le> sqrt ?sxx" |
|
794 by (auto simp: float_plus_up.rep_eq plus_up_def float_round_up.rep_eq intro!: truncate_up_le) |
|
795 also have "\<dots> \<le> ub_sqrt prec ?sxx" |
|
796 using bnds_sqrt'[of ?sxx prec] by auto |
|
797 finally |
|
798 have "sqrt (1 + x*x) \<le> ub_sqrt prec ?sxx" . |
|
799 hence "?R \<le> ?fR" by (auto simp: float_plus_up.rep_eq plus_up_def intro!: truncate_up_le) |
|
800 hence "0 < ?fR" and "0 < real_of_float ?fR" using \<open>0 < ?R\<close> by auto |
|
801 |
|
802 have monotone: "?DIV \<le> x / ?R" |
|
803 proof - |
|
804 have "?DIV \<le> real_of_float x / ?fR" by (rule float_divl) |
|
805 also have "\<dots> \<le> x / ?R" by (rule divide_left_mono[OF \<open>?R \<le> ?fR\<close> \<open>0 \<le> real_of_float x\<close> mult_pos_pos[OF order_less_le_trans[OF divisor_gt0 \<open>?R \<le> real_of_float ?fR\<close>] divisor_gt0]]) |
|
806 finally show ?thesis . |
|
807 qed |
|
808 |
|
809 show ?thesis |
|
810 proof (cases "x \<le> Float 1 1") |
|
811 case True |
|
812 have "x \<le> sqrt (1 + x * x)" |
|
813 using real_sqrt_sum_squares_ge2[where x=1, unfolded numeral_2_eq_2] by auto |
|
814 also note \<open>\<dots> \<le> (ub_sqrt prec ?sxx)\<close> |
|
815 finally have "real_of_float x \<le> ?fR" |
|
816 by (auto simp: float_plus_up.rep_eq plus_up_def intro!: truncate_up_le) |
|
817 moreover have "?DIV \<le> real_of_float x / ?fR" |
|
818 by (rule float_divl) |
|
819 ultimately have "real_of_float ?DIV \<le> 1" |
|
820 unfolding divide_le_eq_1_pos[OF \<open>0 < real_of_float ?fR\<close>, symmetric] by auto |
|
821 |
|
822 have "0 \<le> real_of_float ?DIV" |
|
823 using float_divl_lower_bound[OF \<open>0 \<le> x\<close>] \<open>0 < ?fR\<close> |
|
824 unfolding less_eq_float_def by auto |
|
825 |
|
826 from arctan_0_1_bounds_round[OF \<open>0 \<le> real_of_float (?DIV)\<close> \<open>real_of_float (?DIV) \<le> 1\<close>] |
|
827 have "Float 1 1 * ?lb_horner ?DIV \<le> 2 * arctan ?DIV" |
|
828 by simp |
|
829 also have "\<dots> \<le> 2 * arctan (x / ?R)" |
|
830 using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono arctan_monotone') |
|
831 also have "2 * arctan (x / ?R) = arctan x" |
|
832 using arctan_half[symmetric] unfolding numeral_2_eq_2 power_Suc2 power_0 mult_1_left . |
|
833 finally show ?thesis |
|
834 unfolding lb_arctan.simps Let_def if_not_P[OF \<open>\<not> x < 0\<close>] |
|
835 if_not_P[OF \<open>\<not> x \<le> Float 1 (- 1)\<close>] if_P[OF True] |
|
836 by (auto simp: float_round_down.rep_eq |
|
837 intro!: order_trans[OF mult_left_mono[OF truncate_down]]) |
|
838 next |
|
839 case False |
|
840 hence "2 < real_of_float x" by auto |
|
841 hence "1 \<le> real_of_float x" by auto |
|
842 |
|
843 let "?invx" = "float_divr prec 1 x" |
|
844 have "0 \<le> arctan x" using arctan_monotone'[OF \<open>0 \<le> real_of_float x\<close>] |
|
845 using arctan_tan[of 0, unfolded tan_zero] by auto |
|
846 |
|
847 show ?thesis |
|
848 proof (cases "1 < ?invx") |
|
849 case True |
|
850 show ?thesis |
|
851 unfolding lb_arctan.simps Let_def if_not_P[OF \<open>\<not> x < 0\<close>] |
|
852 if_not_P[OF \<open>\<not> x \<le> Float 1 (- 1)\<close>] if_not_P[OF False] if_P[OF True] |
|
853 using \<open>0 \<le> arctan x\<close> by auto |
|
854 next |
|
855 case False |
|
856 hence "real_of_float ?invx \<le> 1" by auto |
|
857 have "0 \<le> real_of_float ?invx" |
|
858 by (rule order_trans[OF _ float_divr]) (auto simp add: \<open>0 \<le> real_of_float x\<close>) |
|
859 |
|
860 have "1 / x \<noteq> 0" and "0 < 1 / x" |
|
861 using \<open>0 < real_of_float x\<close> by auto |
|
862 |
|
863 have "arctan (1 / x) \<le> arctan ?invx" |
|
864 unfolding one_float.rep_eq[symmetric] by (rule arctan_monotone', rule float_divr) |
|
865 also have "\<dots> \<le> ?ub_horner ?invx" |
|
866 using arctan_0_1_bounds_round[OF \<open>0 \<le> real_of_float ?invx\<close> \<open>real_of_float ?invx \<le> 1\<close>] |
|
867 by (auto intro!: float_round_up_le) |
|
868 also note float_round_up |
|
869 finally have "pi / 2 - float_round_up prec (?ub_horner ?invx) \<le> arctan x" |
|
870 using \<open>0 \<le> arctan x\<close> arctan_inverse[OF \<open>1 / x \<noteq> 0\<close>] |
|
871 unfolding sgn_pos[OF \<open>0 < 1 / real_of_float x\<close>] le_diff_eq by auto |
|
872 moreover |
|
873 have "lb_pi prec * Float 1 (- 1) \<le> pi / 2" |
|
874 unfolding Float_num times_divide_eq_right mult_1_left using pi_boundaries by simp |
|
875 ultimately |
|
876 show ?thesis |
|
877 unfolding lb_arctan.simps Let_def if_not_P[OF \<open>\<not> x < 0\<close>] |
|
878 if_not_P[OF \<open>\<not> x \<le> Float 1 (- 1)\<close>] if_not_P[OF \<open>\<not> x \<le> Float 1 1\<close>] if_not_P[OF False] |
|
879 by (auto intro!: float_plus_down_le) |
|
880 qed |
|
881 qed |
|
882 qed |
|
883 qed |
|
884 |
|
885 lemma ub_arctan_bound': |
|
886 assumes "0 \<le> real_of_float x" |
|
887 shows "arctan x \<le> ub_arctan prec x" |
|
888 proof - |
|
889 have "\<not> x < 0" and "0 \<le> x" |
|
890 using \<open>0 \<le> real_of_float x\<close> by auto |
|
891 |
|
892 let "?ub_horner x" = |
|
893 "float_round_up prec (x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (float_round_down (Suc prec) (x * x)))" |
|
894 let "?lb_horner x" = |
|
895 "float_round_down prec (x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (float_round_up (Suc prec) (x * x)))" |
|
896 |
|
897 show ?thesis |
|
898 proof (cases "x \<le> Float 1 (- 1)") |
|
899 case True |
|
900 hence "real_of_float x \<le> 1" by auto |
|
901 show ?thesis |
|
902 unfolding ub_arctan.simps Let_def if_not_P[OF \<open>\<not> x < 0\<close>] if_P[OF True] |
|
903 using arctan_0_1_bounds_round[OF \<open>0 \<le> real_of_float x\<close> \<open>real_of_float x \<le> 1\<close>] |
|
904 by (auto intro!: float_round_up_le) |
|
905 next |
|
906 case False |
|
907 hence "0 < real_of_float x" by auto |
|
908 let ?R = "1 + sqrt (1 + real_of_float x * real_of_float x)" |
|
909 let ?sxx = "float_plus_down prec 1 (float_round_down prec (x * x))" |
|
910 let ?fR = "float_plus_down (Suc prec) 1 (lb_sqrt prec ?sxx)" |
|
911 let ?DIV = "float_divr prec x ?fR" |
|
912 |
|
913 have sqr_ge0: "0 \<le> 1 + real_of_float x * real_of_float x" |
|
914 using sum_power2_ge_zero[of 1 "real_of_float x", unfolded numeral_2_eq_2] by auto |
|
915 hence "0 \<le> real_of_float (1 + x*x)" by auto |
|
916 |
|
917 hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg) |
|
918 |
|
919 have "lb_sqrt prec ?sxx \<le> sqrt ?sxx" |
|
920 using bnds_sqrt'[of ?sxx] by auto |
|
921 also have "\<dots> \<le> sqrt (1 + x*x)" |
|
922 by (auto simp: float_plus_down.rep_eq plus_down_def float_round_down.rep_eq truncate_down_le) |
|
923 finally have "lb_sqrt prec ?sxx \<le> sqrt (1 + x*x)" . |
|
924 hence "?fR \<le> ?R" |
|
925 by (auto simp: float_plus_down.rep_eq plus_down_def truncate_down_le) |
|
926 have "0 < real_of_float ?fR" |
|
927 by (auto simp: float_plus_down.rep_eq plus_down_def float_round_down.rep_eq |
|
928 intro!: truncate_down_ge1 lb_sqrt_lower_bound order_less_le_trans[OF zero_less_one] |
|
929 truncate_down_nonneg add_nonneg_nonneg) |
|
930 have monotone: "x / ?R \<le> (float_divr prec x ?fR)" |
|
931 proof - |
|
932 from divide_left_mono[OF \<open>?fR \<le> ?R\<close> \<open>0 \<le> real_of_float x\<close> mult_pos_pos[OF divisor_gt0 \<open>0 < real_of_float ?fR\<close>]] |
|
933 have "x / ?R \<le> x / ?fR" . |
|
934 also have "\<dots> \<le> ?DIV" by (rule float_divr) |
|
935 finally show ?thesis . |
|
936 qed |
|
937 |
|
938 show ?thesis |
|
939 proof (cases "x \<le> Float 1 1") |
|
940 case True |
|
941 show ?thesis |
|
942 proof (cases "?DIV > 1") |
|
943 case True |
|
944 have "pi / 2 \<le> ub_pi prec * Float 1 (- 1)" |
|
945 unfolding Float_num times_divide_eq_right mult_1_left using pi_boundaries by auto |
|
946 from order_less_le_trans[OF arctan_ubound this, THEN less_imp_le] |
|
947 show ?thesis |
|
948 unfolding ub_arctan.simps Let_def if_not_P[OF \<open>\<not> x < 0\<close>] |
|
949 if_not_P[OF \<open>\<not> x \<le> Float 1 (- 1)\<close>] if_P[OF \<open>x \<le> Float 1 1\<close>] if_P[OF True] . |
|
950 next |
|
951 case False |
|
952 hence "real_of_float ?DIV \<le> 1" by auto |
|
953 |
|
954 have "0 \<le> x / ?R" |
|
955 using \<open>0 \<le> real_of_float x\<close> \<open>0 < ?R\<close> unfolding zero_le_divide_iff by auto |
|
956 hence "0 \<le> real_of_float ?DIV" |
|
957 using monotone by (rule order_trans) |
|
958 |
|
959 have "arctan x = 2 * arctan (x / ?R)" |
|
960 using arctan_half unfolding numeral_2_eq_2 power_Suc2 power_0 mult_1_left . |
|
961 also have "\<dots> \<le> 2 * arctan (?DIV)" |
|
962 using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono) |
|
963 also have "\<dots> \<le> (Float 1 1 * ?ub_horner ?DIV)" unfolding Float_num |
|
964 using arctan_0_1_bounds_round[OF \<open>0 \<le> real_of_float ?DIV\<close> \<open>real_of_float ?DIV \<le> 1\<close>] |
|
965 by (auto intro!: float_round_up_le) |
|
966 finally show ?thesis |
|
967 unfolding ub_arctan.simps Let_def if_not_P[OF \<open>\<not> x < 0\<close>] |
|
968 if_not_P[OF \<open>\<not> x \<le> Float 1 (- 1)\<close>] if_P[OF \<open>x \<le> Float 1 1\<close>] if_not_P[OF False] . |
|
969 qed |
|
970 next |
|
971 case False |
|
972 hence "2 < real_of_float x" by auto |
|
973 hence "1 \<le> real_of_float x" by auto |
|
974 hence "0 < real_of_float x" by auto |
|
975 hence "0 < x" by auto |
|
976 |
|
977 let "?invx" = "float_divl prec 1 x" |
|
978 have "0 \<le> arctan x" |
|
979 using arctan_monotone'[OF \<open>0 \<le> real_of_float x\<close>] and arctan_tan[of 0, unfolded tan_zero] by auto |
|
980 |
|
981 have "real_of_float ?invx \<le> 1" |
|
982 unfolding less_float_def |
|
983 by (rule order_trans[OF float_divl]) |
|
984 (auto simp add: \<open>1 \<le> real_of_float x\<close> divide_le_eq_1_pos[OF \<open>0 < real_of_float x\<close>]) |
|
985 have "0 \<le> real_of_float ?invx" |
|
986 using \<open>0 < x\<close> by (intro float_divl_lower_bound) auto |
|
987 |
|
988 have "1 / x \<noteq> 0" and "0 < 1 / x" |
|
989 using \<open>0 < real_of_float x\<close> by auto |
|
990 |
|
991 have "(?lb_horner ?invx) \<le> arctan (?invx)" |
|
992 using arctan_0_1_bounds_round[OF \<open>0 \<le> real_of_float ?invx\<close> \<open>real_of_float ?invx \<le> 1\<close>] |
|
993 by (auto intro!: float_round_down_le) |
|
994 also have "\<dots> \<le> arctan (1 / x)" |
|
995 unfolding one_float.rep_eq[symmetric] by (rule arctan_monotone') (rule float_divl) |
|
996 finally have "arctan x \<le> pi / 2 - (?lb_horner ?invx)" |
|
997 using \<open>0 \<le> arctan x\<close> arctan_inverse[OF \<open>1 / x \<noteq> 0\<close>] |
|
998 unfolding sgn_pos[OF \<open>0 < 1 / x\<close>] le_diff_eq by auto |
|
999 moreover |
|
1000 have "pi / 2 \<le> ub_pi prec * Float 1 (- 1)" |
|
1001 unfolding Float_num times_divide_eq_right mult_1_right |
|
1002 using pi_boundaries by auto |
|
1003 ultimately |
|
1004 show ?thesis |
|
1005 unfolding ub_arctan.simps Let_def if_not_P[OF \<open>\<not> x < 0\<close>] |
|
1006 if_not_P[OF \<open>\<not> x \<le> Float 1 (- 1)\<close>] if_not_P[OF False] |
|
1007 by (auto intro!: float_round_up_le float_plus_up_le) |
|
1008 qed |
|
1009 qed |
|
1010 qed |
|
1011 |
|
1012 lemma arctan_boundaries: "arctan x \<in> {(lb_arctan prec x) .. (ub_arctan prec x)}" |
|
1013 proof (cases "0 \<le> x") |
|
1014 case True |
|
1015 hence "0 \<le> real_of_float x" by auto |
|
1016 show ?thesis |
|
1017 using ub_arctan_bound'[OF \<open>0 \<le> real_of_float x\<close>] lb_arctan_bound'[OF \<open>0 \<le> real_of_float x\<close>] |
|
1018 unfolding atLeastAtMost_iff by auto |
|
1019 next |
|
1020 case False |
|
1021 let ?mx = "-x" |
|
1022 from False have "x < 0" and "0 \<le> real_of_float ?mx" |
|
1023 by auto |
|
1024 hence bounds: "lb_arctan prec ?mx \<le> arctan ?mx \<and> arctan ?mx \<le> ub_arctan prec ?mx" |
|
1025 using ub_arctan_bound'[OF \<open>0 \<le> real_of_float ?mx\<close>] lb_arctan_bound'[OF \<open>0 \<le> real_of_float ?mx\<close>] by auto |
|
1026 show ?thesis |
|
1027 unfolding minus_float.rep_eq arctan_minus lb_arctan.simps[where x=x] |
|
1028 ub_arctan.simps[where x=x] Let_def if_P[OF \<open>x < 0\<close>] |
|
1029 unfolding atLeastAtMost_iff using bounds[unfolded minus_float.rep_eq arctan_minus] |
|
1030 by (simp add: arctan_minus) |
|
1031 qed |
|
1032 |
|
1033 lemma bnds_arctan: "\<forall> (x::real) lx ux. (l, u) = (lb_arctan prec lx, ub_arctan prec ux) \<and> x \<in> {lx .. ux} \<longrightarrow> l \<le> arctan x \<and> arctan x \<le> u" |
|
1034 proof (rule allI, rule allI, rule allI, rule impI) |
|
1035 fix x :: real |
|
1036 fix lx ux |
|
1037 assume "(l, u) = (lb_arctan prec lx, ub_arctan prec ux) \<and> x \<in> {lx .. ux}" |
|
1038 hence l: "lb_arctan prec lx = l " |
|
1039 and u: "ub_arctan prec ux = u" |
|
1040 and x: "x \<in> {lx .. ux}" |
|
1041 by auto |
|
1042 show "l \<le> arctan x \<and> arctan x \<le> u" |
|
1043 proof |
|
1044 show "l \<le> arctan x" |
|
1045 proof - |
|
1046 from arctan_boundaries[of lx prec, unfolded l] |
|
1047 have "l \<le> arctan lx" by (auto simp del: lb_arctan.simps) |
|
1048 also have "\<dots> \<le> arctan x" using x by (auto intro: arctan_monotone') |
|
1049 finally show ?thesis . |
|
1050 qed |
|
1051 show "arctan x \<le> u" |
|
1052 proof - |
|
1053 have "arctan x \<le> arctan ux" using x by (auto intro: arctan_monotone') |
|
1054 also have "\<dots> \<le> u" using arctan_boundaries[of ux prec, unfolded u] by (auto simp del: ub_arctan.simps) |
|
1055 finally show ?thesis . |
|
1056 qed |
|
1057 qed |
|
1058 qed |
|
1059 |
|
1060 |
|
1061 section "Sinus and Cosinus" |
|
1062 |
|
1063 subsection "Compute the cosinus and sinus series" |
|
1064 |
|
1065 fun ub_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" |
|
1066 and lb_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where |
|
1067 "ub_sin_cos_aux prec 0 i k x = 0" |
|
1068 | "ub_sin_cos_aux prec (Suc n) i k x = float_plus_up prec |
|
1069 (rapprox_rat prec 1 k) (- |
|
1070 float_round_down prec (x * (lb_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)))" |
|
1071 | "lb_sin_cos_aux prec 0 i k x = 0" |
|
1072 | "lb_sin_cos_aux prec (Suc n) i k x = float_plus_down prec |
|
1073 (lapprox_rat prec 1 k) (- |
|
1074 float_round_up prec (x * (ub_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)))" |
|
1075 |
|
1076 lemma cos_aux: |
|
1077 shows "(lb_sin_cos_aux prec n 1 1 (x * x)) \<le> (\<Sum> i=0..<n. (- 1) ^ i * (1/(fact (2 * i))) * x ^(2 * i))" (is "?lb") |
|
1078 and "(\<Sum> i=0..<n. (- 1) ^ i * (1/(fact (2 * i))) * x^(2 * i)) \<le> (ub_sin_cos_aux prec n 1 1 (x * x))" (is "?ub") |
|
1079 proof - |
|
1080 have "0 \<le> real_of_float (x * x)" by auto |
|
1081 let "?f n" = "fact (2 * n) :: nat" |
|
1082 have f_eq: "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) ^^ n) 1 * (((\<lambda>i. i + 2) ^^ n) 1 + 1)" for n |
|
1083 proof - |
|
1084 have "\<And>m. ((\<lambda>i. i + 2) ^^ n) m = m + 2 * n" by (induct n) auto |
|
1085 then show ?thesis by auto |
|
1086 qed |
|
1087 from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0, |
|
1088 OF \<open>0 \<le> real_of_float (x * x)\<close> f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps] |
|
1089 show ?lb and ?ub |
|
1090 by (auto simp add: power_mult power2_eq_square[of "real_of_float x"]) |
|
1091 qed |
|
1092 |
|
1093 lemma lb_sin_cos_aux_zero_le_one: "lb_sin_cos_aux prec n i j 0 \<le> 1" |
|
1094 by (cases j n rule: nat.exhaust[case_product nat.exhaust]) |
|
1095 (auto intro!: float_plus_down_le order_trans[OF lapprox_rat]) |
|
1096 |
|
1097 lemma one_le_ub_sin_cos_aux: "odd n \<Longrightarrow> 1 \<le> ub_sin_cos_aux prec n i (Suc 0) 0" |
|
1098 by (cases n) (auto intro!: float_plus_up_le order_trans[OF _ rapprox_rat]) |
|
1099 |
|
1100 lemma cos_boundaries: |
|
1101 assumes "0 \<le> real_of_float x" and "x \<le> pi / 2" |
|
1102 shows "cos x \<in> {(lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) .. (ub_sin_cos_aux prec (get_odd n) 1 1 (x * x))}" |
|
1103 proof (cases "real_of_float x = 0") |
|
1104 case False |
|
1105 hence "real_of_float x \<noteq> 0" by auto |
|
1106 hence "0 < x" and "0 < real_of_float x" |
|
1107 using \<open>0 \<le> real_of_float x\<close> by auto |
|
1108 have "0 < x * x" |
|
1109 using \<open>0 < x\<close> by simp |
|
1110 |
|
1111 have morph_to_if_power: "(\<Sum> i=0..<n. (-1::real) ^ i * (1/(fact (2 * i))) * x ^ (2 * i)) = |
|
1112 (\<Sum> i = 0 ..< 2 * n. (if even(i) then ((- 1) ^ (i div 2))/((fact i)) else 0) * x ^ i)" |
|
1113 (is "?sum = ?ifsum") for x n |
|
1114 proof - |
|
1115 have "?sum = ?sum + (\<Sum> j = 0 ..< n. 0)" by auto |
|
1116 also have "\<dots> = |
|
1117 (\<Sum> j = 0 ..< n. (- 1) ^ ((2 * j) div 2) / ((fact (2 * j))) * x ^(2 * j)) + (\<Sum> j = 0 ..< n. 0)" by auto |
|
1118 also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. if even i then (- 1) ^ (i div 2) / ((fact i)) * x ^ i else 0)" |
|
1119 unfolding sum_split_even_odd atLeast0LessThan .. |
|
1120 also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. (if even i then (- 1) ^ (i div 2) / ((fact i)) else 0) * x ^ i)" |
|
1121 by (rule sum.cong) auto |
|
1122 finally show ?thesis . |
|
1123 qed |
|
1124 |
|
1125 { fix n :: nat assume "0 < n" |
|
1126 hence "0 < 2 * n" by auto |
|
1127 obtain t where "0 < t" and "t < real_of_float x" and |
|
1128 cos_eq: "cos x = (\<Sum> i = 0 ..< 2 * n. (if even(i) then ((- 1) ^ (i div 2))/((fact i)) else 0) * (real_of_float x) ^ i) |
|
1129 + (cos (t + 1/2 * (2 * n) * pi) / (fact (2*n))) * (real_of_float x)^(2*n)" |
|
1130 (is "_ = ?SUM + ?rest / ?fact * ?pow") |
|
1131 using Maclaurin_cos_expansion2[OF \<open>0 < real_of_float x\<close> \<open>0 < 2 * n\<close>] |
|
1132 unfolding cos_coeff_def atLeast0LessThan by auto |
|
1133 |
|
1134 have "cos t * (- 1) ^ n = cos t * cos (n * pi) + sin t * sin (n * pi)" by auto |
|
1135 also have "\<dots> = cos (t + n * pi)" by (simp add: cos_add) |
|
1136 also have "\<dots> = ?rest" by auto |
|
1137 finally have "cos t * (- 1) ^ n = ?rest" . |
|
1138 moreover |
|
1139 have "t \<le> pi / 2" using \<open>t < real_of_float x\<close> and \<open>x \<le> pi / 2\<close> by auto |
|
1140 hence "0 \<le> cos t" using \<open>0 < t\<close> and cos_ge_zero by auto |
|
1141 ultimately have even: "even n \<Longrightarrow> 0 \<le> ?rest" and odd: "odd n \<Longrightarrow> 0 \<le> - ?rest " by auto |
|
1142 |
|
1143 have "0 < ?fact" by auto |
|
1144 have "0 < ?pow" using \<open>0 < real_of_float x\<close> by auto |
|
1145 |
|
1146 { |
|
1147 assume "even n" |
|
1148 have "(lb_sin_cos_aux prec n 1 1 (x * x)) \<le> ?SUM" |
|
1149 unfolding morph_to_if_power[symmetric] using cos_aux by auto |
|
1150 also have "\<dots> \<le> cos x" |
|
1151 proof - |
|
1152 from even[OF \<open>even n\<close>] \<open>0 < ?fact\<close> \<open>0 < ?pow\<close> |
|
1153 have "0 \<le> (?rest / ?fact) * ?pow" by simp |
|
1154 thus ?thesis unfolding cos_eq by auto |
|
1155 qed |
|
1156 finally have "(lb_sin_cos_aux prec n 1 1 (x * x)) \<le> cos x" . |
|
1157 } note lb = this |
|
1158 |
|
1159 { |
|
1160 assume "odd n" |
|
1161 have "cos x \<le> ?SUM" |
|
1162 proof - |
|
1163 from \<open>0 < ?fact\<close> and \<open>0 < ?pow\<close> and odd[OF \<open>odd n\<close>] |
|
1164 have "0 \<le> (- ?rest) / ?fact * ?pow" |
|
1165 by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le) |
|
1166 thus ?thesis unfolding cos_eq by auto |
|
1167 qed |
|
1168 also have "\<dots> \<le> (ub_sin_cos_aux prec n 1 1 (x * x))" |
|
1169 unfolding morph_to_if_power[symmetric] using cos_aux by auto |
|
1170 finally have "cos x \<le> (ub_sin_cos_aux prec n 1 1 (x * x))" . |
|
1171 } note ub = this and lb |
|
1172 } note ub = this(1) and lb = this(2) |
|
1173 |
|
1174 have "cos x \<le> (ub_sin_cos_aux prec (get_odd n) 1 1 (x * x))" |
|
1175 using ub[OF odd_pos[OF get_odd] get_odd] . |
|
1176 moreover have "(lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) \<le> cos x" |
|
1177 proof (cases "0 < get_even n") |
|
1178 case True |
|
1179 show ?thesis using lb[OF True get_even] . |
|
1180 next |
|
1181 case False |
|
1182 hence "get_even n = 0" by auto |
|
1183 have "- (pi / 2) \<le> x" |
|
1184 by (rule order_trans[OF _ \<open>0 < real_of_float x\<close>[THEN less_imp_le]]) auto |
|
1185 with \<open>x \<le> pi / 2\<close> show ?thesis |
|
1186 unfolding \<open>get_even n = 0\<close> lb_sin_cos_aux.simps minus_float.rep_eq zero_float.rep_eq |
|
1187 using cos_ge_zero by auto |
|
1188 qed |
|
1189 ultimately show ?thesis by auto |
|
1190 next |
|
1191 case True |
|
1192 hence "x = 0" |
|
1193 by transfer |
|
1194 thus ?thesis |
|
1195 using lb_sin_cos_aux_zero_le_one one_le_ub_sin_cos_aux |
|
1196 by simp |
|
1197 qed |
|
1198 |
|
1199 lemma sin_aux: |
|
1200 assumes "0 \<le> real_of_float x" |
|
1201 shows "(x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> |
|
1202 (\<Sum> i=0..<n. (- 1) ^ i * (1/(fact (2 * i + 1))) * x^(2 * i + 1))" (is "?lb") |
|
1203 and "(\<Sum> i=0..<n. (- 1) ^ i * (1/(fact (2 * i + 1))) * x^(2 * i + 1)) \<le> |
|
1204 (x * ub_sin_cos_aux prec n 2 1 (x * x))" (is "?ub") |
|
1205 proof - |
|
1206 have "0 \<le> real_of_float (x * x)" by auto |
|
1207 let "?f n" = "fact (2 * n + 1) :: nat" |
|
1208 have f_eq: "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) ^^ n) 2 * (((\<lambda>i. i + 2) ^^ n) 2 + 1)" for n |
|
1209 proof - |
|
1210 have F: "\<And>m. ((\<lambda>i. i + 2) ^^ n) m = m + 2 * n" by (induct n) auto |
|
1211 show ?thesis |
|
1212 unfolding F by auto |
|
1213 qed |
|
1214 from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0, |
|
1215 OF \<open>0 \<le> real_of_float (x * x)\<close> f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps] |
|
1216 show "?lb" and "?ub" using \<open>0 \<le> real_of_float x\<close> |
|
1217 apply (simp_all only: power_add power_one_right mult.assoc[symmetric] sum_distrib_right[symmetric]) |
|
1218 apply (simp_all only: mult.commute[where 'a=real] of_nat_fact) |
|
1219 apply (auto intro!: mult_left_mono simp add: power_mult power2_eq_square[of "real_of_float x"]) |
|
1220 done |
|
1221 qed |
|
1222 |
|
1223 lemma sin_boundaries: |
|
1224 assumes "0 \<le> real_of_float x" |
|
1225 and "x \<le> pi / 2" |
|
1226 shows "sin x \<in> {(x * lb_sin_cos_aux prec (get_even n) 2 1 (x * x)) .. (x * ub_sin_cos_aux prec (get_odd n) 2 1 (x * x))}" |
|
1227 proof (cases "real_of_float x = 0") |
|
1228 case False |
|
1229 hence "real_of_float x \<noteq> 0" by auto |
|
1230 hence "0 < x" and "0 < real_of_float x" |
|
1231 using \<open>0 \<le> real_of_float x\<close> by auto |
|
1232 have "0 < x * x" |
|
1233 using \<open>0 < x\<close> by simp |
|
1234 |
|
1235 have sum_morph: "(\<Sum>j = 0 ..< n. (- 1) ^ (((2 * j + 1) - Suc 0) div 2) / ((fact (2 * j + 1))) * x ^(2 * j + 1)) = |
|
1236 (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else ((- 1) ^ ((i - Suc 0) div 2))/((fact i))) * x ^ i)" |
|
1237 (is "?SUM = _") for x :: real and n |
|
1238 proof - |
|
1239 have pow: "!!i. x ^ (2 * i + 1) = x * x ^ (2 * i)" |
|
1240 by auto |
|
1241 have "?SUM = (\<Sum> j = 0 ..< n. 0) + ?SUM" |
|
1242 by auto |
|
1243 also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. if even i then 0 else (- 1) ^ ((i - Suc 0) div 2) / ((fact i)) * x ^ i)" |
|
1244 unfolding sum_split_even_odd atLeast0LessThan .. |
|
1245 also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. (if even i then 0 else (- 1) ^ ((i - Suc 0) div 2) / ((fact i))) * x ^ i)" |
|
1246 by (rule sum.cong) auto |
|
1247 finally show ?thesis . |
|
1248 qed |
|
1249 |
|
1250 { fix n :: nat assume "0 < n" |
|
1251 hence "0 < 2 * n + 1" by auto |
|
1252 obtain t where "0 < t" and "t < real_of_float x" and |
|
1253 sin_eq: "sin x = (\<Sum> i = 0 ..< 2 * n + 1. (if even(i) then 0 else ((- 1) ^ ((i - Suc 0) div 2))/((fact i))) * (real_of_float x) ^ i) |
|
1254 + (sin (t + 1/2 * (2 * n + 1) * pi) / (fact (2*n + 1))) * (real_of_float x)^(2*n + 1)" |
|
1255 (is "_ = ?SUM + ?rest / ?fact * ?pow") |
|
1256 using Maclaurin_sin_expansion3[OF \<open>0 < 2 * n + 1\<close> \<open>0 < real_of_float x\<close>] |
|
1257 unfolding sin_coeff_def atLeast0LessThan by auto |
|
1258 |
|
1259 have "?rest = cos t * (- 1) ^ n" |
|
1260 unfolding sin_add cos_add of_nat_add distrib_right distrib_left by auto |
|
1261 moreover |
|
1262 have "t \<le> pi / 2" |
|
1263 using \<open>t < real_of_float x\<close> and \<open>x \<le> pi / 2\<close> by auto |
|
1264 hence "0 \<le> cos t" |
|
1265 using \<open>0 < t\<close> and cos_ge_zero by auto |
|
1266 ultimately have even: "even n \<Longrightarrow> 0 \<le> ?rest" and odd: "odd n \<Longrightarrow> 0 \<le> - ?rest" |
|
1267 by auto |
|
1268 |
|
1269 have "0 < ?fact" |
|
1270 by (simp del: fact_Suc) |
|
1271 have "0 < ?pow" |
|
1272 using \<open>0 < real_of_float x\<close> by (rule zero_less_power) |
|
1273 |
|
1274 { |
|
1275 assume "even n" |
|
1276 have "(x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> |
|
1277 (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else ((- 1) ^ ((i - Suc 0) div 2))/((fact i))) * (real_of_float x) ^ i)" |
|
1278 using sin_aux[OF \<open>0 \<le> real_of_float x\<close>] unfolding sum_morph[symmetric] by auto |
|
1279 also have "\<dots> \<le> ?SUM" by auto |
|
1280 also have "\<dots> \<le> sin x" |
|
1281 proof - |
|
1282 from even[OF \<open>even n\<close>] \<open>0 < ?fact\<close> \<open>0 < ?pow\<close> |
|
1283 have "0 \<le> (?rest / ?fact) * ?pow" by simp |
|
1284 thus ?thesis unfolding sin_eq by auto |
|
1285 qed |
|
1286 finally have "(x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> sin x" . |
|
1287 } note lb = this |
|
1288 |
|
1289 { |
|
1290 assume "odd n" |
|
1291 have "sin x \<le> ?SUM" |
|
1292 proof - |
|
1293 from \<open>0 < ?fact\<close> and \<open>0 < ?pow\<close> and odd[OF \<open>odd n\<close>] |
|
1294 have "0 \<le> (- ?rest) / ?fact * ?pow" |
|
1295 by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le) |
|
1296 thus ?thesis unfolding sin_eq by auto |
|
1297 qed |
|
1298 also have "\<dots> \<le> (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else ((- 1) ^ ((i - Suc 0) div 2))/((fact i))) * (real_of_float x) ^ i)" |
|
1299 by auto |
|
1300 also have "\<dots> \<le> (x * ub_sin_cos_aux prec n 2 1 (x * x))" |
|
1301 using sin_aux[OF \<open>0 \<le> real_of_float x\<close>] unfolding sum_morph[symmetric] by auto |
|
1302 finally have "sin x \<le> (x * ub_sin_cos_aux prec n 2 1 (x * x))" . |
|
1303 } note ub = this and lb |
|
1304 } note ub = this(1) and lb = this(2) |
|
1305 |
|
1306 have "sin x \<le> (x * ub_sin_cos_aux prec (get_odd n) 2 1 (x * x))" |
|
1307 using ub[OF odd_pos[OF get_odd] get_odd] . |
|
1308 moreover have "(x * lb_sin_cos_aux prec (get_even n) 2 1 (x * x)) \<le> sin x" |
|
1309 proof (cases "0 < get_even n") |
|
1310 case True |
|
1311 show ?thesis |
|
1312 using lb[OF True get_even] . |
|
1313 next |
|
1314 case False |
|
1315 hence "get_even n = 0" by auto |
|
1316 with \<open>x \<le> pi / 2\<close> \<open>0 \<le> real_of_float x\<close> |
|
1317 show ?thesis |
|
1318 unfolding \<open>get_even n = 0\<close> ub_sin_cos_aux.simps minus_float.rep_eq |
|
1319 using sin_ge_zero by auto |
|
1320 qed |
|
1321 ultimately show ?thesis by auto |
|
1322 next |
|
1323 case True |
|
1324 show ?thesis |
|
1325 proof (cases "n = 0") |
|
1326 case True |
|
1327 thus ?thesis |
|
1328 unfolding \<open>n = 0\<close> get_even_def get_odd_def |
|
1329 using \<open>real_of_float x = 0\<close> lapprox_rat[where x="-1" and y=1] by auto |
|
1330 next |
|
1331 case False |
|
1332 with not0_implies_Suc obtain m where "n = Suc m" by blast |
|
1333 thus ?thesis |
|
1334 unfolding \<open>n = Suc m\<close> get_even_def get_odd_def |
|
1335 using \<open>real_of_float x = 0\<close> rapprox_rat[where x=1 and y=1] lapprox_rat[where x=1 and y=1] |
|
1336 by (cases "even (Suc m)") auto |
|
1337 qed |
|
1338 qed |
|
1339 |
|
1340 |
|
1341 subsection "Compute the cosinus in the entire domain" |
|
1342 |
|
1343 definition lb_cos :: "nat \<Rightarrow> float \<Rightarrow> float" where |
|
1344 "lb_cos prec x = (let |
|
1345 horner = \<lambda> x. lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 1 1 (x * x) ; |
|
1346 half = \<lambda> x. if x < 0 then - 1 else float_plus_down prec (Float 1 1 * x * x) (- 1) |
|
1347 in if x < Float 1 (- 1) then horner x |
|
1348 else if x < 1 then half (horner (x * Float 1 (- 1))) |
|
1349 else half (half (horner (x * Float 1 (- 2)))))" |
|
1350 |
|
1351 definition ub_cos :: "nat \<Rightarrow> float \<Rightarrow> float" where |
|
1352 "ub_cos prec x = (let |
|
1353 horner = \<lambda> x. ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x) ; |
|
1354 half = \<lambda> x. float_plus_up prec (Float 1 1 * x * x) (- 1) |
|
1355 in if x < Float 1 (- 1) then horner x |
|
1356 else if x < 1 then half (horner (x * Float 1 (- 1))) |
|
1357 else half (half (horner (x * Float 1 (- 2)))))" |
|
1358 |
|
1359 lemma lb_cos: |
|
1360 assumes "0 \<le> real_of_float x" and "x \<le> pi" |
|
1361 shows "cos x \<in> {(lb_cos prec x) .. (ub_cos prec x)}" (is "?cos x \<in> {(?lb x) .. (?ub x) }") |
|
1362 proof - |
|
1363 have x_half[symmetric]: "cos x = 2 * cos (x / 2) * cos (x / 2) - 1" for x :: real |
|
1364 proof - |
|
1365 have "cos x = cos (x / 2 + x / 2)" |
|
1366 by auto |
|
1367 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" |
|
1368 unfolding cos_add by auto |
|
1369 also have "\<dots> = 2 * cos (x / 2) * cos (x / 2) - 1" |
|
1370 by algebra |
|
1371 finally show ?thesis . |
|
1372 qed |
|
1373 |
|
1374 have "\<not> x < 0" using \<open>0 \<le> real_of_float x\<close> by auto |
|
1375 let "?ub_horner x" = "ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x)" |
|
1376 let "?lb_horner x" = "lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 1 1 (x * x)" |
|
1377 let "?ub_half x" = "float_plus_up prec (Float 1 1 * x * x) (- 1)" |
|
1378 let "?lb_half x" = "if x < 0 then - 1 else float_plus_down prec (Float 1 1 * x * x) (- 1)" |
|
1379 |
|
1380 show ?thesis |
|
1381 proof (cases "x < Float 1 (- 1)") |
|
1382 case True |
|
1383 hence "x \<le> pi / 2" |
|
1384 using pi_ge_two by auto |
|
1385 show ?thesis |
|
1386 unfolding lb_cos_def[where x=x] ub_cos_def[where x=x] |
|
1387 if_not_P[OF \<open>\<not> x < 0\<close>] if_P[OF \<open>x < Float 1 (- 1)\<close>] Let_def |
|
1388 using cos_boundaries[OF \<open>0 \<le> real_of_float x\<close> \<open>x \<le> pi / 2\<close>] . |
|
1389 next |
|
1390 case False |
|
1391 { fix y x :: float let ?x2 = "(x * Float 1 (- 1))" |
|
1392 assume "y \<le> cos ?x2" and "-pi \<le> x" and "x \<le> pi" |
|
1393 hence "- (pi / 2) \<le> ?x2" and "?x2 \<le> pi / 2" |
|
1394 using pi_ge_two unfolding Float_num by auto |
|
1395 hence "0 \<le> cos ?x2" |
|
1396 by (rule cos_ge_zero) |
|
1397 |
|
1398 have "(?lb_half y) \<le> cos x" |
|
1399 proof (cases "y < 0") |
|
1400 case True |
|
1401 show ?thesis |
|
1402 using cos_ge_minus_one unfolding if_P[OF True] by auto |
|
1403 next |
|
1404 case False |
|
1405 hence "0 \<le> real_of_float y" by auto |
|
1406 from mult_mono[OF \<open>y \<le> cos ?x2\<close> \<open>y \<le> cos ?x2\<close> \<open>0 \<le> cos ?x2\<close> this] |
|
1407 have "real_of_float y * real_of_float y \<le> cos ?x2 * cos ?x2" . |
|
1408 hence "2 * real_of_float y * real_of_float y \<le> 2 * cos ?x2 * cos ?x2" |
|
1409 by auto |
|
1410 hence "2 * real_of_float y * real_of_float y - 1 \<le> 2 * cos (x / 2) * cos (x / 2) - 1" |
|
1411 unfolding Float_num by auto |
|
1412 thus ?thesis |
|
1413 unfolding if_not_P[OF False] x_half Float_num |
|
1414 by (auto intro!: float_plus_down_le) |
|
1415 qed |
|
1416 } note lb_half = this |
|
1417 |
|
1418 { fix y x :: float let ?x2 = "(x * Float 1 (- 1))" |
|
1419 assume ub: "cos ?x2 \<le> y" and "- pi \<le> x" and "x \<le> pi" |
|
1420 hence "- (pi / 2) \<le> ?x2" and "?x2 \<le> pi / 2" |
|
1421 using pi_ge_two unfolding Float_num by auto |
|
1422 hence "0 \<le> cos ?x2" by (rule cos_ge_zero) |
|
1423 |
|
1424 have "cos x \<le> (?ub_half y)" |
|
1425 proof - |
|
1426 have "0 \<le> real_of_float y" |
|
1427 using \<open>0 \<le> cos ?x2\<close> ub by (rule order_trans) |
|
1428 from mult_mono[OF ub ub this \<open>0 \<le> cos ?x2\<close>] |
|
1429 have "cos ?x2 * cos ?x2 \<le> real_of_float y * real_of_float y" . |
|
1430 hence "2 * cos ?x2 * cos ?x2 \<le> 2 * real_of_float y * real_of_float y" |
|
1431 by auto |
|
1432 hence "2 * cos (x / 2) * cos (x / 2) - 1 \<le> 2 * real_of_float y * real_of_float y - 1" |
|
1433 unfolding Float_num by auto |
|
1434 thus ?thesis |
|
1435 unfolding x_half Float_num |
|
1436 by (auto intro!: float_plus_up_le) |
|
1437 qed |
|
1438 } note ub_half = this |
|
1439 |
|
1440 let ?x2 = "x * Float 1 (- 1)" |
|
1441 let ?x4 = "x * Float 1 (- 1) * Float 1 (- 1)" |
|
1442 |
|
1443 have "-pi \<le> x" |
|
1444 using pi_ge_zero[THEN le_imp_neg_le, unfolded minus_zero] \<open>0 \<le> real_of_float x\<close> |
|
1445 by (rule order_trans) |
|
1446 |
|
1447 show ?thesis |
|
1448 proof (cases "x < 1") |
|
1449 case True |
|
1450 hence "real_of_float x \<le> 1" by auto |
|
1451 have "0 \<le> real_of_float ?x2" and "?x2 \<le> pi / 2" |
|
1452 using pi_ge_two \<open>0 \<le> real_of_float x\<close> using assms by auto |
|
1453 from cos_boundaries[OF this] |
|
1454 have lb: "(?lb_horner ?x2) \<le> ?cos ?x2" and ub: "?cos ?x2 \<le> (?ub_horner ?x2)" |
|
1455 by auto |
|
1456 |
|
1457 have "(?lb x) \<le> ?cos x" |
|
1458 proof - |
|
1459 from lb_half[OF lb \<open>-pi \<le> x\<close> \<open>x \<le> pi\<close>] |
|
1460 show ?thesis |
|
1461 unfolding lb_cos_def[where x=x] Let_def |
|
1462 using \<open>\<not> x < 0\<close> \<open>\<not> x < Float 1 (- 1)\<close> \<open>x < 1\<close> by auto |
|
1463 qed |
|
1464 moreover have "?cos x \<le> (?ub x)" |
|
1465 proof - |
|
1466 from ub_half[OF ub \<open>-pi \<le> x\<close> \<open>x \<le> pi\<close>] |
|
1467 show ?thesis |
|
1468 unfolding ub_cos_def[where x=x] Let_def |
|
1469 using \<open>\<not> x < 0\<close> \<open>\<not> x < Float 1 (- 1)\<close> \<open>x < 1\<close> by auto |
|
1470 qed |
|
1471 ultimately show ?thesis by auto |
|
1472 next |
|
1473 case False |
|
1474 have "0 \<le> real_of_float ?x4" and "?x4 \<le> pi / 2" |
|
1475 using pi_ge_two \<open>0 \<le> real_of_float x\<close> \<open>x \<le> pi\<close> unfolding Float_num by auto |
|
1476 from cos_boundaries[OF this] |
|
1477 have lb: "(?lb_horner ?x4) \<le> ?cos ?x4" and ub: "?cos ?x4 \<le> (?ub_horner ?x4)" |
|
1478 by auto |
|
1479 |
|
1480 have eq_4: "?x2 * Float 1 (- 1) = x * Float 1 (- 2)" |
|
1481 by transfer simp |
|
1482 |
|
1483 have "(?lb x) \<le> ?cos x" |
|
1484 proof - |
|
1485 have "-pi \<le> ?x2" and "?x2 \<le> pi" |
|
1486 using pi_ge_two \<open>0 \<le> real_of_float x\<close> \<open>x \<le> pi\<close> by auto |
|
1487 from lb_half[OF lb_half[OF lb this] \<open>-pi \<le> x\<close> \<open>x \<le> pi\<close>, unfolded eq_4] |
|
1488 show ?thesis |
|
1489 unfolding lb_cos_def[where x=x] if_not_P[OF \<open>\<not> x < 0\<close>] |
|
1490 if_not_P[OF \<open>\<not> x < Float 1 (- 1)\<close>] if_not_P[OF \<open>\<not> x < 1\<close>] Let_def . |
|
1491 qed |
|
1492 moreover have "?cos x \<le> (?ub x)" |
|
1493 proof - |
|
1494 have "-pi \<le> ?x2" and "?x2 \<le> pi" |
|
1495 using pi_ge_two \<open>0 \<le> real_of_float x\<close> \<open> x \<le> pi\<close> by auto |
|
1496 from ub_half[OF ub_half[OF ub this] \<open>-pi \<le> x\<close> \<open>x \<le> pi\<close>, unfolded eq_4] |
|
1497 show ?thesis |
|
1498 unfolding ub_cos_def[where x=x] if_not_P[OF \<open>\<not> x < 0\<close>] |
|
1499 if_not_P[OF \<open>\<not> x < Float 1 (- 1)\<close>] if_not_P[OF \<open>\<not> x < 1\<close>] Let_def . |
|
1500 qed |
|
1501 ultimately show ?thesis by auto |
|
1502 qed |
|
1503 qed |
|
1504 qed |
|
1505 |
|
1506 lemma lb_cos_minus: |
|
1507 assumes "-pi \<le> x" |
|
1508 and "real_of_float x \<le> 0" |
|
1509 shows "cos (real_of_float(-x)) \<in> {(lb_cos prec (-x)) .. (ub_cos prec (-x))}" |
|
1510 proof - |
|
1511 have "0 \<le> real_of_float (-x)" and "(-x) \<le> pi" |
|
1512 using \<open>-pi \<le> x\<close> \<open>real_of_float x \<le> 0\<close> by auto |
|
1513 from lb_cos[OF this] show ?thesis . |
|
1514 qed |
|
1515 |
|
1516 definition bnds_cos :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where |
|
1517 "bnds_cos prec lx ux = (let |
|
1518 lpi = float_round_down prec (lb_pi prec) ; |
|
1519 upi = float_round_up prec (ub_pi prec) ; |
|
1520 k = floor_fl (float_divr prec (lx + lpi) (2 * lpi)) ; |
|
1521 lx = float_plus_down prec lx (- k * 2 * (if k < 0 then lpi else upi)) ; |
|
1522 ux = float_plus_up prec ux (- k * 2 * (if k < 0 then upi else lpi)) |
|
1523 in if - lpi \<le> lx \<and> ux \<le> 0 then (lb_cos prec (-lx), ub_cos prec (-ux)) |
|
1524 else if 0 \<le> lx \<and> ux \<le> lpi then (lb_cos prec ux, ub_cos prec lx) |
|
1525 else if - lpi \<le> lx \<and> ux \<le> lpi then (min (lb_cos prec (-lx)) (lb_cos prec ux), Float 1 0) |
|
1526 else if 0 \<le> lx \<and> ux \<le> 2 * lpi then (Float (- 1) 0, max (ub_cos prec lx) (ub_cos prec (- (ux - 2 * lpi)))) |
|
1527 else if -2 * lpi \<le> lx \<and> ux \<le> 0 then (Float (- 1) 0, max (ub_cos prec (lx + 2 * lpi)) (ub_cos prec (-ux))) |
|
1528 else (Float (- 1) 0, Float 1 0))" |
|
1529 |
|
1530 lemma floor_int: obtains k :: int where "real_of_int k = (floor_fl f)" |
|
1531 by (simp add: floor_fl_def) |
|
1532 |
|
1533 lemma cos_periodic_nat[simp]: |
|
1534 fixes n :: nat |
|
1535 shows "cos (x + n * (2 * pi)) = cos x" |
|
1536 proof (induct n arbitrary: x) |
|
1537 case 0 |
|
1538 then show ?case by simp |
|
1539 next |
|
1540 case (Suc n) |
|
1541 have split_pi_off: "x + (Suc n) * (2 * pi) = (x + n * (2 * pi)) + 2 * pi" |
|
1542 unfolding Suc_eq_plus1 of_nat_add of_int_1 distrib_right by auto |
|
1543 show ?case |
|
1544 unfolding split_pi_off using Suc by auto |
|
1545 qed |
|
1546 |
|
1547 lemma cos_periodic_int[simp]: |
|
1548 fixes i :: int |
|
1549 shows "cos (x + i * (2 * pi)) = cos x" |
|
1550 proof (cases "0 \<le> i") |
|
1551 case True |
|
1552 hence i_nat: "real_of_int i = nat i" by auto |
|
1553 show ?thesis |
|
1554 unfolding i_nat by auto |
|
1555 next |
|
1556 case False |
|
1557 hence i_nat: "i = - real (nat (-i))" by auto |
|
1558 have "cos x = cos (x + i * (2 * pi) - i * (2 * pi))" |
|
1559 by auto |
|
1560 also have "\<dots> = cos (x + i * (2 * pi))" |
|
1561 unfolding i_nat mult_minus_left diff_minus_eq_add by (rule cos_periodic_nat) |
|
1562 finally show ?thesis by auto |
|
1563 qed |
|
1564 |
|
1565 lemma bnds_cos: "\<forall>(x::real) lx ux. (l, u) = |
|
1566 bnds_cos prec lx ux \<and> x \<in> {lx .. ux} \<longrightarrow> l \<le> cos x \<and> cos x \<le> u" |
|
1567 proof (rule allI | rule impI | erule conjE)+ |
|
1568 fix x :: real |
|
1569 fix lx ux |
|
1570 assume bnds: "(l, u) = bnds_cos prec lx ux" and x: "x \<in> {lx .. ux}" |
|
1571 |
|
1572 let ?lpi = "float_round_down prec (lb_pi prec)" |
|
1573 let ?upi = "float_round_up prec (ub_pi prec)" |
|
1574 let ?k = "floor_fl (float_divr prec (lx + ?lpi) (2 * ?lpi))" |
|
1575 let ?lx2 = "(- ?k * 2 * (if ?k < 0 then ?lpi else ?upi))" |
|
1576 let ?ux2 = "(- ?k * 2 * (if ?k < 0 then ?upi else ?lpi))" |
|
1577 let ?lx = "float_plus_down prec lx ?lx2" |
|
1578 let ?ux = "float_plus_up prec ux ?ux2" |
|
1579 |
|
1580 obtain k :: int where k: "k = real_of_float ?k" |
|
1581 by (rule floor_int) |
|
1582 |
|
1583 have upi: "pi \<le> ?upi" and lpi: "?lpi \<le> pi" |
|
1584 using float_round_up[of "ub_pi prec" prec] pi_boundaries[of prec] |
|
1585 float_round_down[of prec "lb_pi prec"] |
|
1586 by auto |
|
1587 hence "lx + ?lx2 \<le> x - k * (2 * pi) \<and> x - k * (2 * pi) \<le> ux + ?ux2" |
|
1588 using x |
|
1589 by (cases "k = 0") |
|
1590 (auto intro!: add_mono |
|
1591 simp add: k [symmetric] uminus_add_conv_diff [symmetric] |
|
1592 simp del: float_of_numeral uminus_add_conv_diff) |
|
1593 hence "?lx \<le> x - k * (2 * pi) \<and> x - k * (2 * pi) \<le> ?ux" |
|
1594 by (auto intro!: float_plus_down_le float_plus_up_le) |
|
1595 note lx = this[THEN conjunct1] and ux = this[THEN conjunct2] |
|
1596 hence lx_less_ux: "?lx \<le> real_of_float ?ux" by (rule order_trans) |
|
1597 |
|
1598 { assume "- ?lpi \<le> ?lx" and x_le_0: "x - k * (2 * pi) \<le> 0" |
|
1599 with lpi[THEN le_imp_neg_le] lx |
|
1600 have pi_lx: "- pi \<le> ?lx" and lx_0: "real_of_float ?lx \<le> 0" |
|
1601 by simp_all |
|
1602 |
|
1603 have "(lb_cos prec (- ?lx)) \<le> cos (real_of_float (- ?lx))" |
|
1604 using lb_cos_minus[OF pi_lx lx_0] by simp |
|
1605 also have "\<dots> \<le> cos (x + (-k) * (2 * pi))" |
|
1606 using cos_monotone_minus_pi_0'[OF pi_lx lx x_le_0] |
|
1607 by (simp only: uminus_float.rep_eq of_int_minus |
|
1608 cos_minus mult_minus_left) simp |
|
1609 finally have "(lb_cos prec (- ?lx)) \<le> cos x" |
|
1610 unfolding cos_periodic_int . } |
|
1611 note negative_lx = this |
|
1612 |
|
1613 { assume "0 \<le> ?lx" and pi_x: "x - k * (2 * pi) \<le> pi" |
|
1614 with lx |
|
1615 have pi_lx: "?lx \<le> pi" and lx_0: "0 \<le> real_of_float ?lx" |
|
1616 by auto |
|
1617 |
|
1618 have "cos (x + (-k) * (2 * pi)) \<le> cos ?lx" |
|
1619 using cos_monotone_0_pi_le[OF lx_0 lx pi_x] |
|
1620 by (simp only: of_int_minus |
|
1621 cos_minus mult_minus_left) simp |
|
1622 also have "\<dots> \<le> (ub_cos prec ?lx)" |
|
1623 using lb_cos[OF lx_0 pi_lx] by simp |
|
1624 finally have "cos x \<le> (ub_cos prec ?lx)" |
|
1625 unfolding cos_periodic_int . } |
|
1626 note positive_lx = this |
|
1627 |
|
1628 { assume pi_x: "- pi \<le> x - k * (2 * pi)" and "?ux \<le> 0" |
|
1629 with ux |
|
1630 have pi_ux: "- pi \<le> ?ux" and ux_0: "real_of_float ?ux \<le> 0" |
|
1631 by simp_all |
|
1632 |
|
1633 have "cos (x + (-k) * (2 * pi)) \<le> cos (real_of_float (- ?ux))" |
|
1634 using cos_monotone_minus_pi_0'[OF pi_x ux ux_0] |
|
1635 by (simp only: uminus_float.rep_eq of_int_minus |
|
1636 cos_minus mult_minus_left) simp |
|
1637 also have "\<dots> \<le> (ub_cos prec (- ?ux))" |
|
1638 using lb_cos_minus[OF pi_ux ux_0, of prec] by simp |
|
1639 finally have "cos x \<le> (ub_cos prec (- ?ux))" |
|
1640 unfolding cos_periodic_int . } |
|
1641 note negative_ux = this |
|
1642 |
|
1643 { assume "?ux \<le> ?lpi" and x_ge_0: "0 \<le> x - k * (2 * pi)" |
|
1644 with lpi ux |
|
1645 have pi_ux: "?ux \<le> pi" and ux_0: "0 \<le> real_of_float ?ux" |
|
1646 by simp_all |
|
1647 |
|
1648 have "(lb_cos prec ?ux) \<le> cos ?ux" |
|
1649 using lb_cos[OF ux_0 pi_ux] by simp |
|
1650 also have "\<dots> \<le> cos (x + (-k) * (2 * pi))" |
|
1651 using cos_monotone_0_pi_le[OF x_ge_0 ux pi_ux] |
|
1652 by (simp only: of_int_minus |
|
1653 cos_minus mult_minus_left) simp |
|
1654 finally have "(lb_cos prec ?ux) \<le> cos x" |
|
1655 unfolding cos_periodic_int . } |
|
1656 note positive_ux = this |
|
1657 |
|
1658 show "l \<le> cos x \<and> cos x \<le> u" |
|
1659 proof (cases "- ?lpi \<le> ?lx \<and> ?ux \<le> 0") |
|
1660 case True |
|
1661 with bnds have l: "l = lb_cos prec (-?lx)" and u: "u = ub_cos prec (-?ux)" |
|
1662 by (auto simp add: bnds_cos_def Let_def) |
|
1663 from True lpi[THEN le_imp_neg_le] lx ux |
|
1664 have "- pi \<le> x - k * (2 * pi)" and "x - k * (2 * pi) \<le> 0" |
|
1665 by auto |
|
1666 with True negative_ux negative_lx show ?thesis |
|
1667 unfolding l u by simp |
|
1668 next |
|
1669 case 1: False |
|
1670 show ?thesis |
|
1671 proof (cases "0 \<le> ?lx \<and> ?ux \<le> ?lpi") |
|
1672 case True with bnds 1 |
|
1673 have l: "l = lb_cos prec ?ux" |
|
1674 and u: "u = ub_cos prec ?lx" |
|
1675 by (auto simp add: bnds_cos_def Let_def) |
|
1676 from True lpi lx ux |
|
1677 have "0 \<le> x - k * (2 * pi)" and "x - k * (2 * pi) \<le> pi" |
|
1678 by auto |
|
1679 with True positive_ux positive_lx show ?thesis |
|
1680 unfolding l u by simp |
|
1681 next |
|
1682 case 2: False |
|
1683 show ?thesis |
|
1684 proof (cases "- ?lpi \<le> ?lx \<and> ?ux \<le> ?lpi") |
|
1685 case Cond: True |
|
1686 with bnds 1 2 have l: "l = min (lb_cos prec (-?lx)) (lb_cos prec ?ux)" |
|
1687 and u: "u = Float 1 0" |
|
1688 by (auto simp add: bnds_cos_def Let_def) |
|
1689 show ?thesis |
|
1690 unfolding u l using negative_lx positive_ux Cond |
|
1691 by (cases "x - k * (2 * pi) < 0") (auto simp add: real_of_float_min) |
|
1692 next |
|
1693 case 3: False |
|
1694 show ?thesis |
|
1695 proof (cases "0 \<le> ?lx \<and> ?ux \<le> 2 * ?lpi") |
|
1696 case Cond: True |
|
1697 with bnds 1 2 3 |
|
1698 have l: "l = Float (- 1) 0" |
|
1699 and u: "u = max (ub_cos prec ?lx) (ub_cos prec (- (?ux - 2 * ?lpi)))" |
|
1700 by (auto simp add: bnds_cos_def Let_def) |
|
1701 |
|
1702 have "cos x \<le> real_of_float u" |
|
1703 proof (cases "x - k * (2 * pi) < pi") |
|
1704 case True |
|
1705 hence "x - k * (2 * pi) \<le> pi" by simp |
|
1706 from positive_lx[OF Cond[THEN conjunct1] this] show ?thesis |
|
1707 unfolding u by (simp add: real_of_float_max) |
|
1708 next |
|
1709 case False |
|
1710 hence "pi \<le> x - k * (2 * pi)" by simp |
|
1711 hence pi_x: "- pi \<le> x - k * (2 * pi) - 2 * pi" by simp |
|
1712 |
|
1713 have "?ux \<le> 2 * pi" |
|
1714 using Cond lpi by auto |
|
1715 hence "x - k * (2 * pi) - 2 * pi \<le> 0" |
|
1716 using ux by simp |
|
1717 |
|
1718 have ux_0: "real_of_float (?ux - 2 * ?lpi) \<le> 0" |
|
1719 using Cond by auto |
|
1720 |
|
1721 from 2 and Cond have "\<not> ?ux \<le> ?lpi" by auto |
|
1722 hence "- ?lpi \<le> ?ux - 2 * ?lpi" by auto |
|
1723 hence pi_ux: "- pi \<le> (?ux - 2 * ?lpi)" |
|
1724 using lpi[THEN le_imp_neg_le] by auto |
|
1725 |
|
1726 have x_le_ux: "x - k * (2 * pi) - 2 * pi \<le> (?ux - 2 * ?lpi)" |
|
1727 using ux lpi by auto |
|
1728 have "cos x = cos (x + (-k) * (2 * pi) + (-1::int) * (2 * pi))" |
|
1729 unfolding cos_periodic_int .. |
|
1730 also have "\<dots> \<le> cos ((?ux - 2 * ?lpi))" |
|
1731 using cos_monotone_minus_pi_0'[OF pi_x x_le_ux ux_0] |
|
1732 by (simp only: minus_float.rep_eq of_int_minus of_int_1 |
|
1733 mult_minus_left mult_1_left) simp |
|
1734 also have "\<dots> = cos ((- (?ux - 2 * ?lpi)))" |
|
1735 unfolding uminus_float.rep_eq cos_minus .. |
|
1736 also have "\<dots> \<le> (ub_cos prec (- (?ux - 2 * ?lpi)))" |
|
1737 using lb_cos_minus[OF pi_ux ux_0] by simp |
|
1738 finally show ?thesis unfolding u by (simp add: real_of_float_max) |
|
1739 qed |
|
1740 thus ?thesis unfolding l by auto |
|
1741 next |
|
1742 case 4: False |
|
1743 show ?thesis |
|
1744 proof (cases "-2 * ?lpi \<le> ?lx \<and> ?ux \<le> 0") |
|
1745 case Cond: True |
|
1746 with bnds 1 2 3 4 have l: "l = Float (- 1) 0" |
|
1747 and u: "u = max (ub_cos prec (?lx + 2 * ?lpi)) (ub_cos prec (-?ux))" |
|
1748 by (auto simp add: bnds_cos_def Let_def) |
|
1749 |
|
1750 have "cos x \<le> u" |
|
1751 proof (cases "-pi < x - k * (2 * pi)") |
|
1752 case True |
|
1753 hence "-pi \<le> x - k * (2 * pi)" by simp |
|
1754 from negative_ux[OF this Cond[THEN conjunct2]] show ?thesis |
|
1755 unfolding u by (simp add: real_of_float_max) |
|
1756 next |
|
1757 case False |
|
1758 hence "x - k * (2 * pi) \<le> -pi" by simp |
|
1759 hence pi_x: "x - k * (2 * pi) + 2 * pi \<le> pi" by simp |
|
1760 |
|
1761 have "-2 * pi \<le> ?lx" using Cond lpi by auto |
|
1762 |
|
1763 hence "0 \<le> x - k * (2 * pi) + 2 * pi" using lx by simp |
|
1764 |
|
1765 have lx_0: "0 \<le> real_of_float (?lx + 2 * ?lpi)" |
|
1766 using Cond lpi by auto |
|
1767 |
|
1768 from 1 and Cond have "\<not> -?lpi \<le> ?lx" by auto |
|
1769 hence "?lx + 2 * ?lpi \<le> ?lpi" by auto |
|
1770 hence pi_lx: "(?lx + 2 * ?lpi) \<le> pi" |
|
1771 using lpi[THEN le_imp_neg_le] by auto |
|
1772 |
|
1773 have lx_le_x: "(?lx + 2 * ?lpi) \<le> x - k * (2 * pi) + 2 * pi" |
|
1774 using lx lpi by auto |
|
1775 |
|
1776 have "cos x = cos (x + (-k) * (2 * pi) + (1 :: int) * (2 * pi))" |
|
1777 unfolding cos_periodic_int .. |
|
1778 also have "\<dots> \<le> cos ((?lx + 2 * ?lpi))" |
|
1779 using cos_monotone_0_pi_le[OF lx_0 lx_le_x pi_x] |
|
1780 by (simp only: minus_float.rep_eq of_int_minus of_int_1 |
|
1781 mult_minus_left mult_1_left) simp |
|
1782 also have "\<dots> \<le> (ub_cos prec (?lx + 2 * ?lpi))" |
|
1783 using lb_cos[OF lx_0 pi_lx] by simp |
|
1784 finally show ?thesis unfolding u by (simp add: real_of_float_max) |
|
1785 qed |
|
1786 thus ?thesis unfolding l by auto |
|
1787 next |
|
1788 case False |
|
1789 with bnds 1 2 3 4 show ?thesis |
|
1790 by (auto simp add: bnds_cos_def Let_def) |
|
1791 qed |
|
1792 qed |
|
1793 qed |
|
1794 qed |
|
1795 qed |
|
1796 qed |
|
1797 |
|
1798 |
|
1799 section "Exponential function" |
|
1800 |
|
1801 subsection "Compute the series of the exponential function" |
|
1802 |
|
1803 fun ub_exp_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" |
|
1804 and lb_exp_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" |
|
1805 where |
|
1806 "ub_exp_horner prec 0 i k x = 0" | |
|
1807 "ub_exp_horner prec (Suc n) i k x = float_plus_up prec |
|
1808 (rapprox_rat prec 1 (int k)) (float_round_up prec (x * lb_exp_horner prec n (i + 1) (k * i) x))" | |
|
1809 "lb_exp_horner prec 0 i k x = 0" | |
|
1810 "lb_exp_horner prec (Suc n) i k x = float_plus_down prec |
|
1811 (lapprox_rat prec 1 (int k)) (float_round_down prec (x * ub_exp_horner prec n (i + 1) (k * i) x))" |
|
1812 |
|
1813 lemma bnds_exp_horner: |
|
1814 assumes "real_of_float x \<le> 0" |
|
1815 shows "exp x \<in> {lb_exp_horner prec (get_even n) 1 1 x .. ub_exp_horner prec (get_odd n) 1 1 x}" |
|
1816 proof - |
|
1817 have f_eq: "fact (Suc n) = fact n * ((\<lambda>i::nat. i + 1) ^^ n) 1" for n |
|
1818 proof - |
|
1819 have F: "\<And> m. ((\<lambda>i. i + 1) ^^ n) m = n + m" |
|
1820 by (induct n) auto |
|
1821 show ?thesis |
|
1822 unfolding F by auto |
|
1823 qed |
|
1824 |
|
1825 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, |
|
1826 OF assms f_eq lb_exp_horner.simps ub_exp_horner.simps] |
|
1827 |
|
1828 have "lb_exp_horner prec (get_even n) 1 1 x \<le> exp x" |
|
1829 proof - |
|
1830 have "lb_exp_horner prec (get_even n) 1 1 x \<le> (\<Sum>j = 0..<get_even n. 1 / (fact j) * real_of_float x ^ j)" |
|
1831 using bounds(1) by auto |
|
1832 also have "\<dots> \<le> exp x" |
|
1833 proof - |
|
1834 obtain t where "\<bar>t\<bar> \<le> \<bar>real_of_float x\<bar>" and "exp x = (\<Sum>m = 0..<get_even n. real_of_float x ^ m / (fact m)) + exp t / (fact (get_even n)) * (real_of_float x) ^ (get_even n)" |
|
1835 using Maclaurin_exp_le unfolding atLeast0LessThan by blast |
|
1836 moreover have "0 \<le> exp t / (fact (get_even n)) * (real_of_float x) ^ (get_even n)" |
|
1837 by (auto simp: zero_le_even_power) |
|
1838 ultimately show ?thesis using get_odd exp_gt_zero by auto |
|
1839 qed |
|
1840 finally show ?thesis . |
|
1841 qed |
|
1842 moreover |
|
1843 have "exp x \<le> ub_exp_horner prec (get_odd n) 1 1 x" |
|
1844 proof - |
|
1845 have x_less_zero: "real_of_float x ^ get_odd n \<le> 0" |
|
1846 proof (cases "real_of_float x = 0") |
|
1847 case True |
|
1848 have "(get_odd n) \<noteq> 0" using get_odd[THEN odd_pos] by auto |
|
1849 thus ?thesis unfolding True power_0_left by auto |
|
1850 next |
|
1851 case False hence "real_of_float x < 0" using \<open>real_of_float x \<le> 0\<close> by auto |
|
1852 show ?thesis by (rule less_imp_le, auto simp add: \<open>real_of_float x < 0\<close>) |
|
1853 qed |
|
1854 obtain t where "\<bar>t\<bar> \<le> \<bar>real_of_float x\<bar>" |
|
1855 and "exp x = (\<Sum>m = 0..<get_odd n. (real_of_float x) ^ m / (fact m)) + exp t / (fact (get_odd n)) * (real_of_float x) ^ (get_odd n)" |
|
1856 using Maclaurin_exp_le unfolding atLeast0LessThan by blast |
|
1857 moreover have "exp t / (fact (get_odd n)) * (real_of_float x) ^ (get_odd n) \<le> 0" |
|
1858 by (auto intro!: mult_nonneg_nonpos divide_nonpos_pos simp add: x_less_zero) |
|
1859 ultimately have "exp x \<le> (\<Sum>j = 0..<get_odd n. 1 / (fact j) * real_of_float x ^ j)" |
|
1860 using get_odd exp_gt_zero by auto |
|
1861 also have "\<dots> \<le> ub_exp_horner prec (get_odd n) 1 1 x" |
|
1862 using bounds(2) by auto |
|
1863 finally show ?thesis . |
|
1864 qed |
|
1865 ultimately show ?thesis by auto |
|
1866 qed |
|
1867 |
|
1868 lemma ub_exp_horner_nonneg: "real_of_float x \<le> 0 \<Longrightarrow> |
|
1869 0 \<le> real_of_float (ub_exp_horner prec (get_odd n) (Suc 0) (Suc 0) x)" |
|
1870 using bnds_exp_horner[of x prec n] |
|
1871 by (intro order_trans[OF exp_ge_zero]) auto |
|
1872 |
|
1873 |
|
1874 subsection "Compute the exponential function on the entire domain" |
|
1875 |
|
1876 function ub_exp :: "nat \<Rightarrow> float \<Rightarrow> float" and lb_exp :: "nat \<Rightarrow> float \<Rightarrow> float" where |
|
1877 "lb_exp prec x = |
|
1878 (if 0 < x then float_divl prec 1 (ub_exp prec (-x)) |
|
1879 else |
|
1880 let |
|
1881 horner = (\<lambda> x. let y = lb_exp_horner prec (get_even (prec + 2)) 1 1 x in |
|
1882 if y \<le> 0 then Float 1 (- 2) else y) |
|
1883 in |
|
1884 if x < - 1 then |
|
1885 power_down_fl prec (horner (float_divl prec x (- floor_fl x))) (nat (- int_floor_fl x)) |
|
1886 else horner x)" | |
|
1887 "ub_exp prec x = |
|
1888 (if 0 < x then float_divr prec 1 (lb_exp prec (-x)) |
|
1889 else if x < - 1 then |
|
1890 power_up_fl prec |
|
1891 (ub_exp_horner prec (get_odd (prec + 2)) 1 1 |
|
1892 (float_divr prec x (- floor_fl x))) (nat (- int_floor_fl x)) |
|
1893 else ub_exp_horner prec (get_odd (prec + 2)) 1 1 x)" |
|
1894 by pat_completeness auto |
|
1895 termination |
|
1896 by (relation "measure (\<lambda> v. let (prec, x) = case_sum id id v in (if 0 < x then 1 else 0))") auto |
|
1897 |
|
1898 lemma exp_m1_ge_quarter: "(1 / 4 :: real) \<le> exp (- 1)" |
|
1899 proof - |
|
1900 have eq4: "4 = Suc (Suc (Suc (Suc 0)))" by auto |
|
1901 have "1 / 4 = (Float 1 (- 2))" |
|
1902 unfolding Float_num by auto |
|
1903 also have "\<dots> \<le> lb_exp_horner 3 (get_even 3) 1 1 (- 1)" |
|
1904 by (subst less_eq_float.rep_eq [symmetric]) code_simp |
|
1905 also have "\<dots> \<le> exp (- 1 :: float)" |
|
1906 using bnds_exp_horner[where x="- 1"] by auto |
|
1907 finally show ?thesis |
|
1908 by simp |
|
1909 qed |
|
1910 |
|
1911 lemma lb_exp_pos: |
|
1912 assumes "\<not> 0 < x" |
|
1913 shows "0 < lb_exp prec x" |
|
1914 proof - |
|
1915 let "?lb_horner x" = "lb_exp_horner prec (get_even (prec + 2)) 1 1 x" |
|
1916 let "?horner x" = "let y = ?lb_horner x in if y \<le> 0 then Float 1 (- 2) else y" |
|
1917 have pos_horner: "0 < ?horner x" for x |
|
1918 unfolding Let_def by (cases "?lb_horner x \<le> 0") auto |
|
1919 moreover have "0 < real_of_float ((?horner x) ^ num)" for x :: float and num :: nat |
|
1920 proof - |
|
1921 have "0 < real_of_float (?horner x) ^ num" using \<open>0 < ?horner x\<close> by simp |
|
1922 also have "\<dots> = (?horner x) ^ num" by auto |
|
1923 finally show ?thesis . |
|
1924 qed |
|
1925 ultimately show ?thesis |
|
1926 unfolding lb_exp.simps if_not_P[OF \<open>\<not> 0 < x\<close>] Let_def |
|
1927 by (cases "floor_fl x", cases "x < - 1") |
|
1928 (auto simp: real_power_up_fl real_power_down_fl intro!: power_up_less power_down_pos) |
|
1929 qed |
|
1930 |
|
1931 lemma exp_boundaries': |
|
1932 assumes "x \<le> 0" |
|
1933 shows "exp x \<in> { (lb_exp prec x) .. (ub_exp prec x)}" |
|
1934 proof - |
|
1935 let "?lb_exp_horner x" = "lb_exp_horner prec (get_even (prec + 2)) 1 1 x" |
|
1936 let "?ub_exp_horner x" = "ub_exp_horner prec (get_odd (prec + 2)) 1 1 x" |
|
1937 |
|
1938 have "real_of_float x \<le> 0" and "\<not> x > 0" |
|
1939 using \<open>x \<le> 0\<close> by auto |
|
1940 show ?thesis |
|
1941 proof (cases "x < - 1") |
|
1942 case False |
|
1943 hence "- 1 \<le> real_of_float x" by auto |
|
1944 show ?thesis |
|
1945 proof (cases "?lb_exp_horner x \<le> 0") |
|
1946 case True |
|
1947 from \<open>\<not> x < - 1\<close> |
|
1948 have "- 1 \<le> real_of_float x" by auto |
|
1949 hence "exp (- 1) \<le> exp x" |
|
1950 unfolding exp_le_cancel_iff . |
|
1951 from order_trans[OF exp_m1_ge_quarter this] have "Float 1 (- 2) \<le> exp x" |
|
1952 unfolding Float_num . |
|
1953 with True show ?thesis |
|
1954 using bnds_exp_horner \<open>real_of_float x \<le> 0\<close> \<open>\<not> x > 0\<close> \<open>\<not> x < - 1\<close> by auto |
|
1955 next |
|
1956 case False |
|
1957 thus ?thesis |
|
1958 using bnds_exp_horner \<open>real_of_float x \<le> 0\<close> \<open>\<not> x > 0\<close> \<open>\<not> x < - 1\<close> by (auto simp add: Let_def) |
|
1959 qed |
|
1960 next |
|
1961 case True |
|
1962 let ?num = "nat (- int_floor_fl x)" |
|
1963 |
|
1964 have "real_of_int (int_floor_fl x) < - 1" |
|
1965 using int_floor_fl[of x] \<open>x < - 1\<close> by simp |
|
1966 hence "real_of_int (int_floor_fl x) < 0" by simp |
|
1967 hence "int_floor_fl x < 0" by auto |
|
1968 hence "1 \<le> - int_floor_fl x" by auto |
|
1969 hence "0 < nat (- int_floor_fl x)" by auto |
|
1970 hence "0 < ?num" by auto |
|
1971 hence "real ?num \<noteq> 0" by auto |
|
1972 have num_eq: "real ?num = - int_floor_fl x" |
|
1973 using \<open>0 < nat (- int_floor_fl x)\<close> by auto |
|
1974 have "0 < - int_floor_fl x" |
|
1975 using \<open>0 < ?num\<close>[unfolded of_nat_less_iff[symmetric]] by simp |
|
1976 hence "real_of_int (int_floor_fl x) < 0" |
|
1977 unfolding less_float_def by auto |
|
1978 have fl_eq: "real_of_int (- int_floor_fl x) = real_of_float (- floor_fl x)" |
|
1979 by (simp add: floor_fl_def int_floor_fl_def) |
|
1980 from \<open>0 < - int_floor_fl x\<close> have "0 \<le> real_of_float (- floor_fl x)" |
|
1981 by (simp add: floor_fl_def int_floor_fl_def) |
|
1982 from \<open>real_of_int (int_floor_fl x) < 0\<close> have "real_of_float (floor_fl x) < 0" |
|
1983 by (simp add: floor_fl_def int_floor_fl_def) |
|
1984 have "exp x \<le> ub_exp prec x" |
|
1985 proof - |
|
1986 have div_less_zero: "real_of_float (float_divr prec x (- floor_fl x)) \<le> 0" |
|
1987 using float_divr_nonpos_pos_upper_bound[OF \<open>real_of_float x \<le> 0\<close> \<open>0 \<le> real_of_float (- floor_fl x)\<close>] |
|
1988 unfolding less_eq_float_def zero_float.rep_eq . |
|
1989 |
|
1990 have "exp x = exp (?num * (x / ?num))" |
|
1991 using \<open>real ?num \<noteq> 0\<close> by auto |
|
1992 also have "\<dots> = exp (x / ?num) ^ ?num" |
|
1993 unfolding exp_of_nat_mult .. |
|
1994 also have "\<dots> \<le> exp (float_divr prec x (- floor_fl x)) ^ ?num" |
|
1995 unfolding num_eq fl_eq |
|
1996 by (rule power_mono, rule exp_le_cancel_iff[THEN iffD2], rule float_divr) auto |
|
1997 also have "\<dots> \<le> (?ub_exp_horner (float_divr prec x (- floor_fl x))) ^ ?num" |
|
1998 unfolding real_of_float_power |
|
1999 by (rule power_mono, rule bnds_exp_horner[OF div_less_zero, unfolded atLeastAtMost_iff, THEN conjunct2], auto) |
|
2000 also have "\<dots> \<le> real_of_float (power_up_fl prec (?ub_exp_horner (float_divr prec x (- floor_fl x))) ?num)" |
|
2001 by (auto simp add: real_power_up_fl intro!: power_up ub_exp_horner_nonneg div_less_zero) |
|
2002 finally show ?thesis |
|
2003 unfolding ub_exp.simps if_not_P[OF \<open>\<not> 0 < x\<close>] if_P[OF \<open>x < - 1\<close>] floor_fl_def Let_def . |
|
2004 qed |
|
2005 moreover |
|
2006 have "lb_exp prec x \<le> exp x" |
|
2007 proof - |
|
2008 let ?divl = "float_divl prec x (- floor_fl x)" |
|
2009 let ?horner = "?lb_exp_horner ?divl" |
|
2010 |
|
2011 show ?thesis |
|
2012 proof (cases "?horner \<le> 0") |
|
2013 case False |
|
2014 hence "0 \<le> real_of_float ?horner" by auto |
|
2015 |
|
2016 have div_less_zero: "real_of_float (float_divl prec x (- floor_fl x)) \<le> 0" |
|
2017 using \<open>real_of_float (floor_fl x) < 0\<close> \<open>real_of_float x \<le> 0\<close> |
|
2018 by (auto intro!: order_trans[OF float_divl] divide_nonpos_neg) |
|
2019 |
|
2020 have "(?lb_exp_horner (float_divl prec x (- floor_fl x))) ^ ?num \<le> |
|
2021 exp (float_divl prec x (- floor_fl x)) ^ ?num" |
|
2022 using \<open>0 \<le> real_of_float ?horner\<close>[unfolded floor_fl_def[symmetric]] |
|
2023 bnds_exp_horner[OF div_less_zero, unfolded atLeastAtMost_iff, THEN conjunct1] |
|
2024 by (auto intro!: power_mono) |
|
2025 also have "\<dots> \<le> exp (x / ?num) ^ ?num" |
|
2026 unfolding num_eq fl_eq |
|
2027 using float_divl by (auto intro!: power_mono simp del: uminus_float.rep_eq) |
|
2028 also have "\<dots> = exp (?num * (x / ?num))" |
|
2029 unfolding exp_of_nat_mult .. |
|
2030 also have "\<dots> = exp x" |
|
2031 using \<open>real ?num \<noteq> 0\<close> by auto |
|
2032 finally show ?thesis |
|
2033 using False |
|
2034 unfolding lb_exp.simps if_not_P[OF \<open>\<not> 0 < x\<close>] if_P[OF \<open>x < - 1\<close>] |
|
2035 int_floor_fl_def Let_def if_not_P[OF False] |
|
2036 by (auto simp: real_power_down_fl intro!: power_down_le) |
|
2037 next |
|
2038 case True |
|
2039 have "power_down_fl prec (Float 1 (- 2)) ?num \<le> (Float 1 (- 2)) ^ ?num" |
|
2040 by (metis Float_le_zero_iff less_imp_le linorder_not_less |
|
2041 not_numeral_le_zero numeral_One power_down_fl) |
|
2042 then have "power_down_fl prec (Float 1 (- 2)) ?num \<le> real_of_float (Float 1 (- 2)) ^ ?num" |
|
2043 by simp |
|
2044 also |
|
2045 have "real_of_float (floor_fl x) \<noteq> 0" and "real_of_float (floor_fl x) \<le> 0" |
|
2046 using \<open>real_of_float (floor_fl x) < 0\<close> by auto |
|
2047 from divide_right_mono_neg[OF floor_fl[of x] \<open>real_of_float (floor_fl x) \<le> 0\<close>, unfolded divide_self[OF \<open>real_of_float (floor_fl x) \<noteq> 0\<close>]] |
|
2048 have "- 1 \<le> x / (- floor_fl x)" |
|
2049 unfolding minus_float.rep_eq by auto |
|
2050 from order_trans[OF exp_m1_ge_quarter this[unfolded exp_le_cancel_iff[where x="- 1", symmetric]]] |
|
2051 have "Float 1 (- 2) \<le> exp (x / (- floor_fl x))" |
|
2052 unfolding Float_num . |
|
2053 hence "real_of_float (Float 1 (- 2)) ^ ?num \<le> exp (x / (- floor_fl x)) ^ ?num" |
|
2054 by (metis Float_num(5) power_mono zero_le_divide_1_iff zero_le_numeral) |
|
2055 also have "\<dots> = exp x" |
|
2056 unfolding num_eq fl_eq exp_of_nat_mult[symmetric] |
|
2057 using \<open>real_of_float (floor_fl x) \<noteq> 0\<close> by auto |
|
2058 finally show ?thesis |
|
2059 unfolding lb_exp.simps if_not_P[OF \<open>\<not> 0 < x\<close>] if_P[OF \<open>x < - 1\<close>] |
|
2060 int_floor_fl_def Let_def if_P[OF True] real_of_float_power . |
|
2061 qed |
|
2062 qed |
|
2063 ultimately show ?thesis by auto |
|
2064 qed |
|
2065 qed |
|
2066 |
|
2067 lemma exp_boundaries: "exp x \<in> { lb_exp prec x .. ub_exp prec x }" |
|
2068 proof - |
|
2069 show ?thesis |
|
2070 proof (cases "0 < x") |
|
2071 case False |
|
2072 hence "x \<le> 0" by auto |
|
2073 from exp_boundaries'[OF this] show ?thesis . |
|
2074 next |
|
2075 case True |
|
2076 hence "-x \<le> 0" by auto |
|
2077 |
|
2078 have "lb_exp prec x \<le> exp x" |
|
2079 proof - |
|
2080 from exp_boundaries'[OF \<open>-x \<le> 0\<close>] |
|
2081 have ub_exp: "exp (- real_of_float x) \<le> ub_exp prec (-x)" |
|
2082 unfolding atLeastAtMost_iff minus_float.rep_eq by auto |
|
2083 |
|
2084 have "float_divl prec 1 (ub_exp prec (-x)) \<le> 1 / ub_exp prec (-x)" |
|
2085 using float_divl[where x=1] by auto |
|
2086 also have "\<dots> \<le> exp x" |
|
2087 using ub_exp[unfolded inverse_le_iff_le[OF order_less_le_trans[OF exp_gt_zero ub_exp] |
|
2088 exp_gt_zero, symmetric]] |
|
2089 unfolding exp_minus nonzero_inverse_inverse_eq[OF exp_not_eq_zero] inverse_eq_divide |
|
2090 by auto |
|
2091 finally show ?thesis |
|
2092 unfolding lb_exp.simps if_P[OF True] . |
|
2093 qed |
|
2094 moreover |
|
2095 have "exp x \<le> ub_exp prec x" |
|
2096 proof - |
|
2097 have "\<not> 0 < -x" using \<open>0 < x\<close> by auto |
|
2098 |
|
2099 from exp_boundaries'[OF \<open>-x \<le> 0\<close>] |
|
2100 have lb_exp: "lb_exp prec (-x) \<le> exp (- real_of_float x)" |
|
2101 unfolding atLeastAtMost_iff minus_float.rep_eq by auto |
|
2102 |
|
2103 have "exp x \<le> (1 :: float) / lb_exp prec (-x)" |
|
2104 using lb_exp lb_exp_pos[OF \<open>\<not> 0 < -x\<close>, of prec] |
|
2105 by (simp del: lb_exp.simps add: exp_minus field_simps) |
|
2106 also have "\<dots> \<le> float_divr prec 1 (lb_exp prec (-x))" |
|
2107 using float_divr . |
|
2108 finally show ?thesis |
|
2109 unfolding ub_exp.simps if_P[OF True] . |
|
2110 qed |
|
2111 ultimately show ?thesis |
|
2112 by auto |
|
2113 qed |
|
2114 qed |
|
2115 |
|
2116 lemma bnds_exp: "\<forall>(x::real) lx ux. (l, u) = |
|
2117 (lb_exp prec lx, ub_exp prec ux) \<and> x \<in> {lx .. ux} \<longrightarrow> l \<le> exp x \<and> exp x \<le> u" |
|
2118 proof (rule allI, rule allI, rule allI, rule impI) |
|
2119 fix x :: real and lx ux |
|
2120 assume "(l, u) = (lb_exp prec lx, ub_exp prec ux) \<and> x \<in> {lx .. ux}" |
|
2121 hence l: "lb_exp prec lx = l " and u: "ub_exp prec ux = u" and x: "x \<in> {lx .. ux}" |
|
2122 by auto |
|
2123 show "l \<le> exp x \<and> exp x \<le> u" |
|
2124 proof |
|
2125 show "l \<le> exp x" |
|
2126 proof - |
|
2127 from exp_boundaries[of lx prec, unfolded l] |
|
2128 have "l \<le> exp lx" by (auto simp del: lb_exp.simps) |
|
2129 also have "\<dots> \<le> exp x" using x by auto |
|
2130 finally show ?thesis . |
|
2131 qed |
|
2132 show "exp x \<le> u" |
|
2133 proof - |
|
2134 have "exp x \<le> exp ux" using x by auto |
|
2135 also have "\<dots> \<le> u" using exp_boundaries[of ux prec, unfolded u] by (auto simp del: ub_exp.simps) |
|
2136 finally show ?thesis . |
|
2137 qed |
|
2138 qed |
|
2139 qed |
|
2140 |
|
2141 |
|
2142 section "Logarithm" |
|
2143 |
|
2144 subsection "Compute the logarithm series" |
|
2145 |
|
2146 fun ub_ln_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" |
|
2147 and lb_ln_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where |
|
2148 "ub_ln_horner prec 0 i x = 0" | |
|
2149 "ub_ln_horner prec (Suc n) i x = float_plus_up prec |
|
2150 (rapprox_rat prec 1 (int i)) (- float_round_down prec (x * lb_ln_horner prec n (Suc i) x))" | |
|
2151 "lb_ln_horner prec 0 i x = 0" | |
|
2152 "lb_ln_horner prec (Suc n) i x = float_plus_down prec |
|
2153 (lapprox_rat prec 1 (int i)) (- float_round_up prec (x * ub_ln_horner prec n (Suc i) x))" |
|
2154 |
|
2155 lemma ln_bounds: |
|
2156 assumes "0 \<le> x" |
|
2157 and "x < 1" |
|
2158 shows "(\<Sum>i=0..<2*n. (- 1) ^ i * (1 / real (i + 1)) * x ^ (Suc i)) \<le> ln (x + 1)" (is "?lb") |
|
2159 and "ln (x + 1) \<le> (\<Sum>i=0..<2*n + 1. (- 1) ^ i * (1 / real (i + 1)) * x ^ (Suc i))" (is "?ub") |
|
2160 proof - |
|
2161 let "?a n" = "(1/real (n +1)) * x ^ (Suc n)" |
|
2162 |
|
2163 have ln_eq: "(\<Sum> i. (- 1) ^ i * ?a i) = ln (x + 1)" |
|
2164 using ln_series[of "x + 1"] \<open>0 \<le> x\<close> \<open>x < 1\<close> by auto |
|
2165 |
|
2166 have "norm x < 1" using assms by auto |
|
2167 have "?a \<longlonglongrightarrow> 0" unfolding Suc_eq_plus1[symmetric] inverse_eq_divide[symmetric] |
|
2168 using tendsto_mult[OF LIMSEQ_inverse_real_of_nat LIMSEQ_Suc[OF LIMSEQ_power_zero[OF \<open>norm x < 1\<close>]]] by auto |
|
2169 have "0 \<le> ?a n" for n |
|
2170 by (rule mult_nonneg_nonneg) (auto simp: \<open>0 \<le> x\<close>) |
|
2171 have "?a (Suc n) \<le> ?a n" for n |
|
2172 unfolding inverse_eq_divide[symmetric] |
|
2173 proof (rule mult_mono) |
|
2174 show "0 \<le> x ^ Suc (Suc n)" |
|
2175 by (auto simp add: \<open>0 \<le> x\<close>) |
|
2176 have "x ^ Suc (Suc n) \<le> x ^ Suc n * 1" |
|
2177 unfolding power_Suc2 mult.assoc[symmetric] |
|
2178 by (rule mult_left_mono, fact less_imp_le[OF \<open>x < 1\<close>]) (auto simp: \<open>0 \<le> x\<close>) |
|
2179 thus "x ^ Suc (Suc n) \<le> x ^ Suc n" by auto |
|
2180 qed auto |
|
2181 from summable_Leibniz'(2,4)[OF \<open>?a \<longlonglongrightarrow> 0\<close> \<open>\<And>n. 0 \<le> ?a n\<close>, OF \<open>\<And>n. ?a (Suc n) \<le> ?a n\<close>, unfolded ln_eq] |
|
2182 show ?lb and ?ub |
|
2183 unfolding atLeast0LessThan by auto |
|
2184 qed |
|
2185 |
|
2186 lemma ln_float_bounds: |
|
2187 assumes "0 \<le> real_of_float x" |
|
2188 and "real_of_float x < 1" |
|
2189 shows "x * lb_ln_horner prec (get_even n) 1 x \<le> ln (x + 1)" (is "?lb \<le> ?ln") |
|
2190 and "ln (x + 1) \<le> x * ub_ln_horner prec (get_odd n) 1 x" (is "?ln \<le> ?ub") |
|
2191 proof - |
|
2192 obtain ev where ev: "get_even n = 2 * ev" using get_even_double .. |
|
2193 obtain od where od: "get_odd n = 2 * od + 1" using get_odd_double .. |
|
2194 |
|
2195 let "?s n" = "(- 1) ^ n * (1 / real (1 + n)) * (real_of_float x)^(Suc n)" |
|
2196 |
|
2197 have "?lb \<le> sum ?s {0 ..< 2 * ev}" |
|
2198 unfolding power_Suc2 mult.assoc[symmetric] times_float.rep_eq sum_distrib_right[symmetric] |
|
2199 unfolding mult.commute[of "real_of_float x"] ev |
|
2200 using horner_bounds(1)[where G="\<lambda> i k. Suc k" and F="\<lambda>x. x" and f="\<lambda>x. x" |
|
2201 and lb="\<lambda>n i k x. lb_ln_horner prec n k x" |
|
2202 and ub="\<lambda>n i k x. ub_ln_horner prec n k x" and j'=1 and n="2*ev", |
|
2203 OF \<open>0 \<le> real_of_float x\<close> refl lb_ln_horner.simps ub_ln_horner.simps] \<open>0 \<le> real_of_float x\<close> |
|
2204 unfolding real_of_float_power |
|
2205 by (rule mult_right_mono) |
|
2206 also have "\<dots> \<le> ?ln" |
|
2207 using ln_bounds(1)[OF \<open>0 \<le> real_of_float x\<close> \<open>real_of_float x < 1\<close>] by auto |
|
2208 finally show "?lb \<le> ?ln" . |
|
2209 |
|
2210 have "?ln \<le> sum ?s {0 ..< 2 * od + 1}" |
|
2211 using ln_bounds(2)[OF \<open>0 \<le> real_of_float x\<close> \<open>real_of_float x < 1\<close>] by auto |
|
2212 also have "\<dots> \<le> ?ub" |
|
2213 unfolding power_Suc2 mult.assoc[symmetric] times_float.rep_eq sum_distrib_right[symmetric] |
|
2214 unfolding mult.commute[of "real_of_float x"] od |
|
2215 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", |
|
2216 OF \<open>0 \<le> real_of_float x\<close> refl lb_ln_horner.simps ub_ln_horner.simps] \<open>0 \<le> real_of_float x\<close> |
|
2217 unfolding real_of_float_power |
|
2218 by (rule mult_right_mono) |
|
2219 finally show "?ln \<le> ?ub" . |
|
2220 qed |
|
2221 |
|
2222 lemma ln_add: |
|
2223 fixes x :: real |
|
2224 assumes "0 < x" and "0 < y" |
|
2225 shows "ln (x + y) = ln x + ln (1 + y / x)" |
|
2226 proof - |
|
2227 have "x \<noteq> 0" using assms by auto |
|
2228 have "x + y = x * (1 + y / x)" |
|
2229 unfolding distrib_left times_divide_eq_right nonzero_mult_div_cancel_left[OF \<open>x \<noteq> 0\<close>] |
|
2230 by auto |
|
2231 moreover |
|
2232 have "0 < y / x" using assms by auto |
|
2233 hence "0 < 1 + y / x" by auto |
|
2234 ultimately show ?thesis |
|
2235 using ln_mult assms by auto |
|
2236 qed |
|
2237 |
|
2238 |
|
2239 subsection "Compute the logarithm of 2" |
|
2240 |
|
2241 definition ub_ln2 where "ub_ln2 prec = (let third = rapprox_rat (max prec 1) 1 3 |
|
2242 in float_plus_up prec |
|
2243 ((Float 1 (- 1) * ub_ln_horner prec (get_odd prec) 1 (Float 1 (- 1)))) |
|
2244 (float_round_up prec (third * ub_ln_horner prec (get_odd prec) 1 third)))" |
|
2245 definition lb_ln2 where "lb_ln2 prec = (let third = lapprox_rat prec 1 3 |
|
2246 in float_plus_down prec |
|
2247 ((Float 1 (- 1) * lb_ln_horner prec (get_even prec) 1 (Float 1 (- 1)))) |
|
2248 (float_round_down prec (third * lb_ln_horner prec (get_even prec) 1 third)))" |
|
2249 |
|
2250 lemma ub_ln2: "ln 2 \<le> ub_ln2 prec" (is "?ub_ln2") |
|
2251 and lb_ln2: "lb_ln2 prec \<le> ln 2" (is "?lb_ln2") |
|
2252 proof - |
|
2253 let ?uthird = "rapprox_rat (max prec 1) 1 3" |
|
2254 let ?lthird = "lapprox_rat prec 1 3" |
|
2255 |
|
2256 have ln2_sum: "ln 2 = ln (1/2 + 1) + ln (1 / 3 + 1::real)" |
|
2257 using ln_add[of "3 / 2" "1 / 2"] by auto |
|
2258 have lb3: "?lthird \<le> 1 / 3" using lapprox_rat[of prec 1 3] by auto |
|
2259 hence lb3_ub: "real_of_float ?lthird < 1" by auto |
|
2260 have lb3_lb: "0 \<le> real_of_float ?lthird" using lapprox_rat_nonneg[of 1 3] by auto |
|
2261 have ub3: "1 / 3 \<le> ?uthird" using rapprox_rat[of 1 3] by auto |
|
2262 hence ub3_lb: "0 \<le> real_of_float ?uthird" by auto |
|
2263 |
|
2264 have lb2: "0 \<le> real_of_float (Float 1 (- 1))" and ub2: "real_of_float (Float 1 (- 1)) < 1" |
|
2265 unfolding Float_num by auto |
|
2266 |
|
2267 have "0 \<le> (1::int)" and "0 < (3::int)" by auto |
|
2268 have ub3_ub: "real_of_float ?uthird < 1" |
|
2269 by (simp add: Float.compute_rapprox_rat Float.compute_lapprox_rat rapprox_posrat_less1) |
|
2270 |
|
2271 have third_gt0: "(0 :: real) < 1 / 3 + 1" by auto |
|
2272 have uthird_gt0: "0 < real_of_float ?uthird + 1" using ub3_lb by auto |
|
2273 have lthird_gt0: "0 < real_of_float ?lthird + 1" using lb3_lb by auto |
|
2274 |
|
2275 show ?ub_ln2 |
|
2276 unfolding ub_ln2_def Let_def ln2_sum Float_num(4)[symmetric] |
|
2277 proof (rule float_plus_up_le, rule add_mono, fact ln_float_bounds(2)[OF lb2 ub2]) |
|
2278 have "ln (1 / 3 + 1) \<le> ln (real_of_float ?uthird + 1)" |
|
2279 unfolding ln_le_cancel_iff[OF third_gt0 uthird_gt0] using ub3 by auto |
|
2280 also have "\<dots> \<le> ?uthird * ub_ln_horner prec (get_odd prec) 1 ?uthird" |
|
2281 using ln_float_bounds(2)[OF ub3_lb ub3_ub] . |
|
2282 also note float_round_up |
|
2283 finally show "ln (1 / 3 + 1) \<le> float_round_up prec (?uthird * ub_ln_horner prec (get_odd prec) 1 ?uthird)" . |
|
2284 qed |
|
2285 show ?lb_ln2 |
|
2286 unfolding lb_ln2_def Let_def ln2_sum Float_num(4)[symmetric] |
|
2287 proof (rule float_plus_down_le, rule add_mono, fact ln_float_bounds(1)[OF lb2 ub2]) |
|
2288 have "?lthird * lb_ln_horner prec (get_even prec) 1 ?lthird \<le> ln (real_of_float ?lthird + 1)" |
|
2289 using ln_float_bounds(1)[OF lb3_lb lb3_ub] . |
|
2290 note float_round_down_le[OF this] |
|
2291 also have "\<dots> \<le> ln (1 / 3 + 1)" |
|
2292 unfolding ln_le_cancel_iff[OF lthird_gt0 third_gt0] |
|
2293 using lb3 by auto |
|
2294 finally show "float_round_down prec (?lthird * lb_ln_horner prec (get_even prec) 1 ?lthird) \<le> |
|
2295 ln (1 / 3 + 1)" . |
|
2296 qed |
|
2297 qed |
|
2298 |
|
2299 |
|
2300 subsection "Compute the logarithm in the entire domain" |
|
2301 |
|
2302 function ub_ln :: "nat \<Rightarrow> float \<Rightarrow> float option" and lb_ln :: "nat \<Rightarrow> float \<Rightarrow> float option" where |
|
2303 "ub_ln prec x = (if x \<le> 0 then None |
|
2304 else if x < 1 then Some (- the (lb_ln prec (float_divl (max prec 1) 1 x))) |
|
2305 else let horner = \<lambda>x. float_round_up prec (x * ub_ln_horner prec (get_odd prec) 1 x) in |
|
2306 if x \<le> Float 3 (- 1) then Some (horner (x - 1)) |
|
2307 else if x < Float 1 1 then Some (float_round_up prec (horner (Float 1 (- 1)) + horner (x * rapprox_rat prec 2 3 - 1))) |
|
2308 else let l = bitlen (mantissa x) - 1 in |
|
2309 Some (float_plus_up prec (float_round_up prec (ub_ln2 prec * (Float (exponent x + l) 0))) (horner (Float (mantissa x) (- l) - 1))))" | |
|
2310 "lb_ln prec x = (if x \<le> 0 then None |
|
2311 else if x < 1 then Some (- the (ub_ln prec (float_divr prec 1 x))) |
|
2312 else let horner = \<lambda>x. float_round_down prec (x * lb_ln_horner prec (get_even prec) 1 x) in |
|
2313 if x \<le> Float 3 (- 1) then Some (horner (x - 1)) |
|
2314 else if x < Float 1 1 then Some (float_round_down prec (horner (Float 1 (- 1)) + |
|
2315 horner (max (x * lapprox_rat prec 2 3 - 1) 0))) |
|
2316 else let l = bitlen (mantissa x) - 1 in |
|
2317 Some (float_plus_down prec (float_round_down prec (lb_ln2 prec * (Float (exponent x + l) 0))) (horner (Float (mantissa x) (- l) - 1))))" |
|
2318 by pat_completeness auto |
|
2319 |
|
2320 termination |
|
2321 proof (relation "measure (\<lambda> v. let (prec, x) = case_sum id id v in (if x < 1 then 1 else 0))", auto) |
|
2322 fix prec and x :: float |
|
2323 assume "\<not> real_of_float x \<le> 0" and "real_of_float x < 1" and "real_of_float (float_divl (max prec (Suc 0)) 1 x) < 1" |
|
2324 hence "0 < real_of_float x" "1 \<le> max prec (Suc 0)" "real_of_float x < 1" |
|
2325 by auto |
|
2326 from float_divl_pos_less1_bound[OF \<open>0 < real_of_float x\<close> \<open>real_of_float x < 1\<close>[THEN less_imp_le] \<open>1 \<le> max prec (Suc 0)\<close>] |
|
2327 show False |
|
2328 using \<open>real_of_float (float_divl (max prec (Suc 0)) 1 x) < 1\<close> by auto |
|
2329 next |
|
2330 fix prec x |
|
2331 assume "\<not> real_of_float x \<le> 0" and "real_of_float x < 1" and "real_of_float (float_divr prec 1 x) < 1" |
|
2332 hence "0 < x" by auto |
|
2333 from float_divr_pos_less1_lower_bound[OF \<open>0 < x\<close>, of prec] \<open>real_of_float x < 1\<close> show False |
|
2334 using \<open>real_of_float (float_divr prec 1 x) < 1\<close> by auto |
|
2335 qed |
|
2336 |
|
2337 lemma float_pos_eq_mantissa_pos: "x > 0 \<longleftrightarrow> mantissa x > 0" |
|
2338 apply (subst Float_mantissa_exponent[of x, symmetric]) |
|
2339 apply (auto simp add: zero_less_mult_iff zero_float_def dest: less_zeroE) |
|
2340 apply (metis not_le powr_ge_pzero) |
|
2341 done |
|
2342 |
|
2343 lemma Float_pos_eq_mantissa_pos: "Float m e > 0 \<longleftrightarrow> m > 0" |
|
2344 using powr_gt_zero[of 2 "e"] |
|
2345 by (auto simp add: zero_less_mult_iff zero_float_def simp del: powr_gt_zero dest: less_zeroE) |
|
2346 |
|
2347 lemma Float_representation_aux: |
|
2348 fixes m e |
|
2349 defines "x \<equiv> Float m e" |
|
2350 assumes "x > 0" |
|
2351 shows "Float (exponent x + (bitlen (mantissa x) - 1)) 0 = Float (e + (bitlen m - 1)) 0" (is ?th1) |
|
2352 and "Float (mantissa x) (- (bitlen (mantissa x) - 1)) = Float m ( - (bitlen m - 1))" (is ?th2) |
|
2353 proof - |
|
2354 from assms have mantissa_pos: "m > 0" "mantissa x > 0" |
|
2355 using Float_pos_eq_mantissa_pos[of m e] float_pos_eq_mantissa_pos[of x] by simp_all |
|
2356 thus ?th1 |
|
2357 using bitlen_Float[of m e] assms |
|
2358 by (auto simp add: zero_less_mult_iff intro!: arg_cong2[where f=Float]) |
|
2359 have "x \<noteq> float_of 0" |
|
2360 unfolding zero_float_def[symmetric] using \<open>0 < x\<close> by auto |
|
2361 from denormalize_shift[OF assms(1) this] guess i . note i = this |
|
2362 |
|
2363 have "2 powr (1 - (real_of_int (bitlen (mantissa x)) + real_of_int i)) = |
|
2364 2 powr (1 - (real_of_int (bitlen (mantissa x)))) * inverse (2 powr (real i))" |
|
2365 by (simp add: powr_minus[symmetric] powr_add[symmetric] field_simps) |
|
2366 hence "real_of_int (mantissa x) * 2 powr (1 - real_of_int (bitlen (mantissa x))) = |
|
2367 (real_of_int (mantissa x) * 2 ^ i) * 2 powr (1 - real_of_int (bitlen (mantissa x * 2 ^ i)))" |
|
2368 using \<open>mantissa x > 0\<close> by (simp add: powr_realpow) |
|
2369 then show ?th2 |
|
2370 unfolding i by transfer auto |
|
2371 qed |
|
2372 |
|
2373 lemma compute_ln[code]: |
|
2374 fixes m e |
|
2375 defines "x \<equiv> Float m e" |
|
2376 shows "ub_ln prec x = (if x \<le> 0 then None |
|
2377 else if x < 1 then Some (- the (lb_ln prec (float_divl (max prec 1) 1 x))) |
|
2378 else let horner = \<lambda>x. float_round_up prec (x * ub_ln_horner prec (get_odd prec) 1 x) in |
|
2379 if x \<le> Float 3 (- 1) then Some (horner (x - 1)) |
|
2380 else if x < Float 1 1 then Some (float_round_up prec (horner (Float 1 (- 1)) + horner (x * rapprox_rat prec 2 3 - 1))) |
|
2381 else let l = bitlen m - 1 in |
|
2382 Some (float_plus_up prec (float_round_up prec (ub_ln2 prec * (Float (e + l) 0))) (horner (Float m (- l) - 1))))" |
|
2383 (is ?th1) |
|
2384 and "lb_ln prec x = (if x \<le> 0 then None |
|
2385 else if x < 1 then Some (- the (ub_ln prec (float_divr prec 1 x))) |
|
2386 else let horner = \<lambda>x. float_round_down prec (x * lb_ln_horner prec (get_even prec) 1 x) in |
|
2387 if x \<le> Float 3 (- 1) then Some (horner (x - 1)) |
|
2388 else if x < Float 1 1 then Some (float_round_down prec (horner (Float 1 (- 1)) + |
|
2389 horner (max (x * lapprox_rat prec 2 3 - 1) 0))) |
|
2390 else let l = bitlen m - 1 in |
|
2391 Some (float_plus_down prec (float_round_down prec (lb_ln2 prec * (Float (e + l) 0))) (horner (Float m (- l) - 1))))" |
|
2392 (is ?th2) |
|
2393 proof - |
|
2394 from assms Float_pos_eq_mantissa_pos have "x > 0 \<Longrightarrow> m > 0" |
|
2395 by simp |
|
2396 thus ?th1 ?th2 |
|
2397 using Float_representation_aux[of m e] |
|
2398 unfolding x_def[symmetric] |
|
2399 by (auto dest: not_le_imp_less) |
|
2400 qed |
|
2401 |
|
2402 lemma ln_shifted_float: |
|
2403 assumes "0 < m" |
|
2404 shows "ln (Float m e) = ln 2 * (e + (bitlen m - 1)) + ln (Float m (- (bitlen m - 1)))" |
|
2405 proof - |
|
2406 let ?B = "2^nat (bitlen m - 1)" |
|
2407 define bl where "bl = bitlen m - 1" |
|
2408 have "0 < real_of_int m" and "\<And>X. (0 :: real) < 2^X" and "0 < (2 :: real)" and "m \<noteq> 0" |
|
2409 using assms by auto |
|
2410 hence "0 \<le> bl" by (simp add: bitlen_alt_def bl_def) |
|
2411 show ?thesis |
|
2412 proof (cases "0 \<le> e") |
|
2413 case True |
|
2414 thus ?thesis |
|
2415 unfolding bl_def[symmetric] using \<open>0 < real_of_int m\<close> \<open>0 \<le> bl\<close> |
|
2416 apply (simp add: ln_mult) |
|
2417 apply (cases "e=0") |
|
2418 apply (cases "bl = 0", simp_all add: powr_minus ln_inverse ln_powr) |
|
2419 apply (cases "bl = 0", simp_all add: powr_minus ln_inverse ln_powr field_simps) |
|
2420 done |
|
2421 next |
|
2422 case False |
|
2423 hence "0 < -e" by auto |
|
2424 have lne: "ln (2 powr real_of_int e) = ln (inverse (2 powr - e))" |
|
2425 by (simp add: powr_minus) |
|
2426 hence pow_gt0: "(0::real) < 2^nat (-e)" |
|
2427 by auto |
|
2428 hence inv_gt0: "(0::real) < inverse (2^nat (-e))" |
|
2429 by auto |
|
2430 show ?thesis |
|
2431 using False unfolding bl_def[symmetric] |
|
2432 using \<open>0 < real_of_int m\<close> \<open>0 \<le> bl\<close> |
|
2433 by (auto simp add: lne ln_mult ln_powr ln_div field_simps) |
|
2434 qed |
|
2435 qed |
|
2436 |
|
2437 lemma ub_ln_lb_ln_bounds': |
|
2438 assumes "1 \<le> x" |
|
2439 shows "the (lb_ln prec x) \<le> ln x \<and> ln x \<le> the (ub_ln prec x)" |
|
2440 (is "?lb \<le> ?ln \<and> ?ln \<le> ?ub") |
|
2441 proof (cases "x < Float 1 1") |
|
2442 case True |
|
2443 hence "real_of_float (x - 1) < 1" and "real_of_float x < 2" by auto |
|
2444 have "\<not> x \<le> 0" and "\<not> x < 1" using \<open>1 \<le> x\<close> by auto |
|
2445 hence "0 \<le> real_of_float (x - 1)" using \<open>1 \<le> x\<close> by auto |
|
2446 |
|
2447 have [simp]: "(Float 3 (- 1)) = 3 / 2" by simp |
|
2448 |
|
2449 show ?thesis |
|
2450 proof (cases "x \<le> Float 3 (- 1)") |
|
2451 case True |
|
2452 show ?thesis |
|
2453 unfolding lb_ln.simps |
|
2454 unfolding ub_ln.simps Let_def |
|
2455 using ln_float_bounds[OF \<open>0 \<le> real_of_float (x - 1)\<close> \<open>real_of_float (x - 1) < 1\<close>, of prec] |
|
2456 \<open>\<not> x \<le> 0\<close> \<open>\<not> x < 1\<close> True |
|
2457 by (auto intro!: float_round_down_le float_round_up_le) |
|
2458 next |
|
2459 case False |
|
2460 hence *: "3 / 2 < x" by auto |
|
2461 |
|
2462 with ln_add[of "3 / 2" "x - 3 / 2"] |
|
2463 have add: "ln x = ln (3 / 2) + ln (real_of_float x * 2 / 3)" |
|
2464 by (auto simp add: algebra_simps diff_divide_distrib) |
|
2465 |
|
2466 let "?ub_horner x" = "float_round_up prec (x * ub_ln_horner prec (get_odd prec) 1 x)" |
|
2467 let "?lb_horner x" = "float_round_down prec (x * lb_ln_horner prec (get_even prec) 1 x)" |
|
2468 |
|
2469 { have up: "real_of_float (rapprox_rat prec 2 3) \<le> 1" |
|
2470 by (rule rapprox_rat_le1) simp_all |
|
2471 have low: "2 / 3 \<le> rapprox_rat prec 2 3" |
|
2472 by (rule order_trans[OF _ rapprox_rat]) simp |
|
2473 from mult_less_le_imp_less[OF * low] * |
|
2474 have pos: "0 < real_of_float (x * rapprox_rat prec 2 3 - 1)" by auto |
|
2475 |
|
2476 have "ln (real_of_float x * 2/3) |
|
2477 \<le> ln (real_of_float (x * rapprox_rat prec 2 3 - 1) + 1)" |
|
2478 proof (rule ln_le_cancel_iff[symmetric, THEN iffD1]) |
|
2479 show "real_of_float x * 2 / 3 \<le> real_of_float (x * rapprox_rat prec 2 3 - 1) + 1" |
|
2480 using * low by auto |
|
2481 show "0 < real_of_float x * 2 / 3" using * by simp |
|
2482 show "0 < real_of_float (x * rapprox_rat prec 2 3 - 1) + 1" using pos by auto |
|
2483 qed |
|
2484 also have "\<dots> \<le> ?ub_horner (x * rapprox_rat prec 2 3 - 1)" |
|
2485 proof (rule float_round_up_le, rule ln_float_bounds(2)) |
|
2486 from mult_less_le_imp_less[OF \<open>real_of_float x < 2\<close> up] low * |
|
2487 show "real_of_float (x * rapprox_rat prec 2 3 - 1) < 1" by auto |
|
2488 show "0 \<le> real_of_float (x * rapprox_rat prec 2 3 - 1)" using pos by auto |
|
2489 qed |
|
2490 finally have "ln x \<le> ?ub_horner (Float 1 (-1)) |
|
2491 + ?ub_horner ((x * rapprox_rat prec 2 3 - 1))" |
|
2492 using ln_float_bounds(2)[of "Float 1 (- 1)" prec prec] add |
|
2493 by (auto intro!: add_mono float_round_up_le) |
|
2494 note float_round_up_le[OF this, of prec] |
|
2495 } |
|
2496 moreover |
|
2497 { let ?max = "max (x * lapprox_rat prec 2 3 - 1) 0" |
|
2498 |
|
2499 have up: "lapprox_rat prec 2 3 \<le> 2/3" |
|
2500 by (rule order_trans[OF lapprox_rat], simp) |
|
2501 |
|
2502 have low: "0 \<le> real_of_float (lapprox_rat prec 2 3)" |
|
2503 using lapprox_rat_nonneg[of 2 3 prec] by simp |
|
2504 |
|
2505 have "?lb_horner ?max |
|
2506 \<le> ln (real_of_float ?max + 1)" |
|
2507 proof (rule float_round_down_le, rule ln_float_bounds(1)) |
|
2508 from mult_less_le_imp_less[OF \<open>real_of_float x < 2\<close> up] * low |
|
2509 show "real_of_float ?max < 1" by (cases "real_of_float (lapprox_rat prec 2 3) = 0", |
|
2510 auto simp add: real_of_float_max) |
|
2511 show "0 \<le> real_of_float ?max" by (auto simp add: real_of_float_max) |
|
2512 qed |
|
2513 also have "\<dots> \<le> ln (real_of_float x * 2/3)" |
|
2514 proof (rule ln_le_cancel_iff[symmetric, THEN iffD1]) |
|
2515 show "0 < real_of_float ?max + 1" by (auto simp add: real_of_float_max) |
|
2516 show "0 < real_of_float x * 2/3" using * by auto |
|
2517 show "real_of_float ?max + 1 \<le> real_of_float x * 2/3" using * up |
|
2518 by (cases "0 < real_of_float x * real_of_float (lapprox_posrat prec 2 3) - 1", |
|
2519 auto simp add: max_def) |
|
2520 qed |
|
2521 finally have "?lb_horner (Float 1 (- 1)) + ?lb_horner ?max \<le> ln x" |
|
2522 using ln_float_bounds(1)[of "Float 1 (- 1)" prec prec] add |
|
2523 by (auto intro!: add_mono float_round_down_le) |
|
2524 note float_round_down_le[OF this, of prec] |
|
2525 } |
|
2526 ultimately |
|
2527 show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps Let_def |
|
2528 using \<open>\<not> x \<le> 0\<close> \<open>\<not> x < 1\<close> True False by auto |
|
2529 qed |
|
2530 next |
|
2531 case False |
|
2532 hence "\<not> x \<le> 0" and "\<not> x < 1" "0 < x" "\<not> x \<le> Float 3 (- 1)" |
|
2533 using \<open>1 \<le> x\<close> by auto |
|
2534 show ?thesis |
|
2535 proof - |
|
2536 define m where "m = mantissa x" |
|
2537 define e where "e = exponent x" |
|
2538 from Float_mantissa_exponent[of x] have Float: "x = Float m e" |
|
2539 by (simp add: m_def e_def) |
|
2540 let ?s = "Float (e + (bitlen m - 1)) 0" |
|
2541 let ?x = "Float m (- (bitlen m - 1))" |
|
2542 |
|
2543 have "0 < m" and "m \<noteq> 0" using \<open>0 < x\<close> Float powr_gt_zero[of 2 e] |
|
2544 apply (auto simp add: zero_less_mult_iff) |
|
2545 using not_le powr_ge_pzero apply blast |
|
2546 done |
|
2547 define bl where "bl = bitlen m - 1" |
|
2548 hence "bl \<ge> 0" |
|
2549 using \<open>m > 0\<close> by (simp add: bitlen_alt_def) |
|
2550 have "1 \<le> Float m e" |
|
2551 using \<open>1 \<le> x\<close> Float unfolding less_eq_float_def by auto |
|
2552 from bitlen_div[OF \<open>0 < m\<close>] float_gt1_scale[OF \<open>1 \<le> Float m e\<close>] \<open>bl \<ge> 0\<close> |
|
2553 have x_bnds: "0 \<le> real_of_float (?x - 1)" "real_of_float (?x - 1) < 1" |
|
2554 unfolding bl_def[symmetric] |
|
2555 by (auto simp: powr_realpow[symmetric] field_simps) |
|
2556 (auto simp : powr_minus field_simps) |
|
2557 |
|
2558 { |
|
2559 have "float_round_down prec (lb_ln2 prec * ?s) \<le> ln 2 * (e + (bitlen m - 1))" |
|
2560 (is "real_of_float ?lb2 \<le> _") |
|
2561 apply (rule float_round_down_le) |
|
2562 unfolding nat_0 power_0 mult_1_right times_float.rep_eq |
|
2563 using lb_ln2[of prec] |
|
2564 proof (rule mult_mono) |
|
2565 from float_gt1_scale[OF \<open>1 \<le> Float m e\<close>] |
|
2566 show "0 \<le> real_of_float (Float (e + (bitlen m - 1)) 0)" by simp |
|
2567 qed auto |
|
2568 moreover |
|
2569 from ln_float_bounds(1)[OF x_bnds] |
|
2570 have "float_round_down prec ((?x - 1) * lb_ln_horner prec (get_even prec) 1 (?x - 1)) \<le> ln ?x" (is "real_of_float ?lb_horner \<le> _") |
|
2571 by (auto intro!: float_round_down_le) |
|
2572 ultimately have "float_plus_down prec ?lb2 ?lb_horner \<le> ln x" |
|
2573 unfolding Float ln_shifted_float[OF \<open>0 < m\<close>, of e] by (auto intro!: float_plus_down_le) |
|
2574 } |
|
2575 moreover |
|
2576 { |
|
2577 from ln_float_bounds(2)[OF x_bnds] |
|
2578 have "ln ?x \<le> float_round_up prec ((?x - 1) * ub_ln_horner prec (get_odd prec) 1 (?x - 1))" |
|
2579 (is "_ \<le> real_of_float ?ub_horner") |
|
2580 by (auto intro!: float_round_up_le) |
|
2581 moreover |
|
2582 have "ln 2 * (e + (bitlen m - 1)) \<le> float_round_up prec (ub_ln2 prec * ?s)" |
|
2583 (is "_ \<le> real_of_float ?ub2") |
|
2584 apply (rule float_round_up_le) |
|
2585 unfolding nat_0 power_0 mult_1_right times_float.rep_eq |
|
2586 using ub_ln2[of prec] |
|
2587 proof (rule mult_mono) |
|
2588 from float_gt1_scale[OF \<open>1 \<le> Float m e\<close>] |
|
2589 show "0 \<le> real_of_int (e + (bitlen m - 1))" by auto |
|
2590 have "0 \<le> ln (2 :: real)" by simp |
|
2591 thus "0 \<le> real_of_float (ub_ln2 prec)" using ub_ln2[of prec] by arith |
|
2592 qed auto |
|
2593 ultimately have "ln x \<le> float_plus_up prec ?ub2 ?ub_horner" |
|
2594 unfolding Float ln_shifted_float[OF \<open>0 < m\<close>, of e] |
|
2595 by (auto intro!: float_plus_up_le) |
|
2596 } |
|
2597 ultimately show ?thesis |
|
2598 unfolding lb_ln.simps |
|
2599 unfolding ub_ln.simps |
|
2600 unfolding if_not_P[OF \<open>\<not> x \<le> 0\<close>] if_not_P[OF \<open>\<not> x < 1\<close>] |
|
2601 if_not_P[OF False] if_not_P[OF \<open>\<not> x \<le> Float 3 (- 1)\<close>] Let_def |
|
2602 unfolding plus_float.rep_eq e_def[symmetric] m_def[symmetric] |
|
2603 by simp |
|
2604 qed |
|
2605 qed |
|
2606 |
|
2607 lemma ub_ln_lb_ln_bounds: |
|
2608 assumes "0 < x" |
|
2609 shows "the (lb_ln prec x) \<le> ln x \<and> ln x \<le> the (ub_ln prec x)" |
|
2610 (is "?lb \<le> ?ln \<and> ?ln \<le> ?ub") |
|
2611 proof (cases "x < 1") |
|
2612 case False |
|
2613 hence "1 \<le> x" |
|
2614 unfolding less_float_def less_eq_float_def by auto |
|
2615 show ?thesis |
|
2616 using ub_ln_lb_ln_bounds'[OF \<open>1 \<le> x\<close>] . |
|
2617 next |
|
2618 case True |
|
2619 have "\<not> x \<le> 0" using \<open>0 < x\<close> by auto |
|
2620 from True have "real_of_float x \<le> 1" "x \<le> 1" |
|
2621 by simp_all |
|
2622 have "0 < real_of_float x" and "real_of_float x \<noteq> 0" |
|
2623 using \<open>0 < x\<close> by auto |
|
2624 hence A: "0 < 1 / real_of_float x" by auto |
|
2625 |
|
2626 { |
|
2627 let ?divl = "float_divl (max prec 1) 1 x" |
|
2628 have A': "1 \<le> ?divl" using float_divl_pos_less1_bound[OF \<open>0 < real_of_float x\<close> \<open>real_of_float x \<le> 1\<close>] by auto |
|
2629 hence B: "0 < real_of_float ?divl" by auto |
|
2630 |
|
2631 have "ln ?divl \<le> ln (1 / x)" unfolding ln_le_cancel_iff[OF B A] using float_divl[of _ 1 x] by auto |
|
2632 hence "ln x \<le> - ln ?divl" unfolding nonzero_inverse_eq_divide[OF \<open>real_of_float x \<noteq> 0\<close>, symmetric] ln_inverse[OF \<open>0 < real_of_float x\<close>] by auto |
|
2633 from this ub_ln_lb_ln_bounds'[OF A', THEN conjunct1, THEN le_imp_neg_le] |
|
2634 have "?ln \<le> - the (lb_ln prec ?divl)" unfolding uminus_float.rep_eq by (rule order_trans) |
|
2635 } moreover |
|
2636 { |
|
2637 let ?divr = "float_divr prec 1 x" |
|
2638 have A': "1 \<le> ?divr" using float_divr_pos_less1_lower_bound[OF \<open>0 < x\<close> \<open>x \<le> 1\<close>] unfolding less_eq_float_def less_float_def by auto |
|
2639 hence B: "0 < real_of_float ?divr" by auto |
|
2640 |
|
2641 have "ln (1 / x) \<le> ln ?divr" unfolding ln_le_cancel_iff[OF A B] using float_divr[of 1 x] by auto |
|
2642 hence "- ln ?divr \<le> ln x" unfolding nonzero_inverse_eq_divide[OF \<open>real_of_float x \<noteq> 0\<close>, symmetric] ln_inverse[OF \<open>0 < real_of_float x\<close>] by auto |
|
2643 from ub_ln_lb_ln_bounds'[OF A', THEN conjunct2, THEN le_imp_neg_le] this |
|
2644 have "- the (ub_ln prec ?divr) \<le> ?ln" unfolding uminus_float.rep_eq by (rule order_trans) |
|
2645 } |
|
2646 ultimately show ?thesis unfolding lb_ln.simps[where x=x] ub_ln.simps[where x=x] |
|
2647 unfolding if_not_P[OF \<open>\<not> x \<le> 0\<close>] if_P[OF True] by auto |
|
2648 qed |
|
2649 |
|
2650 lemma lb_ln: |
|
2651 assumes "Some y = lb_ln prec x" |
|
2652 shows "y \<le> ln x" and "0 < real_of_float x" |
|
2653 proof - |
|
2654 have "0 < x" |
|
2655 proof (rule ccontr) |
|
2656 assume "\<not> 0 < x" |
|
2657 hence "x \<le> 0" |
|
2658 unfolding less_eq_float_def less_float_def by auto |
|
2659 thus False |
|
2660 using assms by auto |
|
2661 qed |
|
2662 thus "0 < real_of_float x" by auto |
|
2663 have "the (lb_ln prec x) \<le> ln x" |
|
2664 using ub_ln_lb_ln_bounds[OF \<open>0 < x\<close>] .. |
|
2665 thus "y \<le> ln x" |
|
2666 unfolding assms[symmetric] by auto |
|
2667 qed |
|
2668 |
|
2669 lemma ub_ln: |
|
2670 assumes "Some y = ub_ln prec x" |
|
2671 shows "ln x \<le> y" and "0 < real_of_float x" |
|
2672 proof - |
|
2673 have "0 < x" |
|
2674 proof (rule ccontr) |
|
2675 assume "\<not> 0 < x" |
|
2676 hence "x \<le> 0" by auto |
|
2677 thus False |
|
2678 using assms by auto |
|
2679 qed |
|
2680 thus "0 < real_of_float x" by auto |
|
2681 have "ln x \<le> the (ub_ln prec x)" |
|
2682 using ub_ln_lb_ln_bounds[OF \<open>0 < x\<close>] .. |
|
2683 thus "ln x \<le> y" |
|
2684 unfolding assms[symmetric] by auto |
|
2685 qed |
|
2686 |
|
2687 lemma bnds_ln: "\<forall>(x::real) lx ux. (Some l, Some u) = |
|
2688 (lb_ln prec lx, ub_ln prec ux) \<and> x \<in> {lx .. ux} \<longrightarrow> l \<le> ln x \<and> ln x \<le> u" |
|
2689 proof (rule allI, rule allI, rule allI, rule impI) |
|
2690 fix x :: real |
|
2691 fix lx ux |
|
2692 assume "(Some l, Some u) = (lb_ln prec lx, ub_ln prec ux) \<and> x \<in> {lx .. ux}" |
|
2693 hence l: "Some l = lb_ln prec lx " and u: "Some u = ub_ln prec ux" and x: "x \<in> {lx .. ux}" |
|
2694 by auto |
|
2695 |
|
2696 have "ln ux \<le> u" and "0 < real_of_float ux" |
|
2697 using ub_ln u by auto |
|
2698 have "l \<le> ln lx" and "0 < real_of_float lx" and "0 < x" |
|
2699 using lb_ln[OF l] x by auto |
|
2700 |
|
2701 from ln_le_cancel_iff[OF \<open>0 < real_of_float lx\<close> \<open>0 < x\<close>] \<open>l \<le> ln lx\<close> |
|
2702 have "l \<le> ln x" |
|
2703 using x unfolding atLeastAtMost_iff by auto |
|
2704 moreover |
|
2705 from ln_le_cancel_iff[OF \<open>0 < x\<close> \<open>0 < real_of_float ux\<close>] \<open>ln ux \<le> real_of_float u\<close> |
|
2706 have "ln x \<le> u" |
|
2707 using x unfolding atLeastAtMost_iff by auto |
|
2708 ultimately show "l \<le> ln x \<and> ln x \<le> u" .. |
|
2709 qed |
|
2710 |
|
2711 |
|
2712 section \<open>Real power function\<close> |
|
2713 |
|
2714 definition bnds_powr :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> (float \<times> float) option" where |
|
2715 "bnds_powr prec l1 u1 l2 u2 = ( |
|
2716 if l1 = 0 \<and> u1 = 0 then |
|
2717 Some (0, 0) |
|
2718 else if l1 = 0 \<and> l2 \<ge> 1 then |
|
2719 let uln = the (ub_ln prec u1) |
|
2720 in Some (0, ub_exp prec (float_round_up prec (uln * (if uln \<ge> 0 then u2 else l2)))) |
|
2721 else if l1 \<le> 0 then |
|
2722 None |
|
2723 else |
|
2724 Some (map_bnds lb_exp ub_exp prec |
|
2725 (bnds_mult prec (the (lb_ln prec l1)) (the (ub_ln prec u1)) l2 u2)))" |
|
2726 |
|
2727 lemmas [simp del] = lb_exp.simps ub_exp.simps |
|
2728 |
|
2729 lemma mono_exp_real: "mono (exp :: real \<Rightarrow> real)" |
|
2730 by (auto simp: mono_def) |
|
2731 |
|
2732 lemma ub_exp_nonneg: "real_of_float (ub_exp prec x) \<ge> 0" |
|
2733 proof - |
|
2734 have "0 \<le> exp (real_of_float x)" by simp |
|
2735 also from exp_boundaries[of x prec] |
|
2736 have "\<dots> \<le> real_of_float (ub_exp prec x)" by simp |
|
2737 finally show ?thesis . |
|
2738 qed |
|
2739 |
|
2740 lemma bnds_powr: |
|
2741 assumes lu: "Some (l, u) = bnds_powr prec l1 u1 l2 u2" |
|
2742 assumes x: "x \<in> {real_of_float l1..real_of_float u1}" |
|
2743 assumes y: "y \<in> {real_of_float l2..real_of_float u2}" |
|
2744 shows "x powr y \<in> {real_of_float l..real_of_float u}" |
|
2745 proof - |
|
2746 consider "l1 = 0" "u1 = 0" | "l1 = 0" "u1 \<noteq> 0" "l2 \<ge> 1" | |
|
2747 "l1 \<le> 0" "\<not>(l1 = 0 \<and> (u1 = 0 \<or> l2 \<ge> 1))" | "l1 > 0" by force |
|
2748 thus ?thesis |
|
2749 proof cases |
|
2750 assume "l1 = 0" "u1 = 0" |
|
2751 with x lu show ?thesis by (auto simp: bnds_powr_def) |
|
2752 next |
|
2753 assume A: "l1 = 0" "u1 \<noteq> 0" "l2 \<ge> 1" |
|
2754 define uln where "uln = the (ub_ln prec u1)" |
|
2755 show ?thesis |
|
2756 proof (cases "x = 0") |
|
2757 case False |
|
2758 with A x y have "x powr y = exp (ln x * y)" by (simp add: powr_def) |
|
2759 also { |
|
2760 from A x False have "ln x \<le> ln (real_of_float u1)" by simp |
|
2761 also from ub_ln_lb_ln_bounds[of u1 prec] A y x False |
|
2762 have "ln (real_of_float u1) \<le> real_of_float uln" by (simp add: uln_def del: lb_ln.simps) |
|
2763 also from A x y have "\<dots> * y \<le> real_of_float uln * (if uln \<ge> 0 then u2 else l2)" |
|
2764 by (auto intro: mult_left_mono mult_left_mono_neg) |
|
2765 also have "\<dots> \<le> real_of_float (float_round_up prec (uln * (if uln \<ge> 0 then u2 else l2)))" |
|
2766 by (simp add: float_round_up_le) |
|
2767 finally have "ln x * y \<le> \<dots>" using A y by - simp |
|
2768 } |
|
2769 also have "exp (real_of_float (float_round_up prec (uln * (if uln \<ge> 0 then u2 else l2)))) \<le> |
|
2770 real_of_float (ub_exp prec (float_round_up prec |
|
2771 (uln * (if uln \<ge> 0 then u2 else l2))))" |
|
2772 using exp_boundaries by simp |
|
2773 finally show ?thesis using A x y lu |
|
2774 by (simp add: bnds_powr_def uln_def Let_def del: lb_ln.simps ub_ln.simps) |
|
2775 qed (insert x y lu A, simp_all add: bnds_powr_def Let_def ub_exp_nonneg |
|
2776 del: lb_ln.simps ub_ln.simps) |
|
2777 next |
|
2778 assume "l1 \<le> 0" "\<not>(l1 = 0 \<and> (u1 = 0 \<or> l2 \<ge> 1))" |
|
2779 with lu show ?thesis by (simp add: bnds_powr_def split: if_split_asm) |
|
2780 next |
|
2781 assume l1: "l1 > 0" |
|
2782 obtain lm um where lmum: |
|
2783 "(lm, um) = bnds_mult prec (the (lb_ln prec l1)) (the (ub_ln prec u1)) l2 u2" |
|
2784 by (cases "bnds_mult prec (the (lb_ln prec l1)) (the (ub_ln prec u1)) l2 u2") simp |
|
2785 with l1 have "(l, u) = map_bnds lb_exp ub_exp prec (lm, um)" |
|
2786 using lu by (simp add: bnds_powr_def del: lb_ln.simps ub_ln.simps split: if_split_asm) |
|
2787 hence "exp (ln x * y) \<in> {real_of_float l..real_of_float u}" |
|
2788 proof (rule map_bnds[OF _ mono_exp_real], goal_cases) |
|
2789 case 1 |
|
2790 let ?lln = "the (lb_ln prec l1)" and ?uln = "the (ub_ln prec u1)" |
|
2791 from ub_ln_lb_ln_bounds[of l1 prec] ub_ln_lb_ln_bounds[of u1 prec] x l1 |
|
2792 have "real_of_float ?lln \<le> ln (real_of_float l1) \<and> |
|
2793 ln (real_of_float u1) \<le> real_of_float ?uln" |
|
2794 by (auto simp del: lb_ln.simps ub_ln.simps) |
|
2795 moreover from l1 x have "ln (real_of_float l1) \<le> ln x \<and> ln x \<le> ln (real_of_float u1)" |
|
2796 by auto |
|
2797 ultimately have ln: "real_of_float ?lln \<le> ln x \<and> ln x \<le> real_of_float ?uln" by simp |
|
2798 from lmum show ?case |
|
2799 by (rule bnds_mult) (insert y ln, simp_all) |
|
2800 qed (insert exp_boundaries[of lm prec] exp_boundaries[of um prec], simp_all) |
|
2801 with x l1 show ?thesis |
|
2802 by (simp add: powr_def mult_ac) |
|
2803 qed |
|
2804 qed |
|
2805 |
|
2806 end |