1 (* ========================================================================== *)
2 (* FLYSPECK - BOOK FORMALIZATION *)
5 (* Copied from HOL Light jordan directory *)
6 (* Author: Thomas C. Hales *)
8 (* ========================================================================== *)
11 A handheld calculator, implemented in HOL light.
13 It implements + - / * sqrt, ssqrt, abs, atn, acs, pi,
14 The user can supply additional conversons to extend the calculator to other functions.
16 There is an example file float_example.hl.
19 (* based on float.ml in hol/jordan/float.ml *)
22 To do. Allow the user to supply theorems that carry domain constraints,
26 flyspeck_needs "jordan/taylor_atn.hl";;
37 Parse_ext_override_interface.unambiguous_interface();;
38 Parse_ext_override_interface.prioritize_real();;
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.
45 Internal failure: a bug in this code.
47 User input failure: sqrt of negative numbers, false inequalities.
51 exception Insufficient_precision;;
54 exception Mk_comb of (term*hol_type)*(term*hol_type);;
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)));;
59 let add_test,test = new_test_suite();;
63 `twopow x = if (?n. (x = (int_of_num n)))
64 then ((&2) pow (nabs x))
65 else inv((&2) pow (nabs x))`);;
69 `float x n = (real_of_int x)*(twopow n)`);;
73 `interval_eps x f eps = ((abs (x-f)) <= eps)`);;
75 (*--------------------------------------------------------------------*)
78 let mk_interval a b ex =
79 mk_comb(mk_comb (mk_comb (`interval_eps`,a),b),ex);;
81 add_test("mk_interval",
82 mk_interval `#3` `#4` `#1` = `interval_eps #3 #4 #1`);;
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
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));;
95 (*--------------------------------------------------------------------*)
97 let (dest_int:term-> Num.num) =
100 let (op,nat) = dest_comb a in
101 if (fst (dest_const op) = "int_of_num") then (dest_numeral nat)
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
107 Failure _ -> failwith "dest_int ";;
110 let (mk_int:Num.num -> term) =
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));;
119 (mk_int (Int (-1443)) = `int_neg (int_of_num 1443)`) &&
120 (mk_int (Int 37) = `(int_of_num 37)`));;
122 (* ------------------------------------------------------------------ *)
124 let (split_ratio:Num.num -> Num.num*Num.num) =
126 (Ratio r) -> (Big_int (Ratio.numerator_ratio r)),
127 (Big_int (Ratio.denominator_ratio r))|
130 add_test("split_ratio",
131 let (a,b) = split_ratio ((Int 4)//(Int 20)) in
132 (a =/ (Int 1)) && (b =/ (Int 5)));;
134 (* ------------------------------------------------------------------ *)
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) =
139 let intC = (Int c) in
141 if ((Int 0) =/ mod_num a intC) then (divC (a//intC,b+/(Int 1)))
144 if ((Num.is_integer_num r)&& not((Num.sign_num r) = 0)) then
145 divC (r,(Int 0)) else failwith "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))));;
152 (*--------------------------------------------------------------------*)
154 let (dest_float:term -> Num.num) =
156 let (a,b) = dest_binop `float` f in
157 (dest_int a)*/ ((Int 2) **/ (dest_int b));;
159 add_test("dest_float",
160 dest_float `float (int_of_num 3) (int_of_num 17)` = (Int 393216));;
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)`);;
165 (* ------------------------------------------------------------------ *)
166 (* creates float of the form `float a b` with a odd *)
167 let (mk_float:Num.num -> term) =
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
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";;
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)`));;
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))));;
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)`);;
189 (* ------------------------------------------------------------------ *)
190 (* creates decimal of the form `DECIMAL a b` with a prime to 10 *)
191 let (mk_pos_decimal:Num.num -> term) =
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');;
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`));;
216 let (mk_decimal:Num.num->term) =
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;;
222 add_test("mk_decimal",
223 (mk_decimal (Int 3) = `#3`) &&
224 (mk_decimal (Int (-3)) = `real_neg (#3)`));;
228 (*--------------------------------------------------------------------*)
230 let (dest_decimal:term -> Num.num) =
232 let (a,b) = dest_binop `DECIMAL` f in
233 let a1 = dest_numeral a in
234 let b1 = dest_numeral b in
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)`);;
246 (*--------------------------------------------------------------------*)
247 (* Properties of integer powers of 2. *)
248 (* ------------------------------------------------------------------ *)
251 let TWOPOW_POS = prove(`!n. (twopow (int_of_num n) = (&2) pow n)`,
252 (REWRITE_TAC[twopow])
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[]));;
261 let TWOPOW_NEG = prove(`!n. (twopow (--(int_of_num n)) = inv((&2) pow n))`,
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])
273 THEN (ASM_REWRITE_TAC[real_pow;REAL_INV_1]));;
276 let TWOPOW_INV = prove(`!a. (twopow (--: a) = (inv (twopow a)))`,
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])));;
283 let INT_REP3 = prove(`!a .(?n.( (a = int_of_num n) \/ (a = --: (int_of_num (n+1)))))`,
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[]))))
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]
294 THEN ((EXISTS_TAC `PRE n`))
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])
302 let REAL_EQ_INV = prove(`!x y. ((x:real = y) <=> (inv(x) = inv (y)))`,
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[]))]);;
311 prove(`!a. (twopow (a +: (int_of_num 1)) = twopow (a) *. (twopow (int_of_num 1)))`,
314 CHOOSE_TAC (SPEC `a:int` INT_REP3);
315 POP_ASSUM DISJ_CASES_TAC
317 ASM_REWRITE_TAC[TWOPOW_POS;INT_OF_NUM_ADD;REAL_POW_ADD];
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]
329 let TWOPOW_0 = prove(`twopow (int_of_num 0) = (&. 1)`,
330 (REWRITE_TAC[TWOPOW_POS;real_pow]));;
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)))`,
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])));;
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]));;
348 let TWOPOW_ADD_INT = prove(
349 `!a b. (twopow (a +: b) = twopow(a) *. (twopow(b)))`,
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'))`))
359 THENL[ ((REWRITE_TAC[GSYM INT_OF_NUM_ADD]))
360 THEN ((INT_ARITH_TAC));ALL_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]));;
370 let TWOPOW_ABS = prove(`!a. real_abs (twopow a) = (twopow a)`,
372 THEN ((CHOOSE_THEN DISJ_CASES_TAC (SPEC `a:int` INT_REP2)))
374 THEN ((ASM_REWRITE_TAC[TWOPOW_POS;TWOPOW_NEG;REAL_ABS_POW;REAL_ABS_NUM;REAL_ABS_INV])));;
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))
380 THENL[ ((REWRITE_TAC[real_pow;MULT_0]));
382 ((REWRITE_TAC[real_pow]))
383 THEN ((ASM_REWRITE_TAC[ADD1;LEFT_ADD_DISTRIB;REAL_POW_ADD;REAL_MUL_AC;MULT_CLAUSES]))]);;
385 let INT_POW_POW = INT_OF_REAL_THM REAL_POW_POW;;
387 let TWOPOW_POW = prove(
388 `!a n. (twopow a) pow n = twopow (a *: (int_of_num n))`,
390 THEN ((CHOOSE_THEN DISJ_CASES_TAC (SPEC `a:int` INT_REP2)))
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])));;
395 (* ------------------------------------------------------------------ *)
396 (* Arithmetic operations on float *)
397 (* ------------------------------------------------------------------ *)
400 let FLOAT_NEG = prove(`!a m. real_neg (float a m) = float (--: a) m`,
402 REWRITE_TAC[float;GSYM REAL_MUL_LNEG;int_neg_th]);;
406 let FLOAT_MUL = prove(`!a b m n. (float a m) *. (float b n) = (float (a *: b) (m +: n))`,
408 THEN ((REWRITE_TAC[float;GSYM REAL_MUL_ASSOC;TWOPOW_ADD_INT;int_mul_th]))
409 THEN ((MESON_TAC[REAL_MUL_AC])));;
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])));;
419 let FLOAT_ADD_EQ = prove(
420 `!a b m. (float a m) +. (float b m) =
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])));;
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]));;
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)))`,
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));
438 THEN ((ASM_REWRITE_TAC[FLOAT_ADD]))]);;
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))))`,
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]))
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)));
452 THEN ((ASM_REWRITE_TAC[FLOAT_ADD]))]);;
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))))`,
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]))
462 THEN ((REWRITE_TAC[REAL_ADD_AC])));;
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)))))`,
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)]))
473 THEN ((INT_ARITH_TAC));
476 THEN (ASM_REWRITE_TAC[GSYM FLOAT_ADD;REAL_ADD_AC])]);;
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)))))`,
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])))
487 THEN (((REWRITE_TAC[REAL_ADD_AC]))));;
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]);;
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]));;
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]));;
504 `!a b. (a -: b) = (a +: (--: b))`,
505 (REWRITE_TAC[GSYM INT_SUB_RNEG;INT_NEG_NEG]));;
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]));;
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]));;
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]));;
518 let INT_POW_MUL = INT_OF_REAL_THM REAL_POW_MUL;;
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]));;
529 GEN_TAC THEN REWRITE_TAC[num_CONV `1`] THEN
530 REWRITE_TAC[pow; REAL_MUL_RID]);;
533 `!x. x pow 2 = x * x`,
534 GEN_TAC THEN REWRITE_TAC[num_CONV `2`] THEN
535 REWRITE_TAC[pow; POW_1]);;
538 let INT_POW_EVEN_NEG1 = prove(
539 `!x n. (--: (int_of_num x)) **: (NUMERAL (BIT0 n)) = ((int_of_num x) **: (NUMERAL (BIT0 n)))`,
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])));;
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]);;
564 GEN_TAC THEN REWRITE_TAC[num_CONV `1`] THEN
565 REWRITE_TAC[pow; REAL_MUL_RID]);;
567 let INT_POW_ODD_NEG1 = prove(
568 `!x n. (--: (int_of_num x)) **: (NUMERAL (BIT1 n)) = --: ((int_of_num x) **: (NUMERAL (BIT1 n)))`,
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])));;
580 (* subtraction of integers *)
582 let INT_ADD_NEG = prove(
583 `!x y. (x < y ==> ((int_of_num x) +: (--: (int_of_num y)) = --: (int_of_num (y - x))))`,
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`))
589 THENL [(((ASM MESON_TAC)[LE_LT]));((SIMP_TAC[GSYM (INT_OF_REAL_THM REAL_OF_NUM_SUB)]))]
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)`))]);;
595 let INT_ADD_NEGv2 = prove(
596 `!x y. (y <= x ==> ((int_of_num x) +: (--: (int_of_num y)) = (int_of_num (x - y))))`,
599 THEN ((SUBGOAL_MP_TAC `int_of_num (x - y) = (int_of_num x) - (int_of_num y)`))
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]))
607 (* assumes a term of the form int_of_num a +: (--: (int_of_num b)) *)
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
618 (ARITH_SIMP_CONV[thm]) t;; (* (SIMP_CONV[thm;ARITH]) t;; *)
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]);;
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]);;
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]);;
632 (* ------------------------------------------------------------------ *)
633 (* Simplifies an arithmetic expression in floats *)
635 (* ------------------------------------------------------------------ *)
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;
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 *)
649 add_test("FLOAT_CONV1",
651 let (x,y) = dest_eq z in
652 let (u,v) = dest_thm (FLOAT_CONV x) 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))
661 f `real_neg (float (--: (int_of_num 3)) (int_of_num 0)) = float (int_of_num 3) (int_of_num 0)`
664 (* ------------------------------------------------------------------ *)
665 (* Operations on interval_eps. Preliminary stuff to deal with *)
666 (* chains of inequalities. *)
667 (* ------------------------------------------------------------------ *)
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)))`,
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]))));;
679 let REAL_ADD_LE_SUBST_LHS = prove(
680 `!a b c P. (((P(a) <=. b /\ (!x. (P x) = x + (P (&. 0))) /\ (c <=. a)))
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])));;
690 (a::b) -> fun thm ->(SPECL b (SPEC a thm));;
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.
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);;
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.
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);;
748 (* ------------------------------------------------------------------ *)
749 (* INTERVAL ADD, NEG, SUBTRACT *)
750 (* ------------------------------------------------------------------ *)
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)`,
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));
766 ABBREV_TAC `a':real = abs a`;
767 ABBREV_TAC `b':real = abs b`;
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[]]);;
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)`]));;
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))`]));;
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]))
786 THEN (((H_RULER (ONCE_REWRITE_RULE))[THM(INTERVAL_NEG)] (HYP_INT 1)))
787 THEN (((ASM MESON_TAC)[INTERVAL_ADD])));;
789 (* ------------------------------------------------------------------ *)
790 (* INTERVAL MULTIPLICATION *)
791 (* ------------------------------------------------------------------ *)
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]));;
799 let REAL_LE_ABS_MUL = prove(
800 `!x y z w. (( x <=. y) /\ (abs z <=. w)) ==> (x*.w <=. y*.w) `,
802 THEN ((ASSUME_TAC (REAL_ARITH `abs z <=. w ==> (&.0) <=. w`)))
803 THEN (((ASM MESON_TAC)[REAL_LE_RMUL_IMP])));;
805 let ABS_MUL = REAL_ABS_MUL;;
807 let INTERVAL_ABS = REAL_ARITH
808 `!(x:real) z d. (x - d) <= z /\ z <= (x + d) <=> abs(z - x) <= d`;;
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))`,
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))`]))
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)));;
833 (* ------------------------------------------------------------------ *)
834 (* INTERVAL BASIC OPERATIONS *)
835 (* ------------------------------------------------------------------ *)
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]));;
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]))
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[])));;
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]))
854 THEN ((H_VAL2 IWRITE_REAL_LE_RHS (HYP_INT 1) (HYP_INT 0)))
855 THEN ((ASM_REWRITE_TAC[])));;
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);;
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);;
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)
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)
875 let REAL_RINV_2 = prove(
876 `&.2 *. (inv (&.2 )) = &. 1`,
878 MATCH_MP_TAC REAL_MUL_RINV;
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)) ==>
886 ((xmax-.xmin)*.half)`,
888 REWRITE_TAC[LET_DEF;LET_END_DEF];
890 REWRITE_TAC[interval_eps;float;TWOPOW_NEG;INT_NUM_REAL;REAL_POW_1;REAL_MUL_LID];
891 REWRITE_TAC[GSYM INTERVAL_ABS];
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]
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]
907 let INTERVAL_EPS_POS = prove(`!x f ex.
908 (interval_eps x f ex) ==> (&.0 <=. ex)`,
910 REWRITE_TAC[interval_eps];
912 DISCH_THEN(fun x -> (MP_TAC (CONJ (SPEC `x-.f` REAL_ABS_POS) x)));
913 MATCH_ACCEPT_TAC REAL_LE_TRANS]);;
915 let INTERVAL_EPS_0 = prove(
916 `!x f n. (interval_eps x f (float (int_of_num 0) n)) ==> (x = f)`,
918 REWRITE_TAC[interval_eps;float;int_of_num_th;REAL_MUL_LZERO];
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]);;
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]);;
928 (* ------------------------------------------------------------------ *)
929 (* INTERVAL DIVIDE *)
930 (* ------------------------------------------------------------------ *)
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`,
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));
946 SUBGOAL_MP_TAC `~(y= (&.0))`
950 REWRITE_TAC[interval_eps];
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);
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`;
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));
985 (* ------------------------------------------------------------------ *)
986 (* INTERVAL ABS VALUE *)
987 (* ------------------------------------------------------------------ *)
989 let INTERVAL_ABSV = prove(`!x f ex. interval_eps x f ex ==> (interval_eps (abs x) (abs f) ex)`,
991 REWRITE_TAC[interval_eps];
993 ASSUME_TAC (SPECL [`x:real`;`f:real`] REAL_ABS_SUB_ABS);
994 ASM_MESON_TAC[REAL_LE_TRANS]
997 (* ------------------------------------------------------------------ *)
998 (* INTERVAL ATN VALUE *)
999 (* ------------------------------------------------------------------ *)
1001 let taylor_atn = prove_by_refinement(
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`,
1008 REWRITE_TAC[interval_eps;float;TWOPOW_NEG;int_of_num_th;REAL_ARITH `&1 * x = x`];
1010 MATCH_MP_TAC real_taylor_atn_approx;
1012 EXISTS_TAC `inv(&2 pow (6 * n +5))`;
1013 EXISTS_TAC `eps2:real`;
1019 let pi_atn = prove_by_refinement(
1020 `pi = &4 * atn (&1)`,
1029 (* ------------------------------------------------------------------ *)
1031 (* This requires some preliminaries. Extend sqrt by 0 on negatives *)
1032 (* ------------------------------------------------------------------ *)
1034 let ssqrt = new_definition `ssqrt x = if (x <. (&.0)) then (&.0) else sqrt x`;; (*2m*)
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`,
1042 MESON_TAC[REAL_ARITH `&0 < x ==> ~(x < &0)`];
1046 let LET_TAC = REWRITE_TAC[LET_DEF;LET_END_DEF];;
1049 let REAL_SSQRT_NEG = prove(`!x. (x <. (&.0)) ==> (ssqrt x = (&.0))`,
1055 ACCEPT_TAC (REFL `&.0`);
1061 let REAL_SSQRT_NN = prove(`!x. (&.0) <=. x ==> (ssqrt x = (sqrt x))`,
1067 ASM_MESON_TAC[real_lt];
1068 ACCEPT_TAC (REFL `sqrt x`)
1070 ]);; (* 12 min, mostly spent loading *index-shell* *)
1074 let REAL_MK_NN_SSQRT = prove(`!x. (&.0) <=. (ssqrt x)`,
1077 DISJ_CASES_TAC (SPECL[`x:real`;`&.0`] REAL_LTE_TOTAL)
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]
1086 let REAL_SV_SSQRT_0 = prove(`!x. ssqrt (&.0) = (&.0)`,
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 *)
1095 let REAL_SV_SSQRT_n = prove_by_refinement(`!n. ssqrt (&n) = sqrt(&n)`,
1098 MATCH_MP_TAC REAL_SSQRT_NN;
1103 let REAL_SSQRT_EQ_0 = prove(`!(x:real). (ssqrt(x) = (&.0)) ==> (x <=. (&.0))`,
1106 DISJ_CASES_TAC (SPECL[`x:real`;`&.0`] REAL_LTE_TOTAL)
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]
1112 ]);; (* 15 minutes *)
1115 let REAL_SSQRT_MONO = prove(`!x. (x<=. y) ==> (ssqrt x <=. (ssqrt y))`,
1118 DISJ_CASES_TAC (SPECL[`x:real`;`&.0`] REAL_LTE_TOTAL)
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];
1123 ]);; (* 5 minutes *)
1125 let REAL_SSQRT_CHAR = prove(`!x t. (&.0 <=. t /\ (t*t = x)) ==> (t = (ssqrt x))`,
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 *)
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 *)
1135 let REAL_SSQRT_SQUARE' = prove(`!x. (&.0<=. x) ==> (ssqrt (x*.x) = x)`,
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*)
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)`,
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);
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);
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));
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]
1206 (* conversion for interval_eps *)
1208 (* ------------------------------------------------------------------ *)
1209 (* Take a term x of type real. Convert to a thm of the form *)
1210 (* interval_eps x f eps *)
1212 (* ------------------------------------------------------------------ *)
1214 let add_test,test = new_test_suite();;
1216 (* ------------------------------------------------------------------ *)
1218 Take the absolute value of input.
1219 Write it as a*2^k, where 1 <= a < 2, return k.
1222 num_exponent (Int 0) is -1.
1224 let (num_exponent:Num.num -> Num.num) =
1226 let afloat = float_of_num (abs_num a) in
1227 Int ((snd (frexp afloat)) - 1);;
1229 (*test*)let f (u,v) = ((num_exponent u) =(Int v)) in
1230 add_test("num_exponenwt",
1232 [Int 1,0; Int 65,6; Int (-65),6;
1233 Int 0,-1; (Int 3)//(Int 4),-1]);;
1234 (* ------------------------------------------------------------------ *)
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";;
1242 (* ------------------------------------------------------------------ *)
1245 (* finds a nearby (outward-rounded) Int with only prec_b significant bits *)
1246 let (round_outward: int -> Num.num -> Num.num) =
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;;
1254 let (round_inward: int-> Num.num -> Num.num) =
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;;
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);;
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);;
1274 let (round_outward_float: int -> float -> Num.num) =
1276 if (f=0.0) then (Int 0) else
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))))
1285 let (round_inward_float: int -> float -> Num.num) =
1287 if (f=0.0) then (Int 0) else
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))))
1298 (* ------------------------------------------------------------------ *)
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)
1305 add_test("SUBST_TERM",
1306 SUBST_TERM [(`#1`,`a:real`);(`#2`,`b:real`)] (`a +. b +. c`) =
1309 (* ------------------------------------------------------------------ *)
1311 (* take a term of the form `interval_eps x f ex` and clean up the f and ex *)
1313 let INTERVAL_CLEAN_CONV:conv =
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);;
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) &&
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))`
1335 (* ------------------------------------------------------------------ *)
1336 (* GENERAL lemmas *)
1337 (* ------------------------------------------------------------------ *)
1340 (* verifies statement of the form `float a b = float a' b'` *)
1344 let REAL_INV_POS = REAL_LT_INV;;
1346 let TWOPOW_MK_POS = prove(
1347 `!a. (&.0 <. ( twopow a))`,
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 ;
1358 let TWOPOW_NZ = prove(
1359 `!a. ~(twopow a = (&.0))`,
1361 ACCEPT_TAC (MATCH_MP REAL_MK_NZ_POS (SPEC `a:int` TWOPOW_MK_POS)));;
1363 let FLOAT_ZERO = prove(
1364 `!a b. (float a b = (&.0)) <=> (a = (int_of_num 0))`,
1366 REWRITE_TAC[float;REAL_ENTIRE;INT_OF_NUM_DEST];
1367 MESON_TAC[TWOPOW_NZ];
1370 let INT_ZERO = prove(
1371 `!n. ((int_of_num n = (int_of_num 0)) = (n=0))`,REWRITE_TAC[INT_OF_NUM_EQ]);;
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]);;
1377 let FLOAT_NN = Prove_by_refinement.prove_by_refinement(
1378 `!a b. ((&.0) <=. (float a b)) <=> (int_of_num 0 <=: a)`,
1380 REWRITE_TAC[float;INT_OF_NUM_DEST];
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];
1386 let INT_NN = INT_POS;;
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
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]);;
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);;
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);;
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];
1408 MESON_TAC[real_div]]);;
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]);;
1413 let ABS_NUM_P = prove_by_refinement(
1414 `!m n. (m <= n) ==> (abs (&. n -. (&. m)) = &.((m-|n) + (n-|m)))`,
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];
1425 let ABS_NUM = prove_by_refinement (
1426 `!m n. abs (&. n -. (&. m)) = &.((m-|n) + (n-|m))`,
1429 DISJ_CASES_TAC (SPECL [`m:num`;`n:num`] LTE_CASES) ; (* THENL[ *)
1431 MATCH_MP_TAC ABS_NUM_P;
1432 FIRST_X_ASSUM MP_TAC;
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;
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
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");
1454 H_MATCH_MP (THM REAL_LET_TRANS) (H_RULE2 CONJ (HYP "4") (HYP "5"));
1456 H_MATCH_MP (THM REAL_LTE_TRANS) (H_RULE2 CONJ (HYP "6") (HYP "3"));
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))))`,
1463 REWRITE_TAC[interval_eps];
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)`;
1475 ASM_MESON_TAC[REAL_LE_TRANS];
1479 (* end of general lemmas *)
1480 (* ------------------------------------------------------------------ *)
1483 (* ------------------------------------------------------------------ *)
1484 (* Cache of computed constants (abs (r - u) <= k)
1485 or interval r u eps *)
1486 (* ------------------------------------------------------------------ *)
1488 let cache_entry_of_ineq bprec ineq =
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
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"));;
1499 let get_real_constant cache bprec tm =
1500 (fun (_,_,r) -> r) (find (fun (t,p,th) -> (t = tm) & (p>= bprec)) cache);;
1502 let remove_entry t = filter (fun x -> not(x = t));;
1504 (* ------------------------------------------------------------------ *)
1506 (* take |- interval_eps r u e to interval_eps r' u e, where th_r gives r <=> r' *)
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
1519 (* term of the form '&.n'. Assume error checking done already. *)
1521 let INTERVAL_OF_NUM:conv =
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;;
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))`));;
1533 (* ------------------------------------------------------------------ *)
1534 (* conversion for float inequalities *)
1535 (* ------------------------------------------------------------------ *)
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;];;
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;;
1547 add_test("FLOAT_INEQ_CONV",
1549 dest_thm (c x) = ([],x) in
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)));;
1555 let INTERVAL_TO_LESS_CONV =
1556 let rthm = ASSUME `!f g ex ey. (&.0 <. (g -. (ey +. ex +. f)))` in
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;;
1567 add_test("INTERVAL_TO_LESS_CONV",
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)`);;
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)`));;
1579 (* ------------------------------------------------------------------ *)
1581 (* conversion for DEC <= posfloat and posfloat <= DEC *)
1584 `!n m p. ((&.p/(&.m)) <= (&.n)) <=> ((&.p/(&.m)) <= (&.n)/(&.1))`,
1585 MESON_TAC[REAL_DIV_1]);;
1588 `!n m p. ((&.p) <= ((&.n)/(&.m))) <=> ((&.p/(&.1)) <= (&.n)/(&.m))`,
1589 MESON_TAC[REAL_DIV_1]);;
1591 let REAL_LT = REAL_OF_NUM_LT;;
1593 let REAL_LE = REAL_OF_NUM_LE;;
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;
1600 ASM_SIMP_TAC[RAT_LEMMA4;REAL_LT;REAL_OF_NUM_MUL;REAL_LE]]);;
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];;
1607 add_test("DEC_FLOAT",
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)`)));;
1617 (* Conversion for powers *)
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];
1625 let pow_conv = PURE_ONCE_REWRITE_CONV[CONJUNCT1 real_pow;REAL_POW_1;pow_parity];;
1627 (* conversion for negation *)
1630 let th = REAL_ARITH `--x = &0 - x` in
1631 PURE_ONCE_REWRITE_CONV[th];;
1633 (* conversion for sums *)
1636 let sum = `sum:(num->bool)->(num->real)->real` in
1639 let (x,f) = dest_binop sum tm in
1640 let (i,j) = dest_binop rg x in
1644 let sum = `sum:(num->bool)->(num->real)->real` in
1647 let x = mk_binop rg i j in
1651 mk_sum (dest_sum `sum (3..4) f`);;
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
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;;
1671 (* ------------------------------------------------------------------ *)
1673 (* converts a DECIMAL TO A THEOREM *)
1675 let ABS_BOUNDS = REAL_ABS_BOUNDS;; (* new *)
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];
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);;
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`));;
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`)) ());;
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))`,
1712 REWRITE_TAC[interval_eps;REAL_ARITH `abs(x - x) = &0`];
1713 REWRITE_TAC[float;int_of_num_th];
1718 let INTERVAL_FLOAT_CONV flt =
1719 let (a,b) = dest_binop `float` flt in
1720 (ISPECL [a;b] INTERVAL_FLOAT);;
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 (*--------------------------------------------------------------------*)
1730 (* converts exceptions to false *)
1731 let falsify_ex f x = try (f x) with Failure _ -> false;;
1733 let is_bits_format =
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;;
1741 let is_numeral_format =
1743 let (h,t) = dest_comb x in
1744 ((h = `NUMERAL`) && (is_bits_format t)) in
1747 let is_decimal_format =
1749 let (t1,t2) = dest_binop `DECIMAL` x in
1750 ((is_numeral_format t1) && (is_numeral_format t2)) in
1753 let is_pos_int_format =
1755 let (h,t) = dest_comb x in
1756 (h = `int_of_num`) && (is_numeral_format t) in
1759 let is_neg_int_format =
1761 let (h,t) = dest_comb x in
1762 (h = `real_neg`) && (is_pos_int_format t) in
1765 let is_int_format x =
1766 (is_neg_int_format x) or (is_pos_int_format x);;
1768 let is_float_format =
1770 let (t1,t2) = dest_binop `float` x in
1771 (is_int_format t1) && (is_int_format t2) in
1774 let is_interval_format =
1776 let (a,b,c) = dest_interval x in
1777 (is_float_format b) && (is_float_format c) in
1782 let (h,t) = dest_comb x in
1786 let is_real_num_format =
1788 let (h,t) = dest_comb x in
1789 (h=`real_of_num`) && (is_numeral_format t) in
1792 let is_comb_of t u =
1794 t = (fst (dest_comb u)) in
1795 try (fn t u) with Failure _ -> false;;
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 *)
1802 let rec heron_sqrt depth A x eps =
1803 let half = (Int 1)//(Int 2) in
1805 then raise (Failure "heron_sqrt:internal error:sqrt recursion depth exceeded")
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;;
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
1817 (* compute ez parameter for division *)
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
1834 let pureconv (s,t) = (s,PURE_ONCE_REWRITE_CONV[t]);;
1836 (* ------------------------------------------------------------------ *)
1839 bprec gives the bits of precision.
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
1854 let interval_of_binop op opname match_thm =
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
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
1863 with Failure _ -> failwith ("INTERVAL_OF_TERM: "^opname)
1865 (* treat cached values first *)
1866 if (can (get_real_constant cache bprec) tm) then
1869 let int_thm = get_real_constant cache bprec tm in
1870 if (can dest_interval (concl int_thm)) then (cache,int_thm)
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)
1886 with Failure _ -> failwith "INTERVAL_OF_TERM : CONSTANT"
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)
1893 (* treat abs value *)
1894 else if (is_comb_of `real_abs` tm) then
1896 let (_,b) = dest_comb tm in
1897 let (cache,i1) = (INTERVAL_OF_TERM convl cache bprec b) in
1899 let b_int = MATCH_MP INTERVAL_ABSV i1 in
1900 let th = EQ_MP (INTERVAL_CLEAN_CONV (concl b_int)) b_int in
1902 with Failure _ -> failwith "INTERVAL_OF_TERM : ABS"
1905 else if (is_comb_of `twopow` tm) then
1908 let (_,b) = dest_comb tm in
1909 let th = SPEC b INTERVAL_OF_TWOPOW
1912 with Failure _ -> failwith "INTERVAL_OF_TERM : TWOPOW"
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
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
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
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
1944 with Failure _ -> failwith "INTERVAL_OF_TERM :DIV"
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
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
1956 else if can dest_sum tm then
1957 interval_of_conv cache bprec sum_conv tm
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
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
2004 with Failure _ -> failwith "INTERVAL_OF_TERM : SSQRT"
2006 else if (can (dest_unary `sqrt`) tm) then
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
2017 with Failure _ -> failwith "INTERVAL_OF_TERM : SQRT"
2020 else if can (dest_unary `atn`) tm then
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) )
2037 else if can (dest_unary `acs`) tm then
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
2049 else if (tm = `pi`) then
2051 interval_of_conv cache bprec (PURE_ONCE_REWRITE_CONV[pi_atn]) tm
2053 (* user supplied conversions *)
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
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);;
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;;
2080 (* Give it a conjunction of strict inequalities *)
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);;
2089 let REAL_INEQUALITY_CALCULATOR convl tm =
2090 REAL_INEQUALITY_CALCULATOR_PREC (10,70) convl tm;;
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);;
2103 let tmint = map (rec_real_ineq [] 10) tml in
2104 List.fold_left CONJ (hd tmint) (tl tmint);;