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


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

(* Given two terms (\x. f1 x) and (\x. f2 x) and a binary operation,
   returns an equality theorem 
   \x. (\x. f1 x) x (op) (\x. f2 x) x = \x. f1 x (op) f2 x *)
let binary_beta_eq f1_tm f2_tm op_tm =
  let beta_tm1, beta_tm2 = mk_comb (f1_tm, x_var_real), mk_comb (f2_tm, x_var_real) in
  let beta_th1 = if is_abs f1_tm then BETA beta_tm1 else REFL beta_tm1 and
      beta_th2 = if is_abs f2_tm then BETA beta_tm2 else REFL beta_tm2 in
    ABS x_var_real (MK_COMB (AP_TERM op_tm beta_th1, beta_th2));;


let unary_beta_eq f_tm op_tm =
  let beta_tm = mk_comb (f_tm, x_var_real) in
  let beta_th = if is_abs f_tm then BETA beta_tm else REFL beta_tm in
    ABS x_var_real (AP_TERM op_tm beta_th);;



(*****************************************)
(* mk_lin_approx, dest_lin_approx *)


let MK_LIN_APPROX' = (RULE o MATCH_MP EQ_IMP o SYM o SPEC_ALL) lin_approx_eq;;
let DEST_LIN_APPROX' = (RULE o MATCH_MP EQ_IMP o SPEC_ALL) lin_approx_eq;;


let dest_lin_approx approx_tm =
  let lhs, df_bounds = dest_comb approx_tm in
  let lhs2, f_bounds = dest_comb lhs in
  let lhs3, x_tm = dest_comb lhs2 in
  let f_tm = rand lhs3 in
    f_tm, x_tm, f_bounds, df_bounds;;


let lin_approx_components approx_th =
  let f_tm, x_tm, f_bounds, df_bounds = dest_lin_approx (concl approx_th) in
  let th0 = INST[f_tm, f_var_fun; x_tm, x_var_real; 
		 f_bounds, f_bounds_var; df_bounds, df_bounds_var] DEST_LIN_APPROX' in
  let th1 = MY_PROVE_HYP approx_th th0 in
  let [r1; r2; r3] = CONJUNCTS th1 in
    (r1, r2, r3);;


let mk_lin_approx f_bounds_th df_bounds_th diff_th =
  let lhs, f_bounds = dest_comb (concl f_bounds_th) in
  let fx_tm = rand lhs in
  let f_tm, x_tm = dest_comb fx_tm in
  let df_bounds = (rand o concl) df_bounds_th in
    (MY_PROVE_HYP f_bounds_th o MY_PROVE_HYP df_bounds_th o MY_PROVE_HYP diff_th o
       INST[f_tm, f_var_fun; x_tm, x_var_real; f_bounds, f_bounds_var; df_bounds, df_bounds_var])
      MK_LIN_APPROX';;
    
    

(***************************)
(* f(x) = x *)
let LIN_APPROX_X' = RULE lin_approx_x;;

let eval_lin_approx_x x_tm =
  INST[x_tm, x_var_real] LIN_APPROX_X';;


(***************************)
(* f(x) = c *)
let LIN_APPROX_CONST' = RULE lin_approx_const;;


let c_var_real = `c:real`;;


let eval_lin_approx_const c_tm x_tm =
  INST[c_tm, c_var_real; x_tm, x_var_real] LIN_APPROX_CONST';;


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

let MK_LIN_APPROX_ATN' = (UNDISCH_ALL o REWRITE_RULE[derivative_atn] o DISCH_ALL o
			    MY_PROVE_HYP (SPEC_ALL REAL_DIFFERENTIABLE_AT_ATN) o 
			    INST[`atn`, `f:real->real`]) MK_LIN_APPROX';;


let eval_lin_approx_atn pp x_tm =
  let x_th = mk_const_interval x_tm in
  let d_th = float_interval_inv pp (float_interval_add pp one_interval (float_interval_mul pp x_th x_th)) in
  let atn_th = float_interval_atn pp x_th in
  let f_bounds = (rand o concl) atn_th in
  let df_bounds = (rand o concl) d_th in
  let th0 = INST[x_tm, x_var_real; f_bounds, f_bounds_var; df_bounds, df_bounds_var] MK_LIN_APPROX_ATN' in
    MY_PROVE_HYP d_th (MY_PROVE_HYP atn_th th0);;


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

