2 needs "jordan/refinement.hl";;
4 needs "jordan/parse_ext_override_interface.hl";;
5 needs "jordan/real_ext.hl";;
7 needs "jordan/taylor_atn.hl";;
9 needs "../formal_lp/arith/float.hl";;
12 module type Float_atn_sig =
14 val float_interval_pow_simple : int -> int -> thm -> thm
15 val pi_approx_array : thm array
16 val pi2_approx_array : thm array
18 val float_interval_atn : int -> thm -> thm
19 val float_interval_acs : int -> thm -> thm
23 module Float_atn : Float_atn_sig = struct
32 (******************************)
33 let x_var_real = `x:real` and
34 n_var_num = `n:num` and
35 e_var_num = `e:num` and
36 a_var_real = `a:real` and
37 b_var_real = `b:real` and
38 d_var_real = `d:real` and
39 hi_var_real = `hi:real` and
40 lo_var_real = `lo:real`;;
43 let add_op_real = `(+):real->real->real` and
44 sub_op_real = `(-):real->real->real` and
45 mul_op_real = `( * ):real->real->real` and
46 div_op_real = `(/):real->real->real` and
47 neg_op_real = `(--):real->real` and
48 mul_op_num = `( * ):num->num->num` and
49 add_op_num = `(+):num->num->num`;;
53 (******************************)
54 (* halfatn and halfatn4 *)
56 let float_interval_1 = mk_float_interval_small_num 1;;
59 let HALFATN' = (SYM o SPEC_ALL o REWRITE_RULE[REAL_POW_2]) halfatn;;
60 let HALFATN4' = prove(`halfatn(halfatn(halfatn(halfatn x))) = halfatn4 x`,
61 REWRITE_TAC[halfatn4; o_THM]);;
66 let float_interval_halfatn pp x_th =
67 let x_tm = (rand o rator o concl) x_th in
68 let xx_th = float_interval_mul pp x_th x_th in
69 let one_xx_th = float_interval_add pp float_interval_1 xx_th in
70 let sqrt_th = float_interval_sqrt pp one_xx_th in
71 let one_sqrt_th = float_interval_add pp sqrt_th float_interval_1 in
72 let r_th = float_interval_div pp x_th one_sqrt_th in
73 let th0 = INST[x_tm, x_var_real] HALFATN' in
74 let ltm, rtm = dest_comb(concl r_th) in
75 EQ_MP (AP_THM (AP_TERM (rator ltm) th0) rtm) r_th;;
78 let float_interval_halfatn4 pp x_th =
79 let x_tm = (rand o rator o concl) x_th in
80 let r_th = float_interval_halfatn pp
81 (float_interval_halfatn pp
82 (float_interval_halfatn pp (float_interval_halfatn pp x_th))) in
83 let th0 = INST[x_tm, x_var_real] HALFATN4' in
84 let ltm, rtm = dest_comb(concl r_th) in
85 EQ_MP (AP_THM (AP_TERM (rator ltm) th0) rtm) r_th;;
90 let x_th = mk_float_interval_small_num 3;;
91 let x_th2 = float_interval_div 11 float_interval_1 x_th;;
93 float_interval_halfatn4 11 x_th2;;
95 test 10 (float_interval_halfatn4 11) x_th2;;
100 (****************************************)
103 let rec float_interval_calc pp expr x_th =
107 let ltm, r_tm = dest_comb expr in
109 let op, l_tm = dest_comb ltm in
110 let l_th = float_interval_calc pp l_tm x_th in
111 let r_th = float_interval_calc pp r_tm x_th in
112 if op = add_op_real then
113 float_interval_add pp l_th r_th
114 else if op = mul_op_real then
115 float_interval_mul pp l_th r_th
116 else if op = div_op_real then
117 float_interval_div pp l_th r_th
118 else if op = sub_op_real then
119 float_interval_sub pp l_th r_th
121 failwith ("Unknown operation: " ^ (fst o dest_const) op)
123 if ltm = neg_op_real then
124 let r_th = float_interval_calc pp r_tm x_th in
125 float_interval_neg r_th
127 mk_float_interval_num (dest_numeral r_tm);;
130 (*************************************)
131 (* Polynomial functions *)
133 let poly_f = new_definition `poly_f cs x = ITLIST (\c s. c + x * s) cs (&0)`;;
136 let poly_f_even = new_definition `poly_f_even cs x = ITLIST (\c s. c + (x * x) * s) cs (&0)`;;
138 let poly_f_odd = new_definition `poly_f_odd cs x = x * poly_f_even cs x`;;
139 let poly_f_odd' = SPECL[`t:(real)list`; `x:real`] poly_f_odd;;
141 let NUMERALS_TO_NUM = Arith_nat.NUMERALS_TO_NUM;;
143 let POLY_F_EMPTY = (NUMERALS_TO_NUM o prove) (`poly_f [] x = &0`, REWRITE_TAC[poly_f; ITLIST]) and
144 POLY_F_CONS = prove(`poly_f (CONS h t) x = h + x * poly_f t x`, REWRITE_TAC[poly_f; ITLIST]);;
146 let POLY_F_EVEN_EMPTY = (NUMERALS_TO_NUM o prove) (`poly_f_even [] x = &0`, REWRITE_TAC[poly_f_even; ITLIST]) and
147 POLY_F_EVEN_CONS = prove(`poly_f_even (CONS h t) x = h + (x * x) * poly_f_even t x`, REWRITE_TAC[poly_f_even; ITLIST]);;
149 let POLY_F_ODD_EMPTY = (NUMERALS_TO_NUM o prove) (`poly_f_odd [] x = &0`, REWRITE_TAC[poly_f_odd; poly_f_even; ITLIST; REAL_MUL_RZERO]);;
154 let REVERSE_TABLE = define `(REVERSE_TABLE (f:num->A) 0 = []) /\
155 (REVERSE_TABLE f (SUC i) = CONS (f i) ( REVERSE_TABLE f i))`;;
158 let TABLE = new_definition `!(f:num->A) k. TABLE f k = REVERSE (REVERSE_TABLE f k)`;;
161 let LENGTH_REVERSE_TABLE = prove(`!(f:num->A) n. LENGTH (REVERSE_TABLE f n) = n`,
162 GEN_TAC THEN INDUCT_TAC THEN ASM_REWRITE_TAC[REVERSE_TABLE; LENGTH]);;
165 let LENGTH_REVERSE = prove(`!(l:(A)list). LENGTH (REVERSE l) = LENGTH l`,
166 MATCH_MP_TAC list_INDUCT THEN REWRITE_TAC[REVERSE] THEN
167 REPEAT STRIP_TAC THEN
168 ASM_REWRITE_TAC[LENGTH_APPEND; LENGTH] THEN
172 let LENGTH_TABLE = prove(`!(f:num->A) n. LENGTH (TABLE f n) = n`,
173 REWRITE_TAC[TABLE; LENGTH_REVERSE; LENGTH_REVERSE_TABLE]);;
177 let rec reverse_table_conv tm =
178 let ltm, i_tm = dest_comb tm in
180 ONCE_REWRITE_CONV[REVERSE_TABLE] tm
182 let i_suc = num_CONV i_tm in
183 let th1 = ONCE_REWRITE_RULE[REVERSE_TABLE] (AP_TERM ltm i_suc) in
184 let ltm, rtm = dest_comb (rand(concl th1)) in
185 let th2 = reverse_table_conv rtm in
186 TRANS th1 (AP_TERM ltm th2);;
191 let atn_co_table = new_definition `atn_co_table n = TABLE
192 (\k. (if (EVEN k) then &1 else --(&1)) / &(2 * k + 1)) (SUC n)`;;
195 (* Returns a theorem |- atn_co_table n = [...] and
196 a list of interval approximations of the coefficients in the table *)
197 let mk_atn_co_table pp n =
198 let table = SPEC (mk_small_numeral n) atn_co_table in
199 let th = CONV_RULE (DEPTH_CONV NUM_SUC_CONV THENC
200 REWRITE_CONV[TABLE] THENC
201 ONCE_DEPTH_CONV reverse_table_conv THENC
202 REWRITE_CONV[REVERSE; APPEND] THENC
203 NUM_REDUCE_CONV) table in
204 let list = (rand o concl) th in
205 th, map (fun tm -> float_interval_calc pp tm float_interval_1) (dest_list list);;
209 let POLY_F_EVEN_ALT = prove(`poly_f_even cs x = poly_f cs (x * x)`,
210 REWRITE_TAC[poly_f_even; poly_f]);;
213 let POLY_F_APPEND = prove(`!x b a. poly_f (APPEND a b) x = poly_f a x + x pow (LENGTH a) * poly_f b x`,
214 GEN_TAC THEN GEN_TAC THEN
215 MATCH_MP_TAC list_INDUCT THEN
216 REPEAT STRIP_TAC THENL
218 REWRITE_TAC[APPEND; poly_f; ITLIST; LENGTH] THEN
219 REWRITE_TAC[real_pow; REAL_MUL_LID; REAL_ADD_LID];
223 REWRITE_TAC[APPEND; poly_f; ITLIST] THEN
224 ASM_REWRITE_TAC[GSYM poly_f] THEN
225 REWRITE_TAC[LENGTH; real_pow] THEN
229 let POLY_F_EVEN_APPEND = prove(`!x b a. poly_f_even (APPEND a b) x = poly_f_even a x + x pow (2 * LENGTH a) * poly_f_even b x`,
230 REWRITE_TAC[POLY_F_EVEN_ALT; POLY_F_APPEND] THEN
231 REWRITE_TAC[GSYM REAL_POW_2; REAL_POW_POW]);;
234 let POLY_F_ODD_APPEND = prove(`!x b a. poly_f_odd (APPEND a b) x = poly_f_odd a x + x pow (2 * LENGTH a) * poly_f_odd b x`,
236 REWRITE_TAC[poly_f_odd] THEN
237 REWRITE_TAC[POLY_F_EVEN_APPEND] THEN
241 let ATN_SUM_TABLE = prove(`!x n. sum (0..n) (halfatn4_co x) = poly_f_odd (atn_co_table n) (halfatn4 x)`,
242 GEN_TAC THEN INDUCT_TAC THENL
244 REWRITE_TAC[SUM_SING_NUMSEG; atn_co_table; TABLE; REVERSE_TABLE; REVERSE; APPEND] THEN
245 REWRITE_TAC[ARITH_EVEN] THEN
246 REWRITE_TAC[poly_f_odd; poly_f_even; ITLIST] THEN
247 REWRITE_TAC[halfatn4_co; REAL_MUL_RZERO; REAL_ADD_RID] THEN
248 REWRITE_TAC[MULT_CLAUSES; ARITH_ADD; REAL_POW_1; real_pow] THEN
253 REWRITE_TAC[SUM_CLAUSES_NUMSEG; atn_co_table; TABLE; LE_0] THEN
254 ONCE_REWRITE_TAC[REVERSE_TABLE] THEN
255 ONCE_REWRITE_TAC[REVERSE] THEN
256 REWRITE_TAC[GSYM TABLE; GSYM atn_co_table] THEN
257 ASM_REWRITE_TAC[POLY_F_ODD_APPEND] THEN
258 REWRITE_TAC[REAL_EQ_ADD_LCANCEL] THEN
259 REWRITE_TAC[atn_co_table; LENGTH_TABLE; halfatn4_co] THEN
260 REWRITE_TAC[poly_f_odd; poly_f_even; ITLIST; REAL_MUL_RZERO; REAL_ADD_RID] THEN
261 REWRITE_TAC[real_div; REAL_MUL_ASSOC] THEN
262 REWRITE_TAC[REAL_EQ_MUL_RCANCEL] THEN
264 REWRITE_TAC[REAL_MUL_AC] THEN
265 REWRITE_TAC[GSYM real_pow; ARITH_RULE `2 * SUC n + 1 = SUC (2 * SUC n)`] THEN
266 GEN_REWRITE_TAC LAND_CONV [REAL_MUL_AC] THEN
267 REWRITE_TAC[REAL_EQ_MUL_RCANCEL] THEN
269 REWRITE_TAC[REAL_POW_NEG; real_pow; REAL_POW_ONE; REAL_MUL_RID]);;
273 let POLY_F_SING = prove(`poly_f [c] x = c`,
274 REWRITE_TAC[poly_f; ITLIST; REAL_MUL_RZERO; REAL_ADD_RID]);;
277 let c_var_real = `c:real` and
278 cs_var_list = `cs:(real)list` and
279 h_var_real = `h:real` and
280 t_var_list = `t:(real)list`;;
282 let interval_const = `interval_arith`;;
285 let rec float_interval_poly_f pp (cs, l) x_th =
287 failwith "float_interval_poly_f: an empty coefficient list"
289 let ltm, x_bounds = dest_comb (concl x_th) in
290 let x_tm = rand ltm in
292 let ltm, first_bounds = dest_comb (concl first) in
293 let first_tm = rand ltm in
295 let th0 = INST[first_tm, c_var_real; x_tm, x_var_real] POLY_F_SING in
296 EQ_MP (SYM (AP_THM (AP_TERM interval_const th0) first_bounds)) first
298 let ltm, t_tm = dest_comb cs in
299 let h_tm = rand ltm in
301 let th0 = INST[h_tm, h_var_real; t_tm, t_var_list; x_tm, x_var_real] POLY_F_CONS in
302 let r_th = float_interval_poly_f pp (t_tm, tl l) x_th in
303 let th1 = float_interval_add pp first (float_interval_mul pp x_th r_th) in
304 let bounds = rand (concl th1) in
305 EQ_MP (SYM (AP_THM (AP_TERM interval_const th0) bounds)) th1;;
311 let float_interval_poly_f_even pp (cs, l) x_th =
312 let x_tm = (rand o rator o concl) x_th in
313 let xx_th = float_interval_mul pp x_th x_th in
314 let th0 = INST[cs, cs_var_list; x_tm, x_var_real] POLY_F_EVEN_ALT in
315 let th1 = float_interval_poly_f pp (cs, l) xx_th in
316 let bounds = rand(concl th1) in
317 EQ_MP (SYM (AP_THM (AP_TERM interval_const th0) bounds)) th1;;
320 let float_interval_poly_f_odd pp (cs, l) x_th =
321 let x_tm = (rand o rator o concl) x_th in
322 let th0 = INST[cs, t_var_list; x_tm, x_var_real] poly_f_odd' in
323 let even_th = float_interval_poly_f_even pp (cs, l) x_th in
324 let th1 = float_interval_mul pp x_th even_th in
325 let bounds = rand(concl th1) in
326 EQ_MP (SYM (AP_THM (AP_TERM interval_const th0) bounds)) th1;;
331 let poly_f_odd_const = `poly_f_odd`;;
332 let ATN_SUM_TABLE' = SPEC_ALL ATN_SUM_TABLE;;
333 let float_interval_16 = mk_float_interval_small_num 16;;
336 (* Computes an interval for &16 * sum(0..n) (halfatn4_co x) *)
337 let float_interval_atn_sum pp (cs_th, l) x_th =
338 let n_tm = (rand o lhand o concl) cs_th in
339 let cs_tm = rand(concl cs_th) in
340 let halfatn4 = float_interval_halfatn4 pp x_th in
342 let poly_th = float_interval_poly_f_odd pp (cs_tm, l) halfatn4 in
343 let bounds = rand (concl poly_th) in
344 let halfatn4_tm = (rand o rator o concl) halfatn4 in
345 let x_tm = rand halfatn4_tm in
347 let th1 = AP_THM (AP_TERM interval_const (AP_THM (AP_TERM poly_f_odd_const cs_th) halfatn4_tm)) bounds in
348 let poly_atn_th = EQ_MP (SYM th1) poly_th in
349 let bounds = rand (concl poly_atn_th) in
351 let th2 = INST[n_tm, n_var_num; x_tm, x_var_real] ATN_SUM_TABLE' in
352 let th3 = EQ_MP (SYM (AP_THM (AP_TERM interval_const th2) bounds)) poly_atn_th in
353 float_interval_mul pp float_interval_16 th3;;
360 let cs_th, l = mk_atn_co_table pp 4;;
361 let cs = rand (concl cs_th);;
362 let x_th = float_interval_1;;
364 let th = float_interval_poly_f pp (cs, l) x_th;;
365 float_interval_mul pp (mk_float_interval_small_num 4) th;;
368 let th = float_interval_atn_sum pp (cs_th, l) x_th;;
369 float_interval_mul pp (mk_float_interval_small_num 4) th;;
373 (******************************)
375 let bounds_var_pair = `bounds:real#real`;;
378 let FLOAT_INTERVAL_INV = prove(`interval_arith (&1 / x) bounds <=>
379 interval_arith (inv x) bounds`,
380 REWRITE_TAC[real_div; REAL_MUL_LID]);;
384 let float_interval_inv pp x_th =
385 let x_tm = (rand o rator o concl) x_th in
386 let r_th = float_interval_div pp float_interval_1 x_th in
387 let th0 = INST[x_tm, x_var_real; rand(concl r_th), bounds_var_pair] FLOAT_INTERVAL_INV in
391 let REAL_POW_SUC = (SPEC_ALL o CONJUNCT2) real_pow;;
393 let INTERVAL_REAL_POW_0 = prove(mk_comb(`interval_arith (x pow 0)`, (rand o concl) float_interval_1),
394 REWRITE_TAC[real_pow; float_interval_1]);;
396 let INTERVAL_REAL_POW_1 = prove(`interval_arith x bounds <=> interval_arith (x pow 1) bounds`,
397 REWRITE_TAC[REAL_POW_1]);;
402 let rec float_interval_pow_simple pp n x_th =
403 let x_tm = (rand o rator o concl) x_th in
405 INST[x_tm, x_var_real] INTERVAL_REAL_POW_0
408 let bounds = rand(concl x_th) in
409 let th0 = INST[x_tm, x_var_real; bounds, bounds_var_pair] INTERVAL_REAL_POW_1 in
412 let n_tm' = mk_small_numeral n in
413 let n_suc = num_CONV n_tm' in
414 let n_tm = rand(rand(concl n_suc)) in
415 let th0 = INST[x_tm, x_var_real; n_tm, n_var_num] REAL_POW_SUC in
416 let r_th = float_interval_pow_simple pp (n - 1) x_th in
417 let th1 = float_interval_mul pp x_th r_th in
418 let bounds = rand (concl th1) in
419 let th2 = TRANS (AP_TERM (rator(lhand(concl th0))) n_suc) th0 in
420 EQ_MP (SYM (AP_THM (AP_TERM interval_const th2) bounds)) th1;;
424 let float_interval_2 = mk_float_interval_small_num 2 and
430 (* Computes an interval for inv(&2 pow (6 * n + 5)) *)
431 let compute_eps1 pp n =
432 let n_tm = mk_small_numeral n in
433 let n6 = NUM_MULT_CONV (mk_binop mul_op_num six_const n_tm) in
434 let n65_1 = AP_THM (AP_TERM add_op_num n6) five_const in
435 let n65_2 = NUM_ADD_CONV (rand (concl n65_1)) in
436 let n65 = TRANS n65_1 n65_2 in
437 let pow_th = float_interval_pow_simple pp (6 * n + 5) float_interval_2 in
438 let ltm, bounds = dest_comb(concl pow_th) in
439 let pow_tm = (rator o rand) ltm in
440 let th0 = EQ_MP (SYM (AP_THM (AP_TERM interval_const (AP_TERM pow_tm n65)) bounds)) pow_th in
441 float_interval_inv pp th0;;
446 (**********************************)
450 let FLOAT_ATN_LO_HI = prove(`interval_arith (&16 * sum(0..n) (halfatn4_co x)) (a, b) /\
451 interval_arith (inv(&2 pow (6*n + 5))) (c,d) /\
452 b + d <= hi /\ lo <= a - d
453 ==> interval_arith (atn x) (lo, hi)`,
454 REWRITE_TAC[interval_arith] THEN
456 MP_TAC (SPEC_ALL real_taylor_atn_halfatn4) THEN
457 MP_TAC (REAL_ARITH `&0 <= abs(&16)`) THEN
458 REWRITE_TAC[IMP_IMP] THEN
459 DISCH_THEN (MP_TAC o MATCH_MP REAL_LE_LMUL) THEN
460 REWRITE_TAC[GSYM REAL_ABS_MUL; REAL_ARITH `a * (b - c) = a * b - a * c:real`] THEN
461 ONCE_REWRITE_TAC[GSYM atn_halfatn4] THEN
462 REWRITE_TAC[REAL_ARITH `abs (x - v) <= e <=> v - e <= x /\ x <= v + e`] THEN
463 REWRITE_TAC[REAL_ABS_NUM] THEN
464 SUBGOAL_THEN `&16 * inv(&8 pow (2 * n + 3)) = inv(&2 pow (6 * n + 5))` (fun th -> REWRITE_TAC[th]) THENL
466 REWRITE_TAC[GSYM real_div] THEN
467 SUBGOAL_THEN `&16 = &2 pow 4 /\ &8 = &2 pow 3 /\ ~(&2 = &0)` ASSUME_TAC THENL
472 ASM_REWRITE_TAC[REAL_POW_POW] THEN
473 ASM_SIMP_TAC[REAL_DIV_POW2] THEN
474 REWRITE_TAC[ARITH_RULE `~(3 * (2 * n + 3) <= 4)`] THEN
475 REPEAT AP_TERM_TAC THEN
479 REWRITE_TAC[GSYM halfatn4_co] THEN
480 SUBGOAL_THEN `sum (0..n) (\j. halfatn4_co x j) = sum (0..n) (halfatn4_co x)` (fun th -> REWRITE_TAC[th]) THENL
483 REWRITE_TAC[FUN_EQ_THM];
486 REPEAT (POP_ASSUM MP_TAC) THEN
490 let FLOAT_ATN_LO_HI' = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP]) FLOAT_ATN_LO_HI;;
494 let float_interval_atn_0 pp (cs_th, l) eps1_th x_th =
495 let sum_th = float_interval_atn_sum pp (cs_th, l) x_th in
496 let n_tm = (rand o lhand o concl) cs_th in
497 let x_tm = (rand o rator o concl) x_th in
499 let sum_bounds = rand (concl sum_th) in
500 let a_tm, b_tm = dest_pair sum_bounds in
501 let c_tm, d_tm = (dest_pair o rand o concl) eps1_th in
503 let hi_th = float_add_hi pp b_tm d_tm in
504 let lo_th = float_sub_lo pp a_tm d_tm in
505 let hi_tm = rand(concl hi_th) in
506 let lo_tm = lhand(concl lo_th) in
508 let th0 = INST[n_tm, n_var_num; x_tm, x_var_real;
509 a_tm, a_var_real; b_tm, b_var_real;
510 c_tm, c_var_real; d_tm, d_var_real;
511 hi_tm, hi_var_real; lo_tm, lo_var_real] FLOAT_ATN_LO_HI' in
512 MY_PROVE_HYP lo_th (MY_PROVE_HYP hi_th (MY_PROVE_HYP sum_th (MY_PROVE_HYP eps1_th th0)));;
518 let cs_th, l = mk_atn_co_table pp n;;
519 let eps1_th = compute_eps1 pp n;;
520 let x_th = float_interval_2;;
522 float_interval_atn_0 pp (cs_th, l) eps1_th x_th;;
523 (* 10(min_exp = 20): 1.316 *)
524 test 10 (float_interval_atn_0 pp (cs_th, l) eps1_th) x_th;;
529 (* Fill in lookup tables *)
531 (* Computes n such that 2^(-(6n + 5)) <= base^(-(p + 1)) *)
533 let x = (float_of_int (pp + 1) *. log (float_of_int Arith_options.base) /. log (2.0) -. 5.0) /. 6.0 in
534 let n = (int_of_float o ceil) x in
535 if n < 1 then 1 else n;;
538 let atn_co_array = Array.init 21 (fun i -> mk_atn_co_table (i + 1) (n_of_p i));;
539 let eps1_array = Array.init 21 (fun i -> compute_eps1 (i + 1) (n_of_p i));;
542 let float_interval_atn pp x_th =
543 float_interval_atn_0 pp atn_co_array.(pp) eps1_array.(pp) x_th;;
548 (*****************************************)
550 (* pi approximation *)
553 let x_th = float_interval_1;;
554 let th1 = float_interval_atn pp x_th;;
555 let th2 = float_interval_mul pp (mk_float_interval_small_num 4) th1;;
556 let float_interval_pi = REWRITE_RULE[ATN_1; REAL_ARITH `&4 * pi / &4 = pi`] th2;;
557 let float_interval_pi2 = float_interval_div pp float_interval_pi float_interval_2;;
560 let pi_approx_array = Array.init 19 (fun i -> float_interval_round i float_interval_pi);;
561 let pi2_approx_array = Array.init 19 (fun i -> float_interval_round i float_interval_pi2);;
566 (********************************************)
570 let TAN_HALF = prove(`!x. ~(cos x = -- &1) ==> tan (x / &2) = sin x / (&1 + cos x)`,
572 ABBREV_TAC `t = x / &2` THEN
573 SUBGOAL_THEN `x = &2 * t` ASSUME_TAC THENL
575 EXPAND_TAC "t" THEN REAL_ARITH_TAC;
579 ASM_REWRITE_TAC[SIN_DOUBLE; COS_DOUBLE_COS; REAL_ARITH `&1 + a - &1 = a`] THEN
580 REWRITE_TAC[REAL_ARITH `a - &1 = -- &1 <=> a = &0`] THEN
581 REWRITE_TAC[real_div; REAL_INV_MUL; REAL_POW_2] THEN
582 REWRITE_TAC[REAL_ENTIRE; REAL_ARITH `&2 = &0 <=> F`] THEN
584 REWRITE_TAC[REAL_ARITH `(&2 * s * c) * i2 * ic * ic = (&2 * i2) * (c * ic) * s * ic`] THEN
585 ASM_SIMP_TAC[REAL_MUL_RINV; REAL_ARITH `~(&2 = &0)`] THEN
586 REWRITE_TAC[REAL_MUL_LID; tan; real_div]);;
590 let X_EQ_COS_T = prove(`!x. abs x <= &1 ==> ?t. &0 <= t /\ t <= pi /\ x = cos t`,
591 REWRITE_TAC[REAL_ARITH `abs x <= &1 <=> -- &1 <= x /\ x <= &1`] THEN
592 REPEAT STRIP_TAC THEN
593 EXISTS_TAC `acs x` THEN
594 ASM_SIMP_TAC[ACS_BOUNDS; COS_ACS]);;
598 let ACS_ATN_ALT = prove(`!x. -- &1 < x /\ x <= &1 ==>
599 acs x = &2 * atn (sqrt (&1 - x pow 2) / (&1 + x))`,
600 REPEAT STRIP_TAC THEN
601 MP_TAC (SPEC_ALL X_EQ_COS_T) THEN
604 REPEAT (POP_ASSUM MP_TAC) THEN
608 REPEAT STRIP_TAC THEN
609 ASM_REWRITE_TAC[] THEN
610 ASM_SIMP_TAC[ACS_COS] THEN
611 MP_TAC (SPEC `t:real` SIN_COS_SQRT) THEN
614 ASM_SIMP_TAC[SIN_POS_PI_LE];
618 DISCH_THEN (fun th -> REWRITE_TAC[SYM th]) THEN
619 MP_TAC (SPEC `t:real` TAN_HALF) THEN
622 POP_ASSUM (fun th -> REWRITE_TAC[SYM th]) THEN
623 UNDISCH_TAC `-- &1 < x` THEN
628 DISCH_THEN (fun th -> REWRITE_TAC[SYM th]) THEN
629 REWRITE_TAC[REAL_ARITH `t = &2 * a <=> a = t / &2`] THEN
630 MATCH_MP_TAC TAN_ATN THEN
631 REWRITE_TAC[REAL_ARITH `a / &2 < b / &2 <=> a < b`] THEN
632 REWRITE_TAC[REAL_ARITH `--(a / &2) < b / &2 <=> --a < b`] THEN
635 MATCH_MP_TAC REAL_LTE_TRANS THEN
637 ASM_REWRITE_TAC[REAL_NEG_LT0; PI_POS];
638 SUBGOAL_THEN `t = acs x` MP_TAC THENL
640 ASM_SIMP_TAC[ACS_COS];
643 DISCH_THEN (fun th -> REWRITE_TAC[th]) THEN
644 REWRITE_TAC[SYM ACS_NEG_1] THEN
645 MATCH_MP_TAC ACS_MONO_LT THEN
646 ASM_REWRITE_TAC[REAL_LE_REFL]
651 let FLOAT_F_LT = prove(`!n e. &0 < float F n e <=> ~(n = 0)`,
652 REWRITE_TAC[REAL_ARITH `&0 < a <=> &0 <= a /\ ~(a = &0)`] THEN
653 REWRITE_TAC[FLOAT_F_POS; FLOAT_EQ_0]);;
657 let FLOAT_INTERVAL_ACS = prove(`interval_arith (pi / &2 - atn(x / sqrt(&1 - x * x))) bounds /\
658 interval_arith (&1 - x * x) (float F n e, hi) /\
660 ==> interval_arith (acs x) bounds`,
661 REWRITE_TAC[GSYM REAL_POW_2] THEN
663 MP_TAC (SPEC_ALL ACS_ATN) THEN
666 POP_ASSUM MP_TAC THEN POP_ASSUM MP_TAC THEN
667 REWRITE_TAC[interval_arith] THEN
668 REWRITE_TAC[REAL_ARITH `-- &1 < x /\ x < &1 <=> abs x < abs (&1)`] THEN
669 REWRITE_TAC[REAL_LT_SQUARE_ABS] THEN
670 REWRITE_TAC[REAL_ARITH `&1 pow 2 = &1`] THEN
671 REPEAT STRIP_TAC THEN
672 REWRITE_TAC[REAL_ARITH `a < &1 <=> &0 < &1 - a`] THEN
674 MP_TAC (SPEC_ALL FLOAT_F_LT) THEN
675 ASM_REWRITE_TAC[] THEN
676 UNDISCH_TAC `float F n e <= &1 - x pow 2` THEN
681 DISCH_THEN (fun th -> ASM_REWRITE_TAC[th]));;
686 let ZERO_EQ_ZERO_CONST = prove(`0 = _0`, REWRITE_TAC[NUMERAL]);;
689 let FLOAT_INTERVAL_ACS' = (UNDISCH_ALL o
690 PURE_ONCE_REWRITE_RULE[TAUT `~P <=> (P <=> F)`] o
691 REWRITE_RULE[ZERO_EQ_ZERO_CONST; GSYM IMP_IMP]) FLOAT_INTERVAL_ACS;;
696 let float_interval_acs_0 pp (cs_th, l) eps1_th x_th =
697 let int1 = float_interval_sub pp float_interval_1 (float_interval_mul pp x_th x_th) in
698 let int2 = float_interval_div pp x_th (float_interval_sqrt pp int1) in
699 let atn_int = float_interval_atn_0 pp (cs_th, l) eps1_th int2 in
700 let acs_int = float_interval_sub pp pi2_approx_array.(pp + 1) atn_int in
702 let x_tm = (rand o rator o concl) x_th in
703 let bounds = (rand o concl) acs_int in
704 let int1_bounds = (rand o concl) int1 in
705 let lo_tm, hi_tm = dest_pair int1_bounds in
706 let s, n_tm, e_tm = dest_float lo_tm in
708 failwith "float_interval_acs_0: &1 - x pow 2 < &1 is not satisfied"
710 let n_th = Arith_nat.raw_eq0_hash_conv n_tm in
711 let th0 = INST[x_tm, x_var_real; bounds, bounds_var_pair;
712 n_tm, n_var_num; e_tm, e_var_num;
713 hi_tm, hi_var_real] FLOAT_INTERVAL_ACS' in
714 MY_PROVE_HYP acs_int (MY_PROVE_HYP int1 (MY_PROVE_HYP n_th th0));;
717 let float_interval_acs pp x_th =
718 float_interval_acs_0 pp atn_co_array.(pp) eps1_array.(pp) x_th;;
722 (* acs(&1 / &3) <= #1.230959418 *)
725 let cs_th, l = mk_atn_co_table pp n;;
726 let eps1_th = compute_eps1 pp n;;
727 let x_th = float_interval_div pp float_interval_1 (mk_float_interval_small_num 3);;
730 float_interval_acs_0 pp (cs_th, l) eps1_th x_th;;
731 float_interval_acs pp x_th;;
734 test 10 (float_interval_acs_0 pp (cs_th, l) eps1_th) x_th;;
740 (****************************************)
742 (* delta_x, delta_x4, delta_y *)
745 let delta_x' = SPEC_ALL delta_x;;
746 let x1_var_real = `x1:real` and
747 x2_var_real = `x2:real` and
748 x3_var_real = `x3:real` and
749 x4_var_real = `x4:real` and
750 x5_var_real = `x5:real` and
751 x6_var_real = `x6:real`;;
755 let float_interval_delta_x pp x1 x2 x3 x4 x5 x6 =
756 let (+++) = float_interval_add pp in
757 let (---) = float_interval_sub pp in
758 let neg = float_interval_neg in
759 let (|*|) = float_interval_mul pp in
761 let t1 = neg x1 +++ (x2 +++ (x3 --- x4 +++ (x5 +++ x6))) and
762 t2 = x1 --- x2 +++ (x3 +++ (x4 --- x5 +++ x6)) and
763 t3 = x1 +++ (x2 --- x3 +++ (x4 +++ (x5 --- x6))) in
764 let s1 = x1 |*| (x4 |*| t1) and
765 s2 = x2 |*| (x5 |*| t2) and
766 s3 = x3 |*| (x6 |*| t3) and
767 s4 = x2 |*| (x3 |*| x4) and
768 s5 = x1 |*| (x3 |*| x5) and
769 s6 = x1 |*| (x2 |*| x6) and
770 s7 = x4 |*| (x5 |*| x6) in
772 let int_th = s1 +++ (s2 +++ (s3 --- s4 --- s5 --- s6 --- s7)) in
774 let get_tm = rand o rator o concl in
775 let x1_tm = get_tm x1 and
776 x2_tm = get_tm x2 and
777 x3_tm = get_tm x3 and
778 x4_tm = get_tm x4 and
779 x5_tm = get_tm x5 and
782 let delta_th = INST[x1_tm, x1_var_real; x2_tm, x2_var_real; x3_tm, x3_var_real;
783 x4_tm, x4_var_real; x5_tm, x5_var_real; x6_tm, x6_var_real] delta_x' in
784 EQ_MP (SYM (AP_THM (AP_TERM interval_const delta_th) (rand (concl int_th)))) int_th;;
789 let delta_x4' = SPEC_ALL Sphere.delta_x4;;
791 let float_interval_delta_x4 pp x1 x2 x3 x4 x5 x6 =
792 let (+++) = float_interval_add pp in
793 let (---) = float_interval_sub pp in
794 let neg = float_interval_neg in
795 let (|*|) = float_interval_mul pp in
797 let t1 = neg x1 +++ (x2 +++ (x3 --- x4 +++ (x5 +++ x6))) in
798 let s1 = neg x2 |*| x3 and
805 let int_th = s1 --- s2 +++ (s3 +++ (s4 --- s5 +++ s6)) in
807 let get_tm = rand o rator o concl in
808 let x1_tm = get_tm x1 and
809 x2_tm = get_tm x2 and
810 x3_tm = get_tm x3 and
811 x4_tm = get_tm x4 and
812 x5_tm = get_tm x5 and
815 let delta4_th = INST[x1_tm, x1_var_real; x2_tm, x2_var_real; x3_tm, x3_var_real;
816 x4_tm, x4_var_real; x5_tm, x5_var_real; x6_tm, x6_var_real] delta_x4' in
817 EQ_MP (SYM (AP_THM (AP_TERM interval_const delta4_th) (rand (concl int_th)))) int_th;;
823 let delta_y' = SPEC_ALL Sphere.delta_y;;
825 let y1_var_real = `y1:real` and
826 y2_var_real = `y2:real` and
827 y3_var_real = `y3:real` and
828 y4_var_real = `y4:real` and
829 y5_var_real = `y5:real` and
830 y6_var_real = `y6:real`;;
833 let float_interval_delta_y pp y1 y2 y3 y4 y5 y6 =
834 let (|*|) = float_interval_mul pp in
835 let x1 = y1 |*| y1 and
841 let get_tm = rand o rator o concl in
842 let y1_tm = get_tm y1 and
843 y2_tm = get_tm y2 and
844 y3_tm = get_tm y3 and
845 y4_tm = get_tm y4 and
846 y5_tm = get_tm y5 and
849 let int_th = float_interval_delta_x pp x1 x2 x3 x4 x5 x6 in
850 let delta_th = INST[y1_tm, y1_var_real; y2_tm, y2_var_real; y3_tm, y3_var_real;
851 y4_tm, y4_var_real; y5_tm, y5_var_real; y6_tm, y6_var_real] delta_y' in
852 EQ_MP (SYM (AP_THM (AP_TERM interval_const delta_th) (rand (concl int_th)))) int_th;;
858 let x1 = mk_float_interval_small_num 1;;
859 let x2 = mk_float_interval_small_num 2;;
860 let x3 = mk_float_interval_small_num 3;;
861 let x4 = mk_float_interval_small_num 4;;
862 let x5 = mk_float_interval_small_num 5;;
863 let x6 = mk_float_interval_small_num 6;;
865 float_interval_delta_x 5 x1 x2 x3 x4 x5 x6;;
867 test 100 (float_interval_delta_x 5 x1 x2 x3 x4 x5) x6;;
873 let mk_int = mk_float_interval_small_num;;
874 let (///) = float_interval_div pp;;
876 let x1 = mk_int 1 /// mk_int 3 and
877 x2 = float_interval_sqrt pp (mk_int 2) and
878 x3 = mk_int 3 /// mk_int 11 and
879 x4 = mk_int 2 /// mk_int 61 and
880 x5 = mk_int 17 /// mk_int 13 and
881 x6 = pi_approx_array.(pp);;
883 float_interval_delta_x pp x1 x2 x3 x4 x5 x6;;
885 test 10 (float_interval_delta_x pp x1 x2 x3 x4 x5) x6;;
891 let x1 = pi_approx_array.(pp);;
892 float_interval_delta_x pp x1 x1 x1 x1 x1 x1;;
893 float_interval_delta_y pp x1 x1 x1 x1 x1 x1;;
894 float_interval_delta_x4 pp x1 x1 x1 x1 x1 x1;;
898 (**********************************)
903 let dih_x_interval = prove(`&0 < &4 * x1 * delta_x x1 x2 x3 x4 x5 x6 /\
904 interval_arith (pi / &2 + atn (--delta_x4 x1 x2 x3 x4 x5 x6 / sqrt(&4 * x1 * delta_x x1 x2 x3 x4 x5 x6))) bounds
905 ==> interval_arith (dih_x x1 x2 x3 x4 x5 x6) bounds`,
906 REWRITE_TAC[Sphere.dih_x] THEN
907 CONV_TAC (DEPTH_CONV let_CONV) THEN
908 MAP_EVERY ABBREV_TAC [`d = delta_x x1 x2 x3 x4 x5 x6`; `d4 = delta_x4 x1 x2 x3 x4 x5 x6`] THEN
910 SUBGOAL_THEN `atn2 (sqrt (&4 * x1 * d), --d4) = atn(--d4 / sqrt (&4 * x1 * d))` MP_TAC THENL
912 MP_TAC ((CONJUNCT1 o SPECL[`sqrt(&4 * x1 * d)`; `--d4`]) Trigonometry1.ATN2_BREAKDOWN) THEN
915 MATCH_MP_TAC SQRT_POS_LT THEN
923 DISCH_THEN (fun th -> ASM_REWRITE_TAC[th]));;
927 let FLOAT_INTERVAL_GT0 = prove(`interval_arith x (lo, hi) /\ (&0 < lo <=> T)
929 REWRITE_TAC[interval_arith] THEN
932 let FLOAT_INTERVAL_GT0' = (UNDISCH_ALL o PURE_REWRITE_RULE[GSYM IMP_IMP]) FLOAT_INTERVAL_GT0;;
936 let float_interval_gt0 x_th =
937 let x_tm, lo_tm, hi_tm = dest_float_interval (concl x_th) in
938 let gt_th = float_gt0 lo_tm in
939 if (fst o dest_const o rand o concl) gt_th <> "T" then
940 failwith "float_interval_gt0"
942 let th0 = INST[x_tm, x_var_real; lo_tm, lo_var_real; hi_tm, hi_var_real] FLOAT_INTERVAL_GT0' in
943 MY_PROVE_HYP x_th (MY_PROVE_HYP gt_th th0);;
948 let DIH_X_INTERVAL' = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP]) dih_x_interval;;
949 let float_interval_4 = mk_float_interval_small_num 4;;
951 let float_interval_dih_x pp x1 x2 x3 x4 x5 x6 =
952 let (|*|) = float_interval_mul pp in
953 let d_int = float_interval_delta_x pp x1 x2 x3 x4 x5 x6 in
954 let d4_int = float_interval_delta_x4 pp x1 x2 x3 x4 x5 x6 in
955 let int1 = float_interval_4 |*| (x1 |*| d_int) in
956 let x_int = float_interval_sqrt pp int1 in
957 let y_int = float_interval_neg d4_int in
958 let atn_int = float_interval_atn pp (float_interval_div pp y_int x_int) in
959 let int_th = float_interval_add pp pi2_approx_array.(pp) atn_int in
961 let bounds = (rand o concl) int_th in
962 let get_tm = rand o rator o concl in
963 let x1_tm = get_tm x1 and
964 x2_tm = get_tm x2 and
965 x3_tm = get_tm x3 and
966 x4_tm = get_tm x4 and
967 x5_tm = get_tm x5 and
970 let gt_th = float_interval_gt0 int1 in
971 let th0 = INST[bounds, bounds_var_pair;
972 x1_tm, x1_var_real; x2_tm, x2_var_real; x3_tm, x3_var_real;
973 x4_tm, x4_var_real; x5_tm, x5_var_real; x6_tm, x6_var_real] DIH_X_INTERVAL' in
974 MY_PROVE_HYP int_th (MY_PROVE_HYP gt_th th0);;
979 let x1 = pi_approx_array.(pp);;
980 float_interval_dih_x pp x1 x1 x1 x1 x1 x1;;
983 test 10 (float_interval_dih_x pp x1 x1 x1 x1 x1) x1;;
989 let dih_y' = (SYM o SPEC_ALL o CONV_RULE (DEPTH_CONV let_CONV)) Sphere.dih_y;;
992 let float_interval_dih_y pp y1 y2 y3 y4 y5 y6 =
993 let (|*|) = float_interval_mul pp in
994 let x1 = y1 |*| y1 and
1000 let get_tm = rand o rator o concl in
1001 let y1_tm = get_tm y1 and
1002 y2_tm = get_tm y2 and
1003 y3_tm = get_tm y3 and
1004 y4_tm = get_tm y4 and
1005 y5_tm = get_tm y5 and
1006 y6_tm = get_tm y6 in
1008 let int_th = float_interval_dih_x pp x1 x2 x3 x4 x5 x6 in
1009 let y_th = INST[y1_tm, y1_var_real; y2_tm, y2_var_real; y3_tm, y3_var_real;
1010 y4_tm, y4_var_real; y5_tm, y5_var_real; y6_tm, y6_var_real] dih_y' in
1011 EQ_MP (AP_THM (AP_TERM interval_const y_th) (rand (concl int_th))) int_th;;
1017 let x1 = pi_approx_array.(pp);;
1018 float_interval_dih_x pp x1 x1 x1 x1 x1 x1;;
1020 float_interval_dih_y pp x1 x1 x1 x1 x1 x1;;
1021 float_interval_dih_x pp x1 x1 x1 x1 x1 x1;;
1025 (***********************************)