(*****************************)
(* Floating point arithmetic *)
(*****************************)

(* Dependencies *)
needs "../formal_lp/arith/nat.hl";;
needs "../formal_lp/arith/num_exp_theory.hl";;
needs "../formal_lp/arith/float_theory.hl";;
needs "../formal_lp/arith/interval_arith.hl";;


prioritize_real();;


module type Arith_float_sig =
  sig
    val mk_num_exp : term -> term -> term
    val dest_num_exp : term -> (term * term)
    val dest_float : term -> (string * term * term)
    val float_lt0 : term -> thm
    val float_gt0 : term -> thm
    val float_lt : term -> term -> thm
    val float_le0 : term -> thm
    val float_ge0 : term -> thm
    val float_le : term -> term -> thm
    val float_min : term -> term -> thm
    val float_max : term -> term -> thm
    val float_min_max : term -> term -> (thm * thm)
    val float_mul_eq : term -> term -> thm
    val float_mul_lo : int -> term -> term -> thm
    val float_mul_hi : int -> term -> term -> thm
    val float_div_lo : int -> term -> term -> thm
    val float_div_hi : int -> term -> term -> thm
    val float_add_lo : int -> term -> term -> thm
    val float_add_hi : int -> term -> term -> thm
    val float_sub_lo : int -> term -> term -> thm
    val float_sub_hi : int -> term -> term -> thm
    val float_sqrt_lo : int -> term -> thm
    val float_sqrt_hi : int -> term -> thm

    val reset_stat : unit -> unit
    val reset_cache : unit -> unit
    val print_stat : unit -> unit
(*
    val mul_lo_list : (int * term * term) list ref
    val mul_hi_list : (int * term * term) list ref
    val div_lo_list : (int * term * term) list ref
    val div_hi_list : (int * term * term) list ref
    val add_lo_list : (int * term * term) list ref
    val add_hi_list : (int * term * term) list ref
    val sub_lo_list : (int * term * term) list ref
    val sub_hi_list : (int * term * term) list ref
*)

    val dest_float_interval : term -> term * term * term
    val mk_float_interval_small_num : int -> thm
    val mk_float_interval_num : num -> thm

    val float_lo : int -> term -> thm
    val float_hi : int -> term -> thm

    val float_interval_round : int -> thm -> thm

    val float_interval_neg : thm -> thm
    val float_interval_mul : int -> thm -> thm -> thm
    val float_interval_div : int -> thm -> thm -> thm
    val float_interval_add : int -> thm -> thm -> thm
    val float_interval_sub : int -> thm -> thm -> thm
    val float_interval_sqrt : int -> thm -> thm

    val float_abs : term -> thm
	

    val FLOAT_TO_NUM_CONV : term -> thm
end;;



module Arith_float : Arith_float_sig = struct


open Big_int;;

open Arith_options;;
open Arith_misc;;
open Arith_nat;;
open Num_exp_theory;;
open Float_theory;;
open Interval_arith;;



(* interval *)


let APPROX_INTERVAL' = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP]) APPROX_INTERVAL;;


let zero_const = `_0` and
    t_const = `T` and
    f_const = `F`;;

let n_var_num = `n:num` and
    m_var_num = `m:num` and
    k_var_num = `k:num` and
    e_var_num = `e:num` and
    e1_var_num = `e1:num` and
    e2_var_num = `e2:num` and
    r_var_num = `r:num` and
    r1_var_num = `r1:num` and
    r2_var_num = `r2:num` and
    n1_var_num = `n1:num` and
    n2_var_num = `n2:num` and
    m1_var_num = `m1:num` and
    m2_var_num = `m2:num` and
    x_var_num = `x:num` and
    y_var_num = `y:num`;;

let le_op_num = `(<=):num->num->bool` and
    lt_op_num = `(<):num->num->bool` and
    mul_op_num = `( * ):num->num->num` and
    plus_op_num = `(+):num->num->num` and
    minus_op_num = `(-):num->num->num` and
    div_op_real = `(/):real->real->real` and
    div_op_num = `DIV` and
    mul_op_real = `( * ):real->real->real` and
    plus_op_real = `(+):real->real->real` and
    minus_op_real = `(-):real->real->real` and
    neg_op_real = `(--):real->real` and
    le_op_real = `(<=):real->real->bool`;;



let amp_op_real = `(&):num->real` and
    interval_const = `interval_arith` and
    num_exp_const = `num_exp`;;

let x_var_real = `x:real` and
    y_var_real = `y:real` and
    z_var_real = `z:real` and
    a_var_real = `a:real` and
    b_var_real = `b:real` and
    lo_var_real = `lo:real` and
    hi_var_real = `hi:real` and
    f1_var_real = `f1:real` and
    f2_var_real = `f2:real`;;



let b0_const = (fst o dest_comb o lhand o concl) (Arith_hash.def_array.(0));;
let b0_name = (fst o dest_const) b0_const;;
let base_const = mk_small_numeral base;;


let NUM_REMOVE = prove(mk_eq(mk_comb(Arith_hash.num_const, n_var_num), n_var_num), 
		       REWRITE_TAC[Arith_hash.num_def; NUMERAL]);;

(* B0 n = base * n *)
let b0_thm = prove(mk_eq(mk_comb(b0_const, n_var_num),
                         mk_binop mul_op_num base_const n_var_num),
                   REWRITE_TAC[Arith_hash.def_array.(0)] THEN
                     TRY ARITH_TAC THEN
                     GEN_REWRITE_TAC (LAND_CONV o DEPTH_CONV) [BIT0] THEN
                     ARITH_TAC);;

let dest_num_exp tm =
  let ltm, e_tm = dest_comb tm in
    rand ltm, e_tm;;



let num_exp_const = `num_exp`;;
let mk_num_exp n_tm e_tm = mk_binop num_exp_const n_tm e_tm;;


(* float s n e -> "s", n, e *)
let dest_float tm =
  let ltm, e_tm = dest_comb tm in
  let ltm, n_tm = dest_comb ltm in
  let float_tm, s_tm = dest_comb ltm in
    if (fst o dest_const) float_tm <> "float" then
      failwith "dest_float: not float"
    else
      (fst o dest_const) s_tm, n_tm, e_tm;;



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

let NUM_EXP_EXP' = SPEC_ALL NUM_EXP_EXP;;

let NUM_EXP_0' = (SPEC_ALL o REWRITE_RULE[NUMERAL]) NUM_EXP_0;;

let NUM_EXP_LE' = (UNDISCH_ALL o SPEC_ALL) NUM_EXP_LE;;

let NUM_EXP_LT' = (UNDISCH_ALL o SPEC_ALL) NUM_EXP_LT;;


(* NUM_EXP_DIV *)

let NUM_EXP_DIV1 = 
prove(`~(n2 = 0) /\ e2 <= e1 ==>
                           num_exp n1 e1 DIV num_exp n2 e2 = num_exp n1 (e1 - e2) DIV n2`,

   STRIP_TAC THEN
     (*`num_exp n1 e1 = 16 EXP e2 * num_exp n1 (e1 - e2)` MP_TAC THENL*)
     SUBGOAL_THEN (mk_eq(`num_exp n1 e1`, mk_binop mul_op_num (mk_binop `EXP` base_const `e2:num`) `num_exp n1 (e1 - e2)`)) MP_TAC THENL
     [
       REWRITE_TAC[num_exp] THEN
         ONCE_REWRITE_TAC[ARITH_RULE `a * b * c = b * (a * c:num)`] THEN
         REWRITE_TAC[GSYM EXP_ADD] THEN
         ASM_SIMP_TAC[ARITH_RULE `e2 <= e1 ==> e2 + e1 - e2 = e1:num`];
       ALL_TAC
     ] THEN
     DISCH_THEN (fun th -> REWRITE_TAC[th]) THEN
     SUBGOAL_THEN (mk_eq(`num_exp n2 e2`, mk_binop mul_op_num (mk_binop `EXP` base_const `e2:num`) `n2:num`)) MP_TAC THENL
     [
       REWRITE_TAC[num_exp; MULT_AC];
       ALL_TAC
     ] THEN
     DISCH_THEN (fun th -> REWRITE_TAC[th]) THEN
     MATCH_MP_TAC DIV_MULT2 THEN
     ASM_REWRITE_TAC[MULT_EQ_0; DE_MORGAN_THM; EXP_EQ_0] THEN
     ARITH_TAC);;


let NUM_EXP_DIV2 = 
prove(`~(n2 = 0) /\ e1 <= e2 ==>
                           num_exp n1 e1 DIV num_exp n2 e2 = n1 DIV num_exp n2 (e2 - e1)`,

   STRIP_TAC THEN
     (*`num_exp n2 e2 = 16 EXP e1 * num_exp n2 (e2 - e1)` MP_TAC THENL*)
     SUBGOAL_THEN (mk_eq(`num_exp n2 e2`, mk_binop mul_op_num (mk_binop `EXP` base_const `e1:num`) `num_exp n2 (e2 - e1)`)) MP_TAC THENL
     [
       REWRITE_TAC[num_exp] THEN
         ONCE_REWRITE_TAC[ARITH_RULE `a * b * c = b * (a * c:num)`] THEN
         REWRITE_TAC[GSYM EXP_ADD] THEN
         ASM_SIMP_TAC[ARITH_RULE `e1 <= e2 ==> e1 + e2 - e1 = e2:num`];
       ALL_TAC
     ] THEN
     DISCH_THEN (fun th -> REWRITE_TAC[th]) THEN
     SUBGOAL_THEN (mk_eq(`num_exp n1 e1`, mk_binop mul_op_num (mk_binop `EXP` base_const `e1:num`) `n1:num`)) MP_TAC THENL
     [
       REWRITE_TAC[num_exp; MULT_AC];
       ALL_TAC
     ] THEN
     DISCH_THEN (fun th -> REWRITE_TAC[th]) THEN
     MATCH_MP_TAC DIV_MULT2 THEN
     ASM_REWRITE_TAC[num_exp; MULT_EQ_0; DE_MORGAN_THM; EXP_EQ_0] THEN
     ARITH_TAC);;




(* B0 n = num_exp n bits *)
let normal_lemma1 = prove(mk_eq(mk_comb(b0_const, n_var_num), `num_exp n 1`),
   REWRITE_TAC[Arith_hash.def_array.(0); num_exp] THEN
     TRY ARITH_TAC THEN
     GEN_REWRITE_TAC (LAND_CONV o DEPTH_CONV) [BIT0] THEN
     ARITH_TAC);;


let NORMAL_LEMMA1 = NUMERALS_TO_NUM normal_lemma1;;

let normal_lemma2 = prove(mk_eq (mk_comb (b0_const, `num_exp n e`), `num_exp n (SUC e)`),
   REWRITE_TAC[normal_lemma1; NUM_EXP_EXP] THEN ARITH_TAC);;


let rec normalize tm =
  if (is_comb tm) then
    let ltm, rtm = dest_comb tm in
    let lname = (fst o dest_const) ltm in
      if (lname = b0_name) then
        let lth = INST[rtm, n_var_num] NORMAL_LEMMA1 in
        let rth, flag = normalize rtm in
          if flag then
            let ltm, lexp = (dest_comb o snd o dest_eq o concl) lth in
            let ltm, rtm = dest_comb ltm in
            let rn, rexp = (dest_comb o snd o dest_eq o concl) rth in
            let rn = rand rn in
            let th1 = AP_THM (AP_TERM ltm rth) lexp in
            let th2 = INST[rexp, e1_var_num; lexp, e2_var_num; rn, n_var_num] NUM_EXP_EXP' in
            let th3 = TRANS lth (TRANS th1 th2) in
            let ltm, rtm = (dest_comb o snd o dest_eq o concl) th3 in
            let add_th = raw_add_conv_hash rtm in
            let th4 = AP_TERM ltm add_th in
              (TRANS th3 th4, true)
          else
            (lth, true)
      else
        (REFL tm, false)
  else
    (REFL tm, false);;


(* Converts a raw numeral to a num_exp expression *)
let to_num_exp tm =
  let x, flag = normalize tm in
    if flag then x
    else
      INST[tm, n_var_num] NUM_EXP_0';;
	



(*
let tm = `B0 (B0 (B0 (B1 (B2 (B3 _0)))))`;;
(* 8: 0.568 *)
test 10000 normalize tm;; (* 0.568 *)
(* 8: 0.012 *)
test 10000 normalize `B2 (B3 (B4 _0))`;;
*)



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

let SYM_NUM_EXP_0' = SYM NUM_EXP_0';;

let NUM_EXP_n0 = 
prove(`!e. num_exp 0 e = 0`,
 REWRITE_TAC[num_exp; MULT_CLAUSES]);;
let NUM_EXP_n0' = (REWRITE_RULE[NUMERAL] o SPEC_ALL) NUM_EXP_n0;;


let NUM_EXP_DENORM = (UNDISCH_ALL o prove)
  (mk_imp(`e = _0 <=> F`, mk_eq(`num_exp n e`, mk_comb (b0_const, `num_exp n (PRE e)`))),
   REWRITE_TAC[] THEN ONCE_REWRITE_TAC[SYM (REWRITE_CONV[NUMERAL] `0`)] THEN
     REWRITE_TAC[num_exp; b0_thm] THEN
     REWRITE_TAC[ARITH_RULE (mk_eq(mk_binop mul_op_num base_const `n * a:num`,
                                   mk_binop mul_op_num `n:num` (mk_binop mul_op_num base_const `a:num`)))] THEN
     REWRITE_TAC[GSYM EXP] THEN
     SIMP_TAC[ARITH_RULE `~(e = 0) ==> SUC (PRE e) = e`]);;




(* Converts num_exp n e to a numeral by adding e B0's *)
let rec denormalize tm =
  let ltm, etm = dest_comb tm in
  let ntm = rand ltm in
    if (etm = zero_const) then
      INST[ntm, n_var_num] SYM_NUM_EXP_0'
    else
      if ntm = zero_const then
        INST[etm, e_var_num] NUM_EXP_n0'
      else
        let e_th = raw_eq0_hash_conv etm in
        let th0' = INST[etm, e_var_num; ntm, n_var_num] NUM_EXP_DENORM in
        let th0 = MY_PROVE_HYP e_th th0' in
        let b0_tm, rtm = dest_comb(rand(concl th0)) in
        let ltm, pre_tm = dest_comb rtm in
        let pre_th = raw_pre_hash_conv pre_tm in
        let th1 = AP_TERM ltm pre_th in
        let th2 = denormalize (rand(concl th1)) in
          TRANS th0 (AP_TERM b0_tm (TRANS th1 th2));;



(*
let tm = `num_exp (B1 (B3 (B5 _0))) (B15 _0)`;;
denormalize tm;;
(* 4: 0.264 *)
test 1000 denormalize tm;;
*)


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



let rec comb_number tm n =
  if (is_comb tm) then comb_number ((snd o dest_comb) tm) (n + 1) else n;;


let make_lo_thm i =
  let th_concl = mk_binop `(<=):num->num->bool`
    (mk_comb (Arith_hash.const_array.(0), n_var_num))
    (mk_comb (Arith_hash.const_array.(i), n_var_num)) in
    prove(th_concl,
          REWRITE_TAC[Arith_hash.def_array.(i); Arith_hash.def_array.(0)] THEN
            REWRITE_TAC[ARITH_LE; LE_REFL] THEN
            ARITH_TAC);;



let lo_thm_array = Array.init base make_lo_thm;;
let lo_thm_table = Hashtbl.create base;;

for i = 0 to base - 1 do
  Hashtbl.add lo_thm_table Arith_hash.const_array.(i) lo_thm_array.(i);
done;;



let make_lo_thm2 i =
  let th_concl = mk_imp (`n <= m:num`,
                         mk_binop `(<=):num->num->bool`
                          (mk_comb (Arith_hash.const_array.(0), n_var_num))
                          (mk_comb (Arith_hash.const_array.(i), m_var_num))) in
    (UNDISCH_ALL o prove) (th_concl,
          REWRITE_TAC[Arith_hash.def_array.(i); Arith_hash.def_array.(0); ARITH_LE] THEN
                          ARITH_TAC);;


let lo_thm2_array = Array.init base make_lo_thm2;;
let lo_thm2_table = Hashtbl.create base;;

for i = 0 to base - 1 do
  Hashtbl.add lo_thm2_table Arith_hash.const_array.(i) lo_thm2_array.(i);
done;;




let make_hi_thm i =
  let th_concl = mk_imp (`n < m:num`,
                         mk_binop `(<):num->num->bool`
                          (mk_comb (Arith_hash.const_array.(i), n_var_num))
                          (mk_comb (Arith_hash.const_array.(0), m_var_num))) in
    (UNDISCH_ALL o prove) (th_concl,
                           REWRITE_TAC[Arith_hash.def_array.(i); Arith_hash.def_array.(0); ARITH_LT] THEN
                             ARITH_TAC);;



let hi_thm_array = Array.init base make_hi_thm;;
let hi_thm_table = Hashtbl.create base;;

for i = 0 to base - 1 do
  Hashtbl.add hi_thm_table Arith_hash.const_array.(i) hi_thm_array.(i);
done;;



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


let LE_REFL' = SPEC_ALL LE_REFL;;
let LE_TRANS' = (UNDISCH_ALL o SPEC_ALL o REWRITE_RULE[GSYM IMP_IMP]) LE_TRANS;;


let lo_num_conv p tm =
  let n = comb_number tm 0 in
    if (n <= p) then
      INST[tm, n_var_num] LE_REFL'
    else
      let rec lo_bound n tm =
        let btm, rtm = dest_comb tm in
        let th0 = INST[rtm, n_var_num] (Hashtbl.find lo_thm_table btm) in
          if n > 1 then
            let rth = lo_bound (n - 1) rtm in
            let xtm = rand (rator (concl rth)) in
            let th1' = INST[xtm, n_var_num; rtm, m_var_num] (Hashtbl.find lo_thm2_table btm) in
            let th1 = MY_PROVE_HYP rth th1' in
              th1
          else
            th0 in

        lo_bound (n - p) tm;;



let N_LT_SUC = ARITH_RULE `n < SUC n`;;
let LT_IMP_LE' = (UNDISCH_ALL o SPEC_ALL) LT_IMP_LE;;
let N_LT_SUC = ARITH_RULE `n < SUC n`;;
let LT_LE_TRANS = (UNDISCH_ALL o ARITH_RULE) `n < e ==> e <= m ==> n < m:num`;;


(* Generates a theorem |- n <= m such that m contains at most p non-zero digits *)
let hi_num_conv p tm =
  let n = comb_number tm 0 in
    if (n <= p) then
      INST[tm, n_var_num] LE_REFL'
    else
      let k = n - p in

      let rec check_b0s n tm =
        let btm, rtm = dest_comb tm in
          if ((fst o dest_const) btm = b0_name) then
            if n > 1 then check_b0s (n - 1) rtm else true
          else
            false in

        if (check_b0s k tm) then
          INST[tm, n_var_num] LE_REFL'
        else
          let rec hi_bound n tm =
            if n > 0 then
              let btm, rtm = dest_comb tm in
              let r_th = hi_bound (n - 1) rtm in
              let xtm = rand (concl r_th) in
              let th0 = INST[rtm, n_var_num; xtm, m_var_num] (Hashtbl.find hi_thm_table btm) in
                MY_PROVE_HYP r_th th0
            else
              let th0 = INST[tm, n_var_num] N_LT_SUC in
              let ltm, suc_tm = dest_comb (concl th0) in
              let suc_th = raw_suc_conv_hash suc_tm in
                EQ_MP (AP_TERM ltm suc_th) th0 in

          let th = hi_bound k tm in
          let m_tm, l_tm = dest_comb (concl th) in
            MY_PROVE_HYP th (INST[rand m_tm, m_var_num; l_tm, n_var_num] LT_IMP_LE');;




(*
let tm = `B10 (B12 (B11 (B1 (B14 (B15 _0)))))`;;
hi_num_conv 4 tm;;

let tm = `B0 (B0 (B1 (B2 _0)))`;;
hi_num_conv 2 tm;;
*)



(* Generates a theorem |- n < m such that m contains at most p non-zero digits *)
let hi_lt_num_conv p tm =
  let n = comb_number tm 0 in
    if (n <= p) then
      let th0 = INST[tm, n_var_num] N_LT_SUC in
      let ltm, rtm = dest_comb(concl th0) in
      let suc_th = raw_suc_conv_hash rtm in
        EQ_MP (AP_TERM ltm suc_th) th0
    else
      let k = n - p in

      let rec check_b0s n tm =
        let btm, rtm = dest_comb tm in
          if ((fst o dest_const) btm = b0_name) then
            if n > 1 then check_b0s (n - 1) rtm else true
          else
            false in

        if (check_b0s k tm) then
          let th0 = INST[tm, n_var_num] N_LT_SUC in
          let ltm, rtm = dest_comb (concl th0) in
          let suc_th = raw_suc_conv_hash rtm in
          let suc_tm = rand(concl suc_th) in
          let th1 = hi_num_conv p suc_tm in
          let th2 = EQ_MP (AP_TERM ltm suc_th) th0 in
          let th = INST[tm, n_var_num; suc_tm, e_var_num; rand(concl th1), m_var_num] LT_LE_TRANS in
            MY_PROVE_HYP th1 (MY_PROVE_HYP th2 th)

        else
          let rec hi_bound n tm =
            if n > 0 then
              let btm, rtm = dest_comb tm in
              let r_th = hi_bound (n - 1) rtm in
              let xtm = rand (concl r_th) in
              let th0 = INST[rtm, n_var_num; xtm, m_var_num] (Hashtbl.find hi_thm_table btm) in
                MY_PROVE_HYP r_th th0
            else
              let th0 = INST[tm, n_var_num] N_LT_SUC in
              let ltm, suc_tm = dest_comb (concl th0) in
              let suc_th = raw_suc_conv_hash suc_tm in
                EQ_MP (AP_TERM ltm suc_th) th0 in
            hi_bound k tm;;



(*
let tm = `B10 (B12 (B11 (B1 (B14 (B15 _0)))))`;;
hi_num_conv 4 tm;;
hi_lt_num_conv 4 tm;;

let tm = `B0 (B0 (B1 (B2 _0)))`;;
hi_num_conv 3 tm;;
hi_lt_num_conv 3 tm;;
*)


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




let num_exp_lo p tm =
  let ltm, e_tm = dest_comb tm in
  let n_tm = rand ltm in
  let n_th = lo_num_conv p n_tm in
  let m_tm = rand (rator (concl n_th)) in
  let m_norm, flag = normalize m_tm in

  let th0' = INST[m_tm, m_var_num; n_tm, n_var_num; e_tm, e_var_num] NUM_EXP_LE' in
  let th0 = MY_PROVE_HYP n_th th0' in
    if flag then
      let th1 = AP_THM (AP_TERM (rator ltm) m_norm) e_tm in
      let m_tm, me_tm = (dest_comb o rand o concl) m_norm in
      let th2 = INST[me_tm, e1_var_num; e_tm, e2_var_num; rand m_tm, n_var_num] NUM_EXP_EXP' in
      let th3 = TRANS th1 th2 in
      let ltm, rtm = (dest_comb o rand o concl) th3 in
      let th_add = raw_add_conv_hash rtm in
      let th4 = TRANS th3 (AP_TERM ltm th_add) in
        EQ_MP (AP_THM (AP_TERM le_op_num th4) tm) th0


    else
      th0;;



let num_exp_hi p tm =
  let ltm, e_tm = dest_comb tm in
  let n_tm = rand ltm in
  let n_th = hi_num_conv p n_tm in
  let m_tm = rand (concl n_th) in
  let m_norm, flag = normalize m_tm in

  let th0' = INST[m_tm, n_var_num; n_tm, m_var_num; e_tm, e_var_num] NUM_EXP_LE' in
  let th0 = MY_PROVE_HYP n_th th0' in
    if flag then
      let th1 = AP_THM (AP_TERM (rator ltm) m_norm) e_tm in
      let m_tm, me_tm = (dest_comb o rand o concl) m_norm in
      let th2 = INST[me_tm, e1_var_num; e_tm, e2_var_num; rand m_tm, n_var_num] NUM_EXP_EXP' in
      let th3 = TRANS th1 th2 in
      let ltm, rtm = (dest_comb o rand o concl) th3 in
      let th_add = raw_add_conv_hash rtm in
      let th4 = TRANS th3 (AP_TERM ltm th_add) in
        EQ_MP (AP_TERM (rator (concl th0)) th4) th0
    else
      th0;;




let num_exp_hi_lt p tm =
  let ltm, e_tm = dest_comb tm in
  let n_tm = rand ltm in
  let n_th = hi_lt_num_conv p n_tm in
  let m_tm = rand (concl n_th) in
  let m_norm, flag = normalize m_tm in

  let th0' = INST[m_tm, n_var_num; n_tm, m_var_num; e_tm, e_var_num] NUM_EXP_LT' in
  let th0 = MY_PROVE_HYP n_th th0' in
    if flag then
      let th1 = AP_THM (AP_TERM (rator ltm) m_norm) e_tm in
      let m_tm, me_tm = (dest_comb o rand o concl) m_norm in
      let th2 = INST[me_tm, e1_var_num; e_tm, e2_var_num; rand m_tm, n_var_num] NUM_EXP_EXP' in
      let th3 = TRANS th1 th2 in
      let ltm, rtm = (dest_comb o rand o concl) th3 in
      let th_add = raw_add_conv_hash rtm in
      let th4 = TRANS th3 (AP_TERM ltm th_add) in
        EQ_MP (AP_TERM (rator (concl th0)) th4) th0
    else
      th0;;





(*
let tm = `num_exp (B5 (B1 (B2 (B3 _0)))) (B3 _0)`;;
num_exp_lo 2 tm;;
num_exp_hi 2 tm;;
num_exp_hi_lt 2 tm;;
num_exp_hi_lt 2 `num_exp (B0 (B0 (B1 _0))) (B4 _0)`;;
(* 4: 0.848 *)
test 10000 (num_exp_lo 2) tm;;
(* 4: 0.448 *)
test 10000 (num_exp_lo 3) tm;;
(* 4: 0.116 *)
test 10000 (num_exp_lo 5) tm;;
*)


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

(* num_exp_lt, num_exp_le *)

let transform = UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o SPEC_ALL;;

let NUM_EXP_LT1_EQ' = transform NUM_EXP_LT1_EQ and
    NUM_EXP_LT2_EQ' = transform NUM_EXP_LT2_EQ;;


let NUM_EXP_LE1_EQ' = transform NUM_EXP_LE1_EQ and
    NUM_EXP_LE2_EQ' = transform NUM_EXP_LE2_EQ;;


let num_exp_lt tm1 tm2 = 
  let n1_tm, e1_tm = dest_num_exp tm1 in
  let n2_tm, e2_tm = dest_num_exp tm2 in
  let sub_th, le_th = raw_sub_and_le_hash_conv e1_tm e2_tm in
  let r_tm = rand(concl sub_th) in
    if (rand(concl le_th) = e1_tm) then 
      let x_expr = mk_num_exp n1_tm r_tm in
      let x_th = denormalize x_expr in
      let x_tm = rand(concl x_th) in

      let th0 = INST[e2_tm, e2_var_num; e1_tm, e1_var_num;
		     r_tm, r_var_num; x_tm, x_var_num;
		     n1_tm, n1_var_num; n2_tm, n2_var_num] NUM_EXP_LT1_EQ' in
      let th1 = MY_PROVE_HYP x_th (MY_PROVE_HYP sub_th (MY_PROVE_HYP le_th th0)) in
      let lt_th = raw_lt_hash_conv (rand(concl th1)) in
	TRANS th1 lt_th
    else
      let x_expr = mk_num_exp n2_tm r_tm in
      let x_th = denormalize x_expr in
      let x_tm = rand(concl x_th) in

      let th0 = INST[e2_tm, e2_var_num; e1_tm, e1_var_num;
		     r_tm, r_var_num; x_tm, x_var_num;
		     n1_tm, n1_var_num; n2_tm, n2_var_num] NUM_EXP_LT2_EQ' in
      let th1 = MY_PROVE_HYP x_th (MY_PROVE_HYP sub_th (MY_PROVE_HYP le_th th0)) in
      let lt_th = raw_lt_hash_conv (rand(concl th1)) in
	TRANS th1 lt_th;;



let num_exp_le tm1 tm2 = 
  let n1_tm, e1_tm = dest_num_exp tm1 in
  let n2_tm, e2_tm = dest_num_exp tm2 in
  let sub_th, le_th = raw_sub_and_le_hash_conv e1_tm e2_tm in
  let r_tm = rand(concl sub_th) in
    if (rand(concl le_th) = e1_tm) then 
      let x_expr = mk_num_exp n1_tm r_tm in
      let x_th = denormalize x_expr in
      let x_tm = rand(concl x_th) in

      let th0 = INST[e2_tm, e2_var_num; e1_tm, e1_var_num;
		     r_tm, r_var_num; x_tm, x_var_num;
		     n1_tm, n1_var_num; n2_tm, n2_var_num] NUM_EXP_LE1_EQ' in
      let th1 = MY_PROVE_HYP x_th (MY_PROVE_HYP sub_th (MY_PROVE_HYP le_th th0)) in
      let le_th = raw_le_hash_conv (rand(concl th1)) in
	TRANS th1 le_th
    else
      let x_expr = mk_num_exp n2_tm r_tm in
      let x_th = denormalize x_expr in
      let x_tm = rand(concl x_th) in

      let th0 = INST[e2_tm, e2_var_num; e1_tm, e1_var_num;
		     r_tm, r_var_num; x_tm, x_var_num;
		     n1_tm, n1_var_num; n2_tm, n2_var_num] NUM_EXP_LE2_EQ' in
      let th1 = MY_PROVE_HYP x_th (MY_PROVE_HYP sub_th (MY_PROVE_HYP le_th th0)) in
      let le_th = raw_le_hash_conv (rand(concl th1)) in
	TRANS th1 le_th;;




(*
let tm1 = `num_exp (D0 (D4 (D1 _0))) (D2 (D3 _0))`;;
let tm2 = `num_exp (D4 (D6 _0)) (D5 (D3 _0))`;;

num_exp_lt tm1 tm2;;
num_exp_lt tm2 tm1;;
num_exp_le tm1 tm1;;
num_exp_lt tm1 tm1;;

(* 100: 0.164 *)
test 1000 (num_exp_lt tm1) tm2;;
*)



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

(* num_exp_mul *)

let NUM_EXP_MUL' = SPEC_ALL NUM_EXP_MUL;;

let num_exp_mul tm1 tm2 =
  let n1_tm, e1_tm = dest_comb tm1 in
  let n1_tm = rand n1_tm in
  let n2_tm, e2_tm = dest_comb tm2 in
  let n2_tm = rand n2_tm in
  let th0 = INST[n1_tm, n1_var_num; e1_tm, e1_var_num;
                 n2_tm, n2_var_num; e2_tm, e2_var_num] NUM_EXP_MUL' in
  let ltm, tm_add = dest_comb (rand (concl th0)) in
  let tm_mul = rand ltm in
  let th_mul = raw_mul_conv_hash tm_mul in
  let th_add = raw_add_conv_hash tm_add in
    TRANS th0 (MK_COMB (AP_TERM (rator ltm) th_mul, th_add));;


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

(* num_exp_add *)

let NUM_EXP_ADD' = (UNDISCH_ALL o SPEC_ALL) NUM_EXP_ADD;;
let ADD_COMM = ARITH_RULE `m + n = n + m:num`;;
let num_exp_add tm1 tm2 =
  let n1_tm, e1_tm = dest_comb tm1 in
  let n1_tm = rand n1_tm in
  let n2_tm, e2_tm = dest_comb tm2 in
  let n2_tm = rand n2_tm in
  let e_sub, e_le = raw_sub_and_le_hash_conv e1_tm e2_tm in

  let flag = (rand(concl e_le) = e2_tm) in

  let th0' =
    if flag then
      INST[n1_tm, n1_var_num; e1_tm, e1_var_num;
           n2_tm, n2_var_num; e2_tm, e2_var_num] NUM_EXP_ADD'
    else
      INST[n2_tm, n1_var_num; e2_tm, e1_var_num;
           n1_tm, n2_var_num; e1_tm, e2_var_num] NUM_EXP_ADD' in

  let th0 = MY_PROVE_HYP e_le th0' in
  let ltm, e0_tm = dest_comb(rand(concl th0)) in
  let exp_tm, add_tm = dest_comb ltm in
  let ltm, d_tm = dest_comb add_tm in
  let th1 = AP_TERM (rator d_tm) e_sub in
  let th2 = denormalize (rand(concl th1)) in
  let th3 = AP_TERM ltm (TRANS th1 th2) in
  let th4 = raw_add_conv_hash (rand(concl th3)) in
  let th5 = AP_THM (AP_TERM exp_tm (TRANS th3 th4)) e0_tm in
  let th = TRANS th0 th5 in
    if flag then th else
      TRANS (INST[tm1, m_var_num; tm2, n_var_num] ADD_COMM) th;;


(*
let tm1 = `num_exp (B0 (B4 (B1 _0))) (B1 (B3 _0))`;;
let tm2 = `num_exp (B4 (B6 _0)) (B1 (B2 _0))`;;
num_exp_add tm1 tm2;;
num_exp_add tm2 tm1;;
num_exp_mul tm1 tm2;;
(* 4: 0.836 *)
test 10000 (num_exp_mul tm1) tm2;;
(* 4: 4.112 *)
test 10000 (num_exp_add tm1) tm2;;
*)


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

(* num_exp_sub *)

let NUM_EXP_SUB1' = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o SPEC_ALL) NUM_EXP_SUB1 and
    NUM_EXP_SUB2' = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o SPEC_ALL) NUM_EXP_SUB2 and
    NUM_EXP_LE1' = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o SPEC_ALL) NUM_EXP_LE1 and
    NUM_EXP_LE2' = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o SPEC_ALL) NUM_EXP_LE2;;


(* Returns two theorems: |- tm1 - tm2 = tm, |- tm2 <= tm1 or
   |- tm2 - tm1 = tm, |- tm1 <= tm2 *)
let num_exp_sub tm1 tm2 =
  let n1_tm, e1_tm = dest_num_exp tm1 in
  let n2_tm, e2_tm = dest_num_exp tm2 in
  let e_sub, e_le = raw_sub_and_le_hash_conv e1_tm e2_tm in

    if rand(concl e_le) = e1_tm then
      (* e2 <= e1 *)
      let e1_sub_e2 = rand(concl e_sub) in
      let a0 = mk_num_exp n1_tm e1_sub_e2 in
      let b = n2_tm in
      let a_th = denormalize a0 in
      let a = rand(concl a_th) in

      let th_sub, th_le = raw_sub_and_le_hash_conv a b in
	if rand(concl th_le) = a then
	  (* b <= a *)
	  let a_sub_b = TRANS (AP_THM (AP_TERM minus_op_num a_th) b) th_sub in
	  let b_le_a = EQ_MP (SYM (AP_TERM (rator(concl th_le)) a_th)) th_le in
	  let th0 = AP_THM (AP_TERM num_exp_const a_sub_b) e2_tm in

	  let inst = INST[n1_tm, n1_var_num; e1_tm, e1_var_num;
			  n2_tm, n2_var_num; e2_tm, e2_var_num; e1_sub_e2, r_var_num] in
	  let th1_sub = inst NUM_EXP_SUB1' in
	  let th1_le = inst NUM_EXP_LE1' in
	  let th2_sub = MY_PROVE_HYP e_sub (MY_PROVE_HYP e_le th1_sub) in
	  let th2_le = MY_PROVE_HYP e_sub (MY_PROVE_HYP b_le_a (MY_PROVE_HYP e_le th1_le)) in
	    TRANS th2_sub th0, th2_le

	else
	  (* a <= b *)
	  let b_sub_a = TRANS (AP_TERM (rator(lhand(concl th_sub))) a_th) th_sub in
	  let a_le_b = EQ_MP (SYM (AP_THM (AP_TERM le_op_num a_th) b)) th_le in
	  let th0 = AP_THM (AP_TERM num_exp_const b_sub_a) e2_tm in
	  let inst = INST[n2_tm, n1_var_num; e2_tm, e1_var_num;
			  n1_tm, n2_var_num; e1_tm, e2_var_num; e1_sub_e2, r_var_num] in
	  let th1_sub = inst NUM_EXP_SUB2' in
	  let th1_le = inst NUM_EXP_LE2' in
	  let th2_sub = MY_PROVE_HYP e_sub (MY_PROVE_HYP e_le th1_sub) in
	  let th2_le = MY_PROVE_HYP e_sub (MY_PROVE_HYP a_le_b (MY_PROVE_HYP e_le th1_le)) in
	    TRANS th2_sub th0, th2_le

    else
      (* e1 <= e2 *)
      let e2_sub_e1 = rand(concl e_sub) in
      let b0 = mk_num_exp n2_tm e2_sub_e1 in
      let a = n1_tm in
      let b_th = denormalize b0 in
      let b = rand(concl b_th) in

      let th_sub, th_le = raw_sub_and_le_hash_conv a b in
	if rand(concl th_le) = a then
	  (* b <= a *)
	  let a_sub_b = TRANS (AP_TERM (rator(lhand(concl th_sub))) b_th) th_sub in
	  let b_le_a = EQ_MP (SYM (AP_THM (AP_TERM le_op_num b_th) a)) th_le in
	  let th0 = AP_THM (AP_TERM num_exp_const a_sub_b) e1_tm in
	  let inst = INST[n1_tm, n1_var_num; e1_tm, e1_var_num;
			  n2_tm, n2_var_num; e2_tm, e2_var_num; e2_sub_e1, r_var_num] in
	  let th1_sub = inst NUM_EXP_SUB2' in
	  let th1_le = inst NUM_EXP_LE2' in
	  let th2_sub = MY_PROVE_HYP e_sub (MY_PROVE_HYP e_le th1_sub) in
	  let th2_le = MY_PROVE_HYP e_sub (MY_PROVE_HYP b_le_a (MY_PROVE_HYP e_le th1_le)) in
	    TRANS th2_sub th0, th2_le

	else
	  (* a <= b *)
	  let b_sub_a = TRANS (AP_THM (AP_TERM minus_op_num b_th) a) th_sub in
	  let a_le_b = EQ_MP (SYM (AP_TERM (rator(concl th_le)) b_th)) th_le in
	  let th0 = AP_THM (AP_TERM num_exp_const b_sub_a) e1_tm in
	  let inst = INST[n2_tm, n1_var_num; e2_tm, e1_var_num;
			  n1_tm, n2_var_num; e1_tm, e2_var_num; e2_sub_e1, r_var_num] in
	  let th1_sub = inst NUM_EXP_SUB1' in
	  let th1_le = inst NUM_EXP_LE1' in
	  let th2_sub = MY_PROVE_HYP e_sub (MY_PROVE_HYP e_le th1_sub) in
	  let th2_le = MY_PROVE_HYP e_sub (MY_PROVE_HYP a_le_b (MY_PROVE_HYP e_le th1_le)) in
	    TRANS th2_sub th0, th2_le;;



  
(*
let tm1 = `num_exp (D3 (D4 (D1 _0))) (D2 _0)`;;
let tm2 = `num_exp (D1 (D2 (D1 _0))) (D3 _0)`;;
num_exp_sub tm1 tm2;;
(* 10: 0.252 *)
test 1000 (num_exp_sub tm1) tm2;;
*)




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

(* division *)

let NUM_EXP_DIV1' = (UNDISCH_ALL o PURE_REWRITE_RULE[NUMERAL] o
                       PURE_ONCE_REWRITE_RULE[ARITH_RULE `~(x = 0) <=> (x = 0 <=> F)`] o
                       REWRITE_RULE[GSYM IMP_IMP]) NUM_EXP_DIV1;;
let NUM_EXP_DIV2' = (UNDISCH_ALL o PURE_REWRITE_RULE[NUMERAL] o
                       PURE_ONCE_REWRITE_RULE[ARITH_RULE `~(x = 0) <=> (x = 0 <=> F)`] o
                       REWRITE_RULE[GSYM IMP_IMP]) NUM_EXP_DIV2;;