let MK_LIN_APPROX_INV' = (UNDISCH_ALL o prove)
  (`(x = &0 <=> F) ==> 
     interval_arith (inv x) f_bounds ==>
     interval_arith (--(inv x * inv x)) df_bounds ==>
     lin_approx inv x f_bounds df_bounds`,
  REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[lin_approx_eq] THEN
    ASM_SIMP_TAC[derivative_inv; REAL_DIFFERENTIABLE_AT_INV; REAL_INV_MUL]);;

let eval_lin_approx_inv pp x_tm =
  let x_th = mk_const_interval x_tm in
  let inv_th = float_interval_inv pp x_th in
  let d_th =
    let neg = float_interval_neg and
	( * ) = float_interval_mul pp in
      neg (inv_th * inv_th) in
  let f_bounds = (rand o concl) inv_th and
      df_bounds = (rand o concl) d_th in
  let xn0_th = float_eq0 x_tm in
    (MY_PROVE_HYP xn0_th o MY_PROVE_HYP inv_th o MY_PROVE_HYP d_th o
       INST[x_tm, x_var_real; f_bounds, f_bounds_var; df_bounds, df_bounds_var])
      MK_LIN_APPROX_INV';;
      


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

let MK_LIN_APPROX_SQRT' = (UNDISCH_ALL o prove)
  (`(&0 < x <=> T) ==> 
     interval_arith (sqrt x) f_bounds ==>
     interval_arith (inv(&2 * sqrt x)) df_bounds ==>
     lin_approx sqrt x f_bounds df_bounds`,
  REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[lin_approx_eq] THEN
    ASM_SIMP_TAC[derivative_sqrt; REAL_DIFFERENTIABLE_AT_SQRT]);;


let eval_lin_approx_sqrt pp x_tm =
  let x_th = mk_const_interval x_tm in
  let sqrt_th = float_interval_sqrt pp x_th in
  let d_th =
    let ( * ) = float_interval_mul pp and
	inv = float_interval_inv pp in
      inv (two_interval * sqrt_th) in
  let f_bounds = (rand o concl) sqrt_th and
      df_bounds = (rand o concl) d_th in
  let x_pos_th = float_gt0 x_tm in
    (MY_PROVE_HYP x_pos_th o MY_PROVE_HYP sqrt_th o MY_PROVE_HYP d_th o
       INST[x_tm, x_var_real; f_bounds, f_bounds_var; df_bounds, df_bounds_var])
      MK_LIN_APPROX_SQRT';;

(*      
eval_lin_approx_sqrt 5 two_float;;
let x = 2;;
*)


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

let MK_LIN_APPROX_ACS' = (UNDISCH_ALL o GEN_REWRITE_RULE (LAND_CONV o LAND_CONV o RAND_CONV) [(SYM o FLOAT_TO_NUM_CONV) one_float] o prove)
  (`(abs x < &1 <=> T) ==> 
     interval_arith (acs x) f_bounds ==>
     interval_arith (--inv (sqrt (&1 - x * x))) df_bounds ==>
     lin_approx acs x f_bounds df_bounds`,
  REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[lin_approx_eq] THEN
    ASM_SIMP_TAC[derivative_acs; REAL_DIFFERENTIABLE_AT_ACS]);;

let abs_tm = `abs`;;
let check_abs_lt_1 x_tm =
  let abs_eq = float_abs (mk_comb (abs_tm, x_tm)) in
  let lt_th = float_lt (rand (concl abs_eq)) one_float in
    if (fst o dest_const o rand o concl) lt_th <> "T" then
      failwith "check_abs_lt_1: abs < 1 <=> F"
    else
      let th0 = AP_THM (AP_TERM lt_op_real abs_eq) one_float in
	TRANS th0 lt_th;;



