(* File: Stirling_Formula.thy Author: Manuel Eberl A proof of Stirling's formula, i.e. an asymptotic approximation of the Gamma function and the factorial. *) section ‹Stirling's Formula› theory Stirling_Formula imports "HOL-Analysis.Analysis" "Landau_Symbols.Landau_More" "HOL-Real_Asymp.Real_Asymp" begin context begin text ‹ First, we define the $S_n^*$ from Jameson's article: › qualified definition S' :: "nat ⇒ real ⇒ real" where "S' n x = 1/(2*x) + (∑r=1..<n. 1/(of_nat r+x)) + 1 /(2*(n + x))" text ‹ Next, the trapezium (also called $T$ in Jameson's article): › qualified definition T :: "real ⇒ real" where "T x = 1/(2*x) + 1/(2*(x+1))" text ‹ Now we define The difference $\Delta(x)$: › qualified definition D :: "real ⇒ real" where "D x = T x - ln (x + 1) + ln x" qualified lemma S'_telescope_trapezium: assumes "n > 0" shows "S' n x = (∑r<n. T (of_nat r + x))" proof (cases n) case (Suc m) hence m: "Suc m = n" by simp have "(∑r<n. T (of_nat r + x)) = (∑r<Suc m. 1 / (2 * real r + 2 * x)) + (∑r<n. 1 / (2 * real (Suc r) + 2 * x))" unfolding m by (simp add: T_def sum.distrib algebra_simps) also have "(∑r<Suc m. 1 / (2 * real r + 2 * x)) = 1/(2*x) + (∑r<m. 1 / (2 * real (Suc r) + 2 * x))" (is "_ = ?a + ?S") by (subst sum.lessThan_Suc_shift) simp also have "(∑r<n. 1 / (2 * real (Suc r) + 2 * x)) = ?S + 1 / (2*(real m + x + 1))" (is "_ = _ + ?b") by (simp add: Suc) also have "?a + ?S + (?S + ?b) = 2*?S + ?a + ?b" by (simp add: add_ac) also have "2 * ?S = (∑r=0..<m. 1 / (real (Suc r) + x))" unfolding sum_distrib_left by (intro sum.cong) (auto simp add: divide_simps) also have "(∑r=0..<m. 1 / (real (Suc r) + x)) = (∑r=Suc 0..<Suc m. 1 / (real r + x))" by (subst sum.atLeast_Suc_lessThan_Suc_shift) simp_all also have "… = (∑r=1..<n. 1 / (real r + x))" unfolding m by simp also have "… + ?a + ?b = S' n x" by (simp add: S'_def Suc) finally show ?thesis .. qed (insert assms, simp_all) qualified lemma stirling_trapezium: assumes x: "(x::real) > 0" shows "D x ∈ {0 .. 1/(12*x^2) - 1/(12 * (x+1)^2)}" proof - define y where "y = 1 / (2*x + 1)" from x have y: "y > 0" "y < 1" by (simp_all add: divide_simps y_def) from x have "D x = T x - ln ((x + 1) / x)" by (subst ln_div) (simp_all add: D_def) also from x have "(x + 1) / x = 1 + 1 / x" by (simp add: field_simps) finally have D: "D x = T x - ln (1 + 1/x)" . from y have "(λn. y * y^n) sums (y * (1 / (1 - y)))" by (intro geometric_sums sums_mult) simp_all hence "(λn. y ^ Suc n) sums (y / (1 - y))" by simp also from x have "y / (1 - y) = 1 / (2*x)" by (simp add: y_def divide_simps) finally have *: "(λn. y ^ Suc n) sums (1 / (2*x))" . from y have "(λn. (-y) * (-y)^n) sums ((-y) * (1 / (1 - (-y))))" by (intro geometric_sums sums_mult) simp_all hence "(λn. (-y) ^ Suc n) sums (-(y / (1 + y)))" by simp also from x have "y / (1 + y) = 1 / (2*(x+1))" by (simp add: y_def divide_simps) finally have **: "(λn. (-y) ^ Suc n) sums (-(1 / (2*(x+1))))" . from sums_diff[OF * **] have sum1: "(λn. y ^ Suc n - (-y) ^ Suc n) sums T x" by (simp add: T_def) from y have "abs y < 1" "abs (-y) < 1" by simp_all from sums_diff[OF this[THEN ln_series']] have "(λn. y ^ n / real n - (-y) ^ n / real n) sums (ln (1 + y) - ln (1 - y))" by simp also from y have "ln (1 + y) - ln (1 - y) = ln ((1 + y) / (1 - y))" by (simp add: ln_div) also from x have "(1 + y) / (1 - y) = 1 + 1/x" by (simp add: divide_simps y_def) finally have "(λn. y^n / real n - (-y)^n / real n) sums ln (1 + 1/x)" . hence sum2: "(λn. y^Suc n / real (Suc n) - (-y)^Suc n / real (Suc n)) sums ln (1 + 1/x)" by (subst sums_Suc_iff) simp have "ln (1 + 1/x) ≤ T x" proof (rule sums_le [OF _ sum2 sum1]) fix n :: nat show "y ^ Suc n / real (Suc n) - (-y) ^ Suc n / real (Suc n) ≤ y^Suc n - (-y) ^ Suc n" proof (cases "even n") case True hence eq: "A - (-y) ^ Suc n / B = A + y ^ Suc n / B" "A - (-y) ^ Suc n = A + y ^ Suc n" for A B by simp_all from y show ?thesis unfolding eq by (intro add_mono) (auto simp: divide_simps) qed simp_all qed hence "D x ≥ 0" by (simp add: D) define c where "c = (λn. if even n then 2 * (1 - 1 / real (Suc n)) else 0)" note sums_diff[OF sum1 sum2] also have "(λn. y ^ Suc n - (-y) ^ Suc n - (y ^ Suc n / real (Suc n) - (-y) ^ Suc n / real (Suc n))) = (λn. c n * y ^ Suc n)" by (intro ext) (simp add: c_def algebra_simps) finally have sum3: "(λn. c n * y ^ Suc n) sums D x" by (simp add: D) from y have "(λn. y^2 * (of_nat (Suc n) * y^n)) sums (y^2 * (1 / (1 - y)^2))" by (intro sums_mult geometric_deriv_sums) simp_all hence "(λn. of_nat (Suc n) * y^(n+2)) sums (y^2 / (1 - y)^2)" by (simp add: algebra_simps power2_eq_square) also from x have "y^2 / (1 - y)^2 = 1 / (4*x^2)" by (simp add: y_def divide_simps) finally have *: "(λn. real (Suc n) * y ^ (Suc (Suc n))) sums (1 / (4 * x⇧^{2}))" by simp from y have "(λn. y^2 * (of_nat (Suc n) * (-y)^n)) sums (y^2 * (1 / (1 - -y)^2))" by (intro sums_mult geometric_deriv_sums) simp_all hence "(λn. of_nat (Suc n) * (-y)^(n+2)) sums (y^2 / (1 + y)^2)" by (simp add: algebra_simps power2_eq_square) also from x have "y^2 / (1 + y)^2 = 1 / (2^2*(x+1)^2)" unfolding power_mult_distrib [symmetric] by (simp add: y_def divide_simps add_ac) finally have **: "(λn. real (Suc n) * (- y) ^ (Suc (Suc n))) sums (1 / (4 * (x + 1)⇧^{2}))" by simp define d where "d = (λn. if even n then 2 * real n else 0)" note sums_diff[OF * **] also have "(λn. real (Suc n) * y^(Suc (Suc n)) - real (Suc n) * (-y)^(Suc (Suc n))) = (λn. d (Suc n) * y ^ Suc (Suc n))" by (intro ext) (simp_all add: d_def) finally have "(λn. d n * y ^ Suc n) sums (1 / (4 * x⇧^{2}) - 1 / (4 * (x + 1)⇧^{2}))" by (subst (asm) sums_Suc_iff) (simp add: d_def) from sums_mult[OF this, of "1/3"] x have sum4: "(λn. d n / 3 * y ^ Suc n) sums (1 / (12 * x^2) - 1 / (12 * (x + 1)^2))" by (simp add: field_simps) have "D x ≤ (1 / (12 * x^2) - 1 / (12 * (x + 1)^2))" proof (intro sums_le [OF _ sum3 sum4] allI) fix n :: nat define c' :: "nat ⇒ real" where "c' = (λn. if odd n ∨ n = 0 then 0 else if n = 2 then 4/3 else 2)" show "c n * y ^ Suc n ≤ d n / 3 * y ^ Suc n" proof (intro mult_right_mono) have "c n ≤ c' n" by (simp add: c_def c'_def) also consider "n = 0" | "n = 1" | "n = 2" | "n ≥ 3" by force hence "c' n ≤ d n / 3" by cases (simp_all add: c'_def d_def) finally show "c n ≤ d n / 3" . qed (insert y, simp) qed with ‹D x ≥ 0› show ?thesis by simp qed text ‹ The following functions correspond to the $p_n(x)$ from the article. The special case $n = 0$ would not, strictly speaking, be necessary, but it allows some theorems to work even without the precondition $n \neq 0$: › qualified definition p :: "nat ⇒ real ⇒ real" where "p n x = (if n = 0 then 1/x else (∑r<n. D (real r + x)))" text ‹ We can write the Digamma function in terms of @{term S'}: › qualified lemma S'_LIMSEQ_Digamma: assumes x: "x ≠ 0" shows "(λn. ln (real n) - S' n x - 1/(2*x)) ⇢ Digamma x" proof - define c where "c = (λn. ln (real n) - (∑r<n. inverse (x + real r)))" have "eventually (λn. 1 / (2 * (x + real n)) = c n - (ln (real n) - S' n x - 1/(2*x))) at_top" using eventually_gt_at_top[of "0::nat"] proof eventually_elim fix n :: nat assume n: "n > 0" have "c n - (ln (real n) - S' n x - 1/(2*x)) = -(∑r<n. inverse (real r + x)) + (1/x + (∑r=1..<n. inverse (real r + x))) + 1/(2*(real n + x))" using x by (simp add: S'_def c_def field_simps) also have "1/x + (∑r=1..<n. inverse (real r + x)) = (∑r<n. inverse (real r + x))" unfolding lessThan_atLeast0 using n by (subst (2) sum.atLeast_Suc_lessThan) (simp_all add: field_simps) finally show "1 / (2 * (x + real n)) = c n - (ln (real n) - S' n x - 1/(2*x))" by simp qed moreover have "(λn. 1 / (2 * (x + real n))) ⇢ 0" by real_asymp ultimately have "(λn. c n - (ln (real n) - S' n x - 1/(2*x))) ⇢ 0" by (blast intro: Lim_transform_eventually) from tendsto_minus[OF this] have "(λn. (ln (real n) - S' n x - 1/(2*x)) - c n) ⇢ 0" by simp moreover from Digamma_LIMSEQ[OF x] have "c ⇢ Digamma x" by (simp add: c_def) ultimately show "(λn. ln (real n) - S' n x - 1/(2*x)) ⇢ Digamma x" by (rule Lim_transform [rotated]) qed text ‹ Moreover, we can give an expansion of @{term S'} with the @{term p} as variation terms. › qualified lemma S'_approx: "S' n x = ln (real n + x) - ln x + p n x" proof (cases "n = 0") case True thus ?thesis by (simp add: p_def S'_def) next case False hence "S' n x = (∑r<n. T (real r + x))" by (subst S'_telescope_trapezium) simp_all also have "… = (∑r<n. ln (real r + x + 1) - ln (real r + x) + D (real r + x))" by (simp add: D_def) also have "… = (∑r<n. ln (real (Suc r) + x) - ln (real r + x)) + p n x" using False by (simp add: sum.distrib add_ac p_def) also have "(∑r<n. ln (real (Suc r) + x) - ln (real r + x)) = ln (real n + x) - ln x" by (subst sum_lessThan_telescope) simp_all finally show ?thesis . qed text ‹ We define the limit of the @{term p} (simply called $p(x)$ in Jameson's article): › qualified definition P :: "real ⇒ real" where "P x = (∑n. D (real n + x))" qualified lemma D_summable: assumes x: "x > 0" shows "summable (λn. D (real n + x))" proof - have *: "summable (λn. 1 / (12 * (x + real n)⇧^{2}) - 1 / (12 * (x + real (Suc n))⇧^{2}))" by (rule telescope_summable') real_asymp show "summable (λn. D (real n + x))" proof (rule summable_comparison_test[OF _ *], rule exI[of _ 2], safe) fix n :: nat assume "n ≥ 2" show "norm (D (real n + x)) ≤ 1 / (12 * (x + real n)^2) - 1 / (12 * (x + real (Suc n))^2)" using stirling_trapezium[of "real n + x"] x by (auto simp: algebra_simps) qed qed qualified lemma p_LIMSEQ: assumes x: "x > 0" shows "(λn. p n x) ⇢ P x" proof (rule Lim_transform_eventually) from D_summable[OF x] have "(λn. D (real n + x)) sums P x" unfolding P_def by (simp add: sums_iff) then show "(λn. ∑r<n. D (real r + x)) ⇢ P x" by (simp add: sums_def) moreover from eventually_gt_at_top[of 1] show "eventually (λn. (∑r<n. D (real r + x)) = p n x) at_top" by eventually_elim (auto simp: p_def) qed text ‹ This gives us an expansion of the Digamma function: › lemma Digamma_approx: assumes x: "(x :: real) > 0" shows "Digamma x = ln x - 1 / (2 * x) - P x" proof - have "eventually (λn. ln (inverse (1 + x / real n)) + ln x - 1/(2*x) - p n x = ln (real n) - S' n x - 1/(2*x)) at_top" using eventually_gt_at_top[of "1::nat"] proof eventually_elim fix n :: nat assume n: "n > 1" have "ln (real n) - S' n x = ln ((real n) / (real n + x)) + ln x - p n x" using assms n unfolding S'_approx by (subst ln_div) (auto simp: algebra_simps) also from n have "real n / (real n + x) = inverse (1 + x / real n)" by (simp add: field_simps) finally show "ln (inverse (1 + x / real n)) + ln x - 1/(2*x) - p n x = ln (real n) - S' n x - 1/(2*x)" by simp qed moreover have "(λn. ln (inverse (1 + x / real n)) + ln x - 1/(2*x) - p n x) ⇢ ln (inverse (1 + 0)) + ln x - 1/(2*x) - P x" by (rule tendsto_intros p_LIMSEQ x real_tendsto_divide_at_top filterlim_real_sequentially | simp)+ hence "(λn. ln (inverse (1 + x / real n)) + ln x - 1/(2*x) - p n x) ⇢ ln x - 1/(2*x) - P x" by simp ultimately have "(λn. ln (real n) - S' n x - 1 / (2 * x)) ⇢ ln x - 1/(2*x) - P x" by (blast intro: Lim_transform_eventually) moreover from x have "(λn. ln (real n) - S' n x - 1 / (2 * x)) ⇢ Digamma x" by (intro S'_LIMSEQ_Digamma) simp_all ultimately show "Digamma x = ln x - 1 / (2 * x) - P x" by (rule LIMSEQ_unique [rotated]) qed text ‹ Next, we derive some bounds on @{term "P"}: › qualified lemma p_ge_0: "x > 0 ⟹ p n x ≥ 0" using stirling_trapezium[of "real n + x" for n] by (auto simp add: p_def intro!: sum_nonneg) qualified lemma P_ge_0: "x > 0 ⟹ P x ≥ 0" by (rule tendsto_lowerbound[OF p_LIMSEQ]) (insert p_ge_0[of x], simp_all) qualified lemma p_upper_bound: assumes "x > 0" "n > 0" shows "p n x ≤ 1/(12*x^2)" proof - from assms have "p n x = (∑r<n. D (real r + x))" by (simp add: p_def) also have "… ≤ (∑r<n. 1/(12*(real r + x)^2) - 1/(12 * (real (Suc r) + x)^2))" using stirling_trapezium[of "real r + x" for r] assms by (intro sum_mono) (simp add: add_ac) also have "… = 1 / (12 * x⇧^{2}) - 1 / (12 * (real n + x)⇧^{2})" by (subst sum_lessThan_telescope') simp also have "… ≤ 1 / (12 * x^2)" by simp finally show ?thesis . qed qualified lemma P_upper_bound: assumes "x > 0" shows "P x ≤ 1/(12*x^2)" proof (rule tendsto_upperbound) show "eventually (λn. p n x ≤ 1 / (12 * x^2)) at_top" using eventually_gt_at_top[of 0] by eventually_elim (use p_upper_bound[of x] assms in auto) show "(λn. p n x) ⇢ P x" by (simp add: assms p_LIMSEQ) qed auto text ‹ We can now use this approximation of the Digamma function to approximate @{term ln_Gamma} (since the former is the derivative of the latter). We therefore define the function $g$ from Jameson's article, which measures the difference between @{term ln_Gamma} and its approximation: › qualified definition g :: "real ⇒ real" where "g x = ln_Gamma x - (x - 1/2) * ln x + x" qualified lemma DERIV_g: "x > 0 ⟹ (g has_field_derivative -P x) (at x)" unfolding g_def [abs_def] using Digamma_approx[of x] by (auto intro!: derivative_eq_intros simp: field_simps) qualified lemma isCont_P: assumes "x > 0" shows "isCont P x" proof - define D' :: "real ⇒ real" where "D' = (λx. - 1 / (2 * x^2 * (x+1)^2))" have DERIV_D: "(D has_field_derivative D' x) (at x)" if "x > 0" for x unfolding D_def [abs_def] D'_def T_def by (insert that, (rule derivative_eq_intros refl | simp)+) (simp add: power2_eq_square divide_simps, (simp add: algebra_simps)?) note this [THEN DERIV_chain2, derivative_intros] have "(P has_field_derivative (∑n. D' (real n + x))) (at x)" unfolding P_def [abs_def] proof (rule has_field_derivative_series') show "convex {x/2<..}" by simp next fix n :: nat and y :: real assume y: "y ∈ {x/2<..}" with assms have "y > 0" by simp thus "((λa. D (real n + a)) has_real_derivative D' (real n + y)) (at y within {x/2<..})" by (auto intro!: derivative_eq_intros) next from assms D_summable[of x] show "summable (λn. D (real n + x))" by simp next show "uniformly_convergent_on {x/2<..} (λn x. ∑i<n. D' (real i + x))" proof (rule Weierstrass_m_test') fix n :: nat and y :: real assume y: "y ∈ {x/2<..}" with assms have "y > 0" by auto have "norm (D' (real n + y)) = (1 / (2 * (y + real n)^2)) * (1 / (y + real (Suc n))^2)" by (simp add: D'_def add_ac) also from y assms have "… ≤ (1 / (2 * (x/2)^2)) * (1 / (real (Suc n))^2)" by (intro mult_mono divide_left_mono power_mono) simp_all also have "1 / (real (Suc n))^2 = inverse ((real (Suc n))^2)" by (simp add: field_simps) finally show "norm (D' (real n + y)) ≤ (1 / (2 * (x/2)^2)) * inverse ((real (Suc n))^2)" . next show "summable (λn. (1 / (2 * (x/2)^2)) * inverse ((real (Suc n))^2))" by (subst summable_Suc_iff, intro summable_mult inverse_power_summable) simp_all qed qed (insert assms, simp_all add: interior_open) thus ?thesis by (rule DERIV_isCont) qed qualified lemma P_continuous_on [THEN continuous_on_subset]: "continuous_on {0<..} P" by (intro continuous_at_imp_continuous_on ballI isCont_P) auto qualified lemma P_integrable: assumes a: "a > 0" shows "P integrable_on {a..}" proof - define f where "f = (λn x. if x ∈ {a..real n} then P x else 0)" show "P integrable_on {a..}" proof (rule dominated_convergence) fix n :: nat from a have "P integrable_on {a..real n}" by (intro integrable_continuous_real P_continuous_on) auto hence "f n integrable_on {a..real n}" by (rule integrable_eq) (simp add: f_def) thus "f n integrable_on {a..}" by (rule integrable_on_superset) (auto simp: f_def) next fix n :: nat show "norm (f n x) ≤ of_real (1/12) * (1 / x^2)" if "x ∈ {a..}" for x using a P_ge_0 P_upper_bound by (auto simp: f_def) next show "(λx::real. of_real (1/12) * (1 / x^2)) integrable_on {a..}" using has_integral_inverse_power_to_inf[of 2 a] a by (intro integrable_on_cmult_left) auto next show "(λn. f n x) ⇢ P x" if "x∈{a..}" for x proof - have "eventually (λn. real n ≥ x) at_top" using filterlim_real_sequentially by (simp add: filterlim_at_top) with that have "eventually (λn. f n x = P x) at_top" by (auto elim!: eventually_mono simp: f_def) thus "(λn. f n x) ⇢ P x" by (simp add: tendsto_eventually) qed qed qed qualified definition c :: real where "c = integral {1..} (λx. -P x) + g 1" text ‹ We can now give bounds on @{term g}: › qualified lemma g_bounds: assumes x: "x ≥ 1" shows "g x ∈ {c..c + 1/(12*x)}" proof - from assms have int_nonneg: "integral {x..} P ≥ 0" by (intro Henstock_Kurzweil_Integration.integral_nonneg P_integrable) (auto simp: P_ge_0) have int_upper_bound: "integral {x..} P ≤ 1/(12*x)" proof (rule has_integral_le) from x show "(P has_integral integral {x..} P) {x..}" by (intro integrable_integral P_integrable) simp_all from x has_integral_mult_right[OF has_integral_inverse_power_to_inf[of 2 x], of "1/12"] show "((λx. 1/(12*x^2)) has_integral (1/(12*x))) {x..}" by (simp add: field_simps) qed (insert P_upper_bound x, simp_all) note DERIV_g [THEN DERIV_chain2, derivative_intros] from assms have int1: "((λx. -P x) has_integral (g x - g 1)) {1..x}" by (intro fundamental_theorem_of_calculus) (auto simp: has_real_derivative_iff_has_vector_derivative [symmetric] intro!: derivative_eq_intros) from x have int2: "((λx. -P x) has_integral integral {x..} (λx. -P x)) {x..}" by (intro integrable_integral integrable_neg P_integrable) simp_all from has_integral_Un[OF int1 int2] x have "((λx. - P x) has_integral g x - g 1 - integral {x..} P) ({1..x} ∪ {x..})" by (simp add: max_def) also from x have "{1..x} ∪ {x..} = {1..}" by auto finally have "((λx. -P x) has_integral g x - g 1 - integral {x..} P) {1..}" . moreover have "((λx. -P x) has_integral integral {1..} (λx. -P x)) {1..}" by (intro integrable_integral integrable_neg P_integrable) simp_all ultimately have "g x - g 1 - integral {x..} P = integral {1..} (λx. -P x)" by (simp add: has_integral_unique) hence "g x = c + integral {x..} P" by (simp add: c_def algebra_simps) with int_upper_bound int_nonneg show "g x ∈ {c..c + 1/(12*x)}" by simp qed text ‹ Finally, we have bounds on @{term ln_Gamma}, @{term Gamma}, and @{term fact}. › qualified lemma ln_Gamma_bounds_aux: "x ≥ 1 ⟹ ln_Gamma x ≥ c + (x - 1/2) * ln x - x" "x ≥ 1 ⟹ ln_Gamma x ≤ c + (x - 1/2) * ln x - x + 1/(12*x)" using g_bounds[of x] by (simp_all add: g_def) qualified lemma Gamma_bounds_aux: assumes x: "x ≥ 1" shows "Gamma x ≥ exp c * x powr (x - 1/2) / exp x" "Gamma x ≤ exp c * x powr (x - 1/2) / exp x * exp (1/(12*x))" proof - have "exp (ln_Gamma x) ≥ exp (c + (x - 1/2) * ln x - x)" by (subst exp_le_cancel_iff, rule ln_Gamma_bounds_aux) (simp add: x) with x show "Gamma x ≥ exp c * x powr (x - 1/2) / exp x" by (simp add: Gamma_real_pos_exp exp_add exp_diff powr_def del: exp_le_cancel_iff) next have "exp (ln_Gamma x) ≤ exp (c + (x - 1/2) * ln x - x + 1/(12*x))" by (subst exp_le_cancel_iff, rule ln_Gamma_bounds_aux) (simp add: x) with x show "Gamma x ≤ exp c * x powr (x - 1/2) / exp x * exp (1/(12*x))" by (simp add: Gamma_real_pos_exp exp_add exp_diff powr_def del: exp_le_cancel_iff) qed qualified lemma Gamma_asymp_equiv_aux: "Gamma ∼[at_top] (λx. exp c * x powr (x - 1/2) / exp x)" proof (rule asymp_equiv_sandwich) include asymp_equiv_notation show "eventually (λx. exp c * x powr (x - 1/2) / exp x ≤ Gamma x) at_top" "eventually (λx. exp c * x powr (x - 1/2) / exp x * exp (1/(12*x)) ≥ Gamma x) at_top" using eventually_ge_at_top[of "1::real"] by (eventually_elim; use Gamma_bounds_aux in force)+ have "((λx::real. exp (1 / (12 * x))) ⤏ exp 0) at_top" by real_asymp hence "(λx. exp (1 / (12 * x))) ∼ (λx. 1 :: real)" by (intro asymp_equivI') simp_all hence "(λx. exp c * x powr (x - 1 / 2) / exp x * 1) ∼ (λx. exp c * x powr (x - 1 / 2) / exp x * exp (1 / (12 * x)))" by (intro asymp_equiv_mult asymp_equiv_refl) (simp add: asymp_equiv_sym) thus "(λx. exp c * x powr (x - 1 / 2) / exp x) ∼ (λx. exp c * x powr (x - 1 / 2) / exp x * exp (1 / (12 * x)))" by simp qed simp_all qualified lemma exp_1_powr_real [simp]: "exp (1::real) powr x = exp x" by (simp add: powr_def) qualified lemma fact_asymp_equiv_aux: "fact ∼[at_top] (λx. exp c * sqrt (real x) * (real x / exp 1) powr real x)" proof - include asymp_equiv_notation have "fact ∼ (λn. Gamma (real (Suc n)))" by (simp add: Gamma_fact) also have "eventually (λn. Gamma (real (Suc n)) = real n * Gamma (real n)) at_top" using eventually_gt_at_top[of "0::nat"] by eventually_elim (insert Gamma_plus1[of "real n" for n], auto simp: add_ac of_nat_in_nonpos_Ints_iff) also have "(λn. Gamma (real n)) ∼ (λn. exp c * real n powr (real n - 1/2) / exp (real n))" by (rule asymp_equiv_compose'[OF Gamma_asymp_equiv_aux] filterlim_real_sequentially)+ also have "eventually (λn. real n * (exp c * real n powr (real n - 1 / 2) / exp (real n)) = exp c * sqrt (real n) * (real n / exp 1) powr real n) at_top" using eventually_gt_at_top[of "0::nat"] proof eventually_elim fix n :: nat assume n: "n > 0" thus "real n * (exp c * real n powr (real n - 1 / 2) / exp (real n)) = exp c * sqrt (real n) * (real n / exp 1) powr real n" by (subst powr_diff) (simp_all add: powr_divide powr_half_sqrt field_simps) qed finally show ?thesis by - (simp_all add: asymp_equiv_mult) qed text ‹ We cal also bound @{term Digamma} above and below. › lemma Digamma_plus_1_gt_ln: assumes x: "x > (0 :: real)" shows "Digamma (x + 1) > ln x" proof - have "0 < (17 :: real)" by simp also have "17 ≤ 12 * x ^ 2 + 28 * x + 17" using x by auto finally have "0 < (12 * x ^ 2 + 28 * x + 17) / (12 * (x + 1) ^ 2 * (1 + 2 * x))" using x by (intro divide_pos_pos mult_pos_pos zero_less_power) auto also have "… = 2 / (2 * x + 1) - 1 / (2 * (x + 1)) - 1 / (12 * (x + 1) ^ 2)" using x by (simp add: divide_simps) (auto simp: field_simps power2_eq_square add_eq_0_iff) also have "2 / (2 * x + 1) ≤ ln (x + 1) - ln x" using ln_inverse_approx_ge[of x "x + 1"] x by simp also have "… - 1 / (2 * (x + 1)) - 1 / (12 * (x + 1) ^ 2) ≤ ln (x + 1) - ln x - 1 / (2 * (x + 1)) - P (x + 1)" using P_upper_bound[of "x + 1"] x by (intro diff_mono) auto also have "… = Digamma (x + 1) - ln x" by (subst Digamma_approx) (use x in auto) finally show "Digamma (x + 1) > ln x" by simp qed lemma Digamma_less_ln: assumes x: "x > (0 :: real)" shows "Digamma x < ln x" proof - have "Digamma x - ln x = - (1 / (2 * x)) - P x" by (subst Digamma_approx) (use x in auto) also have "… < 0 - P x" using x by (intro diff_strict_right_mono) auto also have "… ≤ 0" using P_ge_0[of x] x by simp finally show "Digamma x < ln x" by simp qed text ‹ We still need to determine the constant term @{term c}, which we do using Wallis' product formula: $$\prod_{n=1}^\infty \frac{4n^2}{4n^2-1} = \frac{\pi}{2}$$ › qualified lemma powr_mult_2: "(x::real) > 0 ⟹ x powr (y * 2) = (x^2) powr y" by (subst mult.commute, subst powr_powr [symmetric]) (simp add: powr_numeral) qualified lemma exp_mult_2: "exp (y * 2 :: real) = exp y * exp y" by (subst exp_add [symmetric]) simp qualified lemma exp_c: "exp c = sqrt (2*pi)" proof - include asymp_equiv_notation define p where "p = (λn. ∏k=1..n. (4*real k^2) / (4*real k^2 - 1))" have p_0 [simp]: "p 0 = 1" by (simp add: p_def) have p_Suc: "p (Suc n) = p n * (4 * real (Suc n)^2) / (4 * real (Suc n)^2 - 1)" for n unfolding p_def by (subst prod.nat_ivl_Suc') simp_all have p: "p = (λn. 16 ^ n * fact n ^ 4 / (fact (2 * n))⇧^{2}/ (2 * real n + 1))" proof fix n :: nat have "p n = (∏k=1..n. (2*real k)^2 / (2*real k - 1) / (2 * real k + 1))" unfolding p_def by (intro prod.cong refl) (simp add: field_simps power2_eq_square) also have "… = (∏k=1..n. (2*real k)^2 / (2*real k - 1)) / (∏k=1..n. (2 * real (Suc k) - 1))" by (simp add: prod_dividef prod.distrib add_ac) also have "(∏k=1..n. (2 * real (Suc k) - 1)) = (∏k=Suc 1..Suc n. (2 * real k - 1))" by (subst prod.atLeast_Suc_atMost_Suc_shift) simp_all also have "… = (∏k=1..n. (2 * real k - 1)) * (2 * real n + 1)" by (induction n) (simp_all add: prod.nat_ivl_Suc') also have "(∏k = 1..n. (2 * real k)⇧^{2}/ (2 * real k - 1)) / … = (∏k = 1..n. (2 * real k)^2 / (2 * real k - 1)^2) / (2 * real n + 1)" unfolding power2_eq_square by (simp add: prod.distrib prod_dividef) also have "(∏k = 1..n. (2 * real k)^2 / (2 * real k - 1)^2) = (∏k = 1..n. (2 * real k)^4 / ((2*real k)*(2 * real k - 1))^2)" by (rule prod.cong) (simp_all add: power2_eq_square eval_nat_numeral) also have "… = 16^n * fact n^4 / (∏k=1..n. (2*real k) * (2*real k - 1))^2" by (simp add: prod.distrib prod_dividef fact_prod prod_power_distrib [symmetric] prod_constant) also have "(∏k=1..n. (2*real k) * (2*real k - 1)) = fact (2*n)" by (induction n) (simp_all add: algebra_simps prod.nat_ivl_Suc') finally show "p n = 16 ^ n * fact n ^ 4 / (fact (2 * n))⇧^{2}/ (2 * real n + 1)" . qed have "p ∼ (λn. 16 ^ n * fact n ^ 4 / (fact (2 * n))⇧^{2}/ (2 * real n + 1))" by (simp add: p) also have "… ∼ (λn. 16^n * (exp c * sqrt (real n) * (real n / exp 1) powr real n)^4 / (exp c * sqrt (real (2*n)) * (real (2*n) / exp 1) powr real (2*n))^2 / (2 * real n + 1))" (is "_ ∼ ?f") by (intro asymp_equiv_mult asymp_equiv_divide asymp_equiv_refl mult_nat_left_at_top fact_asymp_equiv_aux asymp_equiv_power asymp_equiv_compose'[OF fact_asymp_equiv_aux]) simp_all also have "eventually (λn. … n = exp c ^ 2 / (4 + 2/n)) at_top" using eventually_gt_at_top[of "0::nat"] proof eventually_elim fix n :: nat assume n: "n > 0" have [simp]: "16^n = 4^n * (4^n :: real)" by (simp add: power_mult_distrib [symmetric]) from n have "?f n = exp c ^ 2 * (n / (2*(2*n+1)))" by (simp add: power_mult_distrib divide_simps powr_mult real_sqrt_power_even) (simp add: field_simps power2_eq_square eval_nat_numeral powr_mult_2 exp_mult_2 powr_realpow) also from n have "… = exp c ^ 2 / (4 + 2/n)" by (simp add: field_simps) finally show "?f n = …" . qed also have "(λx. 4 + 2 / real x) ∼ (λx. 4)" by (subst asymp_equiv_add_right) auto finally have "p ⇢ exp c ^ 2 / 4" by (rule asymp_equivD_const) (simp_all add: asymp_equiv_divide) moreover have "p ⇢ pi / 2" unfolding p_def by (rule wallis) ultimately have "exp c ^ 2 / 4 = pi / 2" by (rule LIMSEQ_unique) hence "2 * pi = exp c ^ 2" by simp also have "sqrt (exp c ^ 2) = exp c" by simp finally show "exp c = sqrt (2 * pi)" .. qed qualified lemma c: "c = ln (2*pi) / 2" proof - note exp_c [symmetric] also have "ln (exp c) = c" by simp finally show ?thesis by (simp add: ln_sqrt) qed text ‹ This gives us the final bounds: › theorem Gamma_bounds: assumes "x ≥ 1" shows "Gamma x ≥ sqrt (2*pi/x) * (x / exp 1) powr x" (is ?th1) "Gamma x ≤ sqrt (2*pi/x) * (x / exp 1) powr x * exp (1 / (12 * x))" (is ?th2) proof - from assms have "exp c * x powr (x - 1/2) / exp x = sqrt (2*pi/x) * (x / exp 1) powr x" by (subst powr_diff) (simp add: exp_c real_sqrt_divide powr_divide powr_half_sqrt field_simps) with Gamma_bounds_aux[OF assms] show ?th1 ?th2 by simp_all qed theorem ln_Gamma_bounds: assumes "x ≥ 1" shows "ln_Gamma x ≥ ln (2*pi/x) / 2 + x * ln x - x" (is ?th1) "ln_Gamma x ≤ ln (2*pi/x) / 2 + x * ln x - x + 1/(12*x)" (is ?th2) proof - from ln_Gamma_bounds_aux[OF assms] assms show ?th1 ?th2 by (simp_all add: c field_simps ln_div) from assms have "exp c * x powr (x - 1/2) / exp x = sqrt (2*pi/x) * (x / exp 1) powr x" by (subst powr_diff) (simp add: exp_c real_sqrt_divide powr_divide powr_half_sqrt field_simps) qed theorem fact_bounds: assumes "n > 0" shows "(fact n :: real) ≥ sqrt (2*pi*n) * (n / exp 1) ^ n" (is ?th1) "(fact n :: real) ≤ sqrt (2*pi*n) * (n / exp 1) ^ n * exp (1 / (12 * n))" (is ?th2) proof - from assms have n: "real n ≥ 1" by simp from assms Gamma_plus1[of "real n"] have "real n * Gamma (real n) = Gamma (real (Suc n))" by (simp add: of_nat_in_nonpos_Ints_iff add_ac) also have "Gamma (real (Suc n)) = fact n" by (subst Gamma_fact [symmetric]) simp finally have *: "fact n = real n * Gamma (real n)" by simp have "2*pi/n = 2*pi*n / n^2" by (simp add: power2_eq_square) also have "sqrt … = sqrt (2*pi*n) / n" by (subst real_sqrt_divide) simp_all also have "real n * … = sqrt (2*pi*n)" by simp finally have **: "real n * sqrt (2*pi/real n) = sqrt (2*pi*real n)" . note * also note Gamma_bounds(2)[OF n] also have "real n * (sqrt (2 * pi / real n) * (real n / exp 1) powr real n * exp (1 / (12 * real n))) = (real n * sqrt (2*pi/n)) * (n / exp 1) powr n * exp (1 / (12 * n))" by (simp add: algebra_simps) also from n have "(real n / exp 1) powr real n = (real n / exp 1) ^ n" by (subst powr_realpow) simp_all also note ** finally show ?th2 by - (insert assms, simp_all) have "sqrt (2*pi*n) * (n / exp 1) powr n = n * (sqrt (2*pi/n) * (n / exp 1) powr n)" by (subst ** [symmetric]) (simp add: field_simps) also from assms have "… ≤ real n * Gamma (real n)" by (intro mult_left_mono Gamma_bounds(1)) simp_all also from n have "(real n / exp 1) powr real n = (real n / exp 1) ^ n" by (subst powr_realpow) simp_all also note * [symmetric] finally show ?th1 . qed theorem ln_fact_bounds: assumes "n > 0" shows "ln (fact n :: real) ≥ ln (2*pi*n)/2 + n * ln n - n" (is ?th1) "ln (fact n :: real) ≤ ln (2*pi*n)/2 + n * ln n - n + 1/(12*real n)" (is ?th2) proof - have "ln (fact n :: real) ≥ ln (sqrt (2*pi*real n)*(real n/exp 1)^n)" using fact_bounds(1)[OF assms] assms by (subst ln_le_cancel_iff) auto also from assms have "ln (sqrt (2*pi*real n)*(real n/exp 1)^n) = ln (2*pi*n)/2 + n * ln n - n" by (simp add: ln_mult ln_div ln_realpow ln_sqrt field_simps) finally show ?th1 . next have "ln (fact n :: real) ≤ ln (sqrt (2*pi*real n) * (real n/exp 1)^n * exp (1/(12*real n)))" using fact_bounds(2)[OF assms] assms by (subst ln_le_cancel_iff) auto also from assms have "… = ln (2*pi*n)/2 + n * ln n - n + 1/(12*real n)" by (simp add: ln_mult ln_div ln_realpow ln_sqrt field_simps) finally show ?th2 . qed theorem Gamma_asymp_equiv: "Gamma ∼[at_top] (λx. sqrt (2*pi/x) * (x / exp 1) powr x :: real)" proof - note Gamma_asymp_equiv_aux also have "eventually (λx. exp c * x powr (x - 1/2) / exp x = sqrt (2*pi/x) * (x / exp 1) powr x) at_top" using eventually_gt_at_top[of "0::real"] proof eventually_elim fix x :: real assume x: "x > 0" thus "exp c * x powr (x - 1/2) / exp x = sqrt (2*pi/x) * (x / exp 1) powr x" by (subst powr_diff) (simp add: exp_c powr_half_sqrt powr_divide field_simps real_sqrt_divide) qed finally show ?thesis . qed theorem fact_asymp_equiv: "fact ∼[at_top] (λn. sqrt (2*pi*n) * (n / exp 1) ^ n :: real)" proof - note fact_asymp_equiv_aux also have "eventually (λn. exp c * sqrt (real n) = sqrt (2 * pi * real n)) at_top" using eventually_gt_at_top[of "0::nat"] by eventually_elim (simp add: exp_c real_sqrt_mult) also have "eventually (λn. (n / exp 1) powr n = (n / exp 1) ^ n) at_top" using eventually_gt_at_top[of "0::nat"] by eventually_elim (simp add: powr_realpow) finally show ?thesis . qed corollary stirling_tendsto_sqrt_pi: "(λn. fact n / (sqrt (2 * n) * (n / exp 1) ^ n)) ⇢ sqrt pi" proof - have *: "(λn. fact n / (sqrt (2 * pi * n) * (n / exp 1) ^ n)) ⇢ 1" using fact_asymp_equiv by (simp add: asymp_equiv_def) have "(λn. sqrt pi * fact n / (sqrt (2 * pi * real n) * (real n / exp 1) ^ n)) = (λn. fact n / (sqrt (real (2 * n)) * (real n / exp 1) ^ n))" by (force simp add: divide_simps powr_realpow real_sqrt_mult) with tendsto_mult_left[OF *, of "sqrt pi"] show ?thesis by simp qed end end