Session Stirling_Formula

Theory Stirling_Formula

(*
  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"
begin

context
begin

text ‹
  First, we define the $S_n^*$ from Jameson's article:
›
private 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):
›
private definition T :: "real  real" where
  "T x = 1/(2*x) + 1/(2*(x+1))"

text ‹
  Now we define The difference $\Delta(x)$:
›
private definition D :: "real  real" where
  "D x = T x - ln (x + 1) + ln x"

private 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)

private 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 * x2))" 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 * x2) - 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$:
›
private 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'}:
›
private 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 (rule real_tendsto_divide_at_top tendsto_const filterlim_tendsto_pos_mult_at_top
          filterlim_tendsto_add_at_top filterlim_real_sequentially | simp)+
  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.
›
private 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):
›
private definition P :: "real  real" where
  "P x = (n. D (real n + x))"

private 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_tendsto_divide_at_top tendsto_const 
             filterlim_tendsto_pos_mult_at_top filterlim_pow_at_top
             filterlim_tendsto_add_at_top filterlim_real_sequentially | simp)+
  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

private 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"}:
›
private 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)

private lemma P_ge_0: "x > 0  P x  0"
  by (rule tendsto_lowerbound[OF p_LIMSEQ]) 
     (insert p_ge_0[of x], simp_all)

private 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 * x2) - 1 / (12 * (real n + x)2)"
    by (subst sum_lessThan_telescope') simp
  also have "  1 / (12 * x^2)" by simp
  finally show ?thesis .
qed

private 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:
›

private definition g :: "real  real" where
  "g x = ln_Gamma x - (x - 1/2) * ln x + x"

private 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)

private 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

private lemma P_continuous_on [THEN continuous_on_subset]: "continuous_on {0<..} P"
  by (intro continuous_at_imp_continuous_on ballI isCont_P) auto

private 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

private definition c :: real where "c = integral {1..} (λx. -P x) + g 1"

text ‹
  We can now give bounds on @{term g}:
›
private 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_field_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}.
›
private 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)

private 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

private 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 (rule tendsto_intros real_tendsto_divide_at_top filterlim_tendsto_pos_mult_at_top)+
       (simp_all add: filterlim_ident)
  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

private lemma exp_1_powr_real [simp]: "exp (1::real) powr x = exp x"
  by (simp add: powr_def)

private 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 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}$$
›
private 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)

private lemma exp_mult_2: "exp (y * 2 :: real) = exp y * exp y"
  by (subst exp_add [symmetric]) simp

private 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

private 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

end

end

Theory Gamma_Asymptotics

(*
  File:    Gamma_Asymptotics.thy
  Author:  Manuel Eberl

  The complete asymptotics of the real and complex logarithmic Gamma functions.
  Also of the real Polygamma functions (could be extended to the complex ones fairly easily
  if needed).
*)
section ‹Complete asymptotics of the logarithmic Gamma function›
theory Gamma_Asymptotics
imports 
  "HOL-Complex_Analysis.Complex_Analysis"
  "HOL-Real_Asymp.Real_Asymp"
  Bernoulli.Bernoulli_FPS 
  Bernoulli.Periodic_Bernpoly 
  Stirling_Formula
begin

subsection ‹Auxiliary Facts›
  
(* TODO Move *)
lemma arg_of_real [simp]:
  "x > 0  arg (complex_of_real x) = 0"
  "x < 0  arg (complex_of_real x) = pi"
  by (rule arg_unique; simp add: complex_sgn_def scaleR_conv_of_real)+
  
lemma arg_conv_arctan:
  assumes "Re z > 0"
  shows   "arg z = arctan (Im z / Re z)"
proof (rule arg_unique)
  show "sgn z = cis (arctan (Im z / Re z))"
  proof (rule complex_eqI)
    have "Re (cis (arctan (Im z / Re z))) = 1 / sqrt (1 + (Im z)2 / (Re z)2)"
      by (simp add: cos_arctan power_divide)
    also have "1 + Im z ^ 2 / Re z ^ 2 = norm z ^ 2 / Re z ^ 2"
      using assms by (simp add: cmod_def field_simps)
    also have "1 / sqrt  = Re z / norm z"
      using assms by (simp add: real_sqrt_divide)
    finally show "Re (sgn z) = Re (cis (arctan (Im z / Re z)))"
      by simp
  next
    have "Im (cis (arctan (Im z / Re z))) = Im z / (Re z * sqrt (1 + (Im z)2 / (Re z)2))"
      by (simp add: sin_arctan field_simps)
    also have "1 + Im z ^ 2 / Re z ^ 2 = norm z ^ 2 / Re z ^ 2"
      using assms by (simp add: cmod_def field_simps)
    also have "Im z / (Re z * sqrt ) = Im z / norm z"
      using assms by (simp add: real_sqrt_divide)
    finally show "Im (sgn z) = Im (cis (arctan (Im z / Re z)))"
      by simp
  qed
next
  show "arctan (Im z / Re z) > -pi"
    by (rule le_less_trans[OF _ arctan_lbound]) auto
next
  have "arctan (Im z / Re z) < pi / 2"
    by (rule arctan_ubound)
  also have "  pi" by simp
  finally show "arctan (Im z / Re z)  pi"
    by simp
qed

lemma mult_indicator_cong:
  fixes f g :: "_  'a :: semiring_1"
  shows "(x. x  A  f x = g x)  indicator A x * f x = indicator A x * g x"
  by (auto simp: indicator_def)
  
