needs "../formal_lp/arith/arith_hash_int.hl";;


let MY_PROVE_HYP hyp th = EQ_MP (DEDUCT_ANTISYM_RULE hyp th) hyp;;
let REFL_CONV tm = EQT_INTRO (REFL (rand tm));;


(* Constants *)

let div_op_real = `(/):real->real->real` and
    inv_op_real = `inv:real->real`;;


let REPLACE_NUMERALS_CONV = DEPTH_CONV from_numeral;;
let REPLACE_NUMERALS = CONV_RULE REPLACE_NUMERALS_CONV;;
let REPLACE_NUMS = CONV_RULE (DEPTH_CONV to_numeral);;


(******************************)

(* INT_RAT_CONV *)

let int_to_rat = (REPLACE_NUMERALS o prove) (`&n = &n / &1`, 
		       CONV_TAC (DEPTH_CONV to_numeral) THEN REAL_ARITH_TAC) and
    neg_int_to_rat = (REPLACE_NUMERALS o prove) (`-- &n = -- &n / &1`,
			   CONV_TAC (DEPTH_CONV to_numeral) THEN REAL_ARITH_TAC);;


let my_real_int_rat_conv tm =
  if (is_binop div_op_real tm) then
    REFL tm
  else
    let ltm, rtm = dest_comb tm in
      if (ltm = neg_op_real) then
	let amp_tm, n_tm = dest_comb rtm in
	  if (amp_tm = amp_op_real) then
	    INST[n_tm, n_var_num] neg_int_to_rat
	  else
	    failwith "my_real_int_rat_conv: --(&n) expected"
      else
	if (ltm = amp_op_real) then
	  INST[rtm, n_var_num] int_to_rat
	else
	  failwith "my_real_int_rat_conv: &n expected";;


(*
let tm = `-- &343`;;
(* 0.180 *)
test 10000 REAL_INT_RAT_CONV tm;;
(* 0.044 *)
test 10000 my_real_int_rat_conv tm;;
*)



(*********************************)

(* ADD *)

let add_th = (UNDISCH_ALL o REPLACE_NUMERALS o prove) (`0 < n1==> 0 < n2 ==> 0 < n3 ==>
		     (x1 * &n2 + x2 * &n1) * &n3 = x3 * &n1 * &n2 ==>
		     x1 / &n1 + x2 / &n2 = x3 / &n3`,
  REWRITE_TAC[GSYM REAL_OF_NUM_LT] THEN
    MAP_EVERY ABBREV_TAC [`y1 = &n1`; `y2 = &n2`; `y3 = &n3`] THEN
    REPEAT DISCH_TAC THEN
    MP_TAC RAT_LEMMA2 THEN
    ASM_REWRITE_TAC[] THEN
    DISCH_THEN SUBST1_TAC THEN
    REWRITE_TAC[GSYM REAL_INV_MUL; GSYM real_div] THEN
    SUBGOAL_THEN `&0 < y1 * y2 /\ &0 < y3` MP_TAC THENL
     [ASM_REWRITE_TAC[] THEN MATCH_MP_TAC REAL_LT_MUL THEN
      ASM_REWRITE_TAC[];
      DISCH_THEN(fun th -> ASM_REWRITE_TAC[MATCH_MP RAT_LEMMA5 th])]);;

let x1_var = `x1:real` and x2_var = `x2:real` and x3_var = `x3:real` and
    y1_var = `y1:real` and y2_var = `y2:real` and
    n1_var = `n1:num` and n2_var = `n2:num` and n3_var = `n3:num`;;


