needs "arith/interval_arith.hl";;
needs "misc/vars.hl";;

module Constant_approximations = struct

open Interval_arith;;
open Misc_vars;;


let EPS_TO_INTERVAL = 
prove(`abs (f - x) <= e <=> interval_arith f (x - e, x + e)`,

   REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;
let acs3_lo, acs3_hi =
  let verify = MATCH_MP REAL_LT_IMP_LE o fst o M_verifier_main.verify_ineq M_verifier_main.default_params 10 in
    verify `#1.230959417 < acs (&1 / &3)`, verify `acs(&1 / &3) < #1.230959418`;;
  
let mk_decimal a b = 
  let tm = mk_comb(mk_comb(`DECIMAL`, mk_numeral (Num.abs_num a)), mk_numeral b) in
    if (a </ Int 0) then
      mk_comb (neg_op_real, tm)
    else
      tm;;


let approx_interval th precision =
  let th' = CONV_RULE (RAND_CONV (REWRITE_CONV[DECIMAL] THENC REAL_RAT_REDUCE_CONV)) th in
  let lo_tm, hi_tm = dest_pair (rand(concl th')) in
  let lo, hi = rat_of_term lo_tm, rat_of_term hi_tm in
  let m = (Int 10 **/ Int precision) in
  let lo_bound = floor_num (lo */ m) in
  let hi_bound = ceiling_num (hi */ m) in
  let conv = EQT_ELIM o REAL_RAT_LE_CONV in
  let lo_th = conv (mk_binop le_op_real (mk_decimal lo_bound m) lo_tm) in
  let hi_th = conv (mk_binop le_op_real hi_tm (mk_decimal hi_bound m)) in
  let th0 = CONJ (CONJ lo_th hi_th) th' in
    MATCH_MP APPROX_INTERVAL th0;;



(************************************)
(* Square root *)

let INTERVAL_SQRT = 
prove(`interval_arith x (a, b) /\
			    (c * c <= a /\ b <= d * d) ==>
			    interval_arith (sqrt x) (abs c, abs d)`,

   REWRITE_TAC[interval_arith] THEN REPEAT STRIP_TAC THENL
     [
       MATCH_MP_TAC REAL_LE_RSQRT THEN
	 MATCH_MP_TAC REAL_LE_TRANS THEN
	 EXISTS_TAC `a:real` THEN
	 ASM_REWRITE_TAC[REAL_ARITH `abs a pow 2 = a * a`];
       MATCH_MP_TAC REAL_LE_LSQRT THEN
	 ASM_REWRITE_TAC[REAL_ARITH `abs d pow 2 = d * d`; REAL_ABS_POS] THEN
	 CONJ_TAC THENL
	 [
	   MATCH_MP_TAC REAL_LE_TRANS THEN
	     EXISTS_TAC `a:real` THEN
	     ASM_REWRITE_TAC[] THEN
	     MATCH_MP_TAC REAL_LE_TRANS THEN
	     EXISTS_TAC `c * c:real` THEN
	     ASM_REWRITE_TAC[REAL_LE_SQUARE];
	   MATCH_MP_TAC REAL_LE_TRANS THEN
	     EXISTS_TAC `b:real` THEN
	     ASM_REWRITE_TAC[]
	 ]
     ]);;


let interval_sqrt th precision =
  let th' = CONV_RULE (REWRITE_CONV[DECIMAL] THENC REAL_RAT_REDUCE_CONV) th in
  let lo, hi = dest_pair(rand(concl th')) in
  let x_lo, x_hi = float_of_num (rat_of_term lo), float_of_num (rat_of_term hi) in
  let lo_sqrt, hi_sqrt = Pervasives.sqrt x_lo, Pervasives.sqrt x_hi in
  let m = 10.0 ** float_of_int precision in
  let hack n = num_of_string (Int64.to_string (Int64.of_float n)) in
  let sqrt_lo_num, sqrt_hi_num = hack (floor (lo_sqrt *. m)), hack (ceil (hi_sqrt *. m)) in
  let m_num = Int 10 **/ Int precision in
  let x_lo_tm = mk_decimal sqrt_lo_num m_num in
  let x_hi_tm = mk_decimal sqrt_hi_num m_num in
  let conv = EQT_ELIM o REAL_RAT_REDUCE_CONV in
  let lo_th = conv (mk_binop le_op_real (mk_binop mul_op_real x_lo_tm x_lo_tm) lo) in
  let hi_th = conv (mk_binop le_op_real hi (mk_binop mul_op_real x_hi_tm x_hi_tm)) in
  let th0 = CONJ th' (CONJ lo_th hi_th) in
    (CONV_RULE REAL_RAT_REDUCE_CONV) (MATCH_MP INTERVAL_SQRT th0);;



(************************************)
(* Arithmetic of intervals *)

let INTERVAL_ADD = 
prove(`interval_arith x (a, b) /\ interval_arith y (c, d)
			   ==> interval_arith (x + y) (a + c, b + d)`,

   REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;

let INTERVAL_SUB = 
prove(`interval_arith x (a, b) /\ interval_arith y (c, d)
			  ==> interval_arith (x - y) (a - d, b - c)`,

   REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;

let INTERVAL_NEG = 
prove(`interval_arith x (a, b) ==>
			   interval_arith (--x) (--b, --a)`,

   REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;


let INTERVAL_INV = 
prove(`interval_arith x (a, b) /\ (&0 < a \/ b < &0)
			   ==> interval_arith (inv x) (inv b, inv a)`,

   REWRITE_TAC[interval_arith] THEN
     STRIP_TAC THENL
     [
       CONJ_TAC THEN MATCH_MP_TAC REAL_LE_INV2 THEN 
	 ASM_REWRITE_TAC[] THEN
	 REPEAT (POP_ASSUM MP_TAC) THEN
	 REAL_ARITH_TAC;
       ALL_TAC
     ] THEN 
     ONCE_REWRITE_TAC[REAL_ARITH `a <= b <=> --b <= --a`] THEN
     REWRITE_TAC[GSYM REAL_INV_NEG] THEN
     CONJ_TAC THEN MATCH_MP_TAC REAL_LE_INV2 THEN
     REPEAT (POP_ASSUM MP_TAC) THEN
     REAL_ARITH_TAC);;

let INTERVAL_INV_POS = 
prove(`interval_arith x (a,b) /\ &0 < a
			       ==> interval_arith (inv x) (inv b, inv a)`,

			     SIMP_TAC[INTERVAL_INV]);;

let INTERVAL_INV_NEG = 
prove(`interval_arith x (a,b) /\ b < &0
			       ==> interval_arith (inv x) (inv b, inv a)`,

			     SIMP_TAC[INTERVAL_INV]);;


let INTERVAL_MUL_lemma = 
prove(`!x y a b c d. interval_arith x (a, b) /\ interval_arith y (c, d) /\ x <= y
				 ==> x * y <= max (max (a * c) (a * d)) (max (b * c) (b * d))`,

   REPEAT GEN_TAC THEN
     REWRITE_TAC[interval_arith] THEN DISCH_TAC THEN
     ABBREV_TAC `t = max (max (a * c) (a * d)) (max (b * c) (b * d))` THEN
     SUBGOAL_THEN `a * c <= t /\ a * d <= t /\ b * c <= t /\ b * d <= t:real` ASSUME_TAC THENL
     [
       EXPAND_TAC "t" THEN
	 REAL_ARITH_TAC;
       ALL_TAC
     ] THEN
     
     DISJ_CASES_TAC (REAL_ARITH `&0 <= x \/ x <= &0`) THENL
     [
       MATCH_MP_TAC REAL_LE_TRANS THEN
	 EXISTS_TAC `b * d:real` THEN
	 ASM_REWRITE_TAC[] THEN
	 MATCH_MP_TAC REAL_LE_MUL2 THEN
	 ASM_REWRITE_TAC[] THEN
	 MATCH_MP_TAC REAL_LE_TRANS THEN
	 EXISTS_TAC `x:real` THEN
	 ASM_REWRITE_TAC[];
       ALL_TAC
     ] THEN

     DISJ_CASES_TAC (REAL_ARITH `&0 <= b \/ b <= &0`) THENL
     [
       DISJ_CASES_TAC (REAL_ARITH `&0 <= y \/ y <= &0`) THENL
	 [
	   MATCH_MP_TAC REAL_LE_TRANS THEN
	     EXISTS_TAC `&0` THEN
	     CONJ_TAC THENL
	     [
	       ONCE_REWRITE_TAC[REAL_ARITH `&0 = &0 * y`] THEN
		 MATCH_MP_TAC REAL_LE_RMUL THEN
		 ASM_REWRITE_TAC[];
	       ALL_TAC
	     ] THEN
	     
	     MATCH_MP_TAC REAL_LE_TRANS THEN
	     EXISTS_TAC `b * d:real` THEN
	     ASM_REWRITE_TAC[] THEN
	     MATCH_MP_TAC REAL_LE_MUL THEN
	     ASM_REWRITE_TAC[] THEN
	     MATCH_MP_TAC REAL_LE_TRANS THEN
	     EXISTS_TAC `y:real` THEN
	     ASM_REWRITE_TAC[];
	   ALL_TAC
	 ] THEN

	 MATCH_MP_TAC REAL_LE_TRANS THEN
	 EXISTS_TAC `a * c:real` THEN
	 ASM_REWRITE_TAC[] THEN
	 ONCE_REWRITE_TAC[GSYM REAL_NEG_MUL2] THEN
	 MATCH_MP_TAC REAL_LE_MUL2 THEN
	 ASM_REWRITE_TAC[REAL_LE_NEG; REAL_NEG_GE0];
       ALL_TAC
     ] THEN

     DISJ_CASES_TAC (REAL_ARITH `&0 <= c \/ c <= &0`) THENL
     [
       MATCH_MP_TAC REAL_LE_TRANS THEN
	 EXISTS_TAC `b * c:real` THEN
	 ASM_REWRITE_TAC[] THEN
	 ONCE_REWRITE_TAC[REAL_ARITH `x * y <= b * c <=> (--b) * c <= (--x) * y`] THEN
	 MATCH_MP_TAC REAL_LE_MUL2 THEN
	 ASM_REWRITE_TAC[REAL_LE_NEG; REAL_NEG_GE0];
       ALL_TAC
     ] THEN

     DISJ_CASES_TAC (REAL_ARITH `&0 <= y \/ y <= &0`) THENL
     [
       MATCH_MP_TAC REAL_LE_TRANS THEN
	 EXISTS_TAC `&0` THEN
	 CONJ_TAC THENL
	 [
	   ONCE_REWRITE_TAC[REAL_ARITH `&0 = &0 * y`] THEN
	     MATCH_MP_TAC REAL_LE_RMUL THEN
	     ASM_REWRITE_TAC[];
	   ALL_TAC
	 ] THEN
	 
	 MATCH_MP_TAC REAL_LE_TRANS THEN
	 EXISTS_TAC `a * c:real` THEN
	 ASM_REWRITE_TAC[] THEN
	 ONCE_REWRITE_TAC[GSYM REAL_NEG_MUL2] THEN
	 MATCH_MP_TAC REAL_LE_MUL THEN
	 ASM_REWRITE_TAC[REAL_NEG_GE0] THEN
	 MATCH_MP_TAC REAL_LE_TRANS THEN
	 EXISTS_TAC `x:real` THEN
	 ASM_REWRITE_TAC[];
       ALL_TAC
     ] THEN

     MATCH_MP_TAC REAL_LE_TRANS THEN
     EXISTS_TAC `a * c:real` THEN
     ASM_REWRITE_TAC[] THEN
     ONCE_REWRITE_TAC[GSYM REAL_NEG_MUL2] THEN
     MATCH_MP_TAC REAL_LE_MUL2 THEN
     ASM_REWRITE_TAC[REAL_LE_NEG; REAL_NEG_GE0]);;


let INTERVAL_MUL_lemma2 = 
prove(`!x y a b c d. interval_arith x (a,b) /\ interval_arith y (c,d)
				  ==> x * y <= max (max (a * c) (a * d)) (max (b * c) (b * d))`,

   REPEAT STRIP_TAC THEN
     DISJ_CASES_TAC (REAL_ARITH `x <= y \/ y <= x:real`) THENL
     [
       MATCH_MP_TAC INTERVAL_MUL_lemma THEN
	 ASM_REWRITE_TAC[];
       ALL_TAC
     ] THEN

     MP_TAC (SPECL [`y:real`; `x:real`; `c:real`; `d:real`; `a:real`; `b:real`] INTERVAL_MUL_lemma) THEN
     ASM_REWRITE_TAC[] THEN
     REAL_ARITH_TAC);;
     



let INTERVAL_MUL = 
prove(`interval_arith x (a, b) /\ interval_arith y (c, d)
			   ==> interval_arith (x * y)
			   (min (min (a * c) (a * d)) (min (b * c) (b * d)),
			    max (max (a * c) (a * d)) (max (b * c) (b * d)))`,

   DISCH_TAC THEN REWRITE_TAC[interval_arith] THEN
     ASM_SIMP_TAC[INTERVAL_MUL_lemma2] THEN
     MP_TAC (SPECL[`--x:real`; `y:real`; `--b:real`; `--a:real`; `c:real`; `d:real`] INTERVAL_MUL_lemma2) THEN
     ASM_SIMP_TAC[INTERVAL_NEG] THEN
     REAL_ARITH_TAC);;
     
   

(**************************************)
let const_interval tm = SPEC tm CONST_INTERVAL;;

let interval_neg th = MATCH_MP INTERVAL_NEG th;;

let interval_add th1 th2 =
  let th0 = MATCH_MP INTERVAL_ADD (CONJ th1 th2) in
    (CONV_RULE (RAND_CONV REAL_RAT_REDUCE_CONV)) th0;;

let interval_sub th1 th2 =
  let th0 = MATCH_MP INTERVAL_SUB (CONJ th1 th2) in
    (CONV_RULE (RAND_CONV REAL_RAT_REDUCE_CONV)) th0;;

let interval_mul th1 th2 =
  let th0 = MATCH_MP INTERVAL_MUL (CONJ th1 th2) in
    (CONV_RULE (RAND_CONV REAL_RAT_REDUCE_CONV)) th0;;

let interval_inv th =
  let lt_op_real = `(<):real->real->bool` in
  let lo_tm, hi_tm = dest_pair(rand(concl th)) in
  let lo_ineq = REAL_RAT_LT_CONV (mk_binop lt_op_real `&0` lo_tm) in
    if (rand(concl lo_ineq) = `T`) then
      let th0 = CONJ th (EQT_ELIM lo_ineq) in
	(CONV_RULE (RAND_CONV REAL_RAT_REDUCE_CONV)) (MATCH_MP INTERVAL_INV_POS th0)
    else 
      let hi_ineq = REAL_RAT_LT_CONV (mk_binop lt_op_real hi_tm `&0`) in
	if (rand(concl hi_ineq) = `T`) then
	  let th0 = CONJ th (EQT_ELIM hi_ineq) in
	    (CONV_RULE (RAND_CONV REAL_RAT_REDUCE_CONV)) (MATCH_MP INTERVAL_INV_NEG th0)
	else failwith("interval_inv: 0 is inside interval");;

let interval_div th1 th2 =
  let th2' = interval_inv th2 in
    ONCE_REWRITE_RULE[GSYM real_div] (interval_mul th1 th2');;



(*************************)
let acs3_interval = REWRITE_RULE[GSYM interval_arith] (CONJ acs3_lo acs3_hi);;

let pi_interval = 
prove(`interval_arith pi (#3.141592653, #3.141592654)`,

   REWRITE_TAC[interval_arith] THEN
     MP_TAC PI_APPROX_32 THEN REAL_ARITH_TAC);;
let tgt_interval = 
prove(`interval_arith tgt (#1.541, #1.541)`,

   REWRITE_TAC[Tame_defs.tgt; interval_arith; REAL_LE_REFL]);;

let interval_table = Hashtbl.create 10;;
let add_interval th = Hashtbl.add interval_table ((rand o rator o concl) th) th;;


let rec create_interval tm =
  if Hashtbl.mem interval_table tm then
    Hashtbl.find interval_table tm
  else
    if (is_ratconst tm) then
      const_interval tm
    else if (is_binop add_op_real tm) then
      let lhs, rhs = dest_binop add_op_real tm in
	interval_add (create_interval lhs) (create_interval rhs)
    else if (is_binop sub_op_real tm) then
      let lhs, rhs = dest_binop sub_op_real tm in
	interval_sub (create_interval lhs) (create_interval rhs)
    else if (is_binop mul_op_real tm) then
      let lhs, rhs = dest_binop mul_op_real tm in
	interval_mul (create_interval lhs) (create_interval rhs)
    else if (is_binop div_op_real tm) then
      let lhs, rhs = dest_binop div_op_real tm in
	interval_div (create_interval lhs) (create_interval rhs)
    else if (is_comb tm) then
      let ltm, rtm = dest_comb tm in
	if (ltm = inv_op_real) then
	  interval_inv (create_interval rtm)
	else if (ltm = neg_op_real) then
	  interval_neg (create_interval rtm)
	else failwith "create_interval: unknown unary operation"
    else
      failwith "create_interval: unexpected term";;


open Sphere;;

add_interval pi_interval;;
add_interval acs3_interval;;
add_interval tgt_interval;;
add_interval (REWRITE_RULE[GSYM sqrt8] (interval_sqrt (const_interval `&8`) 9));;
add_interval (REWRITE_RULE[GSYM sol0] (create_interval `&3 * acs(&1 / &3) - pi`));;
add_interval (create_interval `sol0 / pi`);;

let rho218 = new_definition `rho218 = rho(#2.18)`;;
let rho218_def = (REWRITE_CONV[rho218; rho; ly; interp; GSYM Tame_general.sol0_over_pi_EQ_const1] THENC
   REAL_RAT_REDUCE_CONV) `rho218`;;

let rho218_interval = REWRITE_RULE[SYM rho218_def] (create_interval(rand(concl rho218_def)));;
add_interval rho218_interval;;


end;;