let num_exp_div tm1 tm2 =
  let n1_tm, e1_tm = dest_comb tm1 in
  let n1_tm = rand n1_tm in
  let n2_tm, e2_tm = dest_comb tm2 in
  let n2_tm = rand n2_tm in
  let e_sub, e_le = raw_sub_and_le_hash_conv e1_tm e2_tm in

  let inst = INST[n1_tm, n1_var_num; e1_tm, e1_var_num;
                  n2_tm, n2_var_num; e2_tm, e2_var_num] in

  let n2_not_0 = raw_eq0_hash_conv n2_tm in
    if ((fst o dest_const o rand o concl) n2_not_0 = "T") then
      failwith "num_exp_div: n2 = 0"
    else
      if (rand(concl e_le) = e1_tm) then
         let th0' = inst NUM_EXP_DIV1' in
        let th0 = MY_PROVE_HYP n2_not_0 (MY_PROVE_HYP e_le th0') in

        let ltm, rtm = dest_comb(rand(concl th0)) in
        let div_tm, rtm2 = dest_comb ltm in
        let num_exp_tm = rator rtm2 in

        let th1 = AP_THM (AP_TERM div_tm (AP_TERM num_exp_tm e_sub)) rtm in
        let ltm, rtm = dest_comb(rand(concl th1)) in
        let tm1 = rand ltm in

        let th2 = AP_THM (AP_TERM div_tm (denormalize tm1)) rtm in
        let th3 = raw_div_hash_conv (rand(concl th2)) in
        let th = TRANS th0 (TRANS th1 (TRANS th2 th3)) in
          TRANS th (INST[rand(concl th), n_var_num] NUM_EXP_0')

      else
        let th0' = inst NUM_EXP_DIV2' in
        let th0 = MY_PROVE_HYP n2_not_0 (MY_PROVE_HYP e_le th0') in

        let ltm, rtm = dest_comb(rand(concl th0)) in
        let num_exp_tm = rator rtm in
        let th1 = AP_TERM ltm (AP_TERM num_exp_tm e_sub) in

        let ltm, rtm = dest_comb(rand(concl th1)) in
        let th2 = AP_TERM ltm (denormalize rtm) in
        let th3 = raw_div_hash_conv (rand(concl th2)) in
        let th = TRANS th0 (TRANS th1 (TRANS th2 th3)) in
          TRANS th (INST[rand(concl th), n_var_num] NUM_EXP_0');;




(*
let tm1 = `num_exp (B0 (B4 (B1 _0))) (B1 (B3 _0))`;;
let tm2 = `num_exp (B4 (B6 _0)) (B1 (B2 _0))`;;
num_exp_div tm1 tm2;;
num_exp_div tm2 tm1;;
(* 4: 1.768 *)
test 1000 (num_exp_div tm1) tm2;;
(* 4: 0.656 *)
test 1000 (num_exp_div tm2) tm1;;
*)

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

(* Computes the lower bound for (op tm1 tm2) *)
let num_exp_op_lo p op tm1 tm2 =
  let op_th = op tm1 tm2 in
  let rtm = rand (concl op_th) in
  let lo_th = num_exp_lo p rtm in
  let ltm = rator (concl lo_th) in
  let th0 = AP_TERM ltm op_th in
    EQ_MP (SYM th0) lo_th;;



let num_exp_op_hi p op tm1 tm2 =
  let op_th = op tm1 tm2 in
  let rtm = rand (concl op_th) in
  let hi_th = num_exp_hi p rtm in
  let tm = rand (concl hi_th) in
  let th0 = AP_THM (AP_TERM le_op_num op_th) tm in
    EQ_MP (SYM th0) hi_th;;



let num_exp_op_hi_lt p op tm1 tm2 =
  let op_th = op tm1 tm2 in
  let rtm = rand (concl op_th) in
  let hi_lt_th = num_exp_hi_lt p rtm in
  let tm = rand (concl hi_lt_th) in
  let th0 = AP_THM (AP_TERM lt_op_num op_th) tm in
    EQ_MP (SYM th0) hi_lt_th;;


(*
num_exp_op_lo 1 num_exp_add tm1 tm2;;
num_exp_op_hi 1 num_exp_div tm1 tm2;;
num_exp_op_hi_lt 1 num_exp_div tm1 tm2;;
*)


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


(* float *)

let mod_plus = new_definition `mod_plus s1 s2 = (~(s1 /\ s2) /\ (s1 \/ s2))`;;



let EXP_INV_lemma = 
prove(`!n e1 e2. ~(n = 0) /\ e2 <= e1 ==> &(n EXP (e1 - e2)) =
                              &(n EXP e1) * inv(&(n EXP e2))`,

   REPEAT STRIP_TAC THEN
     REWRITE_TAC[GSYM REAL_OF_NUM_POW] THEN
     MP_TAC (SPECL [`&n`; `e2:num`; `e1:num`] REAL_POW_SUB) THEN
     ASM_REWRITE_TAC[REAL_OF_NUM_EQ; real_div]);;

let NUM_EXP_SUB_lemma = 
prove(`!n e1 e2. e2 <= e1 ==> &(num_exp n (e1 - e2)) =
                                  &(num_exp n e1) * inv(&(num_exp 1 e2))`,

   REPEAT STRIP_TAC THEN
     REWRITE_TAC[num_exp] THEN
     REWRITE_TAC[GSYM REAL_OF_NUM_MUL] THEN
     MP_TAC (SPECL [base_const; `e1:num`; `e2:num`] EXP_INV_lemma) THEN
     ANTS_TAC THENL
     [
       ASM_REWRITE_TAC[] THEN ARITH_TAC;
       ALL_TAC
     ] THEN
     DISCH_THEN (fun th -> REWRITE_TAC[th]) THEN
     REAL_ARITH_TAC);;

(********************)
(* Float operations *)
(********************)

module Float_ops = struct


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

(* FLOAT_LT *)

let FLOAT_LT_FF = 
prove(`float F n1 e1 < float F n2 e2 <=> num_exp n1 e1 < num_exp n2 e2`,

			REWRITE_TAC[float; GSYM REAL_OF_NUM_LT; REAL_MUL_LID; real_div] THEN
			  MATCH_MP_TAC REAL_LT_RMUL_EQ THEN
			  MATCH_MP_TAC REAL_LT_INV THEN
			  REWRITE_TAC[REAL_OF_NUM_LT; LT_NZ; NUM_EXP_EQ_0] THEN
			  ARITH_TAC);;
let FLOAT_LT_TT = 
prove(`float T n1 e1 < float T n2 e2 <=> num_exp n2 e2 < num_exp n1 e1`,

			REWRITE_TAC[FLOAT_NEG_T; REAL_ARITH `--a < --b <=> b < a`] THEN
			  REWRITE_TAC[FLOAT_LT_FF]);;

let FLOAT_LT_FT = 
prove(`float F n1 e1 < float T n2 e2 <=> F`,

			MP_TAC (SPECL [`n2:num`; `e2:num`] FLOAT_T_NEG) THEN
			  MP_TAC (SPECL [`n1:num`; `e1:num`] FLOAT_F_POS) THEN
			  REAL_ARITH_TAC);;

let FLOAT_LT_TF_00 = (PURE_REWRITE_RULE[NUMERAL] o prove)
  (`float T 0 e1 < float F 0 e2 <=> F`,
   MP_TAC (SPECL [`T`; `0`; `e1:num`] FLOAT_EQ_0) THEN
     MP_TAC (SPECL [`F`; `0`; `e2:num`] FLOAT_EQ_0) THEN
     REWRITE_TAC[] THEN
     REPLICATE_TAC 2 (DISCH_THEN (fun th -> REWRITE_TAC[th])) THEN
     REAL_ARITH_TAC);;


let FLOAT_LT_TF_1 = (UNDISCH_ALL o PURE_REWRITE_RULE[NUMERAL] o prove)
  (`(n1 = 0 <=> F) ==> (float T n1 e1 < float F n2 e2 <=> T)`,
   DISCH_TAC THEN
     MATCH_MP_TAC (REAL_ARITH `a < &0 /\ &0 <= b ==> (a < b <=> T)`) THEN
     REWRITE_TAC[FLOAT_F_POS] THEN
     MATCH_MP_TAC (REAL_ARITH `~(a = &0) /\ a <= &0 ==> a < &0`) THEN
     ASM_REWRITE_TAC[FLOAT_T_NEG; FLOAT_EQ_0]);;


let FLOAT_LT_TF_2 = (UNDISCH_ALL o PURE_REWRITE_RULE[NUMERAL] o prove)
  (`(n2 = 0 <=> F) ==> (float T n1 e1 < float F n2 e2 <=> T)`,
   DISCH_TAC THEN
     MATCH_MP_TAC (REAL_ARITH `a <= &0 /\ &0 < b ==> (a < b <=> T)`) THEN
     REWRITE_TAC[FLOAT_T_NEG] THEN
     MATCH_MP_TAC (REAL_ARITH `~(a = &0) /\ &0 <= a ==> &0 < a`) THEN
     ASM_REWRITE_TAC[FLOAT_F_POS; FLOAT_EQ_0]);;


let FLOAT_F_LT_0 = 
prove(`float F n e < &0 <=> F`,

			 MP_TAC (SPEC_ALL FLOAT_F_POS) THEN
			   REAL_ARITH_TAC);;
let FLOAT_T_LT_0 = (CONV_RULE (RAND_CONV (REWRITE_CONV[NUMERAL])) o prove)
  (`float T n e < &0 <=> (0 < n)`,
   REWRITE_TAC[REAL_ARITH `a < &0 <=> a <= &0 /\ ~(a = &0)`] THEN
     REWRITE_TAC[FLOAT_T_NEG; FLOAT_EQ_0] THEN
     ARITH_TAC);;


let FLOAT_F_GT_0 = (CONV_RULE (RAND_CONV (REWRITE_CONV[NUMERAL])) o prove)
  (`&0 < float F n e <=> 0 < n`,
   REWRITE_TAC[REAL_ARITH `&0 < a <=> &0 <= a /\ ~(a = &0)`] THEN
     REWRITE_TAC[FLOAT_F_POS; FLOAT_EQ_0] THEN
     ARITH_TAC);;

let FLOAT_T_GT_0 = 
prove(`&0 < float T n e <=> F`,

			 MP_TAC (SPEC_ALL FLOAT_T_NEG) THEN
			   REAL_ARITH_TAC);;

(* float_lt0, float_gt0 *)

let float_lt0 f1 = 
  let s, n_tm, e_tm = dest_float f1 in
  let inst = INST[n_tm, n_var_num; e_tm, e_var_num] in
    if s = "F" then
      inst FLOAT_F_LT_0
    else
      let gt_th = raw_gt0_hash_conv n_tm in
	TRANS (inst FLOAT_T_LT_0) gt_th;;


let float_gt0 f1 =
  let s, n_tm, e_tm = dest_float f1 in
  let inst = INST[n_tm, n_var_num; e_tm, e_var_num] in
    if s = "F" then
      let gt_th = raw_gt0_hash_conv n_tm in
	TRANS (inst FLOAT_F_GT_0) gt_th
    else
      inst FLOAT_T_GT_0;;



(* float_lt *)

let float_lt f1 f2 = 
  let s1, n1, e1 = dest_float f1 in
  let s2, n2, e2 = dest_float f2 in
  let inst = INST[n1, n1_var_num; e1, e1_var_num;
		  n2, n2_var_num; e2, e2_var_num] in
    if s1 = "F" then
      if s2 = "F" then
	(* FF *)
	let th0 = inst FLOAT_LT_FF in
	let ltm, tm2 = dest_comb (rand (concl th0)) in
	let lt_th = num_exp_lt (rand ltm) tm2 in
	  TRANS th0 lt_th
      else
	(* FT *)
	inst FLOAT_LT_FT
    else
      if s2 = "F" then
	(* TF *)
	if (is_const n1 && is_const n2) then
	  (* n1 = _0 and n2 = _0 *)
	  inst FLOAT_LT_TF_00
	else
	  let n1_0 = raw_eq0_hash_conv n1 in
	    if (fst o dest_const o rand o concl) n1_0 = "F" then
	      (* n1 <> _0 *)
	      MY_PROVE_HYP n1_0 (inst FLOAT_LT_TF_1)
	    else
	      let n2_0 = raw_eq0_hash_conv n2 in
		if (fst o dest_const o rand o concl) n2_0 = "F" then 
		  (* n2 <> _0 *)
		  MY_PROVE_HYP n2_0 (inst FLOAT_LT_TF_2)
		else
		  failwith "float_lt: D0 _0 exception"
      else
	(* TT *)
	let th0 = inst FLOAT_LT_TT in
	let ltm, tm2 = dest_comb (rand (concl th0)) in
	let lt_th = num_exp_lt (rand ltm) tm2 in
	  TRANS th0 lt_th;;

(*
let f1 = `float T (_0) (D5 _0)`;;
let f2 = `float F (D1 _0) (D3 _0)`;;
float_lt f1 f2;;
*)


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

(* FLOAT_LE *)

let FLOAT_LE_FF = 
prove(`float F n1 e1 <= float F n2 e2 <=> num_exp n1 e1 <= num_exp n2 e2`,

			REWRITE_TAC[float; GSYM REAL_OF_NUM_LE; REAL_MUL_LID; real_div] THEN
			  MATCH_MP_TAC REAL_LE_RMUL_EQ THEN
			  MATCH_MP_TAC REAL_LT_INV THEN
			  REWRITE_TAC[REAL_OF_NUM_LT; LT_NZ; NUM_EXP_EQ_0] THEN
			  ARITH_TAC);;
let FLOAT_LE_TT = 
prove(`float T n1 e1 <= float T n2 e2 <=> num_exp n2 e2 <= num_exp n1 e1`,

			REWRITE_TAC[FLOAT_NEG_T; REAL_ARITH `--a <= --b <=> b <= a`] THEN
			  REWRITE_TAC[FLOAT_LE_FF]);;

let FLOAT_LE_TF = 
prove(`float T n1 e1 <= float F n2 e2 <=> T`,

			MP_TAC (SPECL [`n1:num`; `e1:num`] FLOAT_T_NEG) THEN
			  MP_TAC (SPECL [`n2:num`; `e2:num`] FLOAT_F_POS) THEN
			  REAL_ARITH_TAC);;
			


let FLOAT_LE_FT = 
prove(`float F n1 e1 <= float T n2 e2 <=> n1 = 0 /\ n2 = 0`,

			REWRITE_TAC[REAL_LE_LT; FLOAT_LT_FT] THEN EQ_TAC THENL
			  [
			    DISCH_TAC THEN SUBGOAL_THEN `float F n1 e1 = &0 /\ float T n2 e2 = &0` MP_TAC THENL
			      [
				MP_TAC (SPECL [`n2:num`; `e2:num`] FLOAT_T_NEG) THEN
				  MP_TAC (SPECL [`n1:num`; `e1:num`] FLOAT_F_POS) THEN
				  ASM_REWRITE_TAC[] THEN REAL_ARITH_TAC;
				ALL_TAC
			      ] THEN
			      REWRITE_TAC[FLOAT_EQ_0];
			    DISCH_THEN (fun th -> REWRITE_TAC[th]) THEN
			      REWRITE_TAC[float; NUM_EXP_n0; real_div; REAL_MUL_LZERO; REAL_MUL_RZERO]
			  ]);;
			

let FLOAT_LE_FT_00 = (PURE_REWRITE_RULE[NUMERAL] o prove)
  (`float F 0 e1 <= float T 0 e2 <=> T`, REWRITE_TAC[FLOAT_LE_FT]);;

let FLOAT_LE_FT_1 = (UNDISCH_ALL o PURE_REWRITE_RULE[NUMERAL] o prove)
  (`(n1 = 0 <=> F) ==> (float F n1 e1 <= float T n2 e2 <=> F)`,
   DISCH_TAC THEN ASM_REWRITE_TAC[FLOAT_LE_FT]);;

let FLOAT_LE_FT_2 = (UNDISCH_ALL o PURE_REWRITE_RULE[NUMERAL] o prove)
  (`(n2 = 0 <=> F) ==> (float F n1 e1 <= float T n2 e2 <=> F)`,
   DISCH_TAC THEN ASM_REWRITE_TAC[FLOAT_LE_FT]);;


let FLOAT_F_LE_0 = (CONV_RULE (RAND_CONV (REWRITE_CONV[NUMERAL])) o prove)
  (`float F n e <= &0 <=> n = 0`,
   REWRITE_TAC[GSYM (SPEC `F` FLOAT_EQ_0)] THEN
     MP_TAC (SPEC_ALL FLOAT_F_POS) THEN
     REAL_ARITH_TAC);;

let FLOAT_T_LE_0 = 
prove(`float T n e <= &0 <=> T`,
 REWRITE_TAC[FLOAT_T_NEG]);;

let FLOAT_F_GE_0 = 
prove(`&0 <= float F n e <=> T`,
 REWRITE_TAC[FLOAT_F_POS]);;
let FLOAT_T_GE_0 = (CONV_RULE (RAND_CONV (REWRITE_CONV[NUMERAL])) o prove)
  (`&0 <= float T n e <=> n = 0`,
   REWRITE_TAC[GSYM (SPEC `T` FLOAT_EQ_0)] THEN
     MP_TAC (SPEC_ALL FLOAT_T_NEG) THEN
     REAL_ARITH_TAC);;


(* float_le0, float_ge0 *)

let float_le0 f1 = 
  let s, n_tm, e_tm = dest_float f1 in
  let inst = INST[n_tm, n_var_num; e_tm, e_var_num] in
    if s = "T" then
      inst FLOAT_T_LE_0
    else
      let eq_th = raw_eq0_hash_conv n_tm in
	TRANS (inst FLOAT_F_LE_0) eq_th;;


let float_ge0 f1 =
  let s, n_tm, e_tm = dest_float f1 in
  let inst = INST[n_tm, n_var_num; e_tm, e_var_num] in
    if s = "T" then
      let eq_th = raw_eq0_hash_conv n_tm in
	TRANS (inst FLOAT_T_GE_0) eq_th
    else
      inst FLOAT_F_GE_0;;


(* float_le *)

let float_le f1 f2 = 
  let s1, n1, e1 = dest_float f1 in
  let s2, n2, e2 = dest_float f2 in
  let inst = INST[n1, n1_var_num; e1, e1_var_num;
		  n2, n2_var_num; e2, e2_var_num] in
    if s2 = "F" then
      if s1 = "F" then
	(* FF *)
	let th0 = inst FLOAT_LE_FF in
	let ltm, tm2 = dest_comb (rand (concl th0)) in
	let le_th = num_exp_le (rand ltm) tm2 in
	  TRANS th0 le_th
      else
	(* TF *)
	inst FLOAT_LE_TF
    else
      if s1 = "F" then
	(* FT *)
	if (is_const n1 && is_const n2) then
	  (* n1 = _0 and n2 = _0 *)
	  inst FLOAT_LE_FT_00
	else
	  let n1_0 = raw_eq0_hash_conv n1 in
	    if (fst o dest_const o rand o concl) n1_0 = "F" then
	      (* n1 <> _0 *)
	      MY_PROVE_HYP n1_0 (inst FLOAT_LE_FT_1)
	    else
	      let n2_0 = raw_eq0_hash_conv n2 in
		if (fst o dest_const o rand o concl) n2_0 = "F" then 
		  (* n2 <> _0 *)
		  MY_PROVE_HYP n2_0 (inst FLOAT_LE_FT_2)
		else
		  failwith "float_lt: D0 _0 exception"
      else
	(* TT *)
	let th0 = inst FLOAT_LE_TT in
	let ltm, tm2 = dest_comb (rand (concl th0)) in
	let le_th = num_exp_le (rand ltm) tm2 in
	  TRANS th0 le_th;;


(*
let f1 = `float F (D1 _0) (D5 _0)`;;
let f2 = `float T (D2 _0) (D3 _0)`;;
float_le f1 f2;;
*)


(*************************************)
(* float_max, float_min *)

let FLOAT_MIN_1 = (UNDISCH_ALL o prove)(`(f1 <= f2 <=> T) ==> min f1 f2 = f1`, REAL_ARITH_TAC);;
let FLOAT_MIN_2 = (UNDISCH_ALL o prove)(`(f1 <= f2 <=> F) ==> min f1 f2 = f2`, REAL_ARITH_TAC);;

let FLOAT_MAX_1 = (UNDISCH_ALL o prove)(`(f1 <= f2 <=> T) ==> max f1 f2 = f2`, REAL_ARITH_TAC);;
let FLOAT_MAX_2 = (UNDISCH_ALL o prove)(`(f1 <= f2 <=> F) ==> max f1 f2 = f1`, REAL_ARITH_TAC);;


let float_min f1 f2 =
  let inst = INST[f1, f1_var_real; f2, f2_var_real] in
  let le_th = float_le f1 f2 in
  let th0 =
    if (fst o dest_const o rand o concl) le_th = "T" then
      inst FLOAT_MIN_1
    else
      inst FLOAT_MIN_2 in
    MY_PROVE_HYP le_th th0;;


let float_max f1 f2 =
  let inst = INST[f1, f1_var_real; f2, f2_var_real] in
  let le_th = float_le f1 f2 in
  let th0 =
    if (fst o dest_const o rand o concl) le_th = "T" then
      inst FLOAT_MAX_1
    else
      inst FLOAT_MAX_2 in
    MY_PROVE_HYP le_th th0;;


let float_min_max f1 f2 =
  let inst = INST[f1, f1_var_real; f2, f2_var_real] in
  let le_th = float_le f1 f2 in
  let th_min, th_max =
    if (fst o dest_const o rand o concl) le_th = "T" then
      inst FLOAT_MIN_1, inst FLOAT_MAX_1
    else
      inst FLOAT_MIN_2, inst FLOAT_MAX_2 in
    MY_PROVE_HYP le_th th_min, MY_PROVE_HYP le_th th_max;;
      

(*
let f1 = `float F (D1 _0) (D5 _0)`;;
let f2 = `float T (D2 _0) (D3 _0)`;;
float_min f1 f2;;
float_max f1 f2;;
float_min_max f1 f2;;
*)

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

(* FLOAT_MUL *)

let FLOAT_MUL = 
prove(`!s1 s2. min_exp <= e /\ num_exp n1 e1 * num_exp n2 e2 = num_exp n e
			  ==> float s1 n1 e1 * float s2 n2 e2 =
			      float (mod_plus s1 s2) n (e - min_exp)`,

   REPEAT STRIP_TAC THEN
     REWRITE_TAC[float] THEN
     ONCE_REWRITE_TAC[REAL_ARITH `(a * b / c) * (d * e / f) = (a * d) * (b * e) / c / f`] THEN

     SUBGOAL_THEN `(if s1 then -- &1 else &1) * (if s2 then -- &1 else &1) = if mod_plus s1 s2 then -- &1 else &1` MP_TAC THENL
     [
       REWRITE_TAC[mod_plus] THEN
         COND_CASES_TAC THEN COND_CASES_TAC THEN
         REWRITE_TAC[REAL_ARITH `-- &1 * -- &1 = &1`; REAL_MUL_LID; REAL_MUL_RID];
       ALL_TAC
     ] THEN

     DISCH_THEN (fun th -> REWRITE_TAC[th]) THEN
     REWRITE_TAC[real_div] THEN
     REWRITE_TAC[GSYM REAL_MUL_ASSOC] THEN
     REWRITE_TAC[REAL_EQ_MUL_LCANCEL] THEN
     DISJ2_TAC THEN

     MP_TAC (SPECL[`n:num`; `e:num`; `min_exp`] NUM_EXP_SUB_lemma) THEN
     ASM_REWRITE_TAC[] THEN
     DISCH_THEN (fun th -> REWRITE_TAC[th]) THEN
     REWRITE_TAC[REAL_MUL_ASSOC] THEN
     ASM_REWRITE_TAC[REAL_OF_NUM_MUL]);;


let FLOAT_MUL_FF = 
prove(`min_exp <= e /\ e - min_exp = r /\ num_exp n1 e1 * num_exp n2 e2 = num_exp n e ==>
			     float F n1 e1 * float F n2 e2 = float F n r`,

			 SIMP_TAC[FLOAT_MUL; mod_plus]);;
let FLOAT_MUL_FT = 
prove(`min_exp <= e /\ e - min_exp = r /\ num_exp n1 e1 * num_exp n2 e2 = num_exp n e ==>
			     float F n1 e1 * float T n2 e2 = float T n r`,

			 SIMP_TAC[FLOAT_MUL; mod_plus]);;
let FLOAT_MUL_TF = 
prove(`min_exp <= e /\ e - min_exp = r /\ num_exp n1 e1 * num_exp n2 e2 = num_exp n e ==>
			     float T n1 e1 * float F n2 e2 = float T n r`,

			 SIMP_TAC[FLOAT_MUL; mod_plus]);;
let FLOAT_MUL_TT = 
prove(`min_exp <= e /\ e - min_exp = r /\ num_exp n1 e1 * num_exp n2 e2 = num_exp n e ==>
			     float T n1 e1 * float T n2 e2 = float F n r`,

			 SIMP_TAC[FLOAT_MUL; mod_plus]);;


let FLOAT_MUL_FF_hi, FLOAT_MUL_FF_lo =
  let ff_hi = `min_exp <= e /\ e - min_exp = r /\ num_exp n1 e1 * num_exp n2 e2 <= num_exp n e
                      ==> float F n1 e1 * float F n2 e2 <= float F n r` in
  let ff_lo = `min_exp <= e /\ e - min_exp = r /\ num_exp n e <= num_exp n1 e1 * num_exp n2 e2
                      ==> float F n r <= float F n1 e1 * float F n2 e2` in
  let proof =
    REPEAT STRIP_TAC THEN
      POP_ASSUM MP_TAC THEN
      POP_ASSUM (fun th -> REWRITE_TAC[SYM th]) THEN
      REWRITE_TAC[GSYM REAL_OF_NUM_LE; GSYM REAL_OF_NUM_MUL] THEN
      DISCH_TAC THEN
      MAP_EVERY ABBREV_TAC [`z = &(num_exp n e)`; `x = &(num_exp n1 e1)`; `y = &(num_exp n2 e2)`] THEN
      ASM_REWRITE_TAC[float; REAL_MUL_LID] THEN
      REWRITE_TAC[REAL_ARITH `a / b * c / d = (a * c) / b / d`] THEN
      REWRITE_TAC[real_div] THEN
      REWRITE_TAC[REAL_MUL_ASSOC] THEN
      MATCH_MP_TAC REAL_LE_RMUL THEN
      REWRITE_TAC[REAL_LE_INV_EQ; REAL_POS] THEN

      MP_TAC (SPECL [`n:num`; `e:num`; `min_exp`] NUM_EXP_SUB_lemma) THEN
      ASM_REWRITE_TAC[] THEN
      DISCH_THEN (fun th -> REWRITE_TAC[th]) THEN
      MATCH_MP_TAC REAL_LE_RMUL THEN
      ASM_REWRITE_TAC[REAL_LE_INV_EQ; REAL_POS] in
    prove(ff_hi, proof), prove(ff_lo, proof);;



let FLOAT_MUL_TT_hi, FLOAT_MUL_TT_lo =
  let tt_hi = `min_exp <= e /\ e - min_exp = r /\ num_exp n1 e1 * num_exp n2 e2 <= num_exp n e
                     ==> float T n1 e1 * float T n2 e2 <= float F n r` in
  let tt_lo = `min_exp <= e /\ e - min_exp = r /\ num_exp n e <= num_exp n1 e1 * num_exp n2 e2
                     ==> float F n r <= float T n1 e1 * float T n2 e2` in
  let proof =
    REWRITE_TAC[FLOAT_NEG_T] THEN
      REWRITE_TAC[REAL_ARITH `--a * --b = a * b`] THEN
      REWRITE_TAC[FLOAT_MUL_FF_hi; FLOAT_MUL_FF_lo] in
    prove(tt_hi, proof), prove(tt_lo, proof);;



let FLOAT_MUL_FT_hi, FLOAT_MUL_FT_lo =
  let ft_hi = `min_exp <= e /\ e - min_exp = r /\ num_exp n e <= num_exp n1 e1 * num_exp n2 e2
                     ==> float F n1 e1 * float T n2 e2 <= float T n r` in
  let ft_lo = `min_exp <= e /\ e - min_exp = r /\ num_exp n1 e1 * num_exp n2 e2 <= num_exp n e
                     ==> float T n r <= float F n1 e1 * float T n2 e2` in
  let proof =
    REWRITE_TAC[FLOAT_NEG_T] THEN
      REWRITE_TAC[REAL_ARITH `a * --b <= --c <=> c <= a * b`] THEN
      REWRITE_TAC[REAL_ARITH `--c <= a * --b <=> a * b <= c`] THEN
      REWRITE_TAC[FLOAT_MUL_FF_hi; FLOAT_MUL_FF_lo] in
    prove(ft_hi, proof), prove(ft_lo, proof);;



let FLOAT_MUL_TF_hi, FLOAT_MUL_TF_lo =
  let ft_hi = `min_exp <= e /\ e - min_exp = r /\ num_exp n e <= num_exp n1 e1 * num_exp n2 e2
                     ==> float T n1 e1 * float F n2 e2 <= float T n r` in
  let ft_lo = `min_exp <= e /\ e - min_exp = r /\ num_exp n1 e1 * num_exp n2 e2 <= num_exp n e
                     ==> float T n r <= float T n1 e1 * float F n2 e2` in
  let proof =
    REWRITE_TAC[FLOAT_NEG_T] THEN
      REWRITE_TAC[REAL_ARITH `--a * b <= --c <=> c <= a * b`] THEN
      REWRITE_TAC[REAL_ARITH `--c <= --a * b <=> a * b <= c`] THEN
      REWRITE_TAC[FLOAT_MUL_FF_hi; FLOAT_MUL_FF_lo] in
    prove(ft_hi, proof), prove(ft_lo, proof);;


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

(* float_mul_lo, float_mul_hi *)

let transform = UNDISCH_ALL o NUMERALS_TO_NUM o REWRITE_RULE[min_exp_def; GSYM IMP_IMP];;
let FLOAT_MUL_FF_hi' = transform FLOAT_MUL_FF_hi and
    FLOAT_MUL_FF_lo' = transform FLOAT_MUL_FF_lo and
    FLOAT_MUL_TT_hi' = transform FLOAT_MUL_TT_hi and
    FLOAT_MUL_TT_lo' = transform FLOAT_MUL_TT_lo and
    FLOAT_MUL_FT_hi' = transform FLOAT_MUL_FT_hi and
    FLOAT_MUL_FT_lo' = transform FLOAT_MUL_FT_lo and
    FLOAT_MUL_TF_hi' = transform FLOAT_MUL_TF_hi and
    FLOAT_MUL_TF_lo' = transform FLOAT_MUL_TF_lo;;

let FLOAT_MUL_FF' = transform FLOAT_MUL_FF and
    FLOAT_MUL_TT' = transform FLOAT_MUL_TT and
    FLOAT_MUL_FT' = transform FLOAT_MUL_FT and
    FLOAT_MUL_TF' = transform FLOAT_MUL_TF;;


let float_mul_eq f1 f2 =
  let s1, n1, e1 = dest_float f1 and
      s2, n2, e2 = dest_float f2 in
  let flag = s1 = s2 in
  let num_exp1 = mk_num_exp n1 e1 and
      num_exp2 = mk_num_exp n2 e2 in

  let mul_th = num_exp_mul num_exp1 num_exp2 in
  let n_tm, e_tm = dest_num_exp (rand (concl mul_th)) in

  let sub_th, le_th = raw_sub_and_le_hash_conv e_tm min_exp_num_const in
    if (rand(concl le_th) <> e_tm) then
      failwith "float_mul_lo: underflow"
    else
      let r_tm = rand(concl sub_th) in
      let inst = INST[e_tm, e_var_num; r_tm, r_var_num; n_tm, n_var_num;
		      n1, n1_var_num; e1, e1_var_num; n2, n2_var_num; e2, e2_var_num] in
      let th0 = inst
	(if flag then
	   if s1 = "F" then FLOAT_MUL_FF' else FLOAT_MUL_TT'
	 else
	   if s1 = "F" then FLOAT_MUL_FT' else FLOAT_MUL_TF') in
	MY_PROVE_HYP sub_th (MY_PROVE_HYP mul_th (MY_PROVE_HYP le_th th0));;



let float_mul_lo pp f1 f2 =
  let s1, n1, e1 = dest_float f1 and
      s2, n2, e2 = dest_float f2 in
  let flag = s1 = s2 in
  let num_exp1 = mk_num_exp n1 e1 and
      num_exp2 = mk_num_exp n2 e2 in

  let mul_th, n_tm, e_tm = 
    if flag then
      let th = num_exp_op_lo pp num_exp_mul num_exp1 num_exp2 in
      let n_tm, e_tm = dest_num_exp (lhand (concl th)) in
	th, n_tm, e_tm
    else
      let th = num_exp_op_hi pp num_exp_mul num_exp1 num_exp2 in
      let n_tm, e_tm = dest_num_exp (rand (concl th)) in
	th, n_tm, e_tm in

  let sub_th, le_th = raw_sub_and_le_hash_conv e_tm min_exp_num_const in
    if (rand(concl le_th) <> e_tm) then
      failwith "float_mul_lo: underflow"
    else
      let r_tm = rand(concl sub_th) in
      let inst = INST[e_tm, e_var_num; r_tm, r_var_num; n_tm, n_var_num;
		      n1, n1_var_num; e1, e1_var_num; n2, n2_var_num; e2, e2_var_num] in
      let th0 = inst
	(if flag then
	   if s1 = "F" then FLOAT_MUL_FF_lo' else FLOAT_MUL_TT_lo'
	 else
	   if s1 = "F" then FLOAT_MUL_FT_lo' else FLOAT_MUL_TF_lo') in
	MY_PROVE_HYP sub_th (MY_PROVE_HYP mul_th (MY_PROVE_HYP le_th th0));;



let float_mul_hi pp f1 f2 =
  let s1, n1, e1 = dest_float f1 and
      s2, n2, e2 = dest_float f2 in
  let flag = s1 = s2 in
  let num_exp1 = mk_num_exp n1 e1 and
      num_exp2 = mk_num_exp n2 e2 in

  let mul_th, n_tm, e_tm = 
    if flag then
      let th = num_exp_op_hi pp num_exp_mul num_exp1 num_exp2 in
      let n_tm, e_tm = dest_num_exp (rand (concl th)) in
	th, n_tm, e_tm
    else
      let th = num_exp_op_lo pp num_exp_mul num_exp1 num_exp2 in
      let n_tm, e_tm = dest_num_exp (lhand (concl th)) in
	th, n_tm, e_tm in

  let sub_th, le_th = raw_sub_and_le_hash_conv e_tm min_exp_num_const in
    if (rand(concl le_th) <> e_tm) then
      failwith "float_mul_hi: underflow"
    else
      let r_tm = rand(concl sub_th) in
      let inst = INST[e_tm, e_var_num; r_tm, r_var_num; n_tm, n_var_num;
		      n1, n1_var_num; e1, e1_var_num; n2, n2_var_num; e2, e2_var_num] in
      let th0 = inst
	(if flag then
	   if s1 = "F" then FLOAT_MUL_FF_hi' else FLOAT_MUL_TT_hi'
	 else
	   if s1 = "F" then FLOAT_MUL_FT_hi' else FLOAT_MUL_TF_hi') in
	MY_PROVE_HYP sub_th (MY_PROVE_HYP mul_th (MY_PROVE_HYP le_th th0));;




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

(* FLOAT_DIV *)


let DIV_lemma = 
prove(`!x y. ~(y = 0) ==> &(x DIV y) <= &x / &y /\ &x / &y <= &(x DIV y + 1)`,

   REPEAT GEN_TAC THEN DISCH_TAC THEN
     MP_TAC (SPECL [`y:num`; `x:num`] FLOOR_DIV_DIV) THEN
     ASM_REWRITE_TAC[GSYM REAL_OF_NUM_ADD] THEN
     DISCH_THEN (fun th -> REWRITE_TAC[SYM th]) THEN
     SIMP_TAC[FLOOR; REAL_LT_IMP_LE]);;

let FLOAT_DIV_FF = 
prove(`e2 + k <= min_exp + e + e1 /\ ~(n2 = 0) /\
                           num_exp n1 k DIV num_exp n2 0 = num_exp n e
                             ==> float F n ((min_exp + e + e1) - (e2 + k)) <= float F n1 e1 / float F n2 e2`,

   MAP_EVERY ABBREV_TAC [`z = num_exp n e`; `x = num_exp n1 k`; `y = num_exp n2 0`] THEN
     REPEAT STRIP_TAC THEN
     REWRITE_TAC[float; REAL_MUL_LID] THEN
     REWRITE_TAC[real_div; REAL_INV_MUL; REAL_INV_INV] THEN
     REWRITE_TAC[REAL_ARITH `(a * b) * c * d = (b * d) * (a * c)`] THEN
     SUBGOAL_THEN `~(&(num_exp 1 min_exp) = &0)` ASSUME_TAC THENL
     [
       REWRITE_TAC[num_exp; REAL_OF_NUM_EQ; MULT_CLAUSES; EXP_EQ_0] THEN
         ARITH_TAC;
       ALL_TAC
     ] THEN

     ASM_SIMP_TAC[REAL_MUL_LINV; REAL_MUL_LID] THEN

     ASM_SIMP_TAC[NUM_EXP_SUB_lemma] THEN
     SUBGOAL_THEN `&(num_exp n1 e1) * inv(&(num_exp n2 e2)) = (&x / &y) * &(num_exp 1 e1) * inv(&(num_exp 1 (e2 + k)))` MP_TAC THENL
     [
       EXPAND_TAC "x" THEN EXPAND_TAC "y" THEN
         REWRITE_TAC[real_div] THEN
         REWRITE_TAC[num_exp; GSYM REAL_OF_NUM_MUL; GSYM REAL_OF_NUM_POW] THEN
         REWRITE_TAC[REAL_MUL_LID; REAL_INV_MUL; REAL_INV_1; real_pow; REAL_MUL_RID] THEN
         REWRITE_TAC[REAL_POW_ADD; REAL_INV_MUL] THEN
         REWRITE_TAC[REAL_ARITH `((a * b) * c) * d * e * f = (b * f) * (a * c * d * e)`] THEN

         SUBGOAL_THEN (mk_comb(`(~)`, mk_eq(mk_binop `pow` (mk_comb (`&`, base_const)) `k:num`, `&0`))) ASSUME_TAC THENL
         [
           REWRITE_TAC[REAL_POW_EQ_0] THEN
             REAL_ARITH_TAC;
           ALL_TAC
         ] THEN

         ASM_SIMP_TAC[REAL_MUL_RINV; REAL_MUL_LID] THEN
         REAL_ARITH_TAC;
       ALL_TAC
     ] THEN

     DISCH_THEN (fun th -> REWRITE_TAC[th]) THEN
     ONCE_REWRITE_TAC[REAL_ARITH `(a * b) * c = (a * c) * b`] THEN
     REWRITE_TAC[REAL_MUL_ASSOC] THEN
     MATCH_MP_TAC REAL_LE_RMUL THEN
     REWRITE_TAC[REAL_LE_INV_EQ; REAL_POS] THEN
     ONCE_REWRITE_TAC[NUM_EXP_SUM1] THEN
     REWRITE_TAC[NUM_EXP_SUM] THEN
     REWRITE_TAC[GSYM REAL_OF_NUM_MUL] THEN
     ASM_REWRITE_TAC[REAL_ARITH `(a * b * c) * d = (d * a) * b * c`] THEN
     ASM_SIMP_TAC[REAL_MUL_LINV; REAL_MUL_LID] THEN
     MATCH_MP_TAC REAL_LE_RMUL THEN
     REWRITE_TAC[REAL_LE_INV_EQ; REAL_POS] THEN
     MP_TAC (SPEC_ALL DIV_lemma) THEN
     ANTS_TAC THENL
     [
       EXPAND_TAC "y" THEN
         REWRITE_TAC[num_exp; MULT_EQ_0; DE_MORGAN_THM] THEN
         ASM_REWRITE_TAC[EXP] THEN
         ARITH_TAC;
       ALL_TAC
     ] THEN

     ASM_SIMP_TAC[]);;



let FLOAT_DIV_FF_lo = 
prove(`e2 + k = r1 /\ min_exp + e + e1 = r2 /\ r2 - r1 = r /\
                             r1 <= r2 /\ ~(n2 = 0) /\
                           num_exp n e <= num_exp n1 k DIV num_exp n2 0
                             ==> float F n r <= float F n1 e1 / float F n2 e2`,

   MAP_EVERY ABBREV_TAC [`z = num_exp n e`; `x = num_exp n1 k`; `y = num_exp n2 0`] THEN
     REPEAT STRIP_TAC THEN
     REPLICATE_TAC 3 (POP_ASSUM MP_TAC) THEN
     REPLICATE_TAC 3 (POP_ASSUM (fun th -> REWRITE_TAC[SYM th])) THEN
     REPEAT STRIP_TAC THEN
     REWRITE_TAC[float; REAL_MUL_LID] THEN
     REWRITE_TAC[real_div; REAL_INV_MUL; REAL_INV_INV] THEN
     REWRITE_TAC[REAL_ARITH `(a * b) * c * d = (b * d) * (a * c)`] THEN
     SUBGOAL_THEN `~(&(num_exp 1 min_exp) = &0)` ASSUME_TAC THENL
     [
       REWRITE_TAC[num_exp; REAL_OF_NUM_EQ; MULT_CLAUSES; EXP_EQ_0] THEN
         ARITH_TAC;
       ALL_TAC
     ] THEN

     ASM_SIMP_TAC[REAL_MUL_LINV; REAL_MUL_LID] THEN

     ASM_SIMP_TAC[NUM_EXP_SUB_lemma] THEN
     SUBGOAL_THEN `&(num_exp n1 e1) * inv(&(num_exp n2 e2)) = (&x / &y) * &(num_exp 1 e1) * inv(&(num_exp 1 (e2 + k)))` MP_TAC THENL
     [
       EXPAND_TAC "x" THEN EXPAND_TAC "y" THEN
         REWRITE_TAC[real_div] THEN
         REWRITE_TAC[num_exp; GSYM REAL_OF_NUM_MUL; GSYM REAL_OF_NUM_POW] THEN
         REWRITE_TAC[REAL_MUL_LID; REAL_INV_MUL; REAL_INV_1; real_pow; REAL_MUL_RID] THEN
         REWRITE_TAC[REAL_POW_ADD; REAL_INV_MUL] THEN
         REWRITE_TAC[REAL_ARITH `((a * b) * c) * d * e * f = (b * f) * (a * c * d * e)`] THEN
         SUBGOAL_THEN
         (mk_comb(`(~)`, mk_eq(mk_binop `pow` (mk_comb (`&`, base_const)) `k:num`, `&0`))) ASSUME_TAC THENL
         [
           REWRITE_TAC[REAL_POW_EQ_0] THEN
             REAL_ARITH_TAC;
           ALL_TAC
         ] THEN

         ASM_SIMP_TAC[REAL_MUL_RINV; REAL_MUL_LID] THEN
         REAL_ARITH_TAC;
       ALL_TAC
     ] THEN

     DISCH_THEN (fun th -> REWRITE_TAC[th]) THEN
     ONCE_REWRITE_TAC[REAL_ARITH `(a * b) * c = (a * c) * b`] THEN
     REWRITE_TAC[REAL_MUL_ASSOC] THEN
     MATCH_MP_TAC REAL_LE_RMUL THEN
     REWRITE_TAC[REAL_LE_INV_EQ; REAL_POS] THEN
     ONCE_REWRITE_TAC[NUM_EXP_SUM1] THEN
     REWRITE_TAC[NUM_EXP_SUM] THEN
     REWRITE_TAC[GSYM REAL_OF_NUM_MUL] THEN
     ASM_REWRITE_TAC[REAL_ARITH `(a * b * c) * d = (d * a) * b * c`] THEN
     ASM_SIMP_TAC[REAL_MUL_LINV; REAL_MUL_LID] THEN
     MATCH_MP_TAC REAL_LE_RMUL THEN
     REWRITE_TAC[REAL_LE_INV_EQ; REAL_POS] THEN
     MP_TAC (SPEC_ALL DIV_lemma) THEN
     ANTS_TAC THENL
     [
       EXPAND_TAC "y" THEN
         REWRITE_TAC[num_exp; MULT_EQ_0; DE_MORGAN_THM] THEN
         ASM_REWRITE_TAC[EXP] THEN
         ARITH_TAC;
       ALL_TAC
     ] THEN

     STRIP_TAC THEN
     MATCH_MP_TAC REAL_LE_TRANS THEN
     EXISTS_TAC `&(x DIV y)` THEN
     ASM_REWRITE_TAC[REAL_OF_NUM_LE]);;



let FLOAT_DIV_FF_hi = 
prove(`e2 + k = r1 /\ min_exp + e + e1 = r2 /\ r2 - r1 = r /\
                              r1 <= r2 /\ ~(n2 = 0) /\
                              num_exp n1 k DIV num_exp n2 0 < num_exp n e
                              ==> float F n1 e1 / float F n2 e2 <= float F n r`,

   MAP_EVERY ABBREV_TAC [`z = num_exp n e`; `x = num_exp n1 k`; `y = num_exp n2 0`] THEN
     REPEAT STRIP_TAC THEN
     REPLICATE_TAC 3 (POP_ASSUM MP_TAC) THEN
     REPLICATE_TAC 3 (POP_ASSUM (fun th -> REWRITE_TAC[SYM th])) THEN
     REPEAT STRIP_TAC THEN
     REWRITE_TAC[float; REAL_MUL_LID] THEN
     REWRITE_TAC[real_div; REAL_INV_MUL; REAL_INV_INV] THEN
     REWRITE_TAC[REAL_ARITH `(a * b) * c * d = (b * d) * (a * c)`] THEN
     SUBGOAL_THEN `~(&(num_exp 1 min_exp) = &0)` ASSUME_TAC THENL
     [
       REWRITE_TAC[num_exp; REAL_OF_NUM_EQ; MULT_CLAUSES; EXP_EQ_0] THEN
         ARITH_TAC;
       ALL_TAC
     ] THEN

     ASM_SIMP_TAC[REAL_MUL_LINV; REAL_MUL_LID] THEN

     ASM_SIMP_TAC[NUM_EXP_SUB_lemma] THEN
     SUBGOAL_THEN `&(num_exp n1 e1) * inv(&(num_exp n2 e2)) = (&x / &y) * &(num_exp 1 e1) * inv(&(num_exp 1 (e2 + k)))` MP_TAC THENL
     [
       EXPAND_TAC "x" THEN EXPAND_TAC "y" THEN
         REWRITE_TAC[real_div] THEN
         REWRITE_TAC[num_exp; GSYM REAL_OF_NUM_MUL; GSYM REAL_OF_NUM_POW] THEN
         REWRITE_TAC[REAL_MUL_LID; REAL_INV_MUL; REAL_INV_1; real_pow; REAL_MUL_RID] THEN
         REWRITE_TAC[REAL_POW_ADD; REAL_INV_MUL] THEN
         REWRITE_TAC[REAL_ARITH `((a * b) * c) * d * e * f = (b * f) * (a * c * d * e)`] THEN
         SUBGOAL_THEN
         (mk_comb(`(~)`, mk_eq(mk_binop `pow` (mk_comb (`&`, base_const)) `k:num`, `&0`))) ASSUME_TAC THENL
         [
           REWRITE_TAC[REAL_POW_EQ_0] THEN
             REAL_ARITH_TAC;
           ALL_TAC
         ] THEN

         ASM_SIMP_TAC[REAL_MUL_RINV; REAL_MUL_LID] THEN
         REAL_ARITH_TAC;
       ALL_TAC
     ] THEN

     DISCH_THEN (fun th -> REWRITE_TAC[th]) THEN
     ONCE_REWRITE_TAC[REAL_ARITH `(a * b) * c = (a * c) * b`] THEN
     REWRITE_TAC[REAL_MUL_ASSOC] THEN
     MATCH_MP_TAC REAL_LE_RMUL THEN
     REWRITE_TAC[REAL_LE_INV_EQ; REAL_POS] THEN
     ONCE_REWRITE_TAC[NUM_EXP_SUM1] THEN
     REWRITE_TAC[NUM_EXP_SUM] THEN
     REWRITE_TAC[GSYM REAL_OF_NUM_MUL] THEN
     ASM_REWRITE_TAC[REAL_ARITH `(a * b * c) * d = (d * a) * b * c`] THEN
     ASM_SIMP_TAC[REAL_MUL_LINV; REAL_MUL_LID] THEN
     MATCH_MP_TAC REAL_LE_RMUL THEN
     REWRITE_TAC[REAL_LE_INV_EQ; REAL_POS] THEN
     MP_TAC (SPEC_ALL DIV_lemma) THEN
     ANTS_TAC THENL
     [
       EXPAND_TAC "y" THEN
         REWRITE_TAC[num_exp; MULT_EQ_0; DE_MORGAN_THM] THEN
         ASM_REWRITE_TAC[EXP] THEN
         ARITH_TAC;
       ALL_TAC
     ] THEN

     STRIP_TAC THEN
     MATCH_MP_TAC REAL_LE_TRANS THEN
     EXISTS_TAC `&(x DIV y + 1)` THEN
     ASM_REWRITE_TAC[REAL_OF_NUM_LE] THEN
     UNDISCH_TAC `x DIV y < z` THEN
     ARITH_TAC);;


let FLOAT_DIV_TT_lo = 
prove(`e2 + k = r1 /\ min_exp + e + e1 = r2 /\ r2 - r1 = r /\
                              r1 <= r2 /\ ~(n2 = 0) /\
                              num_exp n e <= num_exp n1 k DIV num_exp n2 0
             ==> float F n r <= float T n1 e1 / float T n2 e2`,

   REWRITE_TAC[FLOAT_NEG_T] THEN
     REWRITE_TAC[real_div; REAL_INV_NEG; REAL_NEG_MUL2] THEN
     REWRITE_TAC[GSYM real_div] THEN
     REWRITE_TAC[FLOAT_DIV_FF_lo]);;
let FLOAT_DIV_TT_hi = 
prove(`e2 + k = r1 /\ min_exp + e + e1 = r2 /\
                              r2 - r1 = r /\ r1 <= r2 /\ ~(n2 = 0) /\
                              num_exp n1 k DIV num_exp n2 0 < num_exp n e
             ==> float T n1 e1 / float T n2 e2 <= float F n r`,

   REWRITE_TAC[FLOAT_NEG_T] THEN
     REWRITE_TAC[real_div; REAL_INV_NEG; REAL_NEG_MUL2] THEN
     REWRITE_TAC[GSYM real_div] THEN
     REWRITE_TAC[FLOAT_DIV_FF_hi]);;

let FLOAT_DIV_FT_lo = 
prove(`e2 + k = r1 /\ min_exp + e + e1 = r2 /\
                              r2 - r1 = r /\ r1 <= r2 /\ ~(n2 = 0) /\
                              num_exp n1 k DIV num_exp n2 0 < num_exp n e
             ==> float T n r <= float F n1 e1 / float T n2 e2`,

   REWRITE_TAC[FLOAT_NEG_T] THEN
     REWRITE_TAC[real_div; REAL_INV_NEG] THEN
     REWRITE_TAC[REAL_ARITH `--a <= b * --c <=> b * c <= a`] THEN
     REWRITE_TAC[GSYM real_div] THEN
     REWRITE_TAC[FLOAT_DIV_FF_hi]);;
let FLOAT_DIV_FT_hi = 
prove(`e2 + k = r1 /\ min_exp + e + e1 = r2 /\
                              r2 - r1 = r /\ r1 <= r2 /\ ~(n2 = 0) /\
                              num_exp n e <= num_exp n1 k DIV num_exp n2 0
             ==> float F n1 e1 / float T n2 e2 <= float T n r`,

   REWRITE_TAC[FLOAT_NEG_T] THEN
     REWRITE_TAC[real_div; REAL_INV_NEG] THEN
     REWRITE_TAC[REAL_ARITH `a * --b <= --c <=> c <= a * b`] THEN
     REWRITE_TAC[GSYM real_div] THEN
     REWRITE_TAC[FLOAT_DIV_FF_lo]);;

let FLOAT_DIV_TF_lo = 
prove(`e2 + k = r1 /\ min_exp + e + e1 = r2 /\
                              r2 - r1 = r /\ r1 <= r2 /\ ~(n2 = 0) /\
                              num_exp n1 k DIV num_exp n2 0 < num_exp n e
             ==> float T n r <= float T n1 e1 / float F n2 e2`,

   REWRITE_TAC[FLOAT_NEG_T] THEN
     REWRITE_TAC[real_div; REAL_INV_NEG] THEN
     REWRITE_TAC[REAL_ARITH `--a <= --b * c <=> b * c <= a`] THEN
     REWRITE_TAC[GSYM real_div] THEN
     REWRITE_TAC[FLOAT_DIV_FF_hi]);;
let FLOAT_DIV_TF_hi = 
prove(`e2 + k = r1 /\ min_exp + e + e1 = r2 /\
                              r2 - r1 = r /\ r1 <= r2 /\ ~(n2 = 0) /\
                              num_exp n e <= num_exp n1 k DIV num_exp n2 0
             ==> float T n1 e1 / float F n2 e2 <= float T n r`,

   REWRITE_TAC[FLOAT_NEG_T] THEN
     REWRITE_TAC[real_div; REAL_INV_NEG] THEN
     REWRITE_TAC[REAL_ARITH `--a * b <= --c <=> c <= a * b`] THEN
     REWRITE_TAC[GSYM real_div] THEN
     REWRITE_TAC[FLOAT_DIV_FF_lo]);;

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

(* float_div_lo, float_div_hi *)

let transform = UNDISCH_ALL o PURE_REWRITE_RULE[TAUT `~P <=> (P <=> F)`] o
  NUMERALS_TO_NUM o REWRITE_RULE[GSYM IMP_IMP; min_exp_def];;

let FLOAT_DIV_FF_hi' = transform FLOAT_DIV_FF_hi and
    FLOAT_DIV_FF_lo' = transform FLOAT_DIV_FF_lo and
    FLOAT_DIV_TT_hi' = transform FLOAT_DIV_TT_hi and
    FLOAT_DIV_TT_lo' = transform FLOAT_DIV_TT_lo and
    FLOAT_DIV_FT_hi' = transform FLOAT_DIV_FT_hi and
    FLOAT_DIV_FT_lo' = transform FLOAT_DIV_FT_lo and
    FLOAT_DIV_TF_hi' = transform FLOAT_DIV_TF_hi and
    FLOAT_DIV_TF_lo' = transform FLOAT_DIV_TF_lo;;



let float_div_lo pp f1 f2 =
  let s1, n1, e1 = dest_float f1 and
      s2, n2, e2 = dest_float f2 in
  let flag = s1 = s2 in

  let k_tm = rand (mk_small_numeral_array (2 * pp)) in
  let num_exp1 = mk_num_exp n1 k_tm and
      num_exp2 = mk_num_exp n2 zero_const in
  let div_th, n_tm, e_tm =
    if flag then
      let th = num_exp_op_lo pp num_exp_div num_exp1 num_exp2 in
      let n_tm, e_tm = dest_num_exp (lhand(concl th)) in
	th, n_tm, e_tm
    else
      let th = num_exp_op_hi_lt pp num_exp_div num_exp1 num_exp2 in
      let n_tm, e_tm = dest_num_exp (rand(concl th)) in
	th, n_tm, e_tm in

  let r1_th = raw_add_conv_hash (mk_binop plus_op_num e2 k_tm) in
  let r1_tm = rand(concl r1_th) in
  let e_plus_e1 = raw_add_conv_hash (mk_binop plus_op_num e_tm e1) in
  let ltm, rtm = dest_comb(concl e_plus_e1) in
  let r2_th' = raw_add_conv_hash (mk_binop plus_op_num min_exp_num_const rtm) in
  let r2_th = TRANS (AP_TERM (mk_comb (plus_op_num, min_exp_num_const)) e_plus_e1) r2_th' in
  let r2_tm = rand(concl r2_th) in
  let sub_th, le_th = raw_sub_and_le_hash_conv r2_tm r1_tm in
    if rand(concl le_th) <> r2_tm then
      failwith "float_div_lo: underflow"
    else
      let r_tm = rand(concl sub_th) in
      let n2_not_zero = raw_eq0_hash_conv n2 in
      let inst = INST[r1_tm, r1_var_num; r2_tm, r2_var_num;
		      n1, n1_var_num; e1, e1_var_num;
		      e_tm, e_var_num; k_tm, k_var_num;
		      n2, n2_var_num; e2, e2_var_num;
		      n_tm, n_var_num; r_tm, r_var_num] in
      let th0 = inst
	(if flag then
	   if s1 = "F" then FLOAT_DIV_FF_lo' else FLOAT_DIV_TT_lo'
	 else
	   if s1 = "F" then FLOAT_DIV_FT_lo' else FLOAT_DIV_TF_lo') in
      let th1 = MY_PROVE_HYP n2_not_zero (MY_PROVE_HYP div_th (MY_PROVE_HYP le_th th0)) in
	MY_PROVE_HYP sub_th (MY_PROVE_HYP r2_th (MY_PROVE_HYP r1_th th1));;





let float_div_hi pp f1 f2 =
  let s1, n1, e1 = dest_float f1 and
      s2, n2, e2 = dest_float f2 in
  let flag = s1 = s2 in

  let k_tm = rand (mk_small_numeral_array (2 * pp)) in
  let num_exp1 = mk_num_exp n1 k_tm and
      num_exp2 = mk_num_exp n2 zero_const in
  let div_th, n_tm, e_tm =
    if flag then
      let th = num_exp_op_hi_lt pp num_exp_div num_exp1 num_exp2 in
      let n_tm, e_tm = dest_num_exp (rand(concl th)) in
	th, n_tm, e_tm
    else
      let th = num_exp_op_lo pp num_exp_div num_exp1 num_exp2 in
      let n_tm, e_tm = dest_num_exp (lhand(concl th)) in
	th, n_tm, e_tm in

  let r1_th = raw_add_conv_hash (mk_binop plus_op_num e2 k_tm) in
  let r1_tm = rand(concl r1_th) in
  let e_plus_e1 = raw_add_conv_hash (mk_binop plus_op_num e_tm e1) in
  let ltm, rtm = dest_comb(concl e_plus_e1) in
  let r2_th' = raw_add_conv_hash (mk_binop plus_op_num min_exp_num_const rtm) in
  let r2_th = TRANS (AP_TERM (mk_comb (plus_op_num, min_exp_num_const)) e_plus_e1) r2_th' in
  let r2_tm = rand(concl r2_th) in
  let sub_th, le_th = raw_sub_and_le_hash_conv r2_tm r1_tm in
    if rand(concl le_th) <> r2_tm then
      failwith "float_div_hi: underflow"
    else
      let r_tm = rand(concl sub_th) in
      let n2_not_zero = raw_eq0_hash_conv n2 in
      let inst = INST[r1_tm, r1_var_num; r2_tm, r2_var_num;
		      n1, n1_var_num; e1, e1_var_num;
		      e_tm, e_var_num; k_tm, k_var_num;
		      n2, n2_var_num; e2, e2_var_num;
		      n_tm, n_var_num; r_tm, r_var_num] in
      let th0 = inst
	(if flag then
	   if s1 = "F" then FLOAT_DIV_FF_hi' else FLOAT_DIV_TT_hi'
	 else
	   if s1 = "F" then FLOAT_DIV_FT_hi' else FLOAT_DIV_TF_hi') in
      let th1 = MY_PROVE_HYP n2_not_zero (MY_PROVE_HYP div_th (MY_PROVE_HYP le_th th0)) in
	MY_PROVE_HYP sub_th (MY_PROVE_HYP r2_th (MY_PROVE_HYP r1_th th1));;





(*
float_div_lo 1 `float F (B1 _0) (B10 _0)` `float F (B3 _0) (B10 _0)`;;
float_div_hi 1 `float F (B1 _0) (B10 _0)` `float F (B3 _0) (B10 _0)`;;
float_div_hi 5 f1 f2;;
*)




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


(* FLOAT_ADD *)

let FLOAT_ADD_FF = 
prove(`num_exp n1 e1 + num_exp n2 e2 = num_exp n e
    ==> float F n1 e1 + float F n2 e2 = float F n e`,

   REPEAT STRIP_TAC THEN
     REWRITE_TAC[float; REAL_MUL_LID] THEN
     REWRITE_TAC[REAL_ARITH `a / b + c / b = (a + c) / b`] THEN
     ASM_REWRITE_TAC[REAL_OF_NUM_ADD]);;

let FLOAT_ADD_TT = 
prove(`num_exp n1 e1 + num_exp n2 e2 = num_exp n e
    ==> float T n1 e1 + float T n2 e2 = float T n e`,

   REWRITE_TAC[FLOAT_NEG_T; REAL_ARITH `--a + --b = --c <=> a + b = c`] THEN
     REWRITE_TAC[FLOAT_ADD_FF]);;


let FLOAT_ADD_FF_lo = 
prove(`num_exp n e <= num_exp n1 e1 + num_exp n2 e2
                              ==> float F n e <= float F n1 e1 + float F n2 e2`,

   REWRITE_TAC[GSYM REAL_OF_NUM_LE; GSYM REAL_OF_NUM_ADD] THEN
     REPEAT STRIP_TAC THEN
     MAP_EVERY ABBREV_TAC [`z = &(num_exp n e)`; `x = &(num_exp n1 e1)`; `y = &(num_exp n2 e2)`] THEN
     ASM_REWRITE_TAC[float; REAL_MUL_LID] THEN
     REWRITE_TAC[REAL_ARITH `a / b + c / b = (a + c) / b`] THEN
     REWRITE_TAC[real_div] THEN
     MATCH_MP_TAC REAL_LE_RMUL THEN
     ASM_REWRITE_TAC[REAL_LE_INV_EQ; REAL_POS]);;




let FLOAT_ADD_FF_hi = 
prove(`num_exp n1 e1 + num_exp n2 e2 <= num_exp n e
                              ==> float F n1 e1 + float F n2 e2 <= float F n e`,

   REWRITE_TAC[GSYM REAL_OF_NUM_LE; GSYM REAL_OF_NUM_ADD] THEN
     REPEAT STRIP_TAC THEN
     MAP_EVERY ABBREV_TAC [`z = &(num_exp n e)`; `x = &(num_exp n1 e1)`; `y = &(num_exp n2 e2)`] THEN
     ASM_REWRITE_TAC[float; REAL_MUL_LID] THEN
     REWRITE_TAC[REAL_ARITH `a / b + c / b = (a + c) / b`] THEN
     REWRITE_TAC[real_div] THEN
     MATCH_MP_TAC REAL_LE_RMUL THEN
     ASM_REWRITE_TAC[REAL_LE_INV_EQ; REAL_POS]);;


let FLOAT_ADD_TT_lo = 
prove(`num_exp n1 e1 + num_exp n2 e2 <= num_exp n e
			      ==> float T n e <= float T n1 e1 + float T n2 e2`,

   REWRITE_TAC[FLOAT_NEG_T; REAL_ARITH `--a <= --b + --c <=> b + c <= a`] THEN
     REWRITE_TAC[FLOAT_ADD_FF_hi]);;


let FLOAT_ADD_TT_hi = 
prove(`num_exp n e <= num_exp n1 e1 + num_exp n2 e2
			      ==> float T n1 e1 + float T n2 e2 <= float T n e`,

   REWRITE_TAC[FLOAT_NEG_T; REAL_ARITH `--b + --c <= --a <=> a <= b + c`] THEN
     REWRITE_TAC[FLOAT_ADD_FF_lo]);;


let FLOAT_ADD_FT_F_lo = 
prove(`num_exp n2 e2 <= num_exp n1 e1 ==>
				num_exp n e <= num_exp n1 e1 - num_exp n2 e2 
				==> float F n e <= float F n1 e1 + float T n2 e2`,

   MAP_EVERY ABBREV_TAC[`z = num_exp n e`; `x = num_exp n1 e1`; `y = num_exp n2 e2`] THEN
     ASM_REWRITE_TAC[FLOAT_NEG_T; float; REAL_MUL_LID] THEN
     DISCH_TAC THEN
     ASM_SIMP_TAC[GSYM REAL_OF_NUM_LE; GSYM REAL_OF_NUM_SUB] THEN
     REWRITE_TAC[num_exp; min_exp_def; MULT_CLAUSES; GSYM REAL_OF_NUM_POW] THEN
     REAL_ARITH_TAC);;

let FLOAT_ADD_FT_T_lo = 
prove(`num_exp n1 e1 <= num_exp n2 e2 ==>
				num_exp n2 e2 - num_exp n1 e1 <= num_exp n e
				==> float T n e <= float F n1 e1 + float T n2 e2`,

   MAP_EVERY ABBREV_TAC[`z = num_exp n e`; `x = num_exp n1 e1`; `y = num_exp n2 e2`] THEN
     ASM_REWRITE_TAC[FLOAT_NEG_T; float; REAL_MUL_LID] THEN
     DISCH_TAC THEN
     ASM_SIMP_TAC[GSYM REAL_OF_NUM_LE; GSYM REAL_OF_NUM_SUB] THEN
     REWRITE_TAC[num_exp; min_exp_def; MULT_CLAUSES; GSYM REAL_OF_NUM_POW] THEN
     REAL_ARITH_TAC);;

let FLOAT_ADD_FT_F_hi = 
prove(`num_exp n2 e2 <= num_exp n1 e1 ==>
				num_exp n1 e1 - num_exp n2 e2 <= num_exp n e
				==> float F n1 e1 + float T n2 e2 <= float F n e`,

   REWRITE_TAC[FLOAT_NEG_T; REAL_ARITH `a + --b <= c <=> --c <= b + --a`] THEN
     REWRITE_TAC[GSYM FLOAT_NEG_T; FLOAT_ADD_FT_T_lo]);;

let FLOAT_ADD_FT_T_hi = 
prove(`num_exp n1 e1 <= num_exp n2 e2 ==>
				num_exp n e <= num_exp n2 e2 - num_exp n1 e1
				==> float F n1 e1 + float T n2 e2 <= float T n e`,

   REWRITE_TAC[FLOAT_NEG_T; REAL_ARITH `a + --b <= --c <=> c <= b + --a`] THEN
     REWRITE_TAC[GSYM FLOAT_NEG_T; FLOAT_ADD_FT_F_lo]);;
     



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

(* float_add_lo, float_add_hi *)

let m_var_real = `m:real` and
    n_var_real = `n:real`;;

let REAL_ADD_COMM = CONJUNCT1 REAL_ADD_AC;;

let transform = UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o NUMERALS_TO_NUM;;
let FLOAT_ADD_FF_hi' = transform FLOAT_ADD_FF_hi and
    FLOAT_ADD_FF_lo' = transform FLOAT_ADD_FF_lo and
    FLOAT_ADD_TT_hi' = transform FLOAT_ADD_TT_hi and
    FLOAT_ADD_TT_lo' = transform FLOAT_ADD_TT_lo and
    FLOAT_ADD_FT_F_lo' = transform FLOAT_ADD_FT_F_lo and
    FLOAT_ADD_FT_T_lo' = transform FLOAT_ADD_FT_T_lo and
    FLOAT_ADD_FT_F_hi' = transform FLOAT_ADD_FT_F_hi and
    FLOAT_ADD_FT_T_hi' = transform FLOAT_ADD_FT_T_hi;;




let float_add_lo pp f1 f2 =
  let s1, n1, e1 = dest_float f1 in
  let s2, n2, e2 = dest_float f2 in
    if s1 = s2 then
      let num_exp1 = mk_num_exp n1 e1 in
      let num_exp2 = mk_num_exp n2 e2 in

	if s1 = "F" then
	  (* F + F *)
	  let add_th = num_exp_op_lo pp num_exp_add num_exp1 num_exp2 in
	  let n_tm, e_tm = dest_num_exp (lhand(concl add_th)) in
	  let th0 = INST[e_tm, e_var_num; n_tm, n_var_num; n1, n1_var_num;
			 e1, e1_var_num; n2, n2_var_num; e2, e2_var_num] FLOAT_ADD_FF_lo' in
	    MY_PROVE_HYP add_th th0
	else
	  (* T + T *)
	  let add_th = num_exp_op_hi pp num_exp_add num_exp1 num_exp2 in
	  let n_tm, e_tm = dest_num_exp (rand(concl add_th)) in
	  let th0 = INST[e_tm, e_var_num; n_tm, n_var_num; n1, n1_var_num;
			 e1, e1_var_num; n2, n2_var_num; e2, e2_var_num] FLOAT_ADD_TT_lo' in
	    MY_PROVE_HYP add_th th0
    else
      (* F + T or T + F *)
      let th0, n1, e1, n2, e2 =
	if s1 = "T" then
	  INST[f2, m_var_real; f1, n_var_real] REAL_ADD_COMM, n2, e2, n1, e1
	else
	  REFL(mk_binop plus_op_real f1 f2), n1, e1, n2, e2 in

      let num_exp1 = mk_num_exp n1 e1 in
      let num_exp2 = mk_num_exp n2 e2 in

      let sub_th, le_th = num_exp_sub num_exp1 num_exp2 in
      let sub_tm = rand(concl sub_th) in

	if rand(concl le_th) = num_exp1 then
	  let lo_th = num_exp_lo pp sub_tm in
	  let n_tm, e_tm = dest_num_exp (lhand(concl lo_th)) in
	  let lo_sub_th = EQ_MP (AP_TERM (rator(concl lo_th)) (SYM sub_th)) lo_th in

	  let th1 = INST[n1, n1_var_num; e1, e1_var_num; n2, n2_var_num; e2, e2_var_num;
			 n_tm, n_var_num; e_tm, e_var_num] FLOAT_ADD_FT_F_lo' in
	  let th2 = MY_PROVE_HYP lo_sub_th (MY_PROVE_HYP le_th th1) in
	    EQ_MP (AP_TERM (rator(concl th2)) th0) th2

	else 
	  let hi_th = num_exp_hi pp sub_tm in
	  let n_tm, e_tm = dest_num_exp(rand(concl hi_th)) in
	  let hi_sub_th = EQ_MP (SYM (AP_THM (AP_TERM le_op_num sub_th) (rand(concl hi_th)))) hi_th in

	  let th1 = INST[n1, n1_var_num; e1, e1_var_num; n2, n2_var_num; e2, e2_var_num;
			 n_tm, n_var_num; e_tm, e_var_num] FLOAT_ADD_FT_T_lo' in
	  let th2 = MY_PROVE_HYP hi_sub_th (MY_PROVE_HYP le_th th1) in
	    EQ_MP (AP_TERM (rator(concl th2)) th0) th2;;






let float_add_hi pp f1 f2 =
  let s1, n1, e1 = dest_float f1 in
  let s2, n2, e2 = dest_float f2 in
    if s1 = s2 then
      let num_exp1 = mk_num_exp n1 e1 in
      let num_exp2 = mk_num_exp n2 e2 in

	if s1 = "F" then
	  (* F + F *)
	  let add_th = num_exp_op_hi pp num_exp_add num_exp1 num_exp2 in
	  let n_tm, e_tm = dest_num_exp (rand(concl add_th)) in
	  let th0 = INST[e_tm, e_var_num; n_tm, n_var_num; n1, n1_var_num;
			 e1, e1_var_num; n2, n2_var_num; e2, e2_var_num] FLOAT_ADD_FF_hi' in
	    MY_PROVE_HYP add_th th0
	else
	  (* T + T *)
	  let add_th = num_exp_op_lo pp num_exp_add num_exp1 num_exp2 in
	  let n_tm, e_tm = dest_num_exp (lhand(concl add_th)) in
	  let th0 = INST[e_tm, e_var_num; n_tm, n_var_num; n1, n1_var_num;
			 e1, e1_var_num; n2, n2_var_num; e2, e2_var_num] FLOAT_ADD_TT_hi' in
	    MY_PROVE_HYP add_th th0
    else
      (* F + T or T + F *)
      let th0, n1, e1, n2, e2 =
	if s1 = "T" then
	  INST[f2, m_var_real; f1, n_var_real] REAL_ADD_COMM, n2, e2, n1, e1
	else
	  REFL(mk_binop plus_op_real f1 f2), n1, e1, n2, e2 in

      let num_exp1 = mk_num_exp n1 e1 in
      let num_exp2 = mk_num_exp n2 e2 in

      let sub_th, le_th = num_exp_sub num_exp1 num_exp2 in
      let sub_tm = rand(concl sub_th) in

	if rand(concl le_th) = num_exp1 then
	  let hi_th = num_exp_hi pp sub_tm in
	  let n_tm, e_tm = dest_num_exp (rand(concl hi_th)) in
	  let hi_sub_th = EQ_MP (SYM (AP_THM (AP_TERM le_op_num sub_th) (rand(concl hi_th)))) hi_th in

	  let th1 = INST[n1, n1_var_num; e1, e1_var_num; n2, n2_var_num; e2, e2_var_num;
			 n_tm, n_var_num; e_tm, e_var_num] FLOAT_ADD_FT_F_hi' in
	  let th2 = MY_PROVE_HYP hi_sub_th (MY_PROVE_HYP le_th th1) in
	    EQ_MP (AP_THM (AP_TERM le_op_real th0) (rand(concl th2))) th2

	else 
	  let lo_th = num_exp_lo pp sub_tm in
	  let n_tm, e_tm = dest_num_exp(lhand(concl lo_th)) in
	  let lo_sub_th = EQ_MP (AP_TERM (rator(concl lo_th)) (SYM sub_th)) lo_th in

	  let th1 = INST[n1, n1_var_num; e1, e1_var_num; n2, n2_var_num; e2, e2_var_num;
			 n_tm, n_var_num; e_tm, e_var_num] FLOAT_ADD_FT_T_hi' in
	  let th2 = MY_PROVE_HYP lo_sub_th (MY_PROVE_HYP le_th th1) in
	    EQ_MP (AP_THM (AP_TERM le_op_real th0) (rand(concl th2))) th2;;



(*
float_add_lo 1 f1 f2;;
float_add_hi 1 f1 f2;;
*)



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

(* float_sub_lo, float_sub_hi *)


let s1_var_bool = `s1:bool` and
    f1_var_real = `f1:real` and
    f2_var_real = `f2:real`;;


let FLOAT_SUB_F_EQ_ADD = (SYM o prove)(`f1 - float F n2 e2 = f1 + float T n2 e2`,
				       REWRITE_TAC[FLOAT_NEG_T] THEN REAL_ARITH_TAC);;

let FLOAT_SUB_T_EQ_ADD = (SYM o prove)(`f1 - float T n2 e2 = f1 + float F n2 e2`,
				       REWRITE_TAC[FLOAT_NEG_T] THEN REAL_ARITH_TAC);;



let float_sub_lo pp f1 f2 =
  let s2, n2, e2 = dest_float f2 in
  let th0 =
    INST[f1, f1_var_real; n2, n2_var_num; e2, e2_var_num] 
      (if s2 = "F" then FLOAT_SUB_F_EQ_ADD else FLOAT_SUB_T_EQ_ADD) in
  let ltm,f2_tm = dest_comb(lhand(concl th0)) in
  let f1_tm = rand ltm in
  let lo_th = float_add_lo pp f1_tm f2_tm in
    EQ_MP (AP_TERM (rator(concl lo_th)) th0) lo_th;;



let float_sub_hi pp f1 f2 =
  let s2, n2, e2 = dest_float f2 in
  let th0 =
    INST[f1, f1_var_real; n2, n2_var_num; e2, e2_var_num]
      (if s2 = "F" then FLOAT_SUB_F_EQ_ADD else FLOAT_SUB_T_EQ_ADD) in
  let ltm, f2_tm = dest_comb(lhand(concl th0)) in
  let f1_tm = rand ltm in
  let hi_th = float_add_hi pp f1_tm f2_tm in
    EQ_MP (AP_THM (AP_TERM le_op_real th0) (rand(concl hi_th))) hi_th;;




(*
let f1 = `float F (D5 (D2 _0)) (D0 (D1 _0))`;;
let f2 = `float F (D3 (D2 _0)) (D9 _0)`;;
float_sub_hi 1 f2 f1;;
*)


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

(* FLOAT_SQRT *)


(* float F m e = float F (B0 m) (PRE e) *)
let FLOAT_PRE_EXP = prove(mk_imp(`~(e = 0) /\ PRE e = e1`,
				 mk_eq(`float F m e`,
				       mk_comb(mk_comb(`float F`, mk_comb(b0_const, m_var_num)), `e1:num`))),
   STRIP_TAC THEN POP_ASSUM (fun th -> REWRITE_TAC[SYM th]) THEN
     REWRITE_TAC[float; REAL_MUL_LID; real_div; REAL_EQ_MUL_RCANCEL] THEN
     DISJ1_TAC THEN
     REWRITE_TAC[num_exp; b0_thm; REAL_OF_NUM_EQ] THEN
     SUBGOAL_THEN `e = SUC (PRE e)` MP_TAC THENL
     [
       POP_ASSUM MP_TAC THEN ARITH_TAC;
       ALL_TAC
     ] THEN
     DISCH_THEN (fun th -> CONV_TAC (LAND_CONV (ONCE_REWRITE_CONV[th]))) THEN
     REWRITE_TAC[EXP] THEN
     ARITH_TAC);;
	



let DIV2_EVEN_lemma = 
prove(`!n. EVEN n ==> 2 * (n DIV 2) = n`,

   GEN_TAC THEN
     REWRITE_TAC[EVEN_EXISTS] THEN
     STRIP_TAC THEN
     ASM_REWRITE_TAC[] THEN
     MATCH_MP_TAC (ARITH_RULE `x = y ==> 2 * x = 2 * y`) THEN
     MATCH_MP_TAC DIV_MULT THEN
     ARITH_TAC);;


let FLOAT_SQRT_EVEN_lo = 
prove(`f1 * f1 = f2 /\ f2 <= x /\
                                num_exp m (2 * p) = x /\ f1 = num_exp n1 e1
 				   /\ EVEN e /\ e DIV 2 = e2 /\
				   e1 + e2 + (min_exp DIV 2) = r /\ 
				   p <= r /\ r - p = r2 
				       ==> float F n1 r2 <= sqrt (float F m e)`,

   STRIP_TAC THEN
     UNDISCH_TAC `f2 <= x:num` THEN
     UNDISCH_TAC `num_exp m (2 * p) = x` THEN
     DISCH_THEN (fun th -> REWRITE_TAC[SYM th]) THEN
     UNDISCH_TAC `f1 * f1 = f2:num` THEN
     DISCH_THEN (fun th -> REWRITE_TAC[SYM th]) THEN
     UNDISCH_TAC `e1 + e2 + min_exp DIV 2 = r` THEN
     UNDISCH_TAC `e DIV 2 = e2` THEN
     DISCH_THEN (fun th -> REWRITE_TAC[SYM th]) THEN
     REPEAT (POP_ASSUM MP_TAC) THEN
     REWRITE_TAC[num_exp; float; REAL_MUL_LID; GSYM REAL_OF_NUM_MUL] THEN
     REPEAT STRIP_TAC THEN
     UNDISCH_TAC `r - p = r2:num` THEN
     DISCH_THEN (fun th -> REWRITE_TAC[SYM th]) THEN

     MATCH_MP_TAC REAL_LE_RSQRT THEN
     REWRITE_TAC[GSYM REAL_OF_NUM_POW; REAL_POW_DIV] THEN
     REWRITE_TAC[REAL_POW_2; real_div; REAL_INV_MUL] THEN
     REWRITE_TAC[REAL_MUL_ASSOC] THEN
     MATCH_MP_TAC REAL_LE_RMUL THEN
     CONJ_TAC THENL
     [
       REWRITE_TAC[REAL_ARITH `(((a * b) * a) * b) * c = (a * a) * (b * b) * c:real`] THEN
	 REWRITE_TAC[GSYM REAL_POW_ADD] THEN
	 REWRITE_TAC[ARITH_RULE `r - p + r - p = 2 * r - 2 * p`] THEN
	 MP_TAC (SPECL[mk_comb(amp_op_real, base_const); `2 * r`; `2 * p`] REAL_DIV_POW2) THEN
	 ANTS_TAC THENL [ REAL_ARITH_TAC; ALL_TAC ] THEN
	 ASM_SIMP_TAC[ARITH_RULE `p <= r ==> 2 * p <= 2 * r`] THEN
	 DISCH_THEN (fun th -> REWRITE_TAC[SYM th; real_div]) THEN

	 SUBGOAL_THEN `2 * r = (e1 + e1) + min_exp + e` (fun th -> REWRITE_TAC[th]) THENL
	 [
	   EXPAND_TAC "r" THEN
	     REWRITE_TAC[ARITH_RULE `2 * (e1 + b + c) = (e1 + e1) + 2 * c + 2 * b`] THEN
	     MATCH_MP_TAC (ARITH_RULE `b1 = b2 /\ c1 = c2 ==> a + b1 + c1 = a + b2 + c2:num`) THEN
	     SUBGOAL_THEN `EVEN min_exp` ASSUME_TAC THENL
	     [
	       REWRITE_TAC[min_exp_def] THEN ARITH_TAC;
	       ALL_TAC
	     ] THEN
	     ASM_SIMP_TAC[DIV2_EVEN_lemma];
	   ALL_TAC
	 ] THEN

	 REWRITE_TAC[REAL_POW_ADD] THEN
	 REWRITE_TAC[REAL_ARITH `(n * n) * (((e * e) * x * y) * z) * u = (n * e) * (n * e) * (x * u) * z * y:real`] THEN
	 SUBGOAL_THEN `~(&(num_exp 1 min_exp) = &0)` MP_TAC THENL
	 [
	   REWRITE_TAC[REAL_OF_NUM_EQ; NUM_EXP_EQ_0] THEN ARITH_TAC;
	   ALL_TAC
	 ] THEN
	 REWRITE_TAC[num_exp; REAL_MUL_LID; GSYM REAL_OF_NUM_POW; GSYM REAL_OF_NUM_MUL] THEN
	 DISCH_THEN (fun th -> SIMP_TAC[th; REAL_MUL_RINV; REAL_MUL_LID]) THEN
	 FIRST_X_ASSUM (MP_TAC o check(fun th -> (fst o dest_var o lhand o concl) th = "f1")) THEN
	 REWRITE_TAC[GSYM REAL_OF_NUM_EQ; GSYM REAL_OF_NUM_MUL; GSYM REAL_OF_NUM_POW] THEN
	 DISCH_THEN (fun th -> REWRITE_TAC[SYM th]) THEN

	 REWRITE_TAC[REAL_MUL_ASSOC] THEN
	 MATCH_MP_TAC REAL_LE_RMUL THEN
	 CONJ_TAC THENL
	 [
	   SUBGOAL_THEN `!x y z. &0 < x /\ y <= z * x ==> y * inv x <= z` MP_TAC THENL
	     [
	       REPEAT STRIP_TAC THEN
		 MP_TAC (SPECL [`y * inv x`; `z:real`; `x:real`] REAL_LE_RMUL_EQ) THEN
		 ASM_REWRITE_TAC[] THEN
		 DISCH_THEN (fun th -> REWRITE_TAC[SYM th; GSYM REAL_MUL_ASSOC]) THEN
		 ASM_SIMP_TAC[REAL_ARITH `&0 < x ==> ~(x = &0)`; REAL_MUL_LINV; REAL_MUL_RID];
	       ALL_TAC
	     ] THEN

	     DISCH_THEN (MP_TAC o SPECL[`&(num_exp 1 (2 * p))`; `&(f1 * f1)`; `&m`]) THEN
	     REWRITE_TAC[num_exp; GSYM REAL_OF_NUM_MUL; GSYM REAL_OF_NUM_POW; REAL_MUL_LID] THEN
	     DISCH_THEN MATCH_MP_TAC THEN
	     ASM_REWRITE_TAC[REAL_OF_NUM_MUL; REAL_OF_NUM_POW; REAL_OF_NUM_LE] THEN
	     REWRITE_TAC[REAL_OF_NUM_LT; EXP_LT_0] THEN
	     ARITH_TAC;
	   ALL_TAC
	 ] THEN

	 MATCH_MP_TAC REAL_POW_LE THEN
	 ARITH_TAC;
       ALL_TAC
     ] THEN

     MATCH_MP_TAC REAL_LE_INV THEN
     MATCH_MP_TAC REAL_POW_LE THEN
     ARITH_TAC);;

let FLOAT_SQRT_EVEN_hi = 
prove(`f1 * f1 = f2 /\ x <= f2 /\
                                num_exp m (2 * p) = x /\ f1 = num_exp n1 e1
 				   /\ EVEN e /\ e DIV 2 = e2 /\
				   e1 + e2 + (min_exp DIV 2) = r /\ 
				   p <= r /\ r - p = r2 
				       ==> sqrt (float F m e) <= float F n1 r2`,

   STRIP_TAC THEN
     UNDISCH_TAC `x <= f2:num` THEN
     UNDISCH_TAC `num_exp m (2 * p) = x` THEN
     DISCH_THEN (fun th -> REWRITE_TAC[SYM th]) THEN
     UNDISCH_TAC `f1 * f1 = f2:num` THEN
     DISCH_THEN (fun th -> REWRITE_TAC[SYM th]) THEN
     UNDISCH_TAC `e1 + e2 + min_exp DIV 2 = r` THEN
     UNDISCH_TAC `e DIV 2 = e2` THEN
     DISCH_THEN (fun th -> REWRITE_TAC[SYM th]) THEN
     REPEAT (POP_ASSUM MP_TAC) THEN
     REWRITE_TAC[num_exp; float; REAL_MUL_LID; GSYM REAL_OF_NUM_MUL] THEN
     REPEAT STRIP_TAC THEN
     UNDISCH_TAC `r - p = r2:num` THEN
     DISCH_THEN (fun th -> REWRITE_TAC[SYM th]) THEN

     MATCH_MP_TAC REAL_LE_LSQRT THEN
     REPEAT CONJ_TAC THENL
     [
       REWRITE_TAC[REAL_OF_NUM_MUL; real_div] THEN
	 MATCH_MP_TAC REAL_LE_MUL THEN
	 REWRITE_TAC[REAL_POS] THEN
	 MATCH_MP_TAC REAL_LE_INV THEN
	 REWRITE_TAC[REAL_POS];
       REWRITE_TAC[REAL_OF_NUM_MUL; real_div] THEN
	 MATCH_MP_TAC REAL_LE_MUL THEN
	 REWRITE_TAC[REAL_POS] THEN
	 MATCH_MP_TAC REAL_LE_INV THEN
	 REWRITE_TAC[REAL_POS];
       ALL_TAC
     ] THEN

     REWRITE_TAC[GSYM REAL_OF_NUM_POW; REAL_POW_DIV] THEN
     REWRITE_TAC[REAL_POW_2; real_div; REAL_INV_MUL] THEN
     REWRITE_TAC[REAL_MUL_ASSOC] THEN
     MATCH_MP_TAC REAL_LE_RMUL THEN
     CONJ_TAC THENL
     [
       REWRITE_TAC[REAL_ARITH `(((a * b) * a) * b) * c = (a * a) * (b * b) * c:real`] THEN
	 REWRITE_TAC[GSYM REAL_POW_ADD] THEN
	 REWRITE_TAC[ARITH_RULE `r - p + r - p = 2 * r - 2 * p`] THEN
	 MP_TAC (SPECL[mk_comb(amp_op_real, base_const); `2 * r`; `2 * p`] REAL_DIV_POW2) THEN
	 ANTS_TAC THENL [ REAL_ARITH_TAC; ALL_TAC ] THEN
	 ASM_SIMP_TAC[ARITH_RULE `p <= r ==> 2 * p <= 2 * r`] THEN
	 DISCH_THEN (fun th -> REWRITE_TAC[SYM th; real_div]) THEN

	 SUBGOAL_THEN `2 * r = (e1 + e1) + min_exp + e` (fun th -> REWRITE_TAC[th]) THENL
	 [
	   EXPAND_TAC "r" THEN
	     REWRITE_TAC[ARITH_RULE `2 * (e1 + b + c) = (e1 + e1) + 2 * c + 2 * b`] THEN
	     MATCH_MP_TAC (ARITH_RULE `b1 = b2 /\ c1 = c2 ==> a + b1 + c1 = a + b2 + c2:num`) THEN
	     SUBGOAL_THEN `EVEN min_exp` ASSUME_TAC THENL
	     [
	       REWRITE_TAC[min_exp_def] THEN ARITH_TAC;
	       ALL_TAC
	     ] THEN
	     ASM_SIMP_TAC[DIV2_EVEN_lemma];
	   ALL_TAC
	 ] THEN

	 REWRITE_TAC[REAL_POW_ADD] THEN
	 REWRITE_TAC[REAL_ARITH `(n * n) * (((e * e) * x * y) * z) * u = (n * e) * (n * e) * (x * u) * z * y:real`] THEN
	 SUBGOAL_THEN `~(&(num_exp 1 min_exp) = &0)` MP_TAC THENL
	 [
	   REWRITE_TAC[REAL_OF_NUM_EQ; NUM_EXP_EQ_0] THEN ARITH_TAC;
	   ALL_TAC
	 ] THEN
	 REWRITE_TAC[num_exp; REAL_MUL_LID; GSYM REAL_OF_NUM_POW; GSYM REAL_OF_NUM_MUL] THEN
	 DISCH_THEN (fun th -> SIMP_TAC[th; REAL_MUL_RINV; REAL_MUL_LID]) THEN
	 FIRST_X_ASSUM (MP_TAC o check(fun th -> (fst o dest_var o lhand o concl) th = "f1")) THEN
	 REWRITE_TAC[GSYM REAL_OF_NUM_EQ; GSYM REAL_OF_NUM_MUL; GSYM REAL_OF_NUM_POW] THEN
	 DISCH_THEN (fun th -> REWRITE_TAC[SYM th]) THEN

	 REWRITE_TAC[REAL_MUL_ASSOC] THEN
	 MATCH_MP_TAC REAL_LE_RMUL THEN
	 CONJ_TAC THENL
	 [
	   SUBGOAL_THEN `!x y z. &0 < x /\ z * x <= y ==> z <= y * inv x` MP_TAC THENL
	     [
	       REPEAT STRIP_TAC THEN
		 MP_TAC (SPECL [`z:real`; `y * inv x`; `x:real`] REAL_LE_RMUL_EQ) THEN
		 ASM_REWRITE_TAC[] THEN
		 DISCH_THEN (fun th -> REWRITE_TAC[SYM th; GSYM REAL_MUL_ASSOC]) THEN
		 ASM_SIMP_TAC[REAL_ARITH `&0 < x ==> ~(x = &0)`; REAL_MUL_LINV; REAL_MUL_RID];
	       ALL_TAC
	     ] THEN

	     DISCH_THEN (MP_TAC o SPECL[`&(num_exp 1 (2 * p))`; `&(f1 * f1)`; `&m`]) THEN
	     REWRITE_TAC[num_exp; GSYM REAL_OF_NUM_MUL; GSYM REAL_OF_NUM_POW; REAL_MUL_LID] THEN
	     DISCH_THEN MATCH_MP_TAC THEN
	     ASM_REWRITE_TAC[REAL_OF_NUM_MUL; REAL_OF_NUM_POW; REAL_OF_NUM_LE] THEN
	     REWRITE_TAC[REAL_OF_NUM_LT; EXP_LT_0] THEN
	     ARITH_TAC;
	   ALL_TAC
	 ] THEN

	 MATCH_MP_TAC REAL_POW_LE THEN
	 ARITH_TAC;
       ALL_TAC
     ] THEN

     MATCH_MP_TAC REAL_LE_INV THEN
     MATCH_MP_TAC REAL_POW_LE THEN
     ARITH_TAC);;




let transform = UNDISCH_ALL o
  PURE_ONCE_REWRITE_RULE[TAUT `EVEN e <=> (EVEN e <=> T)`] o
  NUMERALS_TO_NUM o
  CONV_RULE (DEPTH_CONV NUM_DIV_CONV) o
  REWRITE_RULE[GSYM IMP_IMP; min_exp_def];;


let FLOAT_SQRT_EVEN_lo' = transform FLOAT_SQRT_EVEN_lo and
    FLOAT_SQRT_EVEN_hi' = transform FLOAT_SQRT_EVEN_hi and
    FLOAT_PRE_EXP' = (UNDISCH_ALL o 
			PURE_ONCE_REWRITE_RULE[TAUT `~(e = _0) <=> ((e = _0) <=> F)`] o
			REWRITE_RULE[GSYM IMP_IMP; NUMERAL]) FLOAT_PRE_EXP;;


let even_const = `EVEN` and
    pre_const = `PRE` and
    two_num = rand(mk_small_numeral_array 2) and
    min_exp_div2 = rand(mk_small_numeral_array (min_exp / 2)) and
    f2_var_num = `f2:num` and
    f1_var_num = `f1:num` and
    p_var_num = `p:num`;;


(* Returns the list of digits of the given Big_int n in the base b *)
let rec get_big_int_digits b n =
  let bb = big_int_of_int b in
  if le_big_int n zero_big_int then []
  else
    let q, r = quomod_big_int n bb in
      r :: get_big_int_digits b q;;


(* [1;2;3] -> 123 (base = 10) *)
let rec big_int_from_list b list =
  let rec proc acc list =
    match list with
	[] -> acc
      | h::t -> proc (add_big_int h (mult_int_big_int b acc)) t in
    proc zero_big_int list;;


(* Returns n first elements of the list *)
let rec take n list =
  match list with
      x :: xs -> if n > 0 then x :: take (n - 1) xs else []
    |	[] -> [];;



(* Returns an integer number that contains at most pp "significant" digits
   in the given base b *)
let big_int_round_lo base pp n =
  let digits = rev (get_big_int_digits base n) in
  let n_digits = length digits in
    if n_digits <= pp then
      n
    else
      let m = big_int_from_list base (take pp digits) in
	mult_big_int (power_int_positive_int base (n_digits - pp)) m;;



let big_int_round_hi base pp n =
  let digits = rev (get_big_int_digits base n) in
  let n_digits = length digits in
    if n_digits <= pp then n
    else
      let l1, l2 = chop_list pp digits in
	if forall (eq_big_int zero_big_int) l2 then n
	else
	  let m = succ_big_int (big_int_from_list base l1) in
	    mult_big_int (power_int_positive_int base (n_digits - pp)) m;;
	    
	

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

let rec float_sqrt_lo pp tm =
  let s, m_tm, e_tm = dest_float tm in
  let p_tm = rand (mk_small_numeral_array pp) in
    if s <> "F" then
      failwith "float_sqrt_lo: negative argument"
    else
      let even_th = raw_even_hash_conv (mk_comb (even_const, e_tm)) in
	if (fst o dest_const o rand o concl) even_th <> "T" then
	  (* ODD e *)
	  let pre_e = raw_pre_hash_conv (mk_comb (pre_const, e_tm)) in
	  let e_neq_0 = raw_eq0_hash_conv e_tm in
	  let e1_tm = rand (concl pre_e) in
	  let th0 = INST[e1_tm, e1_var_num; e_tm, e_var_num; m_tm, m_var_num] FLOAT_PRE_EXP' in
	  let th1 = MY_PROVE_HYP pre_e (MY_PROVE_HYP e_neq_0 th0) in
	  let th2 = float_sqrt_lo pp (rand(concl th1)) in
	  let ltm, rtm = dest_comb (concl th2) in
	    EQ_MP (SYM (AP_TERM ltm (AP_TERM (rator rtm) th1))) th2
	else
	  (* EVEN e *)
	  let p2_tm = mk_binop mul_op_num two_num p_tm in
	  let p2_th = raw_mul_conv_hash p2_tm in
	  let f1_1 = AP_TERM (mk_comb(num_exp_const, m_tm)) p2_th in
	  let f1_2 = TRANS f1_1 (denormalize (rand (concl f1_1))) in

	  let x_tm = rand(concl f1_2) in
	  let x = raw_dest_hash x_tm in
	  let f1' = Big_int.sqrt_big_int (big_int_of_num x) in
	  let f1 = num_of_big_int (big_int_round_lo base pp f1') in
	  let f1_tm = rand(mk_numeral_array f1) in
	  let f1_num_exp = to_num_exp f1_tm in

	  let n1_tm, e1_tm = dest_num_exp (rand (concl f1_num_exp)) in
	  let f1f1_eq_f2 = raw_mul_conv_hash (mk_binop mul_op_num f1_tm f1_tm) in
	  let f2_tm = rand(concl f1f1_eq_f2) in
	  let f2_le_x = EQT_ELIM (raw_le_hash_conv (mk_binop le_op_num f2_tm x_tm)) in

	  let e_div2_eq_e2 = raw_div_hash_conv (mk_binop div_op_num e_tm two_num) in
	  let e2_tm = rand(concl e_div2_eq_e2) in
	  let r_th1 = raw_add_conv_hash (mk_binop plus_op_num e2_tm min_exp_div2) in
	  let r_th2 = AP_TERM (mk_comb(plus_op_num, e1_tm)) r_th1 in
	  let r_th = TRANS r_th2 (raw_add_conv_hash (rand (concl r_th2))) in

	  let r_tm = rand(concl r_th) in
	  let r_sub_p, p_le_r = raw_sub_and_le_hash_conv p_tm r_tm in
	  let r2_tm = rand(concl r_sub_p) in
	    if (rand(concl p_le_r) <> r_tm) then
	      failwith "float_sqrt_lo: underflow"
	    else
	      let th0 = INST[f2_tm, f2_var_num; x_tm, x_var_num; p_tm, p_var_num; r_tm, r_var_num;
			     f1_tm, f1_var_num; n1_tm, n1_var_num; e1_tm, e1_var_num; e2_tm, e2_var_num;
			     e_tm, e_var_num; m_tm, m_var_num; r2_tm, r2_var_num] 
            		FLOAT_SQRT_EVEN_lo' in
		MY_PROVE_HYP f1_2 (
		  MY_PROVE_HYP e_div2_eq_e2 (
		    MY_PROVE_HYP r_sub_p (
		      MY_PROVE_HYP r_th (
			MY_PROVE_HYP f1f1_eq_f2 (
			  MY_PROVE_HYP f1_num_exp (
			    MY_PROVE_HYP even_th (
			      MY_PROVE_HYP f2_le_x (
				MY_PROVE_HYP p_le_r th0
			      ))))))));;





let rec float_sqrt_hi pp tm =
  let s, m_tm, e_tm = dest_float tm in
  let p_tm = rand (mk_small_numeral_array pp) in
    if s <> "F" then
      failwith "float_sqrt_lo: negative argument"
    else
      let even_th = raw_even_hash_conv (mk_comb (even_const, e_tm)) in
	if (fst o dest_const o rand o concl) even_th <> "T" then
	  (* ODD e *)
	  let pre_e = raw_pre_hash_conv (mk_comb (pre_const, e_tm)) in
	  let e_neq_0 = raw_eq0_hash_conv e_tm in
	  let e1_tm = rand (concl pre_e) in
	  let th0 = INST[e1_tm, e1_var_num; e_tm, e_var_num; m_tm, m_var_num] FLOAT_PRE_EXP' in
	  let th1 = MY_PROVE_HYP pre_e (MY_PROVE_HYP e_neq_0 th0) in
	  let th2 = float_sqrt_hi pp (rand(concl th1)) in
	  let ltm, rtm = dest_comb (concl th2) in
	  let ltm2, rtm2 = dest_comb ltm in
	  let th3 = AP_THM (AP_TERM ltm2 (AP_TERM (rator rtm2) th1)) rtm in
	    EQ_MP (SYM th3) th2
	else
	  (* EVEN e *)
	  let p2_tm = mk_binop mul_op_num two_num p_tm in
	  let p2_th = raw_mul_conv_hash p2_tm in
	  let f1_1 = AP_TERM (mk_comb(num_exp_const, m_tm)) p2_th in
	  let f1_2 = TRANS f1_1 (denormalize (rand (concl f1_1))) in

	  let x_tm = rand(concl f1_2) in
	  let x = raw_dest_hash x_tm in
	  let x' = big_int_of_num x in
	  let f1' = sqrt_big_int x' in
	  let f1 = (num_of_big_int o big_int_round_hi base pp)
	    (if eq_big_int (mult_big_int f1' f1') x' then f1' else succ_big_int f1') in
	      
	  let f1_tm = rand(mk_numeral_array f1) in
	  let f1_num_exp = to_num_exp f1_tm in

	  let n1_tm, e1_tm = dest_num_exp (rand (concl f1_num_exp)) in
	  let f1f1_eq_f2 = raw_mul_conv_hash (mk_binop mul_op_num f1_tm f1_tm) in
	  let f2_tm = rand(concl f1f1_eq_f2) in
	  let x_le_f2 = EQT_ELIM (raw_le_hash_conv (mk_binop le_op_num x_tm f2_tm)) in

	  let e_div2_eq_e2 = raw_div_hash_conv (mk_binop div_op_num e_tm two_num) in
	  let e2_tm = rand(concl e_div2_eq_e2) in
	  let r_th1 = raw_add_conv_hash (mk_binop plus_op_num e2_tm min_exp_div2) in
	  let r_th2 = AP_TERM (mk_comb(plus_op_num, e1_tm)) r_th1 in
	  let r_th = TRANS r_th2 (raw_add_conv_hash (rand (concl r_th2))) in

	  let r_tm = rand(concl r_th) in
	  let r_sub_p, p_le_r = raw_sub_and_le_hash_conv p_tm r_tm in
	  let r2_tm = rand(concl r_sub_p) in
	    if (rand(concl p_le_r) <> r_tm) then
	      failwith "float_sqrt_lo: underflow"
	    else
	      let th0 = INST[f2_tm, f2_var_num; x_tm, x_var_num; p_tm, p_var_num; r_tm, r_var_num;
			     f1_tm, f1_var_num; n1_tm, n1_var_num; e1_tm, e1_var_num; e2_tm, e2_var_num;
			     e_tm, e_var_num; m_tm, m_var_num; r2_tm, r2_var_num] 
            		FLOAT_SQRT_EVEN_hi' in
		MY_PROVE_HYP f1_2 (
		  MY_PROVE_HYP e_div2_eq_e2 (
		    MY_PROVE_HYP r_sub_p (
		      MY_PROVE_HYP r_th (
			MY_PROVE_HYP f1f1_eq_f2 (
			  MY_PROVE_HYP f1_num_exp (
			    MY_PROVE_HYP even_th (
			      MY_PROVE_HYP x_le_f2 (
				MY_PROVE_HYP p_le_r th0
			      ))))))));;



(*
float_sqrt_lo 3 `float F (D2 _0) (D0 (D2 _0))`;;
float_sqrt_hi 3 `float F (D2 _0) (D0 (D2 _0))`;;

let tm = `float F (D2 _0) (D0 (D2 _0))`;;
(* 10: 1.092 *)
test 1000 (float_sqrt_lo 5) tm;;
(* 10: 1.128 *)
test 1000 (float_sqrt_hi 5) tm;;
*)



end;; (* Float_ops module *)

(************************************)
(* Cached floating point operations *)
(************************************)

(* Counters for collecting stats *)
let lt0_c = ref 0 and
    gt0_c = ref 0 and
    lt_c = ref 0 and
    le0_c = ref 0 and
    ge0_c = ref 0 and
    le_c = ref 0 and
    min_c = ref 0 and
    max_c = ref 0 and
    min_max_c = ref 0 and
    mul_lo_c = ref 0 and
    mul_hi_c = ref 0 and
    div_lo_c = ref 0 and
    div_hi_c = ref 0 and
    add_lo_c = ref 0 and
    add_hi_c = ref 0 and
    sub_lo_c = ref 0 and
    sub_hi_c = ref 0 and
	sqrt_lo_c = ref 0 and
	sqrt_hi_c = ref 0;;

(*
let mul_lo_list = ref [] and
    mul_hi_list = ref [] and
    div_lo_list = ref [] and
    div_hi_list = ref [] and
    add_lo_list = ref [] and
    add_hi_list = ref [] and
    sub_lo_list = ref [] and
    sub_hi_list = ref [];;
*)
	
(* Hash tables *)
let cache_size = if Arith_options.cached then 10000 else 1;;
let max_cache_size = cache_size * 2;;

let my_add h key v =
  if Hashtbl.length h >= max_cache_size then
    let _ = Hashtbl.clear h in
      print_string "Clearing a float hash table"
  else 
    ();
  Hashtbl.add h key v;;
  
  
let mul_table = Hashtbl.create cache_size and
    div_table = Hashtbl.create cache_size and
    add_table = Hashtbl.create cache_size and
    sub_table = Hashtbl.create cache_size and
    sqrt_table = Hashtbl.create cache_size and
    le_table = Hashtbl.create cache_size and
    max_table = Hashtbl.create cache_size;;

let reset_cache () =
  Hashtbl.clear mul_table;
  Hashtbl.clear div_table;
  Hashtbl.clear add_table;
  Hashtbl.clear sub_table;
  Hashtbl.clear sqrt_table;
  Hashtbl.clear le_table;
  Hashtbl.clear max_table;;

let reset_stat () =
  lt0_c := 0;
  gt0_c := 0;
  lt_c := 0;
  le0_c := 0;
  ge0_c := 0;
  le_c := 0;
  min_c := 0;
  max_c := 0;
  min_max_c := 0;
  mul_lo_c := 0;
  mul_hi_c := 0;
  div_lo_c := 0;
  div_hi_c := 0;
  add_lo_c := 0;
  add_hi_c := 0;
  sub_lo_c := 0;
  sub_hi_c := 0;
  sqrt_lo_c := 0;
  sqrt_hi_c := 0;;

let print_stat () =
  let len = Hashtbl.length in
  let cmp_str1 = sprintf "lt0 = %d\ngt0 = %d\nlt = %d\n" !lt0_c !gt0_c !lt_c and
      cmp_str2 = sprintf "le0 = %d\nge0 = %d\n" !le0_c !ge0_c and
      cmp_str3 = sprintf "min = %d\nmin_max = %d\n" !min_c !min_max_c and
      le_str = sprintf "le = %d (le_hash = %d)\n" !le_c (len le_table) and
      max_str = sprintf "max = %d (max_hash = %d)\n" !max_c (len max_table) and
      mul_str = sprintf "mul_lo = %d, mul_hi = %d (mul_hash = %d)\n" !mul_lo_c !mul_hi_c (len mul_table) and
      div_str = sprintf "div_lo = %d, div_hi = %d (div_hash = %d)\n" !div_lo_c !div_hi_c (len div_table) and
      add_str = sprintf "add_lo = %d, add_hi = %d (add_hash = %d)\n" !add_lo_c !add_hi_c (len add_table) and
      sub_str = sprintf "sub_lo = %d, sub_hi = %d (sub_hash = %d)\n" !sub_lo_c !sub_hi_c (len sub_table) and
      sqrt_str = sprintf "sqrt_lo = %d, sqrt_hi = %d (sqrt_hash = %d)\n" !sqrt_lo_c !sqrt_hi_c (len sqrt_table) in
    print_string (cmp_str1 ^ cmp_str2 ^ cmp_str3 ^ 
		    le_str ^ max_str ^ mul_str ^ div_str ^ add_str ^ sub_str ^ sqrt_str);;


(* lt0 *)
let float_lt0 =
  if float_cached then
    fun tm ->
      let _ = lt0_c := !lt0_c + 1 in
	Float_ops.float_lt0 tm
  else
    Float_ops.float_lt0;;

(* gt0 *)
let float_gt0 =
  if float_cached then
    fun tm ->
      let _ = gt0_c := !gt0_c + 1 in
	Float_ops.float_gt0 tm
  else
    Float_ops.float_gt0;;

(* lt *)
let float_lt =
  if float_cached then
    fun tm1 tm2 ->
      let _ = lt_c := !lt_c + 1 in
	Float_ops.float_lt tm1 tm2
  else
    Float_ops.float_lt;;

(* le0 *)
let float_le0 =
  if float_cached then
    fun tm ->
      let _ = le0_c := !le0_c + 1 in
	Float_ops.float_le0 tm
  else
    Float_ops.float_le0;;

(* ge0 *)
let float_ge0 =
  if float_cached then
    fun tm ->
      let _ = ge0_c := !ge0_c + 1 in
	Float_ops.float_ge0 tm
  else
    Float_ops.float_ge0;;

(* min *)
let float_min =
  if float_cached then
    fun tm1 tm2 ->
      let _ = min_c := !min_c + 1 in
	Float_ops.float_min tm1 tm2
  else
    Float_ops.float_min;;

(* min_max *)
let float_min_max =
  if float_cached then
    fun tm1 tm2 ->
      let _ = min_max_c := !min_max_c + 1 in
	Float_ops.float_min_max tm1 tm2
  else
    Float_ops.float_min_max;;


(***************)
let float_hash tm =
  let s, n_tm, e_tm = dest_float tm in
    s ^ (Arith_cache.num_tm_hash n_tm) ^ "e" ^ (Arith_cache.num_tm_hash e_tm);;

let float_op_hash pp tm1 tm2 =
  string_of_int pp ^ float_hash tm1 ^ "x" ^ float_hash tm2;;
  
let float_op_hash1 pp tm =
	string_of_int pp ^ float_hash tm;;

  
(* le *)
let float_le =
  if float_cached then
    fun tm1 tm2 ->
      let _ = le_c := !le_c + 1 in
      let hash = float_op_hash 0 tm1 tm2 in
	try
	  Hashtbl.find le_table hash
	with Not_found ->
	  let result = Float_ops.float_le tm1 tm2 in
	  let _ = my_add le_table hash result in
	    result
  else
    Float_ops.float_le;;

(* max *)
let float_max =
  if float_cached then
    fun tm1 tm2 ->
      let _ = max_c := !max_c + 1 in
      let hash = float_op_hash 0 tm1 tm2 in
	try
	  Hashtbl.find max_table hash
	with Not_found ->
	  let result = Float_ops.float_max tm1 tm2 in
	  let _ = my_add max_table hash result in
	    result
  else
    Float_ops.float_max;;

(* mul_eq *)
let float_mul_eq = Float_ops.float_mul_eq;;

(* mul_lo *)
let float_mul_lo =
  if float_cached then
    fun pp tm1 tm2 ->
      let _ = mul_lo_c := !mul_lo_c + 1 in
(*      let _ = mul_lo_list := (pp,tm1,tm2) :: !mul_lo_list in *)
      let hash = "lo" ^ float_op_hash pp tm1 tm2 in
	try
	  Hashtbl.find mul_table hash
	with Not_found ->
	  let result = Float_ops.float_mul_lo pp tm1 tm2 in
	  let _ = my_add mul_table hash result in
	    result
  else
    Float_ops.float_mul_lo;;

(* mul_hi *)
let float_mul_hi =
  if float_cached then
    fun pp tm1 tm2 ->
      let _ = mul_hi_c := !mul_hi_c + 1 in
(*      let _ = mul_hi_list := (pp,tm1,tm2) :: !mul_hi_list in *)
      let hash = "hi" ^ float_op_hash pp tm1 tm2 in
	try
	  Hashtbl.find mul_table hash
	with Not_found ->
	  let result = Float_ops.float_mul_hi pp tm1 tm2 in
	  let _ = my_add mul_table hash result in
	    result
  else
    Float_ops.float_mul_hi;;

(* div_lo *)
let float_div_lo =
  if float_cached then
    fun pp tm1 tm2 ->
      let _ = div_lo_c := !div_lo_c + 1 in
(*      let _ = div_lo_list := (pp,tm1,tm2) :: !div_lo_list in *)
      let hash = "lo" ^ float_op_hash pp tm1 tm2 in
	try
	  Hashtbl.find div_table hash
	with Not_found ->
	  let result = Float_ops.float_div_lo pp tm1 tm2 in
	  let _ = my_add div_table hash result in
	    result
  else
    Float_ops.float_div_lo;;

(* div_hi *)
let float_div_hi =
  if float_cached then
    fun pp tm1 tm2 ->
      let _ = div_hi_c := !div_hi_c + 1 in
(*      let _ = div_hi_list := (pp,tm1,tm2) :: !div_hi_list in *)
      let hash = "hi" ^ float_op_hash pp tm1 tm2 in
	try
	  Hashtbl.find div_table hash
	with Not_found ->
	  let result = Float_ops.float_div_hi pp tm1 tm2 in
	  let _ = my_add div_table hash result in
	    result
  else
    Float_ops.float_div_hi;;

(* add_lo *)
let float_add_lo =
  if float_cached then
    fun pp tm1 tm2 ->
      let _ = add_lo_c := !add_lo_c + 1 in
(*      let _ = add_lo_list := (pp,tm1,tm2) :: !add_lo_list in *)
      let hash = "lo" ^ float_op_hash pp tm1 tm2 in
	try
	  Hashtbl.find add_table hash
	with Not_found ->
	  let result = Float_ops.float_add_lo pp tm1 tm2 in
	  let _ = my_add add_table hash result in
	    result
  else
    Float_ops.float_add_lo;;

(* add_hi *)
let float_add_hi =
  if float_cached then
    fun pp tm1 tm2 ->
      let _ = add_hi_c := !add_hi_c + 1 in
(*      let _ = add_hi_list := (pp,tm1,tm2) :: !add_hi_list in *)
      let hash = "hi" ^ float_op_hash pp tm1 tm2 in
	try
	  Hashtbl.find add_table hash
	with Not_found ->
	  let result = Float_ops.float_add_hi pp tm1 tm2 in
	  let _ = my_add add_table hash result in
	    result
  else
    Float_ops.float_add_hi;;

(* sub_lo *)
let float_sub_lo =
  if float_cached then
    fun pp tm1 tm2 ->
      let _ = sub_lo_c := !sub_lo_c + 1 in
(*      let _ = sub_lo_list := (pp,tm1,tm2) :: !sub_lo_list in *)
      let hash = "lo" ^ float_op_hash pp tm1 tm2 in
	try
	  Hashtbl.find sub_table hash
	with Not_found ->
	  let result = Float_ops.float_sub_lo pp tm1 tm2 in
	  let _ = my_add sub_table hash result in
	    result
  else
    Float_ops.float_sub_lo;;

(* sub_hi *)
let float_sub_hi =
  if float_cached then
    fun pp tm1 tm2 ->
      let _ = sub_hi_c := !sub_hi_c + 1 in
(*      let _ = sub_hi_list := (pp,tm1,tm2) :: !sub_hi_list in *)
      let hash = "hi" ^ float_op_hash pp tm1 tm2 in
	try
	  Hashtbl.find sub_table hash
	with Not_found ->
	  let result = Float_ops.float_sub_hi pp tm1 tm2 in
	  let _ = my_add sub_table hash result in
	    result
  else
    Float_ops.float_sub_hi;;

	
(* sqrt_lo *)
let float_sqrt_lo =
  if float_cached then
    fun pp tm ->
      let _ = sqrt_lo_c := !sqrt_lo_c + 1 in
      let hash = "lo" ^ float_op_hash1 pp tm in
	try
	  Hashtbl.find sqrt_table hash
	with Not_found ->
	  let result = Float_ops.float_sqrt_lo pp tm in
	  let _ = my_add sqrt_table hash result in
	    result
  else
    Float_ops.float_sqrt_lo;;

(* sqrt_hi *)
let float_sqrt_hi =
  if float_cached then
    fun pp tm ->
      let _ = sqrt_hi_c := !sqrt_hi_c + 1 in
      let hash = "hi" ^ float_op_hash1 pp tm in
	try
	  Hashtbl.find sqrt_table hash
	with Not_found ->
	  let result = Float_ops.float_sqrt_hi pp tm in
	  let _ = my_add sqrt_table hash result in
	    result
  else
    Float_ops.float_sqrt_hi;;
    


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

(* float intervals *)

let FLOAT_OF_NUM' = (SPEC_ALL o REWRITE_RULE[min_exp_def]) FLOAT_OF_NUM;;

let FLOAT_INTERVAL_OF_NUM = (NUMERALS_TO_NUM o REWRITE_RULE[min_exp_def] o prove)(`interval_arith (&n) (float F n min_exp, float F n min_exp)`,
   REWRITE_TAC[FLOAT_OF_NUM; CONST_INTERVAL]);;

let FLOAT_F_bound' = (UNDISCH_ALL o SPEC_ALL) FLOAT_F_bound;;

let FLOAT_T_bound' = (UNDISCH_ALL o SPEC_ALL) FLOAT_T_bound;;






(* interval_arith x (float s1 n1 e1, float s2 n2 e2) -> x, float s1 n1 e1, float s2 n2 e2 *)
let dest_float_interval tm =
  let ltm, rtm = dest_comb tm in
  let f1, f2 = dest_pair rtm in
    rand ltm, f1, f2;;




let mk_float_interval_small_num n =
  let n_tm0 = mk_small_numeral n in
  let n_th = NUMERAL_TO_NUM_CONV n_tm0 in
  let n_tm = rand(rand(concl n_th)) in
  let n_th1 = TRANS n_th (INST[n_tm, n_var_num] NUM_REMOVE) in
  let th1 = AP_TERM amp_op_real n_th1 in
  let int_th = INST[n_tm, n_var_num] FLOAT_INTERVAL_OF_NUM in
  let rtm = rand(concl int_th) in
    EQ_MP (SYM (AP_THM (AP_TERM interval_const th1) rtm)) int_th;;


let mk_float_interval_num n =
  let n_tm0 = mk_numeral n in
  let n_th = NUMERAL_TO_NUM_CONV n_tm0 in
  let n_tm = rand(rand(concl n_th)) in
  let n_th1 = TRANS n_th (INST[n_tm, n_var_num] NUM_REMOVE) in
  let th1 = AP_TERM amp_op_real n_th1 in
  let int_th = INST[n_tm, n_var_num] FLOAT_INTERVAL_OF_NUM in
  let rtm = rand(concl int_th) in
    EQ_MP (SYM (AP_THM (AP_TERM interval_const th1) rtm)) int_th;;






(* Returns the lower bound for the given float *)
let float_lo p tm =
  let s, n_tm, e_tm = dest_float tm in
    if s = "F" then
      let num_exp_tm = mk_num_exp n_tm e_tm in
      let th0 = num_exp_lo p num_exp_tm in
      let ltm, e1_tm = dest_comb(lhand(concl th0)) in
      let n1_tm = rand ltm in
      let th1 = INST[n1_tm, n1_var_num; e1_tm, e1_var_num; n_tm, n2_var_num; e_tm, e2_var_num] FLOAT_F_bound' in
        MY_PROVE_HYP th0 th1
    else
      let num_exp_tm = mk_num_exp n_tm e_tm in
      let th0 = num_exp_hi p num_exp_tm in
      let ltm, e1_tm = dest_comb(rand(concl th0)) in
      let n1_tm = rand ltm in
      let th1 = INST[n_tm, n1_var_num; e_tm, e1_var_num; n1_tm, n2_var_num; e1_tm, e2_var_num] FLOAT_T_bound' in
        MY_PROVE_HYP th0 th1;;



(* Returns the upper bound for the given float *)
let float_hi p tm =
  let s, n_tm, e_tm = dest_float tm in
    if s = "F" then
      let num_exp_tm = mk_num_exp n_tm e_tm in
      let th0 = num_exp_hi p num_exp_tm in
      let ltm, e2_tm = dest_comb(rand(concl th0)) in
      let n2_tm = rand ltm in
      let th1 = INST[n_tm, n1_var_num; e_tm, e1_var_num; n2_tm, n2_var_num; e2_tm, e2_var_num] FLOAT_F_bound' in
        MY_PROVE_HYP th0 th1
    else
      let num_exp_tm = mk_num_exp n_tm e_tm in
      let th0 = num_exp_lo p num_exp_tm in
      let ltm, e1_tm = dest_comb(lhand(concl th0)) in
      let n1_tm = rand ltm in
      let th1 = INST[n1_tm, n1_var_num; e1_tm, e1_var_num; n_tm, n2_var_num; e_tm, e2_var_num] FLOAT_T_bound' in
        MY_PROVE_HYP th0 th1;;






(* Approximates the given interval with p-digit float numbers *)
let float_interval_round p th =
  let x_tm, f1, f2 = dest_float_interval (concl th) in
  let lo_th = float_lo p f1 in
  let hi_th = float_hi p f2 in
  let lo_tm = lhand(concl lo_th) in
  let hi_tm = rand(concl hi_th) in
  let th0 = INST[x_tm, x_var_real; f1, lo_var_real; f2, hi_var_real; lo_tm, a_var_real; hi_tm, b_var_real] APPROX_INTERVAL' in
    MY_PROVE_HYP lo_th (MY_PROVE_HYP hi_th (MY_PROVE_HYP th th0));;




(*
let th = mk_float_interval_num 124235353;;
float_interval_round 5 th;;
(* 4: 0.240 *)
test 1000 (float_interval_round 5) th;;
*)


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

(* float_interval_lt *)

let FLOAT_INTERVAL_LT = 
prove(`interval_arith x (lo1, hi1) /\ interval_arith y (lo2, hi2) /\ hi1 < lo2
				==> x < y`,

			      REWRITE_TAC[interval_arith] THEN
				REAL_ARITH_TAC);;


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

(* float_interval_neg *)


let FLOAT_INTERVAL_NEG = 
prove(`!s1 s2. interval_arith x (float s1 n1 e1, float s2 n2 e2)
                                 ==> interval_arith (--x) (float (~s2) n2 e2, float (~s1) n1 e1)`,

   REPEAT GEN_TAC THEN
     DISCH_THEN (fun th -> MP_TAC (MATCH_MP INTERVAL_NEG th)) THEN
     SIMP_TAC[FLOAT_NEG]);;

let FLOAT_INTERVAL_NEG_FF = (UNDISCH_ALL o REWRITE_RULE[] o SPECL[`F`; `F`]) FLOAT_INTERVAL_NEG;;
let FLOAT_INTERVAL_NEG_FT = (UNDISCH_ALL o REWRITE_RULE[] o SPECL[`F`; `T`]) FLOAT_INTERVAL_NEG;;
let FLOAT_INTERVAL_NEG_TF = (UNDISCH_ALL o REWRITE_RULE[] o SPECL[`T`; `F`]) FLOAT_INTERVAL_NEG;;
let FLOAT_INTERVAL_NEG_TT = (UNDISCH_ALL o REWRITE_RULE[] o SPECL[`T`; `T`]) FLOAT_INTERVAL_NEG;;





(* |- interval x (float s1 n1 e1, float s2 n2 e2) ->
   |- interval (--x) (float ~s2 n2 e2, float ~s1 n1 e1 *)
let float_interval_neg th =
  let x_tm, f1, f2 = dest_float_interval (concl th) in
  let s1, n1_tm, e1_tm = dest_float f1 in
  let s2, n2_tm, e2_tm = dest_float f2 in
  let inst = INST[x_tm, x_var_real; n1_tm, n1_var_num; e1_tm, e1_var_num;
                  n2_tm, n2_var_num; e2_tm, e2_var_num] in
  let th0 =
    if s1 = "F" then
      if s2 = "F" then
        inst FLOAT_INTERVAL_NEG_FF
      else
        inst FLOAT_INTERVAL_NEG_FT
    else
      if s2 = "F" then
        inst FLOAT_INTERVAL_NEG_TF
      else
        inst FLOAT_INTERVAL_NEG_TT in
    MY_PROVE_HYP th th0;;



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

(* float_interval_mul *)


let f1_var_real = `f1:real` and
    f2_var_real = `f2:real` and
    f1_1_var = `f1_1:real` and
    f1_2_var = `f1_2:real` and
    f2_1_var = `f2_1:real` and
    f2_2_var = `f2_2:real`;;



(* FF_FF *)
let FLOAT_INTERVAL_MUL_FF_FF = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o prove)(
  `interval_arith x (float F n1 e1, float F n2 e2) /\
    interval_arith y (float F m1 r1, float F m2 r2) /\
    f1 <= float F n1 e1 * float F m1 r1 /\
    float F n2 e2 * float F m2 r2 <= f2
    ==> interval_arith (x * y) (f1, f2)`,
   MAP_EVERY ABBREV_TAC [`a = float F n1 e1`; `b = float F n2 e2`; `c = float F m1 r1`; `d = float F m2 r2`] THEN
     SUBGOAL_THEN `&0 <= a /\ &0 <= b /\ &0 <= c /\ &0 <= d` MP_TAC THENL
     [
       MAP_EVERY EXPAND_TAC ["a"; "b"; "c"; "d"] THEN
         REWRITE_TAC[FLOAT_F_POS];
       ALL_TAC
     ] THEN
     REPEAT (POP_ASSUM (fun th -> ALL_TAC)) THEN
     REWRITE_TAC[interval_arith] THEN
     REPEAT STRIP_TAC THENL
     [
       MATCH_MP_TAC REAL_LE_TRANS THEN
         EXISTS_TAC `a * c:real` THEN
         ASM_REWRITE_TAC[] THEN
         MATCH_MP_TAC REAL_LE_MUL2 THEN
         ASM_REWRITE_TAC[];
       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
         CONJ_TAC THEN MATCH_MP_TAC REAL_LE_TRANS THENL
         [
           EXISTS_TAC `a:real` THEN ASM_REWRITE_TAC[];
           EXISTS_TAC `c:real` THEN ASM_REWRITE_TAC[]
         ]
     ]);;



(* TT_TT *)
let FLOAT_INTERVAL_MUL_TT_TT = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o prove)(
  `interval_arith x (float T n1 e1, float T n2 e2) /\
    interval_arith y (float T m1 r1, float T m2 r2) /\
    f1 <= float T n2 e2 * float T m2 r2 /\
    float T n1 e1 * float T m1 r1 <= f2
    ==> interval_arith (x * y) (f1, f2)`,
  REWRITE_TAC[FLOAT_NEG_T] THEN
    REWRITE_TAC[interval_arith] THEN
    MAP_EVERY ABBREV_TAC [`a = float F n1 e1`; `b = float F n2 e2`; `c = float F m1 r1`; `d = float F m2 r2`] THEN
    SUBGOAL_THEN `&0 <= a /\ &0 <= b /\ &0 <= c /\ &0 <= d` MP_TAC THENL
    [
      MAP_EVERY EXPAND_TAC ["a"; "b"; "c"; "d"] THEN
	REWRITE_TAC[FLOAT_F_POS];
      ALL_TAC
    ] THEN
    REWRITE_TAC[REAL_NEG_MUL2] THEN
    REPEAT STRIP_TAC THENL
    [
      MATCH_MP_TAC REAL_LE_TRANS THEN
	EXISTS_TAC `b * d:real` THEN
	ASM_REWRITE_TAC[] THEN
	ONCE_REWRITE_TAC[REAL_ARITH `a <= x * y <=> a <= --x * --y`] THEN
	MATCH_MP_TAC REAL_LE_MUL2 THEN
	ONCE_REWRITE_TAC[REAL_ARITH `b <= --x <=> x <= --b`] THEN
	ASM_REWRITE_TAC[];
      MATCH_MP_TAC REAL_LE_TRANS THEN
	EXISTS_TAC `a * c:real` THEN
	ASM_REWRITE_TAC[] THEN
	ONCE_REWRITE_TAC[REAL_ARITH `x * y <= a <=> --x * --y <= a`] THEN
	MATCH_MP_TAC REAL_LE_MUL2 THEN
	ONCE_REWRITE_TAC[REAL_ARITH `--x <= c <=> --c <= x`] THEN
	ASM_REWRITE_TAC[] THEN
	ASSUME_TAC (REAL_ARITH `!b x. &0 <= b /\ x <= --b ==> &0 <= --x`) THEN
	CONJ_TAC THEN FIRST_X_ASSUM MATCH_MP_TAC THENL
	[
	  EXISTS_TAC `b:real` THEN ASM_REWRITE_TAC[];
	  EXISTS_TAC `d:real` THEN ASM_REWRITE_TAC[]
	]
    ]);;



(* FF_TT *)
let FLOAT_INTERVAL_MUL_FF_TT = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o prove)(
  `interval_arith x (float F n1 e1, float F n2 e2) /\
    interval_arith y (float T m1 r1, float T m2 r2) /\
    f1 <= float F n2 e2 * float T m1 r1 /\
    float F n1 e1 * float T m2 r2 <= f2
    ==> interval_arith (x * y) (f1, f2)`,
  REWRITE_TAC[FLOAT_NEG_T] THEN
    REWRITE_TAC[interval_arith] THEN
    MAP_EVERY ABBREV_TAC [`a = float F n1 e1`; `b = float F n2 e2`; `c = float F m1 r1`; `d = float F m2 r2`] THEN
    SUBGOAL_THEN `&0 <= a /\ &0 <= b /\ &0 <= c /\ &0 <= d` MP_TAC THENL
    [
      MAP_EVERY EXPAND_TAC ["a"; "b"; "c"; "d"] THEN
	REWRITE_TAC[FLOAT_F_POS];
      ALL_TAC
    ] THEN

    REPEAT STRIP_TAC THENL
    [
      MATCH_MP_TAC REAL_LE_TRANS THEN
	EXISTS_TAC `b * --c` THEN
	ASM_REWRITE_TAC[] THEN
	ONCE_REWRITE_TAC[REAL_ARITH `b * --c <= x * y <=> x * --y <= b * c`] THEN
	MATCH_MP_TAC REAL_LE_MUL2 THEN
	ONCE_REWRITE_TAC[REAL_ARITH `--y <= c <=> --c <= y`] THEN
	ASM_REWRITE_TAC[] THEN
	CONJ_TAC THEN MATCH_MP_TAC REAL_LE_TRANS THENL
	[
	  EXISTS_TAC `a:real` THEN ASM_REWRITE_TAC[];
	  EXISTS_TAC `d:real` THEN
	    ONCE_REWRITE_TAC[REAL_ARITH `d <= --y <=> y <= --d`] THEN
	    ASM_REWRITE_TAC[]
	];

      MATCH_MP_TAC REAL_LE_TRANS THEN
	EXISTS_TAC `a * --d` THEN
	ASM_REWRITE_TAC[] THEN
	ONCE_REWRITE_TAC[REAL_ARITH `x * y <= a * --d <=> a * d <= x * --y`] THEN
	MATCH_MP_TAC REAL_LE_MUL2 THEN
	ONCE_REWRITE_TAC[REAL_ARITH `d <= --y <=> y <= --d`] THEN
	ASM_REWRITE_TAC[]
    ]);;


(* TT_FF *)
let FLOAT_INTERVAL_MUL_TT_FF = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o prove)(
  `interval_arith x (float T n1 e1, float T n2 e2) /\
    interval_arith y (float F m1 r1, float F m2 r2) /\
    f1 <= float T n1 e1 * float F m2 r2 /\
    float T n2 e2 * float F m1 r1 <= f2
    ==> interval_arith (x * y) (f1, f2)`,
   STRIP_TAC THEN
     MP_TAC ((GEN_ALL o DISCH_ALL) FLOAT_INTERVAL_MUL_FF_TT) THEN
     DISCH_THEN (MP_TAC o SPECL[`n1:num`; `e1:num`; `n2:num`; `e2:num`; `m1:num`; `r1:num`; `m2:num`; `r2:num`]) THEN
     DISCH_THEN (MP_TAC o SPECL[`y:real`; `x:real`; `f1:real`; `f2:real`]) THEN
     REPEAT (POP_ASSUM MP_TAC) THEN
     SIMP_TAC[REAL_MUL_AC]);;



(* TF_FF *)
let FLOAT_INTERVAL_MUL_TF_FF = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o prove)(
  `interval_arith x (float T n1 e1, float F n2 e2) /\
    interval_arith y (float F m1 r1, float F m2 r2) /\
    f1 <= float T n1 e1 * float F m2 r2 /\
    float F n2 e2 * float F m2 r2 <= f2
    ==> interval_arith (x * y) (f1, f2)`,
  REWRITE_TAC[FLOAT_NEG_T] THEN
    REWRITE_TAC[interval_arith] THEN
    MAP_EVERY ABBREV_TAC [`a = float F n1 e1`; `b = float F n2 e2`; `c = float F m1 r1`; `d = float F m2 r2`] THEN
    STRIP_TAC THEN
    SUBGOAL_THEN `&0 <= a /\ &0 <= b /\ &0 <= c /\ &0 <= d` ASSUME_TAC THENL
    [
      MAP_EVERY EXPAND_TAC ["a"; "b"; "c"; "d"] THEN
	REWRITE_TAC[FLOAT_F_POS];
      ALL_TAC
    ] THEN

    SUBGOAL_THEN `&0 <= y` ASSUME_TAC THENL
    [
      MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `c:real` THEN ASM_REWRITE_TAC[];
      ALL_TAC
    ] THEN

    CONJ_TAC THENL
    [
      DISJ_CASES_TAC (REAL_ARITH `&0 <= x \/ &0 <= --x`) THENL
	[
	  MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `&0` THEN
	    CONJ_TAC THENL
	    [
	      MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `--a * d` THEN
		ASM_REWRITE_TAC[REAL_ARITH `--a * d <= &0 <=> &0 <= a * d`] THEN
		MATCH_MP_TAC REAL_LE_MUL THEN ASM_REWRITE_TAC[];
	      MATCH_MP_TAC REAL_LE_MUL THEN ASM_REWRITE_TAC[]
	    ];
	  ALL_TAC
	] THEN

	MATCH_MP_TAC REAL_LE_TRANS THEN	EXISTS_TAC `--a * d` THEN ASM_REWRITE_TAC[] THEN
	ONCE_REWRITE_TAC[REAL_ARITH `--a * d <= x * y <=> --x * y <= a * d`] THEN
	MATCH_MP_TAC REAL_LE_MUL2 THEN
	ONCE_REWRITE_TAC[REAL_ARITH `--x <= a <=> --a <= x`] THEN ASM_REWRITE_TAC[];

      DISJ_CASES_TAC (REAL_ARITH `&0 <= --x \/ &0 <= x`) THENL
	[
	  MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `&0` THEN
	    CONJ_TAC THENL
	    [
	      REWRITE_TAC[REAL_ARITH `x * y <= &0 <=> &0 <= --x * y`] THEN
		MATCH_MP_TAC REAL_LE_MUL THEN ASM_REWRITE_TAC[];
	      MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `b * d` THEN ASM_REWRITE_TAC[] THEN
		MATCH_MP_TAC REAL_LE_MUL THEN ASM_REWRITE_TAC[]
	    ];
	  ALL_TAC
	] THEN

	MATCH_MP_TAC REAL_LE_TRANS THEN	EXISTS_TAC `b * d` THEN ASM_REWRITE_TAC[] THEN
	MATCH_MP_TAC REAL_LE_MUL2 THEN ASM_REWRITE_TAC[]
    ]);;


(* TF_TT *)
let FLOAT_INTERVAL_MUL_TF_TT = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o prove)(
  `interval_arith x (float T n1 e1, float F n2 e2) /\
    interval_arith y (float T m1 r1, float T m2 r2) /\
    f1 <= float F n2 e2 * float T m1 r1 /\
    float T n1 e1 * float T m1 r1 <= f2
    ==> interval_arith (x * y) (f1, f2)`,
  REWRITE_TAC[FLOAT_NEG_T] THEN
    REWRITE_TAC[interval_arith] THEN
    MAP_EVERY ABBREV_TAC [`a = float F n1 e1`; `b = float F n2 e2`; `c = float F m1 r1`; `d = float F m2 r2`] THEN
    STRIP_TAC THEN
    SUBGOAL_THEN `&0 <= a /\ &0 <= b /\ &0 <= c /\ &0 <= d` ASSUME_TAC THENL
    [
      MAP_EVERY EXPAND_TAC ["a"; "b"; "c"; "d"] THEN
	REWRITE_TAC[FLOAT_F_POS];
      ALL_TAC
    ] THEN

    SUBGOAL_THEN `&0 <= --y` ASSUME_TAC THENL
    [
      MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `d:real` THEN 
	ONCE_REWRITE_TAC[REAL_ARITH `d <= --y <=> y <= --d`] THEN ASM_REWRITE_TAC[];
      ALL_TAC
    ] THEN

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

	MATCH_MP_TAC REAL_LE_TRANS THEN	EXISTS_TAC `b * --c` THEN ASM_REWRITE_TAC[] THEN
	ONCE_REWRITE_TAC[REAL_ARITH `b * --c <= x * y <=> x * --y <= b * c`] THEN
	MATCH_MP_TAC REAL_LE_MUL2 THEN ONCE_REWRITE_TAC[REAL_ARITH `--y <= c <=> --c <= y`] THEN
	ASM_REWRITE_TAC[];

      DISJ_CASES_TAC (REAL_ARITH `&0 <= x \/ &0 <= --x`) THENL
	[
	  MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `&0` THEN
	    CONJ_TAC THENL
	    [
	      REWRITE_TAC[REAL_ARITH `x * y <= &0 <=> &0 <= x * --y`] THEN
		MATCH_MP_TAC REAL_LE_MUL THEN ASM_REWRITE_TAC[];
	      MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `--a * --c` THEN ASM_REWRITE_TAC[REAL_NEG_MUL2] THEN
		MATCH_MP_TAC REAL_LE_MUL THEN ASM_REWRITE_TAC[]
	    ];
	  ALL_TAC
	] THEN

	MATCH_MP_TAC REAL_LE_TRANS THEN	EXISTS_TAC `--a * --c` THEN ASM_REWRITE_TAC[REAL_NEG_MUL2] THEN
	ONCE_REWRITE_TAC[REAL_ARITH `x * y <= a * c <=> --x * --y <= a * c`] THEN
	MATCH_MP_TAC REAL_LE_MUL2 THEN ONCE_REWRITE_TAC[REAL_ARITH `--x <= a <=> --a <= x`] THEN
	ASM_REWRITE_TAC[]
    ]);;