let eval_lin_approx_acs pp x_tm =
  let x_th = mk_const_interval x_tm in
  let acs_th = float_interval_acs pp x_th in
  let d_th =
    let ( * ) = float_interval_mul pp and
	inv = float_interval_inv pp and
	sqrt = float_interval_sqrt pp and
	neg = float_interval_neg and
	(-) = float_interval_sub pp in
      neg (inv (sqrt (one_interval - x_th * x_th))) in
  let f_bounds = (rand o concl) acs_th and
      df_bounds = (rand o concl) d_th in
  let x_abs_th = check_abs_lt_1 x_tm in
    (MY_PROVE_HYP x_abs_th o MY_PROVE_HYP acs_th o MY_PROVE_HYP d_th o
       INST[x_tm, x_var_real; f_bounds, f_bounds_var; df_bounds, df_bounds_var])
      MK_LIN_APPROX_ACS';;




(*************************)
(* lin_approx arithmetic *)
(*************************)


(********************************)
(* lin_approx_neg *)

let MK_LIN_APPROX_NEG' = (UNDISCH_ALL o prove)
  (`f real_differentiable atreal x ==>
     interval_arith (--f x) f_bounds ==>
     interval_arith (--derivative f x) df_bounds ==>
     lin_approx (\x. --f x) x f_bounds df_bounds`,
  REPEAT DISCH_TAC THEN ASM_REWRITE_TAC[lin_approx_eq] THEN
    ASM_SIMP_TAC[REAL_DIFFERENTIABLE_NEG; derivative_neg]);;

    
let lin_approx_neg approx1 =
  let diff1_th, f1_th, df1_th = lin_approx_components approx1 in
  let f1_tm = (lhand o concl) diff1_th and
      x_tm = (rand o rand o concl) diff1_th in

  let f_th = float_interval_neg f1_th and
      df_th = float_interval_neg df1_th in

  let f_bounds = (rand o concl) f_th and
      df_bounds = (rand o concl) df_th in

    (MY_PROVE_HYP diff1_th o MY_PROVE_HYP f_th o MY_PROVE_HYP df_th o
       INST[f1_tm, f_var_fun; x_tm, x_var_real;
	    f_bounds, f_bounds_var; df_bounds, df_bounds_var]) MK_LIN_APPROX_NEG';;



(*
let pp = 5;;
let approx1 = eval_lin_approx_atn pp one_float;;
let approx2 = eval_lin_approx_x one_float;;
lin_approx_neg approx1;;
lin_approx_neg approx2;; 
*)



(********************************)
(* lin_approx_add *)

let MK_LIN_APPROX_ADD' = (UNDISCH_ALL o prove)
  (`f real_differentiable atreal x ==> g real_differentiable atreal x ==>
     interval_arith (f x + g x) f_bounds ==>
     interval_arith (derivative f x + derivative g x) df_bounds ==>
     lin_approx (\x. f x + g x) x f_bounds df_bounds`,
  REPEAT DISCH_TAC THEN ASM_REWRITE_TAC[lin_approx_eq] THEN
    ASM_SIMP_TAC[REAL_DIFFERENTIABLE_ADD; derivative_add]);;

    
let lin_approx_add pp approx1 approx2 =
  let diff1_th, f1_th, df1_th = lin_approx_components approx1 and
      diff2_th, f2_th, df2_th = lin_approx_components approx2 in
  let f1_tm = (lhand o concl) diff1_th and
      f2_tm = (lhand o concl) diff2_th and
      x_tm = (rand o rand o concl) diff1_th in

  let f_th = float_interval_add pp f1_th f2_th and
      df_th = float_interval_add pp df1_th df2_th in

  let f_bounds = (rand o concl) f_th and
      df_bounds = (rand o concl) df_th in

    (MY_PROVE_HYP diff1_th o MY_PROVE_HYP diff2_th o MY_PROVE_HYP f_th o MY_PROVE_HYP df_th o
       INST[f1_tm, f_var_fun; f2_tm, g_var_fun; x_tm, x_var_real;
	    f_bounds, f_bounds_var; df_bounds, df_bounds_var]) MK_LIN_APPROX_ADD';;



(*
let pp = 5;;
let approx1 = eval_lin_approx_atn pp one_float;;
let approx2 = eval_lin_approx_x one_float;;
lin_approx_add pp (eval_lin_approx_add pp approx2 approx2) approx1;;
*)



(********************************)
(* lin_approx_sub *)