lemma has_absolute_integral_change_of_variables_1':
  fixes f :: "real  real" and g :: "real  real"
  assumes S: "S  sets lebesgue"
    and der_g: "x. x  S  (g has_field_derivative g' x) (at x within S)"
    and inj: "inj_on g S"
  shows "(λx. ¦g' x¦ *R f(g x)) absolutely_integrable_on S 
           integral S (λx. ¦g' x¦ *R f(g x)) = b
      f absolutely_integrable_on (g ` S)  integral (g ` S) f = b"
proof -
  have "(λx. ¦g' x¦ *R vec (f(g x)) :: real ^ 1) absolutely_integrable_on S 
           integral S (λx. ¦g' x¦ *R vec (f(g x))) = (vec b :: real ^ 1)
          (λx. vec (f x) :: real ^ 1) absolutely_integrable_on (g ` S) 
           integral (g ` S) (λx. vec (f x)) = (vec b :: real ^ 1)"
    using assms unfolding has_field_derivative_iff_has_vector_derivative
    by (intro has_absolute_integral_change_of_variables_1 assms) auto
  thus ?thesis
    by (simp add: absolutely_integrable_on_1_iff integral_on_1_eq)
qed

corollary Ln_times_of_nat:
    "r > 0; z  0  Ln(of_nat r * z :: complex) = ln (of_nat r) + Ln(z)"
  using Ln_times_of_real[of "of_nat r" z] by simp

lemma tendsto_of_real_0_I: 
  "(f  0) G  ((λx. (of_real (f x)))  (0 ::'a::real_normed_div_algebra)) G"
  by (subst (asm) tendsto_of_real_iff [symmetric]) simp  

lemma negligible_atLeastAtMostI: "b  a  negligible {a..(b::real)}"
  by (cases "b < a") auto
    
lemma vector_derivative_cong_eq:
  assumes "eventually (λx. x  A  f x = g x) (nhds x)" "x = y" "A = B" "x  A"
  shows   "vector_derivative f (at x within A) = vector_derivative g (at y within B)"
proof -
  from eventually_nhds_x_imp_x[OF assms(1)] assms(4) have "f x = g x" by blast
  hence "(λD. (f has_vector_derivative D) (at x within A)) = 
           (λD. (g has_vector_derivative D) (at x within A))" using assms
    by (intro ext has_vector_derivative_cong_ev refl assms) simp_all
  thus ?thesis by (simp add: vector_derivative_def assms)
qed

lemma differentiable_of_real [simp]: "of_real differentiable at x within A"
proof -
  have "(of_real has_vector_derivative 1) (at x within A)"
    by (auto intro!: derivative_eq_intros)
  thus ?thesis by (rule differentiableI_vector)
qed
  
lemma higher_deriv_cong_ev:
  assumes "eventually (λx. f x = g x) (nhds x)" "x = y"
  shows   "(deriv ^^ n) f x = (deriv ^^ n) g y"
proof -
  from assms(1) have "eventually (λx. (deriv ^^ n) f x = (deriv ^^ n) g x) (nhds x)"
  proof (induction n arbitrary: f g)
    case (Suc n)
    from Suc.prems have "eventually (λy. eventually (λz. f z = g z) (nhds y)) (nhds x)"
      by (simp add: eventually_eventually)
    hence "eventually (λx. deriv f x = deriv g x) (nhds x)"
      by eventually_elim (rule deriv_cong_ev, simp_all)
    thus ?case by (auto intro!: deriv_cong_ev Suc simp: funpow_Suc_right simp del: funpow.simps)
  qed auto
  from eventually_nhds_x_imp_x[OF this] assms(2) show ?thesis by simp
qed

lemma deriv_of_real [simp]: 
  "at x within A  bot  vector_derivative of_real (at x within A) = 1"
  by (auto intro!: vector_derivative_within derivative_eq_intros)

lemma deriv_Re [simp]: "deriv Re = (λ_. 1)"
  by (auto intro!: DERIV_imp_deriv simp: fun_eq_iff)
    
lemma vector_derivative_of_real_left:
  assumes "f differentiable at x"
  shows   "vector_derivative (λx. of_real (f x)) (at x) = of_real (deriv f x)"
proof -
  have "vector_derivative (of_real  f) (at x) = (of_real (deriv f x))"
    by (subst vector_derivative_chain_at)
       (simp_all add: scaleR_conv_of_real field_derivative_eq_vector_derivative assms)
  thus ?thesis by (simp add: o_def)
qed
  
lemma vector_derivative_of_real_right:
  assumes "f field_differentiable at (of_real x)"
  shows   "vector_derivative (λx. f (of_real x)) (at x) = deriv f (of_real x)"
proof -
  have "vector_derivative (f  of_real) (at x) = deriv f (of_real x)"
    using assms by (subst vector_derivative_chain_at_general) simp_all
  thus ?thesis by (simp add: o_def)
qed
  
lemma Ln_holomorphic [holomorphic_intros]:
  assumes "A  0 = {}"
  shows   "Ln holomorphic_on (A :: complex set)"
proof (intro holomorphic_onI)
  fix z assume "z  A"
  with assms have "(Ln has_field_derivative inverse z) (at z within A)"
    by (auto intro!: derivative_eq_intros)
  thus "Ln field_differentiable at z within A" by (auto simp: field_differentiable_def)
qed

lemma higher_deriv_Polygamma:
  assumes "z  0"
  shows   "(deriv ^^ n) (Polygamma m) z = 
             Polygamma (m + n) (z :: 'a :: {real_normed_field,euclidean_space})"
proof -
  have "eventually (λu. (deriv ^^ n) (Polygamma m) u = Polygamma (m + n) u) (nhds z)"
  proof (induction n)
    case (Suc n)
    from Suc.IH have "eventually (λz. eventually (λu. (deriv ^^ n) (Polygamma m) u = Polygamma (m + n) u) (nhds z)) (nhds z)"
      by (simp add: eventually_eventually)
    hence "eventually (λz. deriv ((deriv ^^ n) (Polygamma m)) z = 
             deriv (Polygamma (m + n)) z) (nhds z)"
      by eventually_elim (intro deriv_cong_ev refl)
    moreover have "eventually (λz. z  UNIV - 0) (nhds z)" using assms
      by (intro eventually_nhds_in_open open_Diff open_UNIV) auto
    ultimately show ?case by eventually_elim (simp_all add: deriv_Polygamma)
  qed simp_all
  thus ?thesis by (rule eventually_nhds_x_imp_x)
qed
  
lemma higher_deriv_cmult:
  assumes "f holomorphic_on A" "x  A" "open A"
  shows   "(deriv ^^ j) (λx. c * f x) x = c * (deriv ^^ j) f x"
  using assms
proof (induction j arbitrary: f x)
  case (Suc j f x)
  have "deriv ((deriv ^^ j) (λx. c * f x)) x = deriv (λx. c * (deriv ^^ j) f x) x"
    using eventually_nhds_in_open[of A x] assms(2,3) Suc.prems
    by (intro deriv_cong_ev refl) (auto elim!: eventually_mono simp: Suc.IH)
  also have " = c * deriv ((deriv ^^ j) f) x" using Suc.prems assms(2,3)
    by (intro deriv_cmult holomorphic_on_imp_differentiable_at holomorphic_higher_deriv) auto
  finally show ?case by simp
qed simp_all

lemma higher_deriv_ln_Gamma_complex:
  assumes "(x::complex)  0"
  shows   "(deriv ^^ j) ln_Gamma x = (if j = 0 then ln_Gamma x else Polygamma (j - 1) x)"
proof (cases j)
  case (Suc j')
  have "(deriv ^^ j') (deriv ln_Gamma) x = (deriv ^^ j') Digamma x"
    using eventually_nhds_in_open[of "UNIV - 0" x] assms
    by (intro higher_deriv_cong_ev refl)
       (auto elim!: eventually_mono simp: open_Diff deriv_ln_Gamma_complex)
  also have " = Polygamma j' x" using assms
    by (subst higher_deriv_Polygamma)
       (auto elim!: nonpos_Ints_cases simp: complex_nonpos_Reals_iff)
  finally show ?thesis using Suc by (simp del: funpow.simps add: funpow_Suc_right)
qed simp_all

lemma higher_deriv_ln_Gamma_real:
  assumes "(x::real) > 0"
  shows   "(deriv ^^ j) ln_Gamma x = (if j = 0 then ln_Gamma x else Polygamma (j - 1) x)"
proof (cases j)
  case (Suc j')
  have "(deriv ^^ j') (deriv ln_Gamma) x = (deriv ^^ j') Digamma x"
    using eventually_nhds_in_open[of "{0<..}" x] assms
    by (intro higher_deriv_cong_ev refl)
       (auto elim!: eventually_mono simp: open_Diff deriv_ln_Gamma_real)
  also have " = Polygamma j' x" using assms
    by (subst higher_deriv_Polygamma)
       (auto elim!: nonpos_Ints_cases simp: complex_nonpos_Reals_iff)
  finally show ?thesis using Suc by (simp del: funpow.simps add: funpow_Suc_right)
qed simp_all
  
lemma higher_deriv_ln_Gamma_complex_of_real:
  assumes "(x :: real) > 0"
  shows   "(deriv ^^ j) ln_Gamma (complex_of_real x) = of_real ((deriv ^^ j) ln_Gamma x)"
    using assms
    by (auto simp: higher_deriv_ln_Gamma_real higher_deriv_ln_Gamma_complex
                   ln_Gamma_complex_of_real Polygamma_of_real)
(* END TODO *)

(* TODO: could be automated with Laurent series expansions in the future *)
lemma stirling_limit_aux1: 
  "((λy. Ln (1 + z * of_real y) / of_real y)  z) (at_right 0)" for z :: complex
proof (cases "z = 0")
  case True
  then show ?thesis by simp
next
  case False
  have "((λy. ln (1 + z * of_real y)) has_vector_derivative 1 * z) (at 0)"
    by (rule has_vector_derivative_real_field) (auto intro!: derivative_eq_intros)
  then have "(λy. (Ln (1 + z * of_real y) - of_real y * z) / of_real ¦y¦) 0 0"
    by (auto simp add: has_vector_derivative_def has_derivative_def netlimit_at 
          scaleR_conv_of_real field_simps)
  then have "((λy. (Ln (1 + z * of_real y) - of_real y * z) / of_real ¦y¦)  0) (at_right 0)"
    by (rule filterlim_mono[OF _ _ at_le]) simp_all
  also have "?this  ((λy. Ln (1 + z * of_real y) / (of_real y) - z)  0) (at_right 0)"
    using eventually_at_right_less[of "0::real"]
    by (intro filterlim_cong refl) (auto elim!: eventually_mono simp: field_simps)
  finally show ?thesis by (simp only: LIM_zero_iff)
qed
  
lemma stirling_limit_aux2: 
  "((λy. y * Ln (1 + z / of_real y))  z) at_top" for z :: complex
  using stirling_limit_aux1[of z] by (subst filterlim_at_top_to_right) (simp add: field_simps)

lemma Union_atLeastAtMost: 
  assumes "N > 0" 
  shows   "(n{0..<N}. {real n..real (n + 1)}) = {0..real N}"
proof (intro equalityI subsetI)
  fix x assume x: "x  {0..real N}"
  thus "x  (n{0..<N}. {real n..real (n + 1)})"
  proof (cases "x = real N")
    case True
    with assms show ?thesis by (auto intro!: bexI[of _ "N - 1"])
  next
    case False
    with x have x: "x  0" "x < real N" by simp_all
    hence "x  real (nat x)" "x  real (nat x + 1)" by linarith+
    moreover from x have "nat x < N" by linarith
    ultimately have "n{0..<N}. x  {real n..real (n + 1)}"
      by (intro bexI[of _ "nat x"]) simp_all
    thus ?thesis by blast
  qed
qed auto


subsection ‹Cones in the complex plane›

definition complex_cone :: "real  real  complex set" where
  "complex_cone a b = {z. y{a..b}. z = rcis (norm z) y}"

abbreviation complex_cone' :: "real  complex set" where
  "complex_cone' a  complex_cone (-a) a"

lemma zero_in_complex_cone [simp, intro]: "a  b  0  complex_cone a b"
  by (auto simp: complex_cone_def)

lemma complex_coneE:
  assumes "z  complex_cone a b"
  obtains r α where "r  0" "α  {a..b}" "z = rcis r α"
proof -
  from assms obtain y where "y  {a..b}" "z = rcis (norm z) y"
    unfolding complex_cone_def by auto
  thus ?thesis using that[of "norm z" y] by auto
qed

lemma arg_cis [simp]:
  assumes "x  {-pi<..pi}"
  shows   "arg (cis x) = x"
  using assms by (intro arg_unique) auto

lemma arg_mult_of_real_left [simp]:
  assumes "r > 0"
  shows   "arg (of_real r * z) = arg z"
proof (cases "z = 0")
  case False
  thus ?thesis
    using arg_bounded[of z] assms
    by (intro arg_unique) (auto simp: sgn_mult sgn_of_real cis_arg)
qed auto

lemma arg_mult_of_real_right [simp]:
  assumes "r > 0"
  shows   "arg (z * of_real r) = arg z"
  by (subst mult.commute, subst arg_mult_of_real_left) (simp_all add: assms)

lemma arg_rcis [simp]:
  assumes "x  {-pi<..pi}" "r > 0"
  shows   "arg (rcis r x) = x"
  using assms by (simp add: rcis_def)

lemma rcis_in_complex_cone [intro]: 
  assumes "α  {a..b}" "r  0"
  shows   "rcis r α  complex_cone a b"
  using assms by (auto simp: complex_cone_def)  

lemma arg_imp_in_complex_cone:
  assumes "arg z  {a..b}"
  shows   "z  complex_cone a b"
proof -
  have "z = rcis (norm z) (arg z)"
    by (simp add: rcis_cmod_arg)
  also have "  complex_cone a b"
    using assms by auto
  finally show ?thesis .
qed

lemma complex_cone_altdef:
  assumes "-pi < a" "a  b" "b  pi"
  shows   "complex_cone a b = insert 0 {z. arg z  {a..b}}"
proof (intro equalityI subsetI)
  fix z assume "z  complex_cone a b"
  then obtain r α where *: "r  0" "α  {a..b}" "z = rcis r α"
    by (auto elim: complex_coneE)
  have "arg z  {a..b}" if [simp]: "z  0"
  proof -
    have "r > 0" using that * by (subst (asm) *) auto
    hence "α  {a..b}"
      using *(1,2) assms by (auto simp: *(1))
    moreover from assms *(2) have "α  {-pi<..pi}"
      by auto
    ultimately show ?thesis using *(3) r > 0
      by (subst *) auto
  qed
  thus "z  insert 0 {z. arg z  {a..b}}"
    by auto
qed (use assms in auto intro: arg_imp_in_complex_cone›)

lemma nonneg_of_real_in_complex_cone [simp, intro]:
  assumes "x  0" "a  0" "0  b"
  shows   "of_real x  complex_cone a b"
proof -
  from assms have "rcis x 0  complex_cone a b"
    by (intro rcis_in_complex_cone) auto
  thus ?thesis by simp
qed

lemma one_in_complex_cone [simp, intro]: "a  0  0  b  1  complex_cone a b"
  using nonneg_of_real_in_complex_cone[of 1] by (simp del: nonneg_of_real_in_complex_cone)

lemma of_nat_in_complex_cone [simp, intro]: "a  0  0  b  of_nat n  complex_cone a b"
  using nonneg_of_real_in_complex_cone[of "real n"] by (simp del: nonneg_of_real_in_complex_cone)


subsection ‹Another integral representation of the Beta function›

lemma complex_cone_inter_nonpos_Reals:
  assumes "-pi < a" "a  b" "b < pi"
  shows   "complex_cone a b  0 = {0}"
proof (safe elim!: nonpos_Reals_cases)
  fix x :: real
  assume "complex_of_real x  complex_cone a b" "x  0"
  hence "¬(x < 0)"
    using assms by (intro notI) (auto simp: complex_cone_altdef)
  with x  0 show "complex_of_real x = 0" by auto  
qed (use assms in auto)

theorem 
  assumes a: "a > 0" and b: "b > (0 :: real)"
  shows has_integral_Beta_real': 
          "((λu. u powr (b - 1) / (1 + u) powr (a + b)) has_integral Beta a b) {0<..}"
    and Beta_conv_nn_integral:
          "Beta a b = (+u. ennreal (indicator {0<..} u * u powr (b - 1) / (1 + u) powr (a + b)) lborel)"
proof -
  define I where 
    "I = (+u. ennreal (indicator {0<..} u * u powr (b - 1) / (1 + u) powr (a + b)) lborel)"
  have "Gamma (a + b) > 0" "Beta a b > 0"
    using assms by (simp_all add: add_pos_pos Beta_def)
  from a b have "ennreal (Gamma a * Gamma b) =
    (+ t. ennreal (indicator {0..} t * t powr (a - 1) / exp t) lborel) *
    (+ t. ennreal (indicator {0..} t * t powr (b - 1) / exp t) lborel)"
    by (subst ennreal_mult') (simp_all add: Gamma_conv_nn_integral_real)
  also have " = (+t. +u. ennreal (indicator {0..} t * t powr (a - 1) / exp t) *
                            ennreal (indicator {0..} u * u powr (b - 1) / exp u) lborel lborel)"
    by (simp add: nn_integral_cmult nn_integral_multc)
  also have " = (+t. indicator {0<..} t * (+u. indicator {0..} u * t powr (a - 1) * u powr (b - 1)
                            / exp (t + u) lborel) lborel)"
    by (intro nn_integral_cong_AE AE_I[of _ _ "{0}"])
       (auto simp: indicator_def divide_ennreal ennreal_mult' [symmetric] exp_add mult_ac)
  also have " = (+t. indicator {0<..} t * (+u. indicator {0..} u * t powr (a - 1) * u powr (b - 1)
                            / exp (t + u) 
                    (density (distr lborel borel ((*) t)) (λx. ennreal ¦t¦))) lborel)"
    by (intro nn_integral_cong mult_indicator_cong, subst lborel_distr_mult' [symmetric]) auto
  also have " = (+(t::real). indicator {0<..} t * (+u. 
                     indicator {0..} (u * t) * t powr a *
                     (u * t) powr (b - 1) / exp (t + t * u) lborel) lborel)"
    by (intro nn_integral_cong mult_indicator_cong)
       (auto simp: nn_integral_density nn_integral_distr algebra_simps powr_diff
             simp flip: ennreal_mult)
  also have " = (+(t::real). +u. indicator ({0<..}×{0..}) (t, u) *
                     t powr a * (u * t) powr (b - 1) / exp (t * (u + 1)) lborel lborel)"
    by (subst nn_integral_cmult [symmetric], simp, intro nn_integral_cong)
       (auto simp: indicator_def zero_le_mult_iff algebra_simps)
  also have " = (+(t::real). +u. indicator ({0<..}×{0..}) (t, u) *
                     t powr (a + b - 1) * u powr (b - 1) / exp (t * (u + 1)) lborel lborel)"
    by (intro nn_integral_cong) (auto simp: powr_add powr_diff indicator_def powr_mult field_simps)
  also have " = (+(u::real). +t. indicator ({0<..}×{0..}) (t, u) *
                     t powr (a + b - 1) * u powr (b - 1) / exp (t * (u + 1)) lborel lborel)"
    by (rule lborel_pair.Fubini') auto
  also have " = (+(u::real). indicator {0..} u * (+t. indicator {0<..} t *
                     t powr (a + b - 1) * u powr (b - 1) / exp (t * (u + 1)) lborel) lborel)"
    by (intro nn_integral_cong mult_indicator_cong) (auto simp: indicator_def)
  also have " = (+(u::real). indicator {0<..} u * (+t. indicator {0<..} t *
                     t powr (a + b - 1) * u powr (b - 1) / exp (t * (u + 1)) lborel) lborel)"
    by (intro nn_integral_cong_AE AE_I[of _ _ "{0}"]) (auto simp: indicator_def)
  also have " = (+(u::real). indicator {0<..} u * (+t. indicator {0<..} t *
                     t powr (a + b - 1) * u powr (b - 1) / exp (t * (u + 1)) 
                    (density (distr lborel borel ((*) (1/(1+u)))) (λx. ennreal ¦1/(1+u)¦))) lborel)"
    by (intro nn_integral_cong mult_indicator_cong, subst lborel_distr_mult' [symmetric]) auto
  also have " = (+(u::real). indicator {0<..} u * 
                    (+t. ennreal (1 / (u + 1)) * ennreal (indicator {0<..} (t / (u + 1)) *
                     (t / (1+u)) powr (a + b - 1) * u powr (b - 1) / exp t)
                    lborel) lborel)"
    by (intro nn_integral_cong mult_indicator_cong)
       (auto simp: nn_integral_distr nn_integral_density add_ac)
  also have " = (+u. +t. indicator ({0<..}×{0<..}) (u, t) * 
                    1/(u+1) * (t / (u+1)) powr (a + b - 1) * u powr (b - 1) / exp t
                    lborel lborel)"
    by (subst nn_integral_cmult [symmetric], simp, intro nn_integral_cong)
       (auto simp: indicator_def field_simps divide_ennreal simp flip: ennreal_mult ennreal_mult')
  also have " = (+u. +t. ennreal (indicator {0<..} u * u powr (b - 1) / (1 + u) powr (a + b)) *
                            ennreal (indicator {0<..} t * t powr (a + b - 1) / exp t)
                    lborel lborel)"
    by (intro nn_integral_cong)
       (auto simp: indicator_def powr_add powr_diff powr_divide powr_minus divide_simps add_ac
             simp flip: ennreal_mult)
  also have " = I * (+t. indicator {0<..} t * t powr (a + b - 1) / exp t lborel)"
    by (simp add: nn_integral_cmult nn_integral_multc I_def)
  also have "(+t. indicator {0<..} t * t powr (a + b - 1) / exp t lborel) =
               ennreal (Gamma (a + b))"
    using assms
    by (subst Gamma_conv_nn_integral_real)
       (auto intro!: nn_integral_cong_AE[OF AE_I[of _ _ "{0}"]] 
             simp: indicator_def split: if_splits)
  finally have "ennreal (Gamma a * Gamma b) = I * ennreal (Gamma (a + b))" .
  hence "ennreal (Gamma a * Gamma b) / ennreal (Gamma (a + b)) =
           I * ennreal (Gamma (a + b)) / ennreal (Gamma (a + b))" by simp
  also have " = I"
    using ‹Gamma (a + b) > 0 by (intro ennreal_mult_divide_eq) (auto simp: )
  also have "ennreal (Gamma a * Gamma b) / ennreal (Gamma (a + b)) =
               ennreal (Gamma a * Gamma b / Gamma (a + b))"
    using assms by (intro divide_ennreal) auto
  also have " = ennreal (Beta a b)"
    by (simp add: Beta_def)
  finally show *: "ennreal (Beta a b) = I" .

  define f where "f = (λu. u powr (b - 1) / (1 + u) powr (a + b))"
  have "((λu. indicator {0<..} u * f u) has_integral Beta a b) UNIV"
    using * ‹Beta a b > 0
    by (subst has_integral_iff_nn_integral_lebesgue)
       (auto simp: f_def measurable_completion nn_integral_completion I_def mult_ac)
  also have "(λu. indicator {0<..} u * f u) = (λu. if u  {0<..} then f u else 0)"
    by (auto simp: fun_eq_iff)
  also have "( has_integral Beta a b) UNIV  (f has_integral Beta a b) {0<..}"
    by (rule has_integral_restrict_UNIV)
  finally show  by (simp add: f_def)
qed

lemma has_integral_Beta2:
  fixes a :: real
  assumes "a < -1/2"
  shows   "((λx. (1 + x ^ 2) powr a) has_integral Beta (- a - 1 / 2) (1 / 2) / 2) {0<..}"
proof -
  define f where "f = (λu. u powr (-1/2) / (1 + u) powr (-a))"
  define C where "C = Beta (-a-1/2) (1/2)"
  have I: "(f has_integral C) {0<..}"
    using has_integral_Beta_real'[of "-a-1/2" "1/2"] assms
    by (simp_all add: diff_divide_distrib f_def C_def)

  define g where "g = (λx. x ^ 2 :: real)"
  have bij: "bij_betw g {0<..} {0<..}"
    by (intro bij_betwI[of _ _ _ sqrt]) (auto simp: g_def)

  have "(f absolutely_integrable_on g ` {0<..}  integral (g ` {0<..}) f = C)"
    using I bij by (simp add: bij_betw_def has_integral_iff absolutely_integrable_on_def f_def)
  also have "?this  ((λx. ¦2 * x¦ *R f (g x)) absolutely_integrable_on {0<..} 
                         integral {0<..} (λx. ¦2 * x¦ *R f (g x)) = C)"
    using bij by (intro has_absolute_integral_change_of_variables_1' [symmetric])
                 (auto intro!: derivative_eq_intros simp: g_def bij_betw_def)
  finally have "((λx. ¦2 * x¦ * f (g x)) has_integral C) {0<..}"
    by (simp add: absolutely_integrable_on_def f_def has_integral_iff)
  also have "?this  ((λx::real. 2 * (1 + x2) powr a) has_integral C) {0<..}"
    by (intro has_integral_cong) (auto simp: f_def g_def powr_def exp_minus ln_realpow field_simps)
  finally have "((λx::real. 1/2 * (2 * (1 + x2) powr a)) has_integral 1/2 * C) {0<..}"
    by (intro has_integral_mult_right)
  thus ?thesis by (simp add: C_def)
qed

lemma has_integral_Beta3:
  fixes a b :: real
  assumes "a < -1/2" and "b > 0"
  shows   "((λx. (b + x ^ 2) powr a) has_integral
             Beta (-a - 1/2) (1/2) / 2 * b powr (a + 1/2)) {0<..}"
