(* ========================================================================= *)
(* Stirling's approximation.                                                 *)
(* ========================================================================= *)

needs "Library/analysis.ml";;
needs "Library/transc.ml";;

override_interface("-->",`(tends_num_real)`);;

(* ------------------------------------------------------------------------- *)
(* This is a handy induction for Wallis's product below.                     *)
(* ------------------------------------------------------------------------- *)

let ODDEVEN_INDUCT = prove
 (`!P. P 0 /\ P 1 /\ (!n. P n ==> P(n + 2)) ==> !n. P n`,
  GEN_TAC THEN STRIP_TAC THEN
  SUBGOAL_THEN `!n. P n /\ P(n + 1)` (fun th -> MESON_TAC[th]) THEN
  INDUCT_TAC THEN ASM_SIMP_TAC[ADD1; GSYM ADD_ASSOC] THEN
  ASM_SIMP_TAC[ARITH]);;

(* ------------------------------------------------------------------------- *)
(* A particular limit we need below.                                         *)
(* ------------------------------------------------------------------------- *)

let LN_LIM_BOUND = prove
 (`!n. ~(n = 0) ==> abs(&n * ln(&1 + &1 / &n) - &1) <= &1 / (&2 * &n)`,
  REPEAT STRIP_TAC THEN MP_TAC(SPECL [`&1 / &n`; `2`] MCLAURIN_LN_POS) THEN
  ASM_SIMP_TAC[SUM_2; REAL_LT_DIV; REAL_OF_NUM_LT; LT_NZ; ARITH] THEN
  DISCH_THEN(X_CHOOSE_THEN `t:real` STRIP_ASSUME_TAC) THEN
  ASM_REWRITE_TAC[] THEN
  REWRITE_TAC[real_div; REAL_INV_0; REAL_MUL_RZERO; REAL_ADD_LID] THEN
  REWRITE_TAC[REAL_POW_1; REAL_POW_2; REAL_MUL_LNEG; REAL_MUL_RNEG;
              REAL_NEG_NEG; REAL_MUL_LID; REAL_INV_1; REAL_POW_NEG;
              REAL_POW_ONE; ARITH; REAL_MUL_RID] THEN
  ASM_SIMP_TAC[REAL_OF_NUM_EQ; REAL_FIELD
   `~(x = &0) ==> x * (inv(x) + a) - &1 = x * a`] THEN
  REWRITE_TAC[REAL_ARITH `n * --((i * i) * a) = --((n * i) * a * i)`] THEN
  ASM_SIMP_TAC[REAL_MUL_RINV; REAL_OF_NUM_EQ; ARITH; REAL_MUL_LID] THEN
  REWRITE_TAC[REAL_ABS_NEG; REAL_ABS_MUL] THEN
  ONCE_REWRITE_TAC[REAL_INV_MUL] THEN REWRITE_TAC[REAL_ABS_MUL] THEN
  CONV_TAC REAL_RAT_REDUCE_CONV THEN REWRITE_TAC[GSYM REAL_MUL_ASSOC] THEN
  MATCH_MP_TAC REAL_LE_LMUL THEN CONV_TAC REAL_RAT_REDUCE_CONV THEN
  GEN_REWRITE_TAC RAND_CONV [GSYM REAL_MUL_LID] THEN
  MATCH_MP_TAC REAL_LE_RMUL THEN
  ASM_SIMP_TAC[real_div; REAL_MUL_LID; REAL_LE_INV_EQ; REAL_POS] THEN
  REWRITE_TAC[REAL_ABS_INV] THEN MATCH_MP_TAC REAL_INV_LE_1 THEN
  REWRITE_TAC[REAL_ABS_MUL] THEN
  GEN_REWRITE_TAC LAND_CONV [GSYM REAL_MUL_LID] THEN
  MATCH_MP_TAC REAL_LE_MUL2 THEN REWRITE_TAC[REAL_POS] THEN
  UNDISCH_TAC `&0 < t` THEN REAL_ARITH_TAC);;

let LN_LIM_LEMMA = prove
 (`(\n. &n * ln(&1 + &1 / &n)) --> &1`,
  GEN_REWRITE_TAC (LAND_CONV o BINDER_CONV)
   [REAL_ARITH `a = (a - &1) + &1`] THEN
  GEN_REWRITE_TAC RAND_CONV [GSYM REAL_ADD_LID] THEN
  MATCH_MP_TAC SEQ_ADD THEN REWRITE_TAC[SEQ_CONST] THEN
  MATCH_MP_TAC SEQ_LE_0 THEN EXISTS_TAC `\n. &1 / &n` THEN
  REWRITE_TAC[SEQ_HARMONIC] THEN
  EXISTS_TAC `1` THEN REWRITE_TAC[ARITH_RULE `n >= 1 <=> ~(n = 0)`] THEN
  REPEAT STRIP_TAC THEN MATCH_MP_TAC REAL_LE_TRANS THEN
  EXISTS_TAC `&1 / (&2 * &n)` THEN ASM_SIMP_TAC[LN_LIM_BOUND] THEN
  REWRITE_TAC[real_div; REAL_MUL_LID; REAL_ABS_INV] THEN
  MATCH_MP_TAC REAL_LE_INV2 THEN UNDISCH_TAC `~(n = 0)` THEN
  REWRITE_TAC[GSYM REAL_OF_NUM_EQ] THEN REAL_ARITH_TAC);;

(* ------------------------------------------------------------------------- *)
(* Lemma for proving inequality via derivative and limit at infinity.        *)
(* ------------------------------------------------------------------------- *)

let POSITIVE_DIFF_LEMMA = prove
 (`!f f'. (!x. &0 < x ==> (f diffl f'(x)) x /\ f'(x) < &0) /\
          (\n. f(&n)) --> &0
          ==> !n. ~(n = 0) ==> &0 < f(&n)`,
  REPEAT STRIP_TAC THEN REWRITE_TAC[GSYM REAL_NOT_LE] THEN DISCH_TAC THEN
  SUBGOAL_THEN `!m p. n <= m /\ m < p ==> (f:real->real)(&p) < f(&m)`
  ASSUME_TAC THENL
   [REPEAT STRIP_TAC THEN
    MP_TAC(SPECL [`f:real->real`; `f':real->real`; `&m`; `&p`] MVT_ALT) THEN
    ANTS_TAC THENL
     [ASM_MESON_TAC[LT_NZ; LTE_TRANS; REAL_OF_NUM_LT; REAL_LTE_TRANS];
      ALL_TAC] THEN
    REWRITE_TAC[REAL_EQ_SUB_RADD] THEN STRIP_TAC THEN ASM_REWRITE_TAC[] THEN
    MATCH_MP_TAC(REAL_ARITH `&0 < z * --y ==> z * y + a < a`) THEN
    MATCH_MP_TAC REAL_LT_MUL THEN
    ASM_REWRITE_TAC[REAL_SUB_LT; REAL_OF_NUM_LT] THEN
    REWRITE_TAC[REAL_ARITH `&0 < --x <=> x < &0`] THEN
    ASM_MESON_TAC[LT_NZ; LTE_TRANS; REAL_OF_NUM_LT; REAL_LT_TRANS];
    ALL_TAC] THEN
  SUBGOAL_THEN `f(&(n + 1)) < &0` ASSUME_TAC THENL
   [FIRST_ASSUM(MP_TAC o SPECL [`n:num`; `n + 1`]) THEN ANTS_TAC THENL
     [ARITH_TAC; UNDISCH_TAC `f(&n) <= &0` THEN REAL_ARITH_TAC];
    ALL_TAC] THEN
  FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE I [SEQ]) THEN
  DISCH_THEN(MP_TAC o SPEC `--f(&(n + 1))`) THEN
  ASM_REWRITE_TAC[REAL_SUB_RZERO; REAL_ARITH `&0 < --x <=> x < &0`] THEN
  DISCH_THEN(X_CHOOSE_THEN `p:num` (MP_TAC o SPEC `n + p + 2`)) THEN
  ANTS_TAC THENL [ARITH_TAC; ALL_TAC] THEN
  MATCH_MP_TAC(REAL_ARITH `y < &0 /\ z < y ==> abs(z) < --y ==> F`) THEN
  ASM_REWRITE_TAC[] THEN FIRST_X_ASSUM MATCH_MP_TAC THEN ARITH_TAC);;

(* ------------------------------------------------------------------------- *)
(* Auxiliary definition.                                                     *)
(* ------------------------------------------------------------------------- *)

let stirling = new_definition
 `stirling n = ln(&(FACT n)) - ((&n + &1 / &2) * ln(&n) - &n)`;;

(* ------------------------------------------------------------------------- *)
(* This difference is a decreasing sequence.                                 *)
(* ------------------------------------------------------------------------- *)

let STIRLING_DIFF = prove
 (`!n. ~(n = 0)
       ==> stirling(n) - stirling(n + 1) =
           (&n + &1 / &2) * ln((&n + &1) / &n) - &1`,
  REPEAT STRIP_TAC THEN REWRITE_TAC[stirling] THEN
  MATCH_MP_TAC(REAL_ARITH
   `(f' - f) + x = (nl' - nl) /\ n' = n + &1
    ==> (f - (nl - n)) - (f' - (nl' - n')) = x - &1`) THEN
  REWRITE_TAC[REAL_OF_NUM_ADD] THEN
  REWRITE_TAC[REWRITE_RULE[ADD1] FACT; GSYM REAL_OF_NUM_MUL] THEN
  SIMP_TAC[LN_MUL; FACT_LT; ADD_EQ_0; ARITH; LT_NZ; REAL_OF_NUM_LT] THEN
  ASM_SIMP_TAC[LN_DIV; REAL_OF_NUM_LT; ADD_EQ_0; ARITH; LT_NZ] THEN
  REWRITE_TAC[GSYM REAL_OF_NUM_ADD] THEN REAL_ARITH_TAC);;

let STIRLING_DELTA_DERIV = prove
 (`!x. &0 < x
       ==> ((\x. ln ((x + &1) / x) - &1 / (x + &1 / &2)) diffl
            (-- &1 / (x * (x + &1) * (&2 * x + &1) pow 2))) x`,
  GEN_TAC THEN DISCH_TAC THEN
  W(fun (asl,w) -> MP_TAC(SPEC(rand w) (DIFF_CONV(lhand(rator w))))) THEN
  REWRITE_TAC[IMP_IMP] THEN ANTS_TAC THENL
   [REPEAT CONJ_TAC THEN TRY(MATCH_MP_TAC REAL_LT_DIV) THEN
    POP_ASSUM MP_TAC THEN CONV_TAC REAL_FIELD;
    ALL_TAC] THEN
  MATCH_MP_TAC EQ_IMP THEN
  AP_THM_TAC THEN AP_TERM_TAC THEN REWRITE_TAC[REAL_POW_2] THEN
  POP_ASSUM MP_TAC THEN CONV_TAC REAL_FIELD);;

let STIRLING_DELTA_LIMIT = prove
 (`(\n. ln ((&n + &1) / &n) - &1 / (&n + &1 / &2)) --> &0`,
  GEN_REWRITE_TAC RAND_CONV [GSYM REAL_SUB_RZERO] THEN
  MATCH_MP_TAC SEQ_SUB THEN CONJ_TAC THEN MATCH_MP_TAC SEQ_LE_0 THEN
  EXISTS_TAC `\n. &1 / &n` THEN REWRITE_TAC[SEQ_HARMONIC] THEN
  EXISTS_TAC `1` THEN X_GEN_TAC `n:num` THEN
  REWRITE_TAC[GE; GSYM REAL_OF_NUM_LE] THEN
  DISCH_TAC THEN MATCH_MP_TAC
   (REAL_ARITH `&0 <= x /\ x <= y ==> abs x <= abs y`)
  THEN CONJ_TAC THENL
   [MATCH_MP_TAC LN_POS THEN
    ASM_SIMP_TAC[REAL_LE_RDIV_EQ; REAL_ARITH `&1 <= x ==> &0 < x`] THEN
    REAL_ARITH_TAC;
    ASM_SIMP_TAC[REAL_FIELD `&1 <= x ==> (x + &1) / x = &1 + &1 / x`] THEN
    MATCH_MP_TAC LN_LE THEN ASM_SIMP_TAC[REAL_LE_DIV; REAL_POS];
    MATCH_MP_TAC REAL_LE_DIV THEN REAL_ARITH_TAC;
    REWRITE_TAC[real_div; REAL_MUL_LID] THEN MATCH_MP_TAC REAL_LE_INV2 THEN
    POP_ASSUM MP_TAC THEN REAL_ARITH_TAC]);;

let STIRLING_DECREASES = prove
 (`!n. ~(n = 0) ==> stirling(n + 1) < stirling(n)`,
  ONCE_REWRITE_TAC[GSYM REAL_SUB_LT] THEN SIMP_TAC[STIRLING_DIFF] THEN
  REWRITE_TAC[REAL_SUB_LT] THEN
  ONCE_REWRITE_TAC[REAL_MUL_SYM] THEN
  SIMP_TAC[GSYM REAL_LT_LDIV_EQ; REAL_ARITH `&0 < &n + &1 / &2`] THEN
  ONCE_REWRITE_TAC[GSYM REAL_SUB_LT] THEN
  MATCH_MP_TAC POSITIVE_DIFF_LEMMA THEN
  EXISTS_TAC `\x. -- &1 / (x * (x + &1) * (&2 * x + &1) pow 2)` THEN
  SIMP_TAC[STIRLING_DELTA_DERIV; STIRLING_DELTA_LIMIT] THEN
  X_GEN_TAC `x:real` THEN DISCH_TAC THEN
  REWRITE_TAC[real_div; REAL_MUL_LNEG; REAL_MUL_LID] THEN
  REWRITE_TAC[REAL_ARITH `--x < &0 <=> &0 < x`; REAL_LT_INV_EQ] THEN
  REWRITE_TAC[REAL_POW_2] THEN
  REPEAT(MATCH_MP_TAC REAL_LT_MUL THEN CONJ_TAC) THEN
  POP_ASSUM MP_TAC THEN REAL_ARITH_TAC);;

(* ------------------------------------------------------------------------- *)
(* However a slight tweak gives an *increasing* sequence.                    *)
(* ------------------------------------------------------------------------- *)

let OTHER_DERIV_LEMMA = prove
 (`!x. &0 < x
       ==> ((\x. &1 / (&12 * x * (x + &1) * (x + &1 / &2))) diffl
            --(&3 * x pow 2 + &3 * x + &1 / &2) /
              (&12 * (x * (x + &1) * (x + &1 / &2)) pow 2)) x`,
  REPEAT STRIP_TAC THEN
  W(fun (asl,w) -> MP_TAC(SPEC(rand w) (DIFF_CONV(lhand(rator w))))) THEN
  REWRITE_TAC[IMP_IMP] THEN ANTS_TAC THENL
   [REWRITE_TAC[REAL_ENTIRE] THEN POP_ASSUM MP_TAC THEN ARITH_TAC;
    ALL_TAC] THEN
  MATCH_MP_TAC EQ_IMP THEN
  AP_THM_TAC THEN AP_TERM_TAC THEN REWRITE_TAC[REAL_POW_2] THEN
  POP_ASSUM MP_TAC THEN CONV_TAC REAL_FIELD);;

let STIRLING_INCREASES = prove
 (`!n. ~(n = 0)
       ==> stirling(n + 1) - &1 / (&12 * (&(n + 1)))
           > stirling(n) - &1 / (&12 * &n)`,
  REWRITE_TAC[GSYM REAL_OF_NUM_EQ; GSYM REAL_OF_NUM_ADD] THEN
  REWRITE_TAC[REAL_ARITH `a - b > c - d <=> c - a < d - b`] THEN
  SIMP_TAC[REAL_FIELD
    `~(&n = &0) ==> &1 / (&12 * &n) - &1 / (&12 * (&n + &1)) =
                    &1 / (&12 * &n * (&n + &1))`] THEN
  SIMP_TAC[REAL_OF_NUM_EQ; STIRLING_DIFF] THEN
  REWRITE_TAC[REAL_ARITH `a * b - &1 < c <=> b * a < c + &1`] THEN
  SIMP_TAC[GSYM REAL_LT_RDIV_EQ; REAL_ARITH `&0 < &n + &1 / &2`] THEN
  REWRITE_TAC[REAL_ARITH `(&1 / x + &1) / y = &1 / x / y + &1 / y`] THEN
  REWRITE_TAC[REAL_ARITH `a < b + c <=> &0 < b - (a - c)`] THEN
  REWRITE_TAC[real_div; GSYM REAL_MUL_ASSOC; GSYM REAL_INV_MUL] THEN
  REWRITE_TAC[GSYM real_div] THEN MATCH_MP_TAC POSITIVE_DIFF_LEMMA THEN
  EXISTS_TAC `\x. &1 / (x * (x + &1) * (&2 * x + &1) pow 2) -
                  (&3 * x pow 2 + &3 * x + &1 / &2) /
                  (&12 * (x * (x + &1) * (x + &1 / &2)) pow 2)` THEN
  REWRITE_TAC[] THEN CONJ_TAC THENL
   [ALL_TAC;
    GEN_REWRITE_TAC RAND_CONV [GSYM REAL_SUB_RZERO] THEN
    MATCH_MP_TAC SEQ_SUB THEN REWRITE_TAC[STIRLING_DELTA_LIMIT] THEN
    REWRITE_TAC[real_div; REAL_INV_MUL; REAL_MUL_LID] THEN
    REWRITE_TAC[REAL_FIELD
     `inv(&12) * x * y * inv(&n + inv(&2)) = x * y * inv(&12 * &n + &6)`] THEN
    GEN_REWRITE_TAC RAND_CONV [SYM(REAL_RAT_REDUCE_CONV `&0 * &0 * &0`)] THEN
    REPEAT(MATCH_MP_TAC SEQ_MUL THEN CONJ_TAC) THEN
    MP_TAC(SPEC `&1` SEQ_HARMONIC) THEN
    REWRITE_TAC[real_div; REAL_MUL_LID] THEN
    DISCH_THEN(MP_TAC o MATCH_MP SEQ_SUBSEQ) THENL
     [DISCH_THEN(MP_TAC o SPECL [`1`; `1`]);
      DISCH_THEN(MP_TAC o SPECL [`12`; `6`])] THEN
    REWRITE_TAC[REAL_OF_NUM_ADD; REAL_OF_NUM_MUL; ARITH; MULT_CLAUSES]] THEN
  X_GEN_TAC `x:real` THEN DISCH_TAC THEN CONJ_TAC THENL
   [GEN_REWRITE_TAC LAND_CONV
     [REAL_ARITH `&1 / x - y / z = --y / z - -- &1 / x`] THEN
    MATCH_MP_TAC DIFF_SUB THEN
    ASM_SIMP_TAC[STIRLING_DELTA_DERIV; OTHER_DERIV_LEMMA];
    ALL_TAC] THEN
  REWRITE_TAC[REAL_ARITH `a - b < &0 <=> a < b`] THEN
  ASM_SIMP_TAC[GSYM REAL_POW_2; REAL_FIELD
   `&0 < x
    ==> &1 / (x * (x + &1) * (&2 * x + &1) pow 2) =
        (&3 * x * (x + &1)) /
        (&12 * (x * (x + &1) * (x + &1 / &2)) *
               (x * (x + &1) * (x + &1 / &2)))`] THEN
  ONCE_REWRITE_TAC[real_div] THEN MATCH_MP_TAC REAL_LT_RMUL THEN
  CONJ_TAC THENL [REAL_ARITH_TAC; ALL_TAC] THEN
  REWRITE_TAC[REAL_LT_INV_EQ; REAL_POW_2] THEN
  REPEAT(MATCH_MP_TAC REAL_LT_MUL THEN CONJ_TAC) THEN
  POP_ASSUM MP_TAC THEN REAL_ARITH_TAC);;

(* ------------------------------------------------------------------------- *)
(* Hence it converges to *something*.                                        *)
(* ------------------------------------------------------------------------- *)

let STIRLING_UPPERBOUND = prove
 (`!n. stirling(SUC n) <= &1`,
  INDUCT_TAC THENL
   [REWRITE_TAC[stirling] THEN CONV_TAC NUM_REDUCE_CONV THEN
    REWRITE_TAC[LN_1] THEN CONV_TAC REAL_RAT_REDUCE_CONV;
    ALL_TAC] THEN
  MATCH_MP_TAC REAL_LE_TRANS THEN
  EXISTS_TAC `stirling(SUC n)` THEN ASM_REWRITE_TAC[] THEN
  REWRITE_TAC[ARITH_RULE `SUC(SUC n) = SUC n + 1`] THEN
  MATCH_MP_TAC REAL_LT_IMP_LE THEN MATCH_MP_TAC STIRLING_DECREASES THEN
  ARITH_TAC);;

let STIRLING_LOWERBOUND = prove
 (`!n. -- &1 <= stirling(SUC n)`,
  GEN_TAC THEN MATCH_MP_TAC REAL_LE_TRANS THEN
  EXISTS_TAC `stirling(SUC n) - &1 / (&12 * &(SUC n))` THEN CONJ_TAC THENL
   [ALL_TAC;
    SIMP_TAC[REAL_ARITH `a - b <= a <=> &0 <= b`] THEN
    REWRITE_TAC[real_div; REAL_MUL_LID; REAL_LE_INV_EQ] THEN
    REWRITE_TAC[REAL_OF_NUM_MUL; REAL_POS]] THEN
  SPEC_TAC(`n:num`,`n:num`) THEN INDUCT_TAC THENL
   [REWRITE_TAC[stirling] THEN CONV_TAC NUM_REDUCE_CONV THEN
    REWRITE_TAC[LN_1] THEN CONV_TAC REAL_RAT_REDUCE_CONV;
    ALL_TAC] THEN
  MATCH_MP_TAC REAL_LE_TRANS THEN
  EXISTS_TAC `stirling(SUC n) - &1 / (&12 * &(SUC n))` THEN
  ASM_REWRITE_TAC[] THEN
  REWRITE_TAC[ARITH_RULE `SUC(SUC n) = SUC n + 1`] THEN
  MATCH_MP_TAC(REAL_ARITH `b > a ==> a <= b`) THEN
  MATCH_MP_TAC STIRLING_INCREASES THEN ARITH_TAC);;

let STIRLING_MONO = prove
 (`!m n. ~(m = 0) /\ m <= n ==> stirling n <= stirling m`,
  REPEAT STRIP_TAC THEN
  FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE I [LE_EXISTS]) THEN
  DISCH_THEN(X_CHOOSE_THEN `d:num` SUBST1_TAC) THEN
  SPEC_TAC(`d:num`,`d:num`) THEN INDUCT_TAC THEN
  REWRITE_TAC[ADD_CLAUSES; REAL_LE_REFL] THEN
  MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `stirling(m + d)` THEN
  ASM_SIMP_TAC[ADD1; REAL_LT_IMP_LE; STIRLING_DECREASES; ADD_EQ_0]);;