let MK_LIN_APPROX_SUB' = (UNDISCH_ALL o prove)
  (`f real_differentiable atreal x ==> g real_differentiable atreal x ==>
     interval_arith (f x - g x) f_bounds ==>
     interval_arith (derivative f x - derivative g x) df_bounds ==>
     lin_approx (\x. f x - g x) x f_bounds df_bounds`,
  REPEAT DISCH_TAC THEN ASM_REWRITE_TAC[lin_approx_eq] THEN
    ASM_SIMP_TAC[REAL_DIFFERENTIABLE_SUB; derivative_sub]);;

    
let lin_approx_sub pp approx1 approx2 =
  let diff1_th, f1_th, df1_th = lin_approx_components approx1 and
      diff2_th, f2_th, df2_th = lin_approx_components approx2 in
  let f1_tm = (lhand o concl) diff1_th and
      f2_tm = (lhand o concl) diff2_th and
      x_tm = (rand o rand o concl) diff1_th in

  let f_th = float_interval_sub pp f1_th f2_th and
      df_th = float_interval_sub pp df1_th df2_th in

  let f_bounds = (rand o concl) f_th and
      df_bounds = (rand o concl) df_th in

    (MY_PROVE_HYP diff1_th o MY_PROVE_HYP diff2_th o MY_PROVE_HYP f_th o MY_PROVE_HYP df_th o
       INST[f1_tm, f_var_fun; f2_tm, g_var_fun; x_tm, x_var_real;
	    f_bounds, f_bounds_var; df_bounds, df_bounds_var]) MK_LIN_APPROX_SUB';;



(*
let pp = 5;;
let approx1 = eval_lin_approx_atn pp one_float;;
let approx2 = eval_lin_approx_x one_float;;
lin_approx_sub pp approx1 approx2;;
*)



(********************************)
(* lin_approx_mul *)

let MK_LIN_APPROX_MUL' = (UNDISCH_ALL o prove)
  (`f real_differentiable atreal x ==> g real_differentiable atreal x ==>
     interval_arith (f x * g x) f_bounds ==>
     interval_arith (f x * derivative g x + derivative f x * g x) df_bounds ==>
     lin_approx (\x. f x * g x) x f_bounds df_bounds`,
  REPEAT DISCH_TAC THEN ASM_REWRITE_TAC[lin_approx_eq] THEN
    ASM_SIMP_TAC[REAL_DIFFERENTIABLE_MUL_ATREAL; derivative_mul]);;

    
let lin_approx_mul pp approx1 approx2 =
  let diff1_th, f1_th, df1_th = lin_approx_components approx1 and
      diff2_th, f2_th, df2_th = lin_approx_components approx2 in
  let f1_tm = (lhand o concl) diff1_th and
      f2_tm = (lhand o concl) diff2_th and
      x_tm = (rand o rand o concl) diff1_th in

  let f_th, df_th =
    let (+) = float_interval_add pp and
	( * ) = float_interval_mul pp in
      f1_th * f2_th, f1_th * df2_th + df1_th * f2_th in

  let f_bounds = (rand o concl) f_th and
      df_bounds = (rand o concl) df_th in

    (MY_PROVE_HYP diff1_th o MY_PROVE_HYP diff2_th o MY_PROVE_HYP f_th o MY_PROVE_HYP df_th o
       INST[f1_tm, f_var_fun; f2_tm, g_var_fun; x_tm, x_var_real;
	    f_bounds, f_bounds_var; df_bounds, df_bounds_var]) MK_LIN_APPROX_MUL';;


(*
let pp = 5;;
let approx1 = eval_lin_approx_atn (pp + 2) one_float;;
let approx2 = eval_lin_approx_x one_float;;
lin_approx_mul pp approx1 approx1;;
*)


(********************************)
(* lin_approx_div *)