proof -
  define C where "C = Beta (- a - 1 / 2) (1 / 2) / 2"
  have int: "nn_integral lborel (λx. indicator {0<..} x * (1 + x ^ 2) powr a) = C"
    using nn_integral_has_integral_lebesgue[OF _ has_integral_Beta2[OF assms(1)]]
    by (auto simp: C_def)
  have "nn_integral lborel (λx. indicator {0<..} x * (b + x ^ 2) powr a) =
        (+x. ennreal (indicat_real {0<..} (x * sqrt b) * (b + (x * sqrt b)2) powr a * sqrt b) lborel)"
    using assms
    by (subst lborel_distr_mult'[of "sqrt b"])
       (auto simp: nn_integral_density nn_integral_distr mult_ac simp flip: ennreal_mult)
  also have " = (+x. ennreal (indicat_real {0<..} x * (b * (1 + x ^ 2)) powr a * sqrt b) lborel)"
    using assms
    by (intro nn_integral_cong) (auto simp: indicator_def field_simps zero_less_mult_iff)
  also have " = (+x. ennreal (indicat_real {0<..} x * b powr (a + 1/2) * (1 + x ^ 2) powr a) lborel)"
    using assms
    by (intro nn_integral_cong) (auto simp: indicator_def powr_add powr_half_sqrt powr_mult)    
  also have " = b powr (a + 1/2) * (+x. ennreal (indicat_real {0<..} x * (1 + x ^ 2) powr a) lborel)"
    using assms by (subst nn_integral_cmult [symmetric]) (simp_all add: mult_ac flip: ennreal_mult)
  also have "(+x. ennreal (indicat_real {0<..} x * (1 + x ^ 2) powr a) lborel) = C"
    using int by simp
  also have "ennreal (b powr (a + 1/2)) * ennreal C = ennreal (C * b powr (a + 1/2))"
    using assms by (subst ennreal_mult) (auto simp: C_def mult_ac Beta_def)
  finally have *: "(+ x. ennreal (indicat_real {0<..} x * (b + x2) powr a) lborel) = " .
  hence "((λx. indicator {0<..} x * (b + x^2) powr a) has_integral C * b powr (a + 1/2)) UNIV"
    using assms
    by (subst has_integral_iff_nn_integral_lebesgue)
       (auto simp: C_def measurable_completion nn_integral_completion Beta_def)
  also have "(λx. indicator {0<..} x * (b + x^2) powr a) =
             (λx. if x  {0<..} then (b + x^2) powr a else 0)"
    by (auto simp: fun_eq_iff)
  finally show ?thesis
    by (subst (asm) has_integral_restrict_UNIV) (auto simp: C_def)
qed

  
subsection ‹Asymptotics of the real $\log\Gamma$ function and its derivatives›

text ‹
  This is the error term that occurs in the expansion of @{term ln_Gamma}. It can be shown to 
  be of order $O(s^{-n})$.
›
definition stirling_integral :: "nat  'a :: {real_normed_div_algebra, banach}  'a" where
  "stirling_integral n s = 
     lim (λN. integral {0..N} (λx. of_real (pbernpoly n x) / (of_real x + s) ^ n))"

context
  fixes s :: complex assumes s: "s  0"
  fixes approx :: "nat  complex"
  defines "approx  (λN. 
    (n = 1..<N. s / of_nat n - ln (1 + s / of_nat n)) - (euler_mascheroni * s + ln s) - ― ‹⟶ ln_Gamma s›
    (ln_Gamma (of_nat N) - ln (2 * pi / of_nat N) / 2 - of_nat N * ln (of_nat N) + of_nat N) - ― ‹⟶ 0›
    s * (harm (N - 1) - ln (of_nat (N - 1)) - euler_mascheroni) + ― ‹⟶ 0›
    s * (ln (of_nat N + s) - ln (of_nat (N - 1))) - ― ‹⟶ 0›
    (1/2) * (ln (of_nat N + s) - ln (of_nat N)) +       ― ‹⟶ 0›
    of_nat N * (ln (of_nat N + s) - ln (of_nat N)) -  ― ‹⟶ s›
    (s - 1/2) * ln s - ln (2 * pi) / 2)"
begin       
  
qualified lemma
  assumes N: "N > 0"
  shows   integrable_pbernpoly_1:
            "(λx. of_real (-pbernpoly 1 x) / (of_real x + s)) integrable_on {0..real N}"
  and     integral_pbernpoly_1_aux:
            "integral {0..real N} (λx. -of_real (pbernpoly 1 x) / (of_real x + s)) = approx N"
  and     has_integral_pbernpoly_1:
            "((λx. pbernpoly 1 x /(x + s)) has_integral 
              (m<N. (of_nat m + 1 / 2 + s) * (ln (of_nat m + s) - 
                        ln (of_nat m + 1 + s)) + 1)) {0..real N}"
proof -
  let ?A = "(λn. {of_nat n..of_nat (n+1)}) ` {0..<N}"
  have has_integral: 
    "((λx. -pbernpoly 1 x / (x + s)) has_integral 
             (of_nat n + 1/2 + s) * (ln (of_nat (n + 1) + s) - ln (of_nat n + s)) - 1) 
           {of_nat n..of_nat (n + 1)}" for n
  proof (rule has_integral_spike)      
    have "((λx. (of_nat n + 1/2 + s) * (1 / (x + s)) - 1) has_integral 
              (of_nat n + 1/2 + s) * (ln (of_real (real (n + 1)) + s) - ln (of_real (real n) + s)) - 1) 
            {of_nat n..of_nat (n + 1)}" 
      using s has_integral_const_real[of 1 "of_nat n" "of_nat (n + 1)"]
      by (intro has_integral_diff has_integral_mult_right fundamental_theorem_of_calculus)
         (auto intro!: derivative_eq_intros has_vector_derivative_real_field
               simp: has_field_derivative_iff_has_vector_derivative [symmetric] field_simps
                     complex_nonpos_Reals_iff)
    thus "((λx. (of_nat n + 1/2 + s) * (1 / (x + s)) - 1) has_integral 
              (of_nat n + 1/2 + s) * (ln (of_nat (n + 1) + s) - ln (of_nat n + s)) - 1) 
            {of_nat n..of_nat (n + 1)}" by simp
             
    show "-pbernpoly 1 x / (x + s) = (of_nat n + 1/2 + s) * (1 / (x + s)) - 1"
         if "x  {of_nat n..of_nat (n + 1)} - {of_nat (n + 1)}" for x
    proof -
      have x: "x  real n" "x < real (n + 1)" using that by simp_all
      hence "floor x = int n" by linarith
      moreover from s x have "complex_of_real x  -s" 
        by (auto simp add: complex_eq_iff complex_nonpos_Reals_iff simp del: of_nat_Suc)
      ultimately show "-pbernpoly 1 x / (x + s) = (of_nat n + 1/2 + s) * (1 / (x + s)) - 1"
        by (auto simp: pbernpoly_def bernpoly_def frac_def divide_simps add_eq_0_iff2)
    qed
  qed simp_all
  hence *: "I. I?A  ((λx. -pbernpoly 1 x / (x + s)) has_integral 
              (Inf I + 1/2 + s) * (ln (Inf I + 1 + s) - ln (Inf I + s)) - 1) I"
    by (auto simp: add_ac)
  have "((λx. - pbernpoly 1 x / (x + s)) has_integral
          (I?A. (Inf I + 1 / 2 + s) * (ln (Inf I + 1 + s) - ln (Inf I + s)) - 1))
          (n{0..<N}. {real n..real (n + 1)})" (is "(_ has_integral ?i) _")
    apply (intro has_integral_Union * finite_imageI)
      apply (force intro!: negligible_atLeastAtMostI pairwiseI)+
    done
  hence has_integral: "((λx. - pbernpoly 1 x / (x + s)) has_integral ?i) {0..real N}"
    by (subst has_integral_spike_set_eq)
       (use Union_atLeastAtMost assms in auto simp: intro!: empty_imp_negligible›)
  hence "(λx. - pbernpoly 1 x / (x + s)) integrable_on {0..real N}"
    and integral:   "integral {0..real N} (λx. - pbernpoly 1 x / (x + s)) = ?i"
    by (simp_all add: has_integral_iff)
  show "(λx. - pbernpoly 1 x / (x + s)) integrable_on {0..real N}" by fact

  note has_integral_neg[OF has_integral]
  also have "-?i = (x<N. (of_nat x + 1 / 2 + s) * (ln (of_nat x + s) - ln (of_nat x + 1 + s)) + 1)" 
    by (subst sum.reindex) 
       (simp_all add: inj_on_def atLeast0LessThan algebra_simps sum_negf [symmetric])
  finally show has_integral: 
    "((λx. of_real (pbernpoly 1 x) / (of_real x + s)) has_integral ) {0..real N}" by simp
      
  note integral
  also have "?i = (n<N. (of_nat n + 1 / 2 + s) * 
                    (ln (of_nat n + 1 + s) - ln (of_nat n + s))) - N" (is "_ = ?S - _")
    by (subst sum.reindex) (simp_all add: inj_on_def sum_subtractf atLeast0LessThan)
  also have "?S = (n<N. of_nat n * (ln (of_nat n + 1 + s) - ln (of_nat n + s))) +
                    (s + 1 / 2) * (n<N. ln (of_nat (Suc n) + s) - ln (of_nat n + s))" 
    (is "_ = ?S1 + _ * ?S2") by (simp add: algebra_simps sum.distrib sum_subtractf sum_distrib_left)
  also have "?S2 = ln (of_nat N + s) - ln s" by (subst sum_lessThan_telescope) simp
  also have "?S1 = (n=1..<N. of_nat n * (ln (of_nat n + 1 + s) - ln (of_nat n + s)))"
    by (intro sum.mono_neutral_right) auto
  also have " = (n=1..<N. of_nat n * ln (of_nat n + 1 + s)) - (n=1..<N. of_nat n * ln (of_nat n + s))"
    by (simp add: algebra_simps sum_subtractf)
  also have "(n=1..<N. of_nat n * ln (of_nat n + 1 + s)) = 
               (n=1..<N. (of_nat n - 1) * ln (of_nat n + s)) + (N - 1) * ln (of_nat N + s)"
    by (induction N) (simp_all add: add_ac of_nat_diff)
  also have " - (n = 1..<N. of_nat n * ln (of_nat n + s)) =
               -(n=1..<N. ln (of_nat n + s)) + (N - 1) * ln (of_nat N + s)"
    by (induction N) (simp_all add: algebra_simps)
  also from s have neq: "s + of_nat x  0" for x
    by (auto simp:  complex_nonpos_Reals_iff complex_eq_iff)
  hence "(n=1..<N. ln (of_nat n + s)) = (n=1..<N. ln (of_nat n) + ln (1 + s/n))"
    by (intro sum.cong refl, subst Ln_times_of_nat [symmetric]) (auto simp: divide_simps add_ac)
  also have " = ln (fact (N - 1)) + (n=1..<N. ln (1 + s/n))"
    by (induction N) (simp_all add: Ln_times_of_nat fact_reduce add_ac)
  also have "(n=1..<N. ln (1 + s/n)) = -(n=1..<N. s / n - ln (1 + s/n)) + s * (n=1..<N. 1 / of_nat n)"
    by (simp add: sum_distrib_left sum_subtractf) 
  also from N have "ln (fact (N - 1)) = ln_Gamma (of_nat N :: complex)" 
    by (simp add: ln_Gamma_complex_conv_fact)
  also have "{1..<N} = {1..N - 1}" by auto
  hence "(n = 1..<N. 1 / of_nat n) = (harm (N - 1) :: complex)" 
    by (simp add: harm_def divide_simps)
  also have "- (ln_Gamma (of_nat N) + (- (n = 1..<N. s / of_nat n - ln (1 + s / of_nat n)) +
                 s * harm (N - 1))) + of_nat (N - 1) * ln (of_nat N + s) +
                (s + 1 / 2) * (ln (of_nat N + s) - ln s) - of_nat N = approx N"
    using N by (simp add: field_simps of_nat_diff ln_div approx_def Ln_of_nat 
                          ln_Gamma_complex_of_real [symmetric])
  finally show "integral {0..of_nat N} (λx. -of_real (pbernpoly 1 x) / (of_real x + s)) = " 
    by simp
qed
  
lemma integrable_ln_Gamma_aux:
  shows   "(λx. of_real (pbernpoly n x) / (of_real x + s) ^ n) integrable_on {0..real N}"
proof (cases "n = 1")
  case True
  with s show ?thesis using integrable_neg[OF integrable_pbernpoly_1[of N]] 
    by (cases "N = 0") (simp_all add: integrable_negligible)
next
  case False
  from s have "of_real x + s  0" if "x  0" for x using that 
    by (auto simp: complex_eq_iff add_eq_0_iff2 complex_nonpos_Reals_iff)
  with False s show ?thesis
    by (auto intro!: integrable_continuous_real continuous_intros)
qed
  
text ‹
  This following proof is based on ``Rudiments of the theory of the gamma function'' 
  by Bruce Berndt~\cite{berndt}.
›
qualified lemma integral_pbernpoly_1: 
  "(λN. integral {0..real N} (λx. pbernpoly 1 x / (x + s)))
      -ln_Gamma s - s + (s - 1 / 2) * ln s + ln (2 * pi) / 2"
