# Theory Gamma_Function

```(*  Title:    HOL/Analysis/Gamma_Function.thy
Author:   Manuel Eberl, TU München
*)

section ‹The Gamma Function›

theory Gamma_Function
imports
Equivalence_Lebesgue_Henstock_Integration
Summation_Tests
Harmonic_Numbers
"HOL-Library.Nonpos_Ints"
"HOL-Library.Periodic_Fun"
begin

text ‹
Several equivalent definitions of the Gamma function and its
most important properties. Also contains the definition and some properties
of the log-Gamma function and the Digamma function and the other Polygamma functions.

Based on the Gamma function, we also prove the Weierstra{\ss} product form of the
sin function and, based on this, the solution of the Basel problem (the
sum over all \<^term>‹1 / (n::nat)^2›.
›

lemma pochhammer_eq_0_imp_nonpos_Int:
"pochhammer (x::'a::field_char_0) n = 0 ⟹ x ∈ ℤ⇩≤⇩0"
by (auto simp: pochhammer_eq_0_iff)

lemma closed_nonpos_Ints [simp]: "closed (ℤ⇩≤⇩0 :: 'a :: real_normed_algebra_1 set)"
proof -
have "ℤ⇩≤⇩0 = (of_int ` {n. n ≤ 0} :: 'a set)"
by (auto elim!: nonpos_Ints_cases intro!: nonpos_Ints_of_int)
also have "closed …" by (rule closed_of_int_image)
finally show ?thesis .
qed

lemma plus_one_in_nonpos_Ints_imp: "z + 1 ∈ ℤ⇩≤⇩0 ⟹ z ∈ ℤ⇩≤⇩0"
using nonpos_Ints_diff_Nats[of "z+1" "1"] by simp_all

lemma of_int_in_nonpos_Ints_iff:
"(of_int n :: 'a :: ring_char_0) ∈ ℤ⇩≤⇩0 ⟷ n ≤ 0"
by (auto simp: nonpos_Ints_def)

lemma one_plus_of_int_in_nonpos_Ints_iff:
"(1 + of_int n :: 'a :: ring_char_0) ∈ ℤ⇩≤⇩0 ⟷ n ≤ -1"
proof -
have "1 + of_int n = (of_int (n + 1) :: 'a)" by simp
also have "… ∈ ℤ⇩≤⇩0 ⟷ n + 1 ≤ 0" by (subst of_int_in_nonpos_Ints_iff) simp_all
also have "… ⟷ n ≤ -1" by presburger
finally show ?thesis .
qed

lemma one_minus_of_nat_in_nonpos_Ints_iff:
"(1 - of_nat n :: 'a :: ring_char_0) ∈ ℤ⇩≤⇩0 ⟷ n > 0"
proof -
have "(1 - of_nat n :: 'a) = of_int (1 - int n)" by simp
also have "… ∈ ℤ⇩≤⇩0 ⟷ n > 0" by (subst of_int_in_nonpos_Ints_iff) presburger
finally show ?thesis .
qed

lemma fraction_not_in_ints:
assumes "¬(n dvd m)" "n ≠ 0"
shows   "of_int m / of_int n ∉ (ℤ :: 'a :: {division_ring,ring_char_0} set)"
proof
assume "of_int m / (of_int n :: 'a) ∈ ℤ"
then obtain k where "of_int m / of_int n = (of_int k :: 'a)" by (elim Ints_cases)
with assms have "of_int m = (of_int (k * n) :: 'a)" by (auto simp add: field_split_simps)
hence "m = k * n" by (subst (asm) of_int_eq_iff)
hence "n dvd m" by simp
with assms(1) show False by contradiction
qed

lemma fraction_not_in_nats:
assumes "¬n dvd m" "n ≠ 0"
shows   "of_int m / of_int n ∉ (ℕ :: 'a :: {division_ring,ring_char_0} set)"
proof
assume "of_int m / of_int n ∈ (ℕ :: 'a set)"
also note Nats_subset_Ints
finally have "of_int m / of_int n ∈ (ℤ :: 'a set)" .
moreover have "of_int m / of_int n ∉ (ℤ :: 'a set)"
using assms by (intro fraction_not_in_ints)
qed

lemma not_in_Ints_imp_not_in_nonpos_Ints: "z ∉ ℤ ⟹ z ∉ ℤ⇩≤⇩0"
by (auto simp: Ints_def nonpos_Ints_def)

