src/HOL/Analysis/Gamma_Function.thy
author wenzelm
Mon Mar 25 17:21:26 2019 +0100 (3 weeks ago)
changeset 69981 3dced198b9ec
parent 69597 ff784d5a5bfb
child 70113 c8deb8ba6d05
permissions -rw-r--r--
more strict AFP properties;
wenzelm@63992
     1
(*  Title:    HOL/Analysis/Gamma_Function.thy
hoelzl@62055
     2
    Author:   Manuel Eberl, TU München
eberlm@62049
     3
*)
hoelzl@62055
     4
hoelzl@62055
     5
section \<open>The Gamma Function\<close>
hoelzl@62055
     6
hoelzl@63594
     7
theory Gamma_Function
paulson@62131
     8
imports
eberlm@66286
     9
  Conformal_Mappings
hoelzl@63594
    10
  Summation_Tests
eberlm@62049
    11
  Harmonic_Numbers
wenzelm@66453
    12
  "HOL-Library.Nonpos_Ints"
wenzelm@66453
    13
  "HOL-Library.Periodic_Fun"
eberlm@62049
    14
begin
eberlm@62049
    15
hoelzl@62055
    16
text \<open>
paulson@62131
    17
  Several equivalent definitions of the Gamma function and its
hoelzl@62055
    18
  most important properties. Also contains the definition and some properties
hoelzl@62055
    19
  of the log-Gamma function and the Digamma function and the other Polygamma functions.
paulson@62131
    20
hoelzl@62055
    21
  Based on the Gamma function, we also prove the Weierstraß product form of the
paulson@62131
    22
  sin function and, based on this, the solution of the Basel problem (the
wenzelm@69597
    23
  sum over all \<^term>\<open>1 / (n::nat)^2\<close>.
paulson@62131
    24
\<close>
hoelzl@62055
    25
paulson@62131
    26
lemma pochhammer_eq_0_imp_nonpos_Int:
eberlm@62049
    27
  "pochhammer (x::'a::field_char_0) n = 0 \<Longrightarrow> x \<in> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
    28
  by (auto simp: pochhammer_eq_0_iff)
paulson@62131
    29
eberlm@62049
    30
lemma closed_nonpos_Ints [simp]: "closed (\<int>\<^sub>\<le>\<^sub>0 :: 'a :: real_normed_algebra_1 set)"
eberlm@62049
    31
proof -
paulson@62131
    32
  have "\<int>\<^sub>\<le>\<^sub>0 = (of_int ` {n. n \<le> 0} :: 'a set)"
eberlm@62049
    33
    by (auto elim!: nonpos_Ints_cases intro!: nonpos_Ints_of_int)
eberlm@62049
    34
  also have "closed \<dots>" by (rule closed_of_int_image)
eberlm@62049
    35
  finally show ?thesis .
eberlm@62049
    36
qed
eberlm@62049
    37
eberlm@62049
    38
lemma plus_one_in_nonpos_Ints_imp: "z + 1 \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> z \<in> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
    39
  using nonpos_Ints_diff_Nats[of "z+1" "1"] by simp_all
eberlm@62049
    40
eberlm@63295
    41
lemma of_int_in_nonpos_Ints_iff:
eberlm@63295
    42
  "(of_int n :: 'a :: ring_char_0) \<in> \<int>\<^sub>\<le>\<^sub>0 \<longleftrightarrow> n \<le> 0"
eberlm@63295
    43
  by (auto simp: nonpos_Ints_def)
eberlm@63295
    44
eberlm@63295
    45
lemma one_plus_of_int_in_nonpos_Ints_iff:
eberlm@63295
    46
  "(1 + of_int n :: 'a :: ring_char_0) \<in> \<int>\<^sub>\<le>\<^sub>0 \<longleftrightarrow> n \<le> -1"
eberlm@63295
    47
proof -
eberlm@63295
    48
  have "1 + of_int n = (of_int (n + 1) :: 'a)" by simp
eberlm@63295
    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
eberlm@63295
    50
  also have "\<dots> \<longleftrightarrow> n \<le> -1" by presburger
eberlm@63295
    51
  finally show ?thesis .
eberlm@63295
    52
qed
eberlm@63295
    53
eberlm@63295
    54
lemma one_minus_of_nat_in_nonpos_Ints_iff:
eberlm@63295
    55
  "(1 - of_nat n :: 'a :: ring_char_0) \<in> \<int>\<^sub>\<le>\<^sub>0 \<longleftrightarrow> n > 0"
eberlm@63295
    56
proof -
eberlm@63295
    57
  have "(1 - of_nat n :: 'a) = of_int (1 - int n)" by simp
eberlm@63295
    58
  also have "\<dots> \<in> \<int>\<^sub>\<le>\<^sub>0 \<longleftrightarrow> n > 0" by (subst of_int_in_nonpos_Ints_iff) presburger
eberlm@63295
    59
  finally show ?thesis .
eberlm@63295
    60
qed
eberlm@63295
    61
eberlm@62049
    62
lemma fraction_not_in_ints:
eberlm@62049
    63
  assumes "\<not>(n dvd m)" "n \<noteq> 0"
eberlm@62049
    64
  shows   "of_int m / of_int n \<notin> (\<int> :: 'a :: {division_ring,ring_char_0} set)"
eberlm@62049
    65
proof
eberlm@62049
    66
  assume "of_int m / (of_int n :: 'a) \<in> \<int>"
eberlm@62049
    67
  then obtain k where "of_int m / of_int n = (of_int k :: 'a)" by (elim Ints_cases)
eberlm@62049
    68
  with assms have "of_int m = (of_int (k * n) :: 'a)" by (auto simp add: divide_simps)
eberlm@62049
    69
  hence "m = k * n" by (subst (asm) of_int_eq_iff)
eberlm@62049
    70
  hence "n dvd m" by simp
eberlm@62049
    71
  with assms(1) show False by contradiction
eberlm@62049
    72
qed
eberlm@62049
    73
eberlm@63317
    74
lemma fraction_not_in_nats:
eberlm@63317
    75
  assumes "\<not>n dvd m" "n \<noteq> 0"
eberlm@63317
    76
  shows   "of_int m / of_int n \<notin> (\<nat> :: 'a :: {division_ring,ring_char_0} set)"
eberlm@63317
    77
proof
eberlm@63317
    78
  assume "of_int m / of_int n \<in> (\<nat> :: 'a set)"
eberlm@63317
    79
  also note Nats_subset_Ints
eberlm@63317
    80
  finally have "of_int m / of_int n \<in> (\<int> :: 'a set)" .
eberlm@63317
    81
  moreover have "of_int m / of_int n \<notin> (\<int> :: 'a set)"
eberlm@63317
    82
    using assms by (intro fraction_not_in_ints)
eberlm@63317
    83
  ultimately show False by contradiction
eberlm@63317
    84
qed
eberlm@63317
    85
eberlm@62049
    86
lemma not_in_Ints_imp_not_in_nonpos_Ints: "z \<notin> \<int> \<Longrightarrow> z \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
    87
  by (auto simp: Ints_def nonpos_Ints_def)
eberlm@62049
    88
eberlm@62049
    89
lemma double_in_nonpos_Ints_imp:
eberlm@62049
    90
  assumes "2 * (z :: 'a :: field_char_0) \<in> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
    91
  shows   "z \<in> \<int>\<^sub>\<le>\<^sub>0 \<or> z + 1/2 \<in> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
    92
proof-
eberlm@62049
    93
  from assms obtain k where k: "2 * z = - of_nat k" by (elim nonpos_Ints_cases')
eberlm@62049
    94
  thus ?thesis by (cases "even k") (auto elim!: evenE oddE simp: field_simps)
eberlm@62049
    95
qed
eberlm@62049
    96
eberlm@62049
    97
eberlm@62049
    98
lemma sin_series: "(\<lambda>n. ((-1)^n / fact (2*n+1)) *\<^sub>R z^(2*n+1)) sums sin z"
eberlm@62049
    99
proof -
eberlm@62049
   100
  from sin_converges[of z] have "(\<lambda>n. sin_coeff n *\<^sub>R z^n) sums sin z" .
paulson@62131
   101
  also have "(\<lambda>n. sin_coeff n *\<^sub>R z^n) sums sin z \<longleftrightarrow>
eberlm@62049
   102
                 (\<lambda>n. ((-1)^n / fact (2*n+1)) *\<^sub>R z^(2*n+1)) sums sin z"
eberlm@62049
   103
    by (subst sums_mono_reindex[of "\<lambda>n. 2*n+1", symmetric])
eberlm@66447
   104
       (auto simp: sin_coeff_def strict_mono_def ac_simps elim!: oddE)
eberlm@62049
   105
  finally show ?thesis .
eberlm@62049
   106
qed
eberlm@62049
   107
eberlm@62049
   108
lemma cos_series: "(\<lambda>n. ((-1)^n / fact (2*n)) *\<^sub>R z^(2*n)) sums cos z"
eberlm@62049
   109
proof -
eberlm@62049
   110
  from cos_converges[of z] have "(\<lambda>n. cos_coeff n *\<^sub>R z^n) sums cos z" .
paulson@62131
   111
  also have "(\<lambda>n. cos_coeff n *\<^sub>R z^n) sums cos z \<longleftrightarrow>
eberlm@62049
   112
                 (\<lambda>n. ((-1)^n / fact (2*n)) *\<^sub>R z^(2*n)) sums cos z"
eberlm@62049
   113
    by (subst sums_mono_reindex[of "\<lambda>n. 2*n", symmetric])
eberlm@66447
   114
       (auto simp: cos_coeff_def strict_mono_def ac_simps elim!: evenE)
eberlm@62049
   115
  finally show ?thesis .
eberlm@62049
   116
qed
eberlm@62049
   117
eberlm@62049
   118
lemma sin_z_over_z_series:
eberlm@62049
   119
  fixes z :: "'a :: {real_normed_field,banach}"
eberlm@62049
   120
  assumes "z \<noteq> 0"
eberlm@62049
   121
  shows   "(\<lambda>n. (-1)^n / fact (2*n+1) * z^(2*n)) sums (sin z / z)"
eberlm@62049
   122
proof -
eberlm@62049
   123
  from sin_series[of z] have "(\<lambda>n. z * ((-1)^n / fact (2*n+1)) * z^(2*n)) sums sin z"
eberlm@62049
   124
    by (simp add: field_simps scaleR_conv_of_real)
eberlm@62049
   125
  from sums_mult[OF this, of "inverse z"] and assms show ?thesis
eberlm@62049
   126
    by (simp add: field_simps)
eberlm@62049
   127
qed
eberlm@62049
   128
eberlm@62049
   129
lemma sin_z_over_z_series':
eberlm@62049
   130
  fixes z :: "'a :: {real_normed_field,banach}"
eberlm@62049
   131
  assumes "z \<noteq> 0"
eberlm@62049
   132
  shows   "(\<lambda>n. sin_coeff (n+1) *\<^sub>R z^n) sums (sin z / z)"
eberlm@62049
   133
proof -
paulson@62131
   134
  from sums_split_initial_segment[OF sin_converges[of z], of 1]
eberlm@62049
   135
    have "(\<lambda>n. z * (sin_coeff (n+1) *\<^sub>R z ^ n)) sums sin z" by simp
eberlm@62049
   136
  from sums_mult[OF this, of "inverse z"] assms show ?thesis by (simp add: field_simps)
eberlm@62049
   137
qed
eberlm@62049
   138
eberlm@62049
   139
lemma has_field_derivative_sin_z_over_z:
eberlm@62049
   140
  fixes A :: "'a :: {real_normed_field,banach} set"
eberlm@62049
   141
  shows "((\<lambda>z. if z = 0 then 1 else sin z / z) has_field_derivative 0) (at 0 within A)"
eberlm@62049
   142
      (is "(?f has_field_derivative ?f') _")
eberlm@62049
   143
proof (rule has_field_derivative_at_within)
paulson@62131
   144
  have "((\<lambda>z::'a. \<Sum>n. of_real (sin_coeff (n+1)) * z^n)
eberlm@62049
   145
            has_field_derivative (\<Sum>n. diffs (\<lambda>n. of_real (sin_coeff (n+1))) n * 0^n)) (at 0)"
eberlm@62049
   146
  proof (rule termdiffs_strong)
eberlm@62049
   147
    from summable_ignore_initial_segment[OF sums_summable[OF sin_converges[of "1::'a"]], of 1]
eberlm@62049
   148
      show "summable (\<lambda>n. of_real (sin_coeff (n+1)) * (1::'a)^n)" by (simp add: of_real_def)
eberlm@62049
   149
  qed simp
eberlm@62049
   150
  also have "(\<lambda>z::'a. \<Sum>n. of_real (sin_coeff (n+1)) * z^n) = ?f"
eberlm@62049
   151
  proof
eberlm@62049
   152
    fix z
eberlm@62049
   153
    show "(\<Sum>n. of_real (sin_coeff (n+1)) * z^n)  = ?f z"
paulson@62131
   154
      by (cases "z = 0") (insert sin_z_over_z_series'[of z],
haftmann@66936
   155
            simp_all add: scaleR_conv_of_real sums_iff sin_coeff_def)
eberlm@62049
   156
  qed
paulson@62131
   157
  also have "(\<Sum>n. diffs (\<lambda>n. of_real (sin_coeff (n + 1))) n * (0::'a) ^ n) =
haftmann@66936
   158
                 diffs (\<lambda>n. of_real (sin_coeff (Suc n))) 0" by simp
eberlm@62049
   159
  also have "\<dots> = 0" by (simp add: sin_coeff_def diffs_def)
eberlm@62049
   160
  finally show "((\<lambda>z::'a. if z = 0 then 1 else sin z / z) has_field_derivative 0) (at 0)" .
eberlm@62049
   161
qed
eberlm@62049
   162
eberlm@62049
   163
lemma round_Re_minimises_norm:
eberlm@62049
   164
  "norm ((z::complex) - of_int m) \<ge> norm (z - of_int (round (Re z)))"
eberlm@62049
   165
proof -
eberlm@62049
   166
  let ?n = "round (Re z)"
eberlm@62049
   167
  have "norm (z - of_int ?n) = sqrt ((Re z - of_int ?n)\<^sup>2 + (Im z)\<^sup>2)"
eberlm@62049
   168
    by (simp add: cmod_def)
eberlm@62049
   169
  also have "\<bar>Re z - of_int ?n\<bar> \<le> \<bar>Re z - of_int m\<bar>" by (rule round_diff_minimal)
eberlm@62049
   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)"
eberlm@62049
   171
    by (intro real_sqrt_le_mono add_mono) (simp_all add: abs_le_square_iff)
eberlm@62049
   172
  also have "\<dots> = norm (z - of_int m)" by (simp add: cmod_def)
eberlm@62049
   173
  finally show ?thesis .
eberlm@62049
   174
qed
eberlm@62049
   175
eberlm@62049
   176
lemma Re_pos_in_ball:
eberlm@62049
   177
  assumes "Re z > 0" "t \<in> ball z (Re z/2)"
eberlm@62049
   178
  shows   "Re t > 0"
eberlm@62049
   179
proof -
eberlm@62049
   180
  have "Re (z - t) \<le> norm (z - t)" by (rule complex_Re_le_cmod)
eberlm@62049
   181
  also from assms have "\<dots> < Re z / 2" by (simp add: dist_complex_def)
eberlm@62049
   182
  finally show "Re t > 0" using assms by simp
eberlm@62049
   183
qed
eberlm@62049
   184
eberlm@62049
   185
lemma no_nonpos_Int_in_ball_complex:
eberlm@62049
   186
  assumes "Re z > 0" "t \<in> ball z (Re z/2)"
eberlm@62049
   187
  shows   "t \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
   188
  using Re_pos_in_ball[OF assms] by (force elim!: nonpos_Ints_cases)
eberlm@62049
   189
paulson@62131
   190
lemma no_nonpos_Int_in_ball:
eberlm@62049
   191
  assumes "t \<in> ball z (dist z (round (Re z)))"
eberlm@62049
   192
  shows   "t \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
   193
proof
eberlm@62049
   194
  assume "t \<in> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
   195
  then obtain n where "t = of_int n" by (auto elim!: nonpos_Ints_cases)
eberlm@62049
   196
  have "dist z (of_int n) \<le> dist z t + dist t (of_int n)" by (rule dist_triangle)
eberlm@62049
   197
  also from assms have "dist z t < dist z (round (Re z))" by simp
eberlm@62049
   198
  also have "\<dots> \<le> dist z (of_int n)"
eberlm@62049
   199
    using round_Re_minimises_norm[of z] by (simp add: dist_complex_def)
eberlm@62049
   200
  finally have "dist t (of_int n) > 0" by simp
wenzelm@62072
   201
  with \<open>t = of_int n\<close> show False by simp
eberlm@62049
   202
qed
eberlm@62049
   203
eberlm@62049
   204
lemma no_nonpos_Int_in_ball':
eberlm@62049
   205
  assumes "(z :: 'a :: {euclidean_space,real_normed_algebra_1}) \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
   206
  obtains d where "d > 0" "\<And>t. t \<in> ball z d \<Longrightarrow> t \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
   207
proof (rule that)
eberlm@62049
   208
  from assms show "setdist {z} \<int>\<^sub>\<le>\<^sub>0 > 0" by (subst setdist_gt_0_compact_closed) auto
eberlm@62049
   209
next
eberlm@62049
   210
  fix t assume "t \<in> ball z (setdist {z} \<int>\<^sub>\<le>\<^sub>0)"
eberlm@62049
   211
  thus "t \<notin> \<int>\<^sub>\<le>\<^sub>0" using setdist_le_dist[of z "{z}" t "\<int>\<^sub>\<le>\<^sub>0"] by force
eberlm@62049
   212
qed
eberlm@62049
   213
paulson@62131
   214
lemma no_nonpos_Real_in_ball:
paulson@62131
   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)"
paulson@62131
   216
  shows   "t \<notin> \<real>\<^sub>\<le>\<^sub>0"
paulson@62131
   217
using z
eberlm@62049
   218
proof (cases "Im z = 0")
eberlm@62049
   219
  assume A: "Im z = 0"
paulson@62131
   220
  with z have "Re z > 0" by (force simp add: complex_nonpos_Reals_iff)
paulson@62131
   221
  with t A Re_pos_in_ball[of z t] show ?thesis by (force simp add: complex_nonpos_Reals_iff)
eberlm@62049
   222
next
eberlm@62049
   223
  assume A: "Im z \<noteq> 0"
eberlm@62049
   224
  have "abs (Im z) - abs (Im t) \<le> abs (Im z - Im t)" by linarith
eberlm@62049
   225
  also have "\<dots> = abs (Im (z - t))" by simp
eberlm@62049
   226
  also have "\<dots> \<le> norm (z - t)" by (rule abs_Im_le_cmod)
paulson@62131
   227
  also from A t have "\<dots> \<le> abs (Im z) / 2" by (simp add: dist_complex_def)
eberlm@62049
   228
  finally have "abs (Im t) > 0" using A by simp
paulson@62131
   229
  thus ?thesis by (force simp add: complex_nonpos_Reals_iff)
eberlm@62049
   230
qed
eberlm@62049
   231
eberlm@62049
   232
eberlm@68624
   233
subsection \<open>The Euler form and the logarithmic Gamma function\<close>
eberlm@62049
   234
eberlm@62049
   235
text \<open>
eberlm@68624
   236
  We define the Gamma function by first defining its multiplicative inverse \<open>rGamma\<close>.
eberlm@68624
   237
  This is more convenient because \<open>rGamma\<close> is entire, which makes proofs of its
eberlm@62049
   238
  properties more convenient because one does not have to watch out for discontinuities.
eberlm@68624
   239
  (e.g. \<open>rGamma\<close> fulfils \<open>rGamma z = z * rGamma (z + 1)\<close> everywhere, whereas the \<open>\<Gamma>\<close> function
eberlm@68624
   240
  does not fulfil the analogous equation on the non-positive integers)
eberlm@68624
   241
eberlm@68624
   242
  We define the \<open>\<Gamma>\<close> function (resp.\ its reciprocale) in the Euler form. This form has the advantage
paulson@62131
   243
  that it is a relatively simple limit that converges everywhere. The limit at the poles is 0
eberlm@68624
   244
  (due to division by 0). The functional equation \<open>Gamma (z + 1) = z * Gamma z\<close> follows
eberlm@62049
   245
  immediately from the definition.
eberlm@62049
   246
\<close>
eberlm@62049
   247
eberlm@68624
   248
definition%important Gamma_series :: "('a :: {banach,real_normed_field}) \<Rightarrow> nat \<Rightarrow> 'a" where
eberlm@62049
   249
  "Gamma_series z n = fact n * exp (z * of_real (ln (of_nat n))) / pochhammer z (n+1)"
eberlm@62049
   250
eberlm@62049
   251
definition Gamma_series' :: "('a :: {banach,real_normed_field}) \<Rightarrow> nat \<Rightarrow> 'a" where
eberlm@62049
   252
  "Gamma_series' z n = fact (n - 1) * exp (z * of_real (ln (of_nat n))) / pochhammer z n"
eberlm@62049
   253
eberlm@62049
   254
definition rGamma_series :: "('a :: {banach,real_normed_field}) \<Rightarrow> nat \<Rightarrow> 'a" where
eberlm@62049
   255
  "rGamma_series z n = pochhammer z (n+1) / (fact n * exp (z * of_real (ln (of_nat n))))"
eberlm@62049
   256
eberlm@62049
   257
lemma Gamma_series_altdef: "Gamma_series z n = inverse (rGamma_series z n)"
eberlm@62049
   258
  and rGamma_series_altdef: "rGamma_series z n = inverse (Gamma_series z n)"
eberlm@62049
   259
  unfolding Gamma_series_def rGamma_series_def by simp_all
paulson@62131
   260
eberlm@62049
   261
lemma rGamma_series_minus_of_nat:
eberlm@62049
   262
  "eventually (\<lambda>n. rGamma_series (- of_nat k) n = 0) sequentially"
eberlm@62049
   263
  using eventually_ge_at_top[of k]
paulson@62131
   264
  by eventually_elim (auto simp: rGamma_series_def pochhammer_of_nat_eq_0_iff)
eberlm@62049
   265
eberlm@62049
   266
lemma Gamma_series_minus_of_nat:
eberlm@62049
   267
  "eventually (\<lambda>n. Gamma_series (- of_nat k) n = 0) sequentially"
eberlm@62049
   268
  using eventually_ge_at_top[of k]
paulson@62131
   269
  by eventually_elim (auto simp: Gamma_series_def pochhammer_of_nat_eq_0_iff)
eberlm@62049
   270
eberlm@62049
   271
lemma Gamma_series'_minus_of_nat:
eberlm@62049
   272
  "eventually (\<lambda>n. Gamma_series' (- of_nat k) n = 0) sequentially"
eberlm@62049
   273
  using eventually_gt_at_top[of k]
