text_formalization/jordan/float.hl

   1 (* ========================================================================== *)
   2 (* FLYSPECK - BOOK FORMALIZATION                                              *)
   3 (*                                                                            *)
   4 (* Chapter: Jordan                                                               *)
   5 (* Copied from HOL Light jordan directory *)
   6 (* Author: Thomas C. Hales                                                    *)
   7 (* Date: 2010-07-08                                                           *)
   8 (* ========================================================================== *)
   9
  10 (*
  11 A handheld calculator, implemented in HOL light.
  12
  13 It implements + - / * sqrt, ssqrt, abs, atn, acs, pi,
  14 The user can supply additional conversons to extend the calculator to other functions.
  15
  16 There is an example file float_example.hl.
  17 *)
  18
  19 (* based on float.ml in hol/jordan/float.ml *)
  20
  21 (*
  22 To do. Allow the user to supply theorems that carry domain constraints,
  23  (such as ACS_ATN).
  24 *)
  25
  26 flyspeck_needs "jordan/taylor_atn.hl";;
  27
  28 module Float = struct
  29
  30 open Lib_ext;;
  31 open Num_ext_nabs;;
  32 open Real_ext;;
  33 open Refinement;;
  34 open Tactics_jordan;;
  35 open Taylor_atn;;
  36
  37 Parse_ext_override_interface.unambiguous_interface();;
  38 Parse_ext_override_interface.prioritize_real();;
  39
  40 (* Three types of error:
  41      Insufficient precision: the inequality might be true but there is not enough
  42      precision to decide.  Divide by zero generally lands here, because there is
  43      no testing for exact equality.
  44
  45      Internal failure: a bug in this code.
  46
  47      User input failure: sqrt of negative numbers, false inequalities.
  48
  49 *)
  50
  51 exception Insufficient_precision;;
  52
  53 (*
  54 exception Mk_comb of (term*hol_type)*(term*hol_type);;
  55
  56 let mk_comb(a,b) = try mk_comb(a,b) with
  57   Failure _ -> raise (Mk_comb((a,type_of a),(b,type_of b)));;
  58 *)
  59 let add_test,test = new_test_suite();;
  60
  61 let twopow =
  62   new_definition(
  63         `twopow x = if (?n. (x = (int_of_num n)))
  64         then ((&2) pow (nabs x))
  65         else inv((&2) pow (nabs x))`);;
  66
  67 let float =
  68   new_definition(
  69                   `float x n = (real_of_int x)*(twopow n)`);;
  70
  71 let interval_eps =
  72   new_definition(
  73                    `interval_eps x f eps = ((abs (x-f)) <= eps)`);;
  74
  75 (*--------------------------------------------------------------------*)
  76
  77
  78 let mk_interval a b ex =
  79    mk_comb(mk_comb (mk_comb (`interval_eps`,a),b),ex);;
  80
  81 add_test("mk_interval",
  82    mk_interval `#3` `#4` `#1` = `interval_eps #3 #4 #1`);;
  83
  84 let dest_interval intv =
  85    let (h1,ex) = dest_comb intv in
  86    let (h2,f) = dest_comb h1 in
  87    let (h3,a) = dest_comb h2 in
  88    let _ = (h3 = `interval_eps`) or (print_term h3; failwith "dest_interval:interval_eps") in
  89    (a,f,ex);;
  90
  91 let _ = add_test("dest_interval",
  92    let a = `#3` and b = `#4` and c = `#1` in
  93    dest_interval (mk_interval a b c) = (a,b,c));;
  94
  95 (*--------------------------------------------------------------------*)
  96
  97 let (dest_int:term-> Num.num) =
  98   fun b ->
  99   let dest_pos_int a =
 100     let (op,nat) = dest_comb a in
 101     if (fst (dest_const op) = "int_of_num") then (dest_numeral nat)
 102       else fail() in
 103     let (op',u) = (dest_comb b) in
 104     try (if (fst (dest_const op') = "int_neg") then
 105            minus_num (dest_pos_int u) else
 106              dest_pos_int b) with
 107         Failure _ -> failwith "dest_int ";;
 108
 109
 110 let (mk_int:Num.num -> term) =
 111   fun a ->
 112     let sgn = Num.sign_num a in
 113     let abv = Num.abs_num a in
 114     let r = mk_comb(`int_of_num`,mk_numeral abv) in
 115     try (if (sgn<0) then mk_comb (`int_neg`,r) else r) with
 116         Failure _ -> failwith ("dest_int "^(string_of_num a));;
 117
 118 add_test("mk_int",
 119    (mk_int (Int (-1443)) = `int_neg (int_of_num 1443)`) &&
 120    (mk_int (Int 37) = `(int_of_num 37)`));;
 121
 122 (* ------------------------------------------------------------------ *)
 123
 124 let (split_ratio:Num.num -> Num.num*Num.num) =
 125   function
 126     (Ratio r) -> (Big_int (Ratio.numerator_ratio r)),
 127          (Big_int (Ratio.denominator_ratio r))|
 128     u -> (u,(Int 1));;
 129
 130 add_test("split_ratio",
 131    let (a,b) = split_ratio ((Int 4)//(Int 20)) in
 132    (a =/ (Int 1)) && (b =/ (Int 5)));;
 133
 134 (* ------------------------------------------------------------------ *)
 135
 136 (* break nonzero int r into a*(C**b) with a prime to C . *)
 137 let (factor_C:int -> Num.num -> Num.num*Num.num) =
 138   function c ->
 139   let intC = (Int c) in
 140   let rec divC (a,b) =
 141     if ((Int 0) =/ mod_num a intC) then (divC (a//intC,b+/(Int 1)))
 142       else (a,b) in
 143   function r->
 144   if ((Num.is_integer_num r)&& not((Num.sign_num r) = 0)) then
 145     divC (r,(Int 0)) else failwith "factor_C";;
 146
 147 add_test("factor_C",
 148    (factor_C 2 (Int (4096+32)) = (Int 129,Int 5)) &&
 149    (factor_C 10 (Int (5000)) = (Int 5,Int 3)) &&
 150    (cannot (factor_C 2) ((Int 50)//(Int 3))));;
 151
 152 (*--------------------------------------------------------------------*)
 153
 154 let (dest_float:term -> Num.num) =
 155   fun f ->
 156     let (a,b) = dest_binop `float` f in
 157     (dest_int a)*/ ((Int 2) **/ (dest_int b));;
 158
 159 add_test("dest_float",
 160    dest_float `float (int_of_num 3) (int_of_num 17)` = (Int 393216));;
 161
 162 add_test("dest_float2", (* must express as numeral first *)
 163    cannot dest_float `float ((int_of_num 3)+:(int_of_num 1)) (int_of_num 17)`);;
 164
 165 (* ------------------------------------------------------------------ *)
 166 (* creates float of the form `float a b` with a odd *)
 167 let (mk_float:Num.num -> term) =
 168   function r ->
 169     let (a,b) = split_ratio r in
 170     let (a',exp_a) = if (a=/(Int 0)) then ((Int 0),(Int 0)) else factor_C 2 a in
 171     let (b',exp_b) = factor_C 2 b in
 172     let c = a'//b' in
 173     if (Num.is_integer_num c) then
 174           mk_binop `float` (mk_int c) (mk_int (exp_a -/ exp_b))
 175           else failwith "mk_float";;
 176
 177 add_test("mk_float",
 178    mk_float (Int (4096+32)) = `float (int_of_num 129) (int_of_num 5)` &&
 179    (mk_float (Int 0) = `float (int_of_num 0) (int_of_num 0)`));;
 180
 181 add_test("mk_float2",  (* throws exception exactly when denom != 2^k *)
 182    let rtest = fun t -> (t =/ dest_float (mk_float t)) in
 183    rtest ((Int 3)//(Int 1024)) &&
 184   (cannot rtest ((Int 1)//(Int 3))));;
 185
 186 add_test("mk_float dest_float",  (* constructs canonical form of float *)
 187   mk_float (dest_float `float (int_of_num 4) (int_of_num 3)`) = `float (int_of_num 1) (int_of_num 5)`);;
 188
 189 (* ------------------------------------------------------------------ *)
 190 (* creates decimal of the form `DECIMAL a b` with a prime to 10 *)
 191 let (mk_pos_decimal:Num.num -> term) =
 192   let zero = `#0` in
 193   function r ->
 194     let _ = assert (r >=/ (Int 0)) in
 195     let (a,b) = split_ratio r in
 196     if (a=/(Int 0)) then zero else
 197     let (a1,exp_a5) = factor_C 5 a in
 198     let (a2,exp_a2) = factor_C 2 a1 in
 199     let (b1,exp_b5) = factor_C 5 b in
 200     let (b2,exp_b2) = factor_C 2 b1 in
 201     let _ = assert(b2 =/ (Int 1)) in
 202     let c = end_itlist Num.max_num [exp_b5-/exp_a5;exp_b2-/exp_a2;(Int 0)] in
 203     let a' = a2*/((Int 2)**/ (c +/ exp_a2 -/ exp_b2))*/
 204              ((Int 5)**/(c +/ exp_a5 -/ exp_b5)) in
 205     let b' = (Int 10) **/ c in
 206     mk_binop `DECIMAL` (mk_numeral a') (mk_numeral b');;
 207
 208 add_test("mk_pos_decimal",
 209    mk_pos_decimal (Int (5000)) = `#5000` &&
 210    (mk_pos_decimal ((Int 30)//(Int 40)) = `#0.75`) &&
 211    (mk_pos_decimal (Int 0) = `#0`) &&
 212    (mk_pos_decimal ((Int 15)//(Int 25)) = `#0.6`) &&
 213    (mk_pos_decimal ((Int 25)//(Int 4)) = `#6.25`) &&
 214    (mk_pos_decimal ((Int 2)//(Int 25)) = `#0.08`));;
 215
 216 let (mk_decimal:Num.num->term) =
 217   function r ->
 218   let a = Num.sign_num r in
 219   let b = mk_pos_decimal (Num.abs_num r) in
 220   if (a < 0) then (mk_comb (`real_neg`, b)) else b;;
 221
 222 add_test("mk_decimal",
 223   (mk_decimal (Int 3) = `#3`) &&
 224   (mk_decimal (Int (-3)) = `real_neg  (#3)`));;
 225
 226
 227
 228 (*--------------------------------------------------------------------*)
 229
 230 let (dest_decimal:term -> Num.num) =
 231   fun f ->
 232     let (a,b) = dest_binop `DECIMAL` f in
 233     let a1 = dest_numeral a in
 234     let b1 = dest_numeral b in
 235         a1//b1;;
 236
 237 add_test("dest_decimal",
 238    dest_decimal `#3.4` =/ ((Int 34)//(Int 10)));;
 239 add_test("dest_decimal2",
 240    cannot dest_decimal `real_neg  (#3.4)`);;
 241
 242
 243
 244
 245
 246 (*--------------------------------------------------------------------*)
 247 (*   Properties of integer powers of 2.                               *)
 248 (* ------------------------------------------------------------------ *)
 249
 250
 251 let TWOPOW_POS = prove(`!n. (twopow (int_of_num n) = (&2) pow n)`,
 252         (REWRITE_TAC[twopow])
 253         THEN GEN_TAC
 254         THEN COND_CASES_TAC
 255         THENL [AP_TERM_TAC;ALL_TAC]
 256         THEN (REWRITE_TAC[NABS_POS])
 257         THEN (UNDISCH_EL_TAC 0)
 258         THEN (TAUT_TAC (` ( A    ) ==> (~ A ==> B)`))
 259         THEN (MESON_TAC[]));;
 260
 261 let TWOPOW_NEG = prove(`!n. (twopow (--(int_of_num n)) = inv((&2) pow n))`,
 262         GEN_TAC
 263         THEN (REWRITE_TAC[twopow])
 264         THEN (COND_CASES_TAC THENL [ALL_TAC;REWRITE_TAC[NABS_NEG]])
 265         THEN (POP_ASSUM CHOOSE_TAC)
 266         THEN (REWRITE_TAC[NABS_NEG])
 267         THEN (UNDISCH_EL_TAC 0)
 268         THEN (REWRITE_TAC[int_eq;int_neg_th;INT_NUM_REAL])
 269         THEN (REWRITE_TAC[prove (`! u y.((--(real_of_num u) = (real_of_num y))=
 270                 ((real_of_num u) +(real_of_num y) = (&0)))`,REAL_ARITH_TAC)])
 271         THEN (REWRITE_TAC[REAL_OF_NUM_ADD;REAL_OF_NUM_EQ;ADD_EQ_0])
 272         THEN (DISCH_TAC)
 273         THEN (ASM_REWRITE_TAC[real_pow;REAL_INV_1]));;
 274
 275
 276 let TWOPOW_INV = prove(`!a. (twopow (--: a) = (inv (twopow a)))`,
 277   (GEN_TAC)
 278   THEN ((ASSUME_TAC (SPEC `a:int` INT_REP2)))
 279   THEN ((POP_ASSUM CHOOSE_TAC))
 280   THEN ((POP_ASSUM DISJ_CASES_TAC))
 281   THEN ((ASM_REWRITE_TAC[TWOPOW_POS;TWOPOW_NEG;REAL_INV_INV;INT_NEG_NEG])));;
 282
 283 let INT_REP3 = prove(`!a .(?n.( (a = int_of_num  n) \/ (a = --: (int_of_num  (n+1)))))`,
 284 (GEN_TAC)
 285 THEN ((ASSUME_TAC (SPEC `a:int` INT_REP2)))
 286 THEN ((POP_ASSUM CHOOSE_TAC))
 287 THEN ((DISJ_CASES_TAC (prove (`((a:int) = (int_of_num  0)) \/ ~((a:int) =(int_of_num  0))`, MESON_TAC[]))))
 288 (* cases *)
 289 THENL[ ((EXISTS_TAC `0`)) THEN ((ASM_REWRITE_TAC[]));ALL_TAC]
 290 THEN ((UNDISCH_EL_TAC 0))
 291 THEN ((POP_ASSUM DISJ_CASES_TAC))
 292 THENL [DISCH_TAC THEN ((ASM MESON_TAC)[]);ALL_TAC]
 293 THEN (DISCH_TAC)
 294 THEN ((EXISTS_TAC `PRE n`))
 295 THEN ((DISJ2_TAC))
 296 THEN ((ASM_REWRITE_TAC[INT_EQ_NEG2]))
 297 (*** Changed by JRH, 2006/03/28 to avoid PRE_ELIM_TAC ***)
 298 THEN (FIRST_X_ASSUM(MP_TAC o check(is_neg o concl)))
 299 THEN (ASM_REWRITE_TAC[INT_NEG_EQ_0; INT_OF_NUM_EQ])
 300 THEN (ARITH_TAC));;
 301
 302 let REAL_EQ_INV = prove(`!x y. ((x:real = y) <=> (inv(x) = inv (y)))`,
 303 ((REPEAT GEN_TAC))
 304 THEN (EQ_TAC)
 305 THENL [((DISCH_TAC THEN (ASM_REWRITE_TAC[])));
 306  (* branch 2*) ((DISCH_TAC))
 307 THEN ((ONCE_REWRITE_TAC [(GSYM REAL_INV_INV)]))
 308 THEN ((ASM_REWRITE_TAC[]))]);;
 309
 310 let TWOPOW_ADD_1 =
 311   prove(`!a. (twopow (a +: (int_of_num 1)) = twopow (a) *. (twopow (int_of_num 1)))`,
 312 EVERY[
 313   GEN_TAC;
 314   CHOOSE_TAC (SPEC `a:int` INT_REP3);
 315   POP_ASSUM DISJ_CASES_TAC
 316   THENL[
 317     ASM_REWRITE_TAC[TWOPOW_POS;INT_OF_NUM_ADD;REAL_POW_ADD];
 318     EVERY[
 319       ASM_REWRITE_TAC[GSYM INT_OF_NUM_ADD;INT_NEG_ADD;GSYM INT_ADD_ASSOC;INT_ADD_LINV;INT_ADD_RID];
 320       REWRITE_TAC[GSYM INT_NEG_ADD;INT_OF_NUM_ADD;TWOPOW_NEG;TWOPOW_POS];
 321       ONCE_REWRITE_TAC[SPEC `(&. 2) pow 1` (GSYM REAL_INV_INV)];
 322       REWRITE_TAC[GSYM REAL_INV_MUL;GSYM REAL_EQ_INV;REAL_POW_ADD;GSYM REAL_MUL_ASSOC;REAL_POW_1];
 323       REWRITE_TAC[MATCH_MP REAL_MUL_RINV (REAL_ARITH `~((&. 2) = (&. 0))`); REAL_MUL_RID]
 324     ]
 325   ]
 326 ]);;
 327
 328
 329 let TWOPOW_0 = prove(`twopow (int_of_num  0) = (&. 1)`,
 330  (REWRITE_TAC[TWOPOW_POS;real_pow]));;
 331
 332 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)))`,
 333 ((INDUCT_TAC))
 334 THENL [REWRITE_TAC[INT_SUB_LZERO;REAL_MUL_LID;TWOPOW_0];ALL_TAC]
 335 THEN ((INDUCT_TAC THEN
 336    ( (ASM_REWRITE_TAC[INT_SUB_RZERO;TWOPOW_0;REAL_MUL_RID;INT_NEG_0;ADD1;GSYM INT_OF_NUM_ADD]))))
 337 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)]))
 338 THEN ((REWRITE_TAC[REAL_INV_MUL]))
 339 THEN ((ABBREV_TAC `a:real = twopow (int_of_num  m)`))
 340 THEN ((ABBREV_TAC `b:real = inv(twopow (int_of_num  n))`))
 341 THEN ((REWRITE_TAC[TWOPOW_POS;REAL_POW_1;GSYM REAL_MUL_ASSOC;prove (`!(x:real). ((&.2)*x = x*(&.2))`,REAL_ARITH_TAC)]))
 342 THEN ((REWRITE_TAC[REAL_INV2;REAL_MUL_RID])));;
 343
 344 let TWOPOW_ADD_NUM = prove(
 345   `!m n. (twopow ((int_of_num m) + (int_of_num n)) = twopow((int_of_num m))*. twopow((int_of_num n)))`,
 346 (REWRITE_TAC[TWOPOW_POS;REAL_POW_ADD;INT_OF_NUM_ADD]));;
 347
 348 let TWOPOW_ADD_INT = prove(
 349   `!a b. (twopow (a +: b) = twopow(a) *. (twopow(b)))`,
 350  ((REPEAT GEN_TAC))
 351 THEN ((ASSUME_TAC (SPEC `a:int` INT_REP)))
 352 THEN ((POP_ASSUM CHOOSE_TAC))
 353 THEN ((POP_ASSUM CHOOSE_TAC))
 354 THEN ((ASSUME_TAC (SPEC `b:int` INT_REP)))
 355 THEN ((REPEAT (POP_ASSUM CHOOSE_TAC)))
 356 THEN ((ASM_REWRITE_TAC[]))
 357 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'))`))
 358 (* branch *)
 359 THENL[ ((REWRITE_TAC[GSYM INT_OF_NUM_ADD]))
 360 THEN ((INT_ARITH_TAC));ALL_TAC]
 361 (* 2nd *)
 362 THEN ((DISCH_TAC))
 363 THEN ((ASM_REWRITE_TAC[TWOPOW_SUB_NUM;TWOPOW_INV;TWOPOW_POS;REAL_POW_ADD;REAL_INV_MUL;GSYM REAL_MUL_ASSOC]))
 364 THEN ((ABBREV_TAC `a':real = inv ((&. 2) pow m)`))
 365 THEN ((ABBREV_TAC `c :real = (&. 2) pow n`))
 366 THEN ((ABBREV_TAC `d :real = (&. 2) pow n'`))
 367 THEN ((ABBREV_TAC `e :real = inv ((&. 2) pow m')`))
 368 THEN (MESON_TAC[REAL_MUL_AC]));;
 369
 370 let TWOPOW_ABS = prove(`!a. real_abs  (twopow a) = (twopow a)`,
 371 (GEN_TAC)
 372 THEN ((CHOOSE_THEN DISJ_CASES_TAC (SPEC `a:int` INT_REP2)))
 373 (* branch *)
 374 THEN ((ASM_REWRITE_TAC[TWOPOW_POS;TWOPOW_NEG;REAL_ABS_POW;REAL_ABS_NUM;REAL_ABS_INV])));;
 375
 376 let REAL_POW_POW = prove(
 377   `!x m n . (x **. m) **. n = x **. (m *| n)`,
 378 ((GEN_TAC THEN GEN_TAC THEN INDUCT_TAC))
 379 (* branch *)
 380 THENL[ ((REWRITE_TAC[real_pow;MULT_0]));
 381 (* second branch *)
 382 ((REWRITE_TAC[real_pow]))
 383 THEN ((ASM_REWRITE_TAC[ADD1;LEFT_ADD_DISTRIB;REAL_POW_ADD;REAL_MUL_AC;MULT_CLAUSES]))]);;
 384
 385 let INT_POW_POW = INT_OF_REAL_THM REAL_POW_POW;;
 386
 387 let TWOPOW_POW = prove(
 388   `!a n. (twopow a) pow n = twopow (a *: (int_of_num  n))`,
 389 ((REPEAT GEN_TAC))
 390 THEN ((CHOOSE_THEN DISJ_CASES_TAC (SPEC `a:int` INT_REP2)))
 391 (* branch *)
 392 THEN ((ASM_REWRITE_TAC[TWOPOW_POS;INT_OF_NUM_MUL;
 393    REAL_POW_POW;TWOPOW_NEG;REAL_POW_INV;INT_OF_NUM_MUL;GSYM INT_NEG_LMUL])));;
 394
 395 (* ------------------------------------------------------------------ *)
 396 (*   Arithmetic operations on float                                   *)
 397 (* ------------------------------------------------------------------ *)
 398
 399
 400 let FLOAT_NEG = prove(`!a m. real_neg  (float a m) = float (--: a) m`,
 401  REPEAT GEN_TAC THEN
 402  REWRITE_TAC[float;GSYM REAL_MUL_LNEG;int_neg_th]);;
 403
 404
 405
 406 let FLOAT_MUL = prove(`!a b m n. (float a m) *. (float b n) = (float (a *: b) (m +: n))`,
 407 ((REPEAT GEN_TAC))
 408 THEN ((REWRITE_TAC[float;GSYM REAL_MUL_ASSOC;TWOPOW_ADD_INT;int_mul_th]))
 409 THEN ((MESON_TAC[REAL_MUL_AC])));;
 410
 411 let FLOAT_ADD = prove(
 412   `!a b c m. (float a (m+: (int_of_num c))) +. (float b m)
 413      = (float ( (int_of_num (2 EXP c))*a +: b) m)`,
 414 ((REWRITE_TAC[float;int_add_th;REAL_ADD_RDISTRIB;int_mul_th;TWOPOW_ADD_INT]))
 415 THEN ((REPEAT GEN_TAC))
 416 THEN ((REWRITE_TAC[TWOPOW_POS;INT_NUM_REAL;GSYM REAL_OF_NUM_POW]))
 417 THEN ((MESON_TAC[REAL_MUL_AC])));;
 418
 419 let FLOAT_ADD_EQ = prove(
 420   `!a b m. (float a  m) +. (float b m) =
 421   (float (a+:b) m)`,
 422  ((REPEAT GEN_TAC))
 423 THEN ((REWRITE_TAC[REWRITE_RULE[INT_ADD_RID] (SPEC `m:int` (SPEC `0` (SPEC `b:int` (SPEC `a:int` FLOAT_ADD))))]))
 424 THEN ((REWRITE_TAC[EXP;INT_MUL_LID])));;
 425
 426 let FLOAT_ADD_NP = prove(
 427   `!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)))`,
 428 (REWRITE_TAC[REAL_ADD_AC]));;
 429
 430 let FLOAT_ADD_PN = prove(
 431   `!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)))`,
 432 ((REPEAT GEN_TAC))
 433 THEN ((SUBGOAL_MP_TAC `int_of_num  m = (--:(int_of_num  n)) + (int_of_num (m+n))`))
 434 THENL[ ((REWRITE_TAC[GSYM INT_OF_NUM_ADD]))
 435 THEN ((INT_ARITH_TAC));
 436 (* branch *)
 437 ((DISCH_TAC))
 438 THEN ((ASM_REWRITE_TAC[FLOAT_ADD]))]);;
 439
 440 let FLOAT_ADD_PP = prove(
 441   `!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))))`,
 442 ((REPEAT GEN_TAC))
 443 THEN (DISCH_TAC)
 444 THEN ((SUBGOAL_MP_TAC `int_of_num  m = (int_of_num  n) + (int_of_num  (m-n))`))
 445 THENL[ ((REWRITE_TAC[INT_OF_NUM_ADD]))
 446 THEN (AP_TERM_TAC)
 447 THEN ((REWRITE_TAC[prove (`!(m:num) n. (n+m-n) = (m-n)+n`,REWRITE_TAC[ADD_AC])]))
 448 THEN ((UNDISCH_EL_TAC 0))
 449 THEN ((MATCH_ACCEPT_TAC(GSYM SUB_ADD)));
 450 (* branch *)
 451 ((DISCH_TAC))
 452 THEN ((ASM_REWRITE_TAC[FLOAT_ADD]))]);;
 453
 454 let FLOAT_ADD_PPv2 = prove(
 455   `!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))))`,
 456 ((REPEAT GEN_TAC))
 457 THEN (DISCH_TAC)
 458 THEN ((H_MATCH_MP (THM (prove(`!m n. m<|n ==> m <=|n`,MESON_TAC[LT_LE]))) (HYP_INT 0)))
 459 THEN ((UNDISCH_EL_TAC 0))
 460 THEN ((SIMP_TAC[GSYM FLOAT_ADD_PP]))
 461 THEN (DISCH_TAC)
 462 THEN ((REWRITE_TAC[REAL_ADD_AC])));;
 463
 464 let FLOAT_ADD_NN = prove(
 465 `!a b m n. ((n <=| m) ==>( (float a (--:(int_of_num  m))) +. (float b (--:(int_of_num  n)))
 466      = (float  ((int_of_num (2 EXP (m -| n))) *b + a) (--:(int_of_num  m)))))`,
 467 ((REPEAT GEN_TAC))
 468 THEN (DISCH_TAC)
 469 THEN ((SUBGOAL_MP_TAC `--: (int_of_num  n) = --: (int_of_num  m) + (int_of_num  (m-n))`))
 470 THENL [((UNDISCH_EL_TAC 0))
 471 THEN ((SIMP_TAC [INT_OF_REAL_THM (GSYM REAL_OF_NUM_SUB)]))
 472 THEN (DISCH_TAC)
 473 THEN ((INT_ARITH_TAC));
 474 (*branch*)
 475 ((DISCH_TAC))
 476 THEN (ASM_REWRITE_TAC[GSYM FLOAT_ADD;REAL_ADD_AC])]);;
 477
 478 let FLOAT_ADD_NNv2 = prove(
 479 `!a b m n. ((m <| n) ==>( (float a (--:(int_of_num  m))) +. (float b (--:(int_of_num  n)))
 480      = (float  ((int_of_num (2 EXP (n -| m))) *a + b) (--:(int_of_num  n)))))`,
 481 ((REPEAT GEN_TAC))
 482 THEN (DISCH_TAC)
 483 THEN (((H_MATCH_MP (THM (prove(`!m n. m<|n ==> m <=|n`,MESON_TAC[LT_LE]))) (HYP_INT 0))))
 484 THEN (((UNDISCH_EL_TAC 0)))
 485 THEN (((SIMP_TAC[GSYM FLOAT_ADD_NN])))
 486 THEN ((DISCH_TAC))
 487 THEN (((REWRITE_TAC[REAL_ADD_AC]))));;
 488
 489 let FLOAT_SUB = prove(
 490   `!a b n m. (float a n) -. (float b m)
 491      = (float a n) +. (float (--: b) m)`,
 492 REWRITE_TAC[float;int_neg_th;real_sub;REAL_NEG_LMUL]);;
 493
 494 let FLOAT_ABS = prove(
 495   `!a n. real_abs  (float a n) = (float (||: a) n)`,
 496 (REWRITE_TAC[float;int_abs_th;REAL_ABS_MUL;TWOPOW_ABS]));;
 497
 498
 499 let FLOAT_POW = prove(
 500   `!a n m. (float a n) **. m = (float (a **: m) (n *: (int_of_num m)))`,
 501 (REWRITE_TAC[float;REAL_POW_MUL;int_pow_th;TWOPOW_POW]));;
 502
 503 let INT_SUB = prove(
 504   `!a b. (a -: b) = (a +: (--: b))`,
 505  (REWRITE_TAC[GSYM INT_SUB_RNEG;INT_NEG_NEG]));;
 506
 507 let INT_ABS_NUM = prove(
 508   `!n. ||: (int_of_num  n) = (int_of_num  n)`,
 509  (REWRITE_TAC[int_eq;int_abs_th;INT_NUM_REAL;REAL_ABS_NUM]));;
 510
 511 let INT_ABS_NEG_NUM = prove(
 512   `!n. ||: (--: (int_of_num  n)) = (int_of_num  n)`,
 513  (REWRITE_TAC[int_eq;int_abs_th;int_neg_th;INT_NUM_REAL;REAL_ABS_NUM;REAL_ABS_NEG]));;
 514
 515 let INT_ADD_NEG_NUM = prove(`!x y. --: (int_of_num  x) +: (int_of_num  y) = (int_of_num  y) +: (--: (int_of_num  x))`,
 516  (REWRITE_TAC[INT_ADD_AC]));;
 517
 518 let INT_POW_MUL = INT_OF_REAL_THM REAL_POW_MUL;;
 519
 520 let INT_POW_NEG1 = prove (
 521   `!x n. (--: (int_of_num  x)) **: n = ((--: (int_of_num  1)) **: n) *: ((int_of_num  x) **: n)`,
 522 (REWRITE_TAC[GSYM INT_POW_MUL; GSYM INT_NEG_MINUS1]));;
 523
 524 let pow = real_pow;;
 525
 526
 527 let POW_1 = prove(
 528   `!x. x pow 1 = x`,
 529   GEN_TAC THEN REWRITE_TAC[num_CONV `1`] THEN
 530   REWRITE_TAC[pow; REAL_MUL_RID]);;
 531
 532 let POW_2 = prove(
 533   `!x. x pow 2 = x * x`,
 534   GEN_TAC THEN REWRITE_TAC[num_CONV `2`] THEN
 535   REWRITE_TAC[pow; POW_1]);;
 536
 537
 538 let INT_POW_EVEN_NEG1 = prove(
 539   `!x n. (--: (int_of_num  x)) **: (NUMERAL (BIT0 n)) =  ((int_of_num  x) **: (NUMERAL (BIT0 n)))`,
 540 ((REPEAT GEN_TAC))
 541 THEN ((ONCE_REWRITE_TAC[INT_POW_NEG1]))
 542 THEN ((ABBREV_TAC `a = int_of_num  1`))
 543 THEN ((ABBREV_TAC `b = (int_of_num  x)**: (NUMERAL (BIT0 n))`))
 544 THEN ((REWRITE_TAC[NUMERAL;BIT0]))
 545 THEN ((REWRITE_TAC[GSYM MULT_2;GSYM INT_POW_POW;INT_OF_REAL_THM REAL_POW_2;INT_NEG_MUL2]))
 546 THEN ((EXPAND_TAC "a"))
 547 THEN ((REWRITE_TAC[INT_MUL_RID;INT_MUL_LID;INT_OF_REAL_THM REAL_POW_ONE])));;
 548
 549
 550
 551 let POW_MINUS1 = prove(
 552   `!n. (--(&1)) pow (2 * n) = &1`,
 553   INDUCT_TAC THEN REWRITE_TAC[MULT_CLAUSES; pow] THEN
 554   REWRITE_TAC[num_CONV `2`; num_CONV `1`; ADD_CLAUSES] THEN
 555   REWRITE_TAC[pow] THEN
 556   REWRITE_TAC[SYM(num_CONV `2`); SYM(num_CONV `1`)] THEN
 557   ASM_REWRITE_TAC[] THEN
 558   REWRITE_TAC[GSYM REAL_NEG_LMUL; GSYM REAL_NEG_RMUL] THEN
 559   REWRITE_TAC[REAL_MUL_LID; REAL_NEGNEG]);;
 560
 561
 562 let POW_1 = prove(
 563   `!x. x pow 1 = x`,
 564   GEN_TAC THEN REWRITE_TAC[num_CONV `1`] THEN
 565   REWRITE_TAC[pow; REAL_MUL_RID]);;
 566
 567 let INT_POW_ODD_NEG1 = prove(
 568   `!x n. (--: (int_of_num  x)) **: (NUMERAL (BIT1 n)) = --: ((int_of_num  x) **: (NUMERAL (BIT1 n)))`,
 569 ((REPEAT GEN_TAC))
 570 THEN ((ONCE_REWRITE_TAC[INT_POW_NEG1]))
 571 THEN (((ABBREV_TAC `a = int_of_num  1`)))
 572 THEN (((ABBREV_TAC `b = (int_of_num  x)**: (NUMERAL (BIT1 n))`)))
 573 THEN ((REWRITE_TAC[NUMERAL;BIT1]))
 574 THEN ((ONCE_REWRITE_TAC[ADD1]))
 575 THEN ((EXPAND_TAC "a"))
 576 THEN ((REWRITE_TAC[GSYM MULT_2]))
 577 THEN ((REWRITE_TAC[INT_OF_REAL_THM POW_MINUS1;INT_OF_REAL_THM REAL_POW_ADD]))
 578 THEN ((REWRITE_TAC[INT_OF_REAL_THM POW_1;INT_MUL_LID;INT_MUL_LNEG])));;
 579
 580 (* subtraction of integers *)
 581
 582 let INT_ADD_NEG = prove(
 583  `!x y. (x < y ==> ((int_of_num  x) +: (--: (int_of_num  y)) = --: (int_of_num  (y - x))))`,
 584 ((REPEAT GEN_TAC))
 585 THEN ((DISCH_TAC))
 586 THEN ((SUBGOAL_MP_TAC `int_of_num  (y-x ) = (int_of_num  y) - (int_of_num  x)`))
 587 THENL [((SUBGOAL_MP_TAC `x <=| y`))
 588          (* branch *)
 589          THENL [(((ASM MESON_TAC)[LE_LT]));((SIMP_TAC[GSYM (INT_OF_REAL_THM REAL_OF_NUM_SUB)]))]
 590 (* branch *)
 591 ; ((DISCH_TAC))
 592 THEN ((ASM_REWRITE_TAC[]))
 593 THEN (ACCEPT_TAC(INT_ARITH `int_of_num  x +: --: (int_of_num  y) = --: (int_of_num  y -: int_of_num  x)`))]);;
 594
 595 let INT_ADD_NEGv2 = prove(
 596  `!x y. (y <= x ==> ((int_of_num  x) +: (--: (int_of_num  y)) = (int_of_num  (x - y))))`,
 597  ((REPEAT GEN_TAC))
 598  THEN ((DISCH_TAC))
 599  THEN ((SUBGOAL_MP_TAC `int_of_num  (x - y) = (int_of_num  x) - (int_of_num  y)`))
 600  THENL[
 601   ((UNDISCH_EL_TAC 0)) THEN ((SIMP_TAC[GSYM (INT_OF_REAL_THM REAL_OF_NUM_SUB)]));
 602   ((DISCH_TAC)) THEN ((ASM_REWRITE_TAC[INT_SUB]))
 603      ]
 604 );;
 605
 606
 607 (* assumes a term of the form int_of_num a +: (--: (int_of_num  b))  *)
 608 let INT_SUB_CONV t =
 609     let a,b = dest_binop `int_add` t in
 610   let (_,a) = dest_comb a in
 611   let (_,b) = dest_comb b in
 612   let (_,b) = dest_comb b in
 613   let a = dest_numeral a in
 614   let b = dest_numeral b in
 615   let thm = if  (b <=/ a) then
 616     INT_ADD_NEGv2
 617   else INT_ADD_NEG in
 618   (ARITH_SIMP_CONV[thm]) t;; (*   (SIMP_CONV[thm;ARITH]) t;; *)
 619
 620 let FLOAT_EQ = prove(
 621   `!a b a' b'.  (float a b = (float a'  b')) <=>
 622         ((float a b) -. (float a' b') = (&.0))`,MESON_TAC[REAL_SUB_0]);;
 623
 624 let FLOAT_LT = prove(
 625   `!a b a' b'. (float a b <. (float a' b')) <=>
 626         ((&.0) <. (float a' b') -. (float a b))`,MESON_TAC[REAL_SUB_LT]);;
 627
 628 let FLOAT_LE = prove(
 629   `!a b a' b'. (float a b <=. (float a' b')) <=>
 630         ((&.0) <=. (float a' b') -. (float a b))`,MESON_TAC[REAL_SUB_LE]);;
 631
 632 (* ------------------------------------------------------------------ *)
 633 (*   Simplifies an arithmetic expression in floats                    *)
 634 (*   A workhorse                                                      *)
 635 (* ------------------------------------------------------------------ *)
 636
 637 let FLOAT_CONV =
 638               (ARITH_SIMP_CONV[FLOAT_LT;FLOAT_LE;FLOAT_MUL;FLOAT_SUB;FLOAT_ABS;FLOAT_POW;
 639               FLOAT_ADD_NN;FLOAT_ADD_NNv2;FLOAT_ADD_PP;FLOAT_ADD_PPv2;
 640               FLOAT_ADD_NP;FLOAT_ADD_PN;FLOAT_NEG;
 641               INT_NEG_NEG;INT_SUB;
 642               INT_ABS_NUM;INT_ABS_NEG_NUM;
 643               INT_MUL_LNEG;INT_MUL_RNEG;INT_NEG_MUL2;INT_OF_NUM_MUL;
 644               INT_OF_NUM_ADD;GSYM INT_NEG_ADD;INT_ADD_NEG_NUM;
 645               INT_OF_NUM_POW;INT_POW_ODD_NEG1;INT_POW_EVEN_NEG1;
 646               INT_ADD_NEG;INT_ADD_NEGv2 (* ; ARITH *)
 647               ]) ;;
 648
 649 add_test("FLOAT_CONV1",
 650   let f z =
 651     let (x,y) =  dest_eq z in
 652     let (u,v) =  dest_thm (FLOAT_CONV x) in
 653     (u=[]) && (z = v) in
 654   f `float (int_of_num 3) (int_of_num 0) = float (int_of_num 3) (int_of_num 0)` &&
 655   f `float (int_of_num 3) (int_of_num 3) = float (int_of_num 3) (int_of_num 3)` &&
 656   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))` &&
 657   f `float (int_of_num 1 + (int_of_num 3)) (int_of_num 4) = float (int_of_num 4) (int_of_num 4)` &&
 658   f `float (int_of_num 3 - (int_of_num 7)) (int_of_num 0) = float (--:(int_of_num 4)) (int_of_num 0)` &&
 659   f `float (int_of_num 3) (int_of_num 4) *. (float (--: (int_of_num 2)) (int_of_num 3)) = float (--: (int_of_num 6))
 660                                                         (int_of_num 7)` &&
 661   f `real_neg  (float (--: (int_of_num 3)) (int_of_num 0)) = float (int_of_num 3) (int_of_num 0)`
 662         );;
 663
 664 (* ------------------------------------------------------------------ *)
 665 (*   Operations on interval_eps. Preliminary stuff to deal with           *)
 666 (*   chains of inequalities.                                          *)
 667 (* ------------------------------------------------------------------ *)
 668
 669
 670 let REAL_ADD_LE_SUBST_RHS = prove(
 671   `!a b c P. ((a <=. ((P b)) /\ (!x. (P x) =  x + (P (&. 0))) /\ (b <=. c)) ==> (a <=. (P c)))`,
 672 (((REPEAT GEN_TAC)))
 673 THEN (((REPEAT (TAUT_TAC `(b ==> a==> c) ==> (a /\ b ==> c)`))))
 674 THEN (((REPEAT DISCH_TAC)))
 675 THEN ((((H_RULER(ONCE_REWRITE_RULE))[HYP_INT 1] (HYP_INT 0))))
 676 THEN ((((ASM ONCE_REWRITE_TAC)[])))
 677 THEN ((((ASM MESON_TAC)[REAL_LE_RADD;REAL_LE_TRANS]))));;
 678
 679 let REAL_ADD_LE_SUBST_LHS = prove(
 680   `!a b c P. (((P(a) <=. b /\ (!x. (P x) =  x + (P (&. 0))) /\ (c <=. a)))
 681      ==> ((P c) <=. b))`,
 682 (REP_GEN_TAC)
 683 THEN (DISCH_ALL_TAC)
 684 THEN (((H_RULER(ONCE_REWRITE_RULE)) [HYP_INT 1] (HYP_INT 0)))
 685 THEN (((ASM ONCE_REWRITE_TAC)[]))
 686 THEN (((ASM MESON_TAC)[REAL_LE_RADD;REAL_LE_TRANS])));;
 687 (*
 688 let rec SPECL =
 689     function [] -> I |
 690     (a::b)  -> fun thm ->(SPECL b (SPEC a thm));;
 691 *)
 692 (*
 693   input:
 694     rel: b <=. c
 695     thm: a <=. P(b).
 696
 697   output: a <=. P(c).
 698
 699   condition: REAL_ARITH must be able to prove !x. P(x) = x+. P(&.0).
 700   condition: the term `a` must appear exactly once the lhs of the thm.
 701   *)
 702
 703 let IWRITE_REAL_LE_RHS rel thm =
 704   let bvar = genvar `:real` in
 705   let (bt,_) = dest_binop `real_le` (concl rel) in
 706   let sub = SUBS_CONV[ASSUME (mk_eq(bt,bvar))] in
 707   let rule = (fun th -> EQ_MP (sub (concl th)) th) in
 708   let (subrel,subthm) = (rule rel,rule thm) in
 709   let (a,p) = dest_binop `real_le` (concl subthm) in
 710   let (_,c) = dest_binop `real_le` (concl subrel) in
 711   let pfn = mk_abs (bvar,p) in
 712   let imp_th = BETA_RULE (SPECL [a;bvar;c;pfn] REAL_ADD_LE_SUBST_RHS) in
 713   let ppart =   REAL_ARITH
 714       (fst(dest_conj(snd(dest_conj(fst(dest_imp(concl imp_th))))))) in
 715   let prethm = MATCH_MP imp_th (CONJ subthm (CONJ ppart subrel)) in
 716   let prethm2 = SPEC bt (GEN bvar (DISCH
 717        (find (fun x -> try(bvar = rhs x) with Failure _ -> false) (hyp prethm)) prethm)) in
 718   MATCH_MP prethm2 (REFL bt);;
 719
 720 (*
 721   input:
 722     rel: c <=. a
 723     thm: P a <=. b
 724
 725   output: P c <=. b
 726
 727   condition: REAL_ARITH must be able to prove !x. P(x) = x+. P(&.0).
 728   condition: the term `a` must appear exactly once the lhs of the thm.
 729   *)
 730
 731 let IWRITE_REAL_LE_LHS rel thm =
 732   let avar = genvar `:real` in
 733   let (_,at) = dest_binop `real_le` (concl rel) in
 734   let sub = SUBS_CONV[ASSUME (mk_eq(at,avar))] in
 735   let rule = (fun th -> EQ_MP (sub (concl th)) th) in
 736   let (subrel,subthm) = (rule rel,rule thm) in
 737   let (p,b) = dest_binop `real_le` (concl subthm) in
 738   let (c,_) = dest_binop `real_le` (concl subrel) in
 739   let pfn = mk_abs (avar,p) in
 740   let imp_th = BETA_RULE (SPECL [avar;b;c;pfn] REAL_ADD_LE_SUBST_LHS) in
 741   let ppart =   REAL_ARITH
 742       (fst(dest_conj(snd(dest_conj(fst(dest_imp(concl imp_th))))))) in
 743   let prethm = MATCH_MP imp_th (CONJ subthm (CONJ ppart subrel)) in
 744   let prethm2 = SPEC at (GEN avar (DISCH
 745        (find (fun x -> try(avar = rhs x) with Failure _ -> false) (hyp prethm)) prethm)) in
 746   MATCH_MP prethm2 (REFL at);;
 747
 748 (* ------------------------------------------------------------------ *)
 749 (*   INTERVAL ADD, NEG, SUBTRACT                                      *)
 750 (* ------------------------------------------------------------------ *)
 751
 752
 753
 754 let INTERVAL_ADD = prove(
 755    `!x f ex y g ey. interval_eps x f ex /\ interval_eps y g ey ==>
 756                          interval_eps (x +. y) (f +. g) (ex +. ey)`,
 757 EVERY[
 758  REPEAT GEN_TAC;
 759  TAUT_TAC `(A==>B==>C)==>(A/\ B ==> C)`;
 760  REWRITE_TAC[interval_eps];
 761  REWRITE_TAC[prove(`(x+.y) -. (f+.g) = (x-.f) +. (y-.g)`,REAL_ARITH_TAC)];
 762  ABBREV_TAC `a = x-.f`;
 763  ABBREV_TAC `b = y-.g`;
 764  ASSUME_TAC (SPEC `b:real` (SPEC `a:real` ABS_TRIANGLE));
 765  UNDISCH_EL_TAC 0;
 766  ABBREV_TAC `a':real = abs a`;
 767  ABBREV_TAC `b':real = abs b`;
 768  REPEAT DISCH_TAC;
 769  (H_VAL2(IWRITE_REAL_LE_RHS)) (HYP_INT 0) (HYP_INT 2);
 770  (H_VAL2(IWRITE_REAL_LE_RHS)) (HYP_INT 2) (HYP_INT 0);
 771  ASM_REWRITE_TAC[]]);;
 772
 773 let INTERVAL_NEG = prove(
 774   `!x f ex. interval_eps x f ex = interval_eps (real_neg  x) (real_neg  f) ex`,
 775 (REWRITE_TAC[interval_eps;REAL_ABS_NEG;REAL_ARITH `!x y. -- x -. (-- y) = real_neg  (x -. y)`]));;
 776
 777 let INTERVAL_NEG2 = prove(
 778   `!x f ex. interval_eps (real_neg  x) f ex = interval_eps x (real_neg  f) ex`,
 779  (REWRITE_TAC[interval_eps;REAL_ABS_NEG;REAL_ARITH `!x y. -- x -. y = real_neg  (x -. (real_neg  y))`]));;
 780
 781
 782 let INTERVAL_SUB = prove(
 783    `!x f ex y g ey. interval_eps x f ex /\ interval_eps y g ey ==> interval_eps (x -. y) (f -. g) (ex +. ey)`,
 784 ((REWRITE_TAC[real_sub]))
 785 THEN (DISCH_ALL_TAC)
 786 THEN (((H_RULER (ONCE_REWRITE_RULE))[THM(INTERVAL_NEG)] (HYP_INT 1)))
 787 THEN (((ASM MESON_TAC)[INTERVAL_ADD])));;
 788
 789 (* ------------------------------------------------------------------ *)
 790 (*   INTERVAL MULTIPLICATION                                          *)
 791 (* ------------------------------------------------------------------ *)
 792
 793
 794
 795 (* renamed from REAL_LE_ABS_RMUL *)
 796 let REAL_PROP_LE_RABS = prove(
 797   `!x y z. (y <=. z) ==> ( y * (abs x) <=. z *(abs x))`,(SIMP_TAC[REAL_LE_RMUL_IMP;ABS_POS]));;
 798
 799 let REAL_LE_ABS_MUL = prove(
 800   `!x y z w. (( x <=. y) /\ (abs z <=. w)) ==> (x*.w <=. y*.w) `,
 801 (DISCH_ALL_TAC)
 802 THEN ((ASSUME_TAC (REAL_ARITH `abs z <=. w ==> (&.0) <=. w`)))
 803 THEN (((ASM MESON_TAC)[REAL_LE_RMUL_IMP])));;
 804
 805 let ABS_MUL = REAL_ABS_MUL;;
 806
 807 let INTERVAL_ABS = REAL_ARITH
 808   `!(x:real) z d. (x - d) <= z /\ z <= (x + d) <=> abs(z - x) <= d`;;
 809
 810 let INTERVAL_MUL = prove(
 811   `!x f ex y g ey. (interval_eps x f ex) /\ (interval_eps y g ey) ==>
 812          (interval_eps (x *. y) (f *. g) (abs(f)*.ey +. ex*. abs(g) +. ex*.ey))`,
 813 (REP_GEN_TAC)
 814 THEN ((REWRITE_TAC[interval_eps]))
 815 THEN ((REWRITE_TAC[REAL_ARITH `(x*. y -. f*. g) = (f *.(y -. g) +. (x -. f)*.g +. (x-.f)*.(y-. g))`]))
 816 THEN (DISCH_ALL_TAC)
 817 THEN ((ASSUME_TAC (SPECL [`f*.(y-g)`;`(x-f)*g +. (x-f)*.(y-g)`] ABS_TRIANGLE)))
 818 THEN ((ASSUME_TAC (SPECL [`(x-f)*.g`;`(x-f)*.(y-g)`] ABS_TRIANGLE)))
 819 THEN (((H_VAL2(IWRITE_REAL_LE_RHS)) (HYP_INT 0) (HYP_INT 1)))
 820 THEN ((H_REWRITE_RULE [THM ABS_MUL] (HYP_INT 0)))
 821 THEN ((H_MATCH_MP (THM (SPECL [`g:real`; `abs (x -. f)`; `ex:real`] REAL_PROP_LE_RABS)) (HYP_INT 4)))
 822 THEN (((H_VAL2(IWRITE_REAL_LE_RHS)) (HYP_INT 0) (HYP_INT 1)))
 823 THEN ((H_MATCH_MP (THM (SPECL [`f:real`; `abs (y -. g)`; `ey:real`] REAL_PROP_LE_LABS)) (HYP_INT 7)))
 824 THEN (((H_VAL2 (IWRITE_REAL_LE_RHS)) (HYP_INT 0) (HYP_INT 1)))
 825 THEN ((H_MATCH_MP (THM (SPECL [`x-.f`; `abs (y -. g)`; `ey:real`] REAL_PROP_LE_LABS)) (HYP_INT 9)))
 826 THEN (((H_VAL2(IWRITE_REAL_LE_RHS)) (HYP_INT 0) (HYP_INT 1)))
 827 THEN ((ASSUME_TAC (SPECL [`abs(x-.f)`;`ex:real`;`y-.g`;`ey:real`] REAL_LE_ABS_MUL)))
 828 THEN ((H_CONJ (HYP_INT 11) (HYP_INT 12)))
 829 THEN ((H_MATCH_MP (HYP_INT 1) (HYP_INT 0)))
 830 THEN (((H_VAL2(IWRITE_REAL_LE_RHS)) (HYP_INT 0) (HYP_INT 3)))
 831 THEN ((POP_ASSUM ACCEPT_TAC)));;
 832
 833 (* ------------------------------------------------------------------ *)
 834 (*   INTERVAL BASIC OPERATIONS                                        *)
 835 (* ------------------------------------------------------------------ *)
 836
 837
 838 let INTERVAL_NUM = prove(
 839   `!n. (interval_eps(&.n) (float(int_of_num n) (int_of_num 0)) (float  (int_of_num  0) (int_of_num 0)))`,
 840 (REWRITE_TAC[interval_eps;float;TWOPOW_POS;pow;REAL_MUL_RID;INT_NUM_REAL;REAL_SUB_REFL;REAL_ABS_0;REAL_LE_REFL]));;
 841
 842 let INTERVAL_CENTER = prove(
 843   `!x f ex g. (interval_eps x f ex) ==> (interval_eps x g (abs(f-g)+.ex))`,
 844 ((REWRITE_TAC[interval_eps]))
 845 THEN (DISCH_ALL_TAC)
 846 THEN ((ASSUME_TAC (REAL_ARITH `abs(x -. g) <=. abs(f-.g) +. abs(x -. f)`)))
 847 THEN ((H_VAL2 IWRITE_REAL_LE_RHS (HYP_INT 1) (HYP_INT 0)))
 848 THEN ((ASM_REWRITE_TAC[])));;
 849
 850 let INTERVAL_WIDTH = prove(
 851   `!x f ex ex'. (ex <=. ex') ==> (interval_eps x f ex) ==> (interval_eps x f ex')`,
 852 ((REWRITE_TAC[interval_eps]))
 853 THEN (DISCH_ALL_TAC)
 854 THEN ((H_VAL2 IWRITE_REAL_LE_RHS (HYP_INT 1) (HYP_INT 0)))
 855 THEN ((ASM_REWRITE_TAC[])));;
 856
 857 let INTERVAL_MAX = prove(
 858   `!x f ex. interval_eps x f ex ==> (x <=. f+.ex)`,
 859 (REWRITE_TAC[interval_eps]) THEN REAL_ARITH_TAC);;
 860
 861 let INTERVAL_MIN = prove(
 862   `!x f ex. interval_eps x f ex ==> (f-. ex <=. x)`,
 863 (REWRITE_TAC[interval_eps]) THEN REAL_ARITH_TAC);;
 864
 865 let INTERVAL_ABS_MIN = prove(
 866   `!x f ex. interval_eps x f ex ==> (abs(f)-. ex <=. abs(x))`,
 867   (REWRITE_TAC[interval_eps] THEN REAL_ARITH_TAC)
 868 );;
 869
 870 let INTERVAL_ABS_MAX = prove(
 871   `!x f ex. interval_eps x f ex ==> (abs(x) <=. abs(f)+. ex)`,
 872   (REWRITE_TAC[interval_eps] THEN REAL_ARITH_TAC)
 873 );;
 874
 875 let REAL_RINV_2 = prove(
 876   `&.2 *. (inv (&.2 )) = &. 1`,
 877 EVERY[
 878   MATCH_MP_TAC REAL_MUL_RINV;
 879   REAL_ARITH_TAC]);;
 880
 881 let INTERVAL_MK = prove(
 882    `let half = float(int_of_num 1)(--:(int_of_num 1)) in
 883     !x xmin xmax. ((xmin <=. x) /\ (x <=. xmax)) ==>
 884       interval_eps x
 885          ((xmin+.xmax)*.half)
 886          ((xmax-.xmin)*.half)`,
 887 EVERY[
 888   REWRITE_TAC[LET_DEF;LET_END_DEF];
 889   DISCH_ALL_TAC;
 890   REWRITE_TAC[interval_eps;float;TWOPOW_NEG;INT_NUM_REAL;REAL_POW_1;REAL_MUL_LID];
 891   REWRITE_TAC[GSYM INTERVAL_ABS];
 892   CONJ_TAC
 893   ]
 894 THENL[
 895   EVERY[
 896     REWRITE_TAC[GSYM REAL_SUB_RDISTRIB];
 897     REWRITE_TAC[REAL_ARITH `(b+.a)-.(a-.b)=b*.(&.2)`;GSYM REAL_MUL_ASSOC];
 898     ASM_REWRITE_TAC[REAL_RINV_2;REAL_MUL_RID]
 899   ];
 900   EVERY[
 901     REWRITE_TAC[GSYM REAL_ADD_RDISTRIB];
 902     REWRITE_TAC[REAL_ARITH `(b+.a)+. a -. b=a*.(&.2)`;GSYM REAL_MUL_ASSOC];
 903     ASM_REWRITE_TAC[REAL_RINV_2;REAL_MUL_RID]
 904   ]
 905 ]);;
 906
 907 let INTERVAL_EPS_POS = prove(`!x f ex.
 908   (interval_eps x f ex) ==> (&.0 <=. ex)`,
 909 EVERY[
 910   REWRITE_TAC[interval_eps];
 911   REPEAT (GEN_TAC);
 912   DISCH_THEN(fun x -> (MP_TAC (CONJ (SPEC `x-.f` REAL_ABS_POS) x)));
 913   MATCH_ACCEPT_TAC REAL_LE_TRANS]);;
 914
 915 let INTERVAL_EPS_0 = prove(
 916   `!x f n. (interval_eps x f (float (int_of_num 0) n)) ==> (x = f)`,
 917 EVERY[
 918   REWRITE_TAC[interval_eps;float;int_of_num_th;REAL_MUL_LZERO];
 919   REAL_ARITH_TAC]);;
 920
 921 let REAL_EQ_RCANCEL_IMP' = prove(`!x y z.(x * z = y * z) ==> (~(z = &0) ==> (x=y))`,
 922   MESON_TAC[REAL_EQ_RCANCEL_IMP]);;
 923
 924 (* renamed from REAL_ABS_POS *)
 925 let REAL_MK_POS_ABS_' = prove (`!x. (~(x=(&.0))) ==> (&.0 < abs(x))`,
 926   MESON_TAC[REAL_PROP_NZ_ABS;ABS_POS;REAL_LT_LE]);;
 927
 928 (* ------------------------------------------------------------------ *)
 929 (*   INTERVAL DIVIDE                                                  *)
 930 (* ------------------------------------------------------------------ *)
 931
 932 let INTERVAL_DIV = prove(`!x f ex y g ey h ez.
 933   (((interval_eps x f ex) /\ (interval_eps y g ey) /\ (ey <. (abs g)) /\
 934   ((ex +. (abs (f -. (h*.g))) +. (abs h)*. ey) <=. (ez*.((abs g) -. ey))))
 935   ==> (interval_eps (x / y) h ez))`,
 936 let lemma1 = prove( `&.0 < u /\ real_abs  z <=. e*. u ==> (&.0) <=. e`,
 937   EVERY[
 938     DISCH_ALL_TAC;
 939     ASSUME_TAC (SPEC `z:real` REAL_MK_NN_ABS);
 940     H_MATCH_MP (THM REAL_LE_TRANS) (H_RULE2 CONJ (HYP_INT 0) (HYP_INT 2));
 941     H_MATCH_MP (THM REAL_PROP_NN_RCANCEL) (H_RULE2 CONJ (HYP_INT 2) (HYP_INT 0));
 942     ASM_REWRITE_TAC[]
 943   ]) in
 944 EVERY[
 945   DISCH_ALL_TAC;
 946   SUBGOAL_MP_TAC `~(y= (&.0))`
 947   THENL[
 948     EVERY[
 949       UNDISCH_LIST[1;2];
 950       REWRITE_TAC[interval_eps];
 951       REAL_ARITH_TAC
 952     ];
 953     EVERY[
 954       REWRITE_TAC[interval_eps];
 955       DISCH_TAC THEN (H I (HYP_INT 0)) THEN (UNDISCH_EL_TAC 0);
 956       DISCH_THEN (fun th -> (MP_TAC(MATCH_MP REAL_MK_POS_ABS_' th)));
 957       MATCH_MP_TAC REAL_MUL_RTIMES_LE;
 958       REWRITE_TAC[GSYM ABS_MUL;REAL_SUB_RDISTRIB;real_div;GSYM REAL_MUL_ASSOC];
 959       ASM_SIMP_TAC[REAL_MUL_LINV;REAL_MUL_RID];
 960       H (REWRITE_RULE[interval_eps]) (HYP_INT 1);
 961       H (REWRITE_RULE[interval_eps]) (HYP_INT 3);
 962       H (MATCH_MP INTERVAL_ABS_MIN) (HYP_INT 4);
 963       POPL_TAC[3;4;5];
 964       H_VAL2 (IWRITE_REAL_LE_LHS) (HYP_INT 2) (HYP_INT 4);
 965       H (REWRITE_RULE[ REAL_ADD_ASSOC]) (HYP_INT 0);
 966       H_VAL2 (IWRITE_REAL_LE_LHS) (THM (SPEC `f-. h*g` (SPEC `x-.f` ABS_TRIANGLE))) (HYP_INT 0);
 967       H (ONCE_REWRITE_RULE[REAL_ABS_SUB]) (HYP_INT 4);
 968       H (MATCH_MP (SPEC `h:real` REAL_PROP_LE_LABS)) (HYP_INT 0);
 969       H (REWRITE_RULE[GSYM ABS_MUL]) (HYP_INT 0);
 970       H_VAL2 (IWRITE_REAL_LE_LHS) (HYP_INT 0) (HYP_INT 3);
 971       H_VAL2 (IWRITE_REAL_LE_LHS) (THM (SPEC `h*.(g-.y)` (SPEC`(x-.f)+(f-. h*g)`  ABS_TRIANGLE))) (HYP_INT 0);
 972       POPL_TAC[1;2;3;4;5;6;7;9;10;12];
 973       H (ONCE_REWRITE_RULE[REAL_ARITH `((x-.f) +. (f -. h*. g)) +. h*.(g-. y) = x -. h*. y `]) (HYP_INT 0);
 974       ABBREV_TAC `z = x -. h*.y`;
 975       H (ONCE_REWRITE_RULE[GSYM REAL_SUB_LT]) (HYP_INT 4);
 976       ABBREV_TAC `u = abs(g) -. ey`;
 977       POPL_TAC[0;2;4;6];
 978       H (MATCH_MP lemma1 ) (H_RULE2 CONJ (HYP_INT 0) (HYP_INT 1));
 979       H (MATCH_MP REAL_PROP_LE_LMUL) (H_RULE2 CONJ (HYP_INT 0) (HYP_INT 3));
 980       H (MATCH_MP REAL_LE_TRANS) (H_RULE2 CONJ (HYP_INT 3) (HYP_INT 0));
 981       ASM_REWRITE_TAC[]
 982   ];
 983   ]]);;
 984
 985 (* ------------------------------------------------------------------ *)
 986 (*   INTERVAL ABS VALUE                                               *)
 987 (* ------------------------------------------------------------------ *)
 988
 989 let INTERVAL_ABSV = prove(`!x f ex. interval_eps x f ex ==> (interval_eps (abs x) (abs f) ex)`,
 990 EVERY[
 991   REWRITE_TAC[interval_eps];
 992   DISCH_ALL_TAC;
 993   ASSUME_TAC (SPECL [`x:real`;`f:real`] REAL_ABS_SUB_ABS);
 994   ASM_MESON_TAC[REAL_LE_TRANS]
 995 ]);;  (* 7 minutes *)
 996
 997 (* ------------------------------------------------------------------ *)
 998 (*   INTERVAL ATN VALUE                                               *)
 999 (* ------------------------------------------------------------------ *)
1000
1001 let taylor_atn = prove_by_refinement(
1002   `!n x v eps2 eps.
1003       interval_eps (&16 * sum (0..n) (halfatn4_co x)) v eps2 /\
1004       float (int_of_num 1) ( -- int_of_num (6 * n + 5)) + eps2 <= eps ==>
1005       interval_eps (atn x) v eps`,
1006   (* {{{ proof *)
1007   [
1008   REWRITE_TAC[interval_eps;float;TWOPOW_NEG;int_of_num_th;REAL_ARITH `&1 * x = x`];
1009   REPEAT STRIP_TAC;
1010   MATCH_MP_TAC real_taylor_atn_approx;
1011   EXISTS_TAC `n:num`;
1012   EXISTS_TAC `inv(&2 pow (6 * n +5))`;
1013   EXISTS_TAC `eps2:real`;
1014   ASM_REWRITE_TAC[];
1015   REAL_ARITH_TAC;
1016   ]);;
1017   (* }}} *)
1018
1019 let pi_atn = prove_by_refinement(
1020   `pi = &4 * atn (&1)`,
1021   (* {{{ proof *)
1022   [
1023   REWRITE_TAC[ATN_1];
1024   REAL_ARITH_TAC;
1025   ]);;
1026   (* }}} *)
1027
1028
1029 (* ------------------------------------------------------------------ *)
1030 (*   INTERVAL SQRT                                                    *)
1031 (*   This requires some preliminaries. Extend sqrt by 0 on negatives  *)
1032 (* ------------------------------------------------------------------ *)
1033
1034 let ssqrt = new_definition `ssqrt x = if (x <. (&.0)) then (&.0) else sqrt x`;; (*2m*)
1035
1036 let ssqrt_sqrt = prove_by_refinement(
1037   `!x u f. &0 < x /\ interval_eps (ssqrt x) u f ==>
1038         interval_eps (sqrt x) u f`,
1039   (* {{{ proof *)
1040   [
1041   REWRITE_TAC[ssqrt];
1042   MESON_TAC[REAL_ARITH `&0 < x ==> ~(x < &0)`];
1043   ]);;
1044   (* }}} *)
1045
1046 let LET_TAC = REWRITE_TAC[LET_DEF;LET_END_DEF];;
1047
1048
1049 let REAL_SSQRT_NEG = prove(`!x. (x <. (&.0)) ==> (ssqrt x = (&.0))`,
1050   EVERY[
1051     DISCH_ALL_TAC;
1052     REWRITE_TAC[ssqrt];
1053     COND_CASES_TAC
1054     THENL[
1055       ACCEPT_TAC (REFL `&.0`);
1056       ASM_MESON_TAC[]
1057     ]
1058   ]);;
1059 (* 5 min*)
1060
1061 let REAL_SSQRT_NN = prove(`!x. (&.0) <=. x ==> (ssqrt x = (sqrt x))`,
1062   EVERY[
1063   DISCH_ALL_TAC;
1064   REWRITE_TAC[ssqrt];
1065   COND_CASES_TAC
1066   THENL[
1067     ASM_MESON_TAC[real_lt];
1068     ACCEPT_TAC (REFL `sqrt x`)
1069   ]
1070   ]);;  (* 12 min, mostly spent loading *index-shell* *)
1071
1072
1073 (*17 minutes*)
1074 let REAL_MK_NN_SSQRT = prove(`!x. (&.0) <=. (ssqrt x)`,
1075   EVERY[
1076     GEN_TAC;
1077     DISJ_CASES_TAC (SPECL[`x:real`;`&.0`] REAL_LTE_TOTAL)
1078     THENL[
1079       POP_ASSUM (fun th -> MP_TAC(MATCH_MP (REAL_SSQRT_NEG) th)) THEN
1080       MESON_TAC[REAL_LE_REFL];
1081       POP_ASSUM (fun th -> ASSUME_TAC(CONJ th (MATCH_MP (REAL_SSQRT_NN) th)))  THEN
1082       ASM_MESON_TAC[REAL_PROP_NN_SQRT]
1083     ]
1084   ]);;
1085
1086 let REAL_SV_SSQRT_0  = prove(`!x. ssqrt (&.0) = (&.0)`,
1087   EVERY[
1088     GEN_TAC;
1089     MP_TAC (SPEC `&.0` REAL_LE_REFL);
1090     DISCH_THEN(fun th -> REWRITE_TAC[MATCH_MP REAL_SSQRT_NN th]);
1091     ACCEPT_TAC REAL_SV_SQRT_0
1092   ]);; (* 6 minutes *)
1093
1094
1095 let REAL_SV_SSQRT_n  = prove_by_refinement(`!n. ssqrt (&n) = sqrt(&n)`,
1096 [
1097     GEN_TAC;
1098     MATCH_MP_TAC REAL_SSQRT_NN;
1099     REAL_ARITH_TAC;
1100   ]);;
1101
1102
1103 let REAL_SSQRT_EQ_0 = prove(`!(x:real). (ssqrt(x) = (&.0)) ==> (x <=. (&.0))`,
1104   EVERY[
1105     GEN_TAC;
1106     DISJ_CASES_TAC (SPECL[`x:real`;`&.0`] REAL_LTE_TOTAL)
1107     THENL[
1108       ASM_MESON_TAC[REAL_LT_IMP_LE];
1109       ASM_SIMP_TAC[REAL_SSQRT_NN] THEN
1110       ASM_MESON_TAC[SQRT_EQ_0;REAL_EQ_IMP_LE]
1111     ]
1112   ]);;  (* 15 minutes *)
1113
1114
1115 let REAL_SSQRT_MONO = prove(`!x. (x<=. y) ==> (ssqrt x <=. (ssqrt y))`,
1116   EVERY[
1117     GEN_TAC;
1118     DISJ_CASES_TAC (SPECL[`x:real`;`&.0`] REAL_LTE_TOTAL)
1119       THENL[
1120         ASM_MESON_TAC[REAL_SSQRT_NEG;REAL_MK_NN_SSQRT];
1121         ASM_MESON_TAC[REAL_LE_TRANS;REAL_SSQRT_NN;REAL_PROP_LE_SQRT];
1122       ]
1123   ]);;  (* 5 minutes *)
1124
1125 let REAL_SSQRT_CHAR = prove(`!x t. (&.0 <=. t /\ (t*t = x)) ==> (t = (ssqrt x))`,
1126   EVERY[
1127     DISCH_ALL_TAC;
1128     H_ASSUME_TAC (H_RULE_LIST REWRITE_RULE[HYP_INT 1] (THM (SPEC `t:real` REAL_MK_NN_SQUARE)));
1129     ASM_MESON_TAC[REAL_SSQRT_NN;SQRT_MUL;POW_2_SQRT_ABS;REAL_POW_2;REAL_ABS_REFL];
1130   ]);;  (* 13 minutes *)
1131
1132 let REAL_SSQRT_SQUARE = prove(`!x. (&.0 <=. x) ==> ((ssqrt x)*.(ssqrt x) = x)`,
1133   MESON_TAC[REAL_SSQRT_NN;POW_2;SQRT_POW_2]);;(* 7min *)
1134
1135 let REAL_SSQRT_SQUARE' = prove(`!x. (&.0<=. x) ==> (ssqrt (x*.x) = x)`,
1136   DISCH_ALL_TAC THEN
1137   REWRITE_TAC[(MATCH_MP REAL_SSQRT_NN (SPEC `x:real` REAL_MK_NN_SQUARE))] THEN
1138   ASM_SIMP_TAC[SQRT_MUL;GSYM POW_2;SQRT_POW_2]);; (*20min*)
1139
1140
1141 (* an alternate proof appears in RCS *)
1142 let INTERVAL_SSQRT = prove(`!x f ex u ey ez v. (interval_eps x f ex) /\ (interval_eps (u*.u) f ey) /\
1143   (ex +.ey <=. ez*.(v+.u)) /\ (v*.v <=. f-.ex) /\ (&.0 <. u) /\ (&.0 <=. v) ==>
1144   (interval_eps (ssqrt x) u ez)`,
1145 EVERY[
1146   DISCH_ALL_TAC;
1147   H (MATCH_MP REAL_LE_TRANS) (H_RULE2 CONJ (THM (SPEC `v:real` REAL_MK_NN_SQUARE)) (HYP_INT 3));
1148   H (MATCH_MP (INTERVAL_MIN)) (HYP_INT 1);
1149   H (MATCH_MP REAL_LE_TRANS)  (H_RULE2 CONJ (HYP_INT 1) (HYP_INT 0));
1150   H (MATCH_MP INTERVAL_EPS_POS) (HYP_INT 3);
1151   H (MATCH_MP INTERVAL_EPS_POS) (HYP_INT 5);
1152   H (MATCH_MP REAL_PROP_NN_ADD2) (H_RULE2 CONJ (HYP_INT 1) (HYP_INT 0));
1153   H (MATCH_MP REAL_PROP_POS_LADD) (H_RULE2 CONJ (HYP_INT 11) (HYP_INT 10));
1154   H (MATCH_MP REAL_PROP_POS_LADD) (H_RULE2 CONJ (THM (SPEC `x:real` REAL_MK_NN_SSQRT)) (HYP_INT 11));
1155   H (MATCH_MP REAL_PROP_POS_INV) (HYP_INT 0);
1156   ASSUME_TAC (REAL_ARITH  `(ssqrt x -. u) = (ssqrt x -. u)*.(&.1)`);
1157   H (MATCH_MP REAL_MK_NZ_POS) (HYP_INT 2);
1158   H (MATCH_MP REAL_MUL_RINV) (HYP_INT 0);
1159   H_REWRITE_RULE[(H_RULE GSYM) (HYP_INT 0)] (HYP_INT 2);
1160   POPL_TAC[1;2;3];
1161   H (REWRITE_RULE[REAL_MUL_ASSOC]) (HYP_INT 0);
1162   H (REWRITE_RULE[ONCE_REWRITE_RULE[REAL_MUL_SYM] REAL_DIFFSQ]) (HYP_INT 0);
1163   POPL_TAC[1;2];
1164   H_SIMP_RULE[HYP_INT 7;THM REAL_SSQRT_SQUARE] (HYP_INT 0);
1165   ASSUME_TAC (REAL_ARITH `abs(x -. u*.u) <=. abs(x -. f) + abs(f-. u*.u)`);
1166   H (REWRITE_RULE[interval_eps]) (HYP_INT 12);
1167   H (ONCE_REWRITE_RULE[interval_eps]) (HYP_INT 14);
1168   H (ONCE_REWRITE_RULE[REAL_ABS_SUB]) (HYP_INT 0);
1169   POPL_TAC[1;5;15;16];
1170   H (MATCH_MP REAL_LE_ADD2) (H_RULE2 CONJ (HYP_INT 1) (HYP_INT 0));
1171   H (MATCH_MP REAL_LE_TRANS) (H_RULE2 CONJ (HYP_INT 3) (HYP_INT 0));
1172   POPL_TAC[1;2;3;4];
1173   H (AP_TERM `real_abs `) (HYP_INT 1);
1174   H (REWRITE_RULE[ABS_MUL]) (HYP_INT 0);
1175   H (MATCH_MP REAL_LT_IMP_LE)  (HYP_INT 4);
1176   H (REWRITE_RULE[GSYM REAL_ABS_REFL]) (HYP_INT 0);
1177   H_REWRITE_RULE [HYP_INT 0] (HYP_INT 2);
1178   H (MATCH_MP REAL_LE_RMUL) (H_RULE2 CONJ (HYP_INT 5) (HYP_INT 2));
1179   H_REWRITE_RULE [H_RULE GSYM (HYP_INT 1)] (HYP_INT 0);
1180   POPL_TAC[1;2;3;5;6;7;8];
1181   H (MATCH_MP REAL_LE_TRANS) (H_RULE2 CONJ (HYP_INT 12) (HYP_INT 9));
1182   H (MATCH_MP REAL_SSQRT_MONO) (HYP_INT 0);
1183   H (MATCH_MP REAL_SSQRT_SQUARE') (HYP_INT 16);
1184   H_REWRITE_RULE [HYP_INT 0] (HYP_INT 1);
1185   H (ONCE_REWRITE_RULE[GSYM (SPECL[`v:real`;`ssqrt x`;`u:real`] REAL_LE_RADD)]) (HYP_INT 0);
1186   H (MATCH_MP REAL_LE_INV2) (H_RULE2 CONJ (HYP_INT 9) (HYP_INT 0));
1187   POPL_TAC[1;2;3;4;5;7;8;9;12;13];
1188   H (MATCH_MP REAL_LE_LMUL) (H_RULE2 CONJ (HYP_INT 3) (HYP_INT 0));
1189   H (MATCH_MP REAL_LE_TRANS) (H_RULE2 CONJ (HYP_INT 2) (HYP_INT 0));
1190   H (MATCH_MP REAL_PROP_POS_INV) (HYP_INT 4);
1191   H (MATCH_MP REAL_LT_IMP_LE) (HYP_INT 0);
1192   H (MATCH_MP REAL_LE_RMUL) (H_RULE2 CONJ (HYP_INT 11) (HYP_INT 0));
1193   H (REWRITE_RULE[GSYM REAL_MUL_ASSOC]) (HYP_INT 0);
1194   H (MATCH_MP REAL_MK_NZ_POS) (HYP_INT 8);
1195   H (MATCH_MP REAL_MUL_RINV) (HYP_INT 0);
1196   H_REWRITE_RULE[HYP_INT 0; THM REAL_MUL_RID] (HYP_INT 2);
1197   H (MATCH_MP REAL_LE_TRANS) (H_RULE2 CONJ (HYP_INT 7) (HYP_INT 0));
1198   ASM_REWRITE_TAC[interval_eps]
1199   ]);;
1200
1201
1202
1203 test();;
1204
1205
1206 (* conversion for interval_eps *)
1207
1208 (* ------------------------------------------------------------------ *)
1209 (*   Take a term x of type real.  Convert to a thm of the form        *)
1210 (*   interval_eps x f eps                                                 *)
1211 (*                                                                    *)
1212 (* ------------------------------------------------------------------ *)
1213
1214 let add_test,test = new_test_suite();;
1215
1216 (* ------------------------------------------------------------------ *)
1217 (* num_exponent
1218    Take the absolute value of input.
1219    Write it as a*2^k, where 1 <= a < 2, return k.
1220
1221    Except:
1222    num_exponent (Int 0) is -1.
1223 *)
1224 let (num_exponent:Num.num -> Num.num) =
1225   fun a ->
1226     let afloat = float_of_num (abs_num a) in
1227     Int ((snd (frexp afloat)) - 1);;
1228
1229 (*test*)let f (u,v) = ((num_exponent u) =(Int v)) in
1230     add_test("num_exponenwt",
1231                 forall f
1232     [Int 1,0; Int 65,6; Int (-65),6;
1233      Int 0,-1; (Int 3)//(Int 4),-1]);;
1234 (* ------------------------------------------------------------------ *)
1235
1236 let dest_unary op tm =
1237   try let xop,r = dest_comb tm in
1238       if xop = op then r else fail()
1239   with Failure _ -> failwith "dest_unary";;
1240
1241
1242 (* ------------------------------------------------------------------ *)
1243
1244
1245 (* finds a nearby (outward-rounded) Int with only prec_b significant bits *)
1246 let (round_outward: int -> Num.num -> Num.num) =
1247   fun prec_b a ->
1248     let b = abs_num a in
1249     let sign = if (a =/ b) then I else minus_num in
1250     let throw_bits = Num.max_num (Int 0) ((num_exponent b)-/ (Int prec_b)) in
1251     let twoexp = power_num (Int 2) throw_bits  in
1252     (sign (ceiling_num (b // twoexp)))*/twoexp;;
1253
1254 let (round_inward: int-> Num.num -> Num.num) =
1255   fun prec_b a ->
1256     let b = abs_num a in
1257     let sign = if (a=/b) then I else minus_num in
1258     let throw_bits = Num.max_num (Int 0) ((num_exponent b)-/ (Int prec_b)) in
1259     let twoexp = power_num (Int 2) throw_bits  in
1260     (sign (floor_num (b // twoexp)))*/twoexp;;
1261
1262 let round_rat bprec n =
1263   let b = abs_num n in
1264   let sign = if (b =/ n) then I else minus_num in
1265   let powt  = ((Int 2) **/ (Int bprec)) in
1266   sign ((round_outward bprec (Num.ceiling_num (b */ powt)))//powt);;
1267
1268 let round_inward_rat bprec n =
1269   let b = abs_num n in
1270   let sign = if (b =/ n) then I else minus_num in
1271   let powt  = ((Int 2) **/ (Int bprec)) in
1272   sign ((round_inward bprec (Num.floor_num (b */ powt)))//powt);;
1273
1274 let (round_outward_float: int -> float -> Num.num) =
1275  fun  bprec f ->
1276   if (f=0.0) then (Int 0) else
1277   begin
1278     let b = abs_float f in
1279     let sign = if (f >= 0.0) then I else minus_num in
1280     let (x,n) = frexp b in
1281     let u = int_of_float( ceil (ldexp x bprec)) in
1282     sign ((Int u)*/ ((Int 2) **/ (Int (n - bprec))))
1283   end;;
1284
1285 let (round_inward_float: int -> float -> Num.num) =
1286  fun  bprec f ->
1287   if (f=0.0) then (Int 0) else
1288   begin
1289     (* avoid overflow on 30 bit integers *)
1290     let bprec = if (bprec > 25) then 25 else bprec in
1291     let b = abs_float f in
1292     let sign = if (f >= 0.0) then I else minus_num in
1293     let (x,n) = frexp b in
1294     let u = int_of_float( floor (ldexp x bprec)) in
1295     sign ((Int u)*/ ((Int 2) **/ (Int (n - bprec))))
1296   end;;
1297
1298 (* ------------------------------------------------------------------ *)
1299
1300 (* This doesn't belong here.  A general term substitution function *)
1301 let SUBST_TERM sublist tm =
1302   rhs (concl ((SPECL (map fst sublist)) (GENL (map snd sublist)
1303                                           (REFL tm))));;
1304
1305 add_test("SUBST_TERM",
1306  SUBST_TERM [(`#1`,`a:real`);(`#2`,`b:real`)] (`a +. b +. c`) =
1307  `#1 + #2 + c`);;
1308
1309 (* ------------------------------------------------------------------ *)
1310
1311 (* take a term of the form `interval_eps x f ex` and clean up the f and ex *)
1312
1313 let INTERVAL_CLEAN_CONV:conv =
1314   fun interv ->
1315     let (ixf,ex) = dest_comb interv in
1316     let (ix,f) = dest_comb ixf in
1317     let fthm = FLOAT_CONV f in
1318     let exthm = FLOAT_CONV ex in
1319     let ixfthm = AP_TERM ix fthm in
1320     MK_COMB (ixfthm, exthm);;
1321
1322 (*test*) add_test("INTERVAL_CLEAN_CONV",
1323   let testval = INTERVAL_CLEAN_CONV `interval_eps ((&.1) +. (&.1))
1324        (float (int_of_num 3) (int_of_num 4) +. (float (int_of_num 2) (--: (int_of_num 3))))
1325        (float (int_of_num 1) (int_of_num 2) *. (float (int_of_num 3) (--: (int_of_num 2))))` in
1326   let hypval = hyp testval in
1327   let concval = concl testval in
1328         (length hypval = 0) &&
1329         concval =
1330      `interval_eps (&1 + &1) (float (int_of_num 3) (int_of_num 4) + float (int_of_num 2) (--: (int_of_num 3)))
1331      (float (int_of_num 1) (int_of_num 2) * float (int_of_num 3) (--: (int_of_num 2))) =
1332      interval_eps (&1 + &1) (float (int_of_num 386) (--: (int_of_num 3))) (float (int_of_num 3) (int_of_num 0))`
1333                   );;
1334
1335 (* ------------------------------------------------------------------ *)
1336 (*   GENERAL lemmas                                                   *)
1337 (* ------------------------------------------------------------------ *)
1338
1339
1340 (* verifies statement of the form `float a b = float a' b'` *)
1341
1342
1343
1344 let REAL_INV_POS = REAL_LT_INV;;
1345
1346 let TWOPOW_MK_POS = prove(
1347   `!a. (&.0 <. ( twopow a))`,
1348 EVERY[
1349   GEN_TAC;
1350   CHOOSE_TAC (SPEC `a:int` INT_REP2);
1351   POP_ASSUM DISJ_CASES_TAC;
1352   ASM_REWRITE_TAC[TWOPOW_POS;TWOPOW_NEG];
1353   TRY (MATCH_MP_TAC REAL_INV_POS);
1354   MATCH_MP_TAC REAL_POW_LT ;
1355   REAL_ARITH_TAC;
1356 ]);;
1357
1358 let TWOPOW_NZ = prove(
1359   `!a. ~(twopow a = (&.0))`,
1360   GEN_TAC THEN
1361   ACCEPT_TAC (MATCH_MP REAL_MK_NZ_POS (SPEC `a:int` TWOPOW_MK_POS)));;
1362
1363 let FLOAT_ZERO = prove(
1364   `!a b. (float a b = (&.0)) <=> (a = (int_of_num 0))`,
1365 EVERY[
1366   REWRITE_TAC[float;REAL_ENTIRE;INT_OF_NUM_DEST];
1367   MESON_TAC[TWOPOW_NZ];
1368 ]);;
1369
1370 let INT_ZERO = prove(
1371   `!n. ((int_of_num n = (int_of_num 0)) = (n=0))`,REWRITE_TAC[INT_OF_NUM_EQ]);;
1372
1373 let INT_ZERO_NEG=prove(
1374   `!n. ((--: (int_of_num n) = (int_of_num 0))) <=> (n=0)`,
1375     REWRITE_TAC[INT_NEG_EQ_0;INT_ZERO]);;
1376
1377 let FLOAT_NN = Prove_by_refinement.prove_by_refinement(
1378   `!a b. ((&.0) <=. (float a b)) <=> (int_of_num 0 <=: a)`,
1379 [
1380   REWRITE_TAC[float;INT_OF_NUM_DEST];
1381   REP_GEN_TAC;
1382   REWRITE_TAC[int_of_num_th;int_le];
1383   MESON_TAC[REAL_PROP_NN_RCANCEL;TWOPOW_MK_POS;REAL_PROP_NN_POS;REAL_PROP_NN_MUL2];
1384 ]);;
1385
1386 let INT_NN = INT_POS;;
1387
1388 let INT_NN_NEG = prove(`!n. ((int_of_num 0) <=: (--:(int_of_num n))) <=> (n = 0)`,
1389   REWRITE_TAC[INT_NEG_GE0;INT_OF_NUM_LE] THEN ARITH_TAC
1390                       );;
1391
1392 let FLOAT_POS = prove(`!a b. ((&.0) <. (float a b)) <=> (int_of_num 0 <: a)`,
1393   MESON_TAC[FLOAT_NN;FLOAT_ZERO;INT_LT_LE;REAL_LT_LE]);;
1394
1395 let INT_POS' = prove(`!n. (int_of_num 0) <: (int_of_num n) <=> (~(n=0) )`,
1396   REWRITE_TAC[INT_OF_NUM_LT] THEN ARITH_TAC);;
1397
1398 let INT_POS_NEG =prove(`!n. ((int_of_num 0) <: (--:(int_of_num n))) <=> F`,
1399   REWRITE_TAC[INT_OF_NUM_LT] THEN ARITH_TAC);;
1400
1401 let RAT_LEMMA1_SUB = prove(`~(y1 = &0) /\ ~(y2 = &0) ==>
1402       ((x1 / y1) - (x2 / y2) = (x1 * y2 - x2 * y1) * inv(y1) * inv(y2))`,
1403   EVERY[REWRITE_TAC[real_div];
1404   REWRITE_TAC[real_sub;GSYM REAL_MUL_LNEG];
1405   REWRITE_TAC[GSYM real_div];
1406   SIMP_TAC[RAT_LEMMA1];
1407   DISCH_TAC;
1408   MESON_TAC[real_div]]);;
1409
1410 let INTERVAL_0 = prove(`! a f ex. (interval_eps a f ex <=> (&.0 <= (ex - (abs (a -. f)))))`,
1411    MESON_TAC[interval_eps;REAL_SUB_LE]);;
1412
1413 let ABS_NUM_P = prove_by_refinement(
1414   `!m n. (m <= n) ==> (abs (&. n -. (&. m)) = &.((m-|n) + (n-|m)))`,
1415 [
1416  REPEAT STRIP_TAC;
1417   ASM_SIMP_TAC[(ARITH_RULE `(m <= n ==> (m  - n = 0) )/\ ( 0 + x = x)`)];
1418   ASM_SIMP_TAC[GSYM REAL_OF_NUM_SUB;REAL_ABS_REFL;];
1419   FIRST_X_ASSUM MP_TAC;
1420   REWRITE_TAC[GSYM REAL_OF_NUM_LE];
1421   REAL_ARITH_TAC;
1422 ])
1423 ;;
1424
1425 let ABS_NUM = prove_by_refinement (
1426 `!m n. abs (&. n -. (&. m)) = &.((m-|n) + (n-|m))`,
1427 [
1428   REPEAT GEN_TAC ;
1429   DISJ_CASES_TAC (SPECL [`m:num`;`n:num`] LTE_CASES) ; (* THENL[ *)
1430   (* first case *)
1431   MATCH_MP_TAC ABS_NUM_P;
1432   FIRST_X_ASSUM MP_TAC;
1433   ARITH_TAC;
1434   (* second case *)
1435   ONCE_REWRITE_TAC[REAL_ARITH `abs(x - y) = abs(y - x)`];
1436   ONCE_REWRITE_TAC[ARITH_RULE `(a:num) + b = b + a`];
1437   MATCH_MP_TAC ABS_NUM_P;
1438   ASM_REWRITE_TAC[];
1439 ]);;
1440
1441
1442 let INTERVAL_TO_LESS = prove(
1443   `!a f ex b g ey. ((interval_eps a f ex) /\ (interval_eps b g ey) /\
1444       (&.0 <. (g -. (ey +. ex +. f)))) ==> (a <. b)`,
1445    let lemma1 = REAL_ARITH `!ex ey f g. (&.0 <.
1446      (g -. (ey +. ex +. f))) ==> ((f +. ex)<. (g -. ey)) ` in
1447    EVERY[
1448    REPEAT GEN_TAC;
1449    DISCH_ALL_TAC;
1450    H_MATCH_MP (THM lemma1) (HYP "2");
1451    H_MATCH_MP (THM INTERVAL_MAX) (HYP "0");
1452    H_MATCH_MP (THM INTERVAL_MIN) (HYP "1");
1453    LABEL_ALL_TAC;
1454    H_MATCH_MP (THM REAL_LET_TRANS) (H_RULE2 CONJ (HYP "4") (HYP "5"));
1455    LABEL_ALL_TAC;
1456    H_MATCH_MP (THM REAL_LTE_TRANS) (H_RULE2 CONJ (HYP "6") (HYP "3"));
1457    ASM_REWRITE_TAC[]
1458    ]);;
1459
1460 let ABS_TO_INTERVAL = prove(
1461   `!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))))`,
1462    EVERY[
1463    REWRITE_TAC[interval_eps];
1464    DISCH_ALL_TAC;
1465    REPEAT GEN_TAC;
1466    DISCH_ALL_TAC;
1467    ONCE_REWRITE_TAC [REAL_ARITH `c -. f = (c-. u) + (u-. f)`];
1468    ONCE_REWRITE_TAC [REAL_ADD_ASSOC];
1469    ASSUME_TAC (SPECL [`c-.u`;`u-.f`] ABS_TRIANGLE);
1470    IMP_RES_THEN ASSUME_TAC (REAL_ARITH `real_abs (k-.g) <=. ey ==> (k <=. (g +. ey))`);
1471    MATCH_MP_TAC (REAL_ARITH `(?a b.((x <=. (a+.b)) /\ (a <=. u) /\ (b <=. v)))  ==> (x <=. (u +. v))`);
1472    EXISTS_TAC `abs (c-.u)`;
1473    EXISTS_TAC `abs(u-.f)`;
1474    ASM_REWRITE_TAC[];
1475    ASM_MESON_TAC[REAL_LE_TRANS];
1476    ]);;
1477
1478
1479 (* end of general lemmas *)
1480 (* ------------------------------------------------------------------ *)
1481
1482
1483 (* ------------------------------------------------------------------ *)
1484 (* Cache of computed constants (abs (r - u) <= k)
1485     or interval r u eps *)
1486 (* ------------------------------------------------------------------ *)
1487
1488 let cache_entry_of_ineq bprec ineq =
1489   try(
1490   let (abst,_) = dest_binop `real_le` (concl ineq) in
1491   let (_,cmu) = dest_comb abst in
1492   let (r,_) = dest_binop `real_sub` cmu in
1493     (r,bprec,ineq))
1494   with Failure _ ->
1495   (try(
1496   let (r,_,_) = dest_interval (concl ineq) in (r,bprec,ineq))
1497   with Failure _ -> (print_term (concl ineq); failwith "calculated_constants format : abs(r - u) <= k"));;
1498
1499 let get_real_constant cache bprec tm =
1500  (fun (_,_,r) -> r) (find (fun (t,p,th) -> (t = tm) & (p>= bprec)) cache);;
1501
1502 let remove_entry t =   filter (fun x -> not(x = t));;
1503
1504 (* ------------------------------------------------------------------ *)
1505
1506 (* take |- interval_eps r u e to interval_eps r' u e, where th_r gives r <=> r' *)
1507
1508 let interval_rule th_r th_int =
1509   let (iru,e) = dest_comb (concl th_int) in
1510   let (ir,u) = dest_comb iru in
1511   let (i,r) = dest_comb ir in
1512   let ir_t = AP_TERM i th_r in
1513   let iru_t = AP_THM ir_t u in
1514   let irue_t = AP_THM iru_t e in
1515   let a = try EQ_MP irue_t th_int
1516      with Failure m -> (print_thm irue_t; print_thm th_int; failwith m) in
1517     a;;
1518
1519 (* term of the form '&.n'. Assume error checking done already. *)
1520
1521 let INTERVAL_OF_NUM:conv =
1522   fun tm ->
1523     let n = snd (dest_comb tm) in
1524     let th_r = (ARITH_REWRITE_CONV[] n) in
1525     let th_int = SPEC (rhs (concl th_r)) INTERVAL_NUM in
1526     interval_rule (SYM (AP_TERM `real_of_num` th_r)) th_int;;
1527
1528
1529 add_test("INTERVAL_OF_NUM",
1530    dest_thm (INTERVAL_OF_NUM `&.3`) = ([],
1531    `interval_eps (&3) (float (int_of_num 3) (int_of_num 0)) (float (int_of_num 0) (int_of_num 0))`));;
1532
1533 (* ------------------------------------------------------------------ *)
1534 (* conversion for float inequalities *)
1535 (* ------------------------------------------------------------------ *)
1536
1537
1538 let FLOAT_LESS_CONV = FLOAT_CONV THENC ARITH_SIMP_CONV[FLOAT_POS;
1539    INT_POS';INT_POS_NEG;FLOAT_NN;INT_NN;INT_NN_NEG;];;
1540
1541 let FLOAT_INEQ_CONV t =
1542   let th =FLOAT_LESS_CONV t in
1543   let (_,rhs) = dest_eq (concl th) in
1544   let _ =  (rhs = `T`) or raise Insufficient_precision in
1545      EQ_MP (SYM th) TRUTH;;
1546
1547 add_test("FLOAT_INEQ_CONV",
1548   let f c x =
1549     dest_thm (c x) = ([],x) in
1550   let t1 =
1551    `(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
1552     ((f FLOAT_INEQ_CONV t1)));;
1553
1554
1555 let INTERVAL_TO_LESS_CONV =
1556   let rthm = ASSUME `!f g ex ey. (&.0 <. (g -. (ey +. ex +. f)))` in
1557    fun thm1 thm2 ->
1558    let (a,f,ex) = dest_interval (concl thm1) in
1559    let (b,g,ey) = dest_interval (concl thm2) in
1560    let rspec = concl (SPECL [f;g;ex;ey] rthm) in
1561    let rspec_th = FLOAT_INEQ_CONV rspec in
1562    let fthm = CONJ thm1 (CONJ thm2 rspec_th) in
1563    MATCH_MP INTERVAL_TO_LESS fthm;;
1564
1565
1566
1567 add_test("INTERVAL_TO_LESS_CONV",
1568   let thm1 = ASSUME
1569    `interval_eps (#0.1) (float (int_of_num 1) (--: (int_of_num 1))) (float (int_of_num 1) (--: (int_of_num 2)))` in
1570   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
1571   let thm3 = INTERVAL_TO_LESS_CONV thm1 thm2 in
1572     concl thm3 = `#0.1 <. (#7)`);;
1573
1574 add_test("INTERVAL_TO_LESS_CONV2",
1575    let (h,c) = dest_thm (INTERVAL_TO_LESS_CONV
1576      (INTERVAL_OF_NUM `&.3`) (INTERVAL_OF_NUM `&.8`)) in
1577      (h=[]) && (c = `&.3 <. (&.8)`));;
1578
1579 (* ------------------------------------------------------------------ *)
1580
1581 (* conversion for DEC <= posfloat and posfloat <= DEC *)
1582
1583 let lemma1 = prove(
1584   `!n m p. ((&.p/(&.m)) <= (&.n)) <=> ((&.p/(&.m)) <= (&.n)/(&.1))`,
1585   MESON_TAC[REAL_DIV_1]);;
1586
1587 let lemma2 = prove(
1588   `!n m p. ((&.p) <= ((&.n)/(&.m))) <=> ((&.p/(&.1)) <= (&.n)/(&.m))`,
1589   MESON_TAC[REAL_DIV_1]);;
1590
1591 let REAL_LT = REAL_OF_NUM_LT;;
1592
1593 let REAL_LE = REAL_OF_NUM_LE;;
1594
1595 let lemma3 = prove(`!a b c d. (
1596    ((0<b) /\ (0<d) /\ (a*d <=| c*b))
1597     ==> (&.a/(&.b) <=. ((&.c)/(&.d))))`,
1598    EVERY[REPEAT GEN_TAC;
1599    DISCH_ALL_TAC;
1600    ASM_SIMP_TAC[RAT_LEMMA4;REAL_LT;REAL_OF_NUM_MUL;REAL_LE]]);;
1601
1602 let DEC_FLOAT = EQT_ELIM o
1603    ARITH_SIMP_CONV[DECIMAL;float;TWOPOW_POS;TWOPOW_NEG;GSYM real_div;
1604        REAL_OF_NUM_POW;INT_NUM_REAL;REAL_OF_NUM_MUL;
1605        lemma1;lemma2;lemma3];;
1606
1607 add_test("DEC_FLOAT",
1608    let f c x =
1609       dest_thm (c x) = ([],x) in
1610    ((f DEC_FLOAT `#10.0 <= (float (int_of_num 3) (int_of_num 2))`) &&
1611     (f DEC_FLOAT `#10 <= (float (int_of_num 3) (int_of_num 2))`) &&
1612     (f DEC_FLOAT `#0.1 <= (float (int_of_num 1) (--: (int_of_num 2)))`) &&
1613     (f DEC_FLOAT `float (int_of_num 3) (int_of_num 2) <= (#13.0)`) &&
1614     (f DEC_FLOAT `float (int_of_num 3) (int_of_num 2) <= (#13)`) &&
1615     (f DEC_FLOAT `float (int_of_num 1) (--: (int_of_num 2)) <= (#0.3)`)));;
1616
1617 (* Conversion for powers *)
1618
1619 let pow_parity = prove(
1620   `!x u. x pow (NUMERAL(BIT0 u)) = (x pow (NUMERAL u) * x pow (NUMERAL u)) /\
1621     x pow (NUMERAL(BIT1 u)) = x* x pow (NUMERAL u) * x pow (NUMERAL u)`,
1622   REWRITE_TAC[GSYM REAL_POW_ADD;NUMERAL;BIT0;BIT1;real_pow];
1623   );;
1624
1625 let pow_conv = PURE_ONCE_REWRITE_CONV[CONJUNCT1 real_pow;REAL_POW_1;pow_parity];;
1626
1627 (* conversion for negation *)
1628
1629 let neg_conv =
1630    let th = REAL_ARITH `--x = &0 - x` in
1631    PURE_ONCE_REWRITE_CONV[th];;
1632
1633 (* conversion for sums *)
1634
1635 let dest_sum  =
1636   let sum = `sum:(num->bool)->(num->real)->real` in
1637   let rg = `(..)` in
1638   fun tm ->
1639   let (x,f) = dest_binop sum tm in
1640   let (i,j) = dest_binop rg x in
1641     (i,j,f);;
1642
1643 let mk_sum =
1644   let sum = `sum:(num->bool)->(num->real)->real` in
1645   let rg = `(..)` in
1646   fun  (i,j,f)  ->
1647     let x = mk_binop rg i j in
1648     mk_binop sum x f;;
1649
1650 (*
1651 mk_sum (dest_sum `sum (3..4) f`);;
1652 *)
1653
1654 let sum_conv  =
1655   fun tm ->
1656   let (i,j,f) = dest_sum tm in
1657   let i',j' = dest_numeral i,dest_numeral j in
1658     if (i'=j') then  (ONCE_REWRITE_CONV[SUM_SING_NUMSEG])  tm
1659     else if (i' >/ j')  then
1660       let th = ARITH_RULE(mk_binop `(<):num->num->bool` j i) in
1661       let th2 = MATCH_MP SUM_TRIV_NUMSEG th in
1662         (PURE_ONCE_REWRITE_CONV[th2]) tm
1663     else
1664       let j_pred = mk_numeral (j'-/ Int 1) in
1665       let j_new = mk_comb (`SUC`, j_pred) in
1666       let th = (ISPECL [i;j_pred] (CONJUNCT2 SUM_CLAUSES_NUMSEG)) in
1667       let rr = ARITH_RULE (mk_eq (j_new,j)) in
1668       let th' = REWRITE_RULE[rr;ARITH_RULE (mk_binop `(<=):num->num->bool` i j)] th in
1669         PURE_ONCE_REWRITE_CONV[th'] tm;;
1670
1671 (* ------------------------------------------------------------------ *)
1672
1673 (* converts a DECIMAL TO A THEOREM *)
1674
1675 let ABS_BOUNDS = REAL_ABS_BOUNDS;; (* new *)
1676
1677 let INTERVAL_MINMAX = prove(`!x f ex.
1678    ((f -. ex) <= x) /\ (x <=. (f +. ex)) ==> (interval_eps x f ex)`,
1679    EVERY[REPEAT GEN_TAC;
1680    REWRITE_TAC[interval_eps;ABS_BOUNDS];
1681    REAL_ARITH_TAC]);;
1682
1683 let INTERVAL_OF_DECIMAL bprec dec =
1684   let a_num = dest_decimal dec in
1685   let f_num = round_rat bprec a_num in
1686   let ex_num = round_rat bprec (Num.abs_num (f_num -/ a_num)) in
1687   let _ = assert (ex_num <=/ f_num) in
1688   let f = mk_float f_num in
1689   let ex= mk_float ex_num in
1690   let fplus_ex = FLOAT_CONV (mk_binop `real_add` f ex) in
1691   let fminus_ex= FLOAT_CONV (mk_binop `real_sub` f ex) in
1692   let fplus_term = rhs (concl fplus_ex) in
1693   let fminus_term = rhs (concl fminus_ex) in
1694   let th1 = DEC_FLOAT (mk_binop `real_le` fminus_term dec) in
1695   let th2 = DEC_FLOAT (mk_binop `real_le` dec fplus_term) in
1696   let intv = mk_interval dec f ex in
1697   EQT_ELIM (SIMP_CONV[INTERVAL_MINMAX;fplus_ex;fminus_ex;th1;th2] intv);;
1698
1699 add_test("INTERVAL_OF_DECIMAL",
1700   let (h,c) = dest_thm (INTERVAL_OF_DECIMAL 4 `#36.1`) in
1701   let (x,f,ex) = dest_interval c in
1702    (h=[]) && (x = `#36.1`));;
1703
1704 add_test("INTERVAL_OF_DECIMAL2",
1705  can (fun() -> INTERVAL_TO_LESS_CONV (INTERVAL_OF_DECIMAL 4 `#33.33`)
1706   (INTERVAL_OF_DECIMAL 4 `#36.1`)) ());;
1707
1708 let INTERVAL_FLOAT = prove_by_refinement(
1709   `!a b. interval_eps (float a b) (float a b) (float (int_of_num 0) (int_of_num 0))`,
1710   (* {{{ proof *)
1711   [
1712   REWRITE_TAC[interval_eps;REAL_ARITH `abs(x - x) = &0`];
1713   REWRITE_TAC[float;int_of_num_th];
1714   REAL_ARITH_TAC;
1715   ]);;
1716   (* }}} *)
1717
1718 let INTERVAL_FLOAT_CONV flt =
1719   let (a,b) = dest_binop `float` flt in
1720     (ISPECL [a;b] INTERVAL_FLOAT);;
1721
1722 (*--------------------------------------------------------------------*)
1723 (*   functions to check format.                                       *)
1724 (*   There are various implicit rules:                                *)
1725 (*   NUMERAL is followed by bits and no other kind of num, etc.       *)
1726 (*   FLOAT a b, both a and b are int_of_num NUMERAL or --:int_of_num NUMERAL, etc.      *)
1727 (*--------------------------------------------------------------------*)
1728
1729
1730 (* converts exceptions to false *)
1731 let falsify_ex f x = try (f x) with Failure _ -> false;;
1732
1733 let is_bits_format  =
1734   let zero = `_0` in
1735     let rec format x =
1736     if (x = zero) then true
1737     else let (h,t) = dest_comb x  in
1738       (((h = `BIT1`) or (h = `BIT0`)) && (format t))
1739     in falsify_ex format;;
1740
1741 let is_numeral_format =
1742     let fn x =
1743     let (h,t) = dest_comb x in
1744       ((h = `NUMERAL`) && (is_bits_format t)) in
1745     falsify_ex fn;;
1746
1747 let is_decimal_format  =
1748     let fn x =
1749       let (t1,t2) = dest_binop `DECIMAL` x in
1750         ((is_numeral_format t1) && (is_numeral_format t2)) in
1751     falsify_ex fn;;
1752
1753 let is_pos_int_format =
1754     let fn x =
1755       let (h,t) = dest_comb x in
1756        (h = `int_of_num`) && (is_numeral_format t) in
1757     falsify_ex fn;;
1758
1759 let is_neg_int_format =
1760     let fn x =
1761       let (h,t) = dest_comb x in
1762         (h = `real_neg`) && (is_pos_int_format t) in
1763       falsify_ex fn;;
1764
1765 let is_int_format x =
1766   (is_neg_int_format x) or (is_pos_int_format x);;
1767
1768 let is_float_format =
1769     let fn x =
1770       let (t1,t2) = dest_binop `float` x in
1771       (is_int_format t1) && (is_int_format t2) in
1772     falsify_ex fn;;
1773
1774 let is_interval_format =
1775   let fn x =
1776     let (a,b,c) = dest_interval x in
1777       (is_float_format b) && (is_float_format c) in
1778     falsify_ex fn;;
1779
1780 let is_neg_real =
1781   let fn x =
1782      let (h,t) = dest_comb x in
1783       (h= `real_neg `) in
1784     falsify_ex fn;;
1785
1786 let is_real_num_format =
1787   let fn x =
1788     let (h,t) = dest_comb x in
1789       (h=`real_of_num`) && (is_numeral_format t) in
1790   falsify_ex fn;;
1791
1792 let is_comb_of t u =
1793   let fn t u =
1794     t = (fst (dest_comb u)) in
1795   try (fn t u) with Failure _ -> false;;
1796
1797 (* ------------------------------------------------------------------ *)
1798 (* Heron's formula for the square root of A
1799    Return a value x that is always at most the actual square root
1800    and such that abs (x  - A/x ) < epsilon *)
1801
1802 let rec heron_sqrt depth A x eps =
1803     let half = (Int 1)//(Int 2) in
1804     if (depth <= 0)
1805     then raise (Failure "heron_sqrt:internal error:sqrt recursion depth exceeded")
1806     else
1807     if (Num.abs_num (x -/ (A//x) ) </ eps) & (x*/ x >=/ A)  then (A//x) else
1808     let x' = half */ (x +/ (A//x)) in
1809     heron_sqrt (depth -1) A x' eps;;
1810
1811 let INTERVAL_OF_TWOPOW = prove(
1812    `!n. interval_eps (twopow n) (float (int_of_num 1) n) (float (int_of_num 0) (int_of_num 0))`,
1813    REWRITE_TAC[interval_eps;float;int_of_num_th] THEN
1814    REAL_ARITH_TAC
1815    );;
1816
1817 (* compute ez parameter for division *)
1818
1819 let  division_ez bprec (f,ex,g,ey) =
1820     let f_num = dest_float f in
1821     let ex_num = dest_float ex in
1822     let g_num = dest_float g in
1823     let ey_num = dest_float ey in
1824     let _ = (ey_num </ abs_num g_num) or raise Insufficient_precision in
1825     let h_num = round_rat bprec (f_num//g_num) in
1826     let h = mk_float h_num in
1827     let ez_rat = (ex_num +/ abs_num (f_num -/ (h_num*/ g_num))
1828         +/ (abs_num h_num */ ey_num))//((abs_num g_num) -/ (ey_num)) in
1829     let ez_num = round_rat bprec (ez_rat) in
1830     let _ = assert((ez_num >=/ (Int 0))) in
1831     let ez = mk_float ez_num
1832     in (h,ez);;
1833
1834 let pureconv (s,t) = (s,PURE_ONCE_REWRITE_CONV[t]);;
1835
1836 (* ------------------------------------------------------------------ *)
1837
1838 (*
1839 bprec gives the bits of precision.
1840 *)
1841
1842
1843 let rec (INTERVAL_OF_TERM:(string*conv)list -> (term *int*thm)list ->int -> term -> (term*int*thm) list * thm) =
1844 fun convl cache  bprec tm ->
1845   let mk c th = (cache_entry_of_ineq bprec th :: c, th) in
1846   (* treat user conversions *)
1847   let interval_of_conv cache bprec conv tm =
1848     let th_r =  conv tm in
1849     let (lhs,rhs) = dest_eq (concl th_r) in
1850     let _ = (not (aconv lhs rhs)) or (print_term lhs; failwith "interval_of_conv: no change"  )in
1851     let (cache,th_rhs) = INTERVAL_OF_TERM convl cache bprec rhs in
1852     let th_lhs = interval_rule (SYM th_r) th_rhs in
1853       mk cache th_lhs in
1854   let interval_of_binop op opname match_thm =
1855     begin
1856     let (a,b) = dest_binop op tm in
1857     let (cache,a_int) = INTERVAL_OF_TERM convl cache  bprec a in
1858     let (cache,b_int) = INTERVAL_OF_TERM convl cache  bprec b in
1859       try(
1860     let c_int = MATCH_MP match_thm (CONJ a_int b_int) in
1861     let th = EQ_MP (INTERVAL_CLEAN_CONV (concl c_int)) c_int in
1862       mk cache th)
1863     with Failure _ -> failwith ("INTERVAL_OF_TERM: "^opname)
1864     end in
1865   (* treat cached values first *)
1866   if (can (get_real_constant cache bprec) tm) then
1867     begin
1868     try(
1869     let int_thm = get_real_constant cache bprec tm in
1870     if (can dest_interval (concl int_thm)) then (cache,int_thm)
1871     else (
1872     let absthm = int_thm in
1873     let (abst,k) = dest_binop `real_le` (concl absthm) in
1874     let (absh,cmu) = dest_comb abst in
1875     let (c,u) = dest_binop `real_sub` cmu in
1876     let (cache,intk) = INTERVAL_OF_TERM convl cache  bprec k in
1877     let (cache,intu) = INTERVAL_OF_TERM convl cache  bprec u in
1878     let thm1 = MATCH_MP ABS_TO_INTERVAL absthm in
1879     let thm2 = MATCH_MP thm1 (CONJ intu intk) in
1880     let (_,f,ex)= dest_interval (concl thm2) in
1881     let fthm = FLOAT_CONV f in
1882     let exthm = FLOAT_CONV ex in
1883     let thm3 = REWRITE_RULE[fthm;exthm] thm2 in
1884      (cache_entry_of_ineq bprec thm3 :: (remove_entry (tm,bprec,int_thm) cache), thm3)
1885     ))
1886     with Failure _ -> failwith "INTERVAL_OF_TERM  : CONSTANT"
1887     end
1888   else if (is_real_num_format tm) then
1889     let th = INTERVAL_OF_NUM tm in (cache,th)
1890   else if (is_decimal_format tm) then
1891     let th =  (INTERVAL_OF_DECIMAL bprec tm)
1892     in mk cache th
1893   (* treat abs value *)
1894   else if (is_comb_of `real_abs` tm) then
1895     begin
1896       let (_,b) = dest_comb tm in
1897       let (cache,i1) = (INTERVAL_OF_TERM convl cache  bprec b) in
1898       try(
1899       let b_int = MATCH_MP INTERVAL_ABSV i1 in
1900       let th = EQ_MP (INTERVAL_CLEAN_CONV (concl b_int)) b_int in
1901         mk cache th)
1902       with Failure _ -> failwith "INTERVAL_OF_TERM  : ABS"
1903     end
1904   (* treat twopow *)
1905   else if (is_comb_of `twopow` tm) then
1906     begin
1907       try(
1908       let (_,b) = dest_comb tm in
1909       let th = SPEC b INTERVAL_OF_TWOPOW
1910       in mk cache th
1911       )
1912       with Failure _ -> failwith "INTERVAL_OF_TERM  : TWOPOW"
1913     end
1914   (* treat addition *)
1915   else if (can (dest_binop `real_add`) tm) then
1916     interval_of_binop `real_add` "ADD" INTERVAL_ADD
1917   (* treat subtraction *)
1918   else if (can (dest_binop `real_sub`) tm) then
1919     interval_of_binop `real_sub` "SUB" INTERVAL_SUB
1920   (* treat multiplication *)
1921   else if (can (dest_binop `real_mul`) tm) then
1922      interval_of_binop `real_mul` "MUL" INTERVAL_MUL
1923   (* treat division : instantiate INTERVAL_DIV *)
1924   else if (can (dest_binop `( / )`) tm) then
1925     begin
1926     let (a,b) = dest_binop `( / )` tm in
1927     let (cache,a_int) = INTERVAL_OF_TERM convl cache  bprec a in
1928     let (cache,b_int) = INTERVAL_OF_TERM convl cache  bprec b in
1929       try(
1930     let (_,f,ex) = dest_interval (concl a_int) in
1931     let (_,g,ey) = dest_interval (concl b_int) in
1932     let (h,ez) = division_ez bprec (f,ex,g,ey) in
1933     let hyp1 = a_int in
1934     let hyp2 = b_int in
1935     let hyp3 = try (FLOAT_INEQ_CONV (mk_binop `real_lt` ey (mk_comb (`real_abs`,g))))
1936       with Failure _ -> raise Insufficient_precision in
1937     let thm = SPECL [a;f;ex;b;g;ey;h;ez] INTERVAL_DIV in
1938     let conj2 x = snd (dest_conj x) in
1939     let hyp4t = (conj2 (conj2 (conj2 (fst(dest_imp (concl thm)))))) in
1940     let hyp4 = FLOAT_INEQ_CONV hyp4t in
1941     let hyp_all = end_itlist CONJ [hyp1;hyp2;hyp3;hyp4] in
1942     let th = MATCH_MP thm hyp_all in
1943       mk cache th)
1944     with Failure _ -> failwith "INTERVAL_OF_TERM  :DIV"
1945     end
1946   (* treat pow *)
1947   else if (can (dest_binary "real_pow") tm) then
1948      interval_of_conv cache bprec pow_conv tm
1949   (* treat negative terms *)
1950   else if (is_neg_real tm) then
1951     interval_of_conv cache bprec neg_conv tm
1952   (* treat beta *)
1953   else if  can (fun t-> let (f,x) = dest_comb t in dest_abs f )  tm then
1954      interval_of_conv cache bprec BETA_CONV tm
1955   (* treat sum *)
1956   else if can dest_sum tm then
1957      interval_of_conv cache bprec sum_conv tm
1958   (* treat float *)
1959   else if can (dest_binop `float`) tm then
1960     (cache,INTERVAL_FLOAT_CONV tm)
1961   (* treat ssqrt : instantiate INTERVAL_SSQRT *)
1962   else if (can (dest_unary `ssqrt`) tm) then
1963     begin
1964     try(
1965     let x = dest_unary `ssqrt` tm in
1966     let (cache,x_int) = INTERVAL_OF_TERM convl cache   bprec x in
1967     let (_,f,ex)  = dest_interval (concl x_int) in
1968     let f_num = dest_float f in
1969     let ex_num = dest_float ex in
1970     let fd_num = f_num -/ ex_num in
1971     let fe_f = Num.float_of_num fd_num in
1972     let apprx_sqrt = Pervasives.sqrt fe_f in
1973     (* put in heron's formula *)
1974     let v_num1 = round_inward_float 25 (apprx_sqrt) in
1975     let v_num = round_inward_rat bprec
1976          (heron_sqrt 10 fd_num v_num1 ((Int 2) **/ (Int (-bprec-4)))) in
1977     let u_num1 = round_inward_float 25
1978         (Pervasives.sqrt (float_of_num f_num)) in
1979     let u_num = round_inward_rat bprec
1980         (heron_sqrt 10 f_num u_num1 ((Int 2) **/ (Int (-bprec-4)))) in
1981     let ey_num = round_rat bprec (abs_num (f_num -/ (u_num */ u_num))) in
1982     let ez_num = round_rat bprec ((ex_num +/ ey_num)//(u_num +/ v_num)) in
1983     let (v,u) = (mk_float v_num,mk_float u_num) in
1984     let (ey,ez) = (mk_float ey_num,mk_float ez_num) in
1985     let thm = SPECL [x;f;ex;u;ey;ez;v] INTERVAL_SSQRT in
1986     let conjhyp = fst (dest_imp (concl thm)) in
1987     let [hyp6;hyp5;hyp4;hyp3;hyp2;hyp1] =
1988       let rec break_conj c acc =
1989         if (not(is_conj c)) then (c::acc) else
1990         let (u,v) = dest_conj c in break_conj v (u::acc) in
1991        (break_conj conjhyp []) in
1992     let thm2 = prove(hyp2,REWRITE_TAC[interval_eps] THEN
1993                        (CONV_TAC FLOAT_INEQ_CONV)) in
1994     let thm3 = FLOAT_INEQ_CONV hyp3 in
1995     let thm4 = FLOAT_INEQ_CONV hyp4 in
1996     let float_tac = REWRITE_TAC[FLOAT_NN;FLOAT_POS;INT_NN;INT_NN_NEG;
1997                        INT_POS';INT_POS_NEG] THEN ARITH_TAC in
1998     let thm5 = prove( hyp5,float_tac) in
1999     let thm6 = prove( hyp6,float_tac) in
2000     let ant  = end_itlist CONJ[x_int;thm2;thm3;thm4;thm5;thm6] in
2001     let th =    MATCH_MP thm ant
2002     in mk cache th
2003     )
2004     with Failure _ -> failwith "INTERVAL_OF_TERM  : SSQRT"
2005     end
2006  else if (can (dest_unary `sqrt`) tm) then
2007     begin
2008     try(
2009     let x = dest_unary `sqrt` tm in
2010     let stm = mk_comb (`ssqrt`, x) in
2011     let (cache,ssqrt_thm) = INTERVAL_OF_TERM convl cache  bprec stm in
2012     let (cache,int1) = INTERVAL_OF_TERM convl cache  bprec`real_of_num 0` in
2013     let (cache,int2)  =INTERVAL_OF_TERM convl cache  bprec x in
2014     let ant = CONJ (INTERVAL_TO_LESS_CONV int1 int2) ssqrt_thm in
2015     let th = MATCH_MP ssqrt_sqrt ant
2016     in mk cache th)
2017     with  Failure _ -> failwith "INTERVAL_OF_TERM : SQRT"
2018     end
2019 (* treat atn *)
2020  else if can (dest_unary `atn`) tm then
2021    begin
2022      let half_co = ("halfatn4_co",ARITH_SIMP_CONV[Taylor_atn.halfatn4_co]) in
2023      let half4 = ("halfatn4",PURE_REWRITE_CONV[Taylor_atn.halfatn4;o_THM]) in
2024      let halfatnconv = pureconv("halfatn",Taylor_atn.halfatn) in
2025      let aconvl =  (halfatnconv::half4::half_co::convl) in
2026      let n = mk_numeral (Int (1 + (bprec - 5) / 6)) in
2027      let x = dest_unary `atn` tm in
2028      let sum_tm = instantiate ([],[(n,`n:num`);(x,`u:real`)],[]) `&16 * sum (0..n) (halfatn4_co u)` in
2029      let (cache,hyp_1) = INTERVAL_OF_TERM (aconvl) cache bprec sum_tm in
2030      let (_,_,eps2) = dest_interval (concl hyp_1) in
2031      let eps_tm = instantiate([],[(n,`n:num`);(eps2,`eps2:real`)],[]) `float (int_of_num 1) ( -- int_of_num (6 * n + 5)) + eps2` in
2032      let eps_thm = FLOAT_CONV eps_tm in
2033      let eps = snd(dest_eq (concl eps_thm)) in
2034      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) )
2035    end
2036 (* treat acs *)
2037  else if can (dest_unary `acs`) tm then
2038    begin
2039      let x = dest_unary `acs` tm in
2040      let (_,int_m1) = INTERVAL_OF_TERM convl cache  bprec`-- real_of_num 1` in
2041      let (_,int_p1) = INTERVAL_OF_TERM convl cache  bprec`real_of_num 1` in
2042      let (cache,int_x)  =INTERVAL_OF_TERM convl cache  bprec x in
2043      let ant1 = INTERVAL_TO_LESS_CONV int_m1 int_x in
2044      let ant2 = INTERVAL_TO_LESS_CONV int_x int_p1 in
2045      let th =  MATCH_MP ACS_ATN (CONJ ant1 ant2) in
2046        interval_of_conv cache bprec (PURE_ONCE_REWRITE_CONV[th]) tm
2047    end
2048 (* treat pi *)
2049  else if (tm = `pi`) then
2050     begin
2051       interval_of_conv cache bprec (PURE_ONCE_REWRITE_CONV[pi_atn]) tm
2052     end
2053 (* user supplied conversions *)
2054   else
2055     begin
2056       let (c,b) = strip_comb tm in
2057       let _ = is_const c or failwith "INTERVAL_OF_TERM: case not installed" in
2058       let (s,_) = dest_const c in
2059       let conv = try (assoc s convl)
2060          with Failure _ -> failwith ("INTERVAL_OF_TERM: function not found "^s) in
2061          interval_of_conv cache bprec conv tm
2062     end;;
2063
2064 let real_ineq convl cache bprec tm =
2065   let (t1,t2) = try (dest_binop `real_lt` tm)
2066     with Failure _ -> failwith "real_ineq: inequality of the form `x < y` expected" in
2067   let (cache,int1) = INTERVAL_OF_TERM convl cache  bprec t1 in
2068   let (cache,int2) = INTERVAL_OF_TERM convl cache  bprec t2 in
2069   (cache,INTERVAL_TO_LESS_CONV int1 int2);;
2070
2071 let rec rec_real_ineq convl cache bprec bprec_max tm =
2072    try (bprec,real_ineq convl cache bprec tm) with
2073    Insufficient_precision ->
2074      let bprec = (16 * bprec) / 10 in
2075      let _ =  (bprec < bprec_max) or failwith ("insufficient precision "^ (string_of_int bprec_max)) in
2076      let _ = Format.print_string ("increasing precision to "^(string_of_int bprec)^" bits");
2077                  Format.print_newline() in
2078      rec_real_ineq convl cache bprec  bprec_max tm;;
2079
2080 (* Give it a conjunction of strict inequalities *)
2081
2082 let REAL_INEQUALITY_CALCULATOR_PREC (bprec_init,bprec_max) convl tm =
2083   let tml = striplist dest_conj tm in
2084   let (b,c,th) = List.fold_left (fun (b,c,thl) tm1 ->
2085                         let (bprec,(cache,th1)) = rec_real_ineq convl c b bprec_max tm1
2086        in (bprec,cache, th1 ::thl)) (bprec_init,[], []) (tml) in
2087     (b, EQ_MP (SYM(REWRITE_CONV th tm)) TRUTH);;
2088
2089 let REAL_INEQUALITY_CALCULATOR convl tm =
2090    REAL_INEQUALITY_CALCULATOR_PREC (10,70) convl tm;;
2091
2092 (*
2093 let REAL_INEQUALITY_CALCULATOR convl tm =
2094   let tml = striplist dest_conj tm in
2095   let bprec_max = 70 in
2096   let (b,c,th) = List.fold_left (fun (b,c,thl) tm1 ->
2097                         let (bprec,(cache,th1)) = rec_real_ineq convl c b bprec_max tm1
2098        in (bprec,cache, th1 ::thl)) (10,[], []) (tml) in
2099     (b, EQ_MP (SYM(REWRITE_CONV th tm)) TRUTH);;
2100 *)
2101
2102 (*
2103   let tmint = map (rec_real_ineq [] 10) tml in
2104    List.fold_left CONJ (hd tmint) (tl tmint);;
2105
2106 *)
2107
2108 end;;