formal_lp/old/formal_interval/m_taylor.hl

   1 (* TODO: move lemmas about TABLE into a separate file *)\r
   2 (* TODO: remove dependencies on Packing3, Arc_properties from theory files as well *)\r
   3 needs "packing/pack3.hl";;\r
   4 \r
   5 needs "../formal_lp/formal_interval/more_float.hl";;\r
   6 needs "../formal_lp/formal_interval/theory/taylor_interval.hl";;\r
   7 needs "../formal_lp/formal_interval/theory/multivariate_taylor.hl";;\r
   8 needs "../formal_lp/formal_interval/second_approx.hl";;\r
   9 needs "../formal_lp/formal_interval/eval_interval.hl";;\r
  10 needs "../formal_lp/list/list_conversions.hl";;\r
  11 needs "../formal_lp/list/list_float.hl";;\r
  12 \r
  13 let max_dim = 8;;\r
  14 \r
  15 let inst_first_type_var ty th =\r
  16   let ty_vars = type_vars_in_term (concl th) in\r
  17     if ty_vars = [] then\r
  18       failwith "inst_first_type: no type variables in the theorem"\r
  19     else\r
  20       INST_TYPE [ty, hd ty_vars] th;;\r
  21 \r
  22 let float0 = mk_float 0 Arith_options.min_exp and\r
  23     interval0 = mk_float_interval_small_num 0;;\r
  24 \r
  25 \r
  26 let real_ty = `:real` and\r
  27     real_list_ty = `:(real)list` and\r
  28     real_pair_ty = `:real#real` and\r
  29     real_pair_list_ty = `:(real#real)list` and\r
  30     nty = `:N`;;\r
  31 \r
  32 let d_bounds_list_var = `d_bounds_list : (real#real)list` and\r
  33     dd_bounds_list_var = `dd_bounds_list : ((real#real)list)list` and\r
  34     error_var = `error : real` and\r
  35     dd_list_var = `dd_list : (real#real)list`;;\r
  36 \r
  37 let has_size_array = Array.init (max_dim + 1) \r
  38   (fun i -> match i with\r
  39      | 0 -> TRUTH\r
  40      | 1 -> HAS_SIZE_1\r
  41      | _ -> define_finite_type i);;\r
  42 \r
  43 let dimindex_array = Array.init (max_dim + 1) \r
  44   (fun i -> if i < 1 then TRUTH else MATCH_MP DIMINDEX_UNIQUE has_size_array.(i));;\r
  45 \r
  46 let n_type_array = Array.init (max_dim + 1)\r
  47   (fun i -> if i < 1 then bool_ty else \r
  48      let dimindex_th = dimindex_array.(i) in\r
  49        (hd o snd o dest_type o snd o dest_const o rand o lhand o concl) dimindex_th);;\r
  50 \r
  51 let n_vector_type_array = Array.init (max_dim + 1)\r
  52   (fun i -> if i < 1 then bool_ty else mk_type ("cart", [real_ty; n_type_array.(i)]));;\r
  53 \r
  54 \r
  55 \r
  56 (************************************)\r
  57 (* m_cell_domain *)\r
  58 \r
  59 let ALL2_ALL_ZIP = prove(`!(P:A->B->bool) l1 l2. LENGTH l1 = LENGTH l2 ==> \r
  60     (ALL2 P l1 l2 <=> ALL (\p. P (FST p) (SND p)) (ZIP l1 l2))`,\r
  61   GEN_TAC THEN LIST_INDUCT_TAC THENL\r
  62     [\r
  63       GEN_TAC THEN REWRITE_TAC[LENGTH; EQ_SYM_EQ; LENGTH_EQ_NIL] THEN \r
  64         DISCH_THEN (fun th -> REWRITE_TAC[th]) THEN\r
  65         REWRITE_TAC[ZIP; ALL2; ALL];\r
  66       ALL_TAC\r
  67     ] THEN\r
  68 \r
  69     LIST_INDUCT_TAC THEN REWRITE_TAC[LENGTH] THENL [ARITH_TAC; ALL_TAC] THEN\r
  70     REWRITE_TAC[eqSS] THEN DISCH_TAC THEN\r
  71     REWRITE_TAC[ALL2; ZIP; ALL] THEN\r
  72     FIRST_X_ASSUM (new_rewrite [] []) THEN ASM_REWRITE_TAC[]);;\r
  73 \r
  74 let EL_ZIP = prove(`!(l1:(A)list) (l2:(B)list) i. LENGTH l1 = LENGTH l2 /\ i < LENGTH l1 ==> \r
  75     EL i (ZIP l1 l2) = (EL i l1, EL i l2)`,\r
  76   LIST_INDUCT_TAC THEN LIST_INDUCT_TAC THEN REWRITE_TAC[ZIP; LENGTH] THEN TRY ARITH_TAC THEN\r
  77     case THEN REWRITE_TAC[EL; HD; TL] THEN GEN_TAC THEN\r
  78     REWRITE_TAC[eqSS; ARITH_RULE `SUC n < SUC x <=> n < x`] THEN STRIP_TAC THEN\r
  79     FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[]);;\r
  80 \r
  81 let LENGTH_ZIP = prove(`!l1 l2. LENGTH l1 = LENGTH l2 ==> LENGTH (ZIP l1 l2) = LENGTH l1`,\r
  82   LIST_INDUCT_TAC THEN LIST_INDUCT_TAC THEN REWRITE_TAC[ZIP; LENGTH] THEN TRY ARITH_TAC THEN\r
  83     REWRITE_TAC[eqSS] THEN DISCH_TAC THEN FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[]);;\r
  84 \r
  85     \r
  86 \r
  87 let VECTOR_COMPONENT = prove(`!l i. i IN 1..dimindex (:N) ==>\r
  88                                (vector l:A^N)$i = EL (i - 1) l`,\r
  89 REWRITE_TAC[IN_NUMSEG] THEN REPEAT GEN_TAC THEN DISCH_TAC THEN REWRITE_TAC[vector] THEN\r
  90   MATCH_MP_TAC LAMBDA_BETA THEN ASM_REWRITE_TAC[]);;\r
  91 \r
  92 \r
  93 let test_domain_xi = new_definition\r
  94   `test_domain_xi xz yw <=> FST xz <= FST yw /\ FST yw <= SND xz /\ \r
  95   FST yw - FST xz <= SND yw /\ SND xz - FST yw <= SND yw`;;\r
  96 \r
  97 let MK_CELL_DOMAIN = prove(`!xz (yw:(real#real)list) x z y w.\r
  98     LENGTH x = dimindex (:N) /\ LENGTH z = dimindex (:N) /\\r
  99     LENGTH y = dimindex (:N) /\ LENGTH w = dimindex (:N) /\\r
 100     ZIP y w = yw /\ ZIP x z = xz /\\r
 101     ALL2 test_domain_xi xz yw ==>\r
 102     m_cell_domain (vector x, vector z:real^N) (vector y) (vector w)`,\r
 103   REPEAT GEN_TAC THEN STRIP_TAC THEN POP_ASSUM MP_TAC THEN\r
 104     SUBGOAL_THEN `LENGTH (xz:(real#real)list) = dimindex (:N) /\ LENGTH (yw:(real#real)list) = dimindex (:N)` ASSUME_TAC THENL\r
 105     [\r
 106       EXPAND_TAC "yw" THEN EXPAND_TAC "xz" THEN\r
 107         REPEAT (new_rewrite [] [] LENGTH_ZIP) THEN ASM_REWRITE_TAC[];\r
 108       ALL_TAC\r
 109     ] THEN\r
 110     rewrite [] [] ALL2_ALL_ZIP THEN ASM_REWRITE_TAC[m_cell_domain; GSYM ALL_EL] THEN DISCH_TAC THEN\r
 111     REWRITE_TAC[m_cell_domain] THEN GEN_TAC THEN DISCH_TAC THEN\r
 112     REPEAT (new_rewrite [] [] VECTOR_COMPONENT) THEN ASM_REWRITE_TAC[] THEN\r
 113     ABBREV_TAC `j = i - 1` THEN\r
 114     SUBGOAL_THEN `j < dimindex (:N)` ASSUME_TAC THENL\r
 115     [\r
 116       POP_ASSUM MP_TAC THEN POP_ASSUM MP_TAC THEN REWRITE_TAC[IN_NUMSEG] THEN ARITH_TAC;\r
 117       ALL_TAC\r
 118     ] THEN\r
 119     FIRST_X_ASSUM (MP_TAC o SPEC `j:num`) THEN REWRITE_TAC[test_domain_xi] THEN\r
 120     rewrite [] [] LENGTH_ZIP THEN ASM_REWRITE_TAC[] THEN\r
 121     rewrite [] [] EL_ZIP THEN ASM_REWRITE_TAC[] THEN\r
 122     EXPAND_TAC "xz" THEN EXPAND_TAC "yw" THEN\r
 123     REPEAT (new_rewrite [] [] EL_ZIP) THEN ASM_REWRITE_TAC[] THEN\r
 124     ARITH_TAC);;\r
 125 \r
 126 \r
 127 \r
 128 \r
 129 let float0_eq = FLOAT_TO_NUM_CONV (mk_float 0 Arith_options.min_exp);;\r
 130      \r
 131 \r
 132 (* array of theorems *)\r
 133 let mk_m_domain_array =\r
 134   let mk_m_domain n =\r
 135     let dimindex_th = dimindex_array.(n) in\r
 136     let n_ty = (hd o snd o dest_type o snd o dest_const o rand o lhand o concl) dimindex_th in\r
 137     let nty = `:N` in\r
 138       (UNDISCH_ALL o REWRITE_RULE[float0_eq] o DISCH_ALL o RULE o \r
 139          REWRITE_RULE[dimindex_th] o INST_TYPE[n_ty, nty]) MK_CELL_DOMAIN in\r
 140     Array.init (max_dim + 1) (fun i -> if i < 1 then TRUTH else mk_m_domain i);;\r
 141 \r
 142 \r
 143 \r
 144 let x_var_real_list = `x:(real)list` and\r
 145     y_var_real_list = `y:(real)list` and\r
 146     z_var_real_list = `z:(real)list` and\r
 147     w_var_real_list = `w:(real)list` and\r
 148     yw_var = `yw:(real#real)list` and\r
 149     xz_var = `xz:(real#real)list` and\r
 150     xz_pair_var = `xz:real#real` and\r
 151     yw_pair_var = `yw:real#real`;;\r
 152 \r
 153 \r
 154 \r
 155 \r
 156 let TEST_DOMAIN_XI' = (EQT_INTRO o RULE o prove)(`xz = (x,z) /\ yw = (y,w) /\\r
 157     x <= y /\ y <= z /\ y - x <= w1 /\ z - y <= w2 /\ w1 <= w /\ w2 <= w ==> test_domain_xi xz yw`,\r
 158   SIMP_TAC[test_domain_xi] THEN REAL_ARITH_TAC);;\r
 159 \r
 160 \r
 161 let eval_test_domain_xi pp test_domain_tm =\r
 162   let ltm, yw = dest_comb test_domain_tm in\r
 163   let xz = rand ltm in\r
 164   let x, z = dest_pair xz and\r
 165       y, w = dest_pair yw in\r
 166   let (<=) = (fun t1 t2 -> EQT_ELIM (float_le t1 t2)) and\r
 167       (-) = float_sub_hi pp in\r
 168   let x_le_y = x <= y and\r
 169       y_le_z = y <= z and\r
 170       yx_le_w1 = y - x and\r
 171       zy_le_w2 = z - y in\r
 172   let w1 = (rand o concl) yx_le_w1 and\r
 173       w2 = (rand o concl) zy_le_w2 in\r
 174   let w1_le_w = w1 <= w and\r
 175       w2_le_w = w2 <= w in\r
 176     (MY_PROVE_HYP (REFL xz) o MY_PROVE_HYP (REFL yw) o\r
 177       MY_PROVE_HYP x_le_y o MY_PROVE_HYP y_le_z o \r
 178        MY_PROVE_HYP yx_le_w1 o MY_PROVE_HYP zy_le_w2 o\r
 179        MY_PROVE_HYP w1_le_w o MY_PROVE_HYP w2_le_w o\r
 180        INST[x, x_var_real; y, y_var_real; z, z_var_real; w, w_var_real;\r
 181             w1, w1_var_real; w2, w2_var_real; \r
 182             xz, xz_pair_var; yw, yw_pair_var]) TEST_DOMAIN_XI';;\r
 183     \r
 184 \r
 185 (* mk_m_center_domain *)\r
 186 let mk_m_center_domain n pp x_list_tm z_list_tm =\r
 187   let x_list = dest_list x_list_tm and\r
 188       z_list = dest_list z_list_tm in\r
 189   let y_list =\r
 190     let ( * ) = (fun t1 t2 -> (rand o concl) (float_mul_eq t1 t2)) and\r
 191         (+) = (fun t1 t2 -> (rand o concl) (float_add_hi pp t1 t2)) in\r
 192       map2 (fun x y -> if x = y then x else float_inv2 * (x + y)) x_list z_list in\r
 193 \r
 194   let w_list =\r
 195     let (-) = (fun t1 t2 -> (rand o concl) (float_sub_hi pp t1 t2)) and\r
 196         max = (fun t1 t2 -> (rand o concl) (float_max t1 t2)) in\r
 197     let w1 = map2 (-) y_list x_list and\r
 198         w2 = map2 (-) z_list y_list in\r
 199       map2 max w1 w2 in\r
 200 \r
 201   let y_list_tm = mk_list (y_list, real_ty) and\r
 202       w_list_tm = mk_list (w_list, real_ty) in\r
 203 \r
 204   let yw_zip_th = eval_zip y_list_tm w_list_tm and\r
 205       xz_zip_th = eval_zip x_list_tm z_list_tm in\r
 206 \r
 207   let yw_list_tm = (rand o concl) yw_zip_th and\r
 208       xz_list_tm = (rand o concl) xz_zip_th in\r
 209 \r
 210   let len_x_th = eval_length x_list_tm and\r
 211       len_z_th = eval_length z_list_tm and\r
 212       len_y_th = eval_length y_list_tm and\r
 213       len_w_th = eval_length w_list_tm in\r
 214   let th0 = (MY_PROVE_HYP len_x_th o MY_PROVE_HYP len_z_th o\r
 215                MY_PROVE_HYP len_y_th o MY_PROVE_HYP len_w_th o\r
 216                MY_PROVE_HYP yw_zip_th o MY_PROVE_HYP xz_zip_th o\r
 217                INST[x_list_tm, x_var_real_list; z_list_tm, z_var_real_list;\r
 218                     y_list_tm, y_var_real_list; w_list_tm, w_var_real_list;\r
 219                     yw_list_tm, yw_var; xz_list_tm, xz_var]) mk_m_domain_array.(n) in\r
 220   let all_th = (EQT_ELIM o all2_conv_univ (eval_test_domain_xi pp) o hd o hyp) th0 in\r
 221     MY_PROVE_HYP all_th th0;;\r
 222 \r
 223 \r
 224 (*\r
 225 let n = 8;;\r
 226 let pp = 5;;\r
 227 \r
 228 let x_list_tm = mk_list (replicate one_float n, real_ty) and\r
 229     z_list_tm = mk_list (replicate two_float n, real_ty);;\r
 230 \r
 231 mk_m_center_domain n pp x_list_tm z_list_tm;;\r
 232 \r
 233 (* 10: 1.572 *)\r
 234 test 100 (mk_m_center_domain n pp x_list_tm) z_list_tm;;\r
 235 *)\r
 236 \r
 237 \r
 238 (***********************)\r
 239 \r
 240 let MK_M_TAYLOR_INTERVAL' = (RULE o MATCH_MP iffRL o SPEC_ALL) m_taylor_interval;;\r
 241 \r
 242 \r
 243 let get_types_and_vars n =\r
 244   let ty = n_type_array.(n) and\r
 245       xty = n_vector_type_array.(n) in\r
 246   let x_var = mk_var ("x", xty) and\r
 247       f_var = mk_var ("f", mk_fun_ty xty real_ty) and\r
 248       y_var = mk_var ("y", xty) and\r
 249       w_var = mk_var ("w", xty) and\r
 250       domain_var = mk_var ("domain", mk_type ("prod", [xty; xty])) in\r
 251     ty, xty, x_var, f_var, y_var, w_var, domain_var;;\r
 252 \r
 253 \r
 254 \r
 255 let dest_m_cell_domain domain_tm =\r
 256   let lhs, w_tm = dest_comb domain_tm in\r
 257   let lhs2, y_tm = dest_comb lhs in\r
 258     rand lhs2, y_tm, w_tm;;\r
 259 \r
 260 \r
 261 (**************************************************)\r
 262 \r
 263 (* Given a variable of the type `:real^N` returns the number N *)\r
 264 let get_dim = int_of_string o fst o dest_type o hd o tl o snd o dest_type o type_of;;\r
 265 \r
 266 \r
 267 (**********************)\r
 268 (* eval_m_taylor_poly *)\r
 269 \r
 270 let partial_pow = prove(`!i f n (y:real^N). lift o f differentiable at y ==>\r
 271                           partial i (\x. f x pow n) y = &n * f y pow (n - 1) * partial i f y`,\r
 272   REPEAT STRIP_TAC THEN\r
 273     SUBGOAL_THEN `(\x:real^N. f x pow n) = (\t. t pow n) o f` (fun th -> REWRITE_TAC[th]) THENL\r
 274     [\r
 275       ONCE_REWRITE_TAC[GSYM eq_ext] THEN REWRITE_TAC[o_THM];\r
 276       ALL_TAC\r
 277     ] THEN\r
 278     new_rewrite [] [] partial_uni_compose THENL\r
 279     [\r
 280       ASM_REWRITE_TAC[] THEN\r
 281         new_rewrite [] [] REAL_DIFFERENTIABLE_POW_ATREAL THEN REWRITE_TAC[REAL_DIFFERENTIABLE_ID];\r
 282       ALL_TAC\r
 283     ] THEN\r
 284     new_rewrite [] [] derivative_pow THEN REWRITE_TAC[REAL_DIFFERENTIABLE_ID; derivative_x] THEN\r
 285     REAL_ARITH_TAC);;\r
 286 \r
 287 \r
 288 let nth_diff2_pow = prove(`!n y. nth_diff_strong 2 (\x. x pow n) y`,\r
 289   REWRITE_TAC[nth_diff_strong2_eq] THEN REPEAT GEN_TAC THEN\r
 290     EXISTS_TAC `(:real)` THEN REWRITE_TAC[REAL_OPEN_UNIV; IN_UNIV] THEN GEN_TAC THEN\r
 291     new_rewrite [] [] REAL_DIFFERENTIABLE_POW_ATREAL THEN REWRITE_TAC[REAL_DIFFERENTIABLE_ID] THEN\r
 292     MATCH_MP_TAC differentiable_local THEN\r
 293     EXISTS_TAC `\x. &n * x pow (n - 1)` THEN EXISTS_TAC `(:real)` THEN\r
 294     REWRITE_TAC[REAL_OPEN_UNIV; IN_UNIV] THEN\r
 295     new_rewrite [] [] REAL_DIFFERENTIABLE_MUL_ATREAL THEN REWRITE_TAC[REAL_DIFFERENTIABLE_CONST] THENL\r
 296     [\r
 297       new_rewrite [] [] REAL_DIFFERENTIABLE_POW_ATREAL THEN REWRITE_TAC[REAL_DIFFERENTIABLE_ID];\r
 298       ALL_TAC\r
 299     ] THEN\r
 300     GEN_TAC THEN new_rewrite [] [] derivative_pow THEN REWRITE_TAC[REAL_DIFFERENTIABLE_ID] THEN\r
 301     REWRITE_TAC[derivative_x; REAL_MUL_RID]);;\r
 302 \r
 303 \r
 304 let diff2c_pow = prove(`!f n (x:real^N). diff2c f x ==> diff2c (\x. f x pow n) x`,\r
 305   REPEAT STRIP_TAC THEN\r
 306     SUBGOAL_THEN `(\x:real^N. f x pow n) = (\t. t pow n) o f` (fun th -> REWRITE_TAC[th]) THENL\r
 307     [\r
 308       ONCE_REWRITE_TAC[GSYM eq_ext] THEN REWRITE_TAC[o_THM];\r
 309       ALL_TAC\r
 310     ] THEN\r
 311     apply_tac diff2c_uni_compose THEN ASM_REWRITE_TAC[nth_diff2_pow] THEN\r
 312     REWRITE_TAC[nth_derivative2] THEN\r
 313     SUBGOAL_THEN `!n. derivative (\t. t pow n) = (\t. &n * t pow (n - 1))` ASSUME_TAC THENL\r
 314     [\r
 315       GEN_TAC THEN REWRITE_TAC[FUN_EQ_THM] THEN GEN_TAC THEN\r
 316         new_rewrite [] [] derivative_pow THEN REWRITE_TAC[REAL_DIFFERENTIABLE_ID] THEN\r
 317         REWRITE_TAC[derivative_x; REAL_MUL_RID];\r
 318       ALL_TAC\r
 319     ] THEN\r
 320     ASM_REWRITE_TAC[] THEN\r
 321     SUBGOAL_THEN `!n. derivative (\t. &n * t pow (n - 1)) = (\t. &n * derivative (\t. t pow (n - 1)) t)` ASSUME_TAC THENL\r
 322     [\r
 323       GEN_TAC THEN REWRITE_TAC[FUN_EQ_THM] THEN GEN_TAC THEN\r
 324         new_rewrite [] [] derivative_scale THEN REWRITE_TAC[] THEN\r
 325         MATCH_MP_TAC REAL_DIFFERENTIABLE_POW_ATREAL THEN REWRITE_TAC[REAL_DIFFERENTIABLE_ID];\r
 326       ALL_TAC\r
 327     ] THEN\r
 328     ASM_REWRITE_TAC[] THEN\r
 329     REPEAT (MATCH_MP_TAC REAL_CONTINUOUS_LMUL) THEN\r
 330     MATCH_MP_TAC REAL_CONTINUOUS_POW THEN\r
 331     REWRITE_TAC[REAL_CONTINUOUS_AT_ID]);;\r
 332 \r
 333 \r
 334 let diff2c_domain_pow = prove(`!f n domain. diff2c_domain domain f ==> \r
 335                                diff2c_domain domain (\x. f x pow n)`,\r
 336   REWRITE_TAC[diff2c_domain] THEN REPEAT STRIP_TAC THEN ASM_SIMP_TAC[diff2c_pow]);;\r
 337 \r
 338 \r
 339 let diff2c_domain_tm = `diff2c_domain domain`;;\r
 340 \r
 341 let rec gen_diff2c_domain_poly poly_tm =\r
 342   let x_var, expr = dest_abs poly_tm in\r
 343   let n = (int_of_string o fst o dest_type o hd o tl o snd o dest_type o type_of) x_var in\r
 344   let diff2c_tm = mk_icomb (diff2c_domain_tm, poly_tm) in\r
 345     if frees expr = [] then\r
 346       (* const *)\r
 347       (SPEC_ALL o ISPEC expr o inst_first_type_var (n_type_array.(n))) diff2c_domain_const\r
 348     else\r
 349       let lhs, r_tm = dest_comb expr in\r
 350         if lhs = neg_op_real then\r
 351           (* -- *)\r
 352           let r_th = gen_diff2c_domain_poly (mk_abs (x_var, r_tm)) in\r
 353             prove(diff2c_tm, MATCH_MP_TAC diff2c_domain_neg THEN REWRITE_TAC[r_th])\r
 354         else\r
 355           let op, l_tm = dest_comb lhs in\r
 356           let name = (fst o dest_const) op in\r
 357             if name = "$" then\r
 358               (* x$k *)\r
 359               let dim_th = dimindex_array.(n) in\r
 360                 prove(diff2c_tm, MATCH_MP_TAC diff2c_domain_x THEN\r
 361                         REWRITE_TAC[IN_NUMSEG; dim_th] THEN ARITH_TAC)\r
 362             else\r
 363               let l_th = gen_diff2c_domain_poly (mk_abs (x_var, l_tm)) in\r
 364                 if name = "real_pow" then\r
 365                   (* f pow n *)\r
 366                   prove(diff2c_tm, MATCH_MP_TAC diff2c_domain_pow THEN REWRITE_TAC[l_th])\r
 367                 else\r
 368                   let r_th = gen_diff2c_domain_poly (mk_abs (x_var, r_tm)) in\r
 369                     prove(diff2c_tm,\r
 370                           MAP_FIRST apply_tac [diff2c_domain_add; diff2c_domain_sub; diff2c_domain_mul] THEN\r
 371                             REWRITE_TAC[l_th; r_th]);;\r
 372 \r
 373 \r
 374 let gen_diff2c_poly =\r
 375   let th_imp = prove(`!f. (!domain. diff2c_domain domain f) ==> !x:real^N. diff2c f x`,\r
 376                      REWRITE_TAC[diff2c_domain] THEN REPEAT STRIP_TAC THEN FIRST_X_ASSUM MATCH_MP_TAC THEN\r
 377                        EXISTS_TAC `(x:real^N, x:real^N)` THEN\r
 378                        REWRITE_TAC[INTERVAL_SING; IN_SING]) in\r
 379     fun poly_tm ->\r
 380       (MATCH_MP th_imp o GEN_ALL o gen_diff2c_domain_poly) poly_tm;;\r
 381 \r
 382 \r
 383 let gen_diff_poly =\r
 384   let th_imp = prove(`!f. (!domain. diff2c_domain domain f) ==> !x:real^N. lift o f differentiable at x`,\r
 385                      REWRITE_TAC[diff2c_domain; diff2c; diff2] THEN REPEAT STRIP_TAC THEN\r
 386                        FIRST_X_ASSUM (MP_TAC o SPECL [`x:real^N, x:real^N`; `x:real^N`]) THEN\r
 387                        REWRITE_TAC[INTERVAL_SING; IN_SING] THEN case THEN REPEAT STRIP_TAC THEN\r
 388                        FIRST_X_ASSUM (MP_TAC o SPEC `x:real^N`) THEN ASM_SIMP_TAC[]) in\r
 389     fun poly_tm ->\r
 390       (MATCH_MP th_imp o GEN_ALL o gen_diff2c_domain_poly) poly_tm;;\r
 391 \r
 392 \r
 393 let in_tm = `IN`;;\r
 394 \r
 395 let gen_partial_poly i poly_tm =\r
 396   let i_tm = mk_small_numeral i in\r
 397   let rec gen_rec poly_tm =\r
 398     let x_var, expr = dest_abs poly_tm in\r
 399     let n = (int_of_string o fst o dest_type o hd o tl o snd o dest_type o type_of) x_var in\r
 400       if frees expr = [] then\r
 401         (* const *)\r
 402         (SPECL [i_tm; expr] o inst_first_type_var (n_type_array.(n))) partial_const\r
 403       else\r
 404         let lhs, r_tm = dest_comb expr in\r
 405           if lhs = neg_op_real then\r
 406             (* -- *)\r
 407             let r_poly = mk_abs (x_var, r_tm) in\r
 408             let r_diff = (SPEC_ALL o gen_diff_poly) r_poly and\r
 409                 r_partial = gen_rec r_poly in\r
 410             let th0 = SPEC i_tm (MATCH_MP partial_neg r_diff) in \r
 411               REWRITE_RULE[r_partial] th0\r
 412             else\r
 413               let op, l_tm = dest_comb lhs in\r
 414               let name = (fst o dest_const) op in\r
 415                 if name = "$" then\r
 416                   (* comp *)\r
 417                   let dim_th = dimindex_array.(n) in\r
 418                   let dim_tm = (lhand o concl) dim_th in\r
 419                   let i_eq_k = NUM_EQ_CONV (mk_eq (i_tm, r_tm)) in\r
 420                   let int_tm = mk_binop `..` `1` dim_tm in\r
 421                   let k_in_dim = prove(mk_comb (mk_icomb(in_tm, r_tm), int_tm),\r
 422                                        REWRITE_TAC[IN_NUMSEG; dim_th] THEN ARITH_TAC) in\r
 423                     (REWRITE_RULE[i_eq_k] o MATCH_MP (SPECL [r_tm; i_tm] partial_x)) k_in_dim\r
 424                 else\r
 425                   let l_poly = mk_abs (x_var, l_tm) in\r
 426                   let l_partial = gen_rec l_poly in\r
 427                   let l_diff = (SPEC_ALL o gen_diff_poly) l_poly in\r
 428                     if name = "real_pow" then\r
 429                       (* f pow n *)\r
 430                       let th0 = SPECL [i_tm; r_tm] (MATCH_MP partial_pow l_diff) in\r
 431                         REWRITE_RULE[l_partial] th0\r
 432                     else\r
 433                       let r_poly = mk_abs (x_var, r_tm) in\r
 434                       let r_partial = gen_rec r_poly in\r
 435                       let r_diff = (SPEC_ALL o gen_diff_poly) r_poly in\r
 436                       let imp_th = assoc op [add_op_real, partial_add; \r
 437                                              sub_op_real, partial_sub;\r
 438                                              mul_op_real, partial_mul] in\r
 439                       let th0 = SPEC i_tm (MATCH_MP (MATCH_MP imp_th l_diff) r_diff) in\r
 440                         REWRITE_RULE[l_partial; r_partial] th0 in\r
 441     \r
 442   let th1 = gen_rec poly_tm in\r
 443   let th2 = ((NUM_REDUCE_CONV THENC REWRITE_CONV[DECIMAL] THENC REAL_POLY_CONV) o rand o concl) th1 in\r
 444     (REWRITE_RULE[ETA_AX] o ONCE_REWRITE_RULE[eq_ext] o GEN_ALL) (TRANS th1 th2);;\r
 445 \r
 446 \r
 447 let gen_partial2_poly i j poly_tm =\r
 448   let partial_j = gen_partial_poly j poly_tm in\r
 449   let partial_ij = gen_partial_poly i (rand (concl partial_j)) in\r
 450   let pi = (rator o lhand o concl) partial_ij in\r
 451     REWRITE_RULE[GSYM partial2] (TRANS (AP_TERM pi partial_j) partial_ij);;\r
 452 \r
 453 \r
 454 (*\r
 455 let poly_tm = `\x:real^2. (x$1 - x$2) pow 3 * (x$1 - x$2) * x$1 pow 2`;;\r
 456 gen_partial2_poly 2 2 poly_tm;;\r
 457 *)\r
 458 \r
 459 let expr_to_vector_fun =\r
 460   let comp_op = `$` in\r
 461     fun expr_tm ->\r
 462       let vars = List.sort Pervasives.compare (frees expr_tm) in\r
 463       let n = length vars in\r
 464       let x_var = mk_var ("x", n_vector_type_array.(n)) in\r
 465       let x_tm = mk_icomb (comp_op, x_var) in\r
 466       let vars2 = map (fun i -> mk_comb (x_tm, mk_small_numeral i)) (1--n) in\r
 467         mk_abs (x_var, subst (zip vars2 vars) expr_tm);;\r
 468 \r
 469 (*\r
 470 let delta_x4_poly = expr_to_vector_fun ((rand o concl o SPEC_ALL) Sphere.delta_x4);;\r
 471 gen_partial2_poly 2 3 delta_x4_poly;;\r
 472 \r
 473 map (fun i -> map (fun j -> gen_partial2_poly i j poly_delta_x4) (1--6)) (1--6);;\r
 474 *)\r
 475 \r
 476 \r
 477 \r
 478 \r
 479 \r
 480 (********************************************)\r
 481 \r
 482 let x_var_names = Array.init (max_dim + 1) (fun i -> "x"^(string_of_int i)) and\r
 483     y_var_names = Array.init (max_dim + 1) (fun i -> "y"^(string_of_int i)) and\r
 484     z_var_names = Array.init (max_dim + 1) (fun i -> "z"^(string_of_int i)) and\r
 485     w_var_names = Array.init (max_dim + 1) (fun i -> "w"^(string_of_int i));;\r
 486 \r
 487 let x_vars_array = Array.init (max_dim + 1) (fun i -> mk_var(x_var_names.(i), real_ty)) and\r
 488     y_vars_array = Array.init (max_dim + 1) (fun i -> mk_var(y_var_names.(i), real_ty)) and\r
 489     z_vars_array = Array.init (max_dim + 1) (fun i -> mk_var(z_var_names.(i), real_ty)) and\r
 490     w_vars_array = Array.init (max_dim + 1) (fun i -> mk_var(w_var_names.(i), real_ty));;\r
 491 \r
 492 let df_vars_array = Array.init (max_dim + 1) (fun i -> mk_var ("df"^(string_of_int i), real_pair_ty));;\r
 493 let dd_vars_array = Array.init (max_dim + 1) (fun i ->\r
 494    Array.init (max_dim + 1) (fun j -> mk_var ("dd"^(string_of_int i)^(string_of_int j), real_pair_ty)));;\r
 495 \r
 496 let dest_vector = dest_list o rand;;\r
 497 \r
 498 let mk_vector list_tm =\r
 499   let n = (length o dest_list) list_tm in\r
 500   let ty = (hd o snd o dest_type o type_of) list_tm in\r
 501   let vec = mk_const ("vector", [ty, aty; n_type_array.(n), nty]) in\r
 502     mk_comb (vec, list_tm);;\r
 503 \r
 504 let mk_vector_list list =\r
 505     mk_vector (mk_list (list, type_of (hd list)));;\r
 506 \r
 507 \r
 508 let el_thms_array =\r
 509   let el_tm = `EL : num->(A)list->A` in\r
 510   let gen0 n =\r
 511     let e_list = mk_list (map (fun i -> mk_var ("e"^(string_of_int i), aty)) (1--n), aty) in\r
 512     let el0_th = REWRITE_CONV[EL; HD] (mk_binop el_tm `0` e_list) in\r
 513       Array.create n el0_th in\r
 514   let array = Array.init (max_dim + 1) gen0 in\r
 515   let gen_i n i =\r
 516     let e_list = (rand o lhand o concl) array.(n).(i) in\r
 517     let prev_thm = array.(n - 1).(i - 1) in\r
 518     let i_tm = mk_small_numeral i in\r
 519     let prev_i = num_CONV i_tm in\r
 520     let el_th = REWRITE_CONV[prev_i; EL; HD; TL; prev_thm] (mk_binop el_tm i_tm e_list) in\r
 521       array.(n).(i) <- el_th in\r
 522   let _ = map (fun n -> map (fun i -> gen_i n i) (1--(n - 1))) (2--max_dim) in\r
 523     array;;\r
 524 \r
 525 \r
 526 let gen_comp_thm n i =\r
 527   let i_tm = mk_small_numeral i and\r
 528       x_list = mk_list (map (fun i -> mk_var("x"^(string_of_int i), aty)) (1--n), aty) in\r
 529   let th0 = (ISPECL [x_list; i_tm] o inst_first_type_var (n_type_array.(n))) VECTOR_COMPONENT in\r
 530   let th1 = (CONV_RULE NUM_REDUCE_CONV o REWRITE_RULE[IN_NUMSEG; dimindex_array.(n)]) th0 in\r
 531     REWRITE_RULE[el_thms_array.(n).(i - 1)] th1;;\r
 532 \r
 533 let comp_thms_array = Array.init (max_dim + 1)\r
 534   (fun n -> Array.init (n + 1)\r
 535        (fun i -> if i < 1 or n < 1 then TRUTH else gen_comp_thm n i));;\r
 536 \r
 537 (*****************************)\r
 538 \r
 539 let eval_diff2_poly diff2_domain_th =\r
 540   fun xx zz ->\r
 541     let domain_tm = mk_pair (xx, zz) in\r
 542       INST[domain_tm, mk_var ("domain", type_of domain_tm)] diff2_domain_th;;\r
 543 \r
 544 \r
 545 (*****************************)\r
 546 \r
 547 \r
 548 let gen_lin_approx_eq_thm n =\r
 549   let ty = n_type_array.(n) in\r
 550   let df_vars = map (fun i -> df_vars_array.(i)) (1--n) in\r
 551   let df_bounds_list = mk_list (df_vars, real_pair_ty) in\r
 552   let th0 = (SPECL[f_bounds_var; df_bounds_list] o inst_first_type_var ty) m_lin_approx in\r
 553   let th1 = (CONV_RULE NUM_REDUCE_CONV o REWRITE_RULE[all_n]) th0 in\r
 554     th1;;\r
 555 \r
 556 \r
 557 let gen_lin_approx_poly_thm poly_tm diff_th partials =\r
 558   let x_var, _ = dest_abs poly_tm in\r
 559   let n = get_dim x_var in\r
 560   let lin_eq = (REWRITE_RULE partials o SPECL [poly_tm]) (gen_lin_approx_eq_thm n) in\r
 561   let x_vec = mk_vector_list (map (fun i -> x_vars_array.(i)) (1--n)) in\r
 562   let th1 = (REWRITE_RULE (Array.to_list comp_thms_array.(n)) o SPEC x_vec o REWRITE_RULE[diff_th]) lin_eq in\r
 563     th1;;\r
 564 \r
 565 \r
 566 let gen_lin_approx_poly_thm0 poly_tm =\r
 567   let x_var, _ = dest_abs poly_tm in\r
 568   let n = get_dim x_var in\r
 569   let partials = map (fun i -> gen_partial_poly i poly_tm) (1--n) in \r
 570   let diff_th = gen_diff_poly poly_tm in\r
 571     gen_lin_approx_poly_thm poly_tm diff_th partials;;\r
 572 \r
 573 \r
 574 let eval_lin_approx pp0 lin_approx_th =\r
 575   let poly_tm, _, _, _ = (dest_lin_approx o lhand o concl) lin_approx_th in\r
 576   let x_var, _ = dest_abs poly_tm in\r
 577   let n = get_dim x_var in\r
 578   let th0 = lin_approx_th in\r
 579   let th1 = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o MATCH_MP iffRL) th0 in\r
 580   let build_eval int_hyp =\r
 581     let expr, b_var = dest_binary "interval_arith" int_hyp in\r
 582       (eval_constants pp0 o build_interval_fun) expr, b_var in\r
 583   let int_fs = map build_eval (hyp th1) in\r
 584     \r
 585   let rec split_rules i_list =\r
 586     match i_list with\r
 587       | [] -> ([], [])\r
 588       | ((i_fun, var_tm) :: es) -> \r
 589           let th_list, i_list' = split_rules es in\r
 590             match i_fun with\r
 591               | Int_const th -> (var_tm, th) :: th_list, i_list'\r
 592               | Int_var v -> (var_tm, INST[v, x_var_real] CONST_INTERVAL') :: th_list, i_list'\r
 593               | _ -> th_list, (var_tm, i_fun) :: i_list' in\r
 594 \r
 595   let const_th_list, i_list0 = split_rules int_fs in\r
 596   let th2 = itlist (fun (var_tm, th) th0 -> \r
 597                       let b_tm = rand (concl th) in\r
 598                         (MY_PROVE_HYP th o INST[b_tm, var_tm]) th0) const_th_list th1 in\r
 599   let v_list, i_list' = unzip i_list0 in\r
 600   let i_list = find_and_replace_all i_list' [] in\r
 601     fun pp vector_tm ->\r
 602       let x_vals = dest_vector vector_tm in\r
 603       if length x_vals <> n then failwith (sprintf "Wrong vector size; expected size: %d" n)\r
 604       else\r
 605         let x_ints = map mk_const_interval x_vals in\r
 606         let vars = map (fun i -> x_vars_array.(i)) (1--n) in\r
 607         let th3 = INST (zip x_vals vars) th2 in\r
 608         let i_vals = eval_interval_fun_list pp i_list (zip vars x_ints) in\r
 609           itlist2 (fun var_tm th th0 ->\r
 610                      let b_tm = rand (concl th) in\r
 611                        (MY_PROVE_HYP th o INST[b_tm, var_tm]) th0) v_list i_vals th3;;\r
 612 \r
 613 \r
 614 let eval_lin_approx_poly0 pp0 poly_tm =\r
 615   eval_lin_approx pp0 (gen_lin_approx_poly_thm0 poly_tm);;\r
 616 \r
 617 \r
 618 (*\r
 619 let pp = 10;;\r
 620 let eval_rd = eval_lin_approx_poly0 pp rd_poly;;\r
 621 let eval_rd_old = eval_lin_approx_poly0_old pp rd_poly;;\r
 622 let pi5 = (fst o dest_pair o rand o concl) pi_approx_array.(5);;\r
 623 let y_list = mk_list (replicate pi5 3, real_ty);;\r
 624 let th1 = eval_rd (mk_vector y_list);;\r
 625 let th2 = eval_rd_old (mk_vector y_list);;\r
 626 th1 = th2;;\r
 627 \r
 628 (* 1.216 (rd) *)\r
 629 test 100 eval_rd (mk_vector y_list);;\r
 630 (* 1.180 *)\r
 631 test 100 eval_rd_old (mk_vector y_list);;\r
 632 \r
 633 \r
 634 let pp = 10;;\r
 635 let eval_rd = eval_lin_approx_poly0 pp delta_y_poly;;\r
 636 let eval_rd_old = eval_lin_approx_poly0_old pp delta_y_poly;;\r
 637 let pi5 = (fst o dest_pair o rand o concl) pi_approx_array.(5);;\r
 638 let y_list = mk_list (replicate pi5 6, real_ty);;\r
 639 let th1 = eval_rd (mk_vector y_list);;\r
 640 let th2 = eval_rd_old (mk_vector y_list);;\r
 641 th1 = th2;;\r
 642 \r
 643 (* 4.888 *)\r
 644 test 10 eval_rd (mk_vector y_list);;\r
 645 (* 6.972 *)\r
 646 test 10 eval_rd_old (mk_vector y_list);;\r
 647 th1 = th2;;\r
 648 *)\r
 649 \r
 650 \r
 651 (*\r
 652 let pp = 5;;\r
 653 let test_poly = `\x:real^2. x$1 * x$2`;;\r
 654 let test_poly = `\x:real^2. &1`;;\r
 655 gen_lin_approx_poly_thm0 test_poly;;\r
 656 let eval_test = eval_lin_approx_poly0 pp test_poly;;\r
 657 let y_list = mk_list (replicate one_float 2, real_ty);;\r
 658 eval_test (mk_vector y_list);;\r
 659 *)\r
 660 \r
 661 \r
 662 (*******************************)\r
 663 \r
 664 (*\r
 665 let pp = 5;;\r
 666 let eval_rd = eval_lin_approx_poly pp rd_poly;;\r
 667 let pi5 = (fst o dest_pair o rand o concl) pi_approx_array.(5);;\r
 668 let y_list = mk_list (replicate pi5 3, real_ty);;\r
 669 eval_rd (mk_vector y_list);;\r
 670 \r
 671 let eval_schwefel = eval_lin_approx_poly pp schwefel_poly;;\r
 672 let y_list = mk_list (replicate pi5 3, real_ty);;\r
 673 eval_schwefel (mk_vector y_list);;\r
 674 \r
 675 \r
 676 let eval_heart = eval_lin_approx_poly pp heart_poly;;\r
 677 let y_list = mk_list (replicate pi5 8, real_ty);;\r
 678 eval_heart (mk_vector y_list);;\r
 679 \r
 680 let eval_magnetism = eval_lin_approx_poly pp magnetism_poly;;\r
 681 let y_list = mk_list (replicate pi5 7, real_ty);;\r
 682 eval_magnetism (mk_vector y_list);;\r
 683 \r
 684 \r
 685 \r
 686 gen_lin_approx_poly_thm magnetism_poly;;\r
 687 gen_lin_approx_poly_thm rd_poly;;\r
 688 *)\r
 689 \r
 690 \r
 691 (*************************************)\r
 692 \r
 693 \r
 694 let le_op_num = `(<=):num->num->bool`;;\r
 695 \r
 696 (* 1 <= i /\ i <= n <=> i = 1 \/ i = 2 \/ ... \/ i = n *)\r
 697 let i_int_array =\r
 698   let i_tm = `i:num` in\r
 699   let i_th0 = prove(`1 <= i /\ i <= SUC n <=> (1 <= i /\ i <= n) \/ i = SUC n`, ARITH_TAC) in\r
 700   let th1 = prove(`1 <= i /\ i <= 1 <=> i = 1`, ARITH_TAC) in\r
 701   let array = Array.create (max_dim + 1) th1 in\r
 702   let prove_next n =\r
 703     let n_tm = mk_small_numeral n in\r
 704     let prev_n = num_CONV n_tm in\r
 705     let tm = mk_conj (`1 <= i`, mk_binop le_op_num i_tm n_tm) in\r
 706     let th = REWRITE_CONV[prev_n; i_th0; array.(n - 1)] tm in\r
 707       array.(n) <- REWRITE_RULE[SYM prev_n; DISJ_ACI] th in\r
 708   let _ = map prove_next (2--max_dim) in\r
 709     array;;\r
 710 \r
 711 \r
 712 \r
 713 (* (!i. 1 <= i /\ i <= n ==> P i) <=> P 1 /\ P 2 /\ ... /\ P n *)\r
 714 let gen_in_interval =\r
 715   let th0 = prove(`(!i:num. (i = k \/ Q i) ==> P i) <=> (P k /\ (!i. Q i ==> P i))`, MESON_TAC[]) in\r
 716   let th1 = prove(`(!i:num. (i = k ==> P i)) <=> P k`, MESON_TAC[]) in\r
 717     fun n ->\r
 718       let n_tm = mk_small_numeral n and\r
 719           i_tm = `i:num` in\r
 720       let lhs1 = mk_conj (`1 <= i`, mk_binop le_op_num i_tm n_tm) in\r
 721       let lhs = mk_forall (i_tm, mk_imp (lhs1, `(P:num->bool) i`)) in\r
 722         REWRITE_CONV[i_int_array.(n); th0; th1] lhs;;\r
 723 \r
 724 \r
 725 let gen_second_bounded_eq_thm n =\r
 726   let ty, _, x_var, _, _, _, domain_var = get_types_and_vars n in\r
 727   let dd_vars = map (fun i -> map (fun j -> dd_vars_array.(j).(i)) (1--i)) (1--n) in\r
 728   let dd_bounds_list = mk_list (map (fun l -> mk_list (l, real_pair_ty)) dd_vars, real_pair_list_ty) in\r
 729   let th0 = (SPECL[domain_var; dd_bounds_list] o inst_first_type_var ty) second_bounded in\r
 730   let th1 = (CONV_RULE NUM_REDUCE_CONV o REWRITE_RULE[all_n]) th0 in\r
 731     th1;;\r
 732 \r
 733 \r
 734 let gen_second_bounded_poly_thm poly_tm partials2 =\r
 735   let x_var, _ = dest_abs poly_tm in\r
 736   let n = get_dim x_var in\r
 737   let partials2' = List.flatten partials2 in\r
 738   let second_th = (REWRITE_RULE partials2' o SPECL [poly_tm]) (gen_second_bounded_eq_thm n) in\r
 739     second_th;;\r
 740 \r
 741 \r
 742 let gen_second_bounded_poly_thm0 poly_tm =\r
 743   let x_var, _ = dest_abs poly_tm in\r
 744   let n = get_dim x_var in\r
 745   let partials = map (fun i -> gen_partial_poly i poly_tm) (1--n) in\r
 746   let get_partial i eq_th = \r
 747     let partial_i = gen_partial_poly i (rand (concl eq_th)) in\r
 748     let pi = (rator o lhand o concl) partial_i in\r
 749       REWRITE_RULE[GSYM partial2] (TRANS (AP_TERM pi eq_th) partial_i) in\r
 750   let partials2 = map (fun th, i -> map (fun j -> get_partial j th) (1--i)) (zip partials (1--n)) in\r
 751     gen_second_bounded_poly_thm poly_tm partials2;;\r
 752 \r
 753 (*  let eq_th = TAUT `(P ==> Q /\ R) <=> ((P ==> Q) /\ (P ==> R))` in\r
 754     REWRITE_RULE[eq_th; FORALL_AND_THM; GSYM m_bounded_on_int] second_th;;*)\r
 755 \r
 756 \r
 757 (* eval_second_bounded *)\r
 758 let eval_second_bounded pp0 second_bounded_th =\r
 759   let poly_tm = (lhand o rator o lhand o concl) second_bounded_th in\r
 760   let th0 = second_bounded_th in\r
 761   let n = (get_dim o fst o dest_abs) poly_tm in\r
 762   let x_vector = mk_vector_list (map (fun i -> x_vars_array.(i)) (1--n)) and\r
 763       z_vector = mk_vector_list (map (fun i -> z_vars_array.(i)) (1--n)) in\r
 764   let _, _, _, _, _, _, domain_var = get_types_and_vars n in\r
 765 \r
 766   let th1 = INST[mk_pair (x_vector, z_vector), domain_var] th0 in\r
 767   let th2 = REWRITE_RULE[IN_INTERVAL; dimindex_array.(n)] th1 in\r
 768   let th3 = REWRITE_RULE[gen_in_interval n; GSYM interval_arith] th2 in\r
 769   let th4 = (REWRITE_RULE[CONJ_ACI] o REWRITE_RULE (Array.to_list comp_thms_array.(n))) th3 in\r
 770   let final_th0 = (UNDISCH_ALL o MATCH_MP iffRL) th4 in\r
 771 \r
 772   let x_var, h_tm = (dest_forall o hd o hyp) final_th0 in\r
 773   let _, h2 = dest_imp h_tm in\r
 774   let concl_ints = striplist dest_conj h2 in\r
 775 \r
 776   let i_funs = map (fun int -> \r
 777                       let expr, var = dest_interval_arith int in\r
 778                         (eval_constants pp0 o build_interval_fun) expr, var) concl_ints in\r
 779 \r
 780   let rec split_rules i_list =\r
 781     match i_list with\r
 782       | [] -> ([], [])\r
 783       | ((i_fun, var_tm) :: es) -> \r
 784           let th_list, i_list' = split_rules es in\r
 785             match i_fun with\r
 786               | Int_const th -> (var_tm, th) :: th_list, i_list'\r
 787 (*            | Int_var v -> (var_tm, INST[v, x_var_real] CONST_INTERVAL') :: th_list, i_list' *)\r
 788               | _ -> th_list, (var_tm, i_fun) :: i_list' in\r
 789 \r
 790   let const_th_list, i_list0 = split_rules i_funs in\r
 791   let th5 = itlist (fun (var_tm, th) th0 -> \r
 792                       let b_tm = rand (concl th) in\r
 793                         (REWRITE_RULE[th] o INST[b_tm, var_tm]) th0) const_th_list (SYM th4) in\r
 794   let final_th = REWRITE_RULE[GSYM IMP_IMP] th5 in\r
 795   let v_list, i_list' = unzip i_list0 in\r
 796   let i_list = find_and_replace_all i_list' [] in\r
 797 \r
 798     fun pp x_vector_tm z_vector_tm ->\r
 799       let x_vals = dest_vector x_vector_tm and\r
 800           z_vals = dest_vector z_vector_tm in\r
 801         if length x_vals <> n or length z_vals <> n then \r
 802           failwith (sprintf "Wrong vector size; expected size: %d" n)\r
 803         else\r
 804           let x_vars = map (fun i -> x_vars_array.(i)) (1--n) and\r
 805               z_vars = map (fun i -> z_vars_array.(i)) (1--n) in\r
 806 \r
 807           let inst_th = (INST (zip x_vals x_vars) o INST (zip z_vals z_vars)) final_th in\r
 808             if (not o is_eq) (concl inst_th) then inst_th \r
 809             else\r
 810               let x_var, lhs = (dest_forall o lhand o concl) inst_th in\r
 811               let hs = (butlast o striplist dest_imp) lhs in\r
 812               let vars = map (rand o rator) hs in\r
 813               let int_vars = zip vars (map ASSUME hs) in\r
 814 \r
 815               let dd_ints = eval_interval_fun_list pp i_list int_vars in\r
 816               let inst_dd = map2 (fun var th -> (rand o concl) th, var) v_list dd_ints in\r
 817               let inst_th2 = INST inst_dd inst_th in\r
 818 \r
 819               let conj_th = end_itlist CONJ dd_ints in\r
 820               let lhs_th = GEN x_var (itlist DISCH hs conj_th) in\r
 821                 EQ_MP inst_th2 lhs_th;;\r
 822 \r
 823 \r
 824 \r
 825 let eval_second_bounded_poly0 pp0 poly_tm =\r
 826   eval_second_bounded pp0 (gen_second_bounded_poly_thm0 poly_tm);;\r
 827 \r
 828 \r
 829 (****************************)\r
 830 \r
 831 (*\r
 832 needs "../formal_lp/formal_interval/m_examples_poly.hl";;\r
 833 \r
 834 let pi5 = (fst o dest_pair o rand o concl) pi_approx_array.(5);;\r
 835 let pp = 5;;\r
 836 \r
 837 (* n = 3 *)\r
 838 \r
 839 let n = 3;;\r
 840 let x_tm = mk_vector_list (replicate one_float n) and\r
 841     z_tm = mk_vector_list (replicate pi5 n);;\r
 842 \r
 843 let eval_poly = eval_second_bounded_poly0 pp schwefel_poly;;\r
 844 eval_poly pp x_tm z_tm;;\r
 845 (* 10: 0.928 (0.688 after find_and_replace_all) *)\r
 846 test 100 (eval_poly pp x_tm) z_tm;;\r
 847 \r
 848 \r
 849 let eval_poly = eval_second_bounded_poly0 pp rd_poly;;\r
 850 eval_poly pp x_tm z_tm;;\r
 851 (* 10: 0.08 *)\r
 852 test 1000 (eval_poly x_tm) z_tm;;\r
 853 \r
 854 \r
 855 (* n = 4 *)\r
 856 let n = 4;;\r
 857 let x_tm = mk_vector_list (replicate one_float n) and\r
 858     z_tm = mk_vector_list (replicate pi5 n);;\r
 859 \r
 860 let eval_poly = eval_second_bounded_poly0 pp lv_poly;;\r
 861 eval_poly 2 x_tm z_tm;;\r
 862 (* 10: 0.460 (0.232 after find_and_replace_all) *)\r
 863 test 100 (eval_poly x_tm) z_tm;;\r
 864 \r
 865 \r
 866 (* n = 6 *)\r
 867 let n = 6;;\r
 868 let x_tm = mk_vector_list (replicate one_float n) and\r
 869     z_tm = mk_vector_list (replicate pi5 n);;\r
 870 \r
 871 let eval_poly = eval_second_bounded_poly0 pp delta_x4_poly;;\r
 872 eval_poly x_tm z_tm;;\r
 873 (* 10: 0.168 *)\r
 874 test 1000 (eval_poly x_tm) z_tm;;\r
 875 \r
 876 (* n = 7 *)\r
 877 let n = 7;;\r
 878 let x_tm = mk_vector_list (replicate one_float n) and\r
 879     z_tm = mk_vector_list (replicate pi5 n);;\r
 880 \r
 881 let eval_poly = eval_second_bounded_poly0 pp magnetism_poly;;\r
 882 eval_poly 3 x_tm z_tm;;\r
 883 (* 10: 0.212 *)\r
 884 test 1000 (eval_poly x_tm) z_tm;;\r
 885 \r
 886 \r
 887 (* n = 8 *)\r
 888 let n = 8;;\r
 889 let x_tm = mk_vector_list (replicate one_float n) and\r
 890     z_tm = mk_vector_list (replicate pi5 n);;\r
 891 \r
 892 let eval_poly = eval_second_bounded_poly0 pp heart_poly;;\r
 893 eval_poly x_tm z_tm;;\r
 894 (* 10: 7.03 (3.272 after find_and_replace_all) *)\r
 895 test 100 (eval_poly x_tm) z_tm;;\r
 896 *)\r
 897 \r
 898 \r
 899 (*************************************)\r
 900 \r
 901 let eval_m_taylor pp0 diff2c_th lin_th second_th =\r
 902   let poly_tm = (rand o concl) diff2c_th in\r
 903   let n = (get_dim o fst o dest_abs) poly_tm in\r
 904   let eval_lin = eval_lin_approx pp0 lin_th and\r
 905       eval_second = eval_second_bounded pp0 second_th in\r
 906 \r
 907   let ty, _, x_var, f_var, y_var, w_var, domain_var = get_types_and_vars n in\r
 908   let th0 = (SPEC_ALL o inst_first_type_var ty) m_taylor_interval in\r
 909   let th1 = INST[poly_tm, f_var] th0 in\r
 910   let th2 = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o MATCH_MP iffRL o REWRITE_RULE[diff2c_th]) th1 in\r
 911 \r
 912     fun p_lin p_second domain_th ->\r
 913       let domain_tm, y_tm, w_tm = dest_m_cell_domain (concl domain_th) in\r
 914       let x_tm, z_tm = dest_pair domain_tm in\r
 915 \r
 916       let lin_th = eval_lin p_lin y_tm and\r
 917           second_th = eval_second p_second x_tm z_tm in\r
 918 \r
 919       let _, _, f_bounds, df_bounds_list = dest_lin_approx (concl lin_th) in\r
 920       let dd_bounds_list = (rand o concl) second_th in\r
 921       let df_var = mk_var ("d_bounds_list", type_of df_bounds_list) and\r
 922           dd_var = mk_var ("dd_bounds_list", type_of dd_bounds_list) in\r
 923 \r
 924         (MY_PROVE_HYP domain_th o MY_PROVE_HYP lin_th o MY_PROVE_HYP second_th o\r
 925            INST[domain_tm, domain_var; y_tm, y_var; w_tm, w_var;\r
 926                 f_bounds, f_bounds_var; df_bounds_list, df_var; dd_bounds_list, dd_var]) th2;;\r
 927 \r
 928 \r
 929 let eval_m_taylor_poly0 pp0 poly_tm =\r
 930   let diff2_th = gen_diff2c_domain_poly poly_tm in\r
 931   let x_var, _ = dest_abs poly_tm in\r
 932   let n = get_dim x_var in\r
 933   let partials = map (fun i -> gen_partial_poly i poly_tm) (1--n) in\r
 934   let get_partial i eq_th = \r
 935     let partial_i = gen_partial_poly i (rand (concl eq_th)) in\r
 936     let pi = (rator o lhand o concl) partial_i in\r
 937       REWRITE_RULE[GSYM partial2] (TRANS (AP_TERM pi eq_th) partial_i) in\r
 938   let partials2 = map2 (fun th i -> map (fun j -> get_partial j th) (1--i)) partials (1--n) in\r
 939   let second_th = gen_second_bounded_poly_thm poly_tm partials2 in\r
 940   let diff_th = gen_diff_poly poly_tm in\r
 941   let lin_th = gen_lin_approx_poly_thm poly_tm diff_th partials in\r
 942     eval_m_taylor pp0 diff2_th lin_th second_th;;\r
 943 \r
 944 \r
 945 (**************************)\r
 946 \r
 947 (*\r
 948 let pp = 5;;    \r
 949 let n = 3;;\r
 950 \r
 951 let x_list = mk_list (replicate one_float n, real_ty) and\r
 952     z_list = mk_list (replicate pi5 n, real_ty);;\r
 953 let domain_th = mk_m_center_domain n pp x_list z_list;;\r
 954 \r
 955 let poly_tm = rd_poly;;\r
 956 let eval_poly = eval_m_taylor_poly0 pp poly_tm;;\r
 957 eval_poly pp pp domain_th;;\r
 958 (* 10: 1.004 *)\r
 959 test 100 eval_poly domain_th;;\r
 960 \r
 961 (* n = 6 *)\r
 962 let n = 6;;\r
 963 \r
 964 let x_list = mk_list (replicate one_float n, real_ty) and\r
 965     z_list = mk_list (replicate pi5 n, real_ty);;\r
 966 let domain_th = mk_m_center_domain n pp x_list z_list;;\r
 967 \r
 968 let poly_tm = delta_x4_poly;;\r
 969 let eval_poly = eval_m_taylor_poly pp poly_tm;;\r
 970 eval_poly domain_th;;\r
 971 (* 10: 2.672 *)\r
 972 test 100 eval_poly domain_th;;\r
 973 \r
 974 \r
 975 \r
 976 (* n = 8 *)\r
 977 let n = 8;;\r
 978 \r
 979 let x_list = mk_list (replicate one_float n, real_ty) and\r
 980     z_list = mk_list (replicate pi5 n, real_ty);;\r
 981 let domain_th = mk_m_center_domain n pp x_list z_list;;\r
 982 \r
 983 let poly_tm = heart_poly;;\r
 984 let eval_poly = eval_m_taylor_poly0 pp poly_tm;;\r
 985 eval_poly pp pp domain_th;;\r
 986 (* 10: 2.208  *)\r
 987 test 10 eval_poly domain_th;;\r
 988 *)\r
 989 \r
 990 \r
 991 (******************************************)\r
 992 \r
 993 (* mk_eval_function *)\r
 994 \r
 995 let mk_eval_function_eq pp0 eq_th =\r
 996   let expr_tm = (rand o concl) eq_th in\r
 997   let tm0 = `!x:real^N. x IN interval [domain] ==> interval_arith (f x) f_bounds` in\r
 998   let n = (get_dim o fst o dest_abs) expr_tm in\r
 999   let x_vector = mk_vector_list (map (fun i -> x_vars_array.(i)) (1--n)) and\r
1000       z_vector = mk_vector_list (map (fun i -> z_vars_array.(i)) (1--n)) in\r
1001   let ty, _, _, _, _, _, domain_var = get_types_and_vars n and\r
1002       f_var = mk_var ("f", type_of expr_tm) in\r
1003   let th1 = (REWRITE_CONV[IN_INTERVAL] o subst[mk_pair(x_vector,z_vector), domain_var] o inst[ty, nty]) tm0 in\r
1004   let th2 = REWRITE_RULE [dimindex_array.(n)] th1 in\r
1005   let th3 = REWRITE_RULE [gen_in_interval n; GSYM interval_arith] th2 in\r
1006   let th4 = (REWRITE_RULE[GSYM IMP_IMP; CONJ_ACI] o REWRITE_RULE (Array.to_list comp_thms_array.(n))) th3 in\r
1007   let final_th0 = (CONV_RULE ((RAND_CONV o ONCE_DEPTH_CONV) BETA_CONV) o INST[expr_tm, f_var]) th4 in\r
1008   let x_var, h_tm = (dest_forall o rand o concl) final_th0 in\r
1009   let f_tm = (fst o dest_interval_arith o last o striplist dest_imp) h_tm in\r
1010   let i_fun = (eval_constants pp0 o build_interval_fun) f_tm in\r
1011   let i_list = find_and_replace_all [i_fun] [] in\r
1012   let final_th = (PURE_REWRITE_RULE[SYM eq_th] o SYM) final_th0 in\r
1013     fun pp x_vector_tm z_vector_tm ->\r
1014       let x_vals = dest_vector x_vector_tm and\r
1015           z_vals = dest_vector z_vector_tm in\r
1016         if length x_vals <> n or length z_vals <> n then \r
1017           failwith (sprintf "Wrong vector size; expected size: %d" n)\r
1018         else\r
1019           let x_vars = map (fun i -> x_vars_array.(i)) (1--n) and\r
1020               z_vars = map (fun i -> z_vars_array.(i)) (1--n) in\r
1021 \r
1022           let inst_th = (INST (zip x_vals x_vars) o INST (zip z_vals z_vars)) final_th in\r
1023           let x_var, lhs = (dest_forall o lhand o concl) inst_th in\r
1024           let hs = (butlast o striplist dest_imp) lhs in\r
1025           let vars = map (rand o rator) hs in\r
1026           let int_vars = zip vars (map ASSUME hs) in\r
1027           let eval_th = hd (eval_interval_fun_list pp i_list int_vars) in\r
1028           let f_bounds = (rand o concl) eval_th in\r
1029           let inst_th2 = INST[f_bounds, f_bounds_var] inst_th in\r
1030           let lhs_th = GEN x_var (itlist DISCH hs eval_th) in\r
1031             EQ_MP inst_th2 lhs_th;;\r
1032 \r
1033 \r
1034 let mk_eval_function pp0 expr_tm = mk_eval_function_eq pp0 (REFL expr_tm);;\r
1035 \r
1036 \r
1037 \r
1038 \r
1039 \r
1040 (********************************)\r
1041 (* m_taylor_error *)\r
1042 \r
1043 \r
1044 (*\r
1045 (* n = 4 *)\r
1046 let n = 4;;\r
1047 let x_list = mk_list (replicate one_float n, real_ty) and\r
1048     z_list = mk_list (replicate pi5 n, real_ty);;\r
1049 let domain_th = mk_m_center_domain n pp x_list z_list;;\r
1050 \r
1051 let eval_poly = eval_m_taylor_poly pp lv_poly;;\r
1052 eval_poly domain_th;;\r
1053 *)\r
1054 \r
1055 (**********)\r
1056 \r
1057 \r
1058 (* Sum of the list elements *)\r
1059 let ITLIST2_EQ_SUM = prove(`!(f:A->B->real) l1 l2. LENGTH l1 <= LENGTH l2 ==>\r
1060                                ITLIST2 (\x y z. f x y + z) l1 l2 (&0) =\r
1061                                sum (1..(LENGTH l1)) (\i. f (EL (i - 1) l1) (EL (i - 1) l2))`,\r
1062    GEN_TAC THEN\r
1063      LIST_INDUCT_TAC THEN LIST_INDUCT_TAC THEN REWRITE_TAC[LENGTH; ITLIST2_DEF] THEN TRY ARITH_TAC THENL\r
1064      [\r
1065        REWRITE_TAC[SUM_CLAUSES_NUMSEG; ARITH];\r
1066        REWRITE_TAC[SUM_CLAUSES_NUMSEG; ARITH];\r
1067        ALL_TAC\r
1068      ] THEN\r
1069      REWRITE_TAC[leqSS] THEN DISCH_TAC THEN\r
1070      FIRST_X_ASSUM (MP_TAC o SPEC `t':(B)list`) THEN ASM_REWRITE_TAC[] THEN\r
1071      DISCH_THEN (fun th -> REWRITE_TAC[TL; th]) THEN\r
1072      REWRITE_TAC[GSYM add1n] THEN\r
1073      new_rewrite [] [] SUM_ADD_SPLIT THEN REWRITE_TAC[ARITH] THEN\r
1074      REWRITE_TAC[TWO; add1n; SUM_SING_NUMSEG; subnn; EL; HD] THEN\r
1075      REWRITE_TAC[GSYM addn1; SUM_OFFSET; REAL_EQ_ADD_LCANCEL] THEN\r
1076      MATCH_MP_TAC SUM_EQ THEN move ["i"] THEN REWRITE_TAC[IN_NUMSEG] THEN DISCH_TAC THEN\r
1077      ASM_SIMP_TAC[ARITH_RULE `1 <= i ==> (i + 1) - 1 = SUC (i - 1)`; EL; TL]);;\r
1078 \r
1079 \r
1080 let interval_arith_abs_le = prove(`!x int y. interval_arith x int ==> iabs int <= y ==> abs x <= y`,\r
1081    GEN_TAC THEN case THEN REWRITE_TAC[interval_arith; IABS'] THEN REAL_ARITH_TAC);;\r
1082 \r
1083 let l_seq = new_definition `l_seq n m = TABLE (\i. n + i) ((m + 1) - n)`;;\r
1084 \r
1085 (* TODO: create a file containing theorems about tables *)\r
1086 \r
1087 let REVERSE_TABLE  = define `(REVERSE_TABLE (f:num->A) 0 = []) /\\r
1088    (REVERSE_TABLE f (SUC i) = CONS (f i)  ( REVERSE_TABLE f i))`;; \r
1089 \r
1090   \r
1091 let TABLE = new_definition `!(f:num->A) k. TABLE f k = REVERSE (REVERSE_TABLE f k)`;;\r
1092 \r
1093 \r
1094 let LENGTH_REVERSE_TABLE = prove(`!(f:num->A) n. LENGTH (REVERSE_TABLE f n) = n`,\r
1095    GEN_TAC THEN INDUCT_TAC THEN ASM_REWRITE_TAC[REVERSE_TABLE; LENGTH]);;\r
1096 \r
1097 \r
1098 let LENGTH_REVERSE = prove(`!(l:(A)list). LENGTH (REVERSE l) = LENGTH l`,\r
1099    MATCH_MP_TAC list_INDUCT THEN REWRITE_TAC[REVERSE] THEN\r
1100      REPEAT STRIP_TAC THEN\r
1101      ASM_REWRITE_TAC[LENGTH_APPEND; LENGTH] THEN\r
1102      ARITH_TAC);;\r
1103 \r
1104 \r
1105 let LENGTH_TABLE = prove(`!(f:num->A) n. LENGTH (TABLE f n) = n`,\r
1106    REWRITE_TAC[TABLE; LENGTH_REVERSE; LENGTH_REVERSE_TABLE]);;\r
1107 \r
1108 \r
1109 let LENGTH_L_SEQ = prove(`LENGTH (l_seq n m) = (m + 1) - n`, REWRITE_TAC[l_seq; LENGTH_TABLE]);;\r
1110 let EL_L_SEQ = prove(`!i m n. i < (m + 1) - n ==> EL i (l_seq n m) = n + i`,\r
1111   REWRITE_TAC[l_seq] THEN REPEAT STRIP_TAC THEN MATCH_MP_TAC Packing3.EL_TABLE THEN ASM_REWRITE_TAC[]);;\r
1112 let L_SEQ_NIL = prove(`!n m. l_seq n m = [] <=> (m < n)`,\r
1113   GEN_TAC THEN GEN_TAC THEN REWRITE_TAC[GSYM LENGTH_EQ_NIL; LENGTH_L_SEQ] THEN ARITH_TAC);;\r
1114 let L_SEQ_NN = prove(`!n. l_seq n n = [n]`,\r
1115   GEN_TAC THEN REWRITE_TAC[l_seq; ARITH_RULE `(n + 1) - n = 1`; ONE; TABLE; REVERSE_TABLE; REVERSE] THEN\r
1116     REWRITE_TAC[APPEND; ADD_0]);;\r
1117 \r
1118 let L_SEQ_CONS = prove(`!n m. n <= m ==> l_seq n m = CONS n (l_seq (n + 1) m)`,\r
1119   REPEAT STRIP_TAC THEN\r
1120     REWRITE_TAC[Packing3.LIST_EL_EQ; LENGTH_L_SEQ; LENGTH] THEN CONJ_TAC THENL\r
1121     [\r
1122       ASM_ARITH_TAC;\r
1123       ALL_TAC\r
1124     ] THEN\r
1125     case THENL\r
1126     [\r
1127       DISCH_TAC THEN new_rewrite [] [] EL_L_SEQ THEN ASM_REWRITE_TAC[EL; HD; addn0];\r
1128       move ["i"; "lt"] THEN\r
1129         new_rewrite [] [] EL_L_SEQ THEN ASM_REWRITE_TAC[EL; TL] THEN\r
1130         new_rewrite [] [] EL_L_SEQ THEN ASM_ARITH_TAC\r
1131     ]);;\r
1132 \r
1133 \r
1134 \r
1135 let ALL_N_ALL2 = prove(`!P (l:(A)list) i0.\r
1136                          (all_n i0 l P <=> if l = [] then T else ALL2 P (l_seq i0 ((i0 + LENGTH l) - 1)) l)`,\r
1137   GEN_TAC THEN LIST_INDUCT_TAC THEN GEN_TAC THEN REWRITE_TAC[all_n; NOT_CONS_NIL] THEN\r
1138     new_rewrite [] [] L_SEQ_CONS THEN REWRITE_TAC[LENGTH; ALL2] THEN TRY ARITH_TAC THEN\r
1139     FIRST_X_ASSUM (new_rewrite [] []) THEN TRY ARITH_TAC THEN\r
1140     REWRITE_TAC[addSn; addnS; addn1] THEN\r
1141     SPEC_TAC (`t:(A)list`, `t:(A)list`) THEN case THEN SIMP_TAC[NOT_CONS_NIL] THEN\r
1142     REWRITE_TAC[LENGTH; addn0] THEN\r
1143     MP_TAC (SPECL [`SUC i0`; `SUC i0 - 1`] L_SEQ_NIL) THEN\r
1144     REWRITE_TAC[ARITH_RULE `SUC i0 - 1 < SUC i0`] THEN DISCH_THEN (fun th -> REWRITE_TAC[th; ALL2]));;\r
1145     \r
1146 \r
1147 let ALL_N_EL = prove(`!P (l:(A)list) i0. all_n i0 l P <=> (!i. i < LENGTH l ==> P (i0 + i) (EL i l))`,\r
1148    REPEAT GEN_TAC THEN REWRITE_TAC[ALL_N_ALL2] THEN\r
1149      SPEC_TAC (`l:(A)list`, `l:(A)list`) THEN case THEN SIMP_TAC[NOT_CONS_NIL; LENGTH; ltn0] THEN\r
1150      REPEAT GEN_TAC THEN\r
1151      new_rewrite [] [] ALL2_ALL_ZIP THENL\r
1152      [\r
1153        REWRITE_TAC[LENGTH_L_SEQ; LENGTH] THEN ARITH_TAC;\r
1154        ALL_TAC\r
1155      ] THEN\r
1156      REWRITE_TAC[GSYM ALL_EL] THEN\r
1157      new_rewrite [] [] LENGTH_ZIP THENL\r
1158      [\r
1159        REWRITE_TAC[LENGTH_L_SEQ; LENGTH] THEN ARITH_TAC;\r
1160        ALL_TAC\r
1161      ] THEN\r
1162      REWRITE_TAC[LENGTH_L_SEQ; ARITH_RULE `((i0 + SUC a) - 1 + 1) - i0 = SUC a`] THEN\r
1163      EQ_TAC THEN REPEAT STRIP_TAC THENL\r
1164      [\r
1165        FIRST_X_ASSUM (MP_TAC o SPEC `i:num`) THEN ASM_REWRITE_TAC[] THEN\r
1166          new_rewrite [] [] EL_ZIP THENL\r
1167          [\r
1168            REWRITE_TAC[LENGTH_L_SEQ; LENGTH] THEN ASM_ARITH_TAC;\r
1169            ALL_TAC\r
1170          ] THEN\r
1171          REWRITE_TAC[] THEN new_rewrite [] [] EL_L_SEQ THEN ASM_ARITH_TAC;\r
1172        ALL_TAC\r
1173      ] THEN\r
1174      new_rewrite [] [] EL_ZIP THENL\r
1175      [\r
1176        REWRITE_TAC[LENGTH_L_SEQ; LENGTH] THEN ASM_ARITH_TAC;\r
1177        ALL_TAC\r
1178      ] THEN\r
1179      REWRITE_TAC[] THEN new_rewrite [] [] EL_L_SEQ THEN TRY ASM_ARITH_TAC THEN\r
1180      FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[]);;\r
1181      \r
1182 \r
1183 \r
1184 let LENGTH_BUTLAST = prove(`!l. LENGTH (BUTLAST l) = LENGTH l - 1`,\r
1185    MATCH_MP_TAC list_INDUCT THEN REWRITE_TAC[BUTLAST; LENGTH; ARITH] THEN REPEAT STRIP_TAC THEN\r
1186      COND_CASES_TAC THEN ASM_REWRITE_TAC[LENGTH; ARITH] THEN\r
1187      POP_ASSUM MP_TAC THEN REWRITE_TAC[GSYM LENGTH_EQ_NIL] THEN\r
1188      ARITH_TAC);;\r
1189 \r
1190 \r
1191 let EL_BUTLAST = prove(`!(l:(A)list) i. i < LENGTH l - 1 ==> EL i (BUTLAST l) = EL i l`,\r
1192    MATCH_MP_TAC list_INDUCT THEN REWRITE_TAC[BUTLAST; LENGTH] THEN REPEAT STRIP_TAC THEN\r
1193      COND_CASES_TAC THENL\r
1194      [\r
1195        UNDISCH_TAC `i < SUC (LENGTH (a1:(A)list)) - 1` THEN\r
1196          ASM_REWRITE_TAC[LENGTH] THEN ARITH_TAC;\r
1197        ALL_TAC\r
1198      ] THEN\r
1199      REWRITE_TAC[EL_CONS] THEN COND_CASES_TAC THEN ASM_REWRITE_TAC[] THEN\r
1200      FIRST_X_ASSUM (MP_TAC o SPEC `i - 1`) THEN\r
1201      ANTS_TAC THENL\r
1202      [\r
1203        ASM_ARITH_TAC;\r
1204        ALL_TAC\r
1205      ] THEN\r
1206      REWRITE_TAC[]);;\r
1207      \r
1208 \r
1209 \r
1210 let M_TAYLOR_ERROR_ITLIST2 = prove(`!f domain y w dd_bounds_list error. \r
1211              m_cell_domain domain y (vector w) ==>\r
1212              diff2c_domain domain f ==>\r
1213              second_bounded (f:real^N->real) domain dd_bounds_list ==>\r
1214              LENGTH w = dimindex (:N) ==>\r
1215                  LENGTH dd_bounds_list = dimindex (:N) ==>\r
1216                  all_n 1 dd_bounds_list (\i list. LENGTH list = i) ==>\r
1217                  ITLIST2 (\list x z. x * (x * iabs (LAST list)\r
1218                                           + &2 * ITLIST2 (\a b c. b * iabs a + c) (BUTLAST list) w (&0)) + z)\r
1219                  dd_bounds_list w (&0) <= error ==>\r
1220                  m_taylor_error f domain (vector w) error`,\r
1221    REPEAT GEN_TAC THEN REWRITE_TAC[second_bounded] THEN\r
1222      set_tac "s" `ITLIST2 _1 _2 _3 _4` THEN\r
1223      move ["domain"; "d2f"; "second"; "lw"; "ldd"; "ldd_all"; "s_le"] THEN\r
1224      ASM_SIMP_TAC[m_taylor_error_eq] THEN move ["x"; "x_in"] THEN\r
1225      SUBGOAL_THEN `!i. i IN 1..dimindex (:N) ==> &0 <= EL (i - 1) w` (LABEL_TAC "w_ge0") THENL\r
1226      [\r
1227        GEN_TAC THEN DISCH_TAC THEN REMOVE_THEN "domain" MP_TAC THEN new_rewrite [] [] pair_eq THEN\r
1228          REWRITE_TAC[m_cell_domain] THEN\r
1229          DISCH_THEN (MP_TAC o SPEC `i:num`) THEN ASM_REWRITE_TAC[] THEN\r
1230          ASM_SIMP_TAC[VECTOR_COMPONENT] THEN REAL_ARITH_TAC;\r
1231        ALL_TAC\r
1232      ] THEN\r
1233      MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `s:real` THEN ASM_REWRITE_TAC[] THEN\r
1234      EXPAND_TAC "s" THEN\r
1235      new_rewrite [] [] ITLIST2_EQ_SUM THEN ASM_REWRITE_TAC[LE_REFL] THEN\r
1236      MATCH_MP_TAC SUM_LE THEN REWRITE_TAC[FINITE_NUMSEG] THEN\r
1237      move ["i"; "i_in"] THEN ASM_SIMP_TAC[VECTOR_COMPONENT] THEN\r
1238      USE_THEN "i_in" (ASSUME_TAC o REWRITE_RULE[IN_NUMSEG]) THEN\r
1239      MATCH_MP_TAC REAL_LE_LMUL THEN ASM_SIMP_TAC[] THEN\r
1240      SUBGOAL_THEN `LENGTH (EL (i - 1) dd_bounds_list:(real#real)list) = i` (LABEL_TAC "len_i") THENL\r
1241      [\r
1242        REMOVE_THEN "ldd_all" MP_TAC THEN\r
1243          REWRITE_TAC[ALL_N_EL] THEN DISCH_THEN (MP_TAC o SPEC `i - 1`) THEN\r
1244          ANTS_TAC THENL\r
1245          [\r
1246            ASM_REWRITE_TAC[] THEN POP_ASSUM MP_TAC THEN ARITH_TAC;\r
1247            ALL_TAC\r
1248          ] THEN\r
1249          DISCH_THEN (fun th -> REWRITE_TAC[th]) THEN POP_ASSUM MP_TAC THEN ARITH_TAC;\r
1250        ALL_TAC\r
1251      ] THEN\r
1252      new_rewrite [] [] ITLIST2_EQ_SUM THEN ASM_REWRITE_TAC[] THENL\r
1253      [\r
1254        REWRITE_TAC[LENGTH_BUTLAST] THEN POP_ASSUM MP_TAC THEN POP_ASSUM MP_TAC THEN ARITH_TAC;\r
1255        ALL_TAC\r
1256      ] THEN\r
1257      MATCH_MP_TAC REAL_LE_ADD2 THEN CONJ_TAC THENL\r
1258      [\r
1259        new_rewrite [] [] LAST_EL THENL\r
1260          [\r
1261            ASM_REWRITE_TAC[GSYM LENGTH_EQ_NIL] THEN \r
1262              REMOVE_THEN "i_in" MP_TAC THEN REWRITE_TAC[IN_NUMSEG] THEN ARITH_TAC;\r
1263            ALL_TAC\r
1264          ] THEN\r
1265          MATCH_MP_TAC REAL_LE_LMUL THEN ASM_SIMP_TAC[] THEN\r
1266          FIRST_X_ASSUM (MP_TAC o SPEC `x:real^N`) THEN ASM_REWRITE_TAC[ALL_N_EL] THEN\r
1267          DISCH_THEN (MP_TAC o SPEC `i - 1`) THEN\r
1268          ANTS_TAC THENL\r
1269          [\r
1270            REMOVE_THEN "i_in" MP_TAC THEN REWRITE_TAC[IN_NUMSEG] THEN ARITH_TAC;\r
1271            ALL_TAC\r
1272          ] THEN\r
1273          DISCH_THEN (MP_TAC o SPEC `i - 1`) THEN ANTS_TAC THENL\r
1274          [\r
1275            UNDISCH_TAC `1 <= i /\ i <= dimindex (:N)` THEN ASM_REWRITE_TAC[] THEN ARITH_TAC;\r
1276            ALL_TAC\r
1277          ] THEN\r
1278          SUBGOAL_THEN `1 + i - 1 = i /\ 1 + i - 1 = i` (fun th -> REWRITE_TAC[th]) THENL\r
1279          [\r
1280            REPLICATE_TAC 3 (POP_ASSUM MP_TAC) THEN ARITH_TAC;\r
1281            ALL_TAC\r
1282          ] THEN\r
1283          REWRITE_TAC[partial2] THEN\r
1284          DISCH_THEN (MP_TAC o MATCH_MP interval_arith_abs_le) THEN\r
1285          DISCH_THEN MATCH_MP_TAC THEN REWRITE_TAC[REAL_LE_REFL];\r
1286        ALL_TAC\r
1287      ] THEN\r
1288      MATCH_MP_TAC REAL_LE_LMUL THEN CONJ_TAC THENL [ REAL_ARITH_TAC; ALL_TAC ] THEN\r
1289      ASM_REWRITE_TAC[LENGTH_BUTLAST] THEN\r
1290      MATCH_MP_TAC SUM_LE THEN REWRITE_TAC[FINITE_NUMSEG] THEN\r
1291      move ["j"; "j_in"] THEN\r
1292      SUBGOAL_THEN `j IN 1..dimindex (:N)` (LABEL_TAC "j_in2") THENL\r
1293      [\r
1294        REMOVE_THEN "i_in" MP_TAC THEN REMOVE_THEN "j_in" MP_TAC THEN\r
1295          REWRITE_TAC[IN_NUMSEG] THEN ARITH_TAC;\r
1296        ALL_TAC\r
1297      ] THEN\r
1298      ASM_SIMP_TAC[VECTOR_COMPONENT] THEN\r
1299      USE_THEN "j_in" (ASSUME_TAC o REWRITE_RULE[IN_NUMSEG]) THEN\r
1300      MATCH_MP_TAC REAL_LE_LMUL THEN ASM_SIMP_TAC[] THEN\r
1301      FIRST_X_ASSUM (MP_TAC o SPEC `x:real^N`) THEN ASM_REWRITE_TAC[ALL_N_EL] THEN\r
1302      DISCH_THEN (MP_TAC o SPEC `i - 1`) THEN ANTS_TAC THENL\r
1303      [\r
1304        REMOVE_THEN "i_in" MP_TAC THEN REWRITE_TAC[IN_NUMSEG] THEN ARITH_TAC;\r
1305        ALL_TAC\r
1306      ] THEN\r
1307      DISCH_THEN (MP_TAC o SPEC `j - 1`) THEN ANTS_TAC THENL\r
1308      [\r
1309        ASM_REWRITE_TAC[] THEN POP_ASSUM MP_TAC THEN ARITH_TAC;\r
1310        ALL_TAC\r
1311      ] THEN\r
1312      SUBGOAL_THEN `1 + j - 1 = j /\ 1 + i - 1 = i` (fun th -> REWRITE_TAC[th]) THENL\r
1313      [\r
1314        REPLICATE_TAC 3 (POP_ASSUM MP_TAC) THEN ARITH_TAC;\r
1315        ALL_TAC\r
1316      ] THEN\r
1317      REWRITE_TAC[partial2] THEN\r
1318      DISCH_THEN (MP_TAC o MATCH_MP interval_arith_abs_le) THEN\r
1319      DISCH_THEN MATCH_MP_TAC THEN\r
1320      new_rewrite [] [] EL_BUTLAST THENL\r
1321      [\r
1322        ASM_REWRITE_TAC[] THEN REMOVE_THEN "j_in" MP_TAC THEN REMOVE_THEN "i_in" MP_TAC THEN\r
1323          REWRITE_TAC[IN_NUMSEG] THEN ARITH_TAC;\r
1324        ALL_TAC\r
1325      ] THEN\r
1326      REWRITE_TAC[REAL_LE_REFL]);;\r
1327 \r
1328 \r
1329 let M_TAYLOR_ERROR_ITLIST2_ALT = prove(`!f domain y w f_lo f_hi d_bounds_list dd_bounds_list error. \r
1330          m_taylor_interval f domain y (vector w:real^N) (f_lo, f_hi) d_bounds_list dd_bounds_list ==>\r
1331          LENGTH w = dimindex (:N) ==>\r
1332              LENGTH dd_bounds_list = dimindex (:N) ==>\r
1333              all_n 1 dd_bounds_list (\i list. LENGTH list = i) ==>\r
1334              ITLIST2 (\list x z. x * (x * iabs (LAST list)\r
1335                                       + &2 * ITLIST2 (\a b c. b * iabs a + c) (BUTLAST list) w (&0)) + z)\r
1336              dd_bounds_list w (&0) <= error ==>\r
1337              m_taylor_error f domain (vector w) error`,\r
1338    REWRITE_TAC[m_taylor_interval] THEN REPEAT STRIP_TAC THEN\r
1339      MP_TAC (SPEC_ALL M_TAYLOR_ERROR_ITLIST2) THEN ASM_REWRITE_TAC[]);;\r
1340 \r
1341 \r
1342 \r
1343 (**********************************)\r
1344 \r
1345 (*\r
1346 let pp = 5;;\r
1347 let eval_poly = eval_m_taylor_poly0 pp lv_poly;;\r
1348 \r
1349 let n = 4;;\r
1350 let xx = mk_vector_list (replicate one_float 4) and\r
1351     zz = mk_vector_list (replicate two_float 4);;\r
1352 \r
1353 let domain4_th = mk_m_center_domain n pp (rand xx) (rand zz);;\r
1354 \r
1355 let m_taylor_th = eval_poly pp pp domain4_th;;\r
1356 *)\r
1357 \r
1358 (****************************)\r
1359 \r
1360 let M_TAYLOR_INTERVAL' = MY_RULE m_taylor_interval;;\r
1361 \r
1362 let dest_m_taylor m_taylor_tm =\r
1363   let ltm1, dd_bounds_list = dest_comb m_taylor_tm in\r
1364   let ltm2, d_bounds_list = dest_comb ltm1 in\r
1365   let ltm3, f_bounds = dest_comb ltm2 in\r
1366   let ltm4, w = dest_comb ltm3 in\r
1367   let ltm5, y = dest_comb ltm4 in\r
1368   let ltm6, domain = dest_comb ltm5 in\r
1369     rand ltm6, domain, y, w, f_bounds, d_bounds_list, dd_bounds_list;;\r
1370 \r
1371 let dest_m_taylor_thms n =\r
1372   let ty, xty, x_var, f_var, y_var, w_var, domain_var = get_types_and_vars n in\r
1373     fun m_taylor_th ->\r
1374       let f, domain, y, w, f_bounds, d_bounds_list, dd_bounds_list = dest_m_taylor (concl m_taylor_th) in\r
1375       let th0 = (INST[f, f_var; domain, domain_var; y, y_var; w, w_var; f_bounds, f_bounds_var;\r
1376                       d_bounds_list, d_bounds_list_var; dd_bounds_list, dd_bounds_list_var] o\r
1377                    inst_first_type_var ty) M_TAYLOR_INTERVAL' in\r
1378       let th1 = EQ_MP th0 m_taylor_th in\r
1379       let [domain_th; d2_th; lin_th; second_th] = CONJUNCTS th1 in\r
1380         domain_th, d2_th, lin_th, second_th;;\r
1381   \r
1382 \r
1383 \r
1384 (*\r
1385 eval_m_taylor_error n pp m_taylor_th;;\r
1386 (* 0.956 *)\r
1387 test 100 (eval_m_taylor_error n pp) m_taylor_th;;\r
1388 *)\r
1389 \r
1390 (**********************)\r
1391 \r
1392 (* bound *)\r
1393 let M_TAYLOR_BOUND' = \r
1394   prove(`m_taylor_interval f domain y (vector w:real^N) (f_lo, f_hi) d_bounds_list dd_bounds_list /\\r
1395           m_taylor_error f domain (vector w) error /\\r
1396           ITLIST2 (\a b c. b * iabs a + c) d_bounds_list w (&0) <= b /\\r
1397           b + inv(&2) * error <= a /\\r
1398           lo <= f_lo - a /\\r
1399           f_hi + a <= hi /\\r
1400           LENGTH w = dimindex (:N) /\ LENGTH d_bounds_list = dimindex (:N) ==>\r
1401               (!x. x IN interval [domain] ==> interval_arith (f x) (lo, hi))`,\r
1402    REWRITE_TAC[GSYM m_bounded_on_int; m_taylor_interval; m_lin_approx; ALL_N_EL] THEN\r
1403      set_tac "s" `ITLIST2 _1 _2 _3 _4` THEN STRIP_TAC THEN\r
1404      SUBGOAL_THEN `diff2_domain domain (f:real^N->real)` ASSUME_TAC THENL\r
1405      [\r
1406        UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN\r
1407          SIMP_TAC[diff2_domain; diff2c_domain; diff2c];\r
1408        ALL_TAC\r
1409      ] THEN\r
1410      apply_tac m_taylor_bounds THEN\r
1411      MAP_EVERY EXISTS_TAC [`y:real^N`; `vector w:real^N`; `error:real`; `f_lo:real`; `f_hi:real`; `a:real`] THEN\r
1412      ASM_REWRITE_TAC[] THEN\r
1413      MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `b + inv (&2) * error` THEN ASM_REWRITE_TAC[] THEN\r
1414      REWRITE_TAC[real_div; REAL_MUL_AC] THEN\r
1415      MATCH_MP_TAC REAL_LE_ADD2 THEN REWRITE_TAC[REAL_LE_REFL] THEN\r
1416      MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `s:real` THEN ASM_REWRITE_TAC[] THEN\r
1417      EXPAND_TAC "s" THEN new_rewrite [] [] ITLIST2_EQ_SUM THEN ASM_REWRITE_TAC[LE_REFL] THEN\r
1418      MATCH_MP_TAC SUM_LE THEN REWRITE_TAC[FINITE_NUMSEG] THEN REPEAT STRIP_TAC THEN\r
1419      MATCH_MP_TAC REAL_LE_MUL2 THEN ASM_SIMP_TAC[VECTOR_COMPONENT; REAL_LE_REFL; REAL_ABS_POS] THEN\r
1420      CONJ_TAC THENL\r
1421      [\r
1422        UNDISCH_TAC `m_cell_domain domain (y:real^N) (vector w)` THEN\r
1423          new_rewrite [] [] pair_eq THEN REWRITE_TAC[m_cell_domain] THEN\r
1424          DISCH_THEN (MP_TAC o SPEC `x:num`) THEN ASM_REWRITE_TAC[] THEN\r
1425          ASM_SIMP_TAC[VECTOR_COMPONENT] THEN REAL_ARITH_TAC;\r
1426        ALL_TAC\r
1427      ] THEN\r
1428      FIRST_X_ASSUM (MP_TAC o SPEC `x - 1`) THEN\r
1429      ANTS_TAC THENL\r
1430      [\r
1431        POP_ASSUM MP_TAC THEN ASM_REWRITE_TAC[IN_NUMSEG] THEN ARITH_TAC;\r
1432        ALL_TAC\r
1433      ] THEN\r
1434      DISCH_THEN (MP_TAC o MATCH_MP interval_arith_abs_le) THEN\r
1435      SUBGOAL_THEN `1 + x - 1 = x` (fun th -> REWRITE_TAC[th]) THENL\r
1436      [\r
1437        POP_ASSUM MP_TAC THEN REWRITE_TAC[IN_NUMSEG] THEN ARITH_TAC;\r
1438        ALL_TAC\r
1439      ] THEN\r
1440      DISCH_THEN MATCH_MP_TAC THEN REWRITE_TAC[REAL_LE_REFL]);;\r
1441 \r
1442 \r
1443 (* upper *)\r
1444 let M_TAYLOR_UPPER_BOUND' = prove(`m_taylor_interval f domain y (vector w) (f_lo, f_hi) \r
1445                                     d_bounds_list dd_bounds_list /\\r
1446                                     m_taylor_error f domain (vector w:real^N) error /\\r
1447                                     ITLIST2 (\a b c. b * iabs a + c) d_bounds_list w (&0) <= b /\\r
1448                                     b + inv(&2) * error <= a /\\r
1449                                     f_hi + a <= hi /\\r
1450                                     LENGTH w = dimindex (:N) /\ LENGTH d_bounds_list = dimindex (:N) ==>\r
1451                                         (!p. p IN interval [domain] ==> f p <= hi)`,\r
1452    STRIP_TAC THEN\r
1453      MP_TAC (INST[`f_lo - a:real`, `lo:real`] M_TAYLOR_BOUND') THEN\r
1454      ASM_SIMP_TAC[interval_arith; REAL_LE_REFL]);;\r
1455 \r
1456 \r
1457 (* lower *)\r
1458 let M_TAYLOR_LOWER_BOUND' = prove(`m_taylor_interval f domain y (vector w:real^N) (f_lo, f_hi) \r
1459                                     d_bounds_list dd_bounds_list /\\r
1460                                     m_taylor_error f domain (vector w) error /\\r
1461                                     ITLIST2 (\a b c. b * iabs a + c) d_bounds_list w (&0) <= b /\\r
1462                                     b + inv(&2) * error <= a /\\r
1463                                     lo <= f_lo - a /\\r
1464                                     LENGTH w = dimindex (:N) /\ LENGTH d_bounds_list = dimindex (:N) ==>\r
1465                                     (!p. p IN interval [domain] ==> lo <= f p)`,\r
1466    STRIP_TAC THEN\r
1467      MP_TAC (INST[`f_hi + a:real`, `hi:real`] M_TAYLOR_BOUND') THEN\r
1468      ASM_SIMP_TAC[interval_arith; REAL_LE_REFL]);;\r
1469 \r
1470 \r
1471 \r
1472 (* arrays *)\r
1473 let gen_taylor_bound_th bound_th n =\r
1474   let th0 = (DISCH_ALL o MY_RULE o REWRITE_RULE[MY_RULE M_TAYLOR_ERROR_ITLIST2_ALT]) bound_th in\r
1475   let ns = 1--n in\r
1476   let mk_list_hd l = mk_list (l, type_of (hd l)) in\r
1477   let w_list = mk_list_hd (map (fun i -> w_vars_array.(i)) ns) in\r
1478   let d_bounds_list = mk_list_hd (map (fun i -> df_vars_array.(i)) ns) in\r
1479   let dd_bounds_list = mk_list_hd (map (fun i -> mk_list_hd (map (fun j -> dd_vars_array.(i).(j)) (1--i))) ns) in\r
1480   let th1 = (INST[w_list, w_var_list; d_bounds_list, d_bounds_list_var; \r
1481                   dd_bounds_list, dd_bounds_list_var] o INST_TYPE[n_type_array.(n), nty]) th0 in\r
1482   let th2 = REWRITE_RULE[LAST; NOT_CONS_NIL; BUTLAST; all_n; ITLIST2_DEF; LENGTH; ARITH; dimindex_array.(n)] th1 in\r
1483   let th3 = REWRITE_RULE[HD; TL; REAL_MUL_RZERO; REAL_ADD_RID; GSYM error_mul_f2] th2 in\r
1484     (MY_RULE o REWRITE_RULE[float_inv2_th; SYM float2_eq]) th3;;\r
1485 \r
1486 \r
1487 let m_taylor_upper_array = Array.init (max_dim + 1)\r
1488   (fun i -> if i < 1 then TRUTH else gen_taylor_bound_th M_TAYLOR_UPPER_BOUND' i);;\r
1489 \r
1490 let m_taylor_lower_array = Array.init (max_dim + 1)\r
1491   (fun i -> if i < 1 then TRUTH else gen_taylor_bound_th M_TAYLOR_LOWER_BOUND' i);;\r
1492 \r
1493 let m_taylor_bound_array = Array.init (max_dim + 1)\r
1494   (fun i -> if i < 1 then TRUTH else gen_taylor_bound_th M_TAYLOR_BOUND' i);;\r
1495 \r
1496 (***************************)\r
1497 (* eval_m_taylor_bounds0 *)\r
1498 let b_var_real = `b:real`;;\r
1499 \r
1500 \r
1501 let eval_m_taylor_bounds0 mode n pp m_taylor_th =\r
1502   let bound_th = \r
1503     if mode = "upper" then \r
1504       m_taylor_upper_array.(n)\r
1505     else if mode = "bound" then\r
1506       m_taylor_bound_array.(n)\r
1507     else\r
1508       m_taylor_lower_array.(n) in\r
1509 \r
1510   let f_tm, domain_tm, y_tm, w_tm, f_bounds, d_bounds_list, dd_bounds_list = dest_m_taylor (concl m_taylor_th) in\r
1511   let f_lo, f_hi = dest_pair f_bounds and\r
1512       ws = dest_list (rand w_tm) and\r
1513       dfs = dest_list d_bounds_list and\r
1514       dds = map dest_list (dest_list dd_bounds_list) in\r
1515   let ns = 1--n in\r
1516   let df_vars = map (fun i -> df_vars_array.(i)) ns and\r
1517       dd_vars = map (fun i -> map (fun j -> dd_vars_array.(i).(j)) (1--i)) ns and\r
1518       w_vars = map (fun i -> w_vars_array.(i)) ns and\r
1519       y_var = mk_var ("y", type_of y_tm) and\r
1520       f_var = mk_var ("f", type_of f_tm) and\r
1521       domain_var = mk_var ("domain", type_of domain_tm) in\r
1522     \r
1523   (* sum of first partials *)\r
1524   let d_th =\r
1525     let mul_wd = map2 (error_mul_f2_le_conv2 pp) ws dfs in\r
1526       end_itlist (add_ineq_hi pp) mul_wd in\r
1527     \r
1528   let b_tm = (rand o concl) d_th in\r
1529     \r
1530   (* sum of second partials *)\r
1531   let dd_th =\r
1532     let ( * ), ( + ) = mul_ineq_pos_const_hi pp, add_ineq_hi pp in\r
1533     let mul_wdd = map2 (fun list i -> my_map2 (error_mul_f2_le_conv2 pp) ws list) dds ns in\r
1534     let sums1 = map (end_itlist ( + ) o butlast) (tl mul_wdd) in\r
1535     let sums2 = (hd o hd) mul_wdd :: map2 (fun list th1 -> last list + two_float * th1) (tl mul_wdd) sums1 in\r
1536     let sums = map2 ( * ) ws sums2 in\r
1537       end_itlist ( + ) sums in\r
1538 \r
1539   let error_tm = (rand o concl) dd_th in\r
1540 \r
1541   (* additional inequalities *)\r
1542   let ineq1_th = \r
1543     let ( * ), ( + ) = float_mul_hi pp, add_ineq_hi pp in\r
1544       mk_refl_ineq b_tm + float_inv2 * error_tm in\r
1545 \r
1546   let a_tm = (rand o concl) ineq1_th in\r
1547 \r
1548   let prove_ineq2, bounds_inst = \r
1549     if mode = "upper" then\r
1550       let ineq2 = float_add_hi pp f_hi a_tm in\r
1551         MY_PROVE_HYP ineq2, [(rand o concl) ineq2, hi_var_real]\r
1552     else if mode = "bound" then\r
1553       let ineq2 = float_add_hi pp f_hi a_tm in\r
1554       let ineq3 = float_sub_lo pp f_lo a_tm in\r
1555         MY_PROVE_HYP ineq2 o MY_PROVE_HYP ineq3, \r
1556     [(rand o concl) ineq2, hi_var_real; (lhand o concl) ineq3, lo_var_real]\r
1557     else\r
1558       let ineq2 = float_sub_lo pp f_lo a_tm in\r
1559         MY_PROVE_HYP ineq2, [(lhand o concl) ineq2, lo_var_real] in\r
1560 \r
1561   (* final step *)\r
1562   let inst_list =\r
1563     let inst1 = zip dfs df_vars in\r
1564     let inst2 = zip (List.flatten dds) (List.flatten dd_vars) in\r
1565     let inst3 = zip ws w_vars in\r
1566       inst1 @ inst2 @ inst3 in\r
1567     \r
1568     (MY_PROVE_HYP m_taylor_th o MY_PROVE_HYP dd_th o MY_PROVE_HYP d_th o \r
1569        MY_PROVE_HYP ineq1_th o prove_ineq2 o\r
1570        INST ([f_hi, f_hi_var; f_lo, f_lo_var; error_tm, error_var; \r
1571               a_tm, a_var_real; b_tm, b_var_real;\r
1572               y_tm, y_var; domain_tm, domain_var;\r
1573               f_tm, f_var] @ bounds_inst @ inst_list)) bound_th;;\r
1574 \r
1575 \r
1576 (* upper *)\r
1577 let eval_m_taylor_upper_bound = eval_m_taylor_bounds0 "upper";;\r
1578 \r
1579 (* lower *)\r
1580 let eval_m_taylor_lower_bound = eval_m_taylor_bounds0 "lower";;\r
1581 \r
1582 (* bound *)\r
1583 let eval_m_taylor_bound = eval_m_taylor_bounds0 "bound";;\r
1584 \r
1585 \r
1586 \r
1587 (*\r
1588 eval_m_taylor_upper_bound n pp m_taylor_th;;\r
1589 (* 10: 1.220 *)\r
1590 test 100 (eval_m_taylor_upper_bound n pp) m_taylor_th;;\r
1591 *)\r
1592 \r
1593 \r
1594 (******************************)\r
1595 (* taylor_upper_partial_bound *)\r
1596 (* taylor_lower_partial_bound *)\r
1597 \r
1598 \r
1599 (* bound *)\r
1600 let M_TAYLOR_PARTIAL_BOUND' = \r
1601   prove(`m_taylor_interval f domain (y:real^N) (vector w) f_bounds d_bounds_list dd_bounds_list /\\r
1602           i IN 1..dimindex (:N) /\\r
1603           EL (i - 1) d_bounds_list = (df_lo, df_hi) /\\r
1604           (!x. x IN interval [domain] ==> all_n 1 dd_list \r
1605              (\j int. interval_arith (if j <= i then partial2 j i f x else partial2 i j f x) int)) /\\r
1606           LENGTH dd_list = dimindex (:N) /\\r
1607           LENGTH d_bounds_list = dimindex (:N) /\ \r
1608           LENGTH w = dimindex (:N) /\\r
1609           ITLIST2 (\a b c. b * iabs a + c) dd_list w (&0) <= error /\\r
1610           df_hi + error <= hi ==>\r
1611           lo <= df_lo - error ==>\r
1612           (!x. x IN interval [domain] ==> interval_arith (partial i f x) (lo, hi))`,\r
1613    REWRITE_TAC[m_taylor_interval; m_lin_approx; ALL_N_EL; GSYM m_bounded_on_int] THEN\r
1614      set_tac "s" `ITLIST2 _1 _2 _3 _4` THEN REPEAT STRIP_TAC THEN\r
1615      SUBGOAL_THEN `1 <= i /\ i <= dimindex (:N)` (LABEL_TAC "i_in") THENL\r
1616      [\r
1617        ASM_REWRITE_TAC[GSYM IN_NUMSEG];\r
1618        ALL_TAC\r
1619      ] THEN\r
1620      SUBGOAL_THEN `diff2_domain domain (f:real^N->real)` ASSUME_TAC THENL\r
1621      [\r
1622        UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN\r
1623          SIMP_TAC[diff2_domain; diff2c_domain; diff2c];\r
1624        ALL_TAC\r
1625      ] THEN\r
1626      REWRITE_TAC[ETA_AX] THEN apply_tac m_taylor_partial_bounds THEN\r
1627      MAP_EVERY EXISTS_TAC [`y:real^N`; `vector w:real^N`; `error:real`; `df_lo:real`; `df_hi:real`] THEN\r
1628      ASM_REWRITE_TAC[] THEN\r
1629      FIRST_X_ASSUM (MP_TAC o SPEC `i - 1`) THEN\r
1630      ASM_SIMP_TAC[ARITH_RULE `1 <= i /\ i <= n ==> i - 1 < n`] THEN\r
1631      ASM_SIMP_TAC[ARITH_RULE `1 <= i ==> 1 + i - 1 = i`; interval_arith] THEN\r
1632      DISCH_THEN (fun th -> ALL_TAC) THEN\r
1633      REWRITE_TAC[m_taylor_partial_error] THEN REPEAT STRIP_TAC THEN\r
1634      MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `s:real` THEN ASM_REWRITE_TAC[] THEN\r
1635      EXPAND_TAC "s" THEN new_rewrite [] [] ITLIST2_EQ_SUM THEN ASM_REWRITE_TAC[LE_REFL] THEN\r
1636      MATCH_MP_TAC SUM_LE THEN REWRITE_TAC[FINITE_NUMSEG] THEN\r
1637      move ["j"; "j_in"] THEN\r
1638      ASM_SIMP_TAC[VECTOR_COMPONENT] THEN\r
1639      MATCH_MP_TAC REAL_LE_LMUL THEN CONJ_TAC THENL\r
1640      [\r
1641        UNDISCH_TAC `m_cell_domain domain (y:real^N) (vector w)` THEN\r
1642          new_rewrite [] [] pair_eq THEN REWRITE_TAC[m_cell_domain] THEN\r
1643          DISCH_THEN (MP_TAC o SPEC `j:num`) THEN ASM_REWRITE_TAC[] THEN\r
1644          ASM_SIMP_TAC[VECTOR_COMPONENT] THEN REAL_ARITH_TAC;\r
1645        ALL_TAC\r
1646      ] THEN\r
1647      FIRST_X_ASSUM (MP_TAC o SPEC `x:real^N`) THEN ASM_REWRITE_TAC[] THEN\r
1648      DISCH_THEN (MP_TAC o SPEC `j - 1`) THEN\r
1649      POP_ASSUM MP_TAC THEN REWRITE_TAC[IN_NUMSEG] THEN DISCH_TAC THEN\r
1650      ASM_SIMP_TAC[ARITH_RULE `1 <= j /\ j <= n ==> j - 1 < n`] THEN\r
1651      ASM_SIMP_TAC[ARITH_RULE `!i. 1 <= i ==> 1 + i - 1 = i`; GSYM partial2] THEN\r
1652      DISCH_THEN (MP_TAC o MATCH_MP interval_arith_abs_le) THEN\r
1653      COND_CASES_TAC THEN TRY (DISCH_THEN MATCH_MP_TAC THEN REWRITE_TAC[REAL_LE_REFL]) THEN\r
1654      new_rewrite [] [] mixed_second_partials THENL\r
1655      [\r
1656        UNDISCH_TAC `diff2c_domain domain (f:real^N->real)` THEN\r
1657          ASM_SIMP_TAC[diff2c_domain];\r
1658        ALL_TAC\r
1659      ] THEN\r
1660      DISCH_THEN MATCH_MP_TAC THEN REWRITE_TAC[REAL_LE_REFL]);;\r
1661 \r
1662 \r
1663 (* upper *)\r
1664 let M_TAYLOR_PARTIAL_UPPER' = \r
1665   prove(`m_taylor_interval f domain (y:real^N) (vector w) f_bounds d_bounds_list dd_bounds_list /\\r
1666           i IN 1..dimindex (:N) /\\r
1667           EL (i - 1) d_bounds_list = (df_lo, df_hi) /\\r
1668           (!x. x IN interval [domain] ==> all_n 1 dd_list \r
1669              (\j int. interval_arith (if j <= i then partial2 j i f x else partial2 i j f x) int)) /\\r
1670           LENGTH dd_list = dimindex (:N) /\\r
1671           LENGTH d_bounds_list = dimindex (:N) /\ \r
1672           LENGTH w = dimindex (:N) /\\r
1673           ITLIST2 (\a b c. b * iabs a + c) dd_list w (&0) <= error /\\r
1674           df_hi + error <= hi ==>\r
1675           (!x. x IN interval [domain] ==> partial i f x <= hi)`,\r
1676    REPEAT STRIP_TAC THEN\r
1677      MP_TAC (INST[`df_lo - error`, `lo:real`] M_TAYLOR_PARTIAL_BOUND') THEN ASM_REWRITE_TAC[REAL_LE_REFL] THEN\r
1678      ASM_SIMP_TAC[interval_arith]);;\r
1679 \r
1680 \r
1681 (* lower *)\r
1682 let M_TAYLOR_PARTIAL_LOWER' = \r
1683   prove(`m_taylor_interval f domain (y:real^N) (vector w) f_bounds d_bounds_list dd_bounds_list /\\r
1684           i IN 1..dimindex (:N) /\\r
1685           EL (i - 1) d_bounds_list = (df_lo, df_hi) /\\r
1686           (!x. x IN interval [domain] ==> all_n 1 dd_list \r
1687              (\j int. interval_arith (if j <= i then partial2 j i f x else partial2 i j f x) int)) /\\r
1688           LENGTH dd_list = dimindex (:N) /\\r
1689           LENGTH d_bounds_list = dimindex (:N) /\ \r
1690           LENGTH w = dimindex (:N) /\\r
1691           ITLIST2 (\a b c. b * iabs a + c) dd_list w (&0) <= error /\\r
1692           lo <= df_lo - error ==>\r
1693           (!x. x IN interval [domain] ==> lo <= partial i f x)`,\r
1694    REPEAT STRIP_TAC THEN\r
1695      MP_TAC (INST[`df_hi + error`, `hi:real`] M_TAYLOR_PARTIAL_BOUND') THEN ASM_REWRITE_TAC[REAL_LE_REFL] THEN\r
1696      ASM_SIMP_TAC[interval_arith]);;\r
1697 \r
1698 \r
1699 \r
1700 \r
1701 (* arrays *)\r
1702 let gen_taylor_partial_bound_th =\r
1703   let imp_and_eq = TAUT `((P ==> Q) /\ (P ==> R)) <=> (P ==> Q /\ R)` in\r
1704   let mk_list_hd l = mk_list (l, type_of (hd l)) in\r
1705   let dd_list_var = `dd_list : (real#real)list` in\r
1706     fun bound_th n i ->\r
1707       let ns = 1--n in\r
1708       let i_tm = mk_small_numeral i in\r
1709       let w_list = mk_list_hd (map (fun i -> w_vars_array.(i)) ns) in\r
1710       let d_bounds_list = mk_list_hd (map (fun i -> df_vars_array.(i)) ns) in\r
1711       let dd_bounds_list = mk_list_hd\r
1712         (map (fun i -> mk_list_hd (map (fun j -> dd_vars_array.(i).(j)) (1--i))) ns) in\r
1713       let dd_list = mk_list_hd \r
1714         (map (fun j -> if j <= i then dd_vars_array.(i).(j) else dd_vars_array.(j).(i)) ns) in\r
1715 \r
1716       let th1 = (INST[w_list, w_var_list; d_bounds_list, d_bounds_list_var; dd_list, dd_list_var; i_tm, `i:num`;\r
1717                       dd_bounds_list, dd_bounds_list_var] o INST_TYPE[n_type_array.(n), nty]) bound_th in\r
1718       let th2 = REWRITE_RULE[REAL_ADD_RID; HD; TL; ITLIST2_DEF; LENGTH; ARITH; dimindex_array.(n)] th1 in\r
1719       let th3 = REWRITE_RULE[IN_NUMSEG; ARITH; el_thms_array.(n).(i - 1)] th2 in\r
1720       let th4 = (REWRITE_RULE[] o INST[`df_lo:real, df_hi:real`, df_vars_array.(i)]) th3 in\r
1721       let th5 = (MY_RULE o REWRITE_RULE[GSYM imp_and_eq; GSYM AND_FORALL_THM; all_n; ARITH]) th4 in\r
1722 \r
1723       let m_taylor_hyp = find (can dest_m_taylor) (hyp th5) in\r
1724       let t_th0 = (REWRITE_RULE[ARITH; all_n; second_bounded; m_taylor_interval] o ASSUME) m_taylor_hyp in\r
1725       let t_th1 = REWRITE_RULE[GSYM imp_and_eq; GSYM AND_FORALL_THM] t_th0 in\r
1726         (MY_RULE_NUM o REWRITE_RULE[GSYM error_mul_f2; t_th1] o DISCH_ALL) th5;;\r
1727 \r
1728 \r
1729 (* The (n, i)-th element is the theorem |- i IN 1..dimindex (:n) *)\r
1730 let i_in_array = Array.init (max_dim + 1)\r
1731   (fun i -> Array.init (i + 1)\r
1732      (fun j ->\r
1733         if j < 1 then TRUTH else\r
1734           let j_tm = mk_small_numeral j in\r
1735           let tm0 = `j IN 1..dimindex (:N)` in\r
1736           let tm1 = (subst [j_tm, `j:num`] o inst [n_type_array.(i), nty]) tm0 in\r
1737             prove(tm1, REWRITE_TAC[dimindex_array.(i); IN_NUMSEG] THEN ARITH_TAC)));;\r
1738 \r
1739 let m_taylor_partial_upper_array, m_taylor_partial_lower_array, m_taylor_partial_bound_array =\r
1740   let gen_array bound_th = Array.init (max_dim + 1)\r
1741     (fun i -> Array.init (i + 1)\r
1742        (fun j -> if j < 1 then TRUTH else gen_taylor_partial_bound_th bound_th i j)) in\r
1743     gen_array M_TAYLOR_PARTIAL_UPPER', gen_array M_TAYLOR_PARTIAL_LOWER', gen_array M_TAYLOR_PARTIAL_BOUND';;\r
1744 \r
1745 \r
1746 (***************************)\r
1747 \r
1748 \r
1749 let eval_m_taylor_partial_bounds0 mode n pp i m_taylor_th =\r
1750   let bound_th = \r
1751     if mode = "upper" then \r
1752       m_taylor_partial_upper_array.(n).(i)\r
1753     else if mode = "bound" then\r
1754       m_taylor_partial_bound_array.(n).(i)\r
1755     else\r
1756       m_taylor_partial_lower_array.(n).(i) in\r
1757 \r
1758   let f_tm, domain_tm, y_tm, w_tm, f_bounds, d_bounds_list, dd_bounds_list = dest_m_taylor (concl m_taylor_th) in\r
1759   let ws = dest_list (rand w_tm) and\r
1760       dfs = dest_list d_bounds_list and\r
1761       dds = map dest_list (dest_list dd_bounds_list) in\r
1762   let ns = 1--n in\r
1763   let df_lo, df_hi = dest_pair (List.nth dfs (i - 1)) and\r
1764       df_vars = map (fun i -> df_vars_array.(i)) ns and\r
1765       dd_vars = map (fun i -> map (fun j -> dd_vars_array.(i).(j)) (1--i)) ns and\r
1766       w_vars = map (fun i -> w_vars_array.(i)) ns and\r
1767       y_var = mk_var ("y", type_of y_tm) and\r
1768       f_var = mk_var ("f", type_of f_tm) and\r
1769       domain_var = mk_var ("domain", type_of domain_tm) in\r
1770 \r
1771   (* sum of second partials *)\r
1772   let dd_list = map \r
1773     (fun j -> if j <= i then List.nth (List.nth dds (i-1)) (j-1) else List.nth (List.nth dds (j-1)) (i-1)) ns in\r
1774 \r
1775   let dd_th =\r
1776     let mul_dd = map2 (error_mul_f2_le_conv2 pp) ws dd_list in\r
1777       end_itlist (add_ineq_hi pp) mul_dd in\r
1778 \r
1779   let error_tm = (rand o concl) dd_th in\r
1780 \r
1781   (* additional inequalities *)\r
1782   let prove_ineq, bounds_inst = \r
1783     if mode = "upper" then\r
1784       let ineq2 = float_add_hi pp df_hi error_tm in\r
1785         MY_PROVE_HYP ineq2, [(rand o concl) ineq2, hi_var_real]\r
1786     else if mode = "bound" then\r
1787       let ineq2 = float_add_hi pp df_hi error_tm in\r
1788       let ineq3 = float_sub_lo pp df_lo error_tm in\r
1789         MY_PROVE_HYP ineq2 o MY_PROVE_HYP ineq3, \r
1790     [(rand o concl) ineq2, hi_var_real; (lhand o concl) ineq3, lo_var_real]\r
1791     else\r
1792       let ineq2 = float_sub_lo pp df_lo error_tm in\r
1793         MY_PROVE_HYP ineq2, [(lhand o concl) ineq2, lo_var_real] in\r
1794 \r
1795   (* final step *)\r
1796   let inst_list =\r
1797     let inst1 = zip dfs df_vars in\r
1798     let inst2 = zip (List.flatten dds) (List.flatten dd_vars) in\r
1799     let inst3 = zip ws w_vars in\r
1800       inst1 @ inst2 @ inst3 in\r
1801     \r
1802     (MY_PROVE_HYP m_taylor_th o MY_PROVE_HYP dd_th o prove_ineq o\r
1803        INST ([df_hi, df_hi_var; df_lo, df_lo_var; error_tm, error_var; \r
1804               y_tm, y_var; domain_tm, domain_var; f_bounds, f_bounds_var;\r
1805               f_tm, f_var] @ bounds_inst @ inst_list)) bound_th;;\r
1806 \r
1807 \r
1808 (* upper *)\r
1809 let eval_m_taylor_partial_upper = eval_m_taylor_partial_bounds0 "upper";;\r
1810 \r
1811 (* lower *)\r
1812 let eval_m_taylor_partial_lower = eval_m_taylor_partial_bounds0 "lower";;\r
1813 \r
1814 (* bound *)\r
1815 let eval_m_taylor_partial_bound = eval_m_taylor_partial_bounds0 "bound";;\r
1816 \r
1817 \r
1818 (*\r
1819 eval_m_taylor_partial_upper n pp 1 m_taylor_th;;\r
1820 eval_m_taylor_partial_upper n pp 4 m_taylor_th;;\r
1821 \r
1822 eval_m_taylor_partial_lower n pp 1 m_taylor_th;;\r
1823 eval_m_taylor_partial_lower n pp 4 m_taylor_th;;\r
1824 \r
1825 (* 10: 0.252 *)\r
1826 test 100 (eval_m_taylor_partial_upper n pp 1) m_taylor_th;;\r
1827 (* 10: 0.296 *)\r
1828 test 100 (eval_m_taylor_partial_upper n pp 4) m_taylor_th;;\r
1829 *)\r
1830 \r
1831 \r
1832 \r
1833 \r
1834 (**********************************)\r
1835 \r
1836 (*\r
1837 let pp = 5;;\r
1838 let eval_butcher = eval_m_taylor_poly0 pp butcher_poly;;\r
1839 \r
1840 let n = 6;;\r
1841 let xx = mk_vector_list (replicate one_float n) and\r
1842     zz = mk_vector_list (replicate two_float n);;\r
1843 \r
1844 let domain6_th = mk_m_center_domain n pp (rand xx) (rand zz);;\r
1845 \r
1846 let m_taylor_th = eval_butcher pp pp domain6_th;;\r
1847 \r
1848 \r
1849 \r
1850 eval_m_taylor_upper_bound n pp m_taylor_th;;\r
1851 eval_m_taylor_partial_upper n pp 1 m_taylor_th;;\r
1852 eval_m_taylor_partial_upper n pp 5 m_taylor_th;;\r
1853 \r
1854 \r
1855 (* 10: 2.524 *)\r
1856 (* 10 (mixed): 1.480 *)\r
1857 test 100 (eval_m_taylor_upper_bound n pp) m_taylor_th;;\r
1858 (* 10: 0.416 *)\r
1859 (* 10 (mixed): 0.348 *)\r
1860 test 100 (eval_m_taylor_partial_upper n pp 1) m_taylor_th;;\r
1861 (* 10: 0.432 *)\r
1862 (* 10 (mixed): 0.300 *)\r
1863 test 100 (eval_m_taylor_partial_upper n pp 5) m_taylor_th;;\r
1864 *)\r
1865 \r
1866 \r
1867 \r