lemma double_in_nonpos_Ints_imp:
assumes "2 * (z :: 'a :: field_char_0) ∈ ℤ⇩≤⇩0"
shows   "z ∈ ℤ⇩≤⇩0 ∨ z + 1/2 ∈ ℤ⇩≤⇩0"
proof-
from assms obtain k where k: "2 * z = - of_nat k" by (elim nonpos_Ints_cases')
thus ?thesis by (cases "even k") (auto elim!: evenE oddE simp: field_simps)
qed

lemma sin_series: "(λn. ((-1)^n / fact (2*n+1)) *⇩R z^(2*n+1)) sums sin z"
proof -
from sin_converges[of z] have "(λn. sin_coeff n *⇩R z^n) sums sin z" .
also have "(λn. sin_coeff n *⇩R z^n) sums sin z ⟷
(λn. ((-1)^n / fact (2*n+1)) *⇩R z^(2*n+1)) sums sin z"
by (subst sums_mono_reindex[of "λn. 2*n+1", symmetric])
(auto simp: sin_coeff_def strict_mono_def ac_simps elim!: oddE)
finally show ?thesis .
qed

lemma cos_series: "(λn. ((-1)^n / fact (2*n)) *⇩R z^(2*n)) sums cos z"
proof -
from cos_converges[of z] have "(λn. cos_coeff n *⇩R z^n) sums cos z" .
also have "(λn. cos_coeff n *⇩R z^n) sums cos z ⟷
(λn. ((-1)^n / fact (2*n)) *⇩R z^(2*n)) sums cos z"
by (subst sums_mono_reindex[of "λn. 2*n", symmetric])
(auto simp: cos_coeff_def strict_mono_def ac_simps elim!: evenE)
finally show ?thesis .
qed

lemma sin_z_over_z_series:
fixes z :: "'a :: {real_normed_field,banach}"
assumes "z ≠ 0"
shows   "(λn. (-1)^n / fact (2*n+1) * z^(2*n)) sums (sin z / z)"
proof -
from sin_series[of z] have "(λn. z * ((-1)^n / fact (2*n+1)) * z^(2*n)) sums sin z"
from sums_mult[OF this, of "inverse z"] and assms show ?thesis
qed

lemma sin_z_over_z_series':
fixes z :: "'a :: {real_normed_field,banach}"
assumes "z ≠ 0"
shows   "(λn. sin_coeff (n+1) *⇩R z^n) sums (sin z / z)"
proof -
from sums_split_initial_segment[OF sin_converges[of z], of 1]
have "(λn. z * (sin_coeff (n+1) *⇩R z ^ n)) sums sin z" by simp
from sums_mult[OF this, of "inverse z"] assms show ?thesis by (simp add: field_simps)
qed

lemma has_field_derivative_sin_z_over_z:
fixes A :: "'a :: {real_normed_field,banach} set"
shows "((λz. if z = 0 then 1 else sin z / z) has_field_derivative 0) (at 0 within A)"
(is "(?f has_field_derivative ?f') _")
proof (rule has_field_derivative_at_within)
have "((λz::'a. ∑n. of_real (sin_coeff (n+1)) * z^n)
has_field_derivative (∑n. diffs (λn. of_real (sin_coeff (n+1))) n * 0^n)) (at 0)"
proof (rule termdiffs_strong)
from summable_ignore_initial_segment[OF sums_summable[OF sin_converges[of "1::'a"]], of 1]
show "summable (λn. of_real (sin_coeff (n+1)) * (1::'a)^n)" by (simp add: of_real_def)
qed simp
also have "(λz::'a. ∑n. of_real (sin_coeff (n+1)) * z^n) = ?f"
proof
fix z
show "(∑n. of_real (sin_coeff (n+1)) * z^n)  = ?f z"
by (cases "z = 0") (insert sin_z_over_z_series'[of z],
qed
also have "(∑n. diffs (λn. of_real (sin_coeff (n + 1))) n * (0::'a) ^ n) =
diffs (λn. of_real (sin_coeff (Suc n))) 0" by simp
also have "… = 0" by (simp add: sin_coeff_def diffs_def)
finally show "((λz::'a. if z = 0 then 1 else sin z / z) has_field_derivative 0) (at 0)" .
qed

lemma round_Re_minimises_norm:
"norm ((z::complex) - of_int m) ≥ norm (z - of_int (round (Re z)))"
proof -
let ?n = "round (Re z)"
have "norm (z - of_int ?n) = sqrt ((Re z - of_int ?n)⇧2 + (Im z)⇧2)"
also have "¦Re z - of_int ?n¦ ≤ ¦Re z - of_int m¦" by (rule round_diff_minimal)
hence "sqrt ((Re z - of_int ?n)⇧2 + (Im z)⇧2) ≤ sqrt ((Re z - of_int m)⇧2 + (Im z)⇧2)"
also have "… = norm (z - of_int m)" by (simp add: cmod_def)
finally show ?thesis .
qed

lemma Re_pos_in_ball:
assumes "Re z > 0" "t ∈ ball z (Re z/2)"
shows   "Re t > 0"
proof -
have "Re (z - t) ≤ norm (z - t)" by (rule complex_Re_le_cmod)
also from assms have "… < Re z / 2" by (simp add: dist_complex_def)
finally show "Re t > 0" using assms by simp
qed

lemma no_nonpos_Int_in_ball_complex:
assumes "Re z > 0" "t ∈ ball z (Re z/2)"
shows   "t ∉ ℤ⇩≤⇩0"
using Re_pos_in_ball[OF assms] by (force elim!: nonpos_Ints_cases)

lemma no_nonpos_Int_in_ball:
assumes "t ∈ ball z (dist z (round (Re z)))"
shows   "t ∉ ℤ⇩≤⇩0"
proof
assume "t ∈ ℤ⇩≤⇩0"
then obtain n where "t = of_int n" by (auto elim!: nonpos_Ints_cases)
have "dist z (of_int n) ≤ dist z t + dist t (of_int n)" by (rule dist_triangle)
also from assms have "dist z t < dist z (round (Re z))" by simp
also have "… ≤ dist z (of_int n)"
using round_Re_minimises_norm[of z] by (simp add: dist_complex_def)
finally have "dist t (of_int n) > 0" by simp
with ‹t = of_int n› show False by simp
qed

lemma no_nonpos_Int_in_ball':
assumes "(z :: 'a :: {euclidean_space,real_normed_algebra_1}) ∉ ℤ⇩≤⇩0"
obtains d where "d > 0" "⋀t. t ∈ ball z d ⟹ t ∉ ℤ⇩≤⇩0"
proof (rule that)
from assms show "setdist {z} ℤ⇩≤⇩0 > 0" by (subst setdist_gt_0_compact_closed) auto
next
fix t assume "t ∈ ball z (setdist {z} ℤ⇩≤⇩0)"
thus "t ∉ ℤ⇩≤⇩0" using setdist_le_dist[of z "{z}" t "ℤ⇩≤⇩0"] by force
qed

lemma no_nonpos_Real_in_ball:
assumes z: "z ∉ ℝ⇩≤⇩0" and t: "t ∈ ball z (if Im z = 0 then Re z / 2 else abs (Im z) / 2)"
shows   "t ∉ ℝ⇩≤⇩0"
using z
proof (cases "Im z = 0")
assume A: "Im z = 0"
with z have "Re z > 0" by (force simp add: complex_nonpos_Reals_iff)
with t A Re_pos_in_ball[of z t] show ?thesis by (force simp add: complex_nonpos_Reals_iff)
next
assume A: "Im z ≠ 0"
have "abs (Im z) - abs (Im t) ≤ abs (Im z - Im t)" by linarith
also have "… = abs (Im (z - t))" by simp
also have "… ≤ norm (z - t)" by (rule abs_Im_le_cmod)
also from A t have "… ≤ abs (Im z) / 2" by (simp add: dist_complex_def)
finally have "abs (Im t) > 0" using A by simp
thus ?thesis by (force simp add: complex_nonpos_Reals_iff)
qed

subsection ‹The Euler form and the logarithmic Gamma function›

text ‹
We define the Gamma function by first defining its multiplicative inverse ‹rGamma›.
This is more convenient because ‹rGamma› is entire, which makes proofs of its
properties more convenient because one does not have to watch out for discontinuities.
(e.g. ‹rGamma› fulfils ‹rGamma z = z * rGamma (z + 1)› everywhere, whereas the ‹Γ› function
does not fulfil the analogous equation on the non-positive integers)

We define the ‹Γ› function (resp.\ its reciprocale) in the Euler form. This form has the advantage
that it is a relatively simple limit that converges everywhere. The limit at the poles is 0
(due to division by 0). The functional equation ‹Gamma (z + 1) = z * Gamma z› follows
immediately from the definition.
›

definition✐‹tag important› Gamma_series :: "('a :: {banach,real_normed_field}) ⇒ nat ⇒ 'a" where
"Gamma_series z n = fact n * exp (z * of_real (ln (of_nat n))) / pochhammer z (n+1)"

definition Gamma_series' :: "('a :: {banach,real_normed_field}) ⇒ nat ⇒ 'a" where
"Gamma_series' z n = fact (n - 1) * exp (z * of_real (ln (of_nat n))) / pochhammer z n"

definition rGamma_series :: "('a :: {banach,real_normed_field}) ⇒ nat ⇒ 'a" where
"rGamma_series z n = pochhammer z (n+1) / (fact n * exp (z * of_real (ln (of_nat n))))"

lemma Gamma_series_altdef: "Gamma_series z n = inverse (rGamma_series z n)"
and rGamma_series_altdef: "rGamma_series z n = inverse (Gamma_series z n)"
unfolding Gamma_series_def rGamma_series_def by simp_all

lemma rGamma_series_minus_of_nat:
"eventually (λn. rGamma_series (- of_nat k) n = 0) sequentially"
using eventually_ge_at_top[of k]
by eventually_elim (auto simp: rGamma_series_def pochhammer_of_nat_eq_0_iff)

lemma Gamma_series_minus_of_nat:
"eventually (λn. Gamma_series (- of_nat k) n = 0) sequentially"
using eventually_ge_at_top[of k]
by eventually_elim (auto simp: Gamma_series_def pochhammer_of_nat_eq_0_iff)

