(* ========================================================================== *)
(* FLYSPECK - BOOK FORMALIZATION                                              *)
(*                                                                            *)
(* Chapter: Jordan                                                               *)
(* Copied from HOL Light jordan directory *)
(* Author: Thomas C. Hales                                                    *)
(* Date: 2010-07-08                                                           *)
(* ========================================================================== *)

(*
A handheld calculator, implemented in HOL light.

It implements + - / * sqrt, ssqrt, abs, atn, acs, pi,
The user can supply additional conversons to extend the calculator to other functions.

There is an example file float_example.hl.
*)

(* based on float.ml in hol/jordan/float.ml *)

(*
To do. Allow the user to supply theorems that carry domain constraints,
 (such as ACS_ATN).
*)

flyspeck_needs "jordan/taylor_atn.hl";;

module Float = struct

open Lib_ext;;
open Num_ext_nabs;;
open Real_ext;;
open Refinement;;
open Tactics_jordan;;
open Taylor_atn;;

Parse_ext_override_interface.unambiguous_interface();;
Parse_ext_override_interface.prioritize_real();;

(* Three types of error:
     Insufficient precision: the inequality might be true but there is not enough
     precision to decide.  Divide by zero generally lands here, because there is
     no testing for exact equality.

     Internal failure: a bug in this code.

     User input failure: sqrt of negative numbers, false inequalities.

*)

exception Insufficient_precision;;

(*
exception Mk_comb of (term*hol_type)*(term*hol_type);;

let mk_comb(a,b) = try mk_comb(a,b) with 
  Failure _ -> raise (Mk_comb((a,type_of a),(b,type_of b)));;
*)
let add_test,test = new_test_suite();;

let twopow =
  new_definition(
        `twopow x = if (?n. (x = (int_of_num n)))
        then ((&2) pow (nabs x))
        else inv((&2) pow (nabs x))`);;
let float =
  new_definition(
                  `float x n = (real_of_int x)*(twopow n)`);;
let interval_eps =
  new_definition(
                   `interval_eps x f eps = ((abs (x-f)) <= eps)`);;