(* FF_TF *)
let FLOAT_INTERVAL_MUL_FF_TF = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o prove)(
  `interval_arith x (float F n1 e1, float F n2 e2) /\
    interval_arith y (float T m1 r1, float F m2 r2) /\
    f1 <= float F n2 e2 * float T m1 r1 /\
    float F n2 e2 * float F m2 r2 <= f2
    ==> interval_arith (x * y) (f1, f2)`,
  STRIP_TAC THEN
    MP_TAC ((SPECL [`n1:num`; `e1:num`; `n2:num`; `e2:num`; `m1:num`; `r1:num`; `m2:num`; `r2:num`; `y:real`; `x:real`; `f1:real`; `f2:real`] o GEN_ALL o DISCH_ALL) FLOAT_INTERVAL_MUL_TF_FF) THEN
    ASM_REWRITE_TAC[REAL_MUL_SYM] THEN DISCH_THEN MATCH_MP_TAC THEN
    POP_ASSUM MP_TAC THEN REAL_ARITH_TAC);;

(* TT_TF *)
let FLOAT_INTERVAL_MUL_TT_TF = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o prove)(
  `interval_arith x (float T n1 e1, float T n2 e2) /\
    interval_arith y (float T m1 r1, float F m2 r2) /\
    f1 <= float T n1 e1 * float F m2 r2 /\
    float T n1 e1 * float T m1 r1 <= f2
    ==> interval_arith (x * y) (f1, f2)`,
  STRIP_TAC THEN
    MP_TAC ((SPECL [`n1:num`; `e1:num`; `n2:num`; `e2:num`; `m1:num`; `r1:num`; `m2:num`; `r2:num`; `y:real`; `x:real`; `f1:real`; `f2:real`] o GEN_ALL o DISCH_ALL) FLOAT_INTERVAL_MUL_TF_TT) THEN
    ASM_REWRITE_TAC[REAL_MUL_SYM] THEN POP_ASSUM MP_TAC THEN POP_ASSUM MP_TAC THEN
    SIMP_TAC[REAL_MUL_SYM]);;