eberlm@62049
   274
  by eventually_elim (auto simp: Gamma_series'_def pochhammer_of_nat_eq_0_iff)
paulson@62131
   275
eberlm@62049
   276
lemma rGamma_series_nonpos_Ints_LIMSEQ: "z \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> rGamma_series z \<longlonglongrightarrow> 0"
eberlm@62049
   277
  by (elim nonpos_Ints_cases', hypsubst, subst tendsto_cong, rule rGamma_series_minus_of_nat, simp)
eberlm@62049
   278
eberlm@62049
   279
lemma Gamma_series_nonpos_Ints_LIMSEQ: "z \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Gamma_series z \<longlonglongrightarrow> 0"
eberlm@62049
   280
  by (elim nonpos_Ints_cases', hypsubst, subst tendsto_cong, rule Gamma_series_minus_of_nat, simp)
eberlm@62049
   281
eberlm@62049
   282
lemma Gamma_series'_nonpos_Ints_LIMSEQ: "z \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Gamma_series' z \<longlonglongrightarrow> 0"
eberlm@62049
   283
  by (elim nonpos_Ints_cases', hypsubst, subst tendsto_cong, rule Gamma_series'_minus_of_nat, simp)
paulson@62131
   284
paulson@62131
   285
lemma Gamma_series_Gamma_series':
eberlm@62049
   286
  assumes z: "z \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
   287
  shows   "(\<lambda>n. Gamma_series' z n / Gamma_series z n) \<longlonglongrightarrow> 1"
eberlm@62049
   288
proof (rule Lim_transform_eventually)
eberlm@62049
   289
  from eventually_gt_at_top[of "0::nat"]
eberlm@62049
   290
    show "eventually (\<lambda>n. z / of_nat n + 1 = Gamma_series' z n / Gamma_series z n) sequentially"
eberlm@62049
   291
  proof eventually_elim
eberlm@62049
   292
    fix n :: nat assume n: "n > 0"
eberlm@62049
   293
    from n z have "Gamma_series' z n / Gamma_series z n = (z + of_nat n) / of_nat n"
eberlm@62049
   294
      by (cases n, simp)
eberlm@62049
   295
         (auto simp add: Gamma_series_def Gamma_series'_def pochhammer_rec'
eberlm@62049
   296
               dest: pochhammer_eq_0_imp_nonpos_Int plus_of_nat_eq_0_imp)
eberlm@62049
   297
    also from n have "\<dots> = z / of_nat n + 1" by (simp add: divide_simps)
eberlm@62049
   298
    finally show "z / of_nat n + 1 = Gamma_series' z n / Gamma_series z n" ..
eberlm@62049
   299
  qed
eberlm@62049
   300
  have "(\<lambda>x. z / of_nat x) \<longlonglongrightarrow> 0"
eberlm@62049
   301
    by (rule tendsto_norm_zero_cancel)
paulson@62131
   302
       (insert tendsto_mult[OF tendsto_const[of "norm z"] lim_inverse_n],
eberlm@62049
   303
        simp add:  norm_divide inverse_eq_divide)
eberlm@62049
   304
  from tendsto_add[OF this tendsto_const[of 1]] show "(\<lambda>n. z / of_nat n + 1) \<longlonglongrightarrow> 1" by simp
eberlm@62049
   305
qed
eberlm@62049
   306
eberlm@62049
   307
text \<open>
eberlm@68624
   308
  We now show that the series that defines the \<open>\<Gamma>\<close> function in the Euler form converges
paulson@62131
   309
  and that the function defined by it is continuous on the complex halfspace with positive
eberlm@62049
   310
  real part.
paulson@62131
   311
paulson@62131
   312
  We do this by showing that the logarithm of the Euler series is continuous and converges
paulson@62131
   313
  locally uniformly, which means that the log-Gamma function defined by its limit is also
eberlm@62049
   314
  continuous.
paulson@62131
   315
paulson@62131
   316
  This will later allow us to lift holomorphicity and continuity from the log-Gamma
eberlm@62085
   317
  function to the inverse of the Gamma function, and from that to the Gamma function itself.
eberlm@62049
   318
\<close>
eberlm@62049
   319
eberlm@68624
   320
definition%important ln_Gamma_series :: "('a :: {banach,real_normed_field,ln}) \<Rightarrow> nat \<Rightarrow> 'a" where
eberlm@62049
   321
  "ln_Gamma_series z n = z * ln (of_nat n) - ln z - (\<Sum>k=1..n. ln (z / of_nat k + 1))"
eberlm@62049
   322
eberlm@68624
   323
definition%unimportant ln_Gamma_series' :: "('a :: {banach,real_normed_field,ln}) \<Rightarrow> nat \<Rightarrow> 'a" where
eberlm@62049
   324
  "ln_Gamma_series' z n =
eberlm@62049
   325
     - euler_mascheroni*z - ln z + (\<Sum>k=1..n. z / of_nat n - ln (z / of_nat k + 1))"
eberlm@62049
   326
eberlm@62049
   327
definition ln_Gamma :: "('a :: {banach,real_normed_field,ln}) \<Rightarrow> 'a" where
eberlm@62049
   328
  "ln_Gamma z = lim (ln_Gamma_series z)"
eberlm@62049
   329
eberlm@62049
   330
eberlm@62049
   331
text \<open>
paulson@62131
   332
  We now show that the log-Gamma series converges locally uniformly for all complex numbers except
eberlm@68624
   333
  the non-positive integers. We do this by proving that the series is locally Cauchy.
eberlm@62049
   334
\<close>
eberlm@62049
   335
eberlm@62049
   336
context
eberlm@62049
   337
begin
eberlm@62049
   338
eberlm@62049
   339
private lemma ln_Gamma_series_complex_converges_aux:
eberlm@62049
   340
  fixes z :: complex and k :: nat
eberlm@62049
   341
  assumes z: "z \<noteq> 0" and k: "of_nat k \<ge> 2*norm z" "k \<ge> 2"
eberlm@62049
   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"
eberlm@62049
   343
proof -
eberlm@62049
   344
  let ?k = "of_nat k :: complex" and ?z = "norm z"
eberlm@62049
   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)"
eberlm@62049
   346
    by (simp add: algebra_simps)
eberlm@62049
   347
  also have "norm ... \<le> ?z * norm (ln (1-1/?k) + 1/?k :: complex) + norm (ln (1+z/?k) - z/?k)"
eberlm@62049
   348
    by (subst norm_mult [symmetric], rule norm_triangle_ineq)
eberlm@62049
   349
  also have "norm (Ln (1 + -1/?k) - (-1/?k)) \<le> (norm (-1/?k))\<^sup>2 / (1 - norm(-1/?k))"
eberlm@62049
   350
    using k by (intro Ln_approx_linear) (simp add: norm_divide)
paulson@62131
   351
  hence "?z * norm (ln (1-1/?k) + 1/?k) \<le> ?z * ((norm (1/?k))^2 / (1 - norm (1/?k)))"
eberlm@62049
   352
    by (intro mult_left_mono) simp_all
eberlm@62049
   353
  also have "... \<le> (?z * (of_nat k / (of_nat k - 1))) / of_nat k^2" using k
eberlm@62049
   354
    by (simp add: field_simps power2_eq_square norm_divide)
eberlm@62049
   355
  also have "... \<le> (?z * 2) / of_nat k^2" using k
eberlm@62049
   356
    by (intro divide_right_mono mult_left_mono) (simp_all add: field_simps)
eberlm@62049
   357
  also have "norm (ln (1+z/?k) - z/?k) \<le> norm (z/?k)^2 / (1 - norm (z/?k))" using k
eberlm@62049
   358
    by (intro Ln_approx_linear) (simp add: norm_divide)
paulson@62131
   359
  hence "norm (ln (1+z/?k) - z/?k) \<le> ?z^2 / of_nat k^2 / (1 - ?z / of_nat k)"
eberlm@62049
   360
    by (simp add: field_simps norm_divide)
eberlm@62049
   361
  also have "... \<le> (?z^2 * (of_nat k / (of_nat k - ?z))) / of_nat k^2" using k
eberlm@62049
   362
    by (simp add: field_simps power2_eq_square)
eberlm@62049
   363
  also have "... \<le> (?z^2 * 2) / of_nat k^2" using k
eberlm@62049
   364
    by (intro divide_right_mono mult_left_mono) (simp_all add: field_simps)
eberlm@62049
   365
  also note add_divide_distrib [symmetric]
eberlm@62049
   366
  finally show ?thesis by (simp only: distrib_left mult.commute)
eberlm@62049
   367
qed
eberlm@62049
   368
eberlm@62049
   369
lemma ln_Gamma_series_complex_converges:
eberlm@62049
   370
  assumes z: "z \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
   371
  assumes d: "d > 0" "\<And>n. n \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> norm (z - of_int n) > d"
eberlm@62049
   372
  shows "uniformly_convergent_on (ball z d) (\<lambda>n z. ln_Gamma_series z n :: complex)"
eberlm@62049
   373
proof (intro Cauchy_uniformly_convergent uniformly_Cauchy_onI')
eberlm@62049
   374
  fix e :: real assume e: "e > 0"
haftmann@69260
   375
  define e'' where "e'' = (SUP t\<in>ball z d. norm t + norm t^2)"
wenzelm@63040
   376
  define e' where "e' = e / (2*e'')"
eberlm@62049
   377
  have "bounded ((\<lambda>t. norm t + norm t^2) ` cball z d)"
eberlm@62049
   378
    by (intro compact_imp_bounded compact_continuous_image) (auto intro!: continuous_intros)
eberlm@62049
   379
  hence "bounded ((\<lambda>t. norm t + norm t^2) ` ball z d)" by (rule bounded_subset) auto
eberlm@62049
   380
  hence bdd: "bdd_above ((\<lambda>t. norm t + norm t^2) ` ball z d)" by (rule bounded_imp_bdd_above)
eberlm@62049
   381
eberlm@62049
   382
  with z d(1) d(2)[of "-1"] have e''_pos: "e'' > 0" unfolding e''_def
eberlm@62049
   383
    by (subst less_cSUP_iff) (auto intro!: add_pos_nonneg bexI[of _ z])
eberlm@62049
   384
  have e'': "norm t + norm t^2 \<le> e''" if "t \<in> ball z d" for t unfolding e''_def using that
eberlm@62049
   385
    by (rule cSUP_upper[OF _ bdd])
paulson@62131
   386
  from e z e''_pos have e': "e' > 0" unfolding e'_def
eberlm@62049
   387
    by (intro divide_pos_pos mult_pos_pos add_pos_pos) (simp_all add: field_simps)
eberlm@62049
   388
eberlm@62049
   389
  have "summable (\<lambda>k. inverse ((real_of_nat k)^2))"
eberlm@62049
   390
    by (rule inverse_power_summable) simp
eberlm@62049
   391
  from summable_partial_sum_bound[OF this e'] guess M . note M = this
eberlm@62049
   392
wenzelm@63040
   393
  define N where "N = max 2 (max (nat \<lceil>2 * (norm z + d)\<rceil>) M)"
eberlm@62049
   394
  {
paulson@62131
   395
    from d have "\<lceil>2 * (cmod z + d)\<rceil> \<ge> \<lceil>0::real\<rceil>"
eberlm@62049
   396
      by (intro ceiling_mono mult_nonneg_nonneg add_nonneg_nonneg) simp_all
eberlm@62049
   397
    hence "2 * (norm z + d) \<le> of_nat (nat \<lceil>2 * (norm z + d)\<rceil>)" unfolding N_def
eberlm@62049
   398
      by (simp_all add: le_of_int_ceiling)
paulson@62131
   399
    also have "... \<le> of_nat N" unfolding N_def
eberlm@62049
   400
      by (subst of_nat_le_iff) (rule max.coboundedI2, rule max.cobounded1)
eberlm@62049
   401
    finally have "of_nat N \<ge> 2 * (norm z + d)" .
eberlm@62049
   402
    moreover have "N \<ge> 2" "N \<ge> M" unfolding N_def by simp_all
eberlm@62049
   403
    moreover have "(\<Sum>k=m..n. 1/(of_nat k)\<^sup>2) < e'" if "m \<ge> N" for m n
wenzelm@62072
   404
      using M[OF order.trans[OF \<open>N \<ge> M\<close> that]] unfolding real_norm_def
nipkow@64267
   405
      by (subst (asm) abs_of_nonneg) (auto intro: sum_nonneg simp: divide_simps)
eberlm@62049
   406
    moreover note calculation
eberlm@62049
   407
  } note N = this
eberlm@62049
   408
eberlm@62049
   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"
eberlm@62049
   410
    unfolding dist_complex_def
eberlm@62049
   411
  proof (intro exI[of _ N] ballI allI impI)
eberlm@62049
   412
    fix t m n assume t: "t \<in> ball z d" and mn: "m \<ge> N" "n > m"
eberlm@62049
   413
    from d(2)[of 0] t have "0 < dist z 0 - dist z t" by (simp add: field_simps dist_complex_def)
eberlm@62049
   414
    also have "dist z 0 - dist z t \<le> dist 0 t" using dist_triangle[of 0 z t]
eberlm@62049
   415
      by (simp add: dist_commute)
eberlm@62049
   416
    finally have t_nz: "t \<noteq> 0" by auto
eberlm@62049
   417
eberlm@62049
   418
    have "norm t \<le> norm z + norm (t - z)" by (rule norm_triangle_sub)
eberlm@62049
   419
    also from t have "norm (t - z) < d" by (simp add: dist_complex_def norm_minus_commute)
eberlm@62049
   420
    also have "2 * (norm z + d) \<le> of_nat N" by (rule N)
eberlm@62049
   421
    also have "N \<le> m" by (rule mn)
eberlm@62049
   422
    finally have norm_t: "2 * norm t < of_nat m" by simp
eberlm@62049
   423
paulson@62131
   424
    have "ln_Gamma_series t m - ln_Gamma_series t n =
eberlm@62049
   425
              (-(t * Ln (of_nat n)) - (-(t * Ln (of_nat m)))) +
eberlm@62049
   426
              ((\<Sum>k=1..n. Ln (t / of_nat k + 1)) - (\<Sum>k=1..m. Ln (t / of_nat k + 1)))"
eberlm@62049
   427
      by (simp add: ln_Gamma_series_def algebra_simps)
eberlm@62049
   428
    also have "(\<Sum>k=1..n. Ln (t / of_nat k + 1)) - (\<Sum>k=1..m. Ln (t / of_nat k + 1)) =
eberlm@62049
   429
                 (\<Sum>k\<in>{1..n}-{1..m}. Ln (t / of_nat k + 1))" using mn
nipkow@64267
   430
      by (simp add: sum_diff)
eberlm@62049
   431
    also from mn have "{1..n}-{1..m} = {Suc m..n}" by fastforce
paulson@62131
   432
    also have "-(t * Ln (of_nat n)) - (-(t * Ln (of_nat m))) =
eberlm@62049
   433
                   (\<Sum>k = Suc m..n. t * Ln (of_nat (k - 1)) - t * Ln (of_nat k))" using mn
nipkow@64267
   434
      by (subst sum_telescope'' [symmetric]) simp_all
eberlm@62049
   435
    also have "... = (\<Sum>k = Suc m..n. t * Ln (of_nat (k - 1) / of_nat k))" using mn N
nipkow@64267
   436
      by (intro sum_cong_Suc)
eberlm@62049
   437
         (simp_all del: of_nat_Suc add: field_simps Ln_of_nat Ln_of_nat_over_of_nat)
eberlm@62049
   438
    also have "of_nat (k - 1) / of_nat k = 1 - 1 / (of_nat k :: complex)" if "k \<in> {Suc m..n}" for k
eberlm@62049
   439
      using that of_nat_eq_0_iff[of "Suc i" for i] by (cases k) (simp_all add: divide_simps)
eberlm@62049
   440
    hence "(\<Sum>k = Suc m..n. t * Ln (of_nat (k - 1) / of_nat k)) =
eberlm@62049
   441
             (\<Sum>k = Suc m..n. t * Ln (1 - 1 / of_nat k))" using mn N
nipkow@64267
   442
      by (intro sum.cong) simp_all
nipkow@64267
   443
    also note sum.distrib [symmetric]
paulson@62131
   444
    also have "norm (\<Sum>k=Suc m..n. t * Ln (1 - 1/of_nat k) + Ln (t/of_nat k + 1)) \<le>
eberlm@62049
   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
nipkow@64267
   446
      by (intro order.trans[OF norm_sum sum_mono[OF ln_Gamma_series_complex_converges_aux]]) simp_all
eberlm@62049
   447
    also have "... \<le> 2 * (norm t + norm t^2) * (\<Sum>k=Suc m..n. 1 / (of_nat k)\<^sup>2)"
nipkow@64267
   448
      by (simp add: sum_distrib_left)
eberlm@62049
   449
    also have "... < 2 * (norm t + norm t^2) * e'" using mn z t_nz
eberlm@62049
   450
      by (intro mult_strict_left_mono N mult_pos_pos add_pos_pos) simp_all
paulson@62131
   451
    also from e''_pos have "... = e * ((cmod t + (cmod t)\<^sup>2) / e'')"
eberlm@62049
   452
      by (simp add: e'_def field_simps power2_eq_square)
paulson@62131
   453
    also from e''[OF t] e''_pos e
eberlm@62049
   454
      have "\<dots> \<le> e * 1" by (intro mult_left_mono) (simp_all add: field_simps)
eberlm@62049
   455
    finally show "norm (ln_Gamma_series t m - ln_Gamma_series t n) < e" by simp
eberlm@62049
   456
  qed
eberlm@62049
   457
qed
eberlm@62049
   458
eberlm@62049
   459
end
eberlm@62049
   460
eberlm@62049
   461
lemma ln_Gamma_series_complex_converges':
eberlm@62049
   462
  assumes z: "(z :: complex) \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
   463
  shows "\<exists>d>0. uniformly_convergent_on (ball z d) (\<lambda>n z. ln_Gamma_series z n)"
eberlm@62049
   464
proof -
wenzelm@63040
   465
  define d' where "d' = Re z"
wenzelm@63040
   466
  define d where "d = (if d' > 0 then d' / 2 else norm (z - of_int (round d')) / 2)"
eberlm@62049
   467
  have "of_int (round d') \<in> \<int>\<^sub>\<le>\<^sub>0" if "d' \<le> 0" using that
eberlm@62049
   468
    by (intro nonpos_Ints_of_int) (simp_all add: round_def)
eberlm@62049
   469
  with assms have d_pos: "d > 0" unfolding d_def by (force simp: not_less)
eberlm@62049
   470
eberlm@62049
   471
  have "d < cmod (z - of_int n)" if "n \<in> \<int>\<^sub>\<le>\<^sub>0" for n
eberlm@62049
   472
  proof (cases "Re z > 0")
eberlm@62049
   473
    case True
eberlm@62049
   474
    from nonpos_Ints_nonpos[OF that] have n: "n \<le> 0" by simp
eberlm@62049
   475
    from True have "d = Re z/2" by (simp add: d_def d'_def)
eberlm@62049
   476
    also from n True have "\<dots> < Re (z - of_int n)" by simp
eberlm@62049
   477
    also have "\<dots> \<le> norm (z - of_int n)" by (rule complex_Re_le_cmod)
eberlm@62049
   478
    finally show ?thesis .
eberlm@62049
   479
  next
eberlm@62049
   480
    case False
paulson@62131
   481
    with assms nonpos_Ints_of_int[of "round (Re z)"]
eberlm@62049
   482
      have "z \<noteq> of_int (round d')" by (auto simp: not_less)
eberlm@62049
   483
    with False have "d < norm (z - of_int (round d'))" by (simp add: d_def d'_def)
eberlm@62049
   484
    also have "\<dots> \<le> norm (z - of_int n)" unfolding d'_def by (rule round_Re_minimises_norm)
eberlm@62049
   485
    finally show ?thesis .
eberlm@62049
   486
  qed
eberlm@62049
   487
  hence conv: "uniformly_convergent_on (ball z d) (\<lambda>n z. ln_Gamma_series z n)"
eberlm@62049
   488
    by (intro ln_Gamma_series_complex_converges d_pos z) simp_all
eberlm@62049
   489
  from d_pos conv show ?thesis by blast
eberlm@62049
   490
qed
eberlm@62049
   491
eberlm@62049
   492
lemma ln_Gamma_series_complex_converges'': "(z :: complex) \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> convergent (ln_Gamma_series z)"
eberlm@62049
   493
  by (drule ln_Gamma_series_complex_converges') (auto intro: uniformly_convergent_imp_convergent)
eberlm@62049
   494
eberlm@68624
   495
theorem ln_Gamma_complex_LIMSEQ: "(z :: complex) \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> ln_Gamma_series z \<longlonglongrightarrow> ln_Gamma z"
eberlm@62049
   496
  using ln_Gamma_series_complex_converges'' by (simp add: convergent_LIMSEQ_iff ln_Gamma_def)
eberlm@62049
   497
eberlm@62049
   498
lemma exp_ln_Gamma_series_complex:
eberlm@62049
   499
  assumes "n > 0" "z \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
   500
  shows   "exp (ln_Gamma_series z n :: complex) = Gamma_series z n"
eberlm@62049
   501
proof -
haftmann@63417
   502
  from assms obtain m where m: "n = Suc m" by (cases n) blast
eberlm@62049
   503
  from assms have "z \<noteq> 0" by (intro notI) auto
paulson@62131
   504
  with assms have "exp (ln_Gamma_series z n) =
eberlm@62049
   505
          (of_nat n) powr z / (z * (\<Prod>k=1..n. exp (Ln (z / of_nat k + 1))))"
nipkow@64267
   506
    unfolding ln_Gamma_series_def powr_def by (simp add: exp_diff exp_sum)
eberlm@62049
   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)"
nipkow@64272
   508
    by (intro prod.cong[OF refl], subst exp_Ln) (auto simp: field_simps plus_of_nat_eq_0_imp)
haftmann@63417
   509
  also have "... = (\<Prod>k=1..n. z + k) / fact n"
nipkow@64272
   510
    by (simp add: fact_prod)
nipkow@64272
   511
    (subst prod_dividef [symmetric], simp_all add: field_simps)
haftmann@63417
   512
  also from m have "z * ... = (\<Prod>k=0..n. z + k) / fact n"
nipkow@64272
   513
    by (simp add: prod.atLeast0_atMost_Suc_shift prod.atLeast_Suc_atMost_Suc_shift)
haftmann@63417
   514
  also have "(\<Prod>k=0..n. z + k) = pochhammer z (Suc n)"
nipkow@64272
   515
    unfolding pochhammer_prod
nipkow@64272
   516
    by (simp add: prod.atLeast0_atMost_Suc atLeastLessThanSuc_atLeastAtMost)
eberlm@62049
   517
  also have "of_nat n powr z / (pochhammer z (Suc n) / fact n) = Gamma_series z n"
haftmann@66936
   518
    unfolding Gamma_series_def using assms by (simp add: divide_simps powr_def)
eberlm@62049
   519
  finally show ?thesis .
eberlm@62049
   520
qed
eberlm@62049
   521
eberlm@62049
   522
eberlm@62049
   523
lemma ln_Gamma_series'_aux:
eberlm@62049
   524
  assumes "(z::complex) \<notin> \<int>\<^sub>\<le>\<^sub>0"
paulson@62131
   525
  shows   "(\<lambda>k. z / of_nat (Suc k) - ln (1 + z / of_nat (Suc k))) sums
eberlm@62049
   526
              (ln_Gamma z + euler_mascheroni * z + ln z)" (is "?f sums ?s")
eberlm@62049
   527
unfolding sums_def
eberlm@62049
   528
proof (rule Lim_transform)
eberlm@62049
   529
  show "(\<lambda>n. ln_Gamma_series z n + of_real (harm n - ln (of_nat n)) * z + ln z) \<longlonglongrightarrow> ?s"
paulson@62131
   530
    (is "?g \<longlonglongrightarrow> _")
eberlm@62049
   531
    by (intro tendsto_intros ln_Gamma_complex_LIMSEQ euler_mascheroni_LIMSEQ_of_real assms)
eberlm@62049
   532
eberlm@62049
   533
  have A: "eventually (\<lambda>n. (\<Sum>k<n. ?f k) - ?g n = 0) sequentially"
eberlm@62049
   534
    using eventually_gt_at_top[of "0::nat"]
eberlm@62049
   535
  proof eventually_elim
eberlm@62049
   536
    fix n :: nat assume n: "n > 0"
eberlm@62049
   537
    have "(\<Sum>k<n. ?f k) = (\<Sum>k=1..n. z / of_nat k - ln (1 + z / of_nat k))"
nipkow@64267
   538
      by (subst atLeast0LessThan [symmetric], subst sum_shift_bounds_Suc_ivl [symmetric],
eberlm@62049
   539
          subst atLeastLessThanSuc_atLeastAtMost) simp_all
eberlm@62049
   540
    also have "\<dots> = z * of_real (harm n) - (\<Sum>k=1..n. ln (1 + z / of_nat k))"
nipkow@64267
   541
      by (simp add: harm_def sum_subtractf sum_distrib_left divide_inverse)
eberlm@62049
   542
    also from n have "\<dots> - ?g n = 0"
nipkow@64267
   543
      by (simp add: ln_Gamma_series_def sum_subtractf algebra_simps Ln_of_nat)
eberlm@62049
   544
    finally show "(\<Sum>k<n. ?f k) - ?g n = 0" .
eberlm@62049
   545
  qed
eberlm@62049
   546
  show "(\<lambda>n. (\<Sum>k<n. ?f k) - ?g n) \<longlonglongrightarrow> 0" by (subst tendsto_cong[OF A]) simp_all
eberlm@62049
   547
qed
eberlm@62049
   548
eberlm@62049
   549
eberlm@62049
   550
lemma uniformly_summable_deriv_ln_Gamma:
eberlm@62049
   551
  assumes z: "(z :: 'a :: {real_normed_field,banach}) \<noteq> 0" and d: "d > 0" "d \<le> norm z/2"
eberlm@62049
   552
  shows "uniformly_convergent_on (ball z d)
eberlm@62049
   553
            (\<lambda>k z. \<Sum>i<k. inverse (of_nat (Suc i)) - inverse (z + of_nat (Suc i)))"
eberlm@62049
   554
           (is "uniformly_convergent_on _ (\<lambda>k z. \<Sum>i<k. ?f i z)")
nipkow@69529
   555
proof (rule Weierstrass_m_test'_ev)
eberlm@62049
   556
  {
eberlm@62049
   557
    fix t assume t: "t \<in> ball z d"
eberlm@62049
   558
    have "norm z = norm (t + (z - t))" by simp
eberlm@62049
   559
    have "norm (t + (z - t)) \<le> norm t + norm (z - t)" by (rule norm_triangle_ineq)
eberlm@62049
   560
    also from t d have "norm (z - t) < norm z / 2" by (simp add: dist_norm)
eberlm@62049
   561
    finally have A: "norm t > norm z / 2" using z by (simp add: field_simps)
paulson@62131
   562
eberlm@62049
   563
    have "norm t = norm (z + (t - z))" by simp
eberlm@62049
   564
    also have "\<dots> \<le> norm z + norm (t - z)" by (rule norm_triangle_ineq)
eberlm@62049
   565
    also from t d have "norm (t - z) \<le> norm z / 2" by (simp add: dist_norm norm_minus_commute)
eberlm@62049
   566
    also from z have "\<dots> < norm z" by simp
eberlm@62049
   567
    finally have B: "norm t < 2 * norm z" by simp
eberlm@62049
   568
    note A B
eberlm@62049
   569
  } note ball = this
eberlm@62049
   570
eberlm@62049
   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"
eberlm@62049
   572
    using eventually_gt_at_top apply eventually_elim
eberlm@62049
   573
  proof safe
eberlm@62049
   574
    fix t :: 'a assume t: "t \<in> ball z d"
eberlm@62049
   575
    from z ball[OF t] have t_nz: "t \<noteq> 0" by auto
eberlm@62049
   576
    fix n :: nat assume n: "n > nat \<lceil>4 * norm z\<rceil>"
eberlm@62049
   577
    from ball[OF t] t_nz have "4 * norm z > 2 * norm t" by simp
eberlm@62049
   578
    also from n have "\<dots>  < of_nat n" by linarith
eberlm@62049
   579
    finally have n: "of_nat n > 2 * norm t" .
eberlm@62049
   580
    hence "of_nat n > norm t" by simp
eberlm@62049
   581
    hence t': "t \<noteq> -of_nat (Suc n)" by (intro notI) (simp del: of_nat_Suc)
eberlm@62049
   582
eberlm@62049
   583
    with t_nz have "?f n t = 1 / (of_nat (Suc n) * (1 + of_nat (Suc n)/t))"
eberlm@62049
   584
      by (simp add: divide_simps eq_neg_iff_add_eq_0 del: of_nat_Suc)
eberlm@62049
   585
    also have "norm \<dots> = inverse (of_nat (Suc n)) * inverse (norm (of_nat (Suc n)/t + 1))"
eberlm@62049
   586
      by (simp add: norm_divide norm_mult divide_simps add_ac del: of_nat_Suc)
eberlm@62049
   587
    also {
eberlm@62049
   588
      from z t_nz ball[OF t] have "of_nat (Suc n) / (4 * norm z) \<le> of_nat (Suc n) / (2 * norm t)"
eberlm@62049
   589
        by (intro divide_left_mono mult_pos_pos) simp_all
paulson@62131
   590
      also have "\<dots> < norm (of_nat (Suc n) / t) - norm (1 :: 'a)"
eberlm@62049
   591
        using t_nz n by (simp add: field_simps norm_divide del: of_nat_Suc)
eberlm@62049
   592
      also have "\<dots> \<le> norm (of_nat (Suc n)/t + 1)" by (rule norm_diff_ineq)
eberlm@62049
   593
      finally have "inverse (norm (of_nat (Suc n)/t + 1)) \<le> 4 * norm z / of_nat (Suc n)"
eberlm@62049
   594
        using z by (simp add: divide_simps norm_divide mult_ac del: of_nat_Suc)
eberlm@62049
   595
    }
eberlm@62049
   596
    also have "inverse (real_of_nat (Suc n)) * (4 * norm z / real_of_nat (Suc n)) =
paulson@62131
   597
                 4 * norm z * inverse (of_nat (Suc n)^2)"
eberlm@62049
   598
                 by (simp add: divide_simps power2_eq_square del: of_nat_Suc)
paulson@62131
   599
    finally show "norm (?f n t) \<le> 4 * norm z * inverse (of_nat (Suc n)^2)"
eberlm@62049
   600
      by (simp del: of_nat_Suc)
eberlm@62049
   601
  qed
eberlm@62049
   602
next
eberlm@62049
   603
  show "summable (\<lambda>n. 4 * norm z * inverse ((of_nat (Suc n))^2))"
eberlm@62049
   604
    by (subst summable_Suc_iff) (simp add: summable_mult inverse_power_summable)
eberlm@62049
   605
qed
eberlm@62049
   606
eberlm@68624
   607
eberlm@68624
   608
subsection \<open>The Polygamma functions\<close>
eberlm@68624
   609
eberlm@62049
   610
lemma summable_deriv_ln_Gamma:
paulson@62131
   611
  "z \<noteq> (0 :: 'a :: {real_normed_field,banach}) \<Longrightarrow>
eberlm@62049
   612
     summable (\<lambda>n. inverse (of_nat (Suc n)) - inverse (z + of_nat (Suc n)))"
eberlm@62049
   613
  unfolding summable_iff_convergent
paulson@62131
   614
  by (rule uniformly_convergent_imp_convergent,
eberlm@62049
   615
      rule uniformly_summable_deriv_ln_Gamma[of z "norm z/2"]) simp_all
eberlm@62049
   616
eberlm@68624
   617
definition%important Polygamma :: "nat \<Rightarrow> ('a :: {real_normed_field,banach}) \<Rightarrow> 'a" where
paulson@62131
   618
  "Polygamma n z = (if n = 0 then
paulson@62131
   619
      (\<Sum>k. inverse (of_nat (Suc k)) - inverse (z + of_nat k)) - euler_mascheroni else
eberlm@62049
   620
      (-1)^Suc n * fact n * (\<Sum>k. inverse ((z + of_nat k)^Suc n)))"
eberlm@62049
   621
eberlm@68624
   622
abbreviation%important Digamma :: "('a :: {real_normed_field,banach}) \<Rightarrow> 'a" where
eberlm@62049
   623
  "Digamma \<equiv> Polygamma 0"
eberlm@62049
   624
paulson@62131
   625
lemma Digamma_def:
eberlm@62049
   626
  "Digamma z = (\<Sum>k. inverse (of_nat (Suc k)) - inverse (z + of_nat k)) - euler_mascheroni"
eberlm@62049
   627
  by (simp add: Polygamma_def)
eberlm@62049
   628
eberlm@62049
   629
paulson@62131
   630
lemma summable_Digamma:
eberlm@62049
   631
  assumes "(z :: 'a :: {real_normed_field,banach}) \<noteq> 0"
eberlm@62049
   632
  shows   "summable (\<lambda>n. inverse (of_nat (Suc n)) - inverse (z + of_nat n))"
eberlm@62049
   633
proof -
paulson@62131
   634
  have sums: "(\<lambda>n. inverse (z + of_nat (Suc n)) - inverse (z + of_nat n)) sums
eberlm@62049
   635
                       (0 - inverse (z + of_nat 0))"
eberlm@62049
   636
    by (intro telescope_sums filterlim_compose[OF tendsto_inverse_0]
eberlm@62049
   637
              tendsto_add_filterlim_at_infinity[OF tendsto_const] tendsto_of_nat)
eberlm@62049
   638
  from summable_add[OF summable_deriv_ln_Gamma[OF assms] sums_summable[OF sums]]
eberlm@62049
   639
    show "summable (\<lambda>n. inverse (of_nat (Suc n)) - inverse (z + of_nat n))" by simp
eberlm@62049
   640
qed
eberlm@62049
   641
eberlm@62049
   642
lemma summable_offset:
eberlm@62049
   643
  assumes "summable (\<lambda>n. f (n + k) :: 'a :: real_normed_vector)"
eberlm@62049
   644
  shows   "summable f"
eberlm@62049
   645
proof -
eberlm@62049
   646
  from assms have "convergent (\<lambda>m. \<Sum>n<m. f (n + k))" by (simp add: summable_iff_convergent)
paulson@62131
   647
  hence "convergent (\<lambda>m. (\<Sum>n<k. f n) + (\<Sum>n<m. f (n + k)))"
eberlm@62049
   648
    by (intro convergent_add convergent_const)
eberlm@62049
   649
  also have "(\<lambda>m. (\<Sum>n<k. f n) + (\<Sum>n<m. f (n + k))) = (\<lambda>m. \<Sum>n<m+k. f n)"
eberlm@62049
   650
  proof
eberlm@62049
   651
    fix m :: nat
eberlm@62049
   652
    have "{..<m+k} = {..<k} \<union> {k..<m+k}" by auto
paulson@62131
   653
    also have "(\<Sum>n\<in>\<dots>. f n) = (\<Sum>n<k. f n) + (\<Sum>n=k..<m+k. f n)"
nipkow@64267
   654
      by (rule sum.union_disjoint) auto
eberlm@62049
   655
    also have "(\<Sum>n=k..<m+k. f n) = (\<Sum>n=0..<m+k-k. f (n + k))"
haftmann@66936
   656
      using sum_shift_bounds_nat_ivl [of f 0 k m] by simp
eberlm@62049
   657
    finally show "(\<Sum>n<k. f n) + (\<Sum>n<m. f (n + k)) = (\<Sum>n<m+k. f n)" by (simp add: atLeast0LessThan)
eberlm@62049
   658
  qed
nipkow@64267
   659
  finally have "(\<lambda>a. sum f {..<a}) \<longlonglongrightarrow> lim (\<lambda>m. sum f {..<m + k})"
eberlm@62049
   660
    by (auto simp: convergent_LIMSEQ_iff dest: LIMSEQ_offset)
eberlm@62049
   661
  thus ?thesis by (auto simp: summable_iff_convergent convergent_def)
eberlm@62049
   662
qed
eberlm@62049
   663
eberlm@62049
   664
lemma Polygamma_converges:
eberlm@62049
   665
  fixes z :: "'a :: {real_normed_field,banach}"
eberlm@62049
   666
  assumes z: "z \<noteq> 0" and n: "n \<ge> 2"
eberlm@62049
   667
  shows "uniformly_convergent_on (ball z d) (\<lambda>k z. \<Sum>i<k. inverse ((z + of_nat i)^n))"
nipkow@69529
   668
proof (rule Weierstrass_m_test'_ev)
wenzelm@63040
   669
  define e where "e = (1 + d / norm z)"
wenzelm@63040
   670
  define m where "m = nat \<lceil>norm z * e\<rceil>"
eberlm@62049
   671
  {
eberlm@62049
   672
    fix t assume t: "t \<in> ball z d"
eberlm@62049
   673
    have "norm t = norm (z + (t - z))" by simp
eberlm@62049
   674
    also have "\<dots> \<le> norm z + norm (t - z)" by (rule norm_triangle_ineq)
eberlm@62049
   675
    also from t have "norm (t - z) < d" by (simp add: dist_norm norm_minus_commute)
eberlm@62049
   676
    finally have "norm t < norm z * e" using z by (simp add: divide_simps e_def)
eberlm@62049
   677
  } note ball = this
paulson@62131
   678
paulson@62131
   679
  show "eventually (\<lambda>k. \<forall>t\<in>ball z d. norm (inverse ((t + of_nat k)^n)) \<le>
eberlm@62049
   680
            inverse (of_nat (k - m)^n)) sequentially"
eberlm@62049
   681
    using eventually_gt_at_top[of m] apply eventually_elim
eberlm@62049
   682
  proof (intro ballI)
eberlm@62049
   683
    fix k :: nat and t :: 'a assume k: "k > m" and t: "t \<in> ball z d"
eberlm@62049
   684
    from k have "real_of_nat (k - m) = of_nat k - of_nat m" by (simp add: of_nat_diff)
paulson@62131
   685
    also have "\<dots> \<le> norm (of_nat k :: 'a) - norm z * e"
eberlm@62049
   686
      unfolding m_def by (subst norm_of_nat) linarith
eberlm@62049
   687
    also from ball[OF t] have "\<dots> \<le> norm (of_nat k :: 'a) - norm t" by simp
eberlm@62049
   688
    also have "\<dots> \<le> norm (of_nat k + t)" by (rule norm_diff_ineq)
eberlm@62049
   689
    finally have "inverse ((norm (t + of_nat k))^n) \<le> inverse (real_of_nat (k - m)^n)" using k n
eberlm@62049
   690
      by (intro le_imp_inverse_le power_mono) (simp_all add: add_ac del: of_nat_Suc)
eberlm@62049
   691
    thus "norm (inverse ((t + of_nat k)^n)) \<le> inverse (of_nat (k - m)^n)"
eberlm@62049
   692
      by (simp add: norm_inverse norm_power power_inverse)
eberlm@62049
   693
  qed
paulson@62131
   694
paulson@62131
   695
  have "summable (\<lambda>k. inverse ((real_of_nat k)^n))"
eberlm@62049
   696
    using inverse_power_summable[of n] n by simp
eberlm@62049
   697
  hence "summable (\<lambda>k. inverse ((real_of_nat (k + m - m))^n))" by simp
eberlm@62049
   698
  thus "summable (\<lambda>k. inverse ((real_of_nat (k - m))^n))" by (rule summable_offset)
eberlm@62049
   699
qed
eberlm@62049
   700
eberlm@62049
   701
lemma Polygamma_converges':
eberlm@62049
   702
  fixes z :: "'a :: {real_normed_field,banach}"
eberlm@62049
   703
  assumes z: "z \<noteq> 0" and n: "n \<ge> 2"
eberlm@62049
   704
  shows "summable (\<lambda>k. inverse ((z + of_nat k)^n))"
eberlm@62049
   705
  using uniformly_convergent_imp_convergent[OF Polygamma_converges[OF assms, of 1], of z]
eberlm@62049
   706
  by (simp add: summable_iff_convergent)
eberlm@62049
   707
eberlm@68624
   708
theorem Digamma_LIMSEQ:
eberlm@63721
   709
  fixes z :: "'a :: {banach,real_normed_field}"
eberlm@63721
   710
  assumes z: "z \<noteq> 0"
eberlm@63721
   711
  shows   "(\<lambda>m. of_real (ln (real m)) - (\<Sum>n<m. inverse (z + of_nat n))) \<longlonglongrightarrow> Digamma z"
eberlm@63721
   712
proof -
eberlm@63721
   713
  have "(\<lambda>n. of_real (ln (real n / (real (Suc n))))) \<longlonglongrightarrow> (of_real (ln 1) :: 'a)"
eberlm@63721
   714
    by (intro tendsto_intros LIMSEQ_n_over_Suc_n) simp_all
eberlm@63721
   715
  hence "(\<lambda>n. of_real (ln (real n / (real n + 1)))) \<longlonglongrightarrow> (0 :: 'a)" by (simp add: add_ac)
eberlm@63721
   716
  hence lim: "(\<lambda>n. of_real (ln (real n)) - of_real (ln (real n + 1))) \<longlonglongrightarrow> (0::'a)"
eberlm@63721
   717
  proof (rule Lim_transform_eventually [rotated])
lp15@66512
   718
    show "eventually (\<lambda>n. of_real (ln (real n / (real n + 1))) =
eberlm@63721
   719
            of_real (ln (real n)) - (of_real (ln (real n + 1)) :: 'a)) at_top"
eberlm@63721
   720
      using eventually_gt_at_top[of "0::nat"] by eventually_elim (simp add: ln_div)
eberlm@63721
   721
  qed
eberlm@63721
   722
eberlm@63721
   723
  from summable_Digamma[OF z]
lp15@66512
   724
    have "(\<lambda>n. inverse (of_nat (n+1)) - inverse (z + of_nat n))
eberlm@63721
   725
              sums (Digamma z + euler_mascheroni)"
eberlm@63721
   726
    by (simp add: Digamma_def summable_sums)
eberlm@63721
   727
  from sums_diff[OF this euler_mascheroni_sum]
eberlm@63721
   728
    have "(\<lambda>n. of_real (ln (real (Suc n) + 1)) - of_real (ln (real n + 1)) - inverse (z + of_nat n))
eberlm@63721
   729
            sums Digamma z" by (simp add: add_ac)
lp15@66512
   730
  hence "(\<lambda>m. (\<Sum>n<m. of_real (ln (real (Suc n) + 1)) - of_real (ln (real n + 1))) -
eberlm@63721
   731
              (\<Sum>n<m. inverse (z + of_nat n))) \<longlonglongrightarrow> Digamma z"
nipkow@64267
   732
    by (simp add: sums_def sum_subtractf)
lp15@66512
   733
  also have "(\<lambda>m. (\<Sum>n<m. of_real (ln (real (Suc n) + 1)) - of_real (ln (real n + 1)))) =
eberlm@63721
   734
                 (\<lambda>m. of_real (ln (m + 1)) :: 'a)"
nipkow@64267
   735
    by (subst sum_lessThan_telescope) simp_all
eberlm@63721
   736
  finally show ?thesis by (rule Lim_transform) (insert lim, simp)
eberlm@63721
   737
qed
eberlm@63721
   738
eberlm@68624
   739
theorem Polygamma_LIMSEQ:
eberlm@68624
   740
  fixes z :: "'a :: {banach,real_normed_field}"
eberlm@68624
   741
  assumes "z \<noteq> 0" and "n > 0"
eberlm@68624
   742
  shows   "(\<lambda>k. inverse ((z + of_nat k)^Suc n)) sums ((-1) ^ Suc n * Polygamma n z / fact n)"
eberlm@68624
   743
  using Polygamma_converges'[OF assms(1), of "Suc n"] assms(2)
eberlm@68624
   744
  by (simp add: sums_iff Polygamma_def)
eberlm@68624
   745
eberlm@68624
   746
theorem has_field_derivative_ln_Gamma_complex [derivative_intros]:
paulson@62131
   747
  fixes z :: complex
paulson@62131
   748
  assumes z: "z \<notin> \<real>\<^sub>\<le>\<^sub>0"
eberlm@62049
   749
  shows   "(ln_Gamma has_field_derivative Digamma z) (at z)"
eberlm@62049
   750
proof -
paulson@62131
   751
  have not_nonpos_Int [simp]: "t \<notin> \<int>\<^sub>\<le>\<^sub>0" if "Re t > 0" for t
eberlm@62049
   752
    using that by (auto elim!: nonpos_Ints_cases')
paulson@62131
   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
paulson@62131
   754
     by blast+
eberlm@62049
   755
  let ?f' = "\<lambda>z k. inverse (of_nat (Suc k)) - inverse (z + of_nat (Suc k))"
eberlm@62049
   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"
wenzelm@63040
   757
  define d where "d = min (norm z/2) (if Im z = 0 then Re z / 2 else abs (Im z) / 2)"
paulson@62131
   758
  from z have d: "d > 0" "norm z/2 \<ge> d" by (auto simp add: complex_nonpos_Reals_iff d_def)
paulson@62131
   759
  have ball: "Im t = 0 \<longrightarrow> Re t > 0" if "dist z t < d" for t
paulson@62131
   760
    using no_nonpos_Real_in_ball[OF z, of t] that unfolding d_def by (force simp add: complex_nonpos_Reals_iff)
paulson@62131
   761
  have sums: "(\<lambda>n. inverse (z + of_nat (Suc n)) - inverse (z + of_nat n)) sums
eberlm@62049
   762
                       (0 - inverse (z + of_nat 0))"
eberlm@62049
   763
    by (intro telescope_sums filterlim_compose[OF tendsto_inverse_0]
eberlm@62049
   764
              tendsto_add_filterlim_at_infinity[OF tendsto_const] tendsto_of_nat)
eberlm@62049
   765
eberlm@62049
   766
  have "((\<lambda>z. \<Sum>n. ?f z n) has_field_derivative ?F' z) (at z)"
eberlm@62049
   767
    using d z ln_Gamma_series'_aux[OF z']
paulson@62131
   768
    apply (intro has_field_derivative_series'(2)[of "ball z d" _ _ z] uniformly_summable_deriv_ln_Gamma)
paulson@62131
   769
    apply (auto intro!: derivative_eq_intros add_pos_pos mult_pos_pos dest!: ball
paulson@62131
   770
             simp: field_simps sums_iff nonpos_Reals_divide_of_nat_iff
eberlm@62049
   771
             simp del: of_nat_Suc)
paulson@62131
   772
    apply (auto simp add: complex_nonpos_Reals_iff)
paulson@62131
   773
    done
paulson@62131
   774
  with z have "((\<lambda>z. (\<Sum>k. ?f z k) - euler_mascheroni * z - Ln z) has_field_derivative
paulson@62131
   775
                   ?F' z - euler_mascheroni - inverse z) (at z)"
eberlm@62049
   776
    by (force intro!: derivative_eq_intros simp: Digamma_def)
eberlm@62049
   777
  also have "?F' z - euler_mascheroni - inverse z = (?F' z + -inverse z) - euler_mascheroni" by simp
eberlm@62049
   778
  also from sums have "-inverse z = (\<Sum>n. inverse (z + of_nat (Suc n)) - inverse (z + of_nat n))"
eberlm@62049
   779
    by (simp add: sums_iff)
paulson@62131
   780
  also from sums summable_deriv_ln_Gamma[OF z'']
eberlm@62049
   781
    have "?F' z + \<dots> =  (\<Sum>n. inverse (of_nat (Suc n)) - inverse (z + of_nat n))"
eberlm@62049
   782
    by (subst suminf_add) (simp_all add: add_ac sums_iff)
eberlm@62049
   783
  also have "\<dots> - euler_mascheroni = Digamma z" by (simp add: Digamma_def)
paulson@62131
   784
  finally have "((\<lambda>z. (\<Sum>k. ?f z k) - euler_mascheroni * z - Ln z)
eberlm@62049
   785
                    has_field_derivative Digamma z) (at z)" .
eberlm@62049
   786
  moreover from eventually_nhds_ball[OF d(1), of z]
eberlm@62049
   787
    have "eventually (\<lambda>z. ln_Gamma z = (\<Sum>k. ?f z k) - euler_mascheroni * z - Ln z) (nhds z)"
eberlm@62049
   788
  proof eventually_elim
eberlm@62049
   789
    fix t assume "t \<in> ball z d"
eberlm@62049
   790
    hence "t \<notin> \<int>\<^sub>\<le>\<^sub>0" by (auto dest!: ball elim!: nonpos_Ints_cases)
paulson@62131
   791
    from ln_Gamma_series'_aux[OF this]
eberlm@62049
   792
      show "ln_Gamma t = (\<Sum>k. ?f t k) - euler_mascheroni * t - Ln t" by (simp add: sums_iff)
eberlm@62049
   793
  qed
eberlm@62049
   794
  ultimately show ?thesis by (subst DERIV_cong_ev[OF refl _ refl])
eberlm@62049
   795
qed
eberlm@62049
   796
eberlm@62049
   797
declare has_field_derivative_ln_Gamma_complex[THEN DERIV_chain2, derivative_intros]
eberlm@62049
   798
eberlm@62049
   799
eberlm@62049
   800
lemma Digamma_1 [simp]: "Digamma (1 :: 'a :: {real_normed_field,banach}) = - euler_mascheroni"
eberlm@62049
   801
  by (simp add: Digamma_def)
paulson@62131
   802
eberlm@62049
   803
lemma Digamma_plus1:
eberlm@62049
   804
  assumes "z \<noteq> 0"
eberlm@62049
   805
  shows   "Digamma (z+1) = Digamma z + 1/z"
eberlm@62049
   806
proof -
paulson@62131
   807
  have sums: "(\<lambda>k. inverse (z + of_nat k) - inverse (z + of_nat (Suc k)))
eberlm@62049
   808
                  sums (inverse (z + of_nat 0) - 0)"
eberlm@62049
   809
    by (intro telescope_sums'[OF filterlim_compose[OF tendsto_inverse_0]]
eberlm@62049
   810
              tendsto_add_filterlim_at_infinity[OF tendsto_const] tendsto_of_nat)
paulson@62131
   811
  have "Digamma (z+1) = (\<Sum>k. inverse (of_nat (Suc k)) - inverse (z + of_nat (Suc k))) -
eberlm@62049
   812
          euler_mascheroni" (is "_ = suminf ?f - _") by (simp add: Digamma_def add_ac)
eberlm@62049
   813
  also have "suminf ?f = (\<Sum>k. inverse (of_nat (Suc k)) - inverse (z + of_nat k)) +
eberlm@62049
   814
                         (\<Sum>k. inverse (z + of_nat k) - inverse (z + of_nat (Suc k)))"
eberlm@62049
   815
    using summable_Digamma[OF assms] sums by (subst suminf_add) (simp_all add: add_ac sums_iff)
eberlm@62049
   816
  also have "(\<Sum>k. inverse (z + of_nat k) - inverse (z + of_nat (Suc k))) = 1/z"
eberlm@62049
   817
    using sums by (simp add: sums_iff inverse_eq_divide)
eberlm@62049
   818
  finally show ?thesis by (simp add: Digamma_def[of z])
eberlm@62049
   819
qed
eberlm@62049
   820
eberlm@68624
   821
theorem Polygamma_plus1:
eberlm@62049
   822
  assumes "z \<noteq> 0"
eberlm@62049
   823
  shows   "Polygamma n (z + 1) = Polygamma n z + (-1)^n * fact n / (z ^ Suc n)"
eberlm@62049
   824
proof (cases "n = 0")
eberlm@62049
   825
  assume n: "n \<noteq> 0"
eberlm@62049
   826
  let ?f = "\<lambda>k. inverse ((z + of_nat k) ^ Suc n)"
eberlm@62049
   827
  have "Polygamma n (z + 1) = (-1) ^ Suc n * fact n * (\<Sum>k. ?f (k+1))"
eberlm@62049
   828
    using n by (simp add: Polygamma_def add_ac)
eberlm@62049
   829
  also have "(\<Sum>k. ?f (k+1)) + (\<Sum>k<1. ?f k) = (\<Sum>k. ?f k)"
eberlm@62049
   830
    using Polygamma_converges'[OF assms, of "Suc n"] n
eberlm@62049
   831
    by (subst suminf_split_initial_segment [symmetric]) simp_all
eberlm@62049
   832
  hence "(\<Sum>k. ?f (k+1)) = (\<Sum>k. ?f k) - inverse (z ^ Suc n)" by (simp add: algebra_simps)
paulson@62131
   833
  also have "(-1) ^ Suc n * fact n * ((\<Sum>k. ?f k) - inverse (z ^ Suc n)) =
eberlm@62049
   834
               Polygamma n z + (-1)^n * fact n / (z ^ Suc n)" using n
eberlm@62049
   835
    by (simp add: inverse_eq_divide algebra_simps Polygamma_def)
eberlm@62049
   836
  finally show ?thesis .
eberlm@62049
   837
qed (insert assms, simp add: Digamma_plus1 inverse_eq_divide)
eberlm@62049
   838
eberlm@68624
   839
theorem Digamma_of_nat:
eberlm@62049
   840
  "Digamma (of_nat (Suc n) :: 'a :: {real_normed_field,banach}) = harm n - euler_mascheroni"
eberlm@62049
   841
proof (induction n)
eberlm@62049
   842
  case (Suc n)
eberlm@62049
   843
  have "Digamma (of_nat (Suc (Suc n)) :: 'a) = Digamma (of_nat (Suc n) + 1)" by simp
paulson@62131
   844
  also have "\<dots> = Digamma (of_nat (Suc n)) + inverse (of_nat (Suc n))"
eberlm@62049
   845
    by (subst Digamma_plus1) (simp_all add: inverse_eq_divide del: of_nat_Suc)
eberlm@62049
   846
  also have "Digamma (of_nat (Suc n) :: 'a) = harm n - euler_mascheroni " by (rule Suc)
eberlm@62049
   847
  also have "\<dots> + inverse (of_nat (Suc n)) = harm (Suc n) - euler_mascheroni"
eberlm@62049
   848
    by (simp add: harm_Suc)
eberlm@62049
   849
  finally show ?case .
eberlm@62049
   850
qed (simp add: harm_def)
eberlm@62049
   851
eberlm@62049
   852
lemma Digamma_numeral: "Digamma (numeral n) = harm (pred_numeral n) - euler_mascheroni"
eberlm@62049
   853
  by (subst of_nat_numeral[symmetric], subst numeral_eq_Suc, subst Digamma_of_nat) (rule refl)
eberlm@62049
   854
eberlm@62049
   855
lemma Polygamma_of_real: "x \<noteq> 0 \<Longrightarrow> Polygamma n (of_real x) = of_real (Polygamma n x)"
eberlm@62049
   856
  unfolding Polygamma_def using summable_Digamma[of x] Polygamma_converges'[of x "Suc n"]
eberlm@62049
   857
  by (simp_all add: suminf_of_real)
eberlm@62049
   858
eberlm@62049
   859
lemma Polygamma_Real: "z \<in> \<real> \<Longrightarrow> z \<noteq> 0 \<Longrightarrow> Polygamma n z \<in> \<real>"
eberlm@62049
   860
  by (elim Reals_cases, hypsubst, subst Polygamma_of_real) simp_all
eberlm@62049
   861
eberlm@62049
   862
lemma Digamma_half_integer:
paulson@62131
   863
  "Digamma (of_nat n + 1/2 :: 'a :: {real_normed_field,banach}) =
eberlm@62049
   864
      (\<Sum>k<n. 2 / (of_nat (2*k+1))) - euler_mascheroni - of_real (2 * ln 2)"
eberlm@62049
   865
proof (induction n)
eberlm@62049
   866
  case 0
eberlm@62049
   867
  have "Digamma (1/2 :: 'a) = of_real (Digamma (1/2))" by (simp add: Polygamma_of_real [symmetric])
paulson@62131
   868
  also have "Digamma (1/2::real) =
paulson@62131
   869
               (\<Sum>k. inverse (of_nat (Suc k)) - inverse (of_nat k + 1/2)) - euler_mascheroni"
eberlm@62049
   870
    by (simp add: Digamma_def add_ac)
eberlm@62049
   871
  also have "(\<Sum>k. inverse (of_nat (Suc k) :: real) - inverse (of_nat k + 1/2)) =
eberlm@62049
   872
             (\<Sum>k. inverse (1/2) * (inverse (2 * of_nat (Suc k)) - inverse (2 * of_nat k + 1)))"
eberlm@62049
   873
    by (simp_all add: add_ac inverse_mult_distrib[symmetric] ring_distribs del: inverse_divide)
eberlm@62049
   874
  also have "\<dots> = - 2 * ln 2" using sums_minus[OF alternating_harmonic_series_sums']
eberlm@62049
   875
    by (subst suminf_mult) (simp_all add: algebra_simps sums_iff)
eberlm@62049
   876
  finally show ?case by simp
eberlm@62049
   877
next
eberlm@62049
   878
  case (Suc n)
eberlm@62049
   879
  have nz: "2 * of_nat n + (1:: 'a) \<noteq> 0"
eberlm@62049
   880
     using of_nat_neq_0[of "2*n"] by (simp only: of_nat_Suc) (simp add: add_ac)
eberlm@62049
   881
  hence nz': "of_nat n + (1/2::'a) \<noteq> 0" by (simp add: field_simps)
eberlm@62049
   882
  have "Digamma (of_nat (Suc n) + 1/2 :: 'a) = Digamma (of_nat n + 1/2 + 1)" by simp
lp15@67976
   883
  also from nz' have "\<dots> = Digamma (of_nat n + 1/2) + 1 / (of_nat n + 1/2)"
eberlm@62049
   884
    by (rule Digamma_plus1)
lp15@67976
   885
  also from nz nz' have "1 / (of_nat n + 1/2 :: 'a) = 2 / (2 * of_nat n + 1)"
eberlm@62049
   886
    by (subst divide_eq_eq) simp_all
eberlm@62049
   887
  also note Suc
eberlm@62049
   888
  finally show ?case by (simp add: add_ac)
eberlm@62049
   889
qed
eberlm@62049
   890
eberlm@62049
   891
lemma Digamma_one_half: "Digamma (1/2) = - euler_mascheroni - of_real (2 * ln 2)"
eberlm@62049
   892
  using Digamma_half_integer[of 0] by simp
eberlm@62049
   893
eberlm@62049
   894
lemma Digamma_real_three_halves_pos: "Digamma (3/2 :: real) > 0"
eberlm@62049
   895
proof -
eberlm@62049
   896
  have "-Digamma (3/2 :: real) = -Digamma (of_nat 1 + 1/2)" by simp
eberlm@62049
   897
  also have "\<dots> = 2 * ln 2 + euler_mascheroni - 2" by (subst Digamma_half_integer) simp
eberlm@62085
   898
  also note euler_mascheroni_less_13_over_22
eberlm@62085
   899
  also note ln2_le_25_over_36
eberlm@62049
   900
  finally show ?thesis by simp
eberlm@62049
   901
qed
eberlm@62049
   902
eberlm@62049
   903
eberlm@68624
   904
theorem has_field_derivative_Polygamma [derivative_intros]:
eberlm@62049
   905
  fixes z :: "'a :: {real_normed_field,euclidean_space}"
eberlm@62049
   906
  assumes z: "z \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
   907
  shows "(Polygamma n has_field_derivative Polygamma (Suc n) z) (at z within A)"
eberlm@62049
   908
proof (rule has_field_derivative_at_within, cases "n = 0")
eberlm@62049
   909
  assume n: "n = 0"
eberlm@62049
   910
  let ?f = "\<lambda>k z. inverse (of_nat (Suc k)) - inverse (z + of_nat k)"
eberlm@62049
   911
  let ?F = "\<lambda>z. \<Sum>k. ?f k z" and ?f' = "\<lambda>k z. inverse ((z + of_nat k)\<^sup>2)"
eberlm@62049
   912
  from no_nonpos_Int_in_ball'[OF z] guess d . note d = this
eberlm@62049
   913
  from z have summable: "summable (\<lambda>k. inverse (of_nat (Suc k)) - inverse (z + of_nat k))"
eberlm@62049
   914
    by (intro summable_Digamma) force
eberlm@62049
   915
  from z have conv: "uniformly_convergent_on (ball z d) (\<lambda>k z. \<Sum>i<k. inverse ((z + of_nat i)\<^sup>2))"
eberlm@62049
   916
    by (intro Polygamma_converges) auto
eberlm@62049
   917
  with d have "summable (\<lambda>k. inverse ((z + of_nat k)\<^sup>2))" unfolding summable_iff_convergent
eberlm@62049
   918
    by (auto dest!: uniformly_convergent_imp_convergent simp: summable_iff_convergent )
eberlm@62049
   919
eberlm@62049
   920
  have "(?F has_field_derivative (\<Sum>k. ?f' k z)) (at z)"
eberlm@62049
   921
  proof (rule has_field_derivative_series'[of "ball z d" _ _ z])
eberlm@62049
   922
    fix k :: nat and t :: 'a assume t: "t \<in> ball z d"
paulson@62131
   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)"
paulson@62131
   924
      by (auto intro!: derivative_eq_intros simp: power2_eq_square simp del: of_nat_Suc
eberlm@62049
   925
               dest!: plus_of_nat_eq_0_imp elim!: nonpos_Ints_cases)
eberlm@62049
   926
  qed (insert d(1) summable conv, (assumption|simp)+)
eberlm@62049
   927
  with z show "(Polygamma n has_field_derivative Polygamma (Suc n) z) (at z)"
eberlm@62049
   928
    unfolding Digamma_def [abs_def] Polygamma_def [abs_def] using n
eberlm@62049
   929
    by (force simp: power2_eq_square intro!: derivative_eq_intros)
eberlm@62049
   930
next
eberlm@62049
   931
  assume n: "n \<noteq> 0"
eberlm@62049
   932
  from z have z': "z \<noteq> 0" by auto
eberlm@62049
   933
  from no_nonpos_Int_in_ball'[OF z] guess d . note d = this
wenzelm@63040
   934
  define n' where "n' = Suc n"
eberlm@62049
   935
  from n have n': "n' \<ge> 2" by (simp add: n'_def)
eberlm@62049
   936
  have "((\<lambda>z. \<Sum>k. inverse ((z + of_nat k) ^ n')) has_field_derivative
eberlm@62049
   937
                (\<Sum>k. - of_nat n' * inverse ((z + of_nat k) ^ (n'+1)))) (at z)"
eberlm@62049
   938
  proof (rule has_field_derivative_series'[of "ball z d" _ _ z])
eberlm@62049
   939
    fix k :: nat and t :: 'a assume t: "t \<in> ball z d"
eberlm@62049
   940
    with d have t': "t \<notin> \<int>\<^sub>\<le>\<^sub>0" "t \<noteq> 0" by auto
eberlm@62049
   941
    show "((\<lambda>a. inverse ((a + of_nat k) ^ n')) has_field_derivative
eberlm@62049
   942
                - of_nat n' * inverse ((t + of_nat k) ^ (n'+1))) (at t within ball z d)" using t'
eberlm@62049
   943
      by (fastforce intro!: derivative_eq_intros simp: divide_simps power_diff dest: plus_of_nat_eq_0_imp)
eberlm@62049
   944
  next
paulson@62131
   945
    have "uniformly_convergent_on (ball z d)
paulson@62131
   946
              (\<lambda>k z. (- of_nat n' :: 'a) * (\<Sum>i<k. inverse ((z + of_nat i) ^ (n'+1))))"
eberlm@62049
   947
      using z' n by (intro uniformly_convergent_mult Polygamma_converges) (simp_all add: n'_def)
paulson@62131
   948
    thus "uniformly_convergent_on (ball z d)
eberlm@62049
   949
              (\<lambda>k z. \<Sum>i<k. - of_nat n' * inverse ((z + of_nat i :: 'a) ^ (n'+1)))"
nipkow@64267
   950
      by (subst (asm) sum_distrib_left) simp
eberlm@62049
   951
  qed (insert Polygamma_converges'[OF z' n'] d, simp_all)
eberlm@62049
   952
  also have "(\<Sum>k. - of_nat n' * inverse ((z + of_nat k) ^ (n' + 1))) =
eberlm@62049
   953
               (- of_nat n') * (\<Sum>k. inverse ((z + of_nat k) ^ (n' + 1)))"
eberlm@62049
   954
    using Polygamma_converges'[OF z', of "n'+1"] n' by (subst suminf_mult) simp_all
eberlm@62049
   955
  finally have "((\<lambda>z. \<Sum>k. inverse ((z + of_nat k) ^ n')) has_field_derivative
eberlm@62049
   956
                    - of_nat n' * (\<Sum>k. inverse ((z + of_nat k) ^ (n' + 1)))) (at z)" .
