needs "../formal_lp/formal_interval/lin_approx.hl";;

(* Explicit formula (float) for inv(&2) *)

let float_inv2_th =
  let one_float_eq_one = FLOAT_TO_NUM_CONV one_float in
  
let inv2_eq_lemma = 
prove(`interval_arith (&2 * x) (&1, &1) ==> inv (&2) = x`,

			    REWRITE_TAC[interval_arith] THEN CONV_TAC REAL_FIELD) in
  let half_tm = (fst o dest_pair o rand o concl) (float_interval_inv 1 two_interval) in
  let half_interval = mk_const_interval half_tm in
  let mul_th = REWRITE_RULE[one_float_eq_one] (float_interval_mul 2 two_interval half_interval) in
    MATCH_MP inv2_eq_lemma mul_th;;

let float_inv2 = rand (concl float_inv2_th);;
let inv2_interval = mk_const_interval float_inv2;;

let INTERVAL_IMP_HI = (RULE o prove)(`interval_arith x (lo, hi) ==> x <= hi`,
				     REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;

let interval_imp_hi int_th =
  let lhs, bounds = dest_comb (concl int_th) in
  let x_tm, (lo_tm, hi_tm) = rand lhs, dest_pair bounds in
  let th0 = INST[x_tm, x_var_real; lo_tm, lo_var_real; hi_tm, hi_var_real] INTERVAL_IMP_HI in
    MY_PROVE_HYP int_th th0;;


(*****************************)
(* cell_domain *)

let dest_cell_domain tm =
  let lhs, w_tm = dest_comb tm in
  let lhs2, z_tm = dest_comb lhs in
  let lhs3, y_tm = dest_comb lhs2 in
    rand lhs3, y_tm, z_tm, w_tm;;


(******************************************)

let CELL_DOMAIN_IMP_INTERVAL_ARITH = 
prove(`cell_domain x y z w ==> interval_arith y (x, z)`,

					   REWRITE_TAC[cell_domain; interval_arith] THEN REAL_ARITH_TAC);;


let MK_CELL_DOMAIN_0' = (UNDISCH_ALL o PURE_REWRITE_RULE[GSYM IMP_IMP] o prove)
  (`(x <= y <=> T) /\ (y <= z <=> T) /\ y - x <= w1 /\ z - y <= w2 /\ (w1 <= w <=> T) /\ (w2 <= w <=> T) ==> 
     cell_domain x y z w`,
   SIMP_TAC[cell_domain] THEN REAL_ARITH_TAC);;


let MK_CELL_DOMAIN' = (UNDISCH_ALL o PURE_REWRITE_RULE[GSYM IMP_IMP] o prove)
  (`(x <= y <=> T) /\ (y <= z <=> T) /\ y - x <= w1 /\ z - y <= w2 /\ max w1 w2 = w ==>
				       cell_domain x y z w`, REWRITE_TAC[cell_domain] THEN REAL_ARITH_TAC);;


let mk_cell_domain0 pp x_tm y_tm z_tm w_tm =
  let xy_th = float_le x_tm y_tm and
      yz_th = float_le y_tm z_tm in
  let yx_w1 = float_sub_hi pp y_tm x_tm and
      zy_w2 = float_sub_hi pp z_tm y_tm in
  let w1_tm = (rand o concl) yx_w1 and
      w2_tm = (rand o concl) zy_w2 in
  let w1w = float_le w1_tm w_tm and
      w2w = float_le w2_tm w_tm in
  let th0 = INST[x_tm, x_var_real; y_tm, y_var_real; z_tm, z_var_real; 
		 w_tm, w_var_real; w1_tm, w1_var_real; w2_tm, w2_var_real] MK_CELL_DOMAIN_0' in
    (MY_PROVE_HYP xy_th o MY_PROVE_HYP yz_th o MY_PROVE_HYP yx_w1 o MY_PROVE_HYP zy_w2 o
       MY_PROVE_HYP w1w o MY_PROVE_HYP w2w) th0;;


let mk_cell_domain pp x_tm y_tm z_tm =
  let w1_th = float_sub_hi pp y_tm x_tm and
      w2_th = float_sub_hi pp z_tm y_tm in
  let w1_tm = (rand o concl) w1_th and
      w2_tm = (rand o concl) w2_th in
  let max_th = float_max w1_tm w2_tm in
  let w_tm = (rand o concl) max_th in
  let x_le_y = float_le x_tm y_tm and
      y_le_z = float_le y_tm z_tm in
  let th0 = INST[x_tm, x_var_real; y_tm, y_var_real; z_tm, z_var_real;
		 w_tm, w_var_real; w1_tm, w1_var_real; w2_tm, w2_var_real] MK_CELL_DOMAIN' in
    (MY_PROVE_HYP x_le_y o MY_PROVE_HYP y_le_z o 
       MY_PROVE_HYP w1_th o MY_PROVE_HYP w2_th o MY_PROVE_HYP max_th) th0;;
  



(* TODO: remove formal computation of y *)
let mk_center_domain pp x_tm z_tm =
  let y0_tm = (rand o concl) (float_add_hi pp x_tm z_tm) in
  let y_tm = (rand o concl) (float_mul_eq float_inv2 y0_tm) in
    mk_cell_domain pp x_tm y_tm z_tm;;


(*
let x_tm = mk_float 1 50;;
let z_tm = mk_float 4 50;;

let domain_th = mk_center_domain 2 x_tm z_tm;;

let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th);;
mk_cell_domain 2 x_tm y_tm z_tm;;
mk_cell_domain0 2 x_tm y_tm z_tm w_tm;;
*)


(***************************************************)
(* Elementary functions and their taylor intervals *)
(***************************************************)


(*****************************)
(* f(x) = x *)

let TAYLOR_X' = RULE taylor_x;;


let eval_taylor_x domain_th =
  let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th) in
  let th0 = INST[x_tm, x_var_real; y_tm, y_var_real; z_tm, z_var_real; w_tm, w_var_real] TAYLOR_X' in
    MY_PROVE_HYP domain_th th0;;

 
(*			  
eval_taylor_x domain_th;;
*)

(********************************)
(* f(x) = const *)

let TAYLOR_CONST' = RULE taylor_const;;
let eval_taylor_const domain_th c_tm =
  let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th) in
    (MY_PROVE_HYP domain_th o
       INST[c_tm, c_var_real; x_tm, x_var_real; y_tm, y_var_real; 
	    z_tm, z_var_real; w_tm, w_var_real]) TAYLOR_CONST';;


let TAYLOR_CONST_INTERVAL' = (RULE o prove)
  (`cell_domain x y z w ==> interval_arith c bounds ==>
     taylor_interval (\x. c) x y z w bounds (&0, &0) (&0, &0)`,
  REPEAT DISCH_TAC THEN
    MP_TAC (SPEC_ALL taylor_const) THEN ASM_SIMP_TAC[taylor_interval; lin_approx_eq]);;

let eval_taylor_const_interval int_th domain_th =
  let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th) in
  let lhs, bounds = dest_comb (concl int_th) in
  let c_tm = rand lhs in
    (MY_PROVE_HYP int_th o MY_PROVE_HYP domain_th o
       INST[c_tm, c_var_real; bounds, bounds_var; x_tm, x_var_real; y_tm, y_var_real;
	    z_tm, z_var_real; w_tm, w_var_real]) TAYLOR_CONST_INTERVAL';;


(*
eval_taylor_const domain_th `&3`;;
*)

(***********************************)
(* f(x) = atn x *)

