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


(************)
(* Examples *)
(************)

(* f(x) = x * x - atn x; f(x) < 0 when 0.01 <= x <= 0.8 *)
prove_ineq1d 5 `x * x - atn x` `#0.01` `#0.8`;;

(* f(x) = x * x * x - atn (x * x); f(x) < 0 when 0.01 <= x <= 0.8 *)
prove_ineq1d 5 `x * x * x - atn (x * x)` `#0.01` `#0.8`;;

(* f(x) = atn x - x; f(x) < 0 when 0.001 <= x <= sqrt(2) *)
prove_ineq1d 10 `atn x - x` `#0.001` `sqrt (&2) + pi`;;


(****************************************)
(* arclength #2.52 #2.52 (&2) >= #0.816 *)
(* arclength (&2) #2.52 (&2) - arclength #2.52 #2.52 (&2) >= #0.073 *)
(* arclength a b c = acs ((a * a + b * b - c * c) / (&2 * a * b)) *)

let arclength_tm = `acs ((a * a + b * b - c * c) / (&2 * a * b))`;;
let tm1 = subst[`#2.52`, `a:real`; `#2.52`, `b:real`; `&2`, `c:real`] arclength_tm and
    tm2 = subst[`&2`, `a:real`; `#2.52`, `b:real`; `&2`, `c:real`] arclength_tm;;

let ineq1, ineq2 =
  let (-) = mk_binop sub_op_real in
    `#0.816` - tm1, `#0.073` - (tm2 - tm1);;

prove_ineq1d 6 ineq1 `&0` `&0`;;
prove_ineq1d 6 ineq2 `&0` `&0`;;


(********************************************)
(* !x. &2 <= x /\ x <= #2.52                                           *)
(*   ==> arclength (&2) x (&2) - arclength (#2.52) x (&2) >= #0.073`   *)

let ineq =
  let tm1 = subst[`&2`, `a:real`; `x:real`, `b:real`; `&2`, `c:real`] arclength_tm and
      tm2 = subst[`#2.52`, `a:real`; `x:real`, `b:real`; `&2`, `c:real`] arclength_tm and
      (-) = mk_binop sub_op_real in
    `#0.073` - (tm1 - tm2);;

prove_ineq1d 7 ineq `&2` `#2.52`;;


(***********************)
(* !x. &2 <= x /\ x <= #2.52  ==> &1 / &4 * --inv (sqrt (&1 - (x / &4) pow 2)) <= -- #0.28 *)

prove_ineq1d 4 `#0.28 + &1 / &4 * --inv (sqrt (&1 - (x * inv(&4)) * (x * inv(&4))))` `&2` `#2.52`;;

(* !x. &2 <= x /\ x <= #2.52 ==>
	inv (sqrt (&1 - ((&3969 / &625 + x * x - &4) / (&2 * #2.52 * x)) pow 2)) <= &2 *)

let ineq = (rand o concl o REWRITE_CONV[REAL_POW_2])
  `inv (sqrt (&1 - ((&3969 / &625 + x * x - &4) / (&2 * #2.52 * x)) pow 2)) - &2`;;

prove_ineq1d 3 ineq `&2` `#2.52`;;



(***********************************************)
(* |atn x - x / (1 + 0.28 * x * x)| < 0.005 when x IN [0, 1] *)
let abs_lemma = prove(`f - eps < &0 ==> -- eps - f < &0 ==> abs f < eps`, REAL_ARITH_TAC);;

let tm = `atn x - x / (&1 + #0.28 * x * x)`;;
let eps = `#0.005`;;
let hi_tm = mk_binop sub_op_real tm eps;;
let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) tm;;
let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 5 hi_tm `&0` `&1`);;
let lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 6 lo_tm `&0` `&1`);;
(GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;



(*************************************)
(* A polynomial approximation of atn *)
let tm = `x * (&1 - (x * x) * (&11184811 / &33554432 - (x * x) * (&13421773 / &67108864))) - atn x`;;
let eps = `#0.0000001`;;
let hi_tm = mk_binop sub_op_real tm eps;;
let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) tm;;
let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 8 hi_tm `&0` `&1 / &30`) and
    lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 8 lo_tm `&0` `&1 / &30`) and
    hi_th' = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 8 hi_tm `-- &1 / &30` `&0`) and
    lo_th' = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 8 lo_tm `-- &1 / &30` `&0`);;


(GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;



(***********************************)
(* A rational approximation of atn *)


let atan1 = `let P00, P01, P02 = #4.26665376811246382, #3.291955407318148, #0.281939359003812 in
  let Q00, Q01 = #4.26665376811354027, #4.714173328637605 in
  let x1 = x * x in
  let x2 = x1 * x1 in
    x * (x1 * P01 + x2 * P02 + P00) / (x1 * Q01 + x2 + Q00) - atn x`;;
let atan1_tm = (rand o concl o REPEATC let_CONV) atan1;;

let eps = `#0.00001`;;
let hi_tm = mk_binop sub_op_real atan1_tm eps;;
let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) atan1_tm;;
let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 hi_tm `&0` `#0.2`);;
let lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 lo_tm `&0` `#0.2`);;
(GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;


(*
(* A higher precision approximation *)
let eps = `#0.0000000001`;;
let hi_tm = mk_binop sub_op_real atan1_tm eps;;
let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) atan1_tm;;
let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 15 hi_tm `&0` `#0.2`);;
let lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 15 lo_tm `&0` `#0.2`);;
(GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;
*)



(****************************************************************)
(* Test 1: f(x) = x * x - atn x; f(x) < 0 when 0.01 <= x <= 0.8 *)
(****************************************************************)

llet eps = `#0.00001`;;
let hi_tm = mk_binop sub_op_real atan1_tm eps;;
let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) atan1_tm;;
let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 hi_tm `&0` `#0.2`);;
let lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 lo_tm `&0` `#0.2`);;
(GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;
let eps = `#0.00001`;;
let hi_tm = mk_binop sub_op_real atan1_tm eps;;
let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) atan1_tm;;
let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 hi_tm `&0` `#0.2`);;
let lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 lo_tm `&0` `#0.2`);;
(GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;
let eps = `#0.00001`;;
let hi_tm = mk_binop sub_op_real atan1_tm eps;;
let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) atan1_tm;;
let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 hi_tm `&0` `#0.2`);;
let lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 lo_tm `&0` `#0.2`);;
(GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;
et x_tm = mk_float 1 48;;
let z_tm = mk_float 8 49;;

let domain_th = mk_center_domain 5 x_tm z_tm;;
let _, y_tm, _, _ = dest_cell_domain (concl domain_th);;


let domain1_th = mk_center_domain 5 x_tm y_tm;;
let _, y1_tm, _, _ = dest_cell_domain (concl domain1_th);;
let domain11_th = mk_center_domain 5 x_tm y1_tm;;
let domain12_th = mk_center_domain 5 y1_tm y_tm;;

let domain2_th = mk_center_domain 5 y_tm z_tm;;


let eval_f pp domain_th =
  let int_x = eval_taylor_x domain_th in
  let int_atn = eval_taylor_atn pp domain_th in
  let int_x2 = taylor_interval_mul pp int_x int_x in
    taylor_interval_sub pp int_x2 int_atn;;


let int_th = eval_f 5 domain_th;;
let int_th11 = eval_f 5 domain11_th;;
let int_th12 = eval_f 5 domain12_th;;
let int_th2 = eval_f 5 domain2_th;;

taylor_cell_pass 5 int_th;;
let pass11 = taylor_cell_pass 5 int_th11;;
let pass12 = taylor_cell_pass 5 int_th12;;
let pass2 = taylor_cell_pass 5 int_th2;;

let pass1 = cell_pass_glue pass11 pass12;;
let pass = cell_pass_glue pass1 pass2;;



(*************************************************************************)
(* Test 2: f(x) = x * x * x - atn (x * x); f(x) < 0 when 0.1 <= x <= 0.8 *)
(*************************************************************************)

let x_tm = mk_float 1 49;;
let z_tm = mk_float 8 49;;

let domain_th = mk_center_domain 5 x_tm z_tm;;
let _, y_tm, _, _ = dest_cell_domain (concl domain_th);;

let domain1_th = mk_center_domain 5 x_tm y_tm;;
let _, y1_tm, _, _ = dest_cell_domain (concl domain1_th);;

let domain11_th = mk_center_domain 5 x_tm y1_tm;;
let _, y11_tm, _, _ = dest_cell_domain (concl domain11_th);;

let domain111_th = mk_center_domain 5 x_tm y11_tm;;
let domain112_th = mk_center_domain 5 y11_tm y1_tm;;

let domain12_th = mk_center_domain 5 y1_tm y_tm;;
let domain2_th = mk_center_domain 5 y_tm z_tm;;

let eval_sqr pp domain_th =
  let x_th = eval_taylor_x domain_th in
    taylor_interval_mul pp x_th x_th;;


let eval_f pp domain_th =
  let int_x = eval_taylor_x domain_th in
  let int_x2 = taylor_interval_mul pp int_x int_x in
  let int_x3 = taylor_interval_mul pp int_x int_x2 in
  let atn_x2 = taylor_interval_compose pp eval_taylor_atn eval_sqr domain_th in
    taylor_interval_sub pp int_x3 atn_x2;;

let int_th = eval_f 5 domain_th;;
let int_th111 = eval_f 5 domain111_th;;
let int_th112 = eval_f 5 domain112_th;;
let int_th12 = eval_f 5 domain12_th;;
let int_th2 = eval_f 5 domain2_th;;

taylor_cell_pass 5 int_th;;
let pass111 = taylor_cell_pass 5 int_th111;;
let pass112 = taylor_cell_pass 5 int_th112;;
let pass12 = taylor_cell_pass 5 int_th12;;
let pass2 = taylor_cell_pass 5 int_th2;;

let pass11 = cell_pass_glue pass111 pass112;;
let pass1 = cell_pass_glue pass11 pass12;;
let pass = cell_pass_glue pass1 pass2;;

let th0 = REWRITE_RULE[cell_pass; interval_arith] pass;;
CONV_RULE (DEPTH_CONV FLOAT_TO_NUM_CONV) th0;;



(********************************************)
(* Test 3:                                                             *)
(* !x. &2 <= x /\ x <= #2.52                                           *)
(*   ==> arclength (&2) x (&2) - arclength (#2.52) x (&2) >= #0.073`   *)

let eval_arc1 pp domain_th =
  let two = eval_taylor_const domain_th two_float and
      ( * ) = taylor_interval_mul pp and
      (+) = taylor_interval_add pp and
      (-) = taylor_interval_sub pp in
  let b = (fun pp domain_th ->
	     let x = eval_taylor_x domain_th in
	       two * two * x) in
  let a = (fun pp domain_th -> 
	     let x = eval_taylor_x domain_th in
	     let rhs = taylor_interval_compose pp eval_taylor_inv b domain_th in
	       (two * two + x * x - two * two) * rhs) in
    taylor_interval_compose pp eval_taylor_acs a domain_th;;


let h_float = mk_float 252 48;;
let eval_arc2 pp domain_th =
  let two = eval_taylor_const domain_th two_float and
      h = eval_taylor_const domain_th h_float and
      ( * ) = taylor_interval_mul pp and
      (+) = taylor_interval_add pp and
      (-) = taylor_interval_sub pp in
  let b = (fun pp domain_th ->
	     let x = eval_taylor_x domain_th in
	       two * h * x) in
  let a = (fun pp domain_th -> 
	     let x = eval_taylor_x domain_th in
	     let rhs = taylor_interval_compose pp eval_taylor_inv b domain_th in
	       (h * h + x * x - two * two) * rhs) in
    taylor_interval_compose pp eval_taylor_acs a domain_th;;


let t_float = mk_float 73 47;;
let eval_f pp domain_th =
  let t = eval_taylor_const domain_th t_float and
      (+) = taylor_interval_add pp and
      (-) = taylor_interval_sub pp and
      arc1 = eval_arc1 pp domain_th and
      arc2 = eval_arc2 pp domain_th in
    t + arc2 - arc1;;

let pp = 6;;  
let domain0 = mk_center_domain pp two_float h_float;;

let _, y1_tm, _, _ = dest_cell_domain (concl domain0);;
let domain1 = mk_center_domain pp two_float y1_tm and
    domain2 = mk_center_domain pp y1_tm h_float;;

let _, y12_tm, _, _ = dest_cell_domain (concl domain2);;
let domain21 = mk_center_domain pp y1_tm y12_tm and
    domain22 = mk_center_domain pp y12_tm h_float;;

let _, y122_tm, _, _ = dest_cell_domain (concl domain22);;
let domain221 = mk_center_domain pp y12_tm y122_tm and
    domain222 = mk_center_domain pp y122_tm h_float;;

let int1 = eval_f pp domain1;;
let int21 = eval_f pp domain21;;
let int221 = eval_f pp domain221;;
let int222 = eval_f pp domain222;;

let pass1 = taylor_cell_pass pp int1;;
let pass21 = taylor_cell_pass pp int21;;

let pass221 = taylor_cell_pass pp int221;;
let pass222 = taylor_cell_pass pp int222;;

let pass22 = cell_pass_glue pass221 pass222;;
let pass2 = cell_pass_glue pass21 pass22;;

let pass = cell_pass_glue pass1 pass2;;

let th0 = REWRITE_RULE[cell_pass; interval_arith] pass;;
CONV_RULE (DEPTH_CONV FLOAT_TO_NUM_CONV) th0;;


let ineq_test pp =
  let domain0 = mk_center_domain pp two_float h_float in
  let _, y1_tm, _, _ = dest_cell_domain (concl domain0) in
  let domain1 = mk_center_domain pp two_float y1_tm and
      domain2 = mk_center_domain pp y1_tm h_float in
  let _, y12_tm, _, _ = dest_cell_domain (concl domain2) in
  let domain21 = mk_center_domain pp y1_tm y12_tm and
      domain22 = mk_center_domain pp y12_tm h_float in
  let _, y122_tm, _, _ = dest_cell_domain (concl domain22) in
  let domain221 = mk_center_domain pp y12_tm y122_tm and
      domain222 = mk_center_domain pp y122_tm h_float in

  let int1 = eval_f pp domain1 and
      int21 = eval_f pp domain21 and
      int221 = eval_f pp domain221 and
      int222 = eval_f pp domain222 in

  let pass1 = taylor_cell_pass pp int1 in
  let pass21 = taylor_cell_pass pp int21 in
  let pass221 = taylor_cell_pass pp int221 in
  let pass222 = taylor_cell_pass pp int222 in
  let pass22 = cell_pass_glue pass221 pass222 in
  let pass2 = cell_pass_glue pass21 pass22 in
    cell_pass_glue pass1 pass2;;


ineq_test 5;;
ineq_test 6;;
ineq_test 7;;
ineq_test 10;;

(* 10: 9.768 *)
test 5 ineq_test 6;;
(* 10: 11.012 *)
test 5 ineq_test 7;;
(* 10: 16.381 *)
test 5 ineq_test 10;;




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

let x_tm = mk_float 20 48;;
let z_tm = mk_float 99 48;;

let domain_th = mk_center_domain 5 x_tm z_tm;;
test 1000 (mk_center_domain 5 x_tm) z_tm;;

let int_x = eval_taylor_x domain_th;;
let int_atn = eval_taylor_atn 5 domain_th;;
let int_x2 = taylor_interval_mul 5 int_x int_x;;
let int_x_atn = taylor_interval_mul 5 int_x int_atn;;

let int_th1 = taylor_interval_sub 5 int_x2 int_x;;
let int_th2 = taylor_interval_sub 5 int_x_atn int_x;;

eval_taylor_f_bounds 5 int_th1;;
eval_taylor_f_bounds 5 int_th2;;

let x_th =
  let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th) in
    (mk_bounded_on_int o ASSUME) (mk_interval x_var_real (mk_pair (x_tm, z_tm)));;


let b_x2 = bounded_on_int_mul 5 x_th x_th;;
bounded_on_int_sub 5 b_x2 x_th;;






(*******************************************)
(* arclength #2.52 #2.52 (&2) >= #0.816 *)
(* arclength (&2) #2.52 (&2) - arclength #2.52 #2.52 (&2) >= #0.073 *)
(* arclength a b c = acs ((a * a + b * b - c * c) / (&2 * a * b)) *)

let two_tm = mk_float 2 50;;
let h_tm = mk_float 252 48;;
let pp = 6;;


let arclength_int pp a b c = 
  let ( * ) = float_interval_mul pp and
      (+) = float_interval_add pp and
      (/) = float_interval_div pp and
      (-) = float_interval_sub pp and
      acs = float_interval_acs pp in
    acs ((a * a + b * b - c * c) / (two_interval * a * b));;

let result =
  let (!) = mk_const_interval and
      (-) = float_interval_sub pp in
    arclength_int pp !h_tm !h_tm !two_tm,
  arclength_int pp !two_tm !h_tm !two_tm - arclength_int pp !h_tm !h_tm !two_tm;;
  

(CONV_RULE (DEPTH_CONV FLOAT_TO_NUM_CONV)) (fst result);;
(CONV_RULE (DEPTH_CONV FLOAT_TO_NUM_CONV)) (snd result);;