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;;