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