proof -  
  have neq: "s + of_real x  0" if "x  0" for x :: real
    using that s by (auto simp: complex_eq_iff complex_nonpos_Reals_iff)
  have "(approx  ln_Gamma s - 0 - 0 + 0 - 0 + s - (s - 1/2) * ln s - ln (2 * pi) / 2) at_top"
    unfolding approx_def
  proof (intro tendsto_add tendsto_diff)
    from s have s': "s  0" by (auto simp: complex_nonpos_Reals_iff elim!: nonpos_Ints_cases)
    have "(λn. i=1..<n. s / of_nat i - ln (1 + s / of_nat i))  
             ln_Gamma s + euler_mascheroni * s + ln s" (is "?f  _")
      using ln_Gamma_series'_aux[OF s'] unfolding sums_def 
      by (subst filterlim_sequentially_Suc [symmetric], subst (asm) sum.atLeast1_atMost_eq [symmetric]) 
         (simp add: atLeastLessThanSuc_atLeastAtMost)
    thus "((λn. ?f n - (euler_mascheroni * s + ln s))  ln_Gamma s) at_top"
      by (auto intro: tendsto_eq_intros)
  next
    show "(λx. complex_of_real (ln_Gamma (real x) - ln (2 * pi / real x) / 2 - 
                 real x * ln (real x) + real x))  0"
    proof (intro tendsto_of_real_0_I 
             filterlim_compose[OF tendsto_sandwich filterlim_real_sequentially])
      show "eventually (λx::real. ln_Gamma x - ln (2 * pi / x) / 2 - x * ln x + x  0) at_top"
        using eventually_ge_at_top[of "1::real"] 
        by eventually_elim (insert ln_Gamma_bounds(1), simp add: algebra_simps)
      show "eventually (λx::real. ln_Gamma x - ln (2 * pi / x) / 2 - x * ln x + x  
              1 / 12 * inverse x) at_top"
        using eventually_ge_at_top[of "1::real"] 
        by eventually_elim (insert ln_Gamma_bounds(2), simp add: field_simps)
      show "((λx::real. 1 / 12 * inverse x)  0) at_top"
        by (intro tendsto_mult_right_zero tendsto_inverse_0_at_top filterlim_ident)
    qed simp_all
  next
    have "(λx. s * of_real (harm (x - 1) - ln (real (x - 1)) - euler_mascheroni))  
            s * of_real (euler_mascheroni - euler_mascheroni)"
      by (subst filterlim_sequentially_Suc [symmetric], intro tendsto_intros) 
         (insert euler_mascheroni_LIMSEQ, simp_all)
    also have "?this  (λx. s * (harm (x - 1) - ln (of_nat (x - 1)) - euler_mascheroni))  0"
      by (intro filterlim_cong refl eventually_mono[OF eventually_gt_at_top[of "1::nat"]]) 
         (auto simp: Ln_of_nat of_real_harm)
    finally show "(λx. s * (harm (x - 1) - ln (of_nat (x - 1)) - euler_mascheroni))  0"  .
  next
    have "((λx. ln (1 + (s + 1) / of_real x))  ln (1 + 0)) at_top" (is ?P)
      by (intro tendsto_intros tendsto_divide_0[OF tendsto_const]) 
         (simp_all add: filterlim_ident filterlim_at_infinity_conv_norm_at_top filterlim_abs_real)
    also have "ln (of_real (x + 1) + s) - ln (complex_of_real x) = ln (1 + (s + 1) / of_real x)" 
      if "x > 1" for x using that s
      using Ln_divide_of_real[of x "of_real (x + 1) + s", symmetric] neq[of "x+1"]
      by (simp add: field_simps Ln_of_real)
    hence "?P  ((λx. ln (of_real (x + 1) + s) - ln (of_real x))  0) at_top"
      by (intro filterlim_cong refl) 
         (auto intro: eventually_mono[OF eventually_gt_at_top[of "1::real"]])
    finally have "((λn. ln (of_real (real n + 1) + s) - ln (of_real (real n)))  0) at_top"
      by (rule filterlim_compose[OF _ filterlim_real_sequentially])
    hence "((λn. ln (of_nat n + s) - ln (of_nat (n - 1)))  0) at_top"
      by (subst filterlim_sequentially_Suc [symmetric]) (simp add: add_ac)
    thus "(λx. s * (ln (of_nat x + s) - ln (of_nat (x - 1))))  0"
      by (rule tendsto_mult_right_zero)
  next
    have "((λx. ln (1 + s / of_real x))  ln (1 + 0)) at_top" (is ?P)
      by (intro tendsto_intros tendsto_divide_0[OF tendsto_const]) 
         (simp_all add: filterlim_ident  filterlim_at_infinity_conv_norm_at_top filterlim_abs_real)
    also have "ln (of_real x + s) - ln (of_real x) = ln (1 + s / of_real x)" if "x > 0" for x
      using Ln_divide_of_real[of x "of_real x + s"] neq[of x] that
      by (auto simp: field_simps Ln_of_real)
    hence "?P  ((λx. ln (of_real x + s) - ln (of_real x))  0) at_top"
      using s by (intro filterlim_cong refl) 
                 (auto intro: eventually_mono [OF eventually_gt_at_top[of "1::real"]])
    finally have "(λx. (1/2) * (ln (of_real (real x) + s) - ln (of_real (real x))))  0"        
      by (rule tendsto_mult_right_zero[OF filterlim_compose[OF _ filterlim_real_sequentially]])
    thus "(λx. (1/2) * (ln (of_nat x + s) - ln (of_nat x)))  0" by simp
  next
    have "((λx. x * (ln (1 + s / of_real x)))  s) at_top" (is ?P) 
      by (rule stirling_limit_aux2)
    also have "ln (1 + s / of_real x) = ln (of_real x + s) - ln (of_real x)" if "x > 1" for x 
      using that s Ln_divide_of_real [of x "of_real x + s", symmetric] neq[of x]
      by (auto simp: Ln_of_real field_simps)
    hence "?P  ((λx. of_real x * (ln (of_real x + s) - ln (of_real x)))  s) at_top"
      by (intro filterlim_cong refl) 
         (auto intro: eventually_mono[OF eventually_gt_at_top[of "1::real"]])
    finally have "(λn. of_real (real n) * (ln (of_real (real n) + s) - ln (of_real (real n))))  s"
      by (rule filterlim_compose[OF _ filterlim_real_sequentially])
    thus "(λn. of_nat n * (ln (of_nat n + s) - ln (of_nat n)))  s" by simp
  qed simp_all
  also have "?this  ((λN. integral {0..real N} (λx. -pbernpoly 1 x / (x + s))) 
                         ln_Gamma s + s - (s - 1/2) * ln s - ln (2 * pi) / 2) at_top"
    using integral_pbernpoly_1_aux
    by (intro filterlim_cong refl) 
       (auto intro: eventually_mono[OF eventually_gt_at_top[of "0::nat"]])
  also have "(λN. integral {0..real N} (λx. -pbernpoly 1 x / (x + s))) =
               (λN. -integral {0..real N} (λx. pbernpoly 1 x / (x + s)))"
    by (simp add: fun_eq_iff)
  finally show ?thesis by (simp add: tendsto_minus_cancel_left [symmetric] algebra_simps)
qed


qualified lemma pbernpoly_integral_conv_pbernpoly_integral_Suc:
  assumes "n  1"
  shows   "integral {0..real N} (λx. pbernpoly n x / (x + s) ^ n) =
             of_real (pbernpoly (Suc n) (real N)) / (of_nat (Suc n) * (s + of_nat N) ^ n) -
             of_real (bernoulli (Suc n)) / (of_nat (Suc n) * s ^ n) + of_nat n / of_nat (Suc n) * 
               integral {0..real N} (λx. of_real (pbernpoly (Suc n) x) / (of_real x + s) ^ Suc n)"
proof - 
  note [derivative_intros] = has_field_derivative_pbernpoly_Suc'
  define I where "I = -of_real (pbernpoly (Suc n) (of_nat N)) / (of_nat (Suc n) * (of_nat N + s) ^ n) +
            of_real (bernoulli (Suc n) / real (Suc n)) / s ^ n +
            integral {0..real N} (λx. of_real (pbernpoly n x) / (of_real x + s) ^ n)"
  have "((λx. (-of_nat n * inverse (of_real x + s) ^ Suc n) * 
          (of_real (pbernpoly (Suc n) x) / (of_nat (Suc n))))
          has_integral -I) {0..real N}"
  proof (rule integration_by_parts_interior_strong[OF bounded_bilinear_mult])
    fix x :: real assume x: "x  {0<..<real N} - real ` {0..N}"
    have "x  "
    proof
      assume "x  "
      then obtain n where "x = of_int n" by (auto elim!: Ints_cases)
      with x have x': "x = of_nat (nat n)" by simp
      from x show False by (auto simp: x')
    qed
    hence "((λx. of_real (pbernpoly (Suc n) x / of_nat (Suc n))) has_vector_derivative
        complex_of_real (pbernpoly n x)) (at x)"
      by (intro has_vector_derivative_of_real) (auto intro!: derivative_eq_intros)
    thus "((λx. of_real (pbernpoly (Suc n) x) / of_nat (Suc n)) has_vector_derivative
            complex_of_real (pbernpoly n x)) (at x)" by simp
    from x s have "complex_of_real x + s  0"
      by (auto simp: complex_eq_iff complex_nonpos_Reals_iff)
    thus "((λx. inverse (of_real x + s) ^ n) has_vector_derivative 
             - of_nat n * inverse (of_real x + s) ^ Suc n) (at x)" using x s assms
      by (auto intro!: derivative_eq_intros has_vector_derivative_real_field simp: divide_simps power_add [symmetric]
               simp del: power_Suc)
  next
    have "complex_of_real x + s  0" if "x  0" for x 
      using that s by (auto simp: complex_eq_iff complex_nonpos_Reals_iff)
    thus "continuous_on {0..real N} (λx. inverse (of_real x + s) ^ n)" 
         "continuous_on {0..real N} (λx. complex_of_real (pbernpoly (Suc n) x) / of_nat (Suc n))"
      using assms s by (auto intro!: continuous_intros simp del: of_nat_Suc)
  next
    have "((λx. inverse (of_real x + s) ^ n * of_real (pbernpoly n x)) has_integral
            pbernpoly (Suc n) (of_nat N) / (of_nat (Suc n) * (of_nat N + s) ^ n) -
            of_real (bernoulli (Suc n) / real (Suc n)) / s ^ n - -I) {0..real N}" 
      using integrable_ln_Gamma_aux[of n N] assms
      by (auto simp: I_def has_integral_integral divide_simps)
    thus "((λx. inverse (of_real x + s) ^ n * of_real (pbernpoly n x)) has_integral
              inverse (of_real (real N) + s) ^ n * (of_real (pbernpoly (Suc n) (real N)) / 
                  of_nat (Suc n)) -
              inverse (of_real 0 + s) ^ n * (of_real (pbernpoly (Suc n) 0) / of_nat (Suc n)) - - I)
            {0..real N}" by (simp_all add: field_simps)
  qed simp_all
  also have "(λx. - of_nat n * inverse (of_real x + s) ^ Suc n * (of_real (pbernpoly (Suc n) x) /
                         of_nat (Suc n))) =
             (λx. - (of_nat n / of_nat (Suc n) * of_real (pbernpoly (Suc n) x) / 
                         (of_real x + s) ^ Suc n))"
    by (simp add: divide_simps fun_eq_iff)
  finally have "((λx. - (of_nat n / of_nat (Suc n) * of_real (pbernpoly (Suc n) x) /
                            (of_real x + s) ^ Suc n)) has_integral - I) {0..real N}" .
  from has_integral_neg[OF this] show ?thesis
    by (auto simp add: I_def has_integral_iff algebra_simps integral_mult_right [symmetric] 
             simp del: power_Suc of_nat_Suc )
qed

lemma pbernpoly_over_power_tendsto_0: 
  assumes "n > 0"
  shows   "(λx. of_real (pbernpoly (Suc n) (real x)) / (of_nat (Suc n) * (s + of_nat x) ^ n))  0"
proof -
  from s have neq: "s + of_nat n  0" for n
    by (auto simp: complex_eq_iff complex_nonpos_Reals_iff)
  from bounded_pbernpoly[of "Suc n"] guess c . note c = this
  have "eventually (λx. real x + Re s > 0) at_top"
    by real_asymp
  hence "eventually (λx. norm (of_real (pbernpoly (Suc n) (real x)) / 
                                    (of_nat (Suc n) * (s + of_nat x) ^ n)) 
                          (c / real (Suc n)) / (real x + Re s) ^ n) at_top"
    using eventually_gt_at_top[of "0::nat"]
  proof eventually_elim
    case (elim x)
    have "norm (of_real (pbernpoly (Suc n) (real x)) / 
                                    (of_nat (Suc n) * (s + of_nat x) ^ n)) 
            (c / real (Suc n)) / norm (s + of_nat x) ^ n" (is "_  ?rhs") using c[of x]
      by (auto simp: norm_divide norm_mult norm_power neq field_simps simp del: of_nat_Suc)
    also have "(real x + Re s)  cmod (s + of_nat x)"
      using complex_Re_le_cmod[of "s + of_nat x"] s by (auto simp add: complex_nonpos_Reals_iff)
    hence "?rhs  (c / real (Suc n)) / (real x + Re s) ^ n" using s elim c[of 0] neq[of x]
      by (intro divide_left_mono power_mono mult_pos_pos divide_nonneg_pos zero_less_power) auto
    finally show ?case .
  qed 
  moreover have "(λx. (c / real (Suc n)) / (real x + Re s) ^ n)  0"
    using n > 0 by real_asymp
  ultimately show ?thesis by (rule Lim_null_comparison)
qed

lemma convergent_stirling_integral:
  assumes "n > 0"
  shows   "convergent (λN. integral {0..real N} 
             (λx. of_real (pbernpoly n x) / (of_real x + s) ^ n))" (is "convergent (?f n)")
proof -
  have "convergent (?f (Suc n))" for n
  proof (induction n)
    case 0
    thus ?case using integral_pbernpoly_1 by (auto intro!: convergentI)
  next
    case (Suc n)
    have "convergent (λN. ?f (Suc n) N -
            of_real (pbernpoly (Suc (Suc n)) (real N)) / 
                (of_nat (Suc (Suc n)) * (s + of_nat N) ^ Suc n) +
            of_real (bernoulli (Suc (Suc n)) / (real (Suc (Suc n)))) / s ^ Suc n)" 
      (is "convergent ?g")
      by (intro convergent_add convergent_diff Suc 
            convergent_const convergentI[OF pbernpoly_over_power_tendsto_0]) simp_all
    also have "?g = (λN. of_nat (Suc n) / of_nat (Suc (Suc n)) * ?f (Suc (Suc n)) N)" using s
      by (subst pbernpoly_integral_conv_pbernpoly_integral_Suc) 
         (auto simp: fun_eq_iff field_simps simp del: of_nat_Suc power_Suc)
    also have "convergent   convergent (?f (Suc (Suc n)))"
      by (intro convergent_mult_const_iff) (simp_all del: of_nat_Suc)
    finally show ?case .
  qed
  from this[of "n - 1"] assms show ?thesis by simp
qed

lemma stirling_integral_conv_stirling_integral_Suc:
  assumes "n > 0"
  shows   "stirling_integral n s =
             of_nat n / of_nat (Suc n) * stirling_integral (Suc n) s -
             of_real (bernoulli (Suc n)) / (of_nat (Suc n) * s ^ n)"
proof -
  have "(λN. of_real (pbernpoly (Suc n) (real N)) / (of_nat (Suc n) * (s + of_nat N) ^ n) -
             of_real (bernoulli (Suc n)) / (real (Suc n) * s ^ n) +
             integral {0..real N} (λx. of_nat n / of_nat (Suc n) * 
                (of_real (pbernpoly (Suc n) x) / (of_real x + s) ^ Suc n)))
            0 - of_real (bernoulli (Suc n)) / (of_nat (Suc n) * s ^ n) +
                   of_nat n / of_nat (Suc n) * stirling_integral (Suc n) s" (is "?f  _")
    unfolding stirling_integral_def integral_mult_right
    using convergent_stirling_integral[of "Suc n"] assms s
    by (intro tendsto_intros pbernpoly_over_power_tendsto_0)
       (auto simp: convergent_LIMSEQ_iff simp del: of_nat_Suc)
  also have "?this  (λN. integral {0..real N} 
               (λx. of_real (pbernpoly n x) / (of_real x + s) ^ n)) 
               of_nat n / of_nat (Suc n) * stirling_integral (Suc n) s -
                 of_real (bernoulli (Suc n)) / (of_nat (Suc n) * s ^ n)" 
    using eventually_gt_at_top[of "0::nat"] pbernpoly_integral_conv_pbernpoly_integral_Suc[of n] 
          assms unfolding integral_mult_right
    by (intro filterlim_cong refl) (auto elim!: eventually_mono simp del: power_Suc)
  finally show ?thesis unfolding stirling_integral_def[of n] by (rule limI)
qed

lemma stirling_integral_1_unfold:
  assumes "m > 0"
  shows   "stirling_integral 1 s = stirling_integral m s / of_nat m - 
             (k=1..<m. of_real (bernoulli (Suc k)) / (of_nat k * of_nat (Suc k) * s ^ k))"
proof -
  have "stirling_integral 1 s = stirling_integral (Suc m) s / of_nat (Suc m) -
          (k=1..<Suc m. of_real (bernoulli (Suc k)) / (of_nat k * of_nat (Suc k) * s ^ k))" for m
  proof (induction m)
    case (Suc m)
    let ?C = "(k = 1..<Suc m. of_real (bernoulli (Suc k)) / (of_nat k * of_nat (Suc k) * s ^ k))"
    note Suc.IH
    also have "stirling_integral (Suc m) s / of_nat (Suc m) = 
                 stirling_integral (Suc (Suc m)) s / of_nat (Suc (Suc m)) -
                 of_real (bernoulli (Suc (Suc m))) / 
                   (of_nat (Suc m) * of_nat (Suc (Suc m)) * s ^ Suc m)"
      (is "_ = ?A - ?B") by (subst stirling_integral_conv_stirling_integral_Suc)
                            (simp_all del: of_nat_Suc power_Suc add: divide_simps)
    also have "?A - ?B - ?C = ?A - (?B + ?C)" by (rule diff_diff_eq)
    also have "?B + ?C = (k = 1..<Suc (Suc m). of_real (bernoulli (Suc k)) /
                             (of_nat k * of_nat (Suc k) * s ^ k))" 
      using s by (simp add: divide_simps)
    finally show ?case .
  qed simp_all
  note this[of "m - 1"]
  also from assms have "Suc (m - 1) = m" by simp
  finally show ?thesis .
qed
  
