src/HOL/Multivariate_Analysis/Gamma.thy
author hoelzl
Fri Feb 19 13:40:50 2016 +0100 (2016-02-19)
changeset 62378 85ed00c1fe7c
parent 62131 1baed43f453e
child 62390 842917225d56
child 62397 5ae24f33d343
permissions -rw-r--r--
generalize more theorems to support enat and ennreal
     1 (*  Title:    HOL/Multivariate_Analysis/Gamma.thy
     2     Author:   Manuel Eberl, TU München
     3 *)
     4 
     5 section \<open>The Gamma Function\<close>
     6 
     7 theory Gamma
     8 imports
     9   Complex_Transcendental
    10   Summation
    11   Harmonic_Numbers
    12   "~~/src/HOL/Library/Nonpos_Ints"
    13   "~~/src/HOL/Library/Periodic_Fun"
    14 begin
    15 
    16 text \<open>
    17   Several equivalent definitions of the Gamma function and its
    18   most important properties. Also contains the definition and some properties
    19   of the log-Gamma function and the Digamma function and the other Polygamma functions.
    20 
    21   Based on the Gamma function, we also prove the Weierstraß product form of the
    22   sin function and, based on this, the solution of the Basel problem (the
    23   sum over all @{term "1 / (n::nat)^2"}.
    24 \<close>
    25 
    26 lemma pochhammer_eq_0_imp_nonpos_Int:
    27   "pochhammer (x::'a::field_char_0) n = 0 \<Longrightarrow> x \<in> \<int>\<^sub>\<le>\<^sub>0"
    28   by (auto simp: pochhammer_eq_0_iff)
    29 
    30 lemma closed_nonpos_Ints [simp]: "closed (\<int>\<^sub>\<le>\<^sub>0 :: 'a :: real_normed_algebra_1 set)"
    31 proof -
    32   have "\<int>\<^sub>\<le>\<^sub>0 = (of_int ` {n. n \<le> 0} :: 'a set)"
    33     by (auto elim!: nonpos_Ints_cases intro!: nonpos_Ints_of_int)
    34   also have "closed \<dots>" by (rule closed_of_int_image)
    35   finally show ?thesis .
    36 qed
    37 
    38 lemma plus_one_in_nonpos_Ints_imp: "z + 1 \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> z \<in> \<int>\<^sub>\<le>\<^sub>0"
    39   using nonpos_Ints_diff_Nats[of "z+1" "1"] by simp_all
    40 
    41 lemma fraction_not_in_ints:
    42   assumes "\<not>(n dvd m)" "n \<noteq> 0"
    43   shows   "of_int m / of_int n \<notin> (\<int> :: 'a :: {division_ring,ring_char_0} set)"
    44 proof
    45   assume "of_int m / (of_int n :: 'a) \<in> \<int>"
    46   then obtain k where "of_int m / of_int n = (of_int k :: 'a)" by (elim Ints_cases)
    47   with assms have "of_int m = (of_int (k * n) :: 'a)" by (auto simp add: divide_simps)
    48   hence "m = k * n" by (subst (asm) of_int_eq_iff)
    49   hence "n dvd m" by simp
    50   with assms(1) show False by contradiction
    51 qed
    52 
    53 lemma not_in_Ints_imp_not_in_nonpos_Ints: "z \<notin> \<int> \<Longrightarrow> z \<notin> \<int>\<^sub>\<le>\<^sub>0"
    54   by (auto simp: Ints_def nonpos_Ints_def)
    55 
    56 lemma double_in_nonpos_Ints_imp:
    57   assumes "2 * (z :: 'a :: field_char_0) \<in> \<int>\<^sub>\<le>\<^sub>0"
    58   shows   "z \<in> \<int>\<^sub>\<le>\<^sub>0 \<or> z + 1/2 \<in> \<int>\<^sub>\<le>\<^sub>0"
    59 proof-
    60   from assms obtain k where k: "2 * z = - of_nat k" by (elim nonpos_Ints_cases')
    61   thus ?thesis by (cases "even k") (auto elim!: evenE oddE simp: field_simps)
    62 qed
    63 
    64 
    65 lemma sin_series: "(\<lambda>n. ((-1)^n / fact (2*n+1)) *\<^sub>R z^(2*n+1)) sums sin z"
    66 proof -
    67   from sin_converges[of z] have "(\<lambda>n. sin_coeff n *\<^sub>R z^n) sums sin z" .
    68   also have "(\<lambda>n. sin_coeff n *\<^sub>R z^n) sums sin z \<longleftrightarrow>
    69                  (\<lambda>n. ((-1)^n / fact (2*n+1)) *\<^sub>R z^(2*n+1)) sums sin z"
    70     by (subst sums_mono_reindex[of "\<lambda>n. 2*n+1", symmetric])
    71        (auto simp: sin_coeff_def subseq_def ac_simps elim!: oddE)
    72   finally show ?thesis .
    73 qed
    74 
    75 lemma cos_series: "(\<lambda>n. ((-1)^n / fact (2*n)) *\<^sub>R z^(2*n)) sums cos z"
    76 proof -
    77   from cos_converges[of z] have "(\<lambda>n. cos_coeff n *\<^sub>R z^n) sums cos z" .
    78   also have "(\<lambda>n. cos_coeff n *\<^sub>R z^n) sums cos z \<longleftrightarrow>
    79                  (\<lambda>n. ((-1)^n / fact (2*n)) *\<^sub>R z^(2*n)) sums cos z"
    80     by (subst sums_mono_reindex[of "\<lambda>n. 2*n", symmetric])
    81        (auto simp: cos_coeff_def subseq_def ac_simps elim!: evenE)
    82   finally show ?thesis .
    83 qed
    84 
    85 lemma sin_z_over_z_series:
    86   fixes z :: "'a :: {real_normed_field,banach}"
    87   assumes "z \<noteq> 0"
    88   shows   "(\<lambda>n. (-1)^n / fact (2*n+1) * z^(2*n)) sums (sin z / z)"
    89 proof -
    90   from sin_series[of z] have "(\<lambda>n. z * ((-1)^n / fact (2*n+1)) * z^(2*n)) sums sin z"
    91     by (simp add: field_simps scaleR_conv_of_real)
    92   from sums_mult[OF this, of "inverse z"] and assms show ?thesis
    93     by (simp add: field_simps)
    94 qed
    95 
    96 lemma sin_z_over_z_series':
    97   fixes z :: "'a :: {real_normed_field,banach}"
    98   assumes "z \<noteq> 0"
    99   shows   "(\<lambda>n. sin_coeff (n+1) *\<^sub>R z^n) sums (sin z / z)"
   100 proof -
   101   from sums_split_initial_segment[OF sin_converges[of z], of 1]
   102     have "(\<lambda>n. z * (sin_coeff (n+1) *\<^sub>R z ^ n)) sums sin z" by simp
   103   from sums_mult[OF this, of "inverse z"] assms show ?thesis by (simp add: field_simps)
   104 qed
   105 
   106 lemma has_field_derivative_sin_z_over_z:
   107   fixes A :: "'a :: {real_normed_field,banach} set"
   108   shows "((\<lambda>z. if z = 0 then 1 else sin z / z) has_field_derivative 0) (at 0 within A)"
   109       (is "(?f has_field_derivative ?f') _")
   110 proof (rule has_field_derivative_at_within)
   111   have "((\<lambda>z::'a. \<Sum>n. of_real (sin_coeff (n+1)) * z^n)
   112             has_field_derivative (\<Sum>n. diffs (\<lambda>n. of_real (sin_coeff (n+1))) n * 0^n)) (at 0)"
   113   proof (rule termdiffs_strong)
   114     from summable_ignore_initial_segment[OF sums_summable[OF sin_converges[of "1::'a"]], of 1]
   115       show "summable (\<lambda>n. of_real (sin_coeff (n+1)) * (1::'a)^n)" by (simp add: of_real_def)
   116   qed simp
   117   also have "(\<lambda>z::'a. \<Sum>n. of_real (sin_coeff (n+1)) * z^n) = ?f"
   118   proof
   119     fix z
   120     show "(\<Sum>n. of_real (sin_coeff (n+1)) * z^n)  = ?f z"
   121       by (cases "z = 0") (insert sin_z_over_z_series'[of z],
   122             simp_all add: scaleR_conv_of_real sums_iff powser_zero sin_coeff_def)
   123   qed
   124   also have "(\<Sum>n. diffs (\<lambda>n. of_real (sin_coeff (n + 1))) n * (0::'a) ^ n) =
   125                  diffs (\<lambda>n. of_real (sin_coeff (Suc n))) 0" by (simp add: powser_zero)
   126   also have "\<dots> = 0" by (simp add: sin_coeff_def diffs_def)
   127   finally show "((\<lambda>z::'a. if z = 0 then 1 else sin z / z) has_field_derivative 0) (at 0)" .
   128 qed
   129 
   130 lemma round_Re_minimises_norm:
   131   "norm ((z::complex) - of_int m) \<ge> norm (z - of_int (round (Re z)))"
   132 proof -
   133   let ?n = "round (Re z)"
   134   have "norm (z - of_int ?n) = sqrt ((Re z - of_int ?n)\<^sup>2 + (Im z)\<^sup>2)"
   135     by (simp add: cmod_def)
   136   also have "\<bar>Re z - of_int ?n\<bar> \<le> \<bar>Re z - of_int m\<bar>" by (rule round_diff_minimal)
   137   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)"
   138     by (intro real_sqrt_le_mono add_mono) (simp_all add: abs_le_square_iff)
   139   also have "\<dots> = norm (z - of_int m)" by (simp add: cmod_def)
   140   finally show ?thesis .
   141 qed
   142 
   143 lemma Re_pos_in_ball:
   144   assumes "Re z > 0" "t \<in> ball z (Re z/2)"
   145   shows   "Re t > 0"
   146 proof -
   147   have "Re (z - t) \<le> norm (z - t)" by (rule complex_Re_le_cmod)
   148   also from assms have "\<dots> < Re z / 2" by (simp add: dist_complex_def)
   149   finally show "Re t > 0" using assms by simp
   150 qed
   151 
   152 lemma no_nonpos_Int_in_ball_complex:
   153   assumes "Re z > 0" "t \<in> ball z (Re z/2)"
   154   shows   "t \<notin> \<int>\<^sub>\<le>\<^sub>0"
   155   using Re_pos_in_ball[OF assms] by (force elim!: nonpos_Ints_cases)
   156 
   157 lemma no_nonpos_Int_in_ball:
   158   assumes "t \<in> ball z (dist z (round (Re z)))"
   159   shows   "t \<notin> \<int>\<^sub>\<le>\<^sub>0"
   160 proof
   161   assume "t \<in> \<int>\<^sub>\<le>\<^sub>0"
   162   then obtain n where "t = of_int n" by (auto elim!: nonpos_Ints_cases)
   163   have "dist z (of_int n) \<le> dist z t + dist t (of_int n)" by (rule dist_triangle)
   164   also from assms have "dist z t < dist z (round (Re z))" by simp
   165   also have "\<dots> \<le> dist z (of_int n)"
   166     using round_Re_minimises_norm[of z] by (simp add: dist_complex_def)
   167   finally have "dist t (of_int n) > 0" by simp
   168   with \<open>t = of_int n\<close> show False by simp
   169 qed
   170 
   171 lemma no_nonpos_Int_in_ball':
   172   assumes "(z :: 'a :: {euclidean_space,real_normed_algebra_1}) \<notin> \<int>\<^sub>\<le>\<^sub>0"
   173   obtains d where "d > 0" "\<And>t. t \<in> ball z d \<Longrightarrow> t \<notin> \<int>\<^sub>\<le>\<^sub>0"
   174 proof (rule that)
   175   from assms show "setdist {z} \<int>\<^sub>\<le>\<^sub>0 > 0" by (subst setdist_gt_0_compact_closed) auto
   176 next
   177   fix t assume "t \<in> ball z (setdist {z} \<int>\<^sub>\<le>\<^sub>0)"
   178   thus "t \<notin> \<int>\<^sub>\<le>\<^sub>0" using setdist_le_dist[of z "{z}" t "\<int>\<^sub>\<le>\<^sub>0"] by force
   179 qed
   180 
   181 lemma no_nonpos_Real_in_ball:
   182   assumes z: "z \<notin> \<real>\<^sub>\<le>\<^sub>0" and t: "t \<in> ball z (if Im z = 0 then Re z / 2 else abs (Im z) / 2)"
   183   shows   "t \<notin> \<real>\<^sub>\<le>\<^sub>0"
   184 using z
   185 proof (cases "Im z = 0")
   186   assume A: "Im z = 0"
   187   with z have "Re z > 0" by (force simp add: complex_nonpos_Reals_iff)
   188   with t A Re_pos_in_ball[of z t] show ?thesis by (force simp add: complex_nonpos_Reals_iff)
   189 next
   190   assume A: "Im z \<noteq> 0"
   191   have "abs (Im z) - abs (Im t) \<le> abs (Im z - Im t)" by linarith
   192   also have "\<dots> = abs (Im (z - t))" by simp
   193   also have "\<dots> \<le> norm (z - t)" by (rule abs_Im_le_cmod)
   194   also from A t have "\<dots> \<le> abs (Im z) / 2" by (simp add: dist_complex_def)
   195   finally have "abs (Im t) > 0" using A by simp
   196   thus ?thesis by (force simp add: complex_nonpos_Reals_iff)
   197 qed
   198 
   199 
   200 subsection \<open>Definitions\<close>
   201 
   202 text \<open>
   203   We define the Gamma function by first defining its multiplicative inverse @{term "Gamma_inv"}.
   204   This is more convenient because @{term "Gamma_inv"} is entire, which makes proofs of its
   205   properties more convenient because one does not have to watch out for discontinuities.
   206   (e.g. @{term "Gamma_inv"} fulfils @{term "rGamma z = z * rGamma (z + 1)"} everywhere,
   207   whereas @{term "Gamma"} does not fulfil the analogous equation on the non-positive integers)
   208 
   209   We define the Gamma function (resp. its inverse) in the Euler form. This form has the advantage
   210   that it is a relatively simple limit that converges everywhere. The limit at the poles is 0
   211   (due to division by 0). The functional equation @{term "Gamma (z + 1) = z * Gamma z"} follows
   212   immediately from the definition.
   213 \<close>
   214 
   215 definition Gamma_series :: "('a :: {banach,real_normed_field}) \<Rightarrow> nat \<Rightarrow> 'a" where
   216   "Gamma_series z n = fact n * exp (z * of_real (ln (of_nat n))) / pochhammer z (n+1)"
   217 
   218 definition Gamma_series' :: "('a :: {banach,real_normed_field}) \<Rightarrow> nat \<Rightarrow> 'a" where
   219   "Gamma_series' z n = fact (n - 1) * exp (z * of_real (ln (of_nat n))) / pochhammer z n"
   220 
   221 definition rGamma_series :: "('a :: {banach,real_normed_field}) \<Rightarrow> nat \<Rightarrow> 'a" where
   222   "rGamma_series z n = pochhammer z (n+1) / (fact n * exp (z * of_real (ln (of_nat n))))"
   223 
   224 lemma Gamma_series_altdef: "Gamma_series z n = inverse (rGamma_series z n)"
   225   and rGamma_series_altdef: "rGamma_series z n = inverse (Gamma_series z n)"
   226   unfolding Gamma_series_def rGamma_series_def by simp_all
   227 
   228 lemma rGamma_series_minus_of_nat:
   229   "eventually (\<lambda>n. rGamma_series (- of_nat k) n = 0) sequentially"
   230   using eventually_ge_at_top[of k]
   231   by eventually_elim (auto simp: rGamma_series_def pochhammer_of_nat_eq_0_iff)
   232 
   233 lemma Gamma_series_minus_of_nat:
   234   "eventually (\<lambda>n. Gamma_series (- of_nat k) n = 0) sequentially"
   235   using eventually_ge_at_top[of k]
   236   by eventually_elim (auto simp: Gamma_series_def pochhammer_of_nat_eq_0_iff)
   237 
   238 lemma Gamma_series'_minus_of_nat:
   239   "eventually (\<lambda>n. Gamma_series' (- of_nat k) n = 0) sequentially"
   240   using eventually_gt_at_top[of k]
   241   by eventually_elim (auto simp: Gamma_series'_def pochhammer_of_nat_eq_0_iff)
   242 
   243 lemma rGamma_series_nonpos_Ints_LIMSEQ: "z \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> rGamma_series z \<longlonglongrightarrow> 0"
   244   by (elim nonpos_Ints_cases', hypsubst, subst tendsto_cong, rule rGamma_series_minus_of_nat, simp)
   245 
   246 lemma Gamma_series_nonpos_Ints_LIMSEQ: "z \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Gamma_series z \<longlonglongrightarrow> 0"
   247   by (elim nonpos_Ints_cases', hypsubst, subst tendsto_cong, rule Gamma_series_minus_of_nat, simp)
   248 
   249 lemma Gamma_series'_nonpos_Ints_LIMSEQ: "z \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Gamma_series' z \<longlonglongrightarrow> 0"
   250   by (elim nonpos_Ints_cases', hypsubst, subst tendsto_cong, rule Gamma_series'_minus_of_nat, simp)
   251 
   252 lemma Gamma_series_Gamma_series':
   253   assumes z: "z \<notin> \<int>\<^sub>\<le>\<^sub>0"
   254   shows   "(\<lambda>n. Gamma_series' z n / Gamma_series z n) \<longlonglongrightarrow> 1"
   255 proof (rule Lim_transform_eventually)
   256   from eventually_gt_at_top[of "0::nat"]
   257     show "eventually (\<lambda>n. z / of_nat n + 1 = Gamma_series' z n / Gamma_series z n) sequentially"
   258   proof eventually_elim
   259     fix n :: nat assume n: "n > 0"
   260     from n z have "Gamma_series' z n / Gamma_series z n = (z + of_nat n) / of_nat n"
   261       by (cases n, simp)
   262          (auto simp add: Gamma_series_def Gamma_series'_def pochhammer_rec'
   263                dest: pochhammer_eq_0_imp_nonpos_Int plus_of_nat_eq_0_imp)
   264     also from n have "\<dots> = z / of_nat n + 1" by (simp add: divide_simps)
   265     finally show "z / of_nat n + 1 = Gamma_series' z n / Gamma_series z n" ..
   266   qed
   267   have "(\<lambda>x. z / of_nat x) \<longlonglongrightarrow> 0"
   268     by (rule tendsto_norm_zero_cancel)
   269        (insert tendsto_mult[OF tendsto_const[of "norm z"] lim_inverse_n],
   270         simp add:  norm_divide inverse_eq_divide)
   271   from tendsto_add[OF this tendsto_const[of 1]] show "(\<lambda>n. z / of_nat n + 1) \<longlonglongrightarrow> 1" by simp
   272 qed
   273 
   274 
   275 subsection \<open>Convergence of the Euler series form\<close>
   276 
   277 text \<open>
   278   We now show that the series that defines the Gamma function in the Euler form converges
   279   and that the function defined by it is continuous on the complex halfspace with positive
   280   real part.
   281 
   282   We do this by showing that the logarithm of the Euler series is continuous and converges
   283   locally uniformly, which means that the log-Gamma function defined by its limit is also
   284   continuous.
   285 
   286   This will later allow us to lift holomorphicity and continuity from the log-Gamma
   287   function to the inverse of the Gamma function, and from that to the Gamma function itself.
   288 \<close>
   289 
   290 definition ln_Gamma_series :: "('a :: {banach,real_normed_field,ln}) \<Rightarrow> nat \<Rightarrow> 'a" where
   291   "ln_Gamma_series z n = z * ln (of_nat n) - ln z - (\<Sum>k=1..n. ln (z / of_nat k + 1))"
   292 
   293 definition ln_Gamma_series' :: "('a :: {banach,real_normed_field,ln}) \<Rightarrow> nat \<Rightarrow> 'a" where
   294   "ln_Gamma_series' z n =
   295      - euler_mascheroni*z - ln z + (\<Sum>k=1..n. z / of_nat n - ln (z / of_nat k + 1))"
   296 
   297 definition ln_Gamma :: "('a :: {banach,real_normed_field,ln}) \<Rightarrow> 'a" where
   298   "ln_Gamma z = lim (ln_Gamma_series z)"
   299 
   300 
   301 text \<open>
   302   We now show that the log-Gamma series converges locally uniformly for all complex numbers except
   303   the non-positive integers. We do this by proving that the series is locally Cauchy, adapting this
   304   proof:
   305   http://math.stackexchange.com/questions/887158/convergence-of-gammaz-lim-n-to-infty-fracnzn-prod-m-0nzm
   306 \<close>
   307 
   308 context
   309 begin
   310 
   311 private lemma ln_Gamma_series_complex_converges_aux:
   312   fixes z :: complex and k :: nat
   313   assumes z: "z \<noteq> 0" and k: "of_nat k \<ge> 2*norm z" "k \<ge> 2"
   314   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"
   315 proof -
   316   let ?k = "of_nat k :: complex" and ?z = "norm z"
   317   have "z *ln (1 - 1/?k) + ln (z/?k+1) = z*(ln (1 - 1/?k :: complex) + 1/?k) + (ln (1+z/?k) - z/?k)"
   318     by (simp add: algebra_simps)
   319   also have "norm ... \<le> ?z * norm (ln (1-1/?k) + 1/?k :: complex) + norm (ln (1+z/?k) - z/?k)"
   320     by (subst norm_mult [symmetric], rule norm_triangle_ineq)
   321   also have "norm (Ln (1 + -1/?k) - (-1/?k)) \<le> (norm (-1/?k))\<^sup>2 / (1 - norm(-1/?k))"
   322     using k by (intro Ln_approx_linear) (simp add: norm_divide)
   323   hence "?z * norm (ln (1-1/?k) + 1/?k) \<le> ?z * ((norm (1/?k))^2 / (1 - norm (1/?k)))"
   324     by (intro mult_left_mono) simp_all
   325   also have "... \<le> (?z * (of_nat k / (of_nat k - 1))) / of_nat k^2" using k
   326     by (simp add: field_simps power2_eq_square norm_divide)
   327   also have "... \<le> (?z * 2) / of_nat k^2" using k
   328     by (intro divide_right_mono mult_left_mono) (simp_all add: field_simps)
   329   also have "norm (ln (1+z/?k) - z/?k) \<le> norm (z/?k)^2 / (1 - norm (z/?k))" using k
   330     by (intro Ln_approx_linear) (simp add: norm_divide)
   331   hence "norm (ln (1+z/?k) - z/?k) \<le> ?z^2 / of_nat k^2 / (1 - ?z / of_nat k)"
   332     by (simp add: field_simps norm_divide)
   333   also have "... \<le> (?z^2 * (of_nat k / (of_nat k - ?z))) / of_nat k^2" using k
   334     by (simp add: field_simps power2_eq_square)
   335   also have "... \<le> (?z^2 * 2) / of_nat k^2" using k
   336     by (intro divide_right_mono mult_left_mono) (simp_all add: field_simps)
   337   also note add_divide_distrib [symmetric]
   338   finally show ?thesis by (simp only: distrib_left mult.commute)
   339 qed
   340 
   341 lemma ln_Gamma_series_complex_converges:
   342   assumes z: "z \<notin> \<int>\<^sub>\<le>\<^sub>0"
   343   assumes d: "d > 0" "\<And>n. n \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> norm (z - of_int n) > d"
   344   shows "uniformly_convergent_on (ball z d) (\<lambda>n z. ln_Gamma_series z n :: complex)"
   345 proof (intro Cauchy_uniformly_convergent uniformly_Cauchy_onI')
   346   fix e :: real assume e: "e > 0"
   347   def e'' \<equiv> "SUP t:ball z d. norm t + norm t^2"
   348   def e' \<equiv> "e / (2*e'')"
   349   have "bounded ((\<lambda>t. norm t + norm t^2) ` cball z d)"
   350     by (intro compact_imp_bounded compact_continuous_image) (auto intro!: continuous_intros)
   351   hence "bounded ((\<lambda>t. norm t + norm t^2) ` ball z d)" by (rule bounded_subset) auto
   352   hence bdd: "bdd_above ((\<lambda>t. norm t + norm t^2) ` ball z d)" by (rule bounded_imp_bdd_above)
   353 
   354   with z d(1) d(2)[of "-1"] have e''_pos: "e'' > 0" unfolding e''_def
   355     by (subst less_cSUP_iff) (auto intro!: add_pos_nonneg bexI[of _ z])
   356   have e'': "norm t + norm t^2 \<le> e''" if "t \<in> ball z d" for t unfolding e''_def using that
   357     by (rule cSUP_upper[OF _ bdd])
   358   from e z e''_pos have e': "e' > 0" unfolding e'_def
   359     by (intro divide_pos_pos mult_pos_pos add_pos_pos) (simp_all add: field_simps)
   360 
   361   have "summable (\<lambda>k. inverse ((real_of_nat k)^2))"
   362     by (rule inverse_power_summable) simp
   363   from summable_partial_sum_bound[OF this e'] guess M . note M = this
   364 
   365   def N \<equiv> "max 2 (max (nat \<lceil>2 * (norm z + d)\<rceil>) M)"
   366   {
   367     from d have "\<lceil>2 * (cmod z + d)\<rceil> \<ge> \<lceil>0::real\<rceil>"
   368       by (intro ceiling_mono mult_nonneg_nonneg add_nonneg_nonneg) simp_all
   369     hence "2 * (norm z + d) \<le> of_nat (nat \<lceil>2 * (norm z + d)\<rceil>)" unfolding N_def
   370       by (simp_all add: le_of_int_ceiling)
   371     also have "... \<le> of_nat N" unfolding N_def
   372       by (subst of_nat_le_iff) (rule max.coboundedI2, rule max.cobounded1)
   373     finally have "of_nat N \<ge> 2 * (norm z + d)" .
   374     moreover have "N \<ge> 2" "N \<ge> M" unfolding N_def by simp_all
   375     moreover have "(\<Sum>k=m..n. 1/(of_nat k)\<^sup>2) < e'" if "m \<ge> N" for m n
   376       using M[OF order.trans[OF \<open>N \<ge> M\<close> that]] unfolding real_norm_def
   377       by (subst (asm) abs_of_nonneg) (auto intro: setsum_nonneg simp: divide_simps)
   378     moreover note calculation
   379   } note N = this
   380 
   381   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"
   382     unfolding dist_complex_def
   383   proof (intro exI[of _ N] ballI allI impI)
   384     fix t m n assume t: "t \<in> ball z d" and mn: "m \<ge> N" "n > m"
   385     from d(2)[of 0] t have "0 < dist z 0 - dist z t" by (simp add: field_simps dist_complex_def)
   386     also have "dist z 0 - dist z t \<le> dist 0 t" using dist_triangle[of 0 z t]
   387       by (simp add: dist_commute)
   388     finally have t_nz: "t \<noteq> 0" by auto
   389 
   390     have "norm t \<le> norm z + norm (t - z)" by (rule norm_triangle_sub)
   391     also from t have "norm (t - z) < d" by (simp add: dist_complex_def norm_minus_commute)
   392     also have "2 * (norm z + d) \<le> of_nat N" by (rule N)
   393     also have "N \<le> m" by (rule mn)
   394     finally have norm_t: "2 * norm t < of_nat m" by simp
   395 
   396     have "ln_Gamma_series t m - ln_Gamma_series t n =
   397               (-(t * Ln (of_nat n)) - (-(t * Ln (of_nat m)))) +
   398               ((\<Sum>k=1..n. Ln (t / of_nat k + 1)) - (\<Sum>k=1..m. Ln (t / of_nat k + 1)))"
   399       by (simp add: ln_Gamma_series_def algebra_simps)
   400     also have "(\<Sum>k=1..n. Ln (t / of_nat k + 1)) - (\<Sum>k=1..m. Ln (t / of_nat k + 1)) =
   401                  (\<Sum>k\<in>{1..n}-{1..m}. Ln (t / of_nat k + 1))" using mn
   402       by (simp add: setsum_diff)
   403     also from mn have "{1..n}-{1..m} = {Suc m..n}" by fastforce
   404     also have "-(t * Ln (of_nat n)) - (-(t * Ln (of_nat m))) =
   405                    (\<Sum>k = Suc m..n. t * Ln (of_nat (k - 1)) - t * Ln (of_nat k))" using mn
   406       by (subst setsum_telescope'' [symmetric]) simp_all
   407     also have "... = (\<Sum>k = Suc m..n. t * Ln (of_nat (k - 1) / of_nat k))" using mn N
   408       by (intro setsum_cong_Suc)
   409          (simp_all del: of_nat_Suc add: field_simps Ln_of_nat Ln_of_nat_over_of_nat)
   410     also have "of_nat (k - 1) / of_nat k = 1 - 1 / (of_nat k :: complex)" if "k \<in> {Suc m..n}" for k
   411       using that of_nat_eq_0_iff[of "Suc i" for i] by (cases k) (simp_all add: divide_simps)
   412     hence "(\<Sum>k = Suc m..n. t * Ln (of_nat (k - 1) / of_nat k)) =
   413              (\<Sum>k = Suc m..n. t * Ln (1 - 1 / of_nat k))" using mn N
   414       by (intro setsum.cong) simp_all
   415     also note setsum.distrib [symmetric]
   416     also have "norm (\<Sum>k=Suc m..n. t * Ln (1 - 1/of_nat k) + Ln (t/of_nat k + 1)) \<le>
   417       (\<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
   418       by (intro order.trans[OF norm_setsum setsum_mono[OF ln_Gamma_series_complex_converges_aux]]) simp_all
   419     also have "... \<le> 2 * (norm t + norm t^2) * (\<Sum>k=Suc m..n. 1 / (of_nat k)\<^sup>2)"
   420       by (simp add: setsum_right_distrib)
   421     also have "... < 2 * (norm t + norm t^2) * e'" using mn z t_nz
   422       by (intro mult_strict_left_mono N mult_pos_pos add_pos_pos) simp_all
   423     also from e''_pos have "... = e * ((cmod t + (cmod t)\<^sup>2) / e'')"
   424       by (simp add: e'_def field_simps power2_eq_square)
   425     also from e''[OF t] e''_pos e
   426       have "\<dots> \<le> e * 1" by (intro mult_left_mono) (simp_all add: field_simps)
   427     finally show "norm (ln_Gamma_series t m - ln_Gamma_series t n) < e" by simp
   428   qed
   429 qed
   430 
   431 end
   432 
   433 lemma ln_Gamma_series_complex_converges':
   434   assumes z: "(z :: complex) \<notin> \<int>\<^sub>\<le>\<^sub>0"
   435   shows "\<exists>d>0. uniformly_convergent_on (ball z d) (\<lambda>n z. ln_Gamma_series z n)"
   436 proof -
   437   def d' \<equiv> "Re z"
   438   def d \<equiv> "if d' > 0 then d' / 2 else norm (z - of_int (round d')) / 2"
   439   have "of_int (round d') \<in> \<int>\<^sub>\<le>\<^sub>0" if "d' \<le> 0" using that
   440     by (intro nonpos_Ints_of_int) (simp_all add: round_def)
   441   with assms have d_pos: "d > 0" unfolding d_def by (force simp: not_less)
   442 
   443   have "d < cmod (z - of_int n)" if "n \<in> \<int>\<^sub>\<le>\<^sub>0" for n
   444   proof (cases "Re z > 0")
   445     case True
   446     from nonpos_Ints_nonpos[OF that] have n: "n \<le> 0" by simp
   447     from True have "d = Re z/2" by (simp add: d_def d'_def)
   448     also from n True have "\<dots> < Re (z - of_int n)" by simp
   449     also have "\<dots> \<le> norm (z - of_int n)" by (rule complex_Re_le_cmod)
   450     finally show ?thesis .
   451   next
   452     case False
   453     with assms nonpos_Ints_of_int[of "round (Re z)"]
   454       have "z \<noteq> of_int (round d')" by (auto simp: not_less)
   455     with False have "d < norm (z - of_int (round d'))" by (simp add: d_def d'_def)
   456     also have "\<dots> \<le> norm (z - of_int n)" unfolding d'_def by (rule round_Re_minimises_norm)
   457     finally show ?thesis .
   458   qed
   459   hence conv: "uniformly_convergent_on (ball z d) (\<lambda>n z. ln_Gamma_series z n)"
   460     by (intro ln_Gamma_series_complex_converges d_pos z) simp_all
   461   from d_pos conv show ?thesis by blast
   462 qed
   463 
   464 lemma ln_Gamma_series_complex_converges'': "(z :: complex) \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> convergent (ln_Gamma_series z)"
   465   by (drule ln_Gamma_series_complex_converges') (auto intro: uniformly_convergent_imp_convergent)
   466 
   467 lemma ln_Gamma_complex_LIMSEQ: "(z :: complex) \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> ln_Gamma_series z \<longlonglongrightarrow> ln_Gamma z"
   468   using ln_Gamma_series_complex_converges'' by (simp add: convergent_LIMSEQ_iff ln_Gamma_def)
   469 
   470 lemma exp_ln_Gamma_series_complex:
   471   assumes "n > 0" "z \<notin> \<int>\<^sub>\<le>\<^sub>0"
   472   shows   "exp (ln_Gamma_series z n :: complex) = Gamma_series z n"
   473 proof -
   474   from assms have "z \<noteq> 0" by (intro notI) auto
   475   with assms have "exp (ln_Gamma_series z n) =
   476           (of_nat n) powr z / (z * (\<Prod>k=1..n. exp (Ln (z / of_nat k + 1))))"
   477     unfolding ln_Gamma_series_def powr_def by (simp add: exp_diff exp_setsum)
   478   also from assms have "(\<Prod>k=1..n. exp (Ln (z / of_nat k + 1))) = (\<Prod>k=1..n. z / of_nat k + 1)"
   479     by (intro setprod.cong[OF refl], subst exp_Ln) (auto simp: field_simps plus_of_nat_eq_0_imp)
   480   also have "... = (\<Prod>k=1..n. z + k) / fact n" unfolding fact_altdef
   481     by (subst setprod_dividef [symmetric]) (simp_all add: field_simps)
   482   also from assms have "z * ... = (\<Prod>k=0..n. z + k) / fact n"
   483     by (cases n) (simp_all add: setprod_nat_ivl_1_Suc)
   484   also have "(\<Prod>k=0..n. z + k) = pochhammer z (Suc n)" unfolding pochhammer_def by simp
   485   also have "of_nat n powr z / (pochhammer z (Suc n) / fact n) = Gamma_series z n"
   486     unfolding Gamma_series_def using assms by (simp add: divide_simps powr_def Ln_of_nat)
   487   finally show ?thesis .
   488 qed
   489 
   490 
   491 lemma ln_Gamma_series'_aux:
   492   assumes "(z::complex) \<notin> \<int>\<^sub>\<le>\<^sub>0"
   493   shows   "(\<lambda>k. z / of_nat (Suc k) - ln (1 + z / of_nat (Suc k))) sums
   494               (ln_Gamma z + euler_mascheroni * z + ln z)" (is "?f sums ?s")
   495 unfolding sums_def
   496 proof (rule Lim_transform)
   497   show "(\<lambda>n. ln_Gamma_series z n + of_real (harm n - ln (of_nat n)) * z + ln z) \<longlonglongrightarrow> ?s"
   498     (is "?g \<longlonglongrightarrow> _")
   499     by (intro tendsto_intros ln_Gamma_complex_LIMSEQ euler_mascheroni_LIMSEQ_of_real assms)
   500 
   501   have A: "eventually (\<lambda>n. (\<Sum>k<n. ?f k) - ?g n = 0) sequentially"
   502     using eventually_gt_at_top[of "0::nat"]
   503   proof eventually_elim
   504     fix n :: nat assume n: "n > 0"
   505     have "(\<Sum>k<n. ?f k) = (\<Sum>k=1..n. z / of_nat k - ln (1 + z / of_nat k))"
   506       by (subst atLeast0LessThan [symmetric], subst setsum_shift_bounds_Suc_ivl [symmetric],
   507           subst atLeastLessThanSuc_atLeastAtMost) simp_all
   508     also have "\<dots> = z * of_real (harm n) - (\<Sum>k=1..n. ln (1 + z / of_nat k))"
   509       by (simp add: harm_def setsum_subtractf setsum_right_distrib divide_inverse)
   510     also from n have "\<dots> - ?g n = 0"
   511       by (simp add: ln_Gamma_series_def setsum_subtractf algebra_simps Ln_of_nat)
   512     finally show "(\<Sum>k<n. ?f k) - ?g n = 0" .
   513   qed
   514   show "(\<lambda>n. (\<Sum>k<n. ?f k) - ?g n) \<longlonglongrightarrow> 0" by (subst tendsto_cong[OF A]) simp_all
   515 qed
   516 
   517 
   518 lemma uniformly_summable_deriv_ln_Gamma:
   519   assumes z: "(z :: 'a :: {real_normed_field,banach}) \<noteq> 0" and d: "d > 0" "d \<le> norm z/2"
   520   shows "uniformly_convergent_on (ball z d)
   521             (\<lambda>k z. \<Sum>i<k. inverse (of_nat (Suc i)) - inverse (z + of_nat (Suc i)))"
   522            (is "uniformly_convergent_on _ (\<lambda>k z. \<Sum>i<k. ?f i z)")
   523 proof (rule weierstrass_m_test'_ev)
   524   {
   525     fix t assume t: "t \<in> ball z d"
   526     have "norm z = norm (t + (z - t))" by simp
   527     have "norm (t + (z - t)) \<le> norm t + norm (z - t)" by (rule norm_triangle_ineq)
   528     also from t d have "norm (z - t) < norm z / 2" by (simp add: dist_norm)
   529     finally have A: "norm t > norm z / 2" using z by (simp add: field_simps)
   530 
   531     have "norm t = norm (z + (t - z))" by simp
   532     also have "\<dots> \<le> norm z + norm (t - z)" by (rule norm_triangle_ineq)
   533     also from t d have "norm (t - z) \<le> norm z / 2" by (simp add: dist_norm norm_minus_commute)
   534     also from z have "\<dots> < norm z" by simp
   535     finally have B: "norm t < 2 * norm z" by simp
   536     note A B
   537   } note ball = this
   538 
   539   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"
   540     using eventually_gt_at_top apply eventually_elim
   541   proof safe
   542     fix t :: 'a assume t: "t \<in> ball z d"
   543     from z ball[OF t] have t_nz: "t \<noteq> 0" by auto
   544     fix n :: nat assume n: "n > nat \<lceil>4 * norm z\<rceil>"
   545     from ball[OF t] t_nz have "4 * norm z > 2 * norm t" by simp
   546     also from n have "\<dots>  < of_nat n" by linarith
   547     finally have n: "of_nat n > 2 * norm t" .
   548     hence "of_nat n > norm t" by simp
   549     hence t': "t \<noteq> -of_nat (Suc n)" by (intro notI) (simp del: of_nat_Suc)
   550 
   551     with t_nz have "?f n t = 1 / (of_nat (Suc n) * (1 + of_nat (Suc n)/t))"
   552       by (simp add: divide_simps eq_neg_iff_add_eq_0 del: of_nat_Suc)
   553     also have "norm \<dots> = inverse (of_nat (Suc n)) * inverse (norm (of_nat (Suc n)/t + 1))"
   554       by (simp add: norm_divide norm_mult divide_simps add_ac del: of_nat_Suc)
   555     also {
   556       from z t_nz ball[OF t] have "of_nat (Suc n) / (4 * norm z) \<le> of_nat (Suc n) / (2 * norm t)"
   557         by (intro divide_left_mono mult_pos_pos) simp_all
   558       also have "\<dots> < norm (of_nat (Suc n) / t) - norm (1 :: 'a)"
   559         using t_nz n by (simp add: field_simps norm_divide del: of_nat_Suc)
   560       also have "\<dots> \<le> norm (of_nat (Suc n)/t + 1)" by (rule norm_diff_ineq)
   561       finally have "inverse (norm (of_nat (Suc n)/t + 1)) \<le> 4 * norm z / of_nat (Suc n)"
   562         using z by (simp add: divide_simps norm_divide mult_ac del: of_nat_Suc)
   563     }
   564     also have "inverse (real_of_nat (Suc n)) * (4 * norm z / real_of_nat (Suc n)) =
   565                  4 * norm z * inverse (of_nat (Suc n)^2)"
   566                  by (simp add: divide_simps power2_eq_square del: of_nat_Suc)
   567     finally show "norm (?f n t) \<le> 4 * norm z * inverse (of_nat (Suc n)^2)"
   568       by (simp del: of_nat_Suc)
   569   qed
   570 next
   571   show "summable (\<lambda>n. 4 * norm z * inverse ((of_nat (Suc n))^2))"
   572     by (subst summable_Suc_iff) (simp add: summable_mult inverse_power_summable)
   573 qed
   574 
   575 lemma summable_deriv_ln_Gamma:
   576   "z \<noteq> (0 :: 'a :: {real_normed_field,banach}) \<Longrightarrow>
   577      summable (\<lambda>n. inverse (of_nat (Suc n)) - inverse (z + of_nat (Suc n)))"
   578   unfolding summable_iff_convergent
   579   by (rule uniformly_convergent_imp_convergent,
   580       rule uniformly_summable_deriv_ln_Gamma[of z "norm z/2"]) simp_all
   581 
   582 
   583 definition Polygamma :: "nat \<Rightarrow> ('a :: {real_normed_field,banach}) \<Rightarrow> 'a" where
   584   "Polygamma n z = (if n = 0 then
   585       (\<Sum>k. inverse (of_nat (Suc k)) - inverse (z + of_nat k)) - euler_mascheroni else
   586       (-1)^Suc n * fact n * (\<Sum>k. inverse ((z + of_nat k)^Suc n)))"
   587 
   588 abbreviation Digamma :: "('a :: {real_normed_field,banach}) \<Rightarrow> 'a" where
   589   "Digamma \<equiv> Polygamma 0"
   590 
   591 lemma Digamma_def:
   592   "Digamma z = (\<Sum>k. inverse (of_nat (Suc k)) - inverse (z + of_nat k)) - euler_mascheroni"
   593   by (simp add: Polygamma_def)
   594 
   595 
   596 lemma summable_Digamma:
   597   assumes "(z :: 'a :: {real_normed_field,banach}) \<noteq> 0"
   598   shows   "summable (\<lambda>n. inverse (of_nat (Suc n)) - inverse (z + of_nat n))"
   599 proof -
   600   have sums: "(\<lambda>n. inverse (z + of_nat (Suc n)) - inverse (z + of_nat n)) sums
   601                        (0 - inverse (z + of_nat 0))"
   602     by (intro telescope_sums filterlim_compose[OF tendsto_inverse_0]
   603               tendsto_add_filterlim_at_infinity[OF tendsto_const] tendsto_of_nat)
   604   from summable_add[OF summable_deriv_ln_Gamma[OF assms] sums_summable[OF sums]]
   605     show "summable (\<lambda>n. inverse (of_nat (Suc n)) - inverse (z + of_nat n))" by simp
   606 qed
   607 
   608 lemma summable_offset:
   609   assumes "summable (\<lambda>n. f (n + k) :: 'a :: real_normed_vector)"
   610   shows   "summable f"
   611 proof -
   612   from assms have "convergent (\<lambda>m. \<Sum>n<m. f (n + k))" by (simp add: summable_iff_convergent)
   613   hence "convergent (\<lambda>m. (\<Sum>n<k. f n) + (\<Sum>n<m. f (n + k)))"
   614     by (intro convergent_add convergent_const)
   615   also have "(\<lambda>m. (\<Sum>n<k. f n) + (\<Sum>n<m. f (n + k))) = (\<lambda>m. \<Sum>n<m+k. f n)"
   616   proof
   617     fix m :: nat
   618     have "{..<m+k} = {..<k} \<union> {k..<m+k}" by auto
   619     also have "(\<Sum>n\<in>\<dots>. f n) = (\<Sum>n<k. f n) + (\<Sum>n=k..<m+k. f n)"
   620       by (rule setsum.union_disjoint) auto
   621     also have "(\<Sum>n=k..<m+k. f n) = (\<Sum>n=0..<m+k-k. f (n + k))"
   622       by (intro setsum.reindex_cong[of "\<lambda>n. n + k"])
   623          (simp, subst image_add_atLeastLessThan, auto)
   624     finally show "(\<Sum>n<k. f n) + (\<Sum>n<m. f (n + k)) = (\<Sum>n<m+k. f n)" by (simp add: atLeast0LessThan)
   625   qed
   626   finally have "(\<lambda>a. setsum f {..<a}) \<longlonglongrightarrow> lim (\<lambda>m. setsum f {..<m + k})"
   627     by (auto simp: convergent_LIMSEQ_iff dest: LIMSEQ_offset)
   628   thus ?thesis by (auto simp: summable_iff_convergent convergent_def)
   629 qed
   630 
   631 lemma Polygamma_converges:
   632   fixes z :: "'a :: {real_normed_field,banach}"
   633   assumes z: "z \<noteq> 0" and n: "n \<ge> 2"
   634   shows "uniformly_convergent_on (ball z d) (\<lambda>k z. \<Sum>i<k. inverse ((z + of_nat i)^n))"
   635 proof (rule weierstrass_m_test'_ev)
   636   def e \<equiv> "(1 + d / norm z)"
   637   def m \<equiv> "nat \<lceil>norm z * e\<rceil>"
   638   {
   639     fix t assume t: "t \<in> ball z d"
   640     have "norm t = norm (z + (t - z))" by simp
   641     also have "\<dots> \<le> norm z + norm (t - z)" by (rule norm_triangle_ineq)
   642     also from t have "norm (t - z) < d" by (simp add: dist_norm norm_minus_commute)
   643     finally have "norm t < norm z * e" using z by (simp add: divide_simps e_def)
   644   } note ball = this
   645 
   646   show "eventually (\<lambda>k. \<forall>t\<in>ball z d. norm (inverse ((t + of_nat k)^n)) \<le>
   647             inverse (of_nat (k - m)^n)) sequentially"
   648     using eventually_gt_at_top[of m] apply eventually_elim
   649   proof (intro ballI)
   650     fix k :: nat and t :: 'a assume k: "k > m" and t: "t \<in> ball z d"
   651     from k have "real_of_nat (k - m) = of_nat k - of_nat m" by (simp add: of_nat_diff)
   652     also have "\<dots> \<le> norm (of_nat k :: 'a) - norm z * e"
   653       unfolding m_def by (subst norm_of_nat) linarith
   654     also from ball[OF t] have "\<dots> \<le> norm (of_nat k :: 'a) - norm t" by simp
   655     also have "\<dots> \<le> norm (of_nat k + t)" by (rule norm_diff_ineq)
   656     finally have "inverse ((norm (t + of_nat k))^n) \<le> inverse (real_of_nat (k - m)^n)" using k n
   657       by (intro le_imp_inverse_le power_mono) (simp_all add: add_ac del: of_nat_Suc)
   658     thus "norm (inverse ((t + of_nat k)^n)) \<le> inverse (of_nat (k - m)^n)"
   659       by (simp add: norm_inverse norm_power power_inverse)
   660   qed
   661 
   662   have "summable (\<lambda>k. inverse ((real_of_nat k)^n))"
   663     using inverse_power_summable[of n] n by simp
   664   hence "summable (\<lambda>k. inverse ((real_of_nat (k + m - m))^n))" by simp
   665   thus "summable (\<lambda>k. inverse ((real_of_nat (k - m))^n))" by (rule summable_offset)
   666 qed
   667 
   668 lemma Polygamma_converges':
   669   fixes z :: "'a :: {real_normed_field,banach}"
   670   assumes z: "z \<noteq> 0" and n: "n \<ge> 2"
   671   shows "summable (\<lambda>k. inverse ((z + of_nat k)^n))"
   672   using uniformly_convergent_imp_convergent[OF Polygamma_converges[OF assms, of 1], of z]
   673   by (simp add: summable_iff_convergent)
   674 
   675 lemma has_field_derivative_ln_Gamma_complex [derivative_intros]:
   676   fixes z :: complex
   677   assumes z: "z \<notin> \<real>\<^sub>\<le>\<^sub>0"
   678   shows   "(ln_Gamma has_field_derivative Digamma z) (at z)"
   679 proof -
   680   have not_nonpos_Int [simp]: "t \<notin> \<int>\<^sub>\<le>\<^sub>0" if "Re t > 0" for t
   681     using that by (auto elim!: nonpos_Ints_cases')
   682   from z have z': "z \<notin> \<int>\<^sub>\<le>\<^sub>0" and z'': "z \<noteq> 0" using nonpos_Ints_subset_nonpos_Reals nonpos_Reals_zero_I
   683      by blast+
   684   let ?f' = "\<lambda>z k. inverse (of_nat (Suc k)) - inverse (z + of_nat (Suc k))"
   685   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"
   686   def d \<equiv> "min (norm z/2) (if Im z = 0 then Re z / 2 else abs (Im z) / 2)"
   687   from z have d: "d > 0" "norm z/2 \<ge> d" by (auto simp add: complex_nonpos_Reals_iff d_def)
   688   have ball: "Im t = 0 \<longrightarrow> Re t > 0" if "dist z t < d" for t
   689     using no_nonpos_Real_in_ball[OF z, of t] that unfolding d_def by (force simp add: complex_nonpos_Reals_iff)
   690   have sums: "(\<lambda>n. inverse (z + of_nat (Suc n)) - inverse (z + of_nat n)) sums
   691                        (0 - inverse (z + of_nat 0))"
   692     by (intro telescope_sums filterlim_compose[OF tendsto_inverse_0]
   693               tendsto_add_filterlim_at_infinity[OF tendsto_const] tendsto_of_nat)
   694 
   695   have "((\<lambda>z. \<Sum>n. ?f z n) has_field_derivative ?F' z) (at z)"
   696     using d z ln_Gamma_series'_aux[OF z']
   697     apply (intro has_field_derivative_series'(2)[of "ball z d" _ _ z] uniformly_summable_deriv_ln_Gamma)
   698     apply (auto intro!: derivative_eq_intros add_pos_pos mult_pos_pos dest!: ball
   699              simp: field_simps sums_iff nonpos_Reals_divide_of_nat_iff
   700              simp del: of_nat_Suc)
   701     apply (auto simp add: complex_nonpos_Reals_iff)
   702     done
   703   with z have "((\<lambda>z. (\<Sum>k. ?f z k) - euler_mascheroni * z - Ln z) has_field_derivative
   704                    ?F' z - euler_mascheroni - inverse z) (at z)"
   705     by (force intro!: derivative_eq_intros simp: Digamma_def)
   706   also have "?F' z - euler_mascheroni - inverse z = (?F' z + -inverse z) - euler_mascheroni" by simp
   707   also from sums have "-inverse z = (\<Sum>n. inverse (z + of_nat (Suc n)) - inverse (z + of_nat n))"
   708     by (simp add: sums_iff)
   709   also from sums summable_deriv_ln_Gamma[OF z'']
   710     have "?F' z + \<dots> =  (\<Sum>n. inverse (of_nat (Suc n)) - inverse (z + of_nat n))"
   711     by (subst suminf_add) (simp_all add: add_ac sums_iff)
   712   also have "\<dots> - euler_mascheroni = Digamma z" by (simp add: Digamma_def)
   713   finally have "((\<lambda>z. (\<Sum>k. ?f z k) - euler_mascheroni * z - Ln z)
   714                     has_field_derivative Digamma z) (at z)" .
   715   moreover from eventually_nhds_ball[OF d(1), of z]
   716     have "eventually (\<lambda>z. ln_Gamma z = (\<Sum>k. ?f z k) - euler_mascheroni * z - Ln z) (nhds z)"
   717   proof eventually_elim
   718     fix t assume "t \<in> ball z d"
   719     hence "t \<notin> \<int>\<^sub>\<le>\<^sub>0" by (auto dest!: ball elim!: nonpos_Ints_cases)
   720     from ln_Gamma_series'_aux[OF this]
   721       show "ln_Gamma t = (\<Sum>k. ?f t k) - euler_mascheroni * t - Ln t" by (simp add: sums_iff)
   722   qed
   723   ultimately show ?thesis by (subst DERIV_cong_ev[OF refl _ refl])
   724 qed
   725 
   726 declare has_field_derivative_ln_Gamma_complex[THEN DERIV_chain2, derivative_intros]
   727 
   728 
   729 lemma Digamma_1 [simp]: "Digamma (1 :: 'a :: {real_normed_field,banach}) = - euler_mascheroni"
   730   by (simp add: Digamma_def)
   731 
   732 lemma Digamma_plus1:
   733   assumes "z \<noteq> 0"
   734   shows   "Digamma (z+1) = Digamma z + 1/z"
   735 proof -
   736   have sums: "(\<lambda>k. inverse (z + of_nat k) - inverse (z + of_nat (Suc k)))
   737                   sums (inverse (z + of_nat 0) - 0)"
   738     by (intro telescope_sums'[OF filterlim_compose[OF tendsto_inverse_0]]
   739               tendsto_add_filterlim_at_infinity[OF tendsto_const] tendsto_of_nat)
   740   have "Digamma (z+1) = (\<Sum>k. inverse (of_nat (Suc k)) - inverse (z + of_nat (Suc k))) -
   741           euler_mascheroni" (is "_ = suminf ?f - _") by (simp add: Digamma_def add_ac)
   742   also have "suminf ?f = (\<Sum>k. inverse (of_nat (Suc k)) - inverse (z + of_nat k)) +
   743                          (\<Sum>k. inverse (z + of_nat k) - inverse (z + of_nat (Suc k)))"
   744     using summable_Digamma[OF assms] sums by (subst suminf_add) (simp_all add: add_ac sums_iff)
   745   also have "(\<Sum>k. inverse (z + of_nat k) - inverse (z + of_nat (Suc k))) = 1/z"
   746     using sums by (simp add: sums_iff inverse_eq_divide)
   747   finally show ?thesis by (simp add: Digamma_def[of z])
   748 qed
   749 
   750 lemma Polygamma_plus1:
   751   assumes "z \<noteq> 0"
   752   shows   "Polygamma n (z + 1) = Polygamma n z + (-1)^n * fact n / (z ^ Suc n)"
   753 proof (cases "n = 0")
   754   assume n: "n \<noteq> 0"
   755   let ?f = "\<lambda>k. inverse ((z + of_nat k) ^ Suc n)"
   756   have "Polygamma n (z + 1) = (-1) ^ Suc n * fact n * (\<Sum>k. ?f (k+1))"
   757     using n by (simp add: Polygamma_def add_ac)
   758   also have "(\<Sum>k. ?f (k+1)) + (\<Sum>k<1. ?f k) = (\<Sum>k. ?f k)"
   759     using Polygamma_converges'[OF assms, of "Suc n"] n
   760     by (subst suminf_split_initial_segment [symmetric]) simp_all
   761   hence "(\<Sum>k. ?f (k+1)) = (\<Sum>k. ?f k) - inverse (z ^ Suc n)" by (simp add: algebra_simps)
   762   also have "(-1) ^ Suc n * fact n * ((\<Sum>k. ?f k) - inverse (z ^ Suc n)) =
   763                Polygamma n z + (-1)^n * fact n / (z ^ Suc n)" using n
   764     by (simp add: inverse_eq_divide algebra_simps Polygamma_def)
   765   finally show ?thesis .
   766 qed (insert assms, simp add: Digamma_plus1 inverse_eq_divide)
   767 
   768 lemma Digamma_of_nat:
   769   "Digamma (of_nat (Suc n) :: 'a :: {real_normed_field,banach}) = harm n - euler_mascheroni"
   770 proof (induction n)
   771   case (Suc n)
   772   have "Digamma (of_nat (Suc (Suc n)) :: 'a) = Digamma (of_nat (Suc n) + 1)" by simp
   773   also have "\<dots> = Digamma (of_nat (Suc n)) + inverse (of_nat (Suc n))"
   774     by (subst Digamma_plus1) (simp_all add: inverse_eq_divide del: of_nat_Suc)
   775   also have "Digamma (of_nat (Suc n) :: 'a) = harm n - euler_mascheroni " by (rule Suc)
   776   also have "\<dots> + inverse (of_nat (Suc n)) = harm (Suc n) - euler_mascheroni"
   777     by (simp add: harm_Suc)
   778   finally show ?case .
   779 qed (simp add: harm_def)
   780 
   781 lemma Digamma_numeral: "Digamma (numeral n) = harm (pred_numeral n) - euler_mascheroni"
   782   by (subst of_nat_numeral[symmetric], subst numeral_eq_Suc, subst Digamma_of_nat) (rule refl)
   783 
   784 lemma Polygamma_of_real: "x \<noteq> 0 \<Longrightarrow> Polygamma n (of_real x) = of_real (Polygamma n x)"
   785   unfolding Polygamma_def using summable_Digamma[of x] Polygamma_converges'[of x "Suc n"]
   786   by (simp_all add: suminf_of_real)
   787 
   788 lemma Polygamma_Real: "z \<in> \<real> \<Longrightarrow> z \<noteq> 0 \<Longrightarrow> Polygamma n z \<in> \<real>"
   789   by (elim Reals_cases, hypsubst, subst Polygamma_of_real) simp_all
   790 
   791 lemma Digamma_half_integer:
   792   "Digamma (of_nat n + 1/2 :: 'a :: {real_normed_field,banach}) =
   793       (\<Sum>k<n. 2 / (of_nat (2*k+1))) - euler_mascheroni - of_real (2 * ln 2)"
   794 proof (induction n)
   795   case 0
   796   have "Digamma (1/2 :: 'a) = of_real (Digamma (1/2))" by (simp add: Polygamma_of_real [symmetric])
   797   also have "Digamma (1/2::real) =
   798                (\<Sum>k. inverse (of_nat (Suc k)) - inverse (of_nat k + 1/2)) - euler_mascheroni"
   799     by (simp add: Digamma_def add_ac)
   800   also have "(\<Sum>k. inverse (of_nat (Suc k) :: real) - inverse (of_nat k + 1/2)) =
   801              (\<Sum>k. inverse (1/2) * (inverse (2 * of_nat (Suc k)) - inverse (2 * of_nat k + 1)))"
   802     by (simp_all add: add_ac inverse_mult_distrib[symmetric] ring_distribs del: inverse_divide)
   803   also have "\<dots> = - 2 * ln 2" using sums_minus[OF alternating_harmonic_series_sums']
   804     by (subst suminf_mult) (simp_all add: algebra_simps sums_iff)
   805   finally show ?case by simp
   806 next
   807   case (Suc n)
   808   have nz: "2 * of_nat n + (1:: 'a) \<noteq> 0"
   809      using of_nat_neq_0[of "2*n"] by (simp only: of_nat_Suc) (simp add: add_ac)
   810   hence nz': "of_nat n + (1/2::'a) \<noteq> 0" by (simp add: field_simps)
   811   have "Digamma (of_nat (Suc n) + 1/2 :: 'a) = Digamma (of_nat n + 1/2 + 1)" by simp
   812   also from nz' have "\<dots> = Digamma (of_nat n + 1 / 2) + 1 / (of_nat n + 1 / 2)"
   813     by (rule Digamma_plus1)
   814   also from nz nz' have "1 / (of_nat n + 1 / 2 :: 'a) = 2 / (2 * of_nat n + 1)"
   815     by (subst divide_eq_eq) simp_all
   816   also note Suc
   817   finally show ?case by (simp add: add_ac)
   818 qed
   819 
   820 lemma Digamma_one_half: "Digamma (1/2) = - euler_mascheroni - of_real (2 * ln 2)"
   821   using Digamma_half_integer[of 0] by simp
   822 
   823 lemma Digamma_real_three_halves_pos: "Digamma (3/2 :: real) > 0"
   824 proof -
   825   have "-Digamma (3/2 :: real) = -Digamma (of_nat 1 + 1/2)" by simp
   826   also have "\<dots> = 2 * ln 2 + euler_mascheroni - 2" by (subst Digamma_half_integer) simp
   827   also note euler_mascheroni_less_13_over_22
   828   also note ln2_le_25_over_36
   829   finally show ?thesis by simp
   830 qed
   831 
   832 
   833 lemma has_field_derivative_Polygamma [derivative_intros]:
   834   fixes z :: "'a :: {real_normed_field,euclidean_space}"
   835   assumes z: "z \<notin> \<int>\<^sub>\<le>\<^sub>0"
   836   shows "(Polygamma n has_field_derivative Polygamma (Suc n) z) (at z within A)"
   837 proof (rule has_field_derivative_at_within, cases "n = 0")
   838   assume n: "n = 0"
   839   let ?f = "\<lambda>k z. inverse (of_nat (Suc k)) - inverse (z + of_nat k)"
   840   let ?F = "\<lambda>z. \<Sum>k. ?f k z" and ?f' = "\<lambda>k z. inverse ((z + of_nat k)\<^sup>2)"
   841   from no_nonpos_Int_in_ball'[OF z] guess d . note d = this
   842   from z have summable: "summable (\<lambda>k. inverse (of_nat (Suc k)) - inverse (z + of_nat k))"
   843     by (intro summable_Digamma) force
   844   from z have conv: "uniformly_convergent_on (ball z d) (\<lambda>k z. \<Sum>i<k. inverse ((z + of_nat i)\<^sup>2))"
   845     by (intro Polygamma_converges) auto
   846   with d have "summable (\<lambda>k. inverse ((z + of_nat k)\<^sup>2))" unfolding summable_iff_convergent
   847     by (auto dest!: uniformly_convergent_imp_convergent simp: summable_iff_convergent )
   848 
   849   have "(?F has_field_derivative (\<Sum>k. ?f' k z)) (at z)"
   850   proof (rule has_field_derivative_series'[of "ball z d" _ _ z])
   851     fix k :: nat and t :: 'a assume t: "t \<in> ball z d"
   852     from t d(2)[of t] show "((\<lambda>z. ?f k z) has_field_derivative ?f' k t) (at t within ball z d)"
   853       by (auto intro!: derivative_eq_intros simp: power2_eq_square simp del: of_nat_Suc
   854                dest!: plus_of_nat_eq_0_imp elim!: nonpos_Ints_cases)
   855   qed (insert d(1) summable conv, (assumption|simp)+)
   856   with z show "(Polygamma n has_field_derivative Polygamma (Suc n) z) (at z)"
   857     unfolding Digamma_def [abs_def] Polygamma_def [abs_def] using n
   858     by (force simp: power2_eq_square intro!: derivative_eq_intros)
   859 next
   860   assume n: "n \<noteq> 0"
   861   from z have z': "z \<noteq> 0" by auto
   862   from no_nonpos_Int_in_ball'[OF z] guess d . note d = this
   863   def n' \<equiv> "Suc n"
   864   from n have n': "n' \<ge> 2" by (simp add: n'_def)
   865   have "((\<lambda>z. \<Sum>k. inverse ((z + of_nat k) ^ n')) has_field_derivative
   866                 (\<Sum>k. - of_nat n' * inverse ((z + of_nat k) ^ (n'+1)))) (at z)"
   867   proof (rule has_field_derivative_series'[of "ball z d" _ _ z])
   868     fix k :: nat and t :: 'a assume t: "t \<in> ball z d"
   869     with d have t': "t \<notin> \<int>\<^sub>\<le>\<^sub>0" "t \<noteq> 0" by auto
   870     show "((\<lambda>a. inverse ((a + of_nat k) ^ n')) has_field_derivative
   871                 - of_nat n' * inverse ((t + of_nat k) ^ (n'+1))) (at t within ball z d)" using t'
   872       by (fastforce intro!: derivative_eq_intros simp: divide_simps power_diff dest: plus_of_nat_eq_0_imp)
   873   next
   874     have "uniformly_convergent_on (ball z d)
   875               (\<lambda>k z. (- of_nat n' :: 'a) * (\<Sum>i<k. inverse ((z + of_nat i) ^ (n'+1))))"
   876       using z' n by (intro uniformly_convergent_mult Polygamma_converges) (simp_all add: n'_def)
   877     thus "uniformly_convergent_on (ball z d)
   878               (\<lambda>k z. \<Sum>i<k. - of_nat n' * inverse ((z + of_nat i :: 'a) ^ (n'+1)))"
   879       by (subst (asm) setsum_right_distrib) simp
   880   qed (insert Polygamma_converges'[OF z' n'] d, simp_all)
   881   also have "(\<Sum>k. - of_nat n' * inverse ((z + of_nat k) ^ (n' + 1))) =
   882                (- of_nat n') * (\<Sum>k. inverse ((z + of_nat k) ^ (n' + 1)))"
   883     using Polygamma_converges'[OF z', of "n'+1"] n' by (subst suminf_mult) simp_all
   884   finally have "((\<lambda>z. \<Sum>k. inverse ((z + of_nat k) ^ n')) has_field_derivative
   885                     - of_nat n' * (\<Sum>k. inverse ((z + of_nat k) ^ (n' + 1)))) (at z)" .
   886   from DERIV_cmult[OF this, of "(-1)^Suc n * fact n :: 'a"]
   887     show "(Polygamma n has_field_derivative Polygamma (Suc n) z) (at z)"
   888     unfolding n'_def Polygamma_def[abs_def] using n by (simp add: algebra_simps)
   889 qed
   890 
   891 declare has_field_derivative_Polygamma[THEN DERIV_chain2, derivative_intros]
   892 
   893 lemma isCont_Polygamma [continuous_intros]:
   894   fixes f :: "_ \<Rightarrow> 'a :: {real_normed_field,euclidean_space}"
   895   shows "isCont f z \<Longrightarrow> f z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> isCont (\<lambda>x. Polygamma n (f x)) z"
   896   by (rule isCont_o2[OF _  DERIV_isCont[OF has_field_derivative_Polygamma]])
   897 
   898 lemma continuous_on_Polygamma:
   899   "A \<inter> \<int>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> continuous_on A (Polygamma n :: _ \<Rightarrow> 'a :: {real_normed_field,euclidean_space})"
   900   by (intro continuous_at_imp_continuous_on isCont_Polygamma[OF continuous_ident] ballI) blast
   901 
   902 lemma isCont_ln_Gamma_complex [continuous_intros]:
   903   fixes f :: "'a::t2_space \<Rightarrow> complex"
   904   shows "isCont f z \<Longrightarrow> f z \<notin> \<real>\<^sub>\<le>\<^sub>0 \<Longrightarrow> isCont (\<lambda>z. ln_Gamma (f z)) z"
   905   by (rule isCont_o2[OF _  DERIV_isCont[OF has_field_derivative_ln_Gamma_complex]])
   906 
   907 lemma continuous_on_ln_Gamma_complex [continuous_intros]:
   908   fixes A :: "complex set"
   909   shows "A \<inter> \<real>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> continuous_on A ln_Gamma"
   910   by (intro continuous_at_imp_continuous_on ballI isCont_ln_Gamma_complex[OF continuous_ident])
   911      fastforce
   912 
   913 text \<open>
   914   We define a type class that captures all the fundamental properties of the inverse of the Gamma function
   915   and defines the Gamma function upon that. This allows us to instantiate the type class both for
   916   the reals and for the complex numbers with a minimal amount of proof duplication.
   917 \<close>
   918 
   919 class Gamma = real_normed_field + complete_space +
   920   fixes rGamma :: "'a \<Rightarrow> 'a"
   921   assumes rGamma_eq_zero_iff_aux: "rGamma z = 0 \<longleftrightarrow> (\<exists>n. z = - of_nat n)"
   922   assumes differentiable_rGamma_aux1:
   923     "(\<And>n. z \<noteq> - of_nat n) \<Longrightarrow>
   924      let d = (THE d. (\<lambda>n. \<Sum>k<n. inverse (of_nat (Suc k)) - inverse (z + of_nat k))
   925                \<longlonglongrightarrow> d) - scaleR euler_mascheroni 1
   926      in  filterlim (\<lambda>y. (rGamma y - rGamma z + rGamma z * d * (y - z)) /\<^sub>R
   927                         norm (y - z)) (nhds 0) (at z)"
   928   assumes differentiable_rGamma_aux2:
   929     "let z = - of_nat n
   930      in  filterlim (\<lambda>y. (rGamma y - rGamma z - (-1)^n * (setprod of_nat {1..n}) * (y - z)) /\<^sub>R
   931                         norm (y - z)) (nhds 0) (at z)"
   932   assumes rGamma_series_aux: "(\<And>n. z \<noteq> - of_nat n) \<Longrightarrow>
   933              let fact' = (\<lambda>n. setprod of_nat {1..n});
   934                  exp = (\<lambda>x. THE e. (\<lambda>n. \<Sum>k<n. x^k /\<^sub>R fact k) \<longlonglongrightarrow> e);
   935                  pochhammer' = (\<lambda>a n. (\<Prod>n = 0..n. a + of_nat n))
   936              in  filterlim (\<lambda>n. pochhammer' z n / (fact' n * exp (z * (ln (of_nat n) *\<^sub>R 1))))
   937                      (nhds (rGamma z)) sequentially"
   938 begin
   939 subclass banach ..
   940 end
   941 
   942 definition "Gamma z = inverse (rGamma z)"
   943 
   944 
   945 subsection \<open>Basic properties\<close>
   946 
   947 lemma Gamma_nonpos_Int: "z \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Gamma z = 0"
   948   and rGamma_nonpos_Int: "z \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> rGamma z = 0"
   949   using rGamma_eq_zero_iff_aux[of z] unfolding Gamma_def by (auto elim!: nonpos_Ints_cases')
   950 
   951 lemma Gamma_nonzero: "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Gamma z \<noteq> 0"
   952   and rGamma_nonzero: "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> rGamma z \<noteq> 0"
   953   using rGamma_eq_zero_iff_aux[of z] unfolding Gamma_def by (auto elim!: nonpos_Ints_cases')
   954 
   955 lemma Gamma_eq_zero_iff: "Gamma z = 0 \<longleftrightarrow> z \<in> \<int>\<^sub>\<le>\<^sub>0"
   956   and rGamma_eq_zero_iff: "rGamma z = 0 \<longleftrightarrow> z \<in> \<int>\<^sub>\<le>\<^sub>0"
   957   using rGamma_eq_zero_iff_aux[of z] unfolding Gamma_def by (auto elim!: nonpos_Ints_cases')
   958 
   959 lemma rGamma_inverse_Gamma: "rGamma z = inverse (Gamma z)"
   960   unfolding Gamma_def by simp
   961 
   962 lemma rGamma_series_LIMSEQ [tendsto_intros]:
   963   "rGamma_series z \<longlonglongrightarrow> rGamma z"
   964 proof (cases "z \<in> \<int>\<^sub>\<le>\<^sub>0")
   965   case False
   966   hence "z \<noteq> - of_nat n" for n by auto
   967   from rGamma_series_aux[OF this] show ?thesis
   968     by (simp add: rGamma_series_def[abs_def] fact_altdef pochhammer_Suc_setprod
   969                   exp_def of_real_def[symmetric] suminf_def sums_def[abs_def])
   970 qed (insert rGamma_eq_zero_iff[of z], simp_all add: rGamma_series_nonpos_Ints_LIMSEQ)
   971 
   972 lemma Gamma_series_LIMSEQ [tendsto_intros]:
   973   "Gamma_series z \<longlonglongrightarrow> Gamma z"
   974 proof (cases "z \<in> \<int>\<^sub>\<le>\<^sub>0")
   975   case False
   976   hence "(\<lambda>n. inverse (rGamma_series z n)) \<longlonglongrightarrow> inverse (rGamma z)"
   977     by (intro tendsto_intros) (simp_all add: rGamma_eq_zero_iff)
   978   also have "(\<lambda>n. inverse (rGamma_series z n)) = Gamma_series z"
   979     by (simp add: rGamma_series_def Gamma_series_def[abs_def])
   980   finally show ?thesis by (simp add: Gamma_def)
   981 qed (insert Gamma_eq_zero_iff[of z], simp_all add: Gamma_series_nonpos_Ints_LIMSEQ)
   982 
   983 lemma Gamma_altdef: "Gamma z = lim (Gamma_series z)"
   984   using Gamma_series_LIMSEQ[of z] by (simp add: limI)
   985 
   986 lemma rGamma_1 [simp]: "rGamma 1 = 1"
   987 proof -
   988   have A: "eventually (\<lambda>n. rGamma_series 1 n = of_nat (Suc n) / of_nat n) sequentially"
   989     using eventually_gt_at_top[of "0::nat"]
   990     by (force elim!: eventually_mono simp: rGamma_series_def exp_of_real pochhammer_fact
   991                     divide_simps pochhammer_rec' dest!: pochhammer_eq_0_imp_nonpos_Int)
   992   have "rGamma_series 1 \<longlonglongrightarrow> 1" by (subst tendsto_cong[OF A]) (rule LIMSEQ_Suc_n_over_n)
   993   moreover have "rGamma_series 1 \<longlonglongrightarrow> rGamma 1" by (rule tendsto_intros)
   994   ultimately show ?thesis by (intro LIMSEQ_unique)
   995 qed
   996 
   997 lemma rGamma_plus1: "z * rGamma (z + 1) = rGamma z"
   998 proof -
   999   let ?f = "\<lambda>n. (z + 1) * inverse (of_nat n) + 1"
  1000   have "eventually (\<lambda>n. ?f n * rGamma_series z n = z * rGamma_series (z + 1) n) sequentially"
  1001     using eventually_gt_at_top[of "0::nat"]
  1002   proof eventually_elim
  1003     fix n :: nat assume n: "n > 0"
  1004     hence "z * rGamma_series (z + 1) n = inverse (of_nat n) *
  1005              pochhammer z (Suc (Suc n)) / (fact n * exp (z * of_real (ln (of_nat n))))"
  1006       by (subst pochhammer_rec) (simp add: rGamma_series_def field_simps exp_add exp_of_real)
  1007     also from n have "\<dots> = ?f n * rGamma_series z n"
  1008       by (subst pochhammer_rec') (simp_all add: divide_simps rGamma_series_def add_ac)
  1009     finally show "?f n * rGamma_series z n = z * rGamma_series (z + 1) n" ..
  1010   qed
  1011   moreover have "(\<lambda>n. ?f n * rGamma_series z n) \<longlonglongrightarrow> ((z+1) * 0 + 1) * rGamma z"
  1012     by (intro tendsto_intros lim_inverse_n)
  1013   hence "(\<lambda>n. ?f n * rGamma_series z n) \<longlonglongrightarrow> rGamma z" by simp
  1014   ultimately have "(\<lambda>n. z * rGamma_series (z + 1) n) \<longlonglongrightarrow> rGamma z"
  1015     by (rule Lim_transform_eventually)
  1016   moreover have "(\<lambda>n. z * rGamma_series (z + 1) n) \<longlonglongrightarrow> z * rGamma (z + 1)"
  1017     by (intro tendsto_intros)
  1018   ultimately show "z * rGamma (z + 1) = rGamma z" using LIMSEQ_unique by blast
  1019 qed
  1020 
  1021 
  1022 lemma pochhammer_rGamma: "rGamma z = pochhammer z n * rGamma (z + of_nat n)"
  1023 proof (induction n arbitrary: z)
  1024   case (Suc n z)
  1025   have "rGamma z = pochhammer z n * rGamma (z + of_nat n)" by (rule Suc.IH)
  1026   also note rGamma_plus1 [symmetric]
  1027   finally show ?case by (simp add: add_ac pochhammer_rec')
  1028 qed simp_all
  1029 
  1030 lemma Gamma_plus1: "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Gamma (z + 1) = z * Gamma z"
  1031   using rGamma_plus1[of z] by (simp add: rGamma_inverse_Gamma field_simps Gamma_eq_zero_iff)
  1032 
  1033 lemma pochhammer_Gamma: "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> pochhammer z n = Gamma (z + of_nat n) / Gamma z"
  1034   using pochhammer_rGamma[of z]
  1035   by (simp add: rGamma_inverse_Gamma Gamma_eq_zero_iff field_simps)
  1036 
  1037 lemma Gamma_0 [simp]: "Gamma 0 = 0"
  1038   and rGamma_0 [simp]: "rGamma 0 = 0"
  1039   and Gamma_neg_1 [simp]: "Gamma (- 1) = 0"
  1040   and rGamma_neg_1 [simp]: "rGamma (- 1) = 0"
  1041   and Gamma_neg_numeral [simp]: "Gamma (- numeral n) = 0"
  1042   and rGamma_neg_numeral [simp]: "rGamma (- numeral n) = 0"
  1043   and Gamma_neg_of_nat [simp]: "Gamma (- of_nat m) = 0"
  1044   and rGamma_neg_of_nat [simp]: "rGamma (- of_nat m) = 0"
  1045   by (simp_all add: rGamma_eq_zero_iff Gamma_eq_zero_iff)
  1046 
  1047 lemma Gamma_1 [simp]: "Gamma 1 = 1" unfolding Gamma_def by simp
  1048 
  1049 lemma Gamma_fact: "Gamma (of_nat (Suc n)) = fact n"
  1050   by (simp add: pochhammer_fact pochhammer_Gamma of_nat_in_nonpos_Ints_iff)
  1051 
  1052 lemma Gamma_numeral: "Gamma (numeral n) = fact (pred_numeral n)"
  1053   by (subst of_nat_numeral[symmetric], subst numeral_eq_Suc, subst Gamma_fact) (rule refl)
  1054 
  1055 lemma Gamma_of_int: "Gamma (of_int n) = (if n > 0 then fact (nat (n - 1)) else 0)"
  1056 proof (cases "n > 0")
  1057   case True
  1058   hence "Gamma (of_int n) = Gamma (of_nat (Suc (nat (n - 1))))" by (subst of_nat_Suc) simp_all
  1059   with True show ?thesis by (subst (asm) Gamma_fact) simp
  1060 qed (simp_all add: Gamma_eq_zero_iff nonpos_Ints_of_int)
  1061 
  1062 lemma rGamma_of_int: "rGamma (of_int n) = (if n > 0 then inverse (fact (nat (n - 1))) else 0)"
  1063   by (simp add: Gamma_of_int rGamma_inverse_Gamma)
  1064 
  1065 lemma Gamma_seriesI:
  1066   assumes "(\<lambda>n. g n / Gamma_series z n) \<longlonglongrightarrow> 1"
  1067   shows   "g \<longlonglongrightarrow> Gamma z"
  1068 proof (rule Lim_transform_eventually)
  1069   have "1/2 > (0::real)" by simp
  1070   from tendstoD[OF assms, OF this]
  1071     show "eventually (\<lambda>n. g n / Gamma_series z n * Gamma_series z n = g n) sequentially"
  1072     by (force elim!: eventually_mono simp: dist_real_def dist_0_norm)
  1073   from assms have "(\<lambda>n. g n / Gamma_series z n * Gamma_series z n) \<longlonglongrightarrow> 1 * Gamma z"
  1074     by (intro tendsto_intros)
  1075   thus "(\<lambda>n. g n / Gamma_series z n * Gamma_series z n) \<longlonglongrightarrow> Gamma z" by simp
  1076 qed
  1077 
  1078 lemma Gamma_seriesI':
  1079   assumes "f \<longlonglongrightarrow> rGamma z"
  1080   assumes "(\<lambda>n. g n * f n) \<longlonglongrightarrow> 1"
  1081   assumes "z \<notin> \<int>\<^sub>\<le>\<^sub>0"
  1082   shows   "g \<longlonglongrightarrow> Gamma z"
  1083 proof (rule Lim_transform_eventually)
  1084   have "1/2 > (0::real)" by simp
  1085   from tendstoD[OF assms(2), OF this] show "eventually (\<lambda>n. g n * f n / f n = g n) sequentially"
  1086     by (force elim!: eventually_mono simp: dist_real_def dist_0_norm)
  1087   from assms have "(\<lambda>n. g n * f n / f n) \<longlonglongrightarrow> 1 / rGamma z"
  1088     by (intro tendsto_divide assms) (simp_all add: rGamma_eq_zero_iff)
  1089   thus "(\<lambda>n. g n * f n / f n) \<longlonglongrightarrow> Gamma z" by (simp add: Gamma_def divide_inverse)
  1090 qed
  1091 
  1092 lemma Gamma_series'_LIMSEQ: "Gamma_series' z \<longlonglongrightarrow> Gamma z"
  1093   by (cases "z \<in> \<int>\<^sub>\<le>\<^sub>0") (simp_all add: Gamma_nonpos_Int Gamma_seriesI[OF Gamma_series_Gamma_series']
  1094                                       Gamma_series'_nonpos_Ints_LIMSEQ[of z])
  1095 
  1096 
  1097 subsection \<open>Differentiability\<close>
  1098 
  1099 lemma has_field_derivative_rGamma_no_nonpos_int:
  1100   assumes "z \<notin> \<int>\<^sub>\<le>\<^sub>0"
  1101   shows   "(rGamma has_field_derivative -rGamma z * Digamma z) (at z within A)"
  1102 proof (rule has_field_derivative_at_within)
  1103   from assms have "z \<noteq> - of_nat n" for n by auto
  1104   from differentiable_rGamma_aux1[OF this]
  1105     show "(rGamma has_field_derivative -rGamma z * Digamma z) (at z)"
  1106          unfolding Digamma_def suminf_def sums_def[abs_def]
  1107                    has_field_derivative_def has_derivative_def netlimit_at
  1108     by (simp add: Let_def bounded_linear_mult_right mult_ac of_real_def [symmetric])
  1109 qed
  1110 
  1111 lemma has_field_derivative_rGamma_nonpos_int:
  1112   "(rGamma has_field_derivative (-1)^n * fact n) (at (- of_nat n) within A)"
  1113   apply (rule has_field_derivative_at_within)
  1114   using differentiable_rGamma_aux2[of n]
  1115   unfolding Let_def has_field_derivative_def has_derivative_def netlimit_at
  1116   by (simp only: bounded_linear_mult_right mult_ac of_real_def [symmetric] fact_altdef)
  1117 
  1118 lemma has_field_derivative_rGamma [derivative_intros]:
  1119   "(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>)
  1120       else -rGamma z * Digamma z)) (at z within A)"
  1121 using has_field_derivative_rGamma_no_nonpos_int[of z A]
  1122       has_field_derivative_rGamma_nonpos_int[of "nat \<lfloor>norm z\<rfloor>" A]
  1123   by (auto elim!: nonpos_Ints_cases')
  1124 
  1125 declare has_field_derivative_rGamma_no_nonpos_int [THEN DERIV_chain2, derivative_intros]
  1126 declare has_field_derivative_rGamma [THEN DERIV_chain2, derivative_intros]
  1127 declare has_field_derivative_rGamma_nonpos_int [derivative_intros]
  1128 declare has_field_derivative_rGamma_no_nonpos_int [derivative_intros]
  1129 declare has_field_derivative_rGamma [derivative_intros]
  1130 
  1131 
  1132 lemma has_field_derivative_Gamma [derivative_intros]:
  1133   "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> (Gamma has_field_derivative Gamma z * Digamma z) (at z within A)"
  1134   unfolding Gamma_def [abs_def]
  1135   by (fastforce intro!: derivative_eq_intros simp: rGamma_eq_zero_iff)
  1136 
  1137 declare has_field_derivative_Gamma[THEN DERIV_chain2, derivative_intros]
  1138 
  1139 (* TODO: Hide ugly facts properly *)
  1140 hide_fact rGamma_eq_zero_iff_aux differentiable_rGamma_aux1 differentiable_rGamma_aux2
  1141           differentiable_rGamma_aux2 rGamma_series_aux Gamma_class.rGamma_eq_zero_iff_aux
  1142 
  1143 
  1144 
  1145 (* TODO: differentiable etc. *)
  1146 
  1147 
  1148 subsection \<open>Continuity\<close>
  1149 
  1150 lemma continuous_on_rGamma [continuous_intros]: "continuous_on A rGamma"
  1151   by (rule DERIV_continuous_on has_field_derivative_rGamma)+
  1152 
  1153 lemma continuous_on_Gamma [continuous_intros]: "A \<inter> \<int>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> continuous_on A Gamma"
  1154   by (rule DERIV_continuous_on has_field_derivative_Gamma)+ blast
  1155 
  1156 lemma isCont_rGamma [continuous_intros]:
  1157   "isCont f z \<Longrightarrow> isCont (\<lambda>x. rGamma (f x)) z"
  1158   by (rule isCont_o2[OF _  DERIV_isCont[OF has_field_derivative_rGamma]])
  1159 
  1160 lemma isCont_Gamma [continuous_intros]:
  1161   "isCont f z \<Longrightarrow> f z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> isCont (\<lambda>x. Gamma (f x)) z"
  1162   by (rule isCont_o2[OF _  DERIV_isCont[OF has_field_derivative_Gamma]])
  1163 
  1164 
  1165 
  1166 text \<open>The complex Gamma function\<close>
  1167 
  1168 instantiation complex :: Gamma
  1169 begin
  1170 
  1171 definition rGamma_complex :: "complex \<Rightarrow> complex" where
  1172   "rGamma_complex z = lim (rGamma_series z)"
  1173 
  1174 lemma rGamma_series_complex_converges:
  1175         "convergent (rGamma_series (z :: complex))" (is "?thesis1")
  1176   and rGamma_complex_altdef:
  1177         "rGamma z = (if z \<in> \<int>\<^sub>\<le>\<^sub>0 then 0 else exp (-ln_Gamma z))" (is "?thesis2")
  1178 proof -
  1179   have "?thesis1 \<and> ?thesis2"
  1180   proof (cases "z \<in> \<int>\<^sub>\<le>\<^sub>0")
  1181     case False
  1182     have "rGamma_series z \<longlonglongrightarrow> exp (- ln_Gamma z)"
  1183     proof (rule Lim_transform_eventually)
  1184       from ln_Gamma_series_complex_converges'[OF False] guess d by (elim exE conjE)
  1185       from this(1) uniformly_convergent_imp_convergent[OF this(2), of z]
  1186         have "ln_Gamma_series z \<longlonglongrightarrow> lim (ln_Gamma_series z)" by (simp add: convergent_LIMSEQ_iff)
  1187       thus "(\<lambda>n. exp (-ln_Gamma_series z n)) \<longlonglongrightarrow> exp (- ln_Gamma z)"
  1188         unfolding convergent_def ln_Gamma_def by (intro tendsto_exp tendsto_minus)
  1189       from eventually_gt_at_top[of "0::nat"] exp_ln_Gamma_series_complex False
  1190         show "eventually (\<lambda>n. exp (-ln_Gamma_series z n) = rGamma_series z n) sequentially"
  1191         by (force elim!: eventually_mono simp: exp_minus Gamma_series_def rGamma_series_def)
  1192     qed
  1193     with False show ?thesis
  1194       by (auto simp: convergent_def rGamma_complex_def intro!: limI)
  1195   next
  1196     case True
  1197     then obtain k where "z = - of_nat k" by (erule nonpos_Ints_cases')
  1198     also have "rGamma_series \<dots> \<longlonglongrightarrow> 0"
  1199       by (subst tendsto_cong[OF rGamma_series_minus_of_nat]) (simp_all add: convergent_const)
  1200     finally show ?thesis using True
  1201       by (auto simp: rGamma_complex_def convergent_def intro!: limI)
  1202   qed
  1203   thus "?thesis1" "?thesis2" by blast+
  1204 qed
  1205 
  1206 context
  1207 begin
  1208 
  1209 (* TODO: duplication *)
  1210 private lemma rGamma_complex_plus1: "z * rGamma (z + 1) = rGamma (z :: complex)"
  1211 proof -
  1212   let ?f = "\<lambda>n. (z + 1) * inverse (of_nat n) + 1"
  1213   have "eventually (\<lambda>n. ?f n * rGamma_series z n = z * rGamma_series (z + 1) n) sequentially"
  1214     using eventually_gt_at_top[of "0::nat"]
  1215   proof eventually_elim
  1216     fix n :: nat assume n: "n > 0"
  1217     hence "z * rGamma_series (z + 1) n = inverse (of_nat n) *
  1218              pochhammer z (Suc (Suc n)) / (fact n * exp (z * of_real (ln (of_nat n))))"
  1219       by (subst pochhammer_rec) (simp add: rGamma_series_def field_simps exp_add exp_of_real)
  1220     also from n have "\<dots> = ?f n * rGamma_series z n"
  1221       by (subst pochhammer_rec') (simp_all add: divide_simps rGamma_series_def add_ac)
  1222     finally show "?f n * rGamma_series z n = z * rGamma_series (z + 1) n" ..
  1223   qed
  1224   moreover have "(\<lambda>n. ?f n * rGamma_series z n) \<longlonglongrightarrow> ((z+1) * 0 + 1) * rGamma z"
  1225     using rGamma_series_complex_converges
  1226     by (intro tendsto_intros lim_inverse_n)
  1227        (simp_all add: convergent_LIMSEQ_iff rGamma_complex_def)
  1228   hence "(\<lambda>n. ?f n * rGamma_series z n) \<longlonglongrightarrow> rGamma z" by simp
  1229   ultimately have "(\<lambda>n. z * rGamma_series (z + 1) n) \<longlonglongrightarrow> rGamma z"
  1230     by (rule Lim_transform_eventually)
  1231   moreover have "(\<lambda>n. z * rGamma_series (z + 1) n) \<longlonglongrightarrow> z * rGamma (z + 1)"
  1232     using rGamma_series_complex_converges
  1233     by (auto intro!: tendsto_mult simp: rGamma_complex_def convergent_LIMSEQ_iff)
  1234   ultimately show "z * rGamma (z + 1) = rGamma z" using LIMSEQ_unique by blast
  1235 qed
  1236 
  1237 private lemma has_field_derivative_rGamma_complex_no_nonpos_Int:
  1238   assumes "(z :: complex) \<notin> \<int>\<^sub>\<le>\<^sub>0"
  1239   shows   "(rGamma has_field_derivative - rGamma z * Digamma z) (at z)"
  1240 proof -
  1241   have diff: "(rGamma has_field_derivative - rGamma z * Digamma z) (at z)" if "Re z > 0" for z
  1242   proof (subst DERIV_cong_ev[OF refl _ refl])
  1243     from that have "eventually (\<lambda>t. t \<in> ball z (Re z/2)) (nhds z)"
  1244       by (intro eventually_nhds_in_nhd) simp_all
  1245     thus "eventually (\<lambda>t. rGamma t = exp (- ln_Gamma t)) (nhds z)"
  1246       using no_nonpos_Int_in_ball_complex[OF that]
  1247       by (auto elim!: eventually_mono simp: rGamma_complex_altdef)
  1248   next
  1249     have "z \<notin> \<real>\<^sub>\<le>\<^sub>0" using that by (simp add: complex_nonpos_Reals_iff)
  1250     with that show "((\<lambda>t. exp (- ln_Gamma t)) has_field_derivative (-rGamma z * Digamma z)) (at z)"
  1251      by (force elim!: nonpos_Ints_cases intro!: derivative_eq_intros simp: rGamma_complex_altdef)
  1252   qed
  1253 
  1254   from assms show "(rGamma has_field_derivative - rGamma z * Digamma z) (at z)"
  1255   proof (induction "nat \<lfloor>1 - Re z\<rfloor>" arbitrary: z)
  1256     case (Suc n z)
  1257     from Suc.prems have z: "z \<noteq> 0" by auto
  1258     from Suc.hyps have "n = nat \<lfloor>- Re z\<rfloor>" by linarith
  1259     hence A: "n = nat \<lfloor>1 - Re (z + 1)\<rfloor>" by simp
  1260     from Suc.prems have B: "z + 1 \<notin> \<int>\<^sub>\<le>\<^sub>0" by (force dest: plus_one_in_nonpos_Ints_imp)
  1261 
  1262     have "((\<lambda>z. z * (rGamma \<circ> (\<lambda>z. z + 1)) z) has_field_derivative
  1263       -rGamma (z + 1) * (Digamma (z + 1) * z - 1)) (at z)"
  1264       by (rule derivative_eq_intros DERIV_chain Suc refl A B)+ (simp add: algebra_simps)
  1265     also have "(\<lambda>z. z * (rGamma \<circ> (\<lambda>z. z + 1 :: complex)) z) = rGamma"
  1266       by (simp add: rGamma_complex_plus1)
  1267     also from z have "Digamma (z + 1) * z - 1 = z * Digamma z"
  1268       by (subst Digamma_plus1) (simp_all add: field_simps)
  1269     also have "-rGamma (z + 1) * (z * Digamma z) = -rGamma z * Digamma z"
  1270       by (simp add: rGamma_complex_plus1[of z, symmetric])
  1271     finally show ?case .
  1272   qed (intro diff, simp)
  1273 qed
  1274 
  1275 private lemma rGamma_complex_1: "rGamma (1 :: complex) = 1"
  1276 proof -
  1277   have A: "eventually (\<lambda>n. rGamma_series 1 n = of_nat (Suc n) / of_nat n) sequentially"
  1278     using eventually_gt_at_top[of "0::nat"]
  1279     by (force elim!: eventually_mono simp: rGamma_series_def exp_of_real pochhammer_fact
  1280                     divide_simps pochhammer_rec' dest!: pochhammer_eq_0_imp_nonpos_Int)
  1281   have "rGamma_series 1 \<longlonglongrightarrow> 1" by (subst tendsto_cong[OF A]) (rule LIMSEQ_Suc_n_over_n)
  1282   thus "rGamma 1 = (1 :: complex)" unfolding rGamma_complex_def by (rule limI)
  1283 qed
  1284 
  1285 private lemma has_field_derivative_rGamma_complex_nonpos_Int:
  1286   "(rGamma has_field_derivative (-1)^n * fact n) (at (- of_nat n :: complex))"
  1287 proof (induction n)
  1288   case 0
  1289   have A: "(0::complex) + 1 \<notin> \<int>\<^sub>\<le>\<^sub>0" by simp
  1290   have "((\<lambda>z. z * (rGamma \<circ> (\<lambda>z. z + 1 :: complex)) z) has_field_derivative 1) (at 0)"
  1291     by (rule derivative_eq_intros DERIV_chain refl
  1292              has_field_derivative_rGamma_complex_no_nonpos_Int A)+ (simp add: rGamma_complex_1)
  1293     thus ?case by (simp add: rGamma_complex_plus1)
  1294 next
  1295   case (Suc n)
  1296   hence A: "(rGamma has_field_derivative (-1)^n * fact n)
  1297                 (at (- of_nat (Suc n) + 1 :: complex))" by simp
  1298    have "((\<lambda>z. z * (rGamma \<circ> (\<lambda>z. z + 1 :: complex)) z) has_field_derivative
  1299              (- 1) ^ Suc n * fact (Suc n)) (at (- of_nat (Suc n)))"
  1300      by (rule derivative_eq_intros refl A DERIV_chain)+
  1301         (simp add: algebra_simps rGamma_complex_altdef)
  1302   thus ?case by (simp add: rGamma_complex_plus1)
  1303 qed
  1304 
  1305 instance proof
  1306   fix z :: complex show "(rGamma z = 0) \<longleftrightarrow> (\<exists>n. z = - of_nat n)"
  1307     by (auto simp: rGamma_complex_altdef elim!: nonpos_Ints_cases')
  1308 next
  1309   fix z :: complex assume "\<And>n. z \<noteq> - of_nat n"
  1310   hence "z \<notin> \<int>\<^sub>\<le>\<^sub>0" by (auto elim!: nonpos_Ints_cases')
  1311   from has_field_derivative_rGamma_complex_no_nonpos_Int[OF this]
  1312     show "let d = (THE d. (\<lambda>n. \<Sum>k<n. inverse (of_nat (Suc k)) - inverse (z + of_nat k))
  1313                        \<longlonglongrightarrow> d) - euler_mascheroni *\<^sub>R 1 in  (\<lambda>y. (rGamma y - rGamma z +
  1314               rGamma z * d * (y - z)) /\<^sub>R  cmod (y - z)) \<midarrow>z\<rightarrow> 0"
  1315       by (simp add: has_field_derivative_def has_derivative_def Digamma_def sums_def [abs_def]
  1316                     netlimit_at of_real_def[symmetric] suminf_def)
  1317 next
  1318   fix n :: nat
  1319   from has_field_derivative_rGamma_complex_nonpos_Int[of n]
  1320   show "let z = - of_nat n in (\<lambda>y. (rGamma y - rGamma z - (- 1) ^ n * setprod of_nat {1..n} *
  1321                   (y - z)) /\<^sub>R cmod (y - z)) \<midarrow>z\<rightarrow> 0"
  1322     by (simp add: has_field_derivative_def has_derivative_def fact_altdef netlimit_at Let_def)
  1323 next
  1324   fix z :: complex
  1325   from rGamma_series_complex_converges[of z] have "rGamma_series z \<longlonglongrightarrow> rGamma z"
  1326     by (simp add: convergent_LIMSEQ_iff rGamma_complex_def)
  1327   thus "let fact' = \<lambda>n. setprod of_nat {1..n};
  1328             exp = \<lambda>x. THE e. (\<lambda>n. \<Sum>k<n. x ^ k /\<^sub>R fact k) \<longlonglongrightarrow> e;
  1329             pochhammer' = \<lambda>a n. \<Prod>n = 0..n. a + of_nat n
  1330         in  (\<lambda>n. pochhammer' z n / (fact' n * exp (z * ln (real_of_nat n) *\<^sub>R 1))) \<longlonglongrightarrow> rGamma z"
  1331     by (simp add: fact_altdef pochhammer_Suc_setprod rGamma_series_def [abs_def] exp_def
  1332                   of_real_def [symmetric] suminf_def sums_def [abs_def])
  1333 qed
  1334 
  1335 end
  1336 end
  1337 
  1338 
  1339 lemma Gamma_complex_altdef:
  1340   "Gamma z = (if z \<in> \<int>\<^sub>\<le>\<^sub>0 then 0 else exp (ln_Gamma (z :: complex)))"
  1341   unfolding Gamma_def rGamma_complex_altdef by (simp add: exp_minus)
  1342 
  1343 lemma cnj_rGamma: "cnj (rGamma z) = rGamma (cnj z)"
  1344 proof -
  1345   have "rGamma_series (cnj z) = (\<lambda>n. cnj (rGamma_series z n))"
  1346     by (intro ext) (simp_all add: rGamma_series_def exp_cnj)
  1347   also have "... \<longlonglongrightarrow> cnj (rGamma z)" by (intro tendsto_cnj tendsto_intros)
  1348   finally show ?thesis unfolding rGamma_complex_def by (intro sym[OF limI])
  1349 qed
  1350 
  1351 lemma cnj_Gamma: "cnj (Gamma z) = Gamma (cnj z)"
  1352   unfolding Gamma_def by (simp add: cnj_rGamma)
  1353 
  1354 lemma Gamma_complex_real:
  1355   "z \<in> \<real> \<Longrightarrow> Gamma z \<in> (\<real> :: complex set)" and rGamma_complex_real: "z \<in> \<real> \<Longrightarrow> rGamma z \<in> \<real>"
  1356   by (simp_all add: Reals_cnj_iff cnj_Gamma cnj_rGamma)
  1357 
  1358 lemma complex_differentiable_rGamma: "rGamma complex_differentiable (at z within A)"
  1359   using has_field_derivative_rGamma[of z] unfolding complex_differentiable_def by blast
  1360 
  1361 lemma holomorphic_on_rGamma: "rGamma holomorphic_on A"
  1362   unfolding holomorphic_on_def by (auto intro!: complex_differentiable_rGamma)
  1363 
  1364 lemma analytic_on_rGamma: "rGamma analytic_on A"
  1365   unfolding analytic_on_def by (auto intro!: exI[of _ 1] holomorphic_on_rGamma)
  1366 
  1367 
  1368 lemma complex_differentiable_Gamma: "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Gamma complex_differentiable (at z within A)"
  1369   using has_field_derivative_Gamma[of z] unfolding complex_differentiable_def by auto
  1370 
  1371 lemma holomorphic_on_Gamma: "A \<inter> \<int>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> Gamma holomorphic_on A"
  1372   unfolding holomorphic_on_def by (auto intro!: complex_differentiable_Gamma)
  1373 
  1374 lemma analytic_on_Gamma: "A \<inter> \<int>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> Gamma analytic_on A"
  1375   by (rule analytic_on_subset[of _ "UNIV - \<int>\<^sub>\<le>\<^sub>0"], subst analytic_on_open)
  1376      (auto intro!: holomorphic_on_Gamma)
  1377 
  1378 lemma has_field_derivative_rGamma_complex' [derivative_intros]:
  1379   "(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
  1380         -rGamma z * Digamma z)) (at z within A)"
  1381   using has_field_derivative_rGamma[of z] by (auto elim!: nonpos_Ints_cases')
  1382 
  1383 declare has_field_derivative_rGamma_complex'[THEN DERIV_chain2, derivative_intros]
  1384 
  1385 
  1386 lemma complex_differentiable_Polygamma:
  1387   "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Polygamma n complex_differentiable (at z within A)"
  1388   using has_field_derivative_Polygamma[of z n] unfolding complex_differentiable_def by auto
  1389 
  1390 lemma holomorphic_on_Polygamma: "A \<inter> \<int>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> Polygamma n holomorphic_on A"
  1391   unfolding holomorphic_on_def by (auto intro!: complex_differentiable_Polygamma)
  1392 
  1393 lemma analytic_on_Polygamma: "A \<inter> \<int>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> Polygamma n analytic_on A"
  1394   by (rule analytic_on_subset[of _ "UNIV - \<int>\<^sub>\<le>\<^sub>0"], subst analytic_on_open)
  1395      (auto intro!: holomorphic_on_Polygamma)
  1396 
  1397 
  1398 
  1399 text \<open>The real Gamma function\<close>
  1400 
  1401 lemma rGamma_series_real:
  1402   "eventually (\<lambda>n. rGamma_series x n = Re (rGamma_series (of_real x) n)) sequentially"
  1403   using eventually_gt_at_top[of "0 :: nat"]
  1404 proof eventually_elim
  1405   fix n :: nat assume n: "n > 0"
  1406   have "Re (rGamma_series (of_real x) n) =
  1407           Re (of_real (pochhammer x (Suc n)) / (fact n * exp (of_real (x * ln (real_of_nat n)))))"
  1408     using n by (simp add: rGamma_series_def powr_def Ln_of_nat pochhammer_of_real)
  1409   also from n have "\<dots> = Re (of_real ((pochhammer x (Suc n)) /
  1410                               (fact n * (exp (x * ln (real_of_nat n))))))"
  1411     by (subst exp_of_real) simp
  1412   also from n have "\<dots> = rGamma_series x n"
  1413     by (subst Re_complex_of_real) (simp add: rGamma_series_def powr_def)
  1414   finally show "rGamma_series x n = Re (rGamma_series (of_real x) n)" ..
  1415 qed
  1416 
  1417 instantiation real :: Gamma
  1418 begin
  1419 
  1420 definition "rGamma_real x = Re (rGamma (of_real x :: complex))"
  1421 
  1422 instance proof
  1423   fix x :: real
  1424   have "rGamma x = Re (rGamma (of_real x))" by (simp add: rGamma_real_def)
  1425   also have "of_real \<dots> = rGamma (of_real x :: complex)"
  1426     by (intro of_real_Re rGamma_complex_real) simp_all
  1427   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)
  1428   also have "\<dots> \<longleftrightarrow> (\<exists>n. x = - of_nat n)" by (auto elim!: nonpos_Ints_cases')
  1429   finally show "(rGamma x) = 0 \<longleftrightarrow> (\<exists>n. x = - real_of_nat n)" by simp
  1430 next
  1431   fix x :: real assume "\<And>n. x \<noteq> - of_nat n"
  1432   hence "complex_of_real x \<notin> \<int>\<^sub>\<le>\<^sub>0"
  1433     by (subst of_real_in_nonpos_Ints_iff) (auto elim!: nonpos_Ints_cases')
  1434   moreover from this have "x \<noteq> 0" by auto
  1435   ultimately have "(rGamma has_field_derivative - rGamma x * Digamma x) (at x)"
  1436     by (fastforce intro!: derivative_eq_intros has_vector_derivative_real_complex
  1437                   simp: Polygamma_of_real rGamma_real_def [abs_def])
  1438   thus "let d = (THE d. (\<lambda>n. \<Sum>k<n. inverse (of_nat (Suc k)) - inverse (x + of_nat k))
  1439                        \<longlonglongrightarrow> d) - euler_mascheroni *\<^sub>R 1 in  (\<lambda>y. (rGamma y - rGamma x +
  1440               rGamma x * d * (y - x)) /\<^sub>R  norm (y - x)) \<midarrow>x\<rightarrow> 0"
  1441       by (simp add: has_field_derivative_def has_derivative_def Digamma_def sums_def [abs_def]
  1442                     netlimit_at of_real_def[symmetric] suminf_def)
  1443 next
  1444   fix n :: nat
  1445   have "(rGamma has_field_derivative (-1)^n * fact n) (at (- of_nat n :: real))"
  1446     by (fastforce intro!: derivative_eq_intros has_vector_derivative_real_complex
  1447                   simp: Polygamma_of_real rGamma_real_def [abs_def])
  1448   thus "let x = - of_nat n in (\<lambda>y. (rGamma y - rGamma x - (- 1) ^ n * setprod of_nat {1..n} *
  1449                   (y - x)) /\<^sub>R norm (y - x)) \<midarrow>x::real\<rightarrow> 0"
  1450     by (simp add: has_field_derivative_def has_derivative_def fact_altdef netlimit_at Let_def)
  1451 next
  1452   fix x :: real
  1453   have "rGamma_series x \<longlonglongrightarrow> rGamma x"
  1454   proof (rule Lim_transform_eventually)
  1455     show "(\<lambda>n. Re (rGamma_series (of_real x) n)) \<longlonglongrightarrow> rGamma x" unfolding rGamma_real_def
  1456       by (intro tendsto_intros)
  1457   qed (insert rGamma_series_real, simp add: eq_commute)
  1458   thus "let fact' = \<lambda>n. setprod of_nat {1..n};
  1459             exp = \<lambda>x. THE e. (\<lambda>n. \<Sum>k<n. x ^ k /\<^sub>R fact k) \<longlonglongrightarrow> e;
  1460             pochhammer' = \<lambda>a n. \<Prod>n = 0..n. a + of_nat n
  1461         in  (\<lambda>n. pochhammer' x n / (fact' n * exp (x * ln (real_of_nat n) *\<^sub>R 1))) \<longlonglongrightarrow> rGamma x"
  1462     by (simp add: fact_altdef pochhammer_Suc_setprod rGamma_series_def [abs_def] exp_def
  1463                   of_real_def [symmetric] suminf_def sums_def [abs_def])
  1464 qed
  1465 
  1466 end
  1467 
  1468 
  1469 lemma rGamma_complex_of_real: "rGamma (complex_of_real x) = complex_of_real (rGamma x)"
  1470   unfolding rGamma_real_def using rGamma_complex_real by simp
  1471 
  1472 lemma Gamma_complex_of_real: "Gamma (complex_of_real x) = complex_of_real (Gamma x)"
  1473   unfolding Gamma_def by (simp add: rGamma_complex_of_real)
  1474 
  1475 lemma rGamma_real_altdef: "rGamma x = lim (rGamma_series (x :: real))"
  1476   by (rule sym, rule limI, rule tendsto_intros)
  1477 
  1478 lemma Gamma_real_altdef1: "Gamma x = lim (Gamma_series (x :: real))"
  1479   by (rule sym, rule limI, rule tendsto_intros)
  1480 
  1481 lemma Gamma_real_altdef2: "Gamma x = Re (Gamma (of_real x))"
  1482   using rGamma_complex_real[OF Reals_of_real[of x]]
  1483   by (elim Reals_cases)
  1484      (simp only: Gamma_def rGamma_real_def of_real_inverse[symmetric] Re_complex_of_real)
  1485 
  1486 lemma ln_Gamma_series_complex_of_real:
  1487   "x > 0 \<Longrightarrow> n > 0 \<Longrightarrow> ln_Gamma_series (complex_of_real x) n = of_real (ln_Gamma_series x n)"
  1488 proof -
  1489   assume xn: "x > 0" "n > 0"
  1490   have "Ln (complex_of_real x / of_nat k + 1) = of_real (ln (x / of_nat k + 1))" if "k \<ge> 1" for k
  1491     using that xn by (subst Ln_of_real [symmetric]) (auto intro!: add_nonneg_pos simp: field_simps)
  1492   with xn show ?thesis by (simp add: ln_Gamma_series_def Ln_of_nat Ln_of_real)
  1493 qed
  1494 
  1495 lemma ln_Gamma_real_converges:
  1496   assumes "(x::real) > 0"
  1497   shows   "convergent (ln_Gamma_series x)"
  1498 proof -
  1499   have "(\<lambda>n. ln_Gamma_series (complex_of_real x) n) \<longlonglongrightarrow> ln_Gamma (of_real x)" using assms
  1500     by (intro ln_Gamma_complex_LIMSEQ) (auto simp: of_real_in_nonpos_Ints_iff)
  1501   moreover from eventually_gt_at_top[of "0::nat"]
  1502     have "eventually (\<lambda>n. complex_of_real (ln_Gamma_series x n) =
  1503             ln_Gamma_series (complex_of_real x) n) sequentially"
  1504     by eventually_elim (simp add: ln_Gamma_series_complex_of_real assms)
  1505   ultimately have "(\<lambda>n. complex_of_real (ln_Gamma_series x n)) \<longlonglongrightarrow> ln_Gamma (of_real x)"
  1506     by (subst tendsto_cong) assumption+
  1507   from tendsto_Re[OF this] show ?thesis by (auto simp: convergent_def)
  1508 qed
  1509 
  1510 lemma ln_Gamma_real_LIMSEQ: "(x::real) > 0 \<Longrightarrow> ln_Gamma_series x \<longlonglongrightarrow> ln_Gamma x"
  1511   using ln_Gamma_real_converges[of x] unfolding ln_Gamma_def by (simp add: convergent_LIMSEQ_iff)
  1512 
  1513 lemma ln_Gamma_complex_of_real: "x > 0 \<Longrightarrow> ln_Gamma (complex_of_real x) = of_real (ln_Gamma x)"
  1514 proof (unfold ln_Gamma_def, rule limI, rule Lim_transform_eventually)
  1515   assume x: "x > 0"
  1516   show "eventually (\<lambda>n. of_real (ln_Gamma_series x n) =
  1517             ln_Gamma_series (complex_of_real x) n) sequentially"
  1518     using eventually_gt_at_top[of "0::nat"]
  1519     by eventually_elim (simp add: ln_Gamma_series_complex_of_real x)
  1520 qed (intro tendsto_of_real, insert ln_Gamma_real_LIMSEQ[of x], simp add: ln_Gamma_def)
  1521 
  1522 lemma Gamma_real_pos_exp: "x > (0 :: real) \<Longrightarrow> Gamma x = exp (ln_Gamma x)"
  1523   by (auto simp: Gamma_real_altdef2 Gamma_complex_altdef of_real_in_nonpos_Ints_iff
  1524                  ln_Gamma_complex_of_real exp_of_real)
  1525 
  1526 lemma ln_Gamma_real_pos: "x > 0 \<Longrightarrow> ln_Gamma x = ln (Gamma x :: real)"
  1527   unfolding Gamma_real_pos_exp by simp
  1528 
  1529 lemma Gamma_real_pos: "x > (0::real) \<Longrightarrow> Gamma x > 0"
  1530   by (simp add: Gamma_real_pos_exp)
  1531 
  1532 lemma has_field_derivative_ln_Gamma_real [derivative_intros]:
  1533   assumes x: "x > (0::real)"
  1534   shows "(ln_Gamma has_field_derivative Digamma x) (at x)"
  1535 proof (subst DERIV_cong_ev[OF refl _ refl])
  1536   from assms show "((Re \<circ> ln_Gamma \<circ> complex_of_real) has_field_derivative Digamma x) (at x)"
  1537     by (auto intro!: derivative_eq_intros has_vector_derivative_real_complex
  1538              simp: Polygamma_of_real o_def)
  1539   from eventually_nhds_in_nhd[of x "{0<..}"] assms
  1540     show "eventually (\<lambda>y. ln_Gamma y = (Re \<circ> ln_Gamma \<circ> of_real) y) (nhds x)"
  1541     by (auto elim!: eventually_mono simp: ln_Gamma_complex_of_real interior_open)
  1542 qed
  1543 
  1544 declare has_field_derivative_ln_Gamma_real[THEN DERIV_chain2, derivative_intros]
  1545 
  1546 
  1547 lemma has_field_derivative_rGamma_real' [derivative_intros]:
  1548   "(rGamma has_field_derivative (if x \<in> \<int>\<^sub>\<le>\<^sub>0 then (-1)^(nat \<lfloor>-x\<rfloor>) * fact (nat \<lfloor>-x\<rfloor>) else
  1549         -rGamma x * Digamma x)) (at x within A)"
  1550   using has_field_derivative_rGamma[of x] by (force elim!: nonpos_Ints_cases')
  1551 
  1552 declare has_field_derivative_rGamma_real'[THEN DERIV_chain2, derivative_intros]
  1553 
  1554 lemma Polygamma_real_odd_pos:
  1555   assumes "(x::real) \<notin> \<int>\<^sub>\<le>\<^sub>0" "odd n"
  1556   shows   "Polygamma n x > 0"
  1557 proof -
  1558   from assms have "x \<noteq> 0" by auto
  1559   with assms show ?thesis
  1560     unfolding Polygamma_def using Polygamma_converges'[of x "Suc n"]
  1561     by (auto simp: zero_less_power_eq simp del: power_Suc
  1562              dest: plus_of_nat_eq_0_imp intro!: mult_pos_pos suminf_pos)
  1563 qed
  1564 
  1565 lemma Polygamma_real_even_neg:
  1566   assumes "(x::real) > 0" "n > 0" "even n"
  1567   shows   "Polygamma n x < 0"
  1568   using assms unfolding Polygamma_def using Polygamma_converges'[of x "Suc n"]
  1569   by (auto intro!: mult_pos_pos suminf_pos)
  1570 
  1571 lemma Polygamma_real_strict_mono:
  1572   assumes "x > 0" "x < (y::real)" "even n"
  1573   shows   "Polygamma n x < Polygamma n y"
  1574 proof -
  1575   have "\<exists>\<xi>. x < \<xi> \<and> \<xi> < y \<and> Polygamma n y - Polygamma n x = (y - x) * Polygamma (Suc n) \<xi>"
  1576     using assms by (intro MVT2 derivative_intros impI allI) (auto elim!: nonpos_Ints_cases)
  1577   then guess \<xi> by (elim exE conjE) note \<xi> = this
  1578   note \<xi>(3)
  1579   also from \<xi>(1,2) assms have "(y - x) * Polygamma (Suc n) \<xi> > 0"
  1580     by (intro mult_pos_pos Polygamma_real_odd_pos) (auto elim!: nonpos_Ints_cases)
  1581   finally show ?thesis by simp
  1582 qed
  1583 
  1584 lemma Polygamma_real_strict_antimono:
  1585   assumes "x > 0" "x < (y::real)" "odd n"
  1586   shows   "Polygamma n x > Polygamma n y"
  1587 proof -
  1588   have "\<exists>\<xi>. x < \<xi> \<and> \<xi> < y \<and> Polygamma n y - Polygamma n x = (y - x) * Polygamma (Suc n) \<xi>"
  1589     using assms by (intro MVT2 derivative_intros impI allI) (auto elim!: nonpos_Ints_cases)
  1590   then guess \<xi> by (elim exE conjE) note \<xi> = this
  1591   note \<xi>(3)
  1592   also from \<xi>(1,2) assms have "(y - x) * Polygamma (Suc n) \<xi> < 0"
  1593     by (intro mult_pos_neg Polygamma_real_even_neg) simp_all
  1594   finally show ?thesis by simp
  1595 qed
  1596 
  1597 lemma Polygamma_real_mono:
  1598   assumes "x > 0" "x \<le> (y::real)" "even n"
  1599   shows   "Polygamma n x \<le> Polygamma n y"
  1600   using Polygamma_real_strict_mono[OF assms(1) _ assms(3), of y] assms(2)
  1601   by (cases "x = y") simp_all
  1602 
  1603 lemma Digamma_real_ge_three_halves_pos:
  1604   assumes "x \<ge> 3/2"
  1605   shows   "Digamma (x :: real) > 0"
  1606 proof -
  1607   have "0 < Digamma (3/2 :: real)" by (fact Digamma_real_three_halves_pos)
  1608   also from assms have "\<dots> \<le> Digamma x" by (intro Polygamma_real_mono) simp_all
  1609   finally show ?thesis .
  1610 qed
  1611 
  1612 lemma ln_Gamma_real_strict_mono:
  1613   assumes "x \<ge> 3/2" "x < y"
  1614   shows   "ln_Gamma (x :: real) < ln_Gamma y"
  1615 proof -
  1616   have "\<exists>\<xi>. x < \<xi> \<and> \<xi> < y \<and> ln_Gamma y - ln_Gamma x = (y - x) * Digamma \<xi>"
  1617     using assms by (intro MVT2 derivative_intros impI allI) (auto elim!: nonpos_Ints_cases)
  1618   then guess \<xi> by (elim exE conjE) note \<xi> = this
  1619   note \<xi>(3)
  1620   also from \<xi>(1,2) assms have "(y - x) * Digamma \<xi> > 0"
  1621     by (intro mult_pos_pos Digamma_real_ge_three_halves_pos) simp_all
  1622   finally show ?thesis by simp
  1623 qed
  1624 
  1625 lemma Gamma_real_strict_mono:
  1626   assumes "x \<ge> 3/2" "x < y"
  1627   shows   "Gamma (x :: real) < Gamma y"
  1628 proof -
  1629   from Gamma_real_pos_exp[of x] assms have "Gamma x = exp (ln_Gamma x)" by simp
  1630   also have "\<dots> < exp (ln_Gamma y)" by (intro exp_less_mono ln_Gamma_real_strict_mono assms)
  1631   also from Gamma_real_pos_exp[of y] assms have "\<dots> = Gamma y" by simp
  1632   finally show ?thesis .
  1633 qed
  1634 
  1635 lemma log_convex_Gamma_real: "convex_on {0<..} (ln \<circ> Gamma :: real \<Rightarrow> real)"
  1636   by (rule convex_on_realI[of _ _ Digamma])
  1637      (auto intro!: derivative_eq_intros Polygamma_real_mono Gamma_real_pos
  1638            simp: o_def Gamma_eq_zero_iff elim!: nonpos_Ints_cases')
  1639 
  1640 
  1641 subsection \<open>Beta function\<close>
  1642 
  1643 definition Beta where "Beta a b = Gamma a * Gamma b / Gamma (a + b)"
  1644 
  1645 lemma Beta_altdef: "Beta a b = Gamma a * Gamma b * rGamma (a + b)"
  1646   by (simp add: inverse_eq_divide Beta_def Gamma_def)
  1647 
  1648 lemma Beta_commute: "Beta a b = Beta b a"
  1649   unfolding Beta_def by (simp add: ac_simps)
  1650 
  1651 lemma has_field_derivative_Beta1 [derivative_intros]:
  1652   assumes "x \<notin> \<int>\<^sub>\<le>\<^sub>0" "x + y \<notin> \<int>\<^sub>\<le>\<^sub>0"
  1653   shows   "((\<lambda>x. Beta x y) has_field_derivative (Beta x y * (Digamma x - Digamma (x + y))))
  1654                (at x within A)" unfolding Beta_altdef
  1655   by (rule DERIV_cong, (rule derivative_intros assms)+) (simp add: algebra_simps)
  1656 
  1657 lemma has_field_derivative_Beta2 [derivative_intros]:
  1658   assumes "y \<notin> \<int>\<^sub>\<le>\<^sub>0" "x + y \<notin> \<int>\<^sub>\<le>\<^sub>0"
  1659   shows   "((\<lambda>y. Beta x y) has_field_derivative (Beta x y * (Digamma y - Digamma (x + y))))
  1660                (at y within A)"
  1661   using has_field_derivative_Beta1[of y x A] assms by (simp add: Beta_commute add_ac)
  1662 
  1663 lemma Beta_plus1_plus1:
  1664   assumes "x \<notin> \<int>\<^sub>\<le>\<^sub>0" "y \<notin> \<int>\<^sub>\<le>\<^sub>0"
  1665   shows   "Beta (x + 1) y + Beta x (y + 1) = Beta x y"
  1666 proof -
  1667   have "Beta (x + 1) y + Beta x (y + 1) =
  1668             (Gamma (x + 1) * Gamma y + Gamma x * Gamma (y + 1)) * rGamma ((x + y) + 1)"
  1669     by (simp add: Beta_altdef add_divide_distrib algebra_simps)
  1670   also have "\<dots> = (Gamma x * Gamma y) * ((x + y) * rGamma ((x + y) + 1))"
  1671     by (subst assms[THEN Gamma_plus1])+ (simp add: algebra_simps)
  1672   also from assms have "\<dots> = Beta x y" unfolding Beta_altdef by (subst rGamma_plus1) simp
  1673   finally show ?thesis .
  1674 qed
  1675 
  1676 lemma Beta_plus1_left:
  1677   assumes "x \<notin> \<int>\<^sub>\<le>\<^sub>0" "y \<notin> \<int>\<^sub>\<le>\<^sub>0"
  1678   shows   "(x + y) * Beta (x + 1) y = x * Beta x y"
  1679 proof -
  1680   have "(x + y) * Beta (x + 1) y = Gamma (x + 1) * Gamma y * ((x + y) * rGamma ((x + y) + 1))"
  1681     unfolding Beta_altdef by (simp only: ac_simps)
  1682   also have "\<dots> = x * Beta x y" unfolding Beta_altdef
  1683      by (subst assms[THEN Gamma_plus1] rGamma_plus1)+ (simp only: ac_simps)
  1684   finally show ?thesis .
  1685 qed
  1686 
  1687 lemma Beta_plus1_right:
  1688   assumes "x \<notin> \<int>\<^sub>\<le>\<^sub>0" "y \<notin> \<int>\<^sub>\<le>\<^sub>0"
  1689   shows   "(x + y) * Beta x (y + 1) = y * Beta x y"
  1690   using Beta_plus1_left[of y x] assms by (simp add: Beta_commute add.commute)
  1691 
  1692 lemma Gamma_Gamma_Beta:
  1693   assumes "x \<notin> \<int>\<^sub>\<le>\<^sub>0" "y \<notin> \<int>\<^sub>\<le>\<^sub>0" "x + y \<notin> \<int>\<^sub>\<le>\<^sub>0"
  1694   shows   "Gamma x * Gamma y = Beta x y * Gamma (x + y)"
  1695   unfolding Beta_altdef using assms Gamma_eq_zero_iff[of "x+y"]
  1696   by (simp add: rGamma_inverse_Gamma)
  1697 
  1698 
  1699 
  1700 subsection \<open>Legendre duplication theorem\<close>
  1701 
  1702 context
  1703 begin
  1704 
  1705 private lemma Gamma_legendre_duplication_aux:
  1706   fixes z :: "'a :: Gamma"
  1707   assumes "z \<notin> \<int>\<^sub>\<le>\<^sub>0" "z + 1/2 \<notin> \<int>\<^sub>\<le>\<^sub>0"
  1708   shows "Gamma z * Gamma (z + 1/2) = exp ((1 - 2*z) * of_real (ln 2)) * Gamma (1/2) * Gamma (2*z)"
  1709 proof -
  1710   let ?powr = "\<lambda>b a. exp (a * of_real (ln (of_nat b)))"
  1711   let ?h = "\<lambda>n. (fact (n-1))\<^sup>2 / fact (2*n-1) * of_nat (2^(2*n)) *
  1712                 exp (1/2 * of_real (ln (real_of_nat n)))"
  1713   {
  1714     fix z :: 'a assume z: "z \<notin> \<int>\<^sub>\<le>\<^sub>0" "z + 1/2 \<notin> \<int>\<^sub>\<le>\<^sub>0"
  1715     let ?g = "\<lambda>n. ?powr 2 (2*z) * Gamma_series' z n * Gamma_series' (z + 1/2) n /
  1716                       Gamma_series' (2*z) (2*n)"
  1717     have "eventually (\<lambda>n. ?g n = ?h n) sequentially" using eventually_gt_at_top
  1718     proof eventually_elim
  1719       fix n :: nat assume n: "n > 0"
  1720       let ?f = "fact (n - 1) :: 'a" and ?f' = "fact (2*n - 1) :: 'a"
  1721       have A: "exp t * exp t = exp (2*t :: 'a)" for t by (subst exp_add [symmetric]) simp
  1722       have A: "Gamma_series' z n * Gamma_series' (z + 1/2) n = ?f^2 * ?powr n (2*z + 1/2) /
  1723                 (pochhammer z n * pochhammer (z + 1/2) n)"
  1724         by (simp add: Gamma_series'_def exp_add ring_distribs power2_eq_square A mult_ac)
  1725       have B: "Gamma_series' (2*z) (2*n) =
  1726                        ?f' * ?powr 2 (2*z) * ?powr n (2*z) /
  1727                        (of_nat (2^(2*n)) * pochhammer z n * pochhammer (z+1/2) n)" using n
  1728         by (simp add: Gamma_series'_def ln_mult exp_add ring_distribs pochhammer_double)
  1729       from z have "pochhammer z n \<noteq> 0" by (auto dest: pochhammer_eq_0_imp_nonpos_Int)
  1730       moreover from z have "pochhammer (z + 1/2) n \<noteq> 0" by (auto dest: pochhammer_eq_0_imp_nonpos_Int)
  1731       ultimately have "?powr 2 (2*z) * (Gamma_series' z n * Gamma_series' (z + 1/2) n) / Gamma_series' (2*z) (2*n) =
  1732          ?f^2 / ?f' * of_nat (2^(2*n)) * (?powr n ((4*z + 1)/2) * ?powr n (-2*z))"
  1733         using n unfolding A B by (simp add: divide_simps exp_minus)
  1734       also have "?powr n ((4*z + 1)/2) * ?powr n (-2*z) = ?powr n (1/2)"
  1735         by (simp add: algebra_simps exp_add[symmetric] add_divide_distrib)
  1736       finally show "?g n = ?h n" by (simp only: mult_ac)
  1737     qed
  1738 
  1739     moreover from z double_in_nonpos_Ints_imp[of z] have "2 * z \<notin> \<int>\<^sub>\<le>\<^sub>0" by auto
  1740     hence "?g \<longlonglongrightarrow> ?powr 2 (2*z) * Gamma z * Gamma (z+1/2) / Gamma (2*z)"
  1741       using lim_subseq[of "op * 2", OF _ Gamma_series'_LIMSEQ, of "2*z"]
  1742       by (intro tendsto_intros Gamma_series'_LIMSEQ)
  1743          (simp_all add: o_def subseq_def Gamma_eq_zero_iff)
  1744     ultimately have "?h \<longlonglongrightarrow> ?powr 2 (2*z) * Gamma z * Gamma (z+1/2) / Gamma (2*z)"
  1745       by (rule Lim_transform_eventually)
  1746   } note lim = this
  1747 
  1748   from assms double_in_nonpos_Ints_imp[of z] have z': "2 * z \<notin> \<int>\<^sub>\<le>\<^sub>0" by auto
  1749   from fraction_not_in_ints[of 2 1] have "(1/2 :: 'a) \<notin> \<int>\<^sub>\<le>\<^sub>0"
  1750     by (intro not_in_Ints_imp_not_in_nonpos_Ints) simp_all
  1751   with lim[of "1/2 :: 'a"] have "?h \<longlonglongrightarrow> 2 * Gamma (1 / 2 :: 'a)" by (simp add: exp_of_real)
  1752   from LIMSEQ_unique[OF this lim[OF assms]] z' show ?thesis
  1753     by (simp add: divide_simps Gamma_eq_zero_iff ring_distribs exp_diff exp_of_real ac_simps)
  1754 qed
  1755 
  1756 (* TODO: perhaps this is unnecessary once we have the fact that a holomorphic function is
  1757    infinitely differentiable *)
  1758 private lemma Gamma_reflection_aux:
  1759   defines "h \<equiv> \<lambda>z::complex. if z \<in> \<int> then 0 else
  1760                  (of_real pi * cot (of_real pi*z) + Digamma z - Digamma (1 - z))"
  1761   defines "a \<equiv> complex_of_real pi"
  1762   obtains h' where "continuous_on UNIV h'" "\<And>z. (h has_field_derivative (h' z)) (at z)"
  1763 proof -
  1764   def f \<equiv> "\<lambda>n. a * of_real (cos_coeff (n+1) - sin_coeff (n+2))"
  1765   def F \<equiv> "\<lambda>z. if z = 0 then 0 else (cos (a*z) - sin (a*z)/(a*z)) / z"
  1766   def g \<equiv> "\<lambda>n. complex_of_real (sin_coeff (n+1))"
  1767   def G \<equiv> "\<lambda>z. if z = 0 then 1 else sin (a*z)/(a*z)"
  1768   have a_nz: "a \<noteq> 0" unfolding a_def by simp
  1769 
  1770   have "(\<lambda>n. f n * (a*z)^n) sums (F z) \<and> (\<lambda>n. g n * (a*z)^n) sums (G z)"
  1771     if "abs (Re z) < 1" for z
  1772   proof (cases "z = 0"; rule conjI)
  1773     assume "z \<noteq> 0"
  1774     note z = this that
  1775 
  1776     from z have sin_nz: "sin (a*z) \<noteq> 0" unfolding a_def by (auto simp: sin_eq_0)
  1777     have "(\<lambda>n. of_real (sin_coeff n) * (a*z)^n) sums (sin (a*z))" using sin_converges[of "a*z"]
  1778       by (simp add: scaleR_conv_of_real)
  1779     from sums_split_initial_segment[OF this, of 1]
  1780       have "(\<lambda>n. (a*z) * of_real (sin_coeff (n+1)) * (a*z)^n) sums (sin (a*z))" by (simp add: mult_ac)
  1781     from sums_mult[OF this, of "inverse (a*z)"] z a_nz
  1782       have A: "(\<lambda>n. g n * (a*z)^n) sums (sin (a*z)/(a*z))"
  1783       by (simp add: field_simps g_def)
  1784     with z show "(\<lambda>n. g n * (a*z)^n) sums (G z)" by (simp add: G_def)
  1785     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)
  1786 
  1787     have [simp]: "sin_coeff (Suc 0) = 1" by (simp add: sin_coeff_def)
  1788     from sums_split_initial_segment[OF sums_diff[OF cos_converges[of "a*z"] A], of 1]
  1789     have "(\<lambda>n. z * f n * (a*z)^n) sums (cos (a*z) - sin (a*z) / (a*z))"
  1790       by (simp add: mult_ac scaleR_conv_of_real ring_distribs f_def g_def)
  1791     from sums_mult[OF this, of "inverse z"] z assms
  1792       show "(\<lambda>n. f n * (a*z)^n) sums (F z)" by (simp add: divide_simps mult_ac f_def F_def)
  1793   next
  1794     assume z: "z = 0"
  1795     have "(\<lambda>n. f n * (a * z) ^ n) sums f 0" using powser_sums_zero[of f] z by simp
  1796     with z show "(\<lambda>n. f n * (a * z) ^ n) sums (F z)"
  1797       by (simp add: f_def F_def sin_coeff_def cos_coeff_def)
  1798     have "(\<lambda>n. g n * (a * z) ^ n) sums g 0" using powser_sums_zero[of g] z by simp
  1799     with z show "(\<lambda>n. g n * (a * z) ^ n) sums (G z)"
  1800       by (simp add: g_def G_def sin_coeff_def cos_coeff_def)
  1801   qed
  1802   note sums = conjunct1[OF this] conjunct2[OF this]
  1803 
  1804   def h2 \<equiv> "\<lambda>z. (\<Sum>n. f n * (a*z)^n) / (\<Sum>n. g n * (a*z)^n) +
  1805             Digamma (1 + z) - Digamma (1 - z)"
  1806   def POWSER \<equiv> "\<lambda>f z. (\<Sum>n. f n * (z^n :: complex))"
  1807   def POWSER' \<equiv> "\<lambda>f z. (\<Sum>n. diffs f n * (z^n :: complex))"
  1808 
  1809   def h2' \<equiv> "\<lambda>z. a * (POWSER g (a*z) * POWSER' f (a*z) - POWSER f (a*z) * POWSER' g (a*z)) /
  1810                      (POWSER g (a*z))^2 + Polygamma 1 (1 + z) + Polygamma 1 (1 - z)"
  1811 
  1812   have h_eq: "h t = h2 t" if "abs (Re t) < 1" for t
  1813   proof -
  1814     from that have t: "t \<in> \<int> \<longleftrightarrow> t = 0" by (auto elim!: Ints_cases simp: dist_0_norm)
  1815     hence "h t = a*cot (a*t) - 1/t + Digamma (1 + t) - Digamma (1 - t)"
  1816       unfolding h_def using Digamma_plus1[of t] by (force simp: field_simps a_def)
  1817     also have "a*cot (a*t) - 1/t = (F t) / (G t)"
  1818       using t by (auto simp add: divide_simps sin_eq_0 cot_def a_def F_def G_def)
  1819     also have "\<dots> = (\<Sum>n. f n * (a*t)^n) / (\<Sum>n. g n * (a*t)^n)"
  1820       using sums[of t] that by (simp add: sums_iff dist_0_norm)
  1821     finally show "h t = h2 t" by (simp only: h2_def)
  1822   qed
  1823 
  1824   let ?A = "{z. abs (Re z) < 1}"
  1825   have "open ({z. Re z < 1} \<inter> {z. Re z > -1})"
  1826     using open_halfspace_Re_gt open_halfspace_Re_lt by auto
  1827   also have "({z. Re z < 1} \<inter> {z. Re z > -1}) = {z. abs (Re z) < 1}" by auto
  1828   finally have open_A: "open ?A" .
  1829   hence [simp]: "interior ?A = ?A" by (simp add: interior_open)
  1830 
  1831   have summable_f: "summable (\<lambda>n. f n * z^n)" for z
  1832     by (rule powser_inside, rule sums_summable, rule sums[of "\<i> * of_real (norm z + 1) / a"])
  1833        (simp_all add: norm_mult a_def del: of_real_add)
  1834   have summable_g: "summable (\<lambda>n. g n * z^n)" for z
  1835     by (rule powser_inside, rule sums_summable, rule sums[of "\<i> * of_real (norm z + 1) / a"])
  1836        (simp_all add: norm_mult a_def del: of_real_add)
  1837   have summable_fg': "summable (\<lambda>n. diffs f n * z^n)" "summable (\<lambda>n. diffs g n * z^n)" for z
  1838     by (intro termdiff_converges_all summable_f summable_g)+
  1839   have "(POWSER f has_field_derivative (POWSER' f z)) (at z)"
  1840                "(POWSER g has_field_derivative (POWSER' g z)) (at z)" for z
  1841     unfolding POWSER_def POWSER'_def
  1842     by (intro termdiffs_strong_converges_everywhere summable_f summable_g)+
  1843   note derivs = this[THEN DERIV_chain2[OF _ DERIV_cmult[OF DERIV_ident]], unfolded POWSER_def]
  1844   have "isCont (POWSER f) z" "isCont (POWSER g) z" "isCont (POWSER' f) z" "isCont (POWSER' g) z"
  1845     for z unfolding POWSER_def POWSER'_def
  1846     by (intro isCont_powser_converges_everywhere summable_f summable_g summable_fg')+
  1847   note cont = this[THEN isCont_o2[rotated], unfolded POWSER_def POWSER'_def]
  1848 
  1849   {
  1850     fix z :: complex assume z: "abs (Re z) < 1"
  1851     def d \<equiv> "\<i> * of_real (norm z + 1)"
  1852     have d: "abs (Re d) < 1" "norm z < norm d" by (simp_all add: d_def norm_mult del: of_real_add)
  1853     have "eventually (\<lambda>z. h z = h2 z) (nhds z)"
  1854       using eventually_nhds_in_nhd[of z ?A] using h_eq z
  1855       by (auto elim!: eventually_mono simp: dist_0_norm)
  1856 
  1857     moreover from sums(2)[OF z] z have nz: "(\<Sum>n. g n * (a * z) ^ n) \<noteq> 0"
  1858       unfolding G_def by (auto simp: sums_iff sin_eq_0 a_def)
  1859     have A: "z \<in> \<int> \<longleftrightarrow> z = 0" using z by (auto elim!: Ints_cases)
  1860     have no_int: "1 + z \<in> \<int> \<longleftrightarrow> z = 0" using z Ints_diff[of "1+z" 1] A
  1861       by (auto elim!: nonpos_Ints_cases)
  1862     have no_int': "1 - z \<in> \<int> \<longleftrightarrow> z = 0" using z Ints_diff[of 1 "1-z"] A
  1863       by (auto elim!: nonpos_Ints_cases)
  1864     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
  1865     have "(h2 has_field_derivative h2' z) (at z)" unfolding h2_def
  1866       by (rule DERIV_cong, (rule derivative_intros refl derivs[unfolded POWSER_def] nz no_int)+)
  1867          (auto simp: h2'_def POWSER_def field_simps power2_eq_square)
  1868     ultimately have deriv: "(h has_field_derivative h2' z) (at z)"
  1869       by (subst DERIV_cong_ev[OF refl _ refl])
  1870 
  1871     from sums(2)[OF z] z have "(\<Sum>n. g n * (a * z) ^ n) \<noteq> 0"
  1872       unfolding G_def by (auto simp: sums_iff a_def sin_eq_0)
  1873     hence "isCont h2' z" using no_int unfolding h2'_def[abs_def] POWSER_def POWSER'_def
  1874       by (intro continuous_intros cont
  1875             continuous_on_compose2[OF _ continuous_on_Polygamma[of "{z. Re z > 0}"]]) auto
  1876     note deriv and this
  1877   } note A = this
  1878 
  1879   interpret h: periodic_fun_simple' h
  1880   proof
  1881     fix z :: complex
  1882     show "h (z + 1) = h z"
  1883     proof (cases "z \<in> \<int>")
  1884       assume z: "z \<notin> \<int>"
  1885       hence A: "z + 1 \<notin> \<int>" "z \<noteq> 0" using Ints_diff[of "z+1" 1] by auto
  1886       hence "Digamma (z + 1) - Digamma (-z) = Digamma z - Digamma (-z + 1)"
  1887         by (subst (1 2) Digamma_plus1) simp_all
  1888       with A z show "h (z + 1) = h z"
  1889         by (simp add: h_def sin_plus_pi cos_plus_pi ring_distribs cot_def)
  1890     qed (simp add: h_def)
  1891   qed
  1892 
  1893   have h2'_eq: "h2' (z - 1) = h2' z" if z: "Re z > 0" "Re z < 1" for z
  1894   proof -
  1895     have "((\<lambda>z. h (z - 1)) has_field_derivative h2' (z - 1)) (at z)"
  1896       by (rule DERIV_cong, rule DERIV_chain'[OF _ A(1)])
  1897          (insert z, auto intro!: derivative_eq_intros)
  1898     hence "(h has_field_derivative h2' (z - 1)) (at z)" by (subst (asm) h.minus_1)
  1899     moreover from z have "(h has_field_derivative h2' z) (at z)" by (intro A) simp_all
  1900     ultimately show "h2' (z - 1) = h2' z" by (rule DERIV_unique)
  1901   qed
  1902 
  1903   def h2'' \<equiv> "\<lambda>z. h2' (z - of_int \<lfloor>Re z\<rfloor>)"
  1904   have deriv: "(h has_field_derivative h2'' z) (at z)" for z
  1905   proof -
  1906     fix z :: complex
  1907     have B: "\<bar>Re z - real_of_int \<lfloor>Re z\<rfloor>\<bar> < 1" by linarith
  1908     have "((\<lambda>t. h (t - of_int \<lfloor>Re z\<rfloor>)) has_field_derivative h2'' z) (at z)"
  1909       unfolding h2''_def by (rule DERIV_cong, rule DERIV_chain'[OF _ A(1)])
  1910                             (insert B, auto intro!: derivative_intros)
  1911     thus "(h has_field_derivative h2'' z) (at z)" by (simp add: h.minus_of_int)
  1912   qed
  1913 
  1914   have cont: "continuous_on UNIV h2''"
  1915   proof (intro continuous_at_imp_continuous_on ballI)
  1916     fix z :: complex
  1917     def r \<equiv> "\<lfloor>Re z\<rfloor>"
  1918     def A \<equiv> "{t. of_int r - 1 < Re t \<and> Re t < of_int r + 1}"
  1919     have "continuous_on A (\<lambda>t. h2' (t - of_int r))" unfolding A_def
  1920       by (intro continuous_at_imp_continuous_on isCont_o2[OF _ A(2)] ballI continuous_intros)
  1921          (simp_all add: abs_real_def)
  1922     moreover have "h2'' t = h2' (t - of_int r)" if t: "t \<in> A" for t
  1923     proof (cases "Re t \<ge> of_int r")
  1924       case True
  1925       from t have "of_int r - 1 < Re t" "Re t < of_int r + 1" by (simp_all add: A_def)
  1926       with True have "\<lfloor>Re t\<rfloor> = \<lfloor>Re z\<rfloor>" unfolding r_def by linarith
  1927       thus ?thesis by (auto simp: r_def h2''_def)
  1928     next
  1929       case False
  1930       from t have t: "of_int r - 1 < Re t" "Re t < of_int r + 1" by (simp_all add: A_def)
  1931       with False have t': "\<lfloor>Re t\<rfloor> = \<lfloor>Re z\<rfloor> - 1" unfolding r_def by linarith
  1932       moreover from t False have "h2' (t - of_int r + 1 - 1) = h2' (t - of_int r + 1)"
  1933         by (intro h2'_eq) simp_all
  1934       ultimately show ?thesis by (auto simp: r_def h2''_def algebra_simps t')
  1935     qed
  1936     ultimately have "continuous_on A h2''" by (subst continuous_on_cong[OF refl])
  1937     moreover {
  1938       have "open ({t. of_int r - 1 < Re t} \<inter> {t. of_int r + 1 > Re t})"
  1939         by (intro open_Int open_halfspace_Re_gt open_halfspace_Re_lt)
  1940       also have "{t. of_int r - 1 < Re t} \<inter> {t. of_int r + 1 > Re t} = A"
  1941         unfolding A_def by blast
  1942       finally have "open A" .
  1943     }
  1944     ultimately have C: "isCont h2'' t" if "t \<in> A" for t using that
  1945       by (subst (asm) continuous_on_eq_continuous_at) auto
  1946     have "of_int r - 1 < Re z" "Re z  < of_int r + 1" unfolding r_def by linarith+
  1947     thus "isCont h2'' z" by (intro C) (simp_all add: A_def)
  1948   qed
  1949 
  1950   from that[OF cont deriv] show ?thesis .
  1951 qed
  1952 
  1953 lemma Gamma_reflection_complex:
  1954   fixes z :: complex
  1955   shows "Gamma z * Gamma (1 - z) = of_real pi / sin (of_real pi * z)"
  1956 proof -
  1957   let ?g = "\<lambda>z::complex. Gamma z * Gamma (1 - z) * sin (of_real pi * z)"
  1958   def g \<equiv> "\<lambda>z::complex. if z \<in> \<int> then of_real pi else ?g z"
  1959   let ?h = "\<lambda>z::complex. (of_real pi * cot (of_real pi*z) + Digamma z - Digamma (1 - z))"
  1960   def h \<equiv> "\<lambda>z::complex. if z \<in> \<int> then 0 else ?h z"
  1961 
  1962   \<comment> \<open>@{term g} is periodic with period 1.\<close>
  1963   interpret g: periodic_fun_simple' g
  1964   proof
  1965     fix z :: complex
  1966     show "g (z + 1) = g z"
  1967     proof (cases "z \<in> \<int>")
  1968       case False
  1969       hence "z * g z = z * Beta z (- z + 1) * sin (of_real pi * z)" by (simp add: g_def Beta_def)
  1970       also have "z * Beta z (- z + 1) = (z + 1 + -z) * Beta (z + 1) (- z + 1)"
  1971         using False Ints_diff[of 1 "1 - z"] nonpos_Ints_subset_Ints
  1972         by (subst Beta_plus1_left [symmetric]) auto
  1973       also have "\<dots> * sin (of_real pi * z) = z * (Beta (z + 1) (-z) * sin (of_real pi * (z + 1)))"
  1974         using False Ints_diff[of "z+1" 1] Ints_minus[of "-z"] nonpos_Ints_subset_Ints
  1975         by (subst Beta_plus1_right) (auto simp: ring_distribs sin_plus_pi)
  1976       also from False have "Beta (z + 1) (-z) * sin (of_real pi * (z + 1)) = g (z + 1)"
  1977         using Ints_diff[of "z+1" 1] by (auto simp: g_def Beta_def)
  1978       finally show "g (z + 1) = g z" using False by (subst (asm) mult_left_cancel) auto
  1979     qed (simp add: g_def)
  1980   qed
  1981 
  1982   \<comment> \<open>@{term g} is entire.\<close>
  1983   have g_g': "(g has_field_derivative (h z * g z)) (at z)" for z :: complex
  1984   proof (cases "z \<in> \<int>")
  1985     let ?h' = "\<lambda>z. Beta z (1 - z) * ((Digamma z - Digamma (1 - z)) * sin (z * of_real pi) +
  1986                      of_real pi * cos (z * of_real pi))"
  1987     case False
  1988     from False have "eventually (\<lambda>t. t \<in> UNIV - \<int>) (nhds z)"
  1989       by (intro eventually_nhds_in_open) (auto simp: open_Diff)
  1990     hence "eventually (\<lambda>t. g t = ?g t) (nhds z)" by eventually_elim (simp add: g_def)
  1991     moreover {
  1992       from False Ints_diff[of 1 "1-z"] have "1 - z \<notin> \<int>" by auto
  1993       hence "(?g has_field_derivative ?h' z) (at z)" using nonpos_Ints_subset_Ints
  1994         by (auto intro!: derivative_eq_intros simp: algebra_simps Beta_def)
  1995       also from False have "sin (of_real pi * z) \<noteq> 0" by (subst sin_eq_0) auto
  1996       hence "?h' z = h z * g z"
  1997         using False unfolding g_def h_def cot_def by (simp add: field_simps Beta_def)
  1998       finally have "(?g has_field_derivative (h z * g z)) (at z)" .
  1999     }
  2000     ultimately show ?thesis by (subst DERIV_cong_ev[OF refl _ refl])
  2001   next
  2002     case True
  2003     then obtain n where z: "z = of_int n" by (auto elim!: Ints_cases)
  2004     let ?t = "(\<lambda>z::complex. if z = 0 then 1 else sin z / z) \<circ> (\<lambda>z. of_real pi * z)"
  2005     have deriv_0: "(g has_field_derivative 0) (at 0)"
  2006     proof (subst DERIV_cong_ev[OF refl _ refl])
  2007       show "eventually (\<lambda>z. g z = of_real pi * Gamma (1 + z) * Gamma (1 - z) * ?t z) (nhds 0)"
  2008         using eventually_nhds_ball[OF zero_less_one, of "0::complex"]
  2009       proof eventually_elim
  2010         fix z :: complex assume z: "z \<in> ball 0 1"
  2011         show "g z = of_real pi * Gamma (1 + z) * Gamma (1 - z) * ?t z"
  2012         proof (cases "z = 0")
  2013           assume z': "z \<noteq> 0"
  2014           with z have z'': "z \<notin> \<int>\<^sub>\<le>\<^sub>0" "z \<notin> \<int>" by (auto elim!: Ints_cases simp: dist_0_norm)
  2015           from Gamma_plus1[OF this(1)] have "Gamma z = Gamma (z + 1) / z" by simp
  2016           with z'' z' show ?thesis by (simp add: g_def ac_simps)
  2017         qed (simp add: g_def)
  2018       qed
  2019       have "(?t has_field_derivative (0 * of_real pi)) (at 0)"
  2020         using has_field_derivative_sin_z_over_z[of "UNIV :: complex set"]
  2021         by (intro DERIV_chain) simp_all
  2022       thus "((\<lambda>z. of_real pi * Gamma (1 + z) * Gamma (1 - z) * ?t z) has_field_derivative 0) (at 0)"
  2023         by (auto intro!: derivative_eq_intros simp: o_def)
  2024     qed
  2025 
  2026     have "((g \<circ> (\<lambda>x. x - of_int n)) has_field_derivative 0 * 1) (at (of_int n))"
  2027       using deriv_0 by (intro DERIV_chain) (auto intro!: derivative_eq_intros)
  2028     also have "g \<circ> (\<lambda>x. x - of_int n) = g" by (intro ext) (simp add: g.minus_of_int)
  2029     finally show "(g has_field_derivative (h z * g z)) (at z)" by (simp add: z h_def)
  2030   qed
  2031 
  2032   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
  2033   proof (cases "z \<in> \<int>")
  2034     case True
  2035     with that have "z = 0 \<or> z = 1" by (force elim!: Ints_cases)
  2036     moreover have "g 0 * g (1/2) = Gamma (1/2)^2 * g 0"
  2037       using fraction_not_in_ints[where 'a = complex, of 2 1] by (simp add: g_def power2_eq_square)
  2038     moreover have "g (1/2) * g 1 = Gamma (1/2)^2 * g 1"
  2039         using fraction_not_in_ints[where 'a = complex, of 2 1]
  2040         by (simp add: g_def power2_eq_square Beta_def algebra_simps)
  2041     ultimately show ?thesis by force
  2042   next
  2043     case False
  2044     hence z: "z/2 \<notin> \<int>" "(z+1)/2 \<notin> \<int>" using Ints_diff[of "z+1" 1] by (auto elim!: Ints_cases)
  2045     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)
  2046     from z have "1-z/2 \<notin> \<int>" "1-((z+1)/2) \<notin> \<int>"
  2047       using Ints_diff[of 1 "1-z/2"] Ints_diff[of 1 "1-((z+1)/2)"] by auto
  2048     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)
  2049     from z have "g (z/2) * g ((z+1)/2) =
  2050       (Gamma (z/2) * Gamma ((z+1)/2)) * (Gamma (1-z/2) * Gamma (1-((z+1)/2))) *
  2051       (sin (of_real pi * z/2) * sin (of_real pi * (z+1)/2))"
  2052       by (simp add: g_def)
  2053     also from z' Gamma_legendre_duplication_aux[of "z/2"]
  2054       have "Gamma (z/2) * Gamma ((z+1)/2) = exp ((1-z) * of_real (ln 2)) * Gamma (1/2) * Gamma z"
  2055       by (simp add: add_divide_distrib)
  2056     also from z'' Gamma_legendre_duplication_aux[of "1-(z+1)/2"]
  2057       have "Gamma (1-z/2) * Gamma (1-(z+1)/2) =
  2058               Gamma (1-z) * Gamma (1/2) * exp (z * of_real (ln 2))"
  2059       by (simp add: add_divide_distrib ac_simps)
  2060     finally have "g (z/2) * g ((z+1)/2) = Gamma (1/2)^2 * (Gamma z * Gamma (1-z) *
  2061                     (2 * (sin (of_real pi*z/2) * sin (of_real pi*(z+1)/2))))"
  2062       by (simp add: add_ac power2_eq_square exp_add ring_distribs exp_diff exp_of_real)
  2063     also have "sin (of_real pi*(z+1)/2) = cos (of_real pi*z/2)"
  2064       using cos_sin_eq[of "- of_real pi * z/2", symmetric]
  2065       by (simp add: ring_distribs add_divide_distrib ac_simps)
  2066     also have "2 * (sin (of_real pi*z/2) * cos (of_real pi*z/2)) = sin (of_real pi * z)"
  2067       by (subst sin_times_cos) (simp add: field_simps)
  2068     also have "Gamma z * Gamma (1 - z) * sin (complex_of_real pi * z) = g z"
  2069       using \<open>z \<notin> \<int>\<close> by (simp add: g_def)
  2070     finally show ?thesis .
  2071   qed
  2072   have g_eq: "g (z/2) * g ((z+1)/2) = Gamma (1/2)^2 * g z" for z
  2073   proof -
  2074     def r \<equiv> "\<lfloor>Re z / 2\<rfloor>"
  2075     have "Gamma (1/2)^2 * g z = Gamma (1/2)^2 * g (z - of_int (2*r))" by (simp only: g.minus_of_int)
  2076     also have "of_int (2*r) = 2 * of_int r" by simp
  2077     also have "Re z - 2 * of_int r > -1" "Re z - 2 * of_int r < 2" unfolding r_def by linarith+
  2078     hence "Gamma (1/2)^2 * g (z - 2 * of_int r) =
  2079                    g ((z - 2 * of_int r)/2) * g ((z - 2 * of_int r + 1)/2)"
  2080       unfolding r_def by (intro g_eq[symmetric]) simp_all
  2081     also have "(z - 2 * of_int r) / 2 = z/2 - of_int r" by simp
  2082     also have "g \<dots> = g (z/2)" by (rule g.minus_of_int)
  2083     also have "(z - 2 * of_int r + 1) / 2 = (z + 1)/2 - of_int r" by simp
  2084     also have "g \<dots> = g ((z+1)/2)" by (rule g.minus_of_int)
  2085     finally show ?thesis ..
  2086   qed
  2087 
  2088   have g_nz [simp]: "g z \<noteq> 0" for z :: complex
  2089   unfolding g_def using Ints_diff[of 1 "1 - z"]
  2090     by (auto simp: Gamma_eq_zero_iff sin_eq_0 dest!: nonpos_Ints_Int)
  2091 
  2092   have h_eq: "h z = (h (z/2) + h ((z+1)/2)) / 2" for z
  2093   proof -
  2094     have "((\<lambda>t. g (t/2) * g ((t+1)/2)) has_field_derivative
  2095                        (g (z/2) * g ((z+1)/2)) * ((h (z/2) + h ((z+1)/2)) / 2)) (at z)"
  2096       by (auto intro!: derivative_eq_intros g_g'[THEN DERIV_chain2] simp: field_simps)
  2097     hence "((\<lambda>t. Gamma (1/2)^2 * g t) has_field_derivative
  2098               Gamma (1/2)^2 * g z * ((h (z/2) + h ((z+1)/2)) / 2)) (at z)"
  2099       by (subst (1 2) g_eq[symmetric]) simp
  2100     from DERIV_cmult[OF this, of "inverse ((Gamma (1/2))^2)"]
  2101       have "(g has_field_derivative (g z * ((h (z/2) + h ((z+1)/2))/2))) (at z)"
  2102       using fraction_not_in_ints[where 'a = complex, of 2 1]
  2103       by (simp add: divide_simps Gamma_eq_zero_iff not_in_Ints_imp_not_in_nonpos_Ints)
  2104     moreover have "(g has_field_derivative (g z * h z)) (at z)"
  2105       using g_g'[of z] by (simp add: ac_simps)
  2106     ultimately have "g z * h z = g z * ((h (z/2) + h ((z+1)/2))/2)"
  2107       by (intro DERIV_unique)
  2108     thus "h z = (h (z/2) + h ((z+1)/2)) / 2" by simp
  2109   qed
  2110 
  2111   obtain h' where h'_cont: "continuous_on UNIV h'" and
  2112                   h_h': "\<And>z. (h has_field_derivative h' z) (at z)"
  2113      unfolding h_def by (erule Gamma_reflection_aux)
  2114 
  2115   have h'_eq: "h' z = (h' (z/2) + h' ((z+1)/2)) / 4" for z
  2116   proof -
  2117     have "((\<lambda>t. (h (t/2) + h ((t+1)/2)) / 2) has_field_derivative
  2118                        ((h' (z/2) + h' ((z+1)/2)) / 4)) (at z)"
  2119       by (fastforce intro!: derivative_eq_intros h_h'[THEN DERIV_chain2])
  2120     hence "(h has_field_derivative ((h' (z/2) + h' ((z+1)/2))/4)) (at z)"
  2121       by (subst (asm) h_eq[symmetric])
  2122     from h_h' and this show "h' z = (h' (z/2) + h' ((z+1)/2)) / 4" by (rule DERIV_unique)
  2123   qed
  2124 
  2125   have h'_zero: "h' z = 0" for z
  2126   proof -
  2127     def m \<equiv> "max 1 \<bar>Re z\<bar>"
  2128     def B \<equiv> "{t. abs (Re t) \<le> m \<and> abs (Im t) \<le> abs (Im z)}"
  2129     have "closed ({t. Re t \<ge> -m} \<inter> {t. Re t \<le> m} \<inter>
  2130                   {t. Im t \<ge> -\<bar>Im z\<bar>} \<inter> {t. Im t \<le> \<bar>Im z\<bar>})"
  2131       (is "closed ?B") by (intro closed_Int closed_halfspace_Re_ge closed_halfspace_Re_le
  2132                                  closed_halfspace_Im_ge closed_halfspace_Im_le)
  2133     also have "?B = B" unfolding B_def by fastforce
  2134     finally have "closed B" .
  2135     moreover have "bounded B" unfolding bounded_iff
  2136     proof (intro ballI exI)
  2137       fix t assume t: "t \<in> B"
  2138       have "norm t \<le> \<bar>Re t\<bar> + \<bar>Im t\<bar>" by (rule cmod_le)
  2139       also from t have "\<bar>Re t\<bar> \<le> m" unfolding B_def by blast
  2140       also from t have "\<bar>Im t\<bar> \<le> \<bar>Im z\<bar>" unfolding B_def by blast
  2141       finally show "norm t \<le> m + \<bar>Im z\<bar>" by - simp
  2142     qed
  2143     ultimately have compact: "compact B" by (subst compact_eq_bounded_closed) blast
  2144 
  2145     def M \<equiv> "SUP z:B. norm (h' z)"
  2146     have "compact (h' ` B)"
  2147       by (intro compact_continuous_image continuous_on_subset[OF h'_cont] compact) blast+
  2148     hence bdd: "bdd_above ((\<lambda>z. norm (h' z)) ` B)"
  2149       using bdd_above_norm[of "h' ` B"] by (simp add: image_comp o_def compact_imp_bounded)
  2150     have "norm (h' z) \<le> M" unfolding M_def by (intro cSUP_upper bdd) (simp_all add: B_def m_def)
  2151     also have "M \<le> M/2"
  2152     proof (subst M_def, subst cSUP_le_iff)
  2153       have "z \<in> B" unfolding B_def m_def by simp
  2154       thus "B \<noteq> {}" by auto
  2155     next
  2156       show "\<forall>z\<in>B. norm (h' z) \<le> M/2"
  2157       proof
  2158         fix t :: complex assume t: "t \<in> B"
  2159         from h'_eq[of t] t have "h' t = (h' (t/2) + h' ((t+1)/2)) / 4" by (simp add: dist_0_norm)
  2160         also have "norm \<dots> = norm (h' (t/2) + h' ((t+1)/2)) / 4" by simp
  2161         also have "norm (h' (t/2) + h' ((t+1)/2)) \<le> norm (h' (t/2)) + norm (h' ((t+1)/2))"
  2162           by (rule norm_triangle_ineq)
  2163         also from t have "abs (Re ((t + 1)/2)) \<le> m" unfolding m_def B_def by auto
  2164         with t have "t/2 \<in> B" "(t+1)/2 \<in> B" unfolding B_def by auto
  2165         hence "norm (h' (t/2)) + norm (h' ((t+1)/2)) \<le> M + M" unfolding M_def
  2166           by (intro add_mono cSUP_upper bdd) (auto simp: B_def)
  2167         also have "(M + M) / 4 = M / 2" by simp
  2168         finally show "norm (h' t) \<le> M/2" by - simp_all
  2169       qed
  2170     qed (insert bdd, auto simp: cball_eq_empty)
  2171     hence "M \<le> 0" by simp
  2172     finally show "h' z = 0" by simp
  2173   qed
  2174   have h_h'_2: "(h has_field_derivative 0) (at z)" for z
  2175     using h_h'[of z] h'_zero[of z] by simp
  2176 
  2177   have g_real: "g z \<in> \<real>" if "z \<in> \<real>" for z
  2178     unfolding g_def using that by (auto intro!: Reals_mult Gamma_complex_real)
  2179   have h_real: "h z \<in> \<real>" if "z \<in> \<real>" for z
  2180     unfolding h_def using that by (auto intro!: Reals_mult Reals_add Reals_diff Polygamma_Real)
  2181   have g_nz: "g z \<noteq> 0" for z unfolding g_def using Ints_diff[of 1 "1-z"]
  2182     by (auto simp: Gamma_eq_zero_iff sin_eq_0)
  2183 
  2184   from h'_zero h_h'_2 have "\<exists>c. \<forall>z\<in>UNIV. h z = c"
  2185     by (intro has_field_derivative_zero_constant) (simp_all add: dist_0_norm)
  2186   then obtain c where c: "\<And>z. h z = c" by auto
  2187   have "\<exists>u. u \<in> closed_segment 0 1 \<and> Re (g 1) - Re (g 0) = Re (h u * g u * (1 - 0))"
  2188     by (intro complex_mvt_line g_g')
  2189     find_theorems name:deriv Reals
  2190   then guess u by (elim exE conjE) note u = this
  2191   from u(1) have u': "u \<in> \<real>" unfolding closed_segment_def
  2192     by (auto simp: scaleR_conv_of_real)
  2193   from u' g_real[of u] g_nz[of u] have "Re (g u) \<noteq> 0" by (auto elim!: Reals_cases)
  2194   with u(2) c[of u] g_real[of u] g_nz[of u] u'
  2195     have "Re c = 0" by (simp add: complex_is_Real_iff g.of_1)
  2196   with h_real[of 0] c[of 0] have "c = 0" by (auto elim!: Reals_cases)
  2197   with c have A: "h z * g z = 0" for z by simp
  2198   hence "(g has_field_derivative 0) (at z)" for z using g_g'[of z] by simp
  2199   hence "\<exists>c'. \<forall>z\<in>UNIV. g z = c'" by (intro has_field_derivative_zero_constant) simp_all
  2200   then obtain c' where c: "\<And>z. g z = c'" by (force simp: dist_0_norm)
  2201   moreover from this[of 0] have "c' = pi" unfolding g_def by simp
  2202   ultimately have "g z = pi" by simp
  2203 
  2204   show ?thesis
  2205   proof (cases "z \<in> \<int>")
  2206     case False
  2207     with \<open>g z = pi\<close> show ?thesis by (auto simp: g_def divide_simps)
  2208   next
  2209     case True
  2210     then obtain n where n: "z = of_int n" by (elim Ints_cases)
  2211     with sin_eq_0[of "of_real pi * z"] have "sin (of_real pi * z) = 0" by force
  2212     moreover have "of_int (1 - n) \<in> \<int>\<^sub>\<le>\<^sub>0" if "n > 0" using that by (intro nonpos_Ints_of_int) simp
  2213     ultimately show ?thesis using n
  2214       by (cases "n \<le> 0") (auto simp: Gamma_eq_zero_iff nonpos_Ints_of_int)
  2215   qed
  2216 qed
  2217 
  2218 lemma rGamma_reflection_complex:
  2219   "rGamma z * rGamma (1 - z :: complex) = sin (of_real pi * z) / of_real pi"
  2220   using Gamma_reflection_complex[of z]
  2221     by (simp add: Gamma_def divide_simps split: split_if_asm)
  2222 
  2223 lemma rGamma_reflection_complex':
  2224   "rGamma z * rGamma (- z :: complex) = -z * sin (of_real pi * z) / of_real pi"
  2225 proof -
  2226   have "rGamma z * rGamma (-z) = -z * (rGamma z * rGamma (1 - z))"
  2227     using rGamma_plus1[of "-z", symmetric] by simp
  2228   also have "rGamma z * rGamma (1 - z) = sin (of_real pi * z) / of_real pi"
  2229     by (rule rGamma_reflection_complex)
  2230   finally show ?thesis by simp
  2231 qed
  2232 
  2233 lemma Gamma_reflection_complex':
  2234   "Gamma z * Gamma (- z :: complex) = - of_real pi / (z * sin (of_real pi * z))"
  2235   using rGamma_reflection_complex'[of z] by (force simp add: Gamma_def divide_simps mult_ac)
  2236 
  2237 
  2238 
  2239 lemma Gamma_one_half_real: "Gamma (1/2 :: real) = sqrt pi"
  2240 proof -
  2241   from Gamma_reflection_complex[of "1/2"] fraction_not_in_ints[where 'a = complex, of 2 1]
  2242     have "Gamma (1/2 :: complex)^2 = of_real pi" by (simp add: power2_eq_square)
  2243   hence "of_real pi = Gamma (complex_of_real (1/2))^2" by simp
  2244   also have "\<dots> = of_real ((Gamma (1/2))^2)" by (subst Gamma_complex_of_real) simp_all
  2245   finally have "Gamma (1/2)^2 = pi" by (subst (asm) of_real_eq_iff) simp_all
  2246   moreover have "Gamma (1/2 :: real) \<ge> 0" using Gamma_real_pos[of "1/2"] by simp
  2247   ultimately show ?thesis by (rule real_sqrt_unique [symmetric])
  2248 qed
  2249 
  2250 lemma Gamma_one_half_complex: "Gamma (1/2 :: complex) = of_real (sqrt pi)"
  2251 proof -
  2252   have "Gamma (1/2 :: complex) = Gamma (of_real (1/2))" by simp
  2253   also have "\<dots> = of_real (sqrt pi)" by (simp only: Gamma_complex_of_real Gamma_one_half_real)
  2254   finally show ?thesis .
  2255 qed
  2256 
  2257 lemma Gamma_legendre_duplication:
  2258   fixes z :: complex
  2259   assumes "z \<notin> \<int>\<^sub>\<le>\<^sub>0" "z + 1/2 \<notin> \<int>\<^sub>\<le>\<^sub>0"
  2260   shows "Gamma z * Gamma (z + 1/2) =
  2261              exp ((1 - 2*z) * of_real (ln 2)) * of_real (sqrt pi) * Gamma (2*z)"
  2262   using Gamma_legendre_duplication_aux[OF assms] by (simp add: Gamma_one_half_complex)
  2263 
  2264 end
  2265 
  2266 
  2267 subsection \<open>Limits and residues\<close>
  2268 
  2269 text \<open>
  2270   The inverse of the Gamma function has simple zeros:
  2271 \<close>
  2272 
  2273 lemma rGamma_zeros:
  2274   "(\<lambda>z. rGamma z / (z + of_nat n)) \<midarrow> (- of_nat n) \<rightarrow> ((-1)^n * fact n :: 'a :: Gamma)"
  2275 proof (subst tendsto_cong)
  2276   let ?f = "\<lambda>z. pochhammer z n * rGamma (z + of_nat (Suc n)) :: 'a"
  2277   from eventually_at_ball'[OF zero_less_one, of "- of_nat n :: 'a" UNIV]
  2278     show "eventually (\<lambda>z. rGamma z / (z + of_nat n) = ?f z) (at (- of_nat n))"
  2279     by (subst pochhammer_rGamma[of _ "Suc n"])
  2280        (auto elim!: eventually_mono simp: divide_simps pochhammer_rec' eq_neg_iff_add_eq_0)
  2281   have "isCont ?f (- of_nat n)" by (intro continuous_intros)
  2282   thus "?f \<midarrow> (- of_nat n) \<rightarrow> (- 1) ^ n * fact n" unfolding isCont_def
  2283     by (simp add: pochhammer_same)
  2284 qed
  2285 
  2286 
  2287 text \<open>
  2288   The simple zeros of the inverse of the Gamma function correspond to simple poles of the Gamma function,
  2289   and their residues can easily be computed from the limit we have just proven:
  2290 \<close>
  2291 
  2292 lemma Gamma_poles: "filterlim Gamma at_infinity (at (- of_nat n :: 'a :: Gamma))"
  2293 proof -
  2294   from eventually_at_ball'[OF zero_less_one, of "- of_nat n :: 'a" UNIV]
  2295     have "eventually (\<lambda>z. rGamma z \<noteq> (0 :: 'a)) (at (- of_nat n))"
  2296     by (auto elim!: eventually_mono nonpos_Ints_cases'
  2297              simp: rGamma_eq_zero_iff dist_of_nat dist_minus)
  2298   with isCont_rGamma[of "- of_nat n :: 'a", OF continuous_ident]
  2299     have "filterlim (\<lambda>z. inverse (rGamma z) :: 'a) at_infinity (at (- of_nat n))"
  2300     unfolding isCont_def by (intro filterlim_compose[OF filterlim_inverse_at_infinity])
  2301                             (simp_all add: filterlim_at)
  2302   moreover have "(\<lambda>z. inverse (rGamma z) :: 'a) = Gamma"
  2303     by (intro ext) (simp add: rGamma_inverse_Gamma)
  2304   ultimately show ?thesis by (simp only: )
  2305 qed
  2306 
  2307 lemma Gamma_residues:
  2308   "(\<lambda>z. Gamma z * (z + of_nat n)) \<midarrow> (- of_nat n) \<rightarrow> ((-1)^n / fact n :: 'a :: Gamma)"
  2309 proof (subst tendsto_cong)
  2310   let ?c = "(- 1) ^ n / fact n :: 'a"
  2311   from eventually_at_ball'[OF zero_less_one, of "- of_nat n :: 'a" UNIV]
  2312     show "eventually (\<lambda>z. Gamma z * (z + of_nat n) = inverse (rGamma z / (z + of_nat n)))
  2313             (at (- of_nat n))"
  2314     by (auto elim!: eventually_mono simp: divide_simps rGamma_inverse_Gamma)
  2315   have "(\<lambda>z. inverse (rGamma z / (z + of_nat n))) \<midarrow> (- of_nat n) \<rightarrow>
  2316           inverse ((- 1) ^ n * fact n :: 'a)"
  2317     by (intro tendsto_intros rGamma_zeros) simp_all
  2318   also have "inverse ((- 1) ^ n * fact n) = ?c"
  2319     by (simp_all add: field_simps power_mult_distrib [symmetric] del: power_mult_distrib)
  2320   finally show "(\<lambda>z. inverse (rGamma z / (z + of_nat n))) \<midarrow> (- of_nat n) \<rightarrow> ?c" .
  2321 qed
  2322 
  2323 
  2324 
  2325 subsection \<open>Alternative definitions\<close>
  2326 
  2327 
  2328 subsubsection \<open>Variant of the Euler form\<close>
  2329 
  2330 
  2331 definition Gamma_series_euler' where
  2332   "Gamma_series_euler' z n =
  2333      inverse z * (\<Prod>k=1..n. exp (z * of_real (ln (1 + inverse (of_nat k)))) / (1 + z / of_nat k))"
  2334 
  2335 context
  2336 begin
  2337 private lemma Gamma_euler'_aux1:
  2338   fixes z :: "'a :: {real_normed_field,banach}"
  2339   assumes n: "n > 0"
  2340   shows "exp (z * of_real (ln (of_nat n + 1))) = (\<Prod>k=1..n. exp (z * of_real (ln (1 + 1 / of_nat k))))"
  2341 proof -
  2342   have "(\<Prod>k=1..n. exp (z * of_real (ln (1 + 1 / of_nat k)))) =
  2343           exp (z * of_real (\<Sum>k = 1..n. ln (1 + 1 / real_of_nat k)))"
  2344     by (subst exp_setsum [symmetric]) (simp_all add: setsum_right_distrib)
  2345   also have "(\<Sum>k=1..n. ln (1 + 1 / of_nat k) :: real) = ln (\<Prod>k=1..n. 1 + 1 / real_of_nat k)"
  2346     by (subst ln_setprod [symmetric]) (auto intro!: add_pos_nonneg)
  2347   also have "(\<Prod>k=1..n. 1 + 1 / of_nat k :: real) = (\<Prod>k=1..n. (of_nat k + 1) / of_nat k)"
  2348     by (intro setprod.cong) (simp_all add: divide_simps)
  2349   also have "(\<Prod>k=1..n. (of_nat k + 1) / of_nat k :: real) = of_nat n + 1"
  2350     by (induction n) (simp_all add: setprod_nat_ivl_Suc' divide_simps)
  2351   finally show ?thesis ..
  2352 qed
  2353 
  2354 lemma Gamma_series_euler':
  2355   assumes z: "(z :: 'a :: Gamma) \<notin> \<int>\<^sub>\<le>\<^sub>0"
  2356   shows "(\<lambda>n. Gamma_series_euler' z n) \<longlonglongrightarrow> Gamma z"
  2357 proof (rule Gamma_seriesI, rule Lim_transform_eventually)
  2358   let ?f = "\<lambda>n. fact n * exp (z * of_real (ln (of_nat n + 1))) / pochhammer z (n + 1)"
  2359   let ?r = "\<lambda>n. ?f n / Gamma_series z n"
  2360   let ?r' = "\<lambda>n. exp (z * of_real (ln (of_nat (Suc n) / of_nat n)))"
  2361   from z have z': "z \<noteq> 0" by auto
  2362 
  2363   have "eventually (\<lambda>n. ?r' n = ?r n) sequentially" using eventually_gt_at_top[of "0::nat"]
  2364     using z by (auto simp: divide_simps Gamma_series_def ring_distribs exp_diff ln_div add_ac
  2365                      elim!: eventually_mono dest: pochhammer_eq_0_imp_nonpos_Int)
  2366   moreover have "?r' \<longlonglongrightarrow> exp (z * of_real (ln 1))"
  2367     by (intro tendsto_intros LIMSEQ_Suc_n_over_n) simp_all
  2368   ultimately show "?r \<longlonglongrightarrow> 1" by (force dest!: Lim_transform_eventually)
  2369 
  2370   from eventually_gt_at_top[of "0::nat"]
  2371     show "eventually (\<lambda>n. ?r n = Gamma_series_euler' z n / Gamma_series z n) sequentially"
  2372   proof eventually_elim
  2373     fix n :: nat assume n: "n > 0"
  2374     from n z' have "Gamma_series_euler' z n =
  2375       exp (z * of_real (ln (of_nat n + 1))) / (z * (\<Prod>k=1..n. (1 + z / of_nat k)))"
  2376       by (subst Gamma_euler'_aux1)
  2377          (simp_all add: Gamma_series_euler'_def setprod.distrib
  2378                         setprod_inversef[symmetric] divide_inverse)
  2379     also have "(\<Prod>k=1..n. (1 + z / of_nat k)) = pochhammer (z + 1) n / fact n"
  2380       by (cases n) (simp_all add: pochhammer_def fact_altdef setprod_shift_bounds_cl_Suc_ivl
  2381                                   setprod_dividef[symmetric] divide_simps add_ac)
  2382     also have "z * \<dots> = pochhammer z (Suc n) / fact n" by (simp add: pochhammer_rec)
  2383     finally show "?r n = Gamma_series_euler' z n / Gamma_series z n" by simp
  2384   qed
  2385 qed
  2386 
  2387 end
  2388 
  2389 
  2390 
  2391 subsubsection \<open>Weierstrass form\<close>
  2392 
  2393 definition Gamma_series_weierstrass :: "'a :: {banach,real_normed_field} \<Rightarrow> nat \<Rightarrow> 'a" where
  2394   "Gamma_series_weierstrass z n =
  2395      exp (-euler_mascheroni * z) / z * (\<Prod>k=1..n. exp (z / of_nat k) / (1 + z / of_nat k))"
  2396 
  2397 definition rGamma_series_weierstrass :: "'a :: {banach,real_normed_field} \<Rightarrow> nat \<Rightarrow> 'a" where
  2398   "rGamma_series_weierstrass z n =
  2399      exp (euler_mascheroni * z) * z * (\<Prod>k=1..n. (1 + z / of_nat k) * exp (-z / of_nat k))"
  2400 
  2401 lemma Gamma_series_weierstrass_nonpos_Ints:
  2402   "eventually (\<lambda>k. Gamma_series_weierstrass (- of_nat n) k = 0) sequentially"
  2403   using eventually_ge_at_top[of n] by eventually_elim (auto simp: Gamma_series_weierstrass_def)
  2404 
  2405 lemma rGamma_series_weierstrass_nonpos_Ints:
  2406   "eventually (\<lambda>k. rGamma_series_weierstrass (- of_nat n) k = 0) sequentially"
  2407   using eventually_ge_at_top[of n] by eventually_elim (auto simp: rGamma_series_weierstrass_def)
  2408 
  2409 lemma Gamma_weierstrass_complex: "Gamma_series_weierstrass z \<longlonglongrightarrow> Gamma (z :: complex)"
  2410 proof (cases "z \<in> \<int>\<^sub>\<le>\<^sub>0")
  2411   case True
  2412   then obtain n where "z = - of_nat n" by (elim nonpos_Ints_cases')
  2413   also from True have "Gamma_series_weierstrass \<dots> \<longlonglongrightarrow> Gamma z"
  2414     by (simp add: tendsto_cong[OF Gamma_series_weierstrass_nonpos_Ints] Gamma_nonpos_Int)
  2415   finally show ?thesis .
  2416 next
  2417   case False
  2418   hence z: "z \<noteq> 0" by auto
  2419   let ?f = "(\<lambda>x. \<Prod>x = Suc 0..x. exp (z / of_nat x) / (1 + z / of_nat x))"
  2420   have A: "exp (ln (1 + z / of_nat n)) = (1 + z / of_nat n)" if "n \<ge> 1" for n :: nat
  2421     using False that by (subst exp_Ln) (auto simp: field_simps dest!: plus_of_nat_eq_0_imp)
  2422   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"
  2423     using ln_Gamma_series'_aux[OF False]
  2424     by (simp only: atLeastLessThanSuc_atLeastAtMost [symmetric] One_nat_def
  2425                    setsum_shift_bounds_Suc_ivl sums_def atLeast0LessThan)
  2426   from tendsto_exp[OF this] False z have "?f \<longlonglongrightarrow> z * exp (euler_mascheroni * z) * Gamma z"
  2427     by (simp add: exp_add exp_setsum exp_diff mult_ac Gamma_complex_altdef A)
  2428   from tendsto_mult[OF tendsto_const[of "exp (-euler_mascheroni * z) / z"] this] z
  2429     show "Gamma_series_weierstrass z \<longlonglongrightarrow> Gamma z"
  2430     by (simp add: exp_minus divide_simps Gamma_series_weierstrass_def [abs_def])
  2431 qed
  2432 
  2433 lemma tendsto_complex_of_real_iff: "((\<lambda>x. complex_of_real (f x)) \<longlongrightarrow> of_real c) F = (f \<longlongrightarrow> c) F"
  2434   by (rule tendsto_of_real_iff)
  2435 
  2436 lemma Gamma_weierstrass_real: "Gamma_series_weierstrass x \<longlonglongrightarrow> Gamma (x :: real)"
  2437   using Gamma_weierstrass_complex[of "of_real x"] unfolding Gamma_series_weierstrass_def[abs_def]
  2438   by (subst tendsto_complex_of_real_iff [symmetric])
  2439      (simp_all add: exp_of_real[symmetric] Gamma_complex_of_real)
  2440 
  2441 lemma rGamma_weierstrass_complex: "rGamma_series_weierstrass z \<longlonglongrightarrow> rGamma (z :: complex)"
  2442 proof (cases "z \<in> \<int>\<^sub>\<le>\<^sub>0")
  2443   case True
  2444   then obtain n where "z = - of_nat n" by (elim nonpos_Ints_cases')
  2445   also from True have "rGamma_series_weierstrass \<dots> \<longlonglongrightarrow> rGamma z"
  2446     by (simp add: tendsto_cong[OF rGamma_series_weierstrass_nonpos_Ints] rGamma_nonpos_Int)
  2447   finally show ?thesis .
  2448 next
  2449   case False
  2450   have "rGamma_series_weierstrass z = (\<lambda>n. inverse (Gamma_series_weierstrass z n))"
  2451     by (simp add: rGamma_series_weierstrass_def[abs_def] Gamma_series_weierstrass_def
  2452                   exp_minus divide_inverse setprod_inversef[symmetric] mult_ac)
  2453   also from False have "\<dots> \<longlonglongrightarrow> inverse (Gamma z)"
  2454     by (intro tendsto_intros Gamma_weierstrass_complex) (simp add: Gamma_eq_zero_iff)
  2455   finally show ?thesis by (simp add: Gamma_def)
  2456 qed
  2457 
  2458 subsubsection \<open>Binomial coefficient form\<close>
  2459 
  2460 lemma Gamma_binomial:
  2461   "(\<lambda>n. ((z + of_nat n) gchoose n) * exp (-z * of_real (ln (of_nat n)))) \<longlonglongrightarrow> rGamma (z+1)"
  2462 proof (cases "z = 0")
  2463   case False
  2464   show ?thesis
  2465   proof (rule Lim_transform_eventually)
  2466     let ?powr = "\<lambda>a b. exp (b * of_real (ln (of_nat a)))"
  2467     show "eventually (\<lambda>n. rGamma_series z n / z =
  2468             ((z + of_nat n) gchoose n) * ?powr n (-z)) sequentially"
  2469     proof (intro always_eventually allI)
  2470       fix n :: nat
  2471       from False have "((z + of_nat n) gchoose n) = pochhammer z (Suc n) / z / fact n"
  2472         by (simp add: gbinomial_pochhammer' pochhammer_rec)
  2473       also have "pochhammer z (Suc n) / z / fact n * ?powr n (-z) = rGamma_series z n / z"
  2474         by (simp add: rGamma_series_def divide_simps exp_minus)
  2475       finally show "rGamma_series z n / z = ((z + of_nat n) gchoose n) * ?powr n (-z)" ..
  2476     qed
  2477 
  2478     from False have "(\<lambda>n. rGamma_series z n / z) \<longlonglongrightarrow> rGamma z / z" by (intro tendsto_intros)
  2479     also from False have "rGamma z / z = rGamma (z + 1)" using rGamma_plus1[of z]
  2480       by (simp add: field_simps)
  2481     finally show "(\<lambda>n. rGamma_series z n / z) \<longlonglongrightarrow> rGamma (z+1)" .
  2482   qed
  2483 qed (simp_all add: binomial_gbinomial [symmetric])
  2484 
  2485 lemma fact_binomial_limit:
  2486   "(\<lambda>n. of_nat ((k + n) choose n) / of_nat (n ^ k) :: 'a :: Gamma) \<longlonglongrightarrow> 1 / fact k"
  2487 proof (rule Lim_transform_eventually)
  2488   have "(\<lambda>n. of_nat ((k + n) choose n) / of_real (exp (of_nat k * ln (real_of_nat n))))
  2489             \<longlonglongrightarrow> 1 / Gamma (of_nat (Suc k) :: 'a)" (is "?f \<longlonglongrightarrow> _")
  2490     using Gamma_binomial[of "of_nat k :: 'a"]
  2491     by (simp add: binomial_gbinomial add_ac Gamma_def divide_simps exp_of_real [symmetric] exp_minus)
  2492   also have "Gamma (of_nat (Suc k)) = fact k" by (rule Gamma_fact)
  2493   finally show "?f \<longlonglongrightarrow> 1 / fact k" .
  2494 
  2495   show "eventually (\<lambda>n. ?f n = of_nat ((k + n) choose n) / of_nat (n ^ k)) sequentially"
  2496     using eventually_gt_at_top[of "0::nat"]
  2497   proof eventually_elim
  2498     fix n :: nat assume n: "n > 0"
  2499     from n have "exp (real_of_nat k * ln (real_of_nat n)) = real_of_nat (n^k)"
  2500       by (simp add: exp_of_nat_mult)
  2501     thus "?f n = of_nat ((k + n) choose n) / of_nat (n ^ k)" by simp
  2502   qed
  2503 qed
  2504 
  2505 lemma binomial_asymptotic:
  2506   "(\<lambda>n. of_nat ((k + n) choose n) / (of_nat (n ^ k) / fact k) :: 'a :: Gamma) \<longlonglongrightarrow> 1"
  2507   using tendsto_mult[OF fact_binomial_limit[of k] tendsto_const[of "fact k :: 'a"]] by simp
  2508 
  2509 
  2510 subsection \<open>The Weierstraß product formula for the sine\<close>
  2511 
  2512 lemma sin_product_formula_complex:
  2513   fixes z :: complex
  2514   shows "(\<lambda>n. of_real pi * z * (\<Prod>k=1..n. 1 - z^2 / of_nat k^2)) \<longlonglongrightarrow> sin (of_real pi * z)"
  2515 proof -
  2516   let ?f = "rGamma_series_weierstrass"
  2517   have "(\<lambda>n. (- of_real pi * inverse z) * (?f z n * ?f (- z) n))
  2518             \<longlonglongrightarrow> (- of_real pi * inverse z) * (rGamma z * rGamma (- z))"
  2519     by (intro tendsto_intros rGamma_weierstrass_complex)
  2520   also have "(\<lambda>n. (- of_real pi * inverse z) * (?f z n * ?f (-z) n)) =
  2521                     (\<lambda>n. of_real pi * z * (\<Prod>k=1..n. 1 - z^2 / of_nat k ^ 2))"
  2522   proof
  2523     fix n :: nat
  2524     have "(- of_real pi * inverse z) * (?f z n * ?f (-z) n) =
  2525               of_real pi * z * (\<Prod>k=1..n. (of_nat k - z) * (of_nat k + z) / of_nat k ^ 2)"
  2526       by (simp add: rGamma_series_weierstrass_def mult_ac exp_minus
  2527                     divide_simps setprod.distrib[symmetric] power2_eq_square)
  2528     also have "(\<Prod>k=1..n. (of_nat k - z) * (of_nat k + z) / of_nat k ^ 2) =
  2529                  (\<Prod>k=1..n. 1 - z^2 / of_nat k ^ 2)"
  2530       by (intro setprod.cong) (simp_all add: power2_eq_square field_simps)
  2531     finally show "(- of_real pi * inverse z) * (?f z n * ?f (-z) n) = of_real pi * z * \<dots>"
  2532       by (simp add: divide_simps)
  2533   qed
  2534   also have "(- of_real pi * inverse z) * (rGamma z * rGamma (- z)) = sin (of_real pi * z)"
  2535     by (subst rGamma_reflection_complex') (simp add: divide_simps)
  2536   finally show ?thesis .
  2537 qed
  2538 
  2539 lemma sin_product_formula_real:
  2540   "(\<lambda>n. pi * (x::real) * (\<Prod>k=1..n. 1 - x^2 / of_nat k^2)) \<longlonglongrightarrow> sin (pi * x)"
  2541 proof -
  2542   from sin_product_formula_complex[of "of_real x"]
  2543     have "(\<lambda>n. of_real pi * of_real x * (\<Prod>k=1..n. 1 - (of_real x)^2 / (of_nat k)^2))
  2544               \<longlonglongrightarrow> sin (of_real pi * of_real x :: complex)" (is "?f \<longlonglongrightarrow> ?y") .
  2545   also have "?f = (\<lambda>n. of_real (pi * x * (\<Prod>k=1..n. 1 - x^2 / (of_nat k^2))))" by simp
  2546   also have "?y = of_real (sin (pi * x))" by (simp only: sin_of_real [symmetric] of_real_mult)
  2547   finally show ?thesis by (subst (asm) tendsto_of_real_iff)
  2548 qed
  2549 
  2550 lemma sin_product_formula_real':
  2551   assumes "x \<noteq> (0::real)"
  2552   shows   "(\<lambda>n. (\<Prod>k=1..n. 1 - x^2 / of_nat k^2)) \<longlonglongrightarrow> sin (pi * x) / (pi * x)"
  2553   using tendsto_divide[OF sin_product_formula_real[of x] tendsto_const[of "pi * x"]] assms
  2554   by simp
  2555 
  2556 
  2557 subsection \<open>The Solution to the Basel problem\<close>
  2558 
  2559 theorem inverse_squares_sums: "(\<lambda>n. 1 / (n + 1)\<^sup>2) sums (pi\<^sup>2 / 6)"
  2560 proof -
  2561   def P \<equiv> "\<lambda>x n. (\<Prod>k=1..n. 1 - x^2 / of_nat k^2 :: real)"
  2562   def K \<equiv> "\<Sum>n. inverse (real_of_nat (Suc n))^2"
  2563   def f \<equiv> "\<lambda>x. \<Sum>n. P x n / of_nat (Suc n)^2"
  2564   def g \<equiv> "\<lambda>x. (1 - sin (pi * x) / (pi * x))"
  2565 
  2566   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
  2567   proof (cases "x = 0")
  2568     assume x: "x = 0"
  2569     have "summable (\<lambda>n. inverse ((real_of_nat (Suc n))\<^sup>2))"
  2570       using inverse_power_summable[of 2] by (subst summable_Suc_iff) simp
  2571     thus ?thesis by (simp add: x g_def P_def K_def inverse_eq_divide power_divide summable_sums)
  2572   next
  2573     assume x: "x \<noteq> 0"
  2574     have "(\<lambda>n. P x n - P x (Suc n)) sums (P x 0 - sin (pi * x) / (pi * x))"
  2575       unfolding P_def using x by (intro telescope_sums' sin_product_formula_real')
  2576     also have "(\<lambda>n. P x n - P x (Suc n)) = (\<lambda>n. (x^2 / of_nat (Suc n)^2) * P x n)"
  2577       unfolding P_def by (simp add: setprod_nat_ivl_Suc' algebra_simps)
  2578     also have "P x 0 = 1" by (simp add: P_def)
  2579     finally have "(\<lambda>n. x\<^sup>2 / (of_nat (Suc n))\<^sup>2 * P x n) sums (1 - sin (pi * x) / (pi * x))" .
  2580     from sums_divide[OF this, of "x^2"] x show ?thesis unfolding g_def by simp
  2581   qed
  2582 
  2583   have "continuous_on (ball 0 1) f"
  2584   proof (rule uniform_limit_theorem; (intro always_eventually allI)?)
  2585     show "uniform_limit (ball 0 1) (\<lambda>n x. \<Sum>k<n. P x k / of_nat (Suc k)^2) f sequentially"
  2586     proof (unfold f_def, rule weierstrass_m_test)
  2587       fix n :: nat and x :: real assume x: "x \<in> ball 0 1"
  2588       {
  2589         fix k :: nat assume k: "k \<ge> 1"
  2590         from x have "x^2 < 1" by (auto simp: dist_0_norm abs_square_less_1)
  2591         also from k have "\<dots> \<le> of_nat k^2" by simp
  2592         finally have "(1 - x^2 / of_nat k^2) \<in> {0..1}" using k
  2593           by (simp_all add: field_simps del: of_nat_Suc)
  2594       }
  2595       hence "(\<Prod>k=1..n. abs (1 - x^2 / of_nat k^2)) \<le> (\<Prod>k=1..n. 1)" by (intro setprod_mono) simp
  2596       thus "norm (P x n / (of_nat (Suc n)^2)) \<le> 1 / of_nat (Suc n)^2"
  2597         unfolding P_def by (simp add: field_simps abs_setprod del: of_nat_Suc)
  2598     qed (subst summable_Suc_iff, insert inverse_power_summable[of 2], simp add: inverse_eq_divide)
  2599   qed (auto simp: P_def intro!: continuous_intros)
  2600   hence "isCont f 0" by (subst (asm) continuous_on_eq_continuous_at) simp_all
  2601   hence "(f \<midarrow> 0 \<rightarrow> f 0)" by (simp add: isCont_def)
  2602   also have "f 0 = K" unfolding f_def P_def K_def by (simp add: inverse_eq_divide power_divide)
  2603   finally have "f \<midarrow> 0 \<rightarrow> K" .
  2604 
  2605   moreover have "f \<midarrow> 0 \<rightarrow> pi^2 / 6"
  2606   proof (rule Lim_transform_eventually)
  2607     def f' \<equiv> "\<lambda>x. \<Sum>n. - sin_coeff (n+3) * pi ^ (n+2) * x^n"
  2608     have "eventually (\<lambda>x. x \<noteq> (0::real)) (at 0)"
  2609       by (auto simp add: eventually_at intro!: exI[of _ 1])
  2610     thus "eventually (\<lambda>x. f' x = f x) (at 0)"
  2611     proof eventually_elim
  2612       fix x :: real assume x: "x \<noteq> 0"
  2613       have "sin_coeff 1 = (1 :: real)" "sin_coeff 2 = (0::real)" by (simp_all add: sin_coeff_def)
  2614       with sums_split_initial_segment[OF sums_minus[OF sin_converges], of 3 "pi*x"]
  2615       have "(\<lambda>n. - (sin_coeff (n+3) * (pi*x)^(n+3))) sums (pi * x - sin (pi*x))"
  2616         by (simp add: eval_nat_numeral)
  2617       from sums_divide[OF this, of "x^3 * pi"] x
  2618         have "(\<lambda>n. - (sin_coeff (n+3) * pi^(n+2) * x^n)) sums ((1 - sin (pi*x) / (pi*x)) / x^2)"
  2619         by (simp add: divide_simps eval_nat_numeral power_mult_distrib mult_ac)
  2620       with x have "(\<lambda>n. - (sin_coeff (n+3) * pi^(n+2) * x^n)) sums (g x / x^2)"
  2621         by (simp add: g_def)
  2622       hence "f' x = g x / x^2" by (simp add: sums_iff f'_def)
  2623       also have "\<dots> = f x" using sums[of x] x by (simp add: sums_iff g_def f_def)
  2624       finally show "f' x = f x" .
  2625     qed
  2626 
  2627     have "isCont f' 0" unfolding f'_def
  2628     proof (intro isCont_powser_converges_everywhere)
  2629       fix x :: real show "summable (\<lambda>n. -sin_coeff (n+3) * pi^(n+2) * x^n)"
  2630       proof (cases "x = 0")
  2631         assume x: "x \<noteq> 0"
  2632         from summable_divide[OF sums_summable[OF sums_split_initial_segment[OF
  2633                sin_converges[of "pi*x"]], of 3], of "-pi*x^3"] x
  2634           show ?thesis by (simp add: mult_ac power_mult_distrib divide_simps eval_nat_numeral)
  2635       qed (simp only: summable_0_powser)
  2636     qed
  2637     hence "f' \<midarrow> 0 \<rightarrow> f' 0" by (simp add: isCont_def)
  2638     also have "f' 0 = pi * pi / fact 3" unfolding f'_def
  2639       by (subst powser_zero) (simp add: sin_coeff_def)
  2640     finally show "f' \<midarrow> 0 \<rightarrow> pi^2 / 6" by (simp add: eval_nat_numeral)
  2641   qed
  2642 
  2643   ultimately have "K = pi^2 / 6" by (rule LIM_unique)
  2644   moreover from inverse_power_summable[of 2]
  2645     have "summable (\<lambda>n. (inverse (real_of_nat (Suc n)))\<^sup>2)"
  2646     by (subst summable_Suc_iff) (simp add: power_inverse)
  2647   ultimately show ?thesis unfolding K_def
  2648     by (auto simp add: sums_iff power_divide inverse_eq_divide)
  2649 qed
  2650 
  2651 
  2652 end