formal_lp/old/arith/arith_hash_rat.hl

   1 needs "../formal_lp/arith/arith_hash_int.hl";;
   2
   3
   4 let MY_PROVE_HYP hyp th = EQ_MP (DEDUCT_ANTISYM_RULE hyp th) hyp;;
   5 let REFL_CONV tm = EQT_INTRO (REFL (rand tm));;
   6
   7
   8 (* Constants *)
   9
  10 let div_op_real = `(/):real->real->real` and
  11     inv_op_real = `inv:real->real`;;
  12
  13
  14 let REPLACE_NUMERALS_CONV = DEPTH_CONV from_numeral;;
  15 let REPLACE_NUMERALS = CONV_RULE REPLACE_NUMERALS_CONV;;
  16 let REPLACE_NUMS = CONV_RULE (DEPTH_CONV to_numeral);;
  17
  18
  19 (******************************)
  20
  21 (* INT_RAT_CONV *)
  22
  23 let int_to_rat = (REPLACE_NUMERALS o prove) (`&n = &n / &1`,
  24                        CONV_TAC (DEPTH_CONV to_numeral) THEN REAL_ARITH_TAC) and
  25     neg_int_to_rat = (REPLACE_NUMERALS o prove) (`-- &n = -- &n / &1`,
  26                            CONV_TAC (DEPTH_CONV to_numeral) THEN REAL_ARITH_TAC);;
  27
  28
  29 let my_real_int_rat_conv tm =
  30   if (is_binop div_op_real tm) then
  31     REFL tm
  32   else
  33     let ltm, rtm = dest_comb tm in
  34       if (ltm = neg_op_real) then
  35         let amp_tm, n_tm = dest_comb rtm in
  36           if (amp_tm = amp_op_real) then
  37             INST[n_tm, n_var_num] neg_int_to_rat
  38           else
  39             failwith "my_real_int_rat_conv: --(&n) expected"
  40       else
  41         if (ltm = amp_op_real) then
  42           INST[rtm, n_var_num] int_to_rat
  43         else
  44           failwith "my_real_int_rat_conv: &n expected";;
  45
  46
  47 (*
  48 let tm = `-- &343`;;
  49 (* 0.180 *)
  50 test 10000 REAL_INT_RAT_CONV tm;;
  51 (* 0.044 *)
  52 test 10000 my_real_int_rat_conv tm;;
  53 *)
  54
  55
  56
  57 (*********************************)
  58
  59 (* ADD *)
  60
  61 let add_th = (UNDISCH_ALL o REPLACE_NUMERALS o prove) (`0 < n1==> 0 < n2 ==> 0 < n3 ==>
  62                      (x1 * &n2 + x2 * &n1) * &n3 = x3 * &n1 * &n2 ==>
  63                      x1 / &n1 + x2 / &n2 = x3 / &n3`,
  64   REWRITE_TAC[GSYM REAL_OF_NUM_LT] THEN
  65     MAP_EVERY ABBREV_TAC [`y1 = &n1`; `y2 = &n2`; `y3 = &n3`] THEN
  66     REPEAT DISCH_TAC THEN
  67     MP_TAC RAT_LEMMA2 THEN
  68     ASM_REWRITE_TAC[] THEN
  69     DISCH_THEN SUBST1_TAC THEN
  70     REWRITE_TAC[GSYM REAL_INV_MUL; GSYM real_div] THEN
  71     SUBGOAL_THEN `&0 < y1 * y2 /\ &0 < y3` MP_TAC THENL
  72      [ASM_REWRITE_TAC[] THEN MATCH_MP_TAC REAL_LT_MUL THEN
  73       ASM_REWRITE_TAC[];
  74       DISCH_THEN(fun th -> ASM_REWRITE_TAC[MATCH_MP RAT_LEMMA5 th])]);;
  75
  76 let x1_var = `x1:real` and x2_var = `x2:real` and x3_var = `x3:real` and
  77     y1_var = `y1:real` and y2_var = `y2:real` and
  78     n1_var = `n1:num` and n2_var = `n2:num` and n3_var = `n3:num`;;
  79
  80
  81 let raw_real_rat_add_conv tm =
  82   let r1, r2 = dest_binop plus_op_real tm in
  83   let x1_tm, y1_tm = dest_binop div_op_real r1 and
  84       x2_tm, y2_tm = dest_binop div_op_real r2 in
  85   let x1n = my_dest_realintconst x1_tm and
  86       y1n = my_dest_realintconst y1_tm and
  87       x2n = my_dest_realintconst x2_tm and
  88       y2n = my_dest_realintconst y2_tm in
  89   let x3n = x1n */ y2n +/ x2n */ y1n and
  90       y3n = y1n */ y2n in
  91   let d = gcd_num x3n y3n in
  92     let x3n' = quo_num x3n d and y3n' = quo_num y3n d in
  93     let x3n'',y3n'' = if y3n' >/ Int 0 then x3n',y3n'
  94                       else minus_num x3n',minus_num y3n' in
  95     let x3_tm = my_mk_realintconst x3n'' and
  96         y3_tm = my_mk_realintconst y3n'' in
  97
  98     let n1_tm = rand y1_tm and n2_tm = rand y2_tm and n3_tm = rand y3_tm in
  99     let n1_pos = EQT_ELIM(num_gt0 n1_tm) and
 100         n2_pos = EQT_ELIM(num_gt0 n2_tm) and
 101         n3_pos = EQT_ELIM(num_gt0 n3_tm) in
 102
 103     let th0 = INST [x1_tm, x1_var; n1_tm, n1_var; x2_tm, x2_var;
 104                     n2_tm, n2_var; x3_tm, x3_var; n3_tm, n3_var] add_th in
 105     let th1 = MY_PROVE_HYP n1_pos (MY_PROVE_HYP n2_pos (MY_PROVE_HYP n3_pos th0)) in
 106     let tm2, tm3 = (dest_eq o hd o hyp) th1 in
 107     let th2 = (LAND_CONV (BINOP_CONV my_real_int_mul_conv THENC
 108                           my_real_int_add_conv) THENC my_real_int_mul_conv) tm2 and
 109         th3 = (RAND_CONV my_real_int_mul_conv THENC my_real_int_mul_conv) tm3 in
 110       MY_PROVE_HYP (TRANS th2 (SYM th3)) th1;;
 111
 112
 113 let my_real_rat_add_conv tm =
 114   (BINOP_CONV my_real_int_rat_conv THENC raw_real_rat_add_conv) tm;;
 115
 116
 117 (*
 118 let tm = `-- &235346 / &146424 + -- &44635 / &3463462`;;
 119 let tm' = replace_numerals tm;;
 120
 121 (* 2.868 *)
 122 test 100 REAL_RAT_ADD_CONV tm;;
 123 (* 0.216 *)
 124 test 100 my_real_rat_add_conv tm';;
 125 *)
 126
 127
 128 (*************************************)
 129
 130 (* MUL *)
 131
 132 let mul_nocancel = prove(`(x1 / y1) * (x2 / y2) = (x1 * x2) / (y1 * y2)`,
 133                          REWRITE_TAC[real_div; REAL_INV_MUL; REAL_MUL_AC]);;
 134 let mul_cancel = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o REPLACE_NUMERALS o prove)
 135                         (`0 < n1 /\ 0 < n2 /\
 136                          (&n1 * u1 = x1) /\ (&n2 * u2 = x2) /\
 137                          (&n2 * v1 = y1) /\ (&n1 * v2 = y2)
 138                        ==> ((x1 / y1) * (x2 / y2) = (u1 * u2) / (v1 * v2))`,
 139    REWRITE_TAC[ARITH_RULE `0 < n <=> ~(n = 0)`] THEN
 140      REWRITE_TAC[GSYM REAL_OF_NUM_EQ] THEN
 141      MAP_EVERY ABBREV_TAC [`d1 = &n1`; `d2 = &n2`] THEN
 142      DISCH_THEN(CONJUNCTS_THEN2 ASSUME_TAC MP_TAC) THEN
 143      DISCH_THEN(CONJUNCTS_THEN2 ASSUME_TAC MP_TAC) THEN
 144      DISCH_THEN(REPEAT_TCL CONJUNCTS_THEN (SUBST1_TAC o SYM)) THEN
 145      ASM_REWRITE_TAC[real_div; REAL_INV_MUL] THEN
 146      ONCE_REWRITE_TAC[AC REAL_MUL_AC
 147         `((d1 * u1) * (id2 * iv1)) * ((d2 * u2) * id1 * iv2) =
 148          (u1 * u2) * (iv1 * iv2) * (id2 * d2) * (id1 * d1)`] THEN
 149      ASM_SIMP_TAC[REAL_MUL_LINV; REAL_MUL_RID]);;
 150
 151 let u1_var = `u1:real` and v1_var = `v1:real` and
 152     u2_var = `u2:real` and v2_var = `v2:real`;;
 153
 154
 155
 156 let raw_real_rat_mul_conv tm =
 157   let r1,r2 = dest_binop mul_op_real tm in
 158   let x1',y1' = dest_binop div_op_real r1 and
 159       x2',y2' = dest_binop div_op_real r2 in
 160
 161   let x1n = my_dest_realintconst x1' and
 162       y1n = my_dest_realintconst y1' and
 163       x2n = my_dest_realintconst x2' and
 164       y2n = my_dest_realintconst y2' in
 165
 166   let d1n = gcd_num x1n y2n and
 167       d2n = gcd_num x2n y1n in
 168     if d1n = num_1 & d2n = num_1 then
 169       let th0 = INST [x1',x1_var; y1',y1_var; x2',x2_var; y2',y2_var] mul_nocancel in
 170       let th1 = BINOP_CONV my_real_int_mul_conv (rand(concl th0)) in
 171         TRANS th0 th1
 172     else
 173       let u1n = quo_num x1n d1n and
 174           u2n = quo_num x2n d2n and
 175           v1n = quo_num y1n d2n and
 176           v2n = quo_num y2n d1n in
 177       let u1' = my_mk_realintconst u1n and
 178           u2' = my_mk_realintconst u2n and
 179           v1' = my_mk_realintconst v1n and
 180           v2' = my_mk_realintconst v2n and
 181           n1' = mk_num d1n and
 182           n2' = mk_num d2n in
 183
 184       let th0 = INST [x1',x1_var; y1',y1_var; x2',x2_var; y2',y2_var;
 185                       u1',u1_var; v1',v1_var; u2',u2_var; v2',v2_var;
 186                       n1',n1_var; n2',n2_var] mul_cancel in
 187       let n1_pos = EQT_ELIM(num_gt0 n1') and
 188           n2_pos = EQT_ELIM(num_gt0 n2') in
 189       let th1 = MY_PROVE_HYP n1_pos (MY_PROVE_HYP n2_pos th0) in
 190       let hyp_ths = map (EQT_ELIM o (LAND_CONV my_real_int_mul_conv THENC REFL_CONV)) (hyp th1) in
 191       let th2 = itlist MY_PROVE_HYP hyp_ths th1 in
 192       let th3 = BINOP_CONV my_real_int_mul_conv (rand(concl th2)) in
 193         TRANS th2 th3;;
 194
 195
 196 let my_real_rat_mul_conv tm =
 197   (BINOP_CONV my_real_int_rat_conv THENC raw_real_rat_mul_conv) tm;;
 198
 199
 200 (*
 201 let tm = `-- &656 / &4567 * &4566 / &3666`;;
 202 let tm' = replace_numerals tm;;
 203
 204 my_real_rat_mul_conv tm';;
 205
 206 (* 6.384 *)
 207 test 1000 REAL_RAT_MUL_CONV tm;;
 208 (* 0.640 *)
 209 test 1000 my_real_rat_mul_conv tm';;
 210 *)
 211
 212
 213
 214 (***********************************)
 215
 216 (* DIV *)
 217
 218
 219 let div_th = prove(`(x1 / y1) / (x2 / y2) = (x1 / y1) * (y2 / x2)`,
 220    REWRITE_TAC[real_div; REAL_INV_MUL; REAL_INV_INV; REAL_MUL_AC]);;
 221
 222
 223 let my_real_rat_div_conv tm =
 224   let th0 = BINOP_CONV my_real_int_rat_conv tm in
 225   let r1, r2 = dest_binop div_op_real (rand(concl th0)) in
 226   let x1, y1 = dest_binop div_op_real r1 and
 227       x2, y2 = dest_binop div_op_real r2 in
 228   let th1 = INST[x1,x1_var; y1,y1_var; x2,x2_var; y2,y2_var] div_th in
 229   let th2 = raw_real_rat_mul_conv (rand(concl th1)) in
 230     TRANS th0 (TRANS th1 th2);;
 231
 232
 233
 234
 235 (***********************************************)
 236
 237 (* Polynomial functions *)
 238
 239
 240 let poly_f = new_definition `poly_f cs x = ITLIST (\c s. c + x * s) cs (&0)`;;
 241
 242 (* Even function *)
 243 let poly_f_even = new_definition `poly_f_even cs x = ITLIST (\c s. c + (x * x) * s) cs (&0)`;;
 244 (* Odd function *)
 245 let poly_f_odd = new_definition `poly_f_odd cs x = x * poly_f_even cs x`;;
 246 let poly_f_odd' = SPECL[`t:(real)list`; `x:real`] poly_f_odd;;
 247
 248 let POLY_F_EMPTY = (REPLACE_NUMERALS o prove) (`poly_f [] x = &0`, REWRITE_TAC[poly_f; ITLIST]) and
 249     POLY_F_CONS = prove(`poly_f (CONS h t) x = h + x * poly_f t x`, REWRITE_TAC[poly_f; ITLIST]);;
 250
 251 let POLY_F_EVEN_EMPTY = (REPLACE_NUMERALS o prove) (`poly_f_even [] x = &0`, REWRITE_TAC[poly_f_even; ITLIST]) and
 252     POLY_F_EVEN_CONS = prove(`poly_f_even (CONS h t) x = h + (x * x) * poly_f_even t x`, REWRITE_TAC[poly_f_even; ITLIST]);;
 253
 254 let POLY_F_ODD_EMPTY = (REPLACE_NUMERALS o prove) (`poly_f_odd [] x = &0`, REWRITE_TAC[poly_f_odd; poly_f_even; ITLIST; REAL_MUL_RZERO]);;
 255
 256
 257 let x_var = `x:real` and
 258     h_var = `h:real` and
 259     t_var = `t:(real)list`;;
 260
 261
 262 (*************************)
 263
 264 (* poly_f_conv, poly_f_even_conv, poly_f_odd_conv *)
 265
 266
 267
 268 let poly_f_conv tm =
 269   let ltm, x_tm = dest_comb tm in
 270   let list_tm = rand ltm in
 271   let inst_t = INST[x_tm, x_var] in
 272   let poly_f_cons, poly_f_empty = inst_t POLY_F_CONS, inst_t POLY_F_EMPTY in
 273
 274   let rec poly_f_conv_raw list_tm =
 275     if (is_comb list_tm) then
 276       let h_tm, t_tm = dest_comb list_tm in
 277       let th0 = INST[rand h_tm, h_var; t_tm, t_var] poly_f_cons in
 278       let h_plus, rtm = dest_comb(rand(concl th0)) in
 279       let x_times = rator rtm in
 280       let th1 = poly_f_conv_raw t_tm in
 281       let th2 = TRANS th0 (AP_TERM h_plus (AP_TERM x_times th1)) in
 282       let th3 = (RAND_CONV my_real_rat_mul_conv THENC my_real_rat_add_conv) (rand(concl th2)) in
 283         TRANS th2 th3
 284     else
 285       poly_f_empty in
 286
 287     poly_f_conv_raw list_tm;;
 288
 289
 290
 291 let poly_f_even_conv tm =
 292   let ltm, x_tm = dest_comb tm in
 293   let list_tm = rand ltm in
 294   let inst_t = INST[x_tm, x_var] in
 295   let poly_f_even_cons = inst_t POLY_F_EVEN_CONS and
 296       poly_f_even_empty = inst_t POLY_F_EVEN_EMPTY in
 297   let x2_times = AP_TERM mul_op_real (my_real_rat_mul_conv (mk_binop mul_op_real x_tm x_tm)) in
 298
 299   let rec poly_f_even_conv_raw list_tm =
 300     if (is_comb list_tm) then
 301       let h_tm, t_tm = dest_comb list_tm in
 302       let th0 = INST[rand h_tm, h_var; t_tm, t_var] poly_f_even_cons in
 303       let h_plus, rtm = dest_comb(rand(concl th0)) in
 304       let th1 = poly_f_even_conv_raw t_tm in
 305       let th2 = TRANS th0 (AP_TERM h_plus (MK_COMB(x2_times, th1))) in
 306       let th3 = (RAND_CONV my_real_rat_mul_conv THENC my_real_rat_add_conv) (rand(concl th2)) in
 307         TRANS th2 th3
 308     else
 309       poly_f_even_empty in
 310
 311     poly_f_even_conv_raw list_tm;;
 312
 313
 314
 315 let poly_f_odd_conv tm =
 316   let ltm, x_tm = dest_comb tm in
 317   let list_tm = rand ltm in
 318   let th0 = INST[list_tm, t_var; x_tm, x_var] poly_f_odd' in
 319   let ltm, rtm = dest_comb(rand(concl th0)) in
 320   let th1 = AP_TERM ltm (poly_f_even_conv rtm) in
 321   let th2 = my_real_rat_mul_conv (rand(concl th1)) in
 322     TRANS th0 (TRANS th1 th2);;
 323
 324
 325
 326 (*
 327 let tm = `poly_f [-- &1436346 / &436346; -- &244664 / &4654235; &3 / &43545] (-- &144545 / &2345345)`;;
 328 let tm2 = (rand o concl) (REWRITE_CONV[poly_f; ITLIST] tm);;
 329 let tm' = replace_numerals tm;;
 330
 331 poly_f_conv tm';;
 332
 333 (* 3.280 *)
 334 test 10 (REWRITE_CONV [poly_f; ITLIST] THENC REAL_RAT_REDUCE_CONV) tm;;
 335 (* 3.188 *)
 336 test 10 REAL_RAT_REDUCE_CONV tm2;;
 337 (* 0.184 *)
 338 test 10 poly_f_conv tm';;
 339 *)
 340
 341
 342 (********************************)
 343
 344
 345 (*
 346 let REVERSE_TABLE  = define `(REVERSE_TABLE (f:num->A) 0 = []) /\
 347    (REVERSE_TABLE f (SUC i) = CONS (f i)  ( REVERSE_TABLE f i))`;;
 348
 349 let TABLE = new_definition `!(f:num->A) k. TABLE f k = REVERSE (REVERSE_TABLE f k)`;;
 350
 351
 352 let rec reverse_table_conv tm =
 353   let ltm, i_tm = dest_comb tm in
 354     if (i_tm = `0`) then
 355       ONCE_REWRITE_CONV[REVERSE_TABLE] tm
 356     else
 357       let i_suc = num_CONV i_tm in
 358       let th1 = ONCE_REWRITE_RULE[REVERSE_TABLE] (AP_TERM ltm i_suc) in
 359       let ltm, rtm = dest_comb (rand(concl th1)) in
 360       let th2 = reverse_table_conv rtm in
 361         TRANS th1 (AP_TERM ltm th2);;
 362
 363
 364 let cos_taylor n =
 365   let tm = mk_comb(`TABLE (\k. (if (EVEN k) then &1 else --(&1)) / &(FACT(2 * k)))`,
 366                    mk_small_numeral n) in
 367     (rand o concl) ((REWRITE_CONV[TABLE] THENC
 368        ONCE_DEPTH_CONV reverse_table_conv THENC
 369        REWRITE_CONV[REVERSE; APPEND] THENC NUM_REDUCE_CONV) tm);;
 370
 371
 372 let cos_poly n x =
 373   let tm = mk_comb(mk_comb(`poly_f_even`, cos_taylor n), x) in
 374   let tm' = replace_numerals tm in
 375     tm';;
 376
 377 let x = rand(concl(REWRITE_CONV[DECIMAL] `#1.230959417`));;
 378
 379 let tm = cos_poly 6 x;;
 380
 381 let tm' = (rand(concl(REPLACE_NUMS (REWRITE_CONV[poly_f_even; ITLIST] tm))));;
 382 (* 5.860 *)
 383 test 1 REAL_RAT_REDUCE_CONV tm';;
 384
 385 (* 0.200 *)
 386 test 1 poly_f_even_conv tm;;
 387 let result = poly_f_even_conv tm;;
 388 REPLACE_NUMS result;;
 389 *)