eberlm@62049
   957
  from DERIV_cmult[OF this, of "(-1)^Suc n * fact n :: 'a"]
eberlm@62049
   958
    show "(Polygamma n has_field_derivative Polygamma (Suc n) z) (at z)"
eberlm@62049
   959
    unfolding n'_def Polygamma_def[abs_def] using n by (simp add: algebra_simps)
eberlm@62049
   960
qed
eberlm@62049
   961
eberlm@62049
   962
declare has_field_derivative_Polygamma[THEN DERIV_chain2, derivative_intros]
eberlm@62049
   963
paulson@62131
   964
lemma isCont_Polygamma [continuous_intros]:
eberlm@62049
   965
  fixes f :: "_ \<Rightarrow> 'a :: {real_normed_field,euclidean_space}"
eberlm@62049
   966
  shows "isCont f z \<Longrightarrow> f z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> isCont (\<lambda>x. Polygamma n (f x)) z"
eberlm@62049
   967
  by (rule isCont_o2[OF _  DERIV_isCont[OF has_field_derivative_Polygamma]])
eberlm@62049
   968
paulson@62131
   969
lemma continuous_on_Polygamma:
eberlm@62049
   970
  "A \<inter> \<int>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> continuous_on A (Polygamma n :: _ \<Rightarrow> 'a :: {real_normed_field,euclidean_space})"
