(* ========================================================================= *)
(* Derivation of Machin's formula and other similar ones.                    *)
(* ========================================================================= *)
needs "Library/transc.ml";;
let REAL_LT_1_POW2 = prove
 (`!n. &1 < &2 pow n <=> ~(n = 0)`,
  GEN_TAC THEN ASM_CASES_TAC `n = 0` THEN ASM_REWRITE_TAC[] THEN
  CONV_TAC REAL_RAT_REDUCE_CONV THEN
  SUBST1_TAC(SYM(REAL_RAT_REDUCE_CONV `&2 pow 0`)) THEN
  MATCH_MP_TAC 
REAL_POW_MONO_LT THEN
  REWRITE_TAC[
REAL_OF_NUM_LT] THEN POP_ASSUM MP_TAC THEN ARITH_TAC);;
 
let REAL_POW2_CLAUSES = prove
 (`(!n. &0 <= &2 pow n) /\
   (!n. &0 < &2 pow n) /\
   (!n. &0 <= inv(&2 pow n)) /\
   (!n. &0 < inv(&2 pow n)) /\
   (!n. inv(&2 pow n) <= &1) /\
   (!n. &1 - inv(&2 pow n) <= &1) /\
   (!n. &1 <= &2 pow n) /\
   (!n. &1 < &2 pow n <=> ~(n = 0)) /\
   (!n. &0 <= &1 - inv(&2 pow n)) /\
   (!n. &0 <= &2 pow n - &1) /\
   (!n. &0 < &1 - inv(&2 pow n) <=> ~(n = 0))`,
 
let REAL_POW2_THM = prove
 (`&0 < &2 pow n /\
   &1 <= &2 pow n /\
   (&1 < &2 pow n <=> ~(n = 0)) /\
   (&2 pow m <= &2 pow n <=> m <= n) /\
   (&2 pow m < &2 pow n <=> m < n) /\
   (inv(&2 pow m) <= inv(&2 pow n) <=> n <= m) /\
   (inv(&2 pow m) < inv(&2 pow n) <=> n < m)`,
 
let REAL_ATN_POWSER_DIFFL = prove
 (`!x. abs(x) < &1
       ==> ((\x. suminf (\n. (if 
EVEN n then &0
                              else --(&1) pow ((n - 1) DIV 2) / &n) * x pow n))
            diffl (inv(&1 + x pow 2))) x`,
 
let REAL_ATN_POWSER = prove
 (`!x. abs(x) < &1
       ==> (\n. (if 
EVEN n then &0
                 else --(&1) pow ((n - 1) DIV 2) / &n) * x pow n)
           sums (atn x)`,
  REPEAT STRIP_TAC THEN
  FIRST_ASSUM(MP_TAC o MATCH_MP 
REAL_ATN_POWSER_SUMMABLE) THEN
  DISCH_THEN(MP_TAC o MATCH_MP 
SUMMABLE_SUM) THEN
  SUBGOAL_THEN
   `suminf (\n. (if 
EVEN n then &0
                 else --(&1) pow ((n - 1) DIV 2) / &n) * x pow n) = atn(x)`
   (fun th -> REWRITE_TAC[th]) THEN
  ONCE_REWRITE_TAC[REAL_ARITH `(a = b) <=> (a - b = &0)`] THEN
  SUBGOAL_THEN
   `suminf (\n. (if 
EVEN n then &0
                 else --(&1) pow ((n - 1) DIV 2) / &n) * &0 pow n) -
    atn(&0) = &0`
  MP_TAC THENL
   [MATCH_MP_TAC(REAL_ARITH `(a = &0) /\ (b = &0) ==> (a - b = &0)`) THEN
    CONJ_TAC THENL
     [CONV_TAC SYM_CONV THEN MATCH_MP_TAC 
SUM_UNIQ THEN
      MP_TAC(SPEC `&0` 
GP) THEN
      REWRITE_TAC[
REAL_ABS_NUM; 
REAL_OF_NUM_LT; ARITH] THEN
      DISCH_THEN(MP_TAC o SPEC `&0` o MATCH_MP 
SER_CMUL) THEN
      REWRITE_TAC[
REAL_MUL_LZERO] THEN
      MATCH_MP_TAC(TAUT `(a = b) ==> a ==> b`) THEN
      AP_THM_TAC THEN AP_TERM_TAC THEN ABS_TAC THEN
      COND_CASES_TAC THEN ASM_REWRITE_TAC[
REAL_MUL_LZERO] THEN
      CONV_TAC SYM_CONV THEN
      REWRITE_TAC[
REAL_ENTIRE; 
REAL_POW_EQ_0] THEN ASM_MESON_TAC[
EVEN];
      GEN_REWRITE_TAC (LAND_CONV o RAND_CONV) [GSYM 
TAN_0] THEN
      MATCH_MP_TAC 
TAN_ATN THEN
      SIMP_TAC[
PI2_BOUNDS; REAL_ARITH `&0 < x ==> --x < &0`]];
    ALL_TAC] THEN
  ASM_CASES_TAC `x = &0` THEN ASM_REWRITE_TAC[] THEN
  DISCH_THEN(fun th -> GEN_REWRITE_TAC RAND_CONV [SYM th]) THEN
  MP_TAC(SPEC `\x. suminf (\n. (if 
EVEN n then &0
                       else --(&1) pow ((n - 1) DIV 2) / &n) * x pow n) -
          atn x` 
DIFF_ISCONST_END_SIMPLE) THEN
  FIRST_ASSUM(DISJ_CASES_TAC o MATCH_MP (REAL_ARITH
    `~(x = &0) ==> &0 < x \/ x < &0`))
  THENL
   [DISCH_THEN(MP_TAC o SPECL [`&0`; `x:real`]);
    CONV_TAC(RAND_CONV SYM_CONV) THEN
    DISCH_THEN(MP_TAC o SPECL [`x:real`; `&0`])] THEN
  (REWRITE_TAC[] THEN DISCH_THEN MATCH_MP_TAC THEN
   ASM_REWRITE_TAC[] THEN
   X_GEN_TAC `u:real` THEN REPEAT STRIP_TAC THEN
   SUBGOAL_THEN `abs(u) < &1` (MP_TAC o MATCH_MP 
REAL_ATN_POWSER_DIFFL) THENL
    [POP_ASSUM_LIST(MP_TAC o end_itlist CONJ) THEN REAL_ARITH_TAC;
     ALL_TAC] THEN
   DISCH_THEN(MP_TAC o C CONJ (SPEC `u:real` 
DIFF_ATN)) THEN
   DISCH_THEN(MP_TAC o MATCH_MP 
DIFF_SUB) THEN
   REWRITE_TAC[
REAL_SUB_REFL]));;
 
let MCLAURIN_ATN_SIMPLE = prove
 (`!x n k. abs(x) <= inv(&2 pow k) /\ ~(k = 0)
           ==> abs(atn x -
                   sum(0,n) (\m. (if 
EVEN m then &0
                                  else --(&1) pow ((m - 1) DIV 2) / &m) *
                                  x pow m))
               <= &2 * abs(x) pow n`,
  REPEAT STRIP_TAC THEN SUBGOAL_THEN `abs(x) < &1` ASSUME_TAC THENL
   [MATCH_MP_TAC 
REAL_LET_TRANS THEN EXISTS_TAC `inv(&2 pow k)` THEN
    ASM_REWRITE_TAC[REAL_ARITH `a < &1 <=> &0 < &1 - a`; 
REAL_POW2_CLAUSES];
    ALL_TAC] THEN
  FIRST_ASSUM(MP_TAC o MATCH_MP 
REAL_ATN_POWSER) THEN
  DISCH_THEN(fun th -> ASSUME_TAC(SYM(MATCH_MP 
SUM_UNIQ th)) THEN
                       MP_TAC(MATCH_MP 
SUM_SUMMABLE th)) THEN
  DISCH_THEN(MP_TAC o MATCH_MP 
SER_OFFSET) THEN
  DISCH_THEN(MP_TAC o SPEC `n:num`) THEN ASM_REWRITE_TAC[] THEN
  DISCH_THEN(MP_TAC o MATCH_MP 
SUM_UNIQ) THEN
  MATCH_MP_TAC(REAL_ARITH
   `abs(r) <= e ==> (f - s = r) ==> abs(f - s) <= e`) THEN
  SUBGOAL_THEN
   `(\m. abs(x) pow (m + n)) sums (abs(x) pow n) * inv(&1 - abs(x))`
  ASSUME_TAC THENL
   [FIRST_ASSUM(MP_TAC o MATCH_MP 
GP o MATCH_MP (REAL_ARITH
      `abs(x) < &1 ==> abs(abs x) < &1`)) THEN
    DISCH_THEN(MP_TAC o SPEC `abs(x) pow n` o MATCH_MP 
SER_CMUL) THEN
    ONCE_REWRITE_TAC[
ADD_SYM] THEN REWRITE_TAC[GSYM 
REAL_POW_ADD];
    ALL_TAC] THEN
  MATCH_MP_TAC 
REAL_LE_TRANS THEN
  EXISTS_TAC `suminf (\m. abs(x) pow (m + n))` THEN CONJ_TAC THENL
   [ALL_TAC;
    FIRST_ASSUM(SUBST1_TAC o SYM o MATCH_MP 
SUM_UNIQ) THEN
    GEN_REWRITE_TAC RAND_CONV [REAL_MUL_SYM] THEN
    MATCH_MP_TAC 
REAL_LE_LMUL THEN SIMP_TAC[
REAL_ABS_POS; 
REAL_POW_LE] THEN
    GEN_REWRITE_TAC RAND_CONV [GSYM 
REAL_INV_INV] THEN
    MATCH_MP_TAC 
REAL_LE_INV2 THEN
    ONCE_REWRITE_TAC[REAL_ARITH `a <= &1 - b <=> b <= &1 - a`] THEN
    CONV_TAC REAL_RAT_REDUCE_CONV THEN
    MATCH_MP_TAC 
REAL_LE_TRANS THEN EXISTS_TAC `inv(&2 pow k)` THEN
    ASM_REWRITE_TAC[] THEN REWRITE_TAC[
real_div; REAL_MUL_LID] THEN
    MATCH_MP_TAC 
REAL_LE_INV2 THEN CONV_TAC REAL_RAT_REDUCE_CONV THEN
    GEN_REWRITE_TAC LAND_CONV [GSYM 
REAL_POW_1] THEN
    ASM_SIMP_TAC[
REAL_POW_MONO; REAL_OF_NUM_LE; ARITH;
                 ARITH_RULE `~(k = 0) ==> 1 <= k`]] THEN
  SUBGOAL_THEN
   `!m. abs((if 
EVEN (m + n) then &0
             else --(&1) pow (((m + n) - 1) DIV 2) / &(m + n)) *
             x pow (m + n))
        <= abs(x) pow (m + n)`
  ASSUME_TAC THENL
   [GEN_TAC THEN COND_CASES_TAC THEN
    SIMP_TAC[
REAL_MUL_LZERO; 
REAL_ABS_NUM; 
REAL_POW_LE; 
REAL_ABS_POS] THEN
    REWRITE_TAC[
REAL_ABS_MUL; 
REAL_ABS_DIV; 
REAL_ABS_POW; 
REAL_ABS_NEG] THEN
    REWRITE_TAC[
REAL_ABS_NUM; 
REAL_POW_ONE; REAL_MUL_LID] THEN
    GEN_REWRITE_TAC RAND_CONV [GSYM REAL_MUL_LID] THEN
    MATCH_MP_TAC 
REAL_LE_RMUL THEN SIMP_TAC[
REAL_POW_LE; 
REAL_ABS_POS] THEN
    REWRITE_TAC[
real_div; REAL_MUL_LID] THEN
    GEN_REWRITE_TAC RAND_CONV [GSYM 
REAL_INV_1] THEN
    MATCH_MP_TAC 
REAL_LE_INV2 THEN CONV_TAC REAL_RAT_REDUCE_CONV THEN
    REWRITE_TAC[REAL_OF_NUM_LE; ARITH_RULE `1 <= n <=> ~(n = 0)`] THEN
    ASM_MESON_TAC[
EVEN]; ALL_TAC] THEN
  MATCH_MP_TAC 
REAL_LE_TRANS THEN
  EXISTS_TAC
   `suminf (\m. abs((if 
EVEN (m + n) then &0
                     else --(&1) pow (((m + n) - 1) DIV 2) / &(m + n)) *
                    x pow (m + n)))` THEN
  CONJ_TAC THENL
   [MATCH_MP_TAC 
SER_ABS THEN MATCH_MP_TAC 
SER_COMPARA THEN
    EXISTS_TAC `\m. abs(x) pow (m + n)` THEN
    ASM_REWRITE_TAC[] THEN ASM_MESON_TAC[
SUM_SUMMABLE]; ALL_TAC] THEN
  MATCH_MP_TAC 
SER_LE THEN ASM_REWRITE_TAC[] THEN CONJ_TAC THENL
   [MATCH_MP_TAC 
SER_COMPARA THEN
    EXISTS_TAC `\m. abs(x) pow (m + n)` THEN
    ASM_REWRITE_TAC[]; ALL_TAC] THEN
  ASM_MESON_TAC[
SUM_SUMMABLE]);;
 
let MCLAURIN_ATN_APPROX = prove
 (`!x n k. abs(x) <= inv(&2 pow k) /\ ~(k = 0)
           ==> abs(atn x -
                   sum(0,n) (\m. (if 
EVEN m then &0
                                  else --(&1) pow ((m - 1) DIV 2) / &m) *
                                  x pow m))
               <= inv(&2 pow (n * k - 1))`,
 
let TAN_ADD_ATN_SIDECOND = prove
 (`!x y. ~(x * y = &1) ==> ~(cos(atn x + atn y) = &0)`,
  REPEAT GEN_TAC THEN DISCH_TAC THEN
  REWRITE_TAC[
COS_ADD; REAL_ARITH `(a - b = &0) <=> (a = b)`] THEN
  DISCH_THEN(MP_TAC o AP_TERM `(*) (inv(cos(atn x)))`) THEN
  SIMP_TAC[REAL_MUL_ASSOC; REAL_MUL_LINV; COS_ATN_NZ; REAL_MUL_LID] THEN
  DISCH_THEN(MP_TAC o AP_TERM `(*) (inv(cos(atn y)))`) THEN
  SIMP_TAC[REAL_MUL_LINV; 
COS_ATN_NZ; REAL_MUL_LID; GSYM REAL_MUL_ASSOC] THEN
  ONCE_REWRITE_TAC[AC 
REAL_MUL_AC `a * b * c * d = (c * b) * (d * a)`] THEN
  ASM_REWRITE_TAC[GSYM tan; GSYM 
real_div; 
ATN_TAN]);;
 
let ATN_ADD = prove
 (`!x y. ~(x * y = &1) /\
         abs(atn x + atn y) < pi / &2
         ==> (atn(x) + atn(y) = atn((x + y) / (&1 - x * y)))`,
  REPEAT STRIP_TAC THEN
  SUBGOAL_THEN `tan(atn(x) + atn(y)) = (x + y) / (&1 - x * y)` MP_TAC THENL
   [ASM_SIMP_TAC[
ATN_TAN; 
TAN_ADD; 
COS_ATN_NZ; 
TAN_ADD_ATN_SIDECOND];
    DISCH_THEN(SUBST1_TAC o SYM) THEN CONV_TAC SYM_CONV THEN
    ASM_SIMP_TAC[
TAN_ATN; REAL_ARITH `abs(x) < e ==> --e < x /\ x < e`]]);;
 
let ATN_ADD_SMALL_LEMMA = prove
 (`!x y. abs(x * y) < &1 ==> abs(atn(x) + atn(y)) < pi / &2`,
  REPEAT STRIP_TAC THEN
  MATCH_MP_TAC(REAL_ARITH
   `--a < x /\ x < a /\ --a < y /\ y < a /\
    (&0 < y ==> x + y < a) /\
    (&0 < --y ==> --x + --y < a)
    ==> abs(x + y) < a`) THEN
  REWRITE_TAC[
ATN_BOUNDS] THEN CONJ_TAC THEN
  REPEAT STRIP_TAC THEN REWRITE_TAC[GSYM 
ATN_NEG] THEN
  MATCH_MP_TAC 
ATN_ADD_SMALL_LEMMA_POS THEN
  ASM_SIMP_TAC[REAL_ARITH `abs(x) < &1 ==> x < &1`;
               REAL_ARITH `--x * -- y = x * y`] THEN
  FIRST_ASSUM(MATCH_MP_TAC o MATCH_MP (REAL_ARITH
   `&0 < y ==> (z <= &0 ==> y <= &0) ==> &0 < z`)) THEN
  MATCH_MP_TAC(REAL_ARITH
   `(y < &0 ==> z < &0) /\ ((y = &0) ==> (z = &0))
    ==> y <= &0 ==> z <= &0`) THEN
  SIMP_TAC[
ATN_0; GSYM 
ATN_NEG] THEN
  GEN_REWRITE_TAC (RAND_CONV o RAND_CONV) [GSYM 
ATN_0] THEN
  SIMP_TAC[
ATN_MONO_LT]);;
 
let pth_base = prove
   (`(&0 * atn(x) = &0) /\
     (&1 * atn(x) = atn(x))`,
    REWRITE_TAC[
REAL_MUL_LZERO; REAL_MUL_LID])
  and 
pth_0,
pth_1 = (CONJ_PAIR o 
prove)
   (`(&(
NUMERAL(
BIT0 n)) * atn(x) =
      &(
NUMERAL n) * atn(x) + &(
NUMERAL n) * atn(x)) /\
     (&(
NUMERAL(
BIT1 n)) * atn(x) =
      atn(x) + &(
NUMERAL n) * atn(x) + &(
NUMERAL n) * atn(x))`,
    REWRITE_TAC[
REAL_MUL_LZERO; REAL_MUL_LID] THEN
    REWRITE_TAC[
NUMERAL; 
BIT0; 
BIT1] THEN
    REWRITE_TAC[GSYM REAL_OF_NUM_ADD; GSYM 
REAL_OF_NUM_SUC] THEN
    REAL_ARITH_TAC) in
  let rewr_base = GEN_REWRITE_CONV I [
pth_base]
  and rewr_0 = GEN_REWRITE_CONV I [
pth_0]
  and rewr_1 = GEN_REWRITE_CONV I [
pth_1] in
  let rec ATN_CMUL_CONV tm =
    if not (is_ratconst(rand(rand tm))) then failwith "ATN_CMUL_CONV" else
    try rewr_base tm with Failure _ ->
    try let th1 = rewr_0 tm in
        let tm1 = rand(concl th1) in
        let th2 = ATN_CMUL_CONV(rand tm1) in
        let th3 = MK_COMB(AP_TERM (rator(rator tm1)) th2,th2) in
        let th4 = TRANS th3 (ATN_ADD_CONV(rand(concl th3))) in
        TRANS th1 th4
    with Failure _ ->
        let th1 = rewr_1 tm in
        let tm1 = rand(rand(concl th1)) in
        let th2 = ATN_CMUL_CONV(rand tm1) in
        let th3 = MK_COMB(AP_TERM (rator(rator tm1)) th2,th2) in
        let th4 = TRANS th3 (ATN_ADD_CONV(rand(concl th3))) in
        let th5 = AP_TERM (rator(rand(concl th1))) th4 in
        let th6 = TRANS th5 (ATN_ADD_CONV(rand(concl th5))) in
        TRANS th1 th6 in
  ATN_CMUL_CONV;;
 
let pth = prove
   (`(atn(x) - atn(y) = atn(x) + atn(--y))`,
    REWRITE_TAC[
real_sub; 
ATN_NEG]) in
  GEN_REWRITE_CONV I [pth] THENC
  RAND_CONV(RAND_CONV REAL_RAT_NEG_CONV) THENC
  ATN_ADD_CONV;;
 
let pth = prove
   (`!p. abs(atn(&0) - &0) <= inv(&2 pow p)`,
    SIMP_TAC[
ATN_0; 
REAL_SUB_REFL; 
REAL_ABS_NUM; 
REAL_LE_INV_EQ;
             
REAL_POW_LE; 
REAL_POS]) in
  let atn_approx_conv p r =
    if r =/ num_0 then num_0 else
    let k = Num.int_of_num(minus_num(log_2(abs_num r))) in
    if k < 1 then failwith "atn_approx_conv: argument too big" else
    let rats = mclaurin_atn_rule k p in
    
POLY rats r
  and ATN_APPROX_CONV p tm =
    let atm,rtm = dest_comb tm in
    if atm <> atn_tm then failwith "ATN_APPROX_CONV" else
    let r = rat_of_term rtm in
    if r =/ num_0 then SPEC (mk_small_numeral p) pth else
    let k = Num.int_of_num(minus_num(log_2(abs_num r))) in
    if k < 1 then failwith "ATN_APPROX_CONV: argument too big" else
    let th1 = MCLAURIN_ATN_RULE k p in
    let th2 = SPEC rtm th1 in
    let th3 = MP th2 (EQT_ELIM(REAL_RAT_REDUCE_CONV(lhand(concl th2)))) in
    CONV_RULE(LAND_CONV(RAND_CONV(RAND_CONV REAL_RAT_REDUCE_CONV))) th3 in
  atn_approx_conv,ATN_APPROX_CONV;;
 
let pth = prove
   (`(q1 = p + 5) /\
     (q2 = p + 6) /\
     abs(atn(&1 / &8) - a1) <= inv(&2 pow q1) /\
     abs(atn(&1 / &57) - a2) <= inv(&2 pow q2) /\
     abs(atn(&1 / &239) - a3) <= inv(&2 pow q2)
     ==> abs(pi - (&24 * a1 + &8 * a2 + &4 * a3)) <= inv(&2 pow p)`,
    DISCH_THEN(CONJUNCTS_THEN2 ASSUME_TAC MP_TAC) THEN
    DISCH_THEN(CONJUNCTS_THEN2 ASSUME_TAC MP_TAC) THEN
    ASM_REWRITE_TAC[] THEN POP_ASSUM_LIST(K ALL_TAC) THEN STRIP_TAC THEN
    MATCH_MP_TAC 
REAL_LE_LCANCEL_IMP THEN EXISTS_TAC `abs(inv(&2 pow 2))` THEN
    SIMP_TAC[
REAL_POW2_CLAUSES; REAL_ARITH `&0 < x ==> &0 < abs(x)`] THEN
    REWRITE_TAC[GSYM 
REAL_ABS_MUL] THEN
    REWRITE_TAC[
REAL_ABS_INV; 
REAL_ABS_POW; 
REAL_ABS_NUM] THEN
    REWRITE_TAC[GSYM 
REAL_INV_MUL; GSYM 
REAL_POW_ADD] THEN
    CONV_TAC REAL_RAT_REDUCE_CONV THEN
    REWRITE_TAC[
REAL_SUB_LDISTRIB; REAL_ADD_LDISTRIB; REAL_MUL_ASSOC] THEN
    CONV_TAC REAL_RAT_REDUCE_CONV THEN
    REWRITE_TAC[
real_div; REAL_MUL_LID] THEN
    GEN_REWRITE_TAC (LAND_CONV o RAND_CONV o LAND_CONV) [REAL_MUL_SYM] THEN
    REWRITE_TAC[GSYM 
real_div; MACHINLIKE_1] THEN
    REWRITE_TAC[REAL_ARITH `(x1 + x2 + x3) - (y1 + y2 + y3) =
                            (x1 - y1) + (x2 - y2) + (x3 - y3)`] THEN
    REWRITE_TAC[GSYM 
REAL_SUB_LDISTRIB] THEN BOUND_SUMPROD_TAC THEN
    GEN_REWRITE_TAC (LAND_CONV o ONCE_DEPTH_CONV) [
ADD_SYM] THEN
    REWRITE_TAC[
REAL_POW_ADD; 
REAL_INV_MUL; REAL_MUL_ASSOC] THEN
    SIMP_TAC[GSYM 
REAL_ADD_RDISTRIB; 
REAL_LE_RMUL_EQ; 
REAL_POW2_CLAUSES] THEN
    CONV_TAC REAL_RAT_REDUCE_CONV) in
  let pi_approx_rule p =
    let q1 = p + 5
    and q2 = p + 6 in
    let a1 = atn_approx_conv q1 const_1_8
    and a2 = atn_approx_conv q2 const_1_57
    and a3 = atn_approx_conv q2 const_1_239 in
    const_24 */ a1 +/ const_8 */ a2 +/ const_4 */ a3
  and PI_APPROX_RULE p =
    let q1 = p + 5
    and q2 = p + 6 in
    let th1 = ATN_APPROX_CONV q1 tm_1_8
    and th2 = ATN_APPROX_CONV q2 tm_1_57
    and th3 = ATN_APPROX_CONV q2 tm_1_239 in
    let th4 = INST [mk_small_numeral p,p_tm;
                    mk_small_numeral q1,q1_tm;
                    mk_small_numeral q2,q2_tm] pth in
    let th5 = EQT_ELIM(NUM_REDUCE_CONV(lhand(lhand(concl th4))))
    and th6 = EQT_ELIM(NUM_REDUCE_CONV(lhand(rand(lhand(concl th4))))) in
    let th7 = MATCH_MP th4 (end_itlist CONJ [th5; th6; th1; th2; th3]) in
    CONV_RULE(LAND_CONV(RAND_CONV(RAND_CONV REAL_RAT_REDUCE_CONV))) th7 in
  pi_approx_rule,PI_APPROX_RULE;;
 
let pth = prove
   (`abs(x - r) <= inv(&2 pow (SUC p))
     ==> !a. abs(&2 pow p * r - a) <= inv(&2)
             ==> abs(x - a / &2 pow p) <= inv(&2 pow p)`,
    REPEAT STRIP_TAC THEN
    FIRST_ASSUM(MATCH_MP_TAC o MATCH_MP (REAL_ARITH
      `abs(x - r) <= q ==> abs(r - r') <= p - q ==> abs(x - r') <= p`)) THEN
    MATCH_MP_TAC 
REAL_LE_LCANCEL_IMP THEN EXISTS_TAC `abs(&2 pow p)` THEN
    SIMP_TAC[REAL_ARITH `&0 < x ==> &0 < abs(x)`; 
REAL_POW2_THM] THEN
    REWRITE_TAC[GSYM 
REAL_ABS_MUL; 
REAL_SUB_LDISTRIB] THEN
    SIMP_TAC[
REAL_ABS_POW; 
REAL_ABS_NUM; GSYM 
real_div;
             
REAL_DIV_LMUL; 
REAL_LT_IMP_NZ; 
REAL_POW2_CLAUSES;
             
REAL_DIV_POW2; REAL_OF_NUM_EQ; 
ARITH_EQ;
             
LE_REFL; ARITH_RULE `~(SUC p <= p)`;
             ARITH_RULE `SUC p - p = 1`; 
SUB_REFL] THEN
    UNDISCH_TAC `abs (&2 pow p * r - a) <= inv (&2)` THEN
    CONV_TAC NUM_REDUCE_CONV THEN CONV_TAC REAL_RAT_REDUCE_CONV)
  and num_2 = Int 2 in
  let pi_approx_binary_rule p =
    let ppow = power_num num_2 (Int p) in
    let r = pi_approx_rule (p + 1) in
    let a = round_num (ppow */ r) in
    a // ppow
  and PI_APPROX_BINARY_RULE p =
    let ppow = power_num num_2 (Int p) in
    let th1 = PI_APPROX_RULE (p + 1) in
    let th2 = CONV_RULE(funpow 3 RAND_CONV num_CONV) th1 in
    let r = rat_of_term(rand(rand(lhand(concl th2)))) in
    let th3 = SPEC (mk_realintconst(round_num(ppow */ r))) (MATCH_MP pth th2) in
    let th4 = MP th3 (EQT_ELIM(REAL_RAT_REDUCE_CONV(lhand(concl th3)))) in
    CONV_RULE(LAND_CONV(RAND_CONV(RAND_CONV REAL_RAT_REDUCE_CONV))) th4 in
  pi_approx_binary_rule,PI_APPROX_BINARY_RULE;;