Examples/machin.ml

   1 (* ========================================================================= *)
   2 (* Derivation of Machin's formula and other similar ones.                    *)
   3 (* ========================================================================= *)
   4
   5 needs "Library/transc.ml";;
   6
   7 let REAL_LE_1_POW2 = prove
   8  (`!n. &1 <= &2 pow n`,
   9   REWRITE_TAC[REAL_OF_NUM_POW; REAL_OF_NUM_LE; ARITH_RULE `1 <= n <=> 0 < n`;
  10               EXP_LT_0; ARITH]);;
  11
  12 let REAL_LT_1_POW2 = prove
  13  (`!n. &1 < &2 pow n <=> ~(n = 0)`,
  14   GEN_TAC THEN ASM_CASES_TAC `n = 0` THEN ASM_REWRITE_TAC[] THEN
  15   CONV_TAC REAL_RAT_REDUCE_CONV THEN
  16   SUBST1_TAC(SYM(REAL_RAT_REDUCE_CONV `&2 pow 0`)) THEN
  17   MATCH_MP_TAC REAL_POW_MONO_LT THEN
  18   REWRITE_TAC[REAL_OF_NUM_LT] THEN POP_ASSUM MP_TAC THEN ARITH_TAC);;
  19
  20 let REAL_POW2_CLAUSES = prove
  21  (`(!n. &0 <= &2 pow n) /\
  22    (!n. &0 < &2 pow n) /\
  23    (!n. &0 <= inv(&2 pow n)) /\
  24    (!n. &0 < inv(&2 pow n)) /\
  25    (!n. inv(&2 pow n) <= &1) /\
  26    (!n. &1 - inv(&2 pow n) <= &1) /\
  27    (!n. &1 <= &2 pow n) /\
  28    (!n. &1 < &2 pow n <=> ~(n = 0)) /\
  29    (!n. &0 <= &1 - inv(&2 pow n)) /\
  30    (!n. &0 <= &2 pow n - &1) /\
  31    (!n. &0 < &1 - inv(&2 pow n) <=> ~(n = 0))`,
  32   SIMP_TAC[REAL_LE_1_POW2; REAL_LT_1_POW2; REAL_SUB_LE; REAL_SUB_LT;
  33            REAL_INV_LE_1] THEN
  34   SIMP_TAC[REAL_LE_INV_EQ; REAL_LT_INV_EQ; REAL_POW_LT; REAL_POW_LE;
  35            REAL_OF_NUM_LE; REAL_OF_NUM_LT; ARITH;
  36            REAL_ARITH `&1 - x <= &1 <=> &0 <= x`] THEN
  37   GEN_TAC THEN ASM_CASES_TAC `n = 0` THEN ASM_REWRITE_TAC[] THEN
  38   CONV_TAC REAL_RAT_REDUCE_CONV THEN
  39   MATCH_MP_TAC REAL_LET_TRANS THEN EXISTS_TAC `inv(&2 pow 1)` THEN
  40   ASM_SIMP_TAC[REAL_LE_INV2; REAL_POW_MONO; REAL_POW_LT; REAL_OF_NUM_LT; ARITH;
  41                REAL_OF_NUM_LE; ARITH_RULE `1 <= n <=> ~(n = 0)`] THEN
  42   CONV_TAC REAL_RAT_REDUCE_CONV);;
  43
  44 let REAL_POW2_THM = prove
  45  (`&0 < &2 pow n /\
  46    &1 <= &2 pow n /\
  47    (&1 < &2 pow n <=> ~(n = 0)) /\
  48    (&2 pow m <= &2 pow n <=> m <= n) /\
  49    (&2 pow m < &2 pow n <=> m < n) /\
  50    (inv(&2 pow m) <= inv(&2 pow n) <=> n <= m) /\
  51    (inv(&2 pow m) < inv(&2 pow n) <=> n < m)`,
  52   REWRITE_TAC[REAL_POW2_CLAUSES] THEN
  53   SUBGOAL_THEN `!m n. &2 pow m <= &2 pow n <=> m <= n` ASSUME_TAC THENL
  54    [REPEAT GEN_TAC THEN EQ_TAC THEN
  55     SIMP_TAC[REAL_POW_MONO; REAL_OF_NUM_LE; ARITH] THEN
  56     CONV_TAC CONTRAPOS_CONV THEN
  57     SIMP_TAC[REAL_NOT_LE; REAL_NOT_LT; REAL_POW_MONO_LT; REAL_OF_NUM_LT;
  58              NOT_LE; ARITH]; ALL_TAC] THEN
  59   ASM_REWRITE_TAC[GSYM REAL_NOT_LE] THEN REWRITE_TAC[GSYM NOT_LE] THEN
  60   SUBGOAL_THEN `!m n. inv(&2 pow m) <= inv(&2 pow n) <=> &2 pow n <= &2 pow m`
  61    (fun th -> ASM_REWRITE_TAC[th]) THEN
  62   REPEAT GEN_TAC THEN EQ_TAC THEN
  63   SIMP_TAC[REAL_LE_INV2; REAL_POW2_CLAUSES] THEN
  64   CONV_TAC CONTRAPOS_CONV THEN
  65   SIMP_TAC[REAL_NOT_LE; REAL_LT_INV2; REAL_POW2_CLAUSES]);;
  66
  67 (* ------------------------------------------------------------------------- *)
  68 (* Compound errors given bounds in assumptions.                              *)
  69 (* ------------------------------------------------------------------------- *)
  70
  71 let BOUND_SUMPROD_RULE =
  72   let pth_add = REAL_ARITH
  73     `abs(x1) <= b1 /\ abs(x2) <= b2 ==> abs(x1 + x2) <= b1 + b2`
  74   and pth_sub = REAL_ARITH
  75     `abs(x1) <= b1 /\ abs(x2) <= b2 ==> abs(x1 - x2) <= b1 + b2`
  76   and pth_mul = prove
  77    (`abs(x1) <= b1 /\ abs(x2) <= b2 ==> abs(x1 * x2) <= b1 * b2`,
  78     REWRITE_TAC[REAL_ABS_MUL] THEN
  79     SIMP_TAC[REAL_LE_MUL2; REAL_ABS_POS])
  80   and pth_neg = REAL_ARITH
  81     `abs(x1) <= b1 ==> abs(--x1) <= b1`
  82   and pth_pow = prove
  83    (`abs(x) <= b1 ==> abs(x pow n) <= b1 pow n`,
  84     REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[REAL_ABS_POW] THEN
  85     MATCH_MP_TAC REAL_POW_LE2 THEN ASM_REWRITE_TAC[REAL_ABS_POS])
  86   and pth_abs = REAL_ARITH `abs(x) <= b ==> abs(abs(x)) <= b`
  87   and pth_triv = REAL_ARITH `abs(x) <= abs(x)`
  88   and n_tm = `n:num` in
  89   let rec BOUND_SUMPROD_RULE (asl,w) =
  90     let tm = rator w in
  91     try tryfind (fun (_,th) -> if rator(concl th) = tm then th
  92                                else fail()) asl
  93     with Failure _ -> try
  94         let pth,th = tryfind
  95           (fun pth -> pth,PART_MATCH (rator o rand) pth tm)
  96           [pth_neg; pth_abs] in
  97         let th1 = BOUND_SUMPROD_RULE (asl,lhand(concl th)) in
  98         MATCH_MP pth th1
  99     with Failure _ -> try
 100         let pth = INST [funpow 3 rand tm,n_tm] pth_pow in
 101         let th = PART_MATCH (rator o rand) pth tm in
 102         let th1 = BOUND_SUMPROD_RULE (asl,lhand(concl th)) in
 103         MATCH_MP (INST [funpow 3 rand tm,n_tm] pth_pow) th1
 104     with Failure _ -> try
 105         let pth,th = tryfind
 106           (fun pth -> pth,PART_MATCH (rator o rand) pth tm)
 107           [pth_add; pth_sub; pth_mul] in
 108         let trm = lhand(concl th) in
 109         let th1 = BOUND_SUMPROD_RULE (asl,lhand trm)
 110         and th2 = BOUND_SUMPROD_RULE (asl,rand trm) in
 111         MATCH_MP pth (CONJ th1 th2)
 112     with Failure _ ->
 113         PART_MATCH rator pth_triv tm in
 114   BOUND_SUMPROD_RULE;;
 115
 116 let BOUND_SUMPROD_TAC =
 117   let tac =
 118     let pth =
 119       REAL_ARITH `x <= a ==> (!b. a <= b ==> x <= b) /\
 120                              (!b. a < b ==> x < b)` in
 121     fun th ->
 122       let th1,th2 = CONJ_PAIR(MATCH_MP pth th) in
 123       MATCH_MP_TAC th1 ORELSE MATCH_MP_TAC th2
 124   and le_tm = `(<=):real->real->bool` in
 125   fun (asl,w as gl) ->
 126     let l,r = dest_comb w in
 127     let gv = genvar(type_of r) in
 128     let tm = mk_comb(mk_comb(le_tm,rand l),gv) in
 129     let th = BOUND_SUMPROD_RULE(asl,tm) in
 130     tac th gl;;
 131
 132 (* ------------------------------------------------------------------------- *)
 133 (* Power series for atn.                                                     *)
 134 (* ------------------------------------------------------------------------- *)
 135
 136 let REAL_ATN_POWSER_SUMMABLE = prove
 137  (`!x. abs(x) < &1
 138        ==> summable (\n. (if EVEN n then &0
 139                           else --(&1) pow ((n - 1) DIV 2) / &n) * x pow n)`,
 140   REPEAT STRIP_TAC THEN MATCH_MP_TAC SER_COMPAR THEN
 141   EXISTS_TAC `\n. abs(x) pow n` THEN CONJ_TAC THENL
 142    [EXISTS_TAC `0` THEN REPEAT STRIP_TAC THEN REWRITE_TAC[] THEN
 143     COND_CASES_TAC THEN
 144     SIMP_TAC[REAL_POW_LE; REAL_MUL_LZERO; REAL_ABS_POS; REAL_ABS_NUM] THEN
 145     REWRITE_TAC[REAL_ABS_MUL; REAL_ABS_DIV; REAL_ABS_NEG; REAL_ABS_POW] THEN
 146     REWRITE_TAC[REAL_ABS_NUM; REAL_POW_ONE; REAL_MUL_LID] THEN
 147     REWRITE_TAC[real_div; REAL_MUL_LID] THEN
 148     ONCE_REWRITE_TAC[REAL_MUL_SYM] THEN REWRITE_TAC[GSYM real_div] THEN
 149     MATCH_MP_TAC REAL_LE_LDIV THEN
 150     CONJ_TAC THENL [ASM_MESON_TAC[REAL_OF_NUM_LT; EVEN; LT_NZ]; ALL_TAC] THEN
 151     GEN_REWRITE_TAC LAND_CONV [GSYM REAL_MUL_RID] THEN
 152     MATCH_MP_TAC REAL_LE_LMUL THEN
 153     SIMP_TAC[REAL_POW_LE; REAL_ABS_POS] THEN
 154     ASM_MESON_TAC[REAL_OF_NUM_LE; EVEN; ARITH_RULE `1 <= n <=> ~(n = 0)`];
 155     ALL_TAC] THEN
 156   REWRITE_TAC[summable] THEN EXISTS_TAC `inv(&1 - abs x)` THEN
 157   MATCH_MP_TAC GP THEN ASM_REWRITE_TAC[REAL_ABS_ABS]);;
 158
 159 let REAL_ATN_POWSER_DIFFS_SUMMABLE = prove
 160  (`!x. abs(x) < &1
 161        ==> summable (\n. diffs (\n. (if EVEN n then &0
 162                                      else --(&1) pow ((n - 1) DIV 2) / &n)) n *
 163                          x pow n)`,
 164   REPEAT STRIP_TAC THEN REWRITE_TAC[diffs] THEN
 165   MATCH_MP_TAC SER_COMPAR THEN
 166   EXISTS_TAC `\n. abs(x) pow n` THEN CONJ_TAC THENL
 167    [EXISTS_TAC `0` THEN REPEAT STRIP_TAC THEN REWRITE_TAC[] THEN
 168     COND_CASES_TAC THEN
 169     SIMP_TAC[REAL_POW_LE; REAL_MUL_LZERO; REAL_MUL_RZERO;
 170              REAL_ABS_POS; REAL_ABS_NUM] THEN
 171     SIMP_TAC[REAL_MUL_ASSOC; REAL_DIV_LMUL; REAL_OF_NUM_EQ; NOT_SUC] THEN
 172     REWRITE_TAC[REAL_ABS_MUL; REAL_ABS_DIV; REAL_ABS_NEG; REAL_ABS_POW] THEN
 173     REWRITE_TAC[REAL_ABS_NUM; REAL_POW_ONE; REAL_MUL_LID; REAL_LE_REFL];
 174     ALL_TAC] THEN
 175   REWRITE_TAC[summable] THEN EXISTS_TAC `inv(&1 - abs x)` THEN
 176   MATCH_MP_TAC GP THEN ASM_REWRITE_TAC[REAL_ABS_ABS]);;
 177
 178 let REAL_ATN_POWSER_DIFFS_SUM = prove
 179  (`!x. abs(x) < &1
 180        ==> (\n. diffs (\n. (if EVEN n then &0
 181                             else --(&1) pow ((n - 1) DIV 2) / &n)) n * x pow n)
 182            sums (inv(&1 + x pow 2))`,
 183   REPEAT STRIP_TAC THEN
 184   FIRST_ASSUM(MP_TAC o MATCH_MP REAL_ATN_POWSER_DIFFS_SUMMABLE) THEN
 185   DISCH_THEN(fun th -> MP_TAC(MATCH_MP SUMMABLE_SUM th) THEN
 186                        MP_TAC(MATCH_MP SER_PAIR th)) THEN
 187   SUBGOAL_THEN
 188    `(\n. sum (2 * n,2) (\n. diffs
 189       (\n. (if EVEN n then &0
 190             else --(&1) pow ((n - 1) DIV 2) / &n)) n * x pow n)) =
 191     (\n. --(x pow 2) pow n)`
 192   SUBST1_TAC THENL
 193    [ABS_TAC THEN
 194     CONV_TAC(LAND_CONV(LAND_CONV(RAND_CONV(TOP_DEPTH_CONV num_CONV)))) THEN
 195     REWRITE_TAC[sum; diffs; ADD_CLAUSES; EVEN_MULT; ARITH_EVEN; EVEN] THEN
 196     REWRITE_TAC[REAL_ADD_LID; REAL_ADD_RID; REAL_MUL_LZERO;
 197                 REAL_MUL_RZERO] THEN
 198     SIMP_TAC[ARITH_RULE `SUC n - 1 = n`; DIV_MULT; ARITH_EQ] THEN
 199     SIMP_TAC[REAL_MUL_ASSOC; REAL_DIV_LMUL; REAL_OF_NUM_EQ; NOT_SUC] THEN
 200     ONCE_REWRITE_TAC[GSYM REAL_POW_POW] THEN
 201     REWRITE_TAC[GSYM REAL_POW_MUL] THEN
 202     REWRITE_TAC[REAL_MUL_LNEG; REAL_MUL_LID]; ALL_TAC] THEN
 203   SUBGOAL_THEN `(\n. --(x pow 2) pow n) sums inv (&1 + x pow 2)` MP_TAC THENL
 204    [ONCE_REWRITE_TAC[REAL_ARITH `&1 + x = &1 - (--x)`] THEN
 205     MATCH_MP_TAC GP THEN
 206     REWRITE_TAC[REAL_ABS_NEG; REAL_ABS_POW] THEN
 207     GEN_REWRITE_TAC RAND_CONV [GSYM REAL_MUL_LID] THEN
 208     ASM_SIMP_TAC[REAL_POW_2; REAL_LT_MUL2; REAL_ABS_POS]; ALL_TAC] THEN
 209   MESON_TAC[SUM_UNIQ]);;
 210
 211 let REAL_ATN_POWSER_DIFFS_DIFFS_SUMMABLE = prove
 212  (`!x. abs(x) < &1
 213        ==> summable
 214              (\n. diffs (diffs
 215                  (\n. (if EVEN n then &0
 216                        else --(&1) pow ((n - 1) DIV 2) / &n))) n * x pow n)`,
 217   REPEAT STRIP_TAC THEN REWRITE_TAC[diffs] THEN
 218   MATCH_MP_TAC SER_COMPAR THEN
 219   EXISTS_TAC `\n. &(SUC n) * abs(x) pow n` THEN CONJ_TAC THENL
 220    [EXISTS_TAC `0` THEN REPEAT STRIP_TAC THEN
 221     REWRITE_TAC[REAL_ABS_NUM; REAL_ABS_MUL; GSYM REAL_MUL_ASSOC] THEN
 222     MATCH_MP_TAC REAL_LE_LMUL THEN REWRITE_TAC[REAL_POS] THEN
 223     COND_CASES_TAC THEN
 224     SIMP_TAC[REAL_POW_LE; REAL_MUL_LZERO; REAL_MUL_RZERO;
 225              REAL_ABS_POS; REAL_ABS_NUM] THEN
 226     REWRITE_TAC[REAL_ABS_DIV; REAL_ABS_NUM; REAL_MUL_ASSOC] THEN
 227     SIMP_TAC[REAL_DIV_LMUL; REAL_OF_NUM_EQ; NOT_SUC] THEN
 228     REWRITE_TAC[REAL_ABS_POW; REAL_ABS_NEG; REAL_POW_ONE; REAL_MUL_LID;
 229                 REAL_ABS_NUM; REAL_LE_REFL]; ALL_TAC] THEN
 230   MATCH_MP_TAC SER_RATIO THEN
 231   SUBGOAL_THEN `?c. abs(x) < c /\ c < &1` STRIP_ASSUME_TAC THENL
 232    [EXISTS_TAC `(&1 + abs(x)) / &2` THEN
 233     SIMP_TAC[REAL_LT_LDIV_EQ; REAL_LT_RDIV_EQ; REAL_OF_NUM_LT; ARITH] THEN
 234     UNDISCH_TAC `abs(x) < &1` THEN REAL_ARITH_TAC; ALL_TAC] THEN
 235   EXISTS_TAC `c:real` THEN ASM_REWRITE_TAC[] THEN
 236   SUBGOAL_THEN `?N. !n. n >= N ==> &(SUC(SUC n)) * abs(x) <= &(SUC n) * c`
 237   STRIP_ASSUME_TAC THENL
 238    [ALL_TAC;
 239     EXISTS_TAC `N:num` THEN REPEAT STRIP_TAC THEN
 240     REWRITE_TAC[real_pow; REAL_ABS_MUL; REAL_MUL_ASSOC] THEN
 241     MATCH_MP_TAC REAL_LE_RMUL THEN REWRITE_TAC[REAL_ABS_POS] THEN
 242     REWRITE_TAC[REAL_ABS_NUM; REAL_ABS_ABS] THEN
 243     GEN_REWRITE_TAC RAND_CONV [REAL_MUL_SYM] THEN ASM_SIMP_TAC[]] THEN
 244   ASM_CASES_TAC `x = &0` THENL
 245    [ASM_REWRITE_TAC[REAL_ABS_NUM; REAL_MUL_RZERO] THEN
 246     EXISTS_TAC `0` THEN REPEAT STRIP_TAC THEN MATCH_MP_TAC REAL_LE_MUL THEN
 247     REWRITE_TAC[REAL_POS] THEN UNDISCH_TAC `abs(x) < c` THEN REAL_ARITH_TAC;
 248     ALL_TAC] THEN
 249   ASM_SIMP_TAC[GSYM REAL_LE_RDIV_EQ; GSYM REAL_ABS_NZ] THEN
 250   REWRITE_TAC[real_div; GSYM REAL_MUL_ASSOC] THEN
 251   REWRITE_TAC[GSYM real_div] THEN
 252   REWRITE_TAC[ADD1; GSYM REAL_OF_NUM_ADD] THEN
 253   ONCE_REWRITE_TAC[REAL_ARITH `x + &1 <= y <=> &1 <= y - x * &1`] THEN
 254   REWRITE_TAC[GSYM REAL_SUB_LDISTRIB] THEN
 255   SUBGOAL_THEN `?N. &1 <= &N * (c / abs x - &1)` STRIP_ASSUME_TAC THENL
 256    [ALL_TAC;
 257     EXISTS_TAC `N:num` THEN REWRITE_TAC[GE] THEN REPEAT STRIP_TAC THEN
 258     FIRST_ASSUM(MATCH_MP_TAC o MATCH_MP (REAL_ARITH
 259      `&1 <= x ==> x <= y ==> &1 <= y`)) THEN
 260     MATCH_MP_TAC REAL_LE_RMUL THEN
 261     ASM_SIMP_TAC[REAL_ARITH `a <= b ==> a <= b + &1`;
 262                  REAL_OF_NUM_LE; REAL_LE_RADD] THEN
 263     REWRITE_TAC[REAL_LE_SUB_LADD; REAL_ADD_LID] THEN
 264     ASM_SIMP_TAC[REAL_LE_RDIV_EQ; GSYM REAL_ABS_NZ; REAL_MUL_LID;
 265                  REAL_LT_IMP_LE]] THEN
 266   ASM_SIMP_TAC[GSYM REAL_LE_LDIV_EQ; REAL_LT_SUB_LADD; REAL_ADD_LID;
 267                REAL_LT_RDIV_EQ; GSYM REAL_ABS_NZ; REAL_MUL_LID;
 268                REAL_ARCH_SIMPLE]);;
 269
 270 let REAL_ATN_POWSER_DIFFL = prove
 271  (`!x. abs(x) < &1
 272        ==> ((\x. suminf (\n. (if EVEN n then &0
 273                               else --(&1) pow ((n - 1) DIV 2) / &n) * x pow n))
 274             diffl (inv(&1 + x pow 2))) x`,
 275   REPEAT STRIP_TAC THEN
 276   FIRST_ASSUM(MP_TAC o MATCH_MP REAL_ATN_POWSER_DIFFS_SUM) THEN
 277   DISCH_THEN(SUBST1_TAC o MATCH_MP SUM_UNIQ) THEN
 278   MATCH_MP_TAC TERMDIFF THEN
 279   SUBGOAL_THEN `?K. abs(x) < abs(K) /\ abs(K) < &1` STRIP_ASSUME_TAC THENL
 280    [EXISTS_TAC `(&1 + abs(x)) / &2` THEN
 281     SIMP_TAC[REAL_LT_LDIV_EQ; REAL_ABS_DIV; REAL_ABS_NUM;
 282              REAL_LT_RDIV_EQ; REAL_OF_NUM_LT; ARITH] THEN
 283     UNDISCH_TAC `abs(x) < &1` THEN REAL_ARITH_TAC; ALL_TAC] THEN
 284   EXISTS_TAC `K:real` THEN ASM_REWRITE_TAC[] THEN
 285   ASM_SIMP_TAC[REAL_ATN_POWSER_SUMMABLE; REAL_ATN_POWSER_DIFFS_SUMMABLE;
 286                REAL_ATN_POWSER_DIFFS_DIFFS_SUMMABLE]);;
 287
 288 let REAL_ATN_POWSER = prove
 289  (`!x. abs(x) < &1
 290        ==> (\n. (if EVEN n then &0
 291                  else --(&1) pow ((n - 1) DIV 2) / &n) * x pow n)
 292            sums (atn x)`,
 293   REPEAT STRIP_TAC THEN
 294   FIRST_ASSUM(MP_TAC o MATCH_MP REAL_ATN_POWSER_SUMMABLE) THEN
 295   DISCH_THEN(MP_TAC o MATCH_MP SUMMABLE_SUM) THEN
 296   SUBGOAL_THEN
 297    `suminf (\n. (if EVEN n then &0
 298                  else --(&1) pow ((n - 1) DIV 2) / &n) * x pow n) = atn(x)`
 299    (fun th -> REWRITE_TAC[th]) THEN
 300   ONCE_REWRITE_TAC[REAL_ARITH `(a = b) <=> (a - b = &0)`] THEN
 301   SUBGOAL_THEN
 302    `suminf (\n. (if EVEN n then &0
 303                  else --(&1) pow ((n - 1) DIV 2) / &n) * &0 pow n) -
 304     atn(&0) = &0`
 305   MP_TAC THENL
 306    [MATCH_MP_TAC(REAL_ARITH `(a = &0) /\ (b = &0) ==> (a - b = &0)`) THEN
 307     CONJ_TAC THENL
 308      [CONV_TAC SYM_CONV THEN MATCH_MP_TAC SUM_UNIQ THEN
 309       MP_TAC(SPEC `&0` GP) THEN
 310       REWRITE_TAC[REAL_ABS_NUM; REAL_OF_NUM_LT; ARITH] THEN
 311       DISCH_THEN(MP_TAC o SPEC `&0` o MATCH_MP SER_CMUL) THEN
 312       REWRITE_TAC[REAL_MUL_LZERO] THEN
 313       MATCH_MP_TAC(TAUT `(a = b) ==> a ==> b`) THEN
 314       AP_THM_TAC THEN AP_TERM_TAC THEN ABS_TAC THEN
 315       COND_CASES_TAC THEN ASM_REWRITE_TAC[REAL_MUL_LZERO] THEN
 316       CONV_TAC SYM_CONV THEN
 317       REWRITE_TAC[REAL_ENTIRE; REAL_POW_EQ_0] THEN ASM_MESON_TAC[EVEN];
 318       GEN_REWRITE_TAC (LAND_CONV o RAND_CONV) [GSYM TAN_0] THEN
 319       MATCH_MP_TAC TAN_ATN THEN
 320       SIMP_TAC[PI2_BOUNDS; REAL_ARITH `&0 < x ==> --x < &0`]];
 321     ALL_TAC] THEN
 322   ASM_CASES_TAC `x = &0` THEN ASM_REWRITE_TAC[] THEN
 323   DISCH_THEN(fun th -> GEN_REWRITE_TAC RAND_CONV [SYM th]) THEN
 324   MP_TAC(SPEC `\x. suminf (\n. (if EVEN n then &0
 325
 326                        else --(&1) pow ((n - 1) DIV 2) / &n) * x pow n) -
 327           atn x` DIFF_ISCONST_END_SIMPLE) THEN
 328   FIRST_ASSUM(DISJ_CASES_TAC o MATCH_MP (REAL_ARITH
 329     `~(x = &0) ==> &0 < x \/ x < &0`))
 330   THENL
 331    [DISCH_THEN(MP_TAC o SPECL [`&0`; `x:real`]);
 332     CONV_TAC(RAND_CONV SYM_CONV) THEN
 333     DISCH_THEN(MP_TAC o SPECL [`x:real`; `&0`])] THEN
 334   (REWRITE_TAC[] THEN DISCH_THEN MATCH_MP_TAC THEN
 335    ASM_REWRITE_TAC[] THEN
 336    X_GEN_TAC `u:real` THEN REPEAT STRIP_TAC THEN
 337    SUBGOAL_THEN `abs(u) < &1` (MP_TAC o MATCH_MP REAL_ATN_POWSER_DIFFL) THENL
 338     [POP_ASSUM_LIST(MP_TAC o end_itlist CONJ) THEN REAL_ARITH_TAC;
 339      ALL_TAC] THEN
 340    DISCH_THEN(MP_TAC o C CONJ (SPEC `u:real` DIFF_ATN)) THEN
 341    DISCH_THEN(MP_TAC o MATCH_MP DIFF_SUB) THEN
 342    REWRITE_TAC[REAL_SUB_REFL]));;
 343
 344 (* ------------------------------------------------------------------------- *)
 345 (* A more Taylor-like version with a simply bounded remainder term.          *)
 346 (* ------------------------------------------------------------------------- *)
 347
 348 let MCLAURIN_ATN_SIMPLE = prove
 349  (`!x n k. abs(x) <= inv(&2 pow k) /\ ~(k = 0)
 350            ==> abs(atn x -
 351                    sum(0,n) (\m. (if EVEN m then &0
 352                                   else --(&1) pow ((m - 1) DIV 2) / &m) *
 353                                   x pow m))
 354                <= &2 * abs(x) pow n`,
 355   REPEAT STRIP_TAC THEN SUBGOAL_THEN `abs(x) < &1` ASSUME_TAC THENL
 356    [MATCH_MP_TAC REAL_LET_TRANS THEN EXISTS_TAC `inv(&2 pow k)` THEN
 357     ASM_REWRITE_TAC[REAL_ARITH `a < &1 <=> &0 < &1 - a`; REAL_POW2_CLAUSES];
 358     ALL_TAC] THEN
 359   FIRST_ASSUM(MP_TAC o MATCH_MP REAL_ATN_POWSER) THEN
 360   DISCH_THEN(fun th -> ASSUME_TAC(SYM(MATCH_MP SUM_UNIQ th)) THEN
 361                        MP_TAC(MATCH_MP SUM_SUMMABLE th)) THEN
 362   DISCH_THEN(MP_TAC o MATCH_MP SER_OFFSET) THEN
 363   DISCH_THEN(MP_TAC o SPEC `n:num`) THEN ASM_REWRITE_TAC[] THEN
 364   DISCH_THEN(MP_TAC o MATCH_MP SUM_UNIQ) THEN
 365   MATCH_MP_TAC(REAL_ARITH
 366    `abs(r) <= e ==> (f - s = r) ==> abs(f - s) <= e`) THEN
 367   SUBGOAL_THEN
 368    `(\m. abs(x) pow (m + n)) sums (abs(x) pow n) * inv(&1 - abs(x))`
 369   ASSUME_TAC THENL
 370    [FIRST_ASSUM(MP_TAC o MATCH_MP GP o MATCH_MP (REAL_ARITH
 371       `abs(x) < &1 ==> abs(abs x) < &1`)) THEN
 372     DISCH_THEN(MP_TAC o SPEC `abs(x) pow n` o MATCH_MP SER_CMUL) THEN
 373     ONCE_REWRITE_TAC[ADD_SYM] THEN REWRITE_TAC[GSYM REAL_POW_ADD];
 374     ALL_TAC] THEN
 375   MATCH_MP_TAC REAL_LE_TRANS THEN
 376   EXISTS_TAC `suminf (\m. abs(x) pow (m + n))` THEN CONJ_TAC THENL
 377    [ALL_TAC;
 378     FIRST_ASSUM(SUBST1_TAC o SYM o MATCH_MP SUM_UNIQ) THEN
 379     GEN_REWRITE_TAC RAND_CONV [REAL_MUL_SYM] THEN
 380     MATCH_MP_TAC REAL_LE_LMUL THEN SIMP_TAC[REAL_ABS_POS; REAL_POW_LE] THEN
 381     GEN_REWRITE_TAC RAND_CONV [GSYM REAL_INV_INV] THEN
 382     MATCH_MP_TAC REAL_LE_INV2 THEN
 383     ONCE_REWRITE_TAC[REAL_ARITH `a <= &1 - b <=> b <= &1 - a`] THEN
 384     CONV_TAC REAL_RAT_REDUCE_CONV THEN
 385     MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `inv(&2 pow k)` THEN
 386     ASM_REWRITE_TAC[] THEN REWRITE_TAC[real_div; REAL_MUL_LID] THEN
 387     MATCH_MP_TAC REAL_LE_INV2 THEN CONV_TAC REAL_RAT_REDUCE_CONV THEN
 388     GEN_REWRITE_TAC LAND_CONV [GSYM REAL_POW_1] THEN
 389     ASM_SIMP_TAC[REAL_POW_MONO; REAL_OF_NUM_LE; ARITH;
 390                  ARITH_RULE `~(k = 0) ==> 1 <= k`]] THEN
 391   SUBGOAL_THEN
 392    `!m. abs((if EVEN (m + n) then &0
 393              else --(&1) pow (((m + n) - 1) DIV 2) / &(m + n)) *
 394              x pow (m + n))
 395         <= abs(x) pow (m + n)`
 396   ASSUME_TAC THENL
 397    [GEN_TAC THEN COND_CASES_TAC THEN
 398     SIMP_TAC[REAL_MUL_LZERO; REAL_ABS_NUM; REAL_POW_LE; REAL_ABS_POS] THEN
 399     REWRITE_TAC[REAL_ABS_MUL; REAL_ABS_DIV; REAL_ABS_POW; REAL_ABS_NEG] THEN
 400     REWRITE_TAC[REAL_ABS_NUM; REAL_POW_ONE; REAL_MUL_LID] THEN
 401     GEN_REWRITE_TAC RAND_CONV [GSYM REAL_MUL_LID] THEN
 402     MATCH_MP_TAC REAL_LE_RMUL THEN SIMP_TAC[REAL_POW_LE; REAL_ABS_POS] THEN
 403     REWRITE_TAC[real_div; REAL_MUL_LID] THEN
 404     GEN_REWRITE_TAC RAND_CONV [GSYM REAL_INV_1] THEN
 405     MATCH_MP_TAC REAL_LE_INV2 THEN CONV_TAC REAL_RAT_REDUCE_CONV THEN
 406     REWRITE_TAC[REAL_OF_NUM_LE; ARITH_RULE `1 <= n <=> ~(n = 0)`] THEN
 407     ASM_MESON_TAC[EVEN]; ALL_TAC] THEN
 408   MATCH_MP_TAC REAL_LE_TRANS THEN
 409   EXISTS_TAC
 410    `suminf (\m. abs((if EVEN (m + n) then &0
 411                      else --(&1) pow (((m + n) - 1) DIV 2) / &(m + n)) *
 412                     x pow (m + n)))` THEN
 413   CONJ_TAC THENL
 414    [MATCH_MP_TAC SER_ABS THEN MATCH_MP_TAC SER_COMPARA THEN
 415     EXISTS_TAC `\m. abs(x) pow (m + n)` THEN
 416     ASM_REWRITE_TAC[] THEN ASM_MESON_TAC[SUM_SUMMABLE]; ALL_TAC] THEN
 417   MATCH_MP_TAC SER_LE THEN ASM_REWRITE_TAC[] THEN CONJ_TAC THENL
 418    [MATCH_MP_TAC SER_COMPARA THEN
 419     EXISTS_TAC `\m. abs(x) pow (m + n)` THEN
 420     ASM_REWRITE_TAC[]; ALL_TAC] THEN
 421   ASM_MESON_TAC[SUM_SUMMABLE]);;
 422
 423 let MCLAURIN_ATN_APPROX = prove
 424  (`!x n k. abs(x) <= inv(&2 pow k) /\ ~(k = 0)
 425            ==> abs(atn x -
 426                    sum(0,n) (\m. (if EVEN m then &0
 427                                   else --(&1) pow ((m - 1) DIV 2) / &m) *
 428                                   x pow m))
 429                <= inv(&2 pow (n * k - 1))`,
 430   REPEAT STRIP_TAC THEN ASM_CASES_TAC `n = 0` THENL
 431    [ASM_REWRITE_TAC[sum; REAL_SUB_RZERO; MULT_CLAUSES; SUB_0] THEN
 432     MP_TAC(SPECL [`x:real`; `2`; `k:num`] MCLAURIN_ATN_SIMPLE) THEN
 433     ASM_REWRITE_TAC[] THEN
 434     MATCH_MP_TAC(REAL_ARITH
 435      `abs(y) + d <= e ==> abs(x - y) <= d ==> abs(x) <= e`) THEN
 436     CONV_TAC(ONCE_DEPTH_CONV REAL_SUM_CONV) THEN
 437     REWRITE_TAC[real_pow; REAL_POW_1] THEN CONV_TAC NUM_REDUCE_CONV THEN
 438     CONV_TAC REAL_RAT_REDUCE_CONV THEN
 439     REWRITE_TAC[real_div; REAL_MUL_LID; REAL_INV_1; REAL_ADD_LID] THEN
 440     SUBGOAL_THEN `abs(x) <= inv(&2)` ASSUME_TAC THENL
 441      [MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `inv(&2 pow k)` THEN
 442       ASM_REWRITE_TAC[] THEN MATCH_MP_TAC REAL_LE_INV2 THEN
 443       CONV_TAC REAL_RAT_REDUCE_CONV THEN
 444       CONV_TAC REAL_RAT_REDUCE_CONV THEN
 445       GEN_REWRITE_TAC LAND_CONV [GSYM REAL_POW_1] THEN
 446       ASM_SIMP_TAC[REAL_POW_MONO; REAL_OF_NUM_LE; ARITH;
 447                    ARITH_RULE `~(k = 0) ==> 1 <= k`];
 448       ALL_TAC] THEN
 449     MATCH_MP_TAC REAL_LE_TRANS THEN
 450     EXISTS_TAC `inv(&2) + &2 * inv(&2) pow 2` THEN
 451     CONV_TAC(RAND_CONV REAL_RAT_REDUCE_CONV) THEN REWRITE_TAC[REAL_POW_1] THEN
 452     MATCH_MP_TAC REAL_LE_ADD2 THEN ASM_REWRITE_TAC[] THEN
 453     ASM_SIMP_TAC[REAL_LE_LMUL_EQ; REAL_OF_NUM_LT; ARITH;
 454                  REAL_POW_LE2; REAL_OF_NUM_LE; REAL_ABS_POS];
 455     ALL_TAC] THEN
 456   MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `&2 * abs(x) pow n` THEN
 457   CONJ_TAC THENL
 458    [MATCH_MP_TAC MCLAURIN_ATN_SIMPLE THEN ASM_MESON_TAC[]; ALL_TAC] THEN
 459   ASM_SIMP_TAC[REAL_POW_SUB; REAL_OF_NUM_EQ; ARITH_EQ;
 460                ARITH_RULE `1 <= x <=> ~(x = 0)`; MULT_EQ_0] THEN
 461   REWRITE_TAC[REAL_INV_DIV; REAL_POW_1] THEN REWRITE_TAC[real_div] THEN
 462   SIMP_TAC[REAL_LE_LMUL_EQ; REAL_OF_NUM_LT; ARITH] THEN
 463   ONCE_REWRITE_TAC[MULT_SYM] THEN REWRITE_TAC[GSYM REAL_POW_POW] THEN
 464   ONCE_REWRITE_TAC[GSYM REAL_POW_INV] THEN MATCH_MP_TAC REAL_POW_LE2 THEN
 465   ASM_REWRITE_TAC[REAL_ABS_POS]);;
 466
 467 (* ------------------------------------------------------------------------- *)
 468 (* Rules to return approximations to atn(x) good to 2^-p given |x| <= 2^-k.  *)
 469 (* ------------------------------------------------------------------------- *)
 470
 471 let mclaurin_atn_rule,MCLAURIN_ATN_RULE =
 472   let x_tm = `x:real`
 473   and n_tm = `n:num`
 474   and k_tm = `k:num`
 475   and inv_tm = `inv`
 476   and le_tm = `(<=):real->real->bool`
 477   and pow2_tm = `(pow) (&2)` in
 478   let pth = SPECL [x_tm; n_tm; k_tm] MCLAURIN_ATN_APPROX
 479   and CLEAN_RULE = REWRITE_RULE[real_pow]
 480   and MATCH_REAL_LE_TRANS = MATCH_MP REAL_LE_TRANS
 481   and num_0 = Int 0
 482   and num_1 = Int 1 in
 483   let mclaurin_atn_rule k0 p0 =
 484     if k0 = 0 then failwith "mclaurin_atn_rule: must have |x| <= 1/2" else
 485     let k = Int k0
 486     and p = Int p0 in
 487     let n = Num.int_of_num(ceiling_num ((p +/ k) // k)) in
 488     let ns = if n mod 2 = 0 then 0--(n - 1) else 0--(n - 2) in
 489     map (fun m -> if m mod 2 = 0 then num_0 else
 490                   (if (m - 1) mod 4 = 0 then I else minus_num)
 491                   (num_1 // Int m)) ns
 492   and MCLAURIN_ATN_RULE k0 p0 =
 493     if k0 = 0 then failwith "MCLAURIN_ATN_RULE: must have |x| <= 1/2" else
 494     let k = Int k0
 495     and p = Int p0 in
 496     let n = ceiling_num ((p +/ k) // k) in
 497     let th1 = INST [mk_numeral k,k_tm; mk_numeral n,n_tm] pth in
 498     let th2 = ASSUME (lhand(lhand(concl th1)))
 499     and th3 = EQF_ELIM(NUM_REDUCE_CONV(rand(rand(lhand(concl th1))))) in
 500     let th4 = MP th1 (CONJ th2 th3) in
 501     let th5 = CONV_RULE(ONCE_DEPTH_CONV REAL_HORNER_SUM_CONV) th4 in
 502     let th6 = CLEAN_RULE th5 in
 503     let th7 = CONV_RULE (NUM_REDUCE_CONV THENC LAND_CONV REAL_RAT_REDUCE_CONV)
 504                         (BETA_RULE th6) in
 505     let tm1 = mk_comb(inv_tm,mk_comb(pow2_tm,mk_numeral p)) in
 506     let tm2 = mk_comb(mk_comb(le_tm,rand(concl th7)),tm1) in
 507     let th8 = EQT_ELIM((NUM_REDUCE_CONV THENC REAL_RAT_REDUCE_CONV) tm2) in
 508     let th9 = MATCH_REAL_LE_TRANS (CONJ th7 th8) in
 509     GEN x_tm (DISCH_ALL th9) in
 510   mclaurin_atn_rule,MCLAURIN_ATN_RULE;;
 511
 512 (* ------------------------------------------------------------------------- *)
 513 (* Lemmas for Machin-type formulas.                                          *)
 514 (* ------------------------------------------------------------------------- *)
 515
 516 let TAN_ADD_ATN_SIDECOND = prove
 517  (`!x y. ~(x * y = &1) ==> ~(cos(atn x + atn y) = &0)`,
 518   REPEAT GEN_TAC THEN DISCH_TAC THEN
 519   REWRITE_TAC[COS_ADD; REAL_ARITH `(a - b = &0) <=> (a = b)`] THEN
 520   DISCH_THEN(MP_TAC o AP_TERM `(*) (inv(cos(atn x)))`) THEN
 521   SIMP_TAC[REAL_MUL_ASSOC; REAL_MUL_LINV; COS_ATN_NZ; REAL_MUL_LID] THEN
 522   DISCH_THEN(MP_TAC o AP_TERM `(*) (inv(cos(atn y)))`) THEN
 523   SIMP_TAC[REAL_MUL_LINV; COS_ATN_NZ; REAL_MUL_LID; GSYM REAL_MUL_ASSOC] THEN
 524   ONCE_REWRITE_TAC[AC REAL_MUL_AC `a * b * c * d = (c * b) * (d * a)`] THEN
 525   ASM_REWRITE_TAC[GSYM tan; GSYM real_div; ATN_TAN]);;
 526
 527 let ATN_ADD = prove
 528  (`!x y. ~(x * y = &1) /\
 529          abs(atn x + atn y) < pi / &2
 530          ==> (atn(x) + atn(y) = atn((x + y) / (&1 - x * y)))`,
 531   REPEAT STRIP_TAC THEN
 532   SUBGOAL_THEN `tan(atn(x) + atn(y)) = (x + y) / (&1 - x * y)` MP_TAC THENL
 533    [ASM_SIMP_TAC[ATN_TAN; TAN_ADD; COS_ATN_NZ; TAN_ADD_ATN_SIDECOND];
 534     DISCH_THEN(SUBST1_TAC o SYM) THEN CONV_TAC SYM_CONV THEN
 535     ASM_SIMP_TAC[TAN_ATN; REAL_ARITH `abs(x) < e ==> --e < x /\ x < e`]]);;
 536
 537 let ATN_ADD_SMALL_LEMMA_POS = prove
 538  (`!x y. &0 < y /\ x * y < &1
 539          ==> atn(x) + atn(y) < pi / &2`,
 540   REPEAT STRIP_TAC THEN REWRITE_TAC[GSYM REAL_LT_SUB_LADD] THEN
 541   SUBGOAL_THEN `pi / &2 - atn y = atn(tan(pi / &2 - atn y))` SUBST1_TAC THENL
 542    [CONV_TAC SYM_CONV THEN MATCH_MP_TAC TAN_ATN THEN
 543     MATCH_MP_TAC(REAL_ARITH
 544      `&0 < x /\ x < a ==> --a < a - x /\ a - x < a`) THEN
 545     REWRITE_TAC[ATN_BOUNDS] THEN
 546     GEN_REWRITE_TAC LAND_CONV [GSYM ATN_0] THEN
 547     ASM_SIMP_TAC[ATN_MONO_LT];
 548     MATCH_MP_TAC ATN_MONO_LT THEN REWRITE_TAC[TAN_COT; ATN_TAN] THEN
 549     GEN_REWRITE_TAC RAND_CONV [GSYM REAL_MUL_LID] THEN
 550     ASM_SIMP_TAC[GSYM real_div; REAL_LT_RDIV_EQ; REAL_LT_IMP_NZ]]);;
 551
 552 let ATN_ADD_SMALL_LEMMA = prove
 553  (`!x y. abs(x * y) < &1 ==> abs(atn(x) + atn(y)) < pi / &2`,
 554   REPEAT STRIP_TAC THEN
 555   MATCH_MP_TAC(REAL_ARITH
 556    `--a < x /\ x < a /\ --a < y /\ y < a /\
 557     (&0 < y ==> x + y < a) /\
 558     (&0 < --y ==> --x + --y < a)
 559     ==> abs(x + y) < a`) THEN
 560   REWRITE_TAC[ATN_BOUNDS] THEN CONJ_TAC THEN
 561   REPEAT STRIP_TAC THEN REWRITE_TAC[GSYM ATN_NEG] THEN
 562   MATCH_MP_TAC ATN_ADD_SMALL_LEMMA_POS THEN
 563   ASM_SIMP_TAC[REAL_ARITH `abs(x) < &1 ==> x < &1`;
 564                REAL_ARITH `--x * -- y = x * y`] THEN
 565   FIRST_ASSUM(MATCH_MP_TAC o MATCH_MP (REAL_ARITH
 566    `&0 < y ==> (z <= &0 ==> y <= &0) ==> &0 < z`)) THEN
 567   MATCH_MP_TAC(REAL_ARITH
 568    `(y < &0 ==> z < &0) /\ ((y = &0) ==> (z = &0))
 569     ==> y <= &0 ==> z <= &0`) THEN
 570   SIMP_TAC[ATN_0; GSYM ATN_NEG] THEN
 571   GEN_REWRITE_TAC (RAND_CONV o RAND_CONV) [GSYM ATN_0] THEN
 572   SIMP_TAC[ATN_MONO_LT]);;
 573
 574 let ATN_ADD_SMALL = prove
 575  (`!x y. abs(x * y) < &1
 576          ==> (atn(x) + atn(y) = atn((x + y) / (&1 - x * y)))`,
 577   REPEAT STRIP_TAC THEN MATCH_MP_TAC ATN_ADD THEN
 578   ASM_SIMP_TAC[ATN_ADD_SMALL_LEMMA] THEN
 579   POP_ASSUM MP_TAC THEN REAL_ARITH_TAC);;
 580
 581 let ATN_ADD_CONV =
 582   let match_fn = PART_MATCH (lhand o rand) ATN_ADD_SMALL in
 583   let overall_fn =
 584     C MP TRUTH o
 585     CONV_RULE
 586        (COMB2_CONV REAL_RAT_REDUCE_CONV
 587          (RAND_CONV REAL_RAT_REDUCE_CONV)) o
 588     match_fn in
 589   fun tm -> if is_ratconst(rand(rand tm)) &
 590                is_ratconst(rand(lhand tm))
 591             then overall_fn tm
 592             else failwith "ATN_ADD_CONV: Atn of nonconstant";;
 593
 594 let ATN_CMUL_CONV =
 595   let pth_base = prove
 596    (`(&0 * atn(x) = &0) /\
 597      (&1 * atn(x) = atn(x))`,
 598     REWRITE_TAC[REAL_MUL_LZERO; REAL_MUL_LID])
 599   and pth_0,pth_1 = (CONJ_PAIR o prove)
 600    (`(&(NUMERAL(BIT0 n)) * atn(x) =
 601       &(NUMERAL n) * atn(x) + &(NUMERAL n) * atn(x)) /\
 602      (&(NUMERAL(BIT1 n)) * atn(x) =
 603       atn(x) + &(NUMERAL n) * atn(x) + &(NUMERAL n) * atn(x))`,
 604     REWRITE_TAC[REAL_MUL_LZERO; REAL_MUL_LID] THEN
 605     REWRITE_TAC[NUMERAL; BIT0; BIT1] THEN
 606     REWRITE_TAC[GSYM REAL_OF_NUM_ADD; GSYM REAL_OF_NUM_SUC] THEN
 607     REAL_ARITH_TAC) in
 608   let rewr_base = GEN_REWRITE_CONV I [pth_base]
 609   and rewr_0 = GEN_REWRITE_CONV I [pth_0]
 610   and rewr_1 = GEN_REWRITE_CONV I [pth_1] in
 611   let rec ATN_CMUL_CONV tm =
 612     if not (is_ratconst(rand(rand tm))) then failwith "ATN_CMUL_CONV" else
 613     try rewr_base tm with Failure _ ->
 614     try let th1 = rewr_0 tm in
 615         let tm1 = rand(concl th1) in
 616         let th2 = ATN_CMUL_CONV(rand tm1) in
 617         let th3 = MK_COMB(AP_TERM (rator(rator tm1)) th2,th2) in
 618         let th4 = TRANS th3 (ATN_ADD_CONV(rand(concl th3))) in
 619         TRANS th1 th4
 620     with Failure _ ->
 621         let th1 = rewr_1 tm in
 622         let tm1 = rand(rand(concl th1)) in
 623         let th2 = ATN_CMUL_CONV(rand tm1) in
 624         let th3 = MK_COMB(AP_TERM (rator(rator tm1)) th2,th2) in
 625         let th4 = TRANS th3 (ATN_ADD_CONV(rand(concl th3))) in
 626         let th5 = AP_TERM (rator(rand(concl th1))) th4 in
 627         let th6 = TRANS th5 (ATN_ADD_CONV(rand(concl th5))) in
 628         TRANS th1 th6 in
 629   ATN_CMUL_CONV;;
 630
 631 let ATN_SUB_CONV =
 632   let pth = prove
 633    (`(atn(x) - atn(y) = atn(x) + atn(--y))`,
 634     REWRITE_TAC[real_sub; ATN_NEG]) in
 635   GEN_REWRITE_CONV I [pth] THENC
 636   RAND_CONV(RAND_CONV REAL_RAT_NEG_CONV) THENC
 637   ATN_ADD_CONV;;
 638
 639 let MACHIN_CONV =
 640   DEPTH_CONV(ATN_ADD_CONV ORELSEC ATN_SUB_CONV ORELSEC ATN_CMUL_CONV);;
 641
 642 let MACHIN_RULE tm = SYM(TRANS (MACHIN_CONV tm) ATN_1);;
 643
 644 let MACHIN_1 = time MACHIN_RULE `&4 * atn(&1 / &5) - atn(&1 / &239)`;;
 645 let MACHIN_2 = time MACHIN_RULE `atn(&1 / &2) + atn(&1 / &3)`;;
 646 let MACHIN_3 = time MACHIN_RULE `&2 * atn(&1 / &2) - atn(&1 / &7)`;;
 647 let MACHIN_4 = time MACHIN_RULE `&2 * atn(&1 / &3) + atn(&1 / &7)`;;
 648
 649 let EULER = time MACHIN_RULE `&5 * atn(&1 / &7) + &2 * atn (&3 / &79)`;;
 650
 651 let GAUSS_MACHIN = time MACHIN_RULE
 652   `&12 * atn(&1 / &18) + &8 * atn (&1 / &57) - &5 * atn(&1 / &239)`;;
 653
 654 let STRASSNITZKY_MACHIN = time MACHIN_RULE
 655   `atn(&1 / &2) + atn (&1 / &5) + atn(&1 / &8)`;;
 656
 657 let MACHINLIKE_1 = time MACHIN_RULE
 658   `&6 * atn(&1 / &8) + &2 * atn(&1 / &57) + atn(&1 / &239)`;;
 659 let MACHINLIKE_2 = time MACHIN_RULE
 660   `&4 * atn(&1 / &5) - &1 * atn(&1 / &70) + atn(&1 / &99)`;;
 661 let MACHINLIKE_3 = time MACHIN_RULE
 662   `&1 * atn(&1 / &2) + &1 * atn(&1 / &5) + atn(&1 / &8)`;;
 663 let MACHINLIKE_4 = time MACHIN_RULE
 664   `&8 * atn(&1 / &10) - &1 * atn(&1 / &239) - &4 * atn(&1 / &515)`;;
 665 let MACHINLIKE_5 = time MACHIN_RULE
 666   `&5 * atn(&1 / &7) + &4 * atn(&1 / &53) + &2 * atn(&1 / &4443)`;;
 667
 668 (***** Hopefully this one would work, but it takes a long time
 669
 670 let HWANG_MACHIN = time MACHIN_RULE
 671   `&183 * atn(&1 / &239) + &32 * atn(&1 / &1023) -
 672    &68 * atn(&1 / &5832) + &12 * atn(&1 / &110443) -
 673    &12 * atn(&1 / &4841182) - &100 * atn(&1 / &6826318)`;;
 674
 675  *****)
 676
 677 (* ------------------------------------------------------------------------- *)
 678 (* Approximate the arctan of a rational number.                              *)
 679 (* ------------------------------------------------------------------------- *)
 680
 681 let rec POLY l x =
 682   if l = [] then num_0
 683   else hd l +/ (x */ POLY (tl l) x);;
 684
 685 let atn_approx_conv,ATN_APPROX_CONV =
 686   let atn_tm = `atn`
 687   and num_0 = Int 0
 688   and num_1 = Int 1
 689   and num_2 = Int 2 in
 690   let rec log_2 x = if x <=/ num_1 then log_2 (num_2 */ x) -/ num_1
 691        else if x >/ num_2 then log_2 (x // num_2) +/ num_1
 692                       else num_1 in
 693   let pth = prove
 694    (`!p. abs(atn(&0) - &0) <= inv(&2 pow p)`,
 695     SIMP_TAC[ATN_0; REAL_SUB_REFL; REAL_ABS_NUM; REAL_LE_INV_EQ;
 696              REAL_POW_LE; REAL_POS]) in
 697   let atn_approx_conv p r =
 698     if r =/ num_0 then num_0 else
 699     let k = Num.int_of_num(minus_num(log_2(abs_num r))) in
 700     if k < 1 then failwith "atn_approx_conv: argument too big" else
 701     let rats = mclaurin_atn_rule k p in
 702     POLY rats r
 703   and ATN_APPROX_CONV p tm =
 704     let atm,rtm = dest_comb tm in
 705     if atm <> atn_tm then failwith "ATN_APPROX_CONV" else
 706     let r = rat_of_term rtm in
 707     if r =/ num_0 then SPEC (mk_small_numeral p) pth else
 708     let k = Num.int_of_num(minus_num(log_2(abs_num r))) in
 709     if k < 1 then failwith "ATN_APPROX_CONV: argument too big" else
 710     let th1 = MCLAURIN_ATN_RULE k p in
 711     let th2 = SPEC rtm th1 in
 712     let th3 = MP th2 (EQT_ELIM(REAL_RAT_REDUCE_CONV(lhand(concl th2)))) in
 713     CONV_RULE(LAND_CONV(RAND_CONV(RAND_CONV REAL_RAT_REDUCE_CONV))) th3 in
 714   atn_approx_conv,ATN_APPROX_CONV;;
 715
 716 (* ------------------------------------------------------------------------- *)
 717 (* Approximate pi using this and a Machin-type formula.                      *)
 718 (* ------------------------------------------------------------------------- *)
 719
 720 let pi_approx_rule,PI_APPROX_RULE =
 721   let const_1_8 = Int 1 // Int 8
 722   and const_1_57 = Int 1 // Int 57
 723   and const_1_239 = Int 1 // Int 239
 724   and const_24 = Int 24
 725   and const_8 = Int 8
 726   and const_4 = Int 4
 727   and tm_1_8 = `atn(&1 / &8)`
 728   and tm_1_57 = `atn(&1 / &57)`
 729   and tm_1_239 = `atn(&1 / &239)`
 730   and q1_tm = `q1:num`
 731   and q2_tm = `q2:num`
 732   and p_tm = `p:num` in
 733   let pth = prove
 734    (`(q1 = p + 5) /\
 735      (q2 = p + 6) /\
 736      abs(atn(&1 / &8) - a1) <= inv(&2 pow q1) /\
 737      abs(atn(&1 / &57) - a2) <= inv(&2 pow q2) /\
 738      abs(atn(&1 / &239) - a3) <= inv(&2 pow q2)
 739      ==> abs(pi - (&24 * a1 + &8 * a2 + &4 * a3)) <= inv(&2 pow p)`,
 740     DISCH_THEN(CONJUNCTS_THEN2 ASSUME_TAC MP_TAC) THEN
 741     DISCH_THEN(CONJUNCTS_THEN2 ASSUME_TAC MP_TAC) THEN
 742     ASM_REWRITE_TAC[] THEN POP_ASSUM_LIST(K ALL_TAC) THEN STRIP_TAC THEN
 743     MATCH_MP_TAC REAL_LE_LCANCEL_IMP THEN EXISTS_TAC `abs(inv(&2 pow 2))` THEN
 744     SIMP_TAC[REAL_POW2_CLAUSES; REAL_ARITH `&0 < x ==> &0 < abs(x)`] THEN
 745     REWRITE_TAC[GSYM REAL_ABS_MUL] THEN
 746     REWRITE_TAC[REAL_ABS_INV; REAL_ABS_POW; REAL_ABS_NUM] THEN
 747     REWRITE_TAC[GSYM REAL_INV_MUL; GSYM REAL_POW_ADD] THEN
 748     CONV_TAC REAL_RAT_REDUCE_CONV THEN
 749     REWRITE_TAC[REAL_SUB_LDISTRIB; REAL_ADD_LDISTRIB; REAL_MUL_ASSOC] THEN
 750     CONV_TAC REAL_RAT_REDUCE_CONV THEN
 751     REWRITE_TAC[real_div; REAL_MUL_LID] THEN
 752     GEN_REWRITE_TAC (LAND_CONV o RAND_CONV o LAND_CONV) [REAL_MUL_SYM] THEN
 753     REWRITE_TAC[GSYM real_div; MACHINLIKE_1] THEN
 754     REWRITE_TAC[REAL_ARITH `(x1 + x2 + x3) - (y1 + y2 + y3) =
 755                             (x1 - y1) + (x2 - y2) + (x3 - y3)`] THEN
 756     REWRITE_TAC[GSYM REAL_SUB_LDISTRIB] THEN BOUND_SUMPROD_TAC THEN
 757     GEN_REWRITE_TAC (LAND_CONV o ONCE_DEPTH_CONV) [ADD_SYM] THEN
 758     REWRITE_TAC[REAL_POW_ADD; REAL_INV_MUL; REAL_MUL_ASSOC] THEN
 759     SIMP_TAC[GSYM REAL_ADD_RDISTRIB; REAL_LE_RMUL_EQ; REAL_POW2_CLAUSES] THEN
 760     CONV_TAC REAL_RAT_REDUCE_CONV) in
 761   let pi_approx_rule p =
 762     let q1 = p + 5
 763     and q2 = p + 6 in
 764     let a1 = atn_approx_conv q1 const_1_8
 765     and a2 = atn_approx_conv q2 const_1_57
 766     and a3 = atn_approx_conv q2 const_1_239 in
 767     const_24 */ a1 +/ const_8 */ a2 +/ const_4 */ a3
 768   and PI_APPROX_RULE p =
 769     let q1 = p + 5
 770     and q2 = p + 6 in
 771     let th1 = ATN_APPROX_CONV q1 tm_1_8
 772     and th2 = ATN_APPROX_CONV q2 tm_1_57
 773     and th3 = ATN_APPROX_CONV q2 tm_1_239 in
 774     let th4 = INST [mk_small_numeral p,p_tm;
 775                     mk_small_numeral q1,q1_tm;
 776                     mk_small_numeral q2,q2_tm] pth in
 777     let th5 = EQT_ELIM(NUM_REDUCE_CONV(lhand(lhand(concl th4))))
 778     and th6 = EQT_ELIM(NUM_REDUCE_CONV(lhand(rand(lhand(concl th4))))) in
 779     let th7 = MATCH_MP th4 (end_itlist CONJ [th5; th6; th1; th2; th3]) in
 780     CONV_RULE(LAND_CONV(RAND_CONV(RAND_CONV REAL_RAT_REDUCE_CONV))) th7 in
 781   pi_approx_rule,PI_APPROX_RULE;;
 782
 783 (* ------------------------------------------------------------------------- *)
 784 (* A version that yields a fraction with power of two denominator.           *)
 785 (* ------------------------------------------------------------------------- *)
 786
 787 let pi_approx_binary_rule,PI_APPROX_BINARY_RULE =
 788   let pth = prove
 789    (`abs(x - r) <= inv(&2 pow (SUC p))
 790      ==> !a. abs(&2 pow p * r - a) <= inv(&2)
 791              ==> abs(x - a / &2 pow p) <= inv(&2 pow p)`,
 792     REPEAT STRIP_TAC THEN
 793     FIRST_ASSUM(MATCH_MP_TAC o MATCH_MP (REAL_ARITH
 794       `abs(x - r) <= q ==> abs(r - r') <= p - q ==> abs(x - r') <= p`)) THEN
 795     MATCH_MP_TAC REAL_LE_LCANCEL_IMP THEN EXISTS_TAC `abs(&2 pow p)` THEN
 796     SIMP_TAC[REAL_ARITH `&0 < x ==> &0 < abs(x)`; REAL_POW2_THM] THEN
 797     REWRITE_TAC[GSYM REAL_ABS_MUL; REAL_SUB_LDISTRIB] THEN
 798     SIMP_TAC[REAL_ABS_POW; REAL_ABS_NUM; GSYM real_div;
 799              REAL_DIV_LMUL; REAL_LT_IMP_NZ; REAL_POW2_CLAUSES;
 800              REAL_DIV_POW2; REAL_OF_NUM_EQ; ARITH_EQ;
 801              LE_REFL; ARITH_RULE `~(SUC p <= p)`;
 802              ARITH_RULE `SUC p - p = 1`; SUB_REFL] THEN
 803     UNDISCH_TAC `abs (&2 pow p * r - a) <= inv (&2)` THEN
 804     CONV_TAC NUM_REDUCE_CONV THEN CONV_TAC REAL_RAT_REDUCE_CONV)
 805   and num_2 = Int 2 in
 806   let pi_approx_binary_rule p =
 807     let ppow = power_num num_2 (Int p) in
 808     let r = pi_approx_rule (p + 1) in
 809     let a = round_num (ppow */ r) in
 810     a // ppow
 811   and PI_APPROX_BINARY_RULE p =
 812     let ppow = power_num num_2 (Int p) in
 813     let th1 = PI_APPROX_RULE (p + 1) in
 814     let th2 = CONV_RULE(funpow 3 RAND_CONV num_CONV) th1 in
 815     let r = rat_of_term(rand(rand(lhand(concl th2)))) in
 816     let th3 = SPEC (mk_realintconst(round_num(ppow */ r))) (MATCH_MP pth th2) in
 817     let th4 = MP th3 (EQT_ELIM(REAL_RAT_REDUCE_CONV(lhand(concl th3)))) in
 818     CONV_RULE(LAND_CONV(RAND_CONV(RAND_CONV REAL_RAT_REDUCE_CONV))) th4 in
 819   pi_approx_binary_rule,PI_APPROX_BINARY_RULE;;
 820
 821 (* ------------------------------------------------------------------------- *)
 822 (* Rule to expand atn(r) for rational r into more easily calculable bits.    *)
 823 (* ------------------------------------------------------------------------- *)
 824
 825 let ATN_EXPAND_CONV =
 826   let num_0 = Int 0
 827   and num_1 = Int 1
 828   and num_2 = Int 2
 829   and eighth = Int 1 // Int 8
 830   and atn_tm = `atn`
 831   and eighth_tm = `&1 / &8`
 832   and mk_mul = mk_binop `(*)`
 833   and mk_add = mk_binop `(+)`
 834   and amp_tm = `&` in
 835   let home_in =
 836     let rec homein n x =
 837       let x' = (x -/ eighth) // (num_1 +/ x */ eighth) in
 838       if x' </ num_0 then
 839         if abs_num(x') </ abs_num(x) then (x',n+1) else (x,n)
 840       else homein (n + 1) x' in
 841     homein 0 in
 842   fun tm ->
 843     let ltm,rtm = dest_comb tm in
 844     if ltm <> atn_tm then failwith "ATN_EXPAND_CONV" else
 845     let r = rat_of_term rtm in
 846     let (x,n) = home_in r in
 847     let xtm = mk_add (mk_mul (mk_comb(amp_tm,mk_small_numeral n))
 848                              (mk_comb(atn_tm,eighth_tm)))
 849                      (mk_comb(atn_tm,term_of_rat x)) in
 850     SYM(MACHIN_CONV xtm);;