eberlm@62049
   971
  by (intro continuous_at_imp_continuous_on isCont_Polygamma[OF continuous_ident] ballI) blast
eberlm@62049
   972
paulson@62131
   973
lemma isCont_ln_Gamma_complex [continuous_intros]:
paulson@62131
   974
  fixes f :: "'a::t2_space \<Rightarrow> complex"
paulson@62131
   975
  shows "isCont f z \<Longrightarrow> f z \<notin> \<real>\<^sub>\<le>\<^sub>0 \<Longrightarrow> isCont (\<lambda>z. ln_Gamma (f z)) z"
eberlm@62049
   976
  by (rule isCont_o2[OF _  DERIV_isCont[OF has_field_derivative_ln_Gamma_complex]])
eberlm@62049
   977
eberlm@62049
   978
lemma continuous_on_ln_Gamma_complex [continuous_intros]:
paulson@62131
   979
  fixes A :: "complex set"
paulson@62131
   980
  shows "A \<inter> \<real>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> continuous_on A ln_Gamma"
paulson@62131
   981
  by (intro continuous_at_imp_continuous_on ballI isCont_ln_Gamma_complex[OF continuous_ident])
eberlm@62049
   982
     fastforce
lp15@66512
   983
eberlm@64969
   984
lemma deriv_Polygamma:
eberlm@64969
   985
  assumes "z \<notin> \<int>\<^sub>\<le>\<^sub>0"
lp15@66512
   986
  shows   "deriv (Polygamma m) z =
eberlm@64969
   987
             Polygamma (Suc m) (z :: 'a :: {real_normed_field,euclidean_space})"
eberlm@64969
   988
  by (intro DERIV_imp_deriv has_field_derivative_Polygamma assms)
eberlm@64969
   989
    thm has_field_derivative_Polygamma
eberlm@64969
   990
eberlm@64969
   991
lemma higher_deriv_Polygamma:
eberlm@64969
   992
  assumes "z \<notin> \<int>\<^sub>\<le>\<^sub>0"
lp15@66512
   993
  shows   "(deriv ^^ n) (Polygamma m) z =
eberlm@64969
   994
             Polygamma (m + n) (z :: 'a :: {real_normed_field,euclidean_space})"
eberlm@64969
   995
proof -
eberlm@64969
   996
  have "eventually (\<lambda>u. (deriv ^^ n) (Polygamma m) u = Polygamma (m + n) u) (nhds z)"
eberlm@64969
   997
  proof (induction n)
eberlm@64969
   998
    case (Suc n)
eberlm@64969
   999
    from Suc.IH have "eventually (\<lambda>z. eventually (\<lambda>u. (deriv ^^ n) (Polygamma m) u = Polygamma (m + n) u) (nhds z)) (nhds z)"
eberlm@64969
  1000
      by (simp add: eventually_eventually)
lp15@66512
  1001
    hence "eventually (\<lambda>z. deriv ((deriv ^^ n) (Polygamma m)) z =
eberlm@64969
  1002
             deriv (Polygamma (m + n)) z) (nhds z)"
eberlm@64969
  1003
      by eventually_elim (intro deriv_cong_ev refl)
eberlm@64969
  1004
    moreover have "eventually (\<lambda>z. z \<in> UNIV - \<int>\<^sub>\<le>\<^sub>0) (nhds z)" using assms
eberlm@64969
  1005
      by (intro eventually_nhds_in_open open_Diff open_UNIV) auto
eberlm@64969
  1006
    ultimately show ?case by eventually_elim (simp_all add: deriv_Polygamma)
eberlm@64969
  1007
  qed simp_all
eberlm@64969
  1008
  thus ?thesis by (rule eventually_nhds_x_imp_x)
eberlm@64969
  1009
qed
eberlm@64969
  1010
eberlm@64969
  1011
lemma deriv_ln_Gamma_complex:
eberlm@64969
  1012
  assumes "z \<notin> \<real>\<^sub>\<le>\<^sub>0"
eberlm@64969
  1013
  shows   "deriv ln_Gamma z = Digamma (z :: complex)"
eberlm@64969
  1014
  by (intro DERIV_imp_deriv has_field_derivative_ln_Gamma_complex assms)
eberlm@64969
  1015
eberlm@64969
  1016
eberlm@62049
  1017
eberlm@62049
  1018
text \<open>
paulson@62131
  1019
  We define a type class that captures all the fundamental properties of the inverse of the Gamma function
paulson@62131
  1020
  and defines the Gamma function upon that. This allows us to instantiate the type class both for
paulson@62131
  1021
  the reals and for the complex numbers with a minimal amount of proof duplication.
eberlm@62049
  1022
\<close>
eberlm@62049
  1023
eberlm@68624
  1024
class%unimportant Gamma = real_normed_field + complete_space +
eberlm@62049
  1025
  fixes rGamma :: "'a \<Rightarrow> 'a"
eberlm@62049
  1026
  assumes rGamma_eq_zero_iff_aux: "rGamma z = 0 \<longleftrightarrow> (\<exists>n. z = - of_nat n)"
paulson@62131
  1027
  assumes differentiable_rGamma_aux1:
eberlm@62049
  1028
    "(\<And>n. z \<noteq> - of_nat n) \<Longrightarrow>
paulson@62131
  1029
     let d = (THE d. (\<lambda>n. \<Sum>k<n. inverse (of_nat (Suc k)) - inverse (z + of_nat k))
eberlm@62049
  1030
               \<longlonglongrightarrow> d) - scaleR euler_mascheroni 1
paulson@62131
  1031
     in  filterlim (\<lambda>y. (rGamma y - rGamma z + rGamma z * d * (y - z)) /\<^sub>R
eberlm@62049
  1032
                        norm (y - z)) (nhds 0) (at z)"
paulson@62131
  1033
  assumes differentiable_rGamma_aux2:
eberlm@62049
  1034
    "let z = - of_nat n
nipkow@64272
  1035
     in  filterlim (\<lambda>y. (rGamma y - rGamma z - (-1)^n * (prod of_nat {1..n}) * (y - z)) /\<^sub>R
eberlm@62049
  1036
                        norm (y - z)) (nhds 0) (at z)"
paulson@62131
  1037
  assumes rGamma_series_aux: "(\<And>n. z \<noteq> - of_nat n) \<Longrightarrow>
nipkow@64272
  1038
             let fact' = (\<lambda>n. prod of_nat {1..n});
eberlm@62049
  1039
                 exp = (\<lambda>x. THE e. (\<lambda>n. \<Sum>k<n. x^k /\<^sub>R fact k) \<longlonglongrightarrow> e);
eberlm@62049
  1040
                 pochhammer' = (\<lambda>a n. (\<Prod>n = 0..n. a + of_nat n))
paulson@62131
  1041
             in  filterlim (\<lambda>n. pochhammer' z n / (fact' n * exp (z * (ln (of_nat n) *\<^sub>R 1))))
eberlm@62049
  1042
                     (nhds (rGamma z)) sequentially"
eberlm@62049
  1043
begin
eberlm@62049
  1044
subclass banach ..
eberlm@62049
  1045
end
eberlm@62049
  1046
eberlm@62049
  1047
definition "Gamma z = inverse (rGamma z)"
eberlm@62049
  1048
eberlm@62049
  1049
eberlm@62049
  1050
subsection \<open>Basic properties\<close>
eberlm@62049
  1051
eberlm@62049
  1052
lemma Gamma_nonpos_Int: "z \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Gamma z = 0"
eberlm@62049
  1053
  and rGamma_nonpos_Int: "z \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> rGamma z = 0"
