formal_lp/old/hypermap/constants_approx.hl

   1 needs "arith/interval_arith.hl";;
   2 needs "misc/vars.hl";;
   3
   4 module Constant_approximations = struct
   5
   6 open Interval_arith;;
   7 open Misc_vars;;
   8
   9
  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);;
  12
  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`;;
  16
  17 let mk_decimal a b =
  18   let tm = mk_comb(mk_comb(`DECIMAL`, mk_numeral (Num.abs_num a)), mk_numeral b) in
  19     if (a </ Int 0) then
  20       mk_comb (neg_op_real, tm)
  21     else
  22       tm;;
  23
  24
  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;;
  37
  38
  39
  40 (************************************)
  41 (* Square root *)
  42
  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
  47      [
  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
  54          CONJ_TAC THENL
  55          [
  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
  64              ASM_REWRITE_TAC[]
  65          ]
  66      ]);;
  67
  68
  69
  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);;
  86
  87
  88
  89 (************************************)
  90 (* Arithmetic of intervals *)
  91
  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);;
  95
  96
  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);;
 100
 101
 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);;
 105
 106
 107
 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
 111      STRIP_TAC THENL
 112      [
 113        CONJ_TAC THEN MATCH_MP_TAC REAL_LE_INV2 THEN
 114          ASM_REWRITE_TAC[] THEN
 115          REPEAT (POP_ASSUM MP_TAC) THEN
 116          REAL_ARITH_TAC;
 117        ALL_TAC
 118      ] 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
 123      REAL_ARITH_TAC);;
 124
 125
 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]);;
 129
 130
 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]);;
 134
 135
 136
 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))`,
 139    REPEAT GEN_TAC THEN
 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
 143      [
 144        EXPAND_TAC "t" THEN
 145          REAL_ARITH_TAC;
 146        ALL_TAC
 147      ] THEN
 148
 149      DISJ_CASES_TAC (REAL_ARITH `&0 <= x \/ x <= &0`) THENL
 150      [
 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
 158          ASM_REWRITE_TAC[];
 159        ALL_TAC
 160      ] THEN
 161
 162      DISJ_CASES_TAC (REAL_ARITH `&0 <= b \/ b <= &0`) THENL
 163      [
 164        DISJ_CASES_TAC (REAL_ARITH `&0 <= y \/ y <= &0`) THENL
 165          [
 166            MATCH_MP_TAC REAL_LE_TRANS THEN
 167              EXISTS_TAC `&0` THEN
 168              CONJ_TAC THENL
 169              [
 170                ONCE_REWRITE_TAC[REAL_ARITH `&0 = &0 * y`] THEN
 171                  MATCH_MP_TAC REAL_LE_RMUL THEN
 172                  ASM_REWRITE_TAC[];
 173                ALL_TAC
 174              ] THEN
 175
 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
 183              ASM_REWRITE_TAC[];
 184            ALL_TAC
 185          ] THEN
 186
 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];
 193        ALL_TAC
 194      ] THEN
 195
 196      DISJ_CASES_TAC (REAL_ARITH `&0 <= c \/ c <= &0`) THENL
 197      [
 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];
 204        ALL_TAC
 205      ] THEN
 206
 207      DISJ_CASES_TAC (REAL_ARITH `&0 <= y \/ y <= &0`) THENL
 208      [
 209        MATCH_MP_TAC REAL_LE_TRANS THEN
 210          EXISTS_TAC `&0` THEN
 211          CONJ_TAC THENL
 212          [
 213            ONCE_REWRITE_TAC[REAL_ARITH `&0 = &0 * y`] THEN
 214              MATCH_MP_TAC REAL_LE_RMUL THEN
 215              ASM_REWRITE_TAC[];
 216            ALL_TAC
 217          ] THEN
 218
 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
 227          ASM_REWRITE_TAC[];
 228        ALL_TAC
 229      ] THEN
 230
 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]);;
 237
 238
 239
 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
 244      [
 245        MATCH_MP_TAC INTERVAL_MUL_lemma THEN
 246          ASM_REWRITE_TAC[];
 247        ALL_TAC
 248      ] THEN
 249
 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
 252      REAL_ARITH_TAC);;
 253
 254
 255
 256
 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
 265      REAL_ARITH_TAC);;
 266
 267
 268
 269 (**************************************)
 270 let const_interval tm = SPEC tm CONST_INTERVAL;;
 271
 272 let interval_neg th = MATCH_MP INTERVAL_NEG th;;
 273
 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;;
 277
 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;;
 281
 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;;
 285
 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)
 293     else
 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");;
 299
 300 let interval_div th1 th2 =
 301   let th2' = interval_inv th2 in
 302     ONCE_REWRITE_RULE[GSYM real_div] (interval_mul th1 th2');;
 303
 304
 305
 306 (*************************)
 307 let acs3_interval = REWRITE_RULE[GSYM interval_arith] (CONJ acs3_lo acs3_hi);;
 308
 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);;
 312
 313 let tgt_interval = prove(`interval_arith tgt (#1.541, #1.541)`,
 314    REWRITE_TAC[Tame_defs.tgt; interval_arith; REAL_LE_REFL]);;
 315
 316
 317 let interval_table = Hashtbl.create 10;;
 318 let add_interval th = Hashtbl.add interval_table ((rand o rator o concl) th) th;;
 319
 320
 321 let rec create_interval tm =
 322   if Hashtbl.mem interval_table tm then
 323     Hashtbl.find interval_table tm
 324   else
 325     if (is_ratconst tm) then
 326       const_interval tm
 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"
 346     else
 347       failwith "create_interval: unexpected term";;
 348
 349
 350 open Sphere;;
 351
 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`);;
 358
 359 let rho218 = new_definition `rho218 = rho(#2.18)`;;
 360
 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`;;
 363
 364 let rho218_interval = REWRITE_RULE[SYM rho218_def] (create_interval(rand(concl rho218_def)));;
 365 add_interval rho218_interval;;
 366
 367
 368 end;;