let raw_real_rat_add_conv tm =
  let r1, r2 = dest_binop plus_op_real tm in
  let x1_tm, y1_tm = dest_binop div_op_real r1 and
      x2_tm, y2_tm = dest_binop div_op_real r2 in
  let x1n = my_dest_realintconst x1_tm and 
      y1n = my_dest_realintconst y1_tm and
      x2n = my_dest_realintconst x2_tm and 
      y2n = my_dest_realintconst y2_tm in
  let x3n = x1n */ y2n +/ x2n */ y1n and
      y3n = y1n */ y2n in
  let d = gcd_num x3n y3n in
    let x3n' = quo_num x3n d and y3n' = quo_num y3n d in
    let x3n'',y3n'' = if y3n' >/ Int 0 then x3n',y3n'
                      else minus_num x3n',minus_num y3n' in
    let x3_tm = my_mk_realintconst x3n'' and 
	y3_tm = my_mk_realintconst y3n'' in

    let n1_tm = rand y1_tm and n2_tm = rand y2_tm and n3_tm = rand y3_tm in
    let n1_pos = EQT_ELIM(num_gt0 n1_tm) and
	n2_pos = EQT_ELIM(num_gt0 n2_tm) and
	n3_pos = EQT_ELIM(num_gt0 n3_tm) in

    let th0 = INST [x1_tm, x1_var; n1_tm, n1_var; x2_tm, x2_var; 
		    n2_tm, n2_var; x3_tm, x3_var; n3_tm, n3_var] add_th in
    let th1 = MY_PROVE_HYP n1_pos (MY_PROVE_HYP n2_pos (MY_PROVE_HYP n3_pos th0)) in
    let tm2, tm3 = (dest_eq o hd o hyp) th1 in
    let th2 = (LAND_CONV (BINOP_CONV my_real_int_mul_conv THENC
                          my_real_int_add_conv) THENC my_real_int_mul_conv) tm2 and
	th3 = (RAND_CONV my_real_int_mul_conv THENC my_real_int_mul_conv) tm3 in
      MY_PROVE_HYP (TRANS th2 (SYM th3)) th1;;


let my_real_rat_add_conv tm =
  (BINOP_CONV my_real_int_rat_conv THENC raw_real_rat_add_conv) tm;;
  

(*
let tm = `-- &235346 / &146424 + -- &44635 / &3463462`;;
let tm' = replace_numerals tm;;

(* 2.868 *)
test 100 REAL_RAT_ADD_CONV tm;;
(* 0.216 *)
test 100 my_real_rat_add_conv tm';;
*)


(*************************************)

(* MUL *)