eberlm@62049
  1054
  using rGamma_eq_zero_iff_aux[of z] unfolding Gamma_def by (auto elim!: nonpos_Ints_cases')
eberlm@62049
  1055
eberlm@62049
  1056
lemma Gamma_nonzero: "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Gamma z \<noteq> 0"
eberlm@62049
  1057
  and rGamma_nonzero: "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> rGamma z \<noteq> 0"
eberlm@62049
  1058
  using rGamma_eq_zero_iff_aux[of z] unfolding Gamma_def by (auto elim!: nonpos_Ints_cases')
eberlm@62049
  1059
paulson@62131
  1060
lemma Gamma_eq_zero_iff: "Gamma z = 0 \<longleftrightarrow> z \<in> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
  1061
  and rGamma_eq_zero_iff: "rGamma z = 0 \<longleftrightarrow> z \<in> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
  1062
  using rGamma_eq_zero_iff_aux[of z] unfolding Gamma_def by (auto elim!: nonpos_Ints_cases')
eberlm@62049
  1063
eberlm@62049
  1064
lemma rGamma_inverse_Gamma: "rGamma z = inverse (Gamma z)"
eberlm@62049
  1065
  unfolding Gamma_def by simp
eberlm@62049
  1066
eberlm@62049
  1067
lemma rGamma_series_LIMSEQ [tendsto_intros]:
eberlm@62049
  1068
  "rGamma_series z \<longlonglongrightarrow> rGamma z"
eberlm@62049
  1069
proof (cases "z \<in> \<int>\<^sub>\<le>\<^sub>0")
eberlm@62049
  1070
  case False
eberlm@62049
  1071
  hence "z \<noteq> - of_nat n" for n by auto
eberlm@62049
  1072
  from rGamma_series_aux[OF this] show ?thesis
nipkow@64272
  1073
    by (simp add: rGamma_series_def[abs_def] fact_prod pochhammer_Suc_prod
haftmann@63367
  1074
                  exp_def of_real_def[symmetric] suminf_def sums_def[abs_def] atLeast0AtMost)
eberlm@62049
  1075
qed (insert rGamma_eq_zero_iff[of z], simp_all add: rGamma_series_nonpos_Ints_LIMSEQ)
eberlm@62049
  1076
eberlm@68624
  1077
theorem Gamma_series_LIMSEQ [tendsto_intros]:
eberlm@62049
  1078
  "Gamma_series z \<longlonglongrightarrow> Gamma z"
eberlm@62049
  1079
proof (cases "z \<in> \<int>\<^sub>\<le>\<^sub>0")
eberlm@62049
  1080
  case False
eberlm@62049
  1081
  hence "(\<lambda>n. inverse (rGamma_series z n)) \<longlonglongrightarrow> inverse (rGamma z)"
eberlm@62049
  1082
    by (intro tendsto_intros) (simp_all add: rGamma_eq_zero_iff)
paulson@62131
  1083
  also have "(\<lambda>n. inverse (rGamma_series z n)) = Gamma_series z"
eberlm@62049
  1084
    by (simp add: rGamma_series_def Gamma_series_def[abs_def])
eberlm@62049
  1085
  finally show ?thesis by (simp add: Gamma_def)
eberlm@62049
  1086
qed (insert Gamma_eq_zero_iff[of z], simp_all add: Gamma_series_nonpos_Ints_LIMSEQ)
eberlm@62049
  1087
eberlm@62049
  1088
lemma Gamma_altdef: "Gamma z = lim (Gamma_series z)"
eberlm@62049
  1089
  using Gamma_series_LIMSEQ[of z] by (simp add: limI)
eberlm@62049
  1090
eberlm@62049
  1091
lemma rGamma_1 [simp]: "rGamma 1 = 1"
eberlm@62049
  1092
proof -
eberlm@62049
  1093
  have A: "eventually (\<lambda>n. rGamma_series 1 n = of_nat (Suc n) / of_nat n) sequentially"
eberlm@62049
  1094
    using eventually_gt_at_top[of "0::nat"]
paulson@62131
  1095
    by (force elim!: eventually_mono simp: rGamma_series_def exp_of_real pochhammer_fact
eberlm@62049
  1096
                    divide_simps pochhammer_rec' dest!: pochhammer_eq_0_imp_nonpos_Int)
eberlm@62049
  1097
  have "rGamma_series 1 \<longlonglongrightarrow> 1" by (subst tendsto_cong[OF A]) (rule LIMSEQ_Suc_n_over_n)
eberlm@62049
  1098
  moreover have "rGamma_series 1 \<longlonglongrightarrow> rGamma 1" by (rule tendsto_intros)
eberlm@62049
  1099
  ultimately show ?thesis by (intro LIMSEQ_unique)
eberlm@62049
  1100
qed
eberlm@62049
  1101
eberlm@62049
  1102
lemma rGamma_plus1: "z * rGamma (z + 1) = rGamma z"
eberlm@62049
  1103
proof -
eberlm@62049
  1104
  let ?f = "\<lambda>n. (z + 1) * inverse (of_nat n) + 1"
eberlm@62049
  1105
  have "eventually (\<lambda>n. ?f n * rGamma_series z n = z * rGamma_series (z + 1) n) sequentially"
eberlm@62049
  1106
    using eventually_gt_at_top[of "0::nat"]
eberlm@62049
  1107
  proof eventually_elim
eberlm@62049
  1108
    fix n :: nat assume n: "n > 0"
eberlm@62049
  1109
    hence "z * rGamma_series (z + 1) n = inverse (of_nat n) *
eberlm@62049
  1110
             pochhammer z (Suc (Suc n)) / (fact n * exp (z * of_real (ln (of_nat n))))"
eberlm@62049
  1111
      by (subst pochhammer_rec) (simp add: rGamma_series_def field_simps exp_add exp_of_real)
eberlm@62049
  1112
    also from n have "\<dots> = ?f n * rGamma_series z n"
eberlm@62049
  1113
      by (subst pochhammer_rec') (simp_all add: divide_simps rGamma_series_def add_ac)
eberlm@62049
  1114
    finally show "?f n * rGamma_series z n = z * rGamma_series (z + 1) n" ..
eberlm@62049
  1115
  qed
eberlm@62049
  1116
  moreover have "(\<lambda>n. ?f n * rGamma_series z n) \<longlonglongrightarrow> ((z+1) * 0 + 1) * rGamma z"
eberlm@62049
  1117
    by (intro tendsto_intros lim_inverse_n)
eberlm@62049
  1118
  hence "(\<lambda>n. ?f n * rGamma_series z n) \<longlonglongrightarrow> rGamma z" by simp
eberlm@62049
  1119
  ultimately have "(\<lambda>n. z * rGamma_series (z + 1) n) \<longlonglongrightarrow> rGamma z"
eberlm@62049
  1120
    by (rule Lim_transform_eventually)
eberlm@62049
  1121
  moreover have "(\<lambda>n. z * rGamma_series (z + 1) n) \<longlonglongrightarrow> z * rGamma (z + 1)"
eberlm@62049
  1122
    by (intro tendsto_intros)
eberlm@62049
  1123
  ultimately show "z * rGamma (z + 1) = rGamma z" using LIMSEQ_unique by blast
eberlm@62049
  1124
qed
eberlm@62049
  1125
eberlm@62049
  1126
eberlm@62049
  1127
lemma pochhammer_rGamma: "rGamma z = pochhammer z n * rGamma (z + of_nat n)"
eberlm@62049
  1128
proof (induction n arbitrary: z)
eberlm@62049
  1129
  case (Suc n z)
eberlm@62049
  1130
  have "rGamma z = pochhammer z n * rGamma (z + of_nat n)" by (rule Suc.IH)
eberlm@62049
  1131
  also note rGamma_plus1 [symmetric]
eberlm@62049
  1132
  finally show ?case by (simp add: add_ac pochhammer_rec')
eberlm@62049
  1133
qed simp_all
eberlm@62049
  1134
eberlm@68624
  1135
theorem Gamma_plus1: "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Gamma (z + 1) = z * Gamma z"
eberlm@62049
  1136
  using rGamma_plus1[of z] by (simp add: rGamma_inverse_Gamma field_simps Gamma_eq_zero_iff)
eberlm@62049
  1137
eberlm@68624
  1138
theorem pochhammer_Gamma: "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> pochhammer z n = Gamma (z + of_nat n) / Gamma z"
paulson@62131
  1139
  using pochhammer_rGamma[of z]
eberlm@62049
  1140
  by (simp add: rGamma_inverse_Gamma Gamma_eq_zero_iff field_simps)
eberlm@62049
  1141
paulson@62131
  1142
lemma Gamma_0 [simp]: "Gamma 0 = 0"
eberlm@62049
  1143
  and rGamma_0 [simp]: "rGamma 0 = 0"
eberlm@62049
  1144
  and Gamma_neg_1 [simp]: "Gamma (- 1) = 0"
eberlm@62049
  1145
  and rGamma_neg_1 [simp]: "rGamma (- 1) = 0"
eberlm@62049
  1146
  and Gamma_neg_numeral [simp]: "Gamma (- numeral n) = 0"
eberlm@62049
  1147
  and rGamma_neg_numeral [simp]: "rGamma (- numeral n) = 0"
eberlm@62049
  1148
  and Gamma_neg_of_nat [simp]: "Gamma (- of_nat m) = 0"
eberlm@62049
  1149
  and rGamma_neg_of_nat [simp]: "rGamma (- of_nat m) = 0"
eberlm@62049
  1150
  by (simp_all add: rGamma_eq_zero_iff Gamma_eq_zero_iff)
eberlm@62049
  1151
eberlm@62049
  1152
lemma Gamma_1 [simp]: "Gamma 1 = 1" unfolding Gamma_def by simp
eberlm@62049
  1153
eberlm@68624
  1154
theorem Gamma_fact: "Gamma (1 + of_nat n) = fact n"
nipkow@68403
  1155
  by (simp add: pochhammer_fact pochhammer_Gamma of_nat_in_nonpos_Ints_iff flip: of_nat_Suc)
eberlm@62049
  1156
eberlm@62049
  1157
lemma Gamma_numeral: "Gamma (numeral n) = fact (pred_numeral n)"
hoelzl@63594
  1158
  by (subst of_nat_numeral[symmetric], subst numeral_eq_Suc,
eberlm@63295
  1159
      subst of_nat_Suc, subst Gamma_fact) (rule refl)
eberlm@62049
  1160
eberlm@62049
  1161
lemma Gamma_of_int: "Gamma (of_int n) = (if n > 0 then fact (nat (n - 1)) else 0)"
eberlm@62049
  1162
proof (cases "n > 0")
eberlm@62049
  1163
  case True
eberlm@62049
  1164
  hence "Gamma (of_int n) = Gamma (of_nat (Suc (nat (n - 1))))" by (subst of_nat_Suc) simp_all
eberlm@63295
  1165
  with True show ?thesis by (subst (asm) of_nat_Suc, subst (asm) Gamma_fact) simp
eberlm@62049
  1166
qed (simp_all add: Gamma_eq_zero_iff nonpos_Ints_of_int)
eberlm@62049
  1167
eberlm@62049
  1168
lemma rGamma_of_int: "rGamma (of_int n) = (if n > 0 then inverse (fact (nat (n - 1))) else 0)"
eberlm@62049
  1169
  by (simp add: Gamma_of_int rGamma_inverse_Gamma)
eberlm@62049
  1170
eberlm@62049
  1171
lemma Gamma_seriesI:
eberlm@62049
  1172
  assumes "(\<lambda>n. g n / Gamma_series z n) \<longlonglongrightarrow> 1"
eberlm@62049
  1173
  shows   "g \<longlonglongrightarrow> Gamma z"
eberlm@62049
  1174
proof (rule Lim_transform_eventually)
eberlm@62049
  1175
  have "1/2 > (0::real)" by simp
paulson@62131
  1176
  from tendstoD[OF assms, OF this]
eberlm@62049
  1177
    show "eventually (\<lambda>n. g n / Gamma_series z n * Gamma_series z n = g n) sequentially"
eberlm@62049
  1178
    by (force elim!: eventually_mono simp: dist_real_def dist_0_norm)
eberlm@62049
  1179
  from assms have "(\<lambda>n. g n / Gamma_series z n * Gamma_series z n) \<longlonglongrightarrow> 1 * Gamma z"
eberlm@62049
  1180
    by (intro tendsto_intros)
eberlm@62049
  1181
  thus "(\<lambda>n. g n / Gamma_series z n * Gamma_series z n) \<longlonglongrightarrow> Gamma z" by simp
eberlm@62049
  1182
qed
eberlm@62049
  1183
eberlm@62049
  1184
lemma Gamma_seriesI':
eberlm@62049
  1185
  assumes "f \<longlonglongrightarrow> rGamma z"
eberlm@62049
  1186
  assumes "(\<lambda>n. g n * f n) \<longlonglongrightarrow> 1"
eberlm@62049
  1187
  assumes "z \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
  1188
  shows   "g \<longlonglongrightarrow> Gamma z"
eberlm@62049
  1189
proof (rule Lim_transform_eventually)
eberlm@62049
  1190
  have "1/2 > (0::real)" by simp
eberlm@62049
  1191
  from tendstoD[OF assms(2), OF this] show "eventually (\<lambda>n. g n * f n / f n = g n) sequentially"
eberlm@62049
  1192
    by (force elim!: eventually_mono simp: dist_real_def dist_0_norm)
eberlm@62049
  1193
  from assms have "(\<lambda>n. g n * f n / f n) \<longlonglongrightarrow> 1 / rGamma z"
eberlm@62049
  1194
    by (intro tendsto_divide assms) (simp_all add: rGamma_eq_zero_iff)
eberlm@62049
  1195
  thus "(\<lambda>n. g n * f n / f n) \<longlonglongrightarrow> Gamma z" by (simp add: Gamma_def divide_inverse)
eberlm@62049
  1196
qed
eberlm@62049
  1197
eberlm@62049
  1198
lemma Gamma_series'_LIMSEQ: "Gamma_series' z \<longlonglongrightarrow> Gamma z"
paulson@62131
  1199
  by (cases "z \<in> \<int>\<^sub>\<le>\<^sub>0") (simp_all add: Gamma_nonpos_Int Gamma_seriesI[OF Gamma_series_Gamma_series']
eberlm@62049
  1200
                                      Gamma_series'_nonpos_Ints_LIMSEQ[of z])
eberlm@62049
  1201
eberlm@62049
  1202
eberlm@62049
  1203
subsection \<open>Differentiability\<close>
eberlm@62049
  1204
paulson@62131
  1205
lemma has_field_derivative_rGamma_no_nonpos_int:
eberlm@62049
  1206
  assumes "z \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
  1207
  shows   "(rGamma has_field_derivative -rGamma z * Digamma z) (at z within A)"
eberlm@62049
  1208
proof (rule has_field_derivative_at_within)
eberlm@62049
  1209
  from assms have "z \<noteq> - of_nat n" for n by auto
paulson@62131
  1210
  from differentiable_rGamma_aux1[OF this]
eberlm@62049
  1211
    show "(rGamma has_field_derivative -rGamma z * Digamma z) (at z)"
eberlm@62049
  1212
         unfolding Digamma_def suminf_def sums_def[abs_def]
eberlm@62049
  1213
                   has_field_derivative_def has_derivative_def netlimit_at
eberlm@62049
  1214
    by (simp add: Let_def bounded_linear_mult_right mult_ac of_real_def [symmetric])
eberlm@62049
  1215
qed
eberlm@62049
  1216
paulson@62131
  1217
lemma has_field_derivative_rGamma_nonpos_int:
paulson@62131
  1218
  "(rGamma has_field_derivative (-1)^n * fact n) (at (- of_nat n) within A)"
eberlm@62049
  1219
  apply (rule has_field_derivative_at_within)
eberlm@62049
  1220
  using differentiable_rGamma_aux2[of n]
eberlm@62049
  1221
  unfolding Let_def has_field_derivative_def has_derivative_def netlimit_at
nipkow@64272
  1222
  by (simp only: bounded_linear_mult_right mult_ac of_real_def [symmetric] fact_prod) simp
eberlm@62049
  1223
paulson@62131
  1224
lemma has_field_derivative_rGamma [derivative_intros]:
eberlm@62049
  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>)
eberlm@62049
  1226
      else -rGamma z * Digamma z)) (at z within A)"
eberlm@62049
  1227
using has_field_derivative_rGamma_no_nonpos_int[of z A]
eberlm@62049
  1228
      has_field_derivative_rGamma_nonpos_int[of "nat \<lfloor>norm z\<rfloor>" A]
eberlm@62049
  1229
  by (auto elim!: nonpos_Ints_cases')
eberlm@62049
  1230
eberlm@62049
  1231
declare has_field_derivative_rGamma_no_nonpos_int [THEN DERIV_chain2, derivative_intros]
eberlm@62049
  1232
declare has_field_derivative_rGamma [THEN DERIV_chain2, derivative_intros]
eberlm@62049
  1233
declare has_field_derivative_rGamma_nonpos_int [derivative_intros]
eberlm@62049
  1234
declare has_field_derivative_rGamma_no_nonpos_int [derivative_intros]
eberlm@62049
  1235
declare has_field_derivative_rGamma [derivative_intros]
eberlm@62049
  1236
eberlm@68624
  1237
theorem has_field_derivative_Gamma [derivative_intros]:
eberlm@62049
  1238
  "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> (Gamma has_field_derivative Gamma z * Digamma z) (at z within A)"
eberlm@62049
  1239
  unfolding Gamma_def [abs_def]
eberlm@62049
  1240
  by (fastforce intro!: derivative_eq_intros simp: rGamma_eq_zero_iff)
eberlm@62049
  1241
eberlm@62049
  1242
declare has_field_derivative_Gamma[THEN DERIV_chain2, derivative_intros]
eberlm@62049
  1243
eberlm@62049
  1244
(* TODO: Hide ugly facts properly *)
eberlm@62049
  1245
hide_fact rGamma_eq_zero_iff_aux differentiable_rGamma_aux1 differentiable_rGamma_aux2
eberlm@62049
  1246
          differentiable_rGamma_aux2 rGamma_series_aux Gamma_class.rGamma_eq_zero_iff_aux
eberlm@62049
  1247
eberlm@62049
  1248
lemma continuous_on_rGamma [continuous_intros]: "continuous_on A rGamma"
eberlm@62049
  1249
  by (rule DERIV_continuous_on has_field_derivative_rGamma)+
eberlm@62049
  1250
eberlm@62049
  1251
lemma continuous_on_Gamma [continuous_intros]: "A \<inter> \<int>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> continuous_on A Gamma"
eberlm@62049
  1252
  by (rule DERIV_continuous_on has_field_derivative_Gamma)+ blast
eberlm@62049
  1253
paulson@62131
  1254
lemma isCont_rGamma [continuous_intros]:
eberlm@62049
  1255
  "isCont f z \<Longrightarrow> isCont (\<lambda>x. rGamma (f x)) z"
eberlm@62049
  1256
  by (rule isCont_o2[OF _  DERIV_isCont[OF has_field_derivative_rGamma]])
eberlm@62049
  1257
paulson@62131
  1258
lemma isCont_Gamma [continuous_intros]:
eberlm@62049
  1259
  "isCont f z \<Longrightarrow> f z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> isCont (\<lambda>x. Gamma (f x)) z"
eberlm@62049
  1260
  by (rule isCont_o2[OF _  DERIV_isCont[OF has_field_derivative_Gamma]])
eberlm@62049
  1261
eberlm@68624
  1262
subsection%unimportant \<open>The complex Gamma function\<close>
eberlm@68624
  1263
eberlm@68624
  1264
instantiation%unimportant complex :: Gamma
eberlm@62049
  1265
begin
eberlm@62049
  1266
eberlm@68624
  1267
definition%unimportant rGamma_complex :: "complex \<Rightarrow> complex" where
eberlm@62049
  1268
  "rGamma_complex z = lim (rGamma_series z)"
eberlm@62049
  1269
paulson@62131
  1270
lemma rGamma_series_complex_converges:
eberlm@62049
  1271
        "convergent (rGamma_series (z :: complex))" (is "?thesis1")
eberlm@62049
  1272
  and rGamma_complex_altdef:
eberlm@62049
  1273
        "rGamma z = (if z \<in> \<int>\<^sub>\<le>\<^sub>0 then 0 else exp (-ln_Gamma z))" (is "?thesis2")
eberlm@62049
  1274
proof -
eberlm@62049
  1275
  have "?thesis1 \<and> ?thesis2"
eberlm@62049
  1276
  proof (cases "z \<in> \<int>\<^sub>\<le>\<^sub>0")
eberlm@62049
  1277
    case False
eberlm@62049
  1278
    have "rGamma_series z \<longlonglongrightarrow> exp (- ln_Gamma z)"
eberlm@62049
  1279
    proof (rule Lim_transform_eventually)
eberlm@62049
  1280
      from ln_Gamma_series_complex_converges'[OF False] guess d by (elim exE conjE)
paulson@62131
  1281
      from this(1) uniformly_convergent_imp_convergent[OF this(2), of z]
eberlm@62049
  1282
        have "ln_Gamma_series z \<longlonglongrightarrow> lim (ln_Gamma_series z)" by (simp add: convergent_LIMSEQ_iff)
eberlm@62049
  1283
      thus "(\<lambda>n. exp (-ln_Gamma_series z n)) \<longlonglongrightarrow> exp (- ln_Gamma z)"
eberlm@62049
  1284
        unfolding convergent_def ln_Gamma_def by (intro tendsto_exp tendsto_minus)
eberlm@62049
  1285
      from eventually_gt_at_top[of "0::nat"] exp_ln_Gamma_series_complex False
eberlm@62049
  1286
        show "eventually (\<lambda>n. exp (-ln_Gamma_series z n) = rGamma_series z n) sequentially"
eberlm@62049
  1287
        by (force elim!: eventually_mono simp: exp_minus Gamma_series_def rGamma_series_def)
eberlm@62049
  1288
    qed
eberlm@62049
  1289
    with False show ?thesis
eberlm@62049
  1290
      by (auto simp: convergent_def rGamma_complex_def intro!: limI)
eberlm@62049
  1291
  next
eberlm@62049
  1292
    case True
eberlm@62049
  1293
    then obtain k where "z = - of_nat k" by (erule nonpos_Ints_cases')
eberlm@62049
  1294
    also have "rGamma_series \<dots> \<longlonglongrightarrow> 0"
eberlm@62049
  1295
      by (subst tendsto_cong[OF rGamma_series_minus_of_nat]) (simp_all add: convergent_const)
eberlm@62049
  1296
    finally show ?thesis using True
eberlm@62049
  1297
      by (auto simp: rGamma_complex_def convergent_def intro!: limI)
eberlm@62049
  1298
  qed
eberlm@62049
  1299
  thus "?thesis1" "?thesis2" by blast+
eberlm@62049
  1300
qed
eberlm@62049
  1301
eberlm@68624
  1302
context%unimportant
eberlm@62049
  1303
begin
eberlm@62049
  1304
eberlm@62049
  1305
(* TODO: duplication *)
eberlm@62049
  1306
private lemma rGamma_complex_plus1: "z * rGamma (z + 1) = rGamma (z :: complex)"
eberlm@62049
  1307
proof -
eberlm@62049
  1308
  let ?f = "\<lambda>n. (z + 1) * inverse (of_nat n) + 1"
eberlm@62049
  1309
  have "eventually (\<lambda>n. ?f n * rGamma_series z n = z * rGamma_series (z + 1) n) sequentially"
eberlm@62049
  1310
    using eventually_gt_at_top[of "0::nat"]
eberlm@62049
  1311
  proof eventually_elim
eberlm@62049
  1312
    fix n :: nat assume n: "n > 0"
eberlm@62049
  1313
    hence "z * rGamma_series (z + 1) n = inverse (of_nat n) *
eberlm@62049
  1314
             pochhammer z (Suc (Suc n)) / (fact n * exp (z * of_real (ln (of_nat n))))"
eberlm@62049
  1315
      by (subst pochhammer_rec) (simp add: rGamma_series_def field_simps exp_add exp_of_real)
eberlm@62049
  1316
    also from n have "\<dots> = ?f n * rGamma_series z n"
eberlm@62049
  1317
      by (subst pochhammer_rec') (simp_all add: divide_simps rGamma_series_def add_ac)
eberlm@62049
  1318
    finally show "?f n * rGamma_series z n = z * rGamma_series (z + 1) n" ..
eberlm@62049
  1319
  qed
eberlm@62049
  1320
  moreover have "(\<lambda>n. ?f n * rGamma_series z n) \<longlonglongrightarrow> ((z+1) * 0 + 1) * rGamma z"
eberlm@62049
  1321
    using rGamma_series_complex_converges
paulson@62131
  1322
    by (intro tendsto_intros lim_inverse_n)
eberlm@62049
  1323
       (simp_all add: convergent_LIMSEQ_iff rGamma_complex_def)
eberlm@62049
  1324
  hence "(\<lambda>n. ?f n * rGamma_series z n) \<longlonglongrightarrow> rGamma z" by simp
eberlm@62049
  1325
  ultimately have "(\<lambda>n. z * rGamma_series (z + 1) n) \<longlonglongrightarrow> rGamma z"
eberlm@62049
  1326
    by (rule Lim_transform_eventually)
eberlm@62049
  1327
  moreover have "(\<lambda>n. z * rGamma_series (z + 1) n) \<longlonglongrightarrow> z * rGamma (z + 1)"
paulson@62131
  1328
    using rGamma_series_complex_converges
eberlm@62049
  1329
    by (auto intro!: tendsto_mult simp: rGamma_complex_def convergent_LIMSEQ_iff)
eberlm@62049
  1330
  ultimately show "z * rGamma (z + 1) = rGamma z" using LIMSEQ_unique by blast
eberlm@62049
  1331
qed
eberlm@62049
  1332
eberlm@62049
  1333
private lemma has_field_derivative_rGamma_complex_no_nonpos_Int:
eberlm@62049
  1334
  assumes "(z :: complex) \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
  1335
  shows   "(rGamma has_field_derivative - rGamma z * Digamma z) (at z)"
eberlm@62049
  1336
proof -
paulson@62131
  1337
  have diff: "(rGamma has_field_derivative - rGamma z * Digamma z) (at z)" if "Re z > 0" for z
eberlm@62049
  1338
  proof (subst DERIV_cong_ev[OF refl _ refl])
paulson@62131
  1339
    from that have "eventually (\<lambda>t. t \<in> ball z (Re z/2)) (nhds z)"
eberlm@62049
  1340
      by (intro eventually_nhds_in_nhd) simp_all
eberlm@62049
  1341
    thus "eventually (\<lambda>t. rGamma t = exp (- ln_Gamma t)) (nhds z)"
eberlm@62049
  1342
      using no_nonpos_Int_in_ball_complex[OF that]
eberlm@62049
  1343
      by (auto elim!: eventually_mono simp: rGamma_complex_altdef)
eberlm@62049
  1344
  next
paulson@62131
  1345
    have "z \<notin> \<real>\<^sub>\<le>\<^sub>0" using that by (simp add: complex_nonpos_Reals_iff)
paulson@62131
  1346
    with that show "((\<lambda>t. exp (- ln_Gamma t)) has_field_derivative (-rGamma z * Digamma z)) (at z)"
paulson@62131
  1347
     by (force elim!: nonpos_Ints_cases intro!: derivative_eq_intros simp: rGamma_complex_altdef)
eberlm@62049
  1348
  qed
paulson@62131
  1349
eberlm@62049
  1350
  from assms show "(rGamma has_field_derivative - rGamma z * Digamma z) (at z)"
eberlm@62049
  1351
  proof (induction "nat \<lfloor>1 - Re z\<rfloor>" arbitrary: z)
eberlm@62049
  1352
    case (Suc n z)
eberlm@62049
  1353
    from Suc.prems have z: "z \<noteq> 0" by auto
eberlm@62049
  1354
    from Suc.hyps have "n = nat \<lfloor>- Re z\<rfloor>" by linarith
eberlm@62049
  1355
    hence A: "n = nat \<lfloor>1 - Re (z + 1)\<rfloor>" by simp
eberlm@62049
  1356
    from Suc.prems have B: "z + 1 \<notin> \<int>\<^sub>\<le>\<^sub>0" by (force dest: plus_one_in_nonpos_Ints_imp)
paulson@62131
  1357
paulson@62131
  1358
    have "((\<lambda>z. z * (rGamma \<circ> (\<lambda>z. z + 1)) z) has_field_derivative
eberlm@62049
  1359
      -rGamma (z + 1) * (Digamma (z + 1) * z - 1)) (at z)"
eberlm@62049
  1360
      by (rule derivative_eq_intros DERIV_chain Suc refl A B)+ (simp add: algebra_simps)
paulson@62131
  1361
    also have "(\<lambda>z. z * (rGamma \<circ> (\<lambda>z. z + 1 :: complex)) z) = rGamma"
eberlm@62049
  1362
      by (simp add: rGamma_complex_plus1)
eberlm@62049
  1363
    also from z have "Digamma (z + 1) * z - 1 = z * Digamma z"
eberlm@62049
  1364
      by (subst Digamma_plus1) (simp_all add: field_simps)
eberlm@62049
  1365
    also have "-rGamma (z + 1) * (z * Digamma z) = -rGamma z * Digamma z"
eberlm@62049
  1366
      by (simp add: rGamma_complex_plus1[of z, symmetric])
eberlm@62049
  1367
    finally show ?case .
eberlm@62049
  1368
  qed (intro diff, simp)
eberlm@62049
  1369
qed
eberlm@62049
  1370
eberlm@62049
  1371
private lemma rGamma_complex_1: "rGamma (1 :: complex) = 1"
eberlm@62049
  1372
proof -
eberlm@62049
  1373
  have A: "eventually (\<lambda>n. rGamma_series 1 n = of_nat (Suc n) / of_nat n) sequentially"
eberlm@62049
  1374
    using eventually_gt_at_top[of "0::nat"]
paulson@62131
  1375
    by (force elim!: eventually_mono simp: rGamma_series_def exp_of_real pochhammer_fact
eberlm@62049
  1376
                    divide_simps pochhammer_rec' dest!: pochhammer_eq_0_imp_nonpos_Int)
eberlm@62049
  1377
  have "rGamma_series 1 \<longlonglongrightarrow> 1" by (subst tendsto_cong[OF A]) (rule LIMSEQ_Suc_n_over_n)
eberlm@62049
  1378
  thus "rGamma 1 = (1 :: complex)" unfolding rGamma_complex_def by (rule limI)
eberlm@62049
  1379
qed
eberlm@62049
  1380
eberlm@62049
  1381
private lemma has_field_derivative_rGamma_complex_nonpos_Int:
eberlm@62049
  1382
  "(rGamma has_field_derivative (-1)^n * fact n) (at (- of_nat n :: complex))"
eberlm@62049
  1383
proof (induction n)
eberlm@62049
  1384
  case 0
eberlm@62049
  1385
  have A: "(0::complex) + 1 \<notin> \<int>\<^sub>\<le>\<^sub>0" by simp
eberlm@62049
  1386
  have "((\<lambda>z. z * (rGamma \<circ> (\<lambda>z. z + 1 :: complex)) z) has_field_derivative 1) (at 0)"
eberlm@62049
  1387
    by (rule derivative_eq_intros DERIV_chain refl
eberlm@62049
  1388
             has_field_derivative_rGamma_complex_no_nonpos_Int A)+ (simp add: rGamma_complex_1)
eberlm@62049
  1389
    thus ?case by (simp add: rGamma_complex_plus1)
eberlm@62049
  1390
next
eberlm@62049
  1391
  case (Suc n)
paulson@62131
  1392
  hence A: "(rGamma has_field_derivative (-1)^n * fact n)
eberlm@62049
  1393
                (at (- of_nat (Suc n) + 1 :: complex))" by simp
paulson@62131
  1394
   have "((\<lambda>z. z * (rGamma \<circ> (\<lambda>z. z + 1 :: complex)) z) has_field_derivative
eberlm@62049
  1395
             (- 1) ^ Suc n * fact (Suc n)) (at (- of_nat (Suc n)))"
paulson@62131
  1396
     by (rule derivative_eq_intros refl A DERIV_chain)+
eberlm@62049
  1397
        (simp add: algebra_simps rGamma_complex_altdef)
eberlm@62049
  1398
  thus ?case by (simp add: rGamma_complex_plus1)
eberlm@62049
  1399
qed
eberlm@62049
  1400
eberlm@62049
  1401
instance proof
eberlm@62049
  1402
  fix z :: complex show "(rGamma z = 0) \<longleftrightarrow> (\<exists>n. z = - of_nat n)"
eberlm@62049
  1403
    by (auto simp: rGamma_complex_altdef elim!: nonpos_Ints_cases')
eberlm@62049
  1404
next
eberlm@62049
  1405
  fix z :: complex assume "\<And>n. z \<noteq> - of_nat n"
eberlm@62049
  1406
  hence "z \<notin> \<int>\<^sub>\<le>\<^sub>0" by (auto elim!: nonpos_Ints_cases')
eberlm@62049
  1407
  from has_field_derivative_rGamma_complex_no_nonpos_Int[OF this]
eberlm@62049
  1408
    show "let d = (THE d. (\<lambda>n. \<Sum>k<n. inverse (of_nat (Suc k)) - inverse (z + of_nat k))
paulson@62131
  1409
                       \<longlonglongrightarrow> d) - euler_mascheroni *\<^sub>R 1 in  (\<lambda>y. (rGamma y - rGamma z +
eberlm@62049
  1410
              rGamma z * d * (y - z)) /\<^sub>R  cmod (y - z)) \<midarrow>z\<rightarrow> 0"
eberlm@62049
  1411
      by (simp add: has_field_derivative_def has_derivative_def Digamma_def sums_def [abs_def]
eberlm@62049
  1412
                    netlimit_at of_real_def[symmetric] suminf_def)
eberlm@62049
  1413
next
eberlm@62049
  1414
  fix n :: nat
eberlm@62049
  1415
  from has_field_derivative_rGamma_complex_nonpos_Int[of n]
nipkow@64272
  1416
  show "let z = - of_nat n in (\<lambda>y. (rGamma y - rGamma z - (- 1) ^ n * prod of_nat {1..n} *
eberlm@62049
  1417
                  (y - z)) /\<^sub>R cmod (y - z)) \<midarrow>z\<rightarrow> 0"
nipkow@64272
  1418
    by (simp add: has_field_derivative_def has_derivative_def fact_prod netlimit_at Let_def)
eberlm@62049
  1419
next
eberlm@62049
  1420
  fix z :: complex
eberlm@62049
  1421
  from rGamma_series_complex_converges[of z] have "rGamma_series z \<longlonglongrightarrow> rGamma z"
eberlm@62049
  1422
    by (simp add: convergent_LIMSEQ_iff rGamma_complex_def)
nipkow@64272
  1423
  thus "let fact' = \<lambda>n. prod of_nat {1..n};
eberlm@62049
  1424
            exp = \<lambda>x. THE e. (\<lambda>n. \<Sum>k<n. x ^ k /\<^sub>R fact k) \<longlonglongrightarrow> e;
eberlm@62049
  1425
            pochhammer' = \<lambda>a n. \<Prod>n = 0..n. a + of_nat n
eberlm@62049
  1426
        in  (\<lambda>n. pochhammer' z n / (fact' n * exp (z * ln (real_of_nat n) *\<^sub>R 1))) \<longlonglongrightarrow> rGamma z"
nipkow@64272
  1427
    by (simp add: fact_prod pochhammer_Suc_prod rGamma_series_def [abs_def] exp_def
haftmann@63367
  1428
                  of_real_def [symmetric] suminf_def sums_def [abs_def] atLeast0AtMost)
eberlm@62049
  1429
qed
eberlm@62049
  1430
eberlm@62049
  1431
end
eberlm@62049
  1432
end
eberlm@62049
  1433
eberlm@62049
  1434
paulson@62131
  1435
lemma Gamma_complex_altdef:
eberlm@62049
  1436
  "Gamma z = (if z \<in> \<int>\<^sub>\<le>\<^sub>0 then 0 else exp (ln_Gamma (z :: complex)))"
eberlm@62049
  1437
  unfolding Gamma_def rGamma_complex_altdef by (simp add: exp_minus)
eberlm@62049
  1438
eberlm@62049
  1439
lemma cnj_rGamma: "cnj (rGamma z) = rGamma (cnj z)"
eberlm@62049
  1440
proof -
eberlm@62049
  1441
  have "rGamma_series (cnj z) = (\<lambda>n. cnj (rGamma_series z n))"
eberlm@62049
  1442
    by (intro ext) (simp_all add: rGamma_series_def exp_cnj)
eberlm@62049
  1443
  also have "... \<longlonglongrightarrow> cnj (rGamma z)" by (intro tendsto_cnj tendsto_intros)
eberlm@62049
  1444
  finally show ?thesis unfolding rGamma_complex_def by (intro sym[OF limI])
eberlm@62049
  1445
qed
eberlm@62049
  1446
eberlm@62049
  1447
lemma cnj_Gamma: "cnj (Gamma z) = Gamma (cnj z)"
eberlm@62049
  1448
  unfolding Gamma_def by (simp add: cnj_rGamma)
paulson@62131
  1449
paulson@62131
  1450
lemma Gamma_complex_real:
eberlm@62049
  1451
  "z \<in> \<real> \<Longrightarrow> Gamma z \<in> (\<real> :: complex set)" and rGamma_complex_real: "z \<in> \<real> \<Longrightarrow> rGamma z \<in> \<real>"
eberlm@62049
  1452
  by (simp_all add: Reals_cnj_iff cnj_Gamma cnj_rGamma)
eberlm@62049
  1453
lp15@62534
  1454
lemma field_differentiable_rGamma: "rGamma field_differentiable (at z within A)"
lp15@62534
  1455
  using has_field_derivative_rGamma[of z] unfolding field_differentiable_def by blast
eberlm@62049
  1456
eberlm@64969
  1457
lemma holomorphic_rGamma [holomorphic_intros]: "rGamma holomorphic_on A"
lp15@62534
  1458
  unfolding holomorphic_on_def by (auto intro!: field_differentiable_rGamma)
eberlm@62049
  1459
eberlm@68721
  1460
lemma holomorphic_rGamma' [holomorphic_intros]: 
eberlm@68721
  1461
  assumes "f holomorphic_on A"
eberlm@68721
  1462
  shows   "(\<lambda>x. rGamma (f x)) holomorphic_on A"
eberlm@68721
  1463
proof -
eberlm@68721
  1464
  have "rGamma \<circ> f holomorphic_on A" using assms
eberlm@68721
  1465
    by (intro holomorphic_on_compose assms holomorphic_rGamma)
eberlm@68721
  1466
  thus ?thesis by (simp only: o_def)
eberlm@68721
  1467
qed
eberlm@68721
  1468
eberlm@64969
  1469
lemma analytic_rGamma: "rGamma analytic_on A"
eberlm@64969
  1470
  unfolding analytic_on_def by (auto intro!: exI[of _ 1] holomorphic_rGamma)
eberlm@62049
  1471
eberlm@62049
  1472
lp15@62534
  1473
lemma field_differentiable_Gamma: "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Gamma field_differentiable (at z within A)"
lp15@62534
  1474
  using has_field_derivative_Gamma[of z] unfolding field_differentiable_def by auto
eberlm@62049
  1475
eberlm@64969
  1476
lemma holomorphic_Gamma [holomorphic_intros]: "A \<inter> \<int>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> Gamma holomorphic_on A"
lp15@62534
  1477
  unfolding holomorphic_on_def by (auto intro!: field_differentiable_Gamma)
eberlm@62049
  1478
eberlm@68721
  1479
lemma holomorphic_Gamma' [holomorphic_intros]: 
eberlm@68721
  1480
  assumes "f holomorphic_on A" and "\<And>x. x \<in> A \<Longrightarrow> f x \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@68721
  1481
  shows   "(\<lambda>x. Gamma (f x)) holomorphic_on A"
eberlm@68721
  1482
proof -
eberlm@68721
  1483
  have "Gamma \<circ> f holomorphic_on A" using assms
eberlm@68721
  1484
    by (intro holomorphic_on_compose assms holomorphic_Gamma) auto
eberlm@68721
  1485
  thus ?thesis by (simp only: o_def)
eberlm@68721
  1486
qed
eberlm@68721
  1487
eberlm@64969
  1488
lemma analytic_Gamma: "A \<inter> \<int>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> Gamma analytic_on A"
eberlm@62049
  1489
  by (rule analytic_on_subset[of _ "UNIV - \<int>\<^sub>\<le>\<^sub>0"], subst analytic_on_open)
eberlm@64969
  1490
     (auto intro!: holomorphic_Gamma)
lp15@66512
  1491
lp15@66512
  1492
lp15@66512
  1493
lemma field_differentiable_ln_Gamma_complex:
eberlm@64969
  1494
  "z \<notin> \<real>\<^sub>\<le>\<^sub>0 \<Longrightarrow> ln_Gamma field_differentiable (at (z::complex) within A)"
eberlm@64969
  1495
  by (rule field_differentiable_within_subset[of _ _ UNIV])
eberlm@64969
  1496
     (force simp: field_differentiable_def intro!: derivative_intros)+
eberlm@64969
  1497
eberlm@64969
  1498
lemma holomorphic_ln_Gamma [holomorphic_intros]: "A \<inter> \<real>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> ln_Gamma holomorphic_on A"
eberlm@64969
  1499
  unfolding holomorphic_on_def by (auto intro!: field_differentiable_ln_Gamma_complex)
eberlm@64969
  1500
eberlm@64969
  1501
lemma analytic_ln_Gamma: "A \<inter> \<real>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> ln_Gamma analytic_on A"
eberlm@64969
  1502
  by (rule analytic_on_subset[of _ "UNIV - \<real>\<^sub>\<le>\<^sub>0"], subst analytic_on_open)
eberlm@64969
  1503
     (auto intro!: holomorphic_ln_Gamma)
eberlm@64969
  1504
eberlm@62049
  1505
eberlm@62049
  1506
lemma has_field_derivative_rGamma_complex' [derivative_intros]:
paulson@62131
  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
eberlm@62049
  1508
        -rGamma z * Digamma z)) (at z within A)"
eberlm@62049
  1509
  using has_field_derivative_rGamma[of z] by (auto elim!: nonpos_Ints_cases')
eberlm@62049
  1510
eberlm@62049
  1511
declare has_field_derivative_rGamma_complex'[THEN DERIV_chain2, derivative_intros]
eberlm@62049
  1512
eberlm@62049
  1513
lp15@62534
  1514
lemma field_differentiable_Polygamma:
eberlm@64969
  1515
  fixes z :: complex
lp15@62533
  1516
  shows
lp15@62534
  1517
  "z \<notin> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Polygamma n field_differentiable (at z within A)"
lp15@62534
  1518
  using has_field_derivative_Polygamma[of z n] unfolding field_differentiable_def by auto
eberlm@62049
  1519
eberlm@64969
  1520
lemma holomorphic_on_Polygamma [holomorphic_intros]: "A \<inter> \<int>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> Polygamma n holomorphic_on A"
lp15@62534
  1521
  unfolding holomorphic_on_def by (auto intro!: field_differentiable_Polygamma)
eberlm@62049
  1522
eberlm@62049
  1523
lemma analytic_on_Polygamma: "A \<inter> \<int>\<^sub>\<le>\<^sub>0 = {} \<Longrightarrow> Polygamma n analytic_on A"
eberlm@62049
  1524
  by (rule analytic_on_subset[of _ "UNIV - \<int>\<^sub>\<le>\<^sub>0"], subst analytic_on_open)
eberlm@62049
  1525
     (auto intro!: holomorphic_on_Polygamma)
eberlm@62049
  1526
eberlm@62049
  1527
eberlm@62049
  1528
eberlm@68624
  1529
subsection%unimportant \<open>The real Gamma function\<close>
eberlm@62049
  1530
eberlm@62049
  1531
lemma rGamma_series_real:
eberlm@62049
  1532
  "eventually (\<lambda>n. rGamma_series x n = Re (rGamma_series (of_real x) n)) sequentially"
eberlm@62049
  1533
  using eventually_gt_at_top[of "0 :: nat"]
eberlm@62049
  1534
proof eventually_elim
eberlm@62049
  1535
  fix n :: nat assume n: "n > 0"
eberlm@62049
  1536
  have "Re (rGamma_series (of_real x) n) =
eberlm@62049
  1537
          Re (of_real (pochhammer x (Suc n)) / (fact n * exp (of_real (x * ln (real_of_nat n)))))"
eberlm@62049
  1538
    using n by (simp add: rGamma_series_def powr_def Ln_of_nat pochhammer_of_real)
eberlm@62049
  1539
  also from n have "\<dots> = Re (of_real ((pochhammer x (Suc n)) /
eberlm@62049
  1540
                              (fact n * (exp (x * ln (real_of_nat n))))))"
eberlm@62049
  1541
    by (subst exp_of_real) simp
eberlm@62049
  1542
  also from n have "\<dots> = rGamma_series x n"
eberlm@62049
  1543
    by (subst Re_complex_of_real) (simp add: rGamma_series_def powr_def)
eberlm@62049
  1544
  finally show "rGamma_series x n = Re (rGamma_series (of_real x) n)" ..
eberlm@62049
  1545
qed
eberlm@62049
  1546
eberlm@68624
  1547
instantiation%unimportant real :: Gamma
eberlm@62049
  1548
begin
eberlm@62049
  1549
eberlm@62049
  1550
definition "rGamma_real x = Re (rGamma (of_real x :: complex))"
eberlm@62049
  1551
eberlm@62049
  1552
instance proof
eberlm@62049
  1553
  fix x :: real
eberlm@62049
  1554
  have "rGamma x = Re (rGamma (of_real x))" by (simp add: rGamma_real_def)
eberlm@62049
  1555
  also have "of_real \<dots> = rGamma (of_real x :: complex)"
eberlm@62049
  1556
    by (intro of_real_Re rGamma_complex_real) simp_all
eberlm@62049
  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)
eberlm@62049
  1558
  also have "\<dots> \<longleftrightarrow> (\<exists>n. x = - of_nat n)" by (auto elim!: nonpos_Ints_cases')
eberlm@62049
  1559
  finally show "(rGamma x) = 0 \<longleftrightarrow> (\<exists>n. x = - real_of_nat n)" by simp
eberlm@62049
  1560
next
eberlm@62049
  1561
  fix x :: real assume "\<And>n. x \<noteq> - of_nat n"
wenzelm@63539
  1562
  hence x: "complex_of_real x \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
  1563
    by (subst of_real_in_nonpos_Ints_iff) (auto elim!: nonpos_Ints_cases')
wenzelm@63539
  1564
  then have "x \<noteq> 0" by auto
wenzelm@63539
  1565
  with x have "(rGamma has_field_derivative - rGamma x * Digamma x) (at x)"
paulson@62131
  1566
    by (fastforce intro!: derivative_eq_intros has_vector_derivative_real_complex
eberlm@62049
  1567
                  simp: Polygamma_of_real rGamma_real_def [abs_def])
eberlm@62049
  1568
  thus "let d = (THE d. (\<lambda>n. \<Sum>k<n. inverse (of_nat (Suc k)) - inverse (x + of_nat k))
paulson@62131
  1569
                       \<longlonglongrightarrow> d) - euler_mascheroni *\<^sub>R 1 in  (\<lambda>y. (rGamma y - rGamma x +
eberlm@62049
  1570
              rGamma x * d * (y - x)) /\<^sub>R  norm (y - x)) \<midarrow>x\<rightarrow> 0"
eberlm@62049
  1571
      by (simp add: has_field_derivative_def has_derivative_def Digamma_def sums_def [abs_def]
eberlm@62049
  1572
                    netlimit_at of_real_def[symmetric] suminf_def)
eberlm@62049
  1573
next
eberlm@62049
  1574
  fix n :: nat
eberlm@62049
  1575
  have "(rGamma has_field_derivative (-1)^n * fact n) (at (- of_nat n :: real))"
paulson@62131
  1576
    by (fastforce intro!: derivative_eq_intros has_vector_derivative_real_complex
eberlm@62049
  1577
                  simp: Polygamma_of_real rGamma_real_def [abs_def])
nipkow@64272
  1578
  thus "let x = - of_nat n in (\<lambda>y. (rGamma y - rGamma x - (- 1) ^ n * prod of_nat {1..n} *
eberlm@62049
  1579
                  (y - x)) /\<^sub>R norm (y - x)) \<midarrow>x::real\<rightarrow> 0"