(* TF_TF *)
let FLOAT_INTERVAL_MUL_TF_TF = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o prove)(
  `interval_arith x (float T n1 e1, float F n2 e2) /\
    interval_arith y (float T m1 r1, float F m2 r2) /\
    f1_1 <= float T n1 e1 * float F m2 r2 /\
    f1_2 <= float F n2 e2 * float T m1 r1 /\
    min f1_1 f1_2 = f1 /\
    float T n1 e1 * float T m1 r1 <= f2_1 /\
    float F n2 e2 * float F m2 r2 <= f2_2 /\
    max f2_1 f2_2 = f2
    ==> interval_arith (x * y) (f1, f2)`,
  REWRITE_TAC[EQ_SYM_EQ; FLOAT_NEG_T] THEN
    REWRITE_TAC[interval_arith] THEN
    MAP_EVERY ABBREV_TAC [`a = float F n1 e1`; `b = float F n2 e2`; `c = float F m1 r1`; `d = float F m2 r2`] THEN
    STRIP_TAC THEN
    SUBGOAL_THEN `&0 <= a /\ &0 <= b /\ &0 <= c /\ &0 <= d` ASSUME_TAC THENL
    [
      MAP_EVERY EXPAND_TAC ["a"; "b"; "c"; "d"] THEN
	REWRITE_TAC[FLOAT_F_POS];
      ALL_TAC
    ] THEN

    DISJ_CASES_TAC (REAL_ARITH `&0 <= x \/ &0 <= --x`) THENL 
    [
      DISJ_CASES_TAC (REAL_ARITH `&0 <= y \/ &0 <= --y`) THENL 
	[
	  CONJ_TAC THENL 
	    [
	      MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `&0` THEN
		CONJ_TAC THENL 
		[
		  MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `f1_1:real` THEN
		    ASM_REWRITE_TAC[REAL_MIN_MIN] THEN
		    MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `--a * d` THEN
		    ASM_REWRITE_TAC[REAL_ARITH `--a * d <= &0 <=> &0 <= a * d`] THEN
		    MATCH_MP_TAC REAL_LE_MUL THEN ASM_REWRITE_TAC[];
		  ALL_TAC
		] THEN
		MATCH_MP_TAC REAL_LE_MUL THEN ASM_REWRITE_TAC[];
	      ALL_TAC
	    ] THEN
	    
	    MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `b * d` THEN
	    CONJ_TAC THENL
	    [
	      MATCH_MP_TAC REAL_LE_MUL2 THEN ASM_REWRITE_TAC[];
	      ALL_TAC
	    ] THEN
	    MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `f2_2:real` THEN
	    ASM_REWRITE_TAC[REAL_MAX_MAX];
	  ALL_TAC
	] THEN

	CONJ_TAC THENL
	[
	  MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `f1_2:real` THEN
	    ASM_REWRITE_TAC[REAL_MIN_MIN] THEN
	    MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `b * --c` THEN
	    ONCE_REWRITE_TAC[REAL_ARITH `b * --c <= x * y <=> x * --y <= b * c`] THEN ASM_REWRITE_TAC[] THEN
	    MATCH_MP_TAC REAL_LE_MUL2 THEN ONCE_REWRITE_TAC[REAL_ARITH `--y <= c <=> --c <= y`] THEN
	    ASM_REWRITE_TAC[];
	  ALL_TAC
	] THEN
	MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `&0` THEN
	CONJ_TAC THENL
	[
	  REWRITE_TAC[REAL_ARITH `x * y <= &0 <=> &0 <= x * --y`] THEN
	    MATCH_MP_TAC REAL_LE_MUL THEN ASM_REWRITE_TAC[];
	  ALL_TAC
	] THEN
	MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `f2_2:real` THEN
	ASM_REWRITE_TAC[REAL_MAX_MAX] THEN
	MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `b * d` THEN ASM_REWRITE_TAC[] THEN
	MATCH_MP_TAC REAL_LE_MUL THEN ASM_REWRITE_TAC[];
      ALL_TAC
    ] THEN

    DISJ_CASES_TAC (REAL_ARITH `&0 <= y \/ &0 <= --y`) THENL
    [
      CONJ_TAC THENL
	[
	  MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `f1_1:real` THEN
	    ASM_REWRITE_TAC[REAL_MIN_MIN] THEN
	    MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `--a * d` THEN
	    ASM_REWRITE_TAC[] THEN ONCE_REWRITE_TAC[REAL_ARITH `--a * d <= x * y <=> --x * y <= a * d`] THEN
	    MATCH_MP_TAC REAL_LE_MUL2 THEN ONCE_REWRITE_TAC[REAL_ARITH `--x <= a <=> --a <= x`] THEN
	    ASM_REWRITE_TAC[];
	  ALL_TAC
	] THEN
	MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `&0` THEN
	CONJ_TAC THENL
	[
	  REWRITE_TAC[REAL_ARITH `x * y <= &0 <=> &0 <= --x * y`] THEN
	    MATCH_MP_TAC REAL_LE_MUL THEN ASM_REWRITE_TAC[];
	  ALL_TAC
	] THEN
	MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `f2_2:real` THEN
	ASM_REWRITE_TAC[REAL_MAX_MAX] THEN
	MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `b * d` THEN ASM_REWRITE_TAC[] THEN
	MATCH_MP_TAC REAL_LE_MUL THEN ASM_REWRITE_TAC[];
      ALL_TAC
    ] THEN

    CONJ_TAC THENL 
    [
      MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `&0` THEN
	CONJ_TAC THENL 
	[
	  MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `f1_1:real` THEN
	    ASM_REWRITE_TAC[REAL_MIN_MIN] THEN
	    MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `--a * d` THEN
	    ASM_REWRITE_TAC[REAL_ARITH `--a * d <= &0 <=> &0 <= a * d`] THEN
	    MATCH_MP_TAC REAL_LE_MUL THEN ASM_REWRITE_TAC[];
	  ALL_TAC
	] THEN
	ONCE_REWRITE_TAC[GSYM REAL_NEG_MUL2] THEN
	MATCH_MP_TAC REAL_LE_MUL THEN ASM_REWRITE_TAC[];
      ALL_TAC
    ] THEN
	    
    MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `--a * --c` THEN
    CONJ_TAC THENL
    [
      ONCE_REWRITE_TAC[GSYM REAL_NEG_MUL2] THEN REWRITE_TAC[REAL_NEG_NEG] THEN
	MATCH_MP_TAC REAL_LE_MUL2 THEN ONCE_REWRITE_TAC[REAL_ARITH `--x <= a <=> --a <= x`] THEN
	ASM_REWRITE_TAC[];
      ALL_TAC
    ] THEN
    MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `f2_1:real` THEN
    ASM_REWRITE_TAC[REAL_MAX_MAX]);;


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

