Library Coqtail.Reals.RStirling

Proof of Stirling equivalence of factorial.

Require Import Reals.
Require Import Rsequence_facts.
Require Import Rseries_def.
Require Import Rseries_facts.
Require Import Rseries_usual.
Require Import Rpser.
Require Import Rintegral.
Require Import RTaylor.
Require Import Fourier.
Require Import Wallis.

Partial result : De Moivre approximation.

Section De_Moivre.

Lemma De_Moivre_equiv : {C | Rseq_fact ¬ (fun n ⇒ C × (INR n ) ^ n × exp (- (INR n )) × sqrt (INR n)) & 0 < C}.
Proof.
destruct Un_cv as [l Hl Hp].
∃ (/ l).
apply Rseq_equiv_sym.
apply Rseq_equiv_1.
intros n; apply not_0_INR.
apply fact_neq_0.
apply Rseq_equiv_sym.
intros eps Heps.
destruct (Hl (eps × / 2 × l)) as [N HN].
  apply Rmult_lt_0_compat; [fourier|assumption].
∃ N; intros n Hn.
unfold Rseq_constant, Rseq_minus, Rseq_div.
rewrite (Rabs_right 1); [|fourier].
rewrite Rmult_1_r.
rewrite Rabs_minus_sym.
apply Rle_trans with (eps / 2); [left|fourier].
apply (Rmult_lt_reg_l l); [assumption|].
pattern l at 1; rewrite <- Rabs_right; [|fourier].
rewrite <- Rabs_mult.
rewrite (Rmult_comm l (eps / 2)).
replace (l × (/ l × INR n ^ n × exp (- INR n) × sqrt (INR n) / Rseq_fact n - 1))
  with (INR n ^ n × exp (- INR n) × sqrt (INR n) / Rseq_fact n - l).
apply HN; assumption.
field; split.
  intro; fourier.
  apply not_0_INR; apply fact_neq_0.
apply Rinv_0_lt_compat; assumption.
Qed.

End De_Moivre.

Final result : Stirling approximation.

Section Stirling.

Local Coercion INR : nat >-> R.

Lemma Stirling_equiv : Rseq_fact ¬ (fun n ⇒ sqrt (2 × PI) × (INR n ) ^ n × exp (- (INR n )) × sqrt (INR n)).
Proof.
destruct De_Moivre_equiv as [l Hl].
assert(l^2/(2×PI ) = 1) as Heq.
eapply Rseq_cv_unique.
apply Wallis_quotient_lim2.
  auto with ×.
  assert (Hrw : (fun n ⇒ l × n ^ n × exp (-n) × sqrt n ) == (fun n ⇒ ((n / exp 1) ^ n × sqrt n × l))).
   intro n.
   rewrite exp_Ropp.
   replace (exp n) with (exp (R1 ×n)) by auto with ×.
   rewrite <- (exp_pow 1 n).
   unfold Rdiv; rewrite Rpow_mult_distr, Rinv_pow.
    field.
    assert (H := exp_pos 1); auto with ×.
  eapply Rseq_equiv_eq_compat.
   apply Rseq_eq_refl.
   apply Hrw.

   trivial.

apply Wallis_quotient_lim1.

replace l with (sqrt(2×PI)) in Hl.
apply Hl.

assert(HPI := PI_RGT_0).
apply sqrt_lem_1.
fourier.
auto with ×.
replace (l×l) with ((l^2 / (2×PI ))* (2×PI )).
rewrite Heq; ring.
unfold Rdiv, Rsqr; field; auto with ×.
Qed.

End Stirling.