nipkow@64272
  1580
    by (simp add: has_field_derivative_def has_derivative_def fact_prod netlimit_at Let_def)
eberlm@62049
  1581
next
eberlm@62049
  1582
  fix x :: real
eberlm@62049
  1583
  have "rGamma_series x \<longlonglongrightarrow> rGamma x"
eberlm@62049
  1584
  proof (rule Lim_transform_eventually)
eberlm@62049
  1585
    show "(\<lambda>n. Re (rGamma_series (of_real x) n)) \<longlonglongrightarrow> rGamma x" unfolding rGamma_real_def
eberlm@62049
  1586
      by (intro tendsto_intros)
eberlm@62049
  1587
  qed (insert rGamma_series_real, simp add: eq_commute)
nipkow@64272
  1588
  thus "let fact' = \<lambda>n. prod of_nat {1..n};
eberlm@62049
  1589
            exp = \<lambda>x. THE e. (\<lambda>n. \<Sum>k<n. x ^ k /\<^sub>R fact k) \<longlonglongrightarrow> e;
eberlm@62049
  1590
            pochhammer' = \<lambda>a n. \<Prod>n = 0..n. a + of_nat n
eberlm@62049
  1591
        in  (\<lambda>n. pochhammer' x n / (fact' n * exp (x * ln (real_of_nat n) *\<^sub>R 1))) \<longlonglongrightarrow> rGamma x"
nipkow@64272
  1592
    by (simp add: fact_prod pochhammer_Suc_prod rGamma_series_def [abs_def] exp_def
haftmann@63367
  1593
                  of_real_def [symmetric] suminf_def sums_def [abs_def] atLeast0AtMost)
eberlm@62049
  1594
qed
eberlm@62049
  1595
eberlm@62049
  1596
end
eberlm@62049
  1597
eberlm@62049
  1598
eberlm@62049
  1599
lemma rGamma_complex_of_real: "rGamma (complex_of_real x) = complex_of_real (rGamma x)"
eberlm@62049
  1600
  unfolding rGamma_real_def using rGamma_complex_real by simp
eberlm@62049
  1601
eberlm@62049
  1602
lemma Gamma_complex_of_real: "Gamma (complex_of_real x) = complex_of_real (Gamma x)"
eberlm@62049
  1603
  unfolding Gamma_def by (simp add: rGamma_complex_of_real)
eberlm@62049
  1604
eberlm@62049
  1605
lemma rGamma_real_altdef: "rGamma x = lim (rGamma_series (x :: real))"
eberlm@62049
  1606
  by (rule sym, rule limI, rule tendsto_intros)
eberlm@62049
  1607
eberlm@62049
  1608
lemma Gamma_real_altdef1: "Gamma x = lim (Gamma_series (x :: real))"
eberlm@62049
  1609
  by (rule sym, rule limI, rule tendsto_intros)
eberlm@62049
  1610
eberlm@62049
  1611
lemma Gamma_real_altdef2: "Gamma x = Re (Gamma (of_real x))"
eberlm@62049
  1612
  using rGamma_complex_real[OF Reals_of_real[of x]]
paulson@62131
  1613
  by (elim Reals_cases)
eberlm@62049
  1614
     (simp only: Gamma_def rGamma_real_def of_real_inverse[symmetric] Re_complex_of_real)
eberlm@62049
  1615
paulson@62131
  1616
lemma ln_Gamma_series_complex_of_real:
eberlm@62049
  1617
  "x > 0 \<Longrightarrow> n > 0 \<Longrightarrow> ln_Gamma_series (complex_of_real x) n = of_real (ln_Gamma_series x n)"
eberlm@62049
  1618
proof -
eberlm@62049
  1619
  assume xn: "x > 0" "n > 0"
eberlm@62049
  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
eberlm@62049
  1621
    using that xn by (subst Ln_of_real [symmetric]) (auto intro!: add_nonneg_pos simp: field_simps)
eberlm@62049
  1622
  with xn show ?thesis by (simp add: ln_Gamma_series_def Ln_of_nat Ln_of_real)
eberlm@62049
  1623
qed
eberlm@62049
  1624
paulson@62131
  1625
lemma ln_Gamma_real_converges:
paulson@62131
  1626
  assumes "(x::real) > 0"
eberlm@62049
  1627
  shows   "convergent (ln_Gamma_series x)"
eberlm@62049
  1628
proof -
eberlm@62049
  1629
  have "(\<lambda>n. ln_Gamma_series (complex_of_real x) n) \<longlonglongrightarrow> ln_Gamma (of_real x)" using assms
eberlm@62049
  1630
    by (intro ln_Gamma_complex_LIMSEQ) (auto simp: of_real_in_nonpos_Ints_iff)
eberlm@62049
  1631
  moreover from eventually_gt_at_top[of "0::nat"]
paulson@62131
  1632
    have "eventually (\<lambda>n. complex_of_real (ln_Gamma_series x n) =
paulson@62131
  1633
            ln_Gamma_series (complex_of_real x) n) sequentially"
eberlm@62049
  1634
    by eventually_elim (simp add: ln_Gamma_series_complex_of_real assms)
eberlm@62049
  1635
  ultimately have "(\<lambda>n. complex_of_real (ln_Gamma_series x n)) \<longlonglongrightarrow> ln_Gamma (of_real x)"
eberlm@62049
  1636
    by (subst tendsto_cong) assumption+
eberlm@62049
  1637
  from tendsto_Re[OF this] show ?thesis by (auto simp: convergent_def)
eberlm@62049
  1638
qed
eberlm@62049
  1639
eberlm@62049
  1640
lemma ln_Gamma_real_LIMSEQ: "(x::real) > 0 \<Longrightarrow> ln_Gamma_series x \<longlonglongrightarrow> ln_Gamma x"
eberlm@62049
  1641
  using ln_Gamma_real_converges[of x] unfolding ln_Gamma_def by (simp add: convergent_LIMSEQ_iff)
eberlm@62049
  1642
eberlm@62049
  1643
lemma ln_Gamma_complex_of_real: "x > 0 \<Longrightarrow> ln_Gamma (complex_of_real x) = of_real (ln_Gamma x)"
eberlm@62049
  1644
