# 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 nC × (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 nsqrt (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 nl × 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.