let float_interval_mul pp th1 th2 =
  let x, l_lo, l_hi = dest_float_interval (concl th1) and
      y, r_lo, r_hi = dest_float_interval (concl th2) in
  let s1, n1, e1 = dest_float l_lo and
      s2, n2, e2 = dest_float l_hi and
      s3, m1, r1 = dest_float r_lo and
      s4, m2, r2 = dest_float r_hi in

    (* Special case *)
    if s1 <> s2 && s3 <> s4 then
      if s1 = "F" or s3 = "F" then
	failwith ("float_interval_mul: FT interval: "^
		    (string_of_term (concl th1))^"; "^(string_of_term (concl th2)))
      else
	let lo1, lo2 = float_mul_lo pp l_lo r_hi, float_mul_lo pp l_hi r_lo and
	    hi1, hi2 = float_mul_hi pp l_lo r_lo, float_mul_hi pp l_hi r_hi in
	let f1_1 = (lhand o concl) lo1 and
	    f1_2 = (lhand o concl) lo2 and
	    f2_1 = (rand o concl) hi1 and
	    f2_2 = (rand o concl) hi2 in
	let min_th = float_min f1_1 f1_2 and
	    max_th = float_max f2_1 f2_2 in
	let f1_tm = (rand o concl) min_th and
	    f2_tm = (rand o concl) max_th in
	let th0 = INST[x, x_var_real; n1, n1_var_num; e1, e1_var_num;
		    y, y_var_real; n2, n2_var_num; e2, e2_var_num;
		    m1, m1_var_num; r1, r1_var_num;
		    m2, m2_var_num; r2, r2_var_num;
		    f1_tm, f1_var_real; f2_tm, f2_var_real;
		    f1_1, f1_1_var; f1_2, f1_2_var; 
		    f2_1, f2_1_var; f2_2, f2_2_var] FLOAT_INTERVAL_MUL_TF_TF in
	  (MY_PROVE_HYP min_th o MY_PROVE_HYP max_th o MY_PROVE_HYP lo1 o MY_PROVE_HYP lo2 o 
	     MY_PROVE_HYP hi1 o MY_PROVE_HYP hi2 o MY_PROVE_HYP th1 o MY_PROVE_HYP th2) th0
    else
      let lo_th, hi_th, th0 =
	if s1 <> s2 then
	  if s1 = "F" then
	    failwith ("float_interval_mul: FT interval: "^
			(string_of_term (concl th1))^"; "^(string_of_term (concl th2)))
	  else
	    if s3 = "F" then
	      float_mul_lo pp l_lo r_hi, float_mul_hi pp l_hi r_hi, FLOAT_INTERVAL_MUL_TF_FF
	    else
	      float_mul_lo pp l_hi r_lo, float_mul_hi pp l_lo r_lo, FLOAT_INTERVAL_MUL_TF_TT
	else
	  if s3 <> s4 then
	    if s3 = "F" then
	      failwith ("float_interval_mul: FT interval: "^
			  (string_of_term (concl th1))^"; "^(string_of_term (concl th2)))
	    else
	      if s1 = "F" then
		float_mul_lo pp l_hi r_lo, float_mul_hi pp l_hi r_hi, FLOAT_INTERVAL_MUL_FF_TF
	      else
		float_mul_lo pp l_lo r_hi, float_mul_hi pp l_lo r_lo, FLOAT_INTERVAL_MUL_TT_TF
	  else
	    if s1 = "F" then
	      if s3 = "F" then
		float_mul_lo pp l_lo r_lo, float_mul_hi pp l_hi r_hi, FLOAT_INTERVAL_MUL_FF_FF
	      else
		float_mul_lo pp l_hi r_lo, float_mul_hi pp l_lo r_hi, FLOAT_INTERVAL_MUL_FF_TT
	    else
	      if s3 = "F" then
		float_mul_lo pp l_lo r_hi, float_mul_hi pp l_hi r_lo, FLOAT_INTERVAL_MUL_TT_FF
	      else
		float_mul_lo pp l_hi r_hi, float_mul_hi pp l_lo r_lo, FLOAT_INTERVAL_MUL_TT_TT in

      let f1_tm = lhand(concl lo_th) and
	  f2_tm = rand(concl hi_th) in

      let th = INST[x, x_var_real; n1, n1_var_num; e1, e1_var_num;
		    y, y_var_real; n2, n2_var_num; e2, e2_var_num;
		    m1, m1_var_num; r1, r1_var_num;
		    m2, m2_var_num; r2, r2_var_num;
		    f1_tm, f1_var_real; f2_tm, f2_var_real] th0 in
	MY_PROVE_HYP lo_th (MY_PROVE_HYP hi_th (MY_PROVE_HYP th1 (MY_PROVE_HYP th2 th)));;