proof (unfold ln_Gamma_def, rule limI, rule Lim_transform_eventually)
eberlm@62049
  1645
  assume x: "x > 0"
paulson@62131
  1646
  show "eventually (\<lambda>n. of_real (ln_Gamma_series x n) =
paulson@62131
  1647
            ln_Gamma_series (complex_of_real x) n) sequentially"
paulson@62131
  1648
    using eventually_gt_at_top[of "0::nat"]
eberlm@62049
  1649
    by eventually_elim (simp add: ln_Gamma_series_complex_of_real x)
eberlm@62049
  1650
qed (intro tendsto_of_real, insert ln_Gamma_real_LIMSEQ[of x], simp add: ln_Gamma_def)
paulson@62131
  1651
eberlm@62049
  1652
lemma Gamma_real_pos_exp: "x > (0 :: real) \<Longrightarrow> Gamma x = exp (ln_Gamma x)"
paulson@62131
  1653
  by (auto simp: Gamma_real_altdef2 Gamma_complex_altdef of_real_in_nonpos_Ints_iff
eberlm@62049
  1654
                 ln_Gamma_complex_of_real exp_of_real)
eberlm@62049
  1655
eberlm@62049
  1656
lemma ln_Gamma_real_pos: "x > 0 \<Longrightarrow> ln_Gamma x = ln (Gamma x :: real)"
eberlm@62049
  1657
  unfolding Gamma_real_pos_exp by simp
eberlm@62049
  1658
eberlm@64969
  1659
lemma ln_Gamma_complex_conv_fact: "n > 0 \<Longrightarrow> ln_Gamma (of_nat n :: complex) = ln (fact (n - 1))"
eberlm@64969
  1660
  using ln_Gamma_complex_of_real[of "real n"] Gamma_fact[of "n - 1", where 'a = real]
lp15@66512
  1661
  by (simp add: ln_Gamma_real_pos of_nat_diff Ln_of_real [symmetric])
lp15@66512
  1662
eberlm@64969
  1663
lemma ln_Gamma_real_conv_fact: "n > 0 \<Longrightarrow> ln_Gamma (real n) = ln (fact (n - 1))"
eberlm@64969
  1664
  using Gamma_fact[of "n - 1", where 'a = real]
eberlm@64969
  1665
  by (simp add: ln_Gamma_real_pos of_nat_diff Ln_of_real [symmetric])
eberlm@64969
  1666
eberlm@67278
  1667
lemma Gamma_real_pos [simp, intro]: "x > (0::real) \<Longrightarrow> Gamma x > 0"
eberlm@67278
  1668
  by (simp add: Gamma_real_pos_exp)
eberlm@67278
  1669
eberlm@67278
  1670
lemma Gamma_real_nonneg [simp, intro]: "x > (0::real) \<Longrightarrow> Gamma x \<ge> 0"
eberlm@62049
  1671
  by (simp add: Gamma_real_pos_exp)
paulson@62131
  1672
eberlm@62049
  1673
lemma has_field_derivative_ln_Gamma_real [derivative_intros]:
eberlm@62049
  1674
  assumes x: "x > (0::real)"
eberlm@62049
  1675
  shows "(ln_Gamma has_field_derivative Digamma x) (at x)"
eberlm@62049
  1676
proof (subst DERIV_cong_ev[OF refl _ refl])
eberlm@62049
  1677
  from assms show "((Re \<circ> ln_Gamma \<circ> complex_of_real) has_field_derivative Digamma x) (at x)"
eberlm@62049
  1678
    by (auto intro!: derivative_eq_intros has_vector_derivative_real_complex
eberlm@62049
  1679
             simp: Polygamma_of_real o_def)
eberlm@62049
  1680
  from eventually_nhds_in_nhd[of x "{0<..}"] assms
eberlm@62049
  1681
    show "eventually (\<lambda>y. ln_Gamma y = (Re \<circ> ln_Gamma \<circ> of_real) y) (nhds x)"
eberlm@62049
  1682
    by (auto elim!: eventually_mono simp: ln_Gamma_complex_of_real interior_open)
eberlm@62049
  1683
qed
lp15@66512
  1684
lp15@66512
  1685
lemma field_differentiable_ln_Gamma_real:
eberlm@64969
  1686
  "x > 0 \<Longrightarrow> ln_Gamma field_differentiable (at (x::real) within A)"
eberlm@64969
  1687
  by (rule field_differentiable_within_subset[of _ _ UNIV])
eberlm@64969
  1688
     (auto simp: field_differentiable_def intro!: derivative_intros)+
eberlm@62049
  1689
eberlm@62049
  1690
declare has_field_derivative_ln_Gamma_real[THEN DERIV_chain2, derivative_intros]
lp15@66512
  1691
eberlm@64969
  1692
lemma deriv_ln_Gamma_real:
eberlm@64969
  1693
  assumes "z > 0"
eberlm@64969
  1694
  shows   "deriv ln_Gamma z = Digamma (z :: real)"
eberlm@64969
  1695
  by (intro DERIV_imp_deriv has_field_derivative_ln_Gamma_real assms)
eberlm@62049
  1696
eberlm@62049
  1697
eberlm@62049
  1698
lemma has_field_derivative_rGamma_real' [derivative_intros]:
paulson@62131
  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
eberlm@62049
  1700
        -rGamma x * Digamma x)) (at x within A)"
