build_tform2.hl

(* ========================================================================== *)
(*      Formal verification of FPTaylor certificates                          *)
(*                                                                            *)
(*      Author: Alexey Solovyev, University of Utah                           *)
(*                                                                            *)
(*      This file is distributed under the terms of the MIT licence           *)
(* ========================================================================== *)

(* -------------------------------------------------------------------------- *)
(* Construction of Taylor forms                                               *)
(* Note: requires the nonlinear inequality verification tool                  *)
(*       https://github.com/monadius/formal_ineqs                             *)
(* -------------------------------------------------------------------------- *)

needs "build_tform.hl";;
needs "bin_float.hl";;
needs "eval_expr.hl";;
needs "rat.hl";;
needs "proof.hl";;

module Build_tform2 = struct

open List;;
open Interval_arith;;
open Arith_float;;
open More_float;;
open Lib;;
open Tform;;
open Rounding;;
open Build_tform;;
open Bin_float;;
open Eval_expr;;
open Proof;;
open Misc_functions;;
open Misc_vars;;

let use_rat_verifier = ref true;;

(* --------------------------------------------- *)
(* Temporary rounding definitions                *)
(* --------------------------------------------- *)

let rnd32_def = new_definition 
  `rnd32(c) = @r. is_rnd_bin (c, ipow2 (-- &24), ipow2 (-- &126)) 
  (real_interval [--ipow2 (&127), ipow2 (&127)]) r`;;

let rnd64_def = new_definition 
  `rnd64(c) = @r. is_rnd_bin (c, ipow2 (-- &53), ipow2 (-- &1022)) 
  (real_interval [--ipow2 (&1023), ipow2 (&1023)]) r`;;

let id_rnd_bin = prove
  (`!c e2 d2 s2. &0 <= e2 /\ &0 <= d2 ==> is_rnd_bin (c,e2,d2) s2 I`,
   REWRITE_TAC[is_rnd_bin; I_THM] THEN REPEAT STRIP_TAC THEN
     MAP_EVERY EXISTS_TAC [`&0`; `&0`] THEN ASM_ARITH_TAC);;

let rnd32_bin = prove
  (`!c. is_rnd_bin (c, ipow2 (-- &24), ipow2 (-- &126)) 
     (real_interval [--ipow2 (&127), ipow2 (&127)]) (rnd32 c)`,
   GEN_TAC THEN REWRITE_TAC[rnd32_def] THEN CONV_TAC SELECT_CONV THEN
     EXISTS_TAC `I:real->real` THEN
     MATCH_MP_TAC id_rnd_bin THEN REWRITE_TAC[ipow2_pos_le]);;

let rnd64_bin = prove
  (`!c. is_rnd_bin (c, ipow2 (-- &53), ipow2 (-- &1022)) 
     (real_interval [--ipow2 (&1023), ipow2 (&1023)]) (rnd64 c)`,
   GEN_TAC THEN REWRITE_TAC[rnd64_def] THEN CONV_TAC SELECT_CONV THEN
     EXISTS_TAC `I:real->real` THEN
     MATCH_MP_TAC id_rnd_bin THEN REWRITE_TAC[ipow2_pos_le]);;

let rnd32 = prove
  (`!c. is_rnd (c, ipow2 (-- &24), ipow2 (-- &126)) 
     (real_interval [--ipow2 (&127), ipow2 (&127)]) (rnd32 c)`,
   GEN_TAC THEN MATCH_MP_TAC is_rnd_bin_is_rnd THEN REWRITE_TAC[rnd32_bin]);;

let rnd64 = prove
  (`!c. is_rnd (c, ipow2 (-- &53), ipow2 (-- &1022)) 
     (real_interval [--ipow2 (&1023), ipow2 (&1023)]) (rnd64 c)`,
   GEN_TAC THEN MATCH_MP_TAC is_rnd_bin_is_rnd THEN REWRITE_TAC[rnd64_bin]);;

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

(* Enumerates elements of a list starting from a given number *)
let gen_enumerate =
  let rec enum n acc = function
    | [] -> rev acc
    | x :: xs -> enum (n + 1) ((n, x) :: acc) xs in
    fun n s ->
      enum n [] s;;

(* Enumerates elements of a list starting from 0 *)
let enumerate s = gen_enumerate 0 s;;

(* Adds a hypothesis to a theorem *)
let add_hyp hyp th =
  let hyp_th = ASSUME hyp in
  let eq = DEDUCT_ANTISYM_RULE hyp_th th in
    EQ_MP eq hyp_th;;

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

let tform_f1_bound = prove
  (`!s h t:(N)tform. approx s h t
   ==> (!x. x IN s ==> 
	  abs (tform_f1 t x) <= sum_list (tform_list t) (\ (f1, e1, e2). abs (f1 x) * e2))`,
   REWRITE_TAC[tform_f1; approx] THEN REPEAT STRIP_TAC THEN
     MATCH_MP_TAC sum_list_abs_le THEN REWRITE_TAC[GSYM ALL_EL] THEN REPEAT STRIP_TAC THEN
     MP_TAC (SPEC `EL i (tform_list (t:(N)tform))` triple_exists) THEN STRIP_TAC THEN
     ASM_REWRITE_TAC[REAL_ABS_MUL] THEN
     MATCH_MP_TAC REAL_LE_MUL2 THEN ASM_REWRITE_TAC[REAL_ABS_POS; REAL_LE_REFL] THEN
     FIRST_X_ASSUM (MP_TAC o SPEC `x:real^N`) THEN ASM_REWRITE_TAC[GSYM ALL_EL] THEN
     STRIP_TAC THEN FIRST_X_ASSUM (MP_TAC o SPEC `i:num`) THEN ASM_REWRITE_TAC[]);;

let opt_exact_th = prove
  (`!s h t e bound total. approx s h t /\ &0 < e /\
     (!x. x IN s ==> sum_list (tform_list t) (\(f1,e1,e2). abs (f1 x) * (e2 / e)) <= bound)
  /\ bound * e <= total
  ==> !x:real^N. x IN s ==> abs (h x - tform_f0 t x) <= total`,
  REPEAT GEN_TAC THEN STRIP_TAC THEN
    MATCH_MP_TAC approx_bound THEN EXISTS_TAC `s:real^N->bool` THEN
    ASM_REWRITE_TAC[SUBSET_REFL] THEN REPEAT STRIP_TAC THEN
    MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `bound * e:real` THEN
    ASM_SIMP_TAC[GSYM REAL_LE_LDIV_EQ] THEN
    REWRITE_TAC[real_div; GSYM sum_list_rmul] THEN
    SUBGOAL_THEN `!y. sum_list (tform_list (t:(N)tform)) 
         (\x. (\ (f1,e1,e2). abs (f1 y) * e2) x * inv e) = 
         sum_list (tform_list t) (\ (f1,e1,e2). abs (f1 y) * e2 / e)` ASSUME_TAC THENL [
      GEN_TAC THEN REPEAT (FIRST [AP_TERM_TAC; AP_THM_TAC]) THEN
        POP_ASSUM (K ALL_TAC) THEN
        REWRITE_TAC[FUN_EQ_THM] THEN GEN_TAC THEN
        STRIP_ASSUME_TAC (SPEC `x:(real^N->real)#(real^N->real)#real` triple_exists) THEN
        ASM_REWRITE_TAC[real_div; REAL_MUL_ASSOC];
      ALL_TAC
    ] THEN
    ASM_SIMP_TAC[]);;

let opt_approx_th = prove
  (`!s h t bounds total. approx s h t /\
     (!x. x IN s ==> ALL2 (\ (f1,e1,e2) b. abs (f1 x) <= b) (tform_list t) bounds) /\
     sum_list (ZIP (tform_list t) bounds) (\ ((f1,e1,e2), b). e2 * b) <= total
     ==> !x:real^N. x IN s ==> abs (h x - tform_f0 t x) <= total`,
   REPEAT GEN_TAC THEN STRIP_TAC THEN
     MATCH_MP_TAC approx_bound THEN EXISTS_TAC `s:real^N->bool` THEN
     ASM_REWRITE_TAC[SUBSET_REFL] THEN REPEAT STRIP_TAC THEN
     MATCH_MP_TAC REAL_LE_TRANS THEN
     EXISTS_TAC `sum_list (ZIP (tform_list (t:(N)tform)) bounds) (\((f1,e1,e2),b). e2 * b)` THEN
     ASM_REWRITE_TAC[] THEN MATCH_MP_TAC sum_list_le_gen THEN
     FIRST_X_ASSUM (MP_TAC o SPEC `x:real^N`) THEN ASM_REWRITE_TAC[] THEN
     REWRITE_TAC[all2_el] THEN REPEAT STRIP_TAC THEN ASM_SIMP_TAC[length_zip; el_zip] THEN
     STRIP_ASSUME_TAC (SPEC `EL i (tform_list (t:(N)tform))` triple_exists) THEN
     FIRST_X_ASSUM (MP_TAC o SPEC `i:num`) THEN ASM_REWRITE_TAC[REAL_MUL_SYM] THEN
     DISCH_TAC THEN MATCH_MP_TAC REAL_LE_LMUL THEN ASM_SIMP_TAC[REAL_LT_IMP_LE] THEN
     UNDISCH_TAC `approx s h (t:(N)tform)` THEN REWRITE_TAC[approx] THEN
     DISCH_THEN (MP_TAC o SPEC `x:real^N`) THEN ASM_REWRITE_TAC[] THEN
     REWRITE_TAC[GSYM ALL_EL] THEN STRIP_TAC THEN
     FIRST_X_ASSUM (MP_TAC o SPEC `i:num`) THEN ASM_REWRITE_TAC[] THEN
     REAL_ARITH_TAC);;

(* --------------------------------------------- *)
(* tform_f1 estimation                           *)
(* --------------------------------------------- *)

let lemma_trans = prove
  (`interval_arith x (lo, hi) /\ a <= x ==> a <= hi`,
   REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;

(* Returns results in the form (x IN s |- tform1 <= b) *)
let eval_tform1_bound pp approx_th dom_vars =
  let th1 = MATCH_MP tform_f1_bound approx_th in
  let th2 = (UNDISCH_ALL o SPEC_ALL) th1 in
  let sum_tm = rand (concl th2) in
  let sum_interval = eval_expr pp sum_tm dom_vars 
    [list_mk; sum_list_cons; sum_list_nil] in
    MATCH_MP lemma_trans (CONJ sum_interval th2);;

(* --------------------------------------------- *)
(* Domain initialization procedures              *)
(* --------------------------------------------- *)

let create_domain =
  let comp_op = `($):real^A->num->real` in
    fun pp list1 list2 ->
      let n = length (dest_list list1) in
      let sub_th0, (f_list1, f_list2) = M_verifier_main.mk_float_domain pp (list1, list2) in
      let dom_tm = lhand (concl sub_th0) in
      let x_tm = mk_var ("x", M_taylor.n_vector_type_array.(n)) in
      let x_comps = 
	let x_app = mk_icomb (comp_op, x_tm) in
	  map (fun i -> mk_comb (x_app, mk_small_numeral i)) (1--n) in
      let x_vec = M_taylor.mk_vector_list x_comps in
      let comp_th = end_itlist CONJ (Array.to_list M_taylor.comp_thms_array.(n)) in
      let x_eq = prove(mk_eq (x_tm, x_vec),
		       REWRITE_TAC[CART_EQ; M_taylor.dimindex_array.(n); 
				   M_taylor.gen_in_interval n; comp_th]) in
      let sub_th = GEN_REWRITE_RULE LAND_CONV [x_eq] 
	(UNDISCH (SPEC x_tm (REWRITE_RULE[SUBSET] sub_th0))) in
      let sub_th1 = REWRITE_RULE[IN_INTERVAL; M_taylor.dimindex_array.(n); 
				 M_taylor.gen_in_interval n] sub_th in
      let interval_ths = REWRITE_RULE[comp_th; GSYM interval_arith] sub_th1 in
      let vars = map (fun th ->
			let tm, _ = dest_interval_arith (concl th) in tm, th) 
	(CONJUNCTS interval_ths) in
	x_tm, dom_tm, vars;;

let build_domain pp vars =
  let names = map (fun v -> v.name) vars and
      ls = map (fun v -> term_of_rat v.low) vars and
      hs = map (fun v -> term_of_rat v.high) vars in
  let name_list = gen_enumerate 1 names in
  let low = mk_list (ls, real_ty) and
      high = mk_list (hs, real_ty) in
  let x_tm, dom_tm, dom_vars = create_domain pp low high in
    x_tm, dom_tm, dom_vars, name_list;;

(* --------------------------------------------- *)
(* Optimization procedures                       *)
(* --------------------------------------------- *)

let gen_vec_eq =
  let comp_op = `($):real^A->num->real` in
    fun var_names ->
      let n = length var_names in
      let x_tm = mk_var ("x", M_taylor.n_vector_type_array.(n)) in
      let x_comps = 
	let x_app = mk_icomb (comp_op, x_tm) in
	  map (fun i -> mk_comb (x_app, mk_small_numeral i)) (1--n) in
      let x_vec = M_taylor.mk_vector_list x_comps in
      let comp_th = end_itlist CONJ (Array.to_list M_taylor.comp_thms_array.(n)) in
      let x_eq = prove(mk_eq (x_tm, x_vec),
		       REWRITE_TAC[CART_EQ; M_taylor.dimindex_array.(n); 
				   M_taylor.gen_in_interval n; comp_th]) in
      let x_list = map (fun i -> mk_var (assoc i var_names, real_ty)) (1--n) in
      let x_vec2 = M_taylor.mk_vector_list x_list in
      let eq_tm = mk_eq (mk_forall (x_tm, mk_icomb (`P:A->bool`, x_tm)), 
			 list_mk_forall (x_list, mk_icomb (`P:A->bool`, x_vec2))) in
	prove(eq_tm,
	      EQ_TAC THEN REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[] THEN
		ONCE_REWRITE_TAC[x_eq] THEN ASM_REWRITE_TAC[]);;

let prove_bound_cond pp tm =
  try
    if !use_rat_verifier then
      let th, _ = Rat.prove_rat_ineq pp bin_float_rat_conv tm in
	th
    else
      let eq_th = bin_float_rat_conv tm in
	if rand (concl eq_th) = t_const then
	  EQT_ELIM eq_th
	else
	  let vars, ineq_tm = strip_forall (rand (concl eq_th)) in
	  let ineq_th, stat = M_verifier_main.verify_ineq M_verifier_main.default_params pp ineq_tm in
	  let ineq_th1 = SPEC_ALL ineq_th in
	  let imp_tm = mk_imp (concl ineq_th1, ineq_tm) in
	  let ineq_th2 = MP (TAUT imp_tm) ineq_th1 in
	    ONCE_REWRITE_RULE[GSYM eq_th] (itlist GEN vars ineq_th2)
  with _ -> error "prove_bound_cond" [tm] [];;

let strip_binop name =
  let rec strip tm =
    match tm with
      | Comb (Comb (Const (op, _), tm1), tm2) when op = name ->
	  tm1 :: strip tm2
      | _ -> [tm] in
    strip;;

let forall_imp_conj = MESON[]
  `(!x. P x ==> A x /\ B x) <=> (!x. P x ==> A x) /\ (!x. P x ==> B x)`;;

let abs_lt_eq =
  REAL_ARITH `!a b. abs a < b <=> a < b /\ --a < b`;;

let abs_le_eq =
  REAL_ARITH `!a b. abs a <= b <=> a <= b /\ --a <= b`;;

let prove_opt_approx pp var_names bounds_tm total_tm approx_th =
  try
    let th0 = REWRITE_RULE[GSYM IMP_IMP] opt_approx_th in
    let th1 = SPECL[bounds_tm; total_tm] (MATCH_MP th0 approx_th) in
    let th2 = REWRITE_RULE[f0_mk; list_mk; ZIP; ALL2; 
			   sum_list_cons; sum_list_nil; REAL_ADD_RID] th1 in
    let n = length var_names in
    let th3 = PURE_ONCE_REWRITE_RULE[gen_vec_eq var_names] th2 in
    let comp_th = end_itlist CONJ (Array.to_list M_taylor.comp_thms_array.(n)) in
    let th4 = PURE_REWRITE_RULE[IN_INTERVAL; M_taylor.dimindex_array.(n); 
				M_taylor.gen_in_interval n; comp_th] th3 in
    let th5 = REWRITE_RULE[forall_imp_conj; FORALL_AND_THM] th4 in
    let bound_conds0 = fst (dest_imp (concl th5)) in
    let bound_eq = REWRITE_CONV[abs_lt_eq; abs_le_eq; forall_imp_conj; 
				FORALL_AND_THM; GSYM CONJ_ASSOC] bound_conds0 in
    let bound_conds = strip_binop "/\\" (rand (concl bound_eq)) in
    let counter = ref 0 in
    let total_bounds = length bound_conds in
    let _ = report (sprintf "\nBounds: %d" total_bounds) in
    let bound_ths = map 
      (fun tm ->
	 let _ = incr counter in
	 let _ = report (sprintf "Verifying: %d/%d" !counter total_bounds) in
	   prove_bound_cond pp tm) bound_conds in
    let bound_th = ONCE_REWRITE_RULE[GSYM bound_eq] (end_itlist CONJ bound_ths) in
    let th6 = MP th5 bound_th in
      CONV_RULE (LAND_CONV bin_float_rat_conv THENC REWRITE_CONV[]) th6
  with Failure _ ->
    error "prove_opt_approx" [bounds_tm; total_tm] [approx_th];;

let prove_opt_exact pp var_names bound_tm e_tm total_tm approx_th =
  try
    let th0 = REWRITE_RULE[GSYM IMP_IMP] opt_exact_th in
    let th1 = SPECL[e_tm; bound_tm; total_tm] (MATCH_MP th0 approx_th) in
    let th2 = REWRITE_RULE[f0_mk; list_mk; sum_list_cons; 
			   sum_list_nil; REAL_ADD_RID] th1 in
    let e_ineq = EQT_ELIM (bin_float_rat_conv (lhand (concl th2))) in
    let th3 = MP th2 e_ineq in
    let th4 = CONV_RULE (LAND_CONV bin_float_rat_conv 
			   THENC REWRITE_CONV[REAL_MUL_RID]) th3 in
    let n = length var_names in
    let th5 = PURE_ONCE_REWRITE_RULE[gen_vec_eq var_names] th4 in
    let comp_th = end_itlist CONJ (Array.to_list M_taylor.comp_thms_array.(n)) in
    let th6 = PURE_REWRITE_RULE[IN_INTERVAL; M_taylor.dimindex_array.(n);
				M_taylor.gen_in_interval n; comp_th] th5 in
    let th7 = REWRITE_RULE[forall_imp_conj; FORALL_AND_THM] th6 in
    let bound_cond = lhand (concl th7) in
    let _ = report "Verifying the exact bound" in
    let rat_flag = !use_rat_verifier in
    let _ = use_rat_verifier := false  in
    let bound_th = prove_bound_cond pp bound_cond in
    let _ = use_rat_verifier := rat_flag in
    let th8 = MP th7 bound_th in
      CONV_RULE (LAND_CONV bin_float_rat_conv THENC REWRITE_CONV[]) th8
  with Failure _ ->
    error "prove_opt_exact" [bound_tm; e_tm; total_tm] [approx_th];;


(* --------------------------------------------- *)
(* Evaluation of tform assumptions               *)
(* --------------------------------------------- *)

let prove_le_expr pp dom_vars h_tm =
  try
    EQT_ELIM (bin_float_rat_conv h_tm)
  with Failure _ ->
    let cond_tm, expr_tm = dest_imp (snd (dest_forall h_tm)) in
    let th0 = EQT_ELIM (eval_expr pp expr_tm dom_vars []) in
    let in_hyp = 
      try
	find (is_binary "IN") (hyp th0)
      with Not_found -> cond_tm in
    let x_tm = lhand in_hyp in
    let th1 = GEN x_tm (DISCH in_hyp th0) in
      if h_tm <> concl th1 then
	error "prove_le_expr" [h_tm] [th1]
      else
	th1;;

let lemma1 = prove
  (`(!x:real^N. x IN s ==> abs (tform_f1 f x) <= r) ==> r <= b
   ==> (!x. x IN s ==> abs (tform_f1 f x) <= b)`,
   REPEAT STRIP_TAC THEN MATCH_MP_TAC REAL_LE_TRANS THEN
     EXISTS_TAC `r:real` THEN ASM_SIMP_TAC[]);;

let prove_bound1 pp dom_vars approx_th h_tm =
  let th0 = eval_tform1_bound pp approx_th dom_vars in
  let in_hyp = find (is_binary "IN") (hyp th0) in
  let x_tm = lhand in_hyp in
  let th1 = GEN x_tm (DISCH in_hyp th0) in
  let th2 = MATCH_MP lemma1 th1 in
  let b_tm = rand (rand (snd (dest_forall h_tm))) in
  let r1 = fst (dest_imp (concl (INST[b_tm, b_var_real] th2))) in
  let r2 = EQT_ELIM ((DEPTH_CONV FLOAT_TO_NUM_CONV THENC bin_float_rat_conv) r1) in
  let th3 = MATCH_MP th2 r2 in
    if h_tm <> concl th3 then
      error "prove_bound1" [h_tm] [th3]
    else
      th3;;

let lemma2 = prove
  (`(!x:real^N. x IN s ==> abs (tform_f1 f1 x) <= r1) ==> 
     (!x. x IN s ==> abs (tform_f1 f2 x) <= r2) ==>
     r1 * r2 <= b
     ==> (!x. x IN s ==> abs (tform_f1 f1 x * tform_f1 f2 x) <= b)`,
   REPEAT STRIP_TAC THEN MATCH_MP_TAC REAL_LE_TRANS THEN
     EXISTS_TAC `r1 * r2` THEN ASM_REWRITE_TAC[REAL_ABS_MUL] THEN
     MATCH_MP_TAC REAL_LE_MUL2 THEN ASM_SIMP_TAC[REAL_ABS_POS]);;

let prove_bound2 pp dom_vars approx1_th approx2_th h_tm = 
  let bound1 = eval_tform1_bound pp approx1_th dom_vars in
  let bound2 = eval_tform1_bound pp approx2_th dom_vars in
  let in_hyp = find (is_binary "IN") (hyp bound1) in
  let x_tm = lhand in_hyp in
  let t1 = GEN x_tm (DISCH in_hyp bound1) in
  let t2 = GEN x_tm (DISCH in_hyp bound2) in
  let th2 = MATCH_MP (MATCH_MP lemma2 t1) t2 in
  let b_tm = rand (rand (snd (dest_forall h_tm))) in
  let r1 = fst (dest_imp (concl (INST[b_tm, b_var_real] th2))) in
  let r2 = EQT_ELIM ((DEPTH_CONV FLOAT_TO_NUM_CONV THENC bin_float_rat_conv) r1) in
  let th3 = MATCH_MP th2 r2 in
    if h_tm <> concl th3 then
      error "prove_bound2" [h_tm] [th3]
    else
      th3;;

let lemma3 = prove
  (`!m1 b. (!x:real^N. x IN s ==> abs (tform_f1 f1 x) <= r1) ==> 
     r1 <= m1 /\ r1 * r1 <= b
     ==> (!x. x IN s ==> abs (tform_f1 f1 x) <= m1 /\ abs (tform_f1 f1 x pow 2) <= b)`,
   REPEAT STRIP_TAC THEN MATCH_MP_TAC REAL_LE_TRANS THENL [
     EXISTS_TAC `r1:real` THEN ASM_SIMP_TAC[];
     EXISTS_TAC `r1 * r1` THEN ASM_REWRITE_TAC[GSYM REAL_POW_2] THEN
       REWRITE_TAC[REAL_ABS_POW; GSYM REAL_LE_SQUARE_ABS] THEN
       POP_ASSUM (fun th -> FIRST_X_ASSUM (MP_TAC o C MATCH_MP th)) THEN
       REAL_ARITH_TAC
   ]);;

let prove_bound3 pp dom_vars approx_th h_tm =
  let bound = eval_tform1_bound pp approx_th dom_vars in
  let in_hyp = find (is_binary "IN") (hyp bound) in
  let x_tm = lhand in_hyp in
  let t = GEN x_tm (DISCH in_hyp bound) in
  let th2 = MATCH_MP lemma3 t in
  let tm1, tm2 = dest_conj (rand (snd (dest_forall h_tm))) in
  let m1_tm = rand tm1 and
      b_tm = rand tm2 in
  let r1 = fst (dest_imp (concl (SPECL[m1_tm; b_tm] th2))) in
  let r2 = EQT_ELIM ((DEPTH_CONV FLOAT_TO_NUM_CONV THENC bin_float_rat_conv) r1) in
  let th3 = MATCH_MP th2 r2 in
    if h_tm <> concl th3 then
      error "prove_bound3" [h_tm] [th3; approx_th]
    else
      th3;;

let interval_abs_lemma = (REWRITE_RULE[GSYM IMP_IMP] o prove)
  (`!y t a b lo. abs y <= t /\ interval_arith t (a, b) /\ lo = --b
     ==> interval_arith y (lo, b)`,
   REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;

let prove_ineq_bound pp dom_vars h_tm =
  let vars, h_tm1 = strip_forall h_tm in
  let x_in, ineqs0 = dest_imp h_tm1 in
  let y_ineq, ineqs1 = dest_imp ineqs0 in
  let y_bound0 = bin_float_interval pp (rand y_ineq) in
  let lo_eq = float_neg (rand (rand (concl y_bound0))) in
  let y_bound = MATCH_MP (MATCH_MP (MATCH_MP interval_abs_lemma (ASSUME y_ineq)) y_bound0) lo_eq in
  let dom_vars2 = (`y:real`, y_bound) :: dom_vars in
  let r0 = EQT_ELIM (eval_expr pp ineqs1 dom_vars2 [IN_REAL_INTERVAL]) in
  let r1 = add_hyp x_in r0 in
  let r2 = itlist GEN vars (DISCH_ALL (DISCH y_ineq r1)) in
    if h_tm <> concl r2 then
      error "prove_inv_bound" [h_tm] [r2]
    else
      r2;;

let prove_rnd_bin_var_hyp var_names tm =
  let n = length var_names in
  let eq_th1 = (PURE_ONCE_REWRITE_CONV[gen_vec_eq var_names] THENC 
		  REWRITE_CONV[IN_REAL_INTERVAL] THENC
		  INT_REDUCE_CONV THENC NUM_REDUCE_CONV THENC
		  bin_float_rat_conv) tm in
  let tm1 = rand (concl eq_th1) in
  let comp_th = end_itlist CONJ (Array.to_list M_taylor.comp_thms_array.(n)) in
  let eq_th2 = PURE_REWRITE_CONV[IN_INTERVAL; M_taylor.dimindex_array.(n); 
				 M_taylor.gen_in_interval n; comp_th] tm1 in
  let r1 = REAL_ARITH (rand (concl eq_th2)) in
  let r2 = EQ_MP (SYM eq_th2) r1 in
    EQ_MP (SYM eq_th1) r2;;

let prove_forall_hyp =
  let table = [
    `!x:real^N. x IN s ==> abs (tform_f1 f x) <= b`,
    (fun pp var_names dom_vars arg_ths h_tm ->
       prove_bound1 pp dom_vars (hd arg_ths) h_tm);
    `!x:real^N. x IN s ==> abs (tform_f1 f1 x * tform_f1 f2 x) <= b`,
    (fun pp var_names dom_vars arg_ths h_tm ->
       prove_bound2 pp dom_vars (nth arg_ths 0) (nth arg_ths 1) h_tm);
    `!x:real^N. x IN s ==> h x IN (:real)`,
    (fun pp var_names dom_vars arg_ths h_tm ->
       prove(h_tm, REWRITE_TAC[IN_UNIV]));
    `!x:real^N. x IN s ==> abs a + abs b <= c`,
    (fun pp var_names dom_vars arg_ths h_tm ->
       prove_le_expr pp dom_vars h_tm);
    `!x:real^N. x IN s ==> abs (tform_f1 f x) <= m1 /\ abs (tform_f1 f x pow 2) <= b`,
    (fun pp var_names dom_vars arg_ths h_tm ->
       prove_bound3 pp dom_vars (hd arg_ths) h_tm);
    `!x:real^N y. x IN s ==> abs y <= m1 ==> P x y`,
    (fun pp var_names dom_vars arg_ths h_tm ->
       prove_ineq_bound pp dom_vars h_tm);
    `!x:real^N. x IN s ==> abs (x$i) <= m1 /\ P x`,
    (fun pp var_names dom_Vars arg_ths h_tm ->
       prove_rnd_bin_var_hyp var_names h_tm);
  ] in
    fun pp var_names dom_vars arg_ths h_tm ->
      let _, f = find (fun (tm, _) -> can (term_match [] tm) h_tm) table in
	f pp var_names dom_vars arg_ths h_tm;;

let prove_other_hyp =
  let table = [
    `HIDDEN (a:A) (x:C) = HIDDEN (b:B) y`,
    (fun pp dom_Vars arg_ths h_tm ->
       let eq_th = PURE_REWRITE_CONV[hidden_def] h_tm in
       let eq = rand (concl eq_th) in
       let tm1, tm2 = dest_eq eq in
       let def_eq = prove_err_def_eq tm1 tm2 in
	 prove(h_tm, REWRITE_TAC[hidden_def; def_eq]))
  ] in
    fun pp dom_vars arg_ths h_tm ->
      let _, f = find (fun (tm, _) -> can (term_match [] tm) h_tm) table in
	f pp dom_vars arg_ths h_tm;;

let prove_rat_hyp tm =
  let conv = 
    REWRITE_CONV[IN_REAL_INTERVAL] THENC
      NUM_REDUCE_CONV THENC INT_REDUCE_CONV THENC
      bin_float_rat_conv in
  EQT_ELIM (conv tm);;

let get_rnd_thm bin_flag rnd =
  let c_bin_tm = mk_bin_float rnd.coefficient in
  let c_tm = rand (concl (bin_float_rat_conv c_bin_tm)) in
  let th0 = 
    match rnd.bits with
      | 32 -> if bin_flag then rnd32_bin else rnd32
      | 64 -> if bin_flag then rnd64_bin else rnd64
      | _ -> failwith "get_rnd_thm: Unsupported rnd mode" in
    SPEC c_tm th0;;

(* --------------------------------------------- *)
(* The main tform construction procedures        *)
(* --------------------------------------------- *)

let build_tform pp proof =
  let _ = reset_index() in
  let x_tm, dom_tm, dom_vars, var_names = build_domain pp proof.proof_vars in
  let steps = sort (fun s1 s2 -> compare s1.step_index s2.step_index) proof.proof_steps in
  let total_steps = length steps in
  let _ = report (sprintf "Proof steps: %d" total_steps) in
  let step_counter = ref 0 in
  let rec replay acc steps =
    match steps with
      | [] -> acc
      | step :: rest ->
	  let _ = incr step_counter in
	  let args = step.proof_args in
	  let arg_ths = map (fun i -> assoc i acc) args.arg_indices in
	  let err_indices = args.err_indices in
	  let th0 =
	    begin
	      match step.proof_op with
		| Proof_var name ->
		    let _ = report (sprintf "%d/%d: Var(%s)" 
				      !step_counter total_steps name) in
		      build_var_tform dom_tm (rev_assoc name var_names)
		| Proof_rnd_bin_var (rnd, name) ->
		    let _ = report (sprintf "%d/%d: Var_rnd(%s, rnd = %d)" 
				      !step_counter total_steps name rnd.bits) in
		    let n_tm = mk_intconst (Int (int_of_float (nth args.bounds 0))) and
			b_tm = mk_bin_float (nth args.bounds 1) in
		    let rnd_th = get_rnd_thm true rnd in
		    let index = rev_assoc name var_names in
		      build_rnd_bin_var_tform dom_tm index rnd_th n_tm b_tm err_indices
		| Proof_const c ->
		    let _ = report (sprintf "%d/%d: Const(%s)" 
				      !step_counter total_steps (string_of_num c)) in
		      build_const_tform dom_tm (term_of_rat c)
		| Proof_rnd_bin_const (rnd, c) ->
		    let _ = report (sprintf "%d/%d: Const_rnd(%s, rnd = %d)"
				      !step_counter total_steps (string_of_num c) rnd.bits) in
		    let n_tm = mk_intconst (Int (int_of_float (nth args.bounds 0))) and
			b_tm = mk_bin_float (nth args.bounds 1) in
		    let rnd_th = get_rnd_thm true rnd in
		      build_rnd_bin_const_tform dom_tm (term_of_rat c) rnd_th n_tm b_tm err_indices
		| Proof_rnd rnd ->
		    let _ = report (sprintf "%d/%d: Rnd(%d, %f)" 
				      !step_counter total_steps rnd.bits rnd.coefficient) in
		    let m2_tm = mk_bin_float (nth args.bounds 0) and
			b_tm = mk_bin_float (nth args.bounds 1) in
		    let arg = nth arg_ths 0 in
		    let rnd_th = get_rnd_thm false rnd in
		      build_rnd_tform rnd_th m2_tm b_tm err_indices arg
		| Proof_neg ->
		    let _ = report (sprintf "%d/%d: Neg" !step_counter total_steps) in
		    let arg = nth arg_ths 0 in
		      build_neg_tform arg
		| Proof_add ->
		    let _ = report (sprintf "%d/%d: Add" !step_counter total_steps) in
		    let arg1 = nth arg_ths 0 and
			arg2 = nth arg_ths 1 in
		      build_add_tform arg1 arg2
		| Proof_sub ->
		    let _ = report (sprintf "%d/%d: Sub" !step_counter total_steps) in
		    let arg1 = nth arg_ths 0 and
			arg2 = nth arg_ths 1 in
		      build_sub_tform arg1 arg2
		| Proof_mul ->
		    let _ = report (sprintf "%d/%d: Mul" !step_counter total_steps) in
		    let m2_tm = mk_bin_float (nth args.bounds 0) and
			e_tm = mk_ipow2 (Int (int_of_float (nth args.bounds 1))) in
		    let arg1 = nth arg_ths 0 and
			arg2 = nth arg_ths 1 in
		      build_mul_tform m2_tm e_tm err_indices arg1 arg2
		| Proof_inv ->
		    let _ = report (sprintf "%d/%d: Inv" !step_counter total_steps) in
		    let m1_tm = mk_bin_float (nth args.bounds 0) and
			m2_tm = mk_bin_float (nth args.bounds 1) and
			e2_tm = mk_ipow2 (Int (int_of_float (nth args.bounds 2))) and
			b_tm = mk_bin_float (nth args.bounds 3) and
			m3_tm = mk_bin_float (nth args.bounds 4) in
		    let arg = nth arg_ths 0 in
		      build_inv_tform m1_tm m2_tm e2_tm b_tm m3_tm err_indices arg
		| Proof_sqrt ->
(*
		    let _ = error "Not implemented" 
		    (map mk_bin_float args.bounds @ map mk_small_numeral err_indices)
		      arg_ths in
*)
		    let _ = report (sprintf "%d/%d: Sqrt" !step_counter total_steps) in
		    let m1_tm = mk_bin_float (nth args.bounds 0) and
			m2_tm = mk_bin_float (nth args.bounds 1) and
			e2_tm = mk_ipow2 (Int (int_of_float (nth args.bounds 2))) and
			b_tm = mk_bin_float (nth args.bounds 3) and
			m3_tm = mk_bin_float (nth args.bounds 4) in
		    let arg = nth arg_ths 0 in
		      build_sqrt_tform m1_tm m2_tm e2_tm b_tm m3_tm err_indices arg
		| Proof_sin ->
		    let _ = report (sprintf "%d/%d: Sin" !step_counter total_steps) in
		    let m1_tm = mk_bin_float (nth args.bounds 0) and
			m2_tm = mk_bin_float (nth args.bounds 1) and
			e2_tm = mk_ipow2 (Int (int_of_float (nth args.bounds 2))) and
			b_tm = mk_bin_float (nth args.bounds 3) and
			m3_tm = mk_bin_float (nth args.bounds 4) in
		    let arg = nth arg_ths 0 in
		      build_sin_tform m1_tm m2_tm e2_tm b_tm m3_tm err_indices arg
		| Proof_cos ->
		    let _ = report (sprintf "%d/%d: Cos" !step_counter total_steps) in
		    let m1_tm = mk_bin_float (nth args.bounds 0) and
			m2_tm = mk_bin_float (nth args.bounds 1) and
			e2_tm = mk_ipow2 (Int (int_of_float (nth args.bounds 2))) and
			b_tm = mk_bin_float (nth args.bounds 3) and
			m3_tm = mk_bin_float (nth args.bounds 4) in
		    let arg = nth arg_ths 0 in
		      build_cos_tform m1_tm m2_tm e2_tm b_tm m3_tm err_indices arg
		| Proof_exp ->
		    let _ = report (sprintf "%d/%d: Exp" !step_counter total_steps) in
		    let m1_tm = mk_bin_float (nth args.bounds 0) and
			m2_tm = mk_bin_float (nth args.bounds 1) and
			e2_tm = mk_ipow2 (Int (int_of_float (nth args.bounds 2))) and
			b_tm = mk_bin_float (nth args.bounds 3) and
			m3_tm = mk_bin_float (nth args.bounds 4) in
		    let arg = nth arg_ths 0 in
		      build_exp_tform m1_tm m2_tm e2_tm b_tm m3_tm err_indices arg
		| Proof_log ->
		    let _ = report (sprintf "%d/%d: Log" !step_counter total_steps) in
		    let m1_tm = mk_bin_float (nth args.bounds 0) and
			m2_tm = mk_bin_float (nth args.bounds 1) and
			e2_tm = mk_ipow2 (Int (int_of_float (nth args.bounds 2))) and
			b_tm = mk_bin_float (nth args.bounds 3) and
			m3_tm = mk_bin_float (nth args.bounds 4) in
		    let arg = nth arg_ths 0 in
		      build_log_tform m1_tm m2_tm e2_tm b_tm m3_tm err_indices arg
		| Proof_atn ->
		    let _ = report (sprintf "%d/%d: Atn" !step_counter total_steps) in
		    let m1_tm = mk_bin_float (nth args.bounds 0) and
			m2_tm = mk_bin_float (nth args.bounds 1) and
			e2_tm = mk_ipow2 (Int (int_of_float (nth args.bounds 2))) and
			b_tm = mk_bin_float (nth args.bounds 3) and
			m3_tm = mk_bin_float (nth args.bounds 4) in
		    let arg = nth arg_ths 0 in
		      build_atn_tform m1_tm m2_tm e2_tm b_tm m3_tm err_indices arg
		| Proof_simpl_eq (i, j) ->
		    let _ = report (sprintf "%d/%d: Simpl_eq" !step_counter total_steps) in
		    let arg = nth arg_ths 0 in
		      build_simpl_eq_tform i j arg
		| Proof_simpl_add (i, j) ->
		    let _ = report (sprintf "%d/%d: Simpl_add" !step_counter total_steps) in
		    let b_tm = mk_bin_float (nth args.bounds 0) and
			e_tm = mk_ipow2 (Int (int_of_float (nth args.bounds 1))) in
		    let arg = nth arg_ths 0 in
		      build_simpl_add_tform i j b_tm e_tm err_indices arg
(*		| _ -> error "Not implemented" 
		    (map mk_bin_float args.bounds @ map mk_small_numeral err_indices)
		      arg_ths *)
	    end in
	  let hs_forall, hs = partition is_forall (hyp th0) in
	  let hs_rat, hs_other = partition (fun tm -> frees tm = []) hs in
	  let th1 = 
	    begin
	      try
		let rat_ths = map prove_rat_hyp hs_rat in
		let forall_ths = map (prove_forall_hyp pp var_names dom_vars arg_ths) hs_forall in
		let other_ths = map (prove_other_hyp pp dom_vars arg_ths) hs_other in
		let r = itlist MY_PROVE_HYP (forall_ths @ other_ths @ rat_ths) th0 in
		  if hyp r <> [] then failwith "non-empty assumptions" else r
	      with 
		| Failure msg -> error ("build_tform: " ^ msg) (hyp th0) (th0 :: arg_ths)
		| Not_found -> error "Not_found" (hyp th0) (th0 :: arg_ths)
	    end in
	    replay ((step.step_index, th1) :: acc) rest
  in
    replay [] steps;;

let get_bounds_tm bounds indices approx_th =
  let _, _, t_tm = dest_approx (concl approx_th) in
  let _, f1_tm = dest_mk_tform t_tm in
  let indices0 = map extract_index (dest_list f1_tm) in
  let indices1 = map (fun i -> index i indices0) indices in
  let bnds = sort (fun (i1, _) (i2, _) -> compare i1 i2) (zip indices1 bounds) in
    (* Check that all bounds are present *)
  let _ = map2 (fun (i1, _) i2 -> 
		  if i1 <> i2 then error "get_bounds_tm" bounds [approx_th])
    bnds (0--(length bounds - 1)) in
    mk_list (map snd bnds, real_ty);;

let prove_bound pp proof approx_th =
  let x_tm, dom_tm, dom_vars, var_names = build_domain pp proof.proof_vars in
  let opt = hd proof.proof_opts in
    match opt.opt_type with
      | Proof_opt_approx ->
	  let bounds = map mk_bin_float opt.opt_bounds in
	  let total_tm = mk_bin_float opt.total_bound in
	  let bounds_tm = get_bounds_tm bounds opt.opt_indices approx_th in
	    prove_opt_approx pp var_names bounds_tm total_tm approx_th
      | Proof_opt_exact ->
	  let bound, e_exp = list_to_pair opt.opt_bounds in
	  let bound_tm = mk_bin_float bound in
	  let e_tm = mk_ipow2 (Int (int_of_float e_exp)) in
	  let total_tm = mk_bin_float opt.total_bound in
	    prove_opt_exact pp var_names bound_tm e_tm total_tm approx_th;;

let validate_proof pp proof =
  let ths = build_tform pp proof in
  let approx_th = snd (hd ths) in
  let bound_th =
    if proof.proof_opts = [] then
      TRUTH
    else
      prove_bound pp proof approx_th in
    approx_th, bound_th;;

let validate_proof_gen pp proof total_bound =
  let ths = build_tform pp proof in
  let approx_th = snd (hd ths) in
  let bound_th =
    if proof.proof_opts = [] then
      TRUTH
    else
      let opt = hd proof.proof_opts in
	match opt.opt_type with
	  | Proof_opt_approx ->
	      prove_bound pp proof approx_th
	  | Proof_opt_exact ->
	      let x_tm, dom_tm, dom_vars, var_names = build_domain pp proof.proof_vars in
	      let _, e_exp = list_to_pair opt.opt_bounds in
	      let bound = ldexp total_bound (-int_of_float e_exp) in
	      let bound_tm = mk_bin_float bound in
	      let e_tm = mk_ipow2 (Int (int_of_float e_exp)) in
	      let total_tm = mk_bin_float total_bound in
		prove_opt_exact pp var_names bound_tm e_tm total_tm approx_th in	          approx_th, bound_th;;

end;;