1 needs "arith/interval_arith.hl";;
4 module Constant_approximations = struct
10 let EPS_TO_INTERVAL = prove(`abs (f - x) <= e <=> interval_arith f (x - e, x + e)`,
11 REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;
13 let acs3_lo, acs3_hi =
14 let verify = MATCH_MP REAL_LT_IMP_LE o fst o M_verifier_main.verify_ineq M_verifier_main.default_params 10 in
15 verify `#1.230959417 < acs (&1 / &3)`, verify `acs(&1 / &3) < #1.230959418`;;
18 let tm = mk_comb(mk_comb(`DECIMAL`, mk_numeral (Num.abs_num a)), mk_numeral b) in
20 mk_comb (neg_op_real, tm)
25 let approx_interval th precision =
26 let th' = CONV_RULE (RAND_CONV (REWRITE_CONV[DECIMAL] THENC REAL_RAT_REDUCE_CONV)) th in
27 let lo_tm, hi_tm = dest_pair (rand(concl th')) in
28 let lo, hi = rat_of_term lo_tm, rat_of_term hi_tm in
29 let m = (Int 10 **/ Int precision) in
30 let lo_bound = floor_num (lo */ m) in
31 let hi_bound = ceiling_num (hi */ m) in
32 let conv = EQT_ELIM o REAL_RAT_LE_CONV in
33 let lo_th = conv (mk_binop le_op_real (mk_decimal lo_bound m) lo_tm) in
34 let hi_th = conv (mk_binop le_op_real hi_tm (mk_decimal hi_bound m)) in
35 let th0 = CONJ (CONJ lo_th hi_th) th' in
36 MATCH_MP APPROX_INTERVAL th0;;
40 (************************************)
43 let INTERVAL_SQRT = prove(`interval_arith x (a, b) /\
44 (c * c <= a /\ b <= d * d) ==>
45 interval_arith (sqrt x) (abs c, abs d)`,
46 REWRITE_TAC[interval_arith] THEN REPEAT STRIP_TAC THENL
48 MATCH_MP_TAC REAL_LE_RSQRT THEN
49 MATCH_MP_TAC REAL_LE_TRANS THEN
50 EXISTS_TAC `a:real` THEN
51 ASM_REWRITE_TAC[REAL_ARITH `abs a pow 2 = a * a`];
52 MATCH_MP_TAC REAL_LE_LSQRT THEN
53 ASM_REWRITE_TAC[REAL_ARITH `abs d pow 2 = d * d`; REAL_ABS_POS] THEN
56 MATCH_MP_TAC REAL_LE_TRANS THEN
57 EXISTS_TAC `a:real` THEN
58 ASM_REWRITE_TAC[] THEN
59 MATCH_MP_TAC REAL_LE_TRANS THEN
60 EXISTS_TAC `c * c:real` THEN
61 ASM_REWRITE_TAC[REAL_LE_SQUARE];
62 MATCH_MP_TAC REAL_LE_TRANS THEN
63 EXISTS_TAC `b:real` THEN
70 let interval_sqrt th precision =
71 let th' = CONV_RULE (REWRITE_CONV[DECIMAL] THENC REAL_RAT_REDUCE_CONV) th in
72 let lo, hi = dest_pair(rand(concl th')) in
73 let x_lo, x_hi = float_of_num (rat_of_term lo), float_of_num (rat_of_term hi) in
74 let lo_sqrt, hi_sqrt = Pervasives.sqrt x_lo, Pervasives.sqrt x_hi in
75 let m = 10.0 ** float_of_int precision in
76 let hack n = num_of_string (Int64.to_string (Int64.of_float n)) in
77 let sqrt_lo_num, sqrt_hi_num = hack (floor (lo_sqrt *. m)), hack (ceil (hi_sqrt *. m)) in
78 let m_num = Int 10 **/ Int precision in
79 let x_lo_tm = mk_decimal sqrt_lo_num m_num in
80 let x_hi_tm = mk_decimal sqrt_hi_num m_num in
81 let conv = EQT_ELIM o REAL_RAT_REDUCE_CONV in
82 let lo_th = conv (mk_binop le_op_real (mk_binop mul_op_real x_lo_tm x_lo_tm) lo) in
83 let hi_th = conv (mk_binop le_op_real hi (mk_binop mul_op_real x_hi_tm x_hi_tm)) in
84 let th0 = CONJ th' (CONJ lo_th hi_th) in
85 (CONV_RULE REAL_RAT_REDUCE_CONV) (MATCH_MP INTERVAL_SQRT th0);;
89 (************************************)
90 (* Arithmetic of intervals *)
92 let INTERVAL_ADD = prove(`interval_arith x (a, b) /\ interval_arith y (c, d)
93 ==> interval_arith (x + y) (a + c, b + d)`,
94 REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;
97 let INTERVAL_SUB = prove(`interval_arith x (a, b) /\ interval_arith y (c, d)
98 ==> interval_arith (x - y) (a - d, b - c)`,
99 REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;
102 let INTERVAL_NEG = prove(`interval_arith x (a, b) ==>
103 interval_arith (--x) (--b, --a)`,
104 REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;
108 let INTERVAL_INV = prove(`interval_arith x (a, b) /\ (&0 < a \/ b < &0)
109 ==> interval_arith (inv x) (inv b, inv a)`,
110 REWRITE_TAC[interval_arith] THEN
113 CONJ_TAC THEN MATCH_MP_TAC REAL_LE_INV2 THEN
114 ASM_REWRITE_TAC[] THEN
115 REPEAT (POP_ASSUM MP_TAC) THEN
119 ONCE_REWRITE_TAC[REAL_ARITH `a <= b <=> --b <= --a`] THEN
120 REWRITE_TAC[GSYM REAL_INV_NEG] THEN
121 CONJ_TAC THEN MATCH_MP_TAC REAL_LE_INV2 THEN
122 REPEAT (POP_ASSUM MP_TAC) THEN
126 let INTERVAL_INV_POS = prove(`interval_arith x (a,b) /\ &0 < a
127 ==> interval_arith (inv x) (inv b, inv a)`,
128 SIMP_TAC[INTERVAL_INV]);;
131 let INTERVAL_INV_NEG = prove(`interval_arith x (a,b) /\ b < &0
132 ==> interval_arith (inv x) (inv b, inv a)`,
133 SIMP_TAC[INTERVAL_INV]);;
137 let INTERVAL_MUL_lemma = prove(`!x y a b c d. interval_arith x (a, b) /\ interval_arith y (c, d) /\ x <= y
138 ==> x * y <= max (max (a * c) (a * d)) (max (b * c) (b * d))`,
140 REWRITE_TAC[interval_arith] THEN DISCH_TAC THEN
141 ABBREV_TAC `t = max (max (a * c) (a * d)) (max (b * c) (b * d))` THEN
142 SUBGOAL_THEN `a * c <= t /\ a * d <= t /\ b * c <= t /\ b * d <= t:real` ASSUME_TAC THENL
149 DISJ_CASES_TAC (REAL_ARITH `&0 <= x \/ x <= &0`) THENL
151 MATCH_MP_TAC REAL_LE_TRANS THEN
152 EXISTS_TAC `b * d:real` THEN
153 ASM_REWRITE_TAC[] THEN
154 MATCH_MP_TAC REAL_LE_MUL2 THEN
155 ASM_REWRITE_TAC[] THEN
156 MATCH_MP_TAC REAL_LE_TRANS THEN
157 EXISTS_TAC `x:real` THEN
162 DISJ_CASES_TAC (REAL_ARITH `&0 <= b \/ b <= &0`) THENL
164 DISJ_CASES_TAC (REAL_ARITH `&0 <= y \/ y <= &0`) THENL
166 MATCH_MP_TAC REAL_LE_TRANS THEN
170 ONCE_REWRITE_TAC[REAL_ARITH `&0 = &0 * y`] THEN
171 MATCH_MP_TAC REAL_LE_RMUL THEN
176 MATCH_MP_TAC REAL_LE_TRANS THEN
177 EXISTS_TAC `b * d:real` THEN
178 ASM_REWRITE_TAC[] THEN
179 MATCH_MP_TAC REAL_LE_MUL THEN
180 ASM_REWRITE_TAC[] THEN
181 MATCH_MP_TAC REAL_LE_TRANS THEN
182 EXISTS_TAC `y:real` THEN
187 MATCH_MP_TAC REAL_LE_TRANS THEN
188 EXISTS_TAC `a * c:real` THEN
189 ASM_REWRITE_TAC[] THEN
190 ONCE_REWRITE_TAC[GSYM REAL_NEG_MUL2] THEN
191 MATCH_MP_TAC REAL_LE_MUL2 THEN
192 ASM_REWRITE_TAC[REAL_LE_NEG; REAL_NEG_GE0];
196 DISJ_CASES_TAC (REAL_ARITH `&0 <= c \/ c <= &0`) THENL
198 MATCH_MP_TAC REAL_LE_TRANS THEN
199 EXISTS_TAC `b * c:real` THEN
200 ASM_REWRITE_TAC[] THEN
201 ONCE_REWRITE_TAC[REAL_ARITH `x * y <= b * c <=> (--b) * c <= (--x) * y`] THEN
202 MATCH_MP_TAC REAL_LE_MUL2 THEN
203 ASM_REWRITE_TAC[REAL_LE_NEG; REAL_NEG_GE0];
207 DISJ_CASES_TAC (REAL_ARITH `&0 <= y \/ y <= &0`) THENL
209 MATCH_MP_TAC REAL_LE_TRANS THEN
213 ONCE_REWRITE_TAC[REAL_ARITH `&0 = &0 * y`] THEN
214 MATCH_MP_TAC REAL_LE_RMUL THEN
219 MATCH_MP_TAC REAL_LE_TRANS THEN
220 EXISTS_TAC `a * c:real` THEN
221 ASM_REWRITE_TAC[] THEN
222 ONCE_REWRITE_TAC[GSYM REAL_NEG_MUL2] THEN
223 MATCH_MP_TAC REAL_LE_MUL THEN
224 ASM_REWRITE_TAC[REAL_NEG_GE0] THEN
225 MATCH_MP_TAC REAL_LE_TRANS THEN
226 EXISTS_TAC `x:real` THEN
231 MATCH_MP_TAC REAL_LE_TRANS THEN
232 EXISTS_TAC `a * c:real` THEN
233 ASM_REWRITE_TAC[] THEN
234 ONCE_REWRITE_TAC[GSYM REAL_NEG_MUL2] THEN
235 MATCH_MP_TAC REAL_LE_MUL2 THEN
236 ASM_REWRITE_TAC[REAL_LE_NEG; REAL_NEG_GE0]);;
240 let INTERVAL_MUL_lemma2 = prove(`!x y a b c d. interval_arith x (a,b) /\ interval_arith y (c,d)
241 ==> x * y <= max (max (a * c) (a * d)) (max (b * c) (b * d))`,
242 REPEAT STRIP_TAC THEN
243 DISJ_CASES_TAC (REAL_ARITH `x <= y \/ y <= x:real`) THENL
245 MATCH_MP_TAC INTERVAL_MUL_lemma THEN
250 MP_TAC (SPECL [`y:real`; `x:real`; `c:real`; `d:real`; `a:real`; `b:real`] INTERVAL_MUL_lemma) THEN
251 ASM_REWRITE_TAC[] THEN
257 let INTERVAL_MUL = prove(`interval_arith x (a, b) /\ interval_arith y (c, d)
258 ==> interval_arith (x * y)
259 (min (min (a * c) (a * d)) (min (b * c) (b * d)),
260 max (max (a * c) (a * d)) (max (b * c) (b * d)))`,
261 DISCH_TAC THEN REWRITE_TAC[interval_arith] THEN
262 ASM_SIMP_TAC[INTERVAL_MUL_lemma2] THEN
263 MP_TAC (SPECL[`--x:real`; `y:real`; `--b:real`; `--a:real`; `c:real`; `d:real`] INTERVAL_MUL_lemma2) THEN
264 ASM_SIMP_TAC[INTERVAL_NEG] THEN
269 (**************************************)
270 let const_interval tm = SPEC tm CONST_INTERVAL;;
272 let interval_neg th = MATCH_MP INTERVAL_NEG th;;
274 let interval_add th1 th2 =
275 let th0 = MATCH_MP INTERVAL_ADD (CONJ th1 th2) in
276 (CONV_RULE (RAND_CONV REAL_RAT_REDUCE_CONV)) th0;;
278 let interval_sub th1 th2 =
279 let th0 = MATCH_MP INTERVAL_SUB (CONJ th1 th2) in
280 (CONV_RULE (RAND_CONV REAL_RAT_REDUCE_CONV)) th0;;
282 let interval_mul th1 th2 =
283 let th0 = MATCH_MP INTERVAL_MUL (CONJ th1 th2) in
284 (CONV_RULE (RAND_CONV REAL_RAT_REDUCE_CONV)) th0;;
286 let interval_inv th =
287 let lt_op_real = `(<):real->real->bool` in
288 let lo_tm, hi_tm = dest_pair(rand(concl th)) in
289 let lo_ineq = REAL_RAT_LT_CONV (mk_binop lt_op_real `&0` lo_tm) in
290 if (rand(concl lo_ineq) = `T`) then
291 let th0 = CONJ th (EQT_ELIM lo_ineq) in
292 (CONV_RULE (RAND_CONV REAL_RAT_REDUCE_CONV)) (MATCH_MP INTERVAL_INV_POS th0)
294 let hi_ineq = REAL_RAT_LT_CONV (mk_binop lt_op_real hi_tm `&0`) in
295 if (rand(concl hi_ineq) = `T`) then
296 let th0 = CONJ th (EQT_ELIM hi_ineq) in
297 (CONV_RULE (RAND_CONV REAL_RAT_REDUCE_CONV)) (MATCH_MP INTERVAL_INV_NEG th0)
298 else failwith("interval_inv: 0 is inside interval");;
300 let interval_div th1 th2 =
301 let th2' = interval_inv th2 in
302 ONCE_REWRITE_RULE[GSYM real_div] (interval_mul th1 th2');;
306 (*************************)
307 let acs3_interval = REWRITE_RULE[GSYM interval_arith] (CONJ acs3_lo acs3_hi);;
309 let pi_interval = prove(`interval_arith pi (#3.141592653, #3.141592654)`,
310 REWRITE_TAC[interval_arith] THEN
311 MP_TAC PI_APPROX_32 THEN REAL_ARITH_TAC);;
313 let tgt_interval = prove(`interval_arith tgt (#1.541, #1.541)`,
314 REWRITE_TAC[Tame_defs.tgt; interval_arith; REAL_LE_REFL]);;
317 let interval_table = Hashtbl.create 10;;
318 let add_interval th = Hashtbl.add interval_table ((rand o rator o concl) th) th;;
321 let rec create_interval tm =
322 if Hashtbl.mem interval_table tm then
323 Hashtbl.find interval_table tm
325 if (is_ratconst tm) then
327 else if (is_binop add_op_real tm) then
328 let lhs, rhs = dest_binop add_op_real tm in
329 interval_add (create_interval lhs) (create_interval rhs)
330 else if (is_binop sub_op_real tm) then
331 let lhs, rhs = dest_binop sub_op_real tm in
332 interval_sub (create_interval lhs) (create_interval rhs)
333 else if (is_binop mul_op_real tm) then
334 let lhs, rhs = dest_binop mul_op_real tm in
335 interval_mul (create_interval lhs) (create_interval rhs)
336 else if (is_binop div_op_real tm) then
337 let lhs, rhs = dest_binop div_op_real tm in
338 interval_div (create_interval lhs) (create_interval rhs)
339 else if (is_comb tm) then
340 let ltm, rtm = dest_comb tm in
341 if (ltm = inv_op_real) then
342 interval_inv (create_interval rtm)
343 else if (ltm = neg_op_real) then
344 interval_neg (create_interval rtm)
345 else failwith "create_interval: unknown unary operation"
347 failwith "create_interval: unexpected term";;
352 add_interval pi_interval;;
353 add_interval acs3_interval;;
354 add_interval tgt_interval;;
355 add_interval (REWRITE_RULE[GSYM sqrt8] (interval_sqrt (const_interval `&8`) 9));;
356 add_interval (REWRITE_RULE[GSYM sol0] (create_interval `&3 * acs(&1 / &3) - pi`));;
357 add_interval (create_interval `sol0 / pi`);;
359 let rho218 = new_definition `rho218 = rho(#2.18)`;;
361 let rho218_def = (REWRITE_CONV[rho218; rho; ly; interp; GSYM Tame_general.sol0_over_pi_EQ_const1] THENC
362 REAL_RAT_REDUCE_CONV) `rho218`;;
364 let rho218_interval = REWRITE_RULE[SYM rho218_def] (create_interval(rand(concl rho218_def)));;
365 add_interval rho218_interval;;