eberlm@62049
  1701
  using has_field_derivative_rGamma[of x] by (force elim!: nonpos_Ints_cases')
eberlm@62049
  1702
eberlm@62049
  1703
declare has_field_derivative_rGamma_real'[THEN DERIV_chain2, derivative_intros]
eberlm@62049
  1704
eberlm@62049
  1705
lemma Polygamma_real_odd_pos:
eberlm@62049
  1706
  assumes "(x::real) \<notin> \<int>\<^sub>\<le>\<^sub>0" "odd n"
eberlm@62049
  1707
  shows   "Polygamma n x > 0"
eberlm@62049
  1708
proof -
eberlm@62049
  1709
  from assms have "x \<noteq> 0" by auto
eberlm@62049
  1710
  with assms show ?thesis
eberlm@62049
  1711
    unfolding Polygamma_def using Polygamma_converges'[of x "Suc n"]
paulson@62131
  1712
    by (auto simp: zero_less_power_eq simp del: power_Suc
eberlm@62049
  1713
             dest: plus_of_nat_eq_0_imp intro!: mult_pos_pos suminf_pos)
eberlm@62049
  1714
qed
eberlm@62049
  1715
eberlm@62049
  1716
lemma Polygamma_real_even_neg:
eberlm@62049
  1717
  assumes "(x::real) > 0" "n > 0" "even n"
eberlm@62049
  1718
  shows   "Polygamma n x < 0"
eberlm@62049
  1719
  using assms unfolding Polygamma_def using Polygamma_converges'[of x "Suc n"]
eberlm@62049
  1720
  by (auto intro!: mult_pos_pos suminf_pos)
paulson@62131
  1721
eberlm@62049
  1722
lemma Polygamma_real_strict_mono:
eberlm@62049
  1723
  assumes "x > 0" "x < (y::real)" "even n"
eberlm@62049
  1724
  shows   "Polygamma n x < Polygamma n y"
eberlm@62049
  1725
proof -
eberlm@62049
  1726
  have "\<exists>\<xi>. x < \<xi> \<and> \<xi> < y \<and> Polygamma n y - Polygamma n x = (y - x) * Polygamma (Suc n) \<xi>"
eberlm@62049
  1727
    using assms by (intro MVT2 derivative_intros impI allI) (auto elim!: nonpos_Ints_cases)
eberlm@62049
  1728
  then guess \<xi> by (elim exE conjE) note \<xi> = this
eberlm@62049
  1729
  note \<xi>(3)
paulson@62131
  1730
  also from \<xi>(1,2) assms have "(y - x) * Polygamma (Suc n) \<xi> > 0"
eberlm@62049
  1731
    by (intro mult_pos_pos Polygamma_real_odd_pos) (auto elim!: nonpos_Ints_cases)
eberlm@62049
  1732
  finally show ?thesis by simp
eberlm@62049
  1733
qed
eberlm@62049
  1734
eberlm@62049
  1735
lemma Polygamma_real_strict_antimono:
eberlm@62049
  1736
  assumes "x > 0" "x < (y::real)" "odd n"
eberlm@62049
  1737
  shows   "Polygamma n x > Polygamma n y"
eberlm@62049
  1738
proof -
eberlm@62049
  1739
  have "\<exists>\<xi>. x < \<xi> \<and> \<xi> < y \<and> Polygamma n y - Polygamma n x = (y - x) * Polygamma (Suc n) \<xi>"
eberlm@62049
  1740
    using assms by (intro MVT2 derivative_intros impI allI) (auto elim!: nonpos_Ints_cases)
eberlm@62049
  1741
  then guess \<xi> by (elim exE conjE) note \<xi> = this
eberlm@62049
  1742
  note \<xi>(3)
paulson@62131
  1743
  also from \<xi>(1,2) assms have "(y - x) * Polygamma (Suc n) \<xi> < 0"
eberlm@62049
  1744
    by (intro mult_pos_neg Polygamma_real_even_neg) simp_all
eberlm@62049
  1745
  finally show ?thesis by simp
eberlm@62049
  1746
qed
eberlm@62049
  1747
eberlm@62049
  1748
lemma Polygamma_real_mono:
eberlm@62049
  1749
  assumes "x > 0" "x \<le> (y::real)" "even n"
eberlm@62049
  1750
  shows   "Polygamma n x \<le> Polygamma n y"
paulson@62131
  1751
  using Polygamma_real_strict_mono[OF assms(1) _ assms(3), of y] assms(2)
eberlm@62049
  1752
  by (cases "x = y") simp_all
eberlm@62049
  1753
eberlm@63721
  1754
lemma Digamma_real_strict_mono: "(0::real) < x \<Longrightarrow> x < y \<Longrightarrow> Digamma x < Digamma y"
eberlm@63721
  1755
  by (rule Polygamma_real_strict_mono) simp_all
eberlm@63721
  1756
eberlm@63721
  1757
lemma Digamma_real_mono: "(0::real) < x \<Longrightarrow> x \<le> y \<Longrightarrow> Digamma x \<le> Digamma y"
eberlm@63721
  1758
  by (rule Polygamma_real_mono) simp_all
eberlm@63721
  1759
eberlm@62049
  1760
lemma Digamma_real_ge_three_halves_pos:
eberlm@62049
  1761
  assumes "x \<ge> 3/2"
eberlm@62049
  1762
  shows   "Digamma (x :: real) > 0"
eberlm@62049
  1763
proof -
eberlm@62049
  1764
  have "0 < Digamma (3/2 :: real)" by (fact Digamma_real_three_halves_pos)
eberlm@62049
  1765
  also from assms have "\<dots> \<le> Digamma x" by (intro Polygamma_real_mono) simp_all
eberlm@62049
  1766
  finally show ?thesis .
eberlm@62049
  1767
qed
eberlm@62049
  1768
eberlm@62049
  1769
lemma ln_Gamma_real_strict_mono:
eberlm@62049
  1770
  assumes "x \<ge> 3/2" "x < y"
eberlm@62049
  1771
  shows   "ln_Gamma (x :: real) < ln_Gamma y"
eberlm@62049
  1772
proof -
eberlm@62049
  1773
  have "\<exists>\<xi>. x < \<xi> \<and> \<xi> < y \<and> ln_Gamma y - ln_Gamma x = (y - x) * Digamma \<xi>"
eberlm@62049
  1774
    using assms by (intro MVT2 derivative_intros impI allI) (auto elim!: nonpos_Ints_cases)
eberlm@62049
  1775
  then guess \<xi> by (elim exE conjE) note \<xi> = this
eberlm@62049
  1776
  note \<xi>(3)
paulson@62131
  1777
  also from \<xi>(1,2) assms have "(y - x) * Digamma \<xi> > 0"
eberlm@62049
  1778
    by (intro mult_pos_pos Digamma_real_ge_three_halves_pos) simp_all
eberlm@62049
  1779
  finally show ?thesis by simp
paulson@62131
  1780
qed
paulson@62131
  1781
eberlm@62049
  1782
lemma Gamma_real_strict_mono:
eberlm@62049
  1783
  assumes "x \<ge> 3/2" "x < y"
eberlm@62049
  1784
  shows   "Gamma (x :: real) < Gamma y"
eberlm@62049
  1785
proof -
eberlm@62049
  1786
  from Gamma_real_pos_exp[of x] assms have "Gamma x = exp (ln_Gamma x)" by simp
eberlm@62049
  1787
  also have "\<dots> < exp (ln_Gamma y)" by (intro exp_less_mono ln_Gamma_real_strict_mono assms)
eberlm@62049
  1788
  also from Gamma_real_pos_exp[of y] assms have "\<dots> = Gamma y" by simp
eberlm@62049
  1789
  finally show ?thesis .
eberlm@62049
  1790
qed
eberlm@62049
  1791
eberlm@68624
  1792
theorem log_convex_Gamma_real: "convex_on {0<..} (ln \<circ> Gamma :: real \<Rightarrow> real)"
eberlm@62049
  1793
  by (rule convex_on_realI[of _ _ Digamma])
paulson@62131
  1794
     (auto intro!: derivative_eq_intros Polygamma_real_mono Gamma_real_pos
eberlm@62049
  1795
           simp: o_def Gamma_eq_zero_iff elim!: nonpos_Ints_cases')
eberlm@62049
  1796
eberlm@62049
  1797
eberlm@63725
  1798
subsection \<open>The uniqueness of the real Gamma function\<close>
eberlm@63725
  1799
eberlm@63725
  1800
text \<open>
lp15@66512
  1801
  The following is a proof of the Bohr--Mollerup theorem, which states that
wenzelm@69566
  1802
  any log-convex function \<open>G\<close> on the positive reals that fulfils \<open>G(1) = 1\<close> and
wenzelm@69566
  1803
  satisfies the functional equation \<open>G(x + 1) = x G(x)\<close> must be equal to the
eberlm@63725
  1804
  Gamma function.
wenzelm@69566
  1805
  In principle, if \<open>G\<close> is a holomorphic complex function, one could then extend
lp15@66512
  1806
  this from the positive reals to the entire complex plane (minus the non-positive
eberlm@63725
  1807
  integers, where the Gamma function is not defined).
eberlm@63725
  1808
\<close>
eberlm@63725
  1809
eberlm@68624
  1810
context%unimportant
eberlm@63725
  1811
  fixes G :: "real \<Rightarrow> real"
eberlm@63725
  1812
  assumes G_1: "G 1 = 1"
eberlm@63725
  1813
  assumes G_plus1: "x > 0 \<Longrightarrow> G (x + 1) = x * G x"
eberlm@63725
  1814
  assumes G_pos: "x > 0 \<Longrightarrow> G x > 0"
eberlm@63725
  1815
  assumes log_convex_G: "convex_on {0<..} (ln \<circ> G)"
eberlm@63725
  1816
begin
eberlm@63725
  1817
eberlm@63725
  1818
private lemma G_fact: "G (of_nat n + 1) = fact n"
eberlm@63725
  1819
  using G_plus1[of "real n + 1" for n]
eberlm@63725
  1820
  by (induction n) (simp_all add: G_1 G_plus1)
eberlm@63725
  1821
eberlm@63725
  1822
private definition S :: "real \<Rightarrow> real \<Rightarrow> real" where
eberlm@63725
  1823
  "S x y = (ln (G y) - ln (G x)) / (y - x)"
eberlm@63725
  1824
lp15@66512
  1825
private lemma S_eq:
eberlm@63725
  1826
  "n \<ge> 2 \<Longrightarrow> S (of_nat n) (of_nat n + x) = (ln (G (real n + x)) - ln (fact (n - 1))) / x"
eberlm@63725
  1827
  by (subst G_fact [symmetric]) (simp add: S_def add_ac of_nat_diff)
eberlm@63725
  1828
eberlm@63725
  1829
private lemma G_lower:
eberlm@63725
  1830
  assumes x: "x > 0" and n: "n \<ge> 1"
eberlm@63725
  1831
  shows  "Gamma_series x n \<le> G x"
eberlm@63725
  1832
proof -
eberlm@63725
  1833
  have "(ln \<circ> G) (real (Suc n)) \<le> ((ln \<circ> G) (real (Suc n) + x) -
eberlm@63725
  1834
          (ln \<circ> G) (real (Suc n) - 1)) / (real (Suc n) + x - (real (Suc n) - 1)) *
eberlm@63725
  1835
           (real (Suc n) - (real (Suc n) - 1)) + (ln \<circ> G) (real (Suc n) - 1)"
eberlm@63725
  1836
    using x n by (intro convex_onD_Icc' convex_on_subset[OF log_convex_G]) auto
lp15@66512
  1837
  hence "S (of_nat n) (of_nat (Suc n)) \<le> S (of_nat (Suc n)) (of_nat (Suc n) + x)"
eberlm@63725
  1838
    unfolding S_def using x by (simp add: field_simps)
eberlm@63725
  1839
  also have "S (of_nat n) (of_nat (Suc n)) = ln (fact n) - ln (fact (n-1))"
lp15@66512
  1840
    unfolding S_def using n
eberlm@63725
  1841
    by (subst (1 2) G_fact [symmetric]) (simp_all add: add_ac of_nat_diff)
eberlm@63725
  1842
  also have "\<dots> = ln (fact n / fact (n-1))" by (subst ln_div) simp_all
eberlm@63725
  1843
  also from n have "fact n / fact (n - 1) = n" by (cases n) simp_all
eberlm@63725
  1844
  finally have "x * ln (real n) + ln (fact n) \<le> ln (G (real (Suc n) + x))"
eberlm@63725
  1845
    using x n by (subst (asm) S_eq) (simp_all add: field_simps)
lp15@66512
  1846
  also have "x * ln (real n) + ln (fact n) = ln (exp (x * ln (real n)) * fact n)"
eberlm@63725
  1847
    using x by (simp add: ln_mult)
eberlm@63725
  1848
  finally have "exp (x * ln (real n)) * fact n \<le> G (real (Suc n) + x)" using x
eberlm@63725
  1849
    by (subst (asm) ln_le_cancel_iff) (simp_all add: G_pos)
eberlm@63725
  1850
  also have "G (real (Suc n) + x) = pochhammer x (Suc n) * G x"
eberlm@63725
  1851
    using G_plus1[of "real (Suc n) + x" for n] G_plus1[of x] x
eberlm@63725
  1852
    by (induction n) (simp_all add: pochhammer_Suc add_ac)
eberlm@63725
  1853
  finally show "Gamma_series x n \<le> G x"
eberlm@63725
  1854
    using x by (simp add: field_simps pochhammer_pos Gamma_series_def)
eberlm@63725
  1855
qed
eberlm@63725
  1856
eberlm@63725
  1857
private lemma G_upper:
eberlm@63725
  1858
  assumes x: "x > 0" "x \<le> 1" and n: "n \<ge> 2"
eberlm@63725
  1859
  shows  "G x \<le> Gamma_series x n * (1 + x / real n)"
eberlm@63725
  1860
proof -
eberlm@63725
  1861
  have "(ln \<circ> G) (real n + x) \<le> ((ln \<circ> G) (real n + 1) -
eberlm@63725
  1862
          (ln \<circ> G) (real n)) / (real n + 1 - (real n)) *
lp15@66512
  1863
           ((real n + x) - real n) + (ln \<circ> G) (real n)"
eberlm@63725
  1864
    using x n by (intro convex_onD_Icc' convex_on_subset[OF log_convex_G]) auto
lp15@66512
  1865
  hence "S (of_nat n) (of_nat n + x) \<le> S (of_nat n) (of_nat n + 1)"
eberlm@63725
  1866
    unfolding S_def using x by (simp add: field_simps)
eberlm@63725
  1867
  also from n have "S (of_nat n) (of_nat n + 1) = ln (fact n) - ln (fact (n-1))"
eberlm@63725
  1868
    by (subst (1 2) G_fact [symmetric]) (simp add: S_def add_ac of_nat_diff)
eberlm@63725
  1869
  also have "\<dots> = ln (fact n / (fact (n-1)))" using n by (subst ln_div) simp_all
eberlm@63725
  1870
  also from n have "fact n / fact (n - 1) = n" by (cases n) simp_all
lp15@66512
  1871
  finally have "ln (G (real n + x)) \<le> x * ln (real n) + ln (fact (n - 1))"
eberlm@63725
  1872
    using x n by (subst (asm) S_eq) (simp_all add: field_simps)
eberlm@63725
  1873
  also have "\<dots> = ln (exp (x * ln (real n)) * fact (n - 1))" using x
eberlm@63725
  1874
    by (simp add: ln_mult)
eberlm@63725
  1875
  finally have "G (real n + x) \<le> exp (x * ln (real n)) * fact (n - 1)" using x
eberlm@63725
  1876
    by (subst (asm) ln_le_cancel_iff) (simp_all add: G_pos)
eberlm@63725
  1877
  also have "G (real n + x) = pochhammer x n * G x"
eberlm@63725
  1878
    using G_plus1[of "real n + x" for n] x
eberlm@63725
  1879
    by (induction n) (simp_all add: pochhammer_Suc add_ac)
eberlm@63725
  1880
  finally have "G x \<le> exp (x * ln (real n)) * fact (n- 1) / pochhammer x n"
eberlm@63725
  1881
    using x by (simp add: field_simps pochhammer_pos)
eberlm@63725
  1882
  also from n have "fact (n - 1) = fact n / n" by (cases n) simp_all
lp15@66512
  1883
  also have "exp (x * ln (real n)) * \<dots> / pochhammer x n =
eberlm@63725
  1884
               Gamma_series x n * (1 + x / real n)" using n x
eberlm@63725
  1885
    by (simp add: Gamma_series_def divide_simps pochhammer_Suc)
eberlm@63725
  1886
  finally show ?thesis .
eberlm@63725
  1887
qed
eberlm@63725
  1888
eberlm@63725
  1889
private lemma G_eq_Gamma_aux:
eberlm@63725
  1890
  assumes x: "x > 0" "x \<le> 1"
eberlm@63725
  1891
  shows   "G x = Gamma x"
eberlm@63725
  1892
proof (rule antisym)
eberlm@63725
  1893
  show "G x \<ge> Gamma x"
lp15@63952
  1894
  proof (rule tendsto_upperbound)
eberlm@63725
  1895
    from G_lower[of x] show "eventually (\<lambda>n. Gamma_series x n \<le> G x) sequentially"
lp15@65578
  1896
      using  x by (auto intro: eventually_mono[OF eventually_ge_at_top[of "1::nat"]])
eberlm@63725
  1897
  qed (simp_all add: Gamma_series_LIMSEQ)
eberlm@63725
  1898
next
eberlm@63725
  1899
  show "G x \<le> Gamma x"
lp15@63952
  1900
  proof (rule tendsto_lowerbound)
eberlm@63725
  1901
    have "(\<lambda>n. Gamma_series x n * (1 + x / real n)) \<longlonglongrightarrow> Gamma x * (1 + 0)"
lp15@66512
  1902
      by (rule tendsto_intros real_tendsto_divide_at_top
eberlm@63725
  1903
               Gamma_series_LIMSEQ filterlim_real_sequentially)+
eberlm@63725
  1904
    thus "(\<lambda>n. Gamma_series x n * (1 + x / real n)) \<longlonglongrightarrow> Gamma x" by simp
eberlm@63725
  1905
  next
eberlm@63725
  1906
    from G_upper[of x] show "eventually (\<lambda>n. Gamma_series x n * (1 + x / real n) \<ge> G x) sequentially"
lp15@65578
  1907
      using x by (auto intro: eventually_mono[OF eventually_ge_at_top[of "2::nat"]])
eberlm@63725
  1908
  qed simp_all
eberlm@63725
  1909
qed
eberlm@63725
  1910
eberlm@63725
  1911
theorem Gamma_pos_real_unique:
eberlm@63725
  1912
  assumes x: "x > 0"
eberlm@63725
  1913
  shows   "G x = Gamma x"
eberlm@63725
  1914
proof -
eberlm@63725
  1915
  have G_eq: "G (real n + x) = Gamma (real n + x)" if "x \<in> {0<..1}" for n x using that
eberlm@63725
  1916
  proof (induction n)
eberlm@63725
  1917
    case (Suc n)
eberlm@63725
  1918
    from Suc have "x + real n > 0" by simp
eberlm@63725
  1919
    hence "x + real n \<notin> \<int>\<^sub>\<le>\<^sub>0" by auto
eberlm@63725
  1920
    with Suc show ?case using G_plus1[of "real n + x"] Gamma_plus1[of "real n + x"]
eberlm@63725
  1921
      by (auto simp: add_ac)
eberlm@63725
  1922
  qed (simp_all add: G_eq_Gamma_aux)
lp15@66512
  1923
eberlm@63725
  1924
  show ?thesis
eberlm@63725
  1925
  proof (cases "frac x = 0")
eberlm@63725
  1926
    case True
eberlm@63725
  1927
    hence "x = of_int (floor x)" by (simp add: frac_def)
eberlm@63725
  1928
    with x have x_eq: "x = of_nat (nat (floor x) - 1) + 1" by simp
eberlm@63725
  1929
    show ?thesis by (subst (1 2) x_eq, rule G_eq) simp_all
eberlm@63725
  1930
  next
eberlm@63725
  1931
    case False
eberlm@63725
  1932
    from assms have x_eq: "x = of_nat (nat (floor x)) + frac x"
eberlm@63725
  1933
      by (simp add: frac_def)
eberlm@63725
  1934
    have frac_le_1: "frac x \<le> 1" unfolding frac_def by linarith
eberlm@63725
  1935
    show ?thesis
eberlm@63725
  1936
      by (subst (1 2) x_eq, rule G_eq, insert False frac_le_1) simp_all
eberlm@63725
  1937
  qed
eberlm@63725
  1938
qed
eberlm@63725
  1939
eberlm@63725
  1940
end
eberlm@63725
  1941
eberlm@63725
  1942
eberlm@68624
  1943
subsection \<open>The Beta function\<close>
eberlm@62049
  1944
eberlm@62049
  1945
definition Beta where "Beta a b = Gamma a * Gamma b / Gamma (a + b)"
eberlm@62049
  1946
eberlm@62049
  1947
lemma Beta_altdef: "Beta a b = Gamma a * Gamma b * rGamma (a + b)"
eberlm@62049
  1948
  by (simp add: inverse_eq_divide Beta_def Gamma_def)
eberlm@62049
  1949
eberlm@62049
  1950
lemma Beta_commute: "Beta a b = Beta b a"
eberlm@62049
  1951
  unfolding Beta_def by (simp add: ac_simps)
eberlm@62049
  1952
eberlm@62049
  1953
lemma has_field_derivative_Beta1 [derivative_intros]:
eberlm@62049
  1954
  assumes "x \<notin> \<int>\<^sub>\<le>\<^sub>0" "x + y \<notin> \<int>\<^sub>\<le>\<^sub>0"
paulson@62131
  1955
  shows   "((\<lambda>x. Beta x y) has_field_derivative (Beta x y * (Digamma x - Digamma (x + y))))
eberlm@62049
  1956
               (at x within A)" unfolding Beta_altdef
eberlm@62049
  1957
  by (rule DERIV_cong, (rule derivative_intros assms)+) (simp add: algebra_simps)
eberlm@62049
  1958
eberlm@63295
  1959
lemma Beta_pole1: "x \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Beta x y = 0"
eberlm@63295
  1960
  by (auto simp add: Beta_def elim!: nonpos_Ints_cases')
eberlm@63295
  1961
eberlm@63295
  1962
lemma Beta_pole2: "y \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Beta x y = 0"
eberlm@63295
  1963
  by (auto simp add: Beta_def elim!: nonpos_Ints_cases')
hoelzl@63594
  1964
eberlm@63295
  1965
lemma Beta_zero: "x + y \<in> \<int>\<^sub>\<le>\<^sub>0 \<Longrightarrow> Beta x y = 0"
eberlm@63295
  1966
  by (auto simp add: Beta_def elim!: nonpos_Ints_cases')
hoelzl@63594
  1967
eberlm@62049
  1968
lemma has_field_derivative_Beta2 [derivative_intros]:
eberlm@62049
  1969
  assumes "y \<notin> \<int>\<^sub>\<le>\<^sub>0" "x + y \<notin> \<int>\<^sub>\<le>\<^sub>0"
paulson@62131
  1970
  shows   "((\<lambda>y. Beta x y) has_field_derivative (Beta x y * (Digamma y - Digamma (x + y))))
eberlm@62049
  1971
               (at y within A)"
eberlm@62049
  1972
  using has_field_derivative_Beta1[of y x A] assms by (simp add: Beta_commute add_ac)
eberlm@62049
  1973
eberlm@68624
  1974
theorem Beta_plus1_plus1:
paulson@62131
  1975
  assumes "x \<notin> \<int>\<^sub>\<le>\<^sub>0" "y \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
  1976
  shows   "Beta (x + 1) y + Beta x (y + 1) = Beta x y"
eberlm@62049
  1977
proof -
paulson@62131
  1978
  have "Beta (x + 1) y + Beta x (y + 1) =
eberlm@62049
  1979
            (Gamma (x + 1) * Gamma y + Gamma x * Gamma (y + 1)) * rGamma ((x + y) + 1)"
eberlm@62049
  1980
    by (simp add: Beta_altdef add_divide_distrib algebra_simps)
eberlm@62049
  1981
  also have "\<dots> = (Gamma x * Gamma y) * ((x + y) * rGamma ((x + y) + 1))"
eberlm@62049
  1982
    by (subst assms[THEN Gamma_plus1])+ (simp add: algebra_simps)
eberlm@62049
  1983
  also from assms have "\<dots> = Beta x y" unfolding Beta_altdef by (subst rGamma_plus1) simp
eberlm@62049
  1984
  finally show ?thesis .
eberlm@62049
  1985
qed
eberlm@62049
  1986
eberlm@68624
  1987
theorem Beta_plus1_left:
eberlm@63295
  1988
  assumes "x \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
  1989
  shows   "(x + y) * Beta (x + 1) y = x * Beta x y"
eberlm@62049
  1990
proof -
eberlm@62049
  1991
  have "(x + y) * Beta (x + 1) y = Gamma (x + 1) * Gamma y * ((x + y) * rGamma ((x + y) + 1))"
eberlm@62049
  1992
    unfolding Beta_altdef by (simp only: ac_simps)
eberlm@62049
  1993
  also have "\<dots> = x * Beta x y" unfolding Beta_altdef
eberlm@62049
  1994
     by (subst assms[THEN Gamma_plus1] rGamma_plus1)+ (simp only: ac_simps)
eberlm@62049
  1995
  finally show ?thesis .
eberlm@62049
  1996
qed
eberlm@62049
  1997
eberlm@68624
  1998
theorem Beta_plus1_right:
eberlm@63295
  1999
  assumes "y \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
  2000
  shows   "(x + y) * Beta x (y + 1) = y * Beta x y"
eberlm@63295
  2001
  using Beta_plus1_left[of y x] assms by (simp_all add: Beta_commute add.commute)
eberlm@62049
  2002
eberlm@62049
  2003
lemma Gamma_Gamma_Beta:
eberlm@63295
  2004
  assumes "x + y \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
  2005
  shows   "Gamma x * Gamma y = Beta x y * Gamma (x + y)"
eberlm@62049
  2006
  unfolding Beta_altdef using assms Gamma_eq_zero_iff[of "x+y"]
eberlm@62049
  2007
  by (simp add: rGamma_inverse_Gamma)
eberlm@62049
  2008
eberlm@62049
  2009
eberlm@62049
  2010
eberlm@62049
  2011
subsection \<open>Legendre duplication theorem\<close>
eberlm@62049
  2012
eberlm@62049
  2013
context
eberlm@62049
  2014
begin
eberlm@62049
  2015
eberlm@62049
  2016
private lemma Gamma_legendre_duplication_aux:
eberlm@62049
  2017
  fixes z :: "'a :: Gamma"
eberlm@62049
  2018
  assumes "z \<notin> \<int>\<^sub>\<le>\<^sub>0" "z + 1/2 \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
  2019
  shows "Gamma z * Gamma (z + 1/2) = exp ((1 - 2*z) * of_real (ln 2)) * Gamma (1/2) * Gamma (2*z)"
eberlm@62049
  2020
proof -
eberlm@62049
  2021
  let ?powr = "\<lambda>b a. exp (a * of_real (ln (of_nat b)))"
paulson@62131
  2022
  let ?h = "\<lambda>n. (fact (n-1))\<^sup>2 / fact (2*n-1) * of_nat (2^(2*n)) *
eberlm@62049
  2023
                exp (1/2 * of_real (ln (real_of_nat n)))"
eberlm@62049
  2024
  {
eberlm@62049
  2025
    fix z :: 'a assume z: "z \<notin> \<int>\<^sub>\<le>\<^sub>0" "z + 1/2 \<notin> \<int>\<^sub>\<le>\<^sub>0"
paulson@62131
  2026
    let ?g = "\<lambda>n. ?powr 2 (2*z) * Gamma_series' z n * Gamma_series' (z + 1/2) n /
eberlm@62049
  2027
                      Gamma_series' (2*z) (2*n)"
eberlm@62049
  2028
    have "eventually (\<lambda>n. ?g n = ?h n) sequentially" using eventually_gt_at_top
eberlm@62049
  2029
    proof eventually_elim
eberlm@62049
  2030
      fix n :: nat assume n: "n > 0"
eberlm@62049
  2031
      let ?f = "fact (n - 1) :: 'a" and ?f' = "fact (2*n - 1) :: 'a"
eberlm@62049
  2032
      have A: "exp t * exp t = exp (2*t :: 'a)" for t by (subst exp_add [symmetric]) simp
eberlm@62049
  2033
      have A: "Gamma_series' z n * Gamma_series' (z + 1/2) n = ?f^2 * ?powr n (2*z + 1/2) /
eberlm@62049
  2034
                (pochhammer z n * pochhammer (z + 1/2) n)"
eberlm@62049
  2035
        by (simp add: Gamma_series'_def exp_add ring_distribs power2_eq_square A mult_ac)
paulson@62131
  2036
      have B: "Gamma_series' (2*z) (2*n) =
paulson@62131
  2037
                       ?f' * ?powr 2 (2*z) * ?powr n (2*z) /
eberlm@62049
  2038
                       (of_nat (2^(2*n)) * pochhammer z n * pochhammer (z+1/2) n)" using n
eberlm@62049
  2039
        by (simp add: Gamma_series'_def ln_mult exp_add ring_distribs pochhammer_double)
eberlm@62049
  2040
      from z have "pochhammer z n \<noteq> 0" by (auto dest: pochhammer_eq_0_imp_nonpos_Int)
eberlm@62049
  2041
      moreover from z have "pochhammer (z + 1/2) n \<noteq> 0" by (auto dest: pochhammer_eq_0_imp_nonpos_Int)
paulson@62131
  2042
      ultimately have "?powr 2 (2*z) * (Gamma_series' z n * Gamma_series' (z + 1/2) n) / Gamma_series' (2*z) (2*n) =
eberlm@62049
  2043
         ?f^2 / ?f' * of_nat (2^(2*n)) * (?powr n ((4*z + 1)/2) * ?powr n (-2*z))"
eberlm@62049
  2044
        using n unfolding A B by (simp add: divide_simps exp_minus)
eberlm@62049
  2045
      also have "?powr n ((4*z + 1)/2) * ?powr n (-2*z) = ?powr n (1/2)"
eberlm@62049
  2046
        by (simp add: algebra_simps exp_add[symmetric] add_divide_distrib)
eberlm@62049
  2047
      finally show "?g n = ?h n" by (simp only: mult_ac)
eberlm@62049
  2048
    qed
paulson@62131
  2049
paulson@62131
  2050
    moreover from z double_in_nonpos_Ints_imp[of z] have "2 * z \<notin> \<int>\<^sub>\<le>\<^sub>0" by auto
eberlm@62049
  2051
    hence "?g \<longlonglongrightarrow> ?powr 2 (2*z) * Gamma z * Gamma (z+1/2) / Gamma (2*z)"
nipkow@69064
  2052
      using LIMSEQ_subseq_LIMSEQ[OF Gamma_series'_LIMSEQ, of "(*)2" "2*z"]
eberlm@62049
  2053
      by (intro tendsto_intros Gamma_series'_LIMSEQ)
eberlm@66447
  2054
         (simp_all add: o_def strict_mono_def Gamma_eq_zero_iff)
eberlm@62049
  2055
    ultimately have "?h \<longlonglongrightarrow> ?powr 2 (2*z) * Gamma z * Gamma (z+1/2) / Gamma (2*z)"
eberlm@62049
  2056
      by (rule Lim_transform_eventually)
eberlm@62049
  2057
  } note lim = this
paulson@62131
  2058
paulson@62131
  2059
  from assms double_in_nonpos_Ints_imp[of z] have z': "2 * z \<notin> \<int>\<^sub>\<le>\<^sub>0" by auto
eberlm@62049
  2060
  from fraction_not_in_ints[of 2 1] have "(1/2 :: 'a) \<notin> \<int>\<^sub>\<le>\<^sub>0"
eberlm@62049
  2061
    by (intro not_in_Ints_imp_not_in_nonpos_Ints) simp_all
lp15@67976
  2062
  with lim[of "1/2 :: 'a"] have "?h \<longlonglongrightarrow> 2 * Gamma (1/2 :: 'a)" by (simp add: exp_of_real)
eberlm@62049
  2063
  from LIMSEQ_unique[OF this lim[OF assms]] z' show ?thesis
eberlm@62049
  2064
    by (simp add: divide_simps Gamma_eq_zero_iff ring_distribs exp_diff exp_of_real ac_simps)
eberlm@62049
  2065
qed
eberlm@62049
  2066
eberlm@68624
  2067
theorem Gamma_reflection_complex:
eberlm@62049
  2068
  fixes z :: complex
eberlm@62049
  2069
  shows "Gamma z * Gamma (1 - z) = of_real pi / sin (of_real pi * z)"
eberlm@62049
  2070
proof -
eberlm@62049
  2071
  let ?g = "\<lambda>z::complex. Gamma z * Gamma (1 - z) * sin (of_real pi * z)"
wenzelm@63040
  2072
  define g where [abs_def]: "g z = (if z \<in> \<int> then of_real pi else ?g z)" for z :: complex
eberlm@62049
  2073
  let ?h = "\<lambda>z::complex. (of_real pi * cot (of_real pi*z) + Digamma z - Digamma (1 - z))"
wenzelm@63040
  2074
  define h where [abs_def]: "h z = (if z \<in> \<int> then 0 else ?h z)" for z :: complex
eberlm@62049
  2075
wenzelm@69597
  2076
  \<comment> \<open>\<^term>\<open>g\<close> is periodic with period 1.\<close>
eberlm@62049
  2077
  interpret g: periodic_fun_simple' g
eberlm@62049
  2078
  proof
eberlm@62049
  2079
    fix z :: complex
eberlm@62049
  2080
    show "g (z + 1) = g z"
eberlm@62049
  2081
    proof (cases "z \<in> \<int>")
eberlm@62049
  2082
      case False
eberlm@62049
  2083
      hence "z * g z = z * Beta z (- z + 1) * sin (of_real pi * z)" by (simp add: g_def Beta_def)
paulson@62131
  2084
      also have "z * Beta z (- z + 1) = (z + 1 + -z) * Beta (z + 1) (- z + 1)"
eberlm@62049
  2085
        using False Ints_diff[of 1 "1 - z"] nonpos_Ints_subset_Ints
eberlm@62049
  2086
        by (subst Beta_plus1_left [symmetric]) auto
eberlm@62049
  2087
      also have "\<dots> * sin (of_real pi * z) = z * (Beta (z + 1) (-z) * sin (of_real pi * (z + 1)))"
eberlm@62049
  2088
        using False Ints_diff[of "z+1" 1] Ints_minus[of "-z"] nonpos_Ints_subset_Ints
eberlm@62049
  2089
        by (subst Beta_plus1_right) (auto simp: ring_distribs sin_plus_pi)
paulson@62131
  2090
      also from False have "Beta (z + 1) (-z) * sin (of_real pi * (z + 1)) = g (z + 1)"
eberlm@62049
  2091
        using Ints_diff[of "z+1" 1] by (auto simp: g_def Beta_def)
eberlm@62049
  2092
      finally show "g (z + 1) = g z" using False by (subst (asm) mult_left_cancel) auto
eberlm@62049
  2093
    qed (simp add: g_def)
eberlm@62049
  2094
  qed
eberlm@62049
  2095
wenzelm@69597
  2096
  \<comment> \<open>\<^term>\<open>g\<close> is entire.\<close>
eberlm@64969
  2097
  have g_g' [derivative_intros]: "(g has_field_derivative (h z * g z)) (at z)" for z :: complex
eberlm@62049
  2098
  proof (cases "z \<in> \<int>")
eberlm@62049
  2099
    let ?h' = "\<lambda>z. Beta z (1 - z) * ((Digamma z - Digamma (1 - z)) * sin (z * of_real pi) +
eberlm@62049
  2100
                     of_real pi * cos (z * of_real pi))"
eberlm@62049
  2101
    case False
eberlm@62049
  2102
    from False have "eventually (\<lambda>t. t \<in> UNIV - \<int>) (nhds z)"
eberlm@62049
  2103
      by (intro eventually_nhds_in_open) (auto simp: open_Diff)
eberlm@62049
  2104
    hence "eventually (\<lambda>t. g t = ?g t) (nhds z)" by eventually_elim (simp add: g_def)
eberlm@62049
  2105
    moreover {
paulson@62131
  2106
      from False Ints_diff[of 1 "1-z"] have "1 - z \<notin> \<int>" by auto
eberlm@62049
  2107
      hence "(?g has_field_derivative ?h' z) (at z)" using nonpos_Ints_subset_Ints
eberlm@62049
  2108
        by (auto intro!: derivative_eq_intros simp: algebra_simps Beta_def)
eberlm@62049
  2109
      also from False have "sin (of_real pi * z) \<noteq> 0" by (subst sin_eq_0) auto
eberlm@62049
  2110
      hence "?h' z = h z * g z"
eberlm@62049
  2111
        using False unfolding g_def h_def cot_def by (simp add: field_simps Beta_def)
eberlm@62049
  2112
      finally have "(?g has_field_derivative (h z * g z)) (at z)" .
eberlm@62049
  2113
    }
eberlm@62049
  2114
    ultimately show ?thesis by (subst DERIV_cong_ev[OF refl _ refl])
eberlm@62049
  2115
  next
eberlm@62049
  2116
    case True
eberlm@62049
  2117
    then obtain n where z: "z = of_int n" by (auto elim!: Ints_cases)
eberlm@62049
  2118
    let ?t = "(\<lambda>z::complex. if z = 0 then 1 else sin z / z) \<circ> (\<lambda>z. of_real pi * z)"
eberlm@62049
  2119
    have deriv_0: "(g has_field_derivative 0) (at 0)"
eberlm@62049
  2120
    proof (subst DERIV_cong_ev[OF refl _ refl])
eberlm@62049
  2121
      show "eventually (\<lambda>z. g z = of_real pi * Gamma (1 + z) * Gamma (1 - z) * ?t z) (nhds 0)"
eberlm@62049
  2122
        using eventually_nhds_ball[OF zero_less_one, of "0::complex"]
eberlm@62049
  2123
      proof eventually_elim
eberlm@62049
  2124
        fix z :: complex assume z: "z \<in> ball 0 1"
eberlm@62049
  2125
        show "g z = of_real pi * Gamma (1 + z) * Gamma (1 - z) * ?t z"
eberlm@62049
  2126
        proof (cases "z = 0")
eberlm@62049
  2127
          assume z': "z \<noteq> 0"
eberlm@62049
  2128
          with z have z'': "z \<notin> \<int>\<^sub>\<le>\<^sub>0" "z \<notin> \<int>" by (auto elim!: Ints_cases simp: dist_0_norm)
eberlm@62049
  2129
          from Gamma_plus1[OF this(1)] have "Gamma z = Gamma (z + 1) / z" by simp
eberlm@62049
  2130
          with z'' z' show ?thesis by (simp add: g_def ac_simps)
eberlm@62049
  2131
        qed (simp add: g_def)
eberlm@62049
  2132
      qed
eberlm@62049
  2133
      have "(?t has_field_derivative (0 * of_real pi)) (at 0)"
paulson@62131
  2134
        using has_field_derivative_sin_z_over_z[of "UNIV :: complex set"]
eberlm@62049
  2135
        by (intro DERIV_chain) simp_all
eberlm@62049
  2136
      thus "((\<lambda>z. of_real pi * Gamma (1 + z) * Gamma (1 - z) * ?t z) has_field_derivative 0) (at 0)"
eberlm@62049
  2137
        by (auto intro!: derivative_eq_intros simp: o_def)
eberlm@62049
  2138
    qed
paulson@62131
  2139
eberlm@62049
  2140
    have "((g \<circ> (\<lambda>x. x - of_int n)) has_field_derivative 0 * 1) (at (of_int n))"
eberlm@62049
  2141
      using deriv_0 by (intro DERIV_chain) (auto intro!: derivative_eq_intros)
eberlm@62049
  2142
    also have "g \<circ> (\<lambda>x. x - of_int n) = g" by (intro ext) (simp add: g.minus_of_int)
eberlm@62049
  2143
    finally show "(g has_field_derivative (h z * g z)) (at z)" by (simp add: z h_def)
eberlm@62049
  2144
  qed
lp15@66512
  2145
eberlm@64969
  2146
  have g_holo [holomorphic_intros]: "g holomorphic_on A" for A
lp15@66512
  2147
    by (rule holomorphic_on_subset[of _ UNIV])
eberlm@64969
  2148
       (force simp: holomorphic_on_open intro!: derivative_intros)+
eberlm@62049
  2149
eberlm@62049
  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
eberlm@62049
  2151
  proof (cases "z \<in> \<int>")
eberlm@62049
  2152
    case True
eberlm@62049
  2153
    with that have "z = 0 \<or> z = 1" by (force elim!: Ints_cases)
eberlm@62049
  2154
    moreover have "g 0 * g (1/2) = Gamma (1/2)^2 * g 0"
eberlm@62049
  2155
      using fraction_not_in_ints[where 'a = complex, of 2 1] by (simp add: g_def power2_eq_square)
eberlm@62049
  2156
    moreover have "g (1/2) * g 1 = Gamma (1/2)^2 * g 1"
eberlm@62049
  2157
        using fraction_not_in_ints[where 'a = complex, of 2 1]
eberlm@62049
  2158
        by (simp add: g_def power2_eq_square Beta_def algebra_simps)
eberlm@62049
  2159
    ultimately show ?thesis by force
eberlm@62049
  2160
  next
eberlm@62049
  2161
    case False
eberlm@62049
  2162
    hence z: "z/2 \<notin> \<int>" "(z+1)/2 \<notin> \<int>" using Ints_diff[of "z+1" 1] by (auto elim!: Ints_cases)
eberlm@62049
  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)
paulson@62131
  2164
    from z have "1-z/2 \<notin> \<int>" "1-((z+1)/2) \<notin> \<int>"
eberlm@62049
  2165
      using Ints_diff[of 1 "1-z/2"] Ints_diff[of 1 "1-((z+1)/2)"] by auto
eberlm@62049
  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)
paulson@62131
  2167
    from z have "g (z/2) * g ((z+1)/2) =
eberlm@62049
  2168
      (Gamma (z/2) * Gamma ((z+1)/2)) * (Gamma (1-z/2) * Gamma (1-((z+1)/2))) *
eberlm@62049
  2169
      (sin (of_real pi * z/2) * sin (of_real pi * (z+1)/2))"
eberlm@62049
  2170
      by (simp add: g_def)
paulson@62131
  2171
    also from z' Gamma_legendre_duplication_aux[of "z/2"]
eberlm@62049
  2172
      have "Gamma (z/2) * Gamma ((z+1)/2) = exp ((1-z) * of_real (ln 2)) * Gamma (1/2) * Gamma z"
eberlm@62049
  2173
      by (simp add: add_divide_distrib)
eberlm@62049
  2174
    also from z'' Gamma_legendre_duplication_aux[of "1-(z+1)/2"]
paulson@62131
  2175
      have "Gamma (1-z/2) * Gamma (1-(z+1)/2) =
eberlm@62049
  2176
              Gamma (1-z) * Gamma (1/2) * exp (z * of_real (ln 2))"
eberlm@62049
  2177
      by (simp add: add_divide_distrib ac_simps)
eberlm@62049
  2178
    finally have "g (z/2) * g ((z+1)/2) = Gamma (1/2)^2 * (Gamma z * Gamma (1-z) *
eberlm@62049
  2179
                    (2 * (sin (of_real pi*z/2) * sin (of_real pi*(z+1)/2))))"
eberlm@62049
  2180
      by (simp add: add_ac power2_eq_square exp_add ring_distribs exp_diff exp_of_real)
eberlm@62049
  2181
    also have "sin (of_real pi*(z+1)/2) = cos (of_real pi*z/2)"
eberlm@62049
  2182
      using cos_sin_eq[of "- of_real pi * z/2", symmetric]
eberlm@62049
  2183
      by (simp add: ring_distribs add_divide_distrib ac_simps)
eberlm@62049
  2184
    also have "2 * (sin (of_real pi*z/2) * cos (of_real pi*z/2)) = sin (of_real pi * z)"
eberlm@62049
  2185
      by (subst sin_times_cos) (simp add: field_simps)
eberlm@62049
  2186
    also have "Gamma z * Gamma (1 - z) * sin (complex_of_real pi * z) = g z"
eberlm@62049
  2187
      using \<open>z \<notin> \<int>\<close> by (simp add: g_def)
eberlm@62049
  2188
    finally show ?thesis .
eberlm@62049
  2189
  qed
eberlm@62049
  2190
  have g_eq: "g (z/2) * g ((z+1)/2) = Gamma (1/2)^2 * g z" for z
eberlm@62049
  2191
  proof -
wenzelm@63040
  2192
    define r where "r = \<lfloor>Re z / 2\<rfloor>"
eberlm@62049
  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)
eberlm@62049
  2194
    also have "of_int (2*r) = 2 * of_int r" by simp