let float_interval_mul_old pp th1 th2 =
  let x, l_lo, l_hi = dest_float_interval (concl th1) and
      y, r_lo, r_hi = dest_float_interval (concl th2) in
  let s1, n1, e1 = dest_float l_lo and
      s2, n2, e2 = dest_float l_hi and
      s3, m1, r1 = dest_float r_lo and
      s4, m2, r2 = dest_float r_hi in

    (* Special case *)
    if s1 <> s2 or s3 <> s4 then
      failwith "float_interval_mul: not implemented"
    else
      let lo_th, hi_th, th0 =
	if s1 = "F" then
	  if s3 = "F" then
	    float_mul_lo pp l_lo r_lo, float_mul_hi pp l_hi r_hi, FLOAT_INTERVAL_MUL_FF_FF
	  else
	    float_mul_lo pp l_hi r_lo, float_mul_hi pp l_lo r_hi, FLOAT_INTERVAL_MUL_FF_TT
	else
	  if s3 = "F" then
	    float_mul_lo pp l_lo r_hi, float_mul_hi pp l_hi r_lo, FLOAT_INTERVAL_MUL_TT_FF
	  else
	    float_mul_lo pp l_hi r_hi, float_mul_hi pp l_lo r_lo, FLOAT_INTERVAL_MUL_TT_TT in

      let f1_tm = lhand(concl lo_th) and
	  f2_tm = rand(concl hi_th) in

      let th = INST[x, x_var_real; n1, n1_var_num; e1, e1_var_num;
		    y, y_var_real; n2, n2_var_num; e2, e2_var_num;
		    m1, m1_var_num; r1, r1_var_num;
		    m2, m2_var_num; r2, r2_var_num;
		    f1_tm, f1_var_real; f2_tm, f2_var_real] th0 in
	MY_PROVE_HYP lo_th (MY_PROVE_HYP hi_th (MY_PROVE_HYP th1 (MY_PROVE_HYP th2 th)));;




