7 theory ComputeFloat |
7 theory ComputeFloat |
8 imports Complex_Main "~~/src/HOL/Library/Lattice_Algebras" |
8 imports Complex_Main "~~/src/HOL/Library/Lattice_Algebras" |
9 uses "~~/src/Tools/float.ML" ("~~/src/HOL/Tools/float_arith.ML") |
9 uses "~~/src/Tools/float.ML" ("~~/src/HOL/Tools/float_arith.ML") |
10 begin |
10 begin |
11 |
11 |
12 definition pow2 :: "int \<Rightarrow> real" |
|
13 where "pow2 a = (if (0 <= a) then (2^(nat a)) else (inverse (2^(nat (-a)))))" |
|
14 |
|
15 definition float :: "int * int \<Rightarrow> real" |
|
16 where "float x = real (fst x) * pow2 (snd x)" |
|
17 |
|
18 lemma pow2_0[simp]: "pow2 0 = 1" |
|
19 by (simp add: pow2_def) |
|
20 |
|
21 lemma pow2_1[simp]: "pow2 1 = 2" |
|
22 by (simp add: pow2_def) |
|
23 |
|
24 lemma pow2_neg: "pow2 x = inverse (pow2 (-x))" |
|
25 by (simp add: pow2_def) |
|
26 |
|
27 lemma pow2_add1: "pow2 (1 + a) = 2 * (pow2 a)" |
|
28 proof - |
|
29 have h: "! n. nat (2 + int n) - Suc 0 = nat (1 + int n)" by arith |
|
30 have g: "! a b. a - -1 = a + (1::int)" by arith |
|
31 have pos: "! n. pow2 (int n + 1) = 2 * pow2 (int n)" |
|
32 apply (auto, induct_tac n) |
|
33 apply (simp_all add: pow2_def) |
|
34 apply (rule_tac m1="2" and n1="nat (2 + int na)" in ssubst[OF realpow_num_eq_if]) |
|
35 by (auto simp add: h) |
|
36 show ?thesis |
|
37 proof (induct a) |
|
38 case (nonneg n) |
|
39 from pos show ?case by (simp add: algebra_simps) |
|
40 next |
|
41 case (neg n) |
|
42 show ?case |
|
43 apply (auto) |
|
44 apply (subst pow2_neg[of "- int n"]) |
|
45 apply (subst pow2_neg[of "-1 - int n"]) |
|
46 apply (auto simp add: g pos) |
|
47 done |
|
48 qed |
|
49 qed |
|
50 |
|
51 lemma pow2_add: "pow2 (a+b) = (pow2 a) * (pow2 b)" |
|
52 proof (induct b) |
|
53 case (nonneg n) |
|
54 show ?case |
|
55 proof (induct n) |
|
56 case 0 |
|
57 show ?case by simp |
|
58 next |
|
59 case (Suc m) |
|
60 show ?case by (auto simp add: algebra_simps pow2_add1 nonneg Suc) |
|
61 qed |
|
62 next |
|
63 case (neg n) |
|
64 show ?case |
|
65 proof (induct n) |
|
66 case 0 |
|
67 show ?case |
|
68 apply (auto) |
|
69 apply (subst pow2_neg[of "a + -1"]) |
|
70 apply (subst pow2_neg[of "-1"]) |
|
71 apply (simp) |
|
72 apply (insert pow2_add1[of "-a"]) |
|
73 apply (simp add: algebra_simps) |
|
74 apply (subst pow2_neg[of "-a"]) |
|
75 apply (simp) |
|
76 done |
|
77 case (Suc m) |
|
78 have a: "int m - (a + -2) = 1 + (int m - a + 1)" by arith |
|
79 have b: "int m - -2 = 1 + (int m + 1)" by arith |
|
80 show ?case |
|
81 apply (auto) |
|
82 apply (subst pow2_neg[of "a + (-2 - int m)"]) |
|
83 apply (subst pow2_neg[of "-2 - int m"]) |
|
84 apply (auto simp add: algebra_simps) |
|
85 apply (subst a) |
|
86 apply (subst b) |
|
87 apply (simp only: pow2_add1) |
|
88 apply (subst pow2_neg[of "int m - a + 1"]) |
|
89 apply (subst pow2_neg[of "int m + 1"]) |
|
90 apply auto |
|
91 apply (insert Suc) |
|
92 apply (auto simp add: algebra_simps) |
|
93 done |
|
94 qed |
|
95 qed |
|
96 |
|
97 lemma "float (a, e) + float (b, e) = float (a + b, e)" |
|
98 by (simp add: float_def algebra_simps) |
|
99 |
|
100 definition int_of_real :: "real \<Rightarrow> int" |
12 definition int_of_real :: "real \<Rightarrow> int" |
101 where "int_of_real x = (SOME y. real y = x)" |
13 where "int_of_real x = (SOME y. real y = x)" |
102 |
14 |
103 definition real_is_int :: "real \<Rightarrow> bool" |
15 definition real_is_int :: "real \<Rightarrow> bool" |
104 where "real_is_int x = (EX (u::int). x = real u)" |
16 where "real_is_int x = (EX (u::int). x = real u)" |
105 |
17 |
106 lemma real_is_int_def2: "real_is_int x = (x = real (int_of_real x))" |
18 lemma real_is_int_def2: "real_is_int x = (x = real (int_of_real x))" |
107 by (auto simp add: real_is_int_def int_of_real_def) |
19 by (auto simp add: real_is_int_def int_of_real_def) |
108 |
|
109 lemma float_transfer: "real_is_int ((real a)*(pow2 c)) \<Longrightarrow> float (a, b) = float (int_of_real ((real a)*(pow2 c)), b - c)" |
|
110 by (simp add: float_def real_is_int_def2 pow2_add[symmetric]) |
|
111 |
|
112 lemma pow2_int: "pow2 (int c) = 2^c" |
|
113 by (simp add: pow2_def) |
|
114 |
|
115 lemma float_transfer_nat: "float (a, b) = float (a * 2^c, b - int c)" |
|
116 by (simp add: float_def pow2_int[symmetric] pow2_add[symmetric]) |
|
117 |
20 |
118 lemma real_is_int_real[simp]: "real_is_int (real (x::int))" |
21 lemma real_is_int_real[simp]: "real_is_int (real (x::int))" |
119 by (auto simp add: real_is_int_def int_of_real_def) |
22 by (auto simp add: real_is_int_def int_of_real_def) |
120 |
23 |
121 lemma int_of_real_real[simp]: "int_of_real (real x) = x" |
24 lemma int_of_real_real[simp]: "int_of_real (real x) = x" |
180 also have "\<dots> = True" by (simp only: real_is_int_real) |
74 also have "\<dots> = True" by (simp only: real_is_int_real) |
181 ultimately show ?thesis by auto |
75 ultimately show ?thesis by auto |
182 qed |
76 qed |
183 |
77 |
184 lemma real_is_int_number_of[simp]: "real_is_int ((number_of \<Colon> int \<Rightarrow> real) x)" |
78 lemma real_is_int_number_of[simp]: "real_is_int ((number_of \<Colon> int \<Rightarrow> real) x)" |
185 proof - |
79 by (auto simp: real_is_int_def intro!: exI[of _ "number_of x"]) |
186 have neg1: "real_is_int (-1::real)" |
|
187 proof - |
|
188 have "real_is_int (-1::real) = real_is_int(real (-1::int))" by auto |
|
189 also have "\<dots> = True" by (simp only: real_is_int_real) |
|
190 ultimately show ?thesis by auto |
|
191 qed |
|
192 |
|
193 { |
|
194 fix x :: int |
|
195 have "real_is_int ((number_of \<Colon> int \<Rightarrow> real) x)" |
|
196 unfolding number_of_eq |
|
197 apply (induct x) |
|
198 apply (induct_tac n) |
|
199 apply (simp) |
|
200 apply (simp) |
|
201 apply (induct_tac n) |
|
202 apply (simp add: neg1) |
|
203 proof - |
|
204 fix n :: nat |
|
205 assume rn: "(real_is_int (of_int (- (int (Suc n)))))" |
|
206 have s: "-(int (Suc (Suc n))) = -1 + - (int (Suc n))" by simp |
|
207 show "real_is_int (of_int (- (int (Suc (Suc n)))))" |
|
208 apply (simp only: s of_int_add) |
|
209 apply (rule real_is_int_add) |
|
210 apply (simp add: neg1) |
|
211 apply (simp only: rn) |
|
212 done |
|
213 qed |
|
214 } |
|
215 note Abs_Bin = this |
|
216 { |
|
217 fix x :: int |
|
218 have "? u. x = u" |
|
219 apply (rule exI[where x = "x"]) |
|
220 apply (simp) |
|
221 done |
|
222 } |
|
223 then obtain u::int where "x = u" by auto |
|
224 with Abs_Bin show ?thesis by auto |
|
225 qed |
|
226 |
80 |
227 lemma int_of_real_0[simp]: "int_of_real (0::real) = (0::int)" |
81 lemma int_of_real_0[simp]: "int_of_real (0::real) = (0::int)" |
228 by (simp add: int_of_real_def) |
82 by (simp add: int_of_real_def) |
229 |
83 |
230 lemma int_of_real_1[simp]: "int_of_real (1::real) = (1::int)" |
84 lemma int_of_real_1[simp]: "int_of_real (1::real) = (1::int)" |
232 have 1: "(1::real) = real (1::int)" by auto |
86 have 1: "(1::real) = real (1::int)" by auto |
233 show ?thesis by (simp only: 1 int_of_real_real) |
87 show ?thesis by (simp only: 1 int_of_real_real) |
234 qed |
88 qed |
235 |
89 |
236 lemma int_of_real_number_of[simp]: "int_of_real (number_of b) = number_of b" |
90 lemma int_of_real_number_of[simp]: "int_of_real (number_of b) = number_of b" |
237 proof - |
91 unfolding int_of_real_def |
238 have "real_is_int (number_of b)" by simp |
92 by (intro some_equality) |
239 then have uu: "?! u::int. number_of b = real u" by (auto simp add: real_is_int_rep) |
93 (auto simp add: real_of_int_inject[symmetric] simp del: real_of_int_inject) |
240 then obtain u::int where u:"number_of b = real u" by auto |
|
241 have "number_of b = real ((number_of b)::int)" |
|
242 by (simp add: number_of_eq real_of_int_def) |
|
243 have ub: "number_of b = real ((number_of b)::int)" |
|
244 by (simp add: number_of_eq real_of_int_def) |
|
245 from uu u ub have unb: "u = number_of b" |
|
246 by blast |
|
247 have "int_of_real (number_of b) = u" by (simp add: u) |
|
248 with unb show ?thesis by simp |
|
249 qed |
|
250 |
|
251 lemma float_transfer_even: "even a \<Longrightarrow> float (a, b) = float (a div 2, b+1)" |
|
252 apply (subst float_transfer[where a="a" and b="b" and c="-1", simplified]) |
|
253 apply (simp_all add: pow2_def even_def real_is_int_def algebra_simps) |
|
254 apply (auto) |
|
255 proof - |
|
256 fix q::int |
|
257 have a:"b - (-1\<Colon>int) = (1\<Colon>int) + b" by arith |
|
258 show "(float (q, (b - (-1\<Colon>int)))) = (float (q, ((1\<Colon>int) + b)))" |
|
259 by (simp add: a) |
|
260 qed |
|
261 |
94 |
262 lemma int_div_zdiv: "int (a div b) = (int a) div (int b)" |
95 lemma int_div_zdiv: "int (a div b) = (int a) div (int b)" |
263 by (rule zdiv_int) |
96 by (rule zdiv_int) |
264 |
97 |
265 lemma int_mod_zmod: "int (a mod b) = (int a) mod (int b)" |
98 lemma int_mod_zmod: "int (a mod b) = (int a) mod (int b)" |
266 by (rule zmod_int) |
99 by (rule zmod_int) |
267 |
100 |
268 lemma abs_div_2_less: "a \<noteq> 0 \<Longrightarrow> a \<noteq> -1 \<Longrightarrow> abs((a::int) div 2) < abs a" |
101 lemma abs_div_2_less: "a \<noteq> 0 \<Longrightarrow> a \<noteq> -1 \<Longrightarrow> abs((a::int) div 2) < abs a" |
269 by arith |
102 by arith |
270 |
|
271 function norm_float :: "int \<Rightarrow> int \<Rightarrow> int \<times> int" where |
|
272 "norm_float a b = (if a \<noteq> 0 \<and> even a then norm_float (a div 2) (b + 1) |
|
273 else if a = 0 then (0, 0) else (a, b))" |
|
274 by auto |
|
275 |
|
276 termination by (relation "measure (nat o abs o fst)") |
|
277 (auto intro: abs_div_2_less) |
|
278 |
|
279 lemma norm_float: "float x = float (split norm_float x)" |
|
280 proof - |
|
281 { |
|
282 fix a b :: int |
|
283 have norm_float_pair: "float (a, b) = float (norm_float a b)" |
|
284 proof (induct a b rule: norm_float.induct) |
|
285 case (1 u v) |
|
286 show ?case |
|
287 proof cases |
|
288 assume u: "u \<noteq> 0 \<and> even u" |
|
289 with 1 have ind: "float (u div 2, v + 1) = float (norm_float (u div 2) (v + 1))" by auto |
|
290 with u have "float (u,v) = float (u div 2, v+1)" by (simp add: float_transfer_even) |
|
291 then show ?thesis |
|
292 apply (subst norm_float.simps) |
|
293 apply (simp add: ind) |
|
294 done |
|
295 next |
|
296 assume nu: "~(u \<noteq> 0 \<and> even u)" |
|
297 show ?thesis |
|
298 by (simp add: nu float_def) |
|
299 qed |
|
300 qed |
|
301 } |
|
302 note helper = this |
|
303 have "? a b. x = (a,b)" by auto |
|
304 then obtain a b where "x = (a, b)" by blast |
|
305 then show ?thesis by (simp add: helper) |
|
306 qed |
|
307 |
|
308 lemma float_add_l0: "float (0, e) + x = x" |
|
309 by (simp add: float_def) |
|
310 |
|
311 lemma float_add_r0: "x + float (0, e) = x" |
|
312 by (simp add: float_def) |
|
313 |
|
314 lemma float_add: |
|
315 "float (a1, e1) + float (a2, e2) = |
|
316 (if e1<=e2 then float (a1+a2*2^(nat(e2-e1)), e1) |
|
317 else float (a1*2^(nat (e1-e2))+a2, e2))" |
|
318 apply (simp add: float_def algebra_simps) |
|
319 apply (auto simp add: pow2_int[symmetric] pow2_add[symmetric]) |
|
320 done |
|
321 |
|
322 lemma float_add_assoc1: |
|
323 "(x + float (y1, e1)) + float (y2, e2) = (float (y1, e1) + float (y2, e2)) + x" |
|
324 by simp |
|
325 |
|
326 lemma float_add_assoc2: |
|
327 "(float (y1, e1) + x) + float (y2, e2) = (float (y1, e1) + float (y2, e2)) + x" |
|
328 by simp |
|
329 |
|
330 lemma float_add_assoc3: |
|
331 "float (y1, e1) + (x + float (y2, e2)) = (float (y1, e1) + float (y2, e2)) + x" |
|
332 by simp |
|
333 |
|
334 lemma float_add_assoc4: |
|
335 "float (y1, e1) + (float (y2, e2) + x) = (float (y1, e1) + float (y2, e2)) + x" |
|
336 by simp |
|
337 |
|
338 lemma float_mult_l0: "float (0, e) * x = float (0, 0)" |
|
339 by (simp add: float_def) |
|
340 |
|
341 lemma float_mult_r0: "x * float (0, e) = float (0, 0)" |
|
342 by (simp add: float_def) |
|
343 |
|
344 definition lbound :: "real \<Rightarrow> real" |
|
345 where "lbound x = min 0 x" |
|
346 |
|
347 definition ubound :: "real \<Rightarrow> real" |
|
348 where "ubound x = max 0 x" |
|
349 |
|
350 lemma lbound: "lbound x \<le> x" |
|
351 by (simp add: lbound_def) |
|
352 |
|
353 lemma ubound: "x \<le> ubound x" |
|
354 by (simp add: ubound_def) |
|
355 |
|
356 lemma float_mult: |
|
357 "float (a1, e1) * float (a2, e2) = |
|
358 (float (a1 * a2, e1 + e2))" |
|
359 by (simp add: float_def pow2_add) |
|
360 |
|
361 lemma float_minus: |
|
362 "- (float (a,b)) = float (-a, b)" |
|
363 by (simp add: float_def) |
|
364 |
|
365 lemma zero_less_pow2: |
|
366 "0 < pow2 x" |
|
367 proof - |
|
368 { |
|
369 fix y |
|
370 have "0 <= y \<Longrightarrow> 0 < pow2 y" |
|
371 by (induct y, induct_tac n, simp_all add: pow2_add) |
|
372 } |
|
373 note helper=this |
|
374 show ?thesis |
|
375 apply (case_tac "0 <= x") |
|
376 apply (simp add: helper) |
|
377 apply (subst pow2_neg) |
|
378 apply (simp add: helper) |
|
379 done |
|
380 qed |
|
381 |
|
382 lemma zero_le_float: |
|
383 "(0 <= float (a,b)) = (0 <= a)" |
|
384 apply (auto simp add: float_def) |
|
385 apply (auto simp add: zero_le_mult_iff zero_less_pow2) |
|
386 apply (insert zero_less_pow2[of b]) |
|
387 apply (simp_all) |
|
388 done |
|
389 |
|
390 lemma float_le_zero: |
|
391 "(float (a,b) <= 0) = (a <= 0)" |
|
392 apply (auto simp add: float_def) |
|
393 apply (auto simp add: mult_le_0_iff) |
|
394 apply (insert zero_less_pow2[of b]) |
|
395 apply auto |
|
396 done |
|
397 |
|
398 lemma float_abs: |
|
399 "abs (float (a,b)) = (if 0 <= a then (float (a,b)) else (float (-a,b)))" |
|
400 apply (auto simp add: abs_if) |
|
401 apply (simp_all add: zero_le_float[symmetric, of a b] float_minus) |
|
402 done |
|
403 |
|
404 lemma float_zero: |
|
405 "float (0, b) = 0" |
|
406 by (simp add: float_def) |
|
407 |
|
408 lemma float_pprt: |
|
409 "pprt (float (a, b)) = (if 0 <= a then (float (a,b)) else (float (0, b)))" |
|
410 by (auto simp add: zero_le_float float_le_zero float_zero) |
|
411 |
|
412 lemma pprt_lbound: "pprt (lbound x) = float (0, 0)" |
|
413 apply (simp add: float_def) |
|
414 apply (rule pprt_eq_0) |
|
415 apply (simp add: lbound_def) |
|
416 done |
|
417 |
|
418 lemma nprt_ubound: "nprt (ubound x) = float (0, 0)" |
|
419 apply (simp add: float_def) |
|
420 apply (rule nprt_eq_0) |
|
421 apply (simp add: ubound_def) |
|
422 done |
|
423 |
|
424 lemma float_nprt: |
|
425 "nprt (float (a, b)) = (if 0 <= a then (float (0,b)) else (float (a, b)))" |
|
426 by (auto simp add: zero_le_float float_le_zero float_zero) |
|
427 |
103 |
428 lemma norm_0_1: "(0::_::number_ring) = Numeral0 & (1::_::number_ring) = Numeral1" |
104 lemma norm_0_1: "(0::_::number_ring) = Numeral0 & (1::_::number_ring) = Numeral1" |
429 by auto |
105 by auto |
430 |
106 |
431 lemma add_left_zero: "0 + a = (a::'a::comm_monoid_add)" |
107 lemma add_left_zero: "0 + a = (a::'a::comm_monoid_add)" |
547 |
223 |
548 lemmas powerarith = nat_number_of zpower_number_of_even |
224 lemmas powerarith = nat_number_of zpower_number_of_even |
549 zpower_number_of_odd[simplified zero_eq_Numeral0_nring one_eq_Numeral1_nring] |
225 zpower_number_of_odd[simplified zero_eq_Numeral0_nring one_eq_Numeral1_nring] |
550 zpower_Pls zpower_Min |
226 zpower_Pls zpower_Min |
551 |
227 |
|
228 definition float :: "(int \<times> int) \<Rightarrow> real" where |
|
229 "float = (\<lambda>(a, b). real a * 2 powr real b)" |
|
230 |
|
231 lemma float_add_l0: "float (0, e) + x = x" |
|
232 by (simp add: float_def) |
|
233 |
|
234 lemma float_add_r0: "x + float (0, e) = x" |
|
235 by (simp add: float_def) |
|
236 |
|
237 lemma float_add: |
|
238 "float (a1, e1) + float (a2, e2) = |
|
239 (if e1<=e2 then float (a1+a2*2^(nat(e2-e1)), e1) else float (a1*2^(nat (e1-e2))+a2, e2))" |
|
240 by (simp add: float_def algebra_simps powr_realpow[symmetric] powr_divide2[symmetric]) |
|
241 |
|
242 lemma float_mult_l0: "float (0, e) * x = float (0, 0)" |
|
243 by (simp add: float_def) |
|
244 |
|
245 lemma float_mult_r0: "x * float (0, e) = float (0, 0)" |
|
246 by (simp add: float_def) |
|
247 |
|
248 lemma float_mult: |
|
249 "float (a1, e1) * float (a2, e2) = (float (a1 * a2, e1 + e2))" |
|
250 by (simp add: float_def powr_add) |
|
251 |
|
252 lemma float_minus: |
|
253 "- (float (a,b)) = float (-a, b)" |
|
254 by (simp add: float_def) |
|
255 |
|
256 lemma zero_le_float: |
|
257 "(0 <= float (a,b)) = (0 <= a)" |
|
258 using powr_gt_zero[of 2 "real b", arith] |
|
259 by (simp add: float_def zero_le_mult_iff) |
|
260 |
|
261 lemma float_le_zero: |
|
262 "(float (a,b) <= 0) = (a <= 0)" |
|
263 using powr_gt_zero[of 2 "real b", arith] |
|
264 by (simp add: float_def mult_le_0_iff) |
|
265 |
|
266 lemma float_abs: |
|
267 "abs (float (a,b)) = (if 0 <= a then (float (a,b)) else (float (-a,b)))" |
|
268 using powr_gt_zero[of 2 "real b", arith] |
|
269 by (simp add: float_def abs_if mult_less_0_iff) |
|
270 |
|
271 lemma float_zero: |
|
272 "float (0, b) = 0" |
|
273 by (simp add: float_def) |
|
274 |
|
275 lemma float_pprt: |
|
276 "pprt (float (a, b)) = (if 0 <= a then (float (a,b)) else (float (0, b)))" |
|
277 by (auto simp add: zero_le_float float_le_zero float_zero) |
|
278 |
|
279 lemma float_nprt: |
|
280 "nprt (float (a, b)) = (if 0 <= a then (float (0,b)) else (float (a, b)))" |
|
281 by (auto simp add: zero_le_float float_le_zero float_zero) |
|
282 |
|
283 definition lbound :: "real \<Rightarrow> real" |
|
284 where "lbound x = min 0 x" |
|
285 |
|
286 definition ubound :: "real \<Rightarrow> real" |
|
287 where "ubound x = max 0 x" |
|
288 |
|
289 lemma lbound: "lbound x \<le> x" |
|
290 by (simp add: lbound_def) |
|
291 |
|
292 lemma ubound: "x \<le> ubound x" |
|
293 by (simp add: ubound_def) |
|
294 |
|
295 lemma pprt_lbound: "pprt (lbound x) = float (0, 0)" |
|
296 by (auto simp: float_def lbound_def) |
|
297 |
|
298 lemma nprt_ubound: "nprt (ubound x) = float (0, 0)" |
|
299 by (auto simp: float_def ubound_def) |
|
300 |
552 lemmas floatarith[simplified norm_0_1] = float_add float_add_l0 float_add_r0 float_mult float_mult_l0 float_mult_r0 |
301 lemmas floatarith[simplified norm_0_1] = float_add float_add_l0 float_add_r0 float_mult float_mult_l0 float_mult_r0 |
553 float_minus float_abs zero_le_float float_pprt float_nprt pprt_lbound nprt_ubound |
302 float_minus float_abs zero_le_float float_pprt float_nprt pprt_lbound nprt_ubound |
554 |
303 |
555 (* for use with the compute oracle *) |
304 (* for use with the compute oracle *) |
556 lemmas arith = binarith intarith intarithrel natarith powerarith floatarith not_false_eq_true not_true_eq_false |
305 lemmas arith = binarith intarith intarithrel natarith powerarith floatarith not_false_eq_true not_true_eq_false |