needs "../formal_lp/arith/float.hl";;
needs "../formal_lp/arith/arith_hash_rat.hl";;



let rec float_interval_test pp expr x_th =
  if is_var expr then
    x_th
  else
    let ltm, r_tm = dest_comb expr in
      if is_comb ltm then
	let op, l_tm = dest_comb ltm in
	let l_th = float_interval_test pp l_tm x_th in
	let r_th = float_interval_test pp r_tm x_th in
	  if op = plus_op_real then
	    float_interval_add pp l_th r_th
	  else if op = mul_op_real then
	    float_interval_mul pp l_th r_th
	  else if op = div_op_real then
	    float_interval_div pp l_th r_th
	  else if op = minus_op_real then
	    float_interval_sub pp l_th r_th
	  else
	    failwith ("Unknown operation: " ^ (fst o dest_const) op)
      else
	if ltm = neg_op_real then
	  let r_th = float_interval_test pp r_tm x_th in
	    float_interval_neg r_th
	else
	  mk_float_interval_num (dest_numeral r_tm);;


let rec rat_arith_test expr x_th =
  if is_var expr then
    x_th
  else
    let ltm, r_tm = dest_comb expr in
      if is_comb ltm then
	let op, l_tm = dest_comb ltm in
	let l_th = rat_arith_test l_tm x_th in
	let r_th = rat_arith_test r_tm x_th in
	let th0 = MK_COMB (AP_TERM op l_th, r_th) in
	let rtm = rand(concl th0) in
	let th1 =
	  if op = plus_op_real then
	    my_real_rat_add_conv rtm
	  else if op = mul_op_real then
	    my_real_rat_mul_conv rtm
	  else if op = div_op_real then
	    my_real_rat_div_conv rtm
	  else
	    failwith ("Unknown operation: " ^ (fst o dest_const) op) in
	  TRANS th0 th1
      else
	AP_TERM ltm (NUMERAL_TO_NUM_CONV r_tm);;



  

let expr1 = `&1 + x + x*x / &2 + x*x*x / &6 + x*x*x*x / &24 + x*x*x*x*x / &120 + x*x*x*x*x*x / &720`;;
let expr2 = `&1 + x*(&1 + x*(&1 / &2 + x*(&1 / &6 + x*(&1 / &24 + x*(&1 / &120 + x / &720)))))`;;
let expr3 = `&1 + x*x*(--(&1 / &2) + x*x*(&1 / &24 + x*x*(--(&1 / &720) + x*x*(&1 / &40320 + x*x*(--(&1 / &3628800) + x*x / &479001600)))))`;;
let eps_expr = `x*x*x*x*x*x*x*x*x*x*x*x*x*x / &87178291200`;;


let f x =
  1. +. x*.x*.(~-. 0.5 +. x*.x*.(1. /. 24. +. x*.x*.(~-.1. /. 720. +. x*.x*.(1. /. 40320. +. x*.x*.(~-. 1. /. 3628800. +. x*.x /. 479001600.)))));;

f 1.230959417;;



let int1 = mk_float_interval_small_num 1 and
    int3 = mk_float_interval_small_num 3 and
    int234451 = mk_float_interval_small_num 234451 and
    int1234567 = mk_float_interval_small_num 1234567;;


let pp = 10;;

let x1_th = int1;;
let x2_th = float_interval_div pp int1 int3;;
let x3_th = float_interval_div pp int234451 int1234567;;

let x_th = 
  let n = mk_float_interval_num (Num.num_of_string "1230959417") in
  let d = mk_float_interval_num (Num.num_of_string "1000000000") in
    float_interval_div pp n d;;

(* pp = 11:
   10, min_exp = 20: 0.180
   100, min_exp = 20: 0.144 *)
(* pp = 6:
   100, min_exp = 20: 0.116 *)
(* pp = 5:
   8, min_exp = 20: 0.060 *)

float_interval_test pp expr3 x_th;;
test 1 (float_interval_test pp expr3) x_th;;

float_interval_test pp eps_expr x_th;;

let th1 = float_interval_test pp expr3 x_th;;
let th2 = float_interval_sub pp th1 x2_th;;


float_interval_test pp expr3 x1_th;;
float_interval_test pp expr3 x2_th;;
float_interval_test pp expr3 x3_th;;


let y1_th = REPLACE_NUMERALS_CONV `&1`;;
let y2_th = REPLACE_NUMERALS_CONV `&1 / &3`;;
let y3_th = REPLACE_NUMERALS_CONV `&234451 / &1234567`;;

rat_arith_test expr1 y1_th;;
float_interval_test pp expr1 x1_th;;

rat_arith_test expr2 y1_th;;
float_interval_test pp expr2 x1_th;;

rat_arith_test expr1 y2_th;;
float_interval_test pp expr1 x2_th;;

rat_arith_test expr2 y2_th;;
float_interval_test pp expr2 x2_th;;

rat_arith_test expr1 y3_th;;
float_interval_test pp expr1 x3_th;;

rat_arith_test expr2 y3_th;;
float_interval_test pp expr2 x3_th;;




(* 4: 0.624 *)
test 100 (rat_arith_test expr1) y1_th;;
(* 4: 2.028 *)
test 100 (float_interval_test pp expr1) x1_th;;

(* 4: 0.496 *)
test 100 (rat_arith_test expr2) y1_th;;
(* 4: 1.756 *)
test 100 (float_interval_test pp expr2) x1_th;;

(* 4: 1.012 *)
test 100 (rat_arith_test expr1) y2_th;;
(* 4: 4.252 *)
test 100 (float_interval_test pp expr1) x2_th;;

(* 4: 0.760 *)
test 100 (rat_arith_test expr2) y2_th;;
(* 4: 2.460 *)
test 100 (float_interval_test pp expr2) x2_th;;

(* 4: 5.732 *)
test 10 (rat_arith_test expr1) y3_th;;
(* 4: 0.468 *)
test 10 (float_interval_test pp expr1) x3_th;;

(* 4: 1.504 *)
test 100 (rat_arith_test expr2) y3_th;;
(* 4: 0.244 *)
test 100 (float_interval_test pp expr2) x3_th;;