lemma ln_Gamma_stirling_complex:
  assumes "m > 0"
  shows   "ln_Gamma s = (s - 1 / 2) * ln s - s + ln (2 * pi) / 2 +
             (k=1..<m. of_real (bernoulli (Suc k)) / (of_nat k * of_nat (Suc k) * s ^ k)) - 
             stirling_integral m s / of_nat m"
proof -
  have "ln_Gamma s = (s - 1 / 2) * ln s - s + ln (2 * pi) / 2 - stirling_integral 1 s"
    using limI[OF integral_pbernpoly_1] by (simp add: stirling_integral_def algebra_simps)
  also have "stirling_integral 1 s = stirling_integral m s / of_nat m -
               (k = 1..<m. of_real (bernoulli (Suc k)) / (of_nat k * of_nat (Suc k) * s ^ k))"
    using assms by (rule stirling_integral_1_unfold)
  finally show ?thesis by simp
qed

lemma LIMSEQ_stirling_integral:
  "n > 0  (λx. integral {0..real x} (λx. of_real (pbernpoly n x) / (of_real x + s) ^ n))
      stirling_integral n s" unfolding stirling_integral_def 
  using convergent_stirling_integral[of n] by (simp only: convergent_LIMSEQ_iff)

end

lemmas has_integral_of_real = has_integral_linear[OF _ bounded_linear_of_real, unfolded o_def]
lemmas integral_of_real = integral_linear[OF _ bounded_linear_of_real, unfolded o_def]

lemma integrable_ln_Gamma_aux_real:
  assumes "0 < s"
  shows   "(λx. pbernpoly n x / (x + s) ^ n) integrable_on {0..real N}"
proof -
  have "(λx. complex_of_real (pbernpoly n x / (x + s) ^ n)) integrable_on {0..real N}"
    using integrable_ln_Gamma_aux[of "of_real s" n N] assms by simp
  from integrable_linear[OF this bounded_linear_Re] show ?thesis 
    by (simp only: o_def Re_complex_of_real) 
qed
  
lemma  
  assumes "x > 0" "n > 0"
  shows   stirling_integral_complex_of_real:
            "stirling_integral n (complex_of_real x) = of_real (stirling_integral n x)"
    and   LIMSEQ_stirling_integral_real:
            "(λN. integral {0..real N} (λt. pbernpoly n t / (t + x) ^ n))
             stirling_integral n x"
    and   stirling_integral_real_convergent:
            "convergent (λN. integral {0..real N} (λt. pbernpoly n t / (t + x) ^ n))"
proof -
  have "(λN. integral {0..real N} (λt. of_real (pbernpoly n t / (t + x) ^ n)))
            stirling_integral n (complex_of_real x)"
    using LIMSEQ_stirling_integral[of "complex_of_real x" n] assms by simp
  hence "(λN. of_real (integral {0..real N} (λt. pbernpoly n t / (t + x) ^ n)))
            stirling_integral n (complex_of_real x)"
    using integrable_ln_Gamma_aux_real[OF assms(1), of n] 
    by (subst (asm) integral_of_real) simp
  from tendsto_Re[OF this] 
    have "(λN. integral {0..real N} (λt. pbernpoly n t / (t + x) ^ n))
            Re (stirling_integral n (complex_of_real x))" by simp
  thus "convergent (λN. integral {0..real N} (λt. pbernpoly n t / (t + x) ^ n))"
    by (rule convergentI)
  thus "(λN. integral {0..real N} (λt. pbernpoly n t / (t + x) ^ n))
           stirling_integral n x" unfolding stirling_integral_def
    by (simp add: convergent_LIMSEQ_iff)
  from tendsto_of_real[OF this, where 'a = complex] 
       integrable_ln_Gamma_aux_real[OF assms(1), of n]
    have "(λxa. integral {0..real xa} 
                    (λxa. complex_of_real (pbernpoly n xa) / (complex_of_real xa + x) ^ n))
              complex_of_real (stirling_integral n x)"
    by (subst (asm) integral_of_real [symmetric]) simp_all
  from LIMSEQ_unique[OF this LIMSEQ_stirling_integral[of "complex_of_real x" n]] assms
    show "stirling_integral n (complex_of_real x) = of_real (stirling_integral n x)" by simp
qed

lemma ln_Gamma_stirling_real:
  assumes "x > (0 :: real)" "m > (0::nat)"
  shows   "ln_Gamma x = (x - 1 / 2) * ln x - x + ln (2 * pi) / 2 +
              (k = 1..<m. bernoulli (Suc k) / (of_nat k * of_nat (Suc k) * x ^ k)) -
              stirling_integral m x / of_nat m"
proof -
  from assms have "complex_of_real (ln_Gamma x) = ln_Gamma (complex_of_real x)"
    by (simp add: ln_Gamma_complex_of_real)
  also have "ln_Gamma (complex_of_real x) = complex_of_real (
                (x - 1 / 2) * ln x - x + ln (2 * pi) / 2 +
                (k = 1..<m. bernoulli (Suc k) / (of_nat k * of_nat (Suc k) * x ^ k)) -
                stirling_integral m x / of_nat m)" using assms
    by (subst ln_Gamma_stirling_complex[of _ m])
       (simp_all add: Ln_of_real stirling_integral_complex_of_real)
  finally show ?thesis by (subst (asm) of_real_eq_iff)
qed


lemma stirling_integral_bound_aux:
  assumes n: "n > (1::nat)"
  obtains c where "s. Re s > 0  norm (stirling_integral n s)   c / Re s ^ (n - 1)"
proof -
  obtain c where c: "norm (pbernpoly n x)  c" for x by (rule bounded_pbernpoly[of n]) blast
  have c': "pbernpoly n x  c" for x using c[of x] by (simp add: abs_real_def split: if_splits)
  from c[of 0] have c_nonneg: "c  0" by simp
  have "norm (stirling_integral n s)  c / (real n - 1) / Re s ^ (n - 1)" if s: "Re s > 0" for s
  proof (rule Lim_norm_ubound[OF _ LIMSEQ_stirling_integral])
    have pos: "x + norm s > 0" if "x  0" for x using s that by (intro add_nonneg_pos) auto
    have nz: "of_real x + s  0" if "x  0" for x using s that by (auto simp: complex_eq_iff)
    let ?bound = "λN. c / (Re s ^ (n - 1) * (real n - 1)) - 
                        c / ((real N + Re s) ^ (n - 1) * (real n - 1))"
    show "eventually (λN. norm (integral {0..real N} 
              (λx. of_real (pbernpoly n x) / (of_real x + s) ^ n))  
            c / (real n - 1) / Re s ^ (n - 1)) at_top"
      using eventually_gt_at_top[of "0::nat"]
    proof eventually_elim
      case (elim N)
      let ?F = "λx. -c / ((x + Re s) ^ (n - 1) * (real n - 1))"
      from n s have "((λx. c / (x + Re s) ^ n) has_integral (?F (real N) - ?F 0)) {0..real N}"
        by (intro fundamental_theorem_of_calculus)
           (auto intro!: derivative_eq_intros simp: divide_simps power_diff add_eq_0_iff2
                   has_field_derivative_iff_has_vector_derivative [symmetric])      
      also have "?F (real N) - ?F 0 = ?bound N" by simp
      finally have *: "((λx. c / (x + Re s) ^ n) has_integral ?bound N) {0..real N}" .
      have "norm (integral {0..real N} (λx. of_real (pbernpoly n x) / (of_real x + s) ^ n)) 
              integral {0..real N} (λx. c / (x + Re s) ^ n)"
      proof (intro integral_norm_bound_integral integrable_ln_Gamma_aux s ballI)
        fix x assume x: "x  {0..real N}"
        have "norm (of_real (pbernpoly n x) / (of_real x + s) ^ n)  c / norm (of_real x + s) ^ n"
          unfolding norm_divide norm_power using c by (intro divide_right_mono) simp_all
        also have "  c / (x + Re s) ^ n" 
          using x c c_nonneg s nz[of x] complex_Re_le_cmod[of "of_real x + s"]
          by (intro divide_left_mono power_mono mult_pos_pos zero_less_power add_nonneg_pos) auto
        finally show "norm (of_real (pbernpoly n x) / (of_real x + s) ^ n)  " .
      qed (insert n s * pos nz c, auto simp: complex_nonpos_Reals_iff)
      also have " = ?bound N" using * by (simp add: has_integral_iff)
      also have "  c / (Re s ^ (n - 1) * (real n - 1))" using c_nonneg elim s n by simp
      also have " = c / (real n - 1) / (Re s ^ (n - 1))" by simp
      finally show "norm (integral {0..real N} (λx. of_real (pbernpoly n x) /
                      (of_real x + s) ^ n))  c / (real n - 1) / Re s ^ (n - 1)" .
    qed
  qed (insert s n, simp_all add: complex_nonpos_Reals_iff)
  thus ?thesis by (rule that)
qed

lemma stirling_integral_bound_aux_integral1:
  fixes a b c :: real and n :: nat
  assumes "a  0" "b > 0" "c  0" "n > 1" "l < a - b" "r > a + b"
  shows "((λx. c / max b ¦x - a¦ ^ n) has_integral
           2*c*(n / (n - 1))/b^(n-1) - c/(n-1) * (1/(a-l)^(n-1) + 1/(r-a)^(n-1))) {l..r}"
proof -
  define x1 x2 where "x1 = a - b" and "x2 = a + b"
  define F1 where "F1 = (λx::real. c / (a - x) ^ (n - 1) / (n - 1))"
  define F3 where "F3 = (λx::real. -c / (x - a) ^ (n - 1) / (n - 1))"
  have deriv: "(F1 has_vector_derivative (c / (a - x) ^ n)) (at x within A)"
              "(F3 has_vector_derivative (c / (x - a) ^ n)) (at x within A)"
    if "x  a" for x :: real and A
    unfolding F1_def F3_def using assms that
    by (auto intro!: derivative_eq_intros simp: divide_simps power_diff add_eq_0_iff2
             simp flip: has_field_derivative_iff_has_vector_derivative)

  from assms have "((λx. c / (a - x) ^ n) has_integral (F1 x1 - F1 l)) {l..x1}"
    by (intro fundamental_theorem_of_calculus deriv) (auto simp: x1_def max_def split: if_splits)
  also have "?this  ((λx. c / max b ¦x - a¦ ^ n) has_integral (F1 x1 - F1 l)) {l..x1}"
    using assms
    by (intro has_integral_spike_finite_eq[of "{l}"]) (auto simp: x1_def max_def split: if_splits)
  finally have I1: "((λx. c / max b ¦x - a¦ ^ n) has_integral (F1 x1 - F1 l)) {l..x1}" .

  have "((λx. c / b ^ n) has_integral (x2 - x1) * c / b ^ n) {x1..x2}"
    using has_integral_const_real[of "c / b ^ n" x1 x2] assms by (simp add: x1_def x2_def)
  also have "?this  ((λx. c / max b ¦x - a¦ ^ n) has_integral ((x2 - x1) * c / b ^ n)) {x1..x2}"
    by (intro has_integral_spike_finite_eq[of "{x1, x2}"])
       (auto simp: x1_def x2_def split: if_splits)
  finally have I2: "((λx. c / max b ¦x - a¦ ^ n) has_integral ((x2 - x1) * c / b ^ n)) {x1..x2}" .

  from assms have I3: "((λx. c / (x - a) ^ n) has_integral (F3 r - F3 x2)) {x2..r}"
    by (intro fundamental_theorem_of_calculus deriv) (auto simp: x2_def min_def split: if_splits)
  also have "?this  ((λx. c / max b ¦x - a¦ ^ n) has_integral (F3 r - F3 x2)) {x2..r}"
    using assms
    by (intro has_integral_spike_finite_eq[of "{r}"]) (auto simp: x2_def min_def split: if_splits)
  finally have I3: "((λx. c / max b ¦x - a¦ ^ n) has_integral (F3 r - F3 x2)) {x2..r}" .

  have "((λx. c / max b ¦x - a¦ ^ n) has_integral (F1 x1 - F1 l) + ((x2 - x1) * c / b ^ n) + (F3 r - F3 x2)) {l..r}"
    using assms
    by (intro has_integral_combine[OF _ _ has_integral_combine[OF _ _ I1 I2] I3])
       (auto simp: x1_def x2_def)
  also have "(F1 x1 - F1 l) + ((x2 - x1) * c / b ^ n) + (F3 r - F3 x2) =
             F1 x1 - F1 l + F3 r - F3 x2 + (x2 - x1) * c / b ^ n"
    by (simp add: algebra_simps)
  also have "x2 - x1 = 2 * b"
    using assms by (simp add: x2_def x1_def min_def max_def)
  also have "2 * b * c / b ^ n = 2 * c / b ^ (n - 1)"
    using assms by (simp add: power_diff field_simps)
  also have "F1 x1 - F1 l + F3 r - F3 x2 =
               c/(n-1) * (2/b^(n-1) - 1/(a-l)^(n-1) - 1/(r-a)^(n-1))"
    using assms by (simp add: x1_def x2_def F1_def F3_def field_simps)
  also have " + 2 * c / b ^ (n - 1) =
             2*c*(1 + 1/(n-1))/b^(n-1) - c/(n-1) * (1/(a-l)^(n-1) + 1/(r-a)^(n-1))"
    using assms by (simp add: field_simps)
  also have "1 + 1 / (n - 1) = n / (n - 1)"
    using assms by (simp add: field_simps)
  finally show ?thesis .
qed

lemma stirling_integral_bound_aux_integral2:
  fixes a b c :: real and n :: nat
  assumes "a  0" "b > 0" "c  0" "n > 1"
  obtains I where "((λx. c / max b ¦x - a¦ ^ n) has_integral I) {l..r}"
                  "I  2 * c * (n / (n - 1)) / b ^ (n-1)"