let taylor_atn = 
prove(`cell_domain x y z w ==>
			 bounded_on_int (\x. (-- &2 * x) / (&1 + x pow 2) pow 2) (x,z) dd_bounds ==>
			 lin_approx atn y f_bounds df_bounds ==>
			 taylor_interval atn x y z w f_bounds df_bounds dd_bounds`,

  REWRITE_TAC[cell_domain; taylor_interval; real_div; REAL_INV_POW] THEN REPEAT DISCH_TAC THEN 
    ASM_REWRITE_TAC[] THEN
    MATCH_MP_TAC second_derivative_atn_bounds THEN ASM_REWRITE_TAC[]);;

let TAYLOR_ATN' = (UNDISCH_ALL o REWRITE_RULE[REAL_POW_2]) taylor_atn;;


let eval_taylor_atn pp domain_th =
  let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th) in
  let atn_lin = eval_lin_approx_atn pp y_tm in
  let _, _, f_bounds, df_bounds = dest_lin_approx (concl atn_lin) in

  let ddf_th = 
    let x_th = ASSUME (mk_interval x_var_real (mk_pair (x_tm, z_tm))) and
	( * ) = float_interval_mul pp and
	(+) = float_interval_add pp and
	(/) = float_interval_div pp in
    let x2_1 = one_interval + x_th * x_th in
      mk_bounded_on_int ((neg_two_interval * x_th) / (x2_1 * x2_1)) in
 
  let dd_bounds = rand (concl ddf_th) in
    (MY_PROVE_HYP domain_th o MY_PROVE_HYP atn_lin o MY_PROVE_HYP ddf_th o
       INST[x_tm, x_var_real; y_tm, y_var_real; z_tm, z_var_real; w_tm, w_var_real;
	    dd_bounds, dd_bounds_var; f_bounds, f_bounds_var; df_bounds, df_bounds_var]) TAYLOR_ATN';;



(*
let domain_th = mk_center_domain 5 (mk_float 1 50) (mk_float 4 50);;
let pp = 5;;
eval_taylor_atn pp domain_th;;
*)



(***********************************)
(* f(x) = inv x *)

let taylor_inv = 
prove(`cell_domain x y z w ==>
			 interval_not_zero (x,z) ==>
			 bounded_on_int (\x. &2 * inv (x pow 3)) (x,z) dd_bounds ==>
			 lin_approx inv y f_bounds df_bounds ==>
			 taylor_interval inv x y z w f_bounds df_bounds dd_bounds`,

  REWRITE_TAC[cell_domain; taylor_interval] THEN REPEAT DISCH_TAC THEN 
    ASM_REWRITE_TAC[] THEN
    MATCH_MP_TAC (REWRITE_RULE[IMP_IMP] second_derivative_inv_bounds) THEN ASM_REWRITE_TAC[]);;

let TAYLOR_INV' = (UNDISCH_ALL o REWRITE_RULE[GSYM real_div; REAL_ARITH `x pow 3 = x * x * x`]) taylor_inv;;


let eval_taylor_inv pp domain_th =
  let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th) in
  let int_tm = mk_pair (x_tm, z_tm) in
  let lin_th = eval_lin_approx_inv pp y_tm in
  let _, _, f_bounds, df_bounds = dest_lin_approx (concl lin_th) in
  let cond_th = check_interval_not_zero int_tm in

  let ddf_th = 
    let x_th = ASSUME (mk_interval x_var_real (mk_pair (x_tm, z_tm))) and
	( * ) = float_interval_mul pp and
	(/) = float_interval_div pp in
      mk_bounded_on_int (two_interval / (x_th * (x_th * x_th))) in
 
  let dd_bounds = rand (concl ddf_th) in
    (MY_PROVE_HYP domain_th o MY_PROVE_HYP lin_th o MY_PROVE_HYP ddf_th o MY_PROVE_HYP cond_th o
       INST[x_tm, x_var_real; y_tm, y_var_real; z_tm, z_var_real; w_tm, w_var_real;
	    dd_bounds, dd_bounds_var; f_bounds, f_bounds_var; df_bounds, df_bounds_var]) TAYLOR_INV';;



(*
let domain_th = mk_center_domain 5 (mk_float 1 50) (mk_float 4 50);;
let pp = 5;;
eval_taylor_inv pp domain_th;;
*)


(***********************************)
(* f(x) = sqrt x *)

let taylor_sqrt = 
prove(`cell_domain x y z w ==>
			 interval_pos (x,z) ==>
			 bounded_on_int (\x. --inv (&4 * sqrt (x pow 3))) (x,z) dd_bounds ==>
			 lin_approx sqrt y f_bounds df_bounds ==>
			 taylor_interval sqrt x y z w f_bounds df_bounds dd_bounds`,

  REWRITE_TAC[cell_domain; taylor_interval] THEN REPEAT DISCH_TAC THEN 
    ASM_REWRITE_TAC[] THEN
    MATCH_MP_TAC (REWRITE_RULE[IMP_IMP] second_derivative_sqrt_bounds) THEN ASM_REWRITE_TAC[]);;

let TAYLOR_SQRT' = (UNDISCH_ALL o REWRITE_RULE[REAL_ARITH `x pow 3 = x * x * x`]) taylor_sqrt;;
let four_interval = mk_float_interval_small_num 4;;


let eval_taylor_sqrt pp domain_th =
  let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th) in
  let int_tm = mk_pair (x_tm, z_tm) in
  let lin_th = eval_lin_approx_sqrt pp y_tm in
  let _, _, f_bounds, df_bounds = dest_lin_approx (concl lin_th) in
  let cond_th = check_interval_pos int_tm in

  let ddf_th = 
    let x_th = ASSUME (mk_interval x_var_real (mk_pair (x_tm, z_tm))) and
	( * ) = float_interval_mul pp and
	sqrt = float_interval_sqrt pp and
	inv = float_interval_inv pp and
	neg = float_interval_neg in
      mk_bounded_on_int (neg (inv (four_interval * sqrt (x_th * (x_th * x_th))))) in
 
  let dd_bounds = rand (concl ddf_th) in
    (MY_PROVE_HYP domain_th o MY_PROVE_HYP lin_th o MY_PROVE_HYP ddf_th o MY_PROVE_HYP cond_th o
       INST[x_tm, x_var_real; y_tm, y_var_real; z_tm, z_var_real; w_tm, w_var_real;
	    dd_bounds, dd_bounds_var; f_bounds, f_bounds_var; df_bounds, df_bounds_var]) TAYLOR_SQRT';;



(*
let domain_th = mk_center_domain 5 (mk_float 1 50) (mk_float 4 50);;
let pp = 5;;
eval_taylor_sqrt pp domain_th;;
*)


(***********************************)
(* f(x) = acs x *)

let iabs_lt1_lemma = 
  (GEN_REWRITE_RULE (RAND_CONV o LAND_CONV o RAND_CONV) [GSYM (FLOAT_TO_NUM_CONV one_float)] o prove)
    (`iabs int < &1 <=> (iabs int < &1 <=> T)`, REWRITE_TAC[]);;

let taylor_acs = 
prove(`cell_domain x y z w ==>
			 iabs (x,z) < &1 ==>
			 bounded_on_int (\x. --(x / sqrt ((&1 - x * x) pow 3))) (x,z) dd_bounds ==>
			 lin_approx acs y f_bounds df_bounds ==>
			 taylor_interval acs x y z w f_bounds df_bounds dd_bounds`,

  REWRITE_TAC[cell_domain; taylor_interval] THEN REPEAT DISCH_TAC THEN 
    ASM_REWRITE_TAC[] THEN
    MATCH_MP_TAC (REWRITE_RULE[IMP_IMP] second_derivative_acs_bounds) THEN ASM_REWRITE_TAC[]);;

let TAYLOR_ACS' = 
  (UNDISCH_ALL o 
     PURE_ONCE_REWRITE_RULE[iabs_lt1_lemma] o
     REWRITE_RULE[REAL_ARITH `x pow 3 = x * x * x`]) taylor_acs;;


