formal_ineqs/taylor/m_taylor_arith.hl

   1 (* =========================================================== *)
   2 (* Formal arithmetic of taylor intervals                       *)
   3 (* Author: Alexey Solovyev                                     *)
   4 (* Date: 2012-10-27                                            *)
   5 (* =========================================================== *)
   6
   7 needs "taylor/m_taylor.hl";;
   8 needs "misc/vars.hl";;
   9
  10 module M_taylor_arith = struct
  11
  12 open Arith_misc;;
  13 open Arith_nat;;
  14 open Arith_float;;
  15 open More_float;;
  16 open Float_atn;;
  17 open Float_theory;;
  18 open M_taylor;;
  19 open Misc_vars;;
  20
  21
  22 let sqrt_tm = `sqrt` and
  23     atn_tm = `atn` and
  24     acs_tm = `acs`;;
  25
  26
  27 (*************************************)
  28
  29 let binary_beta_gen_eq f1_tm f2_tm x_var op_tm =
  30   let beta_tm1, beta_tm2 = mk_comb (f1_tm, x_var), mk_comb (f2_tm, x_var) in
  31   let beta_th1 = if is_abs f1_tm then BETA beta_tm1 else REFL beta_tm1 and
  32       beta_th2 = if is_abs f2_tm then BETA beta_tm2 else REFL beta_tm2 in
  33     ABS x_var (MK_COMB (AP_TERM op_tm beta_th1, beta_th2));;
  34
  35
  36 let unary_beta_gen_eq f_tm x_var op_tm =
  37   let beta_tm = mk_comb (f_tm, x_var) in
  38   let beta_th = if is_abs f_tm then BETA beta_tm else REFL beta_tm in
  39     ABS x_var (AP_TERM op_tm beta_th);;
  40
  41
  42
  43 let m_taylor_interval_norm th eq_th =
  44   let lhs1, d2f = dest_comb (concl th) in
  45   let lhs2, d1f = dest_comb lhs1 in
  46   let lhs3, d0f = dest_comb lhs2 in
  47   let lhs4, w = dest_comb lhs3 in
  48   let lhs5, y = dest_comb lhs4 in
  49   let lhs6, domain = dest_comb lhs5 in
  50   let m_taylor = rator lhs6 in
  51   let th0 = AP_TERM m_taylor eq_th in
  52   let th1 = AP_THM (AP_THM (AP_THM (AP_THM (AP_THM (AP_THM th0 domain) y) w) d0f) d1f) d2f in
  53     EQ_MP th1 th;;
  54
  55
  56
  57
  58 (*****************************************)
  59 (* dest_m_lin_approx *)
  60
  61
  62 let MK_M_LIN_APPROX' = (RULE o MATCH_MP EQ_IMP o SYM o SPEC_ALL) m_lin_approx;;
  63 let DEST_M_LIN_APPROX' = MY_RULE_NUM m_lin_approx;;
  64
  65
  66
  67 let m_lin_approx_components n m_lin_th =
  68   let f_tm, x_tm, f_bounds, d_bounds_list = dest_lin_approx (concl m_lin_th) in
  69   let ty = n_type_array.(n) in
  70   let f_var = mk_var ("f", type_of f_tm) in
  71   let x_var = mk_var ("x", type_of x_tm) in
  72   let th0 = (INST[f_tm, f_var; x_tm, x_var; f_bounds, f_bounds_var;
  73                   d_bounds_list, df_bounds_list_var] o inst_first_type_var ty) DEST_M_LIN_APPROX' in
  74   let th1 = EQ_MP th0 m_lin_th in
  75   let [r1; r2; r3] = CONJUNCTS th1 in
  76     r1, r2, r3;;
  77
  78
  79 (********************************)
  80 (* all_n manipulations *)
  81
  82 let ALL_N_EMPTY' = prove(`all_n n [] (s:num->A->bool)`, REWRITE_TAC[all_n]);;
  83 let ALL_N_CONS_IMP' = (MY_RULE o prove)(`SUC n = m /\ s n (x:A) ==>
  84     (all_n m t s <=> all_n n (CONS x t) s)`, SIMP_TAC[all_n]);;
  85 let ALL_N_CONS_EQ' = (MY_RULE o prove)(`SUC n = m ==>
  86     (all_n n (CONS x t) s <=> (s n (x:A) /\ all_n m t s))`, SIMP_TAC[all_n]);;
  87
  88 let dest_all_n all_n_tm =
  89   let ltm, s_tm = dest_comb all_n_tm in
  90   let ltm2, list_tm = dest_comb ltm in
  91     rand ltm2, list_tm, s_tm;;
  92
  93
  94 (* Splits `|- all_n n list s` into separate components.
  95    Also returns the list of SUC n = m theorems *)
  96 let all_n_components all_n_th =
  97   let n_tm, list_tm, s_tm = dest_all_n (concl all_n_th) in
  98   let list_ty = type_of list_tm in
  99   let ty = (hd o snd o dest_type) list_ty in
 100   let s_var = mk_var ("s", type_of s_tm) and
 101       x_var = mk_var ("x", ty) and
 102       t_var = mk_var ("t", list_ty) in
 103   let all_n_cons_th = (INST[s_tm, s_var] o INST_TYPE[ty, aty]) ALL_N_CONS_EQ' in
 104
 105   let rec get_components n_tm list_tm all_n_th =
 106     if is_const list_tm then [], []
 107     else
 108       let x_tm, t_tm = dest_cons list_tm in
 109       let suc_th = raw_suc_conv_hash (mk_comb (suc_op_num, n_tm)) in
 110       let m_tm = rand (concl suc_th) in
 111       let th0 = INST[n_tm, n_var_num; m_tm, m_var_num; x_tm, x_var; t_tm, t_var] all_n_cons_th in
 112       let th1 = MY_PROVE_HYP suc_th th0 in
 113       let th2 = EQ_MP th1 all_n_th in
 114       let snx_th, all_m_th = CONJUNCT1 th2, CONJUNCT2 th2 in
 115       let comps, suc_list = get_components m_tm t_tm all_m_th in
 116         snx_th :: comps, suc_th :: suc_list in
 117     get_components n_tm list_tm all_n_th;;
 118
 119
 120 (* Builds all_n from the given theorems and SUC n = m results *)
 121 let build_all_n ths suc_ths =
 122   (* The list ths should be not empty *)
 123   let tm0 = (concl o hd) ths in
 124   let lhs, rhs = dest_comb tm0 in
 125   let s_tm = rator lhs in
 126   let ty = type_of rhs in
 127   let list_ty = mk_type ("list", [ty]) in
 128   let s_var = mk_var ("s", type_of s_tm) and
 129       x_var = mk_var ("x", ty) and
 130       t_var = mk_var ("t", list_ty) in
 131   let m_tm = (rand o concl o hd) suc_ths in
 132
 133   let empty_th = (INST[s_tm, s_var; m_tm, n_var_num] o INST_TYPE[ty, aty]) ALL_N_EMPTY' in
 134   let cons_th = (INST[s_tm, s_var] o INST_TYPE[ty, aty]) ALL_N_CONS_IMP' in
 135
 136   let build suc_th s_th th =
 137     let t_tm = (rand o rator o concl) th in
 138     let x_tm = rand (concl s_th) in
 139     let lhs, m_tm = dest_eq (concl suc_th) in
 140     let n_tm = rand lhs in
 141     let th' = INST[n_tm, n_var_num; m_tm, m_var_num; x_tm, x_var; t_tm, t_var] cons_th in
 142       EQ_MP (MY_PROVE_HYP s_th (MY_PROVE_HYP suc_th th')) th in
 143
 144     rev_itlist2 build suc_ths ths empty_th;;
 145
 146
 147 (*************************)
 148
 149 (* Generates |- s D1 a1 /\ ... /\ s D_m a_m <=> all_n D1 [a1; ... ; a_m] s *)
 150 let gen_all_n_th m =
 151   let a_vars = map (fun i -> mk_var ("a"^string_of_int i, aty)) (1--m) in
 152   let list_tm = mk_list (a_vars, aty) in
 153   let all_tm = mk_comb (mk_binop `all_n : num -> (A)list -> (num -> A -> bool) -> bool` `1` list_tm,
 154                           `s : num -> A -> bool`) in
 155     (SYM o MY_RULE_NUM o CONV_RULE NUM_REDUCE_CONV) (REWRITE_CONV[all_n] all_tm);;
 156
 157 let all_n_array = Array.init (max_dim + 1) (fun i -> if i = 0 then TRUTH else gen_all_n_th i);;
 158
 159 (***)
 160
 161 let build2 ths =
 162   let n = length ths in
 163   let th0 = rev_itlist CONJ (tl ths) (hd ths) in
 164   let tm0 = (concl o hd) ths in
 165   let lhs, rhs = dest_comb tm0 in
 166   let a_tms = rev (map (rand o concl) ths) in
 167   let s_tm = rator lhs in
 168   let ty = type_of rhs in
 169   let s_var = mk_var ("s", type_of s_tm) and
 170       a_vars0 = map (fun i -> mk_var ("a"^string_of_int i, ty)) (1--n) in
 171
 172   let th1 = (INST[s_tm, s_var] o INST (zip a_tms a_vars0) o INST_TYPE[ty, aty]) all_n_array.(n) in
 173     EQ_MP th1 th0;;
 174
 175
 176
 177
 178 (************************)
 179
 180
 181 (* Constructs all_n n (map s list1) *)
 182 let eval_all_n all_n1_th beta_flag s =
 183   let ths1', suc_ths = all_n_components all_n1_th in
 184   let ths1 = if beta_flag then map MY_BETA_RULE ths1' else ths1' in
 185   let ths1, suc_ths = List.rev ths1, List.rev suc_ths in
 186   let ths = map s ths1 in
 187 (*    build_all_n ths suc_ths;; *)
 188     build2 ths;;
 189
 190
 191
 192 (* Constructs all_n n (map2 s list1 list2) *)
 193 let eval_all_n2 all_n1_th all_n2_th beta_flag s =
 194   let ths1', suc_ths = all_n_components all_n1_th in
 195   let ths2', _ = all_n_components all_n2_th in
 196   let ths1, ths2 =
 197     if beta_flag then map MY_BETA_RULE ths1', map MY_BETA_RULE ths2' else ths1', ths2' in
 198
 199   let ths1, ths2, suc_ths = List.rev ths1, List.rev ths2, List.rev suc_ths in
 200   let ths = map2 s ths1 ths2 in
 201 (*    build_all_n ths suc_ths;; *)
 202     build2 ths;;
 203
 204
 205
 206 (***************************************)
 207 (* eval_m_taylor_add *)
 208
 209
 210 let SECOND_BOUNDED' = MY_RULE_NUM second_bounded;;
 211
 212 let dest_second_bounded tm =
 213   let ltm, dd = dest_comb tm in
 214   let ltm2, domain = dest_comb ltm in
 215     rand ltm2, domain, dd;;
 216
 217 let second_bounded_components n th =
 218   let f_tm, domain_tm, dd_tm = dest_second_bounded (concl th) in
 219   let x_var = mk_var ("x", n_vector_type_array.(n)) in
 220   let th0 = (INST[f_tm, mk_var ("f", type_of f_tm);
 221                  domain_tm, mk_var ("domain", type_of domain_tm);
 222                  dd_tm, dd_bounds_list_var] o inst_first_type_var n_type_array.(n)) SECOND_BOUNDED' in
 223     UNDISCH (SPEC x_var (EQ_MP th0 th));;
 224
 225
 226
 227
 228 let MK_M_TAYLOR_ADD' = (MY_RULE_NUM o prove)
 229   (`m_cell_domain domain (y:real^N) w ==>
 230      diff2c_domain domain f ==>
 231      diff2c_domain domain g ==>
 232      interval_arith (f y + g y) bounds ==>
 233      all_n 1 d_bounds_list (\i int. interval_arith (partial i f y + partial i g y) int) ==>
 234      (!x. x IN interval [domain] ==> all_n 1 dd_bounds_list (\i list_i. all_n 1 list_i
 235                                (\j int. interval_arith (partial2 j i f x + partial2 j i g x) int))) ==>
 236      m_taylor_interval (\x. f x + g x) domain y w bounds d_bounds_list dd_bounds_list`,
 237    REWRITE_TAC[m_taylor_interval; m_lin_approx; second_bounded] THEN
 238      REPEAT DISCH_TAC THEN
 239      SUBGOAL_THEN `lift o f differentiable at (y:real^N) /\ lift o g differentiable at y` ASSUME_TAC THENL
 240      [
 241        UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN
 242          UNDISCH_TAC `diff2c_domain domain (g:real^N->real)` THEN
 243          REWRITE_TAC[diff2c_domain_alt; diff2_domain] THEN STRIP_TAC THEN STRIP_TAC THEN
 244          REPEAT (new_rewrite [] [] diff2_imp_diff) THEN REWRITE_TAC[] THEN
 245          FIRST_X_ASSUM MATCH_MP_TAC THEN
 246          MATCH_MP_TAC y_in_domain THEN EXISTS_TAC `w:real^N` THEN ASM_REWRITE_TAC[];
 247        ALL_TAC
 248      ] THEN
 249
 250      REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[] THENL
 251      [
 252        new_rewrite [] [] diff2c_domain_add THEN ASM_REWRITE_TAC[];
 253        REWRITE_TAC[f_lift_add] THEN
 254          new_rewrite [] [] DIFFERENTIABLE_ADD THEN
 255          ASM_REWRITE_TAC[ETA_AX];
 256        ASM_SIMP_TAC[partial_add];
 257        ALL_TAC
 258      ] THEN
 259
 260      UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN
 261      UNDISCH_TAC `diff2c_domain domain (g:real^N->real)` THEN
 262      REWRITE_TAC[diff2c_domain_alt; diff2_domain] THEN
 263      REPEAT (DISCH_THEN (MP_TAC o SPEC `x:real^N` o CONJUNCT1) THEN ASM_REWRITE_TAC[] THEN DISCH_TAC) THEN
 264      ASM_SIMP_TAC[second_partial_add]);;
 265
 266
 267 (*************************)
 268
 269
 270 let add_partial_lemma' = prove(`interval_arith (partial i f (x:real^N) + partial i g x) int <=>
 271                                  (\i int. interval_arith (partial i f x + partial i g x) int) i int`,
 272                                REWRITE_TAC[]);;
 273
 274 let add_second_lemma' = prove(`interval_arith (partial2 j i f (x:real^N) + partial2 j i g x) int <=>
 275                          (\j int. interval_arith (partial2 j i f x + partial2 j i g x) int) j int`,
 276                                REWRITE_TAC[]);;
 277
 278 let add_second_lemma'' = (NUMERALS_TO_NUM o prove)(`all_n 1 list
 279                 (\j int. interval_arith (partial2 j i f (x:real^N) + partial2 j i g x) int) <=>
 280                 (\i list. all_n 1 list
 281                    (\j int. interval_arith (partial2 j i f x + partial2 j i g x) int)) i list`,
 282                                REWRITE_TAC[]);;
 283
 284
 285
 286
 287 let eval_m_taylor_add n p_lin p_second taylor1_th taylor2_th =
 288   let domain_th, diff2_f1_th, lin1_th, second1_th = dest_m_taylor_thms n taylor1_th in
 289   let _, diff2_f2_th, lin2_th, second2_th = dest_m_taylor_thms n taylor2_th in
 290   let f1_tm = (rand o concl) diff2_f1_th and
 291       f2_tm = (rand o concl) diff2_f2_th in
 292   let domain_tm, y_tm, w_tm = dest_m_cell_domain (concl domain_th) in
 293   let ty = type_of y_tm in
 294
 295   let x_var = mk_var ("x", ty) and
 296       y_var = mk_var ("y", ty) and
 297       w_var = mk_var ("w", ty) and
 298       f_var = mk_var ("f", type_of f1_tm) and
 299       g_var = mk_var ("g", type_of f2_tm) and
 300       domain_var = mk_var ("domain", type_of domain_tm) in
 301
 302   let _, bounds1_th, df1_th = m_lin_approx_components n lin1_th and
 303       _, bounds2_th, df2_th = m_lin_approx_components n lin2_th in
 304
 305   let bounds_th = float_interval_add p_lin bounds1_th bounds2_th in
 306   let bounds_tm = (rand o concl) bounds_th in
 307
 308   let add_lemma0 = (INST[f1_tm, f_var; f2_tm, g_var; y_tm, x_var] o
 309                       INST_TYPE[n_type_array.(n), nty]) add_partial_lemma' in
 310
 311   let add th1 th2 =
 312     let add_th = float_interval_add p_lin th1 th2 in
 313     let int_tm = rand (concl add_th) and
 314         i_tm = (rand o rator o rator o lhand) (concl th1) in
 315     let th0 = INST[i_tm, i_var_num; int_tm, int_var] add_lemma0 in
 316       EQ_MP th0 add_th in
 317
 318   let df_th = eval_all_n2 df1_th df2_th true add in
 319   let d_bounds_list = (rand o rator o concl) df_th in
 320
 321
 322   let dd1 = second_bounded_components n second1_th in
 323   let dd2 = second_bounded_components n second2_th in
 324
 325   let add_second_lemma0 = (INST[f1_tm, f_var; f2_tm, g_var] o
 326                              INST_TYPE[n_type_array.(n), nty]) add_second_lemma' in
 327
 328   let add_second_lemma1 = (INST[f1_tm, f_var; f2_tm, g_var] o
 329                              INST_TYPE[n_type_array.(n), nty]) add_second_lemma'' in
 330
 331
 332   let add_second2 th1 th2 =
 333     let i_tm = (rand o rator o concl) th1 in
 334     let th1, th2 = MY_BETA_RULE th1, MY_BETA_RULE th2 in
 335     let lemma = INST[i_tm, i_var_num] add_second_lemma0 in
 336     let add_second th1 th2 =
 337       let add_th = float_interval_add p_second th1 th2 in
 338       let int_tm = rand (concl add_th) and
 339           j_tm = (rand o rator o rator o rator o lhand) (concl th1) in
 340       let th0 = INST[j_tm, j_var_num; int_tm, int_var] lemma in
 341         EQ_MP th0 add_th in
 342     let add_th = eval_all_n2 th1 th2 true add_second in
 343     let list_tm = (rand o rator o concl) add_th in
 344     let lemma1 = INST[i_tm, i_var_num; list_tm, list_var_real_pair] add_second_lemma1 in
 345       EQ_MP lemma1 add_th in
 346
 347
 348   let dd_th0 = eval_all_n2 dd1 dd2 false add_second2 in
 349   let dd_list = (rand o rator o concl) dd_th0 in
 350   let dd_th = GEN x_var (DISCH_ALL dd_th0) in
 351
 352   let th = (MY_PROVE_HYP dd_th o MY_PROVE_HYP diff2_f1_th o MY_PROVE_HYP diff2_f2_th o
 353               MY_PROVE_HYP bounds_th o MY_PROVE_HYP df_th o MY_PROVE_HYP domain_th o
 354               INST[f1_tm, f_var; f2_tm, g_var;
 355                    domain_tm, domain_var; y_tm, y_var; w_tm, w_var;
 356                    bounds_tm, bounds_var; d_bounds_list, d_bounds_list_var;
 357                    dd_list, dd_bounds_list_var] o
 358               INST_TYPE[n_type_array.(n), nty]) MK_M_TAYLOR_ADD' in
 359   let eq_th = binary_beta_gen_eq f1_tm f2_tm x_var add_op_real in
 360     m_taylor_interval_norm th eq_th;;
 361
 362
 363
 364 (***************************************)
 365 (* eval_m_taylor_sub *)
 366
 367
 368 let MK_M_TAYLOR_SUB' = (MY_RULE_NUM o prove)
 369   (`m_cell_domain domain (y:real^N) w ==>
 370      diff2c_domain domain f ==>
 371      diff2c_domain domain g ==>
 372      interval_arith (f y - g y) bounds ==>
 373      all_n 1 d_bounds_list (\i int. interval_arith (partial i f y - partial i g y) int) ==>
 374      (!x. x IN interval [domain] ==> all_n 1 dd_bounds_list (\i list_i. all_n 1 list_i
 375                                (\j int. interval_arith (partial2 j i f x - partial2 j i g x) int))) ==>
 376      m_taylor_interval (\x. f x - g x) domain y w bounds d_bounds_list dd_bounds_list`,
 377    REWRITE_TAC[m_taylor_interval; m_lin_approx; second_bounded] THEN
 378      REPEAT DISCH_TAC THEN
 379      SUBGOAL_THEN `lift o f differentiable at (y:real^N) /\ lift o g differentiable at y` ASSUME_TAC THENL
 380      [
 381        UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN
 382          UNDISCH_TAC `diff2c_domain domain (g:real^N->real)` THEN
 383          REWRITE_TAC[diff2c_domain_alt; diff2_domain] THEN STRIP_TAC THEN STRIP_TAC THEN
 384          REPEAT (new_rewrite [] [] diff2_imp_diff) THEN REWRITE_TAC[] THEN
 385          FIRST_X_ASSUM MATCH_MP_TAC THEN
 386          MATCH_MP_TAC y_in_domain THEN EXISTS_TAC `w:real^N` THEN ASM_REWRITE_TAC[];
 387        ALL_TAC
 388      ] THEN
 389
 390      REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[] THENL
 391      [
 392        new_rewrite [] [] diff2c_domain_sub THEN ASM_REWRITE_TAC[];
 393        REWRITE_TAC[f_lift_sub] THEN
 394          new_rewrite [] [] DIFFERENTIABLE_SUB THEN
 395          ASM_REWRITE_TAC[ETA_AX];
 396        ASM_SIMP_TAC[partial_sub];
 397        ALL_TAC
 398      ] THEN
 399
 400      UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN
 401      UNDISCH_TAC `diff2c_domain domain (g:real^N->real)` THEN
 402      REWRITE_TAC[diff2c_domain_alt; diff2_domain] THEN
 403      REPEAT (DISCH_THEN (MP_TAC o SPEC `x:real^N` o CONJUNCT1) THEN ASM_REWRITE_TAC[] THEN DISCH_TAC) THEN
 404      ASM_SIMP_TAC[second_partial_sub]);;
 405
 406
 407 (*************************)
 408
 409
 410 let sub_partial_lemma' = prove(`interval_arith (partial i f (x:real^N) - partial i g x) int <=>
 411                                  (\i int. interval_arith (partial i f x - partial i g x) int) i int`,
 412                                REWRITE_TAC[]);;
 413
 414 let sub_second_lemma' = prove(`interval_arith (partial2 j i f (x:real^N) - partial2 j i g x) int <=>
 415                          (\j int. interval_arith (partial2 j i f x - partial2 j i g x) int) j int`,
 416                                REWRITE_TAC[]);;
 417
 418 let sub_second_lemma'' = (NUMERALS_TO_NUM o prove)(`all_n 1 list
 419                 (\j int. interval_arith (partial2 j i f (x:real^N) - partial2 j i g x) int) <=>
 420                 (\i list. all_n 1 list
 421                    (\j int. interval_arith (partial2 j i f x - partial2 j i g x) int)) i list`,
 422                                REWRITE_TAC[]);;
 423
 424
 425
 426
 427 let eval_m_taylor_sub n p_lin p_second taylor1_th taylor2_th =
 428   let domain_th, diff2_f1_th, lin1_th, second1_th = dest_m_taylor_thms n taylor1_th in
 429   let _, diff2_f2_th, lin2_th, second2_th = dest_m_taylor_thms n taylor2_th in
 430   let f1_tm = (rand o concl) diff2_f1_th and
 431       f2_tm = (rand o concl) diff2_f2_th in
 432   let domain_tm, y_tm, w_tm = dest_m_cell_domain (concl domain_th) in
 433   let ty = type_of y_tm in
 434
 435   let x_var = mk_var ("x", ty) and
 436       y_var = mk_var ("y", ty) and
 437       w_var = mk_var ("w", ty) and
 438       f_var = mk_var ("f", type_of f1_tm) and
 439       g_var = mk_var ("g", type_of f2_tm) and
 440       domain_var = mk_var ("domain", type_of domain_tm) in
 441
 442   let _, bounds1_th, df1_th = m_lin_approx_components n lin1_th and
 443       _, bounds2_th, df2_th = m_lin_approx_components n lin2_th in
 444
 445   let bounds_th = float_interval_sub p_lin bounds1_th bounds2_th in
 446   let bounds_tm = (rand o concl) bounds_th in
 447
 448   let sub_lemma0 = (INST[f1_tm, f_var; f2_tm, g_var; y_tm, x_var] o
 449                       INST_TYPE[n_type_array.(n), nty]) sub_partial_lemma' in
 450
 451   let sub th1 th2 =
 452     let sub_th = float_interval_sub p_lin th1 th2 in
 453     let int_tm = rand (concl sub_th) and
 454         i_tm = (rand o rator o rator o lhand) (concl th1) in
 455     let th0 = INST[i_tm, i_var_num; int_tm, int_var] sub_lemma0 in
 456       EQ_MP th0 sub_th in
 457
 458   let df_th = eval_all_n2 df1_th df2_th true sub in
 459   let d_bounds_list = (rand o rator o concl) df_th in
 460
 461
 462   let dd1 = second_bounded_components n second1_th in
 463   let dd2 = second_bounded_components n second2_th in
 464
 465   let sub_second_lemma0 = (INST[f1_tm, f_var; f2_tm, g_var] o
 466                              INST_TYPE[n_type_array.(n), nty]) sub_second_lemma' in
 467
 468   let sub_second_lemma1 = (INST[f1_tm, f_var; f2_tm, g_var] o
 469                              INST_TYPE[n_type_array.(n), nty]) sub_second_lemma'' in
 470
 471
 472   let sub_second2 th1 th2 =
 473     let i_tm = (rand o rator o concl) th1 in
 474     let th1, th2 = MY_BETA_RULE th1, MY_BETA_RULE th2 in
 475     let lemma = INST[i_tm, i_var_num] sub_second_lemma0 in
 476     let sub_second th1 th2 =
 477       let sub_th = float_interval_sub p_second th1 th2 in
 478       let int_tm = rand (concl sub_th) and
 479           j_tm = (rand o rator o rator o rator o lhand) (concl th1) in
 480       let th0 = INST[j_tm, j_var_num; int_tm, int_var] lemma in
 481         EQ_MP th0 sub_th in
 482     let sub_th = eval_all_n2 th1 th2 true sub_second in
 483     let list_tm = (rand o rator o concl) sub_th in
 484     let lemma1 = INST[i_tm, i_var_num; list_tm, list_var_real_pair] sub_second_lemma1 in
 485       EQ_MP lemma1 sub_th in
 486
 487
 488   let dd_th0 = eval_all_n2 dd1 dd2 false sub_second2 in
 489   let dd_list = (rand o rator o concl) dd_th0 in
 490   let dd_th = GEN x_var (DISCH_ALL dd_th0) in
 491
 492   let th = (MY_PROVE_HYP dd_th o MY_PROVE_HYP diff2_f1_th o MY_PROVE_HYP diff2_f2_th o
 493               MY_PROVE_HYP bounds_th o MY_PROVE_HYP df_th o MY_PROVE_HYP domain_th o
 494               INST[f1_tm, f_var; f2_tm, g_var;
 495                    domain_tm, domain_var; y_tm, y_var; w_tm, w_var;
 496                    bounds_tm, bounds_var; d_bounds_list, d_bounds_list_var;
 497                    dd_list, dd_bounds_list_var] o
 498               INST_TYPE[n_type_array.(n), nty]) MK_M_TAYLOR_SUB' in
 499   let eq_th = binary_beta_gen_eq f1_tm f2_tm x_var sub_op_real in
 500     m_taylor_interval_norm th eq_th;;
 501
 502
 503 (*******************************************************)
 504
 505
 506 (***************************************)
 507 (* eval_m_taylor_mul *)
 508
 509
 510 let MK_M_TAYLOR_MUL' = (MY_RULE_NUM o prove)
 511   (`m_cell_domain domain (y:real^N) w ==>
 512      diff2c_domain domain f ==>
 513      diff2c_domain domain g ==>
 514      interval_arith (f y * g y) bounds ==>
 515      all_n 1 d_bounds_list (\i int. interval_arith (partial i f y * g y + f y * partial i g y) int) ==>
 516      (!x. x IN interval [domain] ==> all_n 1 dd_bounds_list (\i list_i. all_n 1 list_i
 517         (\j int. interval_arith ((partial2 j i f x * g x + partial i f x * partial j g x) +
 518                                 partial j f x * partial i g x + f x * partial2 j i g x) int))) ==>
 519      m_taylor_interval (\x. f x * g x) domain y w bounds d_bounds_list dd_bounds_list`,
 520    REWRITE_TAC[m_taylor_interval; m_lin_approx; second_bounded] THEN
 521      REPEAT DISCH_TAC THEN
 522      SUBGOAL_THEN `lift o f differentiable at (y:real^N) /\ lift o g differentiable at y` ASSUME_TAC THENL
 523      [
 524        UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN
 525          UNDISCH_TAC `diff2c_domain domain (g:real^N->real)` THEN
 526          REWRITE_TAC[diff2c_domain_alt; diff2_domain] THEN STRIP_TAC THEN STRIP_TAC THEN
 527          REPEAT (new_rewrite [] [] diff2_imp_diff) THEN REWRITE_TAC[] THEN
 528          FIRST_X_ASSUM MATCH_MP_TAC THEN
 529          MATCH_MP_TAC y_in_domain THEN EXISTS_TAC `w:real^N` THEN ASM_REWRITE_TAC[];
 530        ALL_TAC
 531      ] THEN
 532
 533      REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[] THENL
 534      [
 535        new_rewrite [] [] diff2c_domain_mul THEN ASM_REWRITE_TAC[];
 536          new_rewrite [] [] differentiable_mul THEN ASM_REWRITE_TAC[];
 537        ASM_SIMP_TAC[partial_mul];
 538        ALL_TAC
 539      ] THEN
 540
 541      UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN
 542      UNDISCH_TAC `diff2c_domain domain (g:real^N->real)` THEN
 543      REWRITE_TAC[diff2c_domain_alt; diff2_domain] THEN
 544      REPEAT (DISCH_THEN (MP_TAC o SPEC `x:real^N` o CONJUNCT1) THEN ASM_REWRITE_TAC[] THEN DISCH_TAC) THEN
 545      ASM_SIMP_TAC[second_partial_mul]);;
 546
 547
 548 (*************************)
 549
 550
 551 let mul_partial_lemma' =
 552   prove(`interval_arith (partial i f (y:real^N) * g y + f y * partial i g y) int <=>
 553           (\i int. interval_arith (partial i f y * g y + f y * partial i g y) int) i int`,
 554         REWRITE_TAC[]);;
 555
 556 let mul_second_lemma' =
 557   prove(`interval_arith ((partial2 j i f x * g (x:real^N) + partial i f x * partial j g x) +
 558                                 partial j f x * partial i g x + f x * partial2 j i g x) int <=>
 559           (\j int. interval_arith ((partial2 j i f x * g x + partial i f x * partial j g x) +
 560                                 partial j f x * partial i g x + f x * partial2 j i g x) int) j int`,
 561         REWRITE_TAC[]);;
 562
 563
 564 let mul_second_lemma'' = (NUMERALS_TO_NUM o prove)
 565   (`all_n 1 list (\j int. interval_arith ((partial2 j i f x * g x + partial i f x * partial j g x) +
 566                 partial j f x * partial i g x + f (x:real^N) * partial2 j i g x) int) <=>
 567      (\i list. all_n 1 list
 568         (\j int. interval_arith ((partial2 j i f x * g x + partial i f x * partial j g x) +
 569                                    partial j f x * partial i g x + f x * partial2 j i g x) int)) i list`,
 570                                REWRITE_TAC[]);;
 571
 572
 573
 574
 575 let eval_m_taylor_mul n p_lin p_second taylor1_th taylor2_th =
 576   let domain_th, diff2_f1_th, lin1_th, second1_th = dest_m_taylor_thms n taylor1_th and
 577       _, diff2_f2_th, lin2_th, second2_th = dest_m_taylor_thms n taylor2_th in
 578   let f1_tm = (rand o concl) diff2_f1_th and
 579       f2_tm = (rand o concl) diff2_f2_th in
 580   let domain_tm, y_tm, w_tm = dest_m_cell_domain (concl domain_th) in
 581   let ty = type_of y_tm in
 582
 583   let x_var = mk_var ("x", ty) and
 584       y_var = mk_var ("y", ty) and
 585       w_var = mk_var ("w", ty) and
 586       f_var = mk_var ("f", type_of f1_tm) and
 587       g_var = mk_var ("g", type_of f2_tm) and
 588       domain_var = mk_var ("domain", type_of domain_tm) in
 589
 590   let _, bounds1_th, df1_th = m_lin_approx_components n lin1_th and
 591       _, bounds2_th, df2_th = m_lin_approx_components n lin2_th in
 592
 593   let bounds_th = float_interval_mul p_lin bounds1_th bounds2_th in
 594   let bounds_tm = (rand o concl) bounds_th in
 595
 596   let mul_lemma0 = (INST[f1_tm, f_var; f2_tm, g_var; y_tm, y_var] o
 597                       INST_TYPE[n_type_array.(n), nty]) mul_partial_lemma' in
 598
 599   let mul th1 th2 =
 600     let mul_th =
 601       let ( * ), ( + ) = float_interval_mul p_lin, float_interval_add p_lin in
 602         th1 * bounds2_th + bounds1_th * th2 in
 603     let int_tm = rand (concl mul_th) in
 604     let i_tm = (rand o rator o rator o lhand) (concl th1) in
 605     let th0 = INST[i_tm, i_var_num; int_tm, int_var] mul_lemma0 in
 606       EQ_MP th0 mul_th in
 607
 608   let df_th = eval_all_n2 df1_th df2_th true mul in
 609   let d_bounds_list = (rand o rator o concl) df_th in
 610
 611
 612   let dd1 = second_bounded_components n second1_th in
 613   let dd2 = second_bounded_components n second2_th in
 614
 615   let mul_second_lemma0 = (INST[f1_tm, f_var; f2_tm, g_var] o
 616                              INST_TYPE[n_type_array.(n), nty]) mul_second_lemma' in
 617
 618   let mul_second_lemma1 = (INST[f1_tm, f_var; f2_tm, g_var] o
 619                              INST_TYPE[n_type_array.(n), nty]) mul_second_lemma'' in
 620
 621   let undisch = UNDISCH o SPEC x_var in
 622
 623   let d1_bounds = map (fun i ->
 624                          let th0 = eval_m_taylor_partial_bound n p_second i taylor1_th in
 625                            undisch th0) (1--n) in
 626
 627   let d2_bounds = map (fun i ->
 628                          let th0 = eval_m_taylor_partial_bound n p_second i taylor2_th in
 629                            undisch th0) (1--n) in
 630
 631   let f1_bound = undisch (eval_m_taylor_bound n p_second taylor1_th) and
 632       f2_bound = undisch (eval_m_taylor_bound n p_second taylor2_th) in
 633
 634   let mul_second2 th1 th2 =
 635     let i_tm = (rand o rator o concl) th1 in
 636     let i_int = (Num.int_of_num o raw_dest_hash) i_tm in
 637     let di1 = List.nth d1_bounds (i_int - 1) and
 638         di2 = List.nth d2_bounds (i_int - 1) in
 639     let th1, th2 = MY_BETA_RULE th1, MY_BETA_RULE th2 in
 640     let lemma = INST[i_tm, i_var_num] mul_second_lemma0 in
 641     let mul_second th1 th2 =
 642       let j_tm = (rand o rator o rator o rator o lhand) (concl th1) in
 643       let j_int = (Num.int_of_num o raw_dest_hash) j_tm in
 644       let dj1 = List.nth d1_bounds (j_int - 1) and
 645           dj2 = List.nth d2_bounds (j_int - 1) in
 646
 647       let mul_th =
 648         let ( * ), ( + ) = float_interval_mul p_second, float_interval_add p_second in
 649           (th1 * f2_bound + di1 * dj2) + (dj1 * di2 + f1_bound * th2) in
 650
 651       let int_tm = rand (concl mul_th) in
 652       let th0 = INST[j_tm, j_var_num; int_tm, int_var] lemma in
 653         EQ_MP th0 mul_th in
 654
 655     let mul_th = eval_all_n2 th1 th2 true mul_second in
 656     let list_tm = (rand o rator o concl) mul_th in
 657     let lemma1 = INST[i_tm, i_var_num; list_tm, list_var_real_pair] mul_second_lemma1 in
 658       EQ_MP lemma1 mul_th in
 659
 660
 661   let dd_th0 = eval_all_n2 dd1 dd2 false mul_second2 in
 662   let dd_list = (rand o rator o concl) dd_th0 in
 663   let dd_th = GEN x_var (DISCH_ALL dd_th0) in
 664
 665   let th = (MY_PROVE_HYP dd_th o MY_PROVE_HYP diff2_f1_th o MY_PROVE_HYP diff2_f2_th o
 666               MY_PROVE_HYP bounds_th o MY_PROVE_HYP df_th o MY_PROVE_HYP domain_th o
 667               INST[f1_tm, f_var; f2_tm, g_var;
 668                    domain_tm, domain_var; y_tm, y_var; w_tm, w_var;
 669                    bounds_tm, bounds_var; d_bounds_list, d_bounds_list_var;
 670                    dd_list, dd_bounds_list_var] o
 671               INST_TYPE[n_type_array.(n), nty]) MK_M_TAYLOR_MUL' in
 672   let eq_th = binary_beta_gen_eq f1_tm f2_tm x_var mul_op_real in
 673     m_taylor_interval_norm th eq_th;;
 674
 675
 676
 677
 678
 679
 680 (*******************************************************)
 681
 682 (* neg, inv, sqrt, atn, acs *)
 683
 684
 685
 686 let partial_uni_compose' =
 687   REWRITE_RULE[SWAP_FORALL_THM; GSYM RIGHT_IMP_FORALL_THM] partial_uni_compose;;
 688 let second_partial_uni_compose' =
 689   REWRITE_RULE[SWAP_FORALL_THM; GSYM RIGHT_IMP_FORALL_THM] second_partial_uni_compose;;
 690
 691
 692 (* neg *)
 693 let MK_M_TAYLOR_NEG' = (MY_RULE_FLOAT o prove)
 694   (`m_cell_domain domain (y:real^N) w ==>
 695      diff2c_domain domain f ==>
 696      interval_arith (-- (f y)) bounds ==>
 697      all_n 1 d_bounds_list (\i int. interval_arith (-- partial i f y) int) ==>
 698      (!x. x IN interval [domain] ==> all_n 1 dd_bounds_list (\i list_i. all_n 1 list_i
 699         (\j int. interval_arith (-- partial2 j i f x) int))) ==>
 700      m_taylor_interval (\x. -- (f x)) domain y w bounds d_bounds_list dd_bounds_list`,
 701    REWRITE_TAC[m_taylor_interval; m_lin_approx; second_bounded] THEN
 702      REPEAT DISCH_TAC THEN
 703      SUBGOAL_THEN `lift o f differentiable at (y:real^N)` ASSUME_TAC THENL
 704      [
 705        UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN
 706          REWRITE_TAC[diff2c_domain_alt; diff2_domain] THEN STRIP_TAC THEN
 707          REPEAT (new_rewrite [] [] diff2_imp_diff) THEN REWRITE_TAC[] THEN
 708          FIRST_X_ASSUM MATCH_MP_TAC THEN
 709          MATCH_MP_TAC y_in_domain THEN EXISTS_TAC `w:real^N` THEN ASM_REWRITE_TAC[];
 710        ALL_TAC
 711      ] THEN
 712
 713      REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[] THENL
 714      [
 715        new_rewrite [] [] diff2c_domain_neg THEN ASM_REWRITE_TAC[];
 716        REWRITE_TAC[f_lift_neg] THEN
 717          new_rewrite [] [] DIFFERENTIABLE_NEG THEN
 718          ASM_REWRITE_TAC[ETA_AX];
 719        ASM_SIMP_TAC[partial_neg];
 720        ALL_TAC
 721      ] THEN
 722
 723      UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN
 724      REWRITE_TAC[diff2c_domain_alt; diff2_domain] THEN
 725      DISCH_THEN (MP_TAC o SPEC `x:real^N` o CONJUNCT1) THEN ASM_REWRITE_TAC[] THEN DISCH_TAC THEN
 726      ASM_SIMP_TAC[second_partial_neg]);;
 727
 728
 729 (* inv *)
 730 let MK_M_TAYLOR_INV' = (UNDISCH_ALL o PURE_REWRITE_RULE[float2_eq] o DISCH_ALL o
 731                           MY_RULE_FLOAT o prove)
 732   (`m_cell_domain domain (y:real^N) w ==>
 733      (!x. x IN interval [domain] ==> interval_arith (f x) f_bounds) ==>
 734      interval_not_zero f_bounds ==>
 735      diff2c_domain domain f ==>
 736      interval_arith (inv (f y)) bounds ==>
 737      all_n 1 d_bounds_list (\i int. interval_arith (--inv (f y * f y) * partial i f y) int) ==>
 738      (!x. x IN interval [domain] ==> all_n 1 dd_bounds_list (\i list_i. all_n 1 list_i
 739         (\j int. interval_arith (((&2 * inv (f x * f x * f x)) * partial j f x) * partial i f x -
 740                                    inv (f x * f x) * partial2 j i f x) int))) ==>
 741      m_taylor_interval (\x. inv (f x)) domain y w bounds d_bounds_list dd_bounds_list`,
 742    REWRITE_TAC[m_taylor_interval; m_lin_approx; second_bounded; ETA_AX] THEN
 743      REPEAT DISCH_TAC THEN
 744      SUBGOAL_THEN `lift o f differentiable at (y:real^N)` ASSUME_TAC THENL
 745      [
 746        UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN
 747          REWRITE_TAC[diff2c_domain_alt; diff2_domain] THEN STRIP_TAC THEN
 748          new_rewrite [] [] diff2_imp_diff THEN REWRITE_TAC[] THEN
 749          FIRST_X_ASSUM MATCH_MP_TAC THEN
 750          MATCH_MP_TAC y_in_domain THEN EXISTS_TAC `w:real^N` THEN ASM_REWRITE_TAC[];
 751        ALL_TAC
 752      ] THEN
 753
 754      SUBGOAL_THEN `!x:real^N. x IN interval [domain] ==> ~(f x = &0)` ASSUME_TAC THENL
 755      [
 756        GEN_TAC THEN DISCH_TAC THEN
 757          apply_tac interval_arith_not_zero THEN
 758          EXISTS_TAC `f_bounds:real#real` THEN
 759          ASM_REWRITE_TAC[] THEN
 760          FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[];
 761        ALL_TAC
 762      ] THEN
 763
 764      SUBGOAL_THEN `~(f (y:real^N) = &0)` ASSUME_TAC THENL
 765      [
 766        FIRST_X_ASSUM MATCH_MP_TAC THEN
 767          MATCH_MP_TAC y_in_domain THEN EXISTS_TAC `w:real^N` THEN
 768          ASM_REWRITE_TAC[];
 769        ALL_TAC
 770      ] THEN
 771
 772      REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[] THENL
 773      [
 774        UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN
 775          REWRITE_TAC[diff2c_domain] THEN REPEAT STRIP_TAC THEN
 776          ONCE_REWRITE_TAC[GSYM o_THM] THEN REWRITE_TAC[ETA_AX] THEN
 777          apply_tac diff2c_inv_compose THEN CONJ_TAC THEN FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[];
 778        new_rewrite [] [`inv _`] (GSYM o_THM) THEN REWRITE_TAC[ETA_AX] THEN
 779          apply_tac diff_uni_compose THEN ASM_REWRITE_TAC[] THEN
 780          MATCH_MP_TAC REAL_DIFFERENTIABLE_AT_INV THEN ASM_REWRITE_TAC[];
 781        new_rewrite [] [`inv _`] (GSYM o_THM) THEN REWRITE_TAC[ETA_AX] THEN
 782          MP_TAC (ISPECL [`y:real^N`; `f:real^N->real`] partial_uni_compose') THEN
 783          ASM_REWRITE_TAC[] THEN
 784          DISCH_THEN (MP_TAC o SPEC `inv`) THEN
 785          ANTS_TAC THENL
 786          [
 787            MATCH_MP_TAC REAL_DIFFERENTIABLE_AT_INV THEN ASM_REWRITE_TAC[];
 788            ALL_TAC
 789          ] THEN
 790          ASM_SIMP_TAC[derivative_inv];
 791        ALL_TAC
 792      ] THEN
 793
 794      UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN
 795      REWRITE_TAC[diff2c_domain_alt; diff2_domain] THEN
 796      DISCH_THEN (MP_TAC o SPEC `x:real^N` o CONJUNCT1) THEN ASM_REWRITE_TAC[] THEN DISCH_TAC THEN
 797      MP_TAC (ISPECL [`x:real^N`; `f:real^N->real`] second_partial_uni_compose') THEN
 798      ASM_REWRITE_TAC[] THEN
 799      DISCH_THEN (MP_TAC o SPEC `inv`) THEN
 800      ANTS_TAC THENL
 801      [
 802        MATCH_MP_TAC diff2_inv THEN
 803          FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[];
 804        ALL_TAC
 805      ] THEN
 806
 807      new_rewrite [] [`inv _`] (GSYM o_THM) THEN REWRITE_TAC[ETA_AX] THEN
 808      ASM_SIMP_TAC[second_derivative_inv; derivative_inv; REAL_MUL_LNEG; GSYM real_sub] THEN
 809      ASM_SIMP_TAC[REAL_ARITH `a pow 3 = a * a * a`]);;
 810
 811
 812 (* sqrt *)
 813 let MK_M_TAYLOR_SQRT' = (UNDISCH_ALL o PURE_REWRITE_RULE[float2_eq; float4_eq] o
 814                            DISCH_ALL o MY_RULE_FLOAT o prove)
 815   (`m_cell_domain domain (y:real^N) w ==>
 816      (!x. x IN interval [domain] ==> interval_arith (f x) f_bounds) ==>
 817      interval_pos f_bounds ==>
 818      diff2c_domain domain f ==>
 819      interval_arith (sqrt (f y)) bounds ==>
 820      all_n 1 d_bounds_list (\i int. interval_arith (inv (&2 * sqrt (f y)) * partial i f y) int) ==>
 821      (!x. x IN interval [domain] ==> all_n 1 dd_bounds_list (\i list_i. all_n 1 list_i
 822         (\j int. interval_arith ((--inv ((&2 * sqrt (f x)) * (&2 * f x)) * partial j f x) * partial i f x +
 823                                    inv (&2 * sqrt (f x)) * partial2 j i f x) int))) ==>
 824      m_taylor_interval (\x. sqrt (f x)) domain y w bounds d_bounds_list dd_bounds_list`,
 825    REWRITE_TAC[m_taylor_interval; m_lin_approx; second_bounded; ETA_AX] THEN
 826      REPEAT DISCH_TAC THEN
 827      SUBGOAL_THEN `lift o f differentiable at (y:real^N)` ASSUME_TAC THENL
 828      [
 829        UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN
 830          REWRITE_TAC[diff2c_domain_alt; diff2_domain] THEN STRIP_TAC THEN
 831          new_rewrite [] [] diff2_imp_diff THEN REWRITE_TAC[] THEN
 832          FIRST_X_ASSUM MATCH_MP_TAC THEN
 833          MATCH_MP_TAC y_in_domain THEN EXISTS_TAC `w:real^N` THEN ASM_REWRITE_TAC[];
 834        ALL_TAC
 835      ] THEN
 836
 837      SUBGOAL_THEN `!x:real^N. x IN interval [domain] ==> &0 < f x` ASSUME_TAC THENL
 838      [
 839        GEN_TAC THEN DISCH_TAC THEN
 840          apply_tac interval_arith_pos THEN
 841          EXISTS_TAC `f_bounds:real#real` THEN
 842          ASM_REWRITE_TAC[] THEN
 843          FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[];
 844        ALL_TAC
 845      ] THEN
 846
 847      SUBGOAL_THEN `&0 < f (y:real^N)` ASSUME_TAC THENL
 848      [
 849        FIRST_X_ASSUM MATCH_MP_TAC THEN
 850          MATCH_MP_TAC y_in_domain THEN EXISTS_TAC `w:real^N` THEN
 851          ASM_REWRITE_TAC[];
 852        ALL_TAC
 853      ] THEN
 854
 855      REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[] THENL
 856      [
 857        UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN
 858          REWRITE_TAC[diff2c_domain] THEN REPEAT STRIP_TAC THEN
 859          ONCE_REWRITE_TAC[GSYM o_THM] THEN REWRITE_TAC[ETA_AX] THEN
 860          apply_tac diff2c_sqrt_compose THEN CONJ_TAC THEN FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[];
 861        new_rewrite [] [`sqrt _`] (GSYM o_THM) THEN REWRITE_TAC[ETA_AX] THEN
 862          apply_tac diff_uni_compose THEN ASM_REWRITE_TAC[] THEN
 863          MATCH_MP_TAC REAL_DIFFERENTIABLE_AT_SQRT THEN ASM_REWRITE_TAC[];
 864        new_rewrite [] [`sqrt _`] (GSYM o_THM) THEN REWRITE_TAC[ETA_AX] THEN
 865          MP_TAC (ISPECL [`y:real^N`; `f:real^N->real`] partial_uni_compose') THEN
 866          ASM_REWRITE_TAC[] THEN
 867          DISCH_THEN (MP_TAC o SPEC `sqrt`) THEN
 868          ANTS_TAC THENL
 869          [
 870            MATCH_MP_TAC REAL_DIFFERENTIABLE_AT_SQRT THEN ASM_REWRITE_TAC[];
 871            ALL_TAC
 872          ] THEN
 873          ASM_SIMP_TAC[derivative_sqrt];
 874        ALL_TAC
 875      ] THEN
 876
 877      UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN
 878      REWRITE_TAC[diff2c_domain_alt; diff2_domain] THEN
 879      DISCH_THEN (MP_TAC o SPEC `x:real^N` o CONJUNCT1) THEN ASM_REWRITE_TAC[] THEN DISCH_TAC THEN
 880      MP_TAC (ISPECL [`x:real^N`; `f:real^N->real`] second_partial_uni_compose') THEN
 881      ASM_REWRITE_TAC[] THEN
 882      DISCH_THEN (MP_TAC o SPEC `sqrt`) THEN
 883      ANTS_TAC THENL
 884      [
 885        MATCH_MP_TAC diff2_sqrt THEN
 886          FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[];
 887        ALL_TAC
 888      ] THEN
 889
 890      new_rewrite [] [`sqrt _`] (GSYM o_THM) THEN REWRITE_TAC[ETA_AX] THEN
 891      ASM_SIMP_TAC[second_derivative_sqrt; derivative_sqrt] THEN DISCH_THEN (fun th -> ALL_TAC) THEN
 892      REWRITE_TAC[REAL_ARITH `a pow 3 = a * a pow 2`] THEN
 893      new_rewrite [] [] SQRT_MUL THENL
 894      [
 895        REWRITE_TAC[REAL_LE_POW_2] THEN
 896          MATCH_MP_TAC REAL_LT_IMP_LE THEN FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[];
 897        ALL_TAC
 898      ] THEN
 899      new_rewrite [] [] POW_2_SQRT THENL
 900      [
 901        MATCH_MP_TAC REAL_LT_IMP_LE THEN FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[];
 902        ALL_TAC
 903      ] THEN
 904      ASM_SIMP_TAC[REAL_ARITH `&4 * a * b = (&2 * a) * (&2 * b)`]);;
 905
 906
 907
 908 (* atn *)
 909 let MK_M_TAYLOR_ATN' = (UNDISCH_ALL o PURE_REWRITE_RULE[float1_eq; float2_eq; num2_eq] o DISCH_ALL o
 910                           MY_RULE_FLOAT o prove)
 911   (`m_cell_domain domain (y:real^N) w ==>
 912      diff2c_domain domain f ==>
 913      interval_arith (atn (f y)) bounds ==>
 914      all_n 1 d_bounds_list (\i int. interval_arith (inv (&1 + f y * f y) * partial i f y) int) ==>
 915      (!x. x IN interval [domain] ==> all_n 1 dd_bounds_list (\i list_i. all_n 1 list_i
 916         (\j int. interval_arith ((((-- &2 * f x) * inv (&1 + f x * f x) pow 2) * partial j f x)
 917                                  * partial i f x + inv (&1 + f x * f x) * partial2 j i f x) int))) ==>
 918      m_taylor_interval (\x. atn (f x)) domain y w bounds d_bounds_list dd_bounds_list`,
 919    REWRITE_TAC[m_taylor_interval; m_lin_approx; second_bounded; ETA_AX] THEN
 920      REPEAT DISCH_TAC THEN
 921      SUBGOAL_THEN `lift o f differentiable at (y:real^N)` ASSUME_TAC THENL
 922      [
 923        UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN
 924          REWRITE_TAC[diff2c_domain_alt; diff2_domain] THEN STRIP_TAC THEN
 925          new_rewrite [] [] diff2_imp_diff THEN REWRITE_TAC[] THEN
 926          FIRST_X_ASSUM MATCH_MP_TAC THEN
 927          MATCH_MP_TAC y_in_domain THEN EXISTS_TAC `w:real^N` THEN ASM_REWRITE_TAC[];
 928        ALL_TAC
 929      ] THEN
 930
 931      REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[] THENL
 932      [
 933        UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN
 934          REWRITE_TAC[diff2c_domain] THEN REPEAT STRIP_TAC THEN
 935          ONCE_REWRITE_TAC[GSYM o_THM] THEN REWRITE_TAC[ETA_AX] THEN
 936          apply_tac diff2c_atn_compose THEN FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[];
 937        new_rewrite [] [`atn _`] (GSYM o_THM) THEN REWRITE_TAC[ETA_AX] THEN
 938          apply_tac diff_uni_compose THEN ASM_REWRITE_TAC[] THEN
 939          REWRITE_TAC[REAL_DIFFERENTIABLE_AT_ATN];
 940        new_rewrite [] [`atn _`] (GSYM o_THM) THEN REWRITE_TAC[ETA_AX] THEN
 941          MP_TAC (ISPECL [`y:real^N`; `f:real^N->real`] partial_uni_compose') THEN
 942          ASM_REWRITE_TAC[] THEN
 943          DISCH_THEN (MP_TAC o SPEC `atn`) THEN REWRITE_TAC[REAL_DIFFERENTIABLE_AT_ATN] THEN
 944          ASM_SIMP_TAC[derivative_atn];
 945        ALL_TAC
 946      ] THEN
 947
 948      UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN
 949      REWRITE_TAC[diff2c_domain_alt; diff2_domain] THEN
 950      DISCH_THEN (MP_TAC o SPEC `x:real^N` o CONJUNCT1) THEN ASM_REWRITE_TAC[] THEN DISCH_TAC THEN
 951      MP_TAC (ISPECL [`x:real^N`; `f:real^N->real`] second_partial_uni_compose') THEN
 952      ASM_REWRITE_TAC[] THEN
 953      DISCH_THEN (MP_TAC o SPEC `atn`) THEN
 954      REWRITE_TAC[diff2_atn] THEN
 955      new_rewrite [] [`atn _`] (GSYM o_THM) THEN REWRITE_TAC[ETA_AX] THEN
 956      ASM_SIMP_TAC[nth_derivative2; second_derivative_atn; derivative_atn] THEN
 957      new_rewrite [] [`f x pow 2`] REAL_POW_2 THEN
 958      ASM_SIMP_TAC[]);;
 959
 960
 961 (* acs *)
 962 let iabs_lemma = GEN_REWRITE_RULE (RAND_CONV o RAND_CONV) [GSYM float1_eq] (REFL `iabs f_bounds < &1`);;
 963
 964 let MK_M_TAYLOR_ACS' = (UNDISCH_ALL o PURE_ONCE_REWRITE_RULE[iabs_lemma] o
 965                           PURE_REWRITE_RULE[float1_eq; num3_eq] o DISCH_ALL o
 966                           MY_RULE_FLOAT o prove)
 967   (`m_cell_domain domain (y:real^N) w ==>
 968      (!x. x IN interval [domain] ==> interval_arith (f x) f_bounds) ==>
 969      iabs f_bounds < &1 ==>
 970      diff2c_domain domain f ==>
 971      interval_arith (acs (f y)) bounds ==>
 972      all_n 1 d_bounds_list
 973      (\i int. interval_arith (--inv (sqrt (&1 - f y * f y)) * partial i f y) int) ==>
 974      (!x. x IN interval [domain] ==> all_n 1 dd_bounds_list (\i list_i. all_n 1 list_i
 975         (\j int. interval_arith ((--(f x / sqrt ((&1 - f x * f x) pow 3)) * partial j f x) * partial i f x -
 976                                    inv (sqrt (&1 - f x * f x)) * partial2 j i f x) int))) ==>
 977      m_taylor_interval (\x. acs (f x)) domain y w bounds d_bounds_list dd_bounds_list`,
 978    REWRITE_TAC[m_taylor_interval; m_lin_approx; second_bounded; ETA_AX] THEN
 979      REPEAT DISCH_TAC THEN
 980      SUBGOAL_THEN `lift o f differentiable at (y:real^N)` ASSUME_TAC THENL
 981      [
 982        UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN
 983          REWRITE_TAC[diff2c_domain_alt; diff2_domain] THEN STRIP_TAC THEN
 984          new_rewrite [] [] diff2_imp_diff THEN REWRITE_TAC[] THEN
 985          FIRST_X_ASSUM MATCH_MP_TAC THEN
 986          MATCH_MP_TAC y_in_domain THEN EXISTS_TAC `w:real^N` THEN ASM_REWRITE_TAC[];
 987        ALL_TAC
 988      ] THEN
 989
 990      SUBGOAL_THEN `!x:real^N. x IN interval [domain] ==> abs (f x) < &1` ASSUME_TAC THENL
 991      [
 992        GEN_TAC THEN DISCH_TAC THEN
 993          apply_tac interval_arith_abs THEN
 994          EXISTS_TAC `f_bounds:real#real` THEN
 995          ASM_REWRITE_TAC[] THEN
 996          FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[];
 997        ALL_TAC
 998      ] THEN
 999
1000      SUBGOAL_THEN `abs (f (y:real^N)) < &1` ASSUME_TAC THENL
1001      [
1002        FIRST_X_ASSUM MATCH_MP_TAC THEN
1003          MATCH_MP_TAC y_in_domain THEN EXISTS_TAC `w:real^N` THEN
1004          ASM_REWRITE_TAC[];
1005        ALL_TAC
1006      ] THEN
1007
1008      REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[] THENL
1009      [
1010        UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN
1011          REWRITE_TAC[diff2c_domain] THEN REPEAT STRIP_TAC THEN
1012          ONCE_REWRITE_TAC[GSYM o_THM] THEN REWRITE_TAC[ETA_AX] THEN
1013          apply_tac diff2c_acs_compose THEN CONJ_TAC THEN FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[];
1014        new_rewrite [] [`acs _`] (GSYM o_THM) THEN REWRITE_TAC[ETA_AX] THEN
1015          apply_tac diff_uni_compose THEN ASM_REWRITE_TAC[] THEN
1016          MATCH_MP_TAC REAL_DIFFERENTIABLE_AT_ACS THEN ASM_REWRITE_TAC[];
1017        new_rewrite [] [`acs _`] (GSYM o_THM) THEN REWRITE_TAC[ETA_AX] THEN
1018          MP_TAC (ISPECL [`y:real^N`; `f:real^N->real`] partial_uni_compose') THEN
1019          ASM_REWRITE_TAC[] THEN
1020          DISCH_THEN (MP_TAC o SPEC `acs`) THEN
1021          ANTS_TAC THENL
1022          [
1023            MATCH_MP_TAC REAL_DIFFERENTIABLE_AT_ACS THEN ASM_REWRITE_TAC[];
1024            ALL_TAC
1025          ] THEN
1026          ASM_SIMP_TAC[derivative_acs];
1027        ALL_TAC
1028      ] THEN
1029
1030      UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN
1031      REWRITE_TAC[diff2c_domain_alt; diff2_domain] THEN
1032      DISCH_THEN (MP_TAC o SPEC `x:real^N` o CONJUNCT1) THEN ASM_REWRITE_TAC[] THEN DISCH_TAC THEN
1033      MP_TAC (ISPECL [`x:real^N`; `f:real^N->real`] second_partial_uni_compose') THEN
1034      ASM_REWRITE_TAC[] THEN
1035      DISCH_THEN (MP_TAC o SPEC `acs`) THEN
1036      ANTS_TAC THENL
1037      [
1038        MATCH_MP_TAC diff2_acs THEN
1039          FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[];
1040        ALL_TAC
1041      ] THEN
1042
1043      new_rewrite [] [`acs _`] (GSYM o_THM) THEN REWRITE_TAC[ETA_AX] THEN
1044      ASM_SIMP_TAC[second_derivative_acs; derivative_acs; REAL_MUL_LNEG; GSYM real_sub] THEN
1045      ASM_SIMP_TAC[GSYM REAL_MUL_LNEG]);;
1046
1047
1048
1049 (*************************)
1050
1051 (***************************************)
1052 (* eval_m_taylor_inv *)
1053
1054
1055 let inv_partial_lemma' =
1056   prove(`interval_arith (--inv (f y * f y) * partial i f (y:real^N)) int <=>
1057           (\i int. interval_arith (--inv (f y * f y) * partial i f y) int) i int`,
1058         REWRITE_TAC[]);;
1059
1060 let inv_second_lemma' =
1061   prove(`interval_arith (((&2 * inv (f x * f x * f x)) * partial j f x) * partial i f (x:real^N) -
1062                                    inv (f x * f x) * partial2 j i f x) int <=>
1063           (\j int. interval_arith (((&2 * inv (f x * f x * f x)) * partial j f x) * partial i f x -
1064                                      inv (f x * f x) * partial2 j i f x) int) j int`,
1065         REWRITE_TAC[]);;
1066
1067
1068 let inv_second_lemma'' = (PURE_REWRITE_RULE[GSYM num1_eq] o prove)
1069   (`all_n 1 list
1070      (\j int. interval_arith (((&2 * inv (f x * f x * f x)) * partial j f x) * partial i f (x:real^N) -
1071                                 inv (f x * f x) * partial2 j i f x) int) <=>
1072      (\i list. all_n 1 list
1073         (\j int. interval_arith (((&2 * inv (f x * f x * f x)) * partial j f x) * partial i f x -
1074                                             inv (f x * f x) * partial2 j i f x) int)) i list`,
1075    REWRITE_TAC[]);;
1076
1077
1078
1079
1080 let eval_m_taylor_inv n p_lin p_second taylor1_th =
1081   let domain_th, diff2_f1_th, lin1_th, second1_th = dest_m_taylor_thms n taylor1_th in
1082   let f1_tm = (rand o concl) diff2_f1_th in
1083   let domain_tm, y_tm, w_tm = dest_m_cell_domain (concl domain_th) in
1084   let ty = type_of y_tm in
1085
1086   let x_var = mk_var ("x", ty) and
1087       y_var = mk_var ("y", ty) and
1088       w_var = mk_var ("w", ty) and
1089       f_var = mk_var ("f", type_of f1_tm) and
1090       domain_var = mk_var ("domain", type_of domain_tm) in
1091
1092   let undisch = UNDISCH o SPEC x_var in
1093
1094   let f1_bound0 = eval_m_taylor_bound n p_second taylor1_th in
1095   let f1_bound = undisch f1_bound0 in
1096   let f_bounds_tm = (rand o concl) f1_bound in
1097
1098   (* cond *)
1099   let cond_th = check_interval_not_zero f_bounds_tm in
1100
1101   let _, bounds1_th, df1_th = m_lin_approx_components n lin1_th in
1102
1103   let bounds_th = float_interval_inv p_lin bounds1_th in
1104   let bounds_tm = (rand o concl) bounds_th in
1105
1106   (* partial_lemma' *)
1107   let u_lemma0 = (INST[f1_tm, f_var; y_tm, y_var] o
1108                     INST_TYPE[n_type_array.(n), nty]) inv_partial_lemma' in
1109
1110   let u_bounds =
1111     let neg, inv, ( * ) = float_interval_neg, float_interval_inv p_lin, float_interval_mul p_lin in
1112       neg (inv (bounds1_th * bounds1_th)) in
1113
1114
1115   let u_lin th1 =
1116     (* partial *)
1117     let u_th =
1118       let ( * ) = float_interval_mul p_lin in
1119         u_bounds * th1 in
1120     let int_tm = rand (concl u_th) in
1121     let i_tm = (rand o rator o rator o lhand) (concl th1) in
1122     let th0 = INST[i_tm, i_var_num; int_tm, int_var] u_lemma0 in
1123       EQ_MP th0 u_th in
1124
1125   let df_th = eval_all_n df1_th true u_lin in
1126   let d_bounds_list = (rand o rator o concl) df_th in
1127
1128
1129   let dd1 = second_bounded_components n second1_th in
1130
1131   (* second_lemma', second_lemma'' *)
1132   let u_second_lemma0 = (INST[f1_tm, f_var] o
1133                            INST_TYPE[n_type_array.(n), nty]) inv_second_lemma' in
1134
1135   let u_second_lemma1 = (INST[f1_tm, f_var] o
1136                            INST_TYPE[n_type_array.(n), nty]) inv_second_lemma'' in
1137
1138
1139   let d1_bounds = map (fun i ->
1140                          let th0 = eval_m_taylor_partial_bound n p_second i taylor1_th in
1141                            undisch th0) (1--n) in
1142
1143   (* u'(f x), u''(f x) *)
1144   let d1_th0, d2_th0 =
1145     let inv, ( * ) = float_interval_inv p_second, float_interval_mul p_second in
1146     let ff = f1_bound * f1_bound in
1147       inv ff,
1148     two_interval * inv (f1_bound * ff) in
1149
1150
1151   let u_second2 th1 =
1152     let i_tm = (rand o rator o concl) th1 in
1153     let i_int = (Num.int_of_num o raw_dest_hash) i_tm in
1154     let di1 = List.nth d1_bounds (i_int - 1) in
1155     let th1 = MY_BETA_RULE th1 in
1156     let lemma = INST[i_tm, i_var_num] u_second_lemma0 in
1157     let u_second th1 =
1158       let j_tm = (rand o rator o rator o rator o lhand) (concl th1) in
1159       let j_int = (Num.int_of_num o raw_dest_hash) j_tm in
1160       let dj1 = List.nth d1_bounds (j_int - 1) in
1161
1162         (* partial2 *)
1163       let u_th =
1164         let ( * ), ( - ) = float_interval_mul p_second, float_interval_sub p_second in
1165           (d2_th0 * dj1) * di1 - d1_th0 * th1 in
1166
1167       let int_tm = rand (concl u_th) in
1168       let th0 = INST[j_tm, j_var_num; int_tm, int_var] lemma in
1169         EQ_MP th0 u_th in
1170
1171     let u_th = eval_all_n th1 true u_second in
1172     let list_tm = (rand o rator o concl) u_th in
1173     let lemma1 = INST[i_tm, i_var_num; list_tm, list_var_real_pair] u_second_lemma1 in
1174       EQ_MP lemma1 u_th in
1175
1176   let dd_th0 = eval_all_n dd1 false u_second2 in
1177   let dd_list = (rand o rator o concl) dd_th0 in
1178   let dd_th = GEN x_var (DISCH_ALL dd_th0) in
1179
1180   let th = (MY_PROVE_HYP dd_th o MY_PROVE_HYP diff2_f1_th o
1181               MY_PROVE_HYP cond_th o MY_PROVE_HYP f1_bound0 o
1182               MY_PROVE_HYP bounds_th o MY_PROVE_HYP df_th o MY_PROVE_HYP domain_th o
1183               INST[f1_tm, f_var; f_bounds_tm, f_bounds_var;
1184                    domain_tm, domain_var; y_tm, y_var; w_tm, w_var;
1185                    bounds_tm, bounds_var; d_bounds_list, d_bounds_list_var;
1186                    dd_list, dd_bounds_list_var] o
1187               INST_TYPE[n_type_array.(n), nty]) MK_M_TAYLOR_INV' in
1188   let eq_th = unary_beta_gen_eq f1_tm x_var inv_op_real in
1189     m_taylor_interval_norm th eq_th;;
1190
1191
1192 (***************************************)
1193 (* eval_m_taylor_sqrt *)
1194
1195 let sqrt_partial_lemma' =
1196   prove(`interval_arith (inv (&2 * sqrt (f y)) * partial i f (y:real^N)) int <=>
1197           (\i int. interval_arith (inv (&2 * sqrt (f y)) * partial i f y) int) i int`,
1198         REWRITE_TAC[]);;
1199
1200
1201 let sqrt_second_lemma' =
1202   prove(`interval_arith ((--inv ((&2 * sqrt (f x)) * (&2 * f x)) * partial j f x) * partial i f (x:real^N) +
1203                            inv (&2 * sqrt (f x)) * partial2 j i f x) int <=>
1204           (\j int. interval_arith ((--inv ((&2 * sqrt (f x))*(&2 * f x)) * partial j f x) * partial i f x +
1205                                      inv (&2 * sqrt (f x)) * partial2 j i f x) int) j int`,
1206         REWRITE_TAC[]);;
1207
1208
1209 let sqrt_second_lemma'' = (PURE_REWRITE_RULE[GSYM num1_eq] o prove)
1210   (`all_n 1 list
1211      (\j int. interval_arith ((--inv ((&2 * sqrt (f x)) * (&2 * f x)) * partial j f x) * partial i f x +
1212                                      inv (&2 * sqrt (f x)) * partial2 j i f (x:real^N)) int) <=>
1213      (\i list. all_n 1 list
1214         (\j int. interval_arith ((--inv ((&2 * sqrt (f x)) * (&2 * f x)) * partial j f x) * partial i f x +
1215                                    inv (&2 * sqrt (f x)) * partial2 j i f (x:real^N)) int)) i list`,
1216    REWRITE_TAC[]);;
1217
1218
1219 let eval_m_taylor_sqrt n p_lin p_second taylor1_th =
1220   let domain_th, diff2_f1_th, lin1_th, second1_th = dest_m_taylor_thms n taylor1_th in
1221   let f1_tm = (rand o concl) diff2_f1_th in
1222   let domain_tm, y_tm, w_tm = dest_m_cell_domain (concl domain_th) in
1223   let ty = type_of y_tm in
1224
1225   let x_var = mk_var ("x", ty) and
1226       y_var = mk_var ("y", ty) and
1227       w_var = mk_var ("w", ty) and
1228       f_var = mk_var ("f", type_of f1_tm) and
1229       domain_var = mk_var ("domain", type_of domain_tm) in
1230
1231   let undisch = UNDISCH o SPEC x_var in
1232
1233   let f1_bound0 = eval_m_taylor_bound n p_second taylor1_th in
1234   let f1_bound = undisch f1_bound0 in
1235   let f_bounds_tm = (rand o concl) f1_bound in
1236
1237   (* cond *)
1238   let cond_th = check_interval_pos f_bounds_tm in
1239
1240   let _, bounds1_th, df1_th = m_lin_approx_components n lin1_th in
1241
1242   let bounds_th = float_interval_sqrt p_lin bounds1_th in
1243   let bounds_tm = (rand o concl) bounds_th in
1244
1245   (* partial_lemma' *)
1246   let u_lemma0 = (INST[f1_tm, f_var; y_tm, y_var] o
1247                     INST_TYPE[n_type_array.(n), nty]) sqrt_partial_lemma' in
1248
1249   let u_bounds =
1250     let inv, ( * ) = float_interval_inv p_lin, float_interval_mul p_lin in
1251       inv (two_interval * bounds_th) in
1252
1253   let u_lin th1 =
1254     (* partial *)
1255     let u_th =
1256       let ( * ) = float_interval_mul p_lin in
1257         u_bounds * th1 in
1258     let int_tm = rand (concl u_th) in
1259     let i_tm = (rand o rator o rator o lhand) (concl th1) in
1260     let th0 = INST[i_tm, i_var_num; int_tm, int_var] u_lemma0 in
1261       EQ_MP th0 u_th in
1262
1263   let df_th = eval_all_n df1_th true u_lin in
1264   let d_bounds_list = (rand o rator o concl) df_th in
1265
1266
1267   let dd1 = second_bounded_components n second1_th in
1268
1269   (* second_lemma', second_lemma'' *)
1270   let u_second_lemma0 = (INST[f1_tm, f_var] o
1271                            INST_TYPE[n_type_array.(n), nty]) sqrt_second_lemma' in
1272
1273   let u_second_lemma1 = (INST[f1_tm, f_var] o
1274                            INST_TYPE[n_type_array.(n), nty]) sqrt_second_lemma'' in
1275
1276   let d1_bounds = map (fun i ->
1277                          let th0 = eval_m_taylor_partial_bound n p_second i taylor1_th in
1278                            undisch th0) (1--n) in
1279
1280   (* u'(f x), u''(f x) *)
1281   let d1_th0, d2_th0 =
1282     let neg, sqrt, inv, ( * ) = float_interval_neg, float_interval_sqrt p_second,
1283       float_interval_inv p_second, float_interval_mul p_second in
1284     let two_sqrt_f = two_interval * sqrt f1_bound in
1285       inv two_sqrt_f,
1286       neg (inv (two_sqrt_f * (two_interval * f1_bound))) in
1287
1288   let u_second2 th1 =
1289     let i_tm = (rand o rator o concl) th1 in
1290     let i_int = (Num.int_of_num o raw_dest_hash) i_tm in
1291     let di1 = List.nth d1_bounds (i_int - 1) in
1292     let th1 = MY_BETA_RULE th1 in
1293     let lemma = INST[i_tm, i_var_num] u_second_lemma0 in
1294     let u_second th1 =
1295       let j_tm = (rand o rator o rator o rator o lhand) (concl th1) in
1296       let j_int = (Num.int_of_num o raw_dest_hash) j_tm in
1297       let dj1 = List.nth d1_bounds (j_int - 1) in
1298
1299         (* partial2 *)
1300       let u_th =
1301         let ( * ), ( + ) = float_interval_mul p_second, float_interval_add p_second in
1302           (d2_th0 * dj1) * di1 + d1_th0 * th1 in
1303
1304       let int_tm = rand (concl u_th) in
1305       let th0 = INST[j_tm, j_var_num; int_tm, int_var] lemma in
1306         EQ_MP th0 u_th in
1307
1308     let u_th = eval_all_n th1 true u_second in
1309     let list_tm = (rand o rator o concl) u_th in
1310     let lemma1 = INST[i_tm, i_var_num; list_tm, list_var_real_pair] u_second_lemma1 in
1311       EQ_MP lemma1 u_th in
1312
1313   let dd_th0 = eval_all_n dd1 false u_second2 in
1314   let dd_list = (rand o rator o concl) dd_th0 in
1315   let dd_th = GEN x_var (DISCH_ALL dd_th0) in
1316
1317   let th = (MY_PROVE_HYP dd_th o MY_PROVE_HYP diff2_f1_th o
1318               MY_PROVE_HYP cond_th o MY_PROVE_HYP f1_bound0 o
1319               MY_PROVE_HYP bounds_th o MY_PROVE_HYP df_th o MY_PROVE_HYP domain_th o
1320               INST[f1_tm, f_var; f_bounds_tm, f_bounds_var;
1321                    domain_tm, domain_var; y_tm, y_var; w_tm, w_var;
1322                    bounds_tm, bounds_var; d_bounds_list, d_bounds_list_var;
1323                    dd_list, dd_bounds_list_var] o
1324               INST_TYPE[n_type_array.(n), nty]) MK_M_TAYLOR_SQRT' in
1325   let eq_th = unary_beta_gen_eq f1_tm x_var sqrt_tm in
1326     m_taylor_interval_norm th eq_th;;
1327
1328
1329
1330 (***************************************)
1331 (* eval_m_taylor_atn *)
1332
1333 let atn_partial_lemma' =
1334   prove(`interval_arith (inv (&1 + f y * f y) * partial i f (y:real^N)) int <=>
1335           (\i int. interval_arith (inv (&1 + f y * f y) * partial i f y) int) i int`,
1336         REWRITE_TAC[]);;
1337
1338
1339 let atn_second_lemma' =
1340   prove(`interval_arith ((((-- &2 * f x) * inv (&1 + f x * f x) pow 2) * partial j f (x:real^N))
1341                          * partial i f x + inv (&1 + f x * f x) * partial2 j i f x) int <=>
1342           (\j int. interval_arith ((((-- &2 * f x) * inv (&1 + f x * f x) pow 2) * partial j f x)
1343                                    * partial i f x + inv (&1 + f x * f x) * partial2 j i f x) int) j int`,
1344         REWRITE_TAC[]);;
1345
1346
1347 let atn_second_lemma'' = (PURE_REWRITE_RULE[float1_eq; float2_eq; num2_eq] o NUMERALS_TO_NUM o
1348                             PURE_REWRITE_RULE[FLOAT_OF_NUM; min_exp_def] o prove)
1349   (`all_n 1 list
1350      (\j int. interval_arith ((((-- &2 * f x) * inv (&1 + f x * f x) pow 2) * partial j f (x:real^N))
1351                               * partial i f x + inv (&1 + f x * f x) * partial2 j i f x) int) <=>
1352      (\i list. all_n 1 list
1353         (\j int. interval_arith ((((-- &2 * f x) * inv (&1 + f x * f x) pow 2) * partial j f x)
1354                                  * partial i f x + inv (&1 + f x * f x) * partial2 j i f x) int)) i list`,
1355    REWRITE_TAC[]);;
1356
1357
1358
1359 let eval_m_taylor_atn n p_lin p_second taylor1_th =
1360   let domain_th, diff2_f1_th, lin1_th, second1_th = dest_m_taylor_thms n taylor1_th in
1361   let f1_tm = (rand o concl) diff2_f1_th in
1362   let domain_tm, y_tm, w_tm = dest_m_cell_domain (concl domain_th) in
1363   let ty = type_of y_tm in
1364
1365   let x_var = mk_var ("x", ty) and
1366       y_var = mk_var ("y", ty) and
1367       w_var = mk_var ("w", ty) and
1368       f_var = mk_var ("f", type_of f1_tm) and
1369       domain_var = mk_var ("domain", type_of domain_tm) in
1370
1371   let undisch = UNDISCH o SPEC x_var in
1372
1373   let f1_bound0 = eval_m_taylor_bound n p_second taylor1_th in
1374   let f1_bound = undisch f1_bound0 in
1375
1376   let _, bounds1_th, df1_th = m_lin_approx_components n lin1_th in
1377
1378   let bounds_th = float_interval_atn p_lin bounds1_th in
1379   let bounds_tm = (rand o concl) bounds_th in
1380
1381   (* partial_lemma' *)
1382   let u_lemma0 = (INST[f1_tm, f_var; y_tm, y_var] o
1383                     INST_TYPE[n_type_array.(n), nty]) atn_partial_lemma' in
1384
1385   let u_bounds =
1386     let inv, ( + ), ( * ) = float_interval_inv p_lin, float_interval_add p_lin, float_interval_mul p_lin in
1387       inv (one_interval + bounds1_th * bounds1_th) in
1388
1389   let u_lin th1 =
1390     (* partial *)
1391     let u_th =
1392       let ( * ) = float_interval_mul p_lin in
1393         u_bounds * th1 in
1394     let int_tm = rand (concl u_th) in
1395     let i_tm = (rand o rator o rator o lhand) (concl th1) in
1396     let th0 = INST[i_tm, i_var_num; int_tm, int_var] u_lemma0 in
1397       EQ_MP th0 u_th in
1398
1399   let df_th = eval_all_n df1_th true u_lin in
1400   let d_bounds_list = (rand o rator o concl) df_th in
1401
1402
1403   let dd1 = second_bounded_components n second1_th in
1404
1405   (* second_lemma', second_lemma'' *)
1406   let u_second_lemma0 = (INST[f1_tm, f_var] o
1407                            INST_TYPE[n_type_array.(n), nty]) atn_second_lemma' in
1408
1409   let u_second_lemma1 = (INST[f1_tm, f_var] o
1410                            INST_TYPE[n_type_array.(n), nty]) atn_second_lemma'' in
1411
1412   let d1_bounds = map (fun i ->
1413                          let th0 = eval_m_taylor_partial_bound n p_second i taylor1_th in
1414                            undisch th0) (1--n) in
1415
1416
1417   (* u'(f x), u''(f x) *)
1418   let d1_th0, d2_th0 =
1419     let neg, inv, ( + ), ( * ), pow2 = float_interval_neg, float_interval_inv p_second,
1420       float_interval_add p_second, float_interval_mul p_second, float_interval_pow_simple p_second 2 in
1421     let inv_one_ff = inv (one_interval + f1_bound * f1_bound) in
1422       inv_one_ff,
1423     (neg_two_interval * f1_bound) * pow2 inv_one_ff in
1424
1425   let u_second2 th1 =
1426     let i_tm = (rand o rator o concl) th1 in
1427     let i_int = (Num.int_of_num o raw_dest_hash) i_tm in
1428     let di1 = List.nth d1_bounds (i_int - 1) in
1429     let th1 = MY_BETA_RULE th1 in
1430     let lemma = INST[i_tm, i_var_num] u_second_lemma0 in
1431     let u_second th1 =
1432       let j_tm = (rand o rator o rator o rator o lhand) (concl th1) in
1433       let j_int = (Num.int_of_num o raw_dest_hash) j_tm in
1434       let dj1 = List.nth d1_bounds (j_int - 1) in
1435
1436       (* partial2 *)
1437       let u_th =
1438         let ( * ), ( + ) = float_interval_mul p_second, float_interval_add p_second in
1439           (d2_th0 * dj1) * di1 + d1_th0 * th1 in
1440
1441       let int_tm = rand (concl u_th) in
1442       let th0 = INST[j_tm, j_var_num; int_tm, int_var] lemma in
1443         EQ_MP th0 u_th in
1444
1445     let u_th = eval_all_n th1 true u_second in
1446     let list_tm = (rand o rator o concl) u_th in
1447     let lemma1 = INST[i_tm, i_var_num; list_tm, list_var_real_pair] u_second_lemma1 in
1448       EQ_MP lemma1 u_th in
1449
1450   let dd_th0 = eval_all_n dd1 false u_second2 in
1451   let dd_list = (rand o rator o concl) dd_th0 in
1452   let dd_th = GEN x_var (DISCH_ALL dd_th0) in
1453
1454   let th = (MY_PROVE_HYP dd_th o MY_PROVE_HYP diff2_f1_th o
1455               MY_PROVE_HYP bounds_th o MY_PROVE_HYP df_th o MY_PROVE_HYP domain_th o
1456               INST[f1_tm, f_var;
1457                    domain_tm, domain_var; y_tm, y_var; w_tm, w_var;
1458                    bounds_tm, bounds_var; d_bounds_list, d_bounds_list_var;
1459                    dd_list, dd_bounds_list_var] o
1460               INST_TYPE[n_type_array.(n), nty]) MK_M_TAYLOR_ATN' in
1461   let eq_th = unary_beta_gen_eq f1_tm x_var atn_tm in
1462     m_taylor_interval_norm th eq_th;;
1463
1464
1465
1466 (***************************************)
1467 (* eval_m_taylor_acs *)
1468
1469 let acs_partial_lemma' =
1470   prove(`interval_arith (--inv (sqrt (&1 - f y * f y)) * partial i f (y:real^N)) int <=>
1471           (\i int. interval_arith (--inv (sqrt (&1 - f y * f y)) * partial i f y) int) i int`,
1472         REWRITE_TAC[]);;
1473
1474
1475 let acs_second_lemma' =
1476   prove(`interval_arith ((--(f x / sqrt ((&1 - f x * f x) pow 3)) * partial j f x) * partial i f (x:real^N) - inv (sqrt (&1 - f x * f x)) * partial2 j i f x) int <=>
1477
1478           (\j int. interval_arith ((--(f x / sqrt ((&1 - f x * f x) pow 3)) * partial j f x) * partial i f (x:real^N) - inv (sqrt (&1 - f x * f x)) * partial2 j i f x) int) j int`,
1479         REWRITE_TAC[]);;
1480
1481
1482 let acs_second_lemma'' = (PURE_REWRITE_RULE[float1_eq; float2_eq; num3_eq; num2_eq] o NUMERALS_TO_NUM o
1483                             PURE_REWRITE_RULE[FLOAT_OF_NUM; min_exp_def] o prove)
1484   (`all_n 1 list
1485      (\j int. interval_arith ((--(f x / sqrt ((&1 - f x * f x) pow 3)) * partial j f x) * partial i f (x:real^N) - inv (sqrt (&1 - f x * f x)) * partial2 j i f x) int) <=>
1486      (\i list. all_n 1 list
1487         (\j int. interval_arith ((--(f x / sqrt ((&1 - f x * f x) pow 3)) * partial j f x) * partial i f (x:real^N) - inv (sqrt (&1 - f x * f x)) * partial2 j i f x) int)) i list`,
1488    REWRITE_TAC[]);;
1489
1490
1491 let eval_m_taylor_acs n p_lin p_second taylor1_th =
1492   let domain_th, diff2_f1_th, lin1_th, second1_th = dest_m_taylor_thms n taylor1_th in
1493   let f1_tm = (rand o concl) diff2_f1_th in
1494   let domain_tm, y_tm, w_tm = dest_m_cell_domain (concl domain_th) in
1495   let ty = type_of y_tm in
1496
1497   let x_var = mk_var ("x", ty) and
1498       y_var = mk_var ("y", ty) and
1499       w_var = mk_var ("w", ty) and
1500       f_var = mk_var ("f", type_of f1_tm) and
1501       domain_var = mk_var ("domain", type_of domain_tm) in
1502
1503   let undisch = UNDISCH o SPEC x_var in
1504
1505   let f1_bound0 = eval_m_taylor_bound n p_second taylor1_th in
1506   let f1_bound = undisch f1_bound0 in
1507   let f_bounds_tm = (rand o concl) f1_bound in
1508
1509   (* cond *)
1510   let cond_th = EQT_ELIM (check_interval_iabs f_bounds_tm one_float) in
1511
1512   let _, bounds1_th, df1_th = m_lin_approx_components n lin1_th in
1513
1514   let bounds_th = float_interval_acs p_lin bounds1_th in
1515   let bounds_tm = (rand o concl) bounds_th in
1516
1517   (* partial_lemma' *)
1518   let u_lemma0 = (INST[f1_tm, f_var; y_tm, y_var] o
1519                     INST_TYPE[n_type_array.(n), nty]) acs_partial_lemma' in
1520
1521   let u_bounds =
1522     let inv, sqrt, neg = float_interval_inv p_lin, float_interval_sqrt p_lin, float_interval_neg in
1523     let ( * ), (-) = float_interval_mul p_lin, float_interval_sub p_lin in
1524       neg (inv (sqrt (one_interval - bounds1_th * bounds1_th))) in
1525
1526   let u_lin th1 =
1527     (* partial *)
1528     let u_th =
1529       let ( * ) = float_interval_mul p_lin in
1530         u_bounds * th1 in
1531     let int_tm = rand (concl u_th) in
1532     let i_tm = (rand o rator o rator o lhand) (concl th1) in
1533     let th0 = INST[i_tm, i_var_num; int_tm, int_var] u_lemma0 in
1534       EQ_MP th0 u_th in
1535
1536   let df_th = eval_all_n df1_th true u_lin in
1537   let d_bounds_list = (rand o rator o concl) df_th in
1538
1539
1540   let dd1 = second_bounded_components n second1_th in
1541
1542   (* second_lemma', second_lemma'' *)
1543   let u_second_lemma0 = (INST[f1_tm, f_var] o
1544                            INST_TYPE[n_type_array.(n), nty]) acs_second_lemma' in
1545
1546   let u_second_lemma1 = (INST[f1_tm, f_var] o
1547                            INST_TYPE[n_type_array.(n), nty]) acs_second_lemma'' in
1548
1549   let d1_bounds = map (fun i ->
1550                          let th0 = eval_m_taylor_partial_bound n p_second i taylor1_th in
1551                            undisch th0) (1--n) in
1552
1553   (* u'(f x), u''(f x) *)
1554   let d1_th0, d2_th0 =
1555     let neg, sqrt, inv = float_interval_neg, float_interval_sqrt p_second, float_interval_inv p_second in
1556     let (-), ( * ), (/) = float_interval_sub p_second, float_interval_mul p_second, float_interval_div p_second in
1557     let pow3 = float_interval_pow_simple p_second 3 in
1558     let ff_1 = one_interval - f1_bound * f1_bound in
1559       inv (sqrt ff_1),
1560     neg (f1_bound / sqrt (pow3 ff_1)) in
1561
1562   let u_second2 th1 =
1563     let i_tm = (rand o rator o concl) th1 in
1564     let i_int = (Num.int_of_num o raw_dest_hash) i_tm in
1565     let di1 = List.nth d1_bounds (i_int - 1) in
1566     let th1 = MY_BETA_RULE th1 in
1567     let lemma = INST[i_tm, i_var_num] u_second_lemma0 in
1568     let u_second th1 =
1569       let j_tm = (rand o rator o rator o rator o lhand) (concl th1) in
1570       let j_int = (Num.int_of_num o raw_dest_hash) j_tm in
1571       let dj1 = List.nth d1_bounds (j_int - 1) in
1572
1573         (* partial2 *)
1574       let u_th =
1575         let ( * ), ( - ) = float_interval_mul p_second, float_interval_sub p_second in
1576           (d2_th0 * dj1) * di1 - d1_th0 * th1 in
1577
1578       let int_tm = rand (concl u_th) in
1579       let th0 = INST[j_tm, j_var_num; int_tm, int_var] lemma in
1580         EQ_MP th0 u_th in
1581
1582     let u_th = eval_all_n th1 true u_second in
1583     let list_tm = (rand o rator o concl) u_th in
1584     let lemma1 = INST[i_tm, i_var_num; list_tm, list_var_real_pair] u_second_lemma1 in
1585       EQ_MP lemma1 u_th in
1586
1587   let dd_th0 = eval_all_n dd1 false u_second2 in
1588   let dd_list = (rand o rator o concl) dd_th0 in
1589   let dd_th = GEN x_var (DISCH_ALL dd_th0) in
1590
1591   let th = (MY_PROVE_HYP dd_th o MY_PROVE_HYP diff2_f1_th o
1592               MY_PROVE_HYP cond_th o MY_PROVE_HYP f1_bound0 o
1593               MY_PROVE_HYP bounds_th o MY_PROVE_HYP df_th o MY_PROVE_HYP domain_th o
1594               INST[f1_tm, f_var; f_bounds_tm, f_bounds_var;
1595                    domain_tm, domain_var; y_tm, y_var; w_tm, w_var;
1596                    bounds_tm, bounds_var; d_bounds_list, d_bounds_list_var;
1597                    dd_list, dd_bounds_list_var] o
1598               INST_TYPE[n_type_array.(n), nty]) MK_M_TAYLOR_ACS' in
1599   let eq_th = unary_beta_gen_eq f1_tm x_var acs_tm in
1600     m_taylor_interval_norm th eq_th;;
1601
1602
1603
1604 end;;