let mul_nocancel = 
prove(`(x1 / y1) * (x2 / y2) = (x1 * x2) / (y1 * y2)`,
REWRITE_TAC[real_div; REAL_INV_MUL; REAL_MUL_AC]);;
let mul_cancel = (UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o REPLACE_NUMERALS o prove) (`0 < n1 /\ 0 < n2 /\ (&n1 * u1 = x1) /\ (&n2 * u2 = x2) /\ (&n2 * v1 = y1) /\ (&n1 * v2 = y2) ==> ((x1 / y1) * (x2 / y2) = (u1 * u2) / (v1 * v2))`, REWRITE_TAC[ARITH_RULE `0 < n <=> ~(n = 0)`] THEN REWRITE_TAC[GSYM REAL_OF_NUM_EQ] THEN MAP_EVERY ABBREV_TAC [`d1 = &n1`; `d2 = &n2`] THEN DISCH_THEN(CONJUNCTS_THEN2 ASSUME_TAC MP_TAC) THEN DISCH_THEN(CONJUNCTS_THEN2 ASSUME_TAC MP_TAC) THEN DISCH_THEN(REPEAT_TCL CONJUNCTS_THEN (SUBST1_TAC o SYM)) THEN ASM_REWRITE_TAC[real_div; REAL_INV_MUL] THEN ONCE_REWRITE_TAC[AC REAL_MUL_AC `((d1 * u1) * (id2 * iv1)) * ((d2 * u2) * id1 * iv2) = (u1 * u2) * (iv1 * iv2) * (id2 * d2) * (id1 * d1)`] THEN ASM_SIMP_TAC[REAL_MUL_LINV; REAL_MUL_RID]);; let u1_var = `u1:real` and v1_var = `v1:real` and u2_var = `u2:real` and v2_var = `v2:real`;; let raw_real_rat_mul_conv tm = let r1,r2 = dest_binop mul_op_real tm in let x1',y1' = dest_binop div_op_real r1 and x2',y2' = dest_binop div_op_real r2 in let x1n = my_dest_realintconst x1' and y1n = my_dest_realintconst y1' and x2n = my_dest_realintconst x2' and y2n = my_dest_realintconst y2' in let d1n = gcd_num x1n y2n and d2n = gcd_num x2n y1n in if d1n = num_1 & d2n = num_1 then let th0 = INST [x1',x1_var; y1',y1_var; x2',x2_var; y2',y2_var] mul_nocancel in let th1 = BINOP_CONV my_real_int_mul_conv (rand(concl th0)) in TRANS th0 th1 else let u1n = quo_num x1n d1n and u2n = quo_num x2n d2n and v1n = quo_num y1n d2n and v2n = quo_num y2n d1n in let u1' = my_mk_realintconst u1n and u2' = my_mk_realintconst u2n and v1' = my_mk_realintconst v1n and v2' = my_mk_realintconst v2n and n1' = mk_num d1n and n2' = mk_num d2n in let th0 = INST [x1',x1_var; y1',y1_var; x2',x2_var; y2',y2_var; u1',u1_var; v1',v1_var; u2',u2_var; v2',v2_var; n1',n1_var; n2',n2_var] mul_cancel in let n1_pos = EQT_ELIM(num_gt0 n1') and n2_pos = EQT_ELIM(num_gt0 n2') in let th1 = MY_PROVE_HYP n1_pos (MY_PROVE_HYP n2_pos th0) in let hyp_ths = map (EQT_ELIM o (LAND_CONV my_real_int_mul_conv THENC REFL_CONV)) (hyp th1) in let th2 = itlist MY_PROVE_HYP hyp_ths th1 in let th3 = BINOP_CONV my_real_int_mul_conv (rand(concl th2)) in TRANS th2 th3;; let my_real_rat_mul_conv tm = (BINOP_CONV my_real_int_rat_conv THENC raw_real_rat_mul_conv) tm;; (* let tm = `-- &656 / &4567 * &4566 / &3666`;; let tm' = replace_numerals tm;; my_real_rat_mul_conv tm';; (* 6.384 *) test 1000 REAL_RAT_MUL_CONV tm;; (* 0.640 *) test 1000 my_real_rat_mul_conv tm';; *) (***********************************) (* DIV *)
let div_th = 
prove(`(x1 / y1) / (x2 / y2) = (x1 / y1) * (y2 / x2)`,
let my_real_rat_div_conv tm = let th0 = BINOP_CONV my_real_int_rat_conv tm in let r1, r2 = dest_binop div_op_real (rand(concl th0)) in let x1, y1 = dest_binop div_op_real r1 and x2, y2 = dest_binop div_op_real r2 in let th1 = INST[x1,x1_var; y1,y1_var; x2,x2_var; y2,y2_var] div_th in let th2 = raw_real_rat_mul_conv (rand(concl th1)) in TRANS th0 (TRANS th1 th2);; (***********************************************) (* Polynomial functions *)
let poly_f = new_definition `poly_f cs x = ITLIST (\c s. c + x * s) cs (&0)`;;
(* Even function *)
let poly_f_even = new_definition `poly_f_even cs x = ITLIST (\c s. c + (x * x) * s) cs (&0)`;;
(* Odd function *)
let poly_f_odd = new_definition `poly_f_odd cs x = x * poly_f_even cs x`;;
let poly_f_odd' = SPECL[`t:(real)list`; `x:real`] poly_f_odd;; let POLY_F_EMPTY = (REPLACE_NUMERALS o prove) (`poly_f [] x = &0`, REWRITE_TAC[poly_f; ITLIST]) and POLY_F_CONS = prove(`poly_f (CONS h t) x = h + x * poly_f t x`, REWRITE_TAC[poly_f; ITLIST]);; let POLY_F_EVEN_EMPTY = (REPLACE_NUMERALS o prove) (`poly_f_even [] x = &0`, REWRITE_TAC[poly_f_even; ITLIST]) and 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]);; 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]);; let x_var = `x:real` and h_var = `h:real` and t_var = `t:(real)list`;; (*************************) (* poly_f_conv, poly_f_even_conv, poly_f_odd_conv *) let poly_f_conv tm = let ltm, x_tm = dest_comb tm in let list_tm = rand ltm in let inst_t = INST[x_tm, x_var] in let poly_f_cons, poly_f_empty = inst_t POLY_F_CONS, inst_t POLY_F_EMPTY in let rec poly_f_conv_raw list_tm = if (is_comb list_tm) then let h_tm, t_tm = dest_comb list_tm in let th0 = INST[rand h_tm, h_var; t_tm, t_var] poly_f_cons in let h_plus, rtm = dest_comb(rand(concl th0)) in let x_times = rator rtm in let th1 = poly_f_conv_raw t_tm in let th2 = TRANS th0 (AP_TERM h_plus (AP_TERM x_times th1)) in let th3 = (RAND_CONV my_real_rat_mul_conv THENC my_real_rat_add_conv) (rand(concl th2)) in TRANS th2 th3 else poly_f_empty in poly_f_conv_raw list_tm;; let poly_f_even_conv tm = let ltm, x_tm = dest_comb tm in let list_tm = rand ltm in let inst_t = INST[x_tm, x_var] in let poly_f_even_cons = inst_t POLY_F_EVEN_CONS and poly_f_even_empty = inst_t POLY_F_EVEN_EMPTY in let x2_times = AP_TERM mul_op_real (my_real_rat_mul_conv (mk_binop mul_op_real x_tm x_tm)) in let rec poly_f_even_conv_raw list_tm = if (is_comb list_tm) then let h_tm, t_tm = dest_comb list_tm in let th0 = INST[rand h_tm, h_var; t_tm, t_var] poly_f_even_cons in let h_plus, rtm = dest_comb(rand(concl th0)) in let th1 = poly_f_even_conv_raw t_tm in let th2 = TRANS th0 (AP_TERM h_plus (MK_COMB(x2_times, th1))) in let th3 = (RAND_CONV my_real_rat_mul_conv THENC my_real_rat_add_conv) (rand(concl th2)) in TRANS th2 th3 else poly_f_even_empty in poly_f_even_conv_raw list_tm;; let poly_f_odd_conv tm = let ltm, x_tm = dest_comb tm in let list_tm = rand ltm in let th0 = INST[list_tm, t_var; x_tm, x_var] poly_f_odd' in let ltm, rtm = dest_comb(rand(concl th0)) in let th1 = AP_TERM ltm (poly_f_even_conv rtm) in let th2 = my_real_rat_mul_conv (rand(concl th1)) in TRANS th0 (TRANS th1 th2);; (* let tm = `poly_f [-- &1436346 / &436346; -- &244664 / &4654235; &3 / &43545] (-- &144545 / &2345345)`;; let tm2 = (rand o concl) (REWRITE_CONV[poly_f; ITLIST] tm);; let tm' = replace_numerals tm;; poly_f_conv tm';; (* 3.280 *) test 10 (REWRITE_CONV [poly_f; ITLIST] THENC REAL_RAT_REDUCE_CONV) tm;; (* 3.188 *) test 10 REAL_RAT_REDUCE_CONV tm2;; (* 0.184 *) test 10 poly_f_conv tm';; *) (********************************) (* let REVERSE_TABLE = define `(REVERSE_TABLE (f:num->A) 0 = []) /\ (REVERSE_TABLE f (SUC i) = CONS (f i) ( REVERSE_TABLE f i))`;; let TABLE = new_definition `!(f:num->A) k. TABLE f k = REVERSE (REVERSE_TABLE f k)`;; let rec reverse_table_conv tm = let ltm, i_tm = dest_comb tm in if (i_tm = `0`) then ONCE_REWRITE_CONV[REVERSE_TABLE] tm else let i_suc = num_CONV i_tm in let th1 = ONCE_REWRITE_RULE[REVERSE_TABLE] (AP_TERM ltm i_suc) in let ltm, rtm = dest_comb (rand(concl th1)) in let th2 = reverse_table_conv rtm in TRANS th1 (AP_TERM ltm th2);; let cos_taylor n = let tm = mk_comb(`TABLE (\k. (if (EVEN k) then &1 else --(&1)) / &(FACT(2 * k)))`, mk_small_numeral n) in (rand o concl) ((REWRITE_CONV[TABLE] THENC ONCE_DEPTH_CONV reverse_table_conv THENC REWRITE_CONV[REVERSE; APPEND] THENC NUM_REDUCE_CONV) tm);; let cos_poly n x = let tm = mk_comb(mk_comb(`poly_f_even`, cos_taylor n), x) in let tm' = replace_numerals tm in tm';; let x = rand(concl(REWRITE_CONV[DECIMAL] `#1.230959417`));; let tm = cos_poly 6 x;; let tm' = (rand(concl(REPLACE_NUMS (REWRITE_CONV[poly_f_even; ITLIST] tm))));; (* 5.860 *) test 1 REAL_RAT_REDUCE_CONV tm';; (* 0.200 *) test 1 poly_f_even_conv tm;; let result = poly_f_even_conv tm;; REPLACE_NUMS result;; *)