proof -
  define l' where "l' = min l (a - b - 1)"
  define r' where "r' = max r (a + b + 1)"

  define A where "A = 2 * c * (n / (n - 1)) / b ^ (n - 1)"
  define B where "B = c / real (n - 1) * (1 / (a - l') ^ (n - 1) + 1 / (r' - a) ^ (n - 1))"

  have has_int: "((λx. c / max b ¦x - a¦ ^ n) has_integral (A - B)) {l'..r'}"
    using assms unfolding A_def B_def
    by (intro stirling_integral_bound_aux_integral1) (auto simp: l'_def r'_def)
  have "(λx. c / max b ¦x - a¦ ^ n) integrable_on {l..r}"
    by (rule integrable_on_subinterval[OF has_integral_integrable[OF has_int]])
       (auto simp: l'_def r'_def)
  then obtain I where has_int': "((λx. c / max b ¦x - a¦ ^ n) has_integral I) {l..r}"
    by (auto simp: integrable_on_def)

  from assms have "I  A - B"
    by (intro has_integral_subset_le[OF _ has_int' has_int]) (auto simp: l'_def r'_def)
  also have "  A"
    using assms by (simp add: B_def l'_def r'_def)
  finally show ?thesis using that[of I] has_int' unfolding A_def by blast
qed

lemma stirling_integral_bound_aux':
  assumes n: "n > (1::nat)" and α: "α  {0<..<pi}"
  obtains c where "s::complex. s  complex_cone' α - {0} 
                     norm (stirling_integral n s)  c / norm s ^ (n - 1)"
proof -
  obtain c where c: "norm (pbernpoly n x)  c" for x by (rule bounded_pbernpoly[of n]) blast
  have c': "pbernpoly n x  c" for x using c[of x] by (simp add: abs_real_def split: if_splits)
  from c[of 0] have c_nonneg: "c  0" by simp

  define D where "D = c * Beta (- (real_of_int (- int n) / 2) - 1 / 2) (1 / 2) / 2"
  define C where "C = max D (2*c*(n/(n-1))/sin α^(n-1))"

  have *: "norm (stirling_integral n s)  C / norm s ^ (n - 1)"
    if s: "s  complex_cone' α - {0}" for s :: complex
  proof (rule Lim_norm_ubound[OF _ LIMSEQ_stirling_integral])
    from s α have arg: "¦arg s¦  α" by (auto simp: complex_cone_altdef)
    have s': "s  0"
      using complex_cone_inter_nonpos_Reals[of "-α" α] α s  by auto
    from s have [simp]: "s  0" by auto

    show "eventually (λN. norm (integral {0..real N}
              (λx. of_real (pbernpoly n x) / (of_real x + s) ^ n))  
            C / norm s ^ (n - 1)) at_top"
      using eventually_gt_at_top[of "0::nat"]
    proof eventually_elim
      case (elim N)
      show ?case
      proof (cases "Re s > 0")
        case True
        have int: "((λx. c * (x^2 + norm s^2) powr (-n / 2)) has_integral
                  D * (norm s ^ 2) powr (-n / 2 + 1 / 2)) {0<..}"
          using has_integral_mult_left[OF has_integral_Beta3[of "-n/2" "norm s ^ 2"], of c] assms True
          unfolding D_def by (simp add: algebra_simps)
        hence int': "((λx. c * (x^2 + norm s^2) powr (-n / 2)) has_integral
                  D * (norm s ^ 2) powr (-n / 2 + 1 / 2)) {0..}"
          by (subst has_integral_interior [symmetric]) simp_all
        hence integrable: "(λx. c * (x^2 + norm s^2) powr (-n / 2)) integrable_on {0..}"
          by (simp add: has_integral_iff)

        have "norm (integral {0..real N} (λx. of_real (pbernpoly n x) / (of_real x + s) ^ n)) 
                integral {0..real N} (λx. c * (x^2 + norm s^2) powr (-n / 2))"
        proof (intro integral_norm_bound_integral s ballI integrable_ln_Gamma_aux)
          have [simp]: "{0<..} - {0::real..} = {}" "{0..} - {0<..} = {0::real}"
            by auto
          have "(λx. c * (x2 + (cmod s)2) powr (real_of_int (- int n) / 2)) integrable_on {0<..}"
            using int by (simp add: has_integral_iff)
          also have "?this  (λx. c * (x2 + (cmod s)2) powr (real_of_int (- int n) / 2)) integrable_on {0..}"
            by (intro integrable_spike_set_eq) auto
          finally show "(λx. c * (x2 + (cmod s)2) powr (real_of_int (- int n) / 2)) integrable_on
                   {0..real N}" by (rule integrable_on_subinterval) auto
        next
          fix x assume x: "x  {0..real N}"
          have nz: "complex_of_real x + s  0"
            using True x by (auto simp: complex_eq_iff)
          have "norm (of_real (pbernpoly n x) / (of_real x + s) ^ n)  c / norm (of_real x + s) ^ n"
            unfolding norm_divide norm_power using c by (intro divide_right_mono) simp_all
          also have "  c / sqrt (x ^ 2 + norm s ^ 2) ^ n"
          proof (intro divide_left_mono mult_pos_pos zero_less_power power_mono)
            show "sqrt (x2 + (cmod s)2)  cmod (complex_of_real x + s)"
              using x True by (simp add: cmod_def algebra_simps power2_eq_square)
          qed (use x True c_nonneg assms nz in auto simp: add_nonneg_pos›)
          also have "sqrt (x ^ 2 + norm s ^ 2) ^ n = (x ^ 2 + norm s ^ 2) powr (1/2 * n)"
            by (subst powr_powr [symmetric], subst powr_realpow)
               (auto simp: powr_half_sqrt add_nonneg_pos)
          also have "c /  = c * (x^2 + norm s^2) powr (-n / 2)"
            by (simp add: powr_minus field_simps)
          finally show "norm (complex_of_real (pbernpoly n x) / (complex_of_real x + s) ^ n)  " .
        qed fact+
        also have "  integral {0..} (λx. c * (x^2 + norm s^2) powr (-n / 2))"
          using c_nonneg
          by (intro integral_subset_le integrable integrable_on_subinterval[OF integrable]) auto
        also have " = D * (norm s ^ 2) powr (-n / 2 + 1 / 2)"
          using int' by (simp add: has_integral_iff)
        also have "(norm s ^ 2) powr (-n / 2 + 1 / 2) = norm s powr (2 * (-n / 2 + 1 / 2))"
          by (subst powr_powr [symmetric]) auto
        also have " = norm s powr (-real (n - 1))"
          using assms by (simp add: of_nat_diff)
        also have "D *  = D / norm s ^ (n - 1)"
          by (auto simp: powr_minus powr_realpow field_simps)
        also have "  C / norm s ^ (n - 1)"
          by (intro divide_right_mono) (auto simp: C_def)
        finally show "norm (integral {0..real N} (λx. of_real (pbernpoly n x) / (of_real x + s) ^ n))  " .

      next

        case False
        have "cos ¦arg s¦ = cos (arg s)"
          by (simp add: abs_if)
        also have "cos (arg s) = Re (rcis (norm s) (arg s)) / norm s"
          by (subst Re_rcis) auto
        also have " = Re s / norm s"
          by (subst rcis_cmod_arg) auto
        also have "  cos (pi / 2)"
          using False by (auto simp: field_simps)
        finally have "¦arg s¦  pi / 2"
          using arg α by (subst (asm) cos_mono_le_eq) auto

        have "sin α * norm s = sin (pi - α) * norm s"
          by simp
        also have "  sin (pi - ¦arg s¦) * norm s"
          using α arg ¦arg s¦  pi / 2
          by (intro mult_right_mono sin_monotone_2pi_le) auto
        also have "sin ¦arg s¦  0"
          using arg_bounded[of s] by (intro sin_ge_zero) auto
        hence "sin (pi - ¦arg s¦) = ¦sin ¦arg s¦¦"
          by simp 
        also have " = ¦sin (arg s)¦"
          by (simp add: abs_if)
        also have " * norm s = ¦Im (rcis (norm s) (arg s))¦"
          by (simp add: abs_mult)
        also have " = ¦Im s¦"
          by (subst rcis_cmod_arg) auto
        finally have abs_Im_ge: "¦Im s¦  sin α * norm s" .

        have [simp]: "Im s  0" "s  0"
          using s s  0 False
          by (auto simp: cmod_def zero_le_mult_iff complex_nonpos_Reals_iff)
        have "sin α > 0"
          using assms by (intro sin_gt_zero) auto
  
        obtain I where I: "((λx. c / max ¦Im s¦ ¦x + Re s¦ ^ n) has_integral I) {0..real N}"
                          "I  2*c*(n/(n-1)) / ¦Im s¦ ^ (n - 1)"
          using s c_nonneg assms False 
                stirling_integral_bound_aux_integral2[of "-Re s" "¦Im s¦" c n 0 "real N"] by auto
  
        have "norm (integral {0..real N} (λx. of_real (pbernpoly n x) / (of_real x + s) ^ n)) 
                integral {0..real N} (λx. c / max ¦Im s¦ ¦x + Re s¦ ^ n)"
        proof (intro integral_norm_bound_integral integrable_ln_Gamma_aux s ballI)
          show "(λx. c / max ¦Im s¦ ¦x + Re s¦ ^ n) integrable_on {0..real N}"
            using I(1) by (simp add: has_integral_iff)
        next
          fix x assume x: "x  {0..real N}"
          have nz: "complex_of_real x + s  0"
            by (auto simp: complex_eq_iff)
          have "norm (complex_of_real (pbernpoly n x) / (complex_of_real x + s) ^ n) 
                  c / norm (complex_of_real x + s) ^ n"
            unfolding norm_divide norm_power using c[of x] by (intro divide_right_mono) simp_all
          also have "  c / max ¦Im s¦ ¦x + Re s¦ ^ n"
            using c_nonneg nz abs_Re_le_cmod[of "of_real x + s"] abs_Im_le_cmod[of "of_real x + s"]
            by (intro divide_left_mono power_mono mult_pos_pos zero_less_power)
               (auto simp: less_max_iff_disj)
          finally show "norm (complex_of_real (pbernpoly n x) / (complex_of_real x + s) ^ n)  " .
        qed (auto simp: complex_nonpos_Reals_iff)
        also have "  2*c*(n/(n-1)) / ¦Im s¦ ^ (n - 1)"
          using I by (simp add: has_integral_iff)
        also have "  2*c*(n/(n-1)) / (sin α * norm s) ^ (n - 1)"
          using ‹sin α > 0 s c_nonneg abs_Im_ge
          by (intro divide_left_mono mult_pos_pos zero_less_power power_mono mult_nonneg_nonneg) auto
        also have " = 2*c*(n/(n-1))/sin α^(n-1) / norm s ^ (n - 1)"
          by (simp add: field_simps)
        also have "  C / norm s ^ (n - 1)"
          by (intro divide_right_mono) (auto simp: C_def)
        finally show ?thesis .
      qed
    qed
  qed (use that assms complex_cone_inter_nonpos_Reals[of "-α" α] α in auto)
  thus ?thesis by (rule that)
qed

lemma stirling_integral_bound:
  assumes "n > 0"
  obtains c where 
    "s. Re s > 0  norm (stirling_integral n s)  c / Re s ^ n"
proof -
  let ?f = "λs. of_nat n / of_nat (Suc n) * stirling_integral (Suc n) s -
                  of_real (bernoulli (Suc n)) / (of_nat (Suc n) * s ^ n)"
  from stirling_integral_bound_aux[of "Suc n"] assms obtain c where 
    c: "s. Re s > 0  norm (stirling_integral (Suc n) s)  c / Re s ^ n" by auto
  define c1 where "c1 = real n / real (Suc n) * c"
  define c2 where "c2 = ¦bernoulli (Suc n)¦ / real (Suc n)"
  have c2_nonneg: "c2  0" by (simp add: c2_def)
  show ?thesis
  proof (rule that)
    fix s :: complex assume s: "Re s > 0"
    hence s': "s  0" by (auto simp: complex_nonpos_Reals_iff)
    have "stirling_integral n s = ?f s" using s' assms 
      by (rule stirling_integral_conv_stirling_integral_Suc)
    also have "norm   norm (of_nat n / of_nat (Suc n) * stirling_integral (Suc n) s) +
                           norm (of_real (bernoulli (Suc n)) / (of_nat (Suc n) * s ^ n))"
      by (rule norm_triangle_ineq4)
    also have " = real n / real (Suc n) * norm (stirling_integral (Suc n) s) +
                      c2 / norm s ^ n" (is "_ = ?A + ?B")
      by (simp add: norm_divide norm_mult norm_power c2_def field_simps del: of_nat_Suc)
    also have "?A  real n / real (Suc n) * (c / Re s ^ n)"
      by (intro mult_left_mono c s) simp_all
    also have " = c1 / Re s ^ n" by (simp add: c1_def)
    also have "c2 / norm s ^ n  c2 / Re s ^ n" using s c2_nonneg
      by (intro divide_left_mono power_mono complex_Re_le_cmod mult_pos_pos zero_less_power) auto
    also have "c1 / Re s ^ n + c2 / Re s ^ n = (c1 + c2) / Re s ^ n" 
      using s by (simp add: field_simps)
    finally show "norm (stirling_integral n s)  (c1 + c2) / Re s ^ n" by - simp_all
  qed
qed

lemma stirling_integral_bound':
  assumes "n > 0" and "α  {0<..<pi}"
  obtains c where 
    "s::complex. s  complex_cone' α - {0}  norm (stirling_integral n s)  c / norm s ^ n"
proof -
  let ?f = "λs. of_nat n / of_nat (Suc n) * stirling_integral (Suc n) s -
                  of_real (bernoulli (Suc n)) / (of_nat (Suc n) * s ^ n)"
  from stirling_integral_bound_aux'[of "Suc n"] assms obtain c where 
    c: "s::complex. s  complex_cone' α - {0} 
            norm (stirling_integral (Suc n) s)  c / norm s ^ n" by auto
  define c1 where "c1 = real n / real (Suc n) * c"
  define c2 where "c2 = ¦bernoulli (Suc n)¦ / real (Suc n)"
  have c2_nonneg: "c2  0" by (simp add: c2_def)
  show ?thesis
  proof (rule that)
    fix s :: complex assume s: "s  complex_cone' α - {0}"
    have s': "s  0"
      using complex_cone_inter_nonpos_Reals[of "-α" α] assms s by auto
      
    have "stirling_integral n s = ?f s" using s' assms 
      by (intro stirling_integral_conv_stirling_integral_Suc) auto
    also have "norm   norm (of_nat n / of_nat (Suc n) * stirling_integral (Suc n) s) +
                           norm (of_real (bernoulli (Suc n)) / (of_nat (Suc n) * s ^ n))"
      by (rule norm_triangle_ineq4)
    also have " = real n / real (Suc n) * norm (stirling_integral (Suc n) s) +
                      c2 / norm s ^ n" (is "_ = ?A + ?B")
      by (simp add: norm_divide norm_mult norm_power c2_def field_simps del: of_nat_Suc)
    also have "?A  real n / real (Suc n) * (c / norm s ^ n)"
      by (intro mult_left_mono c s) simp_all
    also have " = c1 / norm s ^ n" by (simp add: c1_def)
    also have "c1 / norm s ^ n + c2 / norm s ^ n = (c1 + c2) / norm s ^ n" 
      using s by (simp add: divide_simps)
    finally show "norm (stirling_integral n s)  (c1 + c2) / norm s ^ n" by - simp_all
  qed
qed


lemma stirling_integral_holomorphic [holomorphic_intros]:
  assumes m: "m > 0" and "A  0 = {}"
  shows   "stirling_integral m holomorphic_on A"  
proof -
  from assms have [simp]: "z  0" if "z  A" for z
    using that by auto
  let ?f = "λs::complex. of_nat m * ((s - 1 / 2) * Ln s - s + of_real (ln (2 * pi) / 2) +
          (k=1..<m. of_real (bernoulli (Suc k)) / (of_nat k * of_nat (Suc k) * s ^ k)) - 
          ln_Gamma s)"
  have "?f holomorphic_on A" using assms
    by (auto intro!: holomorphic_intros simp del: of_nat_Suc elim!: nonpos_Reals_cases)
  also have "?this  stirling_integral m holomorphic_on A" 
    using assms by (intro holomorphic_cong refl) 
                   (simp_all add: field_simps ln_Gamma_stirling_complex)
  finally show "stirling_integral m holomorphic_on A" .
qed

lemma stirling_integral_continuous_on_complex [continuous_intros]:
  assumes m: "m > 0" and "A  0 = {}"
  shows   "continuous_on A (stirling_integral m :: _  complex)"
  by (intro holomorphic_on_imp_continuous_on stirling_integral_holomorphic assms)
    
lemma has_field_derivative_stirling_integral_complex:
  fixes x :: complex
  assumes "x  0" "n > 0"
  shows   "(stirling_integral n has_field_derivative deriv (stirling_integral n) x) (at x)"
  using assms
  by (intro holomorphic_derivI[OF stirling_integral_holomorphic, of n  "-0"]) auto


 
lemma
  assumes n: "n > 0" and "x > 0"
  shows   deriv_stirling_integral_complex_of_real:
            "(deriv ^^ j) (stirling_integral n) (complex_of_real x) =
               complex_of_real ((deriv ^^ j) (stirling_integral n) x)" (is "?lhs x = ?rhs x")
  and     differentiable_stirling_integral_real:
            "(deriv ^^ j) (stirling_integral n) field_differentiable at x" (is ?thesis2)