let eval_taylor_acs pp domain_th =
  let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th) in
  let int_tm = mk_pair (x_tm, z_tm) in
  let lin_th = eval_lin_approx_acs pp y_tm in
  let _, _, f_bounds, df_bounds = dest_lin_approx (concl lin_th) in
  let cond_th = check_interval_iabs int_tm one_float in

  let ddf_th = 
    let x_th = ASSUME (mk_interval x_var_real (mk_pair (x_tm, z_tm))) and
	( * ) = float_interval_mul pp and
	(-) = float_interval_sub pp and
	(/) = float_interval_div pp and
	sqrt = float_interval_sqrt pp and
	neg = float_interval_neg in
    let x2_1 = one_interval - x_th * x_th in
      mk_bounded_on_int (neg (x_th / sqrt (x2_1 * (x2_1 * x2_1)))) in
 
  let dd_bounds = rand (concl ddf_th) in
    (MY_PROVE_HYP domain_th o MY_PROVE_HYP lin_th o MY_PROVE_HYP ddf_th o MY_PROVE_HYP cond_th o
       INST[x_tm, x_var_real; y_tm, y_var_real; z_tm, z_var_real; w_tm, w_var_real;
	    dd_bounds, dd_bounds_var; f_bounds, f_bounds_var; df_bounds, df_bounds_var]) TAYLOR_ACS';;



(*
let domain_th1 = mk_center_domain 5 (mk_float 1 49) (mk_float 8 49);;
let domain_th2 = mk_center_domain 5 (mk_float (-5) 49) (mk_float 1 49);;
let pp = 5;;
eval_taylor_acs pp domain_th1;;
eval_taylor_atn pp domain_th1;;
eval_taylor_inv pp domain_th1;;
eval_taylor_sqrt pp domain_th1;;

eval_taylor_acs pp domain_th2;;
eval_taylor_atn pp domain_th2;;
eval_taylor_inv pp domain_th2;;
eval_taylor_sqrt pp domain_th2;;
*)


(***********************************************************)


(***********************************)
(* Taylor bounds *)


let dest_taylor_interval tm =
  let lhs, dd_bounds = dest_comb tm in
  let lhs, df_bounds = dest_comb lhs in
  let lhs, f_bounds = dest_comb lhs in
  let lhs, w_tm = dest_comb lhs in
  let lhs, z_tm = dest_comb lhs in
  let lhs, y_tm = dest_comb lhs in
  let lhs, x_tm = dest_comb lhs in
    rand lhs, x_tm, y_tm, z_tm, w_tm, f_bounds, df_bounds, dd_bounds;;


(****************************)
(* f_bounds *)

let TAYLOR_F_BOUNDS' = (RULE o REWRITE_RULE[float_inv2_th]) taylor_f_bounds;;


let eval_taylor_f_bounds pp th =
  let f_tm, x_tm, y_tm, z_tm, w_tm, f_bounds, df_bounds, dd_bounds = dest_taylor_interval (concl th) in
  let f_lo, f_hi = dest_pair f_bounds in
  let iabs_df_th = float_iabs df_bounds and
      iabs_dd_th = float_iabs dd_bounds in
  let df_tm = (rand o concl) iabs_df_th and
      dd_tm = (rand o concl) iabs_dd_th in

  let error_th =
    let df = mk_refl_ineq df_tm and
	inv2 = mk_refl_ineq float_inv2 and
	( * ) = mul_ineq_pos_const_hi pp and
	( + ) = add_ineq_hi pp in
      w_tm * (df + w_tm * (dd_tm * inv2)) in

  let t_tm = (rand o concl) error_th in
  let lo_th = float_sub_lo pp f_lo t_tm and
      hi_th = float_add_hi pp f_hi t_tm in
  let lo_tm = (lhand o concl) lo_th and
      hi_tm = (rand o concl) hi_th in

    (MY_PROVE_HYP lo_th o MY_PROVE_HYP hi_th o MY_PROVE_HYP error_th o 
       MY_PROVE_HYP iabs_df_th o MY_PROVE_HYP iabs_dd_th o MY_PROVE_HYP th o
       INST[dd_bounds, dd_bounds_var; dd_tm, dd_var_real; df_bounds, df_bounds_var;
	    df_tm, df_var_real; t_tm, t_var_real;
	    lo_tm, lo_var_real; f_lo, f_lo_var;
	    w_tm, w_var_real; f_hi, f_hi_var; hi_tm, hi_var_real;
	    f_tm, f_var_fun; x_tm, x_var_real; 
	    y_tm, y_var_real; z_tm, z_var_real]) TAYLOR_F_BOUNDS';;


(*
let pp = 3;;
let domain_th = mk_center_domain pp one_float two_float;;
let th = eval_taylor_atn pp domain_th;;
eval_taylor_f_bounds pp th;;
(* 10: 1.620 *)
test 1000 (eval_taylor_f_bounds pp) th;;
*)


(****************************)
(* df_bounds *)


let TAYLOR_DF_BOUNDS' = RULE taylor_df_bounds;;

let eval_taylor_df_bounds pp th =
  let f_tm, x_tm, y_tm, z_tm, w_tm, f_bounds, df_bounds, dd_bounds = dest_taylor_interval (concl th) in
  let df_lo, df_hi = dest_pair df_bounds in
  let iabs_dd_th = float_iabs dd_bounds in
  let dd_tm = (rand o concl) iabs_dd_th in
  let w_dd_hi = float_mul_hi pp w_tm dd_tm in

  let lo_th = sub_ineq_lo pp (mk_refl_ineq df_lo) w_dd_hi and
      hi_th = add_ineq_hi pp (mk_refl_ineq df_hi) w_dd_hi in
  let lo_tm = (lhand o concl) lo_th and
      hi_tm = (rand o concl) hi_th in

    (MY_PROVE_HYP lo_th o MY_PROVE_HYP hi_th o MY_PROVE_HYP th o MY_PROVE_HYP iabs_dd_th o
       INST[dd_bounds, dd_bounds_var; dd_tm, dd_var_real; f_bounds, f_bounds_var;
	    lo_tm, lo_var_real; df_lo, df_lo_var;
	    w_tm, w_var_real; df_hi, df_hi_var; hi_tm, hi_var_real;
	    f_tm, f_var_fun; x_tm, x_var_real; 
	    y_tm, y_var_real; z_tm, z_var_real]) TAYLOR_DF_BOUNDS';;


(*
let pp = 3;;
let th = eval_taylor_atn pp domain_th;;
eval_taylor_df_bounds pp th;;
*)

(**************************************************)
(* mk_taylor_interval_th, dest_taylor_interval_th *)

