formal_lp/old/formal_interval/test.hl

   1 needs "../formal_lp/formal_interval/verifier.hl";;
   2
   3
   4 (************)
   5 (* Examples *)
   6 (************)
   7
   8 (* f(x) = x * x - atn x; f(x) < 0 when 0.01 <= x <= 0.8 *)
   9 prove_ineq1d 5 `x * x - atn x` `#0.01` `#0.8`;;
  10
  11 (* f(x) = x * x * x - atn (x * x); f(x) < 0 when 0.01 <= x <= 0.8 *)
  12 prove_ineq1d 5 `x * x * x - atn (x * x)` `#0.01` `#0.8`;;
  13
  14 (* f(x) = atn x - x; f(x) < 0 when 0.001 <= x <= sqrt(2) *)
  15 prove_ineq1d 10 `atn x - x` `#0.001` `sqrt (&2) + pi`;;
  16
  17
  18 (****************************************)
  19 (* arclength #2.52 #2.52 (&2) >= #0.816 *)
  20 (* arclength (&2) #2.52 (&2) - arclength #2.52 #2.52 (&2) >= #0.073 *)
  21 (* arclength a b c = acs ((a * a + b * b - c * c) / (&2 * a * b)) *)
  22
  23 let arclength_tm = `acs ((a * a + b * b - c * c) / (&2 * a * b))`;;
  24 let tm1 = subst[`#2.52`, `a:real`; `#2.52`, `b:real`; `&2`, `c:real`] arclength_tm and
  25     tm2 = subst[`&2`, `a:real`; `#2.52`, `b:real`; `&2`, `c:real`] arclength_tm;;
  26
  27 let ineq1, ineq2 =
  28   let (-) = mk_binop sub_op_real in
  29     `#0.816` - tm1, `#0.073` - (tm2 - tm1);;
  30
  31 prove_ineq1d 6 ineq1 `&0` `&0`;;
  32 prove_ineq1d 6 ineq2 `&0` `&0`;;
  33
  34
  35 (********************************************)
  36 (* !x. &2 <= x /\ x <= #2.52                                           *)
  37 (*   ==> arclength (&2) x (&2) - arclength (#2.52) x (&2) >= #0.073`   *)
  38
  39 let ineq =
  40   let tm1 = subst[`&2`, `a:real`; `x:real`, `b:real`; `&2`, `c:real`] arclength_tm and
  41       tm2 = subst[`#2.52`, `a:real`; `x:real`, `b:real`; `&2`, `c:real`] arclength_tm and
  42       (-) = mk_binop sub_op_real in
  43     `#0.073` - (tm1 - tm2);;
  44
  45 prove_ineq1d 7 ineq `&2` `#2.52`;;
  46
  47
  48 (***********************)
  49 (* !x. &2 <= x /\ x <= #2.52  ==> &1 / &4 * --inv (sqrt (&1 - (x / &4) pow 2)) <= -- #0.28 *)
  50
  51 prove_ineq1d 4 `#0.28 + &1 / &4 * --inv (sqrt (&1 - (x * inv(&4)) * (x * inv(&4))))` `&2` `#2.52`;;
  52
  53 (* !x. &2 <= x /\ x <= #2.52 ==>
  54         inv (sqrt (&1 - ((&3969 / &625 + x * x - &4) / (&2 * #2.52 * x)) pow 2)) <= &2 *)
  55
  56 let ineq = (rand o concl o REWRITE_CONV[REAL_POW_2])
  57   `inv (sqrt (&1 - ((&3969 / &625 + x * x - &4) / (&2 * #2.52 * x)) pow 2)) - &2`;;
  58
  59 prove_ineq1d 3 ineq `&2` `#2.52`;;
  60
  61
  62
  63 (***********************************************)
  64 (* |atn x - x / (1 + 0.28 * x * x)| < 0.005 when x IN [0, 1] *)
  65 let abs_lemma = prove(`f - eps < &0 ==> -- eps - f < &0 ==> abs f < eps`, REAL_ARITH_TAC);;
  66
  67 let tm = `atn x - x / (&1 + #0.28 * x * x)`;;
  68 let eps = `#0.005`;;
  69 let hi_tm = mk_binop sub_op_real tm eps;;
  70 let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) tm;;
  71 let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 5 hi_tm `&0` `&1`);;
  72 let lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 6 lo_tm `&0` `&1`);;
  73 (GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;
  74
  75
  76
  77 (*************************************)
  78 (* A polynomial approximation of atn *)
  79 let tm = `x * (&1 - (x * x) * (&11184811 / &33554432 - (x * x) * (&13421773 / &67108864))) - atn x`;;
  80 let eps = `#0.0000001`;;
  81 let hi_tm = mk_binop sub_op_real tm eps;;
  82 let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) tm;;
  83 let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 8 hi_tm `&0` `&1 / &30`) and
  84     lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 8 lo_tm `&0` `&1 / &30`) and
  85     hi_th' = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 8 hi_tm `-- &1 / &30` `&0`) and
  86     lo_th' = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 8 lo_tm `-- &1 / &30` `&0`);;
  87
  88
  89 (GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;
  90
  91
  92
  93 (***********************************)
  94 (* A rational approximation of atn *)
  95
  96
  97 let atan1 = `let P00, P01, P02 = #4.26665376811246382, #3.291955407318148, #0.281939359003812 in
  98   let Q00, Q01 = #4.26665376811354027, #4.714173328637605 in
  99   let x1 = x * x in
 100   let x2 = x1 * x1 in
 101     x * (x1 * P01 + x2 * P02 + P00) / (x1 * Q01 + x2 + Q00) - atn x`;;
 102 let atan1_tm = (rand o concl o REPEATC let_CONV) atan1;;
 103
 104 let eps = `#0.00001`;;
 105 let hi_tm = mk_binop sub_op_real atan1_tm eps;;
 106 let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) atan1_tm;;
 107 let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 hi_tm `&0` `#0.2`);;
 108 let lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 lo_tm `&0` `#0.2`);;
 109 (GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;
 110
 111
 112 (*
 113 (* A higher precision approximation *)
 114 let eps = `#0.0000000001`;;
 115 let hi_tm = mk_binop sub_op_real atan1_tm eps;;
 116 let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) atan1_tm;;
 117 let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 15 hi_tm `&0` `#0.2`);;
 118 let lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 15 lo_tm `&0` `#0.2`);;
 119 (GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;
 120 *)
 121
 122
 123
 124 (****************************************************************)
 125 (* Test 1: f(x) = x * x - atn x; f(x) < 0 when 0.01 <= x <= 0.8 *)
 126 (****************************************************************)
 127
 128 llet eps = `#0.00001`;;
 129 let hi_tm = mk_binop sub_op_real atan1_tm eps;;
 130 let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) atan1_tm;;
 131 let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 hi_tm `&0` `#0.2`);;
 132 let lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 lo_tm `&0` `#0.2`);;
 133 (GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;
 134 let eps = `#0.00001`;;
 135 let hi_tm = mk_binop sub_op_real atan1_tm eps;;
 136 let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) atan1_tm;;
 137 let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 hi_tm `&0` `#0.2`);;
 138 let lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 lo_tm `&0` `#0.2`);;
 139 (GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;
 140 let eps = `#0.00001`;;
 141 let hi_tm = mk_binop sub_op_real atan1_tm eps;;
 142 let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) atan1_tm;;
 143 let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 hi_tm `&0` `#0.2`);;
 144 let lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 lo_tm `&0` `#0.2`);;
 145 (GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;
 146 et x_tm = mk_float 1 48;;
 147 let z_tm = mk_float 8 49;;
 148
 149 let domain_th = mk_center_domain 5 x_tm z_tm;;
 150 let _, y_tm, _, _ = dest_cell_domain (concl domain_th);;
 151
 152
 153 let domain1_th = mk_center_domain 5 x_tm y_tm;;
 154 let _, y1_tm, _, _ = dest_cell_domain (concl domain1_th);;
 155 let domain11_th = mk_center_domain 5 x_tm y1_tm;;
 156 let domain12_th = mk_center_domain 5 y1_tm y_tm;;
 157
 158 let domain2_th = mk_center_domain 5 y_tm z_tm;;
 159
 160
 161 let eval_f pp domain_th =
 162   let int_x = eval_taylor_x domain_th in
 163   let int_atn = eval_taylor_atn pp domain_th in
 164   let int_x2 = taylor_interval_mul pp int_x int_x in
 165     taylor_interval_sub pp int_x2 int_atn;;
 166
 167
 168 let int_th = eval_f 5 domain_th;;
 169 let int_th11 = eval_f 5 domain11_th;;
 170 let int_th12 = eval_f 5 domain12_th;;
 171 let int_th2 = eval_f 5 domain2_th;;
 172
 173 taylor_cell_pass 5 int_th;;
 174 let pass11 = taylor_cell_pass 5 int_th11;;
 175 let pass12 = taylor_cell_pass 5 int_th12;;
 176 let pass2 = taylor_cell_pass 5 int_th2;;
 177
 178 let pass1 = cell_pass_glue pass11 pass12;;
 179 let pass = cell_pass_glue pass1 pass2;;
 180
 181
 182
 183 (*************************************************************************)
 184 (* Test 2: f(x) = x * x * x - atn (x * x); f(x) < 0 when 0.1 <= x <= 0.8 *)
 185 (*************************************************************************)
 186
 187 let x_tm = mk_float 1 49;;
 188 let z_tm = mk_float 8 49;;
 189
 190 let domain_th = mk_center_domain 5 x_tm z_tm;;
 191 let _, y_tm, _, _ = dest_cell_domain (concl domain_th);;
 192
 193 let domain1_th = mk_center_domain 5 x_tm y_tm;;
 194 let _, y1_tm, _, _ = dest_cell_domain (concl domain1_th);;
 195
 196 let domain11_th = mk_center_domain 5 x_tm y1_tm;;
 197 let _, y11_tm, _, _ = dest_cell_domain (concl domain11_th);;
 198
 199 let domain111_th = mk_center_domain 5 x_tm y11_tm;;
 200 let domain112_th = mk_center_domain 5 y11_tm y1_tm;;
 201
 202 let domain12_th = mk_center_domain 5 y1_tm y_tm;;
 203 let domain2_th = mk_center_domain 5 y_tm z_tm;;
 204
 205 let eval_sqr pp domain_th =
 206   let x_th = eval_taylor_x domain_th in
 207     taylor_interval_mul pp x_th x_th;;
 208
 209
 210 let eval_f pp domain_th =
 211   let int_x = eval_taylor_x domain_th in
 212   let int_x2 = taylor_interval_mul pp int_x int_x in
 213   let int_x3 = taylor_interval_mul pp int_x int_x2 in
 214   let atn_x2 = taylor_interval_compose pp eval_taylor_atn eval_sqr domain_th in
 215     taylor_interval_sub pp int_x3 atn_x2;;
 216
 217 let int_th = eval_f 5 domain_th;;
 218 let int_th111 = eval_f 5 domain111_th;;
 219 let int_th112 = eval_f 5 domain112_th;;
 220 let int_th12 = eval_f 5 domain12_th;;
 221 let int_th2 = eval_f 5 domain2_th;;
 222
 223 taylor_cell_pass 5 int_th;;
 224 let pass111 = taylor_cell_pass 5 int_th111;;
 225 let pass112 = taylor_cell_pass 5 int_th112;;
 226 let pass12 = taylor_cell_pass 5 int_th12;;
 227 let pass2 = taylor_cell_pass 5 int_th2;;
 228
 229 let pass11 = cell_pass_glue pass111 pass112;;
 230 let pass1 = cell_pass_glue pass11 pass12;;
 231 let pass = cell_pass_glue pass1 pass2;;
 232
 233 let th0 = REWRITE_RULE[cell_pass; interval_arith] pass;;
 234 CONV_RULE (DEPTH_CONV FLOAT_TO_NUM_CONV) th0;;
 235
 236
 237
 238 (********************************************)
 239 (* Test 3:                                                             *)
 240 (* !x. &2 <= x /\ x <= #2.52                                           *)
 241 (*   ==> arclength (&2) x (&2) - arclength (#2.52) x (&2) >= #0.073`   *)
 242
 243 let eval_arc1 pp domain_th =
 244   let two = eval_taylor_const domain_th two_float and
 245       ( * ) = taylor_interval_mul pp and
 246       (+) = taylor_interval_add pp and
 247       (-) = taylor_interval_sub pp in
 248   let b = (fun pp domain_th ->
 249              let x = eval_taylor_x domain_th in
 250                two * two * x) in
 251   let a = (fun pp domain_th ->
 252              let x = eval_taylor_x domain_th in
 253              let rhs = taylor_interval_compose pp eval_taylor_inv b domain_th in
 254                (two * two + x * x - two * two) * rhs) in
 255     taylor_interval_compose pp eval_taylor_acs a domain_th;;
 256
 257
 258 let h_float = mk_float 252 48;;
 259 let eval_arc2 pp domain_th =
 260   let two = eval_taylor_const domain_th two_float and
 261       h = eval_taylor_const domain_th h_float and
 262       ( * ) = taylor_interval_mul pp and
 263       (+) = taylor_interval_add pp and
 264       (-) = taylor_interval_sub pp in
 265   let b = (fun pp domain_th ->
 266              let x = eval_taylor_x domain_th in
 267                two * h * x) in
 268   let a = (fun pp domain_th ->
 269              let x = eval_taylor_x domain_th in
 270              let rhs = taylor_interval_compose pp eval_taylor_inv b domain_th in
 271                (h * h + x * x - two * two) * rhs) in
 272     taylor_interval_compose pp eval_taylor_acs a domain_th;;
 273
 274
 275 let t_float = mk_float 73 47;;
 276 let eval_f pp domain_th =
 277   let t = eval_taylor_const domain_th t_float and
 278       (+) = taylor_interval_add pp and
 279       (-) = taylor_interval_sub pp and
 280       arc1 = eval_arc1 pp domain_th and
 281       arc2 = eval_arc2 pp domain_th in
 282     t + arc2 - arc1;;
 283
 284 let pp = 6;;
 285 let domain0 = mk_center_domain pp two_float h_float;;
 286
 287 let _, y1_tm, _, _ = dest_cell_domain (concl domain0);;
 288 let domain1 = mk_center_domain pp two_float y1_tm and
 289     domain2 = mk_center_domain pp y1_tm h_float;;
 290
 291 let _, y12_tm, _, _ = dest_cell_domain (concl domain2);;
 292 let domain21 = mk_center_domain pp y1_tm y12_tm and
 293     domain22 = mk_center_domain pp y12_tm h_float;;
 294
 295 let _, y122_tm, _, _ = dest_cell_domain (concl domain22);;
 296 let domain221 = mk_center_domain pp y12_tm y122_tm and
 297     domain222 = mk_center_domain pp y122_tm h_float;;
 298
 299 let int1 = eval_f pp domain1;;
 300 let int21 = eval_f pp domain21;;
 301 let int221 = eval_f pp domain221;;
 302 let int222 = eval_f pp domain222;;
 303
 304 let pass1 = taylor_cell_pass pp int1;;
 305 let pass21 = taylor_cell_pass pp int21;;
 306
 307 let pass221 = taylor_cell_pass pp int221;;
 308 let pass222 = taylor_cell_pass pp int222;;
 309
 310 let pass22 = cell_pass_glue pass221 pass222;;
 311 let pass2 = cell_pass_glue pass21 pass22;;
 312
 313 let pass = cell_pass_glue pass1 pass2;;
 314
 315 let th0 = REWRITE_RULE[cell_pass; interval_arith] pass;;
 316 CONV_RULE (DEPTH_CONV FLOAT_TO_NUM_CONV) th0;;
 317
 318
 319 let ineq_test pp =
 320   let domain0 = mk_center_domain pp two_float h_float in
 321   let _, y1_tm, _, _ = dest_cell_domain (concl domain0) in
 322   let domain1 = mk_center_domain pp two_float y1_tm and
 323       domain2 = mk_center_domain pp y1_tm h_float in
 324   let _, y12_tm, _, _ = dest_cell_domain (concl domain2) in
 325   let domain21 = mk_center_domain pp y1_tm y12_tm and
 326       domain22 = mk_center_domain pp y12_tm h_float in
 327   let _, y122_tm, _, _ = dest_cell_domain (concl domain22) in
 328   let domain221 = mk_center_domain pp y12_tm y122_tm and
 329       domain222 = mk_center_domain pp y122_tm h_float in
 330
 331   let int1 = eval_f pp domain1 and
 332       int21 = eval_f pp domain21 and
 333       int221 = eval_f pp domain221 and
 334       int222 = eval_f pp domain222 in
 335
 336   let pass1 = taylor_cell_pass pp int1 in
 337   let pass21 = taylor_cell_pass pp int21 in
 338   let pass221 = taylor_cell_pass pp int221 in
 339   let pass222 = taylor_cell_pass pp int222 in
 340   let pass22 = cell_pass_glue pass221 pass222 in
 341   let pass2 = cell_pass_glue pass21 pass22 in
 342     cell_pass_glue pass1 pass2;;
 343
 344
 345 ineq_test 5;;
 346 ineq_test 6;;
 347 ineq_test 7;;
 348 ineq_test 10;;
 349
 350 (* 10: 9.768 *)
 351 test 5 ineq_test 6;;
 352 (* 10: 11.012 *)
 353 test 5 ineq_test 7;;
 354 (* 10: 16.381 *)
 355 test 5 ineq_test 10;;
 356
 357
 358
 359
 360 (***************************)
 361
 362 let x_tm = mk_float 20 48;;
 363 let z_tm = mk_float 99 48;;
 364
 365 let domain_th = mk_center_domain 5 x_tm z_tm;;
 366 test 1000 (mk_center_domain 5 x_tm) z_tm;;
 367
 368 let int_x = eval_taylor_x domain_th;;
 369 let int_atn = eval_taylor_atn 5 domain_th;;
 370 let int_x2 = taylor_interval_mul 5 int_x int_x;;
 371 let int_x_atn = taylor_interval_mul 5 int_x int_atn;;
 372
 373 let int_th1 = taylor_interval_sub 5 int_x2 int_x;;
 374 let int_th2 = taylor_interval_sub 5 int_x_atn int_x;;
 375
 376 eval_taylor_f_bounds 5 int_th1;;
 377 eval_taylor_f_bounds 5 int_th2;;
 378
 379 let x_th =
 380   let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th) in
 381     (mk_bounded_on_int o ASSUME) (mk_interval x_var_real (mk_pair (x_tm, z_tm)));;
 382
 383
 384 let b_x2 = bounded_on_int_mul 5 x_th x_th;;
 385 bounded_on_int_sub 5 b_x2 x_th;;
 386
 387
 388
 389
 390
 391
 392 (*******************************************)
 393 (* arclength #2.52 #2.52 (&2) >= #0.816 *)
 394 (* arclength (&2) #2.52 (&2) - arclength #2.52 #2.52 (&2) >= #0.073 *)
 395 (* arclength a b c = acs ((a * a + b * b - c * c) / (&2 * a * b)) *)
 396
 397 let two_tm = mk_float 2 50;;
 398 let h_tm = mk_float 252 48;;
 399 let pp = 6;;
 400
 401
 402 let arclength_int pp a b c =
 403   let ( * ) = float_interval_mul pp and
 404       (+) = float_interval_add pp and
 405       (/) = float_interval_div pp and
 406       (-) = float_interval_sub pp and
 407       acs = float_interval_acs pp in
 408     acs ((a * a + b * b - c * c) / (two_interval * a * b));;
 409
 410 let result =
 411   let (!) = mk_const_interval and
 412       (-) = float_interval_sub pp in
 413     arclength_int pp !h_tm !h_tm !two_tm,
 414   arclength_int pp !two_tm !h_tm !two_tm - arclength_int pp !h_tm !h_tm !two_tm;;
 415
 416
 417 (CONV_RULE (DEPTH_CONV FLOAT_TO_NUM_CONV)) (fst result);;
 418 (CONV_RULE (DEPTH_CONV FLOAT_TO_NUM_CONV)) (snd result);;