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