text_formalization/jordan/compute_pi.hl

   1
   2
   3 (* split off from misc_defs_and_lemmas.ml *)
   4
   5
   6 (* ------------------------------------------------------------------ *)
   7 (* COMPUTING PI *)
   8 (* ------------------------------------------------------------------ *)
   9
  10 Parse_ext_override_interface.unambiguous_interface();;
  11 Parse_ext_override_interface.prioritize_real();;
  12
  13 (* ------------------------------------------------------------------ *)
  14 (* general series approximations *)
  15 (* ------------------------------------------------------------------ *)
  16
  17 let SER_APPROX1 = prove_by_refinement(
  18   `!s f g.  (f sums s) /\ (summable g) ==>
  19     (!k. ((!n. (||. (f (n+k)) <=. (g (n+k)))) ==>
  20     ( (s - (sum(0,k) f)) <=. (suminf (\n. (g (n +| k)))))))`,
  21   (* {{{ proof *)
  22   [
  23   REPEAT GEN_TAC;
  24   DISCH_ALL_TAC;
  25   GEN_TAC;
  26   DISCH_TAC;
  27   IMP_RES_THEN ASSUME_TAC SUM_SUMMABLE;
  28   IMP_RES_THEN (fun th -> (ASSUME_TAC (SPEC `k:num` th))) SER_OFFSET;
  29   IMP_RES_THEN ASSUME_TAC SUM_UNIQ;
  30   SUBGOAL_THEN `(\n. (f (n+ k))) sums (s - (sum(0,k) f))` ASSUME_TAC;
  31   ASM_MESON_TAC[];
  32   SUBGOAL_THEN `summable (\n. (f (n+k))) /\ (suminf (\n. (f (n+k))) <=. (suminf (\n. (g (n+k)))))` ASSUME_TAC;
  33   MATCH_MP_TAC SER_LE2;
  34   BETA_TAC;
  35   ASM_REWRITE_TAC[];
  36   IMP_RES_THEN ASSUME_TAC SER_OFFSET;
  37   FIRST_X_ASSUM (fun th -> ACCEPT_TAC (MATCH_MP SUM_SUMMABLE (((SPEC `k:num`) th))));
  38   ASM_MESON_TAC[SUM_UNIQ]
  39   ]);;
  40   (* }}} *)
  41
  42 let SER_APPROX = prove_by_refinement(
  43   `!s f g.  (f sums s) /\ (!n. (||. (f n) <=. (g n))) /\
  44        (summable g) ==>
  45     (!k. (abs (s - (sum(0,k) f)) <=. (suminf (\n. (g (n +| k))))))`,
  46   (* {{{ proof *)
  47   [
  48   REPEAT GEN_TAC;
  49   DISCH_ALL_TAC;
  50   GEN_TAC;
  51   REWRITE_TAC[REAL_ABS_BOUNDS];
  52   CONJ_TAC;
  53   SUBGOAL_THEN `(!k. ((!n. (||. ((\p. (--. (f p))) (n+k))) <=. (g (n+k)))) ==> ((--.s) - (sum(0,k) (\p. (--. (f p)))) <=. (suminf (\n. (g (n +k))))))` ASSUME_TAC;
  54   MATCH_MP_TAC SER_APPROX1;
  55   ASM_REWRITE_TAC[];
  56   MATCH_MP_TAC SER_NEG ;
  57   ASM_REWRITE_TAC[];
  58   MATCH_MP_TAC (REAL_ARITH (`(--. s -. (--. u) <=. x) ==> (--. x <=. (s -. u))`));
  59   ONCE_REWRITE_TAC[GSYM SUM_NEG];
  60   FIRST_X_ASSUM (fun th -> (MATCH_MP_TAC th));
  61   BETA_TAC;
  62   ASM_REWRITE_TAC[REAL_ABS_NEG];
  63   H_VAL2 CONJ (HYP "0") (HYP "2");
  64   IMP_RES_THEN MATCH_MP_TAC SER_APPROX1 ;
  65   GEN_TAC;
  66   ASM_MESON_TAC[];
  67   ]);;
  68   (* }}} *)
  69
  70 (* ------------------------------------------------------------------ *)
  71 (* now for pi calculation stuff *)
  72 (* ------------------------------------------------------------------ *)
  73
  74
  75 let local_def = local_definition "trig";;
  76
  77
  78 let PI_EST = prove_by_refinement(
  79                `!n. (1 <=| n) ==> (abs(&4 / &(8 * n + 1) -
  80             &2 / &(8 * n + 4) -
  81             &1 / &(8 * n + 5) -
  82             &1 / &(8 * n + 6)) <= &.622/(&.819))`,
  83   (* {{{ proof *)
  84    [
  85    GEN_TAC THEN DISCH_ALL_TAC;
  86    REWRITE_TAC[real_div];
  87    MATCH_MP_TAC (REWRITE_RULE[real_div] (REWRITE_RULE[REAL_RAT_REDUCE_CONV `(&.4/(&.9) +(&.2/(&.12)) + (&.1/(&.13))+ (&.1/(&.14)))`] (REAL_ARITH `(abs((&.4)*.u)<=. (&.4)/(&.9)) /\ (abs((&.2)*.v)<=. (&.2)/(&.12)) /\ (abs((&.1)*w) <=. (&.1)/(&.13)) /\ (abs((&.1)*x) <=. (&.1)/(&.14)) ==> (abs((&.4)*u -(&.2)*v - (&.1)*w - (&.1)*x) <= (&.4/(&.9) +(&.2/(&.12)) + (&.1/(&.13))+ (&.1/(&.14))))`)));
  88    IMP_RES_THEN ASSUME_TAC (ARITH_RULE `1 <=| n ==> (0 < n)`);
  89    FIRST_X_ASSUM (fun th -> ASSUME_TAC (REWRITE_RULE[GSYM REAL_OF_NUM_LT] th));
  90    ASSUME_TAC (prove(`(a<=.b) ==> (&.n*a <=. (&.n)*b)`,MESON_TAC[REAL_PROP_LE_LMUL;REAL_POS]));
  91    REWRITE_TAC[REAL_ABS_MUL;REAL_ABS_INV;prove(`||.(&.n) = (&.n)`,MESON_TAC[REAL_POS;REAL_ABS_REFL])];
  92    REPEAT CONJ_TAC THEN (POP_ASSUM (fun th -> MATCH_MP_TAC th)) THEN (MATCH_MP_TAC (prove(`((&.0 <. (&.n)) /\ (&.n <=. a)) ==> (inv(a)<=. (inv(&.n)))`,MESON_TAC[REAL_ABS_REFL;REAL_ABS_INV;REAL_LE_INV2]))) THEN
  93    REWRITE_TAC[REAL_LT;REAL_LE] THEN (H_UNDISCH_TAC (HYP"0")) THEN
  94    ARITH_TAC]);;
  95   (* }}} *)
  96
  97 let pi_fun = local_def `pi_fun n = inv (&.16 **. n) *.
  98           (&.4 / &.(8 *| n +| 1) -.
  99            &.2 / &.(8 *| n +| 4) -.
 100            &.1 / &.(8 *| n +| 5) -.
 101            &.1 / &.(8 *| n +| 6))`;;
 102
 103 let pi_bound_fun = local_def `pi_bound_fun n = if (n=0) then (&.8) else
 104     (((&.15)/(&.16))*(inv(&.16 **. n))) `;;
 105
 106 let PI_EST2 = prove_by_refinement(
 107     `!k. abs(pi_fun k) <=. (pi_bound_fun k)`,
 108   (* {{{ proof *)
 109    [
 110    GEN_TAC;
 111    REWRITE_TAC[pi_fun;pi_bound_fun];
 112    COND_CASES_TAC;
 113    ASM_REWRITE_TAC[];
 114    CONV_TAC (NUM_REDUCE_CONV);
 115    (CONV_TAC (REAL_RAT_REDUCE_CONV));
 116    CONV_TAC (RAND_CONV (REWR_CONV (REAL_ARITH `a*b = b*.a`)));
 117    REWRITE_TAC[REAL_ABS_MUL;REAL_ABS_INV;REAL_ABS_POW;prove(`||.(&.n) = (&.n)`,MESON_TAC[REAL_POS;REAL_ABS_REFL])];
 118    MATCH_MP_TAC (prove(`!x y z. (&.0 <. z /\ (y <=. x) ==> (z*y <=. (z*x)))`,MESON_TAC[REAL_LE_LMUL_EQ]));
 119    ASSUME_TAC (REWRITE_RULE[] (REAL_RAT_REDUCE_CONV `(&.622)/(&.819) <=. (&.15)/(&.16)`));
 120    IMP_RES_THEN ASSUME_TAC (ARITH_RULE `~(k=0) ==> (1<=| k)`);
 121    IMP_RES_THEN ASSUME_TAC (PI_EST);
 122    CONJ_TAC;
 123    SIMP_TAC[REAL_POW_LT;REAL_LT_INV;ARITH_RULE `&.0 < (&.16)`];
 124    ASM_MESON_TAC[REAL_LE_TRANS];
 125    ]);;
 126   (* }}} *)
 127
 128 let GP16 = prove_by_refinement(
 129   `!k. (\n. inv (&16 pow k) * inv (&16 pow n)) sums
 130          inv (&16 pow k) * &16 / &15`,
 131   (* {{{ proof *)
 132   [
 133   GEN_TAC;
 134   ASSUME_TAC (REWRITE_RULE[] (REAL_RAT_REDUCE_CONV `abs (&.1 / (&. 16)) <. (&.1)`));
 135   IMP_RES_THEN (fun th -> ASSUME_TAC (CONV_RULE REAL_RAT_REDUCE_CONV th)) GP;
 136   MATCH_MP_TAC SER_CMUL;
 137   ASM_REWRITE_TAC[GSYM REAL_POW_INV;REAL_INV_1OVER];
 138   ]);;
 139   (* }}} *)
 140
 141 let GP16a = prove_by_refinement(
 142    `!k. (0<|k) ==> (\n. (pi_bound_fun (n+k))) sums (inv(&.16 **. k))`,
 143   (* {{{ proof *)
 144    [
 145    GEN_TAC;
 146    DISCH_TAC;
 147    SUBGOAL_THEN `(\n. pi_bound_fun (n+k)) = (\n. ((&.15/(&.16))* (inv(&.16)**. k) *. inv(&.16 **. n)))` (fun th-> REWRITE_TAC[th]);
 148    MATCH_MP_TAC EQ_EXT;
 149    X_GEN_TAC `n:num` THEN BETA_TAC;
 150    REWRITE_TAC[pi_bound_fun];
 151    COND_CASES_TAC;
 152    ASM_MESON_TAC[ARITH_RULE `0<| k ==> (~(n+k = 0))`];
 153    REWRITE_TAC[GSYM REAL_MUL_ASSOC];
 154    AP_TERM_TAC;
 155    REWRITE_TAC[REAL_INV_MUL;REAL_POW_ADD;REAL_POW_INV;REAL_MUL_AC];
 156    SUBGOAL_THEN `(\n. (&.15/(&.16)) *. ((inv(&.16)**. k)*. inv(&.16 **. n))) sums ((&.15/(&.16)) *.(inv(&.16**. k)*. ((&.16)/(&.15))))` ASSUME_TAC;
 157    MATCH_MP_TAC SER_CMUL;
 158    REWRITE_TAC[REAL_POW_INV];
 159    ACCEPT_TAC (SPEC `k:num` GP16);
 160    FIRST_X_ASSUM MP_TAC;
 161    REWRITE_TAC[REAL_MUL_ASSOC];
 162    MATCH_MP_TAC (prove (`(x=y) ==> ((a sums x) ==> (a sums y))`,MESON_TAC[]));
 163    MATCH_MP_TAC (REAL_ARITH `(b*(a*c) = (b*(&.1))) ==> ((a*b)*c = b)`);
 164    AP_TERM_TAC;
 165    CONV_TAC (REAL_RAT_REDUCE_CONV);
 166    ]);;
 167   (* }}} *)
 168
 169 let PI_SER = prove_by_refinement(
 170   `!k. (0<|k) ==> (abs(pi - (sum(0,k) pi_fun)) <=. (inv(&.16 **. (k))))`,
 171   (* {{{ proof *)
 172    [
 173    GEN_TAC THEN DISCH_TAC;
 174    ASSUME_TAC (ONCE_REWRITE_RULE[ETA_AX] (REWRITE_RULE[GSYM pi_fun] POLYLOG_THM));
 175    ASSUME_TAC PI_EST2;
 176    IMP_RES_THEN (ASSUME_TAC) GP16a;
 177    IMP_RES_THEN (ASSUME_TAC) SUM_SUMMABLE;
 178    IMP_RES_THEN (ASSUME_TAC) SER_OFFSET_REV;
 179    IMP_RES_THEN (ASSUME_TAC) SUM_SUMMABLE;
 180    MP_TAC (SPECL [`pi`;`pi_fun`;`pi_bound_fun` ] SER_APPROX);
 181    ASM_REWRITE_TAC[];
 182    DISCH_THEN (fun th -> MP_TAC (SPEC `k:num` th));
 183    SUBGOAL_THEN `suminf (\n. pi_bound_fun (n + k)) = inv (&.16 **. k)` (fun th -> (MESON_TAC[th]));
 184    ASM_MESON_TAC[SUM_UNIQ];
 185    ]);;
 186   (* }}} *)
 187
 188 (* replace 3 by SUC (SUC (SUC 0)) *)
 189 let SUC_EXPAND_CONV tm =
 190    let count = dest_numeral tm in
 191    let rec add_suc i r =
 192      if (i <=/ (Int 0)) then r
 193      else add_suc (i -/ (Int 1)) (mk_comb (`SUC`,r)) in
 194    let tm' = add_suc count `0` in
 195    REWRITE_RULE[] (ARITH_REWRITE_CONV[] (mk_eq (tm,tm')));;
 196
 197 let inv_twopow = prove(
 198   `!n. inv (&.16 **. n) = (twopow (--: (&:(4*n)))) `,
 199     REWRITE_TAC[TWOPOW_NEG;GSYM (NUM_RED_CONV `2 EXP 4`);
 200     REAL_OF_NUM_POW;EXP_MULT]);;
 201
 202 let PI_SERn n =
 203    let SUM_EXPAND_CONV =
 204            (ARITH_REWRITE_CONV[]) THENC
 205            (TOP_DEPTH_CONV SUC_EXPAND_CONV) THENC
 206            (REWRITE_CONV[sum]) THENC
 207            (ARITH_REWRITE_CONV[REAL_ADD_LID;GSYM REAL_ADD_ASSOC]) in
 208    let sum_thm = SUM_EXPAND_CONV (vsubst [n,`i:num`] `sum(0,i) f`) in
 209    let gt_thm = ARITH_RULE (vsubst [n,`i:num`] `0 <| i`) in
 210    ((* CONV_RULE REAL_RAT_REDUCE_CONV *)(CONV_RULE (ARITH_REWRITE_CONV[]) (BETA_RULE (REWRITE_RULE[sum_thm;pi_fun;inv_twopow] (MATCH_MP PI_SER gt_thm)))));;
 211
 212 (* abs(pi - u ) < e *)
 213 let recompute_pi bprec =
 214    let n = (bprec /4) in
 215    let pi_ser = PI_SERn (mk_numeral (Int n)) in
 216    let _ = remove_real_constant `pi` in
 217    (add_real_constant pi_ser; INTERVAL_OF_TERM bprec `pi`);;
 218
 219 (* ------------------------------------------------------------------ *)
 220 (* restore defaults *)
 221 (* ------------------------------------------------------------------ *)
 222
 223 reduce_local_interface("trig");;