proof -
  let ?A = "{s. Re s > 0}"
  let ?f = "λj x. (deriv ^^ j) (stirling_integral n) (complex_of_real x)"
  let ?f' = "λj x. complex_of_real ((deriv ^^ j) (stirling_integral n) x)"
    
  have [simp]: "open ?A" by (simp add: open_halfspace_Re_gt)      

  have "?lhs x = ?rhs x  (deriv ^^ j) (stirling_integral n) field_differentiable at x" 
    if "x > 0" for x using that
  proof (induction j arbitrary: x)
    case 0
    have "((λx. Re (stirling_integral n (of_real x))) has_field_derivative 
                  Re (deriv (λx. stirling_integral n x) (of_real x))) (at x)" using 0 n
      by (auto intro!: derivative_intros has_vector_derivative_real_field
                 field_differentiable_derivI holomorphic_on_imp_differentiable_at[of _ ?A]
                 stirling_integral_holomorphic simp: complex_nonpos_Reals_iff)
    also have "?this  (stirling_integral n has_field_derivative 
             Re (deriv (λx. stirling_integral n x) (of_real x))) (at x)"
      using eventually_nhds_in_open[of "{0<..}" x] 0 n
      by (intro has_field_derivative_cong_ev refl) 
         (auto elim!: eventually_mono simp: stirling_integral_complex_of_real)
    finally have "stirling_integral n field_differentiable at x"
      by (auto simp: field_differentiable_def)
    with 0 n show ?case by (auto simp: stirling_integral_complex_of_real)
  next
    case (Suc j x)
    note IH = conjunct1[OF Suc.IH] conjunct2[OF Suc.IH]
    have *: "(deriv ^^ Suc j) (stirling_integral n) (complex_of_real x) =
                 of_real ((deriv ^^ Suc j) (stirling_integral n) x)" if x: "x > 0" for x
    proof -
      have "deriv ((deriv ^^ j) (stirling_integral n)) (complex_of_real x) =
              vector_derivative (λx. (deriv ^^ j) (stirling_integral n) (of_real x)) (at x)"
        using n x
        by (intro vector_derivative_of_real_right [symmetric] 
                   holomorphic_on_imp_differentiable_at[of _ ?A] holomorphic_higher_deriv
                   stirling_integral_holomorphic) (auto simp: complex_nonpos_Reals_iff)
      also have " = vector_derivative (λx. of_real ((deriv ^^ j) (stirling_integral n) x)) (at x)"
        using eventually_nhds_in_open[of "{0<..}" x] x
        by (intro vector_derivative_cong_eq) (auto elim!: eventually_mono simp: IH(1))
      also have " = of_real (deriv ((deriv ^^ j) (stirling_integral n)) x)"
        by (intro vector_derivative_of_real_left holomorphic_on_imp_differentiable_at[of _ ?A]
              field_differentiable_imp_differentiable IH(2) x)
      finally show ?thesis by simp
    qed
    have "((λx. Re ((deriv ^^ Suc j) (stirling_integral n) (of_real x))) has_field_derivative 
             Re (deriv ((deriv ^^ Suc j) (stirling_integral n)) (of_real x))) (at x)"
      using Suc.prems n
      by (intro derivative_intros has_vector_derivative_real_field field_differentiable_derivI
                holomorphic_on_imp_differentiable_at[of _ ?A] stirling_integral_holomorphic
                holomorphic_higher_deriv) (auto simp: complex_nonpos_Reals_iff)
    also have "?this  ((deriv ^^ Suc j) (stirling_integral n) has_field_derivative 
                  Re (deriv ((deriv ^^ Suc j) (stirling_integral n)) (of_real x))) (at x)"  
      using eventually_nhds_in_open[of "{0<..}" x] Suc.prems *
      by (intro has_field_derivative_cong_ev refl) (auto elim!: eventually_mono)
    finally have "(deriv ^^ Suc j) (stirling_integral n) field_differentiable at x"
      by (auto simp: field_differentiable_def)
    with *[OF Suc.prems] show ?case by blast
  qed
  from this[OF assms(2)] show "?lhs x = ?rhs x" ?thesis2 by blast+
qed

text ‹
  Unfortunately, asymptotic power series cannot, in general, be differentiated. However, since 
  @{term ln_Gamma} is holomorphic on the entire positive real half-space, we can differentiate 
  its asymptotic expansion after all.

  To do this, we use an ad-hoc version of the more general approach outlined in Erdelyi's
  ``Asymptotic Expansions'' for holomorphic functions: We bound the value of the $j$-th derivative 
  of the remainder term at some value $x$ by applying Cauchy's integral formula along a circle 
  centred at $x$ with radius $\frac{1}{2} x$.
›
lemma deriv_stirling_integral_real_bound:
  assumes m: "m > 0"
  shows   "(deriv ^^ j) (stirling_integral m)  O(λx::real. 1 / x ^ (m + j))"
