formal_ineqs/arith/float_atn.hl

   1 (* =========================================================== *)
   2 (* Formal arctangent and arccosine                             *)
   3 (* Author: Alexey Solovyev                                     *)
   4 (* Date: 2012-10-27                                            *)
   5 (* =========================================================== *)
   6
   7 (* Dependencies *)
   8 needs "Multivariate/realanalysis.ml";;
   9
  10 needs "jordan/refinement.hl";;
  11 open Refinement;;
  12 needs "jordan/parse_ext_override_interface.hl";;
  13 needs "jordan/real_ext.hl";;
  14 prioritize_real();;
  15 needs "jordan/taylor_atn.hl";;
  16
  17 needs "arith/float.hl";;
  18 needs "list/more_list.hl";;
  19
  20 module type Float_atn_sig =
  21   sig
  22     val float_interval_pow_simple : int -> int -> thm -> thm
  23     val pi_approx_array : thm array
  24     val pi2_approx_array : thm array
  25
  26     val float_interval_atn : int -> thm -> thm
  27     val float_interval_acs : int -> thm -> thm
  28 end;;
  29
  30 module Float_atn : Float_atn_sig = struct
  31
  32 open Arith_misc;;
  33 open Interval_arith;;
  34 open Float_theory;;
  35 open Arith_float;;
  36 open Taylor_atn;;
  37 open More_list;;
  38
  39 (******************************)
  40 let x_var_real = `x:real` and
  41     n_var_num = `n:num` and
  42     e_var_num = `e:num` and
  43     a_var_real = `a:real` and
  44     b_var_real = `b:real` and
  45     d_var_real = `d:real` and
  46     hi_var_real = `hi:real` and
  47     lo_var_real = `lo:real`;;
  48
  49 let add_op_real = `(+):real->real->real` and
  50     sub_op_real = `(-):real->real->real` and
  51     mul_op_real = `( * ):real->real->real` and
  52     div_op_real = `(/):real->real->real` and
  53     neg_op_real = `(--):real->real` and
  54     mul_op_num = `( * ):num->num->num` and
  55     add_op_num = `(+):num->num->num`;;
  56
  57 (******************************)
  58 (* halfatn and halfatn4 *)
  59
  60 let float_interval_1 = mk_float_interval_small_num 1;;
  61
  62 let HALFATN' = (SYM o SPEC_ALL o REWRITE_RULE[REAL_POW_2]) halfatn;;
  63 let HALFATN4' = prove(`halfatn(halfatn(halfatn(halfatn x))) = halfatn4 x`,
  64                       REWRITE_TAC[halfatn4; o_THM]);;
  65
  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;;
  76
  77 let float_interval_halfatn4 pp x_th =
  78   let x_tm = (rand o rator o concl) x_th in
  79   let r_th = float_interval_halfatn pp
  80     (float_interval_halfatn pp
  81        (float_interval_halfatn pp (float_interval_halfatn pp x_th))) in
  82   let th0 = INST[x_tm, x_var_real] HALFATN4' in
  83   let ltm, rtm = dest_comb(concl r_th) in
  84     EQ_MP (AP_THM (AP_TERM (rator ltm) th0) rtm) r_th;;
  85
  86
  87 (****************************************)
  88 let rec float_interval_calc pp expr x_th =
  89   if is_var expr then
  90     x_th
  91   else
  92     let ltm, r_tm = dest_comb expr in
  93       if is_comb ltm then
  94         let op, l_tm = dest_comb ltm in
  95         let l_th = float_interval_calc pp l_tm x_th in
  96         let r_th = float_interval_calc pp r_tm x_th in
  97           if op = add_op_real then
  98             float_interval_add pp l_th r_th
  99           else if op = mul_op_real then
 100             float_interval_mul pp l_th r_th
 101           else if op = div_op_real then
 102             float_interval_div pp l_th r_th
 103           else if op = sub_op_real then
 104             float_interval_sub pp l_th r_th
 105           else
 106             failwith ("Unknown operation: " ^ (fst o dest_const) op)
 107       else
 108         if ltm = neg_op_real then
 109           let r_th = float_interval_calc pp r_tm x_th in
 110             float_interval_neg r_th
 111         else
 112           mk_float_interval_num (dest_numeral r_tm);;
 113
 114
 115 (*************************************)
 116 (* Polynomial functions *)
 117 let poly_f = new_definition `poly_f cs x = ITLIST (\c s. c + x * s) cs (&0)`;;
 118
 119 (* Even function *)
 120 let poly_f_even = new_definition `poly_f_even cs x = ITLIST (\c s. c + (x * x) * s) cs (&0)`;;
 121 (* Odd function *)
 122 let poly_f_odd = new_definition `poly_f_odd cs x = x * poly_f_even cs x`;;
 123 let poly_f_odd' = SPECL[`t:(real)list`; `x:real`] poly_f_odd;;
 124
 125 let NUMERALS_TO_NUM = Arith_nat.NUMERALS_TO_NUM;;
 126
 127 let POLY_F_EMPTY = (NUMERALS_TO_NUM o prove) (`poly_f [] x = &0`, REWRITE_TAC[poly_f; ITLIST]) and
 128     POLY_F_CONS = prove(`poly_f (CONS h t) x = h + x * poly_f t x`, REWRITE_TAC[poly_f; ITLIST]);;
 129
 130 let POLY_F_EVEN_EMPTY = (NUMERALS_TO_NUM o prove) (`poly_f_even [] x = &0`, REWRITE_TAC[poly_f_even; ITLIST]) and
 131     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]);;
 132
 133 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]);;
 134
 135 (* TABLE *)
 136 let rec reverse_table_conv tm =
 137   let ltm, i_tm = dest_comb tm in
 138     if (i_tm = `0`) then
 139       ONCE_REWRITE_CONV[REVERSE_TABLE] tm
 140     else
 141       let i_suc = num_CONV i_tm in
 142       let th1 = ONCE_REWRITE_RULE[REVERSE_TABLE] (AP_TERM ltm i_suc) in
 143       let ltm, rtm = dest_comb (rand(concl th1)) in
 144       let th2 = reverse_table_conv rtm in
 145         TRANS th1 (AP_TERM ltm th2);;
 146
 147 let atn_co_table = new_definition `atn_co_table n = TABLE
 148   (\k. (if (EVEN k) then &1 else --(&1)) / &(2 * k + 1)) (SUC n)`;;
 149
 150 (* Returns a theorem |- atn_co_table n = [...] and
 151    a list of interval approximations of the coefficients in the table *)
 152 let mk_atn_co_table pp n =
 153   let table = SPEC (mk_small_numeral n) atn_co_table in
 154   let th = CONV_RULE (DEPTH_CONV NUM_SUC_CONV THENC
 155                         REWRITE_CONV[TABLE] THENC
 156                         ONCE_DEPTH_CONV reverse_table_conv THENC
 157                         REWRITE_CONV[REVERSE; APPEND] THENC
 158                         NUM_REDUCE_CONV) table in
 159   let list = (rand o concl) th in
 160     th, map (fun tm -> float_interval_calc pp tm float_interval_1) (dest_list list);;
 161
 162 let POLY_F_EVEN_ALT = prove(`poly_f_even cs x = poly_f cs (x * x)`,
 163                             REWRITE_TAC[poly_f_even; poly_f]);;
 164
 165 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`,
 166    GEN_TAC THEN GEN_TAC THEN
 167      MATCH_MP_TAC list_INDUCT THEN
 168      REPEAT STRIP_TAC THENL
 169      [
 170        REWRITE_TAC[APPEND; poly_f; ITLIST; LENGTH] THEN
 171          REWRITE_TAC[real_pow; REAL_MUL_LID; REAL_ADD_LID];
 172        ALL_TAC
 173      ] THEN
 174
 175      REWRITE_TAC[APPEND; poly_f; ITLIST] THEN
 176      ASM_REWRITE_TAC[GSYM poly_f] THEN
 177      REWRITE_TAC[LENGTH; real_pow] THEN
 178      REAL_ARITH_TAC);;
 179
 180 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`,
 181    REWRITE_TAC[POLY_F_EVEN_ALT; POLY_F_APPEND] THEN
 182      REWRITE_TAC[GSYM REAL_POW_2; REAL_POW_POW]);;
 183
 184 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`,
 185    REPEAT GEN_TAC THEN
 186      REWRITE_TAC[poly_f_odd] THEN
 187      REWRITE_TAC[POLY_F_EVEN_APPEND] THEN
 188      REAL_ARITH_TAC);;
 189
 190 let ATN_SUM_TABLE = prove(`!x n. sum (0..n) (halfatn4_co x) = poly_f_odd (atn_co_table n) (halfatn4 x)`,
 191    GEN_TAC THEN INDUCT_TAC THENL
 192      [
 193        REWRITE_TAC[SUM_SING_NUMSEG; atn_co_table; TABLE; REVERSE_TABLE; REVERSE; APPEND] THEN
 194          REWRITE_TAC[ARITH_EVEN] THEN
 195          REWRITE_TAC[poly_f_odd; poly_f_even; ITLIST] THEN
 196          REWRITE_TAC[halfatn4_co; REAL_MUL_RZERO; REAL_ADD_RID] THEN
 197          REWRITE_TAC[MULT_CLAUSES; ARITH_ADD; REAL_POW_1; real_pow] THEN
 198          REAL_ARITH_TAC;
 199        ALL_TAC
 200      ] THEN
 201
 202      REWRITE_TAC[SUM_CLAUSES_NUMSEG; atn_co_table; TABLE; LE_0] THEN
 203      ONCE_REWRITE_TAC[REVERSE_TABLE] THEN
 204      ONCE_REWRITE_TAC[REVERSE] THEN
 205      REWRITE_TAC[GSYM TABLE; GSYM atn_co_table] THEN
 206      ASM_REWRITE_TAC[POLY_F_ODD_APPEND] THEN
 207      REWRITE_TAC[REAL_EQ_ADD_LCANCEL] THEN
 208      REWRITE_TAC[atn_co_table; LENGTH_TABLE; halfatn4_co] THEN
 209      REWRITE_TAC[poly_f_odd; poly_f_even; ITLIST; REAL_MUL_RZERO; REAL_ADD_RID] THEN
 210      REWRITE_TAC[real_div; REAL_MUL_ASSOC] THEN
 211      REWRITE_TAC[REAL_EQ_MUL_RCANCEL] THEN
 212      DISJ1_TAC THEN
 213      REWRITE_TAC[REAL_MUL_AC] THEN
 214      REWRITE_TAC[GSYM real_pow; ARITH_RULE `2 * SUC n + 1 = SUC (2 * SUC n)`] THEN
 215      GEN_REWRITE_TAC LAND_CONV [REAL_MUL_AC] THEN
 216      REWRITE_TAC[REAL_EQ_MUL_RCANCEL] THEN
 217      DISJ1_TAC THEN
 218      REWRITE_TAC[REAL_POW_NEG; real_pow; REAL_POW_ONE; REAL_MUL_RID]);;
 219
 220 let POLY_F_SING = prove(`poly_f [c] x = c`,
 221    REWRITE_TAC[poly_f; ITLIST; REAL_MUL_RZERO; REAL_ADD_RID]);;
 222
 223 (********************)
 224 let c_var_real = `c:real` and
 225     cs_var_list = `cs:(real)list` and
 226     h_var_real = `h:real` and
 227     t_var_list = `t:(real)list`;;
 228
 229 let interval_const = `interval_arith`;;
 230
 231 let rec float_interval_poly_f pp (cs, l) x_th =
 232   if length l = 0 then
 233     failwith "float_interval_poly_f: an empty coefficient list"
 234   else
 235     let ltm, x_bounds = dest_comb (concl x_th) in
 236     let x_tm = rand ltm in
 237     let first = hd l in
 238     let ltm, first_bounds = dest_comb (concl first) in
 239     let first_tm = rand ltm in
 240       if length l = 1 then
 241         let th0 = INST[first_tm, c_var_real; x_tm, x_var_real] POLY_F_SING in
 242           EQ_MP (SYM (AP_THM (AP_TERM interval_const th0) first_bounds)) first
 243       else
 244         let ltm, t_tm = dest_comb cs in
 245         let h_tm = rand ltm in
 246
 247         let th0 = INST[h_tm, h_var_real; t_tm, t_var_list; x_tm, x_var_real] POLY_F_CONS in
 248         let r_th = float_interval_poly_f pp (t_tm, tl l) x_th in
 249         let th1 = float_interval_add pp first (float_interval_mul pp x_th r_th) in
 250         let bounds = rand (concl th1) in
 251           EQ_MP (SYM (AP_THM (AP_TERM interval_const th0) bounds)) th1;;
 252
 253 let float_interval_poly_f_even pp (cs, l) x_th =
 254   let x_tm = (rand o rator o concl) x_th in
 255   let xx_th = float_interval_mul pp x_th x_th in
 256   let th0 = INST[cs, cs_var_list; x_tm, x_var_real] POLY_F_EVEN_ALT in
 257   let th1 = float_interval_poly_f pp (cs, l) xx_th in
 258   let bounds = rand(concl th1) in
 259     EQ_MP (SYM (AP_THM (AP_TERM interval_const th0) bounds)) th1;;
 260
 261 let float_interval_poly_f_odd pp (cs, l) x_th =
 262   let x_tm = (rand o rator o concl) x_th in
 263   let th0 = INST[cs, t_var_list; x_tm, x_var_real] poly_f_odd' in
 264   let even_th = float_interval_poly_f_even pp (cs, l) x_th in
 265   let th1 = float_interval_mul pp x_th even_th in
 266   let bounds = rand(concl th1) in
 267     EQ_MP (SYM (AP_THM (AP_TERM interval_const th0) bounds)) th1;;
 268
 269
 270 let poly_f_odd_const = `poly_f_odd`;;
 271 let ATN_SUM_TABLE' = SPEC_ALL ATN_SUM_TABLE;;
 272 let float_interval_16 = mk_float_interval_small_num 16;;
 273
 274 (* Computes an interval for &16 * sum(0..n) (halfatn4_co x) *)
 275 let float_interval_atn_sum pp (cs_th, l) x_th =
 276   let n_tm = (rand o lhand o concl) cs_th in
 277   let cs_tm = rand(concl cs_th) in
 278   let halfatn4 = float_interval_halfatn4 pp x_th in
 279
 280   let poly_th = float_interval_poly_f_odd pp (cs_tm, l) halfatn4 in
 281   let bounds = rand (concl poly_th) in
 282   let halfatn4_tm = (rand o rator o concl) halfatn4 in
 283   let x_tm = rand halfatn4_tm in
 284
 285   let th1 = AP_THM (AP_TERM interval_const (AP_THM (AP_TERM poly_f_odd_const cs_th) halfatn4_tm)) bounds in
 286   let poly_atn_th = EQ_MP (SYM th1) poly_th in
 287   let bounds = rand (concl poly_atn_th) in
 288
 289   let th2 = INST[n_tm, n_var_num; x_tm, x_var_real] ATN_SUM_TABLE' in
 290   let th3 = EQ_MP (SYM (AP_THM (AP_TERM interval_const th2) bounds)) poly_atn_th in
 291     float_interval_mul pp float_interval_16 th3;;
 292
 293
 294 (******************************)
 295 let bounds_var_pair = `bounds:real#real`;;
 296
 297 let FLOAT_INTERVAL_INV = prove(`interval_arith (&1 / x) bounds <=>
 298                                  interval_arith (inv x) bounds`,
 299    REWRITE_TAC[real_div; REAL_MUL_LID]);;
 300
 301 let float_interval_inv pp x_th =
 302   let x_tm = (rand o rator o concl) x_th in
 303   let r_th = float_interval_div pp float_interval_1 x_th in
 304   let th0 = INST[x_tm, x_var_real; rand(concl r_th), bounds_var_pair] FLOAT_INTERVAL_INV in
 305     EQ_MP th0 r_th;;
 306
 307 let REAL_POW_SUC = (SPEC_ALL o CONJUNCT2) real_pow;;
 308
 309 let INTERVAL_REAL_POW_0 = prove(mk_comb(`interval_arith (x pow 0)`, (rand o concl) float_interval_1),
 310                                 REWRITE_TAC[real_pow; float_interval_1]);;
 311
 312 let INTERVAL_REAL_POW_1 = prove(`interval_arith x bounds <=> interval_arith (x pow 1) bounds`,
 313                                 REWRITE_TAC[REAL_POW_1]);;
 314
 315
 316 let rec float_interval_pow_simple pp n x_th =
 317   let x_tm = (rand o rator o concl) x_th in
 318     if n = 0 then
 319       INST[x_tm, x_var_real] INTERVAL_REAL_POW_0
 320     else
 321       if n = 1 then
 322         let bounds = rand(concl x_th) in
 323         let th0 = INST[x_tm, x_var_real; bounds, bounds_var_pair] INTERVAL_REAL_POW_1 in
 324           EQ_MP th0 x_th
 325       else
 326         let n_tm' = mk_small_numeral n in
 327         let n_suc = num_CONV n_tm' in
 328         let n_tm = rand(rand(concl n_suc)) in
 329         let th0 = INST[x_tm, x_var_real; n_tm, n_var_num] REAL_POW_SUC in
 330         let r_th = float_interval_pow_simple pp (n - 1) x_th in
 331         let th1 = float_interval_mul pp x_th r_th in
 332         let bounds = rand (concl th1) in
 333         let th2 = TRANS (AP_TERM (rator(lhand(concl th0))) n_suc) th0 in
 334           EQ_MP (SYM (AP_THM (AP_TERM interval_const th2) bounds)) th1;;
 335
 336 let float_interval_2 = mk_float_interval_small_num 2 and
 337     six_const = `6` and
 338     five_const = `5`;;
 339
 340 (* Computes an interval for inv(&2 pow (6 * n + 5)) *)
 341 let compute_eps1 pp n =
 342   let n_tm = mk_small_numeral n in
 343   let n6 = NUM_MULT_CONV (mk_binop mul_op_num six_const n_tm) in
 344   let n65_1 = AP_THM (AP_TERM add_op_num n6) five_const in
 345   let n65_2 = NUM_ADD_CONV (rand (concl n65_1)) in
 346   let n65 = TRANS n65_1 n65_2 in
 347   let pow_th = float_interval_pow_simple pp (6 * n + 5) float_interval_2 in
 348   let ltm, bounds = dest_comb(concl pow_th) in
 349   let pow_tm = (rator o rand) ltm in
 350   let th0 = EQ_MP (SYM (AP_THM (AP_TERM interval_const (AP_TERM pow_tm n65)) bounds)) pow_th in
 351     float_interval_inv pp th0;;
 352
 353
 354 (**********************************)
 355 let FLOAT_ATN_LO_HI = prove(`interval_arith (&16 * sum(0..n) (halfatn4_co x)) (a, b) /\
 356                               interval_arith (inv(&2 pow (6*n + 5))) (c,d) /\
 357                               b + d <= hi /\ lo <= a - d
 358                               ==> interval_arith (atn x) (lo, hi)`,
 359    REWRITE_TAC[interval_arith] THEN
 360      STRIP_TAC THEN
 361      MP_TAC (SPEC_ALL real_taylor_atn_halfatn4) THEN
 362      MP_TAC (REAL_ARITH `&0 <= abs(&16)`) THEN
 363      REWRITE_TAC[IMP_IMP] THEN
 364      DISCH_THEN (MP_TAC o MATCH_MP REAL_LE_LMUL) THEN
 365      REWRITE_TAC[GSYM REAL_ABS_MUL; REAL_ARITH `a * (b - c) = a * b - a * c:real`] THEN
 366      ONCE_REWRITE_TAC[GSYM atn_halfatn4] THEN
 367      REWRITE_TAC[REAL_ARITH `abs (x - v) <= e <=> v - e <= x /\ x <= v + e`] THEN
 368      REWRITE_TAC[REAL_ABS_NUM] THEN
 369      SUBGOAL_THEN `&16 * inv(&8 pow (2 * n + 3)) = inv(&2 pow (6 * n + 5))` (fun th -> REWRITE_TAC[th]) THENL
 370      [
 371        REWRITE_TAC[GSYM real_div] THEN
 372          SUBGOAL_THEN `&16 = &2 pow 4 /\ &8 = &2 pow 3 /\ ~(&2 = &0)` ASSUME_TAC THENL
 373          [
 374            REAL_ARITH_TAC;
 375            ALL_TAC
 376          ] THEN
 377          ASM_REWRITE_TAC[REAL_POW_POW] THEN
 378          ASM_SIMP_TAC[REAL_DIV_POW2] THEN
 379          REWRITE_TAC[ARITH_RULE `~(3 * (2 * n + 3) <= 4)`] THEN
 380          REPEAT AP_TERM_TAC THEN
 381          ARITH_TAC;
 382        ALL_TAC
 383      ] THEN
 384      REWRITE_TAC[GSYM halfatn4_co] THEN
 385      SUBGOAL_THEN `sum (0..n) (\j. halfatn4_co x j) = sum (0..n) (halfatn4_co x)` (fun th -> REWRITE_TAC[th]) THENL
 386      [
 387        AP_TERM_TAC THEN
 388          REWRITE_TAC[FUN_EQ_THM];
 389        ALL_TAC
 390      ] THEN
 391      REPEAT (POP_ASSUM MP_TAC) THEN
 392      REAL_ARITH_TAC);;
 393
 394 let FLOAT_ATN_LO_HI' = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP]) FLOAT_ATN_LO_HI;;
 395
 396 let float_interval_atn_0 pp (cs_th, l) eps1_th x_th =
 397   let sum_th = float_interval_atn_sum pp (cs_th, l) x_th in
 398   let n_tm = (rand o lhand o concl) cs_th in
 399   let x_tm = (rand o rator o concl) x_th in
 400
 401   let sum_bounds = rand (concl sum_th) in
 402   let a_tm, b_tm = dest_pair sum_bounds in
 403   let c_tm, d_tm = (dest_pair o rand o concl) eps1_th in
 404
 405   let hi_th = float_add_hi pp b_tm d_tm in
 406   let lo_th = float_sub_lo pp a_tm d_tm in
 407   let hi_tm = rand(concl hi_th) in
 408   let lo_tm = lhand(concl lo_th) in
 409
 410   let th0 = INST[n_tm, n_var_num; x_tm, x_var_real;
 411                  a_tm, a_var_real; b_tm, b_var_real;
 412                  c_tm, c_var_real; d_tm, d_var_real;
 413                  hi_tm, hi_var_real; lo_tm, lo_var_real] FLOAT_ATN_LO_HI' in
 414     MY_PROVE_HYP lo_th (MY_PROVE_HYP hi_th (MY_PROVE_HYP sum_th (MY_PROVE_HYP eps1_th th0)));;
 415
 416 (* Fill in lookup tables *)
 417
 418 (* Computes n such that 2^(-(6n + 5)) <= base^(-(p + 1)) *)
 419 let n_of_p pp =
 420   let x = (float_of_int (pp + 1) *. log (float_of_int Arith_hash.arith_base) /. log (2.0) -. 5.0) /. 6.0 in
 421   let n = (int_of_float o ceil) x in
 422     if n < 1 then 1 else n;;
 423
 424 let atn_co_array = Array.init 21 (fun i -> mk_atn_co_table (i + 1) (n_of_p i));;
 425 let eps1_array = Array.init 21 (fun i -> compute_eps1 (i + 1) (n_of_p i));;
 426
 427 let float_interval_atn pp x_th =
 428   float_interval_atn_0 pp atn_co_array.(pp) eps1_array.(pp) x_th;;
 429
 430 (*****************************************)
 431 (* pi approximation *)
 432
 433 let pp = 20;;
 434 let x_th = float_interval_1;;
 435 let th1 = float_interval_atn pp x_th;;
 436 let th2 = float_interval_mul pp (mk_float_interval_small_num 4) th1;;
 437 let float_interval_pi = REWRITE_RULE[ATN_1; REAL_ARITH `&4 * pi / &4 = pi`] th2;;
 438 let float_interval_pi2 = float_interval_div pp float_interval_pi float_interval_2;;
 439
 440
 441 let pi_approx_array = Array.init 19 (fun i -> float_interval_round i float_interval_pi);;
 442 let pi2_approx_array = Array.init 19 (fun i -> float_interval_round i float_interval_pi2);;
 443
 444 (********************************************)
 445 (* acs *)
 446
 447 let TAN_HALF = prove(`!x. ~(cos x = -- &1) ==> tan (x / &2) = sin x / (&1 + cos x)`,
 448    GEN_TAC THEN
 449      ABBREV_TAC `t = x / &2` THEN
 450      SUBGOAL_THEN `x = &2 * t` ASSUME_TAC THENL
 451      [
 452        EXPAND_TAC "t" THEN REAL_ARITH_TAC;
 453        ALL_TAC
 454      ] THEN
 455
 456      ASM_REWRITE_TAC[SIN_DOUBLE; COS_DOUBLE_COS; REAL_ARITH `&1 + a - &1 = a`] THEN
 457      REWRITE_TAC[REAL_ARITH `a - &1 = -- &1 <=> a = &0`] THEN
 458      REWRITE_TAC[real_div; REAL_INV_MUL; REAL_POW_2] THEN
 459      REWRITE_TAC[REAL_ENTIRE; REAL_ARITH `&2 = &0 <=> F`] THEN
 460      DISCH_TAC THEN
 461      REWRITE_TAC[REAL_ARITH `(&2 * s * c) * i2 * ic * ic = (&2 * i2) * (c * ic) * s * ic`] THEN
 462      ASM_SIMP_TAC[REAL_MUL_RINV; REAL_ARITH `~(&2 = &0)`] THEN
 463      REWRITE_TAC[REAL_MUL_LID; tan; real_div]);;
 464
 465 let X_EQ_COS_T = prove(`!x. abs x <= &1 ==> ?t. &0 <= t /\ t <= pi /\ x = cos t`,
 466    REWRITE_TAC[REAL_ARITH `abs x <= &1 <=> -- &1 <= x /\ x <= &1`] THEN
 467      REPEAT STRIP_TAC THEN
 468      EXISTS_TAC `acs x` THEN
 469      ASM_SIMP_TAC[ACS_BOUNDS; COS_ACS]);;
 470
 471 let ACS_ATN_ALT = prove(`!x. -- &1 < x /\ x <= &1 ==>
 472                           acs x = &2 * atn (sqrt (&1 - x pow 2) / (&1 + x))`,
 473    REPEAT STRIP_TAC THEN
 474      MP_TAC (SPEC_ALL X_EQ_COS_T) THEN
 475      ANTS_TAC THENL
 476      [
 477        REPEAT (POP_ASSUM MP_TAC) THEN
 478          REAL_ARITH_TAC;
 479        ALL_TAC
 480      ] THEN
 481      REPEAT STRIP_TAC THEN
 482      ASM_REWRITE_TAC[] THEN
 483      ASM_SIMP_TAC[ACS_COS] THEN
 484      MP_TAC (SPEC `t:real` SIN_COS_SQRT) THEN
 485      ANTS_TAC THENL
 486      [
 487        ASM_SIMP_TAC[SIN_POS_PI_LE];
 488        ALL_TAC
 489      ] THEN
 490
 491      DISCH_THEN (fun th -> REWRITE_TAC[SYM th]) THEN
 492      MP_TAC (SPEC `t:real` TAN_HALF) THEN
 493      ANTS_TAC THENL
 494      [
 495        POP_ASSUM (fun th -> REWRITE_TAC[SYM th]) THEN
 496          UNDISCH_TAC `-- &1 < x` THEN
 497          REAL_ARITH_TAC;
 498        ALL_TAC
 499      ] THEN
 500
 501      DISCH_THEN (fun th -> REWRITE_TAC[SYM th]) THEN
 502      REWRITE_TAC[REAL_ARITH `t = &2 * a <=> a = t / &2`] THEN
 503      MATCH_MP_TAC TAN_ATN THEN
 504      REWRITE_TAC[REAL_ARITH `a / &2 < b / &2 <=> a < b`] THEN
 505      REWRITE_TAC[REAL_ARITH `--(a / &2) < b / &2 <=> --a < b`] THEN
 506      CONJ_TAC THENL
 507      [
 508        MATCH_MP_TAC REAL_LTE_TRANS THEN
 509          EXISTS_TAC `&0` THEN
 510          ASM_REWRITE_TAC[REAL_NEG_LT0; PI_POS];
 511        SUBGOAL_THEN `t = acs x` MP_TAC THENL
 512          [
 513            ASM_SIMP_TAC[ACS_COS];
 514            ALL_TAC
 515          ] THEN
 516          DISCH_THEN (fun th -> REWRITE_TAC[th]) THEN
 517          REWRITE_TAC[SYM ACS_NEG_1] THEN
 518          MATCH_MP_TAC ACS_MONO_LT THEN
 519          ASM_REWRITE_TAC[REAL_LE_REFL]
 520      ]);;
 521
 522 let FLOAT_F_LT = prove(`!n e. &0 < float_num F n e <=> ~(n = 0)`,
 523                        REWRITE_TAC[REAL_ARITH `&0 < a <=> &0 <= a /\ ~(a = &0)`] THEN
 524                          REWRITE_TAC[FLOAT_F_POS; FLOAT_EQ_0]);;
 525
 526 let FLOAT_INTERVAL_ACS = prove(`interval_arith (pi / &2 - atn(x / sqrt(&1 - x * x))) bounds /\
 527                                  interval_arith (&1 - x * x) (float_num F n e, hi) /\
 528                                  ~(n = 0)
 529                                  ==> interval_arith (acs x) bounds`,
 530    REWRITE_TAC[GSYM REAL_POW_2] THEN
 531      STRIP_TAC THEN
 532      MP_TAC (SPEC_ALL ACS_ATN) THEN
 533      ANTS_TAC THENL
 534      [
 535        POP_ASSUM MP_TAC THEN POP_ASSUM MP_TAC THEN
 536          REWRITE_TAC[interval_arith] THEN
 537          REWRITE_TAC[REAL_ARITH `-- &1 < x /\ x < &1 <=> abs x < abs (&1)`] THEN
 538          REWRITE_TAC[REAL_LT_SQUARE_ABS] THEN
 539          REWRITE_TAC[REAL_ARITH `&1 pow 2 = &1`] THEN
 540          REPEAT STRIP_TAC THEN
 541          REWRITE_TAC[REAL_ARITH `a < &1 <=> &0 < &1 - a`] THEN
 542
 543          MP_TAC (SPEC_ALL FLOAT_F_LT) THEN
 544          ASM_REWRITE_TAC[] THEN
 545          UNDISCH_TAC `float_num F n e <= &1 - x pow 2` THEN
 546          REAL_ARITH_TAC;
 547        ALL_TAC
 548      ] THEN
 549
 550      DISCH_THEN (fun th -> ASM_REWRITE_TAC[th]));;
 551
 552 let ZERO_EQ_ZERO_CONST = prove(`0 = _0`, REWRITE_TAC[NUMERAL]);;
 553
 554 let FLOAT_INTERVAL_ACS' = (UNDISCH_ALL o
 555                              PURE_ONCE_REWRITE_RULE[TAUT `~P <=> (P <=> F)`] o
 556                              REWRITE_RULE[ZERO_EQ_ZERO_CONST; GSYM IMP_IMP]) FLOAT_INTERVAL_ACS;;
 557
 558 let float_interval_acs_0 pp (cs_th, l) eps1_th x_th =
 559   let int1 = float_interval_sub pp float_interval_1 (float_interval_mul pp x_th x_th) in
 560   let int2 = float_interval_div pp x_th (float_interval_sqrt pp int1) in
 561   let atn_int = float_interval_atn_0 pp (cs_th, l) eps1_th int2 in
 562   let acs_int = float_interval_sub pp pi2_approx_array.(pp + 1) atn_int in
 563
 564   let x_tm = (rand o rator o concl) x_th in
 565   let bounds = (rand o concl) acs_int in
 566   let int1_bounds = (rand o concl) int1 in
 567   let lo_tm, hi_tm = dest_pair int1_bounds in
 568   let s, n_tm, e_tm = dest_float lo_tm in
 569     if s <> "F" then
 570       failwith "float_interval_acs_0: &1 - x pow 2 < &1 is not satisfied"
 571     else
 572       let n_th = Arith_nat.raw_eq0_hash_conv n_tm in
 573       let th0 = INST[x_tm, x_var_real; bounds, bounds_var_pair;
 574                      n_tm, n_var_num; e_tm, e_var_num;
 575                      hi_tm, hi_var_real] FLOAT_INTERVAL_ACS' in
 576         MY_PROVE_HYP acs_int (MY_PROVE_HYP int1 (MY_PROVE_HYP n_th th0));;
 577
 578 let float_interval_acs pp x_th =
 579   float_interval_acs_0 pp atn_co_array.(pp) eps1_array.(pp) x_th;;
 580
 581 (****************************************)
 582 end;;