let MK_LIN_APPROX_DIV' = (UNDISCH_ALL o prove)
  (`interval_arith (g x) g_bounds ==> interval_not_zero g_bounds ==>
     f real_differentiable atreal x ==> g real_differentiable atreal x ==>
     interval_arith (f x / g x) f_bounds ==>
     interval_arith ((derivative f x * g x - f x * derivative g x) / (g x * g x)) df_bounds ==>
     lin_approx (\x. f x / g x) x f_bounds df_bounds`,
  REPEAT DISCH_TAC THEN ASM_REWRITE_TAC[lin_approx_eq] THEN
    SUBGOAL_THEN `~((g:real->real) x = &0)` ASSUME_TAC THENL
    [
      MATCH_MP_TAC (REWRITE_RULE[IMP_IMP] interval_arith_not_zero) THEN
	EXISTS_TAC `g_bounds:real#real` THEN ASM_REWRITE_TAC[];
      ALL_TAC
    ] THEN
    ASM_SIMP_TAC[REAL_DIFFERENTIABLE_DIV_ATREAL; derivative_div]);;

    
let lin_approx_div pp approx1 approx2 =
  let diff1_th, f1_th, df1_th = lin_approx_components approx1 and
      diff2_th, f2_th, df2_th = lin_approx_components approx2 in
  let f1_tm = (lhand o concl) diff1_th and
      f2_tm = (lhand o concl) diff2_th and
      x_tm = (rand o rand o concl) diff1_th and
      g_bounds = (rand o concl) f2_th in

  let not_zero_th = check_interval_not_zero g_bounds in
  let f_th, df_th =
    let (-) = float_interval_sub pp and
	( * ) = float_interval_mul pp and
	(/) = float_interval_div pp in
      f1_th / f2_th, (df1_th * f2_th - f1_th * df2_th) / (f2_th * f2_th) in

  let f_bounds = (rand o concl) f_th and
      df_bounds = (rand o concl) df_th in

    (MY_PROVE_HYP not_zero_th o MY_PROVE_HYP f2_th o 
       MY_PROVE_HYP diff1_th o MY_PROVE_HYP diff2_th o 
       MY_PROVE_HYP f_th o MY_PROVE_HYP df_th o
       INST[f1_tm, f_var_fun; f2_tm, g_var_fun; x_tm, x_var_real;
	    f_bounds, f_bounds_var; df_bounds, df_bounds_var; 
	    g_bounds, g_bounds_var]) MK_LIN_APPROX_DIV';;


(*
let pp = 5;;
let approx1 = eval_lin_approx_atn (pp + 2) one_float;;
let approx2 = eval_lin_approx_x one_float;;
let approx3 = eval_lin_approx_atn 4 (mk_float 0 50);;
lin_approx_div pp approx1 approx2;;
lin_approx_div pp approx2 approx1;;
*)



(********************************)
(* lin_approx_atn *)

let atn_tm = `atn`;;

let MK_LIN_APPROX_COMPOSE_ATN' = (UNDISCH_ALL o prove)
  (`f real_differentiable atreal x ==>
     interval_arith (atn (f x)) f_bounds ==>
     interval_arith (derivative f x / (&1 + f x * f x)) df_bounds ==>
     lin_approx (\x. atn (f x)) x f_bounds df_bounds`,
  REPEAT DISCH_TAC THEN ASM_REWRITE_TAC[lin_approx_eq] THEN
    ASM_SIMP_TAC[derivative_compose_atn]);;


let lin_approx_atn pp approx1 =
  let diff1_th, f1_th, df1_th = lin_approx_components approx1 in
  let f1_tm = (lhand o concl) diff1_th and
      x_tm = (rand o rand o concl) diff1_th in

  let f_th, df_th =
    let (+) = float_interval_add pp and
	( * ) = float_interval_mul pp and
	(/) = float_interval_div pp and
	atn = float_interval_atn pp in
      atn f1_th, df1_th / (one_interval + f1_th * f1_th) in

  let f_bounds = (rand o concl) f_th and
      df_bounds = (rand o concl) df_th in

    (MY_PROVE_HYP diff1_th o MY_PROVE_HYP f_th o MY_PROVE_HYP df_th o
       INST[f1_tm, f_var_fun; x_tm, x_var_real;
	    f_bounds, f_bounds_var; df_bounds, df_bounds_var]) MK_LIN_APPROX_COMPOSE_ATN';;



(*
let pp = 5;;
let approx1 = eval_lin_approx_atn (pp + 2) one_float;;
let approx2 = eval_lin_approx_x one_float;;
let approx3 = eval_lin_approx_atn 4 (mk_float 0 50);;
lin_approx_atn pp approx1;;
lin_approx_atn pp approx2;;
lin_approx_atn pp approx3;;
*)



