formal_lp/old/arith/float_atn.hl

   1 (* Dependencies *)
   2 needs "jordan/refinement.hl";;
   3 open Refinement;;
   4 needs "jordan/parse_ext_override_interface.hl";;
   5 needs "jordan/real_ext.hl";;
   6 prioritize_real();;
   7 needs "jordan/taylor_atn.hl";;
   8
   9 needs "../formal_lp/arith/float.hl";;
  10
  11
  12 module type Float_atn_sig =
  13   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
  17
  18     val float_interval_atn : int -> thm -> thm
  19     val float_interval_acs : int -> thm -> thm
  20 end;;
  21
  22
  23 module Float_atn : Float_atn_sig = struct
  24
  25 open Arith_misc;;
  26 open Interval_arith;;
  27 open Float_theory;;
  28 open Arith_float;;
  29 open Taylor_atn;;
  30
  31
  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`;;
  41
  42
  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`;;
  50
  51
  52
  53 (******************************)
  54 (* halfatn and halfatn4 *)
  55
  56 let float_interval_1 = mk_float_interval_small_num 1;;
  57
  58
  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]);;
  62
  63
  64
  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
  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;;
  86
  87
  88
  89 (*
  90 let x_th = mk_float_interval_small_num 3;;
  91 let x_th2 = float_interval_div 11 float_interval_1 x_th;;
  92
  93 float_interval_halfatn4 11 x_th2;;
  94 (* 10: 1.180 *)
  95 test 10 (float_interval_halfatn4 11) x_th2;;
  96 *)
  97
  98
  99
 100 (****************************************)
 101
 102
 103 let rec float_interval_calc pp expr x_th =
 104   if is_var expr then
 105     x_th
 106   else
 107     let ltm, r_tm = dest_comb expr in
 108       if is_comb ltm then
 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
 120           else
 121             failwith ("Unknown operation: " ^ (fst o dest_const) op)
 122       else
 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
 126         else
 127           mk_float_interval_num (dest_numeral r_tm);;
 128
 129
 130 (*************************************)
 131 (* Polynomial functions *)
 132
 133 let poly_f = new_definition `poly_f cs x = ITLIST (\c s. c + x * s) cs (&0)`;;
 134
 135 (* Even function *)
 136 let poly_f_even = new_definition `poly_f_even cs x = ITLIST (\c s. c + (x * x) * s) cs (&0)`;;
 137 (* Odd function *)
 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;;
 140
 141 let NUMERALS_TO_NUM = Arith_nat.NUMERALS_TO_NUM;;
 142
 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]);;
 145
 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]);;
 148
 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]);;
 150
 151
 152 (* TABLE *)
 153
 154 let REVERSE_TABLE  = define `(REVERSE_TABLE (f:num->A) 0 = []) /\
 155    (REVERSE_TABLE f (SUC i) = CONS (f i)  ( REVERSE_TABLE f i))`;;
 156
 157
 158 let TABLE = new_definition `!(f:num->A) k. TABLE f k = REVERSE (REVERSE_TABLE f k)`;;
 159
 160
 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]);;
 163
 164
 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
 169      ARITH_TAC);;
 170
 171
 172 let LENGTH_TABLE = prove(`!(f:num->A) n. LENGTH (TABLE f n) = n`,
 173    REWRITE_TAC[TABLE; LENGTH_REVERSE; LENGTH_REVERSE_TABLE]);;
 174
 175
 176
 177 let rec reverse_table_conv tm =
 178   let ltm, i_tm = dest_comb tm in
 179     if (i_tm = `0`) then
 180       ONCE_REWRITE_CONV[REVERSE_TABLE] tm
 181     else
 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);;
 187
 188
 189
 190
 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)`;;
 193
 194
 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);;
 206
 207
 208
 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]);;
 211
 212
 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
 217      [
 218        REWRITE_TAC[APPEND; poly_f; ITLIST; LENGTH] THEN
 219          REWRITE_TAC[real_pow; REAL_MUL_LID; REAL_ADD_LID];
 220        ALL_TAC
 221      ] THEN
 222
 223      REWRITE_TAC[APPEND; poly_f; ITLIST] THEN
 224      ASM_REWRITE_TAC[GSYM poly_f] THEN
 225      REWRITE_TAC[LENGTH; real_pow] THEN
 226      REAL_ARITH_TAC);;
 227
 228
 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]);;
 232
 233
 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`,
 235    REPEAT GEN_TAC THEN
 236      REWRITE_TAC[poly_f_odd] THEN
 237      REWRITE_TAC[POLY_F_EVEN_APPEND] THEN
 238      REAL_ARITH_TAC);;
 239
 240
 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
 243      [
 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
 249          REAL_ARITH_TAC;
 250        ALL_TAC
 251      ] THEN
 252
 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
 263      DISJ1_TAC 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
 268      DISJ1_TAC THEN
 269      REWRITE_TAC[REAL_POW_NEG; real_pow; REAL_POW_ONE; REAL_MUL_RID]);;
 270
 271
 272
 273 let POLY_F_SING = prove(`poly_f [c] x = c`,
 274    REWRITE_TAC[poly_f; ITLIST; REAL_MUL_RZERO; REAL_ADD_RID]);;
 275
 276
 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`;;
 281
 282 let interval_const = `interval_arith`;;
 283
 284
 285 let rec float_interval_poly_f pp (cs, l) x_th =
 286   if length l = 0 then
 287     failwith "float_interval_poly_f: an empty coefficient list"
 288   else
 289     let ltm, x_bounds = dest_comb (concl x_th) in
 290     let x_tm = rand ltm in
 291     let first = hd l in
 292     let ltm, first_bounds = dest_comb (concl first) in
 293     let first_tm = rand ltm in
 294       if length l = 1 then
 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
 297       else
 298         let ltm, t_tm = dest_comb cs in
 299         let h_tm = rand ltm in
 300
 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;;
 306
 307
 308
 309
 310
 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;;
 318
 319
 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;;
 327
 328
 329
 330
 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;;
 334
 335
 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
 341
 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
 346
 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
 350
 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;;
 354
 355
 356
 357
 358 (*
 359 let pp = 10;;
 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;;
 363
 364 let th = float_interval_poly_f pp (cs, l) x_th;;
 365 float_interval_mul pp (mk_float_interval_small_num 4) th;;
 366
 367
 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;;
 370 *)
 371
 372
 373 (******************************)
 374
 375 let bounds_var_pair = `bounds:real#real`;;
 376
 377
 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]);;
 381
 382
 383
 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
 388     EQ_MP th0 r_th;;
 389
 390
 391 let REAL_POW_SUC = (SPEC_ALL o CONJUNCT2) real_pow;;
 392
 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]);;
 395
 396 let INTERVAL_REAL_POW_1 = prove(`interval_arith x bounds <=> interval_arith (x pow 1) bounds`,
 397                                 REWRITE_TAC[REAL_POW_1]);;
 398
 399
 400
 401
 402 let rec float_interval_pow_simple pp n x_th =
 403   let x_tm = (rand o rator o concl) x_th in
 404     if n = 0 then
 405       INST[x_tm, x_var_real] INTERVAL_REAL_POW_0
 406     else
 407       if n = 1 then
 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
 410           EQ_MP th0 x_th
 411       else
 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;;
 421
 422
 423
 424 let float_interval_2 = mk_float_interval_small_num 2 and
 425     six_const = `6` and
 426     five_const = `5`;;
 427
 428
 429
 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;;
 442
 443
 444
 445
 446 (**********************************)
 447
 448
 449
 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
 455      STRIP_TAC 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
 465      [
 466        REWRITE_TAC[GSYM real_div] THEN
 467          SUBGOAL_THEN `&16 = &2 pow 4 /\ &8 = &2 pow 3 /\ ~(&2 = &0)` ASSUME_TAC THENL
 468          [
 469            REAL_ARITH_TAC;
 470            ALL_TAC
 471          ] THEN
 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
 476          ARITH_TAC;
 477        ALL_TAC
 478      ] 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
 481      [
 482        AP_TERM_TAC THEN
 483          REWRITE_TAC[FUN_EQ_THM];
 484        ALL_TAC
 485      ] THEN
 486      REPEAT (POP_ASSUM MP_TAC) THEN
 487      REAL_ARITH_TAC);;
 488
 489
 490 let FLOAT_ATN_LO_HI' = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP]) FLOAT_ATN_LO_HI;;
 491
 492
 493
 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
 498
 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
 502
 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
 507
 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)));;
 513
 514
 515 (*
 516 let pp = 10;;
 517 let n = 5;;
 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;;
 521
 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;;
 525 *)
 526
 527
 528
 529 (* Fill in lookup tables *)
 530
 531 (* Computes n such that 2^(-(6n + 5)) <= base^(-(p + 1)) *)
 532 let n_of_p pp =
 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;;
 536
 537
 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));;
 540
 541
 542 let float_interval_atn pp x_th =
 543   float_interval_atn_0 pp atn_co_array.(pp) eps1_array.(pp) x_th;;
 544
 545
 546
 547
 548 (*****************************************)
 549
 550 (* pi approximation *)
 551
 552 let pp = 20;;
 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;;
 558
 559
 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);;
 562
 563
 564
 565
 566 (********************************************)
 567
 568 (* acs *)
 569
 570 let TAN_HALF = prove(`!x. ~(cos x = -- &1) ==> tan (x / &2) = sin x / (&1 + cos x)`,
 571    GEN_TAC THEN
 572      ABBREV_TAC `t = x / &2` THEN
 573      SUBGOAL_THEN `x = &2 * t` ASSUME_TAC THENL
 574      [
 575        EXPAND_TAC "t" THEN REAL_ARITH_TAC;
 576        ALL_TAC
 577      ] THEN
 578
 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
 583      DISCH_TAC 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]);;
 587
 588
 589
 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]);;
 595
 596
 597
 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
 602      ANTS_TAC THENL
 603      [
 604        REPEAT (POP_ASSUM MP_TAC) THEN
 605          REAL_ARITH_TAC;
 606        ALL_TAC
 607      ] 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
 612      ANTS_TAC THENL
 613      [
 614        ASM_SIMP_TAC[SIN_POS_PI_LE];
 615        ALL_TAC
 616      ] THEN
 617
 618      DISCH_THEN (fun th -> REWRITE_TAC[SYM th]) THEN
 619      MP_TAC (SPEC `t:real` TAN_HALF) THEN
 620      ANTS_TAC THENL
 621      [
 622        POP_ASSUM (fun th -> REWRITE_TAC[SYM th]) THEN
 623          UNDISCH_TAC `-- &1 < x` THEN
 624          REAL_ARITH_TAC;
 625        ALL_TAC
 626      ] THEN
 627
 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
 633      CONJ_TAC THENL
 634      [
 635        MATCH_MP_TAC REAL_LTE_TRANS THEN
 636          EXISTS_TAC `&0` THEN
 637          ASM_REWRITE_TAC[REAL_NEG_LT0; PI_POS];
 638        SUBGOAL_THEN `t = acs x` MP_TAC THENL
 639          [
 640            ASM_SIMP_TAC[ACS_COS];
 641            ALL_TAC
 642          ] THEN
 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]
 647      ]);;
 648
 649
 650
 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]);;
 654
 655
 656
 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) /\
 659                                  ~(n = 0)
 660                                  ==> interval_arith (acs x) bounds`,
 661    REWRITE_TAC[GSYM REAL_POW_2] THEN
 662      STRIP_TAC THEN
 663      MP_TAC (SPEC_ALL ACS_ATN) THEN
 664      ANTS_TAC THENL
 665      [
 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
 673
 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
 677          REAL_ARITH_TAC;
 678        ALL_TAC
 679      ] THEN
 680
 681      DISCH_THEN (fun th -> ASM_REWRITE_TAC[th]));;
 682
 683
 684
 685
 686 let ZERO_EQ_ZERO_CONST = prove(`0 = _0`, REWRITE_TAC[NUMERAL]);;
 687
 688
 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;;
 692
 693
 694
 695
 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
 701
 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
 707     if s <> "F" then
 708       failwith "float_interval_acs_0: &1 - x pow 2 < &1 is not satisfied"
 709     else
 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));;
 715
 716
 717 let float_interval_acs pp x_th =
 718   float_interval_acs_0 pp atn_co_array.(pp) eps1_array.(pp) x_th;;
 719
 720
 721 (*
 722 (* acs(&1 / &3) <= #1.230959418 *)
 723 let pp = 11;;
 724 let n = 5;;
 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);;
 728
 729
 730 float_interval_acs_0 pp (cs_th, l) eps1_th x_th;;
 731 float_interval_acs pp x_th;;
 732
 733 (* 10: 1.908 *)
 734 test 10 (float_interval_acs_0 pp (cs_th, l) eps1_th) x_th;;
 735
 736 *)
 737
 738
 739
 740 (****************************************)
 741
 742 (* delta_x, delta_x4, delta_y *)
 743
 744
 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`;;
 752
 753
 754
 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
 760
 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
 771
 772   let int_th = s1 +++ (s2 +++ (s3 --- s4 --- s5 --- s6 --- s7)) in
 773
 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
 780       x6_tm = get_tm x6 in
 781
 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;;
 785
 786
 787 (* delta_x4 *)
 788
 789 let delta_x4' = SPEC_ALL Sphere.delta_x4;;
 790
 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
 796
 797   let t1 = neg x1 +++ (x2 +++ (x3 --- x4 +++ (x5 +++ x6))) in
 798   let s1 = neg x2 |*| x3 and
 799       s2 = x1 |*| x4 and
 800       s3 = x2 |*| x5 and
 801       s4 = x3 |*| x6 and
 802       s5 = x5 |*| x6 and
 803       s6 = x1 |*| t1 in
 804
 805   let int_th = s1 --- s2 +++ (s3 +++ (s4 --- s5 +++ s6)) in
 806
 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
 813       x6_tm = get_tm x6 in
 814
 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;;
 818
 819
 820
 821 (* delta_y *)
 822
 823 let delta_y' = SPEC_ALL Sphere.delta_y;;
 824
 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`;;
 831
 832
 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
 836       x2 = y2 |*| y2 and
 837       x3 = y3 |*| y3 and
 838       x4 = y4 |*| y4 and
 839       x5 = y5 |*| y5 and
 840       x6 = y6 |*| y6 in
 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
 847       y6_tm = get_tm y6 in
 848
 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;;
 853
 854
 855
 856
 857 (*
 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;;
 864
 865 float_interval_delta_x 5 x1 x2 x3 x4 x5 x6;;
 866 (* 100: 1.248 *)
 867 test 100 (float_interval_delta_x 5 x1 x2 x3 x4 x5) x6;;
 868 *)
 869
 870
 871 (*
 872 let pp = 7;;
 873 let mk_int = mk_float_interval_small_num;;
 874 let (///) = float_interval_div pp;;
 875
 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);;
 882
 883 float_interval_delta_x pp x1 x2 x3 x4 x5 x6;;
 884 (* 100: 0.632 *)
 885 test 10 (float_interval_delta_x pp x1 x2 x3 x4 x5) x6;;
 886 *)
 887
 888
 889 (*
 890 let pp = 10;;
 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;;
 895 *)
 896
 897
 898 (**********************************)
 899
 900 (* dih_x *)
 901
 902
 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
 909      STRIP_TAC THEN
 910      SUBGOAL_THEN `atn2 (sqrt (&4 * x1 * d), --d4) = atn(--d4 / sqrt (&4 * x1 * d))` MP_TAC THENL
 911      [
 912        MP_TAC ((CONJUNCT1 o SPECL[`sqrt(&4 * x1 * d)`; `--d4`]) Trigonometry1.ATN2_BREAKDOWN) THEN
 913          ANTS_TAC THENL
 914          [
 915            MATCH_MP_TAC SQRT_POS_LT THEN
 916              ASM_REWRITE_TAC[];
 917            ALL_TAC
 918          ] THEN
 919          REWRITE_TAC[];
 920        ALL_TAC
 921      ] THEN
 922
 923      DISCH_THEN (fun th -> ASM_REWRITE_TAC[th]));;
 924
 925
 926
 927 let FLOAT_INTERVAL_GT0 = prove(`interval_arith x (lo, hi) /\ (&0 < lo <=> T)
 928                                  ==> &0 < x`,
 929                                REWRITE_TAC[interval_arith] THEN
 930                                  REAL_ARITH_TAC);;
 931
 932 let FLOAT_INTERVAL_GT0' = (UNDISCH_ALL o PURE_REWRITE_RULE[GSYM IMP_IMP]) FLOAT_INTERVAL_GT0;;
 933
 934
 935
 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"
 941   else
 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);;
 944
 945
 946
 947
 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;;
 950
 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
 960
 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
 968       x6_tm = get_tm x6 in
 969
 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);;
 975
 976
 977 (*
 978 let pp = 7;;
 979 let x1 = pi_approx_array.(pp);;
 980 float_interval_dih_x pp x1 x1 x1 x1 x1 x1;;
 981
 982 (* 10: 1.628 *)
 983 test 10 (float_interval_dih_x pp x1 x1 x1 x1 x1) x1;;
 984 *)
 985
 986
 987 (* dih_y *)
 988
 989 let dih_y' = (SYM o SPEC_ALL o CONV_RULE (DEPTH_CONV let_CONV)) Sphere.dih_y;;
 990
 991
 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
 995       x2 = y2 |*| y2 and
 996       x3 = y3 |*| y3 and
 997       x4 = y4 |*| y4 and
 998       x5 = y5 |*| y5 and
 999       x6 = y6 |*| y6 in
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
1007
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;;
1012
1013
1014 (*
1015
1016 let pp = 15;;
1017 let x1 = pi_approx_array.(pp);;
1018 float_interval_dih_x pp x1 x1 x1 x1 x1 x1;;
1019
1020 float_interval_dih_y pp x1 x1 x1 x1 x1 x1;;
1021 float_interval_dih_x pp x1 x1 x1 x1 x1 x1;;
1022
1023 *)
1024
1025 (***********************************)
1026
1027
1028 end;;