proof -
  from stirling_integral_bound[OF m] guess c . note c = this
  have "0  cmod (stirling_integral m 1)" by simp
  also have "  c" using c[of 1] by simp
  finally have c_nonneg: "c  0" .
  define B where "B = c * 2 ^ (m + Suc j)"
  define B' where "B' = B * fact j / 2"

  have "eventually (λx::real. norm ((deriv ^^ j) (stirling_integral m) x)  
          B' * norm (1 / x ^ (m+ j))) at_top"
    using eventually_gt_at_top[of "0::real"]
  proof eventually_elim
    case (elim x)
    have "s  0" if "s  cball (of_real x) (x/2)" for s :: complex
    proof -
      have "x - Re s  norm (of_real x - s)" using complex_Re_le_cmod[of "of_real x - s"] by simp
      also from that have "  x/2" by (simp add: dist_complex_def)
      finally show ?thesis using elim by (auto simp: complex_nonpos_Reals_iff)
    qed
    hence "((λu. stirling_integral m u / (u - of_real x) ^ Suc j) has_contour_integral
            complex_of_real (2 * pi) * 𝗂 / fact j * 
              (deriv ^^ j) (stirling_integral m) (of_real x)) (circlepath (of_real x) (x/2))"
      using m elim
      by (intro Cauchy_has_contour_integral_higher_derivative_circlepath 
                stirling_integral_continuous_on_complex stirling_integral_holomorphic) auto
    hence "norm (of_real (2 * pi) * 𝗂 / fact j * (deriv ^^ j) (stirling_integral m) (of_real x)) 
            B / x ^ (m + Suc j) * (2 * pi * (x / 2))"
    proof (rule has_contour_integral_bound_circlepath)
      fix u :: complex assume dist: "norm (u - of_real x) = x / 2"
      have "Re (of_real x - u)  norm (of_real x - u)" by (rule complex_Re_le_cmod)
      also have " = x / 2" using dist by (simp add: norm_minus_commute)
      finally have Re_u: "Re u  x/2" using elim by simp
      have "norm (stirling_integral m u / (u - of_real x) ^ Suc j)  
              c / Re u ^ m / (x / 2) ^ Suc j" using Re_u elim
        unfolding norm_divide norm_power dist
        by (intro divide_right_mono zero_le_power c) simp_all
      also have "  c / (x/2) ^ m / (x / 2) ^ Suc j" using c_nonneg elim Re_u
        by (intro divide_right_mono divide_left_mono power_mono) simp_all
      also have " = B / x ^ (m + Suc j)" using elim by (simp add: B_def field_simps power_add)
      finally show "norm (stirling_integral m u / (u - of_real x) ^ Suc j)  B / x ^ (m + Suc j)" .
    qed (insert elim c_nonneg, auto simp: B_def simp del: power_Suc)
    hence "cmod ((deriv ^^ j) (stirling_integral m) (of_real x))  B' / x ^ (j + m)"
      using elim by (simp add: field_simps norm_divide norm_mult norm_power B'_def)
    with elim m show ?case by (simp_all add: add_ac deriv_stirling_integral_complex_of_real)
  qed
  thus ?thesis by (rule bigoI)
qed

definition stirling_sum where
  "stirling_sum j m x = 
     (-1) ^ j * (k = 1..<m. (of_real (bernoulli (Suc k)) * pochhammer (of_nat k) j / (of_nat k *
                                 of_nat (Suc k))) * inverse x ^ (k + j))"
  
definition stirling_sum' where
  "stirling_sum' j m x = 
     (-1) ^ (Suc j) * (km. (of_real (bernoulli' k) * 
       pochhammer (of_nat (Suc k)) (j - 1) * inverse x ^ (k + j)))"

lemma stirling_sum_complex_of_real:
  "stirling_sum j m (complex_of_real x) = complex_of_real (stirling_sum j m x)"
  by (simp add: stirling_sum_def pochhammer_of_real [symmetric] del: of_nat_Suc)

lemma stirling_sum'_complex_of_real:
  "stirling_sum' j m (complex_of_real x) = complex_of_real (stirling_sum' j m x)"
  by (simp add: stirling_sum'_def pochhammer_of_real [symmetric] del: of_nat_Suc)

lemma has_field_derivative_stirling_sum_complex [derivative_intros]:
  "Re x > 0  (stirling_sum j m has_field_derivative stirling_sum (Suc j) m x) (at x)"
    unfolding stirling_sum_def [abs_def] sum_distrib_left
    by (rule DERIV_sum) (auto intro!: derivative_eq_intros simp del: of_nat_Suc 
                              simp: pochhammer_Suc power_diff)

lemma has_field_derivative_stirling_sum_real [derivative_intros]:
  "x > (0::real)  (stirling_sum j m has_field_derivative stirling_sum (Suc j) m x) (at x)"
    unfolding stirling_sum_def [abs_def] sum_distrib_left
    by (rule DERIV_sum) (auto intro!: derivative_eq_intros simp del: of_nat_Suc 
                              simp: pochhammer_Suc power_diff)

lemma has_field_derivative_stirling_sum'_complex [derivative_intros]:
  assumes "j > 0" "Re x > 0"
  shows   "(stirling_sum' j m has_field_derivative stirling_sum' (Suc j) m x) (at x)"
proof (cases j)
  case (Suc j')
  from assms have [simp]: "x  0" by auto
  define c where "c = (λn. (-1) ^ Suc j * complex_of_real (bernoulli' n) * 
                          pochhammer (of_nat (Suc n)) j')"
  define T where "T = (λn x. c n * inverse x ^ (j + n))"
  define T' where "T' = (λn x. - (of_nat (j + n)) * c n * inverse x ^ (Suc (j + n)))"
  have "((λx. km. T k x) has_field_derivative (km. T' k x)) (at x)" using assms Suc
    by (intro DERIV_sum)
       (auto simp: T_def T'_def intro!: derivative_eq_intros 
             simp: field_simps power_add [symmetric]  simp del: of_nat_Suc power_Suc of_nat_add)
  also have "(λx. (km. T k x)) = stirling_sum' j m"
    by (simp add: Suc T_def c_def stirling_sum'_def fun_eq_iff add_ac mult.assoc sum_distrib_left)
  also have "(km. T' k x) = stirling_sum' (Suc j) m x"
    by (simp add: T'_def c_def Suc stirling_sum'_def sum_distrib_left 
          sum_distrib_right algebra_simps pochhammer_Suc)
  finally show ?thesis .
qed (insert assms, simp_all)
  
lemma has_field_derivative_stirling_sum'_real [derivative_intros]:
  assumes "j > 0" "x > (0::real)"
  shows   "(stirling_sum' j m has_field_derivative stirling_sum' (Suc j) m x) (at x)"
proof (cases j)
  case (Suc j')
  from assms have [simp]: "x  0" by auto
  define c where "c = (λn. (-1) ^ Suc j * (bernoulli' n) * pochhammer (of_nat (Suc n)) j')"
  define T where "T = (λn x. c n * inverse x ^ (j + n))"
  define T' where "T' = (λn x. - (of_nat (j + n)) * c n * inverse x ^ (Suc (j + n)))"
  have "((λx. km. T k x) has_field_derivative (km. T' k x)) (at x)" using assms Suc
    by (intro DERIV_sum)
       (auto simp: T_def T'_def intro!: derivative_eq_intros 
             simp: field_simps power_add [symmetric]  simp del: of_nat_Suc power_Suc of_nat_add)
  also have "(λx. (km. T k x)) = stirling_sum' j m"
    by (simp add: Suc T_def c_def stirling_sum'_def fun_eq_iff add_ac mult.assoc sum_distrib_left)
  also have "(km. T' k x) = stirling_sum' (Suc j) m x"
    by (simp add: T'_def c_def Suc stirling_sum'_def sum_distrib_left 
          sum_distrib_right algebra_simps pochhammer_Suc)
  finally show ?thesis .
qed (insert assms, simp_all)

lemma higher_deriv_stirling_sum_complex:
  "Re x > 0  (deriv ^^ i) (stirling_sum j m) x = stirling_sum (i + j) m x"
proof (induction i arbitrary: x)
  case (Suc i)
  have "deriv ((deriv ^^ i) (stirling_sum j m)) x = deriv (stirling_sum (i + j) m) x"
    using eventually_nhds_in_open[of "{x. Re x > 0}" x] Suc.prems
    by (intro deriv_cong_ev refl) (auto elim!: eventually_mono simp: open_halfspace_Re_gt Suc.IH)
  also from Suc.prems have " = stirling_sum (Suc (i + j)) m x"
    by (intro DERIV_imp_deriv has_field_derivative_stirling_sum_complex)
  finally show ?case by simp
qed simp_all
  

definition Polygamma_approx :: "nat  nat  'a  'a :: {real_normed_field, ln}" where
  "Polygamma_approx j m = 
     (deriv ^^ j) (λx. (x - 1 / 2) * ln x - x + of_real (ln (2 * pi)) / 2 + stirling_sum 0 m x)"
  
lemma Polygamma_approx_Suc: "Polygamma_approx (Suc j) m = deriv (Polygamma_approx j m)"
  by (simp add: Polygamma_approx_def)  

lemma Polygamma_approx_0: 
  "Polygamma_approx 0 m x = (x - 1/2) * ln x - x + of_real (ln (2*pi)) / 2 + stirling_sum 0 m x"
  by (simp add: Polygamma_approx_def)
    
lemma Polygamma_approx_1_complex: 
  "Re x > 0  
     Polygamma_approx (Suc 0) m x = ln x - 1 / (2*x) + stirling_sum (Suc 0) m x"
  unfolding Polygamma_approx_Suc Polygamma_approx_0
  by (intro DERIV_imp_deriv) 
     (auto intro!: derivative_eq_intros elim!: nonpos_Reals_cases simp: field_simps)
     
lemma Polygamma_approx_1_real: 
  "x > (0 :: real)  
     Polygamma_approx (Suc 0) m x = ln x - 1 / (2*x) + stirling_sum (Suc 0) m x"
  unfolding Polygamma_approx_Suc Polygamma_approx_0
  by (intro DERIV_imp_deriv) 
     (auto intro!: derivative_eq_intros elim!: nonpos_Reals_cases simp: field_simps)
 
lemma stirling_sum_2_conv_stirling_sum'_1:
  fixes x :: "'a :: {real_div_algebra, field_char_0}"
  assumes "m > 0" "x  0"
  shows   "stirling_sum' 1 m x = 1 / x + 1 / (2 * x^2) + stirling_sum 2 m x"
proof -
  have pochhammer_2: "pochhammer (of_nat k) 2 = of_nat k * of_nat (Suc k)" for k 
    by (simp add: pochhammer_Suc eval_nat_numeral add_ac)
  have "stirling_sum 2 m x = 
          (k = Suc 0..<m. of_real (bernoulli' (Suc k)) * inverse x ^ Suc (Suc k))"
    unfolding stirling_sum_def pochhammer_2 power2_minus power_one mult_1_left
    by (intro sum.cong refl)
       (simp_all add: stirling_sum_def pochhammer_2 power2_eq_square divide_simps bernoulli'_def
                 del: of_nat_Suc power_Suc)
  also have "1 / (2 * x^2) +  = 
               (k=0..<m. of_real (bernoulli' (Suc k)) * inverse x ^ Suc (Suc k))" using assms
    by (subst (2) sum.atLeast_Suc_lessThan) (simp_all add: power2_eq_square field_simps)
  also have "1 / x +  = (k=0..<Suc m. of_real (bernoulli' k) * inverse x ^ Suc k)"
    by (subst sum.atLeast0_lessThan_Suc_shift) (simp_all add: bernoulli'_def divide_simps)
  also have " = (km. of_real (bernoulli' k) * inverse x ^ Suc k)"
    by (intro sum.cong) auto
  also have " = stirling_sum' 1 m x" by (simp add: stirling_sum'_def)
  finally show ?thesis by (simp add: add_ac)
qed

lemma Polygamma_approx_2_real: 
  assumes "x > (0::real)" "m > 0"
  shows   "Polygamma_approx (Suc (Suc 0)) m x = stirling_sum' 1 m x"
proof -
  have "Polygamma_approx (Suc (Suc 0)) m x = deriv (Polygamma_approx (Suc 0) m) x" 
    by (simp add: Polygamma_approx_Suc)
  also have " = deriv (λx. ln x - 1 / (2*x) + stirling_sum (Suc 0) m x) x"
    using eventually_nhds_in_open[of "{0<..}" x] assms
    by (intro deriv_cong_ev) (auto elim!: eventually_mono simp: Polygamma_approx_1_real)
  also have " = 1 / x + 1 / (2*x^2) + stirling_sum (Suc (Suc 0)) m x" using assms
    by (intro DERIV_imp_deriv) (auto intro!: derivative_eq_intros 
           elim!: nonpos_Reals_cases simp: field_simps power2_eq_square)
  also have " = stirling_sum' 1 m x" using stirling_sum_2_conv_stirling_sum'_1[of m x] assms
    by (simp add: eval_nat_numeral)
  finally show ?thesis .
qed
  
lemma Polygamma_approx_2_complex: 
  assumes "Re x > 0" "m > 0"
  shows   "Polygamma_approx (Suc (Suc 0)) m x = stirling_sum' 1 m x"
proof -
  have "Polygamma_approx (Suc (Suc 0)) m x = deriv (Polygamma_approx (Suc 0) m) x" 
    by (simp add: Polygamma_approx_Suc)
  also have " = deriv (λx. ln x - 1 / (2*x) + stirling_sum (Suc 0) m x) x"
    using eventually_nhds_in_open[of "{s. Re s > 0}" x] assms
    by (intro deriv_cong_ev)
       (auto simp: open_halfspace_Re_gt elim!: eventually_mono simp: Polygamma_approx_1_complex)
  also have " = 1 / x + 1 / (2*x^2) + stirling_sum (Suc (Suc 0)) m x" using assms
    by (intro DERIV_imp_deriv) (auto intro!: derivative_eq_intros 
           elim!: nonpos_Reals_cases simp: field_simps power2_eq_square)
  also have " = stirling_sum' 1 m x" using stirling_sum_2_conv_stirling_sum'_1[of m x] assms
      by (subst stirling_sum_2_conv_stirling_sum'_1) (auto simp: eval_nat_numeral)
  finally show ?thesis .
qed
  
lemma Polygamma_approx_ge_2_real: 
  assumes "x > (0::real)" "m > 0"
  shows   "Polygamma_approx (Suc (Suc j)) m x = stirling_sum' (Suc j) m x"
using assms(1)
proof (induction j arbitrary: x)
  case (0 x)
  with assms show ?case by (simp add: Polygamma_approx_2_real)
next
  case (Suc j x)
  have "Polygamma_approx (Suc (Suc (Suc j))) m x = deriv (Polygamma_approx (Suc (Suc j)) m) x"
    by (simp add: Polygamma_approx_Suc)
  also have " = deriv (stirling_sum' (Suc j) m) x"
    using eventually_nhds_in_open[of "{0<..}" x] Suc.prems
    by (intro deriv_cong_ev refl) (auto elim!: eventually_mono simp: Suc.IH)
  also have " = stirling_sum' (Suc (Suc j)) m x" using Suc.prems
    by (intro DERIV_imp_deriv derivative_intros) simp_all
  finally show ?case .
qed
  
lemma Polygamma_approx_ge_2_complex:
  assumes "Re x > 0" "m > 0"
  shows   "Polygamma_approx (Suc (Suc j)) m x = stirling_sum' (Suc j) m x"
using assms(1)
proof (induction j arbitrary: x)
  case (0 x)
  with assms show ?case by (simp add: Polygamma_approx_2_complex)
next
  case (Suc j x)
  have "Polygamma_approx (Suc (Suc (Suc j))) m x = deriv (Polygamma_approx (Suc (Suc j)) m) x"
    by (simp add: Polygamma_approx_Suc)
  also have " = deriv (stirling_sum' (Suc j) m) x"
    using eventually_nhds_in_open[of "{x. Re x > 0}" x] Suc.prems
    by (intro deriv_cong_ev refl) (auto elim!: eventually_mono simp: Suc.IH open_halfspace_Re_gt)
  also have " = stirling_sum' (Suc (Suc j)) m x" using Suc.prems
    by (intro DERIV_imp_deriv derivative_intros) simp_all
  finally show ?case .
qed

lemma Polygamma_approx_complex_of_real:
  assumes "x > 0" "m > 0"
  shows   "Polygamma_approx j m (complex_of_real x) = of_real (Polygamma_approx j m x)"
proof (cases j)
  case 0
  with assms show ?thesis by (simp add: Polygamma_approx_0 Ln_of_real stirling_sum_complex_of_real)
next
  case [simp]: (Suc j')
  thus ?thesis
  proof (cases j')
    case 0
    with assms show ?thesis 
      by (simp add: Polygamma_approx_1_complex 
                    Polygamma_approx_1_real stirling_sum_complex_of_real Ln_of_real)
  next
    case (Suc j'')
    with assms show ?thesis
      by (simp add: Polygamma_approx_ge_2_complex Polygamma_approx_ge_2_real 
                    stirling_sum'_complex_of_real)
  qed
qed
  
lemma higher_deriv_Polygamma_approx [simp]: 
  "(deriv ^^ j) (Polygamma_approx i m) = Polygamma_approx (j + i) m"
  by (simp add: Polygamma_approx_def funpow_add)
  
lemma stirling_sum_holomorphic [holomorphic_intros]:
  "0  A  stirling_sum j m holomorphic_on A"
  unfolding stirling_sum_def by (intro holomorphic_intros) auto
  
lemma Polygamma_approx_holomorphic [holomorphic_intros]:
  "Polygamma_approx j m holomorphic_on {s. Re s > 0}"
  unfolding Polygamma_approx_def
  by (intro holomorphic_intros) (auto simp: open_halfspace_Re_gt elim!: nonpos_Reals_cases)
  
lemma higher_deriv_lnGamma_stirling:
  assumes m: "m > 0"
  shows   "(λx::real. (deriv ^^ j) ln_Gamma x - Polygamma_approx j m x)  O(λx. 1 / x ^ (m + j))"
proof -
  have "eventually (λx. ¦(deriv ^^ j) ln_Gamma x - Polygamma_approx j m x¦ =
                          inverse (real m) * ¦(deriv ^^ j) (stirling_integral m) x¦) at_top"
    using eventually_gt_at_top[of "0::real"]
  proof eventually_elim
    case (elim x)
    note x = this
    have "F y in nhds (complex_of_real x). y  - 0"
      using elim by (intro eventually_nhds_in_open) auto
    hence "(deriv ^^ j) (λx. ln_Gamma x - Polygamma_approx 0 m x) (complex_of_real x) =
            (deriv ^^ j) (λx. (-inverse (of_nat m)) * stirling_integral m x) (complex_of_real x)"
      using x m
      by (intro higher_deriv_cong_ev refl)
         (auto elim!: eventually_mono simp: ln_Gamma_stirling_complex Polygamma_approx_def 
            field_simps open_halfspace_Re_gt stirling_sum_def)
    also have " = - inverse (of_nat m) * (deriv ^^ j) (stirling_integral m) (of_real x)" using x m
      by (intro higher_deriv_cmult[of _ "-0"] stirling_integral_holomorphic)
         (auto simp: open_halfspace_Re_gt)
    also have "(deriv ^^ j) (λx. ln_Gamma x - Polygamma_approx 0 m x) (complex_of_real x) = 
                 (deriv ^^ j) ln_Gamma (of_real x) - (deriv ^^ j) (Polygamma_approx 0 m) (of_real x)"
      using x 
      by (intro higher_deriv_diff[of _ "{s. Re s > 0}"])
         (auto intro!: holomorphic_intros elim!: nonpos_Reals_cases simp: open_halfspace_Re_gt)
    also have "(deriv ^^ j) (Polygamma_approx 0 m) (complex_of_real x) =
                 of_real (Polygamma_approx j m x)" using x m
      by (simp add: Polygamma_approx_complex_of_real)
    also have "norm (- inverse (of_nat m) * (deriv ^^ j) (stirling_integral m) (complex_of_real x)) = 
                 inverse (real m) * ¦(deriv ^^ j) (stirling_integral m) x¦"
      using x m by (simp add: norm_mult norm_inverse deriv_stirling_integral_complex_of_real)
    also have "(deriv ^^ j) ln_Gamma (complex_of_real x) = of_real ((deriv ^^ j) ln_Gamma x)" using x
      by (simp add: higher_deriv_ln_Gamma_complex_of_real)
    also have "norm ( - of_real (Polygamma_approx j m x)) = 
                 ¦(deriv ^^ j) ln_Gamma x - Polygamma_approx j m x¦"
      by (simp only: of_real_diff [symmetric] norm_of_real)
    finally show ?case .
  qed
  from bigthetaI_cong[OF this] m
    have "(λx::real. (deriv ^^ j) ln_Gamma x - Polygamma_approx j m x)  
             Θ(λx. (deriv ^^ j) (stirling_integral m) x)" by simp
  also have "(λx::real. (deriv ^^ j) (stirling_integral m) x)  O(λx. 1 / x ^ (m + j))" using m
      by (rule deriv_stirling_integral_real_bound)
  finally show ?thesis .
qed

lemma Polygamma_approx_1_real':
  assumes x: "(x::real) > 0" and m: "m > 0"
  shows   "Polygamma_approx 1 m x = ln x - (k = Suc 0..m. bernoulli' k * inverse x ^ k / real k)"
proof -
  have "Polygamma_approx 1 m x = ln x - (1 / (2 * x) + 
          (k=Suc 0..<m. bernoulli (Suc k) * inverse x ^ Suc k / real (Suc k)))"
    (is "_ = _ - (_ + ?S)") using x by (simp add: Polygamma_approx_1_real stirling_sum_def)
  also have "?S = (k=Suc 0..<m. bernoulli' (Suc k) * inverse x ^ Suc k / real (Suc k))"
    by (intro sum.cong refl) (simp_all add: bernoulli'_def)
  also have "1 / (2 * x) +  = 
               (k=0..<m. bernoulli' (Suc k) * inverse x ^ Suc k / real (Suc k))" using m
    by (subst (2) sum.atLeast_Suc_lessThan) (simp_all add: field_simps)
  also have " = (k = Suc 0..m. bernoulli' k * inverse x ^ k / real k)" using assms
    by (subst sum.shift_bounds_Suc_ivl [symmetric]) (simp add: atLeastLessThanSuc_atLeastAtMost)
  finally show ?thesis .
qed

theorem
  assumes m: "m > 0"
  shows   ln_Gamma_real_asymptotics:
            "(λx. ln_Gamma x - ((x - 1 / 2) * ln x - x + ln (2 * pi) / 2 +
                     (k = 1..<m. bernoulli (Suc k) / (real k * real (Suc k)) / x^k)))
                 O(λx. 1 / x ^ m)" (is ?th1)
    and   Digamma_real_asymptotics:
            "(λx. Digamma x - (ln x - (k=1..m. bernoulli' k / real k / x ^ k)))
                 O(λx. 1 / (x ^ Suc m))" (is ?th2)
    and   Polygamma_real_asymptotics: "j > 0  
             (λx. Polygamma j x - (- 1) ^ Suc j * (km. bernoulli' k *
                     pochhammer (real (Suc k)) (j - 1) / x ^ (k + j))) 
                 O(λx. 1 / x ^ (m+j+1))" (is "_  ?th3")
proof -
  define G :: "nat  real  real" where 
    "G = (λm. if m = 0 then ln_Gamma else Polygamma (m - 1))"
  have *: "(λx. G j x - h x)  O(λx. 1 / x ^ (m + j))"
    if "x::real. x > 0  Polygamma_approx j m x = h x" for j h
  proof -
    have "(λx. G j x - h x)  
            Θ(λx. (deriv ^^ j) ln_Gamma x - Polygamma_approx j m x)" (is "_  Θ(?f)")
      using that
      by (intro bigthetaI_cong) (auto intro: eventually_mono[OF eventually_gt_at_top[of "0::real"]]
            simp del: funpow.simps simp: higher_deriv_ln_Gamma_real G_def)
    also have "?f  O(λx::real. 1 / x ^ (m + j))" using m
      by (rule higher_deriv_lnGamma_stirling)
    finally show ?thesis .
  qed  
    
  note [[simproc del: simplify_landau_sum]]
  from *[OF Polygamma_approx_0] assms show ?th1 
    by (simp add: G_def Polygamma_approx_0 stirling_sum_def field_simps)
  from *[OF Polygamma_approx_1_real'] assms show ?th2 by (simp add: G_def field_simps)
      
  assume j: "j > 0"
  from *[OF Polygamma_approx_ge_2_real, of "j - 1"] assms j show ?th3
    by (simp add: G_def stirling_sum'_def power_add power_diff field_simps)
qed




subsection ‹Asymptotics of the complex Gamma function›

text ‹
  The m›-th order remainder of Stirling's formula for $\log\Gamma$ is $O(s^{-m})$ uniformly over
  any complex cone $\text{arg}(z) \leq \alpha$, $z\neq 0$ for any angle
  $\alpha\in(0, \pi)$. This means that there is bounded by $c z^{-m}$ for some constant $c$ for
  all $z$ in this cone.
›
context
  fixes F and α
  assumes α: "α  {0<..<pi}"
  defines "F  principal (complex_cone' α - {0})"
begin

lemma stirling_integral_bigo:
  fixes m :: nat
  assumes m: "m > 0"
  shows   "stirling_integral m  O[F](λs. 1 / s ^ m)"
proof -
  obtain c where c: "s. s  complex_cone' α - {0}  norm (stirling_integral m s)  c / norm s ^ m"
    using stirling_integral_bound'[OF m > 0 α] by blast
  have "0  norm (stirling_integral m 1 :: complex)"
    by simp
  also have "  c"
    using c[of 1] α by simp
  finally have "c  0" .

  have "eventually (λs. s  complex_cone' α - {0}) F"
    unfolding F_def by (auto simp: eventually_principal)
  hence "eventually (λs. norm (stirling_integral m s) 
                     c * norm (1 / s ^ m)) F"
    by eventually_elim (use c in simp add: norm_divide norm_power›)
  thus "stirling_integral m  O[F](λs. 1 / s ^ m)"
    by (intro bigoI[of _ c]) auto
qed

end

text ‹
  The following is a more explicit statement of this:
›
theorem ln_Gamma_complex_asymptotics_explicit:
  fixes m :: nat and α :: real
  assumes "m > 0" and "α  {0<..<pi}"
  obtains C :: real and R :: "complex  complex"
  where "s::complex. s  0 
               ln_Gamma s = (s - 1/2) * ln s - s + ln (2 * pi) / 2 +
                            (k=1..<m. bernoulli (k+1) / (k * (k+1) * s ^ k)) - R s"
    and "s. s  0  ¦arg s¦  α  norm (R s)  C / norm s ^ m"    
proof -
  obtain c where c: "s. s  complex_cone' α - {0}  norm (stirling_integral m s)  c / norm s ^ m"
    using stirling_integral_bound'[OF assms] by blast
  have "0  norm (stirling_integral m 1 :: complex)"
    by simp
  also have "  c"
    using c[of 1] assms by simp
  finally have "c  0" .
  define R where "R = (λs::complex. stirling_integral m s / of_nat m)"
  show ?thesis
  proof (rule that)
    from ln_Gamma_stirling_complex[of _ m] assms show
           "s::complex. s  0 
               ln_Gamma s = (s - 1 / 2) * ln