(* ------------------------------------------------------------------ *)
(* Author and Copyright: Thomas C. Hales                              *)
(* License: GPL http://www.gnu.org/copyleft/gpl.html                  *)
(* Project: FLYSPECK http://www.math.pitt.edu/~thales/flyspeck/       *)
(* ------------------------------------------------------------------ *)



prioritize_real();;

let add_test,test = new_test_suite();;


let twopow =
  new_definition(
        `twopow x = if (?n. (x = (int_of_num n)))
        then ((&2) pow (nabs x))
        else inv((&2) pow (nabs x))`);;


let float =
  new_definition(
                  `float x n = (real_of_int x)*(twopow n)`);;


let interval =
  new_definition(
                   `interval x f eps = ((abs (x-f)) <= eps)`);;

(*--------------------------------------------------------------------*)

let mk_interval a b ex =
   mk_comb(mk_comb (mk_comb (`interval`,a),b),ex);;

add_test("mk_interval",
   mk_interval `#3` `#4` `#1` = `interval #3 #4 #1`);;

let dest_interval intv =
   let (h1,ex) = dest_comb intv in
   let (h2,f) = dest_comb h1 in
   let (h3,a) = dest_comb h2 in
   let _ = assert(h3 = `interval`) in
   (a,f,ex);;

add_test("dest_interval",
   let a = `#3` and b = `#4` and c = `#1` in
   dest_interval (mk_interval a b c) = (a,b,c));;

(*--------------------------------------------------------------------*)

let (dest_int:term-> Num.num) =
  fun b ->
  let dest_pos_int a =
    let (op,nat) = dest_comb a in
    if (fst (dest_const op) = "int_of_num") then (dest_numeral nat)
      else fail() in
    let (op',u) = (dest_comb b) in
    try (if (fst (dest_const op') = "int_neg") then
           minus_num (dest_pos_int u) else
             dest_pos_int b) with
        Failure _ -> failwith "dest_int ";;


let (mk_int:Num.num -> term) =
  fun a ->
    let sgn = Num.sign_num a in
    let abv = Num.abs_num a in
    let r = mk_comb(`&:`,mk_numeral abv) in
    try (if (sgn<0) then mk_comb (`--:`,r) else r) with
        Failure _ -> failwith ("dest_int "^(string_of_num a));;

add_test("mk_int",
   (mk_int (Int (-1443)) = `--: (&:1443)`) &&
   (mk_int (Int 37) = `(&:37)`));;

(* ------------------------------------------------------------------ *)

let (split_ratio:Num.num -> Num.num*Num.num) =
  function
    (Ratio r) -> (Big_int (Ratio.numerator_ratio r)),
         (Big_int (Ratio.denominator_ratio r))|
    u -> (u,(Int 1));;

add_test("split_ratio",
   let (a,b) = split_ratio ((Int 4)//(Int 20)) in
   (a =/ (Int 1)) && (b =/ (Int 5)));;

(* ------------------------------------------------------------------ *)

(* break nonzero int r into a*(C**b) with a prime to C . *)
let (factor_C:int -> Num.num -> Num.num*Num.num) =
  function c ->
  let intC = (Int c) in
  let rec divC (a,b) =
    if ((Int 0) =/ mod_num a intC) then (divC (a//intC,b+/(Int 1)))
      else (a,b) in
  function r->
  if ((Num.is_integer_num r)&& not((Num.sign_num r) = 0)) then
    divC (r,(Int 0)) else failwith "factor_C";;

add_test("factor_C",
   (factor_C 2 (Int (4096+32)) = (Int 129,Int 5)) &&
   (factor_C 10 (Int (5000)) = (Int 5,Int 3)) &&
   (cannot (factor_C 2) ((Int 50)//(Int 3))));;

(*--------------------------------------------------------------------*)

let (dest_float:term -> Num.num) =
  fun f ->
    let (a,b) = dest_binop `float` f in
    (dest_int a)*/ ((Int 2) **/ (dest_int b));;

add_test("dest_float",
   dest_float `float (&:3) (&:17)` = (Int 393216));;

add_test("dest_float2", (* must express as numeral first *)
   cannot dest_float `float ((&:3)+:(&:1)) (&:17)`);;

(* ------------------------------------------------------------------ *)
(* creates float of the form `float a b` with a odd *)
let (mk_float:Num.num -> term) =
  function r ->
    let (a,b) = split_ratio r in
    let (a',exp_a) = if (a=/(Int 0)) then ((Int 0),(Int 0)) else factor_C 2 a in
    let (b',exp_b) = factor_C 2 b in
    let c = a'//b' in
    if (Num.is_integer_num c) then
          mk_binop `float` (mk_int c) (mk_int (exp_a -/ exp_b))
          else failwith "mk_float";;

add_test("mk_float",
   mk_float (Int (4096+32)) = `float (&:129) (&:5)` &&
   (mk_float (Int 0) = `float (&:0) (&:0)`));;

add_test("mk_float2",  (* throws exception exactly when denom != 2^k *)
   let rtest = fun t -> (t =/ dest_float (mk_float t)) in
   rtest ((Int 3)//(Int 1024)) &&
  (cannot rtest ((Int 1)//(Int 3))));;

add_test("mk_float dest_float",  (* constructs canonical form of float *)
  mk_float (dest_float `float (&:4) (&:3)`) = `float (&:1) (&:5)`);;

(* ------------------------------------------------------------------ *)
(* creates decimal of the form `DECIMAL a b` with a prime to 10 *)
let (mk_pos_decimal:Num.num -> term) =
  function r ->
    let _ = assert (r >=/ (Int 0)) in
    let (a,b) = split_ratio r in
    if (a=/(Int 0)) then `#0` else
    let (a1,exp_a5) = factor_C 5 a in
    let (a2,exp_a2) = factor_C 2 a1 in
    let (b1,exp_b5) = factor_C 5 b in
    let (b2,exp_b2) = factor_C 2 b1 in
    let _ = assert(b2 =/ (Int 1)) in
    let c = end_itlist Num.max_num [exp_b5-/exp_a5;exp_b2-/exp_a2;(Int 0)] in
    let a' = a2*/((Int 2)**/ (c +/ exp_a2 -/ exp_b2))*/
             ((Int 5)**/(c +/ exp_a5 -/ exp_b5)) in
    let b' = (Int 10) **/ c in
    mk_binop `DECIMAL` (mk_numeral a') (mk_numeral b');;

add_test("mk_pos_decimal",
   mk_pos_decimal (Int (5000)) = `#5000` &&
   (mk_pos_decimal ((Int 30)//(Int 40)) = `#0.75`) &&
   (mk_pos_decimal (Int 0) = `#0`) &&
   (mk_pos_decimal ((Int 15)//(Int 25)) = `#0.6`) &&
   (mk_pos_decimal ((Int 25)//(Int 4)) = `#6.25`) &&
   (mk_pos_decimal ((Int 2)//(Int 25)) = `#0.08`));;

let (mk_decimal:Num.num->term) =
  function r ->
  let a = Num.sign_num r in
  let b = mk_pos_decimal (Num.abs_num r) in
  if (a < 0) then (mk_comb (`--.`, b)) else b;;

add_test("mk_decimal",
  (mk_decimal (Int 3) = `#3`) &&
  (mk_decimal (Int (-3)) = `--. (#3)`));;



(*--------------------------------------------------------------------*)

let (dest_decimal:term -> Num.num) =
  fun f ->
    let (a,b) = dest_binop `DECIMAL` f in
    let a1 = dest_numeral a in
    let b1 = dest_numeral b in
        a1//b1;;

add_test("dest_decimal",
   dest_decimal `#3.4` =/ ((Int 34)//(Int 10)));;
add_test("dest_decimal2",
   cannot dest_decimal `--. (#3.4)`);;





(*--------------------------------------------------------------------*)
(*   Properties of integer powers of 2.                               *)
(* ------------------------------------------------------------------ *)



let TWOPOW_POS = 
prove(`!n. (twopow (int_of_num n) = (&2) pow n)`,

        (REWRITE_TAC[twopow])
        THEN GEN_TAC
        THEN COND_CASES_TAC
        THENL [AP_TERM_TAC;ALL_TAC]
        THEN (REWRITE_TAC[NABS_POS])
        THEN (UNDISCH_EL_TAC 0)
        THEN (TAUT_TAC (` ( A    ) ==> (~ A ==> B)`))
        THEN (MESON_TAC[]));;


let TWOPOW_NEG = 
prove(`!n. (twopow (--(int_of_num n)) = inv((&2) pow n))`,

        GEN_TAC
        THEN (REWRITE_TAC[twopow])
        THEN (COND_CASES_TAC THENL [ALL_TAC;REWRITE_TAC[NABS_NEG]])
        THEN (POP_ASSUM CHOOSE_TAC)
        THEN (REWRITE_TAC[NABS_NEG])
        THEN (UNDISCH_EL_TAC 0)
        THEN (REWRITE_TAC[int_eq;int_neg_th;INT_NUM_REAL])
        THEN (REWRITE_TAC[prove (`! u y.((--(real_of_num u) = (real_of_num y))=
                ((real_of_num u) +(real_of_num y) = (&0)))`,REAL_ARITH_TAC)])
        THEN (REWRITE_TAC[REAL_OF_NUM_ADD;REAL_OF_NUM_EQ;ADD_EQ_0])
        THEN (DISCH_TAC)
        THEN (ASM_REWRITE_TAC[real_pow;REAL_INV_1]));;



let TWOPOW_INV = 
prove(`!a. (twopow (--: a) = (inv (twopow a)))`,

  (GEN_TAC)
  THEN ((ASSUME_TAC (SPEC `a:int` INT_REP2)))
  THEN ((POP_ASSUM CHOOSE_TAC))
  THEN ((POP_ASSUM DISJ_CASES_TAC))
  THEN ((ASM_REWRITE_TAC[TWOPOW_POS;TWOPOW_NEG;REAL_INV_INV;INT_NEG_NEG])));;


let INT_REP3 = 
prove(`!a .(?n.( (a = &: n) \/ (a = --: (&: (n+1)))))`,

(GEN_TAC)
THEN ((ASSUME_TAC (SPEC `a:int` INT_REP2)))
THEN ((POP_ASSUM CHOOSE_TAC))
THEN ((DISJ_CASES_TAC (prove (`((a:int) = (&: 0)) \/ ~((a:int) =(&: 0))`, MESON_TAC[]))))
(* cases *)
THENL[ ((EXISTS_TAC `0`)) THEN ((ASM_REWRITE_TAC[]));ALL_TAC]
THEN ((UNDISCH_EL_TAC 0))
THEN ((POP_ASSUM DISJ_CASES_TAC))
THENL [DISCH_TAC THEN ((ASM MESON_TAC)[]);ALL_TAC]
THEN (DISCH_TAC)
THEN ((EXISTS_TAC `PRE n`))
THEN ((DISJ2_TAC))
THEN ((ASM_REWRITE_TAC[INT_EQ_NEG2]))
(*** Changed by JRH, 2006/03/28 to avoid PRE_ELIM_TAC ***)
THEN (FIRST_X_ASSUM(MP_TAC o check(is_neg o concl)))
THEN (ASM_REWRITE_TAC[INT_NEG_EQ_0; INT_OF_NUM_EQ])
THEN (ARITH_TAC));;


let REAL_EQ_INV = 
prove(`!x y. ((x:real = y) <=> (inv(x) = inv (y)))`,

((REPEAT GEN_TAC))
THEN (EQ_TAC)
THENL [((DISCH_TAC THEN (ASM_REWRITE_TAC[])));
 (* branch 2*) ((DISCH_TAC))
THEN ((ONCE_REWRITE_TAC [(GSYM REAL_INV_INV)]))
THEN ((ASM_REWRITE_TAC[]))]);;


let TWOPOW_ADD_1 = 
prove(`!a. (twopow (a +: (&:1)) = twopow (a) *. (twopow (&:1)))`,

EVERY[
  GEN_TAC;
  CHOOSE_TAC (SPEC `a:int` INT_REP3);
  POP_ASSUM DISJ_CASES_TAC
  THENL[
    ASM_REWRITE_TAC[TWOPOW_POS;INT_OF_NUM_ADD;REAL_POW_ADD];
    EVERY[
      ASM_REWRITE_TAC[GSYM INT_OF_NUM_ADD;INT_NEG_ADD;GSYM INT_ADD_ASSOC;INT_ADD_LINV;INT_ADD_RID];
      REWRITE_TAC[GSYM INT_NEG_ADD;INT_OF_NUM_ADD;TWOPOW_NEG;TWOPOW_POS];
      ONCE_REWRITE_TAC[SPEC `(&. 2) pow 1` (GSYM REAL_INV_INV)];
      REWRITE_TAC[GSYM REAL_INV_MUL;GSYM REAL_EQ_INV;REAL_POW_ADD;GSYM REAL_MUL_ASSOC;REAL_POW_1];
      REWRITE_TAC[MATCH_MP REAL_MUL_RINV (REAL_ARITH `~((&. 2) = (&. 0))`); REAL_MUL_RID]
    ]
  ]
]);;


let REAL_INV2 = 
prove(
  `(inv(&. 2)*(&. 2) = (&.1)) /\ ((&. 2)*inv(&. 2) = (&.1))`,

  SUBGOAL_TAC `~((&.2) = (&.0))`
THENL[
  REAL_ARITH_TAC;
  SIMP_TAC[REAL_MUL_RINV;REAL_MUL_LINV]]);;


let TWOPOW_0 = 
prove(`twopow (&: 0) = (&. 1)`,

 (REWRITE_TAC[TWOPOW_POS;real_pow]));;


let TWOPOW_SUB_NUM = 
prove(`!m n.( twopow((&:m) - (&: n)) = twopow((&:m))*. twopow(--: (&:n)))`,

((INDUCT_TAC))
THENL [REWRITE_TAC[INT_SUB_LZERO;REAL_MUL_LID;TWOPOW_0];ALL_TAC]
THEN ((INDUCT_TAC THEN
   ( (ASM_REWRITE_TAC[INT_SUB_RZERO;TWOPOW_0;REAL_MUL_RID;INT_NEG_0;ADD1;GSYM INT_OF_NUM_ADD]))))
THEN ((ASM_REWRITE_TAC [TWOPOW_ADD_1;TWOPOW_INV;prove (`((&:m)+(&:1)) -: ((&:n) +: (&:1)) = ((&:m)-: (&:n))`,INT_ARITH_TAC)]))
THEN ((REWRITE_TAC[REAL_INV_MUL]))
THEN ((ABBREV_TAC `a:real = twopow (&: m)`))
THEN ((ABBREV_TAC `b:real = inv(twopow (&: n))`))
THEN ((REWRITE_TAC[TWOPOW_POS;REAL_POW_1;GSYM REAL_MUL_ASSOC;prove (`!(x:real). ((&.2)*x = x*(&.2))`,REAL_ARITH_TAC)]))
THEN ((REWRITE_TAC[REAL_INV2;REAL_MUL_RID])));;


let TWOPOW_ADD_NUM = 
prove(
  `!m n. (twopow ((&:m) + (&:n)) = twopow((&:m))*. twopow((&:n)))`,


let TWOPOW_ADD_INT = 
prove(
  `!a b. (twopow (a +: b) = twopow(a) *. (twopow(b)))`,

 ((REPEAT GEN_TAC))
THEN ((ASSUME_TAC (SPEC `a:int` INT_REP)))
THEN ((POP_ASSUM CHOOSE_TAC))
THEN ((POP_ASSUM CHOOSE_TAC))
THEN ((ASSUME_TAC (SPEC `b:int` INT_REP)))
THEN ((REPEAT (POP_ASSUM CHOOSE_TAC)))
THEN ((ASM_REWRITE_TAC[]))
THEN ((SUBGOAL_TAC `&: n -: &: m +: &: n' -: &: m' = (&: (n+n')) -: (&: (m+m'))`))
(* branch *)
THENL[ ((REWRITE_TAC[GSYM INT_OF_NUM_ADD]))
THEN ((INT_ARITH_TAC));ALL_TAC]
(* 2nd *)
THEN ((DISCH_TAC))
THEN ((ASM_REWRITE_TAC[TWOPOW_SUB_NUM;TWOPOW_INV;TWOPOW_POS;REAL_POW_ADD;REAL_INV_MUL;GSYM REAL_MUL_ASSOC]))
THEN ((ABBREV_TAC `a':real = inv ((&. 2) pow m)`))
THEN ((ABBREV_TAC `c :real = (&. 2) pow n`))
THEN ((ABBREV_TAC `d :real = (&. 2) pow n'`))
THEN ((ABBREV_TAC `e :real = inv ((&. 2) pow m')`))
THEN (MESON_TAC[REAL_MUL_AC]));;


let TWOPOW_ABS = 
prove(`!a. ||. (twopow a) = (twopow a)`,

(GEN_TAC)
THEN ((CHOOSE_THEN DISJ_CASES_TAC (SPEC `a:int` INT_REP2)))
(* branch *)
THEN ((ASM_REWRITE_TAC[TWOPOW_POS;TWOPOW_NEG;REAL_ABS_POW;REAL_ABS_NUM;REAL_ABS_INV])));;


let REAL_POW_POW = 
prove(
  `!x m n . (x **. m) **. n = x **. (m *| n)`,

((GEN_TAC THEN GEN_TAC THEN INDUCT_TAC))
(* branch *)
THENL[ ((REWRITE_TAC[real_pow;MULT_0]));
(* second branch *)
((REWRITE_TAC[real_pow]))
THEN ((ASM_REWRITE_TAC[ADD1;LEFT_ADD_DISTRIB;REAL_POW_ADD;REAL_MUL_AC;MULT_CLAUSES]))]);;


let INT_POW_POW = INT_OF_REAL_THM REAL_POW_POW;;


let TWOPOW_POW = 
prove(
  `!a n. (twopow a) pow n = twopow (a *: (&: n))`,

((REPEAT GEN_TAC))
THEN ((CHOOSE_THEN DISJ_CASES_TAC (SPEC `a:int` INT_REP2)))
(* branch *)
THEN ((ASM_REWRITE_TAC[TWOPOW_POS;INT_OF_NUM_MUL;
   REAL_POW_POW;TWOPOW_NEG;REAL_POW_INV;INT_OF_NUM_MUL;GSYM INT_NEG_LMUL])));;

(* ------------------------------------------------------------------ *)
(*   Arithmetic operations on float                                   *)
(* ------------------------------------------------------------------ *)

let FLOAT_NEG = 
prove(`!a m. --. (float a m) = float (--: a) m`,

 REPEAT GEN_TAC THEN
 REWRITE_TAC[float;GSYM REAL_MUL_LNEG;int_neg_th]);;




let FLOAT_MUL = 
prove(`!a b m n. (float a m) *. (float b n) = (float (a *: b) (m +: n))`,

((REPEAT GEN_TAC))
THEN ((REWRITE_TAC[float;GSYM REAL_MUL_ASSOC;TWOPOW_ADD_INT;int_mul_th]))
THEN ((MESON_TAC[REAL_MUL_AC])));;


let FLOAT_ADD = 
prove(
  `!a b c m. (float a (m+: (&:c))) +. (float b m)
     = (float ( (&:(2 EXP c))*a +: b) m)`,

((REWRITE_TAC[float;int_add_th;REAL_ADD_RDISTRIB;int_mul_th;TWOPOW_ADD_INT]))
THEN ((REPEAT GEN_TAC))
THEN ((REWRITE_TAC[TWOPOW_POS;INT_NUM_REAL;GSYM REAL_OF_NUM_POW]))
THEN ((MESON_TAC[REAL_MUL_AC])));;


let FLOAT_ADD_EQ = 
prove(
  `!a b m. (float a  m) +. (float b m) =
  (float (a+:b) m)`,

 ((REPEAT GEN_TAC))
THEN ((REWRITE_TAC[REWRITE_RULE[INT_ADD_RID] (SPEC `m:int` (SPEC `0` (SPEC `b:int` (SPEC `a:int` FLOAT_ADD))))]))
THEN ((REWRITE_TAC[EXP;INT_MUL_LID])));;


let FLOAT_ADD_NP = 
prove(
  `!a b m n.  (float b (--:(&: n))) +. (float a (&: m)) = (float a (&: m)) +. (float b (--:(&: n)))`,

(REWRITE_TAC[REAL_ADD_AC]));;


let FLOAT_ADD_PN = 
prove(
  `!a b m n. (float a (&: m)) +. (float b (--(&: n))) = (float ( (&:(2 EXP (m+| n)))*a + b) (--:(&: n)))`,

((REPEAT GEN_TAC))
THEN ((SUBGOAL_TAC `&: m = (--:(&: n)) + (&:(m+n))`))
THENL[ ((REWRITE_TAC[GSYM INT_OF_NUM_ADD]))
THEN ((INT_ARITH_TAC));
(* branch *)
((DISCH_TAC))
THEN ((ASM_REWRITE_TAC[FLOAT_ADD]))]);;


let FLOAT_ADD_PP = 
prove(
  `!a b m n. ((n <=| m) ==>( (float a (&: m)) +. (float b (&: n)) = (float  ((&:(2 EXP (m -| n))) *a + b) (&: n))))`,

((REPEAT GEN_TAC))
THEN (DISCH_TAC)
THEN ((SUBGOAL_TAC `&: m = (&: n) + (&: (m-n))`))
THENL[ ((REWRITE_TAC[INT_OF_NUM_ADD]))
THEN (AP_TERM_TAC)
THEN ((REWRITE_TAC[prove (`!(m:num) n. (n+m-n) = (m-n)+n`,REWRITE_TAC[ADD_AC])]))
THEN ((UNDISCH_EL_TAC 0))
THEN ((MATCH_ACCEPT_TAC(GSYM SUB_ADD)));
(* branch *)
((DISCH_TAC))
THEN ((ASM_REWRITE_TAC[FLOAT_ADD]))]);;


let FLOAT_ADD_PPv2 = 
prove(
  `!a b m n. ((m <| n) ==>( (float a (&: m)) +. (float b (&: n)) = (float  ((&:(2 EXP (n -| m))) *b + a) (&: m))))`,

((REPEAT GEN_TAC))
THEN (DISCH_TAC)
THEN ((H_MATCH_MP (THM (prove(`!m n. m<|n ==> m <=|n`,MESON_TAC[LT_LE]))) (HYP_INT 0)))
THEN ((UNDISCH_EL_TAC 0))
THEN ((SIMP_TAC[GSYM FLOAT_ADD_PP]))
THEN (DISCH_TAC)
THEN ((REWRITE_TAC[REAL_ADD_AC])));;


let FLOAT_ADD_NN = 
prove(
`!a b m n. ((n <=| m) ==>( (float a (--:(&: m))) +. (float b (--:(&: n)))
     = (float  ((&:(2 EXP (m -| n))) *b + a) (--:(&: m)))))`,

((REPEAT GEN_TAC))
THEN (DISCH_TAC)
THEN ((SUBGOAL_TAC `--: (&: n) = --: (&: m) + (&: (m-n))`))
THENL [((UNDISCH_EL_TAC 0))
THEN ((SIMP_TAC [INT_OF_REAL_THM (GSYM REAL_OF_NUM_SUB)]))
THEN (DISCH_TAC)
THEN ((INT_ARITH_TAC));
(*branch*)
((DISCH_TAC))
THEN (ASM_REWRITE_TAC[GSYM FLOAT_ADD;REAL_ADD_AC])]);;


let FLOAT_ADD_NNv2 = 
prove(
`!a b m n. ((m <| n) ==>( (float a (--:(&: m))) +. (float b (--:(&: n)))
     = (float  ((&:(2 EXP (n -| m))) *a + b) (--:(&: n)))))`,

((REPEAT GEN_TAC))
THEN (DISCH_TAC)
THEN (((H_MATCH_MP (THM (prove(`!m n. m<|n ==> m <=|n`,MESON_TAC[LT_LE]))) (HYP_INT 0))))
THEN (((UNDISCH_EL_TAC 0)))
THEN (((SIMP_TAC[GSYM FLOAT_ADD_NN])))
THEN ((DISCH_TAC))
THEN (((REWRITE_TAC[REAL_ADD_AC]))));;


let FLOAT_SUB = 
prove(
  `!a b n m. (float a n) -. (float b m)
     = (float a n) +. (float (--: b) m)`,

REWRITE_TAC[float;int_neg_th;real_sub;REAL_NEG_LMUL]);;


let FLOAT_ABS = 
prove(
  `!a n. ||. (float a n) = (float (||: a) n)`,

(REWRITE_TAC[float;int_abs_th;REAL_ABS_MUL;TWOPOW_ABS]));;



let FLOAT_POW = 
prove(
  `!a n m. (float a n) **. m = (float (a **: m) (n *: (&:m)))`,

(REWRITE_TAC[float;REAL_POW_MUL;int_pow_th;TWOPOW_POW]));;


let INT_SUB = 
prove(
  `!a b. (a -: b) = (a +: (--: b))`,

 (REWRITE_TAC[GSYM INT_SUB_RNEG;INT_NEG_NEG]));;


let INT_ABS_NUM = 
prove(
  `!n. ||: (&: n) = (&: n)`,


let INT_ABS_NEG_NUM = 
prove(
  `!n. ||: (--: (&: n)) = (&: n)`,


let INT_ADD_NEG_NUM = 
prove(`!x y. --: (&: x) +: (&: y) = (&: y) +: (--: (&: x))`,

 (REWRITE_TAC[INT_ADD_AC]));;


let INT_POW_MUL = INT_OF_REAL_THM REAL_POW_MUL;;


let INT_POW_NEG1 = 
prove (
  `!x n. (--: (&: x)) **: n = ((--: (&: 1)) **: n) *: ((&: x) **: n)`,

(REWRITE_TAC[GSYM INT_POW_MUL; GSYM INT_NEG_MINUS1]));;




let INT_POW_EVEN_NEG1 = 
prove(
  `!x n. (--: (&: x)) **: (NUMERAL (BIT0 n)) =  ((&: x) **: (NUMERAL (BIT0 n)))`,

((REPEAT GEN_TAC))
THEN ((ONCE_REWRITE_TAC[INT_POW_NEG1]))
THEN ((ABBREV_TAC `a = &: 1`))
THEN ((ABBREV_TAC `b = (&: x)**: (NUMERAL (BIT0 n))`))
THEN ((REWRITE_TAC[NUMERAL;BIT0]))
THEN ((REWRITE_TAC[GSYM MULT_2;GSYM INT_POW_POW;INT_OF_REAL_THM REAL_POW_2;INT_NEG_MUL2]))
THEN ((EXPAND_TAC "a"))
THEN ((REWRITE_TAC[INT_MUL_RID;INT_MUL_LID;INT_OF_REAL_THM REAL_POW_ONE])));;


let INT_POW_ODD_NEG1 = 
prove(
  `!x n. (--: (&: x)) **: (NUMERAL (BIT1 n)) = --: ((&: x) **: (NUMERAL (BIT1 n)))`,

((REPEAT GEN_TAC))
THEN ((ONCE_REWRITE_TAC[INT_POW_NEG1]))
THEN (((ABBREV_TAC `a = &: 1`)))
THEN (((ABBREV_TAC `b = (&: x)**: (NUMERAL (BIT1 n))`)))
THEN ((REWRITE_TAC[NUMERAL;BIT1]))
THEN ((ONCE_REWRITE_TAC[ADD1]))
THEN ((EXPAND_TAC "a"))
THEN ((REWRITE_TAC[GSYM MULT_2]))
THEN ((REWRITE_TAC[INT_OF_REAL_THM POW_MINUS1;INT_OF_REAL_THM REAL_POW_ADD]))
THEN ((REWRITE_TAC[INT_OF_REAL_THM POW_1;INT_MUL_LID;INT_MUL_LNEG])));;

(* subtraction of integers *)


let INT_ADD_NEG = 
prove(
 `!x y. (x < y ==> ((&: x) +: (--: (&: y)) = --: (&: (y - x))))`,

((REPEAT GEN_TAC))
THEN ((DISCH_TAC))
THEN ((SUBGOAL_TAC `&: (y-x ) = (&: y) - (&: x)`))
THENL [((SUBGOAL_TAC `x <=| y`))
         (* branch *)
         THENL [(((ASM MESON_TAC)[LE_LT]));((SIMP_TAC[GSYM (INT_OF_REAL_THM REAL_OF_NUM_SUB)]))]
(* branch *)
; ((DISCH_TAC))
THEN ((ASM_REWRITE_TAC[]))
THEN (ACCEPT_TAC(INT_ARITH `&: x +: --: (&: y) = --: (&: y -: &: x)`))]);;


let INT_ADD_NEGv2 = 
prove(
 `!x y. (y <= x ==> ((&: x) +: (--: (&: y)) = (&: (x - y))))`,

 ((REPEAT GEN_TAC))
 THEN ((DISCH_TAC))
 THEN ((SUBGOAL_TAC `&: (x - y) = (&: x) - (&: y)`))
 THENL[
  ((UNDISCH_EL_TAC 0)) THEN ((SIMP_TAC[GSYM (INT_OF_REAL_THM REAL_OF_NUM_SUB)]));
  ((DISCH_TAC)) THEN ((ASM_REWRITE_TAC[INT_SUB]))
     ]
);;

(* assumes a term of the form &:a +: (--: (&: b))  *)
let INT_SUB_CONV t =
    let a,b = dest_binop `(+:)` t in
  let (_,a) = dest_comb a in
  let (_,b) = dest_comb b in
  let (_,b) = dest_comb b in
  let a = dest_numeral a in
  let b = dest_numeral b in
  let thm = if  (b <=/ a) then
    INT_ADD_NEGv2
  else INT_ADD_NEG in
  (ARITH_SIMP_CONV[thm]) t;; (*   (SIMP_CONV[thm;ARITH]) t;; *)


(* ------------------------------------------------------------------ *)
(*   Simplifies an arithmetic expression in floats                    *)
(*   A workhorse                                                      *)
(* ------------------------------------------------------------------ *)

let FLOAT_CONV =
              (ARITH_SIMP_CONV[FLOAT_MUL;FLOAT_SUB;FLOAT_ABS;FLOAT_POW;
              FLOAT_ADD_NN;FLOAT_ADD_NNv2;FLOAT_ADD_PP;FLOAT_ADD_PPv2;
              FLOAT_ADD_NP;FLOAT_ADD_PN;FLOAT_NEG;
              INT_NEG_NEG;INT_SUB;
              INT_ABS_NUM;INT_ABS_NEG_NUM;
              INT_MUL_LNEG;INT_MUL_RNEG;INT_NEG_MUL2;INT_OF_NUM_MUL;
              INT_OF_NUM_ADD;GSYM INT_NEG_ADD;INT_ADD_NEG_NUM;
              INT_OF_NUM_POW;INT_POW_ODD_NEG1;INT_POW_EVEN_NEG1;
              INT_ADD_NEG;INT_ADD_NEGv2 (* ; ARITH *)
              ]) ;;

add_test("FLOAT_CONV1",
  let f z =
    let (x,y) =  dest_eq z in
    let (u,v) =  dest_thm (FLOAT_CONV x) in
    (u=[]) && (z = v) in
  f `float (&:3) (&:0) = float (&:3) (&:0)` &&
  f `float (&:3) (&:3) = float (&:3) (&:3)` &&
  f `float (&:3) (&:0) +. (float (&:4) (&:0)) = (float (&:7) (&:0))` &&
  f `float (&:1 + (&:3)) (&:4) = float (&:4) (&:4)` &&
  f `float (&:3 - (&:7)) (&:0) = float (--:(&:4)) (&:0)` &&
  f `float (&:3) (&:4) *. (float (--: (&:2)) (&:3)) = float (--: (&:6))
                                                        (&:7)` &&
  f `--. (float (--: (&:3)) (&:0)) = float (&:3) (&:0)`
        );;

(* ------------------------------------------------------------------ *)
(*   Operations on interval. Preliminary stuff to deal with           *)
(*   chains of inequalities.                                          *)
(* ------------------------------------------------------------------ *)



let REAL_ADD_LE_SUBST_RHS = 
prove(
  `!a b c P. ((a <=. ((P b)) /\ (!x. (P x) =  x + (P (&. 0))) /\ (b <=. c)) ==> (a <=. (P c)))`,

(((REPEAT GEN_TAC)))
THEN (((REPEAT (TAUT_TAC `(b ==> a==> c) ==> (a /\ b ==> c)`))))
THEN (((REPEAT DISCH_TAC)))
THEN ((((H_RULER(ONCE_REWRITE_RULE))[HYP_INT 1] (HYP_INT 0))))
THEN ((((ASM ONCE_REWRITE_TAC)[])))
THEN ((((ASM MESON_TAC)[REAL_LE_RADD;REAL_LE_TRANS]))));;


let REAL_ADD_LE_SUBST_LHS = 
prove(
  `!a b c P. (((P(a) <=. b /\ (!x. (P x) =  x + (P (&. 0))) /\ (c <=. a)))
     ==> ((P c) <=. b))`,

(REP_GEN_TAC)
THEN (DISCH_ALL_TAC)
THEN (((H_RULER(ONCE_REWRITE_RULE)) [HYP_INT 1] (HYP_INT 0)))
THEN (((ASM ONCE_REWRITE_TAC)[]))
THEN (((ASM MESON_TAC)[REAL_LE_RADD;REAL_LE_TRANS])));;
(*
let rec SPECL =
    function [] -> I |
    (a::b)  -> fun thm ->(SPECL b (SPEC a thm));;
*)
(*
  input:
    rel: b <=. c
    thm: a <=. P(b).

  output: a <=. P(c).

  condition: REAL_ARITH must be able to prove !x. P(x) = x+. P(&.0).
  condition: the term `a` must appear exactly once the lhs of the thm.
  *)

let IWRITE_REAL_LE_RHS rel thm =
  let bvar = genvar `:real` in
  let (bt,_) = dest_binop `(<=.)` (concl rel) in
  let sub = SUBS_CONV[ASSUME (mk_eq(bt,bvar))] in
  let rule = (fun th -> EQ_MP (sub (concl th)) th) in
  let (subrel,subthm) = (rule rel,rule thm) in
  let (a,p) = dest_binop `(<=.)` (concl subthm) in
  let (_,c) = dest_binop `(<=.)` (concl subrel) in
  let pfn = mk_abs (bvar,p) in
  let imp_th = BETA_RULE (SPECL [a;bvar;c;pfn] REAL_ADD_LE_SUBST_RHS) in
  let ppart =   REAL_ARITH
      (fst(dest_conj(snd(dest_conj(fst(dest_imp(concl imp_th))))))) in
  let prethm = MATCH_MP imp_th (CONJ subthm (CONJ ppart subrel)) in
  let prethm2 = SPEC bt (GEN bvar (DISCH
       (find (fun x -> try(bvar = rhs x) with failure -> false) (hyp prethm)) prethm)) in
  MATCH_MP prethm2 (REFL bt);;

(*
  input:
    rel: c <=. a
    thm: P a <=. b

  output: P c <=. b

  condition: REAL_ARITH must be able to prove !x. P(x) = x+. P(&.0).
  condition: the term `a` must appear exactly once the lhs of the thm.
  *)

let IWRITE_REAL_LE_LHS rel thm =
  let avar = genvar `:real` in
  let (_,at) = dest_binop `(<=.)` (concl rel) in
  let sub = SUBS_CONV[ASSUME (mk_eq(at,avar))] in
  let rule = (fun th -> EQ_MP (sub (concl th)) th) in
  let (subrel,subthm) = (rule rel,rule thm) in
  let (p,b) = dest_binop `(<=.)` (concl subthm) in
  let (c,_) = dest_binop `(<=.)` (concl subrel) in
  let pfn = mk_abs (avar,p) in
  let imp_th = BETA_RULE (SPECL [avar;b;c;pfn] REAL_ADD_LE_SUBST_LHS) in
  let ppart =   REAL_ARITH
      (fst(dest_conj(snd(dest_conj(fst(dest_imp(concl imp_th))))))) in
  let prethm = MATCH_MP imp_th (CONJ subthm (CONJ ppart subrel)) in
  let prethm2 = SPEC at (GEN avar (DISCH
       (find (fun x -> try(avar = rhs x) with failure -> false) (hyp prethm)) prethm)) in
  MATCH_MP prethm2 (REFL at);;

(* ------------------------------------------------------------------ *)
(*   INTERVAL ADD, NEG, SUBTRACT                                      *)
(* ------------------------------------------------------------------ *)



let INTERVAL_ADD = 
prove(
   `!x f ex y g ey. interval x f ex /\ interval y g ey ==>
                         interval (x +. y) (f +. g) (ex +. ey)`,

EVERY[
 REPEAT GEN_TAC;
 TAUT_TAC `(A==>B==>C)==>(A/\ B ==> C)`;
 REWRITE_TAC[interval];
 REWRITE_TAC[prove(`(x+.y) -. (f+.g) = (x-.f) +. (y-.g)`,REAL_ARITH_TAC)];
 ABBREV_TAC `a = x-.f`;
 ABBREV_TAC `b = y-.g`;
 ASSUME_TAC (SPEC `b:real` (SPEC `a:real` ABS_TRIANGLE));
 UNDISCH_EL_TAC 0;
 ABBREV_TAC `a':real = abs a`;
 ABBREV_TAC `b':real = abs b`;
 REPEAT DISCH_TAC;
 (H_VAL2(IWRITE_REAL_LE_RHS)) (HYP_INT 0) (HYP_INT 2);
 (H_VAL2(IWRITE_REAL_LE_RHS)) (HYP_INT 2) (HYP_INT 0);
 ASM_REWRITE_TAC[]]);;


let INTERVAL_NEG = 
prove(
  `!x f ex. interval x f ex = interval (--. x) (--. f) ex`,

(REWRITE_TAC[interval;REAL_ABS_NEG;REAL_ARITH `!x y. -- x -. (-- y) = --. (x -. y)`]));;


let INTERVAL_NEG2 = 
prove(
  `!x f ex. interval (--. x) f ex = interval x (--. f) ex`,

 (REWRITE_TAC[interval;REAL_ABS_NEG;REAL_ARITH `!x y. -- x -. y = --. (x -. (--. y))`]));;



let INTERVAL_SUB = 
prove(
   `!x f ex y g ey. interval x f ex /\ interval y g ey ==> interval (x -. y) (f -. g) (ex +. ey)`,

((REWRITE_TAC[real_sub]))
THEN (DISCH_ALL_TAC)
THEN (((H_RULER (ONCE_REWRITE_RULE))[THM(INTERVAL_NEG)] (HYP_INT 1)))
THEN (((ASM MESON_TAC)[INTERVAL_ADD])));;

(* ------------------------------------------------------------------ *)
(*   INTERVAL MULTIPLICATION                                          *)
(* ------------------------------------------------------------------ *)



let REAL_PROP_LE_LABS = 
prove(
  `!x y z. (y <=. z) ==> ((abs x)* y <=. (abs x) *z)`,
(SIMP_TAC[REAL_LE_LMUL_IMP;ABS_POS]));;

(* renamed from REAL_LE_ABS_RMUL *)

let REAL_PROP_LE_RABS = 
prove(
  `!x y z. (y <=. z) ==> ( y * (abs x) <=. z *(abs x))`,
(SIMP_TAC[REAL_LE_RMUL_IMP;ABS_POS]));;


let REAL_LE_ABS_MUL = 
prove(
  `!x y z w. (( x <=. y) /\ (abs z <=. w)) ==> (x*.w <=. y*.w) `,

(DISCH_ALL_TAC)
THEN ((ASSUME_TAC (REAL_ARITH `abs z <=. w ==> (&.0) <=. w`)))
THEN (((ASM MESON_TAC)[REAL_LE_RMUL_IMP])));;


let INTERVAL_MUL = 
prove(
  `!x f ex y g ey. (interval x f ex) /\ (interval y g ey) ==>
         (interval (x *. y) (f *. g) (abs(f)*.ey +. ex*. abs(g) +. ex*.ey))`,

(REP_GEN_TAC)
THEN ((REWRITE_TAC[interval]))
THEN ((REWRITE_TAC[REAL_ARITH `(x*. y -. f*. g) = (f *.(y -. g) +. (x -. f)*.g +. (x-.f)*.(y-. g))`]))
THEN (DISCH_ALL_TAC)
THEN ((ASSUME_TAC (SPECL [`f*.(y-g)`;`(x-f)*g +. (x-f)*.(y-g)`] ABS_TRIANGLE)))
THEN ((ASSUME_TAC (SPECL [`(x-f)*.g`;`(x-f)*.(y-g)`] ABS_TRIANGLE)))
THEN (((H_VAL2(IWRITE_REAL_LE_RHS)) (HYP_INT 0) (HYP_INT 1)))
THEN ((H_REWRITE_RULE [THM ABS_MUL] (HYP_INT 0)))
THEN ((H_MATCH_MP (THM (SPECL [`g:real`; `abs (x -. f)`; `ex:real`] REAL_PROP_LE_RABS)) (HYP_INT 4)))
THEN (((H_VAL2(IWRITE_REAL_LE_RHS)) (HYP_INT 0) (HYP_INT 1)))
THEN ((H_MATCH_MP (THM (SPECL [`f:real`; `abs (y -. g)`; `ey:real`] REAL_PROP_LE_LABS)) (HYP_INT 7)))
THEN (((H_VAL2 (IWRITE_REAL_LE_RHS)) (HYP_INT 0) (HYP_INT 1)))
THEN ((H_MATCH_MP (THM (SPECL [`x-.f`; `abs (y -. g)`; `ey:real`] REAL_PROP_LE_LABS)) (HYP_INT 9)))
THEN (((H_VAL2(IWRITE_REAL_LE_RHS)) (HYP_INT 0) (HYP_INT 1)))
THEN ((ASSUME_TAC (SPECL [`abs(x-.f)`;`ex:real`;`y-.g`;`ey:real`] REAL_LE_ABS_MUL)))
THEN ((H_CONJ (HYP_INT 11) (HYP_INT 12)))
THEN ((H_MATCH_MP (HYP_INT 1) (HYP_INT 0)))
THEN (((H_VAL2(IWRITE_REAL_LE_RHS)) (HYP_INT 0) (HYP_INT 3)))
THEN ((POP_ASSUM ACCEPT_TAC)));;

(* ------------------------------------------------------------------ *)
(*   INTERVAL BASIC OPERATIONS                                        *)
(* ------------------------------------------------------------------ *)



let INTERVAL_NUM = 
prove(
  `!n. (interval(&.n) (float(&:n) (&:0)) (float  (&: 0) (&:0)))`,


let INTERVAL_CENTER = 
prove(
  `!x f ex g. (interval x f ex) ==> (interval x g (abs(f-g)+.ex))`,

((REWRITE_TAC[interval]))
THEN (DISCH_ALL_TAC)
THEN ((ASSUME_TAC (REAL_ARITH `abs(x -. g) <=. abs(f-.g) +. abs(x -. f)`)))
THEN ((H_VAL2 IWRITE_REAL_LE_RHS (HYP_INT 1) (HYP_INT 0)))
THEN ((ASM_REWRITE_TAC[])));;


let INTERVAL_WIDTH = 
prove(
  `!x f ex ex'. (ex <=. ex') ==> (interval x f ex) ==> (interval x f ex')`,

((REWRITE_TAC[interval]))
THEN (DISCH_ALL_TAC)
THEN ((H_VAL2 IWRITE_REAL_LE_RHS (HYP_INT 1) (HYP_INT 0)))
THEN ((ASM_REWRITE_TAC[])));;


let INTERVAL_MAX = 
prove(
  `!x f ex. interval x f ex ==> (x <=. f+.ex)`,

(REWRITE_TAC[interval]) THEN REAL_ARITH_TAC);;


let INTERVAL_MIN = 
prove(
  `!x f ex. interval x f ex ==> (f-. ex <=. x)`,

(REWRITE_TAC[interval]) THEN REAL_ARITH_TAC);;


let INTERVAL_ABS_MIN = 
prove(
  `!x f ex. interval x f ex ==> (abs(f)-. ex <=. abs(x))`,

  (REWRITE_TAC[interval] THEN REAL_ARITH_TAC)
);;


let INTERVAL_ABS_MAX = 
prove(
  `!x f ex. interval x f ex ==> (abs(x) <=. abs(f)+. ex)`,

  (REWRITE_TAC[interval] THEN REAL_ARITH_TAC)
);;


let REAL_RINV_2 = 
prove(
  `&.2 *. (inv (&.2 )) = &. 1`,

EVERY[
  MATCH_MP_TAC REAL_MUL_RINV;
  REAL_ARITH_TAC]);;


let INTERVAL_MK = 
prove(
   `let half = float(&:1)(--:(&:1)) in
    !x xmin xmax. ((xmin <=. x) /\ (x <=. xmax)) ==>
      interval x
         ((xmin+.xmax)*.half)
         ((xmax-.xmin)*.half)`,

EVERY[
  REWRITE_TAC[LET_DEF;LET_END_DEF];
  DISCH_ALL_TAC;
  REWRITE_TAC[interval;float;TWOPOW_NEG;INT_NUM_REAL;REAL_POW_1;REAL_MUL_LID];
  REWRITE_TAC[GSYM INTERVAL_ABS];
  CONJ_TAC
  ]
THENL[
  EVERY[
    REWRITE_TAC[GSYM REAL_SUB_RDISTRIB];
    REWRITE_TAC[REAL_ARITH `(b+.a)-.(a-.b)=b*.(&.2)`;GSYM REAL_MUL_ASSOC];
    ASM_REWRITE_TAC[REAL_RINV_2;REAL_MUL_RID]
  ];
  EVERY[
    REWRITE_TAC[GSYM REAL_ADD_RDISTRIB];
    REWRITE_TAC[REAL_ARITH `(b+.a)+. a -. b=a*.(&.2)`;GSYM REAL_MUL_ASSOC];
    ASM_REWRITE_TAC[REAL_RINV_2;REAL_MUL_RID]
  ]
]);;


let INTERVAL_EPS_POS = 
prove(`!x f ex.
  (interval x f ex) ==> (&.0 <=. ex)`,

EVERY[
  REWRITE_TAC[interval];
  REPEAT (GEN_TAC);
  DISCH_THEN(fun x -> (MP_TAC (CONJ (SPEC `x-.f` REAL_ABS_POS) x)));
  MATCH_ACCEPT_TAC REAL_LE_TRANS]);;


let INTERVAL_EPS_0 = 
prove(
  `!x f n. (interval x f (float (&:0) n)) ==> (x = f)`,

EVERY[
  REWRITE_TAC[interval;float;int_of_num_th;REAL_MUL_LZERO];
  REAL_ARITH_TAC]);;




let REAL_EQ_RCANCEL_IMP' = 
prove(`!x y z.(x * z = y * z) ==> (~(z = &0) ==> (x=y))`,

  MESON_TAC[REAL_EQ_RCANCEL_IMP]);;

(* renamed from REAL_ABS_POS *)

let REAL_MK_POS_ABS_' = 
prove (`!x. (~(x=(&.0))) ==> (&.0 < abs(x))`,

(* ------------------------------------------------------------------ *)
(*   INTERVAL DIVIDE                                                  *)
(* ------------------------------------------------------------------ *)


let INTERVAL_DIV = 
prove(`!x f ex y g ey h ez.
  (((interval x f ex) /\ (interval y g ey) /\ (ey <. (abs g)) /\
  ((ex +. (abs (f -. (h*.g))) +. (abs h)*. ey) <=. (ez*.((abs g) -. ey))))
  ==> (interval (x / y) h ez))`,



let lemma1 = prove( `&.0 < u /\ ||. z <=. e*. u ==> (&.0) <=. e`,
  EVERY[
    DISCH_ALL_TAC;
    ASSUME_TAC (SPEC `z:real` REAL_MK_NN_ABS);
    H_MATCH_MP (THM REAL_LE_TRANS) (H_RULE2 CONJ (HYP_INT 0) (HYP_INT 2));
    H_MATCH_MP (THM REAL_PROP_NN_RCANCEL) (H_RULE2 CONJ (HYP_INT 2) (HYP_INT 0));
    ASM_REWRITE_TAC[]
  ]) in
EVERY[
  DISCH_ALL_TAC;
  SUBGOAL_TAC `~(y= (&.0))`
  THENL[
    EVERY[
      UNDISCH_LIST[1;2];
      REWRITE_TAC[interval];
      REAL_ARITH_TAC
    ];
    EVERY[
      REWRITE_TAC[interval];
      DISCH_TAC THEN (H I (HYP_INT 0)) THEN (UNDISCH_EL_TAC 0);
      DISCH_THEN (fun th -> (MP_TAC(MATCH_MP REAL_MK_POS_ABS_' th)));
      MATCH_MP_TAC REAL_MUL_RTIMES_LE;
      REWRITE_TAC[GSYM ABS_MUL;REAL_SUB_RDISTRIB;real_div;GSYM REAL_MUL_ASSOC];
      ASM_SIMP_TAC[REAL_MUL_LINV;REAL_MUL_RID];
      H (REWRITE_RULE[interval]) (HYP_INT 1);
      H (REWRITE_RULE[interval]) (HYP_INT 3);
      H (MATCH_MP INTERVAL_ABS_MIN) (HYP_INT 4);
      POPL_TAC[3;4;5];
      H_VAL2 (IWRITE_REAL_LE_LHS) (HYP_INT 2) (HYP_INT 4);
      H (REWRITE_RULE[ REAL_ADD_ASSOC]) (HYP_INT 0);
      H_VAL2 (IWRITE_REAL_LE_LHS) (THM (SPEC `f-. h*g` (SPEC `x-.f` ABS_TRIANGLE))) (HYP_INT 0);
      H (ONCE_REWRITE_RULE[REAL_ABS_SUB]) (HYP_INT 4);
      H (MATCH_MP (SPEC `h:real` REAL_PROP_LE_LABS)) (HYP_INT 0);
      H (REWRITE_RULE[GSYM ABS_MUL]) (HYP_INT 0);
      H_VAL2 (IWRITE_REAL_LE_LHS) (HYP_INT 0) (HYP_INT 3);
      H_VAL2 (IWRITE_REAL_LE_LHS) (THM (SPEC `h*.(g-.y)` (SPEC`(x-.f)+(f-. h*g)`  ABS_TRIANGLE))) (HYP_INT 0);
      POPL_TAC[1;2;3;4;5;6;7;9;10;12];
      H (ONCE_REWRITE_RULE[REAL_ARITH `((x-.f) +. (f -. h*. g)) +. h*.(g-. y) = x -. h*. y `]) (HYP_INT 0);
      ABBREV_TAC `z = x -. h*.y`;
      H (ONCE_REWRITE_RULE[GSYM REAL_SUB_LT]) (HYP_INT 4);
      ABBREV_TAC `u = abs(g) -. ey`;
      POPL_TAC[0;2;4;6];
      H (MATCH_MP lemma1 ) (H_RULE2 CONJ (HYP_INT 0) (HYP_INT 1));
      H (MATCH_MP REAL_PROP_LE_LMUL) (H_RULE2 CONJ (HYP_INT 0) (HYP_INT 3));
      H (MATCH_MP REAL_LE_TRANS) (H_RULE2 CONJ (HYP_INT 3) (HYP_INT 0));
      ASM_REWRITE_TAC[]
  ];
  ]]);;

(* ------------------------------------------------------------------ *)
(*   INTERVAL ABS VALUE                                               *)
(* ------------------------------------------------------------------ *)


let INTERVAL_ABSV = 
prove(`!x f ex. interval x f ex ==> (interval (abs x) (abs f) ex)`,

EVERY[
  REWRITE_TAC[interval];
  DISCH_ALL_TAC;
  ASSUME_TAC (SPECL [`x:real`;`f:real`] REAL_ABS_SUB_ABS);
  ASM_MESON_TAC[REAL_LE_TRANS]
]);;
  (* 7 minutes *)

(* ------------------------------------------------------------------ *)
(*   INTERVAL SQRT                                                    *)
(*   This requires some preliminaries. Extend sqrt by 0 on negatives  *)
(* ------------------------------------------------------------------ *)


let ssqrt = new_definition `ssqrt x = if (x <. (&.0)) then (&.0) else sqrt x`;;
 (*2m*)

let LET_TAC = REWRITE_TAC[LET_DEF;LET_END_DEF];;



let REAL_SSQRT_NEG = 
prove(`!x. (x <. (&.0)) ==> (ssqrt x = (&.0))`,

  EVERY[
    DISCH_ALL_TAC;
    REWRITE_TAC[ssqrt];
    COND_CASES_TAC
    THENL[
      ACCEPT_TAC (REFL `&.0`);
      ASM_MESON_TAC[]
    ]
  ]);;
(* 5 min*)


let REAL_SSQRT_NN = 
prove(`!x. (&.0) <=. x ==> (ssqrt x = (sqrt x))`,

  EVERY[
  DISCH_ALL_TAC;
  REWRITE_TAC[ssqrt];
  COND_CASES_TAC
  THENL[
    ASM_MESON_TAC[real_lt];
    ACCEPT_TAC (REFL `sqrt x`)
  ]
  ]);;
  (* 12 min, mostly spent loading *index-shell* *)


(*17 minutes*)

let REAL_MK_NN_SSQRT = 
prove(`!x. (&.0) <=. (ssqrt x)`,

  EVERY[
    GEN_TAC;
    DISJ_CASES_TAC (SPECL[`x:real`;`&.0`] REAL_LTE_TOTAL)
    THENL[
      POP_ASSUM (fun th -> MP_TAC(MATCH_MP (REAL_SSQRT_NEG) th)) THEN
      MESON_TAC[REAL_LE_REFL];
      POP_ASSUM (fun th -> ASSUME_TAC(CONJ th (MATCH_MP (REAL_SSQRT_NN) th)))  THEN
      ASM_MESON_TAC[REAL_PROP_NN_SQRT]
    ]
  ]);;


let REAL_SV_SSQRT_0  = 
prove(`!x. ssqrt (&.0) = (&.0)`,

  EVERY[
    GEN_TAC;
    MP_TAC (SPEC `&.0` REAL_LE_REFL);
    DISCH_THEN(fun th -> REWRITE_TAC[MATCH_MP REAL_SSQRT_NN th]);
    ACCEPT_TAC REAL_SV_SQRT_0
  ]);;
 (* 6 minutes *)



let REAL_SSQRT_EQ_0 = 
prove(`!(x:real). (ssqrt(x) = (&.0)) ==> (x <=. (&.0))`,

  EVERY[
    GEN_TAC;
    DISJ_CASES_TAC (SPECL[`x:real`;`&.0`] REAL_LTE_TOTAL)
    THENL[
      ASM_MESON_TAC[REAL_LT_IMP_LE];
      ASM_SIMP_TAC[REAL_SSQRT_NN] THEN
      ASM_MESON_TAC[SQRT_EQ_0;REAL_EQ_IMP_LE]
    ]
  ]);;
  (* 15 minutes *)



let REAL_SSQRT_MONO = 
prove(`!x. (x<=. y) ==> (ssqrt x <=. (ssqrt y))`,

  EVERY[
    GEN_TAC;
    DISJ_CASES_TAC (SPECL[`x:real`;`&.0`] REAL_LTE_TOTAL)
      THENL[
        ASM_MESON_TAC[REAL_SSQRT_NEG;REAL_MK_NN_SSQRT];
        ASM_MESON_TAC[REAL_LE_TRANS;REAL_SSQRT_NN;REAL_PROP_LE_SQRT];
      ]
  ]);;
  (* 5 minutes *)


let REAL_SSQRT_CHAR = 
prove(`!x t. (&.0 <=. t /\ (t*t = x)) ==> (t = (ssqrt x))`,

  EVERY[
    DISCH_ALL_TAC;
    H_ASSUME_TAC (H_RULE_LIST REWRITE_RULE[HYP_INT 1] (THM (SPEC `t:real` REAL_MK_NN_SQUARE)));
    ASM_MESON_TAC[REAL_SSQRT_NN;SQRT_MUL;POW_2_SQRT_ABS;REAL_POW_2;REAL_ABS_REFL];
  ]);;
  (* 13 minutes *)


let REAL_SSQRT_SQUARE = 
prove(`!x. (&.0 <=. x) ==> ((ssqrt x)*.(ssqrt x) = x)`,
(* 7min *)


let REAL_SSQRT_SQUARE' = 
prove(`!x. (&.0<=. x) ==> (ssqrt (x*.x) = x)`,

  DISCH_ALL_TAC THEN
  REWRITE_TAC[(MATCH_MP REAL_SSQRT_NN (SPEC `x:real` REAL_MK_NN_SQUARE))] THEN
  ASM_SIMP_TAC[SQRT_MUL;GSYM POW_2;SQRT_POW_2]);;
 (*20min*)


(* an alternate proof appears in RCS *)

let INTERVAL_SSQRT = 
prove(`!x f ex u ey ez v. (interval x f ex) /\ (interval (u*.u) f ey) /\
  (ex +.ey <=. ez*.(v+.u)) /\ (v*.v <=. f-.ex) /\ (&.0 <. u) /\ (&.0 <=. v) ==>
  (interval (ssqrt x) u ez)`,

EVERY[
  DISCH_ALL_TAC;
  H (MATCH_MP REAL_LE_TRANS) (H_RULE2 CONJ (THM (SPEC `v:real` REAL_MK_NN_SQUARE)) (HYP_INT 3));
  H (MATCH_MP (INTERVAL_MIN)) (HYP_INT 1);
  H (MATCH_MP REAL_LE_TRANS)  (H_RULE2 CONJ (HYP_INT 1) (HYP_INT 0));
  H (MATCH_MP INTERVAL_EPS_POS) (HYP_INT 3);
  H (MATCH_MP INTERVAL_EPS_POS) (HYP_INT 5);
  H (MATCH_MP REAL_PROP_NN_ADD2) (H_RULE2 CONJ (HYP_INT 1) (HYP_INT 0));
  H (MATCH_MP REAL_PROP_POS_LADD) (H_RULE2 CONJ (HYP_INT 11) (HYP_INT 10));
  H (MATCH_MP REAL_PROP_POS_LADD) (H_RULE2 CONJ (THM (SPEC `x:real` REAL_MK_NN_SSQRT)) (HYP_INT 11));
  H (MATCH_MP REAL_PROP_POS_INV) (HYP_INT 0);
  ASSUME_TAC (REAL_ARITH  `(ssqrt x -. u) = (ssqrt x -. u)*.(&.1)`);
  H (MATCH_MP REAL_MK_NZ_POS) (HYP_INT 2);
  H (MATCH_MP REAL_MUL_RINV) (HYP_INT 0);
  H_REWRITE_RULE[(H_RULE GSYM) (HYP_INT 0)] (HYP_INT 2);
  POPL_TAC[1;2;3];
  H (REWRITE_RULE[REAL_MUL_ASSOC]) (HYP_INT 0);
  H (REWRITE_RULE[ONCE_REWRITE_RULE[REAL_MUL_SYM] REAL_DIFFSQ]) (HYP_INT 0);
  POPL_TAC[1;2];
  H_SIMP_RULE[HYP_INT 7;THM REAL_SSQRT_SQUARE] (HYP_INT 0);
  ASSUME_TAC (REAL_ARITH `abs(x -. u*.u) <=. abs(x -. f) + abs(f-. u*.u)`);
  H (REWRITE_RULE[interval]) (HYP_INT 12);
  H (ONCE_REWRITE_RULE[interval]) (HYP_INT 14);
  H (ONCE_REWRITE_RULE[REAL_ABS_SUB]) (HYP_INT 0);
  POPL_TAC[1;5;15;16];
  H (MATCH_MP REAL_LE_ADD2) (H_RULE2 CONJ (HYP_INT 1) (HYP_INT 0));
  H (MATCH_MP REAL_LE_TRANS) (H_RULE2 CONJ (HYP_INT 3) (HYP_INT 0));
  POPL_TAC[1;2;3;4];
  H (AP_TERM `||.`) (HYP_INT 1);
  H (REWRITE_RULE[ABS_MUL]) (HYP_INT 0);
  H (MATCH_MP REAL_LT_IMP_LE)  (HYP_INT 4);
  H (REWRITE_RULE[GSYM REAL_ABS_REFL]) (HYP_INT 0);
  H_REWRITE_RULE [HYP_INT 0] (HYP_INT 2);
  H (MATCH_MP REAL_LE_RMUL) (H_RULE2 CONJ (HYP_INT 5) (HYP_INT 2));
  H_REWRITE_RULE [H_RULE GSYM (HYP_INT 1)] (HYP_INT 0);
  POPL_TAC[1;2;3;5;6;7;8];
  H (MATCH_MP REAL_LE_TRANS) (H_RULE2 CONJ (HYP_INT 12) (HYP_INT 9));
  H (MATCH_MP REAL_SSQRT_MONO) (HYP_INT 0);
  H (MATCH_MP REAL_SSQRT_SQUARE') (HYP_INT 16);
  H_REWRITE_RULE [HYP_INT 0] (HYP_INT 1);
  H (ONCE_REWRITE_RULE[GSYM (SPECL[`v:real`;`ssqrt x`;`u:real`] REAL_LE_RADD)]) (HYP_INT 0);
  H (MATCH_MP REAL_LE_INV2) (H_RULE2 CONJ (HYP_INT 9) (HYP_INT 0));
  POPL_TAC[1;2;3;4;5;7;8;9;12;13];
  H (MATCH_MP REAL_LE_LMUL) (H_RULE2 CONJ (HYP_INT 3) (HYP_INT 0));
  H (MATCH_MP REAL_LE_TRANS) (H_RULE2 CONJ (HYP_INT 2) (HYP_INT 0));
  H (MATCH_MP REAL_PROP_POS_INV) (HYP_INT 4);
  H (MATCH_MP REAL_LT_IMP_LE) (HYP_INT 0);
  H (MATCH_MP REAL_LE_RMUL) (H_RULE2 CONJ (HYP_INT 11) (HYP_INT 0));
  H (REWRITE_RULE[GSYM REAL_MUL_ASSOC]) (HYP_INT 0);
  H (MATCH_MP REAL_MK_NZ_POS) (HYP_INT 8);
  H (MATCH_MP REAL_MUL_RINV) (HYP_INT 0);
  H_REWRITE_RULE[HYP_INT 0; THM REAL_MUL_RID] (HYP_INT 2);
  H (MATCH_MP REAL_LE_TRANS) (H_RULE2 CONJ (HYP_INT 7) (HYP_INT 0));
  ASM_REWRITE_TAC[interval]
  ]);;



test();;


(* conversion for interval *)

(* ------------------------------------------------------------------ *)
(*   Take a term x of type real.  Convert to a thm of the form        *)
(*   interval x f eps                                                 *)
(*                                                                    *)
(* ------------------------------------------------------------------ *)

let DOUBLE_CONV_FILE=true;;

let add_test,test = new_test_suite();;

(* Num package docs at http://caml.inria.fr/ocaml/htmlman/libref/Num.html *)

(* ------------------------------------------------------------------ *)
(* num_exponent
   Take the absolute value of input.
   Write it as a*2^k, where 1 <= a < 2, return k.

   Except:
   num_exponent (Int 0) is -1.
*)
let (num_exponent:Num.num -> Num.num) =
  fun a ->
    let afloat = float_of_num (abs_num a) in
    Int ((snd (frexp afloat)) - 1);;

(*test*)let f (u,v) = ((num_exponent u) =(Int v)) in
    add_test("num_exponenwt",
                forall f
    [Int 1,0; Int 65,6; Int (-65),6;
     Int 0,-1; (Int 3)//(Int 4),-1]);;
(* ------------------------------------------------------------------ *)

let dest_unary op tm =
  try let xop,r = dest_comb tm in
      if xop = op then r else fail()
  with Failure _ -> failwith "dest_unary";;


(* ------------------------------------------------------------------ *)


(* finds a nearby (outward-rounded) Int with only prec_b significant bits *)
let (round_outward: int -> Num.num -> Num.num) =
  fun prec_b a ->
    let b = abs_num a in
    let sign = if (a =/ b) then I else minus_num in
    let throw_bits = Num.max_num (Int 0) ((num_exponent b)-/ (Int prec_b)) in
    let twoexp = power_num (Int 2) throw_bits  in
    (sign (ceiling_num (b // twoexp)))*/twoexp;;

let (round_inward: int-> Num.num -> Num.num) =
  fun prec_b a ->
    let b = abs_num a in
    let sign = if (a=/b) then I else minus_num in
    let throw_bits = Num.max_num (Int 0) ((num_exponent b)-/ (Int prec_b)) in
    let twoexp = power_num (Int 2) throw_bits  in
    (sign (floor_num (b // twoexp)))*/twoexp;;

let round_rat bprec n =
  let b = abs_num n in
  let sign = if (b =/ n) then I else minus_num in
  let powt  = ((Int 2) **/ (Int bprec)) in
  sign ((round_outward bprec (Num.ceiling_num (b */ powt)))//powt);;

let round_inward_rat bprec n =
  let b = abs_num n in
  let sign = if (b =/ n) then I else minus_num in
  let powt  = ((Int 2) **/ (Int bprec)) in
  sign ((round_inward bprec (Num.floor_num (b */ powt)))//powt);;

let (round_outward_float: int -> float -> Num.num) =
 fun  bprec f ->
  if (f=0.0) then (Int 0) else
  begin
    let b = abs_float f in
    let sign = if (f >= 0.0) then I else minus_num in
    let (x,n) = frexp b in
    let u = int_of_float( ceil (ldexp x bprec)) in
    sign ((Int u)*/ ((Int 2) **/ (Int (n - bprec))))
  end;;

let (round_inward_float: int -> float -> Num.num) =
 fun  bprec f ->
  if (f=0.0) then (Int 0) else
  begin
    (* avoid overflow on 30 bit integers *)
    let bprec = if (bprec > 25) then 25 else bprec in
    let b = abs_float f in
    let sign = if (f >= 0.0) then I else minus_num in
    let (x,n) = frexp b in
    let u = int_of_float( floor (ldexp x bprec)) in
    sign ((Int u)*/ ((Int 2) **/ (Int (n - bprec))))
  end;;

(* ------------------------------------------------------------------ *)

(* This doesn't belong here.  A general term substitution function *)
let SUBST_TERM sublist tm =
  rhs (concl ((SPECL (map fst sublist)) (GENL (map snd sublist)
                                          (REFL tm))));;

add_test("SUBST_TERM",
 SUBST_TERM [(`#1`,`a:real`);(`#2`,`b:real`)] (`a +. b +. c`) =
 `#1 + #2 + c`);;

(* ------------------------------------------------------------------ *)

(* take a term of the form `interval x f ex` and clean up the f and ex *)

let INTERVAL_CLEAN_CONV:conv =
  fun interv ->
    let (ixf,ex) = dest_comb interv in
    let (ix,f) = dest_comb ixf in
    let fthm = FLOAT_CONV f in
    let exthm = FLOAT_CONV ex in
    let ixfthm = AP_TERM ix fthm in
    MK_COMB (ixfthm, exthm);;

(*test*) add_test("INTERVAL_CLEAN_CONV",
  let testval = INTERVAL_CLEAN_CONV `interval ((&.1) +. (&.1))
       (float (&:3) (&:4) +. (float (&:2) (--: (&:3))))
       (float (&:1) (&:2) *. (float (&:3) (--: (&:2))))` in
  let hypval = hyp testval in
  let concval = concl testval in
        (length hypval = 0) &&
        concval =
     `interval (&1 + &1) (float (&:3) (&:4) + float (&:2) (--: (&:3)))
     (float (&:1) (&:2) * float (&:3) (--: (&:2))) =
     interval (&1 + &1) (float (&:386) (--: (&:3))) (float (&:3) (&:0))`
                  );;

(* ------------------------------------------------------------------ *)
(*   GENERAL lemmas                                                   *)
(* ------------------------------------------------------------------ *)


(* verifies statement of the form `float a b = float a' b'` *)


let FLOAT_EQ = 
prove(
  `!a b a' b'.  (float a b = (float a'  b')) <=>
        ((float a b) -. (float a' b') = (&.0))`,
MESON_TAC[REAL_SUB_0]);;


let FLOAT_LT = 
prove(
  `!a b a' b'. (float a b <. (float a' b')) <=>
        ((&.0) <. (float a' b') -. (float a b))`,
MESON_TAC[REAL_SUB_LT]);;


let FLOAT_LE = 
prove(
  `!a b a' b'. (float a b <=. (float a' b')) <=>
        ((&.0) <=. (float a' b') -. (float a b))`,
MESON_TAC[REAL_SUB_LE]);;


let TWOPOW_MK_POS = 
prove(
  `!a. (&.0 <. ( twopow a))`,

EVERY[
  GEN_TAC;
  CHOOSE_TAC (SPEC `a:int` INT_REP2);
  POP_ASSUM DISJ_CASES_TAC;
  ASM_REWRITE_TAC[TWOPOW_POS;TWOPOW_NEG];
  TRY (MATCH_MP_TAC REAL_INV_POS);
  MATCH_MP_TAC REAL_POW_LT ;
  REAL_ARITH_TAC;
]);;


let TWOPOW_NZ = 
prove(
  `!a. ~(twopow a = (&.0))`,

  GEN_TAC THEN
  ACCEPT_TAC (MATCH_MP REAL_MK_NZ_POS (SPEC `a:int` TWOPOW_MK_POS)));;


let FLOAT_ZERO = 
prove(
  `!a b. (float a b = (&.0)) <=> (a = (&:0))`,

EVERY[
  REWRITE_TAC[float;REAL_ENTIRE;INT_OF_NUM_DEST];
  MESON_TAC[TWOPOW_NZ];
]);;


let INT_ZERO = 
prove(
  `!n. ((&:n = (&:0)) = (n=0))`,
REWRITE_TAC[INT_OF_NUM_EQ]);;


let INT_ZERO_NEG=
prove(
  `!n. ((--: (&:n) = (&:0))) <=> (n=0)`,

    REWRITE_TAC[INT_NEG_EQ_0;INT_ZERO]);;


let FLOAT_NN = 
prove(
  `!a b. ((&.0) <=. (float a b)) <=> (&:0 <=: a)`,

EVERY[
  REWRITE_TAC[float;INT_OF_NUM_DEST];
  REP_GEN_TAC;
  EQ_TAC THENL[EVERY[
  DISCH_ALL_TAC;
  INPUT_COMBO[THM REAL_PROP_NN_RCANCEL;THM (SPEC `b:int` TWOPOW_MK_POS) &&& (HYP"0")];
  ASM_MESON_TAC[int_le;int_of_num_th]];
  EVERY[
  DISCH_ALL_TAC;
  INPUT_COMBO[THM REAL_PROP_NN_POS;THM(SPEC`b:int`TWOPOW_MK_POS)];
  INPUT_COMBO[THM int_of_num_th   ; THM int_le ;(HYP"0")];
  INPUT_COMBO[THM REAL_PROP_NN_MUL2; (HYP"2")&&&(HYP"1")];
  ASM_REWRITE_TAC[]]]
]);;

let INT_NN = INT_POS;;


let INT_NN_NEG = 
prove(`!n. ((&:0) <=: (--:(&:n))) <=> (n = 0)`,

  REWRITE_TAC[INT_NEG_GE0;INT_OF_NUM_LE] THEN ARITH_TAC
                      );;


let FLOAT_POS = 
prove(`!a b. ((&.0) <. (float a b)) <=> (&:0 <: a)`,


let INT_POS' = 
prove(`!n. (&:0) <: (&:n) <=> (~(n=0) )`,

  REWRITE_TAC[INT_OF_NUM_LT] THEN ARITH_TAC);;


let INT_POS_NEG =
prove(`!n. ((&:0) <: (--:(&:n))) <=> F`,

  REWRITE_TAC[INT_OF_NUM_LT] THEN ARITH_TAC);;


let RAT_LEMMA1_SUB = 
prove(`~(y1 = &0) /\ ~(y2 = &0) ==>
      ((x1 / y1) - (x2 / y2) = (x1 * y2 - x2 * y1) * inv(y1) * inv(y2))`,

  EVERY[REWRITE_TAC[real_div];
  REWRITE_TAC[real_sub;GSYM REAL_MUL_LNEG];
  REWRITE_TAC[GSYM real_div];
  SIMP_TAC[RAT_LEMMA1];
  DISCH_TAC;
  MESON_TAC[real_div]]);;


let INTERVAL_0 = 
prove(`! a f ex. (interval a f ex <=> (&.0 <= (ex - (abs (a -. f)))))`,

   MESON_TAC[interval;REAL_SUB_LE]);;




let ABS_NUM = 
prove (`!m n. abs (&. n -. (&. m)) = &.((m-|n) + (n-|m))`,

  REPEAT GEN_TAC THEN
  DISJ_CASES_TAC (SPECL [`m:num`;`n:num`] LTE_CASES) THENL[
  (* first case *)
  EVERY[ LABEL_ALL_TAC;
  H_REWRITE_RULE [THM (GSYM REAL_OF_NUM_LT)] (HYP "0");
  LABEL_ALL_TAC;
  H_ONCE_REWRITE_RULE[THM (GSYM REAL_SUB_LT)] (HYP "1");
  LABEL_ALL_TAC;
  H_MATCH_MP (THM REAL_LT_IMP_LE) (HYP "2");
  LABEL_ALL_TAC;
  H_REWRITE_RULE [THM (GSYM ABS_REFL)] (HYP "3");
  ASM_REWRITE_TAC[];
  H_MATCH_MP (THM LT_IMP_LE) (HYP "0");
  ASM_SIMP_TAC[REAL_OF_NUM_SUB];
  REWRITE_TAC[REAL_OF_NUM_EQ];
  ONCE_REWRITE_TAC[ARITH_RULE `!x:num y:num. (x = y) = (y  = x)`];
  REWRITE_TAC[EQ_ADD_RCANCEL_0];
  ASM_REWRITE_TAC[SUB_EQ_0]];
  (* second case *)
  EVERY[LABEL_ALL_TAC;
  H_REWRITE_RULE [THM (GSYM REAL_OF_NUM_LE)] (HYP "0");
  LABEL_ALL_TAC;
  H_ONCE_REWRITE_RULE[THM (GSYM REAL_SUB_LE)] (HYP "1");
  LABEL_ALL_TAC;
  H_REWRITE_RULE [THM (GSYM ABS_REFL)] (HYP "2");
  ONCE_REWRITE_TAC[GSYM REAL_ABS_NEG];
  REWRITE_TAC[REAL_ARITH `!x y. --.(x -. y) = (y-x)`];
  ASM_REWRITE_TAC[];
  ASM_SIMP_TAC[REAL_OF_NUM_SUB];
  REWRITE_TAC[REAL_OF_NUM_EQ];
  ONCE_REWRITE_TAC[ARITH_RULE `!x:num y:num. (x = y) <=> (y  = x)`];
  REWRITE_TAC[EQ_ADD_LCANCEL_0];
  ASM_REWRITE_TAC[SUB_EQ_0]]]);;


let INTERVAL_TO_LESS = 
prove(
  `!a f ex b g ey. ((interval a f ex) /\ (interval b g ey) /\
      (&.0 <. (g -. (ey +. ex +. f)))) ==> (a <. b)`,

   let lemma1 = REAL_ARITH `!ex ey f g. (&.0 <.
     (g -. (ey +. ex +. f))) ==> ((f +. ex)<. (g -. ey)) ` in
   EVERY[
   REPEAT GEN_TAC;
   DISCH_ALL_TAC;
   H_MATCH_MP (THM lemma1) (HYP "2");
   H_MATCH_MP (THM INTERVAL_MAX) (HYP "0");
   H_MATCH_MP (THM INTERVAL_MIN) (HYP "1");
   LABEL_ALL_TAC;
   H_MATCH_MP (THM REAL_LET_TRANS) (H_RULE2 CONJ (HYP "4") (HYP "5"));
   LABEL_ALL_TAC;
   H_MATCH_MP (THM REAL_LTE_TRANS) (H_RULE2 CONJ (HYP "6") (HYP "3"));
   ASM_REWRITE_TAC[]
   ]);;


let ABS_TO_INTERVAL = 
prove(
  `!c u k. (abs (c - u) <=. k) ==> (!f g ex ey.((interval u f ex) /\ (interval k g ey) ==> (interval c f (g+.ey+.ex))))`,

   EVERY[
   REWRITE_TAC[interval];
   DISCH_ALL_TAC;
   REPEAT GEN_TAC;
   DISCH_ALL_TAC;
   ONCE_REWRITE_TAC [REAL_ARITH `c -. f = (c-. u) + (u-. f)`];
   ONCE_REWRITE_TAC [REAL_ADD_ASSOC];
   ASSUME_TAC (SPECL [`c-.u`;`u-.f`] ABS_TRIANGLE);
   IMP_RES_THEN ASSUME_TAC (REAL_ARITH `||.(k-.g) <=. ey ==> (k <=. (g +. ey))`);
   MATCH_MP_TAC (REAL_ARITH `(?a b.((x <=. (a+.b)) /\ (a <=. u) /\ (b <=. v)))  ==> (x <=. (u +. v))`);
   EXISTS_TAC `abs (c-.u)`;
   EXISTS_TAC `abs(u-.f)`;
   ASM_REWRITE_TAC[];
   ASM_MESON_TAC[REAL_LE_TRANS];
   ]);;


(* end of general lemmas *)
(* ------------------------------------------------------------------ *)


(* ------------------------------------------------------------------ *)
(* Cache of computed constants (abs (c - u) <= k)  *)
(* ------------------------------------------------------------------ *)

let calculated_constants = ref ([]:(term*thm) list);;

let add_real_constant ineq =
  try(
  let (abst,k) = dest_binop `(<=.)` (concl ineq) in
  let (absh,cmu) = dest_comb abst in
  let (c,u) = dest_binop `(-.)` cmu in
  calculated_constants := (c,ineq)::(!calculated_constants))
  with _ ->
  (try(
  let (c,f,ex) = dest_interval (concl ineq) in
  calculated_constants :=  (c,ineq)::(!calculated_constants))
  with _ -> failwith "calculated_constants format : abs(c - u) <= k");;

let get_real_constant tm =
  assoc tm !calculated_constants;;

let remove_real_constant tm =
  calculated_constants :=
    filter (fun t -> not ((fst t) = tm)) !calculated_constants;;



(* ------------------------------------------------------------------ *)

(* term of the form '&.n'. Assume error checking done already. *)
let INTERVAL_OF_NUM:conv =
  fun tm ->
    let tm1 = snd (dest_comb tm) in
    let th1 = (ARITH_REWRITE_CONV[] tm1) in
    ONCE_REWRITE_RULE[AP_TERM `&.` (GSYM th1)]
      (SPEC (rhs (concl th1)) INTERVAL_NUM);;

add_test("INTERVAL_OF_NUM",
   dest_thm (INTERVAL_OF_NUM `&.3`) = ([],
   `interval (&3) (float (&:3) (&:0)) (float (&:0) (&:0))`));;

(* term of the form `--. (&.n)`.  Assume format checking already done. *)
let INTERVAL_OF_NEG:conv =
  fun tm ->
    let (sign,u) = dest_comb tm in
    let _ = assert(sign = `--.`) in
    let (amp,tm1) = (dest_comb u) in
    let _ = assert(amp = `&.`) in
    let th1 = (ARITH_REWRITE_CONV[] tm1) in
    ONCE_REWRITE_RULE[FLOAT_NEG] (
    ONCE_REWRITE_RULE[INTERVAL_NEG] (
    ONCE_REWRITE_RULE[AP_TERM `&.` (GSYM th1)] (
       (SPEC (rhs (concl th1)) INTERVAL_NUM))));;

add_test("INTERVAL_OF_NEG",
   dest_thm (INTERVAL_OF_NEG `--.(&. (3+4))`) =
   ([],`interval( --.(&.(3 + 4)) )
      (float (--: (&:7)) (&:0)) (float (&:0) (&:0))`));;

(* ------------------------------------------------------------------ *)

let INTERVAL_TO_LESS_CONV = fun thm1 thm2 ->
   let (a,f,ex) = dest_interval (concl thm1) in
   let (b,g,ey) = dest_interval (concl thm2) in
   let rthm = ASSUME `!f g ex ey. (&.0 <. (g -. (ey +. ex +. f)))` in
   let rspec = concl (SPECL [f;g;ex;ey] rthm) in
   let rspec_simp = FLOAT_CONV (snd (dest_binop `(<.)` rspec)) in
   
let rthm2 = prove (rspec,REWRITE_TAC[rspec_simp;FLOAT_POS;INT_POS';
                                        INT_POS_NEG] THEN ARITH_TAC) in
   let fthm = CONJ thm1 (CONJ thm2 rthm2) in
   MATCH_MP INTERVAL_TO_LESS fthm;;

add_test("INTERVAL_TO_LESS_CONV",
  let thm1 = ASSUME
   `interval (#0.1) (float (&:1) (--: (&:1))) (float (&:1) (--: (&:2)))` in
  let thm2 = ASSUME `interval (#7) (float (&:4) (&:1)) (float (&:1) (&:0))` in
  let thm3 = INTERVAL_TO_LESS_CONV thm1 thm2 in
    concl thm3 = `#0.1 <. (#7)`);;

add_test("INTERVAL_TO_LESS_CONV2",
   let (h,c) = dest_thm (INTERVAL_TO_LESS_CONV
     (INTERVAL_OF_NUM `&.3`) (INTERVAL_OF_NUM `&.8`)) in
     (h=[]) && (c = `&.3 <. (&.8)`));;

(* ------------------------------------------------------------------ *)

(* conversion for DEC <= posfloat and posfloat <= DEC *)


let lemma1 = 
prove(
  `!n m p. ((&.p/(&.m)) <= (&.n)) <=> ((&.p/(&.m)) <= (&.n)/(&.1))`,

  MESON_TAC[REAL_DIV_1]);;


let lemma2 = 
prove(
  `!n m p. ((&.p) <= ((&.n)/(&.m))) <=> ((&.p/(&.1)) <= (&.n)/(&.m))`,

  MESON_TAC[REAL_DIV_1]);;


let lemma3 = 
prove(`!a b c d. (
   ((0<b) /\ (0<d) /\ (a*d <=| c*b))
    ==> (&.a/(&.b) <=. ((&.c)/(&.d))))`,

   EVERY[REPEAT GEN_TAC;
   DISCH_ALL_TAC;
   ASM_SIMP_TAC[RAT_LEMMA4;REAL_LT;REAL_OF_NUM_MUL;REAL_LE]]);;

let DEC_FLOAT = EQT_ELIM o
   ARITH_SIMP_CONV[DECIMAL;float;TWOPOW_POS;TWOPOW_NEG;GSYM real_div;
       REAL_OF_NUM_POW;INT_NUM_REAL;REAL_OF_NUM_MUL;
       lemma1;lemma2;lemma3];;

add_test("DEC_FLOAT",
   let f c x =
      dest_thm (c x) = ([],x) in
   ((f DEC_FLOAT `#10.0 <= (float (&:3) (&:2))`) &&
    (f DEC_FLOAT `#10 <= (float (&:3) (&:2))`) &&
    (f DEC_FLOAT `#0.1 <= (float (&:1) (--: (&:2)))`) &&
    (f DEC_FLOAT `float (&:3) (&:2) <= (#13.0)`) &&
    (f DEC_FLOAT `float (&:3) (&:2) <= (#13)`) &&
    (f DEC_FLOAT `float (&:1) (--: (&:2)) <= (#0.3)`)));;
(* ------------------------------------------------------------------ *)
(* conversion for float inequalities *)

let FLOAT_INEQ_CONV t =
  let thm1=  (ONCE_REWRITE_CONV[GSYM REAL_SUB_LT;GSYM REAL_SUB_LE] t) in
  let rhsx= rhs (concl thm1) in
  
let thm2= prove(rhsx,REWRITE_TAC[FLOAT_CONV (snd (dest_comb rhsx))] THEN
                    REWRITE_TAC[FLOAT_NN;FLOAT_POS;INT_NN;INT_NN_NEG;
                       INT_POS';
INT_POS_NEG] THEN ARITH_TAC) in
  REWRITE_RULE[GSYM thm1] thm2;;

let t1 = `(float (&:3) (&:0)) +. (float (&:4) (&:0)) <. (float (&:8) (&:1))`;;


add_test("FLOAT_INEQ_CONV",
  let f c x =
    dest_thm (c x) = ([],x) in
  let t1 =
   `(float (&:3) (&:0)) +. (float (&:4) (&:0)) <. (float (&:8) (&:1))` in
    ((f FLOAT_INEQ_CONV t1)));;




(* ------------------------------------------------------------------ *)

(* converts a DECIMAL TO A THEOREM *)


let INTERVAL_MINMAX = 
prove(`!x f ex.
   ((f -. ex) <= x) /\ (x <=. (f +. ex)) ==> (interval x f ex)`,

   EVERY[REPEAT GEN_TAC;
   REWRITE_TAC[interval;ABS_BOUNDS];
   REAL_ARITH_TAC]);;


let INTERVAL_OF_DECIMAL bprec dec =
  let a_num = dest_decimal dec in
  let f_num = round_rat bprec a_num in
  let ex_num = round_rat bprec (Num.abs_num (f_num -/ a_num)) in
  let _ = assert (ex_num <=/ f_num) in
  let f = mk_float f_num in
  let ex= mk_float ex_num in
  let fplus_ex = FLOAT_CONV (mk_binop `(+.)` f ex) in
  let fminus_ex= FLOAT_CONV (mk_binop `(-.)` f ex) in
  let fplus_term = rhs (concl fplus_ex) in
  let fminus_term = rhs (concl fminus_ex) in
  let th1 = DEC_FLOAT (mk_binop `(<=.)` fminus_term dec) in
  let th2 = DEC_FLOAT (mk_binop `(<=.)` dec fplus_term) in
  let intv = mk_interval dec f ex in
  EQT_ELIM (SIMP_CONV[INTERVAL_MINMAX;fplus_ex;fminus_ex;th1;th2] intv);;

add_test("INTERVAL_OF_DECIMAL",
  let (h,c) = dest_thm (INTERVAL_OF_DECIMAL 4 `#36.1`) in
  let (x,f,ex) = dest_interval c in
   (h=[]) && (x = `#36.1`));;

add_test("INTERVAL_OF_DECIMAL2",
 can (fun() -> INTERVAL_TO_LESS_CONV (INTERVAL_OF_DECIMAL 4 `#33.33`)
  (INTERVAL_OF_DECIMAL 4 `#36.1`)) ());;

(*--------------------------------------------------------------------*)
(*   functions to check format.                                       *)
(*   There are various implicit rules:                                *)
(*   NUMERAL is followed by bits and no other kind of num, etc.       *)
(*   FLOAT a b, both a and b are &:NUMERAL or --:&:NUMERAL, etc.      *)
(*--------------------------------------------------------------------*)


(* converts exceptions to false *)
let falsify_ex f x = try (f x) with _ -> false;;

let is_bits_format  =
    let rec format x =
    if (x = `_0`) then true
    else let (h,t) = dest_comb x  in
      (((h = `BIT1`) or (h = `BIT0`)) && (format t))
    in falsify_ex format;;

let is_numeral_format =
    let fn x =
    let (h,t) = dest_comb x in
      ((h = `NUMERAL`) && (is_bits_format t)) in
    falsify_ex fn;;

let is_decimal_format  =
    let fn x =
      let (t1,t2) = dest_binop `DECIMAL` x in
        ((is_numeral_format t1) && (is_numeral_format t2)) in
    falsify_ex fn;;

let is_pos_int_format =
    let fn x =
      let (h,t) = dest_comb x in
       (h = `&:`) && (is_numeral_format t) in
    falsify_ex fn;;

let is_neg_int_format =
    let fn x =
      let (h,t) = dest_comb x in
        (h = `--:`) && (is_pos_int_format t) in
      falsify_ex fn;;

let is_int_format x =
  (is_neg_int_format x) or (is_pos_int_format x);;

let is_float_format =
    let fn x =
      let (t1,t2) = dest_binop `float` x in
      (is_int_format t1) && (is_int_format t2) in
    falsify_ex fn;;

let is_interval_format =
  let fn x =
    let (a,b,c) = dest_interval x in
      (is_float_format b) && (is_float_format c) in
    falsify_ex fn;;

let is_neg_real =
  let fn x =
     let (h,t) = dest_comb x in
      (h= `--.`) in
    falsify_ex fn;;

let is_real_num_format =
  let fn x =
    let (h,t) = dest_comb x in
      (h=`&.`) && (is_numeral_format t) in
  falsify_ex fn;;

let is_comb_of t u =
  let fn t u =
    t = (fst (dest_comb u)) in
  try (fn t u) with failure -> false;;

(* ------------------------------------------------------------------ *)
(* Heron's formula for the square root of A
   Return a value x that is always at most the actual square root
   and such that abs (x  - A/x ) < epsilon *)

let rec heron_sqrt depth A x eps =
    let half = (Int 1)//(Int 2) in
    if (depth <= 0) then raise (Failure "sqrt recursion depth exceeded") else
    if (Num.abs_num (x -/ (A//x) ) </ eps) & (x*/ x >=/ A)  then (A//x) else
    let x' = half */ (x +/ (A//x)) in
    heron_sqrt (depth -1) A x' eps;;


let INTERVAL_OF_TWOPOW = 
prove(
   `!n. interval (twopow n) (float (&:1) n) (float (&:0) (&:0))`,

   REWRITE_TAC[interval;float;int_of_num_th] THEN
   REAL_ARITH_TAC
   );;

(* ------------------------------------------------------------------ *)

let rec INTERVAL_OF_TERM bprec tm =
  (* treat cached values first *)
  if (can get_real_constant tm) then
    begin
    try(
    let int_thm = get_real_constant tm in
    if (can dest_interval (concl int_thm)) then int_thm
    else (
    let absthm = get_real_constant tm in
    let (abst,k) = dest_binop `(<=.)` (concl absthm) in
    let (absh,cmu) = dest_comb abst in
    let (c,u) = dest_binop `(-.)` cmu in
    let intk = INTERVAL_OF_TERM bprec k in
    let intu = INTERVAL_OF_TERM bprec u in
    let thm1 = MATCH_MP ABS_TO_INTERVAL absthm in
    let thm2 = MATCH_MP thm1 (CONJ intu intk) in
    let (_,f,ex)= dest_interval (concl thm2) in
    let fthm = FLOAT_CONV f in
    let exthm = FLOAT_CONV ex in
    let thm3 = REWRITE_RULE[fthm;exthm] thm2 in
    (add_real_constant thm3; thm3)
    ))
    with _ -> failwith "INTERVAL_OF_TERM : CONSTANT"
    end
  else if (is_real_num_format tm) then (INTERVAL_OF_NUM tm)
  else if (is_decimal_format tm) then (INTERVAL_OF_DECIMAL bprec tm)
  (* treat negative terms *)
  else if (is_neg_real tm) then
    begin
    try(
    let (_,t) = dest_comb tm in
    let int1 = INTERVAL_OF_TERM bprec t in
    let (_,b,_) = dest_interval (concl int1) in
    let thm1  = FLOAT_CONV (mk_comb (`--.`, b)) in
    REWRITE_RULE[thm1] (ONCE_REWRITE_RULE[INTERVAL_NEG] int1))
    with _ -> failwith "INTERVAL_OF_TERM : NEG"
    end
  (* treat abs value *)
  else if (is_comb_of `||.` tm) then
    begin
      try(
      let (_,b) = dest_comb tm in
      let b_int = MATCH_MP INTERVAL_ABSV (INTERVAL_OF_TERM bprec b) in
      let (_,f,_) = dest_interval (concl b_int) in
      let thm1 = FLOAT_CONV f in
      REWRITE_RULE[thm1] b_int)
      with _ -> failwith "INTERVAL_OF_TERM : ABS"
    end
  (* treat twopow *)
  else if (is_comb_of `twopow` tm) then
    begin
      try(
      let (_,b) = dest_comb tm in
      SPEC b INTERVAL_OF_TWOPOW
      )
      with _ -> failwith "INTERVAL_OF_TERM : TWOPOW"
    end
  (* treat addition *)
  else if (can (dest_binop `(+.)`) tm) then
    begin
    try(
    let (a,b) = dest_binop `(+.)` tm in
    let a_int = INTERVAL_OF_TERM bprec a in
    let b_int = INTERVAL_OF_TERM bprec b in
    let c_int = MATCH_MP INTERVAL_ADD (CONJ a_int b_int) in
    let (_,f,ex) = dest_interval (concl c_int) in
    let thm1 = FLOAT_CONV f and thm2 = FLOAT_CONV ex in
    REWRITE_RULE[thm1;thm2] c_int)
    with _ -> failwith "INTERVAL_OF_TERM : ADD"
    end
  (* treat subtraction *)
  else if (can (dest_binop `(-.)`) tm) then
    begin
    try(
    let (a,b) = dest_binop `(-.)` tm in
    let a_int = INTERVAL_OF_TERM bprec a in
    let b_int = INTERVAL_OF_TERM bprec b in
    let c_int = MATCH_MP INTERVAL_SUB (CONJ a_int b_int) in
    let (_,f,ex) = dest_interval (concl c_int) in
    let thm1 = FLOAT_CONV f and thm2 = FLOAT_CONV ex in
    REWRITE_RULE[thm1;thm2] c_int)
    with _ -> failwith "INTERVAL_OF_TERM : SUB"
    end
  (* treat multiplication *)
  else if (can (dest_binop `( *. )`) tm) then
    begin
    try(
    let (a,b) = dest_binop `( *. )` tm in
    let a_int = INTERVAL_OF_TERM bprec a in
    let b_int = INTERVAL_OF_TERM bprec b in
    let c_int = MATCH_MP INTERVAL_MUL (CONJ a_int b_int) in
    let (_,f,ex) = dest_interval (concl c_int) in
    let thm1 = FLOAT_CONV f and thm2 = FLOAT_CONV ex in
    REWRITE_RULE[thm1;thm2] c_int)
    with _ -> failwith "INTERVAL_OF_TERM : MUL"
    end
  (* treat division : instantiate INTERVAL_DIV *)
  else if (can (dest_binop `( / )`) tm) then
    begin
    try(
    let (a,b) = dest_binop `( / )` tm in
    let a_int = INTERVAL_OF_TERM bprec a in
    let b_int = INTERVAL_OF_TERM bprec b in
    let (_,f,ex) = dest_interval (concl a_int) in
    let (_,g,ey) = dest_interval (concl b_int) in
    let f_num = dest_float f in
    let ex_num = dest_float ex in
    let g_num = dest_float g in
    let ey_num = dest_float ey in
    let h_num = round_rat bprec (f_num//g_num) in
    let h = mk_float h_num in
    let ez_rat = (ex_num +/ abs_num (f_num -/ (h_num*/ g_num))
        +/ (abs_num h_num */ ey_num))//((abs_num g_num) -/ (ey_num)) in
    let ez_num = round_rat bprec (ez_rat) in
    let _ = assert((ez_num >=/ (Int 0))) in
    let ez = mk_float ez_num in
    let hyp1 = a_int in
    let hyp2 = b_int in
    let hyp3 = FLOAT_INEQ_CONV (mk_binop `(<.)` ey (mk_comb (`||.`,g))) in
    let thm = SPECL [a;f;ex;b;g;ey;h;ez] INTERVAL_DIV in
    let conj2 x = snd (dest_conj x) in
    let hyp4t = (conj2 (conj2 (conj2 (fst(dest_imp (concl thm)))))) in
    let hyp4 = FLOAT_INEQ_CONV hyp4t in
    let hyp_all = end_itlist CONJ [hyp1;hyp2;hyp3;hyp4] in
    MATCH_MP thm hyp_all)
    with _ -> failwith "INTERVAL_OF_TERM :DIV"
    end
  (* treat sqrt : instantiate INTERVAL_SSQRT *)
  else if (can (dest_unary `ssqrt`) tm) then
    begin
    try(
    let x = dest_unary `ssqrt` tm in
    let x_int = INTERVAL_OF_TERM bprec x in
    let (_,f,ex)  = dest_interval (concl x_int) in
    let f_num = dest_float f in
    let ex_num = dest_float ex in
    let fd_num = f_num -/ ex_num in
    let fe_f = Num.float_of_num fd_num in
    let apprx_sqrt = Pervasives.sqrt fe_f in
    (* put in heron's formula *)
    let v_num1 = round_inward_float 25 (apprx_sqrt) in
    let v_num = round_inward_rat bprec
         (heron_sqrt 10 fd_num v_num1 ((Int 2) **/ (Int (-bprec-4)))) in
    let u_num1 = round_inward_float 25
        (Pervasives.sqrt (float_of_num f_num)) in
    let u_num = round_inward_rat bprec
        (heron_sqrt 10 f_num u_num1 ((Int 2) **/ (Int (-bprec-4)))) in
    let ey_num = round_rat bprec (abs_num (f_num -/ (u_num */ u_num))) in
    let ez_num = round_rat bprec ((ex_num +/ ey_num)//(u_num +/ v_num)) in
    let (v,u) = (mk_float v_num,mk_float u_num) in
    let (ey,ez) = (mk_float ey_num,mk_float ez_num) in
    let thm = SPECL [x;f;ex;u;ey;ez;v] INTERVAL_SSQRT in
    let conjhyp = fst (dest_imp (concl thm)) in
    let [hyp6;hyp5;hyp4;hyp3;hyp2;hyp1] =
      let rec break_conj c acc =
        if (not(is_conj c)) then (c::acc) else
        let (u,v) = dest_conj c in break_conj v (u::acc) in
       (break_conj conjhyp []) in
    
let thm2 = prove(hyp2,REWRITE_TAC[interval] THEN
                       (CONV_TAC FLOAT_INEQ_CONV)) in
    let thm3 = FLOAT_INEQ_CONV hyp3 in
    let thm4 = FLOAT_INEQ_CONV hyp4 in
    let float_tac = REWRITE_TAC[FLOAT_NN;FLOAT_POS;INT_NN;INT_NN_NEG;
                       INT_POS';
INT_POS_NEG] THEN ARITH_TAC in
    
let thm5 = prove( hyp5,float_tac) in
    let thm6 = prove( hyp6,float_tac) in
    let ant  = end_itlist CONJ[x_int;thm2;thm3;thm4;thm5;thm6] in
    MATCH_MP thm ant
    )
    with _ -> failwith "INTERVAL_OF_TERM : SSQRT"
    end
  else failwith "INTERVAL_OF_TERM : case not installed";
;


let real_ineq bprec tm =
  let (t1,t2) = dest_binop `(<.)` tm in
  let int1 = INTERVAL_OF_TERM bprec t1 in
  let int2 = INTERVAL_OF_TERM bprec t2 in
  INTERVAL_TO_LESS_CONV int1 int2;;

pop_priority();;