(*--------------------------------------------------------------------*) let mk_interval a b ex = mk_comb(mk_comb (mk_comb (`interval_eps`,a),b),ex);; add_test("mk_interval", mk_interval `#3` `#4` `#1` = `interval_eps #3 #4 #1`);; let dest_interval intv = let (h1,ex) = dest_comb intv in let (h2,f) = dest_comb h1 in let (h3,a) = dest_comb h2 in let _ = (h3 = `interval_eps`) or (print_term h3; failwith "dest_interval:interval_eps") in (a,f,ex);; let _ = add_test("dest_interval", let a = `#3` and b = `#4` and c = `#1` in dest_interval (mk_interval a b c) = (a,b,c));; (*--------------------------------------------------------------------*) let (dest_int:term-> Num.num) = fun b -> let dest_pos_int a = let (op,nat) = dest_comb a in if (fst (dest_const op) = "int_of_num") then (dest_numeral nat) else fail() in let (op',u) = (dest_comb b) in try (if (fst (dest_const op') = "int_neg") then minus_num (dest_pos_int u) else dest_pos_int b) with Failure _ -> failwith "dest_int ";; let (mk_int:Num.num -> term) = fun a -> let sgn = Num.sign_num a in let abv = Num.abs_num a in let r = mk_comb(`int_of_num`,mk_numeral abv) in try (if (sgn<0) then mk_comb (`int_neg`,r) else r) with Failure _ -> failwith ("dest_int "^(string_of_num a));; add_test("mk_int", (mk_int (Int (-1443)) = `int_neg (int_of_num 1443)`) && (mk_int (Int 37) = `(int_of_num 37)`));; (* ------------------------------------------------------------------ *) let (split_ratio:Num.num -> Num.num*Num.num) = function (Ratio r) -> (Big_int (Ratio.numerator_ratio r)), (Big_int (Ratio.denominator_ratio r))| u -> (u,(Int 1));; add_test("split_ratio", let (a,b) = split_ratio ((Int 4)//(Int 20)) in (a =/ (Int 1)) && (b =/ (Int 5)));; (* ------------------------------------------------------------------ *) (* break nonzero int r into a*(C**b) with a prime to C . *) let (factor_C:int -> Num.num -> Num.num*Num.num) = function c -> let intC = (Int c) in let rec divC (a,b) = if ((Int 0) =/ mod_num a intC) then (divC (a//intC,b+/(Int 1))) else (a,b) in function r-> if ((Num.is_integer_num r)&& not((Num.sign_num r) = 0)) then divC (r,(Int 0)) else failwith "factor_C";; add_test("factor_C", (factor_C 2 (Int (4096+32)) = (Int 129,Int 5)) && (factor_C 10 (Int (5000)) = (Int 5,Int 3)) && (cannot (factor_C 2) ((Int 50)//(Int 3))));; (*--------------------------------------------------------------------*) let (dest_float:term -> Num.num) = fun f -> let (a,b) = dest_binop `float` f in (dest_int a)*/ ((Int 2) **/ (dest_int b));; add_test("dest_float", dest_float `float (int_of_num 3) (int_of_num 17)` = (Int 393216));; add_test("dest_float2", (* must express as numeral first *) cannot dest_float `float ((int_of_num 3)+:(int_of_num 1)) (int_of_num 17)`);; (* ------------------------------------------------------------------ *) (* creates float of the form `float a b` with a odd *) let (mk_float:Num.num -> term) = function r -> let (a,b) = split_ratio r in let (a',exp_a) = if (a=/(Int 0)) then ((Int 0),(Int 0)) else factor_C 2 a in let (b',exp_b) = factor_C 2 b in let c = a'//b' in if (Num.is_integer_num c) then mk_binop `float` (mk_int c) (mk_int (exp_a -/ exp_b)) else failwith "mk_float";; add_test("mk_float", mk_float (Int (4096+32)) = `float (int_of_num 129) (int_of_num 5)` && (mk_float (Int 0) = `float (int_of_num 0) (int_of_num 0)`));; add_test("mk_float2", (* throws exception exactly when denom != 2^k *) let rtest = fun t -> (t =/ dest_float (mk_float t)) in rtest ((Int 3)//(Int 1024)) && (cannot rtest ((Int 1)//(Int 3))));; add_test("mk_float dest_float", (* constructs canonical form of float *) mk_float (dest_float `float (int_of_num 4) (int_of_num 3)`) = `float (int_of_num 1) (int_of_num 5)`);; (* ------------------------------------------------------------------ *) (* creates decimal of the form `DECIMAL a b` with a prime to 10 *) let (mk_pos_decimal:Num.num -> term) = let zero = `#0` in function r -> let _ = assert (r >=/ (Int 0)) in let (a,b) = split_ratio r in if (a=/(Int 0)) then zero else let (a1,exp_a5) = factor_C 5 a in let (a2,exp_a2) = factor_C 2 a1 in let (b1,exp_b5) = factor_C 5 b in let (b2,exp_b2) = factor_C 2 b1 in let _ = assert(b2 =/ (Int 1)) in let c = end_itlist Num.max_num [exp_b5-/exp_a5;exp_b2-/exp_a2;(Int 0)] in let a' = a2*/((Int 2)**/ (c +/ exp_a2 -/ exp_b2))*/ ((Int 5)**/(c +/ exp_a5 -/ exp_b5)) in let b' = (Int 10) **/ c in mk_binop `DECIMAL` (mk_numeral a') (mk_numeral b');; add_test("mk_pos_decimal", mk_pos_decimal (Int (5000)) = `#5000` && (mk_pos_decimal ((Int 30)//(Int 40)) = `#0.75`) && (mk_pos_decimal (Int 0) = `#0`) && (mk_pos_decimal ((Int 15)//(Int 25)) = `#0.6`) && (mk_pos_decimal ((Int 25)//(Int 4)) = `#6.25`) && (mk_pos_decimal ((Int 2)//(Int 25)) = `#0.08`));; let (mk_decimal:Num.num->term) = function r -> let a = Num.sign_num r in let b = mk_pos_decimal (Num.abs_num r) in if (a < 0) then (mk_comb (`real_neg`, b)) else b;; add_test("mk_decimal", (mk_decimal (Int 3) = `#3`) && (mk_decimal (Int (-3)) = `real_neg (#3)`));; (*--------------------------------------------------------------------*) let (dest_decimal:term -> Num.num) = fun f -> let (a,b) = dest_binop `DECIMAL` f in let a1 = dest_numeral a in let b1 = dest_numeral b in a1//b1;; add_test("dest_decimal", dest_decimal `#3.4` =/ ((Int 34)//(Int 10)));; add_test("dest_decimal2", cannot dest_decimal `real_neg (#3.4)`);; (*--------------------------------------------------------------------*) (* Properties of integer powers of 2. *) (* ------------------------------------------------------------------ *)
let TWOPOW_POS = 
prove(`!n. (twopow (int_of_num n) = (&2) pow n)`,
(REWRITE_TAC[twopow]) THEN GEN_TAC THEN COND_CASES_TAC THENL [AP_TERM_TAC;ALL_TAC] THEN (REWRITE_TAC[NABS_POS]) THEN (UNDISCH_EL_TAC 0) THEN (TAUT_TAC (` ( A ) ==> (~ A ==> B)`)) THEN (MESON_TAC[]));;
let TWOPOW_NEG = 
prove(`!n. (twopow (--(int_of_num n)) = inv((&2) pow n))`,
GEN_TAC THEN (REWRITE_TAC[twopow]) THEN (COND_CASES_TAC THENL [ALL_TAC;REWRITE_TAC[NABS_NEG]]) THEN (POP_ASSUM CHOOSE_TAC) THEN (REWRITE_TAC[NABS_NEG]) THEN (UNDISCH_EL_TAC 0) THEN (REWRITE_TAC[int_eq;int_neg_th;INT_NUM_REAL]) THEN (REWRITE_TAC[prove (`! u y.((--(real_of_num u) = (real_of_num y))= ((real_of_num u) +(real_of_num y) = (&0)))`,REAL_ARITH_TAC)]) THEN (REWRITE_TAC[REAL_OF_NUM_ADD;REAL_OF_NUM_EQ;ADD_EQ_0]) THEN (DISCH_TAC) THEN (ASM_REWRITE_TAC[real_pow;REAL_INV_1]));;
let TWOPOW_INV = 
prove(`!a. (twopow (--: a) = (inv (twopow a)))`,
(GEN_TAC) THEN ((ASSUME_TAC (SPEC `a:int` INT_REP2))) THEN ((POP_ASSUM CHOOSE_TAC)) THEN ((POP_ASSUM DISJ_CASES_TAC)) THEN ((ASM_REWRITE_TAC[TWOPOW_POS;TWOPOW_NEG;REAL_INV_INV;INT_NEG_NEG])));;
let INT_REP3 = 
prove(`!a .(?n.( (a = int_of_num n) \/ (a = --: (int_of_num (n+1)))))`,
(GEN_TAC) THEN ((ASSUME_TAC (SPEC `a:int` INT_REP2))) THEN ((POP_ASSUM CHOOSE_TAC)) THEN ((DISJ_CASES_TAC (prove (`((a:int) = (int_of_num 0)) \/ ~((a:int) =(int_of_num 0))`, MESON_TAC[])))) (* cases *) THENL[ ((EXISTS_TAC `0`)) THEN ((ASM_REWRITE_TAC[]));ALL_TAC] THEN ((UNDISCH_EL_TAC 0)) THEN ((POP_ASSUM DISJ_CASES_TAC)) THENL [DISCH_TAC THEN ((ASM MESON_TAC)[]);ALL_TAC] THEN (DISCH_TAC) THEN ((EXISTS_TAC `PRE n`)) THEN ((DISJ2_TAC)) THEN ((ASM_REWRITE_TAC[INT_EQ_NEG2])) (*** Changed by JRH, 2006/03/28 to avoid PRE_ELIM_TAC ***) THEN (FIRST_X_ASSUM(MP_TAC o check(is_neg o concl))) THEN (ASM_REWRITE_TAC[INT_NEG_EQ_0; INT_OF_NUM_EQ]) THEN (ARITH_TAC));;
let REAL_EQ_INV = 
prove(`!x y. ((x:real = y) <=> (inv(x) = inv (y)))`,
((REPEAT GEN_TAC)) THEN (EQ_TAC) THENL [((DISCH_TAC THEN (ASM_REWRITE_TAC[]))); (* branch 2*) ((DISCH_TAC)) THEN ((ONCE_REWRITE_TAC [(GSYM REAL_INV_INV)])) THEN ((ASM_REWRITE_TAC[]))]);;
let TWOPOW_ADD_1 = 
prove(`!a. (twopow (a +: (int_of_num 1)) = twopow (a) *. (twopow (int_of_num 1)))`,
EVERY[ GEN_TAC; CHOOSE_TAC (SPEC `a:int` INT_REP3); POP_ASSUM DISJ_CASES_TAC THENL[ ASM_REWRITE_TAC[TWOPOW_POS;INT_OF_NUM_ADD;REAL_POW_ADD]; EVERY[ ASM_REWRITE_TAC[GSYM INT_OF_NUM_ADD;INT_NEG_ADD;GSYM INT_ADD_ASSOC;INT_ADD_LINV;INT_ADD_RID]; REWRITE_TAC[GSYM INT_NEG_ADD;INT_OF_NUM_ADD;TWOPOW_NEG;TWOPOW_POS]; ONCE_REWRITE_TAC[SPEC `(&. 2) pow 1` (GSYM REAL_INV_INV)]; REWRITE_TAC[GSYM REAL_INV_MUL;GSYM REAL_EQ_INV;REAL_POW_ADD;GSYM REAL_MUL_ASSOC;REAL_POW_1]; REWRITE_TAC[MATCH_MP REAL_MUL_RINV (REAL_ARITH `~((&. 2) = (&. 0))`); REAL_MUL_RID] ] ] ]);;
let TWOPOW_0 = 
prove(`twopow (int_of_num 0) = (&. 1)`,
(REWRITE_TAC[TWOPOW_POS;real_pow]));;
let TWOPOW_SUB_NUM = 
prove(`!m n.( twopow((int_of_num m) - (int_of_num n)) = twopow((int_of_num m))*. twopow(--: (int_of_num n)))`,
((INDUCT_TAC)) THENL [REWRITE_TAC[INT_SUB_LZERO;REAL_MUL_LID;TWOPOW_0];ALL_TAC] THEN ((INDUCT_TAC THEN ( (ASM_REWRITE_TAC[INT_SUB_RZERO;TWOPOW_0;REAL_MUL_RID;INT_NEG_0;ADD1;GSYM INT_OF_NUM_ADD])))) THEN ((ASM_REWRITE_TAC [TWOPOW_ADD_1;TWOPOW_INV;prove (`((int_of_num m)+(int_of_num 1)) -: ((int_of_num n) +: (int_of_num 1)) = ((int_of_num m)-: (int_of_num n))`,INT_ARITH_TAC)])) THEN ((REWRITE_TAC[REAL_INV_MUL])) THEN ((ABBREV_TAC `a:real = twopow (int_of_num m)`)) THEN ((ABBREV_TAC `b:real = inv(twopow (int_of_num n))`)) THEN ((REWRITE_TAC[TWOPOW_POS;REAL_POW_1;GSYM REAL_MUL_ASSOC;prove (`!(x:real). ((&.2)*x = x*(&.2))`,REAL_ARITH_TAC)])) THEN ((REWRITE_TAC[REAL_INV2;REAL_MUL_RID])));;
let TWOPOW_ADD_NUM = 
prove( `!m n. (twopow ((int_of_num m) + (int_of_num n)) = twopow((int_of_num m))*. twopow((int_of_num n)))`,
let TWOPOW_ADD_INT = 
prove( `!a b. (twopow (a +: b) = twopow(a) *. (twopow(b)))`,
((REPEAT GEN_TAC)) THEN ((ASSUME_TAC (SPEC `a:int` INT_REP))) THEN ((POP_ASSUM CHOOSE_TAC)) THEN ((POP_ASSUM CHOOSE_TAC)) THEN ((ASSUME_TAC (SPEC `b:int` INT_REP))) THEN ((REPEAT (POP_ASSUM CHOOSE_TAC))) THEN ((ASM_REWRITE_TAC[])) THEN ((SUBGOAL_MP_TAC `int_of_num n -: int_of_num m +: int_of_num n' -: int_of_num m' = (int_of_num (n+n')) -: (int_of_num (m+m'))`)) (* branch *) THENL[ ((REWRITE_TAC[GSYM INT_OF_NUM_ADD])) THEN ((INT_ARITH_TAC));ALL_TAC] (* 2nd *) THEN ((DISCH_TAC)) THEN ((ASM_REWRITE_TAC[TWOPOW_SUB_NUM;TWOPOW_INV;TWOPOW_POS;REAL_POW_ADD;REAL_INV_MUL;GSYM REAL_MUL_ASSOC])) THEN ((ABBREV_TAC `a':real = inv ((&. 2) pow m)`)) THEN ((ABBREV_TAC `c :real = (&. 2) pow n`)) THEN ((ABBREV_TAC `d :real = (&. 2) pow n'`)) THEN ((ABBREV_TAC `e :real = inv ((&. 2) pow m')`)) THEN (MESON_TAC[REAL_MUL_AC]));;
let TWOPOW_ABS = 
prove(`!a. real_abs (twopow a) = (twopow a)`,
(GEN_TAC) THEN ((CHOOSE_THEN DISJ_CASES_TAC (SPEC `a:int` INT_REP2))) (* branch *) THEN ((ASM_REWRITE_TAC[TWOPOW_POS;TWOPOW_NEG;REAL_ABS_POW;REAL_ABS_NUM;REAL_ABS_INV])));;
let REAL_POW_POW = 
prove( `!x m n . (x **. m) **. n = x **. (m *| n)`,
((GEN_TAC THEN GEN_TAC THEN INDUCT_TAC)) (* branch *) THENL[ ((REWRITE_TAC[real_pow;MULT_0])); (* second branch *) ((REWRITE_TAC[real_pow])) THEN ((ASM_REWRITE_TAC[ADD1;LEFT_ADD_DISTRIB;REAL_POW_ADD;REAL_MUL_AC;MULT_CLAUSES]))]);;
let INT_POW_POW = INT_OF_REAL_THM REAL_POW_POW;;
let TWOPOW_POW = 
prove( `!a n. (twopow a) pow n = twopow (a *: (int_of_num n))`,
((REPEAT GEN_TAC)) THEN ((CHOOSE_THEN DISJ_CASES_TAC (SPEC `a:int` INT_REP2))) (* branch *) THEN ((ASM_REWRITE_TAC[TWOPOW_POS;INT_OF_NUM_MUL; REAL_POW_POW;TWOPOW_NEG;REAL_POW_INV;INT_OF_NUM_MUL;GSYM INT_NEG_LMUL])));;
(* ------------------------------------------------------------------ *) (* Arithmetic operations on float *) (* ------------------------------------------------------------------ *)
let FLOAT_NEG = 
prove(`!a m. real_neg (float a m) = float (--: a) m`,
REPEAT GEN_TAC THEN REWRITE_TAC[float;GSYM REAL_MUL_LNEG;int_neg_th]);;
let FLOAT_MUL = 
prove(`!a b m n. (float a m) *. (float b n) = (float (a *: b) (m +: n))`,
((REPEAT GEN_TAC)) THEN ((REWRITE_TAC[float;GSYM REAL_MUL_ASSOC;TWOPOW_ADD_INT;int_mul_th])) THEN ((MESON_TAC[REAL_MUL_AC])));;
let FLOAT_ADD = 
prove( `!a b c m. (float a (m+: (int_of_num c))) +. (float b m) = (float ( (int_of_num (2 EXP c))*a +: b) m)`,
((REWRITE_TAC[float;int_add_th;REAL_ADD_RDISTRIB;int_mul_th;TWOPOW_ADD_INT])) THEN ((REPEAT GEN_TAC)) THEN ((REWRITE_TAC[TWOPOW_POS;INT_NUM_REAL;GSYM REAL_OF_NUM_POW])) THEN ((MESON_TAC[REAL_MUL_AC])));;
let FLOAT_ADD_EQ = 
prove( `!a b m. (float a m) +. (float b m) = (float (a+:b) m)`,
((REPEAT GEN_TAC)) THEN ((REWRITE_TAC[REWRITE_RULE[INT_ADD_RID] (SPEC `m:int` (SPEC `0` (SPEC `b:int` (SPEC `a:int` FLOAT_ADD))))])) THEN ((REWRITE_TAC[EXP;INT_MUL_LID])));;
let FLOAT_ADD_NP = 
prove( `!a b m n. (float b (--:(int_of_num n))) +. (float a (int_of_num m)) = (float a (int_of_num m)) +. (float b (--:(int_of_num n)))`,
(REWRITE_TAC[REAL_ADD_AC]));;
let FLOAT_ADD_PN = 
prove( `!a b m n. (float a (int_of_num m)) +. (float b (--(int_of_num n))) = (float ( (int_of_num (2 EXP (m+| n)))*a + b) (--:(int_of_num n)))`,
((REPEAT GEN_TAC)) THEN ((SUBGOAL_MP_TAC `int_of_num m = (--:(int_of_num n)) + (int_of_num (m+n))`)) THENL[ ((REWRITE_TAC[GSYM INT_OF_NUM_ADD])) THEN ((INT_ARITH_TAC)); (* branch *) ((DISCH_TAC)) THEN ((ASM_REWRITE_TAC[FLOAT_ADD]))]);;
let FLOAT_ADD_PP = 
prove( `!a b m n. ((n <=| m) ==>( (float a (int_of_num m)) +. (float b (int_of_num n)) = (float ((int_of_num (2 EXP (m -| n))) *a + b) (int_of_num n))))`,
((REPEAT GEN_TAC)) THEN (DISCH_TAC) THEN ((SUBGOAL_MP_TAC `int_of_num m = (int_of_num n) + (int_of_num (m-n))`)) THENL[ ((REWRITE_TAC[INT_OF_NUM_ADD])) THEN (AP_TERM_TAC) THEN ((REWRITE_TAC[prove (`!(m:num) n. (n+m-n) = (m-n)+n`,REWRITE_TAC[ADD_AC])])) THEN ((UNDISCH_EL_TAC 0)) THEN ((MATCH_ACCEPT_TAC(GSYM SUB_ADD))); (* branch *) ((DISCH_TAC)) THEN ((ASM_REWRITE_TAC[FLOAT_ADD]))]);;
let FLOAT_ADD_PPv2 = 
prove( `!a b m n. ((m <| n) ==>( (float a (int_of_num m)) +. (float b (int_of_num n)) = (float ((int_of_num (2 EXP (n -| m))) *b + a) (int_of_num m))))`,
((REPEAT GEN_TAC)) THEN (DISCH_TAC) THEN ((H_MATCH_MP (THM (prove(`!m n. m<|n ==> m <=|n`,MESON_TAC[LT_LE]))) (HYP_INT 0))) THEN ((UNDISCH_EL_TAC 0)) THEN ((SIMP_TAC[GSYM FLOAT_ADD_PP])) THEN (DISCH_TAC) THEN ((REWRITE_TAC[REAL_ADD_AC])));;
let FLOAT_ADD_NN = 
prove( `!a b m n. ((n <=| m) ==>( (float a (--:(int_of_num m))) +. (float b (--:(int_of_num n))) = (float ((int_of_num (2 EXP (m -| n))) *b + a) (--:(int_of_num m)))))`,
((REPEAT GEN_TAC)) THEN (DISCH_TAC) THEN ((SUBGOAL_MP_TAC `--: (int_of_num n) = --: (int_of_num m) + (int_of_num (m-n))`)) THENL [((UNDISCH_EL_TAC 0)) THEN ((SIMP_TAC [INT_OF_REAL_THM (GSYM REAL_OF_NUM_SUB)])) THEN (DISCH_TAC) THEN ((INT_ARITH_TAC)); (*branch*) ((DISCH_TAC)) THEN (ASM_REWRITE_TAC[GSYM FLOAT_ADD;REAL_ADD_AC])]);;
let FLOAT_ADD_NNv2 = 
prove( `!a b m n. ((m <| n) ==>( (float a (--:(int_of_num m))) +. (float b (--:(int_of_num n))) = (float ((int_of_num (2 EXP (n -| m))) *a + b) (--:(int_of_num n)))))`,
((REPEAT GEN_TAC)) THEN (DISCH_TAC) THEN (((H_MATCH_MP (THM (prove(`!m n. m<|n ==> m <=|n`,MESON_TAC[LT_LE]))) (HYP_INT 0)))) THEN (((UNDISCH_EL_TAC 0))) THEN (((SIMP_TAC[GSYM FLOAT_ADD_NN]))) THEN ((DISCH_TAC)) THEN (((REWRITE_TAC[REAL_ADD_AC]))));;
let FLOAT_SUB = 
prove( `!a b n m. (float a n) -. (float b m) = (float a n) +. (float (--: b) m)`,
let FLOAT_ABS = 
prove( `!a n. real_abs (float a n) = (float (||: a) n)`,
let FLOAT_POW = 
prove( `!a n m. (float a n) **. m = (float (a **: m) (n *: (int_of_num m)))`,
let INT_SUB = 
prove( `!a b. (a -: b) = (a +: (--: b))`,
(REWRITE_TAC[GSYM INT_SUB_RNEG;INT_NEG_NEG]));;
let INT_ABS_NUM = 
prove( `!n. ||: (int_of_num n) = (int_of_num n)`,
let INT_ABS_NEG_NUM = 
prove( `!n. ||: (--: (int_of_num n)) = (int_of_num n)`,
let INT_ADD_NEG_NUM = 
prove(`!x y. --: (int_of_num x) +: (int_of_num y) = (int_of_num y) +: (--: (int_of_num x))`,
(REWRITE_TAC[INT_ADD_AC]));;
let INT_POW_MUL = INT_OF_REAL_THM REAL_POW_MUL;;
let INT_POW_NEG1 = 
prove ( `!x n. (--: (int_of_num x)) **: n = ((--: (int_of_num 1)) **: n) *: ((int_of_num x) **: n)`,
(REWRITE_TAC[GSYM INT_POW_MUL; GSYM INT_NEG_MINUS1]));;
let pow = real_pow;;
let POW_1 = 
prove( `!x. x pow 1 = x`,
GEN_TAC THEN REWRITE_TAC[num_CONV `1`] THEN REWRITE_TAC[pow; REAL_MUL_RID]);;
let POW_2 = 
prove( `!x. x pow 2 = x * x`,
GEN_TAC THEN REWRITE_TAC[num_CONV `2`] THEN REWRITE_TAC[pow; POW_1]);;
let INT_POW_EVEN_NEG1 = 
prove( `!x n. (--: (int_of_num x)) **: (NUMERAL (BIT0 n)) = ((int_of_num x) **: (NUMERAL (BIT0 n)))`,
((REPEAT GEN_TAC)) THEN ((ONCE_REWRITE_TAC[INT_POW_NEG1])) THEN ((ABBREV_TAC `a = int_of_num 1`)) THEN ((ABBREV_TAC `b = (int_of_num x)**: (NUMERAL (BIT0 n))`)) THEN ((REWRITE_TAC[NUMERAL;BIT0])) THEN ((REWRITE_TAC[GSYM MULT_2;GSYM INT_POW_POW;INT_OF_REAL_THM REAL_POW_2;INT_NEG_MUL2])) THEN ((EXPAND_TAC "a")) THEN ((REWRITE_TAC[INT_MUL_RID;INT_MUL_LID;INT_OF_REAL_THM REAL_POW_ONE])));;
let POW_MINUS1 = 
prove( `!n. (--(&1)) pow (2 * n) = &1`,
INDUCT_TAC THEN REWRITE_TAC[MULT_CLAUSES; pow] THEN REWRITE_TAC[num_CONV `2`; num_CONV `1`; ADD_CLAUSES] THEN REWRITE_TAC[pow] THEN REWRITE_TAC[SYM(num_CONV `2`); SYM(num_CONV `1`)] THEN ASM_REWRITE_TAC[] THEN REWRITE_TAC[GSYM REAL_NEG_LMUL; GSYM REAL_NEG_RMUL] THEN REWRITE_TAC[REAL_MUL_LID; REAL_NEGNEG]);;
let POW_1 = 
prove( `!x. x pow 1 = x`,
GEN_TAC THEN REWRITE_TAC[num_CONV `1`] THEN REWRITE_TAC[pow; REAL_MUL_RID]);;
let INT_POW_ODD_NEG1 = 
prove( `!x n. (--: (int_of_num x)) **: (NUMERAL (BIT1 n)) = --: ((int_of_num x) **: (NUMERAL (BIT1 n)))`,
((REPEAT GEN_TAC)) THEN ((ONCE_REWRITE_TAC[INT_POW_NEG1])) THEN (((ABBREV_TAC `a = int_of_num 1`))) THEN (((ABBREV_TAC `b = (int_of_num x)**: (NUMERAL (BIT1 n))`))) THEN ((REWRITE_TAC[NUMERAL;BIT1])) THEN ((ONCE_REWRITE_TAC[ADD1])) THEN ((EXPAND_TAC "a")) THEN ((REWRITE_TAC[GSYM MULT_2])) THEN ((REWRITE_TAC[INT_OF_REAL_THM POW_MINUS1;INT_OF_REAL_THM REAL_POW_ADD])) THEN ((REWRITE_TAC[INT_OF_REAL_THM POW_1;INT_MUL_LID;INT_MUL_LNEG])));;
(* subtraction of integers *)
let INT_ADD_NEG = 
prove( `!x y. (x < y ==> ((int_of_num x) +: (--: (int_of_num y)) = --: (int_of_num (y - x))))`,
((REPEAT GEN_TAC)) THEN ((DISCH_TAC)) THEN ((SUBGOAL_MP_TAC `int_of_num (y-x ) = (int_of_num y) - (int_of_num x)`)) THENL [((SUBGOAL_MP_TAC `x <=| y`)) (* branch *) THENL [(((ASM MESON_TAC)[LE_LT]));((SIMP_TAC[GSYM (INT_OF_REAL_THM REAL_OF_NUM_SUB)]))] (* branch *) ; ((DISCH_TAC)) THEN ((ASM_REWRITE_TAC[])) THEN (ACCEPT_TAC(INT_ARITH `int_of_num x +: --: (int_of_num y) = --: (int_of_num y -: int_of_num x)`))]);;
let INT_ADD_NEGv2 = 
prove( `!x y. (y <= x ==> ((int_of_num x) +: (--: (int_of_num y)) = (int_of_num (x - y))))`,
((REPEAT GEN_TAC)) THEN ((DISCH_TAC)) THEN ((SUBGOAL_MP_TAC `int_of_num (x - y) = (int_of_num x) - (int_of_num y)`)) THENL[ ((UNDISCH_EL_TAC 0)) THEN ((SIMP_TAC[GSYM (INT_OF_REAL_THM REAL_OF_NUM_SUB)])); ((DISCH_TAC)) THEN ((ASM_REWRITE_TAC[INT_SUB])) ] );;
(* assumes a term of the form int_of_num a +: (--: (int_of_num b)) *) let INT_SUB_CONV t = let a,b = dest_binop `int_add` t in let (_,a) = dest_comb a in let (_,b) = dest_comb b in let (_,b) = dest_comb b in let a = dest_numeral a in let b = dest_numeral b in let thm = if (b <=/ a) then INT_ADD_NEGv2 else INT_ADD_NEG in (ARITH_SIMP_CONV[thm]) t;; (* (SIMP_CONV[thm;ARITH]) t;; *)
let FLOAT_EQ = 
prove( `!a b a' b'. (float a b = (float a' b')) <=> ((float a b) -. (float a' b') = (&.0))`,
MESON_TAC[REAL_SUB_0]);;
let FLOAT_LT = 
prove( `!a b a' b'. (float a b <. (float a' b')) <=> ((&.0) <. (float a' b') -. (float a b))`,
MESON_TAC[REAL_SUB_LT]);;
let FLOAT_LE = 
prove( `!a b a' b'. (float a b <=. (float a' b')) <=> ((&.0) <=. (float a' b') -. (float a b))`,
MESON_TAC[REAL_SUB_LE]);;
(* ------------------------------------------------------------------ *) (* Simplifies an arithmetic expression in floats *) (* A workhorse *) (* ------------------------------------------------------------------ *) let FLOAT_CONV = (ARITH_SIMP_CONV[FLOAT_LT;FLOAT_LE;FLOAT_MUL;FLOAT_SUB;FLOAT_ABS;FLOAT_POW; FLOAT_ADD_NN;FLOAT_ADD_NNv2;FLOAT_ADD_PP;FLOAT_ADD_PPv2; FLOAT_ADD_NP;FLOAT_ADD_PN;FLOAT_NEG; INT_NEG_NEG;INT_SUB; INT_ABS_NUM;INT_ABS_NEG_NUM; INT_MUL_LNEG;INT_MUL_RNEG;INT_NEG_MUL2;INT_OF_NUM_MUL; INT_OF_NUM_ADD;GSYM INT_NEG_ADD;INT_ADD_NEG_NUM; INT_OF_NUM_POW;INT_POW_ODD_NEG1;INT_POW_EVEN_NEG1; INT_ADD_NEG;INT_ADD_NEGv2 (* ; ARITH *) ]) ;; add_test("FLOAT_CONV1", let f z = let (x,y) = dest_eq z in let (u,v) = dest_thm (FLOAT_CONV x) in (u=[]) && (z = v) in f `float (int_of_num 3) (int_of_num 0) = float (int_of_num 3) (int_of_num 0)` && f `float (int_of_num 3) (int_of_num 3) = float (int_of_num 3) (int_of_num 3)` && f `float (int_of_num 3) (int_of_num 0) +. (float (int_of_num 4) (int_of_num 0)) = (float (int_of_num 7) (int_of_num 0))` && f `float (int_of_num 1 + (int_of_num 3)) (int_of_num 4) = float (int_of_num 4) (int_of_num 4)` && f `float (int_of_num 3 - (int_of_num 7)) (int_of_num 0) = float (--:(int_of_num 4)) (int_of_num 0)` && f `float (int_of_num 3) (int_of_num 4) *. (float (--: (int_of_num 2)) (int_of_num 3)) = float (--: (int_of_num 6)) (int_of_num 7)` && f `real_neg (float (--: (int_of_num 3)) (int_of_num 0)) = float (int_of_num 3) (int_of_num 0)` );; (* ------------------------------------------------------------------ *) (* Operations on interval_eps. Preliminary stuff to deal with *) (* chains of inequalities. *) (* ------------------------------------------------------------------ *)
let REAL_ADD_LE_SUBST_RHS = 
prove( `!a b c P. ((a <=. ((P b)) /\ (!x. (P x) = x + (P (&. 0))) /\ (b <=. c)) ==> (a <=. (P c)))`,
(((REPEAT GEN_TAC))) THEN (((REPEAT (TAUT_TAC `(b ==> a==> c) ==> (a /\ b ==> c)`)))) THEN (((REPEAT DISCH_TAC))) THEN ((((H_RULER(ONCE_REWRITE_RULE))[HYP_INT 1] (HYP_INT 0)))) THEN ((((ASM ONCE_REWRITE_TAC)[]))) THEN ((((ASM MESON_TAC)[REAL_LE_RADD;REAL_LE_TRANS]))));;
let REAL_ADD_LE_SUBST_LHS = 
prove( `!a b c P. (((P(a) <=. b /\ (!x. (P x) = x + (P (&. 0))) /\ (c <=. a))) ==> ((P c) <=. b))`,
(REP_GEN_TAC) THEN (DISCH_ALL_TAC) THEN (((H_RULER(ONCE_REWRITE_RULE)) [HYP_INT 1] (HYP_INT 0))) THEN (((ASM ONCE_REWRITE_TAC)[])) THEN (((ASM MESON_TAC)[REAL_LE_RADD;REAL_LE_TRANS])));;
(* let rec SPECL = function [] -> I | (a::b) -> fun thm ->(SPECL b (SPEC a thm));; *) (* input: rel: b <=. c thm: a <=. P(b). output: a <=. P(c). condition: REAL_ARITH must be able to prove !x. P(x) = x+. P(&.0). condition: the term `a` must appear exactly once the lhs of the thm. *) let IWRITE_REAL_LE_RHS rel thm = let bvar = genvar `:real` in let (bt,_) = dest_binop `real_le` (concl rel) in let sub = SUBS_CONV[ASSUME (mk_eq(bt,bvar))] in let rule = (fun th -> EQ_MP (sub (concl th)) th) in let (subrel,subthm) = (rule rel,rule thm) in let (a,p) = dest_binop `real_le` (concl subthm) in let (_,c) = dest_binop `real_le` (concl subrel) in let pfn = mk_abs (bvar,p) in let imp_th = BETA_RULE (SPECL [a;bvar;c;pfn] REAL_ADD_LE_SUBST_RHS) in let ppart = REAL_ARITH (fst(dest_conj(snd(dest_conj(fst(dest_imp(concl imp_th))))))) in let prethm = MATCH_MP imp_th (CONJ subthm (CONJ ppart subrel)) in let prethm2 = SPEC bt (GEN bvar (DISCH (find (fun x -> try(bvar = rhs x) with Failure _ -> false) (hyp prethm)) prethm)) in MATCH_MP prethm2 (REFL bt);; (* input: rel: c <=. a thm: P a <=. b output: P c <=. b condition: REAL_ARITH must be able to prove !x. P(x) = x+. P(&.0). condition: the term `a` must appear exactly once the lhs of the thm. *) let IWRITE_REAL_LE_LHS rel thm = let avar = genvar `:real` in let (_,at) = dest_binop `real_le` (concl rel) in let sub = SUBS_CONV[ASSUME (mk_eq(at,avar))] in let rule = (fun th -> EQ_MP (sub (concl th)) th) in let (subrel,subthm) = (rule rel,rule thm) in let (p,b) = dest_binop `real_le` (concl subthm) in let (c,_) = dest_binop `real_le` (concl subrel) in let pfn = mk_abs (avar,p) in let imp_th = BETA_RULE (SPECL [avar;b;c;pfn] REAL_ADD_LE_SUBST_LHS) in let ppart = REAL_ARITH (fst(dest_conj(snd(dest_conj(fst(dest_imp(concl imp_th))))))) in let prethm = MATCH_MP imp_th (CONJ subthm (CONJ ppart subrel)) in let prethm2 = SPEC at (GEN avar (DISCH (find (fun x -> try(avar = rhs x) with Failure _ -> false) (hyp prethm)) prethm)) in MATCH_MP prethm2 (REFL at);; (* ------------------------------------------------------------------ *) (* INTERVAL ADD, NEG, SUBTRACT *) (* ------------------------------------------------------------------ *)
let INTERVAL_ADD = 
prove( `!x f ex y g ey. interval_eps x f ex /\ interval_eps y g ey ==> interval_eps (x +. y) (f +. g) (ex +. ey)`,
EVERY[ REPEAT GEN_TAC; TAUT_TAC `(A==>B==>C)==>(A/\ B ==> C)`; REWRITE_TAC[interval_eps]; REWRITE_TAC[prove(`(x+.y) -. (f+.g) = (x-.f) +. (y-.g)`,REAL_ARITH_TAC)]; ABBREV_TAC `a = x-.f`; ABBREV_TAC `b = y-.g`; ASSUME_TAC (SPEC `b:real` (SPEC `a:real` ABS_TRIANGLE)); UNDISCH_EL_TAC 0; ABBREV_TAC `a':real = abs a`; ABBREV_TAC `b':real = abs b`; REPEAT DISCH_TAC; (H_VAL2(IWRITE_REAL_LE_RHS)) (HYP_INT 0) (HYP_INT 2); (H_VAL2(IWRITE_REAL_LE_RHS)) (HYP_INT 2) (HYP_INT 0); ASM_REWRITE_TAC[]]);;
let INTERVAL_NEG = 
prove( `!x f ex. interval_eps x f ex = interval_eps (real_neg x) (real_neg f) ex`,
(REWRITE_TAC[interval_eps;REAL_ABS_NEG;REAL_ARITH `!x y. -- x -. (-- y) = real_neg (x -. y)`]));;
let INTERVAL_NEG2 = 
prove( `!x f ex. interval_eps (real_neg x) f ex = interval_eps x (real_neg f) ex`,
(REWRITE_TAC[interval_eps;REAL_ABS_NEG;REAL_ARITH `!x y. -- x -. y = real_neg (x -. (real_neg y))`]));;
let INTERVAL_SUB = 
prove( `!x f ex y g ey. interval_eps x f ex /\ interval_eps y g ey ==> interval_eps (x -. y) (f -. g) (ex +. ey)`,
((REWRITE_TAC[real_sub])) THEN (DISCH_ALL_TAC) THEN (((H_RULER (ONCE_REWRITE_RULE))[THM(INTERVAL_NEG)] (HYP_INT 1))) THEN (((ASM MESON_TAC)[INTERVAL_ADD])));;
(* ------------------------------------------------------------------ *) (* INTERVAL MULTIPLICATION *) (* ------------------------------------------------------------------ *) (* renamed from REAL_LE_ABS_RMUL *)
let REAL_PROP_LE_RABS = 
prove( `!x y z. (y <=. z) ==> ( y * (abs x) <=. z *(abs x))`,
(SIMP_TAC[REAL_LE_RMUL_IMP;ABS_POS]));;
let REAL_LE_ABS_MUL = 
prove( `!x y z w. (( x <=. y) /\ (abs z <=. w)) ==> (x*.w <=. y*.w) `,
(DISCH_ALL_TAC) THEN ((ASSUME_TAC (REAL_ARITH `abs z <=. w ==> (&.0) <=. w`))) THEN (((ASM MESON_TAC)[REAL_LE_RMUL_IMP])));;
let ABS_MUL = REAL_ABS_MUL;; let INTERVAL_ABS = REAL_ARITH `!(x:real) z d. (x - d) <= z /\ z <= (x + d) <=> abs(z - x) <= d`;;
let INTERVAL_MUL = 
prove( `!x f ex y g ey. (interval_eps x f ex) /\ (interval_eps y g ey) ==> (interval_eps (x *. y) (f *. g) (abs(f)*.ey +. ex*. abs(g) +. ex*.ey))`,
(REP_GEN_TAC) THEN ((REWRITE_TAC[interval_eps])) THEN ((REWRITE_TAC[REAL_ARITH `(x*. y -. f*. g) = (f *.(y -. g) +. (x -. f)*.g +. (x-.f)*.(y-. g))`])) THEN (DISCH_ALL_TAC) THEN ((ASSUME_TAC (SPECL [`f*.(y-g)`;`(x-f)*g +. (x-f)*.(y-g)`] ABS_TRIANGLE))) THEN ((ASSUME_TAC (SPECL [`(x-f)*.g`;`(x-f)*.(y-g)`] ABS_TRIANGLE))) THEN (((H_VAL2(IWRITE_REAL_LE_RHS)) (HYP_INT 0) (HYP_INT 1))) THEN ((H_REWRITE_RULE [THM ABS_MUL] (HYP_INT 0))) THEN ((H_MATCH_MP (THM (SPECL [`g:real`; `abs (x -. f)`; `ex:real`] REAL_PROP_LE_RABS)) (HYP_INT 4))) THEN (((H_VAL2(IWRITE_REAL_LE_RHS)) (HYP_INT 0) (HYP_INT 1))) THEN ((H_MATCH_MP (THM (SPECL [`f:real`; `abs (y -. g)`; `ey:real`] REAL_PROP_LE_LABS)) (HYP_INT 7))) THEN (((H_VAL2 (IWRITE_REAL_LE_RHS)) (HYP_INT 0) (HYP_INT 1))) THEN ((H_MATCH_MP (THM (SPECL [`x-.f`; `abs (y -. g)`; `ey:real`] REAL_PROP_LE_LABS)) (HYP_INT 9))) THEN (((H_VAL2(IWRITE_REAL_LE_RHS)) (HYP_INT 0) (HYP_INT 1))) THEN ((ASSUME_TAC (SPECL [`abs(x-.f)`;`ex:real`;`y-.g`;`ey:real`] REAL_LE_ABS_MUL))) THEN ((H_CONJ (HYP_INT 11) (HYP_INT 12))) THEN ((H_MATCH_MP (HYP_INT 1) (HYP_INT 0))) THEN (((H_VAL2(IWRITE_REAL_LE_RHS)) (HYP_INT 0) (HYP_INT 3))) THEN ((POP_ASSUM ACCEPT_TAC)));;
(* ------------------------------------------------------------------ *) (* INTERVAL BASIC OPERATIONS *) (* ------------------------------------------------------------------ *)
let INTERVAL_NUM = 
prove( `!n. (interval_eps(&.n) (float(int_of_num n) (int_of_num 0)) (float (int_of_num 0) (int_of_num 0)))`,
let INTERVAL_CENTER = 
prove( `!x f ex g. (interval_eps x f ex) ==> (interval_eps x g (abs(f-g)+.ex))`,
((REWRITE_TAC[interval_eps])) THEN (DISCH_ALL_TAC) THEN ((ASSUME_TAC (REAL_ARITH `abs(x -. g) <=. abs(f-.g) +. abs(x -. f)`))) THEN ((H_VAL2 IWRITE_REAL_LE_RHS (HYP_INT 1) (HYP_INT 0))) THEN ((ASM_REWRITE_TAC[])));;
let INTERVAL_WIDTH = 
prove( `!x f ex ex'. (ex <=. ex') ==> (interval_eps x f ex) ==> (interval_eps x f ex')`,
((REWRITE_TAC[interval_eps])) THEN (DISCH_ALL_TAC) THEN ((H_VAL2 IWRITE_REAL_LE_RHS (HYP_INT 1) (HYP_INT 0))) THEN ((ASM_REWRITE_TAC[])));;
let INTERVAL_MAX = 
prove( `!x f ex. interval_eps x f ex ==> (x <=. f+.ex)`,
(REWRITE_TAC[interval_eps]) THEN REAL_ARITH_TAC);;
let INTERVAL_MIN = 
prove( `!x f ex. interval_eps x f ex ==> (f-. ex <=. x)`,
(REWRITE_TAC[interval_eps]) THEN REAL_ARITH_TAC);;
let INTERVAL_ABS_MIN = 
prove( `!x f ex. interval_eps x f ex ==> (abs(f)-. ex <=. abs(x))`,
(REWRITE_TAC[interval_eps] THEN REAL_ARITH_TAC) );;
let INTERVAL_ABS_MAX = 
prove( `!x f ex. interval_eps x f ex ==> (abs(x) <=. abs(f)+. ex)`,
(REWRITE_TAC[interval_eps] THEN REAL_ARITH_TAC) );;
let REAL_RINV_2 = 
prove( `&.2 *. (inv (&.2 )) = &. 1`,
EVERY[ MATCH_MP_TAC REAL_MUL_RINV; REAL_ARITH_TAC]);;
let INTERVAL_MK = 
prove( `let half = float(int_of_num 1)(--:(int_of_num 1)) in !x xmin xmax. ((xmin <=. x) /\ (x <=. xmax)) ==> interval_eps x ((xmin+.xmax)*.half) ((xmax-.xmin)*.half)`,
EVERY[ REWRITE_TAC[LET_DEF;LET_END_DEF]; DISCH_ALL_TAC; REWRITE_TAC[interval_eps;float;TWOPOW_NEG;INT_NUM_REAL;REAL_POW_1;REAL_MUL_LID]; REWRITE_TAC[GSYM INTERVAL_ABS]; CONJ_TAC ] THENL[ EVERY[ REWRITE_TAC[GSYM REAL_SUB_RDISTRIB]; REWRITE_TAC[REAL_ARITH `(b+.a)-.(a-.b)=b*.(&.2)`;GSYM REAL_MUL_ASSOC]; ASM_REWRITE_TAC[REAL_RINV_2;REAL_MUL_RID] ]; EVERY[ REWRITE_TAC[GSYM REAL_ADD_RDISTRIB]; REWRITE_TAC[REAL_ARITH `(b+.a)+. a -. b=a*.(&.2)`;GSYM REAL_MUL_ASSOC]; ASM_REWRITE_TAC[REAL_RINV_2;REAL_MUL_RID] ] ]);;
let INTERVAL_EPS_POS = 
prove(`!x f ex. (interval_eps x f ex) ==> (&.0 <=. ex)`,
EVERY[ REWRITE_TAC[interval_eps]; REPEAT (GEN_TAC); DISCH_THEN(fun x -> (MP_TAC (CONJ (SPEC `x-.f` REAL_ABS_POS) x))); MATCH_ACCEPT_TAC REAL_LE_TRANS]);;
let INTERVAL_EPS_0 = 
prove( `!x f n. (interval_eps x f (float (int_of_num 0) n)) ==> (x = f)`,
EVERY[ REWRITE_TAC[interval_eps;float;int_of_num_th;REAL_MUL_LZERO]; REAL_ARITH_TAC]);;
let REAL_EQ_RCANCEL_IMP' = 
prove(`!x y z.(x * z = y * z) ==> (~(z = &0) ==> (x=y))`,
MESON_TAC[REAL_EQ_RCANCEL_IMP]);;
(* renamed from REAL_ABS_POS *)
let REAL_MK_POS_ABS_' = 
prove (`!x. (~(x=(&.0))) ==> (&.0 < abs(x))`,
(* ------------------------------------------------------------------ *) (* INTERVAL DIVIDE *) (* ------------------------------------------------------------------ *)
let INTERVAL_DIV = 
prove(`!x f ex y g ey h ez. (((interval_eps x f ex) /\ (interval_eps y g ey) /\ (ey <. (abs g)) /\ ((ex +. (abs (f -. (h*.g))) +. (abs h)*. ey) <=. (ez*.((abs g) -. ey)))) ==> (interval_eps (x / y) h ez))`,
let lemma1 = prove( `&.0 < u /\ real_abs  z <=. e*. u ==> (&.0) <=. e`,
  EVERY[
    DISCH_ALL_TAC;
    ASSUME_TAC (SPEC `z:real` REAL_MK_NN_ABS);
    H_MATCH_MP (THM REAL_LE_TRANS) (H_RULE2 CONJ (HYP_INT 0) (HYP_INT 2));
    H_MATCH_MP (THM REAL_PROP_NN_RCANCEL) (H_RULE2 CONJ (HYP_INT 2) (HYP_INT 0));
    ASM_REWRITE_TAC[]
  ]) in
EVERY[
  DISCH_ALL_TAC;
  SUBGOAL_MP_TAC `~(y= (&.0))`
  THENL[
    EVERY[
      UNDISCH_LIST[1;2];
      REWRITE_TAC[interval_eps];
      REAL_ARITH_TAC
    ];
    EVERY[
      REWRITE_TAC[interval_eps];
      DISCH_TAC THEN (H I (HYP_INT 0)) THEN (UNDISCH_EL_TAC 0);
      DISCH_THEN (fun th -> (MP_TAC(MATCH_MP REAL_MK_POS_ABS_' th)));
      MATCH_MP_TAC REAL_MUL_RTIMES_LE;
      REWRITE_TAC[GSYM ABS_MUL;REAL_SUB_RDISTRIB;real_div;GSYM REAL_MUL_ASSOC];
      ASM_SIMP_TAC[REAL_MUL_LINV;REAL_MUL_RID];
      H (REWRITE_RULE[interval_eps]) (HYP_INT 1);
      H (REWRITE_RULE[interval_eps]) (HYP_INT 3);
      H (MATCH_MP INTERVAL_ABS_MIN) (HYP_INT 4);
      POPL_TAC[3;4;5];
      H_VAL2 (IWRITE_REAL_LE_LHS) (HYP_INT 2) (HYP_INT 4);
      H (REWRITE_RULE[ REAL_ADD_ASSOC]) (HYP_INT 0);
      H_VAL2 (IWRITE_REAL_LE_LHS) (THM (SPEC `f-. h*g` (SPEC `x-.f` ABS_TRIANGLE))) (HYP_INT 0);
      H (ONCE_REWRITE_RULE[REAL_ABS_SUB]) (HYP_INT 4);
      H (MATCH_MP (SPEC `h:real` REAL_PROP_LE_LABS)) (HYP_INT 0);
      H (REWRITE_RULE[GSYM ABS_MUL]) (HYP_INT 0);
      H_VAL2 (IWRITE_REAL_LE_LHS) (HYP_INT 0) (HYP_INT 3);
      H_VAL2 (IWRITE_REAL_LE_LHS) (THM (SPEC `h*.(g-.y)` (SPEC`(x-.f)+(f-. h*g)`  ABS_TRIANGLE))) (HYP_INT 0);
      POPL_TAC[1;2;3;4;5;6;7;9;10;12];
      H (ONCE_REWRITE_RULE[REAL_ARITH `((x-.f) +. (f -. h*. g)) +. h*.(g-. y) = x -. h*. y `]) (HYP_INT 0);
      ABBREV_TAC `z = x -. h*.y`;
      H (ONCE_REWRITE_RULE[GSYM REAL_SUB_LT]) (HYP_INT 4);
      ABBREV_TAC `u = abs(g) -. ey`;
      POPL_TAC[0;2;4;6];
      H (MATCH_MP lemma1 ) (H_RULE2 CONJ (HYP_INT 0) (HYP_INT 1));
      H (MATCH_MP REAL_PROP_LE_LMUL) (H_RULE2 CONJ (HYP_INT 0) (HYP_INT 3));
      H (MATCH_MP REAL_LE_TRANS) (H_RULE2 CONJ (HYP_INT 3) (HYP_INT 0));
      ASM_REWRITE_TAC[]
  ];
  ]]);;
(* ------------------------------------------------------------------ *) (* INTERVAL ABS VALUE *) (* ------------------------------------------------------------------ *)
let INTERVAL_ABSV = 
prove(`!x f ex. interval_eps x f ex ==> (interval_eps (abs x) (abs f) ex)`,
EVERY[ REWRITE_TAC[interval_eps]; DISCH_ALL_TAC; ASSUME_TAC (SPECL [`x:real`;`f:real`] REAL_ABS_SUB_ABS); ASM_MESON_TAC[REAL_LE_TRANS] ]);;
(* 7 minutes *) (* ------------------------------------------------------------------ *) (* INTERVAL ATN VALUE *) (* ------------------------------------------------------------------ *)
let taylor_atn = 
prove_by_refinement( `!n x v eps2 eps. interval_eps (&16 * sum (0..n) (halfatn4_co x)) v eps2 /\ float (int_of_num 1) ( -- int_of_num (6 * n + 5)) + eps2 <= eps ==> interval_eps (atn x) v eps`,
(* {{{ proof *) [ REWRITE_TAC[interval_eps;float;TWOPOW_NEG;int_of_num_th;REAL_ARITH `&1 * x = x`]; REPEAT STRIP_TAC; MATCH_MP_TAC real_taylor_atn_approx; EXISTS_TAC `n:num`; EXISTS_TAC `inv(&2 pow (6 * n +5))`; EXISTS_TAC `eps2:real`; ASM_REWRITE_TAC[]; REAL_ARITH_TAC; ]);;
(* }}} *)
let pi_atn = 
prove_by_refinement( `pi = &4 * atn (&1)`,
(* {{{ proof *) [ REWRITE_TAC[ATN_1]; REAL_ARITH_TAC; ]);;
(* }}} *) (* ------------------------------------------------------------------ *) (* INTERVAL SQRT *) (* This requires some preliminaries. Extend sqrt by 0 on negatives *) (* ------------------------------------------------------------------ *)
let ssqrt = new_definition `ssqrt x = if (x <. (&.0)) then (&.0) else sqrt x`;;
(*2m*)
let ssqrt_sqrt = 
prove_by_refinement( `!x u f. &0 < x /\ interval_eps (ssqrt x) u f ==> interval_eps (sqrt x) u f`,
(* {{{ proof *) [ REWRITE_TAC[ssqrt]; MESON_TAC[REAL_ARITH `&0 < x ==> ~(x < &0)`]; ]);;
(* }}} *) let LET_TAC = REWRITE_TAC[LET_DEF;LET_END_DEF];;
let REAL_SSQRT_NEG = 
prove(`!x. (x <. (&.0)) ==> (ssqrt x = (&.0))`,
EVERY[ DISCH_ALL_TAC; REWRITE_TAC[ssqrt]; COND_CASES_TAC THENL[ ACCEPT_TAC (REFL `&.0`); ASM_MESON_TAC[] ] ]);;
(* 5 min*)
let REAL_SSQRT_NN = 
prove(`!x. (&.0) <=. x ==> (ssqrt x = (sqrt x))`,
EVERY[ DISCH_ALL_TAC; REWRITE_TAC[ssqrt]; COND_CASES_TAC THENL[ ASM_MESON_TAC[real_lt]; ACCEPT_TAC (REFL `sqrt x`) ] ]);;
(* 12 min, mostly spent loading *index-shell* *) (*17 minutes*)
let REAL_MK_NN_SSQRT = 
prove(`!x. (&.0) <=. (ssqrt x)`,
EVERY[ GEN_TAC; DISJ_CASES_TAC (SPECL[`x:real`;`&.0`] REAL_LTE_TOTAL) THENL[ POP_ASSUM (fun th -> MP_TAC(MATCH_MP (REAL_SSQRT_NEG) th)) THEN MESON_TAC[REAL_LE_REFL]; POP_ASSUM (fun th -> ASSUME_TAC(CONJ th (MATCH_MP (REAL_SSQRT_NN) th))) THEN ASM_MESON_TAC[REAL_PROP_NN_SQRT] ] ]);;
let REAL_SV_SSQRT_0  = 
prove(`!x. ssqrt (&.0) = (&.0)`,
EVERY[ GEN_TAC; MP_TAC (SPEC `&.0` REAL_LE_REFL); DISCH_THEN(fun th -> REWRITE_TAC[MATCH_MP REAL_SSQRT_NN th]); ACCEPT_TAC REAL_SV_SQRT_0 ]);;
(* 6 minutes *)
let REAL_SV_SSQRT_n  = 
prove_by_refinement(`!n. ssqrt (&n) = sqrt(&n)`,
[ GEN_TAC; MATCH_MP_TAC REAL_SSQRT_NN; REAL_ARITH_TAC; ]);;
let REAL_SSQRT_EQ_0 = 
prove(`!(x:real). (ssqrt(x) = (&.0)) ==> (x <=. (&.0))`,
EVERY[ GEN_TAC; DISJ_CASES_TAC (SPECL[`x:real`;`&.0`] REAL_LTE_TOTAL) THENL[ ASM_MESON_TAC[REAL_LT_IMP_LE]; ASM_SIMP_TAC[REAL_SSQRT_NN] THEN ASM_MESON_TAC[SQRT_EQ_0;REAL_EQ_IMP_LE] ] ]);;
(* 15 minutes *)
let REAL_SSQRT_MONO = 
prove(`!x. (x<=. y) ==> (ssqrt x <=. (ssqrt y))`,
EVERY[ GEN_TAC; DISJ_CASES_TAC (SPECL[`x:real`;`&.0`] REAL_LTE_TOTAL) THENL[ ASM_MESON_TAC[REAL_SSQRT_NEG;REAL_MK_NN_SSQRT]; ASM_MESON_TAC[REAL_LE_TRANS;REAL_SSQRT_NN;REAL_PROP_LE_SQRT]; ] ]);;
(* 5 minutes *)
let REAL_SSQRT_CHAR = 
prove(`!x t. (&.0 <=. t /\ (t*t = x)) ==> (t = (ssqrt x))`,
EVERY[ DISCH_ALL_TAC; H_ASSUME_TAC (H_RULE_LIST REWRITE_RULE[HYP_INT 1] (THM (SPEC `t:real` REAL_MK_NN_SQUARE))); ASM_MESON_TAC[REAL_SSQRT_NN;SQRT_MUL;POW_2_SQRT_ABS;REAL_POW_2;REAL_ABS_REFL]; ]);;
(* 13 minutes *)
let REAL_SSQRT_SQUARE = 
prove(`!x. (&.0 <=. x) ==> ((ssqrt x)*.(ssqrt x) = x)`,
(* 7min *)
let REAL_SSQRT_SQUARE' = 
prove(`!x. (&.0<=. x) ==> (ssqrt (x*.x) = x)`,
DISCH_ALL_TAC THEN REWRITE_TAC[(MATCH_MP REAL_SSQRT_NN (SPEC `x:real` REAL_MK_NN_SQUARE))] THEN ASM_SIMP_TAC[SQRT_MUL;GSYM POW_2;SQRT_POW_2]);;
(*20min*) (* an alternate proof appears in RCS *)
let INTERVAL_SSQRT = 
prove(`!x f ex u ey ez v. (interval_eps x f ex) /\ (interval_eps (u*.u) f ey) /\ (ex +.ey <=. ez*.(v+.u)) /\ (v*.v <=. f-.ex) /\ (&.0 <. u) /\ (&.0 <=. v) ==> (interval_eps (ssqrt x) u ez)`,
EVERY[ DISCH_ALL_TAC; H (MATCH_MP REAL_LE_TRANS) (H_RULE2 CONJ (THM (SPEC `v:real` REAL_MK_NN_SQUARE)) (HYP_INT 3)); H (MATCH_MP (INTERVAL_MIN)) (HYP_INT 1); H (MATCH_MP REAL_LE_TRANS) (H_RULE2 CONJ (HYP_INT 1) (HYP_INT 0)); H (MATCH_MP INTERVAL_EPS_POS) (HYP_INT 3); H (MATCH_MP INTERVAL_EPS_POS) (HYP_INT 5); H (MATCH_MP REAL_PROP_NN_ADD2) (H_RULE2 CONJ (HYP_INT 1) (HYP_INT 0)); H (MATCH_MP REAL_PROP_POS_LADD) (H_RULE2 CONJ (HYP_INT 11) (HYP_INT 10)); H (MATCH_MP REAL_PROP_POS_LADD) (H_RULE2 CONJ (THM (SPEC `x:real` REAL_MK_NN_SSQRT)) (HYP_INT 11)); H (MATCH_MP REAL_PROP_POS_INV) (HYP_INT 0); ASSUME_TAC (REAL_ARITH `(ssqrt x -. u) = (ssqrt x -. u)*.(&.1)`); H (MATCH_MP REAL_MK_NZ_POS) (HYP_INT 2); H (MATCH_MP REAL_MUL_RINV) (HYP_INT 0); H_REWRITE_RULE[(H_RULE GSYM) (HYP_INT 0)] (HYP_INT 2); POPL_TAC[1;2;3]; H (REWRITE_RULE[REAL_MUL_ASSOC]) (HYP_INT 0); H (REWRITE_RULE[ONCE_REWRITE_RULE[REAL_MUL_SYM] REAL_DIFFSQ]) (HYP_INT 0); POPL_TAC[1;2]; H_SIMP_RULE[HYP_INT 7;THM REAL_SSQRT_SQUARE] (HYP_INT 0); ASSUME_TAC (REAL_ARITH `abs(x -. u*.u) <=. abs(x -. f) + abs(f-. u*.u)`); H (REWRITE_RULE[interval_eps]) (HYP_INT 12); H (ONCE_REWRITE_RULE[interval_eps]) (HYP_INT 14); H (ONCE_REWRITE_RULE[REAL_ABS_SUB]) (HYP_INT 0); POPL_TAC[1;5;15;16]; H (MATCH_MP REAL_LE_ADD2) (H_RULE2 CONJ (HYP_INT 1) (HYP_INT 0)); H (MATCH_MP REAL_LE_TRANS) (H_RULE2 CONJ (HYP_INT 3) (HYP_INT 0)); POPL_TAC[1;2;3;4]; H (AP_TERM `real_abs `) (HYP_INT 1); H (REWRITE_RULE[ABS_MUL]) (HYP_INT 0); H (MATCH_MP REAL_LT_IMP_LE) (HYP_INT 4); H (REWRITE_RULE[GSYM REAL_ABS_REFL]) (HYP_INT 0); H_REWRITE_RULE [HYP_INT 0] (HYP_INT 2); H (MATCH_MP REAL_LE_RMUL) (H_RULE2 CONJ (HYP_INT 5) (HYP_INT 2)); H_REWRITE_RULE [H_RULE GSYM (HYP_INT 1)] (HYP_INT 0); POPL_TAC[1;2;3;5;6;7;8]; H (MATCH_MP REAL_LE_TRANS) (H_RULE2 CONJ (HYP_INT 12) (HYP_INT 9)); H (MATCH_MP REAL_SSQRT_MONO) (HYP_INT 0); H (MATCH_MP REAL_SSQRT_SQUARE') (HYP_INT 16); H_REWRITE_RULE [HYP_INT 0] (HYP_INT 1); H (ONCE_REWRITE_RULE[GSYM (SPECL[`v:real`;`ssqrt x`;`u:real`] REAL_LE_RADD)]) (HYP_INT 0); H (MATCH_MP REAL_LE_INV2) (H_RULE2 CONJ (HYP_INT 9) (HYP_INT 0)); POPL_TAC[1;2;3;4;5;7;8;9;12;13]; H (MATCH_MP REAL_LE_LMUL) (H_RULE2 CONJ (HYP_INT 3) (HYP_INT 0)); H (MATCH_MP REAL_LE_TRANS) (H_RULE2 CONJ (HYP_INT 2) (HYP_INT 0)); H (MATCH_MP REAL_PROP_POS_INV) (HYP_INT 4); H (MATCH_MP REAL_LT_IMP_LE) (HYP_INT 0); H (MATCH_MP REAL_LE_RMUL) (H_RULE2 CONJ (HYP_INT 11) (HYP_INT 0)); H (REWRITE_RULE[GSYM REAL_MUL_ASSOC]) (HYP_INT 0); H (MATCH_MP REAL_MK_NZ_POS) (HYP_INT 8); H (MATCH_MP REAL_MUL_RINV) (HYP_INT 0); H_REWRITE_RULE[HYP_INT 0; THM REAL_MUL_RID] (HYP_INT 2); H (MATCH_MP REAL_LE_TRANS) (H_RULE2 CONJ (HYP_INT 7) (HYP_INT 0)); ASM_REWRITE_TAC[interval_eps] ]);;
test();; (* conversion for interval_eps *) (* ------------------------------------------------------------------ *) (* Take a term x of type real. Convert to a thm of the form *) (* interval_eps x f eps *) (* *) (* ------------------------------------------------------------------ *) let add_test,test = new_test_suite();; (* ------------------------------------------------------------------ *) (* num_exponent Take the absolute value of input. Write it as a*2^k, where 1 <= a < 2, return k. Except: num_exponent (Int 0) is -1. *) let (num_exponent:Num.num -> Num.num) = fun a -> let afloat = float_of_num (abs_num a) in Int ((snd (frexp afloat)) - 1);; (*test*)let f (u,v) = ((num_exponent u) =(Int v)) in add_test("num_exponenwt", forall f [Int 1,0; Int 65,6; Int (-65),6; Int 0,-1; (Int 3)//(Int 4),-1]);; (* ------------------------------------------------------------------ *) let dest_unary op tm = try let xop,r = dest_comb tm in if xop = op then r else fail() with Failure _ -> failwith "dest_unary";; (* ------------------------------------------------------------------ *) (* finds a nearby (outward-rounded) Int with only prec_b significant bits *) let (round_outward: int -> Num.num -> Num.num) = fun prec_b a -> let b = abs_num a in let sign = if (a =/ b) then I else minus_num in let throw_bits = Num.max_num (Int 0) ((num_exponent b)-/ (Int prec_b)) in let twoexp = power_num (Int 2) throw_bits in (sign (ceiling_num (b // twoexp)))*/twoexp;; let (round_inward: int-> Num.num -> Num.num) = fun prec_b a -> let b = abs_num a in let sign = if (a=/b) then I else minus_num in let throw_bits = Num.max_num (Int 0) ((num_exponent b)-/ (Int prec_b)) in let twoexp = power_num (Int 2) throw_bits in (sign (floor_num (b // twoexp)))*/twoexp;; let round_rat bprec n = let b = abs_num n in let sign = if (b =/ n) then I else minus_num in let powt = ((Int 2) **/ (Int bprec)) in sign ((round_outward bprec (Num.ceiling_num (b */ powt)))//powt);; let round_inward_rat bprec n = let b = abs_num n in let sign = if (b =/ n) then I else minus_num in let powt = ((Int 2) **/ (Int bprec)) in sign ((round_inward bprec (Num.floor_num (b */ powt)))//powt);; let (round_outward_float: int -> float -> Num.num) = fun bprec f -> if (f=0.0) then (Int 0) else begin let b = abs_float f in let sign = if (f >= 0.0) then I else minus_num in let (x,n) = frexp b in let u = int_of_float( ceil (ldexp x bprec)) in sign ((Int u)*/ ((Int 2) **/ (Int (n - bprec)))) end;; let (round_inward_float: int -> float -> Num.num) = fun bprec f -> if (f=0.0) then (Int 0) else begin (* avoid overflow on 30 bit integers *) let bprec = if (bprec > 25) then 25 else bprec in let b = abs_float f in let sign = if (f >= 0.0) then I else minus_num in let (x,n) = frexp b in let u = int_of_float( floor (ldexp x bprec)) in sign ((Int u)*/ ((Int 2) **/ (Int (n - bprec)))) end;; (* ------------------------------------------------------------------ *) (* This doesn't belong here. A general term substitution function *) let SUBST_TERM sublist tm = rhs (concl ((SPECL (map fst sublist)) (GENL (map snd sublist) (REFL tm))));; add_test("SUBST_TERM", SUBST_TERM [(`#1`,`a:real`);(`#2`,`b:real`)] (`a +. b +. c`) = `#1 + #2 + c`);; (* ------------------------------------------------------------------ *) (* take a term of the form `interval_eps x f ex` and clean up the f and ex *) let INTERVAL_CLEAN_CONV:conv = fun interv -> let (ixf,ex) = dest_comb interv in let (ix,f) = dest_comb ixf in let fthm = FLOAT_CONV f in let exthm = FLOAT_CONV ex in let ixfthm = AP_TERM ix fthm in MK_COMB (ixfthm, exthm);; (*test*) add_test("INTERVAL_CLEAN_CONV", let testval = INTERVAL_CLEAN_CONV `interval_eps ((&.1) +. (&.1)) (float (int_of_num 3) (int_of_num 4) +. (float (int_of_num 2) (--: (int_of_num 3)))) (float (int_of_num 1) (int_of_num 2) *. (float (int_of_num 3) (--: (int_of_num 2))))` in let hypval = hyp testval in let concval = concl testval in (length hypval = 0) && concval = `interval_eps (&1 + &1) (float (int_of_num 3) (int_of_num 4) + float (int_of_num 2) (--: (int_of_num 3))) (float (int_of_num 1) (int_of_num 2) * float (int_of_num 3) (--: (int_of_num 2))) = interval_eps (&1 + &1) (float (int_of_num 386) (--: (int_of_num 3))) (float (int_of_num 3) (int_of_num 0))` );; (* ------------------------------------------------------------------ *) (* GENERAL lemmas *) (* ------------------------------------------------------------------ *) (* verifies statement of the form `float a b = float a' b'` *) let REAL_INV_POS = REAL_LT_INV;;
let TWOPOW_MK_POS = 
prove( `!a. (&.0 <. ( twopow a))`,
EVERY[ GEN_TAC; CHOOSE_TAC (SPEC `a:int` INT_REP2); POP_ASSUM DISJ_CASES_TAC; ASM_REWRITE_TAC[TWOPOW_POS;TWOPOW_NEG]; TRY (MATCH_MP_TAC REAL_INV_POS); MATCH_MP_TAC REAL_POW_LT ; REAL_ARITH_TAC; ]);;
let TWOPOW_NZ = 
prove( `!a. ~(twopow a = (&.0))`,
GEN_TAC THEN ACCEPT_TAC (MATCH_MP REAL_MK_NZ_POS (SPEC `a:int` TWOPOW_MK_POS)));;
let FLOAT_ZERO = 
prove( `!a b. (float a b = (&.0)) <=> (a = (int_of_num 0))`,
EVERY[ REWRITE_TAC[float;REAL_ENTIRE;INT_OF_NUM_DEST]; MESON_TAC[TWOPOW_NZ]; ]);;
let INT_ZERO = 
prove( `!n. ((int_of_num n = (int_of_num 0)) = (n=0))`,
REWRITE_TAC[INT_OF_NUM_EQ]);;
let INT_ZERO_NEG=
prove( `!n. ((--: (int_of_num n) = (int_of_num 0))) <=> (n=0)`,
REWRITE_TAC[INT_NEG_EQ_0;INT_ZERO]);;
let FLOAT_NN = Prove_by_refinement.prove_by_refinement( `!a b. ((&.0) <=. (float a b)) <=> (int_of_num 0 <=: a)`, [ REWRITE_TAC[float;INT_OF_NUM_DEST]; REP_GEN_TAC; REWRITE_TAC[int_of_num_th;int_le]; MESON_TAC[REAL_PROP_NN_RCANCEL;TWOPOW_MK_POS;REAL_PROP_NN_POS;REAL_PROP_NN_MUL2]; ]);; let INT_NN = INT_POS;;
let INT_NN_NEG = 
prove(`!n. ((int_of_num 0) <=: (--:(int_of_num n))) <=> (n = 0)`,
REWRITE_TAC[INT_NEG_GE0;INT_OF_NUM_LE] THEN ARITH_TAC );;
let FLOAT_POS = 
prove(`!a b. ((&.0) <. (float a b)) <=> (int_of_num 0 <: a)`,
let INT_POS' = 
prove(`!n. (int_of_num 0) <: (int_of_num n) <=> (~(n=0) )`,
REWRITE_TAC[INT_OF_NUM_LT] THEN ARITH_TAC);;
let INT_POS_NEG =
prove(`!n. ((int_of_num 0) <: (--:(int_of_num n))) <=> F`,
REWRITE_TAC[INT_OF_NUM_LT] THEN ARITH_TAC);;
let RAT_LEMMA1_SUB = 
prove(`~(y1 = &0) /\ ~(y2 = &0) ==> ((x1 / y1) - (x2 / y2) = (x1 * y2 - x2 * y1) * inv(y1) * inv(y2))`,
EVERY[REWRITE_TAC[real_div]; REWRITE_TAC[real_sub;GSYM REAL_MUL_LNEG]; REWRITE_TAC[GSYM real_div]; SIMP_TAC[RAT_LEMMA1]; DISCH_TAC; MESON_TAC[real_div]]);;
let INTERVAL_0 = 
prove(`! a f ex. (interval_eps a f ex <=> (&.0 <= (ex - (abs (a -. f)))))`,
MESON_TAC[interval_eps;REAL_SUB_LE]);;
let ABS_NUM_P = 
prove_by_refinement( `!m n. (m <= n) ==> (abs (&. n -. (&. m)) = &.((m-|n) + (n-|m)))`,
[ REPEAT STRIP_TAC; ASM_SIMP_TAC[(ARITH_RULE `(m <= n ==> (m - n = 0) )/\ ( 0 + x = x)`)]; ASM_SIMP_TAC[GSYM REAL_OF_NUM_SUB;REAL_ABS_REFL;]; FIRST_X_ASSUM MP_TAC; REWRITE_TAC[GSYM REAL_OF_NUM_LE]; REAL_ARITH_TAC; ]) ;;
let ABS_NUM = 
prove_by_refinement ( `!m n. abs (&. n -. (&. m)) = &.((m-|n) + (n-|m))`,
[ REPEAT GEN_TAC ; DISJ_CASES_TAC (SPECL [`m:num`;`n:num`] LTE_CASES) ; (* THENL[ *) (* first case *) MATCH_MP_TAC ABS_NUM_P; FIRST_X_ASSUM MP_TAC; ARITH_TAC; (* second case *) ONCE_REWRITE_TAC[REAL_ARITH `abs(x - y) = abs(y - x)`]; ONCE_REWRITE_TAC[ARITH_RULE `(a:num) + b = b + a`]; MATCH_MP_TAC ABS_NUM_P; ASM_REWRITE_TAC[]; ]);;
let INTERVAL_TO_LESS = 
prove( `!a f ex b g ey. ((interval_eps a f ex) /\ (interval_eps b g ey) /\ (&.0 <. (g -. (ey +. ex +. f)))) ==> (a <. b)`,
let lemma1 = REAL_ARITH `!ex ey f g. (&.0 <. (g -. (ey +. ex +. f))) ==> ((f +. ex)<. (g -. ey)) ` in EVERY[ REPEAT GEN_TAC; DISCH_ALL_TAC; H_MATCH_MP (THM lemma1) (HYP "2"); H_MATCH_MP (THM INTERVAL_MAX) (HYP "0"); H_MATCH_MP (THM INTERVAL_MIN) (HYP "1"); LABEL_ALL_TAC; H_MATCH_MP (THM REAL_LET_TRANS) (H_RULE2 CONJ (HYP "4") (HYP "5")); LABEL_ALL_TAC; H_MATCH_MP (THM REAL_LTE_TRANS) (H_RULE2 CONJ (HYP "6") (HYP "3")); ASM_REWRITE_TAC[] ]);;
let ABS_TO_INTERVAL = 
prove( `!c u k. (abs (c - u) <=. k) ==> (!f g ex ey.((interval_eps u f ex) /\ (interval_eps k g ey) ==> (interval_eps c f (g+.ey+.ex))))`,
EVERY[ REWRITE_TAC[interval_eps]; DISCH_ALL_TAC; REPEAT GEN_TAC; DISCH_ALL_TAC; ONCE_REWRITE_TAC [REAL_ARITH `c -. f = (c-. u) + (u-. f)`]; ONCE_REWRITE_TAC [REAL_ADD_ASSOC]; ASSUME_TAC (SPECL [`c-.u`;`u-.f`] ABS_TRIANGLE); IMP_RES_THEN ASSUME_TAC (REAL_ARITH `real_abs (k-.g) <=. ey ==> (k <=. (g +. ey))`); MATCH_MP_TAC (REAL_ARITH `(?a b.((x <=. (a+.b)) /\ (a <=. u) /\ (b <=. v))) ==> (x <=. (u +. v))`); EXISTS_TAC `abs (c-.u)`; EXISTS_TAC `abs(u-.f)`; ASM_REWRITE_TAC[]; ASM_MESON_TAC[REAL_LE_TRANS]; ]);;
(* end of general lemmas *) (* ------------------------------------------------------------------ *) (* ------------------------------------------------------------------ *) (* Cache of computed constants (abs (r - u) <= k) or interval r u eps *) (* ------------------------------------------------------------------ *) let cache_entry_of_ineq bprec ineq = try( let (abst,_) = dest_binop `real_le` (concl ineq) in let (_,cmu) = dest_comb abst in let (r,_) = dest_binop `real_sub` cmu in (r,bprec,ineq)) with Failure _ -> (try( let (r,_,_) = dest_interval (concl ineq) in (r,bprec,ineq)) with Failure _ -> (print_term (concl ineq); failwith "calculated_constants format : abs(r - u) <= k"));; let get_real_constant cache bprec tm = (fun (_,_,r) -> r) (find (fun (t,p,th) -> (t = tm) & (p>= bprec)) cache);; let remove_entry t = filter (fun x -> not(x = t));; (* ------------------------------------------------------------------ *) (* take |- interval_eps r u e to interval_eps r' u e, where th_r gives r <=> r' *) let interval_rule th_r th_int = let (iru,e) = dest_comb (concl th_int) in let (ir,u) = dest_comb iru in let (i,r) = dest_comb ir in let ir_t = AP_TERM i th_r in let iru_t = AP_THM ir_t u in let irue_t = AP_THM iru_t e in let a = try EQ_MP irue_t th_int with Failure m -> (print_thm irue_t; print_thm th_int; failwith m) in a;; (* term of the form '&.n'. Assume error checking done already. *) let INTERVAL_OF_NUM:conv = fun tm -> let n = snd (dest_comb tm) in let th_r = (ARITH_REWRITE_CONV[] n) in let th_int = SPEC (rhs (concl th_r)) INTERVAL_NUM in interval_rule (SYM (AP_TERM `real_of_num` th_r)) th_int;; add_test("INTERVAL_OF_NUM", dest_thm (INTERVAL_OF_NUM `&.3`) = ([], `interval_eps (&3) (float (int_of_num 3) (int_of_num 0)) (float (int_of_num 0) (int_of_num 0))`));; (* ------------------------------------------------------------------ *) (* conversion for float inequalities *) (* ------------------------------------------------------------------ *) let FLOAT_LESS_CONV = FLOAT_CONV THENC ARITH_SIMP_CONV[FLOAT_POS; INT_POS';INT_POS_NEG;FLOAT_NN;INT_NN;INT_NN_NEG;];; let FLOAT_INEQ_CONV t = let th =FLOAT_LESS_CONV t in let (_,rhs) = dest_eq (concl th) in let _ = (rhs = `T`) or raise Insufficient_precision in EQ_MP (SYM th) TRUTH;; add_test("FLOAT_INEQ_CONV", let f c x = dest_thm (c x) = ([],x) in let t1 = `(float (int_of_num 3) (int_of_num 0)) +. (float (int_of_num 4) (int_of_num 0)) <. (float (int_of_num 8) (int_of_num 1))` in ((f FLOAT_INEQ_CONV t1)));; let INTERVAL_TO_LESS_CONV = let rthm = ASSUME `!f g ex ey. (&.0 <. (g -. (ey +. ex +. f)))` in fun thm1 thm2 -> let (a,f,ex) = dest_interval (concl thm1) in let (b,g,ey) = dest_interval (concl thm2) in let rspec = concl (SPECL [f;g;ex;ey] rthm) in let rspec_th = FLOAT_INEQ_CONV rspec in let fthm = CONJ thm1 (CONJ thm2 rspec_th) in MATCH_MP INTERVAL_TO_LESS fthm;; add_test("INTERVAL_TO_LESS_CONV", let thm1 = ASSUME `interval_eps (#0.1) (float (int_of_num 1) (--: (int_of_num 1))) (float (int_of_num 1) (--: (int_of_num 2)))` in let thm2 = ASSUME `interval_eps (#7) (float (int_of_num 4) (int_of_num 1)) (float (int_of_num 1) (int_of_num 0))` in let thm3 = INTERVAL_TO_LESS_CONV thm1 thm2 in concl thm3 = `#0.1 <. (#7)`);; add_test("INTERVAL_TO_LESS_CONV2", let (h,c) = dest_thm (INTERVAL_TO_LESS_CONV (INTERVAL_OF_NUM `&.3`) (INTERVAL_OF_NUM `&.8`)) in (h=[]) && (c = `&.3 <. (&.8)`));; (* ------------------------------------------------------------------ *) (* conversion for DEC <= posfloat and posfloat <= DEC *)
let lemma1 = 
prove( `!n m p. ((&.p/(&.m)) <= (&.n)) <=> ((&.p/(&.m)) <= (&.n)/(&.1))`,
MESON_TAC[REAL_DIV_1]);;
let lemma2 = 
prove( `!n m p. ((&.p) <= ((&.n)/(&.m))) <=> ((&.p/(&.1)) <= (&.n)/(&.m))`,
MESON_TAC[REAL_DIV_1]);;
let REAL_LT = REAL_OF_NUM_LT;; let REAL_LE = REAL_OF_NUM_LE;;
let lemma3 = 
prove(`!a b c d. ( ((0<b) /\ (0<d) /\ (a*d <=| c*b)) ==> (&.a/(&.b) <=. ((&.c)/(&.d))))`,
EVERY[REPEAT GEN_TAC; DISCH_ALL_TAC; ASM_SIMP_TAC[RAT_LEMMA4;REAL_LT;REAL_OF_NUM_MUL;REAL_LE]]);;
let DEC_FLOAT = EQT_ELIM o ARITH_SIMP_CONV[DECIMAL;float;TWOPOW_POS;TWOPOW_NEG;GSYM real_div; REAL_OF_NUM_POW;INT_NUM_REAL;REAL_OF_NUM_MUL; lemma1;lemma2;lemma3];; add_test("DEC_FLOAT", let f c x = dest_thm (c x) = ([],x) in ((f DEC_FLOAT `#10.0 <= (float (int_of_num 3) (int_of_num 2))`) && (f DEC_FLOAT `#10 <= (float (int_of_num 3) (int_of_num 2))`) && (f DEC_FLOAT `#0.1 <= (float (int_of_num 1) (--: (int_of_num 2)))`) && (f DEC_FLOAT `float (int_of_num 3) (int_of_num 2) <= (#13.0)`) && (f DEC_FLOAT `float (int_of_num 3) (int_of_num 2) <= (#13)`) && (f DEC_FLOAT `float (int_of_num 1) (--: (int_of_num 2)) <= (#0.3)`)));; (* Conversion for powers *)
let pow_parity = 
prove( `!x u. x pow (NUMERAL(BIT0 u)) = (x pow (NUMERAL u) * x pow (NUMERAL u)) /\ x pow (NUMERAL(BIT1 u)) = x* x pow (NUMERAL u) * x pow (NUMERAL u)`,
REWRITE_TAC[GSYM REAL_POW_ADD;NUMERAL;BIT0;BIT1;real_pow]; );;
let pow_conv = PURE_ONCE_REWRITE_CONV[CONJUNCT1 real_pow;REAL_POW_1;pow_parity];; (* conversion for negation *) let neg_conv = let th = REAL_ARITH `--x = &0 - x` in PURE_ONCE_REWRITE_CONV[th];; (* conversion for sums *) let dest_sum = let sum = `sum:(num->bool)->(num->real)->real` in let rg = `(..)` in fun tm -> let (x,f) = dest_binop sum tm in let (i,j) = dest_binop rg x in (i,j,f);; let mk_sum = let sum = `sum:(num->bool)->(num->real)->real` in let rg = `(..)` in fun (i,j,f) -> let x = mk_binop rg i j in mk_binop sum x f;; (* mk_sum (dest_sum `sum (3..4) f`);; *) let sum_conv = fun tm -> let (i,j,f) = dest_sum tm in let i',j' = dest_numeral i,dest_numeral j in if (i'=j') then (ONCE_REWRITE_CONV[SUM_SING_NUMSEG]) tm else if (i' >/ j') then
let th = ARITH_RULE(mk_binop `(<):num->num->bool` j i) in
      let th2 = MATCH_MP SUM_TRIV_NUMSEG th in
	(PURE_ONCE_REWRITE_CONV[th2]) tm
    else 
      let j_pred = mk_numeral (j'-/ Int 1) in
      let j_new = mk_comb (`SUC`, j_pred) in
      let th = (ISPECL [i;j_pred] (CONJUNCT2 SUM_CLAUSES_NUMSEG)) in
      let rr = ARITH_RULE (mk_eq (j_new,j)) in
      let th' = REWRITE_RULE[rr;ARITH_RULE (mk_binop `(<=):num->num->bool` i j)] th in
	PURE_ONCE_REWRITE_CONV[th'] tm;;
(* ------------------------------------------------------------------ *) (* converts a DECIMAL TO A THEOREM *) let ABS_BOUNDS = REAL_ABS_BOUNDS;; (* new *)
let INTERVAL_MINMAX = 
prove(`!x f ex. ((f -. ex) <= x) /\ (x <=. (f +. ex)) ==> (interval_eps x f ex)`,
EVERY[REPEAT GEN_TAC; REWRITE_TAC[interval_eps;ABS_BOUNDS]; REAL_ARITH_TAC]);;
let INTERVAL_OF_DECIMAL bprec dec = let a_num = dest_decimal dec in let f_num = round_rat bprec a_num in let ex_num = round_rat bprec (Num.abs_num (f_num -/ a_num)) in let _ = assert (ex_num <=/ f_num) in let f = mk_float f_num in let ex= mk_float ex_num in let fplus_ex = FLOAT_CONV (mk_binop `real_add` f ex) in let fminus_ex= FLOAT_CONV (mk_binop `real_sub` f ex) in let fplus_term = rhs (concl fplus_ex) in let fminus_term = rhs (concl fminus_ex) in let th1 = DEC_FLOAT (mk_binop `real_le` fminus_term dec) in let th2 = DEC_FLOAT (mk_binop `real_le` dec fplus_term) in let intv = mk_interval dec f ex in EQT_ELIM (SIMP_CONV[INTERVAL_MINMAX;fplus_ex;fminus_ex;th1;th2] intv);; add_test("INTERVAL_OF_DECIMAL", let (h,c) = dest_thm (INTERVAL_OF_DECIMAL 4 `#36.1`) in let (x,f,ex) = dest_interval c in (h=[]) && (x = `#36.1`));; add_test("INTERVAL_OF_DECIMAL2", can (fun() -> INTERVAL_TO_LESS_CONV (INTERVAL_OF_DECIMAL 4 `#33.33`) (INTERVAL_OF_DECIMAL 4 `#36.1`)) ());;
let INTERVAL_FLOAT = 
prove_by_refinement( `!a b. interval_eps (float a b) (float a b) (float (int_of_num 0) (int_of_num 0))`,
(* {{{ proof *) [ REWRITE_TAC[interval_eps;REAL_ARITH `abs(x - x) = &0`]; REWRITE_TAC[float;int_of_num_th]; REAL_ARITH_TAC; ]);;
(* }}} *) let INTERVAL_FLOAT_CONV flt = let (a,b) = dest_binop `float` flt in (ISPECL [a;b] INTERVAL_FLOAT);; (*--------------------------------------------------------------------*) (* functions to check format. *) (* There are various implicit rules: *) (* NUMERAL is followed by bits and no other kind of num, etc. *) (* FLOAT a b, both a and b are int_of_num NUMERAL or --:int_of_num NUMERAL, etc. *) (*--------------------------------------------------------------------*) (* converts exceptions to false *) let falsify_ex f x = try (f x) with Failure _ -> false;; let is_bits_format = let zero = `_0` in let rec format x = if (x = zero) then true else let (h,t) = dest_comb x in (((h = `BIT1`) or (h = `BIT0`)) && (format t)) in falsify_ex format;; let is_numeral_format = let fn x = let (h,t) = dest_comb x in ((h = `NUMERAL`) && (is_bits_format t)) in falsify_ex fn;; let is_decimal_format = let fn x = let (t1,t2) = dest_binop `DECIMAL` x in ((is_numeral_format t1) && (is_numeral_format t2)) in falsify_ex fn;; let is_pos_int_format = let fn x = let (h,t) = dest_comb x in (h = `int_of_num`) && (is_numeral_format t) in falsify_ex fn;; let is_neg_int_format = let fn x = let (h,t) = dest_comb x in (h = `real_neg`) && (is_pos_int_format t) in falsify_ex fn;; let is_int_format x = (is_neg_int_format x) or (is_pos_int_format x);; let is_float_format = let fn x = let (t1,t2) = dest_binop `float` x in (is_int_format t1) && (is_int_format t2) in falsify_ex fn;; let is_interval_format = let fn x = let (a,b,c) = dest_interval x in (is_float_format b) && (is_float_format c) in falsify_ex fn;; let is_neg_real = let fn x = let (h,t) = dest_comb x in (h= `real_neg `) in falsify_ex fn;; let is_real_num_format = let fn x = let (h,t) = dest_comb x in (h=`real_of_num`) && (is_numeral_format t) in falsify_ex fn;; let is_comb_of t u = let fn t u = t = (fst (dest_comb u)) in try (fn t u) with Failure _ -> false;; (* ------------------------------------------------------------------ *) (* Heron's formula for the square root of A Return a value x that is always at most the actual square root and such that abs (x - A/x ) < epsilon *) let rec heron_sqrt depth A x eps = let half = (Int 1)//(Int 2) in if (depth <= 0) then raise (Failure "heron_sqrt:internal error:sqrt recursion depth exceeded") else if (Num.abs_num (x -/ (A//x) ) </ eps) & (x*/ x >=/ A) then (A//x) else let x' = half */ (x +/ (A//x)) in heron_sqrt (depth -1) A x' eps;;
let INTERVAL_OF_TWOPOW = 
prove( `!n. interval_eps (twopow n) (float (int_of_num 1) n) (float (int_of_num 0) (int_of_num 0))`,
REWRITE_TAC[interval_eps;float;int_of_num_th] THEN REAL_ARITH_TAC );;
(* compute ez parameter for division *) let division_ez bprec (f,ex,g,ey) = let f_num = dest_float f in let ex_num = dest_float ex in let g_num = dest_float g in let ey_num = dest_float ey in let _ = (ey_num </ abs_num g_num) or raise Insufficient_precision in let h_num = round_rat bprec (f_num//g_num) in let h = mk_float h_num in let ez_rat = (ex_num +/ abs_num (f_num -/ (h_num*/ g_num)) +/ (abs_num h_num */ ey_num))//((abs_num g_num) -/ (ey_num)) in let ez_num = round_rat bprec (ez_rat) in let _ = assert((ez_num >=/ (Int 0))) in let ez = mk_float ez_num in (h,ez);; let pureconv (s,t) = (s,PURE_ONCE_REWRITE_CONV[t]);; (* ------------------------------------------------------------------ *) (* bprec gives the bits of precision. *) let rec (INTERVAL_OF_TERM:(string*conv)list -> (term *int*thm)list ->int -> term -> (term*int*thm) list * thm) = fun convl cache bprec tm -> let mk c th = (cache_entry_of_ineq bprec th :: c, th) in (* treat user conversions *) let interval_of_conv cache bprec conv tm = let th_r = conv tm in let (lhs,rhs) = dest_eq (concl th_r) in let _ = (not (aconv lhs rhs)) or (print_term lhs; failwith "interval_of_conv: no change" )in let (cache,th_rhs) = INTERVAL_OF_TERM convl cache bprec rhs in let th_lhs = interval_rule (SYM th_r) th_rhs in mk cache th_lhs in let interval_of_binop op opname match_thm = begin let (a,b) = dest_binop op tm in let (cache,a_int) = INTERVAL_OF_TERM convl cache bprec a in let (cache,b_int) = INTERVAL_OF_TERM convl cache bprec b in try( let c_int = MATCH_MP match_thm (CONJ a_int b_int) in let th = EQ_MP (INTERVAL_CLEAN_CONV (concl c_int)) c_int in mk cache th) with Failure _ -> failwith ("INTERVAL_OF_TERM: "^opname) end in (* treat cached values first *) if (can (get_real_constant cache bprec) tm) then begin try( let int_thm = get_real_constant cache bprec tm in if (can dest_interval (concl int_thm)) then (cache,int_thm) else ( let absthm = int_thm in let (abst,k) = dest_binop `real_le` (concl absthm) in let (absh,cmu) = dest_comb abst in let (c,u) = dest_binop `real_sub` cmu in let (cache,intk) = INTERVAL_OF_TERM convl cache bprec k in let (cache,intu) = INTERVAL_OF_TERM convl cache bprec u in let thm1 = MATCH_MP ABS_TO_INTERVAL absthm in let thm2 = MATCH_MP thm1 (CONJ intu intk) in let (_,f,ex)= dest_interval (concl thm2) in let fthm = FLOAT_CONV f in let exthm = FLOAT_CONV ex in let thm3 = REWRITE_RULE[fthm;exthm] thm2 in (cache_entry_of_ineq bprec thm3 :: (remove_entry (tm,bprec,int_thm) cache), thm3) )) with Failure _ -> failwith "INTERVAL_OF_TERM : CONSTANT" end else if (is_real_num_format tm) then let th = INTERVAL_OF_NUM tm in (cache,th) else if (is_decimal_format tm) then let th = (INTERVAL_OF_DECIMAL bprec tm) in mk cache th (* treat abs value *) else if (is_comb_of `real_abs` tm) then begin let (_,b) = dest_comb tm in let (cache,i1) = (INTERVAL_OF_TERM convl cache bprec b) in try( let b_int = MATCH_MP INTERVAL_ABSV i1 in let th = EQ_MP (INTERVAL_CLEAN_CONV (concl b_int)) b_int in mk cache th) with Failure _ -> failwith "INTERVAL_OF_TERM : ABS" end (* treat twopow *) else if (is_comb_of `twopow` tm) then begin try( let (_,b) = dest_comb tm in let th = SPEC b INTERVAL_OF_TWOPOW in mk cache th ) with Failure _ -> failwith "INTERVAL_OF_TERM : TWOPOW" end (* treat addition *) else if (can (dest_binop `real_add`) tm) then interval_of_binop `real_add` "ADD" INTERVAL_ADD (* treat subtraction *) else if (can (dest_binop `real_sub`) tm) then interval_of_binop `real_sub` "SUB" INTERVAL_SUB (* treat multiplication *) else if (can (dest_binop `real_mul`) tm) then interval_of_binop `real_mul` "MUL" INTERVAL_MUL (* treat division : instantiate INTERVAL_DIV *) else if (can (dest_binop `( / )`) tm) then begin let (a,b) = dest_binop `( / )` tm in let (cache,a_int) = INTERVAL_OF_TERM convl cache bprec a in let (cache,b_int) = INTERVAL_OF_TERM convl cache bprec b in try( let (_,f,ex) = dest_interval (concl a_int) in let (_,g,ey) = dest_interval (concl b_int) in let (h,ez) = division_ez bprec (f,ex,g,ey) in let hyp1 = a_int in let hyp2 = b_int in let hyp3 = try (FLOAT_INEQ_CONV (mk_binop `real_lt` ey (mk_comb (`real_abs`,g)))) with Failure _ -> raise Insufficient_precision in let thm = SPECL [a;f;ex;b;g;ey;h;ez] INTERVAL_DIV in let conj2 x = snd (dest_conj x) in let hyp4t = (conj2 (conj2 (conj2 (fst(dest_imp (concl thm)))))) in let hyp4 = FLOAT_INEQ_CONV hyp4t in let hyp_all = end_itlist CONJ [hyp1;hyp2;hyp3;hyp4] in let th = MATCH_MP thm hyp_all in mk cache th) with Failure _ -> failwith "INTERVAL_OF_TERM :DIV" end (* treat pow *) else if (can (dest_binary "real_pow") tm) then interval_of_conv cache bprec pow_conv tm (* treat negative terms *) else if (is_neg_real tm) then interval_of_conv cache bprec neg_conv tm (* treat beta *) else if can (fun t-> let (f,x) = dest_comb t in dest_abs f ) tm then interval_of_conv cache bprec BETA_CONV tm (* treat sum *) else if can dest_sum tm then interval_of_conv cache bprec sum_conv tm (* treat float *) else if can (dest_binop `float`) tm then (cache,INTERVAL_FLOAT_CONV tm) (* treat ssqrt : instantiate INTERVAL_SSQRT *) else if (can (dest_unary `ssqrt`) tm) then begin try( let x = dest_unary `ssqrt` tm in let (cache,x_int) = INTERVAL_OF_TERM convl cache bprec x in let (_,f,ex) = dest_interval (concl x_int) in let f_num = dest_float f in let ex_num = dest_float ex in let fd_num = f_num -/ ex_num in let fe_f = Num.float_of_num fd_num in let apprx_sqrt = Pervasives.sqrt fe_f in (* put in heron's formula *) let v_num1 = round_inward_float 25 (apprx_sqrt) in let v_num = round_inward_rat bprec (heron_sqrt 10 fd_num v_num1 ((Int 2) **/ (Int (-bprec-4)))) in let u_num1 = round_inward_float 25 (Pervasives.sqrt (float_of_num f_num)) in let u_num = round_inward_rat bprec (heron_sqrt 10 f_num u_num1 ((Int 2) **/ (Int (-bprec-4)))) in let ey_num = round_rat bprec (abs_num (f_num -/ (u_num */ u_num))) in let ez_num = round_rat bprec ((ex_num +/ ey_num)//(u_num +/ v_num)) in let (v,u) = (mk_float v_num,mk_float u_num) in let (ey,ez) = (mk_float ey_num,mk_float ez_num) in let thm = SPECL [x;f;ex;u;ey;ez;v] INTERVAL_SSQRT in let conjhyp = fst (dest_imp (concl thm)) in let [hyp6;hyp5;hyp4;hyp3;hyp2;hyp1] = let rec break_conj c acc = if (not(is_conj c)) then (c::acc) else let (u,v) = dest_conj c in break_conj v (u::acc) in (break_conj conjhyp []) in
let thm2 = prove(hyp2,REWRITE_TAC[interval_eps] THEN
                       (CONV_TAC FLOAT_INEQ_CONV)) in
    let thm3 = FLOAT_INEQ_CONV hyp3 in
    let thm4 = FLOAT_INEQ_CONV hyp4 in
    let float_tac = REWRITE_TAC[FLOAT_NN;FLOAT_POS;INT_NN;INT_NN_NEG;
                       INT_POS';
INT_POS_NEG] THEN ARITH_TAC in
let thm5 = prove( hyp5,float_tac) in
    let thm6 = prove( hyp6,float_tac) in
    let ant  = end_itlist CONJ[x_int;thm2;thm3;thm4;thm5;thm6] in
    let th =    MATCH_MP thm ant
    in mk cache th
    )
    with Failure _ -> failwith "INTERVAL_OF_TERM  : SSQRT"
    end
 else if (can (dest_unary `sqrt`) tm) then
    begin
    try(
    let x = dest_unary `sqrt` tm in
    let stm = mk_comb (`ssqrt`, x) in
    let (cache,ssqrt_thm) = INTERVAL_OF_TERM convl cache  bprec stm in
    let (cache,int1) = INTERVAL_OF_TERM convl cache  bprec`real_of_num 0` in
    let (cache,int2)  =INTERVAL_OF_TERM convl cache  bprec x in
    let ant = CONJ (INTERVAL_TO_LESS_CONV int1 int2) ssqrt_thm in
    let th = MATCH_MP ssqrt_sqrt ant
    in mk cache th)
    with  Failure _ -> failwith "INTERVAL_OF_TERM : SQRT"
    end
(* treat atn *)
 else if can (dest_unary `atn`) tm then
   begin
     let half_co = ("halfatn4_co",ARITH_SIMP_CONV[Taylor_atn.halfatn4_co]) in
     let half4 = ("halfatn4",PURE_REWRITE_CONV[Taylor_atn.halfatn4;o_THM]) in
     let halfatnconv = pureconv("halfatn",Taylor_atn.halfatn) in
     let aconvl =  (halfatnconv::half4::half_co::convl) in
     let n = mk_numeral (Int (1 + (bprec - 5) / 6)) in
     let x = dest_unary `atn` tm in
     let sum_tm = instantiate ([],[(n,`n:num`);(x,`u:real`)],[]) `&16 * sum (0..n) (halfatn4_co u)` in
     let (cache,hyp_1) = INTERVAL_OF_TERM (aconvl) cache bprec sum_tm in
     let (_,_,eps2) = dest_interval (concl hyp_1) in
     let eps_tm = instantiate([],[(n,`n:num`);(eps2,`eps2:real`)],[]) `float (int_of_num 1) ( -- int_of_num (6 * n + 5)) + eps2` in
     let eps_thm = FLOAT_CONV eps_tm in
     let eps = snd(dest_eq (concl eps_thm)) in
     let hyp_2  = FLOAT_INEQ_CONV (mk_binop `(<=):real->real->bool` eps_tm eps) in (cache, MATCH_MP  taylor_atn   (CONJ hyp_1 hyp_2) )
   end 
(* treat acs *)
 else if can (dest_unary `acs`) tm then
   begin
     let x = dest_unary `acs` tm in
     let (_,int_m1) = INTERVAL_OF_TERM convl cache  bprec`-- real_of_num 1` in
     let (_,int_p1) = INTERVAL_OF_TERM convl cache  bprec`real_of_num 1` in
     let (cache,int_x)  =INTERVAL_OF_TERM convl cache  bprec x in
     let ant1 = INTERVAL_TO_LESS_CONV int_m1 int_x in
     let ant2 = INTERVAL_TO_LESS_CONV int_x int_p1 in
     let th =  MATCH_MP ACS_ATN (CONJ ant1 ant2) in
       interval_of_conv cache bprec (PURE_ONCE_REWRITE_CONV[th]) tm
   end
(* treat pi *)
 else if (tm = `pi`) then
    begin
      interval_of_conv cache bprec (PURE_ONCE_REWRITE_CONV[pi_atn]) tm
    end
(* user supplied conversions *)
  else 
    begin
      let (c,b) = strip_comb tm in
      let _ = is_const c or failwith "INTERVAL_OF_TERM: case not installed" in
      let (s,_) = dest_const c in
      let conv = try (assoc s convl) 
         with Failure _ -> failwith ("INTERVAL_OF_TERM: function not found "^s) in
         interval_of_conv cache bprec conv tm
    end;;
let real_ineq convl cache bprec tm = let (t1,t2) = try (dest_binop `real_lt` tm) with Failure _ -> failwith "real_ineq: inequality of the form `x < y` expected" in let (cache,int1) = INTERVAL_OF_TERM convl cache bprec t1 in let (cache,int2) = INTERVAL_OF_TERM convl cache bprec t2 in (cache,INTERVAL_TO_LESS_CONV int1 int2);; let rec rec_real_ineq convl cache bprec bprec_max tm = try (bprec,real_ineq convl cache bprec tm) with Insufficient_precision -> let bprec = (16 * bprec) / 10 in let _ = (bprec < bprec_max) or failwith ("insufficient precision "^ (string_of_int bprec_max)) in let _ = Format.print_string ("increasing precision to "^(string_of_int bprec)^" bits"); Format.print_newline() in rec_real_ineq convl cache bprec bprec_max tm;; (* Give it a conjunction of strict inequalities *) let REAL_INEQUALITY_CALCULATOR_PREC (bprec_init,bprec_max) convl tm = let tml = striplist dest_conj tm in let (b,c,th) = List.fold_left (fun (b,c,thl) tm1 -> let (bprec,(cache,th1)) = rec_real_ineq convl c b bprec_max tm1 in (bprec,cache, th1 ::thl)) (bprec_init,[], []) (tml) in (b, EQ_MP (SYM(REWRITE_CONV th tm)) TRUTH);; let REAL_INEQUALITY_CALCULATOR convl tm = REAL_INEQUALITY_CALCULATOR_PREC (10,70) convl tm;; (* let REAL_INEQUALITY_CALCULATOR convl tm = let tml = striplist dest_conj tm in let bprec_max = 70 in let (b,c,th) = List.fold_left (fun (b,c,thl) tm1 -> let (bprec,(cache,th1)) = rec_real_ineq convl c b bprec_max tm1 in (bprec,cache, th1 ::thl)) (10,[], []) (tml) in (b, EQ_MP (SYM(REWRITE_CONV th tm)) TRUTH);; *) (* let tmint = map (rec_real_ineq [] 10) tml in List.fold_left CONJ (hd tmint) (tl tmint);; *) end;;