formal_lp/old/formal_interval/m_taylor_arith.hl

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