Jordan/float.ml

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