let TAYLOR_INTERVAL_EQ' = (RULE o prove)(`taylor_interval f x y z w f_bounds df_bounds dd_bounds <=>
					     cell_domain x y z w /\ lin_approx f y f_bounds df_bounds /\
					     has_bounded_second_derivative f (x,z) dd_bounds`,
					   SIMP_TAC[cell_domain; taylor_interval; CONJ_ACI]);;

let MK_TAYLOR_INTERVAL' = (RULE o prove)(`cell_domain x y z w /\ 
					  lin_approx f y f_bounds df_bounds /\
					  has_bounded_second_derivative f (x,z) dd_bounds ==>
					  taylor_interval f x y z w f_bounds df_bounds dd_bounds`,
					SIMP_TAC[cell_domain; taylor_interval]);;


let dest_taylor_interval_th th =
  let f_tm, x_tm, y_tm, z_tm, w_tm, f_bounds, df_bounds, dd_bounds = dest_taylor_interval (concl th) in
  let th0 = INST[f_tm, f_var_fun; x_tm, x_var_real; y_tm, y_var_real; z_tm, z_var_real;
		 w_tm, w_var_real; f_bounds, f_bounds_var;
		 df_bounds, df_bounds_var; dd_bounds, dd_bounds_var] TAYLOR_INTERVAL_EQ' in
  let [domain; approx; second] = CONJUNCTS (EQ_MP th0 th) in
    domain, approx, second;;



let mk_taylor_interval_th domain_th approx_th second_th =
  let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th) in
  let f_tm, _, f_bounds, df_bounds = dest_lin_approx (concl approx_th) in
  let dd_bounds = rand (concl second_th) in
  let th0 = INST[f_tm, f_var_fun; x_tm, x_var_real; y_tm, y_var_real; z_tm, z_var_real;
		 w_tm, w_var_real; f_bounds, f_bounds_var;
		 df_bounds, df_bounds_var; dd_bounds, dd_bounds_var] MK_TAYLOR_INTERVAL' in
    (MY_PROVE_HYP approx_th o MY_PROVE_HYP domain_th o MY_PROVE_HYP second_th) th0;;
  


(*
let th = eval_taylor_atn 5 domain_th;;
let domain_th, approx_th, second_th = dest_taylor_interval_th th;;
mk_taylor_interval_th domain_th approx_th second_th;;
*)


(**************************************)
(* has_bounded_second_derivative_dest *)

let HAS_BOUNDED_SECOND_DERIVATIVE_EQ' = SPEC_ALL has_bounded_second_derivative;;


let dest_has_bounded_second_derivative tm =
  let lhs, dd_bounds = dest_comb tm in
  let lhs2, int_tm = dest_comb lhs in
    rand lhs2, int_tm, dd_bounds;;


let dest_has_bounded_second_derivative_th th =
  let f_tm, int_tm, dd_bounds = dest_has_bounded_second_derivative (concl th) in
  let th0 = INST[f_tm, f_var_fun; int_tm, int_var; 
		 dd_bounds, dd_bounds_var] HAS_BOUNDED_SECOND_DERIVATIVE_EQ' in
  let [diff_th; bounded_th] = CONJUNCTS (EQ_MP th0 th) in
    diff_th, bounded_th;;


(*
let diff_th, bounded_th = dest_has_bounded_second_derivative_th second_th;;
*)

(********************************************)


(******************************)
(* Taylor interval arithmetic *)
(******************************)


(**************************************)
let taylor_interval_tm = `taylor_interval`;;

let taylor_interval_norm th eq_th =
  let lhs1, d2f = dest_comb (concl th) in
  let lhs2, d1f = dest_comb lhs1 in
  let lhs3, d0f = dest_comb lhs2 in
  let lhs4, w = dest_comb lhs3 in
  let lhs5, z = dest_comb lhs4 in
  let lhs6, y = dest_comb lhs5 in
  let x = rand lhs6 in
  let th0 = AP_TERM taylor_interval_tm eq_th in
  let th1 = (AP_THM (AP_THM (AP_THM (AP_THM (AP_THM (AP_THM (AP_THM th0 x) y) z) w) d0f) d1f) d2f) in
    EQ_MP th1 th;;


(*************************************)
(* taylor_interval_neg *)
let TAYLOR_INTERVAL_NEG' = (UNDISCH_ALL o prove)
  (`cell_domain x y z w ==>
     lin_approx (\x. -- (f x)) y bounds d_bounds ==>
     nth_diff_strong_int 2 (x,z) f ==>
     bounded_on_int (\x. -- nth_derivative 2 f x) (x,z) dd_bounds ==>
     taylor_interval (\x. -- (f x)) x y z w bounds d_bounds dd_bounds`,
   SIMP_TAC[taylor_interval; cell_domain] THEN ONCE_REWRITE_TAC[REAL_ARITH `--a = (-- (&1)) * a`] THEN
     REPEAT STRIP_TAC THEN
     ASM_SIMP_TAC[second_derivative_scale_bounds]);;


let taylor_interval_neg th1 =
  let domain_th, approx1_th, second1_th = dest_taylor_interval_th th1 in
  let diff1_th, dd1_th = dest_has_bounded_second_derivative_th second1_th in

  let approx_th = lin_approx_neg approx1_th in
  let dd_th = bounded_on_int_neg dd1_th in

  let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th) and
      f_tm = (rand o concl) diff1_th and
      _, _, bounds, d_bounds = dest_lin_approx (concl approx_th) and
      dd_bounds = (rand o concl) dd_th in

  let th0 = (MY_PROVE_HYP domain_th o MY_PROVE_HYP approx_th o 
	       MY_PROVE_HYP diff1_th o MY_PROVE_HYP dd_th o
	       INST[x_tm, x_var_real; y_tm, y_var_real; z_tm, z_var_real; w_tm, w_var_real;
		    f_tm, f_var_fun;
		    bounds, bounds_var; d_bounds, d_bounds_var; dd_bounds, dd_bounds_var])
    TAYLOR_INTERVAL_NEG' in
  let eq_th = unary_beta_eq f_tm neg_op_real in
    taylor_interval_norm th0 eq_th;;
     

(*
let th1 = eval_taylor_x domain_th and
    th2 = eval_taylor_sqrt pp domain_th;;
taylor_interval_neg th1;;
taylor_interval_neg th2;;
*)



(*************************************)
(* taylor_interval_add *)

let TAYLOR_INTERVAL_ADD' = (UNDISCH_ALL o prove)
  (`cell_domain x y z w ==>
     lin_approx (\x. f x + g x) y bounds d_bounds ==>
     nth_diff_strong_int 2 (x,z) f ==>
     nth_diff_strong_int 2 (x,z) g ==>
     bounded_on_int (\x. nth_derivative 2 f x + nth_derivative 2 g x) (x,z) dd_bounds ==>
     taylor_interval (\x. f x + g x) x y z w bounds d_bounds dd_bounds`,
   SIMP_TAC[taylor_interval; cell_domain] THEN REPEAT STRIP_TAC THEN
     ASM_SIMP_TAC[second_derivative_add_bounds]);;


let taylor_interval_add pp th1 th2 =
  let domain_th, approx1_th, second1_th = dest_taylor_interval_th th1 and
      _, approx2_th, second2_th = dest_taylor_interval_th th2 in
  let diff1_th, dd1_th = dest_has_bounded_second_derivative_th second1_th and
      diff2_th, dd2_th = dest_has_bounded_second_derivative_th second2_th in

  let approx_th = lin_approx_add pp approx1_th approx2_th in
  let dd_th = 
    let (+) = bounded_on_int_add pp in
      dd1_th + dd2_th in

  let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th) and
      f_tm = (rand o concl) diff1_th and
      g_tm = (rand o concl) diff2_th and
      _, _, bounds, d_bounds = dest_lin_approx (concl approx_th) and
      dd_bounds = (rand o concl) dd_th in

  let th0 = (MY_PROVE_HYP domain_th o MY_PROVE_HYP approx_th o 
	       MY_PROVE_HYP diff1_th o MY_PROVE_HYP diff2_th o MY_PROVE_HYP dd_th o
	       INST[x_tm, x_var_real; y_tm, y_var_real; z_tm, z_var_real; w_tm, w_var_real;
		    f_tm, f_var_fun; g_tm, g_var_fun;
		    bounds, bounds_var; d_bounds, d_bounds_var; dd_bounds, dd_bounds_var])
    TAYLOR_INTERVAL_ADD' in
  let eq_th = binary_beta_eq f_tm g_tm add_op_real in
    taylor_interval_norm th0 eq_th;;



(*
let th1 = eval_taylor_x domain_th and
    th2 = eval_taylor_sqrt pp domain_th;;
taylor_interval_add pp th1 th2;;
*)


(*************************************)
(* taylor_interval_sub *)

let TAYLOR_INTERVAL_SUB' = (UNDISCH_ALL o prove)
  (`cell_domain x y z w ==>
     lin_approx (\x. f x - g x) y bounds d_bounds ==>
     nth_diff_strong_int 2 (x,z) f ==>
     nth_diff_strong_int 2 (x,z) g ==>
     bounded_on_int (\x. nth_derivative 2 f x - nth_derivative 2 g x) (x,z) dd_bounds ==>
     taylor_interval (\x. f x - g x) x y z w bounds d_bounds dd_bounds`,
   SIMP_TAC[taylor_interval; cell_domain] THEN REPEAT STRIP_TAC THEN
     ASM_SIMP_TAC[second_derivative_sub_bounds]);;


let taylor_interval_sub pp th1 th2 =
  let domain_th, approx1_th, second1_th = dest_taylor_interval_th th1 and
      _, approx2_th, second2_th = dest_taylor_interval_th th2 in
  let diff1_th, dd1_th = dest_has_bounded_second_derivative_th second1_th and
      diff2_th, dd2_th = dest_has_bounded_second_derivative_th second2_th in

  let approx_th = lin_approx_sub pp approx1_th approx2_th in
  let dd_th = 
    let (-) = bounded_on_int_sub pp in
      dd1_th - dd2_th in

  let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th) and
      f_tm = (rand o concl) diff1_th and
      g_tm = (rand o concl) diff2_th and
      _, _, bounds, d_bounds = dest_lin_approx (concl approx_th) and
      dd_bounds = (rand o concl) dd_th in

  let th0 = (MY_PROVE_HYP domain_th o MY_PROVE_HYP approx_th o 
	       MY_PROVE_HYP diff1_th o MY_PROVE_HYP diff2_th o MY_PROVE_HYP dd_th o
	       INST[x_tm, x_var_real; y_tm, y_var_real; z_tm, z_var_real; w_tm, w_var_real;
		    f_tm, f_var_fun; g_tm, g_var_fun;
		    bounds, bounds_var; d_bounds, d_bounds_var; dd_bounds, dd_bounds_var])
    TAYLOR_INTERVAL_SUB' in
  let eq_th = binary_beta_eq f_tm g_tm sub_op_real in
    taylor_interval_norm th0 eq_th;;



(*
let th1 = eval_taylor_x domain_th and
    th2 = eval_taylor_sqrt pp domain_th;;
taylor_interval_sub pp th1 th2;;
*)


(*************************************)
(* taylor_interval_mul *)

let TAYLOR_INTERVAL_MUL' = (UNDISCH_ALL o prove)
  (`cell_domain x y z w ==>
     lin_approx (\x. f x * g x) y bounds d_bounds ==>
     nth_diff_strong_int 2 (x,z) f ==>
     nth_diff_strong_int 2 (x,z) g ==>
     bounded_on_int (\x. f x * nth_derivative 2 g x + &2 * derivative f x * derivative g x +
		     nth_derivative 2 f x * g x) (x,z) dd_bounds ==>
     taylor_interval (\x. f x * g x) x y z w bounds d_bounds dd_bounds`,
   SIMP_TAC[taylor_interval; cell_domain] THEN REPEAT STRIP_TAC THEN
     ASM_SIMP_TAC[second_derivative_mul_bounds]);;


let taylor_interval_mul pp th1 th2 =
  let domain_th, approx1_th, second1_th = dest_taylor_interval_th th1 and
      _, approx2_th, second2_th = dest_taylor_interval_th th2 in
  let diff1_th, dd1_th = dest_has_bounded_second_derivative_th second1_th and
      diff2_th, dd2_th = dest_has_bounded_second_derivative_th second2_th in

  let approx_th = lin_approx_mul pp approx1_th approx2_th in
  let dd_th =
    let dest = dest_bounded_on_int_raw in
    let f = (dest o eval_taylor_f_bounds pp) th1 and
	g = (dest o eval_taylor_f_bounds pp) th2 and
	df = (dest o eval_taylor_df_bounds pp) th1 and
	dg = (dest o eval_taylor_df_bounds pp) th2 and
	ddf = dest dd1_th and
	ddg = dest dd2_th and
	(+) = float_interval_add pp and
	( * ) = float_interval_mul pp in
      mk_bounded_on_int (f * ddg + (two_interval * (df * dg) + ddf * g)) in

  let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th) and
      f_tm = (rand o concl) diff1_th and
      g_tm = (rand o concl) diff2_th and
      _, _, bounds, d_bounds = dest_lin_approx (concl approx_th) and
      dd_bounds = (rand o concl) dd_th in

  let th0 = (MY_PROVE_HYP domain_th o MY_PROVE_HYP approx_th o 
	       MY_PROVE_HYP diff1_th o MY_PROVE_HYP diff2_th o MY_PROVE_HYP dd_th o
	       INST[x_tm, x_var_real; y_tm, y_var_real; z_tm, z_var_real; w_tm, w_var_real;
		    f_tm, f_var_fun; g_tm, g_var_fun;
		    bounds, bounds_var; d_bounds, d_bounds_var; dd_bounds, dd_bounds_var])
    TAYLOR_INTERVAL_MUL' in
  let eq_th = binary_beta_eq f_tm g_tm mul_op_real in
    taylor_interval_norm th0 eq_th;;



(*
let th1 = eval_taylor_x domain_th and
    th2 = eval_taylor_sqrt pp domain_th;;
let th3 = taylor_interval_sub pp th1 th2;;
taylor_interval_mul pp th1 th3;;
*)



(*************************************)
(* taylor_interval_div *)

let TAYLOR_INTERVAL_DIV' = (UNDISCH_ALL o prove)
 (`cell_domain x y z w ==>
    bounded_on_int g (x,z) g_bounds ==>
    interval_not_zero g_bounds ==>
    lin_approx (\x. f x / g x) y bounds d_bounds ==>
    nth_diff_strong_int 2 (x,z) f ==>
    nth_diff_strong_int 2 (x,z) g ==>
    bounded_on_int (\x. ((nth_derivative 2 f x * g x - f x * nth_derivative 2 g x) * g x -
			 &2 * derivative g x * (derivative f x * g x - f x * derivative g x)) / 
		      (g x * g x * g x)) (x,z) dd_bounds ==>
    taylor_interval (\x. f x / g x) x y z w bounds d_bounds dd_bounds`,
  SIMP_TAC[taylor_interval; cell_domain] THEN REPEAT STRIP_TAC THEN
    MATCH_MP_TAC ((SPEC_ALL o REWRITE_RULE[IMP_IMP]) second_derivative_div_bounds) THEN
    ASM_REWRITE_TAC[REAL_ARITH `g x pow 3 = g x * g x * (g:real->real) x`]);;


let taylor_interval_div pp th1 th2 =
  let domain_th, approx1_th, second1_th = dest_taylor_interval_th th1 and
      _, approx2_th, second2_th = dest_taylor_interval_th th2 in
  let diff1_th, dd1_th = dest_has_bounded_second_derivative_th second1_th and
      diff2_th, dd2_th = dest_has_bounded_second_derivative_th second2_th in

  let g_bounds_th = eval_taylor_f_bounds pp th2 and
      approx_th = lin_approx_div pp approx1_th approx2_th in
  let g_bounds = (rand o concl) g_bounds_th in
  let cond_th = check_interval_not_zero g_bounds in

  let dd_th =
    let dest = dest_bounded_on_int_raw in
    let f = (dest o eval_taylor_f_bounds pp) th1 and
	g = dest g_bounds_th and
	df = (dest o eval_taylor_df_bounds pp) th1 and
	dg = (dest o eval_taylor_df_bounds pp) th2 and
	ddf = dest dd1_th and
	ddg = dest dd2_th and
	(-) = float_interval_sub pp and
	( * ) = float_interval_mul pp and
	(/) = float_interval_div pp in
      mk_bounded_on_int (((ddf * g - f * ddg) * g - two_interval * (dg * (df * g - f * dg))) / (g * (g * g))) in

  let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th) and
      f_tm = (rand o concl) diff1_th and
      g_tm = (rand o concl) diff2_th and
      _, _, bounds, d_bounds = dest_lin_approx (concl approx_th) and
      dd_bounds = (rand o concl) dd_th in

  let th0 = (MY_PROVE_HYP domain_th o MY_PROVE_HYP approx_th o 
	       MY_PROVE_HYP g_bounds_th o MY_PROVE_HYP cond_th o
	       MY_PROVE_HYP diff1_th o MY_PROVE_HYP diff2_th o MY_PROVE_HYP dd_th o
	       INST[x_tm, x_var_real; y_tm, y_var_real; z_tm, z_var_real; w_tm, w_var_real;
		    f_tm, f_var_fun; g_tm, g_var_fun; g_bounds, g_bounds_var;
		    bounds, bounds_var; d_bounds, d_bounds_var; dd_bounds, dd_bounds_var])
    TAYLOR_INTERVAL_DIV' in
  let eq_th = binary_beta_eq f_tm g_tm div_op_real in
    taylor_interval_norm th0 eq_th;;


(*
let pp = 5;;
let domain_th = mk_center_domain pp one_float two_float;;
let th1 = eval_taylor_x domain_th and
    th2 = eval_taylor_sqrt pp domain_th;;
let th3 = taylor_interval_sub pp th1 th2;;
taylor_interval_div pp th1 th2;;
taylor_interval_div pp th3 th2;;
(* 10: 3.372 *)
test 100 (taylor_interval_div pp th1) th2;;
*)


(*************************************)
(* taylor_interval_atn *)

let TAYLOR_INTERVAL_ATN' = (UNDISCH_ALL o REWRITE_RULE[REAL_POW_2] o prove)
 (`cell_domain x y z w ==>
    lin_approx (\x. atn (f x)) y bounds d_bounds ==>
    nth_diff_strong_int 2 (x,z) f ==>
    bounded_on_int (\x. (nth_derivative 2 f x * (&1 + f x * f x) - &2 * f x * derivative f x pow 2) /
		      (&1 + f x * f x) pow 2) (x,z) dd_bounds ==>
    taylor_interval (\x. atn (f x)) x y z w bounds d_bounds dd_bounds`,
  SIMP_TAC[taylor_interval; cell_domain] THEN REPEAT STRIP_TAC THEN
    MATCH_MP_TAC ((SPEC_ALL o REWRITE_RULE[IMP_IMP]) second_derivative_compose_atn_bounds) THEN
    ASM_REWRITE_TAC[]);;


let taylor_interval_atn pp th1 =
  let domain_th, approx1_th, second1_th = dest_taylor_interval_th th1 in
  let diff1_th, dd1_th = dest_has_bounded_second_derivative_th second1_th in

  let approx_th = lin_approx_atn pp approx1_th in
  let dd_th =
    let dest = dest_bounded_on_int_raw in
    let f = (dest o eval_taylor_f_bounds pp) th1 and
	df = (dest o eval_taylor_df_bounds pp) th1 and
	ddf = dest dd1_th and
	(+) = float_interval_add pp and
	(-) = float_interval_sub pp and
	( * ) = float_interval_mul pp and
	(/) = float_interval_div pp in
    let ff1 = one_interval + f * f in
      mk_bounded_on_int ((ddf * ff1 - two_interval * (f * (df * df))) / (ff1 * ff1)) in

  let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th) and
      f_tm = (rand o concl) diff1_th and
      _, _, bounds, d_bounds = dest_lin_approx (concl approx_th) and
      dd_bounds = (rand o concl) dd_th in

  let th0 = (MY_PROVE_HYP domain_th o MY_PROVE_HYP approx_th o 
	       MY_PROVE_HYP diff1_th o MY_PROVE_HYP dd_th o
	       INST[x_tm, x_var_real; y_tm, y_var_real; z_tm, z_var_real; w_tm, w_var_real;
		    f_tm, f_var_fun;
		    bounds, bounds_var; d_bounds, d_bounds_var; dd_bounds, dd_bounds_var])
    TAYLOR_INTERVAL_ATN' in
  let eq_th = unary_beta_eq f_tm atn_tm in
    taylor_interval_norm th0 eq_th;;


(*
let pp = 5;;
let domain_th = mk_center_domain pp one_float two_float;;
let th1 = eval_taylor_x domain_th and
    th2 = eval_taylor_sqrt pp domain_th;;
let th3 = taylor_interval_add pp th1 th2;;
taylor_interval_atn pp th1;;
taylor_interval_atn pp th2;;
taylor_interval_atn pp th3;;
*)




(*************************************)
(* taylor_interval_inv *)

let TAYLOR_INTERVAL_INV' = (UNDISCH_ALL o REWRITE_RULE[REAL_POW_2; REAL_ARITH `a pow 3 = a * a * a`] o prove)
 (`cell_domain x y z w ==>
    bounded_on_int f (x,z) f_bounds ==>
    interval_not_zero f_bounds ==>
    lin_approx (\x. inv (f x)) y bounds d_bounds ==>
    nth_diff_strong_int 2 (x,z) f ==>
    bounded_on_int (\x. (&2 * derivative f x pow 2 - nth_derivative 2 f x * f x) / (f x pow 3)) (x,z) dd_bounds
    ==> taylor_interval (\x. inv (f x)) x y z w bounds d_bounds dd_bounds`,
  SIMP_TAC[taylor_interval; cell_domain] THEN REPEAT STRIP_TAC THEN
    MATCH_MP_TAC ((SPEC_ALL o REWRITE_RULE[IMP_IMP]) second_derivative_compose_inv_bounds) THEN
    ASM_REWRITE_TAC[]);;


let taylor_interval_inv pp th1 =
  let domain_th, approx1_th, second1_th = dest_taylor_interval_th th1 in
  let diff1_th, dd1_th = dest_has_bounded_second_derivative_th second1_th in

  let f_bounds_th = eval_taylor_f_bounds pp th1 and
      approx_th = lin_approx_inv pp approx1_th in
  let f_bounds = (rand o concl) f_bounds_th in
  let cond_th = check_interval_not_zero f_bounds in

  let dd_th =
    let dest = dest_bounded_on_int_raw in
    let f = dest f_bounds_th and
	df = (dest o eval_taylor_df_bounds pp) th1 and
	ddf = dest dd1_th and
	(-) = float_interval_sub pp and
	( * ) = float_interval_mul pp and
	(/) = float_interval_div pp in
      mk_bounded_on_int ((two_interval * (df * df) - ddf * f) / (f * (f * f))) in

  let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th) and
      f_tm = (rand o concl) diff1_th and
      _, _, bounds, d_bounds = dest_lin_approx (concl approx_th) and
      dd_bounds = (rand o concl) dd_th in

  let th0 = (MY_PROVE_HYP domain_th o MY_PROVE_HYP approx_th o 
	       MY_PROVE_HYP f_bounds_th o MY_PROVE_HYP cond_th o
	       MY_PROVE_HYP diff1_th o MY_PROVE_HYP dd_th o
	       INST[x_tm, x_var_real; y_tm, y_var_real; z_tm, z_var_real; w_tm, w_var_real;
		    f_tm, f_var_fun; f_bounds, f_bounds_var;
		    bounds, bounds_var; d_bounds, d_bounds_var; dd_bounds, dd_bounds_var])
    TAYLOR_INTERVAL_INV' in
  let eq_th = unary_beta_eq f_tm inv_tm in
    taylor_interval_norm th0 eq_th;;


(*
let pp = 5;;
let domain_th = mk_center_domain pp one_float two_float;;
let th1 = eval_taylor_x domain_th and
    th2 = eval_taylor_sqrt pp domain_th;;
let th3 = taylor_interval_add pp th1 th2;;
taylor_interval_inv pp th1;;
taylor_interval_inv pp th2;;
taylor_interval_inv pp th3;;
*)



(*************************************)
(* taylor_interval_sqrt *)

let TAYLOR_INTERVAL_SQRT' = (UNDISCH_ALL o REWRITE_RULE[REAL_POW_2; REAL_ARITH `a pow 3 = a * a * a`] o prove)
 (`cell_domain x y z w ==>
    bounded_on_int f (x,z) f_bounds ==>
    interval_pos f_bounds ==>
    lin_approx (\x. sqrt (f x)) y bounds d_bounds ==>
    nth_diff_strong_int 2 (x,z) f ==>
    bounded_on_int (\x. (&2 * nth_derivative 2 f x * f x - derivative f x pow 2) /
		      (&4 * sqrt (f x pow 3))) (x,z) dd_bounds ==>
    taylor_interval (\x. sqrt (f x)) x y z w bounds d_bounds dd_bounds`,
  SIMP_TAC[taylor_interval; cell_domain] THEN REPEAT STRIP_TAC THEN
    MATCH_MP_TAC ((SPEC_ALL o REWRITE_RULE[IMP_IMP]) second_derivative_compose_sqrt_bounds) THEN
    ASM_REWRITE_TAC[]);;


let taylor_interval_sqrt pp th1 =
  let domain_th, approx1_th, second1_th = dest_taylor_interval_th th1 in
  let diff1_th, dd1_th = dest_has_bounded_second_derivative_th second1_th in

  let f_bounds_th = eval_taylor_f_bounds pp th1 and
      approx_th = lin_approx_sqrt pp approx1_th in
  let f_bounds = (rand o concl) f_bounds_th in
  let cond_th = check_interval_pos f_bounds in

  let dd_th =
    let dest = dest_bounded_on_int_raw in
    let f = dest f_bounds_th and
	df = (dest o eval_taylor_df_bounds pp) th1 and
	ddf = dest dd1_th and
	(-) = float_interval_sub pp and
	( * ) = float_interval_mul pp and
	(/) = float_interval_div pp and
	sqrt = float_interval_sqrt pp in
      mk_bounded_on_int ((two_interval * (ddf * f) - df * df) / (four_interval * sqrt (f * (f * f)))) in

  let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th) and
      f_tm = (rand o concl) diff1_th and
      _, _, bounds, d_bounds = dest_lin_approx (concl approx_th) and
      dd_bounds = (rand o concl) dd_th in

  let th0 = (MY_PROVE_HYP domain_th o MY_PROVE_HYP approx_th o 
	       MY_PROVE_HYP f_bounds_th o MY_PROVE_HYP cond_th o
	       MY_PROVE_HYP diff1_th o MY_PROVE_HYP dd_th o
	       INST[x_tm, x_var_real; y_tm, y_var_real; z_tm, z_var_real; w_tm, w_var_real;
		    f_tm, f_var_fun; f_bounds, f_bounds_var;
		    bounds, bounds_var; d_bounds, d_bounds_var; dd_bounds, dd_bounds_var])
    TAYLOR_INTERVAL_SQRT' in
  let eq_th = unary_beta_eq f_tm sqrt_tm in
    taylor_interval_norm th0 eq_th;;


(*
let pp = 5;;
let domain_th = mk_center_domain pp one_float two_float;;
let th1 = eval_taylor_x domain_th and
    th2 = eval_taylor_sqrt pp domain_th;;
let th3 = taylor_interval_add pp th1 th2;;
taylor_interval_sqrt pp th1;;
taylor_interval_sqrt pp th2;;
taylor_interval_sqrt pp th3;;
*)


(*************************************)
(* taylor_interval_acs *)

let TAYLOR_INTERVAL_ACS' = 
  (UNDISCH_ALL o PURE_ONCE_REWRITE_RULE[iabs_lt1_lemma] o
     REWRITE_RULE[REAL_POW_2; REAL_ARITH `a pow 3 = a * a * a`] o prove)
 (`cell_domain x y z w ==>
    bounded_on_int f (x,z) f_bounds ==>
    (iabs f_bounds < &1 <=> T)==>
    lin_approx (\x. acs (f x)) y bounds d_bounds ==>
    nth_diff_strong_int 2 (x,z) f ==>
    bounded_on_int (\x. --((nth_derivative 2 f x * (&1 - f x * f x) + f x * derivative f x pow 2) /
			     sqrt ((&1 - f x * f x) pow 3))) (x,z) dd_bounds ==>
    taylor_interval (\x. acs (f x)) x y z w bounds d_bounds dd_bounds`,
  SIMP_TAC[taylor_interval; cell_domain] THEN REPEAT STRIP_TAC THEN
    MATCH_MP_TAC ((SPEC_ALL o REWRITE_RULE[IMP_IMP]) second_derivative_compose_acs_bounds) THEN
    ASM_REWRITE_TAC[]);;


let taylor_interval_acs pp th1 =
  let domain_th, approx1_th, second1_th = dest_taylor_interval_th th1 in
  let diff1_th, dd1_th = dest_has_bounded_second_derivative_th second1_th in

  let f_bounds_th = eval_taylor_f_bounds pp th1 and
      approx_th = lin_approx_acs pp approx1_th in
  let f_bounds = (rand o concl) f_bounds_th in
  let cond_th = check_interval_iabs f_bounds one_float in

  let dd_th =
    let dest = dest_bounded_on_int_raw in
    let f = dest f_bounds_th and
	df = (dest o eval_taylor_df_bounds pp) th1 and
	ddf = dest dd1_th and
	(+) = float_interval_add pp and
	(-) = float_interval_sub pp and
	( * ) = float_interval_mul pp and
	(/) = float_interval_div pp and
	sqrt = float_interval_sqrt pp and
	neg = float_interval_neg in
    let ff1 = one_interval - f * f in
      mk_bounded_on_int (neg ((ddf * ff1 + f * (df * df)) / sqrt (ff1 * (ff1 * ff1)))) in

  let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th) and
      f_tm = (rand o concl) diff1_th and
      _, _, bounds, d_bounds = dest_lin_approx (concl approx_th) and
      dd_bounds = (rand o concl) dd_th in

  let th0 = (MY_PROVE_HYP domain_th o MY_PROVE_HYP approx_th o 
	       MY_PROVE_HYP f_bounds_th o MY_PROVE_HYP cond_th o
	       MY_PROVE_HYP diff1_th o MY_PROVE_HYP dd_th o
	       INST[x_tm, x_var_real; y_tm, y_var_real; z_tm, z_var_real; w_tm, w_var_real;
		    f_tm, f_var_fun; f_bounds, f_bounds_var;
		    bounds, bounds_var; d_bounds, d_bounds_var; dd_bounds, dd_bounds_var])
    TAYLOR_INTERVAL_ACS' in
  let eq_th = unary_beta_eq f_tm acs_tm in
    taylor_interval_norm th0 eq_th;;


(*
let pp = 5;;
let domain_th = mk_center_domain pp (mk_float (-5) 49) (mk_float 3 49);;
let th1 = eval_taylor_x domain_th and
    th2 = eval_taylor_atn pp domain_th;;
let th3 = taylor_interval_mul pp th1 th2;;

taylor_interval_acs pp th1;;
taylor_interval_acs pp th2;;
taylor_interval_acs pp th3;;
*)







(****************************************)
(* Composition *)

let dest_bounded_on_int_tm tm =
  let lhs, f_bounds = dest_comb tm in
  let lhs2, int_tm = dest_comb lhs in
    rand lhs2, int_tm, f_bounds;;

let dest_interval_arith tm =
  let lhs, int_tm = dest_comb tm in
    rand lhs, int_tm;;


let BOUNDED_SECOND_DERIVATIVE_COMPOSE' = (UNDISCH_ALL o REWRITE_RULE[REAL_POW_2] o SPEC_ALL) 
  second_derivative_compose_bounds;;

let BOUNDED_ON_INT_COMPOSE' = (UNDISCH_ALL o SPEC_ALL) bounded_on_int_compose;;

let TAYLOR_INTERVAL_NARROW' = (UNDISCH_ALL o PURE_ONCE_REWRITE_RULE[REAL_ARITH `a <= b <=> (a <= b <=> T)`] o
				 SPEC_ALL) taylor_interval_narrow;;

let BOUNDED_ON_INT_IMP_INTERVAL' = (UNDISCH_ALL o SPEC_ALL) bounded_on_int_imp_interval_arith;;

let LIN_APPROX_COMPOSE' = (UNDISCH_ALL o SPEC_ALL) lin_approx_compose;;



let taylor_interval_compose pp eval_f eval_g domain_th =
  let g_taylor = eval_g pp domain_th in
  let _, lin_g, second_g =  dest_taylor_interval_th g_taylor in
  let dg_th, g2_th = dest_has_bounded_second_derivative_th second_g in
  let g0_th = eval_taylor_f_bounds pp g_taylor in
  let g1_th = eval_taylor_df_bounds pp g_taylor in
  let g_tm, y_tm, g_bounds0, _ = dest_lin_approx (concl lin_g) in

  (* narrow domain *)
  let g0_lo, g0_hi = dest_pair g_bounds0 in
  let fu_domain = mk_center_domain (pp + 1) g0_lo g0_hi in
  let w0_tm = (rand o concl) fu_domain in

  (* wide domain *)
  let _, u_tm, _, _ = dest_cell_domain (concl fu_domain) in
  let g_lo, g_hi = (dest_pair o rand o concl) g0_th in
  let fw_domain = mk_cell_domain pp g_lo u_tm g_hi in
  let w_tm = (rand o concl) fw_domain in

  (* taylor_interval f *)
  let f_taylor = eval_f pp fw_domain in
  let f1_th = eval_taylor_df_bounds pp f_taylor in
  let _, lin_f, second_f = dest_taylor_interval_th f_taylor in
  let f_tm, _, f_bounds0, df_bounds0 = dest_lin_approx (concl lin_f) in
  let df_th, f2_th = dest_has_bounded_second_derivative_th second_f in

  (* has_bounded_second_derivative (\x. f (g x)) *)
  let _, int_tm, g_bounds = dest_bounded_on_int_tm (concl g0_th) in
  let df_tm, _, df_bounds = dest_bounded_on_int_tm (concl f1_th) in
  let ddf_tm, _, ddf_bounds = dest_bounded_on_int_tm (concl f2_th) in

  let f2_g_bounded_th = (MY_PROVE_HYP g0_th o MY_PROVE_HYP f2_th o
			   INST[ddf_tm, f_var_fun; g_tm, g_var_fun; g_bounds, g_bounds_var; 
				ddf_bounds, f_bounds_var; int_tm, int_var]) BOUNDED_ON_INT_COMPOSE' in
  let f1_g_bounded_th = (MY_PROVE_HYP g0_th o MY_PROVE_HYP f1_th o
			   INST[df_tm, f_var_fun; g_tm, g_var_fun; g_bounds, g_bounds_var;
				df_bounds, f_bounds_var; int_tm, int_var]) BOUNDED_ON_INT_COMPOSE' in

  let fg2_bounded_th0 = 
    let ( * ) = bounded_on_int_mul pp and
	(+) = bounded_on_int_add pp in
      (f2_g_bounded_th * (g1_th * g1_th) + f1_g_bounded_th * g2_th) in

  let dd_bounds = (rand o concl) fg2_bounded_th0 in
  let fg2_bounded_th = (MY_PROVE_HYP fg2_bounded_th0 o MY_PROVE_HYP g0_th o 
			  MY_PROVE_HYP dg_th o MY_PROVE_HYP df_th o
			  INST[f_tm, f_var_fun; g_tm, g_var_fun; int_tm, int_var; 
			       dd_bounds, dd_bounds_var; g_bounds, g_bounds_var]) 
    BOUNDED_SECOND_DERIVATIVE_COMPOSE' in

  (* lin_approx (\x. f (g x)) *)
  let fu_taylor = (MY_PROVE_HYP (float_le g_lo g0_lo) o MY_PROVE_HYP (float_le g0_hi g_hi) o
		     MY_PROVE_HYP fu_domain o MY_PROVE_HYP f_taylor o
		     INST[g0_lo, x0_var_real; u_tm, y_var_real; g0_hi, z0_var_real; w0_tm, w0_var_real;
			  f_tm, f_var_fun; g_lo, x_var_real; g_hi, z_var_real; w_tm, w_var_real;
			  f_bounds0, f_bounds_var; df_bounds0, df_bounds_var; ddf_bounds, dd_bounds_var])
    TAYLOR_INTERVAL_NARROW' in

  let df0_th = 
    let _, _, f_second = dest_taylor_interval_th fu_taylor in
      (fst o dest_has_bounded_second_derivative_th) f_second in
  let f00_th = eval_taylor_f_bounds pp fu_taylor in
  let f01_th = eval_taylor_df_bounds pp fu_taylor in
  let fg_bounds0 = (rand o concl) f00_th in
  let dfg_bounds0 = (rand o concl) f01_th in

  let dg0_th, int_g_th, int_dg_th = lin_approx_components lin_g in
  let gy_tm, _ = dest_interval_arith (concl int_g_th) in

  let int_df_g_th = (MY_PROVE_HYP int_g_th o MY_PROVE_HYP f01_th o
		       INST[df_tm, f_var_fun; gy_tm, y_var_real; g_bounds0, int_var; dfg_bounds0, f_bounds_var])
    BOUNDED_ON_INT_IMP_INTERVAL' in

  let int_dfg_th = float_interval_mul pp int_dg_th int_df_g_th in
  let dfg_bounds00 = rand (concl int_dfg_th) in

  let fg_lin_th = (MY_PROVE_HYP int_dfg_th o MY_PROVE_HYP int_g_th o MY_PROVE_HYP f00_th o 
		     MY_PROVE_HYP dg0_th o MY_PROVE_HYP df0_th o
		     INST[f_tm, f_var_fun; g_tm, g_var_fun; y_tm, y_var_real; dfg_bounds00, d_bounds_var;
			  g_bounds0, g_bounds_var; fg_bounds0, f_bounds_var]) LIN_APPROX_COMPOSE' in

  (* beta reduction *)
  let gx_tm = mk_comb (g_tm, x_var_real) in
  let gx_th = if is_abs g_tm then BETA gx_tm else REFL gx_tm in

  let fgx_th0 = MK_COMB (REFL f_tm, gx_th) in
  let fgx_eq_th = ABS x_var_real 
    (if is_abs f_tm then TRANS fgx_th0 (BETA_CONV (rand (concl fgx_th0))) else fgx_th0) in

  let fg_lin_th1 = norm_lin_approx fg_lin_th fgx_eq_th and
      fg2_bounded_th1 = norm_second_derivative fg2_bounded_th fgx_eq_th in

    mk_taylor_interval_th domain_th fg_lin_th1 fg2_bounded_th1;;