(*
let th1 = mk_float_interval_small_num 135;;
let th2 = mk_float_interval_small_num 34;;
let th1' = float_interval_neg th1;;
let th2' = float_interval_neg th2;;

float_interval_mul_old 2 th1 th2;;
float_interval_mul_old 2 th1 th2';;
float_interval_mul_old 2 th1' th2;;
float_interval_mul_old 2 th1' th2';;

float_interval_mul 2 th1 th2;;
float_interval_mul 2 th1 th2';;
float_interval_mul 2 th1' th2;;
float_interval_mul 2 th1' th2';;

let th3 = (NUMERALS_TO_NUM o REWRITE_RULE[FLOAT_OF_NUM; min_exp_def; GSYM FLOAT_NEG_T] o prove)
  (`interval_arith (&0) (-- &5, &10)`, REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;
let th4 = (NUMERALS_TO_NUM o REWRITE_RULE[FLOAT_OF_NUM; min_exp_def; GSYM FLOAT_NEG_T] o prove)
  (`interval_arith (&0) (-- &10, &5)`, REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;

float_interval_mul_old 2 th3 th1;;
float_interval_mul 2 th3 th1;;
float_interval_mul 2 th1 th3;;
float_interval_mul 2 th3 th1';;
float_interval_mul 2 th1' th3;;
float_interval_mul 2 th3 th3;;
float_interval_mul 2 th3 th4;;


(* 10: 0.708 *)
test 1000 (float_interval_mul_old 2 th1) th2;;
(* 10: 0.712 *)
test 1000 (float_interval_mul 2 th1) th2;;

(* 10: 0.504 *)
test 1000 (float_interval_mul 2 th3) th2;;
(* 10: 1.096 *)
test 1000 (float_interval_mul 2 th3) th4;;


let th1 = mk_float_interval_small_num 123456781;;
let th2 = mk_float_interval_small_num 514553443;;
(* 4: 0.340 *)
test 100 (float_interval_mul_old 5 th1) th2;;
(* 4: 0.368 *)
test 100 (float_interval_mul_old 1 th1) th2;;


let tm1 = mk_small_numeral_array 123456781;;
let tm2 = mk_small_numeral_array 514553443;;
let tm = mk_binop mul_op_num tm1 tm2;;

(* 4: 0.148 *)
test 100 NUM_MULT_HASH_CONV tm;;


let th1' = float_interval_round 5 th1;;
let th2' = float_interval_round 5 th2;;

(* 4: 0.176 *)
test 100 (float_interval_mul_old 5 th1') th2';;
*)


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

(* float_interval_div *)

(* FF_FF *)
let FLOAT_INTERVAL_DIV_FF_FF = 
prove(
  `~(m1 = 0) /\
    interval_arith x (float F n1 e1, float F n2 e2) /\
    interval_arith y (float F m1 r1, float F m2 r2) /\
    f1 <= float F n1 e1 / float F m2 r2 /\
    float F n2 e2 / float F m1 r1 <= f2
    ==> interval_arith (x / y) (f1, f2)`,

   MAP_EVERY ABBREV_TAC [`a = float F n1 e1`; `b = float F n2 e2`; `c = float F m1 r1`; `d = float F m2 r2`] THEN
     REWRITE_TAC[real_div] THEN
     STRIP_TAC THEN
     SUBGOAL_THEN `&0 <= a /\ &0 <= b /\ &0 <= c /\ &0 <= d` MP_TAC THENL
     [
       MAP_EVERY EXPAND_TAC ["a";
 "b"; "c"; "d"] THEN
         REWRITE_TAC[FLOAT_F_POS];
       ALL_TAC
     ] THEN
     SUBGOAL_THEN `~(c = &0)` ASSUME_TAC THENL
     [
       EXPAND_TAC "c" THEN ASM_REWRITE_TAC[FLOAT_EQ_0];
       ALL_TAC
     ] THEN
     STRIP_TAC THEN
     SUBGOAL_THEN `~(d = &0)` MP_TAC THENL
     [
       MATCH_MP_TAC (REAL_ARITH `~(c = &0) /\ &0 <= c /\ c <= d ==> ~(d = &0)`) THEN
	 ASM_REWRITE_TAC[] THEN
	 UNDISCH_TAC `interval_arith y (c,d)` THEN
	 REWRITE_TAC[interval_arith] THEN
	 REAL_ARITH_TAC;
       ALL_TAC
     ] THEN

     REPLICATE_TAC 10 (POP_ASSUM MP_TAC) THEN
     REPEAT (POP_ASSUM (fun th -> ALL_TAC)) THEN
     REWRITE_TAC[interval_arith] THEN
     REPEAT STRIP_TAC THENL
     [
       MATCH_MP_TAC REAL_LE_TRANS THEN
         EXISTS_TAC `a * inv d` THEN
         ASM_REWRITE_TAC[] THEN
         MATCH_MP_TAC REAL_LE_MUL2 THEN
         ASM_REWRITE_TAC[REAL_LE_INV_EQ] THEN
         MATCH_MP_TAC REAL_LE_INV2 THEN
         ASM_REWRITE_TAC[] THEN
         MATCH_MP_TAC REAL_LTE_TRANS THEN
         EXISTS_TAC `c:real` THEN
         ASM_REWRITE_TAC[REAL_ARITH `&0 < c <=> ~(c = &0) /\ &0 <= c`];

       MATCH_MP_TAC REAL_LE_TRANS THEN
         EXISTS_TAC `b * inv c` THEN
         ASM_REWRITE_TAC[] THEN
         MATCH_MP_TAC REAL_LE_MUL2 THEN
         ASM_REWRITE_TAC[REAL_LE_INV_EQ] THEN
         REPEAT CONJ_TAC THENL
         [
           MATCH_MP_TAC REAL_LE_TRANS THEN
             EXISTS_TAC `a:real` THEN ASM_REWRITE_TAC[];
           MATCH_MP_TAC REAL_LE_TRANS THEN
             EXISTS_TAC `c:real` THEN ASM_REWRITE_TAC[];
           ALL_TAC
         ] THEN
         MATCH_MP_TAC REAL_LE_INV2 THEN
         ASM_REWRITE_TAC[REAL_ARITH `&0 < c <=> ~(c = &0) /\ &0 <= c`]
     ]);;


(* TT_TT *)
let FLOAT_INTERVAL_DIV_TT_TT = 
prove(
  `~(m2 = 0) /\
    interval_arith x (float T n1 e1, float T n2 e2) /\
    interval_arith y (float T m1 r1, float T m2 r2) /\
    f1 <= float T n2 e2 / float T m1 r1 /\
    float T n1 e1 / float T m2 r2 <= f2
    ==> interval_arith (x / y) (f1,f2)`,

   REWRITE_TAC[FLOAT_NEG_T] THEN
     REWRITE_TAC[real_div; REAL_INV_NEG; REAL_NEG_MUL2] THEN
     GEN_REWRITE_TAC (LAND_CONV o DEPTH_CONV)[interval_arith] THEN
     REWRITE_TAC[REAL_ARITH `--a <= x /\ x <= --b <=> b <= --x /\ --x <= a`] THEN
     GEN_REWRITE_TAC (LAND_CONV o DEPTH_CONV)[GSYM interval_arith] THEN
     STRIP_TAC THEN
     MP_TAC ((SPECL[`n2:num`; `e2:num`; `m1:num`; `r1:num`; `n1:num`; `e1:num`; `m2:num`; `r2:num`; `--x`; `--y`] o GEN_ALL) FLOAT_INTERVAL_DIV_FF_FF) THEN
     REWRITE_TAC[real_div; REAL_INV_NEG; REAL_NEG_MUL2] THEN
     DISCH_THEN MATCH_MP_TAC THEN
     ASM_REWRITE_TAC[]);;
     


(* FF_TT *)
let FLOAT_INTERVAL_DIV_FF_TT = 
prove(
  `~(m2 = 0) /\
    interval_arith x (float F n1 e1, float F n2 e2) /\
    interval_arith y (float T m1 r1, float T m2 r2) /\
    f1 <= float F n2 e2 / float T m2 r2 /\
    float F n1 e1 / float T m1 r1 <= f2
    ==> interval_arith (x / y) (f1,f2)`,

   REWRITE_TAC[FLOAT_NEG_T] THEN
     REWRITE_TAC[real_div; REAL_INV_NEG] THEN
     REWRITE_TAC[REAL_ARITH `a * --b <= c <=> --c <= a * b`] THEN
     REWRITE_TAC[REAL_ARITH `c <= a * --b <=> a * b <= --c`] THEN
     GEN_REWRITE_TAC (LAND_CONV o DEPTH_CONV)[interval_arith] THEN
     REWRITE_TAC[REAL_ARITH `--a <= x /\ x <= --b <=> b <= --x /\ --x <= a`] THEN
     GEN_REWRITE_TAC (LAND_CONV o DEPTH_CONV)[GSYM interval_arith] THEN
     STRIP_TAC THEN
     MP_TAC ((SPECL[`n1:num`; `e1:num`; `m1:num`; `r1:num`; `n2:num`; `e2:num`; `m2:num`; `r2:num`; `x:real`; `--y`; `--f2`; `--f1`] o GEN_ALL) FLOAT_INTERVAL_DIV_FF_FF) THEN
     ANTS_TAC THENL
     [
       ASM_REWRITE_TAC[real_div];
       ALL_TAC
     ] THEN
     REWRITE_TAC[real_div; REAL_INV_NEG; interval_arith] THEN
     REAL_ARITH_TAC);;

(* TT_FF *)
let FLOAT_INTERVAL_DIV_TT_FF = 
prove(
  `~(m1 = 0) /\
    interval_arith x (float T n1 e1, float T n2 e2) /\
    interval_arith y (float F m1 r1, float F m2 r2) /\
    f1 <= float T n1 e1 / float F m1 r1 /\
    float T n2 e2 / float F m2 r2 <= f2
    ==> interval_arith (x / y) (f1,f2)`,

   REWRITE_TAC[FLOAT_NEG_T] THEN
     REWRITE_TAC[real_div; REAL_INV_NEG] THEN
     REWRITE_TAC[REAL_ARITH `--a * b <= c <=> --c <= a * b`] THEN
     REWRITE_TAC[REAL_ARITH `c <= --a * b <=> a * b <= --c`] THEN
     GEN_REWRITE_TAC (LAND_CONV o DEPTH_CONV)[interval_arith] THEN
     REWRITE_TAC[REAL_ARITH `--a <= x /\ x <= --b <=> b <= --x /\ --x <= a`] THEN
     GEN_REWRITE_TAC (LAND_CONV o DEPTH_CONV)[GSYM interval_arith] THEN
     STRIP_TAC THEN
     MP_TAC ((SPECL[`n2:num`; `e2:num`; `m2:num`; `r2:num`; `n1:num`; `e1:num`; `m1:num`; `r1:num`; `--x:real`; `y:real`; `--f2`; `--f1`] o GEN_ALL) FLOAT_INTERVAL_DIV_FF_FF) THEN
     ANTS_TAC THENL
     [
       ASM_REWRITE_TAC[real_div];
       ALL_TAC
     ] THEN
     REWRITE_TAC[real_div; REAL_INV_NEG; interval_arith] THEN
     REAL_ARITH_TAC);;

let FLOAT_0 = 
prove(`!s e. float s 0 e = &0`,
 REWRITE_TAC[FLOAT_EQ_0]);;

(* TF_FF *)
let FLOAT_INTERVAL_DIV_TF_FF = 
prove(
  `~(m1 = 0) /\
    interval_arith x (float T n1 e1, float F n2 e2) /\
    interval_arith y (float F m1 r1, float F m2 r2) /\
    f1 <= float T n1 e1 / float F m1 r1 /\
    float F n2 e2 / float F m1 r1 <= f2
    ==> interval_arith (x / y) (f1,f2)`,

   REWRITE_TAC[FLOAT_NEG_T] THEN STRIP_TAC THEN
     MAP_EVERY ABBREV_TAC [`a = float F n1 e1`; `b = float F n2 e2`; `c = float F m1 r1`; `d = float F m2 r2`] THEN
     SUBGOAL_THEN `&0 <= a /\ &0 <= b /\ &0 <= c /\ &0 <= d` ASSUME_TAC THENL
     [
       MAP_EVERY EXPAND_TAC ["a";
 "b"; "c"; "d"] THEN REWRITE_TAC[FLOAT_F_POS];
       ALL_TAC
     ] THEN

     DISJ_CASES_TAC (REAL_ARITH `&0 <= x \/ x <= &0`) THENL
     [
       SUBGOAL_THEN `interval_arith x (float F 0 0, b)` ASSUME_TAC THENL
	 [
	   UNDISCH_TAC `interval_arith x (--a, b)` THEN POP_ASSUM MP_TAC THEN
	     REWRITE_TAC[interval_arith; FLOAT_0] THEN REAL_ARITH_TAC;
	   ALL_TAC
	 ] THEN

	 MP_TAC ((SPEC_ALL o SPECL [`0`; `0`] o GEN_ALL) FLOAT_INTERVAL_DIV_FF_FF) THEN
	 ASM_REWRITE_TAC[] THEN DISCH_THEN MATCH_MP_TAC THEN
	 MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `&0` THEN
	 REWRITE_TAC[real_div] THEN CONJ_TAC THENL
	 [
	   MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `--a / c` THEN
	     ASM_REWRITE_TAC[real_div; ARITH_RULE `--a * b <= &0 <=> &0 <= a * b`] THEN
	     MATCH_MP_TAC REAL_LE_MUL THEN ASM_REWRITE_TAC[REAL_LE_INV_EQ];
	   ALL_TAC
	 ] THEN

	 MATCH_MP_TAC REAL_LE_MUL THEN ASM_REWRITE_TAC[REAL_LE_INV_EQ; FLOAT_F_POS];
       ALL_TAC
     ] THEN

     SUBGOAL_THEN `interval_arith x (--a, float T 0 0)` ASSUME_TAC THENL
     [
       UNDISCH_TAC `interval_arith x (--a, b)` THEN POP_ASSUM MP_TAC THEN
	 REWRITE_TAC[interval_arith; FLOAT_0] THEN REAL_ARITH_TAC;
       ALL_TAC
     ] THEN

     MP_TAC ((INST[`0`, `n2:num`; `0`, `e2:num`]) FLOAT_INTERVAL_DIV_TT_FF) THEN
     ASM_REWRITE_TAC[FLOAT_NEG_T] THEN DISCH_THEN MATCH_MP_TAC THEN
     ASM_REWRITE_TAC[GSYM FLOAT_NEG_T] THEN MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `&0` THEN
     REWRITE_TAC[real_div] THEN CONJ_TAC THENL
     [
       REWRITE_TAC[FLOAT_0; REAL_MUL_LZERO; REAL_LE_REFL];
       ALL_TAC
     ] THEN

     MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `b / c` THEN ASM_REWRITE_TAC[] THEN
     MATCH_MP_TAC REAL_LE_DIV THEN ASM_REWRITE_TAC[]);;


(* TF_TT *)
let FLOAT_INTERVAL_DIV_TF_TT = 
prove(
  `~(m2 = 0) /\
    interval_arith x (float T n1 e1, float F n2 e2) /\
    interval_arith y (float T m1 r1, float T m2 r2) /\
    f1 <= float F n2 e2 / float T m2 r2 /\
    float T n1 e1 / float T m2 r2 <= f2
    ==> interval_arith (x / y) (f1,f2)`,

   REWRITE_TAC[FLOAT_NEG_T] THEN STRIP_TAC THEN
     MAP_EVERY ABBREV_TAC [`a = float F n1 e1`; `b = float F n2 e2`; `c = float F m1 r1`; `d = float F m2 r2`] THEN
     SUBGOAL_THEN `&0 <= a /\ &0 <= b /\ &0 <= c /\ &0 <= d` ASSUME_TAC THENL
     [
       MAP_EVERY EXPAND_TAC ["a";
 "b"; "c"; "d"] THEN REWRITE_TAC[FLOAT_F_POS];
       ALL_TAC
     ] THEN

     DISJ_CASES_TAC (REAL_ARITH `x <= &0 \/ &0 <= x`) THENL
     [
       SUBGOAL_THEN `interval_arith x (--a, float T 0 0)` ASSUME_TAC THENL
	 [
	   UNDISCH_TAC `interval_arith x (--a, b)` THEN POP_ASSUM MP_TAC THEN
	     REWRITE_TAC[interval_arith; FLOAT_0] THEN REAL_ARITH_TAC;
	   ALL_TAC
	 ] THEN

	 MP_TAC ((SPEC_ALL o SPECL [`0`; `0`] o GEN_ALL) FLOAT_INTERVAL_DIV_TT_TT) THEN
	 ASM_REWRITE_TAC[] THEN DISCH_THEN MATCH_MP_TAC THEN ASM_REWRITE_TAC[FLOAT_NEG_T] THEN
	 ASM_REWRITE_TAC[GSYM FLOAT_NEG_T] THEN
	 MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `&0` THEN
	 REWRITE_TAC[real_div] THEN CONJ_TAC THENL
	 [
	   MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `b / --d` THEN
	     ASM_REWRITE_TAC[real_div; REAL_INV_NEG; REAL_ARITH `b * --d <= &0 <=> &0 <= b * d`] THEN
	     MATCH_MP_TAC REAL_LE_MUL THEN ASM_REWRITE_TAC[REAL_LE_INV_EQ];
	   ALL_TAC
	 ] THEN

	 REWRITE_TAC[FLOAT_0; REAL_MUL_LZERO; REAL_LE_REFL];
       ALL_TAC
     ] THEN

     SUBGOAL_THEN `interval_arith x (float F 0 0, b)` ASSUME_TAC THENL
     [
       UNDISCH_TAC `interval_arith x (--a, b)` THEN POP_ASSUM MP_TAC THEN
	 REWRITE_TAC[interval_arith; FLOAT_0] THEN REAL_ARITH_TAC;
       ALL_TAC
     ] THEN

     MP_TAC ((INST[`0`, `n1:num`; `0`, `e1:num`]) FLOAT_INTERVAL_DIV_FF_TT) THEN
     ASM_REWRITE_TAC[FLOAT_NEG_T] THEN DISCH_THEN MATCH_MP_TAC THEN
     MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `&0` THEN
     REWRITE_TAC[real_div] THEN CONJ_TAC THENL
     [
       REWRITE_TAC[FLOAT_0; REAL_MUL_LZERO; REAL_LE_REFL];
       ALL_TAC
     ] THEN

     MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `--a / --d` THEN ASM_REWRITE_TAC[] THEN
     REWRITE_TAC[real_div; REAL_INV_NEG; REAL_NEG_MUL2] THEN MATCH_MP_TAC REAL_LE_MUL THEN
     ASM_REWRITE_TAC[REAL_LE_INV_EQ]);;




let transform = UNDISCH_ALL o
  PURE_REWRITE_RULE[TAUT `~P <=> (P <=> F)`] o
  REWRITE_RULE[GSYM IMP_IMP] o
  NUMERALS_TO_NUM;;


let FLOAT_INTERVAL_DIV_FF_FF' = transform FLOAT_INTERVAL_DIV_FF_FF and
    FLOAT_INTERVAL_DIV_TT_TT' = transform FLOAT_INTERVAL_DIV_TT_TT and
    FLOAT_INTERVAL_DIV_FF_TT' = transform FLOAT_INTERVAL_DIV_FF_TT and
    FLOAT_INTERVAL_DIV_TT_FF' = transform FLOAT_INTERVAL_DIV_TT_FF and
    FLOAT_INTERVAL_DIV_TF_FF' = transform FLOAT_INTERVAL_DIV_TF_FF and
    FLOAT_INTERVAL_DIV_TF_TT' = transform FLOAT_INTERVAL_DIV_TF_TT;;



let float_interval_div pp th1 th2 =
  let x, l_lo, l_hi = dest_float_interval (concl th1) and
      y, r_lo, r_hi = dest_float_interval (concl th2) in
  let s1, n1, e1 = dest_float l_lo and
      s2, n2, e2 = dest_float l_hi and
      s3, m1, r1 = dest_float r_lo and
      s4, m2, r2 = dest_float r_hi in

    if s3 <> s4 then
      failwith "float_interval_div: division by an interval containing 0"
    else
      let lo_th, hi_th, th0, zero_th =
	if s1 = s2 then
	  if s1 = "F" then
	    if s3 = "F" then
	      float_div_lo pp l_lo r_hi, float_div_hi pp l_hi r_lo, 
	FLOAT_INTERVAL_DIV_FF_FF', raw_eq0_hash_conv m1
	    else
	      float_div_lo pp l_hi r_hi, float_div_hi pp l_lo r_lo, 
	FLOAT_INTERVAL_DIV_FF_TT', raw_eq0_hash_conv m2
	  else
	    if s3 = "F" then
	      float_div_lo pp l_lo r_lo, float_div_hi pp l_hi r_hi, 
	FLOAT_INTERVAL_DIV_TT_FF', raw_eq0_hash_conv m1
	    else
	      float_div_lo pp l_hi r_lo, float_div_hi pp l_lo r_hi, 
	FLOAT_INTERVAL_DIV_TT_TT', raw_eq0_hash_conv m2
	else
	  if s1 = "F" then
	    failwith "float_interval_div: interval (F,T)"
	  else
	    if s3 = "F" then
	      float_div_lo pp l_lo r_lo, float_div_hi pp l_hi r_lo,
	FLOAT_INTERVAL_DIV_TF_FF', raw_eq0_hash_conv m1
	    else
	      float_div_lo pp l_hi r_hi, float_div_hi pp l_lo r_hi,
	FLOAT_INTERVAL_DIV_TF_TT', raw_eq0_hash_conv m2 in

      let f1_tm = lhand(concl lo_th) and
	  f2_tm = rand(concl hi_th) in

      let th = INST[x, x_var_real; n1, n1_var_num; e1, e1_var_num;
		    y, y_var_real; n2, n2_var_num; e2, e2_var_num;
		    m1, m1_var_num; r1, r1_var_num;
		    m2, m2_var_num; r2, r2_var_num;
		    f1_tm, f1_var_real; f2_tm, f2_var_real] th0 in
	(MY_PROVE_HYP lo_th o MY_PROVE_HYP hi_th o 
	   MY_PROVE_HYP th1 o MY_PROVE_HYP th2 o MY_PROVE_HYP zero_th) th;;





(*
let th1 = mk_float_interval_small_num 123456781;;
let th2 = mk_float_interval_small_num 514;;
let th1' = float_interval_neg th1;;
let th2' = float_interval_neg th2;;

let th3 = (NUMERALS_TO_NUM o REWRITE_RULE[FLOAT_OF_NUM; min_exp_def; GSYM FLOAT_NEG_T] o prove)
  (`interval_arith (&0) (-- &5, &10)`, REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;
let th4 = (NUMERALS_TO_NUM o REWRITE_RULE[FLOAT_OF_NUM; min_exp_def; GSYM FLOAT_NEG_T] o prove)
  (`interval_arith (&0) (-- &10, &5)`, REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;

float_interval_div 5 th1 th2;;
float_interval_div 5 th1' th2;;
float_interval_div 5 th1 th2';;
float_interval_div 5 th1' th2';;
float_interval_div 1 th1 th2;;

float_interval_div 5 th3 th1;;
float_interval_div 5 th3 th1';;
float_interval_div 5 th3 th4;;

(* 4: 0.392 *)
test 100 (float_interval_div 5 th1) th2;;
(* 4: 0.224 *)
test 100 (float_interval_div 1 th1) th2;;


let tm1 = mk_small_numeral_array 123456781;;
let tm2 = mk_small_numeral_array 514;;
let tm = mk_binop `DIV` tm1 tm2;;

(* 4: 0.056 *)
test 100 NUM_DIV_HASH_CONV tm;;
NUM_DIV_HASH_CONV tm;;


let th1' = float_interval_round 3 th1;;
let th2' = float_interval_round 3 th2;;

(* 4: 0.256 *)
test 100 (float_interval_div 3 th1') th2';;




let th1 = mk_float_interval_num 1;;
let th2 = mk_float_interval_num 17;;

float_interval_div 5 th1 th2;;
*)



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

(* float_interval_add, float_interval_sub *)

let n1_var_real = `n1:real` and
    n2_var_real = `n2:real` and
    m1_var_real = `m1:real` and
    m2_var_real = `m2:real` and
    n_var_real = `n:real` and
    m_var_real = `m:real`;;


let INTERVAL_ADD = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o prove)(
  `interval_arith x (n1, m1) /\
    interval_arith y (n2, m2) /\
    n <= n1 + n2 /\ m1 + m2 <= m
    ==> interval_arith (x + y) (n, m)`,
   REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;



let INTERVAL_SUB = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o prove)(
  `interval_arith x (n1, m1) /\
    interval_arith y (n2, m2) /\
    n <= n1 - m2 /\ m1 - n2 <= m
    ==> interval_arith (x - y) (n, m)`,
  REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;




let float_interval_add pp th1 th2 =
  let x, n1, m1 = dest_float_interval (concl th1) in
  let y, n2, m2 = dest_float_interval (concl th2) in
  let lo_th = float_add_lo pp n1 n2 in
  let hi_th = float_add_hi pp m1 m2 in
  let n_tm = lhand (concl lo_th) in
  let m_tm = rand (concl hi_th) in
  let th0 = INST[x, x_var_real; n1, n1_var_real; m1, m1_var_real;
                 y, y_var_real; n2, n2_var_real; m2, m2_var_real;
                 n_tm, n_var_real; m_tm, m_var_real] INTERVAL_ADD in
    MY_PROVE_HYP lo_th (MY_PROVE_HYP hi_th (MY_PROVE_HYP th2 (MY_PROVE_HYP th1 th0)));;



let float_interval_sub pp th1 th2 =
  let x, n1, m1 = dest_float_interval (concl th1) in
  let y, n2, m2 = dest_float_interval (concl th2) in
  let lo_th = float_sub_lo pp n1 m2 in
  let hi_th = float_sub_hi pp m1 n2 in
  let n_tm = lhand(concl lo_th) in
  let m_tm = rand(concl hi_th) in
  let th0 = INST[x, x_var_real; n1, n1_var_real; m1, m1_var_real;
		 y, y_var_real; n2, n2_var_real; m2, m2_var_real;
		 n_tm, n_var_real; m_tm, m_var_real] INTERVAL_SUB in
    MY_PROVE_HYP lo_th (MY_PROVE_HYP hi_th (MY_PROVE_HYP th2 (MY_PROVE_HYP th1 th0)));;



(*
let th1 = mk_float_interval_num (Num.num_of_int 123456781);;
let th2 = mk_float_interval_num (Num.num_of_int 514);;
let th3 = float_interval_neg th1;;
float_interval_add 5 th1 th2;;
float_interval_sub 5 th2 th1;;
float_interval_add 1 th1 th2;;

float_interval_add 5 th3 th2;;
float_interval_add 1 th3 th2;;

(* 4: 0.424 *)
test 1000 (float_interval_add 5 th1) th2;;
(* 4: 0.648 *)
test 1000 (float_interval_add 1 th1) th2;;


let tm1 = mk_small_numeral_array 123456781;;
let tm2 = mk_small_numeral_array 514;;
let tm = mk_binop plus_op_num tm1 tm2;;

(* 4: 0.036 *)
test 1000 NUM_ADD_HASH_CONV tm;;
NUM_ADD_HASH_CONV tm;;


let th1' = float_interval_round 3 th1;;
let th2' = float_interval_round 3 th2;;
float_interval_add 3 th1' th2';;

(* 4: 0.724 *)
test 1000 (float_interval_add 3 th1') th2';;
*)

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

(* FLOAT_ABS *)

let s_var_bool = `s:bool`;;

let FLOAT_ABS = 
prove(`abs (float s n e) = float F n e`,

   BOOL_CASES_TAC `s:bool` THEN
     REWRITE_TAC[FLOAT_NEG_T; REAL_ABS_NEG; REAL_ABS_REFL; FLOAT_F_POS]);;

let float_abs tm =
  let ltm, rtm = dest_comb tm in
    if ((fst o dest_const) ltm <> "real_abs") then
      failwith "float_abs: no abs"
    else
      let ltm, e_tm = dest_comb rtm in
      let ltm, n_tm = dest_comb ltm in
      let s_tm = rand ltm in
	INST[s_tm, s_var_bool; n_tm, n_var_num; e_tm, e_var_num] FLOAT_ABS;;
  





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

(* float_interval_sqrt *)

let FLOAT_INTERVAL_SQRT = 
prove(`interval_arith x (float F n1 e1, hi) /\
				  f1 <= sqrt (float F n1 e1) /\ sqrt hi <= f2
				  ==> interval_arith (sqrt x) (f1, f2)`,

   ABBREV_TAC `lo = float F n1 e1` THEN
     REWRITE_TAC[interval_arith] THEN
     STRIP_TAC THEN
     SUBGOAL_THEN `&0 <= lo /\ &0 <= hi` ASSUME_TAC THENL
     [
       EXPAND_TAC "lo" THEN
	 REWRITE_TAC[FLOAT_F_POS] THEN
	 MATCH_MP_TAC REAL_LE_TRANS THEN
	 EXISTS_TAC `x:real` THEN
	 ASM_REWRITE_TAC[] THEN
	 MATCH_MP_TAC REAL_LE_TRANS THEN
	 EXISTS_TAC `lo:real` THEN
	 ASM_REWRITE_TAC[] THEN
	 EXPAND_TAC "lo" THEN
	 REWRITE_TAC[FLOAT_F_POS];
       ALL_TAC
     ] THEN
     SUBGOAL_THEN `&0 <= x` ASSUME_TAC THENL
     [
       MATCH_MP_TAC REAL_LE_TRANS THEN
	 EXISTS_TAC `lo:real` THEN
	 ASM_REWRITE_TAC[];
       ALL_TAC
     ] THEN
     CONJ_TAC THENL
     [
       MATCH_MP_TAC REAL_LE_TRANS THEN
	 EXISTS_TAC `sqrt lo` THEN
	 ASM_REWRITE_TAC[] THEN
	 MATCH_MP_TAC SQRT_MONO_LE THEN
	 ASM_REWRITE_TAC[];
       MATCH_MP_TAC REAL_LE_TRANS THEN
	 EXISTS_TAC `sqrt hi` THEN
	 ASM_REWRITE_TAC[] THEN
	 MATCH_MP_TAC SQRT_MONO_LE THEN
	 ASM_REWRITE_TAC[]
     ]);;

let FLOAT_INTERVAL_SQRT' = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP]) FLOAT_INTERVAL_SQRT;;


let float_interval_sqrt pp th =
  let x_tm, lo_tm, hi_tm = dest_float_interval (concl th) in
  let s1, n1_tm, e1_tm = dest_float lo_tm in
    if s1 <> "F" then
      failwith "float_interval_sqrt: negative low bound"
    else
      let lo_th = float_sqrt_lo pp lo_tm in
      let hi_th = float_sqrt_hi pp hi_tm in
      let f1_tm = lhand (concl lo_th) in
      let f2_tm = rand (concl hi_th) in
      let th0 = INST[x_tm, x_var_real; n1_tm, n1_var_num; e1_tm, e1_var_num;
		     hi_tm, hi_var_real; f1_tm, f1_var_real; 
		     f2_tm, f2_var_real] FLOAT_INTERVAL_SQRT' in
	MY_PROVE_HYP lo_th (MY_PROVE_HYP hi_th (MY_PROVE_HYP th th0));;




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

(* FLOAT_TO_NUM_CONV *)


let FLOAT_TO_NUM_CONV tm =
  let ltm, e_tm = dest_comb tm in
  let f_tm, n_tm = dest_comb ltm in
    if (fst o dest_const o rator) f_tm <> "float" then
      failwith "FLOAT_TO_NUM_CONV"
    else
      let n_th' = SYM (INST[n_tm, n_var_num] Arith_hash.NUM_THM) in
      let e_th' = SYM (INST[e_tm, n_var_num] Arith_hash.NUM_THM) in
      let n_th = TRANS n_th' (NUM_TO_NUMERAL_CONV (mk_comb(Arith_hash.num_const, n_tm))) in
      let e_th = TRANS e_th' (NUM_TO_NUMERAL_CONV (mk_comb(Arith_hash.num_const, e_tm))) in
      let th0 = MK_COMB (AP_TERM f_tm n_th, e_th) in
      let tm0 = rand(concl th0) in
      let th1 = REWRITE_CONV[float; num_exp; REAL_MUL_LID; GSYM REAL_OF_NUM_MUL; GSYM REAL_OF_NUM_POW; min_exp_def] tm0 in
      let th2 = REAL_RAT_REDUCE_CONV (rand(concl th1)) in
	TRANS th0 (TRANS th1 th2);;


end;;


(**************************************)
(* Printer for floating-point numbers *)
(**************************************)

let print_float tm=
  try
    let s, m_tm, e_tm = Arith_float.dest_float tm in
    let m = Arith_hash.raw_dest_hash m_tm and
	e = Arith_hash.raw_dest_hash e_tm -/ Num.num_of_int Arith_options.min_exp in
    let s_str = if s = "T" then "-" else "" in
    let m_str = Num.string_of_num m in
    let e_str = if e = num_0 then "" 
    else "*" ^ string_of_int Arith_options.base ^ "^" ^ Num.string_of_num e in
    let str = "##" ^ s_str ^ m_str ^ e_str in
      print_string str
  with _ -> failwith "print_float";;


install_user_printer ("float", print_float);;