(********************************)
(* lin_approx_inv *)

let inv_tm = `inv`;;

let MK_LIN_APPROX_COMPOSE_INV' = (UNDISCH_ALL o prove)
  (`interval_arith (f x) int ==> interval_not_zero int ==>
     f real_differentiable atreal x ==>
     interval_arith (inv (f x)) f_bounds ==>
     interval_arith (--(derivative f x * inv(f x) * inv(f x))) df_bounds ==>
     lin_approx (\x. inv (f x)) x f_bounds df_bounds`,
  REWRITE_TAC[GSYM REAL_INV_MUL] THEN REWRITE_TAC[GSYM REAL_MUL_RNEG] THEN ONCE_REWRITE_TAC[REAL_MUL_SYM] THEN
    REPEAT DISCH_TAC THEN ASM_REWRITE_TAC[lin_approx_eq] THEN
    SUBGOAL_THEN `~((f:real->real) x = &0)` ASSUME_TAC THENL
    [
      MATCH_MP_TAC (REWRITE_RULE[IMP_IMP] interval_arith_not_zero) THEN
	EXISTS_TAC `int:real#real` THEN ASM_REWRITE_TAC[];
      ALL_TAC
    ] THEN
    ASM_SIMP_TAC[derivative_compose_inv]);;


let lin_approx_inv pp approx1 =
  let diff1_th, f1_th, df1_th = lin_approx_components approx1 in
  let f1_tm = (lhand o concl) diff1_th and
      x_tm = (rand o rand o concl) diff1_th and
      int_tm = (rand o concl) f1_th in

  let not_zero_th = check_interval_not_zero int_tm in
  let f_th = float_interval_inv pp f1_th in
  let df_th =
    let ( * ) = float_interval_mul pp and
	neg = float_interval_neg in
      neg (df1_th * (f_th * f_th)) in

  let f_bounds = (rand o concl) f_th and
      df_bounds = (rand o concl) df_th in

    (MY_PROVE_HYP diff1_th o MY_PROVE_HYP f_th o MY_PROVE_HYP df_th o
       MY_PROVE_HYP f1_th o MY_PROVE_HYP not_zero_th o
       INST[f1_tm, f_var_fun; x_tm, x_var_real;
	    f_bounds, f_bounds_var; df_bounds, df_bounds_var;
	    int_tm, int_var]) MK_LIN_APPROX_COMPOSE_INV';;



(*
let pp = 5;;
let approx1 = eval_lin_approx_atn (pp + 2) one_float;;
let approx2 = eval_lin_approx_x one_float;;
let approx3 = eval_lin_approx_atn 4 (mk_float 0 50);;
lin_approx_inv pp approx1;;
lin_approx_inv pp approx2;;
lin_approx_inv pp approx3;;
*)


(********************************)
(* lin_approx_sqrt *)

let sqrt_tm = `sqrt`;;

let MK_LIN_APPROX_COMPOSE_SQRT' = (UNDISCH_ALL o prove)
  (`interval_arith (f x) int ==> interval_pos int ==>
     f real_differentiable atreal x ==>
     interval_arith (sqrt (f x)) f_bounds ==>
     interval_arith (derivative f x / (&2 * sqrt (f x))) df_bounds ==>
     lin_approx (\x. sqrt (f x)) x f_bounds df_bounds`,
   REPEAT DISCH_TAC THEN ASM_REWRITE_TAC[lin_approx_eq] THEN
    SUBGOAL_THEN `&0 < ((f:real->real) x)` ASSUME_TAC THENL
    [
      MATCH_MP_TAC (REWRITE_RULE[IMP_IMP] interval_arith_pos) THEN
	EXISTS_TAC `int:real#real` THEN ASM_REWRITE_TAC[];
      ALL_TAC
    ] THEN
    ASM_SIMP_TAC[derivative_compose_sqrt]);;


