--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Multivariate_Analysis/Gamma.thy Mon Jan 04 17:45:36 2016 +0100
@@ -0,0 +1,2640 @@
+(*
+ Title: HOL/Multivariate_Analysis/Gamma.thy
+ Author: Manuel Eberl, TU München
+
+ 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ß 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"}.
+*)
+theory Gamma
+imports
+ Complex_Transcendental
+ Summation
+ Harmonic_Numbers
+ "~~/src/HOL/Library/Nonpos_Ints"
+ "~~/src/HOL/Library/Periodic_Fun"
+begin
+
+lemma pochhammer_eq_0_imp_nonpos_Int:
+ "pochhammer (x::'a::field_char_0) n = 0 \<Longrightarrow> x \<in> \<int>\<^sub>\<le>\<^sub>0"
+ by (auto simp: pochhammer_eq_0_iff)
+
+lemma closed_nonpos_Ints [simp]: "closed (\<int>\<^sub>\<le>\<^sub>0 :: 'a :: real_normed_algebra_1 set)"
+proof -
+ have "\<int>\<^sub>\<le>\<^sub>0 = (of_int ` {n. n \<le> 0} :: 'a set)"
+ by (auto elim!: nonpos_Ints_cases intro!: nonpos_Ints_of_int)
+ also have "closed \<dots>" by (rule closed_of_int_image)
+ finally show ?thesis .
+qed
+
+lemma plus_one_in_nonpos_Ints_imp: "z + 1 \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> z \<in> \<int>\<^sub>\<le>\<^sub>0"
+ using nonpos_Ints_diff_Nats[of "z+1" "1"] by simp_all
+
+lemma fraction_not_in_ints:
+ assumes "\<not>(n dvd m)" "n \<noteq> 0"
+ shows "of_int m / of_int n \<notin> (\<int> :: 'a :: {division_ring,ring_char_0} set)"
+proof
+ assume "of_int m / (of_int n :: 'a) \<in> \<int>"
+ 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: divide_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 not_in_Ints_imp_not_in_nonpos_Ints: "z \<notin> \<int> \<Longrightarrow> z \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ by (auto simp: Ints_def nonpos_Ints_def)
+
+lemma double_in_nonpos_Ints_imp:
+ assumes "2 * (z :: 'a :: field_char_0) \<in> \<int>\<^sub>\<le>\<^sub>0"
+ shows "z \<in> \<int>\<^sub>\<le>\<^sub>0 \<or> z + 1/2 \<in> \<int>\<^sub>\<le>\<^sub>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: "(\<lambda>n. ((-1)^n / fact (2*n+1)) *\<^sub>R z^(2*n+1)) sums sin z"
+proof -
+ from sin_converges[of z] have "(\<lambda>n. sin_coeff n *\<^sub>R z^n) sums sin z" .
+ also have "(\<lambda>n. sin_coeff n *\<^sub>R z^n) sums sin z \<longleftrightarrow>
+ (\<lambda>n. ((-1)^n / fact (2*n+1)) *\<^sub>R z^(2*n+1)) sums sin z"
+ by (subst sums_mono_reindex[of "\<lambda>n. 2*n+1", symmetric])
+ (auto simp: sin_coeff_def subseq_def ac_simps elim!: oddE)
+ finally show ?thesis .
+qed
+
+lemma cos_series: "(\<lambda>n. ((-1)^n / fact (2*n)) *\<^sub>R z^(2*n)) sums cos z"
+proof -
+ from cos_converges[of z] have "(\<lambda>n. cos_coeff n *\<^sub>R z^n) sums cos z" .
+ also have "(\<lambda>n. cos_coeff n *\<^sub>R z^n) sums cos z \<longleftrightarrow>
+ (\<lambda>n. ((-1)^n / fact (2*n)) *\<^sub>R z^(2*n)) sums cos z"
+ by (subst sums_mono_reindex[of "\<lambda>n. 2*n", symmetric])
+ (auto simp: cos_coeff_def subseq_def ac_simps elim!: evenE)
+ finally show ?thesis .
+qed
+
+lemma sin_z_over_z_series:
+ fixes z :: "'a :: {real_normed_field,banach}"
+ assumes "z \<noteq> 0"
+ shows "(\<lambda>n. (-1)^n / fact (2*n+1) * z^(2*n)) sums (sin z / z)"
+proof -
+ from sin_series[of z] have "(\<lambda>n. z * ((-1)^n / fact (2*n+1)) * z^(2*n)) sums sin z"
+ by (simp add: field_simps scaleR_conv_of_real)
+ from sums_mult[OF this, of "inverse z"] and assms show ?thesis
+ by (simp add: field_simps)
+qed
+
+lemma sin_z_over_z_series':
+ fixes z :: "'a :: {real_normed_field,banach}"
+ assumes "z \<noteq> 0"
+ shows "(\<lambda>n. sin_coeff (n+1) *\<^sub>R z^n) sums (sin z / z)"
+proof -
+ from sums_split_initial_segment[OF sin_converges[of z], of 1]
+ have "(\<lambda>n. z * (sin_coeff (n+1) *\<^sub>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 "((\<lambda>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 "((\<lambda>z::'a. \<Sum>n. of_real (sin_coeff (n+1)) * z^n)
+ has_field_derivative (\<Sum>n. diffs (\<lambda>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 (\<lambda>n. of_real (sin_coeff (n+1)) * (1::'a)^n)" by (simp add: of_real_def)
+ qed simp
+ also have "(\<lambda>z::'a. \<Sum>n. of_real (sin_coeff (n+1)) * z^n) = ?f"
+ proof
+ fix z
+ show "(\<Sum>n. of_real (sin_coeff (n+1)) * z^n) = ?f z"
+ by (cases "z = 0") (insert sin_z_over_z_series'[of z],
+ simp_all add: scaleR_conv_of_real sums_iff powser_zero sin_coeff_def)
+ qed
+ also have "(\<Sum>n. diffs (\<lambda>n. of_real (sin_coeff (n + 1))) n * (0::'a) ^ n) =
+ diffs (\<lambda>n. of_real (sin_coeff (Suc n))) 0" by (simp add: powser_zero)
+ also have "\<dots> = 0" by (simp add: sin_coeff_def diffs_def)
+ finally show "((\<lambda>z::'a. if z = 0 then 1 else sin z / z) has_field_derivative 0) (at 0)" .
+qed
+
+
+lemma tendsto_complex_of_real_iff: "((\<lambda>x. complex_of_real (f x)) \<longlongrightarrow> of_real c) F = (f \<longlongrightarrow> c) F"
+ by (rule tendsto_of_real_iff)
+
+lemma round_Re_minimises_norm:
+ "norm ((z::complex) - of_int m) \<ge> 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)\<^sup>2 + (Im z)\<^sup>2)"
+ by (simp add: cmod_def)
+ also have "\<bar>Re z - of_int ?n\<bar> \<le> \<bar>Re z - of_int m\<bar>" by (rule round_diff_minimal)
+ hence "sqrt ((Re z - of_int ?n)\<^sup>2 + (Im z)\<^sup>2) \<le> sqrt ((Re z - of_int m)\<^sup>2 + (Im z)\<^sup>2)"
+ by (intro real_sqrt_le_mono add_mono) (simp_all add: abs_le_square_iff)
+ also have "\<dots> = norm (z - of_int m)" by (simp add: cmod_def)
+ finally show ?thesis .
+qed
+
+lemma Re_pos_in_ball:
+ assumes "Re z > 0" "t \<in> ball z (Re z/2)"
+ shows "Re t > 0"
+proof -
+ have "Re (z - t) \<le> norm (z - t)" by (rule complex_Re_le_cmod)
+ also from assms have "\<dots> < 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 \<in> ball z (Re z/2)"
+ shows "t \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ using Re_pos_in_ball[OF assms] by (force elim!: nonpos_Ints_cases)
+
+lemma no_nonpos_Int_in_ball:
+ assumes "t \<in> ball z (dist z (round (Re z)))"
+ shows "t \<notin> \<int>\<^sub>\<le>\<^sub>0"
+proof
+ assume "t \<in> \<int>\<^sub>\<le>\<^sub>0"
+ then obtain n where "t = of_int n" by (auto elim!: nonpos_Ints_cases)
+ have "dist z (of_int n) \<le> 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 "\<dots> \<le> 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}) \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ obtains d where "d > 0" "\<And>t. t \<in> ball z d \<Longrightarrow> t \<notin> \<int>\<^sub>\<le>\<^sub>0"
+proof (rule that)
+ from assms show "setdist {z} \<int>\<^sub>\<le>\<^sub>0 > 0" by (subst setdist_gt_0_compact_closed) auto
+next
+ fix t assume "t \<in> ball z (setdist {z} \<int>\<^sub>\<le>\<^sub>0)"
+ thus "t \<notin> \<int>\<^sub>\<le>\<^sub>0" using setdist_le_dist[of z "{z}" t "\<int>\<^sub>\<le>\<^sub>0"] by force
+qed
+
+lemma Re_pos_or_Im_nz_in_ball:
+ assumes "Re z > 0 \<or> Im z \<noteq> 0" "t \<in> ball z (if Im z = 0 then Re z / 2 else abs (Im z) / 2)"
+ shows "Re t > 0 \<or> Im t \<noteq> 0"
+using assms(1)
+proof (cases "Im z = 0")
+ assume A: "Im z = 0"
+ with assms(1) have "Re z > 0" by blast
+ with assms(2) A Re_pos_in_ball[of z t] show ?thesis by simp
+next
+ assume A: "Im z \<noteq> 0"
+ have "abs (Im z) - abs (Im t) \<le> abs (Im z - Im t)" by linarith
+ also have "\<dots> = abs (Im (z - t))" by simp
+ also have "\<dots> \<le> norm (z - t)" by (rule abs_Im_le_cmod)
+ also from A assms(2) have "\<dots> \<le> abs (Im z) / 2" by (simp add: dist_complex_def)
+ finally have "abs (Im t) > 0" using A by simp
+ thus ?thesis by force
+qed
+
+
+subsection \<open>Definitions\<close>
+
+text \<open>
+ We define the Gamma function by first defining its multiplicative inverse @{term "Gamma_inv"}.
+ This is more convenient because @{term "Gamma_inv"} is entire, which makes proofs of its
+ properties more convenient because one does not have to watch out for discontinuities.
+ (e.g. @{term "Gamma_inv"} fulfils @{term "rGamma z = z * rGamma (z + 1)"} everywhere,
+ whereas @{term "Gamma"} does not fulfil the analogous equation on the non-positive integers)
+
+ We define the Gamma function (resp. its inverse) 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 @{term "Gamma (z + 1) = z * Gamma z"} follows
+ immediately from the definition.
+\<close>
+
+definition Gamma_series :: "('a :: {banach,real_normed_field}) \<Rightarrow> nat \<Rightarrow> '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}) \<Rightarrow> nat \<Rightarrow> '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}) \<Rightarrow> nat \<Rightarrow> '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 (\<lambda>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 (\<lambda>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 (\<lambda>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 \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> rGamma_series z \<longlonglongrightarrow> 0"
+ by (elim nonpos_Ints_cases', hypsubst, subst tendsto_cong, rule rGamma_series_minus_of_nat, simp)
+
+lemma Gamma_series_nonpos_Ints_LIMSEQ: "z \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Gamma_series z \<longlonglongrightarrow> 0"
+ by (elim nonpos_Ints_cases', hypsubst, subst tendsto_cong, rule Gamma_series_minus_of_nat, simp)
+
+lemma Gamma_series'_nonpos_Ints_LIMSEQ: "z \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Gamma_series' z \<longlonglongrightarrow> 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 \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ shows "(\<lambda>n. Gamma_series' z n / Gamma_series z n) \<longlonglongrightarrow> 1"
+proof (rule Lim_transform_eventually)
+ from eventually_gt_at_top[of "0::nat"]
+ show "eventually (\<lambda>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 "\<dots> = z / of_nat n + 1" by (simp add: divide_simps)
+ finally show "z / of_nat n + 1 = Gamma_series' z n / Gamma_series z n" ..
+ qed
+ have "(\<lambda>x. z / of_nat x) \<longlonglongrightarrow> 0"
+ by (rule tendsto_norm_zero_cancel)
+ (insert tendsto_mult[OF tendsto_const[of "norm z"] lim_inverse_n],
+ simp add: norm_divide inverse_eq_divide)
+ from tendsto_add[OF this tendsto_const[of 1]] show "(\<lambda>n. z / of_nat n + 1) \<longlonglongrightarrow> 1" by simp
+qed
+
+
+subsection \<open>Convergence of the Euler series form\<close>
+
+text \<open>
+ We now show that the series that defines the Gamma 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 Gamma function, and from that to the Gamma function itself.
+\<close>
+
+definition ln_Gamma_series :: "('a :: {banach,real_normed_field,ln}) \<Rightarrow> nat \<Rightarrow> 'a" where
+ "ln_Gamma_series z n = z * ln (of_nat n) - ln z - (\<Sum>k=1..n. ln (z / of_nat k + 1))"
+
+definition ln_Gamma_series' :: "('a :: {banach,real_normed_field,ln}) \<Rightarrow> nat \<Rightarrow> 'a" where
+ "ln_Gamma_series' z n =
+ - euler_mascheroni*z - ln z + (\<Sum>k=1..n. z / of_nat n - ln (z / of_nat k + 1))"
+
+definition ln_Gamma :: "('a :: {banach,real_normed_field,ln}) \<Rightarrow> 'a" where
+ "ln_Gamma z = lim (ln_Gamma_series z)"
+
+
+text \<open>
+ 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, adapting this
+ proof:
+ http://math.stackexchange.com/questions/887158/convergence-of-gammaz-lim-n-to-infty-fracnzn-prod-m-0nzm
+\<close>
+
+context
+begin
+
+private lemma ln_Gamma_series_complex_converges_aux:
+ fixes z :: complex and k :: nat
+ assumes z: "z \<noteq> 0" and k: "of_nat k \<ge> 2*norm z" "k \<ge> 2"
+ shows "norm (z * ln (1 - 1/of_nat k) + ln (z/of_nat k + 1)) \<le> 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)"
+ by (simp add: algebra_simps)
+ also have "norm ... \<le> ?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)) \<le> (norm (-1/?k))\<^sup>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) \<le> ?z * ((norm (1/?k))^2 / (1 - norm (1/?k)))"
+ by (intro mult_left_mono) simp_all
+ also have "... \<le> (?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 "... \<le> (?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) \<le> 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) \<le> ?z^2 / of_nat k^2 / (1 - ?z / of_nat k)"
+ by (simp add: field_simps norm_divide)
+ also have "... \<le> (?z^2 * (of_nat k / (of_nat k - ?z))) / of_nat k^2" using k
+ by (simp add: field_simps power2_eq_square)
+ also have "... \<le> (?z^2 * 2) / of_nat k^2" using k
+ by (intro divide_right_mono mult_left_mono) (simp_all add: field_simps)
+ also note add_divide_distrib [symmetric]
+ finally show ?thesis by (simp only: distrib_left mult.commute)
+qed
+
+lemma ln_Gamma_series_complex_converges:
+ assumes z: "z \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ assumes d: "d > 0" "\<And>n. n \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> norm (z - of_int n) > d"
+ shows "uniformly_convergent_on (ball z d) (\<lambda>n z. ln_Gamma_series z n :: complex)"
+proof (intro Cauchy_uniformly_convergent uniformly_Cauchy_onI')
+ fix e :: real assume e: "e > 0"
+ def e'' \<equiv> "SUP t:ball z d. norm t + norm t^2"
+ def e' \<equiv> "e / (2*e'')"
+ have "bounded ((\<lambda>t. norm t + norm t^2) ` cball z d)"
+ by (intro compact_imp_bounded compact_continuous_image) (auto intro!: continuous_intros)
+ hence "bounded ((\<lambda>t. norm t + norm t^2) ` ball z d)" by (rule bounded_subset) auto
+ hence bdd: "bdd_above ((\<lambda>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 \<le> e''" if "t \<in> 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
+ by (intro divide_pos_pos mult_pos_pos add_pos_pos) (simp_all add: field_simps)
+
+ have "summable (\<lambda>k. inverse ((real_of_nat k)^2))"
+ by (rule inverse_power_summable) simp
+ from summable_partial_sum_bound[OF this e'] guess M . note M = this
+
+ def N \<equiv> "max 2 (max (nat \<lceil>2 * (norm z + d)\<rceil>) M)"
+ {
+ from d have "\<lceil>2 * (cmod z + d)\<rceil> \<ge> \<lceil>0::real\<rceil>"
+ by (intro ceiling_mono mult_nonneg_nonneg add_nonneg_nonneg) simp_all
+ hence "2 * (norm z + d) \<le> of_nat (nat \<lceil>2 * (norm z + d)\<rceil>)" unfolding N_def
+ by (simp_all add: le_of_int_ceiling)
+ also have "... \<le> of_nat N" unfolding N_def
+ by (subst of_nat_le_iff) (rule max.coboundedI2, rule max.cobounded1)
+ finally have "of_nat N \<ge> 2 * (norm z + d)" .
+ moreover have "N \<ge> 2" "N \<ge> M" unfolding N_def by simp_all
+ moreover have "(\<Sum>k=m..n. 1/(of_nat k)\<^sup>2) < e'" if "m \<ge> N" for m n
+ using M[OF order.trans[OF `N \<ge> M` that]] unfolding real_norm_def
+ by (subst (asm) abs_of_nonneg) (auto intro: setsum_nonneg simp: divide_simps)
+ moreover note calculation
+ } note N = this
+
+ show "\<exists>M. \<forall>t\<in>ball z d. \<forall>m\<ge>M. \<forall>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 \<in> ball z d" and mn: "m \<ge> 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 \<le> dist 0 t" using dist_triangle[of 0 z t]
+ by (simp add: dist_commute)
+ finally have t_nz: "t \<noteq> 0" by auto
+
+ have "norm t \<le> 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) \<le> of_nat N" by (rule N)
+ also have "N \<le> 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)))) +
+ ((\<Sum>k=1..n. Ln (t / of_nat k + 1)) - (\<Sum>k=1..m. Ln (t / of_nat k + 1)))"
+ by (simp add: ln_Gamma_series_def algebra_simps)
+ also have "(\<Sum>k=1..n. Ln (t / of_nat k + 1)) - (\<Sum>k=1..m. Ln (t / of_nat k + 1)) =
+ (\<Sum>k\<in>{1..n}-{1..m}. Ln (t / of_nat k + 1))" using mn
+ by (simp add: setsum_diff)
+ 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))) =
+ (\<Sum>k = Suc m..n. t * Ln (of_nat (k - 1)) - t * Ln (of_nat k))" using mn
+ by (subst setsum_telescope'' [symmetric]) simp_all
+ also have "... = (\<Sum>k = Suc m..n. t * Ln (of_nat (k - 1) / of_nat k))" using mn N
+ by (intro setsum_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 \<in> {Suc m..n}" for k
+ using that of_nat_eq_0_iff[of "Suc i" for i] by (cases k) (simp_all add: divide_simps)
+ hence "(\<Sum>k = Suc m..n. t * Ln (of_nat (k - 1) / of_nat k)) =
+ (\<Sum>k = Suc m..n. t * Ln (1 - 1 / of_nat k))" using mn N
+ by (intro setsum.cong) simp_all
+ also note setsum.distrib [symmetric]
+ also have "norm (\<Sum>k=Suc m..n. t * Ln (1 - 1/of_nat k) + Ln (t/of_nat k + 1)) \<le>
+ (\<Sum>k=Suc m..n. 2 * (norm t + (norm t)\<^sup>2) / (real_of_nat k)\<^sup>2)" using t_nz N(2) mn norm_t
+ by (intro order.trans[OF norm_setsum setsum_mono[OF ln_Gamma_series_complex_converges_aux]]) simp_all
+ also have "... \<le> 2 * (norm t + norm t^2) * (\<Sum>k=Suc m..n. 1 / (of_nat k)\<^sup>2)"
+ by (simp add: setsum_right_distrib)
+ 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)\<^sup>2) / e'')"
+ by (simp add: e'_def field_simps power2_eq_square)
+ also from e''[OF t] e''_pos e
+ have "\<dots> \<le> 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) \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ shows "\<exists>d>0. uniformly_convergent_on (ball z d) (\<lambda>n z. ln_Gamma_series z n)"
+proof -
+ def d' \<equiv> "Re z"
+ def d \<equiv> "if d' > 0 then d' / 2 else norm (z - of_int (round d')) / 2"
+ have "of_int (round d') \<in> \<int>\<^sub>\<le>\<^sub>0" if "d' \<le> 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 \<in> \<int>\<^sub>\<le>\<^sub>0" for n
+ proof (cases "Re z > 0")
+ case True
+ from nonpos_Ints_nonpos[OF that] have n: "n \<le> 0" by simp
+ from True have "d = Re z/2" by (simp add: d_def d'_def)
+ also from n True have "\<dots> < Re (z - of_int n)" by simp
+ also have "\<dots> \<le> 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 \<noteq> 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 "\<dots> \<le> 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) (\<lambda>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) \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> convergent (ln_Gamma_series z)"
+ by (drule ln_Gamma_series_complex_converges') (auto intro: uniformly_convergent_imp_convergent)
+
+lemma ln_Gamma_complex_LIMSEQ: "(z :: complex) \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> ln_Gamma_series z \<longlonglongrightarrow> 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 \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ shows "exp (ln_Gamma_series z n :: complex) = Gamma_series z n"
+proof -
+ from assms have "z \<noteq> 0" by (intro notI) auto
+ with assms have "exp (ln_Gamma_series z n) =
+ (of_nat n) powr z / (z * (\<Prod>k=1..n. exp (Ln (z / of_nat k + 1))))"
+ unfolding ln_Gamma_series_def powr_def by (simp add: exp_diff exp_setsum)
+ also from assms have "(\<Prod>k=1..n. exp (Ln (z / of_nat k + 1))) = (\<Prod>k=1..n. z / of_nat k + 1)"
+ by (intro setprod.cong[OF refl], subst exp_Ln) (auto simp: field_simps plus_of_nat_eq_0_imp)
+ also have "... = (\<Prod>k=1..n. z + k) / fact n" unfolding fact_altdef
+ by (subst setprod_dividef [symmetric]) (simp_all add: field_simps)
+ also from assms have "z * ... = (\<Prod>k=0..n. z + k) / fact n"
+ by (cases n) (simp_all add: setprod_nat_ivl_1_Suc)
+ also have "(\<Prod>k=0..n. z + k) = pochhammer z (Suc n)" unfolding pochhammer_def by simp
+ 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: divide_simps powr_def Ln_of_nat)
+ finally show ?thesis .
+qed
+
+
+lemma ln_Gamma_series'_aux:
+ assumes "(z::complex) \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ shows "(\<lambda>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 "(\<lambda>n. ln_Gamma_series z n + of_real (harm n - ln (of_nat n)) * z + ln z) \<longlonglongrightarrow> ?s"
+ (is "?g \<longlonglongrightarrow> _")
+ by (intro tendsto_intros ln_Gamma_complex_LIMSEQ euler_mascheroni_LIMSEQ_of_real assms)
+
+ have A: "eventually (\<lambda>n. (\<Sum>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 "(\<Sum>k<n. ?f k) = (\<Sum>k=1..n. z / of_nat k - ln (1 + z / of_nat k))"
+ by (subst atLeast0LessThan [symmetric], subst setsum_shift_bounds_Suc_ivl [symmetric],
+ subst atLeastLessThanSuc_atLeastAtMost) simp_all
+ also have "\<dots> = z * of_real (harm n) - (\<Sum>k=1..n. ln (1 + z / of_nat k))"
+ by (simp add: harm_def setsum_subtractf setsum_right_distrib divide_inverse)
+ also from n have "\<dots> - ?g n = 0"
+ by (simp add: ln_Gamma_series_def setsum_subtractf algebra_simps Ln_of_nat)
+ finally show "(\<Sum>k<n. ?f k) - ?g n = 0" .
+ qed
+ show "(\<lambda>n. (\<Sum>k<n. ?f k) - ?g n) \<longlonglongrightarrow> 0" by (subst tendsto_cong[OF A]) simp_all
+qed
+
+
+lemma uniformly_summable_deriv_ln_Gamma:
+ assumes z: "(z :: 'a :: {real_normed_field,banach}) \<noteq> 0" and d: "d > 0" "d \<le> norm z/2"
+ shows "uniformly_convergent_on (ball z d)
+ (\<lambda>k z. \<Sum>i<k. inverse (of_nat (Suc i)) - inverse (z + of_nat (Suc i)))"
+ (is "uniformly_convergent_on _ (\<lambda>k z. \<Sum>i<k. ?f i z)")
+proof (rule weierstrass_m_test'_ev)
+ {
+ fix t assume t: "t \<in> ball z d"
+ have "norm z = norm (t + (z - t))" by simp
+ have "norm (t + (z - t)) \<le> 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 "\<dots> \<le> norm z + norm (t - z)" by (rule norm_triangle_ineq)
+ also from t d have "norm (t - z) \<le> norm z / 2" by (simp add: dist_norm norm_minus_commute)
+ also from z have "\<dots> < norm z" by simp
+ finally have B: "norm t < 2 * norm z" by simp
+ note A B
+ } note ball = this
+
+ show "eventually (\<lambda>n. \<forall>t\<in>ball z d. norm (?f n t) \<le> 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 \<in> ball z d"
+ from z ball[OF t] have t_nz: "t \<noteq> 0" by auto
+ fix n :: nat assume n: "n > nat \<lceil>4 * norm z\<rceil>"
+ from ball[OF t] t_nz have "4 * norm z > 2 * norm t" by simp
+ also from n have "\<dots> < 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 \<noteq> -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))"
+ by (simp add: divide_simps eq_neg_iff_add_eq_0 del: of_nat_Suc)
+ also have "norm \<dots> = inverse (of_nat (Suc n)) * inverse (norm (of_nat (Suc n)/t + 1))"
+ by (simp add: norm_divide norm_mult divide_simps add_ac del: of_nat_Suc)
+ also {
+ from z t_nz ball[OF t] have "of_nat (Suc n) / (4 * norm z) \<le> of_nat (Suc n) / (2 * norm t)"
+ by (intro divide_left_mono mult_pos_pos) simp_all
+ also have "\<dots> < 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 "\<dots> \<le> norm (of_nat (Suc n)/t + 1)" by (rule norm_diff_ineq)
+ finally have "inverse (norm (of_nat (Suc n)/t + 1)) \<le> 4 * norm z / of_nat (Suc n)"
+ using z by (simp add: divide_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: divide_simps power2_eq_square del: of_nat_Suc)
+ finally show "norm (?f n t) \<le> 4 * norm z * inverse (of_nat (Suc n)^2)"
+ by (simp del: of_nat_Suc)
+ qed
+next
+ show "summable (\<lambda>n. 4 * norm z * inverse ((of_nat (Suc n))^2))"
+ by (subst summable_Suc_iff) (simp add: summable_mult inverse_power_summable)
+qed
+
+lemma summable_deriv_ln_Gamma:
+ "z \<noteq> (0 :: 'a :: {real_normed_field,banach}) \<Longrightarrow>
+ summable (\<lambda>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 Polygamma :: "nat \<Rightarrow> ('a :: {real_normed_field,banach}) \<Rightarrow> 'a" where
+ "Polygamma n z = (if n = 0 then
+ (\<Sum>k. inverse (of_nat (Suc k)) - inverse (z + of_nat k)) - euler_mascheroni else
+ (-1)^Suc n * fact n * (\<Sum>k. inverse ((z + of_nat k)^Suc n)))"
+
+abbreviation Digamma :: "('a :: {real_normed_field,banach}) \<Rightarrow> 'a" where
+ "Digamma \<equiv> Polygamma 0"
+
+lemma Digamma_def:
+ "Digamma z = (\<Sum>k. inverse (of_nat (Suc k)) - inverse (z + of_nat k)) - euler_mascheroni"
+ by (simp add: Polygamma_def)
+
+
+lemma summable_Digamma:
+ assumes "(z :: 'a :: {real_normed_field,banach}) \<noteq> 0"
+ shows "summable (\<lambda>n. inverse (of_nat (Suc n)) - inverse (z + of_nat n))"
+proof -
+ have sums: "(\<lambda>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]
+ tendsto_add_filterlim_at_infinity[OF tendsto_const] tendsto_of_nat)
+ from summable_add[OF summable_deriv_ln_Gamma[OF assms] sums_summable[OF sums]]
+ show "summable (\<lambda>n. inverse (of_nat (Suc n)) - inverse (z + of_nat n))" by simp
+qed
+
+lemma summable_offset:
+ assumes "summable (\<lambda>n. f (n + k) :: 'a :: real_normed_vector)"
+ shows "summable f"
+proof -
+ from assms have "convergent (\<lambda>m. \<Sum>n<m. f (n + k))" by (simp add: summable_iff_convergent)
+ hence "convergent (\<lambda>m. (\<Sum>n<k. f n) + (\<Sum>n<m. f (n + k)))"
+ by (intro convergent_add convergent_const)
+ also have "(\<lambda>m. (\<Sum>n<k. f n) + (\<Sum>n<m. f (n + k))) = (\<lambda>m. \<Sum>n<m+k. f n)"
+ proof
+ fix m :: nat
+ have "{..<m+k} = {..<k} \<union> {k..<m+k}" by auto
+ also have "(\<Sum>n\<in>\<dots>. f n) = (\<Sum>n<k. f n) + (\<Sum>n=k..<m+k. f n)"
+ by (rule setsum.union_disjoint) auto
+ also have "(\<Sum>n=k..<m+k. f n) = (\<Sum>n=0..<m+k-k. f (n + k))"
+ by (intro setsum.reindex_cong[of "\<lambda>n. n + k"])
+ (simp, subst image_add_atLeastLessThan, auto)
+ finally show "(\<Sum>n<k. f n) + (\<Sum>n<m. f (n + k)) = (\<Sum>n<m+k. f n)" by (simp add: atLeast0LessThan)
+ qed
+ finally have "(\<lambda>a. setsum f {..<a}) \<longlonglongrightarrow> lim (\<lambda>m. setsum 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 \<noteq> 0" and n: "n \<ge> 2"
+ shows "uniformly_convergent_on (ball z d) (\<lambda>k z. \<Sum>i<k. inverse ((z + of_nat i)^n))"
+proof (rule weierstrass_m_test'_ev)
+ def e \<equiv> "(1 + d / norm z)"
+ def m \<equiv> "nat \<lceil>norm z * e\<rceil>"
+ {
+ fix t assume t: "t \<in> ball z d"
+ have "norm t = norm (z + (t - z))" by simp
+ also have "\<dots> \<le> 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 (\<lambda>k. \<forall>t\<in>ball z d. norm (inverse ((t + of_nat k)^n)) \<le>
+ 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 \<in> 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 "\<dots> \<le> norm (of_nat k :: 'a) - norm z * e"
+ unfolding m_def by (subst norm_of_nat) linarith
+ also from ball[OF t] have "\<dots> \<le> norm (of_nat k :: 'a) - norm t" by simp
+ also have "\<dots> \<le> norm (of_nat k + t)" by (rule norm_diff_ineq)
+ finally have "inverse ((norm (t + of_nat k))^n) \<le> inverse (real_of_nat (k - m)^n)" using k n
+ by (intro le_imp_inverse_le power_mono) (simp_all add: add_ac del: of_nat_Suc)
+ thus "norm (inverse ((t + of_nat k)^n)) \<le> inverse (of_nat (k - m)^n)"
+ by (simp add: norm_inverse norm_power power_inverse)
+ qed
+
+ have "summable (\<lambda>k. inverse ((real_of_nat k)^n))"
+ using inverse_power_summable[of n] n by simp
+ hence "summable (\<lambda>k. inverse ((real_of_nat (k + m - m))^n))" by simp
+ thus "summable (\<lambda>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 \<noteq> 0" and n: "n \<ge> 2"
+ shows "summable (\<lambda>k. inverse ((z + of_nat k)^n))"
+ using uniformly_convergent_imp_convergent[OF Polygamma_converges[OF assms, of 1], of z]
+ by (simp add: summable_iff_convergent)
+
+lemma has_field_derivative_ln_Gamma_complex [derivative_intros]:
+ assumes z: "Re z > 0 \<or> Im z \<noteq> 0"
+ shows "(ln_Gamma has_field_derivative Digamma z) (at z)"
+proof -
+ have not_nonpos_Int [simp]: "t \<notin> \<int>\<^sub>\<le>\<^sub>0" if "Re t > 0" for t
+ using that by (auto elim!: nonpos_Ints_cases')
+ from z have z': "z \<notin> \<int>\<^sub>\<le>\<^sub>0" and z'': "z \<noteq> 0" by (auto elim!: nonpos_Ints_cases)
+ let ?f' = "\<lambda>z k. inverse (of_nat (Suc k)) - inverse (z + of_nat (Suc k))"
+ let ?f = "\<lambda>z k. z / of_nat (Suc k) - ln (1 + z / of_nat (Suc k))" and ?F' = "\<lambda>z. \<Sum>n. ?f' z n"
+ def d \<equiv> "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 \<ge> d" by (auto simp add: d_def)
+ have ball: "Im t = 0 \<longrightarrow> Re t > 0" if "dist z t < d" for t
+ using Re_pos_or_Im_nz_in_ball[OF z, of t] that unfolding d_def by auto
+ have sums: "(\<lambda>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]
+ tendsto_add_filterlim_at_infinity[OF tendsto_const] tendsto_of_nat)
+
+ have "((\<lambda>z. \<Sum>n. ?f z n) has_field_derivative ?F' z) (at z)"
+ using d z ln_Gamma_series'_aux[OF z']
+ by (intro has_field_derivative_series'(2)[of "ball z d" _ _ z] uniformly_summable_deriv_ln_Gamma)
+ (auto intro!: derivative_eq_intros add_pos_pos mult_pos_pos dest!: ball
+ simp: field_simps Re_complex_div_gt_0 sums_iff Im_divide_of_nat
+ simp del: of_nat_Suc)
+ with z have "((\<lambda>z. (\<Sum>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 = (\<Sum>n. inverse (z + of_nat (Suc n)) - inverse (z + of_nat n))"
+ by (simp add: sums_iff)
+ also from sums summable_deriv_ln_Gamma[OF z'']
+ have "?F' z + \<dots> = (\<Sum>n. inverse (of_nat (Suc n)) - inverse (z + of_nat n))"
+ by (subst suminf_add) (simp_all add: add_ac sums_iff)
+ also have "\<dots> - euler_mascheroni = Digamma z" by (simp add: Digamma_def)
+ finally have "((\<lambda>z. (\<Sum>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 (\<lambda>z. ln_Gamma z = (\<Sum>k. ?f z k) - euler_mascheroni * z - Ln z) (nhds z)"
+ proof eventually_elim
+ fix t assume "t \<in> ball z d"
+ hence "t \<notin> \<int>\<^sub>\<le>\<^sub>0" by (auto dest!: ball elim!: nonpos_Ints_cases)
+ from ln_Gamma_series'_aux[OF this]
+ show "ln_Gamma t = (\<Sum>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"
+ by (simp add: Digamma_def)
+
+lemma Digamma_plus1:
+ assumes "z \<noteq> 0"
+ shows "Digamma (z+1) = Digamma z + 1/z"
+proof -
+ have sums: "(\<lambda>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]]
+ tendsto_add_filterlim_at_infinity[OF tendsto_const] tendsto_of_nat)
+ have "Digamma (z+1) = (\<Sum>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 = (\<Sum>k. inverse (of_nat (Suc k)) - inverse (z + of_nat k)) +
+ (\<Sum>k. inverse (z + of_nat k) - inverse (z + of_nat (Suc k)))"
+ using summable_Digamma[OF assms] sums by (subst suminf_add) (simp_all add: add_ac sums_iff)
+ also have "(\<Sum>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
+
+lemma Polygamma_plus1:
+ assumes "z \<noteq> 0"
+ shows "Polygamma n (z + 1) = Polygamma n z + (-1)^n * fact n / (z ^ Suc n)"
+proof (cases "n = 0")
+ assume n: "n \<noteq> 0"
+ let ?f = "\<lambda>k. inverse ((z + of_nat k) ^ Suc n)"
+ have "Polygamma n (z + 1) = (-1) ^ Suc n * fact n * (\<Sum>k. ?f (k+1))"
+ using n by (simp add: Polygamma_def add_ac)
+ also have "(\<Sum>k. ?f (k+1)) + (\<Sum>k<1. ?f k) = (\<Sum>k. ?f k)"
+ using Polygamma_converges'[OF assms, of "Suc n"] n
+ by (subst suminf_split_initial_segment [symmetric]) simp_all
+ hence "(\<Sum>k. ?f (k+1)) = (\<Sum>k. ?f k) - inverse (z ^ Suc n)" by (simp add: algebra_simps)
+ also have "(-1) ^ Suc n * fact n * ((\<Sum>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)
+
+lemma 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 "\<dots> = 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 "\<dots> + inverse (of_nat (Suc n)) = harm (Suc n) - euler_mascheroni"
+ by (simp add: harm_Suc)
+ finally show ?case .
+qed (simp add: harm_def)
+
+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 \<noteq> 0 \<Longrightarrow> Polygamma n (of_real x) = of_real (Polygamma n x)"
+ unfolding Polygamma_def using summable_Digamma[of x] Polygamma_converges'[of x "Suc n"]
+ by (simp_all add: suminf_of_real)
+
+lemma Polygamma_Real: "z \<in> \<real> \<Longrightarrow> z \<noteq> 0 \<Longrightarrow> Polygamma n z \<in> \<real>"
+ 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}) =
+ (\<Sum>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) =
+ (\<Sum>k. inverse (of_nat (Suc k)) - inverse (of_nat k + 1/2)) - euler_mascheroni"
+ by (simp add: Digamma_def add_ac)
+ also have "(\<Sum>k. inverse (of_nat (Suc k) :: real) - inverse (of_nat k + 1/2)) =
+ (\<Sum>k. inverse (1/2) * (inverse (2 * of_nat (Suc k)) - inverse (2 * of_nat k + 1)))"
+ by (simp_all add: add_ac inverse_mult_distrib[symmetric] ring_distribs del: inverse_divide)
+ also have "\<dots> = - 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) \<noteq> 0"
+ using of_nat_neq_0[of "2*n"] by (simp only: of_nat_Suc) (simp add: add_ac)
+ hence nz': "of_nat n + (1/2::'a) \<noteq> 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 "\<dots> = 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
+ finally show ?case by (simp add: add_ac)
+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 "\<dots> = 2 * ln 2 + euler_mascheroni - 2" by (subst Digamma_half_integer) simp
+ also from euler_mascheroni_approx have "euler_mascheroni \<le> (0.58::real)"
+ by (simp add: abs_real_def split: split_if_asm)
+ also from ln_2_bounds have "ln 2 < (0.7 :: real)" by simp
+ finally show ?thesis by simp
+qed
+
+
+lemma has_field_derivative_Polygamma [derivative_intros]:
+ fixes z :: "'a :: {real_normed_field,euclidean_space}"
+ assumes z: "z \<notin> \<int>\<^sub>\<le>\<^sub>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 = "\<lambda>k z. inverse (of_nat (Suc k)) - inverse (z + of_nat k)"
+ let ?F = "\<lambda>z. \<Sum>k. ?f k z" and ?f' = "\<lambda>k z. inverse ((z + of_nat k)\<^sup>2)"
+ from no_nonpos_Int_in_ball'[OF z] guess d . note d = this
+ from z have summable: "summable (\<lambda>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) (\<lambda>k z. \<Sum>i<k. inverse ((z + of_nat i)\<^sup>2))"
+ by (intro Polygamma_converges) auto
+ with d have "summable (\<lambda>k. inverse ((z + of_nat k)\<^sup>2))" unfolding summable_iff_convergent
+ by (auto dest!: uniformly_convergent_imp_convergent simp: summable_iff_convergent )
+
+ have "(?F has_field_derivative (\<Sum>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 \<in> ball z d"
+ from t d(2)[of t] show "((\<lambda>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 \<noteq> 0"
+ from z have z': "z \<noteq> 0" by auto
+ from no_nonpos_Int_in_ball'[OF z] guess d . note d = this
+ def n' \<equiv> "Suc n"
+ from n have n': "n' \<ge> 2" by (simp add: n'_def)
+ have "((\<lambda>z. \<Sum>k. inverse ((z + of_nat k) ^ n')) has_field_derivative
+ (\<Sum>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 \<in> ball z d"
+ with d have t': "t \<notin> \<int>\<^sub>\<le>\<^sub>0" "t \<noteq> 0" by auto
+ show "((\<lambda>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)
+ (\<lambda>k z. (- of_nat n' :: 'a) * (\<Sum>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)
+ (\<lambda>k z. \<Sum>i<k. - of_nat n' * inverse ((z + of_nat i :: 'a) ^ (n'+1)))"
+ by (subst (asm) setsum_right_distrib) simp
+ qed (insert Polygamma_converges'[OF z' n'] d, simp_all)
+ also have "(\<Sum>k. - of_nat n' * inverse ((z + of_nat k) ^ (n' + 1))) =
+ (- of_nat n') * (\<Sum>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 "((\<lambda>z. \<Sum>k. inverse ((z + of_nat k) ^ n')) has_field_derivative
+ - of_nat n' * (\<Sum>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 :: "_ \<Rightarrow> 'a :: {real_normed_field,euclidean_space}"
+ shows "isCont f z \<Longrightarrow> f z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> isCont (\<lambda>x. Polygamma n (f x)) z"
+ by (rule isCont_o2[OF _ DERIV_isCont[OF has_field_derivative_Polygamma]])
+
+lemma continuous_on_Polygamma:
+ "A \<inter> \<int>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> continuous_on A (Polygamma n :: _ \<Rightarrow> '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]:
+ "isCont f z \<Longrightarrow> Re (f z) > 0 \<or> Im (f z) \<noteq> 0 \<Longrightarrow> isCont (\<lambda>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]:
+ "A \<inter> {z. Re z \<le> 0 \<and> Im z = 0} = {} \<Longrightarrow> continuous_on A ln_Gamma"
+ by (intro continuous_at_imp_continuous_on ballI isCont_ln_Gamma_complex[OF continuous_ident])
+ fastforce
+
+
+text \<open>
+ We define a type class that captures all the fundamental properties of the inverse 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.
+\<close>
+
+class Gamma = real_normed_field + complete_space +
+ fixes rGamma :: "'a \<Rightarrow> 'a"
+ assumes rGamma_eq_zero_iff_aux: "rGamma z = 0 \<longleftrightarrow> (\<exists>n. z = - of_nat n)"
+ assumes differentiable_rGamma_aux1:
+ "(\<And>n. z \<noteq> - of_nat n) \<Longrightarrow>
+ let d = (THE d. (\<lambda>n. \<Sum>k<n. inverse (of_nat (Suc k)) - inverse (z + of_nat k))
+ \<longlonglongrightarrow> d) - scaleR euler_mascheroni 1
+ in filterlim (\<lambda>y. (rGamma y - rGamma z + rGamma z * d * (y - z)) /\<^sub>R
+ norm (y - z)) (nhds 0) (at z)"
+ assumes differentiable_rGamma_aux2:
+ "let z = - of_nat n
+ in filterlim (\<lambda>y. (rGamma y - rGamma z - (-1)^n * (setprod of_nat {1..n}) * (y - z)) /\<^sub>R
+ norm (y - z)) (nhds 0) (at z)"
+ assumes rGamma_series_aux: "(\<And>n. z \<noteq> - of_nat n) \<Longrightarrow>
+ let fact' = (\<lambda>n. setprod of_nat {1..n});
+ exp = (\<lambda>x. THE e. (\<lambda>n. \<Sum>k<n. x^k /\<^sub>R fact k) \<longlonglongrightarrow> e);
+ pochhammer' = (\<lambda>a n. (\<Prod>n = 0..n. a + of_nat n))
+ in filterlim (\<lambda>n. pochhammer' z n / (fact' n * exp (z * (ln (of_nat n) *\<^sub>R 1))))
+ (nhds (rGamma z)) sequentially"
+begin
+subclass banach ..
+end
+
+definition "Gamma z = inverse (rGamma z)"
+
+
+subsection \<open>Basic properties\<close>
+
+lemma Gamma_nonpos_Int: "z \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Gamma z = 0"
+ and rGamma_nonpos_Int: "z \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> rGamma z = 0"
+ using rGamma_eq_zero_iff_aux[of z] unfolding Gamma_def by (auto elim!: nonpos_Ints_cases')
+
+lemma Gamma_nonzero: "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Gamma z \<noteq> 0"
+ and rGamma_nonzero: "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> rGamma z \<noteq> 0"
+ using rGamma_eq_zero_iff_aux[of z] unfolding Gamma_def by (auto elim!: nonpos_Ints_cases')
+
+lemma Gamma_eq_zero_iff: "Gamma z = 0 \<longleftrightarrow> z \<in> \<int>\<^sub>\<le>\<^sub>0"
+ and rGamma_eq_zero_iff: "rGamma z = 0 \<longleftrightarrow> z \<in> \<int>\<^sub>\<le>\<^sub>0"
+ using rGamma_eq_zero_iff_aux[of z] unfolding Gamma_def by (auto elim!: nonpos_Ints_cases')
+
+lemma rGamma_inverse_Gamma: "rGamma z = inverse (Gamma z)"
+ unfolding Gamma_def by simp
+
+lemma rGamma_series_LIMSEQ [tendsto_intros]:
+ "rGamma_series z \<longlonglongrightarrow> rGamma z"
+proof (cases "z \<in> \<int>\<^sub>\<le>\<^sub>0")
+ case False
+ hence "z \<noteq> - of_nat n" for n by auto
+ from rGamma_series_aux[OF this] show ?thesis
+ by (simp add: rGamma_series_def[abs_def] fact_altdef pochhammer_Suc_setprod
+ exp_def of_real_def[symmetric] suminf_def sums_def[abs_def])
+qed (insert rGamma_eq_zero_iff[of z], simp_all add: rGamma_series_nonpos_Ints_LIMSEQ)
+
+lemma Gamma_series_LIMSEQ [tendsto_intros]:
+ "Gamma_series z \<longlonglongrightarrow> Gamma z"
+proof (cases "z \<in> \<int>\<^sub>\<le>\<^sub>0")
+ case False
+ hence "(\<lambda>n. inverse (rGamma_series z n)) \<longlonglongrightarrow> inverse (rGamma z)"
+ by (intro tendsto_intros) (simp_all add: rGamma_eq_zero_iff)
+ also have "(\<lambda>n. inverse (rGamma_series z n)) = Gamma_series z"
+ by (simp add: rGamma_series_def Gamma_series_def[abs_def])
+ finally show ?thesis by (simp add: Gamma_def)
+qed (insert Gamma_eq_zero_iff[of z], simp_all add: Gamma_series_nonpos_Ints_LIMSEQ)
+
+lemma Gamma_altdef: "Gamma z = lim (Gamma_series z)"
+ using Gamma_series_LIMSEQ[of z] by (simp add: limI)
+
+lemma rGamma_1 [simp]: "rGamma 1 = 1"
+proof -
+ have A: "eventually (\<lambda>n. rGamma_series 1 n = of_nat (Suc n) / of_nat n) sequentially"
+ using eventually_gt_at_top[of "0::nat"]
+ by (force elim!: eventually_mono simp: rGamma_series_def exp_of_real pochhammer_fact
+ divide_simps pochhammer_rec' dest!: pochhammer_eq_0_imp_nonpos_Int)
+ have "rGamma_series 1 \<longlonglongrightarrow> 1" by (subst tendsto_cong[OF A]) (rule LIMSEQ_Suc_n_over_n)
+ moreover have "rGamma_series 1 \<longlonglongrightarrow> rGamma 1" by (rule tendsto_intros)
+ ultimately show ?thesis by (intro LIMSEQ_unique)
+qed
+
+lemma rGamma_plus1: "z * rGamma (z + 1) = rGamma z"
+proof -
+ let ?f = "\<lambda>n. (z + 1) * inverse (of_nat n) + 1"
+ have "eventually (\<lambda>n. ?f n * rGamma_series z n = z * rGamma_series (z + 1) n) sequentially"
+ using eventually_gt_at_top[of "0::nat"]
+ proof eventually_elim
+ fix n :: nat assume n: "n > 0"
+ hence "z * rGamma_series (z + 1) n = inverse (of_nat n) *
+ pochhammer z (Suc (Suc n)) / (fact n * exp (z * of_real (ln (of_nat n))))"
+ by (subst pochhammer_rec) (simp add: rGamma_series_def field_simps exp_add exp_of_real)
+ also from n have "\<dots> = ?f n * rGamma_series z n"
+ by (subst pochhammer_rec') (simp_all add: divide_simps rGamma_series_def add_ac)
+ finally show "?f n * rGamma_series z n = z * rGamma_series (z + 1) n" ..
+ qed
+ moreover have "(\<lambda>n. ?f n * rGamma_series z n) \<longlonglongrightarrow> ((z+1) * 0 + 1) * rGamma z"
+ by (intro tendsto_intros lim_inverse_n)
+ hence "(\<lambda>n. ?f n * rGamma_series z n) \<longlonglongrightarrow> rGamma z" by simp
+ ultimately have "(\<lambda>n. z * rGamma_series (z + 1) n) \<longlonglongrightarrow> rGamma z"
+ by (rule Lim_transform_eventually)
+ moreover have "(\<lambda>n. z * rGamma_series (z + 1) n) \<longlonglongrightarrow> z * rGamma (z + 1)"
+ by (intro tendsto_intros)
+ ultimately show "z * rGamma (z + 1) = rGamma z" using LIMSEQ_unique by blast
+qed
+
+
+lemma pochhammer_rGamma: "rGamma z = pochhammer z n * rGamma (z + of_nat n)"
+proof (induction n arbitrary: z)
+ case (Suc n z)
+ have "rGamma z = pochhammer z n * rGamma (z + of_nat n)" by (rule Suc.IH)
+ also note rGamma_plus1 [symmetric]
+ finally show ?case by (simp add: add_ac pochhammer_rec')
+qed simp_all
+
+lemma Gamma_plus1: "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Gamma (z + 1) = z * Gamma z"
+ using rGamma_plus1[of z] by (simp add: rGamma_inverse_Gamma field_simps Gamma_eq_zero_iff)
+
+lemma pochhammer_Gamma: "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> pochhammer z n = Gamma (z + of_nat n) / Gamma z"
+ using pochhammer_rGamma[of z]
+ by (simp add: rGamma_inverse_Gamma Gamma_eq_zero_iff field_simps)
+
+lemma Gamma_0 [simp]: "Gamma 0 = 0"
+ and rGamma_0 [simp]: "rGamma 0 = 0"
+ and Gamma_neg_1 [simp]: "Gamma (- 1) = 0"
+ and rGamma_neg_1 [simp]: "rGamma (- 1) = 0"
+ and Gamma_neg_numeral [simp]: "Gamma (- numeral n) = 0"
+ and rGamma_neg_numeral [simp]: "rGamma (- numeral n) = 0"
+ and Gamma_neg_of_nat [simp]: "Gamma (- of_nat m) = 0"
+ and rGamma_neg_of_nat [simp]: "rGamma (- of_nat m) = 0"
+ by (simp_all add: rGamma_eq_zero_iff Gamma_eq_zero_iff)
+
+lemma Gamma_1 [simp]: "Gamma 1 = 1" unfolding Gamma_def by simp
+
+lemma Gamma_fact: "Gamma (of_nat (Suc n)) = fact n"
+ by (simp add: pochhammer_fact pochhammer_Gamma of_nat_in_nonpos_Ints_iff)
+
+lemma Gamma_numeral: "Gamma (numeral n) = fact (pred_numeral n)"
+ by (subst of_nat_numeral[symmetric], subst numeral_eq_Suc, subst Gamma_fact) (rule refl)
+
+lemma Gamma_of_int: "Gamma (of_int n) = (if n > 0 then fact (nat (n - 1)) else 0)"
+proof (cases "n > 0")
+ case True
+ hence "Gamma (of_int n) = Gamma (of_nat (Suc (nat (n - 1))))" by (subst of_nat_Suc) simp_all
+ with True show ?thesis by (subst (asm) Gamma_fact) simp
+qed (simp_all add: Gamma_eq_zero_iff nonpos_Ints_of_int)
+
+lemma rGamma_of_int: "rGamma (of_int n) = (if n > 0 then inverse (fact (nat (n - 1))) else 0)"
+ by (simp add: Gamma_of_int rGamma_inverse_Gamma)
+
+lemma Gamma_seriesI:
+ assumes "(\<lambda>n. g n / Gamma_series z n) \<longlonglongrightarrow> 1"
+ shows "g \<longlonglongrightarrow> Gamma z"
+proof (rule Lim_transform_eventually)
+ have "1/2 > (0::real)" by simp
+ from tendstoD[OF assms, OF this]
+ show "eventually (\<lambda>n. g n / Gamma_series z n * Gamma_series z n = g n) sequentially"
+ by (force elim!: eventually_mono simp: dist_real_def dist_0_norm)
+ from assms have "(\<lambda>n. g n / Gamma_series z n * Gamma_series z n) \<longlonglongrightarrow> 1 * Gamma z"
+ by (intro tendsto_intros)
+ thus "(\<lambda>n. g n / Gamma_series z n * Gamma_series z n) \<longlonglongrightarrow> Gamma z" by simp
+qed
+
+lemma Gamma_seriesI':
+ assumes "f \<longlonglongrightarrow> rGamma z"
+ assumes "(\<lambda>n. g n * f n) \<longlonglongrightarrow> 1"
+ assumes "z \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ shows "g \<longlonglongrightarrow> Gamma z"
+proof (rule Lim_transform_eventually)
+ have "1/2 > (0::real)" by simp
+ from tendstoD[OF assms(2), OF this] show "eventually (\<lambda>n. g n * f n / f n = g n) sequentially"
+ by (force elim!: eventually_mono simp: dist_real_def dist_0_norm)
+ from assms have "(\<lambda>n. g n * f n / f n) \<longlonglongrightarrow> 1 / rGamma z"
+ by (intro tendsto_divide assms) (simp_all add: rGamma_eq_zero_iff)
+ thus "(\<lambda>n. g n * f n / f n) \<longlonglongrightarrow> Gamma z" by (simp add: Gamma_def divide_inverse)
+qed
+
+lemma Gamma_series'_LIMSEQ: "Gamma_series' z \<longlonglongrightarrow> Gamma z"
+ by (cases "z \<in> \<int>\<^sub>\<le>\<^sub>0") (simp_all add: Gamma_nonpos_Int Gamma_seriesI[OF Gamma_series_Gamma_series']
+ Gamma_series'_nonpos_Ints_LIMSEQ[of z])
+
+
+subsection \<open>Differentiability\<close>
+
+lemma has_field_derivative_rGamma_no_nonpos_int:
+ assumes "z \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ shows "(rGamma has_field_derivative -rGamma z * Digamma z) (at z within A)"
+proof (rule has_field_derivative_at_within)
+ from assms have "z \<noteq> - of_nat n" for n by auto
+ from differentiable_rGamma_aux1[OF this]
+ show "(rGamma has_field_derivative -rGamma z * Digamma z) (at z)"
+ unfolding Digamma_def suminf_def sums_def[abs_def]
+ has_field_derivative_def has_derivative_def netlimit_at
+ by (simp add: Let_def bounded_linear_mult_right mult_ac of_real_def [symmetric])
+qed
+
+lemma has_field_derivative_rGamma_nonpos_int:
+ "(rGamma has_field_derivative (-1)^n * fact n) (at (- of_nat n) within A)"
+ apply (rule has_field_derivative_at_within)
+ using differentiable_rGamma_aux2[of n]
+ unfolding Let_def has_field_derivative_def has_derivative_def netlimit_at
+ by (simp only: bounded_linear_mult_right mult_ac of_real_def [symmetric] fact_altdef)
+
+lemma has_field_derivative_rGamma [derivative_intros]:
+ "(rGamma has_field_derivative (if z \<in> \<int>\<^sub>\<le>\<^sub>0 then (-1)^(nat \<lfloor>norm z\<rfloor>) * fact (nat \<lfloor>norm z\<rfloor>)
+ else -rGamma z * Digamma z)) (at z within A)"
+using has_field_derivative_rGamma_no_nonpos_int[of z A]
+ has_field_derivative_rGamma_nonpos_int[of "nat \<lfloor>norm z\<rfloor>" A]
+ by (auto elim!: nonpos_Ints_cases')
+
+declare has_field_derivative_rGamma_no_nonpos_int [THEN DERIV_chain2, derivative_intros]
+declare has_field_derivative_rGamma [THEN DERIV_chain2, derivative_intros]
+declare has_field_derivative_rGamma_nonpos_int [derivative_intros]
+declare has_field_derivative_rGamma_no_nonpos_int [derivative_intros]
+declare has_field_derivative_rGamma [derivative_intros]
+
+
+lemma has_field_derivative_Gamma [derivative_intros]:
+ "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> (Gamma has_field_derivative Gamma z * Digamma z) (at z within A)"
+ unfolding Gamma_def [abs_def]
+ by (fastforce intro!: derivative_eq_intros simp: rGamma_eq_zero_iff)
+
+declare has_field_derivative_Gamma[THEN DERIV_chain2, derivative_intros]
+
+(* TODO: Hide ugly facts properly *)
+hide_fact rGamma_eq_zero_iff_aux differentiable_rGamma_aux1 differentiable_rGamma_aux2
+ differentiable_rGamma_aux2 rGamma_series_aux Gamma_class.rGamma_eq_zero_iff_aux
+
+
+
+(* TODO: differentiable etc. *)
+
+
+subsection \<open>Continuity\<close>
+
+lemma continuous_on_rGamma [continuous_intros]: "continuous_on A rGamma"
+ by (rule DERIV_continuous_on has_field_derivative_rGamma)+
+
+lemma continuous_on_Gamma [continuous_intros]: "A \<inter> \<int>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> continuous_on A Gamma"
+ by (rule DERIV_continuous_on has_field_derivative_Gamma)+ blast
+
+lemma isCont_rGamma [continuous_intros]:
+ "isCont f z \<Longrightarrow> isCont (\<lambda>x. rGamma (f x)) z"
+ by (rule isCont_o2[OF _ DERIV_isCont[OF has_field_derivative_rGamma]])
+
+lemma isCont_Gamma [continuous_intros]:
+ "isCont f z \<Longrightarrow> f z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> isCont (\<lambda>x. Gamma (f x)) z"
+ by (rule isCont_o2[OF _ DERIV_isCont[OF has_field_derivative_Gamma]])
+
+
+
+text \<open>The complex Gamma function\<close>
+
+instantiation complex :: Gamma
+begin
+
+definition rGamma_complex :: "complex \<Rightarrow> complex" where
+ "rGamma_complex z = lim (rGamma_series z)"
+
+lemma rGamma_series_complex_converges:
+ "convergent (rGamma_series (z :: complex))" (is "?thesis1")
+ and rGamma_complex_altdef:
+ "rGamma z = (if z \<in> \<int>\<^sub>\<le>\<^sub>0 then 0 else exp (-ln_Gamma z))" (is "?thesis2")
+proof -
+ have "?thesis1 \<and> ?thesis2"
+ proof (cases "z \<in> \<int>\<^sub>\<le>\<^sub>0")
+ case False
+ have "rGamma_series z \<longlonglongrightarrow> exp (- ln_Gamma z)"
+ proof (rule Lim_transform_eventually)
+ from ln_Gamma_series_complex_converges'[OF False] guess d by (elim exE conjE)
+ from this(1) uniformly_convergent_imp_convergent[OF this(2), of z]
+ have "ln_Gamma_series z \<longlonglongrightarrow> lim (ln_Gamma_series z)" by (simp add: convergent_LIMSEQ_iff)
+ thus "(\<lambda>n. exp (-ln_Gamma_series z n)) \<longlonglongrightarrow> exp (- ln_Gamma z)"
+ unfolding convergent_def ln_Gamma_def by (intro tendsto_exp tendsto_minus)
+ from eventually_gt_at_top[of "0::nat"] exp_ln_Gamma_series_complex False
+ show "eventually (\<lambda>n. exp (-ln_Gamma_series z n) = rGamma_series z n) sequentially"
+ by (force elim!: eventually_mono simp: exp_minus Gamma_series_def rGamma_series_def)
+ qed
+ with False show ?thesis
+ by (auto simp: convergent_def rGamma_complex_def intro!: limI)
+ next
+ case True
+ then obtain k where "z = - of_nat k" by (erule nonpos_Ints_cases')
+ also have "rGamma_series \<dots> \<longlonglongrightarrow> 0"
+ by (subst tendsto_cong[OF rGamma_series_minus_of_nat]) (simp_all add: convergent_const)
+ finally show ?thesis using True
+ by (auto simp: rGamma_complex_def convergent_def intro!: limI)
+ qed
+ thus "?thesis1" "?thesis2" by blast+
+qed
+
+context
+begin
+
+(* TODO: duplication *)
+private lemma rGamma_complex_plus1: "z * rGamma (z + 1) = rGamma (z :: complex)"
+proof -
+ let ?f = "\<lambda>n. (z + 1) * inverse (of_nat n) + 1"
+ have "eventually (\<lambda>n. ?f n * rGamma_series z n = z * rGamma_series (z + 1) n) sequentially"
+ using eventually_gt_at_top[of "0::nat"]
+ proof eventually_elim
+ fix n :: nat assume n: "n > 0"
+ hence "z * rGamma_series (z + 1) n = inverse (of_nat n) *
+ pochhammer z (Suc (Suc n)) / (fact n * exp (z * of_real (ln (of_nat n))))"
+ by (subst pochhammer_rec) (simp add: rGamma_series_def field_simps exp_add exp_of_real)
+ also from n have "\<dots> = ?f n * rGamma_series z n"
+ by (subst pochhammer_rec') (simp_all add: divide_simps rGamma_series_def add_ac)
+ finally show "?f n * rGamma_series z n = z * rGamma_series (z + 1) n" ..
+ qed
+ moreover have "(\<lambda>n. ?f n * rGamma_series z n) \<longlonglongrightarrow> ((z+1) * 0 + 1) * rGamma z"
+ using rGamma_series_complex_converges
+ by (intro tendsto_intros lim_inverse_n)
+ (simp_all add: convergent_LIMSEQ_iff rGamma_complex_def)
+ hence "(\<lambda>n. ?f n * rGamma_series z n) \<longlonglongrightarrow> rGamma z" by simp
+ ultimately have "(\<lambda>n. z * rGamma_series (z + 1) n) \<longlonglongrightarrow> rGamma z"
+ by (rule Lim_transform_eventually)
+ moreover have "(\<lambda>n. z * rGamma_series (z + 1) n) \<longlonglongrightarrow> z * rGamma (z + 1)"
+ using rGamma_series_complex_converges
+ by (auto intro!: tendsto_mult simp: rGamma_complex_def convergent_LIMSEQ_iff)
+ ultimately show "z * rGamma (z + 1) = rGamma z" using LIMSEQ_unique by blast
+qed
+
+private lemma has_field_derivative_rGamma_complex_no_nonpos_Int:
+ assumes "(z :: complex) \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ shows "(rGamma has_field_derivative - rGamma z * Digamma z) (at z)"
+proof -
+ have diff: "(rGamma has_field_derivative - rGamma z * Digamma z) (at z)" if "Re z > 0" for z
+ proof (subst DERIV_cong_ev[OF refl _ refl])
+ from that have "eventually (\<lambda>t. t \<in> ball z (Re z/2)) (nhds z)"
+ by (intro eventually_nhds_in_nhd) simp_all
+ thus "eventually (\<lambda>t. rGamma t = exp (- ln_Gamma t)) (nhds z)"
+ using no_nonpos_Int_in_ball_complex[OF that]
+ by (auto elim!: eventually_mono simp: rGamma_complex_altdef)
+ next
+ show "((\<lambda>t. exp (- ln_Gamma t)) has_field_derivative (-rGamma z * Digamma z)) (at z)"
+ using that by (force elim!: nonpos_Ints_cases intro!: derivative_eq_intros
+ simp: rGamma_complex_altdef)
+ qed
+
+ from assms show "(rGamma has_field_derivative - rGamma z * Digamma z) (at z)"
+ proof (induction "nat \<lfloor>1 - Re z\<rfloor>" arbitrary: z)
+ case (Suc n z)
+ from Suc.prems have z: "z \<noteq> 0" by auto
+ from Suc.hyps have "n = nat \<lfloor>- Re z\<rfloor>" by linarith
+ hence A: "n = nat \<lfloor>1 - Re (z + 1)\<rfloor>" by simp
+ from Suc.prems have B: "z + 1 \<notin> \<int>\<^sub>\<le>\<^sub>0" by (force dest: plus_one_in_nonpos_Ints_imp)
+
+ have "((\<lambda>z. z * (rGamma \<circ> (\<lambda>z. z + 1)) z) has_field_derivative
+ -rGamma (z + 1) * (Digamma (z + 1) * z - 1)) (at z)"
+ by (rule derivative_eq_intros DERIV_chain Suc refl A B)+ (simp add: algebra_simps)
+ also have "(\<lambda>z. z * (rGamma \<circ> (\<lambda>z. z + 1 :: complex)) z) = rGamma"
+ by (simp add: rGamma_complex_plus1)
+ also from z have "Digamma (z + 1) * z - 1 = z * Digamma z"
+ by (subst Digamma_plus1) (simp_all add: field_simps)
+ also have "-rGamma (z + 1) * (z * Digamma z) = -rGamma z * Digamma z"
+ by (simp add: rGamma_complex_plus1[of z, symmetric])
+ finally show ?case .
+ qed (intro diff, simp)
+qed
+
+private lemma rGamma_complex_1: "rGamma (1 :: complex) = 1"
+proof -
+ have A: "eventually (\<lambda>n. rGamma_series 1 n = of_nat (Suc n) / of_nat n) sequentially"
+ using eventually_gt_at_top[of "0::nat"]
+ by (force elim!: eventually_mono simp: rGamma_series_def exp_of_real pochhammer_fact
+ divide_simps pochhammer_rec' dest!: pochhammer_eq_0_imp_nonpos_Int)
+ have "rGamma_series 1 \<longlonglongrightarrow> 1" by (subst tendsto_cong[OF A]) (rule LIMSEQ_Suc_n_over_n)
+ thus "rGamma 1 = (1 :: complex)" unfolding rGamma_complex_def by (rule limI)
+qed
+
+private lemma has_field_derivative_rGamma_complex_nonpos_Int:
+ "(rGamma has_field_derivative (-1)^n * fact n) (at (- of_nat n :: complex))"
+proof (induction n)
+ case 0
+ have A: "(0::complex) + 1 \<notin> \<int>\<^sub>\<le>\<^sub>0" by simp
+ have "((\<lambda>z. z * (rGamma \<circ> (\<lambda>z. z + 1 :: complex)) z) has_field_derivative 1) (at 0)"
+ by (rule derivative_eq_intros DERIV_chain refl
+ has_field_derivative_rGamma_complex_no_nonpos_Int A)+ (simp add: rGamma_complex_1)
+ thus ?case by (simp add: rGamma_complex_plus1)
+next
+ case (Suc n)
+ hence A: "(rGamma has_field_derivative (-1)^n * fact n)
+ (at (- of_nat (Suc n) + 1 :: complex))" by simp
+ have "((\<lambda>z. z * (rGamma \<circ> (\<lambda>z. z + 1 :: complex)) z) has_field_derivative
+ (- 1) ^ Suc n * fact (Suc n)) (at (- of_nat (Suc n)))"
+ by (rule derivative_eq_intros refl A DERIV_chain)+
+ (simp add: algebra_simps rGamma_complex_altdef)
+ thus ?case by (simp add: rGamma_complex_plus1)
+qed
+
+instance proof
+ fix z :: complex show "(rGamma z = 0) \<longleftrightarrow> (\<exists>n. z = - of_nat n)"
+ by (auto simp: rGamma_complex_altdef elim!: nonpos_Ints_cases')
+next
+ fix z :: complex assume "\<And>n. z \<noteq> - of_nat n"
+ hence "z \<notin> \<int>\<^sub>\<le>\<^sub>0" by (auto elim!: nonpos_Ints_cases')
+ from has_field_derivative_rGamma_complex_no_nonpos_Int[OF this]
+ show "let d = (THE d. (\<lambda>n. \<Sum>k<n. inverse (of_nat (Suc k)) - inverse (z + of_nat k))
+ \<longlonglongrightarrow> d) - euler_mascheroni *\<^sub>R 1 in (\<lambda>y. (rGamma y - rGamma z +
+ rGamma z * d * (y - z)) /\<^sub>R cmod (y - z)) \<midarrow>z\<rightarrow> 0"
+ by (simp add: has_field_derivative_def has_derivative_def Digamma_def sums_def [abs_def]
+ netlimit_at of_real_def[symmetric] suminf_def)
+next
+ fix n :: nat
+ from has_field_derivative_rGamma_complex_nonpos_Int[of n]
+ show "let z = - of_nat n in (\<lambda>y. (rGamma y - rGamma z - (- 1) ^ n * setprod of_nat {1..n} *
+ (y - z)) /\<^sub>R cmod (y - z)) \<midarrow>z\<rightarrow> 0"
+ by (simp add: has_field_derivative_def has_derivative_def fact_altdef netlimit_at Let_def)
+next
+ fix z :: complex
+ from rGamma_series_complex_converges[of z] have "rGamma_series z \<longlonglongrightarrow> rGamma z"
+ by (simp add: convergent_LIMSEQ_iff rGamma_complex_def)
+ thus "let fact' = \<lambda>n. setprod of_nat {1..n};
+ exp = \<lambda>x. THE e. (\<lambda>n. \<Sum>k<n. x ^ k /\<^sub>R fact k) \<longlonglongrightarrow> e;
+ pochhammer' = \<lambda>a n. \<Prod>n = 0..n. a + of_nat n
+ in (\<lambda>n. pochhammer' z n / (fact' n * exp (z * ln (real_of_nat n) *\<^sub>R 1))) \<longlonglongrightarrow> rGamma z"
+ by (simp add: fact_altdef pochhammer_Suc_setprod rGamma_series_def [abs_def] exp_def
+ of_real_def [symmetric] suminf_def sums_def [abs_def])
+qed
+
+end
+end
+
+
+lemma Gamma_complex_altdef:
+ "Gamma z = (if z \<in> \<int>\<^sub>\<le>\<^sub>0 then 0 else exp (ln_Gamma (z :: complex)))"
+ unfolding Gamma_def rGamma_complex_altdef by (simp add: exp_minus)
+
+lemma cnj_rGamma: "cnj (rGamma z) = rGamma (cnj z)"
+proof -
+ have "rGamma_series (cnj z) = (\<lambda>n. cnj (rGamma_series z n))"
+ by (intro ext) (simp_all add: rGamma_series_def exp_cnj)
+ also have "... \<longlonglongrightarrow> cnj (rGamma z)" by (intro tendsto_cnj tendsto_intros)
+ finally show ?thesis unfolding rGamma_complex_def by (intro sym[OF limI])
+qed
+
+lemma cnj_Gamma: "cnj (Gamma z) = Gamma (cnj z)"
+ unfolding Gamma_def by (simp add: cnj_rGamma)
+
+lemma Gamma_complex_real:
+ "z \<in> \<real> \<Longrightarrow> Gamma z \<in> (\<real> :: complex set)" and rGamma_complex_real: "z \<in> \<real> \<Longrightarrow> rGamma z \<in> \<real>"
+ by (simp_all add: Reals_cnj_iff cnj_Gamma cnj_rGamma)
+
+lemma complex_differentiable_rGamma: "rGamma complex_differentiable (at z within A)"
+ using has_field_derivative_rGamma[of z] unfolding complex_differentiable_def by blast
+
+lemma holomorphic_on_rGamma: "rGamma holomorphic_on A"
+ unfolding holomorphic_on_def by (auto intro!: complex_differentiable_rGamma)
+
+lemma analytic_on_rGamma: "rGamma analytic_on A"
+ unfolding analytic_on_def by (auto intro!: exI[of _ 1] holomorphic_on_rGamma)
+
+
+lemma complex_differentiable_Gamma: "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Gamma complex_differentiable (at z within A)"
+ using has_field_derivative_Gamma[of z] unfolding complex_differentiable_def by auto
+
+lemma holomorphic_on_Gamma: "A \<inter> \<int>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> Gamma holomorphic_on A"
+ unfolding holomorphic_on_def by (auto intro!: complex_differentiable_Gamma)
+
+lemma analytic_on_Gamma: "A \<inter> \<int>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> Gamma analytic_on A"
+ by (rule analytic_on_subset[of _ "UNIV - \<int>\<^sub>\<le>\<^sub>0"], subst analytic_on_open)
+ (auto intro!: holomorphic_on_Gamma)
+
+lemma has_field_derivative_rGamma_complex' [derivative_intros]:
+ "(rGamma has_field_derivative (if z \<in> \<int>\<^sub>\<le>\<^sub>0 then (-1)^(nat \<lfloor>-Re z\<rfloor>) * fact (nat \<lfloor>-Re z\<rfloor>) else
+ -rGamma z * Digamma z)) (at z within A)"
+ using has_field_derivative_rGamma[of z] by (auto elim!: nonpos_Ints_cases')
+
+declare has_field_derivative_rGamma_complex'[THEN DERIV_chain2, derivative_intros]
+
+
+lemma complex_differentiable_Polygamma:
+ "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Polygamma n complex_differentiable (at z within A)"
+ using has_field_derivative_Polygamma[of z n] unfolding complex_differentiable_def by auto
+
+lemma holomorphic_on_Polygamma: "A \<inter> \<int>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> Polygamma n holomorphic_on A"
+ unfolding holomorphic_on_def by (auto intro!: complex_differentiable_Polygamma)
+
+lemma analytic_on_Polygamma: "A \<inter> \<int>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> Polygamma n analytic_on A"
+ by (rule analytic_on_subset[of _ "UNIV - \<int>\<^sub>\<le>\<^sub>0"], subst analytic_on_open)
+ (auto intro!: holomorphic_on_Polygamma)
+
+
+
+text \<open>The real Gamma function\<close>
+
+lemma rGamma_series_real:
+ "eventually (\<lambda>n. rGamma_series x n = Re (rGamma_series (of_real x) n)) sequentially"
+ using eventually_gt_at_top[of "0 :: nat"]
+proof eventually_elim
+ fix n :: nat assume n: "n > 0"
+ have "Re (rGamma_series (of_real x) n) =
+ Re (of_real (pochhammer x (Suc n)) / (fact n * exp (of_real (x * ln (real_of_nat n)))))"
+ using n by (simp add: rGamma_series_def powr_def Ln_of_nat pochhammer_of_real)
+ also from n have "\<dots> = Re (of_real ((pochhammer x (Suc n)) /
+ (fact n * (exp (x * ln (real_of_nat n))))))"
+ by (subst exp_of_real) simp
+ also from n have "\<dots> = rGamma_series x n"
+ by (subst Re_complex_of_real) (simp add: rGamma_series_def powr_def)
+ finally show "rGamma_series x n = Re (rGamma_series (of_real x) n)" ..
+qed
+
+instantiation real :: Gamma
+begin
+
+definition "rGamma_real x = Re (rGamma (of_real x :: complex))"
+
+instance proof
+ fix x :: real
+ have "rGamma x = Re (rGamma (of_real x))" by (simp add: rGamma_real_def)
+ also have "of_real \<dots> = rGamma (of_real x :: complex)"
+ by (intro of_real_Re rGamma_complex_real) simp_all
+ also have "\<dots> = 0 \<longleftrightarrow> x \<in> \<int>\<^sub>\<le>\<^sub>0" by (simp add: rGamma_eq_zero_iff of_real_in_nonpos_Ints_iff)
+ also have "\<dots> \<longleftrightarrow> (\<exists>n. x = - of_nat n)" by (auto elim!: nonpos_Ints_cases')
+ finally show "(rGamma x) = 0 \<longleftrightarrow> (\<exists>n. x = - real_of_nat n)" by simp
+next
+ fix x :: real assume "\<And>n. x \<noteq> - of_nat n"
+ hence "complex_of_real x \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ by (subst of_real_in_nonpos_Ints_iff) (auto elim!: nonpos_Ints_cases')
+ moreover from this have "x \<noteq> 0" by auto
+ ultimately have "(rGamma has_field_derivative - rGamma x * Digamma x) (at x)"
+ by (fastforce intro!: derivative_eq_intros has_vector_derivative_real_complex
+ simp: Polygamma_of_real rGamma_real_def [abs_def])
+ thus "let d = (THE d. (\<lambda>n. \<Sum>k<n. inverse (of_nat (Suc k)) - inverse (x + of_nat k))
+ \<longlonglongrightarrow> d) - euler_mascheroni *\<^sub>R 1 in (\<lambda>y. (rGamma y - rGamma x +
+ rGamma x * d * (y - x)) /\<^sub>R norm (y - x)) \<midarrow>x\<rightarrow> 0"
+ by (simp add: has_field_derivative_def has_derivative_def Digamma_def sums_def [abs_def]
+ netlimit_at of_real_def[symmetric] suminf_def)
+next
+ fix n :: nat
+ have "(rGamma has_field_derivative (-1)^n * fact n) (at (- of_nat n :: real))"
+ by (fastforce intro!: derivative_eq_intros has_vector_derivative_real_complex
+ simp: Polygamma_of_real rGamma_real_def [abs_def])
+ thus "let x = - of_nat n in (\<lambda>y. (rGamma y - rGamma x - (- 1) ^ n * setprod of_nat {1..n} *
+ (y - x)) /\<^sub>R norm (y - x)) \<midarrow>x::real\<rightarrow> 0"
+ by (simp add: has_field_derivative_def has_derivative_def fact_altdef netlimit_at Let_def)
+next
+ fix x :: real
+ have "rGamma_series x \<longlonglongrightarrow> rGamma x"
+ proof (rule Lim_transform_eventually)
+ show "(\<lambda>n. Re (rGamma_series (of_real x) n)) \<longlonglongrightarrow> rGamma x" unfolding rGamma_real_def
+ by (intro tendsto_intros)
+ qed (insert rGamma_series_real, simp add: eq_commute)
+ thus "let fact' = \<lambda>n. setprod of_nat {1..n};
+ exp = \<lambda>x. THE e. (\<lambda>n. \<Sum>k<n. x ^ k /\<^sub>R fact k) \<longlonglongrightarrow> e;
+ pochhammer' = \<lambda>a n. \<Prod>n = 0..n. a + of_nat n
+ in (\<lambda>n. pochhammer' x n / (fact' n * exp (x * ln (real_of_nat n) *\<^sub>R 1))) \<longlonglongrightarrow> rGamma x"
+ by (simp add: fact_altdef pochhammer_Suc_setprod rGamma_series_def [abs_def] exp_def
+ of_real_def [symmetric] suminf_def sums_def [abs_def])
+qed
+
+end
+
+
+lemma rGamma_complex_of_real: "rGamma (complex_of_real x) = complex_of_real (rGamma x)"
+ unfolding rGamma_real_def using rGamma_complex_real by simp
+
+lemma Gamma_complex_of_real: "Gamma (complex_of_real x) = complex_of_real (Gamma x)"
+ unfolding Gamma_def by (simp add: rGamma_complex_of_real)
+
+lemma rGamma_real_altdef: "rGamma x = lim (rGamma_series (x :: real))"
+ by (rule sym, rule limI, rule tendsto_intros)
+
+lemma Gamma_real_altdef1: "Gamma x = lim (Gamma_series (x :: real))"
+ by (rule sym, rule limI, rule tendsto_intros)
+
+lemma Gamma_real_altdef2: "Gamma x = Re (Gamma (of_real x))"
+ using rGamma_complex_real[OF Reals_of_real[of x]]
+ by (elim Reals_cases)
+ (simp only: Gamma_def rGamma_real_def of_real_inverse[symmetric] Re_complex_of_real)
+
+lemma ln_Gamma_series_complex_of_real:
+ "x > 0 \<Longrightarrow> n > 0 \<Longrightarrow> ln_Gamma_series (complex_of_real x) n = of_real (ln_Gamma_series x n)"
+proof -
+ assume xn: "x > 0" "n > 0"
+ have "Ln (complex_of_real x / of_nat k + 1) = of_real (ln (x / of_nat k + 1))" if "k \<ge> 1" for k
+ using that xn by (subst Ln_of_real [symmetric]) (auto intro!: add_nonneg_pos simp: field_simps)
+ with xn show ?thesis by (simp add: ln_Gamma_series_def Ln_of_nat Ln_of_real)
+qed
+
+lemma ln_Gamma_real_converges:
+ assumes "(x::real) > 0"
+ shows "convergent (ln_Gamma_series x)"
+proof -
+ have "(\<lambda>n. ln_Gamma_series (complex_of_real x) n) \<longlonglongrightarrow> ln_Gamma (of_real x)" using assms
+ by (intro ln_Gamma_complex_LIMSEQ) (auto simp: of_real_in_nonpos_Ints_iff)
+ moreover from eventually_gt_at_top[of "0::nat"]
+ have "eventually (\<lambda>n. complex_of_real (ln_Gamma_series x n) =
+ ln_Gamma_series (complex_of_real x) n) sequentially"
+ by eventually_elim (simp add: ln_Gamma_series_complex_of_real assms)
+ ultimately have "(\<lambda>n. complex_of_real (ln_Gamma_series x n)) \<longlonglongrightarrow> ln_Gamma (of_real x)"
+ by (subst tendsto_cong) assumption+
+ from tendsto_Re[OF this] show ?thesis by (auto simp: convergent_def)
+qed
+
+lemma ln_Gamma_real_LIMSEQ: "(x::real) > 0 \<Longrightarrow> ln_Gamma_series x \<longlonglongrightarrow> ln_Gamma x"
+ using ln_Gamma_real_converges[of x] unfolding ln_Gamma_def by (simp add: convergent_LIMSEQ_iff)
+
+lemma ln_Gamma_complex_of_real: "x > 0 \<Longrightarrow> ln_Gamma (complex_of_real x) = of_real (ln_Gamma x)"
+proof (unfold ln_Gamma_def, rule limI, rule Lim_transform_eventually)
+ assume x: "x > 0"
+ show "eventually (\<lambda>n. of_real (ln_Gamma_series x n) =
+ ln_Gamma_series (complex_of_real x) n) sequentially"
+ using eventually_gt_at_top[of "0::nat"]
+ by eventually_elim (simp add: ln_Gamma_series_complex_of_real x)
+qed (intro tendsto_of_real, insert ln_Gamma_real_LIMSEQ[of x], simp add: ln_Gamma_def)
+
+lemma Gamma_real_pos_exp: "x > (0 :: real) \<Longrightarrow> Gamma x = exp (ln_Gamma x)"
+ by (auto simp: Gamma_real_altdef2 Gamma_complex_altdef of_real_in_nonpos_Ints_iff
+ ln_Gamma_complex_of_real exp_of_real)
+
+lemma ln_Gamma_real_pos: "x > 0 \<Longrightarrow> ln_Gamma x = ln (Gamma x :: real)"
+ unfolding Gamma_real_pos_exp by simp
+
+lemma Gamma_real_pos: "x > (0::real) \<Longrightarrow> Gamma x > 0"
+ by (simp add: Gamma_real_pos_exp)
+
+lemma has_field_derivative_ln_Gamma_real [derivative_intros]:
+ assumes x: "x > (0::real)"
+ shows "(ln_Gamma has_field_derivative Digamma x) (at x)"
+proof (subst DERIV_cong_ev[OF refl _ refl])
+ from assms show "((Re \<circ> ln_Gamma \<circ> complex_of_real) has_field_derivative Digamma x) (at x)"
+ by (auto intro!: derivative_eq_intros has_vector_derivative_real_complex
+ simp: Polygamma_of_real o_def)
+ from eventually_nhds_in_nhd[of x "{0<..}"] assms
+ show "eventually (\<lambda>y. ln_Gamma y = (Re \<circ> ln_Gamma \<circ> of_real) y) (nhds x)"
+ by (auto elim!: eventually_mono simp: ln_Gamma_complex_of_real interior_open)
+qed
+
+declare has_field_derivative_ln_Gamma_real[THEN DERIV_chain2, derivative_intros]
+
+
+lemma has_field_derivative_rGamma_real' [derivative_intros]:
+ "(rGamma has_field_derivative (if x \<in> \<int>\<^sub>\<le>\<^sub>0 then (-1)^(nat \<lfloor>-x\<rfloor>) * fact (nat \<lfloor>-x\<rfloor>) else
+ -rGamma x * Digamma x)) (at x within A)"
+ using has_field_derivative_rGamma[of x] by (force elim!: nonpos_Ints_cases')
+
+declare has_field_derivative_rGamma_real'[THEN DERIV_chain2, derivative_intros]
+
+lemma Polygamma_real_odd_pos:
+ assumes "(x::real) \<notin> \<int>\<^sub>\<le>\<^sub>0" "odd n"
+ shows "Polygamma n x > 0"
+proof -
+ from assms have "x \<noteq> 0" by auto
+ with assms show ?thesis
+ unfolding Polygamma_def using Polygamma_converges'[of x "Suc n"]
+ by (auto simp: zero_less_power_eq simp del: power_Suc
+ dest: plus_of_nat_eq_0_imp intro!: mult_pos_pos suminf_pos)
+qed
+
+lemma Polygamma_real_even_neg:
+ assumes "(x::real) > 0" "n > 0" "even n"
+ shows "Polygamma n x < 0"
+ using assms unfolding Polygamma_def using Polygamma_converges'[of x "Suc n"]
+ by (auto intro!: mult_pos_pos suminf_pos)
+
+lemma Polygamma_real_strict_mono:
+ assumes "x > 0" "x < (y::real)" "even n"
+ shows "Polygamma n x < Polygamma n y"
+proof -
+ have "\<exists>\<xi>. x < \<xi> \<and> \<xi> < y \<and> Polygamma n y - Polygamma n x = (y - x) * Polygamma (Suc n) \<xi>"
+ using assms by (intro MVT2 derivative_intros impI allI) (auto elim!: nonpos_Ints_cases)
+ then guess \<xi> by (elim exE conjE) note \<xi> = this
+ note \<xi>(3)
+ also from \<xi>(1,2) assms have "(y - x) * Polygamma (Suc n) \<xi> > 0"
+ by (intro mult_pos_pos Polygamma_real_odd_pos) (auto elim!: nonpos_Ints_cases)
+ finally show ?thesis by simp
+qed
+
+lemma Polygamma_real_strict_antimono:
+ assumes "x > 0" "x < (y::real)" "odd n"
+ shows "Polygamma n x > Polygamma n y"
+proof -
+ have "\<exists>\<xi>. x < \<xi> \<and> \<xi> < y \<and> Polygamma n y - Polygamma n x = (y - x) * Polygamma (Suc n) \<xi>"
+ using assms by (intro MVT2 derivative_intros impI allI) (auto elim!: nonpos_Ints_cases)
+ then guess \<xi> by (elim exE conjE) note \<xi> = this
+ note \<xi>(3)
+ also from \<xi>(1,2) assms have "(y - x) * Polygamma (Suc n) \<xi> < 0"
+ by (intro mult_pos_neg Polygamma_real_even_neg) simp_all
+ finally show ?thesis by simp
+qed
+
+lemma Polygamma_real_mono:
+ assumes "x > 0" "x \<le> (y::real)" "even n"
+ shows "Polygamma n x \<le> Polygamma n y"
+ using Polygamma_real_strict_mono[OF assms(1) _ assms(3), of y] assms(2)
+ by (cases "x = y") simp_all
+
+lemma Digamma_real_ge_three_halves_pos:
+ assumes "x \<ge> 3/2"
+ shows "Digamma (x :: real) > 0"
+proof -
+ have "0 < Digamma (3/2 :: real)" by (fact Digamma_real_three_halves_pos)
+ also from assms have "\<dots> \<le> Digamma x" by (intro Polygamma_real_mono) simp_all
+ finally show ?thesis .
+qed
+
+lemma ln_Gamma_real_strict_mono:
+ assumes "x \<ge> 3/2" "x < y"
+ shows "ln_Gamma (x :: real) < ln_Gamma y"
+proof -
+ have "\<exists>\<xi>. x < \<xi> \<and> \<xi> < y \<and> ln_Gamma y - ln_Gamma x = (y - x) * Digamma \<xi>"
+ using assms by (intro MVT2 derivative_intros impI allI) (auto elim!: nonpos_Ints_cases)
+ then guess \<xi> by (elim exE conjE) note \<xi> = this
+ note \<xi>(3)
+ also from \<xi>(1,2) assms have "(y - x) * Digamma \<xi> > 0"
+ by (intro mult_pos_pos Digamma_real_ge_three_halves_pos) simp_all
+ finally show ?thesis by simp
+qed
+
+lemma Gamma_real_strict_mono:
+ assumes "x \<ge> 3/2" "x < y"
+ shows "Gamma (x :: real) < Gamma y"
+proof -
+ from Gamma_real_pos_exp[of x] assms have "Gamma x = exp (ln_Gamma x)" by simp
+ also have "\<dots> < exp (ln_Gamma y)" by (intro exp_less_mono ln_Gamma_real_strict_mono assms)
+ also from Gamma_real_pos_exp[of y] assms have "\<dots> = Gamma y" by simp
+ finally show ?thesis .
+qed
+
+lemma log_convex_Gamma_real: "convex_on {0<..} (ln \<circ> Gamma :: real \<Rightarrow> real)"
+ by (rule convex_on_realI[of _ _ Digamma])
+ (auto intro!: derivative_eq_intros Polygamma_real_mono Gamma_real_pos
+ simp: o_def Gamma_eq_zero_iff elim!: nonpos_Ints_cases')
+
+
+subsection \<open>Beta function\<close>
+
+definition Beta where "Beta a b = Gamma a * Gamma b / Gamma (a + b)"
+
+lemma Beta_altdef: "Beta a b = Gamma a * Gamma b * rGamma (a + b)"
+ by (simp add: inverse_eq_divide Beta_def Gamma_def)
+
+lemma Beta_commute: "Beta a b = Beta b a"
+ unfolding Beta_def by (simp add: ac_simps)
+
+lemma has_field_derivative_Beta1 [derivative_intros]:
+ assumes "x \<notin> \<int>\<^sub>\<le>\<^sub>0" "x + y \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ shows "((\<lambda>x. Beta x y) has_field_derivative (Beta x y * (Digamma x - Digamma (x + y))))
+ (at x within A)" unfolding Beta_altdef
+ by (rule DERIV_cong, (rule derivative_intros assms)+) (simp add: algebra_simps)
+
+lemma has_field_derivative_Beta2 [derivative_intros]:
+ assumes "y \<notin> \<int>\<^sub>\<le>\<^sub>0" "x + y \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ shows "((\<lambda>y. Beta x y) has_field_derivative (Beta x y * (Digamma y - Digamma (x + y))))
+ (at y within A)"
+ using has_field_derivative_Beta1[of y x A] assms by (simp add: Beta_commute add_ac)
+
+lemma Beta_plus1_plus1:
+ assumes "x \<notin> \<int>\<^sub>\<le>\<^sub>0" "y \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ shows "Beta (x + 1) y + Beta x (y + 1) = Beta x y"
+proof -
+ have "Beta (x + 1) y + Beta x (y + 1) =
+ (Gamma (x + 1) * Gamma y + Gamma x * Gamma (y + 1)) * rGamma ((x + y) + 1)"
+ by (simp add: Beta_altdef add_divide_distrib algebra_simps)
+ also have "\<dots> = (Gamma x * Gamma y) * ((x + y) * rGamma ((x + y) + 1))"
+ by (subst assms[THEN Gamma_plus1])+ (simp add: algebra_simps)
+ also from assms have "\<dots> = Beta x y" unfolding Beta_altdef by (subst rGamma_plus1) simp
+ finally show ?thesis .
+qed
+
+lemma Beta_plus1_left:
+ assumes "x \<notin> \<int>\<^sub>\<le>\<^sub>0" "y \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ shows "(x + y) * Beta (x + 1) y = x * Beta x y"
+proof -
+ have "(x + y) * Beta (x + 1) y = Gamma (x + 1) * Gamma y * ((x + y) * rGamma ((x + y) + 1))"
+ unfolding Beta_altdef by (simp only: ac_simps)
+ also have "\<dots> = x * Beta x y" unfolding Beta_altdef
+ by (subst assms[THEN Gamma_plus1] rGamma_plus1)+ (simp only: ac_simps)
+ finally show ?thesis .
+qed
+
+lemma Beta_plus1_right:
+ assumes "x \<notin> \<int>\<^sub>\<le>\<^sub>0" "y \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ shows "(x + y) * Beta x (y + 1) = y * Beta x y"
+ using Beta_plus1_left[of y x] assms by (simp add: Beta_commute add.commute)
+
+lemma Gamma_Gamma_Beta:
+ assumes "x \<notin> \<int>\<^sub>\<le>\<^sub>0" "y \<notin> \<int>\<^sub>\<le>\<^sub>0" "x + y \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ shows "Gamma x * Gamma y = Beta x y * Gamma (x + y)"
+ unfolding Beta_altdef using assms Gamma_eq_zero_iff[of "x+y"]
+ by (simp add: rGamma_inverse_Gamma)
+
+
+
+subsection \<open>Legendre duplication theorem\<close>
+
+context
+begin
+
+private lemma Gamma_legendre_duplication_aux:
+ fixes z :: "'a :: Gamma"
+ assumes "z \<notin> \<int>\<^sub>\<le>\<^sub>0" "z + 1/2 \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ shows "Gamma z * Gamma (z + 1/2) = exp ((1 - 2*z) * of_real (ln 2)) * Gamma (1/2) * Gamma (2*z)"
+proof -
+ let ?powr = "\<lambda>b a. exp (a * of_real (ln (of_nat b)))"
+ let ?h = "\<lambda>n. (fact (n-1))\<^sup>2 / fact (2*n-1) * of_nat (2^(2*n)) *
+ exp (1/2 * of_real (ln (real_of_nat n)))"
+ {
+ fix z :: 'a assume z: "z \<notin> \<int>\<^sub>\<le>\<^sub>0" "z + 1/2 \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ let ?g = "\<lambda>n. ?powr 2 (2*z) * Gamma_series' z n * Gamma_series' (z + 1/2) n /
+ Gamma_series' (2*z) (2*n)"
+ have "eventually (\<lambda>n. ?g n = ?h n) sequentially" using eventually_gt_at_top
+ proof eventually_elim
+ fix n :: nat assume n: "n > 0"
+ let ?f = "fact (n - 1) :: 'a" and ?f' = "fact (2*n - 1) :: 'a"
+ have A: "exp t * exp t = exp (2*t :: 'a)" for t by (subst exp_add [symmetric]) simp
+ have A: "Gamma_series' z n * Gamma_series' (z + 1/2) n = ?f^2 * ?powr n (2*z + 1/2) /
+ (pochhammer z n * pochhammer (z + 1/2) n)"
+ by (simp add: Gamma_series'_def exp_add ring_distribs power2_eq_square A mult_ac)
+ have B: "Gamma_series' (2*z) (2*n) =
+ ?f' * ?powr 2 (2*z) * ?powr n (2*z) /
+ (of_nat (2^(2*n)) * pochhammer z n * pochhammer (z+1/2) n)" using n
+ by (simp add: Gamma_series'_def ln_mult exp_add ring_distribs pochhammer_double)
+ from z have "pochhammer z n \<noteq> 0" by (auto dest: pochhammer_eq_0_imp_nonpos_Int)
+ moreover from z have "pochhammer (z + 1/2) n \<noteq> 0" by (auto dest: pochhammer_eq_0_imp_nonpos_Int)
+ ultimately have "?powr 2 (2*z) * (Gamma_series' z n * Gamma_series' (z + 1/2) n) / Gamma_series' (2*z) (2*n) =
+ ?f^2 / ?f' * of_nat (2^(2*n)) * (?powr n ((4*z + 1)/2) * ?powr n (-2*z))"
+ using n unfolding A B by (simp add: divide_simps exp_minus)
+ also have "?powr n ((4*z + 1)/2) * ?powr n (-2*z) = ?powr n (1/2)"
+ by (simp add: algebra_simps exp_add[symmetric] add_divide_distrib)
+ finally show "?g n = ?h n" by (simp only: mult_ac)
+ qed
+
+ moreover from z double_in_nonpos_Ints_imp[of z] have "2 * z \<notin> \<int>\<^sub>\<le>\<^sub>0" by auto
+ hence "?g \<longlonglongrightarrow> ?powr 2 (2*z) * Gamma z * Gamma (z+1/2) / Gamma (2*z)"
+ using lim_subseq[of "op * 2", OF _ Gamma_series'_LIMSEQ, of "2*z"]
+ by (intro tendsto_intros Gamma_series'_LIMSEQ)
+ (simp_all add: o_def subseq_def Gamma_eq_zero_iff)
+ ultimately have "?h \<longlonglongrightarrow> ?powr 2 (2*z) * Gamma z * Gamma (z+1/2) / Gamma (2*z)"
+ by (rule Lim_transform_eventually)
+ } note lim = this
+
+ from assms double_in_nonpos_Ints_imp[of z] have z': "2 * z \<notin> \<int>\<^sub>\<le>\<^sub>0" by auto
+ from fraction_not_in_ints[of 2 1] have "(1/2 :: 'a) \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ by (intro not_in_Ints_imp_not_in_nonpos_Ints) simp_all
+ with lim[of "1/2 :: 'a"] have "?h \<longlonglongrightarrow> 2 * Gamma (1 / 2 :: 'a)" by (simp add: exp_of_real)
+ from LIMSEQ_unique[OF this lim[OF assms]] z' show ?thesis
+ by (simp add: divide_simps Gamma_eq_zero_iff ring_distribs exp_diff exp_of_real ac_simps)
+qed
+
+(* TODO: perhaps this is unnecessary once we have the fact that a holomorphic function is
+ infinitely differentiable *)
+private lemma Gamma_reflection_aux:
+ defines "h \<equiv> \<lambda>z::complex. if z \<in> \<int> then 0 else
+ (of_real pi * cot (of_real pi*z) + Digamma z - Digamma (1 - z))"
+ defines "a \<equiv> complex_of_real pi"
+ obtains h' where "continuous_on UNIV h'" "\<And>z. (h has_field_derivative (h' z)) (at z)"
+proof -
+ def f \<equiv> "\<lambda>n. a * of_real (cos_coeff (n+1) - sin_coeff (n+2))"
+ def F \<equiv> "\<lambda>z. if z = 0 then 0 else (cos (a*z) - sin (a*z)/(a*z)) / z"
+ def g \<equiv> "\<lambda>n. complex_of_real (sin_coeff (n+1))"
+ def G \<equiv> "\<lambda>z. if z = 0 then 1 else sin (a*z)/(a*z)"
+ have a_nz: "a \<noteq> 0" unfolding a_def by simp
+
+ have "(\<lambda>n. f n * (a*z)^n) sums (F z) \<and> (\<lambda>n. g n * (a*z)^n) sums (G z)"
+ if "abs (Re z) < 1" for z
+ proof (cases "z = 0"; rule conjI)
+ assume "z \<noteq> 0"
+ note z = this that
+
+ from z have sin_nz: "sin (a*z) \<noteq> 0" unfolding a_def by (auto simp: sin_eq_0)
+ have "(\<lambda>n. of_real (sin_coeff n) * (a*z)^n) sums (sin (a*z))" using sin_converges[of "a*z"]
+ by (simp add: scaleR_conv_of_real)
+ from sums_split_initial_segment[OF this, of 1]
+ have "(\<lambda>n. (a*z) * of_real (sin_coeff (n+1)) * (a*z)^n) sums (sin (a*z))" by (simp add: mult_ac)
+ from sums_mult[OF this, of "inverse (a*z)"] z a_nz
+ have A: "(\<lambda>n. g n * (a*z)^n) sums (sin (a*z)/(a*z))"
+ by (simp add: field_simps g_def)
+ with z show "(\<lambda>n. g n * (a*z)^n) sums (G z)" by (simp add: G_def)
+ from A z a_nz sin_nz have g_nz: "(\<Sum>n. g n * (a*z)^n) \<noteq> 0" by (simp add: sums_iff g_def)
+
+ have [simp]: "sin_coeff (Suc 0) = 1" by (simp add: sin_coeff_def)
+ from sums_split_initial_segment[OF sums_diff[OF cos_converges[of "a*z"] A], of 1]
+ have "(\<lambda>n. z * f n * (a*z)^n) sums (cos (a*z) - sin (a*z) / (a*z))"
+ by (simp add: mult_ac scaleR_conv_of_real ring_distribs f_def g_def)
+ from sums_mult[OF this, of "inverse z"] z assms
+ show "(\<lambda>n. f n * (a*z)^n) sums (F z)" by (simp add: divide_simps mult_ac f_def F_def)
+ next
+ assume z: "z = 0"
+ have "(\<lambda>n. f n * (a * z) ^ n) sums f 0" using powser_sums_zero[of f] z by simp
+ with z show "(\<lambda>n. f n * (a * z) ^ n) sums (F z)"
+ by (simp add: f_def F_def sin_coeff_def cos_coeff_def)
+ have "(\<lambda>n. g n * (a * z) ^ n) sums g 0" using powser_sums_zero[of g] z by simp
+ with z show "(\<lambda>n. g n * (a * z) ^ n) sums (G z)"
+ by (simp add: g_def G_def sin_coeff_def cos_coeff_def)
+ qed
+ note sums = conjunct1[OF this] conjunct2[OF this]
+
+ def h2 \<equiv> "\<lambda>z. (\<Sum>n. f n * (a*z)^n) / (\<Sum>n. g n * (a*z)^n) +
+ Digamma (1 + z) - Digamma (1 - z)"
+ def POWSER \<equiv> "\<lambda>f z. (\<Sum>n. f n * (z^n :: complex))"
+ def POWSER' \<equiv> "\<lambda>f z. (\<Sum>n. diffs f n * (z^n :: complex))"
+
+ def h2' \<equiv> "\<lambda>z. a * (POWSER g (a*z) * POWSER' f (a*z) - POWSER f (a*z) * POWSER' g (a*z)) /
+ (POWSER g (a*z))^2 + Polygamma 1 (1 + z) + Polygamma 1 (1 - z)"
+
+ have h_eq: "h t = h2 t" if "abs (Re t) < 1" for t
+ proof -
+ from that have t: "t \<in> \<int> \<longleftrightarrow> t = 0" by (auto elim!: Ints_cases simp: dist_0_norm)
+ hence "h t = a*cot (a*t) - 1/t + Digamma (1 + t) - Digamma (1 - t)"
+ unfolding h_def using Digamma_plus1[of t] by (force simp: field_simps a_def)
+ also have "a*cot (a*t) - 1/t = (F t) / (G t)"
+ using t by (auto simp add: divide_simps sin_eq_0 cot_def a_def F_def G_def)
+ also have "\<dots> = (\<Sum>n. f n * (a*t)^n) / (\<Sum>n. g n * (a*t)^n)"
+ using sums[of t] that by (simp add: sums_iff dist_0_norm)
+ finally show "h t = h2 t" by (simp only: h2_def)
+ qed
+
+ let ?A = "{z. abs (Re z) < 1}"
+ have "open ({z. Re z < 1} \<inter> {z. Re z > -1})"
+ using open_halfspace_Re_gt open_halfspace_Re_lt by auto
+ also have "({z. Re z < 1} \<inter> {z. Re z > -1}) = {z. abs (Re z) < 1}" by auto
+ finally have open_A: "open ?A" .
+ hence [simp]: "interior ?A = ?A" by (simp add: interior_open)
+
+ have summable_f: "summable (\<lambda>n. f n * z^n)" for z
+ by (rule powser_inside, rule sums_summable, rule sums[of "\<i> * of_real (norm z + 1) / a"])
+ (simp_all add: norm_mult a_def del: of_real_add)
+ have summable_g: "summable (\<lambda>n. g n * z^n)" for z
+ by (rule powser_inside, rule sums_summable, rule sums[of "\<i> * of_real (norm z + 1) / a"])
+ (simp_all add: norm_mult a_def del: of_real_add)
+ have summable_fg': "summable (\<lambda>n. diffs f n * z^n)" "summable (\<lambda>n. diffs g n * z^n)" for z
+ by (intro termdiff_converges_all summable_f summable_g)+
+ have "(POWSER f has_field_derivative (POWSER' f z)) (at z)"
+ "(POWSER g has_field_derivative (POWSER' g z)) (at z)" for z
+ unfolding POWSER_def POWSER'_def
+ by (intro termdiffs_strong_converges_everywhere summable_f summable_g)+
+ note derivs = this[THEN DERIV_chain2[OF _ DERIV_cmult[OF DERIV_ident]], unfolded POWSER_def]
+ have "isCont (POWSER f) z" "isCont (POWSER g) z" "isCont (POWSER' f) z" "isCont (POWSER' g) z"
+ for z unfolding POWSER_def POWSER'_def
+ by (intro isCont_powser_converges_everywhere summable_f summable_g summable_fg')+
+ note cont = this[THEN isCont_o2[rotated], unfolded POWSER_def POWSER'_def]
+
+ {
+ fix z :: complex assume z: "abs (Re z) < 1"
+ def d \<equiv> "\<i> * of_real (norm z + 1)"
+ have d: "abs (Re d) < 1" "norm z < norm d" by (simp_all add: d_def norm_mult del: of_real_add)
+ have "eventually (\<lambda>z. h z = h2 z) (nhds z)"
+ using eventually_nhds_in_nhd[of z ?A] using h_eq z
+ by (auto elim!: eventually_mono simp: dist_0_norm)
+
+ moreover from sums(2)[OF z] z have nz: "(\<Sum>n. g n * (a * z) ^ n) \<noteq> 0"
+ unfolding G_def by (auto simp: sums_iff sin_eq_0 a_def)
+ have A: "z \<in> \<int> \<longleftrightarrow> z = 0" using z by (auto elim!: Ints_cases)
+ have no_int: "1 + z \<in> \<int> \<longleftrightarrow> z = 0" using z Ints_diff[of "1+z" 1] A
+ by (auto elim!: nonpos_Ints_cases)
+ have no_int': "1 - z \<in> \<int> \<longleftrightarrow> z = 0" using z Ints_diff[of 1 "1-z"] A
+ by (auto elim!: nonpos_Ints_cases)
+ from no_int no_int' have no_int: "1 - z \<notin> \<int>\<^sub>\<le>\<^sub>0" "1 + z \<notin> \<int>\<^sub>\<le>\<^sub>0" by auto
+ have "(h2 has_field_derivative h2' z) (at z)" unfolding h2_def
+ by (rule DERIV_cong, (rule derivative_intros refl derivs[unfolded POWSER_def] nz no_int)+)
+ (auto simp: h2'_def POWSER_def field_simps power2_eq_square)
+ ultimately have deriv: "(h has_field_derivative h2' z) (at z)"
+ by (subst DERIV_cong_ev[OF refl _ refl])
+
+ from sums(2)[OF z] z have "(\<Sum>n. g n * (a * z) ^ n) \<noteq> 0"
+ unfolding G_def by (auto simp: sums_iff a_def sin_eq_0)
+ hence "isCont h2' z" using no_int unfolding h2'_def[abs_def] POWSER_def POWSER'_def
+ by (intro continuous_intros cont
+ continuous_on_compose2[OF _ continuous_on_Polygamma[of "{z. Re z > 0}"]]) auto
+ note deriv and this
+ } note A = this
+
+ interpret h: periodic_fun_simple' h
+ proof
+ fix z :: complex
+ show "h (z + 1) = h z"
+ proof (cases "z \<in> \<int>")
+ assume z: "z \<notin> \<int>"
+ hence A: "z + 1 \<notin> \<int>" "z \<noteq> 0" using Ints_diff[of "z+1" 1] by auto
+ hence "Digamma (z + 1) - Digamma (-z) = Digamma z - Digamma (-z + 1)"
+ by (subst (1 2) Digamma_plus1) simp_all
+ with A z show "h (z + 1) = h z"
+ by (simp add: h_def sin_plus_pi cos_plus_pi ring_distribs cot_def)
+ qed (simp add: h_def)
+ qed
+
+ have h2'_eq: "h2' (z - 1) = h2' z" if z: "Re z > 0" "Re z < 1" for z
+ proof -
+ have "((\<lambda>z. h (z - 1)) has_field_derivative h2' (z - 1)) (at z)"
+ by (rule DERIV_cong, rule DERIV_chain'[OF _ A(1)])
+ (insert z, auto intro!: derivative_eq_intros)
+ hence "(h has_field_derivative h2' (z - 1)) (at z)" by (subst (asm) h.minus_1)
+ moreover from z have "(h has_field_derivative h2' z) (at z)" by (intro A) simp_all
+ ultimately show "h2' (z - 1) = h2' z" by (rule DERIV_unique)
+ qed
+
+ def h2'' \<equiv> "\<lambda>z. h2' (z - of_int \<lfloor>Re z\<rfloor>)"
+ have deriv: "(h has_field_derivative h2'' z) (at z)" for z
+ proof -
+ fix z :: complex
+ have B: "\<bar>Re z - real_of_int \<lfloor>Re z\<rfloor>\<bar> < 1" by linarith
+ have "((\<lambda>t. h (t - of_int \<lfloor>Re z\<rfloor>)) has_field_derivative h2'' z) (at z)"
+ unfolding h2''_def by (rule DERIV_cong, rule DERIV_chain'[OF _ A(1)])
+ (insert B, auto intro!: derivative_intros)
+ thus "(h has_field_derivative h2'' z) (at z)" by (simp add: h.minus_of_int)
+ qed
+
+ have cont: "continuous_on UNIV h2''"
+ proof (intro continuous_at_imp_continuous_on ballI)
+ fix z :: complex
+ def r \<equiv> "\<lfloor>Re z\<rfloor>"
+ def A \<equiv> "{t. of_int r - 1 < Re t \<and> Re t < of_int r + 1}"
+ have "continuous_on A (\<lambda>t. h2' (t - of_int r))" unfolding A_def
+ by (intro continuous_at_imp_continuous_on isCont_o2[OF _ A(2)] ballI continuous_intros)
+ (simp_all add: abs_real_def)
+ moreover have "h2'' t = h2' (t - of_int r)" if t: "t \<in> A" for t
+ proof (cases "Re t \<ge> of_int r")
+ case True
+ from t have "of_int r - 1 < Re t" "Re t < of_int r + 1" by (simp_all add: A_def)
+ with True have "\<lfloor>Re t\<rfloor> = \<lfloor>Re z\<rfloor>" unfolding r_def by linarith
+ thus ?thesis by (auto simp: r_def h2''_def)
+ next
+ case False
+ from t have t: "of_int r - 1 < Re t" "Re t < of_int r + 1" by (simp_all add: A_def)
+ with False have t': "\<lfloor>Re t\<rfloor> = \<lfloor>Re z\<rfloor> - 1" unfolding r_def by linarith
+ moreover from t False have "h2' (t - of_int r + 1 - 1) = h2' (t - of_int r + 1)"
+ by (intro h2'_eq) simp_all
+ ultimately show ?thesis by (auto simp: r_def h2''_def algebra_simps t')
+ qed
+ ultimately have "continuous_on A h2''" by (subst continuous_on_cong[OF refl])
+ moreover {
+ have "open ({t. of_int r - 1 < Re t} \<inter> {t. of_int r + 1 > Re t})"
+ by (intro open_Int open_halfspace_Re_gt open_halfspace_Re_lt)
+ also have "{t. of_int r - 1 < Re t} \<inter> {t. of_int r + 1 > Re t} = A"
+ unfolding A_def by blast
+ finally have "open A" .
+ }
+ ultimately have C: "isCont h2'' t" if "t \<in> A" for t using that
+ by (subst (asm) continuous_on_eq_continuous_at) auto
+ have "of_int r - 1 < Re z" "Re z < of_int r + 1" unfolding r_def by linarith+
+ thus "isCont h2'' z" by (intro C) (simp_all add: A_def)
+ qed
+
+ from that[OF cont deriv] show ?thesis .
+qed
+
+lemma Gamma_reflection_complex:
+ fixes z :: complex
+ shows "Gamma z * Gamma (1 - z) = of_real pi / sin (of_real pi * z)"
+proof -
+ let ?g = "\<lambda>z::complex. Gamma z * Gamma (1 - z) * sin (of_real pi * z)"
+ def g \<equiv> "\<lambda>z::complex. if z \<in> \<int> then of_real pi else ?g z"
+ let ?h = "\<lambda>z::complex. (of_real pi * cot (of_real pi*z) + Digamma z - Digamma (1 - z))"
+ def h \<equiv> "\<lambda>z::complex. if z \<in> \<int> then 0 else ?h z"
+
+ -- \<open>@{term g} is periodic with period 1.\<close>
+ interpret g: periodic_fun_simple' g
+ proof
+ fix z :: complex
+ show "g (z + 1) = g z"
+ proof (cases "z \<in> \<int>")
+ case False
+ hence "z * g z = z * Beta z (- z + 1) * sin (of_real pi * z)" by (simp add: g_def Beta_def)
+ also have "z * Beta z (- z + 1) = (z + 1 + -z) * Beta (z + 1) (- z + 1)"
+ using False Ints_diff[of 1 "1 - z"] nonpos_Ints_subset_Ints
+ by (subst Beta_plus1_left [symmetric]) auto
+ also have "\<dots> * sin (of_real pi * z) = z * (Beta (z + 1) (-z) * sin (of_real pi * (z + 1)))"
+ using False Ints_diff[of "z+1" 1] Ints_minus[of "-z"] nonpos_Ints_subset_Ints
+ by (subst Beta_plus1_right) (auto simp: ring_distribs sin_plus_pi)
+ also from False have "Beta (z + 1) (-z) * sin (of_real pi * (z + 1)) = g (z + 1)"
+ using Ints_diff[of "z+1" 1] by (auto simp: g_def Beta_def)
+ finally show "g (z + 1) = g z" using False by (subst (asm) mult_left_cancel) auto
+ qed (simp add: g_def)
+ qed
+
+ -- \<open>@{term g} is entire.\<close>
+ have g_g': "(g has_field_derivative (h z * g z)) (at z)" for z :: complex
+ proof (cases "z \<in> \<int>")
+ let ?h' = "\<lambda>z. Beta z (1 - z) * ((Digamma z - Digamma (1 - z)) * sin (z * of_real pi) +
+ of_real pi * cos (z * of_real pi))"
+ case False
+ from False have "eventually (\<lambda>t. t \<in> UNIV - \<int>) (nhds z)"
+ by (intro eventually_nhds_in_open) (auto simp: open_Diff)
+ hence "eventually (\<lambda>t. g t = ?g t) (nhds z)" by eventually_elim (simp add: g_def)
+ moreover {
+ from False Ints_diff[of 1 "1-z"] have "1 - z \<notin> \<int>" by auto
+ hence "(?g has_field_derivative ?h' z) (at z)" using nonpos_Ints_subset_Ints
+ by (auto intro!: derivative_eq_intros simp: algebra_simps Beta_def)
+ also from False have "sin (of_real pi * z) \<noteq> 0" by (subst sin_eq_0) auto
+ hence "?h' z = h z * g z"
+ using False unfolding g_def h_def cot_def by (simp add: field_simps Beta_def)
+ finally have "(?g has_field_derivative (h z * g z)) (at z)" .
+ }
+ ultimately show ?thesis by (subst DERIV_cong_ev[OF refl _ refl])
+ next
+ case True
+ then obtain n where z: "z = of_int n" by (auto elim!: Ints_cases)
+ let ?t = "(\<lambda>z::complex. if z = 0 then 1 else sin z / z) \<circ> (\<lambda>z. of_real pi * z)"
+ have deriv_0: "(g has_field_derivative 0) (at 0)"
+ proof (subst DERIV_cong_ev[OF refl _ refl])
+ show "eventually (\<lambda>z. g z = of_real pi * Gamma (1 + z) * Gamma (1 - z) * ?t z) (nhds 0)"
+ using eventually_nhds_ball[OF zero_less_one, of "0::complex"]
+ proof eventually_elim
+ fix z :: complex assume z: "z \<in> ball 0 1"
+ show "g z = of_real pi * Gamma (1 + z) * Gamma (1 - z) * ?t z"
+ proof (cases "z = 0")
+ assume z': "z \<noteq> 0"
+ with z have z'': "z \<notin> \<int>\<^sub>\<le>\<^sub>0" "z \<notin> \<int>" by (auto elim!: Ints_cases simp: dist_0_norm)
+ from Gamma_plus1[OF this(1)] have "Gamma z = Gamma (z + 1) / z" by simp
+ with z'' z' show ?thesis by (simp add: g_def ac_simps)
+ qed (simp add: g_def)
+ qed
+ have "(?t has_field_derivative (0 * of_real pi)) (at 0)"
+ using has_field_derivative_sin_z_over_z[of "UNIV :: complex set"]
+ by (intro DERIV_chain) simp_all
+ thus "((\<lambda>z. of_real pi * Gamma (1 + z) * Gamma (1 - z) * ?t z) has_field_derivative 0) (at 0)"
+ by (auto intro!: derivative_eq_intros simp: o_def)
+ qed
+
+ have "((g \<circ> (\<lambda>x. x - of_int n)) has_field_derivative 0 * 1) (at (of_int n))"
+ using deriv_0 by (intro DERIV_chain) (auto intro!: derivative_eq_intros)
+ also have "g \<circ> (\<lambda>x. x - of_int n) = g" by (intro ext) (simp add: g.minus_of_int)
+ finally show "(g has_field_derivative (h z * g z)) (at z)" by (simp add: z h_def)
+ qed
+
+ have g_eq: "g (z/2) * g ((z+1)/2) = Gamma (1/2)^2 * g z" if "Re z > -1" "Re z < 2" for z
+ proof (cases "z \<in> \<int>")
+ case True
+ with that have "z = 0 \<or> z = 1" by (force elim!: Ints_cases)
+ moreover have "g 0 * g (1/2) = Gamma (1/2)^2 * g 0"
+ using fraction_not_in_ints[where 'a = complex, of 2 1] by (simp add: g_def power2_eq_square)
+ moreover have "g (1/2) * g 1 = Gamma (1/2)^2 * g 1"
+ using fraction_not_in_ints[where 'a = complex, of 2 1]
+ by (simp add: g_def power2_eq_square Beta_def algebra_simps)
+ ultimately show ?thesis by force
+ next
+ case False
+ hence z: "z/2 \<notin> \<int>" "(z+1)/2 \<notin> \<int>" using Ints_diff[of "z+1" 1] by (auto elim!: Ints_cases)
+ hence z': "z/2 \<notin> \<int>\<^sub>\<le>\<^sub>0" "(z+1)/2 \<notin> \<int>\<^sub>\<le>\<^sub>0" by (auto elim!: nonpos_Ints_cases)
+ from z have "1-z/2 \<notin> \<int>" "1-((z+1)/2) \<notin> \<int>"
+ using Ints_diff[of 1 "1-z/2"] Ints_diff[of 1 "1-((z+1)/2)"] by auto
+ hence z'': "1-z/2 \<notin> \<int>\<^sub>\<le>\<^sub>0" "1-((z+1)/2) \<notin> \<int>\<^sub>\<le>\<^sub>0" by (auto elim!: nonpos_Ints_cases)
+ from z have "g (z/2) * g ((z+1)/2) =
+ (Gamma (z/2) * Gamma ((z+1)/2)) * (Gamma (1-z/2) * Gamma (1-((z+1)/2))) *
+ (sin (of_real pi * z/2) * sin (of_real pi * (z+1)/2))"
+ by (simp add: g_def)
+ also from z' Gamma_legendre_duplication_aux[of "z/2"]
+ have "Gamma (z/2) * Gamma ((z+1)/2) = exp ((1-z) * of_real (ln 2)) * Gamma (1/2) * Gamma z"
+ by (simp add: add_divide_distrib)
+ also from z'' Gamma_legendre_duplication_aux[of "1-(z+1)/2"]
+ have "Gamma (1-z/2) * Gamma (1-(z+1)/2) =
+ Gamma (1-z) * Gamma (1/2) * exp (z * of_real (ln 2))"
+ by (simp add: add_divide_distrib ac_simps)
+ finally have "g (z/2) * g ((z+1)/2) = Gamma (1/2)^2 * (Gamma z * Gamma (1-z) *
+ (2 * (sin (of_real pi*z/2) * sin (of_real pi*(z+1)/2))))"
+ by (simp add: add_ac power2_eq_square exp_add ring_distribs exp_diff exp_of_real)
+ also have "sin (of_real pi*(z+1)/2) = cos (of_real pi*z/2)"
+ using cos_sin_eq[of "- of_real pi * z/2", symmetric]
+ by (simp add: ring_distribs add_divide_distrib ac_simps)
+ also have "2 * (sin (of_real pi*z/2) * cos (of_real pi*z/2)) = sin (of_real pi * z)"
+ by (subst sin_times_cos) (simp add: field_simps)
+ also have "Gamma z * Gamma (1 - z) * sin (complex_of_real pi * z) = g z"
+ using \<open>z \<notin> \<int>\<close> by (simp add: g_def)
+ finally show ?thesis .
+ qed
+ have g_eq: "g (z/2) * g ((z+1)/2) = Gamma (1/2)^2 * g z" for z
+ proof -
+ def r \<equiv> "\<lfloor>Re z / 2\<rfloor>"
+ have "Gamma (1/2)^2 * g z = Gamma (1/2)^2 * g (z - of_int (2*r))" by (simp only: g.minus_of_int)
+ also have "of_int (2*r) = 2 * of_int r" by simp
+ also have "Re z - 2 * of_int r > -1" "Re z - 2 * of_int r < 2" unfolding r_def by linarith+
+ hence "Gamma (1/2)^2 * g (z - 2 * of_int r) =
+ g ((z - 2 * of_int r)/2) * g ((z - 2 * of_int r + 1)/2)"
+ unfolding r_def by (intro g_eq[symmetric]) simp_all
+ also have "(z - 2 * of_int r) / 2 = z/2 - of_int r" by simp
+ also have "g \<dots> = g (z/2)" by (rule g.minus_of_int)
+ also have "(z - 2 * of_int r + 1) / 2 = (z + 1)/2 - of_int r" by simp
+ also have "g \<dots> = g ((z+1)/2)" by (rule g.minus_of_int)
+ finally show ?thesis ..
+ qed
+
+ have g_nz [simp]: "g z \<noteq> 0" for z :: complex
+ unfolding g_def using Ints_diff[of 1 "1 - z"]
+ by (auto simp: Gamma_eq_zero_iff sin_eq_0 dest!: nonpos_Ints_Int)
+
+ have h_eq: "h z = (h (z/2) + h ((z+1)/2)) / 2" for z
+ proof -
+ have "((\<lambda>t. g (t/2) * g ((t+1)/2)) has_field_derivative
+ (g (z/2) * g ((z+1)/2)) * ((h (z/2) + h ((z+1)/2)) / 2)) (at z)"
+ by (auto intro!: derivative_eq_intros g_g'[THEN DERIV_chain2] simp: field_simps)
+ hence "((\<lambda>t. Gamma (1/2)^2 * g t) has_field_derivative
+ Gamma (1/2)^2 * g z * ((h (z/2) + h ((z+1)/2)) / 2)) (at z)"
+ by (subst (1 2) g_eq[symmetric]) simp
+ from DERIV_cmult[OF this, of "inverse ((Gamma (1/2))^2)"]
+ have "(g has_field_derivative (g z * ((h (z/2) + h ((z+1)/2))/2))) (at z)"
+ using fraction_not_in_ints[where 'a = complex, of 2 1]
+ by (simp add: divide_simps Gamma_eq_zero_iff not_in_Ints_imp_not_in_nonpos_Ints)
+ moreover have "(g has_field_derivative (g z * h z)) (at z)"
+ using g_g'[of z] by (simp add: ac_simps)
+ ultimately have "g z * h z = g z * ((h (z/2) + h ((z+1)/2))/2)"
+ by (intro DERIV_unique)
+ thus "h z = (h (z/2) + h ((z+1)/2)) / 2" by simp
+ qed
+
+ obtain h' where h'_cont: "continuous_on UNIV h'" and
+ h_h': "\<And>z. (h has_field_derivative h' z) (at z)"
+ unfolding h_def by (erule Gamma_reflection_aux)
+
+ have h'_eq: "h' z = (h' (z/2) + h' ((z+1)/2)) / 4" for z
+ proof -
+ have "((\<lambda>t. (h (t/2) + h ((t+1)/2)) / 2) has_field_derivative
+ ((h' (z/2) + h' ((z+1)/2)) / 4)) (at z)"
+ by (fastforce intro!: derivative_eq_intros h_h'[THEN DERIV_chain2])
+ hence "(h has_field_derivative ((h' (z/2) + h' ((z+1)/2))/4)) (at z)"
+ by (subst (asm) h_eq[symmetric])
+ from h_h' and this show "h' z = (h' (z/2) + h' ((z+1)/2)) / 4" by (rule DERIV_unique)
+ qed
+
+ have h'_zero: "h' z = 0" for z
+ proof -
+ def m \<equiv> "max 1 \<bar>Re z\<bar>"
+ def B \<equiv> "{t. abs (Re t) \<le> m \<and> abs (Im t) \<le> abs (Im z)}"
+ have "closed ({t. Re t \<ge> -m} \<inter> {t. Re t \<le> m} \<inter>
+ {t. Im t \<ge> -\<bar>Im z\<bar>} \<inter> {t. Im t \<le> \<bar>Im z\<bar>})"
+ (is "closed ?B") by (intro closed_Int closed_halfspace_Re_ge closed_halfspace_Re_le
+ closed_halfspace_Im_ge closed_halfspace_Im_le)
+ also have "?B = B" unfolding B_def by fastforce
+ finally have "closed B" .
+ moreover have "bounded B" unfolding bounded_iff
+ proof (intro ballI exI)
+ fix t assume t: "t \<in> B"
+ have "norm t \<le> \<bar>Re t\<bar> + \<bar>Im t\<bar>" by (rule cmod_le)
+ also from t have "\<bar>Re t\<bar> \<le> m" unfolding B_def by blast
+ also from t have "\<bar>Im t\<bar> \<le> \<bar>Im z\<bar>" unfolding B_def by blast
+ finally show "norm t \<le> m + \<bar>Im z\<bar>" by - simp
+ qed
+ ultimately have compact: "compact B" by (subst compact_eq_bounded_closed) blast
+
+ def M \<equiv> "SUP z:B. norm (h' z)"
+ have "compact (h' ` B)"
+ by (intro compact_continuous_image continuous_on_subset[OF h'_cont] compact) blast+
+ hence bdd: "bdd_above ((\<lambda>z. norm (h' z)) ` B)"
+ using bdd_above_norm[of "h' ` B"] by (simp add: image_comp o_def compact_imp_bounded)
+ have "norm (h' z) \<le> M" unfolding M_def by (intro cSUP_upper bdd) (simp_all add: B_def m_def)
+ also have "M \<le> M/2"
+ proof (subst M_def, subst cSUP_le_iff)
+ have "z \<in> B" unfolding B_def m_def by simp
+ thus "B \<noteq> {}" by auto
+ next
+ show "\<forall>z\<in>B. norm (h' z) \<le> M/2"
+ proof
+ fix t :: complex assume t: "t \<in> B"
+ from h'_eq[of t] t have "h' t = (h' (t/2) + h' ((t+1)/2)) / 4" by (simp add: dist_0_norm)
+ also have "norm \<dots> = norm (h' (t/2) + h' ((t+1)/2)) / 4" by simp
+ also have "norm (h' (t/2) + h' ((t+1)/2)) \<le> norm (h' (t/2)) + norm (h' ((t+1)/2))"
+ by (rule norm_triangle_ineq)
+ also from t have "abs (Re ((t + 1)/2)) \<le> m" unfolding m_def B_def by auto
+ with t have "t/2 \<in> B" "(t+1)/2 \<in> B" unfolding B_def by auto
+ hence "norm (h' (t/2)) + norm (h' ((t+1)/2)) \<le> M + M" unfolding M_def
+ by (intro add_mono cSUP_upper bdd) (auto simp: B_def)
+ also have "(M + M) / 4 = M / 2" by simp
+ finally show "norm (h' t) \<le> M/2" by - simp_all
+ qed
+ qed (insert bdd, auto simp: cball_eq_empty)
+ hence "M \<le> 0" by simp
+ finally show "h' z = 0" by simp
+ qed
+ have h_h'_2: "(h has_field_derivative 0) (at z)" for z
+ using h_h'[of z] h'_zero[of z] by simp
+
+ have g_real: "g z \<in> \<real>" if "z \<in> \<real>" for z
+ unfolding g_def using that by (auto intro!: Reals_mult Gamma_complex_real)
+ have h_real: "h z \<in> \<real>" if "z \<in> \<real>" for z
+ unfolding h_def using that by (auto intro!: Reals_mult Reals_add Reals_diff Polygamma_Real)
+ have g_nz: "g z \<noteq> 0" for z unfolding g_def using Ints_diff[of 1 "1-z"]
+ by (auto simp: Gamma_eq_zero_iff sin_eq_0)
+
+ from h'_zero h_h'_2 have "\<exists>c. \<forall>z\<in>UNIV. h z = c"
+ by (intro has_field_derivative_zero_constant) (simp_all add: dist_0_norm)
+ then obtain c where c: "\<And>z. h z = c" by auto
+ have "\<exists>u. u \<in> closed_segment 0 1 \<and> Re (g 1) - Re (g 0) = Re (h u * g u * (1 - 0))"
+ by (intro complex_mvt_line g_g')
+ find_theorems name:deriv Reals
+ then guess u by (elim exE conjE) note u = this
+ from u(1) have u': "u \<in> \<real>" unfolding closed_segment_def
+ by (auto simp: scaleR_conv_of_real)
+ from u' g_real[of u] g_nz[of u] have "Re (g u) \<noteq> 0" by (auto elim!: Reals_cases)
+ with u(2) c[of u] g_real[of u] g_nz[of u] u'
+ have "Re c = 0" by (simp add: complex_is_Real_iff g.of_1)
+ with h_real[of 0] c[of 0] have "c = 0" by (auto elim!: Reals_cases)
+ with c have A: "h z * g z = 0" for z by simp
+ hence "(g has_field_derivative 0) (at z)" for z using g_g'[of z] by simp
+ hence "\<exists>c'. \<forall>z\<in>UNIV. g z = c'" by (intro has_field_derivative_zero_constant) simp_all
+ then obtain c' where c: "\<And>z. g z = c'" by (force simp: dist_0_norm)
+ moreover from this[of 0] have "c' = pi" unfolding g_def by simp
+ ultimately have "g z = pi" by simp
+
+ show ?thesis
+ proof (cases "z \<in> \<int>")
+ case False
+ with `g z = pi` show ?thesis by (auto simp: g_def divide_simps)
+ next
+ case True
+ then obtain n where n: "z = of_int n" by (elim Ints_cases)
+ with sin_eq_0[of "of_real pi * z"] have "sin (of_real pi * z) = 0" by force
+ moreover have "of_int (1 - n) \<in> \<int>\<^sub>\<le>\<^sub>0" if "n > 0" using that by (intro nonpos_Ints_of_int) simp
+ ultimately show ?thesis using n
+ by (cases "n \<le> 0") (auto simp: Gamma_eq_zero_iff nonpos_Ints_of_int)
+ qed
+qed
+
+lemma rGamma_reflection_complex:
+ "rGamma z * rGamma (1 - z :: complex) = sin (of_real pi * z) / of_real pi"
+ using Gamma_reflection_complex[of z]
+ by (simp add: Gamma_def divide_simps split: split_if_asm)
+
+lemma rGamma_reflection_complex':
+ "rGamma z * rGamma (- z :: complex) = -z * sin (of_real pi * z) / of_real pi"
+proof -
+ have "rGamma z * rGamma (-z) = -z * (rGamma z * rGamma (1 - z))"
+ using rGamma_plus1[of "-z", symmetric] by simp
+ also have "rGamma z * rGamma (1 - z) = sin (of_real pi * z) / of_real pi"
+ by (rule rGamma_reflection_complex)
+ finally show ?thesis by simp
+qed
+
+lemma Gamma_reflection_complex':
+ "Gamma z * Gamma (- z :: complex) = - of_real pi / (z * sin (of_real pi * z))"
+ using rGamma_reflection_complex'[of z] by (force simp add: Gamma_def divide_simps mult_ac)
+
+
+
+lemma Gamma_one_half_real: "Gamma (1/2 :: real) = sqrt pi"
+proof -
+ from Gamma_reflection_complex[of "1/2"] fraction_not_in_ints[where 'a = complex, of 2 1]
+ have "Gamma (1/2 :: complex)^2 = of_real pi" by (simp add: power2_eq_square)
+ hence "of_real pi = Gamma (complex_of_real (1/2))^2" by simp
+ also have "\<dots> = of_real ((Gamma (1/2))^2)" by (subst Gamma_complex_of_real) simp_all
+ finally have "Gamma (1/2)^2 = pi" by (subst (asm) of_real_eq_iff) simp_all
+ moreover have "Gamma (1/2 :: real) \<ge> 0" using Gamma_real_pos[of "1/2"] by simp
+ ultimately show ?thesis by (rule real_sqrt_unique [symmetric])
+qed
+
+lemma Gamma_one_half_complex: "Gamma (1/2 :: complex) = of_real (sqrt pi)"
+proof -
+ have "Gamma (1/2 :: complex) = Gamma (of_real (1/2))" by simp
+ also have "\<dots> = of_real (sqrt pi)" by (simp only: Gamma_complex_of_real Gamma_one_half_real)
+ finally show ?thesis .
+qed
+
+lemma Gamma_legendre_duplication:
+ fixes z :: complex
+ assumes "z \<notin> \<int>\<^sub>\<le>\<^sub>0" "z + 1/2 \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ shows "Gamma z * Gamma (z + 1/2) =
+ exp ((1 - 2*z) * of_real (ln 2)) * of_real (sqrt pi) * Gamma (2*z)"
+ using Gamma_legendre_duplication_aux[OF assms] by (simp add: Gamma_one_half_complex)
+
+end
+
+
+subsection \<open>Limits and residues\<close>
+
+text \<open>
+ The inverse Gamma function has simple zeros:
+\<close>
+
+lemma rGamma_zeros:
+ "(\<lambda>z. rGamma z / (z + of_nat n)) \<midarrow> (- of_nat n) \<rightarrow> ((-1)^n * fact n :: 'a :: Gamma)"
+proof (subst tendsto_cong)
+ let ?f = "\<lambda>z. pochhammer z n * rGamma (z + of_nat (Suc n)) :: 'a"
+ from eventually_at_ball'[OF zero_less_one, of "- of_nat n :: 'a" UNIV]
+ show "eventually (\<lambda>z. rGamma z / (z + of_nat n) = ?f z) (at (- of_nat n))"
+ by (subst pochhammer_rGamma[of _ "Suc n"])
+ (auto elim!: eventually_mono simp: divide_simps pochhammer_rec' eq_neg_iff_add_eq_0)
+ have "isCont ?f (- of_nat n)" by (intro continuous_intros)
+ thus "?f \<midarrow> (- of_nat n) \<rightarrow> (- 1) ^ n * fact n" unfolding isCont_def
+ by (simp add: pochhammer_same)
+qed
+
+
+text \<open>
+ The simple zeros of the inverse Gamma function correspond to simple poles of the Gamma function,
+ and their residues can easily be computed from the limit we have just proven:
+\<close>
+
+lemma Gamma_poles: "filterlim Gamma at_infinity (at (- of_nat n :: 'a :: Gamma))"
+proof -
+ from eventually_at_ball'[OF zero_less_one, of "- of_nat n :: 'a" UNIV]
+ have "eventually (\<lambda>z. rGamma z \<noteq> (0 :: 'a)) (at (- of_nat n))"
+ by (auto elim!: eventually_mono nonpos_Ints_cases'
+ simp: rGamma_eq_zero_iff dist_of_nat dist_minus)
+ with isCont_rGamma[of "- of_nat n :: 'a", OF continuous_ident]
+ have "filterlim (\<lambda>z. inverse (rGamma z) :: 'a) at_infinity (at (- of_nat n))"
+ unfolding isCont_def by (intro filterlim_compose[OF filterlim_inverse_at_infinity])
+ (simp_all add: filterlim_at)
+ moreover have "(\<lambda>z. inverse (rGamma z) :: 'a) = Gamma"
+ by (intro ext) (simp add: rGamma_inverse_Gamma)
+ ultimately show ?thesis by (simp only: )
+qed
+
+lemma Gamma_residues:
+ "(\<lambda>z. Gamma z * (z + of_nat n)) \<midarrow> (- of_nat n) \<rightarrow> ((-1)^n / fact n :: 'a :: Gamma)"
+proof (subst tendsto_cong)
+ let ?c = "(- 1) ^ n / fact n :: 'a"
+ from eventually_at_ball'[OF zero_less_one, of "- of_nat n :: 'a" UNIV]
+ show "eventually (\<lambda>z. Gamma z * (z + of_nat n) = inverse (rGamma z / (z + of_nat n)))
+ (at (- of_nat n))"
+ by (auto elim!: eventually_mono simp: divide_simps rGamma_inverse_Gamma)
+ have "(\<lambda>z. inverse (rGamma z / (z + of_nat n))) \<midarrow> (- of_nat n) \<rightarrow>
+ inverse ((- 1) ^ n * fact n :: 'a)"
+ by (intro tendsto_intros rGamma_zeros) simp_all
+ also have "inverse ((- 1) ^ n * fact n) = ?c"
+ by (simp_all add: field_simps power_mult_distrib [symmetric] del: power_mult_distrib)
+ finally show "(\<lambda>z. inverse (rGamma z / (z + of_nat n))) \<midarrow> (- of_nat n) \<rightarrow> ?c" .
+qed
+
+
+
+subsection \<open>Alternative definitions\<close>
+
+
+subsubsection \<open>Variant of the Euler form\<close>
+
+
+definition Gamma_series_euler' where
+ "Gamma_series_euler' z n =
+ inverse z * (\<Prod>k=1..n. exp (z * of_real (ln (1 + inverse (of_nat k)))) / (1 + z / of_nat k))"
+
+context
+begin
+private lemma Gamma_euler'_aux1:
+ fixes z :: "'a :: {real_normed_field,banach}"
+ assumes n: "n > 0"
+ shows "exp (z * of_real (ln (of_nat n + 1))) = (\<Prod>k=1..n. exp (z * of_real (ln (1 + 1 / of_nat k))))"
+proof -
+ have "(\<Prod>k=1..n. exp (z * of_real (ln (1 + 1 / of_nat k)))) =
+ exp (z * of_real (\<Sum>k = 1..n. ln (1 + 1 / real_of_nat k)))"
+ by (subst exp_setsum [symmetric]) (simp_all add: setsum_right_distrib)
+ also have "(\<Sum>k=1..n. ln (1 + 1 / of_nat k) :: real) = ln (\<Prod>k=1..n. 1 + 1 / real_of_nat k)"
+ by (subst ln_setprod [symmetric]) (auto intro!: add_pos_nonneg)
+ also have "(\<Prod>k=1..n. 1 + 1 / of_nat k :: real) = (\<Prod>k=1..n. (of_nat k + 1) / of_nat k)"
+ by (intro setprod.cong) (simp_all add: divide_simps)
+ also have "(\<Prod>k=1..n. (of_nat k + 1) / of_nat k :: real) = of_nat n + 1"
+ by (induction n) (simp_all add: setprod_nat_ivl_Suc' divide_simps)
+ finally show ?thesis ..
+qed
+
+lemma Gamma_series_euler':
+ assumes z: "(z :: 'a :: Gamma) \<notin> \<int>\<^sub>\<le>\<^sub>0"
+ shows "(\<lambda>n. Gamma_series_euler' z n) \<longlonglongrightarrow> Gamma z"
+proof (rule Gamma_seriesI, rule Lim_transform_eventually)
+ let ?f = "\<lambda>n. fact n * exp (z * of_real (ln (of_nat n + 1))) / pochhammer z (n + 1)"
+ let ?r = "\<lambda>n. ?f n / Gamma_series z n"
+ let ?r' = "\<lambda>n. exp (z * of_real (ln (of_nat (Suc n) / of_nat n)))"
+ from z have z': "z \<noteq> 0" by auto
+
+ have "eventually (\<lambda>n. ?r' n = ?r n) sequentially" using eventually_gt_at_top[of "0::nat"]
+ using z by (auto simp: divide_simps Gamma_series_def ring_distribs exp_diff ln_div add_ac
+ elim!: eventually_mono dest: pochhammer_eq_0_imp_nonpos_Int)
+ moreover have "?r' \<longlonglongrightarrow> exp (z * of_real (ln 1))"
+ by (intro tendsto_intros LIMSEQ_Suc_n_over_n) simp_all
+ ultimately show "?r \<longlonglongrightarrow> 1" by (force dest!: Lim_transform_eventually)
+
+ from eventually_gt_at_top[of "0::nat"]
+ show "eventually (\<lambda>n. ?r n = Gamma_series_euler' z n / Gamma_series z n) sequentially"
+ proof eventually_elim
+ fix n :: nat assume n: "n > 0"
+ from n z' have "Gamma_series_euler' z n =
+ exp (z * of_real (ln (of_nat n + 1))) / (z * (\<Prod>k=1..n. (1 + z / of_nat k)))"
+ by (subst Gamma_euler'_aux1)
+ (simp_all add: Gamma_series_euler'_def setprod.distrib
+ setprod_inversef[symmetric] divide_inverse)
+ also have "(\<Prod>k=1..n. (1 + z / of_nat k)) = pochhammer (z + 1) n / fact n"
+ by (cases n) (simp_all add: pochhammer_def fact_altdef setprod_shift_bounds_cl_Suc_ivl
+ setprod_dividef[symmetric] divide_simps add_ac)
+ also have "z * \<dots> = pochhammer z (Suc n) / fact n" by (simp add: pochhammer_rec)
+ finally show "?r n = Gamma_series_euler' z n / Gamma_series z n" by simp
+ qed
+qed
+
+end
+
+
+
+subsubsection \<open>Weierstrass form\<close>
+
+definition Gamma_series_weierstrass :: "'a :: {banach,real_normed_field} \<Rightarrow> nat \<Rightarrow> 'a" where
+ "Gamma_series_weierstrass z n =
+ exp (-euler_mascheroni * z) / z * (\<Prod>k=1..n. exp (z / of_nat k) / (1 + z / of_nat k))"
+
+definition rGamma_series_weierstrass :: "'a :: {banach,real_normed_field} \<Rightarrow> nat \<Rightarrow> 'a" where
+ "rGamma_series_weierstrass z n =
+ exp (euler_mascheroni * z) * z * (\<Prod>k=1..n. (1 + z / of_nat k) * exp (-z / of_nat k))"
+
+lemma Gamma_series_weierstrass_nonpos_Ints:
+ "eventually (\<lambda>k. Gamma_series_weierstrass (- of_nat n) k = 0) sequentially"
+ using eventually_ge_at_top[of n] by eventually_elim (auto simp: Gamma_series_weierstrass_def)
+
+lemma rGamma_series_weierstrass_nonpos_Ints:
+ "eventually (\<lambda>k. rGamma_series_weierstrass (- of_nat n) k = 0) sequentially"
+ using eventually_ge_at_top[of n] by eventually_elim (auto simp: rGamma_series_weierstrass_def)
+
+lemma Gamma_weierstrass_complex: "Gamma_series_weierstrass z \<longlonglongrightarrow> Gamma (z :: complex)"
+proof (cases "z \<in> \<int>\<^sub>\<le>\<^sub>0")
+ case True
+ then obtain n where "z = - of_nat n" by (elim nonpos_Ints_cases')
+ also from True have "Gamma_series_weierstrass \<dots> \<longlonglongrightarrow> Gamma z"
+ by (simp add: tendsto_cong[OF Gamma_series_weierstrass_nonpos_Ints] Gamma_nonpos_Int)
+ finally show ?thesis .
+next
+ case False
+ hence z: "z \<noteq> 0" by auto
+ let ?f = "(\<lambda>x. \<Prod>x = Suc 0..x. exp (z / of_nat x) / (1 + z / of_nat x))"
+ have A: "exp (ln (1 + z / of_nat n)) = (1 + z / of_nat n)" if "n \<ge> 1" for n :: nat
+ using False that by (subst exp_Ln) (auto simp: field_simps dest!: plus_of_nat_eq_0_imp)
+ have "(\<lambda>n. \<Sum>k=1..n. z / of_nat k - ln (1 + z / of_nat k)) \<longlonglongrightarrow> ln_Gamma z + euler_mascheroni * z + ln z"
+ using ln_Gamma_series'_aux[OF False]
+ by (simp only: atLeastLessThanSuc_atLeastAtMost [symmetric] One_nat_def
+ setsum_shift_bounds_Suc_ivl sums_def atLeast0LessThan)
+ from tendsto_exp[OF this] False z have "?f \<longlonglongrightarrow> z * exp (euler_mascheroni * z) * Gamma z"
+ by (simp add: exp_add exp_setsum exp_diff mult_ac Gamma_complex_altdef A)
+ from tendsto_mult[OF tendsto_const[of "exp (-euler_mascheroni * z) / z"] this] z
+ show "Gamma_series_weierstrass z \<longlonglongrightarrow> Gamma z"
+ by (simp add: exp_minus divide_simps Gamma_series_weierstrass_def [abs_def])
+qed
+
+lemma Gamma_weierstrass_real: "Gamma_series_weierstrass x \<longlonglongrightarrow> Gamma (x :: real)"
+ using Gamma_weierstrass_complex[of "of_real x"] unfolding Gamma_series_weierstrass_def[abs_def]
+ by (subst tendsto_complex_of_real_iff [symmetric])
+ (simp_all add: exp_of_real[symmetric] Gamma_complex_of_real)
+
+lemma rGamma_weierstrass_complex: "rGamma_series_weierstrass z \<longlonglongrightarrow> rGamma (z :: complex)"
+proof (cases "z \<in> \<int>\<^sub>\<le>\<^sub>0")
+ case True
+ then obtain n where "z = - of_nat n" by (elim nonpos_Ints_cases')
+ also from True have "rGamma_series_weierstrass \<dots> \<longlonglongrightarrow> rGamma z"
+ by (simp add: tendsto_cong[OF rGamma_series_weierstrass_nonpos_Ints] rGamma_nonpos_Int)
+ finally show ?thesis .
+next
+ case False
+ have "rGamma_series_weierstrass z = (\<lambda>n. inverse (Gamma_series_weierstrass z n))"
+ by (simp add: rGamma_series_weierstrass_def[abs_def] Gamma_series_weierstrass_def
+ exp_minus divide_inverse setprod_inversef[symmetric] mult_ac)
+ also from False have "\<dots> \<longlonglongrightarrow> inverse (Gamma z)"
+ by (intro tendsto_intros Gamma_weierstrass_complex) (simp add: Gamma_eq_zero_iff)
+ finally show ?thesis by (simp add: Gamma_def)
+qed
+
+lemma sin_product_formula_complex:
+ fixes z :: complex
+ shows "(\<lambda>n. of_real pi * z * (\<Prod>k=1..n. 1 - z^2 / of_nat k^2)) \<longlonglongrightarrow> sin (of_real pi * z)"
+proof -
+ let ?f = "rGamma_series_weierstrass"
+ have "(\<lambda>n. (- of_real pi * inverse z) * (?f z n * ?f (- z) n))
+ \<longlonglongrightarrow> (- of_real pi * inverse z) * (rGamma z * rGamma (- z))"
+ by (intro tendsto_intros rGamma_weierstrass_complex)
+ also have "(\<lambda>n. (- of_real pi * inverse z) * (?f z n * ?f (-z) n)) =
+ (\<lambda>n. of_real pi * z * (\<Prod>k=1..n. 1 - z^2 / of_nat k ^ 2))"
+ proof
+ fix n :: nat
+ have "(- of_real pi * inverse z) * (?f z n * ?f (-z) n) =
+ of_real pi * z * (\<Prod>k=1..n. (of_nat k - z) * (of_nat k + z) / of_nat k ^ 2)"
+ by (simp add: rGamma_series_weierstrass_def mult_ac exp_minus
+ divide_simps setprod.distrib[symmetric] power2_eq_square)
+ also have "(\<Prod>k=1..n. (of_nat k - z) * (of_nat k + z) / of_nat k ^ 2) =
+ (\<Prod>k=1..n. 1 - z^2 / of_nat k ^ 2)"
+ by (intro setprod.cong) (simp_all add: power2_eq_square field_simps)
+ finally show "(- of_real pi * inverse z) * (?f z n * ?f (-z) n) = of_real pi * z * \<dots>"
+ by (simp add: divide_simps)
+ qed
+ also have "(- of_real pi * inverse z) * (rGamma z * rGamma (- z)) = sin (of_real pi * z)"
+ by (subst rGamma_reflection_complex') (simp add: divide_simps)
+ finally show ?thesis .
+qed
+
+lemma sin_product_formula_real:
+ "(\<lambda>n. pi * (x::real) * (\<Prod>k=1..n. 1 - x^2 / of_nat k^2)) \<longlonglongrightarrow> sin (pi * x)"
+proof -
+ from sin_product_formula_complex[of "of_real x"]
+ have "(\<lambda>n. of_real pi * of_real x * (\<Prod>k=1..n. 1 - (of_real x)^2 / (of_nat k)^2))
+ \<longlonglongrightarrow> sin (of_real pi * of_real x :: complex)" (is "?f \<longlonglongrightarrow> ?y") .
+ also have "?f = (\<lambda>n. of_real (pi * x * (\<Prod>k=1..n. 1 - x^2 / (of_nat k^2))))" by simp
+ also have "?y = of_real (sin (pi * x))" by (simp only: sin_of_real [symmetric] of_real_mult)
+ finally show ?thesis by (subst (asm) tendsto_of_real_iff)
+qed
+
+lemma sin_product_formula_real':
+ assumes "x \<noteq> (0::real)"
+ shows "(\<lambda>n. (\<Prod>k=1..n. 1 - x^2 / of_nat k^2)) \<longlonglongrightarrow> sin (pi * x) / (pi * x)"
+ using tendsto_divide[OF sin_product_formula_real[of x] tendsto_const[of "pi * x"]] assms
+ by simp
+
+theorem inverse_squares_sums: "(\<lambda>n. 1 / (n + 1)\<^sup>2) sums (pi\<^sup>2 / 6)"
+proof -
+ def P \<equiv> "\<lambda>x n. (\<Prod>k=1..n. 1 - x^2 / of_nat k^2 :: real)"
+ def K \<equiv> "\<Sum>n. inverse (real_of_nat (Suc n))^2"
+ def f \<equiv> "\<lambda>x. \<Sum>n. P x n / of_nat (Suc n)^2"
+ def g \<equiv> "\<lambda>x. (1 - sin (pi * x) / (pi * x))"
+
+ have sums: "(\<lambda>n. P x n / of_nat (Suc n)^2) sums (if x = 0 then K else g x / x^2)" for x
+ proof (cases "x = 0")
+ assume x: "x = 0"
+ have "summable (\<lambda>n. inverse ((real_of_nat (Suc n))\<^sup>2))"
+ using inverse_power_summable[of 2] by (subst summable_Suc_iff) simp
+ thus ?thesis by (simp add: x g_def P_def K_def inverse_eq_divide power_divide summable_sums)
+ next
+ assume x: "x \<noteq> 0"
+ have "(\<lambda>n. P x n - P x (Suc n)) sums (P x 0 - sin (pi * x) / (pi * x))"
+ unfolding P_def using x by (intro telescope_sums' sin_product_formula_real')
+ also have "(\<lambda>n. P x n - P x (Suc n)) = (\<lambda>n. (x^2 / of_nat (Suc n)^2) * P x n)"
+ unfolding P_def by (simp add: setprod_nat_ivl_Suc' algebra_simps)
+ also have "P x 0 = 1" by (simp add: P_def)
+ finally have "(\<lambda>n. x\<^sup>2 / (of_nat (Suc n))\<^sup>2 * P x n) sums (1 - sin (pi * x) / (pi * x))" .
+ from sums_divide[OF this, of "x^2"] x show ?thesis unfolding g_def by simp
+ qed
+
+ have "continuous_on (ball 0 1) f"
+ proof (rule uniform_limit_theorem; (intro always_eventually allI)?)
+ show "uniform_limit (ball 0 1) (\<lambda>n x. \<Sum>k<n. P x k / of_nat (Suc k)^2) f sequentially"
+ proof (unfold f_def, rule weierstrass_m_test)
+ fix n :: nat and x :: real assume x: "x \<in> ball 0 1"
+ {
+ fix k :: nat assume k: "k \<ge> 1"
+ from x have "x^2 < 1" by (auto simp: dist_0_norm abs_square_less_1)
+ also from k have "\<dots> \<le> of_nat k^2" by simp
+ finally have "(1 - x^2 / of_nat k^2) \<in> {0..1}" using k
+ by (simp_all add: field_simps del: of_nat_Suc)
+ }
+ hence "(\<Prod>k=1..n. abs (1 - x^2 / of_nat k^2)) \<le> (\<Prod>k=1..n. 1)" by (intro setprod_mono) simp
+ thus "norm (P x n / (of_nat (Suc n)^2)) \<le> 1 / of_nat (Suc n)^2"
+ unfolding P_def by (simp add: field_simps abs_setprod del: of_nat_Suc)
+ qed (subst summable_Suc_iff, insert inverse_power_summable[of 2], simp add: inverse_eq_divide)
+ qed (auto simp: P_def intro!: continuous_intros)
+ hence "isCont f 0" by (subst (asm) continuous_on_eq_continuous_at) simp_all
+ hence "(f \<midarrow> 0 \<rightarrow> f 0)" by (simp add: isCont_def)
+ also have "f 0 = K" unfolding f_def P_def K_def by (simp add: inverse_eq_divide power_divide)
+ finally have "f \<midarrow> 0 \<rightarrow> K" .
+
+ moreover have "f \<midarrow> 0 \<rightarrow> pi^2 / 6"
+ proof (rule Lim_transform_eventually)
+ def f' \<equiv> "\<lambda>x. \<Sum>n. - sin_coeff (n+3) * pi ^ (n+2) * x^n"
+ have "eventually (\<lambda>x. x \<noteq> (0::real)) (at 0)"
+ by (auto simp add: eventually_at intro!: exI[of _ 1])
+ thus "eventually (\<lambda>x. f' x = f x) (at 0)"
+ proof eventually_elim
+ fix x :: real assume x: "x \<noteq> 0"
+ have "sin_coeff 1 = (1 :: real)" "sin_coeff 2 = (0::real)" by (simp_all add: sin_coeff_def)
+ with sums_split_initial_segment[OF sums_minus[OF sin_converges], of 3 "pi*x"]
+ have "(\<lambda>n. - (sin_coeff (n+3) * (pi*x)^(n+3))) sums (pi * x - sin (pi*x))"
+ by (simp add: eval_nat_numeral)
+ from sums_divide[OF this, of "x^3 * pi"] x
+ have "(\<lambda>n. - (sin_coeff (n+3) * pi^(n+2) * x^n)) sums ((1 - sin (pi*x) / (pi*x)) / x^2)"
+ by (simp add: divide_simps eval_nat_numeral power_mult_distrib mult_ac)
+ with x have "(\<lambda>n. - (sin_coeff (n+3) * pi^(n+2) * x^n)) sums (g x / x^2)"
+ by (simp add: g_def)
+ hence "f' x = g x / x^2" by (simp add: sums_iff f'_def)
+ also have "\<dots> = f x" using sums[of x] x by (simp add: sums_iff g_def f_def)
+ finally show "f' x = f x" .
+ qed
+
+ have "isCont f' 0" unfolding f'_def
+ proof (intro isCont_powser_converges_everywhere)
+ fix x :: real show "summable (\<lambda>n. -sin_coeff (n+3) * pi^(n+2) * x^n)"
+ proof (cases "x = 0")
+ assume x: "x \<noteq> 0"
+ from summable_divide[OF sums_summable[OF sums_split_initial_segment[OF
+ sin_converges[of "pi*x"]], of 3], of "-pi*x^3"] x
+ show ?thesis by (simp add: mult_ac power_mult_distrib divide_simps eval_nat_numeral)
+ qed (simp only: summable_0_powser)
+ qed
+ hence "f' \<midarrow> 0 \<rightarrow> f' 0" by (simp add: isCont_def)
+ also have "f' 0 = pi * pi / fact 3" unfolding f'_def
+ by (subst powser_zero) (simp add: sin_coeff_def)
+ finally show "f' \<midarrow> 0 \<rightarrow> pi^2 / 6" by (simp add: eval_nat_numeral)
+ qed
+
+ ultimately have "K = pi^2 / 6" by (rule LIM_unique)
+ moreover from inverse_power_summable[of 2]
+ have "summable (\<lambda>n. (inverse (real_of_nat (Suc n)))\<^sup>2)"
+ by (subst summable_Suc_iff) (simp add: power_inverse)
+ ultimately show ?thesis unfolding K_def
+ by (auto simp add: sums_iff power_divide inverse_eq_divide)
+qed
+
+
+
+subsection \<open>Binomial coefficient form\<close>
+
+lemma Gamma_binomial:
+ "(\<lambda>n. ((z + of_nat n) gchoose n) * exp (-z * of_real (ln (of_nat n)))) \<longlonglongrightarrow> rGamma (z+1)"
+proof (cases "z = 0")
+ case False
+ show ?thesis
+ proof (rule Lim_transform_eventually)
+ let ?powr = "\<lambda>a b. exp (b * of_real (ln (of_nat a)))"
+ show "eventually (\<lambda>n. rGamma_series z n / z =
+ ((z + of_nat n) gchoose n) * ?powr n (-z)) sequentially"
+ proof (intro always_eventually allI)
+ fix n :: nat
+ from False have "((z + of_nat n) gchoose n) = pochhammer z (Suc n) / z / fact n"
+ by (simp add: gbinomial_pochhammer' pochhammer_rec)
+ also have "pochhammer z (Suc n) / z / fact n * ?powr n (-z) = rGamma_series z n / z"
+ by (simp add: rGamma_series_def divide_simps exp_minus)
+ finally show "rGamma_series z n / z = ((z + of_nat n) gchoose n) * ?powr n (-z)" ..
+ qed
+
+ from False have "(\<lambda>n. rGamma_series z n / z) \<longlonglongrightarrow> rGamma z / z" by (intro tendsto_intros)
+ also from False have "rGamma z / z = rGamma (z + 1)" using rGamma_plus1[of z]
+ by (simp add: field_simps)
+ finally show "(\<lambda>n. rGamma_series z n / z) \<longlonglongrightarrow> rGamma (z+1)" .
+ qed
+qed (simp_all add: binomial_gbinomial [symmetric])
+
+lemma fact_binomial_limit:
+ "(\<lambda>n. of_nat ((k + n) choose n) / of_nat (n ^ k) :: 'a :: Gamma) \<longlonglongrightarrow> 1 / fact k"
+proof (rule Lim_transform_eventually)
+ have "(\<lambda>n. of_nat ((k + n) choose n) / of_real (exp (of_nat k * ln (real_of_nat n))))
+ \<longlonglongrightarrow> 1 / Gamma (of_nat (Suc k) :: 'a)" (is "?f \<longlonglongrightarrow> _")
+ using Gamma_binomial[of "of_nat k :: 'a"]
+ by (simp add: binomial_gbinomial add_ac Gamma_def divide_simps exp_of_real [symmetric] exp_minus)
+ also have "Gamma (of_nat (Suc k)) = fact k" by (rule Gamma_fact)
+ finally show "?f \<longlonglongrightarrow> 1 / fact k" .
+
+ show "eventually (\<lambda>n. ?f n = of_nat ((k + n) choose n) / of_nat (n ^ k)) sequentially"
+ using eventually_gt_at_top[of "0::nat"]
+ proof eventually_elim
+ fix n :: nat assume n: "n > 0"
+ from n have "exp (real_of_nat k * ln (real_of_nat n)) = real_of_nat (n^k)"
+ by (simp add: exp_of_nat_mult)
+ thus "?f n = of_nat ((k + n) choose n) / of_nat (n ^ k)" by simp
+ qed
+qed
+
+lemma binomial_asymptotic:
+ "(\<lambda>n. of_nat ((k + n) choose n) / (of_nat (n ^ k) / fact k) :: 'a :: Gamma) \<longlonglongrightarrow> 1"
+ using tendsto_mult[OF fact_binomial_limit[of k] tendsto_const[of "fact k :: 'a"]] by simp
+
+end