let STIRLING_CONVERGES = prove
 (`?c. stirling --> c`,
  ONCE_REWRITE_TAC[SEQ_SUC] THEN
  REWRITE_TAC[GSYM convergent] THEN MATCH_MP_TAC SEQ_BCONV THEN CONJ_TAC THENL
   [ALL_TAC;
    REWRITE_TAC[mono; real_ge] THEN DISJ2_TAC THEN REPEAT GEN_TAC THEN
    DISCH_TAC THEN MATCH_MP_TAC STIRLING_MONO THEN
    POP_ASSUM MP_TAC THEN ARITH_TAC] THEN
  REWRITE_TAC[MR1_BOUNDED; GE; LE_REFL] THEN
  MAP_EVERY EXISTS_TAC [`&2`; `0`] THEN REWRITE_TAC[LE_0] THEN
  SIMP_TAC[REAL_ARITH `-- &1 <= x /\ x <= &1 ==> abs(x) < &2`;
           STIRLING_UPPERBOUND; STIRLING_LOWERBOUND]);;

(* ------------------------------------------------------------------------- *)
(* Now derive Wallis's infinite product.                                     *)
(* ------------------------------------------------------------------------- *)

let [PI2_LT; PI2_LE; PI2_NZ] = (CONJUNCTS o prove)
 (`&0 < pi / &2 /\ &0 <= pi / &2 /\ ~(pi / &2 = &0)`,
  MP_TAC PI_POS THEN REAL_ARITH_TAC);;

let WALLIS_PARTS = prove
 (`!n. (&n + &2) * integral(&0,pi / &2) (\x. sin(x) pow (n + 2)) =
       (&n + &1) * integral(&0,pi / &2) (\x. sin(x) pow n)`,
  GEN_TAC THEN
  MP_TAC(SPECL [`\x. sin(x) pow (n + 1)`; `\x. --cos(x)`;
                `\x. (&n + &1) * sin(x) pow n * cos(x)`;
                `\x. sin(x)`; `&0`; `pi / &2`] INTEGRAL_BY_PARTS) THEN
  REWRITE_TAC[] THEN ANTS_TAC THENL
   [SIMP_TAC[REAL_LT_IMP_LE; PI_POS; REAL_LT_DIV; REAL_OF_NUM_LT; ARITH] THEN
    CONV_TAC(ONCE_DEPTH_CONV INTEGRABLE_CONV) THEN REWRITE_TAC[] THEN
    CONJ_TAC THEN GEN_TAC THEN STRIP_TAC THEN DIFF_TAC THEN
    REWRITE_TAC[REAL_OF_NUM_ADD; ADD_SUB] THEN REAL_ARITH_TAC;
    ALL_TAC] THEN
  REWRITE_TAC[SIN_PI2; COS_PI2; SIN_0; COS_0] THEN
  REWRITE_TAC[REAL_ARITH `s pow k * s = s * s pow k`] THEN
  REWRITE_TAC[GSYM(CONJUNCT2 real_pow); ARITH_RULE `SUC(n + 1) = n + 2`] THEN
  REWRITE_TAC[GSYM ADD1; real_pow; REAL_MUL_LZERO; REAL_MUL_RZERO;
              REAL_NEG_0; REAL_SUB_LZERO] THEN
  REWRITE_TAC[C MATCH_MP (SPEC_ALL SIN_CIRCLE) (REAL_RING
    `sin(x) pow 2 + cos(x) pow 2 = &1
     ==> (n * sn * cos(x)) * --cos(x) = (n * sn) * (sin(x) pow 2 - &1)`)] THEN
  REWRITE_TAC[REAL_SUB_LDISTRIB; GSYM REAL_MUL_ASSOC; GSYM REAL_POW_ADD] THEN
  REWRITE_TAC[REAL_MUL_RID] THEN
  SUBGOAL_THEN
   `integral(&0,pi / &2)
        (\x. (&n + &1) * sin x pow (n + 2) - (&n + &1) * sin x pow n) =
    (&n + &1) * (integral(&0,pi / &2) (\x. sin(x) pow (n + 2)) -
                 integral(&0,pi / &2) (\x. sin(x) pow n))`
   (fun th -> SUBST1_TAC th THEN REAL_ARITH_TAC) THEN
  REWRITE_TAC[GSYM REAL_SUB_LDISTRIB] THEN MATCH_MP_TAC EQ_TRANS THEN
  EXISTS_TAC
   `(&n + &1) *
    integral(&0,pi / &2) (\x. sin x pow (n + 2) - sin x pow n)` THEN
  CONJ_TAC THENL
   [MATCH_MP_TAC INTEGRAL_CMUL;
    AP_TERM_TAC THEN MATCH_MP_TAC INTEGRAL_SUB] THEN
  CONV_TAC(ONCE_DEPTH_CONV INTEGRABLE_CONV) THEN SIMP_TAC[PI2_LE]);;

let WALLIS_PARTS' = prove
 (`!n. integral(&0,pi / &2) (\x. sin(x) pow (n + 2)) =
       (&n + &1) / (&n + &2) * integral(&0,pi / &2) (\x. sin(x) pow n)`,
  MP_TAC WALLIS_PARTS THEN MATCH_MP_TAC MONO_FORALL THEN
  CONV_TAC REAL_FIELD);;

let WALLIS_0 = prove
 (`integral(&0,pi / &2) (\x. sin(x) pow 0) = pi / &2`,
  REWRITE_TAC[real_pow; REAL_DIV_1; REAL_MUL_LID] THEN
  SIMP_TAC[INTEGRAL_CONST; REAL_LT_IMP_LE; PI_POS; REAL_LT_DIV;
           REAL_OF_NUM_LT; ARITH; REAL_MUL_LID; REAL_SUB_RZERO]);;

let WALLIS_1 = prove
 (`integral(&0,pi / &2) (\x. sin(x) pow 1) = &1`,
  MATCH_MP_TAC DEFINT_INTEGRAL THEN REWRITE_TAC[PI2_LE; REAL_POW_1] THEN
  MP_TAC(SPECL [`\x. --cos(x)`; `\x. sin x`; `&0`; `pi / &2`] FTC1) THEN
  REWRITE_TAC[COS_0; COS_PI2] THEN CONV_TAC REAL_RAT_REDUCE_CONV THEN
  DISCH_THEN MATCH_MP_TAC THEN REWRITE_TAC[PI2_LE] THEN
  REPEAT STRIP_TAC THEN DIFF_TAC THEN REAL_ARITH_TAC);;

let WALLIS_EVEN = prove
 (`!n. integral(&0,pi / &2) (\x. sin(x) pow (2 * n)) =
         (&(FACT(2 * n)) / (&2 pow n * &(FACT n)) pow 2) * pi / &2`,
  INDUCT_TAC THENL
   [CONV_TAC NUM_REDUCE_CONV THEN CONV_TAC REAL_RAT_REDUCE_CONV THEN
    REWRITE_TAC[WALLIS_0] THEN REAL_ARITH_TAC;
    ALL_TAC] THEN
  ASM_REWRITE_TAC[ARITH_RULE `2 * SUC n = 2 * n + 2`; WALLIS_PARTS'] THEN
  REWRITE_TAC[REAL_MUL_ASSOC] THEN AP_THM_TAC THEN AP_TERM_TAC THEN
  REWRITE_TAC[FACT; real_pow; GSYM REAL_OF_NUM_MUL] THEN
  ONCE_REWRITE_TAC[REAL_ARITH `(&2 * x) * y * z = (&2 * y) * (x * z)`] THEN
  GEN_REWRITE_TAC (RAND_CONV o ONCE_DEPTH_CONV) [REAL_POW_MUL] THEN
  REWRITE_TAC[real_div; REAL_INV_MUL; REAL_MUL_ASSOC] THEN
  AP_THM_TAC THEN AP_TERM_TAC THEN
  REWRITE_TAC[ARITH_RULE `2 * n + 2 = SUC(SUC(2 * n))`] THEN
  REWRITE_TAC[GSYM REAL_OF_NUM_SUC; REAL_POW_2; FACT;
              GSYM REAL_OF_NUM_MUL] THEN
  CONV_TAC REAL_FIELD);;

let WALLIS_ODD = prove
 (`!n. integral(&0,pi / &2) (\x. sin(x) pow (2 * n + 1)) =
         (&2 pow n * &(FACT n)) pow 2 / &(FACT(2 * n + 1))`,
  INDUCT_TAC THENL
   [CONV_TAC NUM_REDUCE_CONV THEN CONV_TAC REAL_RAT_REDUCE_CONV THEN
    REWRITE_TAC[WALLIS_1] THEN CONV_TAC REAL_RAT_REDUCE_CONV;
    ALL_TAC] THEN
  ASM_REWRITE_TAC[ARITH_RULE `2 * SUC n + 1 = (2 * n + 1) + 2`;
                  WALLIS_PARTS'] THEN
  GEN_REWRITE_TAC LAND_CONV [REAL_MUL_SYM] THEN
  REWRITE_TAC[FACT; real_pow; GSYM REAL_OF_NUM_MUL] THEN
  ONCE_REWRITE_TAC[REAL_ARITH `(&2 * x) * y * z = (x * z) * (&2 * y)`] THEN
  GEN_REWRITE_TAC (RAND_CONV o ONCE_DEPTH_CONV) [REAL_POW_MUL] THEN
  REWRITE_TAC[real_div; GSYM REAL_MUL_ASSOC] THEN AP_TERM_TAC THEN
  GEN_REWRITE_TAC RAND_CONV [REAL_MUL_SYM] THEN
  REWRITE_TAC[ARITH_RULE `n + 2 = SUC(SUC n)`] THEN
  REWRITE_TAC[GSYM REAL_OF_NUM_SUC; REAL_POW_2; FACT;
              GSYM REAL_OF_NUM_ADD; GSYM REAL_OF_NUM_MUL] THEN
  MP_TAC(SPEC `2 * n + 1` FACT_LT) THEN REWRITE_TAC[GSYM REAL_OF_NUM_LT] THEN
  CONV_TAC REAL_FIELD);;

let WALLIS_QUOTIENT = prove
 (`!n. integral(&0,pi / &2) (\x. sin(x) pow (2 * n)) /
       integral(&0,pi / &2) (\x. sin(x) pow (2 * n + 1)) =
        (&(FACT(2 * n)) * &(FACT(2 * n + 1))) / (&2 pow n * &(FACT n)) pow 4 *
        pi / &2`,
  GEN_TAC THEN REWRITE_TAC[WALLIS_EVEN; WALLIS_ODD] THEN
  REWRITE_TAC[real_div; REAL_INV_MUL; GSYM REAL_POW_INV; REAL_INV_INV] THEN
  REAL_ARITH_TAC);;

let WALLIS_QUOTIENT' = prove
 (`!n. integral(&0,pi / &2) (\x. sin(x) pow (2 * n)) /
       integral(&0,pi / &2) (\x. sin(x) pow (2 * n + 1)) * &2 / pi =
         (&(FACT(2 * n)) * &(FACT(2 * n + 1))) / (&2 pow n * &(FACT n)) pow 4`,
  GEN_TAC THEN REWRITE_TAC[WALLIS_QUOTIENT] THEN
  GEN_REWRITE_TAC (LAND_CONV o LAND_CONV o RAND_CONV) [GSYM REAL_INV_DIV] THEN
  REWRITE_TAC[GSYM real_div] THEN MATCH_MP_TAC REAL_DIV_RMUL THEN
  MP_TAC PI2_NZ THEN CONV_TAC REAL_FIELD);;

let WALLIS_MONO = prove
 (`!m n. m <= n
         ==> integral(&0,pi / &2) (\x. sin(x) pow n)
                <= integral(&0,pi / &2) (\x. sin(x) pow m)`,
  REWRITE_TAC[LE_EXISTS] THEN REPEAT STRIP_TAC THEN
  ASM_REWRITE_TAC[] THEN MATCH_MP_TAC INTEGRAL_LE THEN
  CONV_TAC(ONCE_DEPTH_CONV INTEGRABLE_CONV) THEN REWRITE_TAC[PI2_LE] THEN
  REPEAT STRIP_TAC THEN REWRITE_TAC[REAL_POW_ADD] THEN
  GEN_REWRITE_TAC RAND_CONV [GSYM REAL_MUL_RID] THEN
  MATCH_MP_TAC REAL_LE_LMUL THEN CONJ_TAC THENL
   [MATCH_MP_TAC REAL_POW_LE; MATCH_MP_TAC REAL_POW_1_LE] THEN
  REWRITE_TAC[SIN_BOUNDS] THEN
  (MP_TAC(SPEC `x:real` SIN_POS_PI_LE) THEN
   ANTS_TAC THENL [ALL_TAC; REAL_ARITH_TAC] THEN
   REPEAT(POP_ASSUM MP_TAC) THEN MP_TAC PI2_LT THEN REAL_ARITH_TAC));;

let WALLIS_LT = prove
 (`!n. &0 < integral(&0,pi / &2) (\x. sin(x) pow n)`,
  MATCH_MP_TAC ODDEVEN_INDUCT THEN
  REWRITE_TAC[WALLIS_0; WALLIS_1; PI2_LT; REAL_OF_NUM_LT; ARITH] THEN
  REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[WALLIS_PARTS'] THEN
  MATCH_MP_TAC REAL_LT_MUL THEN ASM_REWRITE_TAC[] THEN
  MATCH_MP_TAC REAL_LT_DIV THEN REAL_ARITH_TAC);;

let WALLIS_NZ = prove
 (`!n. ~(integral(&0,pi / &2) (\x. sin(x) pow n) = &0)`,
  MP_TAC WALLIS_LT THEN MATCH_MP_TAC MONO_FORALL THEN REAL_ARITH_TAC);;

let WALLIS_BOUNDS = prove
 (`!n. integral(&0,pi / &2) (\x. sin(x) pow (n + 1))
        <= integral(&0,pi / &2) (\x. sin(x) pow n) /\
       integral(&0,pi / &2) (\x. sin(x) pow n) <=
        (&n + &2) / (&n + &1) * integral(&0,pi / &2) (\x. sin(x) pow (n + 1))`,
  GEN_TAC THEN SIMP_TAC[WALLIS_MONO; LE_ADD] THEN
  MATCH_MP_TAC REAL_LE_TRANS THEN
  EXISTS_TAC `(&n + &2) / (&n + &1) *
              integral (&0,pi / &2) (\x. sin x pow (n + 2))` THEN
  CONJ_TAC THENL
   [REWRITE_TAC[WALLIS_PARTS'] THEN
    MATCH_MP_TAC REAL_EQ_IMP_LE THEN CONV_TAC REAL_FIELD;
    ALL_TAC] THEN
  MATCH_MP_TAC REAL_LE_LMUL THEN
  SIMP_TAC[WALLIS_MONO; ARITH_RULE `n + 1 <= n + 2`] THEN
  MATCH_MP_TAC REAL_LE_DIV THEN REAL_ARITH_TAC);;

let WALLIS_RATIO_BOUNDS = prove
 (`!n. &1 <= integral(&0,pi / &2) (\x. sin(x) pow n) /
            integral(&0,pi / &2) (\x. sin(x) pow (n + 1)) /\
      integral(&0,pi / &2) (\x. sin(x) pow n) /
      integral(&0,pi / &2) (\x. sin(x) pow (n + 1)) <= (&n + &2) / (&n + &1)`,
  GEN_TAC THEN CONJ_TAC THENL
   [SIMP_TAC[REAL_LE_RDIV_EQ; WALLIS_LT; REAL_MUL_LID; WALLIS_BOUNDS];
    SIMP_TAC[REAL_LE_LDIV_EQ; WALLIS_LT; WALLIS_BOUNDS]]);;

let WALLIS = prove
 (`(\n. (&2 pow n * &(FACT n)) pow 4 / (&(FACT(2 * n)) * &(FACT(2 * n + 1))))
   --> pi / &2`,
  ONCE_REWRITE_TAC[GSYM REAL_INV_DIV] THEN MATCH_MP_TAC SEQ_INV THEN
  CONJ_TAC THENL [ALL_TAC; MP_TAC PI2_NZ THEN CONV_TAC REAL_FIELD] THEN
  REWRITE_TAC[GSYM WALLIS_QUOTIENT'] THEN ONCE_REWRITE_TAC[REAL_MUL_SYM] THEN
  GEN_REWRITE_TAC RAND_CONV [GSYM REAL_MUL_RID] THEN
  MATCH_MP_TAC SEQ_MUL THEN REWRITE_TAC[SEQ_CONST] THEN
  GEN_REWRITE_TAC (LAND_CONV o ABS_CONV) [REAL_ARITH `x = (x - &1) + &1`] THEN
  GEN_REWRITE_TAC RAND_CONV [GSYM REAL_ADD_LID] THEN
  MATCH_MP_TAC SEQ_ADD THEN REWRITE_TAC[SEQ_CONST] THEN
  MATCH_MP_TAC SEQ_LE_0 THEN EXISTS_TAC `\n. &1 / &n` THEN
  REWRITE_TAC[SEQ_HARMONIC] THEN EXISTS_TAC `1` THEN
  REWRITE_TAC[GE] THEN REPEAT STRIP_TAC THEN
  MATCH_MP_TAC(REAL_ARITH
   `!d. &1 <= x /\ x <= d /\ d - &1 <= e ==> abs(x - &1) <= e`) THEN
  EXISTS_TAC `(&(2 * n) + &2) / (&(2 * n) + &1)` THEN
  REWRITE_TAC[WALLIS_RATIO_BOUNDS] THEN
  REWRITE_TAC[GSYM REAL_OF_NUM_MUL] THEN
  REWRITE_TAC[REAL_FIELD
   `(&2 * &n + &2) / (&2 * &n + &1) - &1 = &1 / (&2 * &n + &1)`] THEN
  REWRITE_TAC[real_div; REAL_MUL_LID; REAL_ABS_INV; REAL_ABS_NUM] THEN
  MATCH_MP_TAC REAL_LE_INV2 THEN POP_ASSUM MP_TAC THEN
  REWRITE_TAC[GSYM REAL_OF_NUM_LE] THEN REAL_ARITH_TAC);;

(* ------------------------------------------------------------------------- *)
(* Hence determine the actual value of the limit.                            *)
(* ------------------------------------------------------------------------- *)

let LN_WALLIS = prove
 (`(\n. &4 * &n * ln(&2) + &4 * ln(&(FACT n)) -
        (ln(&(FACT(2 * n))) + ln(&(FACT(2 * n + 1))))) --> ln(pi / &2)`,
  REWRITE_TAC[REAL_ARITH `&4 * x + &4 * y - z = &4 * (x + y) - z`] THEN
  SUBGOAL_THEN `!n. &0 < &2 pow n`
   (fun th -> SIMP_TAC[th; GSYM LN_POW; REAL_OF_NUM_LT; ARITH; GSYM LN_MUL;
                       FACT_LT; REAL_POW_LT; REAL_LT_MUL; GSYM LN_DIV]) THEN
  SIMP_TAC[REAL_POW_LT; REAL_OF_NUM_LT; ARITH] THEN
  MATCH_MP_TAC CONTL_SEQ THEN REWRITE_TAC[WALLIS] THEN
  MATCH_MP_TAC DIFF_CONT THEN EXISTS_TAC `inv(pi / &2)` THEN
  MP_TAC(SPEC `pi / &2` (DIFF_CONV `\x. ln(x)`)) THEN
  SIMP_TAC[ETA_AX; PI2_LT; REAL_LT_DIV; REAL_OF_NUM_LT; ARITH] THEN
  REWRITE_TAC[REAL_MUL_RID]);;

let STIRLING = prove
 (`(\n. ln(&(FACT n)) - ((&n + &1 / &2) * ln(&n) - &n + ln(&2 * pi) / &2))
   --> &0`,
  REWRITE_TAC[REAL_ARITH `a - (b - c + d) = (a - (b - c)) - d`] THEN
  SUBST1_TAC(SYM(SPEC `ln(&2 * pi) / &2` REAL_SUB_REFL)) THEN
  MATCH_MP_TAC SEQ_SUB THEN REWRITE_TAC[SEQ_CONST] THEN
  X_CHOOSE_THEN `C:real` MP_TAC STIRLING_CONVERGES THEN
  GEN_REWRITE_TAC (funpow 2 LAND_CONV) [GSYM ETA_AX] THEN
  REWRITE_TAC[stirling] THEN DISCH_TAC THEN
  FIRST_ASSUM(MP_TAC o SPECL [`2`; `1`] o MATCH_MP SEQ_SUBSEQ) THEN
  FIRST_ASSUM(MP_TAC o SPECL [`2`; `0`] o MATCH_MP SEQ_SUBSEQ) THEN
  REWRITE_TAC[ARITH; ADD_CLAUSES; IMP_IMP] THEN
  DISCH_THEN(MP_TAC o MATCH_MP SEQ_ADD) THEN
  FIRST_ASSUM(MP_TAC o MATCH_MP SEQ_MUL o CONJ (SPEC `&4` SEQ_CONST)) THEN
  REWRITE_TAC[IMP_IMP] THEN DISCH_THEN(MP_TAC o MATCH_MP SEQ_SUB) THEN
  MP_TAC LN_WALLIS THEN REWRITE_TAC[IMP_IMP] THEN
  DISCH_THEN(MP_TAC o MATCH_MP SEQ_SUB) THEN
  REWRITE_TAC[ARITH_RULE
   `(a + &4 * x - (y + z)) - (&4 * (x - b) - ((y - c) + (z - d))) =
    (a + &4 * b) - (c + d)`] THEN
  DISCH_THEN(fun th -> POP_ASSUM MP_TAC THEN ASSUME_TAC th) THEN
  SUBGOAL_THEN `C = ln(&2 * pi) / &2` (fun th -> REWRITE_TAC[th]) THEN
  POP_ASSUM(MP_TAC o CONJ (SPEC `&2 * ln(&2)` SEQ_CONST)) THEN
  DISCH_THEN(MP_TAC o MATCH_MP SEQ_ADD) THEN
  SIMP_TAC[LN_DIV; PI_POS; REAL_OF_NUM_LT; ARITH; LN_MUL] THEN
  REWRITE_TAC[REAL_ARITH `&2 * l + p - l - (&4 * C - (C + C)) =
                          (l + p) - &2 * C`] THEN
  SIMP_TAC[REAL_ARITH `C = (l + p) / &2 <=> (l + p) - &2 * C = &0`] THEN
  MATCH_MP_TAC(REWRITE_RULE[TAUT `a /\ b ==> c <=> b ==> a ==> c`]
    SEQ_UNIQ) THEN
  REWRITE_TAC[GSYM REAL_OF_NUM_ADD; GSYM REAL_OF_NUM_MUL] THEN
  REWRITE_TAC[REAL_ARITH
   `a + (b + &4 * (c - x)) - ((d - &2 * x) + (e - (&2 * x + &1))) =
    (a + b + &4 * c + &1) - (d + e)`] THEN
  REWRITE_TAC[REAL_ARITH `&2 * l + &4 * n * l + &4 * (n + &1 / &2) * x + &1 =
                          (&4 * n + &2) * (l + x) + &1`] THEN
  ONCE_REWRITE_TAC[SEQ_SUC] THEN
  SIMP_TAC[GSYM LN_MUL; REAL_OF_NUM_LT; ARITH; LT_0] THEN
  REWRITE_TAC[GSYM SEQ_SUC] THEN
  CONV_TAC(LAND_CONV(GEN_ALPHA_CONV `n:num`)) THEN
  REWRITE_TAC[REAL_ARITH
   `((&4 * n + &2) * l + &1) - ((&2 * n + &1 / &2) * l + z) =
    (&2 * n + &3 / &2) * l + &1 - z`] THEN
  REWRITE_TAC[REAL_ARITH
   `(&2 * n + &3 / &2) * l + &1 - ((&2 * n + &1) + &1 / &2) * l' =
    (&2 * n + &3 / &2) * (l - l') + &1`] THEN
  GEN_REWRITE_TAC RAND_CONV [REAL_ARITH `&0 = -- &1 + &1`] THEN
  MATCH_MP_TAC SEQ_ADD THEN REWRITE_TAC[SEQ_CONST] THEN
  ONCE_REWRITE_TAC[REAL_ARITH `a * (b - c) = --(a * (c - b))`] THEN
  REWRITE_TAC[GSYM SEQ_NEG] THEN
  ONCE_REWRITE_TAC[SEQ_SUC] THEN
  SIMP_TAC[GSYM LN_DIV; REAL_LT_MUL; REAL_OF_NUM_LT; LT_0; ARITH;
           REAL_ARITH `&0 < &2 * &n + &1`] THEN
  SIMP_TAC[REAL_OF_NUM_LT; LT_0; REAL_FIELD
   `&0 < x ==> (&2 * x + &1) / (&2 * x) = &1 + &1 / (&2 * x)`] THEN
  REWRITE_TAC[GSYM SEQ_SUC] THEN
  CONV_TAC(LAND_CONV(GEN_ALPHA_CONV `n:num`)) THEN
  MP_TAC SEQ_SUBSEQ THEN REWRITE_TAC[RIGHT_IMP_FORALL_THM] THEN
  DISCH_THEN(MP_TAC o GENL [`f:num->real`; `l:real`] o
   SPECL [`f:num->real`; `l:real`; `2`; `0`]) THEN
  REWRITE_TAC[ADD_CLAUSES; ARITH; REAL_OF_NUM_MUL] THEN
  DISCH_THEN MATCH_MP_TAC THEN CONV_TAC(LAND_CONV(GEN_ALPHA_CONV `n:num`)) THEN
  REWRITE_TAC[REAL_ADD_RDISTRIB] THEN
  GEN_REWRITE_TAC RAND_CONV [REAL_ARITH `&1 = &1 + &3 / &2 * &0`] THEN
  MATCH_MP_TAC SEQ_ADD THEN REWRITE_TAC[LN_LIM_LEMMA] THEN
  MATCH_MP_TAC SEQ_MUL THEN REWRITE_TAC[SEQ_CONST] THEN
  MP_TAC LN_LIM_LEMMA THEN MP_TAC(SPEC `&1` SEQ_HARMONIC) THEN
  REWRITE_TAC[IMP_IMP] THEN DISCH_THEN(MP_TAC o MATCH_MP SEQ_MUL) THEN
  REWRITE_TAC[] THEN ONCE_REWRITE_TAC[SEQ_SUC] THEN
  SIMP_TAC[real_div; REAL_MUL_LID; REAL_MUL_ASSOC; REAL_MUL_LINV;
           REAL_MUL_RID; REAL_OF_NUM_EQ; NOT_SUC]);;