lemma Gamma_series'_minus_of_nat:
"eventually (λn. Gamma_series' (- of_nat k) n = 0) sequentially"
using eventually_gt_at_top[of k]
by eventually_elim (auto simp: Gamma_series'_def pochhammer_of_nat_eq_0_iff)

lemma rGamma_series_nonpos_Ints_LIMSEQ: "z ∈ ℤ⇩≤⇩0 ⟹ rGamma_series z ⇢ 0"
by (elim nonpos_Ints_cases', hypsubst, subst tendsto_cong, rule rGamma_series_minus_of_nat, simp)

lemma Gamma_series_nonpos_Ints_LIMSEQ: "z ∈ ℤ⇩≤⇩0 ⟹ Gamma_series z ⇢ 0"
by (elim nonpos_Ints_cases', hypsubst, subst tendsto_cong, rule Gamma_series_minus_of_nat, simp)

lemma Gamma_series'_nonpos_Ints_LIMSEQ: "z ∈ ℤ⇩≤⇩0 ⟹ Gamma_series' z ⇢ 0"
by (elim nonpos_Ints_cases', hypsubst, subst tendsto_cong, rule Gamma_series'_minus_of_nat, simp)

lemma Gamma_series_Gamma_series':
assumes z: "z ∉ ℤ⇩≤⇩0"
shows   "(λn. Gamma_series' z n / Gamma_series z n) ⇢ 1"
proof (rule Lim_transform_eventually)
from eventually_gt_at_top[of "0::nat"]
show "eventually (λn. z / of_nat n + 1 = Gamma_series' z n / Gamma_series z n) sequentially"
proof eventually_elim
fix n :: nat assume n: "n > 0"
from n z have "Gamma_series' z n / Gamma_series z n = (z + of_nat n) / of_nat n"
by (cases n, simp)
(auto simp add: Gamma_series_def Gamma_series'_def pochhammer_rec'
dest: pochhammer_eq_0_imp_nonpos_Int plus_of_nat_eq_0_imp)
also from n have "… = z / of_nat n + 1" by (simp add: field_split_simps)
finally show "z / of_nat n + 1 = Gamma_series' z n / Gamma_series z n" ..
qed
have "(λx. z / of_nat x) ⇢ 0"
by (rule tendsto_norm_zero_cancel)
(insert tendsto_mult[OF tendsto_const[of "norm z"] lim_inverse_n],
from tendsto_add[OF this tendsto_const[of 1]] show "(λn. z / of_nat n + 1) ⇢ 1" by simp
qed

text ‹
We now show that the series that defines the ‹Γ› function in the Euler form converges
and that the function defined by it is continuous on the complex halfspace with positive
real part.

We do this by showing that the logarithm of the Euler series is continuous and converges
locally uniformly, which means that the log-Gamma function defined by its limit is also
continuous.

This will later allow us to lift holomorphicity and continuity from the log-Gamma
function to the inverse of the Gamma function, and from that to the Gamma function itself.
›

definition✐‹tag important› ln_Gamma_series :: "('a :: {banach,real_normed_field,ln}) ⇒ nat ⇒ 'a" where
"ln_Gamma_series z n = z * ln (of_nat n) - ln z - (∑k=1..n. ln (z / of_nat k + 1))"

definition✐‹tag unimportant› ln_Gamma_series' :: "('a :: {banach,real_normed_field,ln}) ⇒ nat ⇒ 'a" where
"ln_Gamma_series' z n =
- euler_mascheroni*z - ln z + (∑k=1..n. z / of_nat n - ln (z / of_nat k + 1))"

definition ln_Gamma :: "('a :: {banach,real_normed_field,ln}) ⇒ 'a" where
"ln_Gamma z = lim (ln_Gamma_series z)"

text ‹
We now show that the log-Gamma series converges locally uniformly for all complex numbers except
the non-positive integers. We do this by proving that the series is locally Cauchy.
›

context
begin

private lemma ln_Gamma_series_complex_converges_aux:
fixes z :: complex and k :: nat
assumes z: "z ≠ 0" and k: "of_nat k ≥ 2*norm z" "k ≥ 2"
shows "norm (z * ln (1 - 1/of_nat k) + ln (z/of_nat k + 1)) ≤ 2*(norm z + norm z^2) / of_nat k^2"
proof -
let ?k = "of_nat k :: complex" and ?z = "norm z"
have "z *ln (1 - 1/?k) + ln (z/?k+1) = z*(ln (1 - 1/?k :: complex) + 1/?k) + (ln (1+z/?k) - z/?k)"
also have "norm ... ≤ ?z * norm (ln (1-1/?k) + 1/?k :: complex) + norm (ln (1+z/?k) - z/?k)"
by (subst norm_mult [symmetric], rule norm_triangle_ineq)
also have "norm (Ln (1 + -1/?k) - (-1/?k)) ≤ (norm (-1/?k))⇧2 / (1 - norm(-1/?k))"
using k by (intro Ln_approx_linear) (simp add: norm_divide)
hence "?z * norm (ln (1-1/?k) + 1/?k) ≤ ?z * ((norm (1/?k))^2 / (1 - norm (1/?k)))"
by (intro mult_left_mono) simp_all
also have "... ≤ (?z * (of_nat k / (of_nat k - 1))) / of_nat k^2" using k
by (simp add: field_simps power2_eq_square norm_divide)
also have "... ≤ (?z * 2) / of_nat k^2" using k
by (intro divide_right_mono mult_left_mono) (simp_all add: field_simps)
also have "norm (ln (1+z/?k) - z/?k) ≤ norm (z/?k)^2 / (1 - norm (z/?k))" using k
by (intro Ln_approx_linear) (simp add: norm_divide)
hence "norm (ln (1+z/?k) - z/?k) ≤ ?z^2 / of_nat k^2 / (1 - ?z / of_nat k)"
also have "... ≤ (?z^2 * (of_nat k / (of_nat k - ?z))) / of_nat k^2" using k
also have "... ≤ (?z^2 * 2) / of_nat k^2" using k
by (intro divide_right_mono mult_left_mono) (simp_all add: field_simps)
finally show ?thesis by (simp only: distrib_left mult.commute)
qed

lemma ln_Gamma_series_complex_converges:
assumes z: "z ∉ ℤ⇩≤⇩0"
assumes d: "d > 0" "⋀n. n ∈ ℤ⇩≤⇩0 ⟹ norm (z - of_int n) > d"
shows "uniformly_convergent_on (ball z d) (λn z. ln_Gamma_series z n :: complex)"
proof (intro Cauchy_uniformly_convergent uniformly_Cauchy_onI')
fix e :: real assume e: "e > 0"
define e'' where "e'' = (SUP t∈ball z d. norm t + norm t^2)"
define e' where "e' = e / (2*e'')"
have "bounded ((λt. norm t + norm t^2) ` cball z d)"
by (intro compact_imp_bounded compact_continuous_image) (auto intro!: continuous_intros)
hence "bounded ((λt. norm t + norm t^2) ` ball z d)" by (rule bounded_subset) auto
hence bdd: "bdd_above ((λt. norm t + norm t^2) ` ball z d)" by (rule bounded_imp_bdd_above)

with z d(1) d(2)[of "-1"] have e''_pos: "e'' > 0" unfolding e''_def
by (subst less_cSUP_iff) (auto intro!: add_pos_nonneg bexI[of _ z])
have e'': "norm t + norm t^2 ≤ e''" if "t ∈ ball z d" for t unfolding e''_def using that
by (rule cSUP_upper[OF _ bdd])
from e z e''_pos have e': "e' > 0" unfolding e'_def

have "summable (λk. inverse ((real_of_nat k)^2))"
by (rule inverse_power_summable) simp
from summable_partial_sum_bound[OF this e']
obtain M where M: "⋀m n. M ≤ m ⟹ norm (∑k = m..n. inverse ((real k)⇧2)) < e'"
by auto

define N where "N = max 2 (max (nat ⌈2 * (norm z + d)⌉) M)"
{
from d have "⌈2 * (cmod z + d)⌉ ≥ ⌈0::real⌉"
by (intro ceiling_mono mult_nonneg_nonneg add_nonneg_nonneg) simp_all
hence "2 * (norm z + d) ≤ of_nat (nat ⌈2 * (norm z + d)⌉)" unfolding N_def
by (simp_all)
also have "... ≤ of_nat N" unfolding N_def
by (subst of_nat_le_iff) (rule max.coboundedI2, rule max.cobounded1)
finally have "of_nat N ≥ 2 * (norm z + d)" .
moreover have "N ≥ 2" "N ≥ M" unfolding N_def by simp_all
moreover have "(∑k=m..n. 1/(of_nat k)⇧2) < e'" if "m ≥ N" for m n
using M[OF order.trans[OF ‹N ≥ M› that]] unfolding real_norm_def
by (subst (asm) abs_of_nonneg) (auto intro: sum_nonneg simp: field_split_simps)
moreover note calculation
} note N = this

show "∃M. ∀t∈ball z d. ∀m≥M. ∀n>m. dist (ln_Gamma_series t m) (ln_Gamma_series t n) < e"
unfolding dist_complex_def
proof (intro exI[of _ N] ballI allI impI)
fix t m n assume t: "t ∈ ball z d" and mn: "m ≥ N" "n > m"
from d(2)[of 0] t have "0 < dist z 0 - dist z t" by (simp add: field_simps dist_complex_def)
also have "dist z 0 - dist z t ≤ dist 0 t" using dist_triangle[of 0 z t]
finally have t_nz: "t ≠ 0" by auto

have "norm t ≤ norm z + norm (t - z)" by (rule norm_triangle_sub)
also from t have "norm (t - z) < d" by (simp add: dist_complex_def norm_minus_commute)
also have "2 * (norm z + d) ≤ of_nat N" by (rule N)
also have "N ≤ m" by (rule mn)
finally have norm_t: "2 * norm t < of_nat m" by simp

have "ln_Gamma_series t m - ln_Gamma_series t n =
(-(t * Ln (of_nat n)) - (-(t * Ln (of_nat m)))) +
((∑k=1..n. Ln (t / of_nat k + 1)) - (∑k=1..m. Ln (t / of_nat k + 1)))"
also have "(∑k=1..n. Ln (t / of_nat k + 1)) - (∑k=1..m. Ln (t / of_nat k + 1)) =
(∑k∈{1..n}-{1..m}. Ln (t / of_nat k + 1))" using mn
also from mn have "{1..n}-{1..m} = {Suc m..n}" by fastforce
also have "-(t * Ln (of_nat n)) - (-(t * Ln (of_nat m))) =
(∑k = Suc m..n. t * Ln (of_nat (k - 1)) - t * Ln (of_nat k))" using mn
by (subst sum_telescope'' [symmetric]) simp_all
also have "... = (∑k = Suc m..n. t * Ln (of_nat (k - 1) / of_nat k))" using mn N
by (intro sum_cong_Suc)
(simp_all del: of_nat_Suc add: field_simps Ln_of_nat Ln_of_nat_over_of_nat)
also have "of_nat (k - 1) / of_nat k = 1 - 1 / (of_nat k :: complex)" if "k ∈ {Suc m..n}" for k
using that of_nat_eq_0_iff[of "Suc i" for i] by (cases k) (simp_all add: field_split_simps)
hence "(∑k = Suc m..n. t * Ln (of_nat (k - 1) / of_nat k)) =
(∑k = Suc m..n. t * Ln (1 - 1 / of_nat k))" using mn N
by (intro sum.cong) simp_all
also note sum.distrib [symmetric]
also have "norm (∑k=Suc m..n. t * Ln (1 - 1/of_nat k) + Ln (t/of_nat k + 1)) ≤
(∑k=Suc m..n. 2 * (norm t + (norm t)⇧2) / (real_of_nat k)⇧2)" using t_nz N(2) mn norm_t
by (intro order.trans[OF norm_sum sum_mono[OF ln_Gamma_series_complex_converges_aux]]) simp_all
also have "... ≤ 2 * (norm t + norm t^2) * (∑k=Suc m..n. 1 / (of_nat k)⇧2)"
also have "... < 2 * (norm t + norm t^2) * e'" using mn z t_nz
by (intro mult_strict_left_mono N mult_pos_pos add_pos_pos) simp_all
also from e''_pos have "... = e * ((cmod t + (cmod t)⇧2) / e'')"
by (simp add: e'_def field_simps power2_eq_square)
also from e''[OF t] e''_pos e
have "… ≤ e * 1" by (intro mult_left_mono) (simp_all add: field_simps)
finally show "norm (ln_Gamma_series t m - ln_Gamma_series t n) < e" by simp
qed
qed

end

lemma ln_Gamma_series_complex_converges':
assumes z: "(z :: complex) ∉ ℤ⇩≤⇩0"
shows "∃d>0. uniformly_convergent_on (ball z d) (λn z. ln_Gamma_series z n)"
proof -
define d' where "d' = Re z"
define d where "d = (if d' > 0 then d' / 2 else norm (z - of_int (round d')) / 2)"
have "of_int (round d') ∈ ℤ⇩≤⇩0" if "d' ≤ 0" using that
by (intro nonpos_Ints_of_int) (simp_all add: round_def)
with assms have d_pos: "d > 0" unfolding d_def by (force simp: not_less)

have "d < cmod (z - of_int n)" if "n ∈ ℤ⇩≤⇩0" for n
proof (cases "Re z > 0")
case True
from nonpos_Ints_nonpos[OF that] have n: "n ≤ 0" by simp
from True have "d = Re z/2" by (simp add: d_def d'_def)
also from n True have "… < Re (z - of_int n)" by simp
also have "… ≤ norm (z - of_int n)" by (rule complex_Re_le_cmod)
finally show ?thesis .
next
case False
with assms nonpos_Ints_of_int[of "round (Re z)"]
have "z ≠ of_int (round d')" by (auto simp: not_less)
with False have "d < norm (z - of_int (round d'))" by (simp add: d_def d'_def)
also have "… ≤ norm (z - of_int n)" unfolding d'_def by (rule round_Re_minimises_norm)
finally show ?thesis .
qed
hence conv: "uniformly_convergent_on (ball z d) (λn z. ln_Gamma_series z n)"
by (intro ln_Gamma_series_complex_converges d_pos z) simp_all
from d_pos conv show ?thesis by blast
qed

lemma ln_Gamma_series_complex_converges'': "(z :: complex) ∉ ℤ⇩≤⇩0 ⟹ convergent (ln_Gamma_series z)"
by (drule ln_Gamma_series_complex_converges') (auto intro: uniformly_convergent_imp_convergent)

theorem ln_Gamma_complex_LIMSEQ: "(z :: complex) ∉ ℤ⇩≤⇩0 ⟹ ln_Gamma_series z ⇢ ln_Gamma z"
using ln_Gamma_series_complex_converges'' by (simp add: convergent_LIMSEQ_iff ln_Gamma_def)

lemma exp_ln_Gamma_series_complex:
assumes "n > 0" "z ∉ ℤ⇩≤⇩0"
shows   "exp (ln_Gamma_series z n :: complex) = Gamma_series z n"
proof -
from assms obtain m where m: "n = Suc m" by (cases n) blast
from assms have "z ≠ 0" by (intro notI) auto
with assms have "exp (ln_Gamma_series z n) =
(of_nat n) powr z / (z * (∏k=1..n. exp (Ln (z / of_nat k + 1))))"
unfolding ln_Gamma_series_def powr_def by (simp add: exp_diff exp_sum)
also from assms have "(∏k=1..n. exp (Ln (z / of_nat k + 1))) = (∏k=1..n. z / of_nat k + 1)"
by (intro prod.cong[OF refl], subst exp_Ln) (auto simp: field_simps plus_of_nat_eq_0_imp)
also have "... = (∏k=1..n. z + k) / fact n"
(subst prod_dividef [symmetric], simp_all add: field_simps)
also from m have "z * ... = (∏k=0..n. z + k) / fact n"
by (simp add: prod.atLeast0_atMost_Suc_shift prod.atLeast_Suc_atMost_Suc_shift del: prod.cl_ivl_Suc)
also have "(∏k=0..n. z + k) = pochhammer z (Suc n)"
unfolding pochhammer_prod
also have "of_nat n powr z / (pochhammer z (Suc n) / fact n) = Gamma_series z n"
unfolding Gamma_series_def using assms by (simp add: field_split_simps powr_def)
finally show ?thesis .
qed

lemma ln_Gamma_series'_aux:
assumes "(z::complex) ∉ ℤ⇩≤⇩0"
shows   "(λk. z / of_nat (Suc k) - ln (1 + z / of_nat (Suc k))) sums
(ln_Gamma z + euler_mascheroni * z + ln z)" (is "?f sums ?s")
unfolding sums_def
proof (rule Lim_transform)
show "(λn. ln_Gamma_series z n + of_real (harm n - ln (of_nat n)) * z + ln z) ⇢ ?s"
(is "?g ⇢ _")
by (intro tendsto_intros ln_Gamma_complex_LIMSEQ euler_mascheroni_LIMSEQ_of_real assms)

have A: "eventually (λn. (∑k<n. ?f k) - ?g n = 0) sequentially"
using eventually_gt_at_top[of "0::nat"]
proof eventually_elim
fix n :: nat assume n: "n > 0"
have "(∑k<n. ?f k) = (∑k=1..n. z / of_nat k - ln (1 + z / of_nat k))"
by (subst atLeast0LessThan [symmetric], subst sum.shift_bounds_Suc_ivl [symmetric],
subst atLeastLessThanSuc_atLeastAtMost) simp_all
also have "… = z * of_real (harm n) - (∑k=1..n. ln (1 + z / of_nat k))"
by (simp add: harm_def sum_subtractf sum_distrib_left divide_inverse)
also from n have "… - ?g n = 0"
by (simp add: ln_Gamma_series_def sum_subtractf algebra_simps)
finally show "(∑k<n. ?f k) - ?g n = 0" .
qed
show "(λn. (∑k<n. ?f k) - ?g n) ⇢ 0" by (subst tendsto_cong[OF A]) simp_all
qed

lemma uniformly_summable_deriv_ln_Gamma:
assumes z: "(z :: 'a :: {real_normed_field,banach}) ≠ 0" and d: "d > 0" "d ≤ norm z/2"
shows "uniformly_convergent_on (ball z d)
(λk z. ∑i<k. inverse (of_nat (Suc i)) - inverse (z + of_nat (Suc i)))"
(is "uniformly_convergent_on _ (λk z. ∑i<k. ?f i z)")
proof (rule Weierstrass_m_test'_ev)
{
fix t assume t: "t ∈ ball z d"
have "norm z = norm (t + (z - t))" by simp
have "norm (t + (z - t)) ≤ norm t + norm (z - t)" by (rule norm_triangle_ineq)
also from t d have "norm (z - t) < norm z / 2" by (simp add: dist_norm)
finally have A: "norm t > norm z / 2" using z by (simp add: field_simps)

have "norm t = norm (z + (t - z))" by simp
also have "… ≤ norm z + norm (t - z)" by (rule norm_triangle_ineq)
also from t d have "norm (t - z) ≤ norm z / 2" by (simp add: dist_norm norm_minus_commute)
also from z have "… < norm z" by simp
finally have B: "norm t < 2 * norm z" by simp
note A B
} note ball = this

show "eventually (λn. ∀t∈ball z d. norm (?f n t) ≤ 4 * norm z * inverse (of_nat (Suc n)^2)) sequentially"
using eventually_gt_at_top apply eventually_elim
proof safe
fix t :: 'a assume t: "t ∈ ball z d"
from z ball[OF t] have t_nz: "t ≠ 0" by auto
fix n :: nat assume n: "n > nat ⌈4 * norm z⌉"
from ball[OF t] t_nz have "4 * norm z > 2 * norm t" by simp
also from n have "…  < of_nat n" by linarith
finally have n: "of_nat n > 2 * norm t" .
hence "of_nat n > norm t" by simp
hence t': "t ≠ -of_nat (Suc n)" by (intro notI) (simp del: of_nat_Suc)

with t_nz have "?f n t = 1 / (of_nat (Suc n) * (1 + of_nat (Suc n)/t))"
also have "norm … = inverse (of_nat (Suc n)) * inverse (norm (of_nat (Suc n)/t + 1))"
by (simp add: norm_divide norm_mult field_split_simps del: of_nat_Suc)
also {
from z t_nz ball[OF t] have "of_nat (Suc n) / (4 * norm z) ≤ of_nat (Suc n) / (2 * norm t)"
by (intro divide_left_mono mult_pos_pos) simp_all
also have "… < norm (of_nat (Suc n) / t) - norm (1 :: 'a)"
using t_nz n by (simp add: field_simps norm_divide del: of_nat_Suc)
also have "… ≤ norm (of_nat (Suc n)/t + 1)" by (rule norm_diff_ineq)
finally have "inverse (norm (of_nat (Suc n)/t + 1)) ≤ 4 * norm z / of_nat (Suc n)"
using z by (simp add: field_split_simps norm_divide mult_ac del: of_nat_Suc)
}
also have "inverse (real_of_nat (Suc n)) * (4 * norm z / real_of_nat (Suc n)) =
4 * norm z * inverse (of_nat (Suc n)^2)"
by (simp add: field_split_simps power2_eq_square del: of_nat_Suc)
finally show "norm (?f n t) ≤ 4 * norm z * inverse (of_nat (Suc n)^2)"
by (simp del: of_nat_Suc)
qed
next
show "summable (λn. 4 * norm z * inverse ((of_nat (Suc n))^2))"
by (subst summable_Suc_iff) (simp add: summable_mult inverse_power_summable)
qed

subsection ‹The Polygamma functions›

lemma summable_deriv_ln_Gamma:
"z ≠ (0 :: 'a :: {real_normed_field,banach}) ⟹
summable (λn. inverse (of_nat (Suc n)) - inverse (z + of_nat (Suc n)))"
unfolding summable_iff_convergent
by (rule uniformly_convergent_imp_convergent,
rule uniformly_summable_deriv_ln_Gamma[of z "norm z/2"]) simp_all

definition✐‹tag important› Polygamma :: "nat ⇒ ('a :: {real_normed_field,banach}) ⇒ 'a" where
"Polygamma n z = (if n = 0 then
(∑k. inverse (of_nat (Suc k)) - inverse (z + of_nat k)) - euler_mascheroni else
(-1)^Suc n * fact n * (∑k. inverse ((z + of_nat k)^Suc n)))"

abbreviation✐‹tag important› Digamma :: "('a :: {real_normed_field,banach}) ⇒ 'a" where
"Digamma ≡ Polygamma 0"

lemma Digamma_def:
"Digamma z = (∑k. inverse (of_nat (Suc k)) - inverse (z + of_nat k)) - euler_mascheroni"

lemma summable_Digamma:
assumes "(z :: 'a :: {real_normed_field,banach}) ≠ 0"
shows   "summable (λn. inverse (of_nat (Suc n)) - inverse (z + of_nat n))"
proof -
have sums: "(λn. inverse (z + of_nat (Suc n)) - inverse (z + of_nat n)) sums
(0 - inverse (z + of_nat 0))"
by (intro telescope_sums filterlim_compose[OF tendsto_inverse_0]
from summable_add[OF summable_deriv_ln_Gamma[OF assms] sums_summable[OF sums]]
show "summable (λn. inverse (of_nat (Suc n)) - inverse (z + of_nat n))" by simp
qed

lemma summable_offset:
assumes "summable (λn. f (n + k) :: 'a :: real_normed_vector)"
shows   "summable f"
proof -
from assms have "convergent (λm. ∑n<m. f (n + k))"
using summable_iff_convergent by blast
hence "convergent (λm. (∑n<k. f n) + (∑n<m. f (n + k)))"
also have "(λm. (∑n<k. f n) + (∑n<m. f (n + k))) = (λm. ∑n<m+k. f n)"
proof
fix m :: nat
have "{..<m+k} = {..<k} ∪ {k..<m+k}" by auto
also have "(∑n∈…. f n) = (∑n<k. f n) + (∑n=k..<m+k. f n)"
by (rule sum.union_disjoint) auto
also have "(∑n=k..<m+k. f n) = (∑n=0..<m+k-k. f (n + k))"
using sum.shift_bounds_nat_ivl [of f 0 k m] by simp
finally show "(∑n<k. f n) + (∑n<m. f (n + k)) = (∑n<m+k. f n)" by (simp add: atLeast0LessThan)
qed
finally have "(λa. sum f {..<a}) ⇢ lim (λm. sum f {..<m + k})"
by (auto simp: convergent_LIMSEQ_iff dest: LIMSEQ_offset)
thus ?thesis by (auto simp: summable_iff_convergent convergent_def)
qed

lemma Polygamma_converges:
fixes z :: "'a :: {real_normed_field,banach}"
assumes z: "z ≠ 0" and n: "n ≥ 2"
shows "uniformly_convergent_on (ball z d) (λk z. ∑i<k. inverse ((z + of_nat i)^n))"
proof (rule Weierstrass_m_test'_ev)
define e where "e = (1 + d / norm z)"
define m where "m = nat ⌈norm z * e⌉"
{
fix t assume t: "t ∈ ball z d"
have "norm t = norm (z + (t - z))" by simp
also have "… ≤ norm z + norm (t - z)" by (rule norm_triangle_ineq)
also from t have "norm (t - z) < d" by (simp add: dist_norm norm_minus_commute)
finally have "norm t < norm z * e" using z by (simp add: divide_simps e_def)
} note ball = this

show "eventually (λk. ∀t∈ball z d. norm (inverse ((t + of_nat k)^n)) ≤
inverse (of_nat (k - m)^n)) sequentially"
using eventually_gt_at_top[of m] apply eventually_elim
proof (intro ballI)
fix k :: nat and t :: 'a assume k: "k > m" and t: "t ∈ ball z d"
from k have "real_of_nat (k - m) = of_nat k - of_nat m" by (simp add: of_nat_diff)
also have "… ≤ norm (of_nat k :: 'a) - norm z * e"
unfolding m_def by (subst norm_of_nat) linarith
also from ball[OF t] have "… ≤ norm (of_nat k :: 'a) - norm t" by simp
also have "… ≤ norm (of_nat k + t)" by (rule norm_diff_ineq)
finally have "inverse ((norm (t + of_nat k))^n) ≤ inverse (real_of_nat (k - m)^n)" using k n
thus "norm (inverse ((t + of_nat k)^n)) ≤ inverse (of_nat (k - m)^n)"
by (simp add: norm_inverse norm_power power_inverse)
qed

have "summable (λk. inverse ((real_of_nat k)^n))"
using inverse_power_summable[of n] n by simp
hence "summable (λk. inverse ((real_of_nat (k + m - m))^n))" by simp
thus "summable (λk. inverse ((real_of_nat (k - m))^n))" by (rule summable_offset)
qed

lemma Polygamma_converges':
fixes z :: "'a :: {real_normed_field,banach}"
assumes z: "z ≠ 0" and n: "n ≥ 2"
shows "summable (λk. inverse ((z + of_nat k)^n))"
using uniformly_convergent_imp_convergent[OF Polygamma_converges[OF assms, of 1], of z]

theorem Digamma_LIMSEQ:
fixes z :: "'a :: {banach,real_normed_field}"
assumes z: "z ≠ 0"
shows   "(λm. of_real (ln (real m)) - (∑n<m. inverse (z + of_nat n))) ⇢ Digamma z"
proof -
have "(λn. of_real (ln (real n / (real (Suc n))))) ⇢ (of_real (ln 1) :: 'a)"
by (intro tendsto_intros LIMSEQ_n_over_Suc_n) simp_all
hence "(λn. of_real (ln (real n / (real n + 1)))) ⇢ (0 :: 'a)" by (simp add: add_ac)
hence lim: "(λn. of_real (ln (real n)) - of_real (ln (real n + 1))) ⇢ (0::'a)"
proof (rule Lim_transform_eventually)
show "eventually (λn. of_real (ln (real n / (real n + 1))) =
of_real (ln (real n)) - (of_real (ln (real n + 1)) :: 'a)) at_top"
using eventually_gt_at_top[of "0::nat"] by eventually_elim (simp add: ln_div)
qed

from summable_Digamma[OF z]
have "(λn. inverse (of_nat (n+1)) - inverse (z + of_nat n))
sums (Digamma z + euler_mascheroni)"
from sums_diff[OF this euler_mascheroni_sum]
have "(λn. of_real (ln (real (Suc n) + 1)) - of_real (ln (real n + 1)) - inverse (z + of_nat n))
hence "(λm. (∑n<m. of_real (ln (real (Suc n) + 1)) - of_real (ln (real n + 1))) -
(∑n<m. inverse (z + of_nat n))) ⇢ Digamma z"
also have "(λm. (∑n<m. of_real (ln (real (Suc n) + 1)) - of_real (ln (real n + 1)))) =
(λm. of_real (ln (m + 1)) :: 'a)"
by (subst sum_lessThan_telescope) simp_all
finally show ?thesis by (rule Lim_transform) (insert lim, simp)
qed

theorem Polygamma_LIMSEQ:
fixes z :: "'a :: {banach,real_normed_field}"
assumes "z ≠ 0" and "n > 0"
shows   "(λk. inverse ((z + of_nat k)^Suc n)) sums ((-1) ^ Suc n * Polygamma n z / fact n)"
using Polygamma_converges'[OF assms(1), of "Suc n"] assms(2)

theorem has_field_derivative_ln_Gamma_complex [derivative_intros]:
fixes z :: complex
assumes z: "z ∉ ℝ⇩≤⇩0"
shows   "(ln_Gamma has_field_derivative Digamma z) (at z)"
proof -
have not_nonpos_Int [simp]: "t ∉ ℤ⇩≤⇩0" if "Re t > 0" for t
using that by (auto elim!: nonpos_Ints_cases')
from z have z': "z ∉ ℤ⇩≤⇩0" and z'': "z ≠ 0" using nonpos_Ints_subset_nonpos_Reals nonpos_Reals_zero_I
by blast+
let ?f' = "λz k. inverse (of_nat (Suc k)) - inverse (z + of_nat (Suc k))"
let ?f = "λz k. z / of_nat (Suc k) - ln (1 + z / of_nat (Suc k))" and ?F' = "λz. ∑n. ?f' z n"
define d where "d = min (norm z/2) (if Im z = 0 then Re z / 2 else abs (Im z) / 2)"
from z have d: "d > 0" "norm z/2 ≥ d" by (auto simp add: complex_nonpos_Reals_iff d_def)
have ball: "Im t = 0 ⟶ Re t > 0" if "dist z t < d" for t
using no_nonpos_Real_in_ball[OF z, of t] that unfolding d_def by (force simp add: complex_nonpos_Reals_iff)
have sums: "(λn. inverse (z + of_nat (Suc n)) - inverse (z + of_nat n)) sums
(0 - inverse (z + of_nat 0))"
by (intro telescope_sums filterlim_compose[OF tendsto_inverse_0]

have "((λz. ∑n. ?f z n) has_field_derivative ?F' z) (at z)"
using d z ln_Gamma_series'_aux[OF z']
apply (intro has_field_derivative_series'(2)[of "ball z d" _ _ z] uniformly_summable_deriv_ln_Gamma)
apply (auto intro!: derivative_eq_intros add_pos_pos mult_pos_pos dest!: ball
simp: field_simps sums_iff nonpos_Reals_divide_of_nat_iff
simp del: of_nat_Suc)
done
with z have "((λz. (∑k. ?f z k) - euler_mascheroni * z - Ln z) has_field_derivative
?F' z - euler_mascheroni - inverse z) (at z)"
by (force intro!: derivative_eq_intros simp: Digamma_def)
also have "?F' z - euler_mascheroni - inverse z = (?F' z + -inverse z) - euler_mascheroni" by simp
also from sums have "-inverse z = (∑n. inverse (z + of_nat (Suc n)) - inverse (z + of_nat n))"
also from sums summable_deriv_ln_Gamma[OF z'']
have "?F' z + … =  (∑n. inverse (of_nat (Suc n)) - inverse (z + of_nat n))"
also have "… - euler_mascheroni = Digamma z" by (simp add: Digamma_def)
finally have "((λz. (∑k. ?f z k) - euler_mascheroni * z - Ln z)
has_field_derivative Digamma z) (at z)" .
moreover from eventually_nhds_ball[OF d(1), of z]
have "eventually (λz. ln_Gamma z = (∑k. ?f z k) - euler_mascheroni * z - Ln z) (nhds z)"
proof eventually_elim
fix t assume "t ∈ ball z d"
hence "t ∉ ℤ⇩≤⇩0" by (auto dest!: ball elim!: nonpos_Ints_cases)
from ln_Gamma_series'_aux[OF this]
show "ln_Gamma t = (∑k. ?f t k) - euler_mascheroni * t - Ln t" by (simp add: sums_iff)
qed
ultimately show ?thesis by (subst DERIV_cong_ev[OF refl _ refl])
qed

declare has_field_derivative_ln_Gamma_complex[THEN DERIV_chain2, derivative_intros]

lemma Digamma_1 [simp]: "Digamma (1 :: 'a :: {real_normed_field,banach}) = - euler_mascheroni"

lemma Digamma_plus1:
assumes "z ≠ 0"
shows   "Digamma (z+1) = Digamma z + 1/z"
proof -
have sums: "(λk. inverse (z + of_nat k) - inverse (z + of_nat (Suc k)))
sums (inverse (z + of_nat 0) - 0)"
by (intro telescope_sums'[OF filterlim_compose[OF tendsto_inverse_0]]
have "Digamma (z+1) = (∑k. inverse (of_nat (Suc k)) - inverse (z + of_nat (Suc k))) -
euler_mascheroni" (is "_ = suminf ?f - _") by (simp add: Digamma_def add_ac)
also have "suminf ?f = (∑k. inverse (of_nat (Suc k)) - inverse (z + of_nat k)) +
(∑k. inverse (z + of_nat k) - inverse (z + of_nat (Suc k)))"
also have "(∑k. inverse (z + of_nat k) - inverse (z + of_nat (Suc k))) = 1/z"
using sums by (simp add: sums_iff inverse_eq_divide)
finally show ?thesis by (simp add: Digamma_def[of z])
qed

theorem Polygamma_plus1:
assumes "z ≠ 0"
shows   "Polygamma n (z + 1) = Polygamma n z + (-1)^n * fact n / (z ^ Suc n)"
proof (cases "n = 0")
assume n: "n ≠ 0"
let ?f = "λk. inverse ((z + of_nat k) ^ Suc n)"
have "Polygamma n (z + 1) = (-1) ^ Suc n * fact n * (∑k. ?f (k+1))"
also have "(∑k. ?f (k+1)) + (∑k<1. ?f k) = (∑k. ?f k)"
using Polygamma_converges'[OF assms, of "Suc n"] n
by (subst suminf_split_initial_segment [symmetric]) simp_all
hence "(∑k. ?f (k+1)) = (∑k. ?f k) - inverse (z ^ Suc n)" by (simp add: algebra_simps)
also have "(-1) ^ Suc n * fact n * ((∑k. ?f k) - inverse (z ^ Suc n)) =
Polygamma n z + (-1)^n * fact n / (z ^ Suc n)" using n
by (simp add: inverse_eq_divide algebra_simps Polygamma_def)
finally show ?thesis .
qed (insert assms, simp add: Digamma_plus1 inverse_eq_divide)

theorem Digamma_of_nat:
"Digamma (of_nat (Suc n) :: 'a :: {real_normed_field,banach}) = harm n - euler_mascheroni"
proof (induction n)
case (Suc n)
have "Digamma (of_nat (Suc (Suc n)) :: 'a) = Digamma (of_nat (Suc n) + 1)" by simp
also have "… = Digamma (of_nat (Suc n)) + inverse (of_nat (Suc n))"
by (subst Digamma_plus1) (simp_all add: inverse_eq_divide del: of_nat_Suc)
also have "Digamma (of_nat (Suc n) :: 'a) = harm n - euler_mascheroni " by (rule Suc)
also have "… + inverse (of_nat (Suc n)) = harm (Suc n) - euler_mascheroni"
finally show ?case .

lemma Digamma_numeral: "Digamma (numeral n) = harm (pred_numeral n) - euler_mascheroni"
by (subst of_nat_numeral[symmetric], subst numeral_eq_Suc, subst Digamma_of_nat) (rule refl)

lemma Polygamma_of_real: "x ≠ 0 ⟹ Polygamma n (of_real x) = of_real (Polygamma n x)"
unfolding Polygamma_def using summable_Digamma[of x] Polygamma_converges'[of x "Suc n"]

lemma Polygamma_Real: "z ∈ ℝ ⟹ z ≠ 0 ⟹ Polygamma n z ∈ ℝ"
by (elim Reals_cases, hypsubst, subst Polygamma_of_real) simp_all

lemma Digamma_half_integer:
"Digamma (of_nat n + 1/2 :: 'a :: {real_normed_field,banach}) =
(∑k<n. 2 / (of_nat (2*k+1))) - euler_mascheroni - of_real (2 * ln 2)"
proof (induction n)
case 0
have "Digamma (1/2 :: 'a) = of_real (Digamma (1/2))" by (simp add: Polygamma_of_real [symmetric])
also have "Digamma (1/2::real) =
(∑k. inverse (of_nat (Suc k)) - inverse (of_nat k + 1/2)) - euler_mascheroni"
also have "(∑k. inverse (of_nat (Suc k) :: real) - inverse (of_nat k + 1/2)) =
(∑k. inverse (1/2) * (inverse (2 * of_nat (Suc k)) - inverse (2 * of_nat k + 1)))"
also have "… = - 2 * ln 2" using sums_minus[OF alternating_harmonic_series_sums']
by (subst suminf_mult) (simp_all add: algebra_simps sums_iff)
finally show ?case by simp
next
case (Suc n)
have nz: "2 * of_nat n + (1:: 'a) ≠ 0"
hence nz': "of_nat n + (1/2::'a) ≠ 0" by (simp add: field_simps)
have "Digamma (of_nat (Suc n) + 1/2 :: 'a) = Digamma (of_nat n + 1/2 + 1)" by simp
also from nz' have "… = Digamma (of_nat n + 1/2) + 1 / (of_nat n + 1/2)"
by (rule Digamma_plus1)
also from nz nz' have "1 / (of_nat n + 1/2 :: 'a) = 2 / (2 * of_nat n + 1)"
by (subst divide_eq_eq) simp_all
also note Suc
qed

lemma Digamma_one_half: "Digamma (1/2) = - euler_mascheroni - of_real (2 * ln 2)"
using Digamma_half_integer[of 0] by simp

lemma Digamma_real_three_halves_pos: "Digamma (3/2 :: real) > 0"
proof -
have "-Digamma (3/2 :: real) = -Digamma (of_nat 1 + 1/2)" by simp
also have "… = 2 * ln 2 + euler_mascheroni - 2" by (subst Digamma_half_integer) simp
also note euler_mascheroni_less_13_over_22
also note ln2_le_25_over_36
finally show ?thesis by simp
qed

theorem has_field_derivative_Polygamma [derivative_intros]:
fixes z :: "'a :: {real_normed_field,euclidean_space}"
assumes z: "z ∉ ℤ⇩≤⇩0"
shows "(Polygamma n has_field_derivative Polygamma (Suc n) z) (at z within A)"
proof (rule has_field_derivative_at_within, cases "n = 0")
assume n: "n = 0"
let ?f = "λk z. inverse (of_nat (Suc k)) - inverse (z + of_nat k)"
let ?F = "λz. ∑k. ?f k z" and ?f' = "λk z. inverse ((z + of_nat k)⇧2)"
from no_nonpos_Int_in_ball'[OF z] obtain d where d: "0 < d" "⋀t. t ∈ ball z d ⟹ t ∉ ℤ⇩≤⇩0"
by auto
from z have summable: "summable (λk. inverse (of_nat (Suc k)) - inverse (z + of_nat k))"
by (intro summable_Digamma) force
from z have conv: "uniformly_convergent_on (ball z d) (λk z. ∑i<k. inverse ((z + of_nat i)⇧2))"
by (intro Polygamma_converges) auto
with d have "summable (λk. inverse ((z + of_nat k)⇧2))" unfolding summable_iff_convergent
by (auto dest!: uniformly_convergent_imp_convergent simp: summable_iff_convergent )

have "(?F has_field_derivative (∑k. ?f' k z)) (at z)"
proof (rule has_field_derivative_series'[of "ball z d" _ _ z])
fix k :: nat and t :: 'a assume t: "t ∈ ball z d"
from t d(2)[of t] show "((λz. ?f k z) has_field_derivative ?f' k t) (at t within ball z d)"
by (auto intro!: derivative_eq_intros simp: power2_eq_square simp del: of_nat_Suc
dest!: plus_of_nat_eq_0_imp elim!: nonpos_Ints_cases)
qed (insert d(1) summable conv, (assumption|simp)+)
with z show "(Polygamma n has_field_derivative Polygamma (Suc n) z) (at z)"
unfolding Digamma_def [abs_def] Polygamma_def [abs_def] using n
by (force simp: power2_eq_square intro!: derivative_eq_intros)
next
assume n: "n ≠ 0"
from z have z': "z ≠ 0" by auto
from no_nonpos_Int_in_ball'[OF z] obtain d where d: "0 < d" "⋀t. t ∈ ball z d ⟹ t ∉ ℤ⇩≤⇩0"
by auto
define n' where "n' = Suc n"
from n have n': "n' ≥ 2" by (simp add: n'_def)
have "((λz. ∑k. inverse ((z + of_nat k) ^ n')) has_field_derivative
(∑k. - of_nat n' * inverse ((z + of_nat k) ^ (n'+1)))) (at z)"
proof (rule has_field_derivative_series'[of "ball z d" _ _ z])
fix k :: nat and t :: 'a assume t: "t ∈ ball z d"
with d have t': "t ∉ ℤ⇩≤⇩0" "t ≠ 0" by auto
show "((λa. inverse ((a + of_nat k) ^ n')) has_field_derivative
- of_nat n' * inverse ((t + of_nat k) ^ (n'+1))) (at t within ball z d)" using t'
by (fastforce intro!: derivative_eq_intros simp: divide_simps power_diff dest: plus_of_nat_eq_0_imp)
next
have "uniformly_convergent_on (ball z d)
(λk z. (- of_nat n' :: 'a) * (∑i<k. inverse ((z + of_nat i) ^ (n'+1))))"
using z' n by (intro uniformly_convergent_mult Polygamma_converges) (simp_all add: n'_def)
thus "uniformly_convergent_on (ball z d)
(λk z. ∑i<k. - of_nat n' * inverse ((z + of_nat i :: 'a) ^ (n'+1)))"
by (subst (asm) sum_distrib_left) simp
qed (insert Polygamma_converges'[OF z' n'] d, simp_all)
also have "(∑k. - of_nat n' * inverse ((z + of_nat k) ^ (n' + 1))) =
(- of_nat n') * (∑k. inverse ((z + of_nat k) ^ (n' + 1)))"
using Polygamma_converges'[OF z', of "n'+1"] n' by (subst suminf_mult) simp_all
finally have "((λz. ∑k. inverse ((z + of_nat k) ^ n')) has_field_derivative
- of_nat n' * (∑k. inverse ((z + of_nat k) ^ (n' + 1)))) (at z)" .
from DERIV_cmult[OF this, of "(-1)^Suc n * fact n :: 'a"]
show "(Polygamma n has_field_derivative Polygamma (Suc n) z) (at z)"
unfolding n'_def Polygamma_def[abs_def] using n by (simp add: algebra_simps)
qed

declare has_field_derivative_Polygamma[THEN DERIV_chain2, derivative_intros]

lemma isCont_Polygamma [continuous_intros]:
fixes f :: "_ ⇒ 'a :: {real_normed_field,euclidean_space}"
shows "isCont f z ⟹ f z ∉ ℤ⇩≤⇩0 ⟹ isCont (λx. Polygamma n (f x)) z"
by (rule isCont_o2[OF _  DERIV_isCont[OF has_field_derivative_Polygamma]])

lemma continuous_on_Polygamma:
"A ∩ ℤ⇩≤⇩0 = {} ⟹ continuous_on A (Polygamma n :: _ ⇒ 'a :: {real_normed_field,euclidean_space})"
by (intro continuous_at_imp_continuous_on isCont_Polygamma[OF continuous_ident] ballI) blast

lemma isCont_ln_Gamma_complex [continuous_intros]:
fixes f :: "'a::t2_space ⇒ complex"
shows "isCont f z ⟹ f z ∉ ℝ⇩≤⇩0 ⟹ isCont (λz. ln_Gamma (f z)) z"
by (rule isCont_o2[OF _  DERIV_isCont[OF has_field_derivative_ln_Gamma_complex]])

lemma continuous_on_ln_Gamma_complex [continuous_intros]:
fixes A :: "complex set"
shows "A ∩ ℝ⇩≤⇩0 = {} ⟹ continuous_on A ln_Gamma"
by (intro continuous_at_imp_continuous_on ballI isCont_ln_Gamma_complex[OF continuous_ident])
fastforce

lemma deriv_Polygamma:
assumes "z ∉ ℤ⇩≤⇩0"
shows   "deriv (Polygamma m) z =
Polygamma (Suc m) (z :: 'a :: {real_normed_field,euclidean_space})"
by (intro DERIV_imp_deriv has_field_derivative_Polygamma assms)
thm has_field_derivative_Polygamma

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)"
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 deriv_ln_Gamma_complex:
assumes "z ∉ ℝ⇩≤⇩0"
shows   "deriv ln_Gamma z = Digamma (z :: complex)"
by (intro DERIV_imp_deriv has_field_derivative_ln_Gamma_complex assms)

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

text ‹
We define a type class that captures all the fundamental properties of the inverse of the Gamma function
and defines the Gamma function upon that. This allows us to instantiate the type class both for
the reals and for the complex numbers with a minimal amount of proof duplication.
›

class✐‹tag unimportant› Gamma = real_normed_field + complete_space +
fixes rGamma :: "'a ⇒ 'a"
assumes rGamma_eq_zero_iff_aux: "rGamma z = 0 ⟷ (∃n. z = - of_nat n)"
assumes differentiable_rGamma_aux1:
"(⋀n. z ≠ - of_nat n) ⟹
let d = (THE d. (λn. ∑k<n. inverse (of_nat (Suc k)) - inverse (z + of_nat k))
⇢ d) - scaleR euler_mascheroni 1
in  filterlim (λy. (rGamma y - rGamma z + rGamma z * d * (y - z)) /⇩R
norm (y - z)) (nhds 0) (at z)"
assumes differentiable_rGamma_aux2:
"let z = - of_na```