src/HOL/Analysis/Harmonic_Numbers.thy
changeset 63627 6ddb43c6b711
parent 63594 bd218a9320b5
child 63721 492bb53c3420
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/HOL/Analysis/Harmonic_Numbers.thy	Mon Aug 08 14:13:14 2016 +0200
     1.3 @@ -0,0 +1,533 @@
     1.4 +(*  Title:    HOL/Analysis/Harmonic_Numbers.thy
     1.5 +    Author:   Manuel Eberl, TU M√ľnchen
     1.6 +*)
     1.7 +
     1.8 +section \<open>Harmonic Numbers\<close>
     1.9 +
    1.10 +theory Harmonic_Numbers
    1.11 +imports
    1.12 +  Complex_Transcendental
    1.13 +  Summation_Tests
    1.14 +  Integral_Test
    1.15 +begin
    1.16 +
    1.17 +text \<open>
    1.18 +  The definition of the Harmonic Numbers and the Euler-Mascheroni constant.
    1.19 +  Also provides a reasonably accurate approximation of @{term "ln 2 :: real"}
    1.20 +  and the Euler-Mascheroni constant.
    1.21 +\<close>
    1.22 +
    1.23 +lemma ln_2_less_1: "ln 2 < (1::real)"
    1.24 +proof -
    1.25 +  have "2 < 5/(2::real)" by simp
    1.26 +  also have "5/2 \<le> exp (1::real)" using exp_lower_taylor_quadratic[of 1, simplified] by simp
    1.27 +  finally have "exp (ln 2) < exp (1::real)" by simp
    1.28 +  thus "ln 2 < (1::real)" by (subst (asm) exp_less_cancel_iff) simp
    1.29 +qed
    1.30 +
    1.31 +lemma setsum_Suc_diff':
    1.32 +  fixes f :: "nat \<Rightarrow> 'a::ab_group_add"
    1.33 +  assumes "m \<le> n"
    1.34 +  shows "(\<Sum>i = m..<n. f (Suc i) - f i) = f n - f m"
    1.35 +using assms by (induct n) (auto simp: le_Suc_eq)
    1.36 +
    1.37 +
    1.38 +subsection \<open>The Harmonic numbers\<close>
    1.39 +
    1.40 +definition harm :: "nat \<Rightarrow> 'a :: real_normed_field" where
    1.41 +  "harm n = (\<Sum>k=1..n. inverse (of_nat k))"
    1.42 +
    1.43 +lemma harm_altdef: "harm n = (\<Sum>k<n. inverse (of_nat (Suc k)))"
    1.44 +  unfolding harm_def by (induction n) simp_all
    1.45 +
    1.46 +lemma harm_Suc: "harm (Suc n) = harm n + inverse (of_nat (Suc n))"
    1.47 +  by (simp add: harm_def)
    1.48 +
    1.49 +lemma harm_nonneg: "harm n \<ge> (0 :: 'a :: {real_normed_field,linordered_field})"
    1.50 +  unfolding harm_def by (intro setsum_nonneg) simp_all
    1.51 +
    1.52 +lemma harm_pos: "n > 0 \<Longrightarrow> harm n > (0 :: 'a :: {real_normed_field,linordered_field})"
    1.53 +  unfolding harm_def by (intro setsum_pos) simp_all
    1.54 +
    1.55 +lemma of_real_harm: "of_real (harm n) = harm n"
    1.56 +  unfolding harm_def by simp
    1.57 +
    1.58 +lemma norm_harm: "norm (harm n) = harm n"
    1.59 +  by (subst of_real_harm [symmetric]) (simp add: harm_nonneg)
    1.60 +
    1.61 +lemma harm_expand:
    1.62 +  "harm 0 = 0"
    1.63 +  "harm (Suc 0) = 1"
    1.64 +  "harm (numeral n) = harm (pred_numeral n) + inverse (numeral n)"
    1.65 +proof -
    1.66 +  have "numeral n = Suc (pred_numeral n)" by simp
    1.67 +  also have "harm \<dots> = harm (pred_numeral n) + inverse (numeral n)"
    1.68 +    by (subst harm_Suc, subst numeral_eq_Suc[symmetric]) simp
    1.69 +  finally show "harm (numeral n) = harm (pred_numeral n) + inverse (numeral n)" .
    1.70 +qed (simp_all add: harm_def)
    1.71 +
    1.72 +lemma not_convergent_harm: "\<not>convergent (harm :: nat \<Rightarrow> 'a :: real_normed_field)"
    1.73 +proof -
    1.74 +  have "convergent (\<lambda>n. norm (harm n :: 'a)) \<longleftrightarrow>
    1.75 +            convergent (harm :: nat \<Rightarrow> real)" by (simp add: norm_harm)
    1.76 +  also have "\<dots> \<longleftrightarrow> convergent (\<lambda>n. \<Sum>k=Suc 0..Suc n. inverse (of_nat k) :: real)"
    1.77 +    unfolding harm_def[abs_def] by (subst convergent_Suc_iff) simp_all
    1.78 +  also have "... \<longleftrightarrow> convergent (\<lambda>n. \<Sum>k\<le>n. inverse (of_nat (Suc k)) :: real)"
    1.79 +    by (subst setsum_shift_bounds_cl_Suc_ivl) (simp add: atLeast0AtMost)
    1.80 +  also have "... \<longleftrightarrow> summable (\<lambda>n. inverse (of_nat n) :: real)"
    1.81 +    by (subst summable_Suc_iff [symmetric]) (simp add: summable_iff_convergent')
    1.82 +  also have "\<not>..." by (rule not_summable_harmonic)
    1.83 +  finally show ?thesis by (blast dest: convergent_norm)
    1.84 +qed
    1.85 +
    1.86 +lemma harm_pos_iff [simp]: "harm n > (0 :: 'a :: {real_normed_field,linordered_field}) \<longleftrightarrow> n > 0"
    1.87 +  by (rule iffI, cases n, simp add: harm_expand, simp, rule harm_pos)
    1.88 +
    1.89 +lemma ln_diff_le_inverse:
    1.90 +  assumes "x \<ge> (1::real)"
    1.91 +  shows   "ln (x + 1) - ln x < 1 / x"
    1.92 +proof -
    1.93 +  from assms have "\<exists>z>x. z < x + 1 \<and> ln (x + 1) - ln x = (x + 1 - x) * inverse z"
    1.94 +    by (intro MVT2) (auto intro!: derivative_eq_intros simp: field_simps)
    1.95 +  then obtain z where z: "z > x" "z < x + 1" "ln (x + 1) - ln x = inverse z" by auto
    1.96 +  have "ln (x + 1) - ln x = inverse z" by fact
    1.97 +  also from z(1,2) assms have "\<dots> < 1 / x" by (simp add: field_simps)
    1.98 +  finally show ?thesis .
    1.99 +qed
   1.100 +
   1.101 +lemma ln_le_harm: "ln (real n + 1) \<le> (harm n :: real)"
   1.102 +proof (induction n)
   1.103 +  fix n assume IH: "ln (real n + 1) \<le> harm n"
   1.104 +  have "ln (real (Suc n) + 1) = ln (real n + 1) + (ln (real n + 2) - ln (real n + 1))" by simp
   1.105 +  also have "(ln (real n + 2) - ln (real n + 1)) \<le> 1 / real (Suc n)"
   1.106 +    using ln_diff_le_inverse[of "real n + 1"] by (simp add: add_ac)
   1.107 +  also note IH
   1.108 +  also have "harm n + 1 / real (Suc n) = harm (Suc n)" by (simp add: harm_Suc field_simps)
   1.109 +  finally show "ln (real (Suc n) + 1) \<le> harm (Suc n)" by - simp
   1.110 +qed (simp_all add: harm_def)
   1.111 +
   1.112 +
   1.113 +subsection \<open>The Euler--Mascheroni constant\<close>
   1.114 +
   1.115 +text \<open>
   1.116 +  The limit of the difference between the partial harmonic sum and the natural logarithm
   1.117 +  (approximately 0.577216). This value occurs e.g. in the definition of the Gamma function.
   1.118 + \<close>
   1.119 +definition euler_mascheroni :: "'a :: real_normed_algebra_1" where
   1.120 +  "euler_mascheroni = of_real (lim (\<lambda>n. harm n - ln (of_nat n)))"
   1.121 +
   1.122 +lemma of_real_euler_mascheroni [simp]: "of_real euler_mascheroni = euler_mascheroni"
   1.123 +  by (simp add: euler_mascheroni_def)
   1.124 +
   1.125 +interpretation euler_mascheroni: antimono_fun_sum_integral_diff "\<lambda>x. inverse (x + 1)"
   1.126 +  by unfold_locales (auto intro!: continuous_intros)
   1.127 +
   1.128 +lemma euler_mascheroni_sum_integral_diff_series:
   1.129 +  "euler_mascheroni.sum_integral_diff_series n = harm (Suc n) - ln (of_nat (Suc n))"
   1.130 +proof -
   1.131 +  have "harm (Suc n) = (\<Sum>k=0..n. inverse (of_nat k + 1) :: real)" unfolding harm_def
   1.132 +    unfolding One_nat_def by (subst setsum_shift_bounds_cl_Suc_ivl) (simp add: add_ac)
   1.133 +  moreover have "((\<lambda>x. inverse (x + 1) :: real) has_integral ln (of_nat n + 1) - ln (0 + 1))
   1.134 +                   {0..of_nat n}"
   1.135 +    by (intro fundamental_theorem_of_calculus)
   1.136 +       (auto intro!: derivative_eq_intros simp: divide_inverse
   1.137 +           has_field_derivative_iff_has_vector_derivative[symmetric])
   1.138 +  hence "integral {0..of_nat n} (\<lambda>x. inverse (x + 1) :: real) = ln (of_nat (Suc n))"
   1.139 +    by (auto dest!: integral_unique)
   1.140 +  ultimately show ?thesis
   1.141 +    by (simp add: euler_mascheroni.sum_integral_diff_series_def atLeast0AtMost)
   1.142 +qed
   1.143 +
   1.144 +lemma euler_mascheroni_sequence_decreasing:
   1.145 +  "m > 0 \<Longrightarrow> m \<le> n \<Longrightarrow> harm n - ln (of_nat n) \<le> harm m - ln (of_nat m :: real)"
   1.146 +  by (cases m, simp, cases n, simp, hypsubst,
   1.147 +      subst (1 2) euler_mascheroni_sum_integral_diff_series [symmetric],
   1.148 +      rule euler_mascheroni.sum_integral_diff_series_antimono, simp)
   1.149 +
   1.150 +lemma euler_mascheroni_sequence_nonneg:
   1.151 +  "n > 0 \<Longrightarrow> harm n - ln (of_nat n) \<ge> (0::real)"
   1.152 +  by (cases n, simp, hypsubst, subst euler_mascheroni_sum_integral_diff_series [symmetric],
   1.153 +      rule euler_mascheroni.sum_integral_diff_series_nonneg)
   1.154 +
   1.155 +lemma euler_mascheroni_convergent: "convergent (\<lambda>n. harm n - ln (of_nat n) :: real)"
   1.156 +proof -
   1.157 +  have A: "(\<lambda>n. harm (Suc n) - ln (of_nat (Suc n))) =
   1.158 +             euler_mascheroni.sum_integral_diff_series"
   1.159 +    by (subst euler_mascheroni_sum_integral_diff_series [symmetric]) (rule refl)
   1.160 +  have "convergent (\<lambda>n. harm (Suc n) - ln (of_nat (Suc n) :: real))"
   1.161 +    by (subst A) (fact euler_mascheroni.sum_integral_diff_series_convergent)
   1.162 +  thus ?thesis by (subst (asm) convergent_Suc_iff)
   1.163 +qed
   1.164 +
   1.165 +lemma euler_mascheroni_LIMSEQ:
   1.166 +  "(\<lambda>n. harm n - ln (of_nat n) :: real) \<longlonglongrightarrow> euler_mascheroni"
   1.167 +  unfolding euler_mascheroni_def
   1.168 +  by (simp add: convergent_LIMSEQ_iff [symmetric] euler_mascheroni_convergent)
   1.169 +
   1.170 +lemma euler_mascheroni_LIMSEQ_of_real:
   1.171 +  "(\<lambda>n. of_real (harm n - ln (of_nat n))) \<longlonglongrightarrow>
   1.172 +      (euler_mascheroni :: 'a :: {real_normed_algebra_1, topological_space})"
   1.173 +proof -
   1.174 +  have "(\<lambda>n. of_real (harm n - ln (of_nat n))) \<longlonglongrightarrow> (of_real (euler_mascheroni) :: 'a)"
   1.175 +    by (intro tendsto_of_real euler_mascheroni_LIMSEQ)
   1.176 +  thus ?thesis by simp
   1.177 +qed
   1.178 +
   1.179 +lemma euler_mascheroni_sum:
   1.180 +  "(\<lambda>n. inverse (of_nat (n+1)) + ln (of_nat (n+1)) - ln (of_nat (n+2)) :: real)
   1.181 +       sums euler_mascheroni"
   1.182 + using sums_add[OF telescope_sums[OF LIMSEQ_Suc[OF euler_mascheroni_LIMSEQ]]
   1.183 +                   telescope_sums'[OF LIMSEQ_inverse_real_of_nat]]
   1.184 +  by (simp_all add: harm_def algebra_simps)
   1.185 +
   1.186 +lemma alternating_harmonic_series_sums: "(\<lambda>k. (-1)^k / real_of_nat (Suc k)) sums ln 2"
   1.187 +proof -
   1.188 +  let ?f = "\<lambda>n. harm n - ln (real_of_nat n)"
   1.189 +  let ?g = "\<lambda>n. if even n then 0 else (2::real)"
   1.190 +  let ?em = "\<lambda>n. harm n - ln (real_of_nat n)"
   1.191 +  have "eventually (\<lambda>n. ?em (2*n) - ?em n + ln 2 = (\<Sum>k<2*n. (-1)^k / real_of_nat (Suc k))) at_top"
   1.192 +    using eventually_gt_at_top[of "0::nat"]
   1.193 +  proof eventually_elim
   1.194 +    fix n :: nat assume n: "n > 0"
   1.195 +    have "(\<Sum>k<2*n. (-1)^k / real_of_nat (Suc k)) =
   1.196 +              (\<Sum>k<2*n. ((-1)^k + ?g k) / of_nat (Suc k)) - (\<Sum>k<2*n. ?g k / of_nat (Suc k))"
   1.197 +      by (simp add: setsum.distrib algebra_simps divide_inverse)
   1.198 +    also have "(\<Sum>k<2*n. ((-1)^k + ?g k) / real_of_nat (Suc k)) = harm (2*n)"
   1.199 +      unfolding harm_altdef by (intro setsum.cong) (auto simp: field_simps)
   1.200 +    also have "(\<Sum>k<2*n. ?g k / real_of_nat (Suc k)) = (\<Sum>k|k<2*n \<and> odd k. ?g k / of_nat (Suc k))"
   1.201 +      by (intro setsum.mono_neutral_right) auto
   1.202 +    also have "\<dots> = (\<Sum>k|k<2*n \<and> odd k. 2 / (real_of_nat (Suc k)))"
   1.203 +      by (intro setsum.cong) auto
   1.204 +    also have "(\<Sum>k|k<2*n \<and> odd k. 2 / (real_of_nat (Suc k))) = harm n"
   1.205 +      unfolding harm_altdef
   1.206 +      by (intro setsum.reindex_cong[of "\<lambda>n. 2*n+1"]) (auto simp: inj_on_def field_simps elim!: oddE)
   1.207 +    also have "harm (2*n) - harm n = ?em (2*n) - ?em n + ln 2" using n
   1.208 +      by (simp_all add: algebra_simps ln_mult)
   1.209 +    finally show "?em (2*n) - ?em n + ln 2 = (\<Sum>k<2*n. (-1)^k / real_of_nat (Suc k))" ..
   1.210 +  qed
   1.211 +  moreover have "(\<lambda>n. ?em (2*n) - ?em n + ln (2::real))
   1.212 +                     \<longlonglongrightarrow> euler_mascheroni - euler_mascheroni + ln 2"
   1.213 +    by (intro tendsto_intros euler_mascheroni_LIMSEQ filterlim_compose[OF euler_mascheroni_LIMSEQ]
   1.214 +              filterlim_subseq) (auto simp: subseq_def)
   1.215 +  hence "(\<lambda>n. ?em (2*n) - ?em n + ln (2::real)) \<longlonglongrightarrow> ln 2" by simp
   1.216 +  ultimately have "(\<lambda>n. (\<Sum>k<2*n. (-1)^k / real_of_nat (Suc k))) \<longlonglongrightarrow> ln 2"
   1.217 +    by (rule Lim_transform_eventually)
   1.218 +
   1.219 +  moreover have "summable (\<lambda>k. (-1)^k * inverse (real_of_nat (Suc k)))"
   1.220 +    using LIMSEQ_inverse_real_of_nat
   1.221 +    by (intro summable_Leibniz(1) decseq_imp_monoseq decseq_SucI) simp_all
   1.222 +  hence A: "(\<lambda>n. \<Sum>k<n. (-1)^k / real_of_nat (Suc k)) \<longlonglongrightarrow> (\<Sum>k. (-1)^k / real_of_nat (Suc k))"
   1.223 +    by (simp add: summable_sums_iff divide_inverse sums_def)
   1.224 +  from filterlim_compose[OF this filterlim_subseq[of "op * (2::nat)"]]
   1.225 +    have "(\<lambda>n. \<Sum>k<2*n. (-1)^k / real_of_nat (Suc k)) \<longlonglongrightarrow> (\<Sum>k. (-1)^k / real_of_nat (Suc k))"
   1.226 +    by (simp add: subseq_def)
   1.227 +  ultimately have "(\<Sum>k. (- 1) ^ k / real_of_nat (Suc k)) = ln 2" by (intro LIMSEQ_unique)
   1.228 +  with A show ?thesis by (simp add: sums_def)
   1.229 +qed
   1.230 +
   1.231 +lemma alternating_harmonic_series_sums':
   1.232 +  "(\<lambda>k. inverse (real_of_nat (2*k+1)) - inverse (real_of_nat (2*k+2))) sums ln 2"
   1.233 +unfolding sums_def
   1.234 +proof (rule Lim_transform_eventually)
   1.235 +  show "(\<lambda>n. \<Sum>k<2*n. (-1)^k / (real_of_nat (Suc k))) \<longlonglongrightarrow> ln 2"
   1.236 +    using alternating_harmonic_series_sums unfolding sums_def
   1.237 +    by (rule filterlim_compose) (rule mult_nat_left_at_top, simp)
   1.238 +  show "eventually (\<lambda>n. (\<Sum>k<2*n. (-1)^k / (real_of_nat (Suc k))) =
   1.239 +            (\<Sum>k<n. inverse (real_of_nat (2*k+1)) - inverse (real_of_nat (2*k+2)))) sequentially"
   1.240 +  proof (intro always_eventually allI)
   1.241 +    fix n :: nat
   1.242 +    show "(\<Sum>k<2*n. (-1)^k / (real_of_nat (Suc k))) =
   1.243 +              (\<Sum>k<n. inverse (real_of_nat (2*k+1)) - inverse (real_of_nat (2*k+2)))"
   1.244 +      by (induction n) (simp_all add: inverse_eq_divide)
   1.245 +  qed
   1.246 +qed
   1.247 +
   1.248 +
   1.249 +subsection \<open>Bounds on the Euler--Mascheroni constant\<close>
   1.250 +
   1.251 +(* TODO: Move? *)
   1.252 +lemma ln_inverse_approx_le:
   1.253 +  assumes "(x::real) > 0" "a > 0"
   1.254 +  shows   "ln (x + a) - ln x \<le> a * (inverse x + inverse (x + a))/2" (is "_ \<le> ?A")
   1.255 +proof -
   1.256 +  define f' where "f' = (inverse (x + a) - inverse x)/a"
   1.257 +  have f'_nonpos: "f' \<le> 0" using assms by (simp add: f'_def divide_simps)
   1.258 +  let ?f = "\<lambda>t. (t - x) * f' + inverse x"
   1.259 +  let ?F = "\<lambda>t. (t - x)^2 * f' / 2 + t * inverse x"
   1.260 +  have diff: "\<forall>t\<in>{x..x+a}. (?F has_vector_derivative ?f t)
   1.261 +                               (at t within {x..x+a})" using assms
   1.262 +    by (auto intro!: derivative_eq_intros
   1.263 +             simp: has_field_derivative_iff_has_vector_derivative[symmetric])
   1.264 +  from assms have "(?f has_integral (?F (x+a) - ?F x)) {x..x+a}"
   1.265 +    by (intro fundamental_theorem_of_calculus[OF _ diff])
   1.266 +       (auto simp: has_field_derivative_iff_has_vector_derivative[symmetric] field_simps
   1.267 +             intro!: derivative_eq_intros)
   1.268 +  also have "?F (x+a) - ?F x = (a*2 + f'*a\<^sup>2*x) / (2*x)" using assms by (simp add: field_simps)
   1.269 +  also have "f'*a^2 = - (a^2) / (x*(x + a))" using assms
   1.270 +    by (simp add: divide_simps f'_def power2_eq_square)
   1.271 +  also have "(a*2 + - a\<^sup>2/(x*(x+a))*x) / (2*x) = ?A" using assms
   1.272 +    by (simp add: divide_simps power2_eq_square) (simp add: algebra_simps)
   1.273 +  finally have int1: "((\<lambda>t. (t - x) * f' + inverse x) has_integral ?A) {x..x + a}" .
   1.274 +
   1.275 +  from assms have int2: "(inverse has_integral (ln (x + a) - ln x)) {x..x+a}"
   1.276 +    by (intro fundamental_theorem_of_calculus)
   1.277 +       (auto simp: has_field_derivative_iff_has_vector_derivative[symmetric] divide_simps
   1.278 +             intro!: derivative_eq_intros)
   1.279 +  hence "ln (x + a) - ln x = integral {x..x+a} inverse" by (simp add: integral_unique)
   1.280 +  also have ineq: "\<forall>xa\<in>{x..x + a}. inverse xa \<le> (xa - x) * f' + inverse x"
   1.281 +  proof
   1.282 +    fix t assume t': "t \<in> {x..x+a}"
   1.283 +    with assms have t: "0 \<le> (t - x) / a" "(t - x) / a \<le> 1" by simp_all
   1.284 +    have "inverse t = inverse ((1 - (t - x) / a) *\<^sub>R x + ((t - x) / a) *\<^sub>R (x + a))" (is "_ = ?A")
   1.285 +      using assms t' by (simp add: field_simps)
   1.286 +    also from assms have "convex_on {x..x+a} inverse" by (intro convex_on_inverse) auto
   1.287 +    from convex_onD_Icc[OF this _ t] assms
   1.288 +      have "?A \<le> (1 - (t - x) / a) * inverse x + (t - x) / a * inverse (x + a)" by simp
   1.289 +    also have "\<dots> = (t - x) * f' + inverse x" using assms
   1.290 +      by (simp add: f'_def divide_simps) (simp add: f'_def field_simps)
   1.291 +    finally show "inverse t \<le> (t - x) * f' + inverse x" .
   1.292 +  qed
   1.293 +  hence "integral {x..x+a} inverse \<le> integral {x..x+a} ?f" using f'_nonpos assms
   1.294 +    by (intro integral_le has_integral_integrable[OF int1] has_integral_integrable[OF int2] ineq)
   1.295 +  also have "\<dots> = ?A" using int1 by (rule integral_unique)
   1.296 +  finally show ?thesis .
   1.297 +qed
   1.298 +
   1.299 +lemma ln_inverse_approx_ge:
   1.300 +  assumes "(x::real) > 0" "x < y"
   1.301 +  shows   "ln y - ln x \<ge> 2 * (y - x) / (x + y)" (is "_ \<ge> ?A")
   1.302 +proof -
   1.303 +  define m where "m = (x+y)/2"
   1.304 +  define f' where "f' = -inverse (m^2)"
   1.305 +  from assms have m: "m > 0" by (simp add: m_def)
   1.306 +  let ?F = "\<lambda>t. (t - m)^2 * f' / 2 + t / m"
   1.307 +  from assms have "((\<lambda>t. (t - m) * f' + inverse m) has_integral (?F y - ?F x)) {x..y}"
   1.308 +    by (intro fundamental_theorem_of_calculus)
   1.309 +       (auto simp: has_field_derivative_iff_has_vector_derivative[symmetric] divide_simps
   1.310 +             intro!: derivative_eq_intros)
   1.311 +  also from m have "?F y - ?F x = ((y - m)^2 - (x - m)^2) * f' / 2 + (y - x) / m"
   1.312 +    by (simp add: field_simps)
   1.313 +  also have "((y - m)^2 - (x - m)^2) = 0" by (simp add: m_def power2_eq_square field_simps)
   1.314 +  also have "0 * f' / 2 + (y - x) / m = ?A" by (simp add: m_def)
   1.315 +  finally have int1: "((\<lambda>t. (t - m) * f' + inverse m) has_integral ?A) {x..y}" .
   1.316 +
   1.317 +  from assms have int2: "(inverse has_integral (ln y - ln x)) {x..y}"
   1.318 +    by (intro fundamental_theorem_of_calculus)
   1.319 +       (auto simp: has_field_derivative_iff_has_vector_derivative[symmetric] divide_simps
   1.320 +             intro!: derivative_eq_intros)
   1.321 +  hence "ln y - ln x = integral {x..y} inverse" by (simp add: integral_unique)
   1.322 +  also have ineq: "\<forall>xa\<in>{x..y}. inverse xa \<ge> (xa - m) * f' + inverse m"
   1.323 +  proof
   1.324 +    fix t assume t: "t \<in> {x..y}"
   1.325 +    from t assms have "inverse t - inverse m \<ge> f' * (t - m)"
   1.326 +      by (intro convex_on_imp_above_tangent[of "{0<..}"] convex_on_inverse)
   1.327 +         (auto simp: m_def interior_open f'_def power2_eq_square intro!: derivative_eq_intros)
   1.328 +    thus "(t - m) * f' + inverse m \<le> inverse t" by (simp add: algebra_simps)
   1.329 +  qed
   1.330 +  hence "integral {x..y} inverse \<ge> integral {x..y} (\<lambda>t. (t - m) * f' + inverse m)"
   1.331 +    using int1 int2 by (intro integral_le has_integral_integrable)
   1.332 +  also have "integral {x..y} (\<lambda>t. (t - m) * f' + inverse m) = ?A"
   1.333 +    using integral_unique[OF int1] by simp
   1.334 +  finally show ?thesis .
   1.335 +qed
   1.336 +
   1.337 +
   1.338 +lemma euler_mascheroni_lower:
   1.339 +        "euler_mascheroni \<ge> harm (Suc n) - ln (real_of_nat (n + 2)) + 1/real_of_nat (2 * (n + 2))"
   1.340 +  and euler_mascheroni_upper:
   1.341 +        "euler_mascheroni \<le> harm (Suc n) - ln (real_of_nat (n + 2)) + 1/real_of_nat (2 * (n + 1))"
   1.342 +proof -
   1.343 +  define D :: "_ \<Rightarrow> real"
   1.344 +    where "D n = inverse (of_nat (n+1)) + ln (of_nat (n+1)) - ln (of_nat (n+2))" for n
   1.345 +  let ?g = "\<lambda>n. ln (of_nat (n+2)) - ln (of_nat (n+1)) - inverse (of_nat (n+1)) :: real"
   1.346 +  define inv where [abs_def]: "inv n = inverse (real_of_nat n)" for n
   1.347 +  fix n :: nat
   1.348 +  note summable = sums_summable[OF euler_mascheroni_sum, folded D_def]
   1.349 +  have sums: "(\<lambda>k. (inv (Suc (k + (n+1))) - inv (Suc (Suc k + (n+1))))/2) sums ((inv (Suc (0 + (n+1))) - 0)/2)"
   1.350 +    unfolding inv_def
   1.351 +    by (intro sums_divide telescope_sums' LIMSEQ_ignore_initial_segment LIMSEQ_inverse_real_of_nat)
   1.352 +  have sums': "(\<lambda>k. (inv (Suc (k + n)) - inv (Suc (Suc k + n)))/2) sums ((inv (Suc (0 + n)) - 0)/2)"
   1.353 +    unfolding inv_def
   1.354 +    by (intro sums_divide telescope_sums' LIMSEQ_ignore_initial_segment LIMSEQ_inverse_real_of_nat)
   1.355 +  from euler_mascheroni_sum have "euler_mascheroni = (\<Sum>k. D k)"
   1.356 +    by (simp add: sums_iff D_def)
   1.357 +  also have "\<dots> = (\<Sum>k. D (k + Suc n)) + (\<Sum>k\<le>n. D k)"
   1.358 +    by (subst suminf_split_initial_segment[OF summable, of "Suc n"], subst lessThan_Suc_atMost) simp
   1.359 +  finally have sum: "(\<Sum>k\<le>n. D k) - euler_mascheroni = -(\<Sum>k. D (k + Suc n))" by simp
   1.360 +
   1.361 +  note sum
   1.362 +  also have "\<dots> \<le> -(\<Sum>k. (inv (k + Suc n + 1) - inv (k + Suc n + 2)) / 2)"
   1.363 +  proof (intro le_imp_neg_le suminf_le allI summable_ignore_initial_segment[OF summable])
   1.364 +    fix k' :: nat
   1.365 +    define k where "k = k' + Suc n"
   1.366 +    hence k: "k > 0" by (simp add: k_def)
   1.367 +    have "real_of_nat (k+1) > 0" by (simp add: k_def)
   1.368 +    with ln_inverse_approx_le[OF this zero_less_one]
   1.369 +      have "ln (of_nat k + 2) - ln (of_nat k + 1) \<le> (inv (k+1) + inv (k+2))/2"
   1.370 +      by (simp add: inv_def add_ac)
   1.371 +    hence "(inv (k+1) - inv (k+2))/2 \<le> inv (k+1) + ln (of_nat (k+1)) - ln (of_nat (k+2))"
   1.372 +      by (simp add: field_simps)
   1.373 +    also have "\<dots> = D k" unfolding D_def inv_def ..
   1.374 +    finally show "D (k' + Suc n) \<ge> (inv (k' + Suc n + 1) - inv (k' + Suc n + 2)) / 2"
   1.375 +      by (simp add: k_def)
   1.376 +    from sums_summable[OF sums]
   1.377 +      show "summable (\<lambda>k. (inv (k + Suc n + 1) - inv (k + Suc n + 2))/2)" by simp
   1.378 +  qed
   1.379 +  also from sums have "\<dots> = -inv (n+2) / 2" by (simp add: sums_iff)
   1.380 +  finally have "euler_mascheroni \<ge> (\<Sum>k\<le>n. D k) + 1 / (of_nat (2 * (n+2)))"
   1.381 +    by (simp add: inv_def field_simps)
   1.382 +  also have "(\<Sum>k\<le>n. D k) = harm (Suc n) - (\<Sum>k\<le>n. ln (real_of_nat (Suc k+1)) - ln (of_nat (k+1)))"
   1.383 +    unfolding harm_altdef D_def by (subst lessThan_Suc_atMost) (simp add:  setsum.distrib setsum_subtractf)
   1.384 +  also have "(\<Sum>k\<le>n. ln (real_of_nat (Suc k+1)) - ln (of_nat (k+1))) = ln (of_nat (n+2))"
   1.385 +    by (subst atLeast0AtMost [symmetric], subst setsum_Suc_diff) simp_all
   1.386 +  finally show "euler_mascheroni \<ge> harm (Suc n) - ln (real_of_nat (n + 2)) + 1/real_of_nat (2 * (n + 2))"
   1.387 +    by simp
   1.388 +
   1.389 +  note sum
   1.390 +  also have "-(\<Sum>k. D (k + Suc n)) \<ge> -(\<Sum>k. (inv (Suc (k + n)) - inv (Suc (Suc k + n)))/2)"
   1.391 +  proof (intro le_imp_neg_le suminf_le allI summable_ignore_initial_segment[OF summable])
   1.392 +    fix k' :: nat
   1.393 +    define k where "k = k' + Suc n"
   1.394 +    hence k: "k > 0" by (simp add: k_def)
   1.395 +    have "real_of_nat (k+1) > 0" by (simp add: k_def)
   1.396 +    from ln_inverse_approx_ge[of "of_nat k + 1" "of_nat k + 2"]
   1.397 +      have "2 / (2 * real_of_nat k + 3) \<le> ln (of_nat (k+2)) - ln (real_of_nat (k+1))"
   1.398 +      by (simp add: add_ac)
   1.399 +    hence "D k \<le> 1 / real_of_nat (k+1) - 2 / (2 * real_of_nat k + 3)"
   1.400 +      by (simp add: D_def inverse_eq_divide inv_def)
   1.401 +    also have "\<dots> = inv ((k+1)*(2*k+3))" unfolding inv_def by (simp add: field_simps)
   1.402 +    also have "\<dots> \<le> inv (2*k*(k+1))" unfolding inv_def using k
   1.403 +      by (intro le_imp_inverse_le)
   1.404 +         (simp add: algebra_simps, simp del: of_nat_add)
   1.405 +    also have "\<dots> = (inv k - inv (k+1))/2" unfolding inv_def using k
   1.406 +      by (simp add: divide_simps del: of_nat_mult) (simp add: algebra_simps)
   1.407 +    finally show "D k \<le> (inv (Suc (k' + n)) - inv (Suc (Suc k' + n)))/2" unfolding k_def by simp
   1.408 +  next
   1.409 +    from sums_summable[OF sums']
   1.410 +      show "summable (\<lambda>k. (inv (Suc (k + n)) - inv (Suc (Suc k + n)))/2)" by simp
   1.411 +  qed
   1.412 +  also from sums' have "(\<Sum>k. (inv (Suc (k + n)) - inv (Suc (Suc k + n)))/2) = inv (n+1)/2"
   1.413 +    by (simp add: sums_iff)
   1.414 +  finally have "euler_mascheroni \<le> (\<Sum>k\<le>n. D k) + 1 / of_nat (2 * (n+1))"
   1.415 +    by (simp add: inv_def field_simps)
   1.416 +  also have "(\<Sum>k\<le>n. D k) = harm (Suc n) - (\<Sum>k\<le>n. ln (real_of_nat (Suc k+1)) - ln (of_nat (k+1)))"
   1.417 +    unfolding harm_altdef D_def by (subst lessThan_Suc_atMost) (simp add:  setsum.distrib setsum_subtractf)
   1.418 +  also have "(\<Sum>k\<le>n. ln (real_of_nat (Suc k+1)) - ln (of_nat (k+1))) = ln (of_nat (n+2))"
   1.419 +    by (subst atLeast0AtMost [symmetric], subst setsum_Suc_diff) simp_all
   1.420 +  finally show "euler_mascheroni \<le> harm (Suc n) - ln (real_of_nat (n + 2)) + 1/real_of_nat (2 * (n + 1))"
   1.421 +    by simp
   1.422 +qed
   1.423 +
   1.424 +lemma euler_mascheroni_pos: "euler_mascheroni > (0::real)"
   1.425 +  using euler_mascheroni_lower[of 0] ln_2_less_1 by (simp add: harm_def)
   1.426 +
   1.427 +context
   1.428 +begin
   1.429 +
   1.430 +private lemma ln_approx_aux:
   1.431 +  fixes n :: nat and x :: real
   1.432 +  defines "y \<equiv> (x-1)/(x+1)"
   1.433 +  assumes x: "x > 0" "x \<noteq> 1"
   1.434 +  shows "inverse (2*y^(2*n+1)) * (ln x - (\<Sum>k<n. 2*y^(2*k+1) / of_nat (2*k+1))) \<in>
   1.435 +            {0..(1 / (1 - y^2) / of_nat (2*n+1))}"
   1.436 +proof -
   1.437 +  from x have norm_y: "norm y < 1" unfolding y_def by simp
   1.438 +  from power_strict_mono[OF this, of 2] have norm_y': "norm y^2 < 1" by simp
   1.439 +
   1.440 +  let ?f = "\<lambda>k. 2 * y ^ (2*k+1) / of_nat (2*k+1)"
   1.441 +  note sums = ln_series_quadratic[OF x(1)]
   1.442 +  define c where "c = inverse (2*y^(2*n+1))"
   1.443 +  let ?d = "c * (ln x - (\<Sum>k<n. ?f k))"
   1.444 +  have "\<forall>k. y\<^sup>2^k / of_nat (2*(k+n)+1) \<le> y\<^sup>2 ^ k / of_nat (2*n+1)"
   1.445 +    by (intro allI divide_left_mono mult_right_mono mult_pos_pos zero_le_power[of "y^2"]) simp_all
   1.446 +  moreover {
   1.447 +    have "(\<lambda>k. ?f (k + n)) sums (ln x - (\<Sum>k<n. ?f k))"
   1.448 +      using sums_split_initial_segment[OF sums] by (simp add: y_def)
   1.449 +    hence "(\<lambda>k. c * ?f (k + n)) sums ?d" by (rule sums_mult)
   1.450 +    also have "(\<lambda>k. c * (2*y^(2*(k+n)+1) / of_nat (2*(k+n)+1))) =
   1.451 +                   (\<lambda>k. (c * (2*y^(2*n+1))) * ((y^2)^k / of_nat (2*(k+n)+1)))"
   1.452 +      by (simp only: ring_distribs power_add power_mult) (simp add: mult_ac)
   1.453 +    also from x have "c * (2*y^(2*n+1)) = 1" by (simp add: c_def y_def)
   1.454 +    finally have "(\<lambda>k. (y^2)^k / of_nat (2*(k+n)+1)) sums ?d" by simp
   1.455 +  } note sums' = this
   1.456 +  moreover from norm_y' have "(\<lambda>k. (y^2)^k / of_nat (2*n+1)) sums (1 / (1 - y^2) / of_nat (2*n+1))"
   1.457 +    by (intro sums_divide geometric_sums) (simp_all add: norm_power)
   1.458 +  ultimately have "?d \<le> (1 / (1 - y^2) / of_nat (2*n+1))" by (rule sums_le)
   1.459 +  moreover have "c * (ln x - (\<Sum>k<n. 2 * y ^ (2 * k + 1) / real_of_nat (2 * k + 1))) \<ge> 0"
   1.460 +    by (intro sums_le[OF _ sums_zero sums']) simp_all
   1.461 +  ultimately show ?thesis unfolding c_def by simp
   1.462 +qed
   1.463 +
   1.464 +lemma
   1.465 +  fixes n :: nat and x :: real
   1.466 +  defines "y \<equiv> (x-1)/(x+1)"
   1.467 +  defines "approx \<equiv> (\<Sum>k<n. 2*y^(2*k+1) / of_nat (2*k+1))"
   1.468 +  defines "d \<equiv> y^(2*n+1) / (1 - y^2) / of_nat (2*n+1)"
   1.469 +  assumes x: "x > 1"
   1.470 +  shows   ln_approx_bounds: "ln x \<in> {approx..approx + 2*d}"
   1.471 +  and     ln_approx_abs:    "abs (ln x - (approx + d)) \<le> d"
   1.472 +proof -
   1.473 +  define c where "c = 2*y^(2*n+1)"
   1.474 +  from x have c_pos: "c > 0" unfolding c_def y_def
   1.475 +    by (intro mult_pos_pos zero_less_power) simp_all
   1.476 +  have A: "inverse c * (ln x - (\<Sum>k<n. 2*y^(2*k+1) / of_nat (2*k+1))) \<in>
   1.477 +              {0.. (1 / (1 - y^2) / of_nat (2*n+1))}" using assms unfolding y_def c_def
   1.478 +    by (intro ln_approx_aux) simp_all
   1.479 +  hence "inverse c * (ln x - (\<Sum>k<n. 2*y^(2*k+1)/of_nat (2*k+1))) \<le> (1 / (1-y^2) / of_nat (2*n+1))"
   1.480 +    by simp
   1.481 +  hence "(ln x - (\<Sum>k<n. 2*y^(2*k+1) / of_nat (2*k+1))) / c \<le> (1 / (1 - y^2) / of_nat (2*n+1))"
   1.482 +    by (auto simp add: divide_simps)
   1.483 +  with c_pos have "ln x \<le> c / (1 - y^2) / of_nat (2*n+1) + approx"
   1.484 +    by (subst (asm) pos_divide_le_eq) (simp_all add: mult_ac approx_def)
   1.485 +  moreover {
   1.486 +    from A c_pos have "0 \<le> c * (inverse c * (ln x - (\<Sum>k<n. 2*y^(2*k+1) / of_nat (2*k+1))))"
   1.487 +      by (intro mult_nonneg_nonneg[of c]) simp_all
   1.488 +    also have "\<dots> = (c * inverse c) * (ln x - (\<Sum>k<n. 2*y^(2*k+1) / of_nat (2*k+1)))"
   1.489 +      by (simp add: mult_ac)
   1.490 +    also from c_pos have "c * inverse c = 1" by simp
   1.491 +    finally have "ln x \<ge> approx" by (simp add: approx_def)
   1.492 +  }
   1.493 +  ultimately show "ln x \<in> {approx..approx + 2*d}" by (simp add: c_def d_def)
   1.494 +  thus "abs (ln x - (approx + d)) \<le> d" by auto
   1.495 +qed
   1.496 +
   1.497 +end
   1.498 +
   1.499 +lemma euler_mascheroni_bounds:
   1.500 +  fixes n :: nat assumes "n \<ge> 1" defines "t \<equiv> harm n - ln (of_nat (Suc n)) :: real"
   1.501 +  shows "euler_mascheroni \<in> {t + inverse (of_nat (2*(n+1)))..t + inverse (of_nat (2*n))}"
   1.502 +  using assms euler_mascheroni_upper[of "n-1"] euler_mascheroni_lower[of "n-1"]
   1.503 +  unfolding t_def by (cases n) (simp_all add: harm_Suc t_def inverse_eq_divide)
   1.504 +
   1.505 +lemma euler_mascheroni_bounds':
   1.506 +  fixes n :: nat assumes "n \<ge> 1" "ln (real_of_nat (Suc n)) \<in> {l<..<u}"
   1.507 +  shows "euler_mascheroni \<in>
   1.508 +           {harm n - u + inverse (of_nat (2*(n+1)))<..<harm n - l + inverse (of_nat (2*n))}"
   1.509 +  using euler_mascheroni_bounds[OF assms(1)] assms(2) by auto
   1.510 +
   1.511 +
   1.512 +text \<open>
   1.513 +  Approximation of @{term "ln 2"}. The lower bound is accurate to about 0.03; the upper
   1.514 +  bound is accurate to about 0.0015.
   1.515 +\<close>
   1.516 +lemma ln2_ge_two_thirds: "2/3 \<le> ln (2::real)"
   1.517 +  and ln2_le_25_over_36: "ln (2::real) \<le> 25/36"
   1.518 +  using ln_approx_bounds[of 2 1, simplified, simplified eval_nat_numeral, simplified] by simp_all
   1.519 +
   1.520 +
   1.521 +text \<open>
   1.522 +  Approximation of the Euler--Mascheroni constant. The lower bound is accurate to about 0.0015;
   1.523 +  the upper bound is accurate to about 0.015.
   1.524 +\<close>
   1.525 +lemma euler_mascheroni_gt_19_over_33: "(euler_mascheroni :: real) > 19/33" (is ?th1)
   1.526 +  and euler_mascheroni_less_13_over_22: "(euler_mascheroni :: real) < 13/22" (is ?th2)
   1.527 +proof -
   1.528 +  have "ln (real (Suc 7)) = 3 * ln 2" by (simp add: ln_powr [symmetric] powr_numeral)
   1.529 +  also from ln_approx_bounds[of 2 3] have "\<dots> \<in> {3*307/443<..<3*4615/6658}"
   1.530 +    by (simp add: eval_nat_numeral)
   1.531 +  finally have "ln (real (Suc 7)) \<in> \<dots>" .
   1.532 +  from euler_mascheroni_bounds'[OF _ this] have "?th1 \<and> ?th2" by (simp_all add: harm_expand)
   1.533 +  thus ?th1 ?th2 by blast+
   1.534 +qed
   1.535 +
   1.536 +end