let lin_approx_sqrt pp approx1 =
  let diff1_th, f1_th, df1_th = lin_approx_components approx1 in
  let f1_tm = (lhand o concl) diff1_th and
      x_tm = (rand o rand o concl) diff1_th and
      int_tm = (rand o concl) f1_th in

  let pos_th = check_interval_pos int_tm in
  let f_th = float_interval_sqrt pp f1_th in
  let df_th =
    let ( * ) = float_interval_mul pp and
	(/) = float_interval_div pp in
      (df1_th / (two_interval * f_th)) in

  let f_bounds = (rand o concl) f_th and
      df_bounds = (rand o concl) df_th in

    (MY_PROVE_HYP diff1_th o MY_PROVE_HYP f_th o MY_PROVE_HYP df_th o
       MY_PROVE_HYP f1_th o MY_PROVE_HYP pos_th o
       INST[f1_tm, f_var_fun; x_tm, x_var_real;
	    f_bounds, f_bounds_var; df_bounds, df_bounds_var;
	    int_tm, int_var]) MK_LIN_APPROX_COMPOSE_SQRT';;



(*
let pp = 5;;
let approx1 = eval_lin_approx_atn (pp + 2) one_float;;
let approx2 = eval_lin_approx_x one_float;;
let approx3 = eval_lin_approx_atn 4 (mk_float 0 50);;
lin_approx_sqrt pp approx1;;
lin_approx_sqrt pp approx2;;
lin_approx_sqrt pp approx3;;
*)



(********************************)
(* lin_approx_acs *)

let acs_tm = `acs`;;

let MK_LIN_APPROX_COMPOSE_ACS' = 
  (UNDISCH_ALL o 
     GEN_REWRITE_RULE (RAND_CONV o LAND_CONV o LAND_CONV o RAND_CONV) [GSYM (FLOAT_TO_NUM_CONV one_float)] o 
     prove)
  (`interval_arith (f x) int ==> (iabs int < &1 <=> T) ==>
     f real_differentiable atreal x ==>
     interval_arith (acs (f x)) f_bounds ==>
     interval_arith (--(derivative f x / sqrt (&1 - f x * f x))) df_bounds ==>
     lin_approx (\x. acs (f x)) x f_bounds df_bounds`,
   REPEAT DISCH_TAC THEN ASM_REWRITE_TAC[lin_approx_eq] THEN
    SUBGOAL_THEN `abs ((f:real->real) x) < &1` ASSUME_TAC THENL
    [
      MATCH_MP_TAC (REWRITE_RULE[IMP_IMP] interval_arith_abs) THEN
	EXISTS_TAC `int:real#real` THEN ASM_REWRITE_TAC[];
      ALL_TAC
    ] THEN
    ASM_SIMP_TAC[derivative_compose_acs]);;


let lin_approx_acs pp approx1 =
  let diff1_th, f1_th, df1_th = lin_approx_components approx1 in
  let f1_tm = (lhand o concl) diff1_th and
      x_tm = (rand o rand o concl) diff1_th and
      int_tm = (rand o concl) f1_th in

  let abs_th = check_interval_iabs int_tm one_float in
  let f_th = float_interval_acs pp f1_th in
  let df_th =
    let (-) = float_interval_sub pp and
	( * ) = float_interval_mul pp and
	(/) = float_interval_div pp and
	neg = float_interval_neg and
	sqrt = float_interval_sqrt pp in
      neg (df1_th / sqrt (one_interval - f1_th * f1_th)) in

  let f_bounds = (rand o concl) f_th and
      df_bounds = (rand o concl) df_th in

    (MY_PROVE_HYP diff1_th o MY_PROVE_HYP f_th o MY_PROVE_HYP df_th o
       MY_PROVE_HYP f1_th o MY_PROVE_HYP abs_th o
       INST[f1_tm, f_var_fun; x_tm, x_var_real;
	    f_bounds, f_bounds_var; df_bounds, df_bounds_var;
	    int_tm, int_var]) MK_LIN_APPROX_COMPOSE_ACS';;



(*
let pp = 5;;
let approx1 = eval_lin_approx_atn (pp + 2) one_float;;
let approx2 = eval_lin_approx_x one_float;;
let approx3 = eval_lin_approx_atn 4 (mk_float 0 50);;
lin_approx_acs pp approx1;;
lin_approx_acs pp approx2;;
lin_approx_acs pp approx3;;
*)