formal_lp/ineqs/constants_approx.hl

   1 (* Basic interval arithmetic for approximation of constants *)
   2 needs "tame/tame_defs.hl";;
   3 needs "tame/TameGeneral.hl";;
   4
   5 module Constants_approx = struct
   6
   7 let interval_arith = new_definition `interval_arith (x:real) (lo, hi) <=> lo <= x /\ x <= hi`;;
   8
   9
  10 let CONST_INTERVAL = prove(`!x. interval_arith x (x,x)`, REWRITE_TAC[interval_arith; REAL_LE_REFL]);;
  11
  12
  13 let EPS_TO_INTERVAL = prove(`abs (f - x) <= e <=> interval_arith f (x - e, x + e)`,
  14    REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;
  15
  16
  17 let APPROX_INTERVAL = prove(`(a <= lo /\ hi <= b) /\ interval_arith x (lo, hi)
  18                               ==> interval_arith x (a,b)`,
  19    REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;
  20
  21
  22
  23 (* Numerical approximations for cos and acos *)
  24
  25
  26 let COS_EQ_NEG_SIN = prove(`!x. cos (x + pi / &2) = --sin x`,
  27      REWRITE_TAC[COS_SIN; REAL_ARITH `a - (b + a) = --b`; SIN_NEG]);;
  28
  29
  30
  31
  32 let COS_DERIVATIVES = prove(`!x n. ((\x. cos (x + &n * pi / &2)) has_real_derivative cos (x + &(n + 1) * pi / &2)) (atreal x)`,
  33    REPEAT GEN_TAC THEN REWRITE_TAC[] THEN
  34      MP_TAC (REAL_DIFF_CONV `((\x. cos (x + &n * pi / &2)) has_real_derivative f) (atreal x)`) THEN
  35      SUBGOAL_THEN `(&1 + &0) * --sin (x + &n * pi / &2) = cos (x + &(n + 1) * pi / &2)` (fun th -> REWRITE_TAC[th]) THEN
  36      REWRITE_TAC[REAL_ARITH `(&1 + &0) * --a = --a`] THEN
  37      REWRITE_TAC[GSYM REAL_OF_NUM_ADD] THEN
  38      REWRITE_TAC[REAL_ARITH `x + (a + &1) * t = (x + a * t) + t`] THEN
  39      REWRITE_TAC[COS_EQ_NEG_SIN]);;
  40
  41
  42
  43
  44 let REAL_TAYLOR_COS_RAW = prove(`!x n. abs (cos x - sum (0..n) (\k. if (EVEN k) then ((-- &1) pow (k DIV 2) * x pow k) / &(FACT k) else &0)) <=
  45                                   abs x pow (n + 1) / &(FACT (n + 1))`,
  46    REPEAT GEN_TAC THEN
  47      MP_TAC (SPECL [`(\i x. cos (x + &i * pi / &2))`; `n:num`; `UNIV:real->bool`; `&1`] REAL_TAYLOR) THEN
  48      ANTS_TAC THENL
  49      [
  50        REWRITE_TAC[is_realinterval; IN_UNIV; WITHINREAL_UNIV; COS_DERIVATIVES; COS_BOUND];
  51        ALL_TAC
  52      ] THEN
  53      REWRITE_TAC[IN_UNIV] THEN
  54      DISCH_THEN (MP_TAC o SPECL [`&0`; `x:real`]) THEN
  55      REWRITE_TAC[REAL_MUL_LZERO; REAL_ADD_RID; REAL_ADD_LID; REAL_SUB_RZERO; REAL_MUL_LID] THEN
  56      SUBGOAL_THEN `!i. cos (&i * pi / &2) * x pow i / &(FACT i) = if EVEN i then (-- &1 pow (i DIV 2) * x pow i) / &(FACT i) else &0` (fun th -> SIMP_TAC[th]) THEN
  57      GEN_TAC THEN
  58      ASM_CASES_TAC `EVEN i` THEN ASM_REWRITE_TAC[] THENL
  59      [
  60        POP_ASSUM MP_TAC THEN
  61          REWRITE_TAC[EVEN_EXISTS] THEN
  62          STRIP_TAC THEN ASM_REWRITE_TAC[] THEN
  63          SUBGOAL_THEN `(2 * m) DIV 2 = m` (fun th -> REWRITE_TAC[th]) THENL
  64          [
  65            MATCH_MP_TAC DIV_MULT THEN
  66              ARITH_TAC;
  67            ALL_TAC
  68          ] THEN
  69          REWRITE_TAC[GSYM REAL_OF_NUM_MUL] THEN
  70          REWRITE_TAC[REAL_ARITH `(&2 * a) * b / &2 = a * b`] THEN
  71          REWRITE_TAC[COS_NPI] THEN
  72          REAL_ARITH_TAC;
  73        ALL_TAC
  74      ] THEN
  75      SUBGOAL_THEN `cos (&i * pi / &2) = &0` (fun th -> REWRITE_TAC[th]) THENL
  76      [
  77        REWRITE_TAC[COS_ZERO] THEN
  78          DISJ1_TAC THEN EXISTS_TAC `i:num` THEN
  79          ASM_REWRITE_TAC[];
  80        ALL_TAC
  81      ] THEN
  82      REWRITE_TAC[REAL_MUL_LZERO]);;
  83
  84
  85 let SUM_PAIR_0 = prove(`!f n. sum (0..2 * n + 1) f = sum(0..n) (\i. f (2 * i) + f (2 * i + 1))`,
  86    REPEAT GEN_TAC THEN
  87      MP_TAC (SPECL [`f:num->real`; `0`; `n:num`] SUM_PAIR) THEN
  88      REWRITE_TAC[ARITH_RULE `2 * 0 = 0`]);;
  89
  90
  91 let REAL_TAYLOR_COS = prove(`!x n. abs (cos x - sum (0..n) (\i. (-- &1) pow i * x pow (2 * i) / &(FACT (2 * i)))) <= abs x pow (2*n + 2) / &(FACT (2*n + 2))`,
  92    REPEAT GEN_TAC THEN
  93      MP_TAC (SPECL [`x:real`; `2 * n + 1`] REAL_TAYLOR_COS_RAW) THEN
  94      REWRITE_TAC[SUM_PAIR_0; EVEN_DOUBLE; ARITH_RULE `(2 * n + 1) + 1 = 2 *n + 2`] THEN
  95      SUBGOAL_THEN `!i. ~(EVEN (2 * i + 1))` MP_TAC THENL
  96      [
  97        GEN_TAC THEN REWRITE_TAC[NOT_EVEN] THEN
  98          REWRITE_TAC[ARITH_ODD; ODD_ADD; ODD_MULT];
  99        ALL_TAC
 100      ] THEN
 101      DISCH_THEN (fun th -> SIMP_TAC[th]) THEN
 102      SUBGOAL_THEN `!i. (2 * i) DIV 2 = i` MP_TAC THENL
 103      [
 104        GEN_TAC THEN
 105          MATCH_MP_TAC DIV_MULT THEN
 106          REWRITE_TAC[ARITH];
 107        ALL_TAC
 108      ] THEN
 109      DISCH_THEN (fun th -> REWRITE_TAC[th; REAL_ADD_RID]) THEN
 110      REWRITE_TAC[REAL_ARITH `(a * b) / c = a * b / c`]);;
 111
 112
 113
 114
 115
 116 let gen_sum_thm n =
 117   let SUM_lemma n =
 118     let tm = list_mk_comb (`sum:(num->bool)->(num->real)->real`, [mk_comb (`(..) 0`, mk_small_numeral n); `f:num->real`]) in
 119     let suc_th = REWRITE_RULE[EQ_CLAUSES] (REWRITE_CONV[ARITH] (mk_eq (mk_small_numeral n, mk_comb (`SUC`, mk_small_numeral (n - 1))))) in
 120     let th1 = REWRITE_CONV[suc_th] tm in
 121       REWRITE_RULE[SUM_CLAUSES_NUMSEG; ARITH] th1 in
 122   let rec rewriter th n =
 123     if n > 0 then
 124       rewriter (REWRITE_RULE[SUM_lemma n; GSYM REAL_ADD_ASSOC] th) (n - 1)
 125     else
 126       th in
 127     GEN_ALL (rewriter (SUM_lemma n) (n - 1));;
 128
 129
 130 (* Evaluates cos at a given point using n terms from the cosine Taylor series *)
 131 let cos_eval x n =
 132   let th1 = (SPECL [x; mk_small_numeral n] REAL_TAYLOR_COS) in
 133   let th2 = REWRITE_RULE[gen_sum_thm n] th1 in
 134   let th4 = CONV_RULE (NUM_REDUCE_CONV) th2 in
 135   let th5 = CONV_RULE (DEPTH_CONV REAL_INT_POW_CONV) th4 in
 136     CONV_RULE (REAL_RAT_REDUCE_CONV) th5;;
 137
 138
 139
 140 let acs3_lo = prove(`#1.230959417 <= acs (&1 / &3)`,
 141    SUBGOAL_THEN `#1.230959417 = acs (cos(#1.230959417))` (fun th -> ONCE_REWRITE_TAC[th]) THENL
 142    [
 143      MATCH_MP_TAC (GSYM ACS_COS) THEN
 144        MP_TAC PI_APPROX_32 THEN
 145        REAL_ARITH_TAC;
 146      ALL_TAC
 147    ] THEN
 148    MATCH_MP_TAC ACS_MONO_LE THEN
 149    REWRITE_TAC[COS_BOUNDS] THEN
 150    MP_TAC (cos_eval `#1.230959417` 6) THEN
 151    REAL_ARITH_TAC);;
 152
 153
 154 (* 1.23095941734077 *)
 155 let acs3_hi = prove(`acs(&1 / &3) <= #1.230959418`,
 156    SUBGOAL_THEN `#1.230959418 = acs(cos(#1.230959418))` (fun th -> ONCE_REWRITE_TAC[th]) THENL
 157    [
 158      MATCH_MP_TAC (GSYM ACS_COS) THEN
 159        MP_TAC PI_APPROX_32 THEN
 160        REAL_ARITH_TAC;
 161      ALL_TAC
 162    ] THEN
 163    MATCH_MP_TAC ACS_MONO_LE THEN
 164    REWRITE_TAC[COS_BOUNDS] THEN
 165    MP_TAC (cos_eval `#1.230959418` 6) THEN
 166    REAL_ARITH_TAC);;
 167
 168
 169 let le_op_real = `(<=):real->real->bool` and
 170     minus_op_real = `(-):real->real->real` and
 171     plus_op_real = `(+):real->real->real` and
 172     mul_op_real = `( * ):real->real->real` and
 173     div_op_real = `(/):real->real->real` and
 174     inv_op_real = `inv:real->real` and
 175     neg_op_real = `(--):real->real`;;
 176
 177 let mk_decimal a b =
 178   let tm = mk_comb(mk_comb(`DECIMAL`, mk_numeral (Num.abs_num a)), mk_numeral b) in
 179     if (a </ Int 0) then
 180       mk_comb (neg_op_real, tm)
 181     else
 182       tm;;
 183
 184
 185
 186 let approx_interval th precision =
 187   let th' = CONV_RULE (RAND_CONV (REWRITE_CONV[DECIMAL] THENC REAL_RAT_REDUCE_CONV)) th in
 188   let lo_tm, hi_tm = dest_pair (rand(concl th')) in
 189   let lo, hi = rat_of_term lo_tm, rat_of_term hi_tm in
 190   let m = (Int 10 **/ Int precision) in
 191   let lo_bound = floor_num (lo */ m) in
 192   let hi_bound = ceiling_num (hi */ m) in
 193   let conv = EQT_ELIM o REAL_RAT_LE_CONV in
 194   let lo_th = conv (mk_binop le_op_real (mk_decimal lo_bound m) lo_tm) in
 195   let hi_th = conv (mk_binop le_op_real hi_tm (mk_decimal hi_bound m)) in
 196   let th0 = CONJ (CONJ lo_th hi_th) th' in
 197     MATCH_MP APPROX_INTERVAL th0;;
 198
 199
 200
 201 (*
 202 approx_interval th1 9;;
 203
 204 let th = cos_eval `#0.61547970867` 5;;
 205 let th1 = (CONV_RULE REAL_RAT_REDUCE_CONV) (REWRITE_RULE[EPS_TO_INTERVAL] th);;
 206 float_of_num (rat_of_term (rand(concl th)));;
 207
 208 approx_interval (concl th1) 10;;
 209 *)
 210
 211
 212
 213 (************************************)
 214
 215 (* Square root *)
 216
 217
 218 let INTERVAL_SQRT = prove(`interval_arith x (a, b) /\
 219                             (c * c <= a /\ b <= d * d) ==>
 220                             interval_arith (sqrt x) (abs c, abs d)`,
 221    REWRITE_TAC[interval_arith] THEN REPEAT STRIP_TAC THENL
 222      [
 223        MATCH_MP_TAC REAL_LE_RSQRT THEN
 224          MATCH_MP_TAC REAL_LE_TRANS THEN
 225          EXISTS_TAC `a:real` THEN
 226          ASM_REWRITE_TAC[REAL_ARITH `abs a pow 2 = a * a`];
 227        MATCH_MP_TAC REAL_LE_LSQRT THEN
 228          ASM_REWRITE_TAC[REAL_ARITH `abs d pow 2 = d * d`; REAL_ABS_POS] THEN
 229          CONJ_TAC THENL
 230          [
 231            MATCH_MP_TAC REAL_LE_TRANS THEN
 232              EXISTS_TAC `a:real` THEN
 233              ASM_REWRITE_TAC[] THEN
 234              MATCH_MP_TAC REAL_LE_TRANS THEN
 235              EXISTS_TAC `c * c:real` THEN
 236              ASM_REWRITE_TAC[REAL_LE_SQUARE];
 237            MATCH_MP_TAC REAL_LE_TRANS THEN
 238              EXISTS_TAC `b:real` THEN
 239              ASM_REWRITE_TAC[]
 240          ]
 241      ]);;
 242
 243
 244
 245 let interval_sqrt th precision =
 246   let th' = CONV_RULE (REWRITE_CONV[DECIMAL] THENC REAL_RAT_REDUCE_CONV) th in
 247   let lo, hi = dest_pair(rand(concl th')) in
 248   let x_lo, x_hi = float_of_num (rat_of_term lo), float_of_num (rat_of_term hi) in
 249   let lo_sqrt, hi_sqrt = Pervasives.sqrt x_lo, Pervasives.sqrt x_hi in
 250   let m = 10.0 ** float_of_int precision in
 251   let hack n = num_of_string (Int64.to_string (Int64.of_float n)) in
 252   let sqrt_lo_num, sqrt_hi_num = hack (floor (lo_sqrt *. m)), hack (ceil (hi_sqrt *. m)) in
 253   let m_num = Int 10 **/ Int precision in
 254   let x_lo_tm = mk_decimal sqrt_lo_num m_num in
 255   let x_hi_tm = mk_decimal sqrt_hi_num m_num in
 256   let conv = EQT_ELIM o REAL_RAT_REDUCE_CONV in
 257   let lo_th = conv (mk_binop le_op_real (mk_binop mul_op_real x_lo_tm x_lo_tm) lo) in
 258   let hi_th = conv (mk_binop le_op_real hi (mk_binop mul_op_real x_hi_tm x_hi_tm)) in
 259   let th0 = CONJ th' (CONJ lo_th hi_th) in
 260     (CONV_RULE REAL_RAT_REDUCE_CONV) (MATCH_MP INTERVAL_SQRT th0);;
 261
 262
 263
 264 (************************************)
 265
 266 (* Arithmetic of intervals *)
 267
 268 let INTERVAL_ADD = prove(`interval_arith x (a, b) /\ interval_arith y (c, d)
 269                            ==> interval_arith (x + y) (a + c, b + d)`,
 270    REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;
 271
 272
 273 let INTERVAL_SUB = prove(`interval_arith x (a, b) /\ interval_arith y (c, d)
 274                           ==> interval_arith (x - y) (a - d, b - c)`,
 275    REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;
 276
 277
 278 let INTERVAL_NEG = prove(`interval_arith x (a, b) ==>
 279                            interval_arith (--x) (--b, --a)`,
 280    REWRITE_TAC[interval_arith] THEN REAL_ARITH_TAC);;
 281
 282
 283
 284 let INTERVAL_INV = prove(`interval_arith x (a, b) /\ (&0 < a \/ b < &0)
 285                            ==> interval_arith (inv x) (inv b, inv a)`,
 286    REWRITE_TAC[interval_arith] THEN
 287      STRIP_TAC THENL
 288      [
 289        CONJ_TAC THEN MATCH_MP_TAC REAL_LE_INV2 THEN
 290          ASM_REWRITE_TAC[] THEN
 291          REPEAT (POP_ASSUM MP_TAC) THEN
 292          REAL_ARITH_TAC;
 293        ALL_TAC
 294      ] THEN
 295      ONCE_REWRITE_TAC[REAL_ARITH `a <= b <=> --b <= --a`] THEN
 296      REWRITE_TAC[GSYM REAL_INV_NEG] THEN
 297      CONJ_TAC THEN MATCH_MP_TAC REAL_LE_INV2 THEN
 298      REPEAT (POP_ASSUM MP_TAC) THEN
 299      REAL_ARITH_TAC);;
 300
 301
 302 let INTERVAL_INV_POS = prove(`interval_arith x (a,b) /\ &0 < a
 303                                ==> interval_arith (inv x) (inv b, inv a)`,
 304                              SIMP_TAC[INTERVAL_INV]);;
 305
 306
 307 let INTERVAL_INV_NEG = prove(`interval_arith x (a,b) /\ b < &0
 308                                ==> interval_arith (inv x) (inv b, inv a)`,
 309                              SIMP_TAC[INTERVAL_INV]);;
 310
 311
 312
 313 let INTERVAL_MUL_lemma = prove(`!x y a b c d. interval_arith x (a, b) /\ interval_arith y (c, d) /\ x <= y
 314                                  ==> x * y <= max (max (a * c) (a * d)) (max (b * c) (b * d))`,
 315    REPEAT GEN_TAC THEN
 316      REWRITE_TAC[interval_arith] THEN DISCH_TAC THEN
 317      ABBREV_TAC `t = max (max (a * c) (a * d)) (max (b * c) (b * d))` THEN
 318      SUBGOAL_THEN `a * c <= t /\ a * d <= t /\ b * c <= t /\ b * d <= t:real` ASSUME_TAC THENL
 319      [
 320        EXPAND_TAC "t" THEN
 321          REAL_ARITH_TAC;
 322        ALL_TAC
 323      ] THEN
 324
 325      DISJ_CASES_TAC (REAL_ARITH `&0 <= x \/ x <= &0`) THENL
 326      [
 327        MATCH_MP_TAC REAL_LE_TRANS THEN
 328          EXISTS_TAC `b * d:real` THEN
 329          ASM_REWRITE_TAC[] THEN
 330          MATCH_MP_TAC REAL_LE_MUL2 THEN
 331          ASM_REWRITE_TAC[] THEN
 332          MATCH_MP_TAC REAL_LE_TRANS THEN
 333          EXISTS_TAC `x:real` THEN
 334          ASM_REWRITE_TAC[];
 335        ALL_TAC
 336      ] THEN
 337
 338      DISJ_CASES_TAC (REAL_ARITH `&0 <= b \/ b <= &0`) THENL
 339      [
 340        DISJ_CASES_TAC (REAL_ARITH `&0 <= y \/ y <= &0`) THENL
 341          [
 342            MATCH_MP_TAC REAL_LE_TRANS THEN
 343              EXISTS_TAC `&0` THEN
 344              CONJ_TAC THENL
 345              [
 346                ONCE_REWRITE_TAC[REAL_ARITH `&0 = &0 * y`] THEN
 347                  MATCH_MP_TAC REAL_LE_RMUL THEN
 348                  ASM_REWRITE_TAC[];
 349                ALL_TAC
 350              ] THEN
 351
 352              MATCH_MP_TAC REAL_LE_TRANS THEN
 353              EXISTS_TAC `b * d:real` THEN
 354              ASM_REWRITE_TAC[] THEN
 355              MATCH_MP_TAC REAL_LE_MUL THEN
 356              ASM_REWRITE_TAC[] THEN
 357              MATCH_MP_TAC REAL_LE_TRANS THEN
 358              EXISTS_TAC `y:real` THEN
 359              ASM_REWRITE_TAC[];
 360            ALL_TAC
 361          ] THEN
 362
 363          MATCH_MP_TAC REAL_LE_TRANS THEN
 364          EXISTS_TAC `a * c:real` THEN
 365          ASM_REWRITE_TAC[] THEN
 366          ONCE_REWRITE_TAC[GSYM REAL_NEG_MUL2] THEN
 367          MATCH_MP_TAC REAL_LE_MUL2 THEN
 368          ASM_REWRITE_TAC[REAL_LE_NEG; REAL_NEG_GE0];
 369        ALL_TAC
 370      ] THEN
 371
 372      DISJ_CASES_TAC (REAL_ARITH `&0 <= c \/ c <= &0`) THENL
 373      [
 374        MATCH_MP_TAC REAL_LE_TRANS THEN
 375          EXISTS_TAC `b * c:real` THEN
 376          ASM_REWRITE_TAC[] THEN
 377          ONCE_REWRITE_TAC[REAL_ARITH `x * y <= b * c <=> (--b) * c <= (--x) * y`] THEN
 378          MATCH_MP_TAC REAL_LE_MUL2 THEN
 379          ASM_REWRITE_TAC[REAL_LE_NEG; REAL_NEG_GE0];
 380        ALL_TAC
 381      ] THEN
 382
 383      DISJ_CASES_TAC (REAL_ARITH `&0 <= y \/ y <= &0`) THENL
 384      [
 385        MATCH_MP_TAC REAL_LE_TRANS THEN
 386          EXISTS_TAC `&0` THEN
 387          CONJ_TAC THENL
 388          [
 389            ONCE_REWRITE_TAC[REAL_ARITH `&0 = &0 * y`] THEN
 390              MATCH_MP_TAC REAL_LE_RMUL THEN
 391              ASM_REWRITE_TAC[];
 392            ALL_TAC
 393          ] THEN
 394
 395          MATCH_MP_TAC REAL_LE_TRANS THEN
 396          EXISTS_TAC `a * c:real` THEN
 397          ASM_REWRITE_TAC[] THEN
 398          ONCE_REWRITE_TAC[GSYM REAL_NEG_MUL2] THEN
 399          MATCH_MP_TAC REAL_LE_MUL THEN
 400          ASM_REWRITE_TAC[REAL_NEG_GE0] THEN
 401          MATCH_MP_TAC REAL_LE_TRANS THEN
 402          EXISTS_TAC `x:real` THEN
 403          ASM_REWRITE_TAC[];
 404        ALL_TAC
 405      ] THEN
 406
 407      MATCH_MP_TAC REAL_LE_TRANS THEN
 408      EXISTS_TAC `a * c:real` THEN
 409      ASM_REWRITE_TAC[] THEN
 410      ONCE_REWRITE_TAC[GSYM REAL_NEG_MUL2] THEN
 411      MATCH_MP_TAC REAL_LE_MUL2 THEN
 412      ASM_REWRITE_TAC[REAL_LE_NEG; REAL_NEG_GE0]);;
 413
 414
 415
 416 let INTERVAL_MUL_lemma2 = prove(`!x y a b c d. interval_arith x (a,b) /\ interval_arith y (c,d)
 417                                   ==> x * y <= max (max (a * c) (a * d)) (max (b * c) (b * d))`,
 418    REPEAT STRIP_TAC THEN
 419      DISJ_CASES_TAC (REAL_ARITH `x <= y \/ y <= x:real`) THENL
 420      [
 421        MATCH_MP_TAC INTERVAL_MUL_lemma THEN
 422          ASM_REWRITE_TAC[];
 423        ALL_TAC
 424      ] THEN
 425
 426      MP_TAC (SPECL [`y:real`; `x:real`; `c:real`; `d:real`; `a:real`; `b:real`] INTERVAL_MUL_lemma) THEN
 427      ASM_REWRITE_TAC[] THEN
 428      REAL_ARITH_TAC);;
 429
 430
 431
 432
 433 let INTERVAL_MUL = prove(`interval_arith x (a, b) /\ interval_arith y (c, d)
 434                            ==> interval_arith (x * y)
 435                            (min (min (a * c) (a * d)) (min (b * c) (b * d)),
 436                             max (max (a * c) (a * d)) (max (b * c) (b * d)))`,
 437    DISCH_TAC THEN REWRITE_TAC[interval_arith] THEN
 438      ASM_SIMP_TAC[INTERVAL_MUL_lemma2] THEN
 439      MP_TAC (SPECL[`--x:real`; `y:real`; `--b:real`; `--a:real`; `c:real`; `d:real`] INTERVAL_MUL_lemma2) THEN
 440      ASM_SIMP_TAC[INTERVAL_NEG] THEN
 441      REAL_ARITH_TAC);;
 442
 443
 444
 445
 446
 447
 448 (**************************************)
 449
 450
 451 let const_interval tm = SPEC tm CONST_INTERVAL;;
 452
 453
 454 let interval_neg th = MATCH_MP INTERVAL_NEG th;;
 455
 456
 457 let interval_add th1 th2 =
 458   let th0 = MATCH_MP INTERVAL_ADD (CONJ th1 th2) in
 459     (CONV_RULE (RAND_CONV REAL_RAT_REDUCE_CONV)) th0;;
 460
 461 let interval_sub th1 th2 =
 462   let th0 = MATCH_MP INTERVAL_SUB (CONJ th1 th2) in
 463     (CONV_RULE (RAND_CONV REAL_RAT_REDUCE_CONV)) th0;;
 464
 465
 466 let interval_mul th1 th2 =
 467   let th0 = MATCH_MP INTERVAL_MUL (CONJ th1 th2) in
 468     (CONV_RULE (RAND_CONV REAL_RAT_REDUCE_CONV)) th0;;
 469
 470
 471 let interval_inv th =
 472   let lt_op_real = `(<):real->real->bool` in
 473   let lo_tm, hi_tm = dest_pair(rand(concl th)) in
 474   let lo_ineq = REAL_RAT_LT_CONV (mk_binop lt_op_real `&0` lo_tm) in
 475     if (rand(concl lo_ineq) = `T`) then
 476       let th0 = CONJ th (EQT_ELIM lo_ineq) in
 477         (CONV_RULE (RAND_CONV REAL_RAT_REDUCE_CONV)) (MATCH_MP INTERVAL_INV_POS th0)
 478     else
 479       let hi_ineq = REAL_RAT_LT_CONV (mk_binop lt_op_real hi_tm `&0`) in
 480         if (rand(concl hi_ineq) = `T`) then
 481           let th0 = CONJ th (EQT_ELIM hi_ineq) in
 482             (CONV_RULE (RAND_CONV REAL_RAT_REDUCE_CONV)) (MATCH_MP INTERVAL_INV_NEG th0)
 483         else failwith("interval_inv: 0 is inside interval");;
 484
 485 let interval_div th1 th2 =
 486   let th2' = interval_inv th2 in
 487     ONCE_REWRITE_RULE[GSYM real_div] (interval_mul th1 th2');;
 488
 489
 490
 491 (*************************)
 492
 493
 494 let acs3_interval = REWRITE_RULE[GSYM interval_arith] (CONJ acs3_lo acs3_hi);;
 495
 496
 497 let pi_interval = prove(`interval_arith pi (#3.141592653, #3.141592654)`,
 498    REWRITE_TAC[interval_arith] THEN
 499      MP_TAC PI_APPROX_32 THEN REAL_ARITH_TAC);;
 500
 501
 502 let tgt_interval = prove(`interval_arith tgt (#1.541, #1.541)`,
 503    REWRITE_TAC[Tame_defs.tgt; interval_arith; REAL_LE_REFL]);;
 504
 505
 506
 507 let interval_table = Hashtbl.create 10;;
 508 let add_interval th = Hashtbl.add interval_table ((rand o rator o concl) th) th;;
 509
 510
 511
 512 let rec create_interval tm =
 513   if Hashtbl.mem interval_table tm then
 514     Hashtbl.find interval_table tm
 515   else
 516     if (is_ratconst tm) then
 517       const_interval tm
 518     else if (is_binop plus_op_real tm) then
 519       let lhs, rhs = dest_binop plus_op_real tm in
 520         interval_add (create_interval lhs) (create_interval rhs)
 521     else if (is_binop minus_op_real tm) then
 522       let lhs, rhs = dest_binop minus_op_real tm in
 523         interval_sub (create_interval lhs) (create_interval rhs)
 524     else if (is_binop mul_op_real tm) then
 525       let lhs, rhs = dest_binop mul_op_real tm in
 526         interval_mul (create_interval lhs) (create_interval rhs)
 527     else if (is_binop div_op_real tm) then
 528       let lhs, rhs = dest_binop div_op_real tm in
 529         interval_div (create_interval lhs) (create_interval rhs)
 530     else if (is_comb tm) then
 531       let ltm, rtm = dest_comb tm in
 532         if (ltm = inv_op_real) then
 533           interval_inv (create_interval rtm)
 534         else if (ltm = neg_op_real) then
 535           interval_neg (create_interval rtm)
 536         else failwith "create_interval: unknown unary operation"
 537     else
 538       failwith "create_interval: unexpected term";;
 539
 540
 541
 542 open Sphere;;
 543
 544
 545 add_interval pi_interval;;
 546 add_interval acs3_interval;;
 547 add_interval tgt_interval;;
 548 add_interval (REWRITE_RULE[GSYM sqrt8] (interval_sqrt (const_interval `&8`) 9));;
 549 add_interval (REWRITE_RULE[GSYM sol0] (create_interval `&3 * acs(&1 / &3) - pi`));;
 550 add_interval (create_interval `sol0 / pi`);;
 551
 552 let rho218 = new_definition `rho218 = rho(#2.18)`;;
 553
 554 let rho218_def = (REWRITE_CONV[rho218; rho; ly; interp; GSYM Tame_general.sol0_over_pi_EQ_const1] THENC
 555    REAL_RAT_REDUCE_CONV) `rho218`;;
 556
 557 let rho218_interval = REWRITE_RULE[SYM rho218_def] (create_interval(rand(concl rho218_def)));;
 558
 559 add_interval rho218_interval;;
 560
 561
 562 (*
 563 approx_interval (create_interval `rho (#2.18)`) 8;;
 564
 565 approx_interval (create_interval `&2 * (pi + sol0)`) 6;;
 566
 567 approx_interval (create_interval `sqrt8 - pi / sol0 + rho(#2.18) * &3`) 6;;
 568
 569 approx_interval (create_interval `pi * pi`) 6;;
 570 *)
 571
 572 end;;