src/HOL/Multivariate_Analysis/Harmonic_Numbers.thy
author eberlm
Mon Jan 04 17:45:36 2016 +0100 (2016-01-04)
changeset 62049 b0f941e207cf
child 62055 755fda743c49
permissions -rw-r--r--
Added lots of material on infinite sums, convergence radii, harmonic numbers, Gamma function
eberlm@62049
     1
(*
eberlm@62049
     2
  Title:    HOL/Multivariate_Analysis/Harmonic_Numbers.thy
eberlm@62049
     3
  Author:   Manuel Eberl, TU München
eberlm@62049
     4
  
eberlm@62049
     5
  The definition of the Harmonic Numbers and the Euler–Mascheroni constant.
eberlm@62049
     6
  Also provides a reasonably accurate approximation of @{term "ln 2 :: real"} 
eberlm@62049
     7
  and the Euler–Mascheroni constant.
eberlm@62049
     8
*)
eberlm@62049
     9
theory Harmonic_Numbers
eberlm@62049
    10
imports 
eberlm@62049
    11
  Complex_Transcendental
eberlm@62049
    12
  Summation
eberlm@62049
    13
  Integral_Test
eberlm@62049
    14
begin
eberlm@62049
    15
eberlm@62049
    16
lemma ln_2_less_1: "ln 2 < (1::real)"
eberlm@62049
    17
proof -
eberlm@62049
    18
  have "2 < 5/(2::real)" by simp
eberlm@62049
    19
  also have "5/2 \<le> exp (1::real)" using exp_lower_taylor_quadratic[of 1, simplified] by simp
eberlm@62049
    20
  finally have "exp (ln 2) < exp (1::real)" by simp
eberlm@62049
    21
  thus "ln 2 < (1::real)" by (subst (asm) exp_less_cancel_iff) simp
eberlm@62049
    22
qed
eberlm@62049
    23
eberlm@62049
    24
lemma setsum_Suc_diff':
eberlm@62049
    25
  fixes f :: "nat \<Rightarrow> 'a::ab_group_add"
eberlm@62049
    26
  assumes "m \<le> n"
eberlm@62049
    27
  shows "(\<Sum>i = m..<n. f(Suc i) - f i) = f n - f m"
eberlm@62049
    28
using assms by (induct n) (auto simp: le_Suc_eq)
eberlm@62049
    29
eberlm@62049
    30
lemma eval_fact:
eberlm@62049
    31
  "fact 0 = 1"
eberlm@62049
    32
  "fact (Suc 0) = 1"
eberlm@62049
    33
  "fact (numeral n) = numeral n * fact (pred_numeral n)"
eberlm@62049
    34
  by (simp, simp, simp_all only: numeral_eq_Suc fact_Suc,
eberlm@62049
    35
      simp only: numeral_eq_Suc [symmetric] of_nat_numeral)
eberlm@62049
    36
eberlm@62049
    37
lemma setsum_poly_horner_expand:
eberlm@62049
    38
  "(\<Sum>k<(numeral n::nat). f k * x^k) = f 0 + (\<Sum>k<pred_numeral n. f (k+1) * x^k) * x"
eberlm@62049
    39
  "(\<Sum>k<Suc 0. f k * x^k) = (f 0 :: 'a :: semiring_1)"
eberlm@62049
    40
  "(\<Sum>k<(0::nat). f k * x^k) = 0"
eberlm@62049
    41
proof -
eberlm@62049
    42
  {
eberlm@62049
    43
    fix m :: nat
eberlm@62049
    44
    have "(\<Sum>k<Suc m. f k * x^k) = f 0 + (\<Sum>k=Suc 0..<Suc m. f k * x^k)"
eberlm@62049
    45
      by (subst atLeast0LessThan [symmetric], subst setsum_head_upt_Suc) simp_all
eberlm@62049
    46
    also have "(\<Sum>k=Suc 0..<Suc m. f k * x^k) = (\<Sum>k<m. f (k+1) * x^k) * x"
eberlm@62049
    47
      by (subst setsum_shift_bounds_Suc_ivl)
eberlm@62049
    48
         (simp add: setsum_left_distrib algebra_simps atLeast0LessThan power_commutes)
eberlm@62049
    49
    finally have "(\<Sum>k<Suc m. f k * x ^ k) = f 0 + (\<Sum>k<m. f (k + 1) * x ^ k) * x" .
eberlm@62049
    50
  }
eberlm@62049
    51
  from this[of "pred_numeral n"] 
eberlm@62049
    52
    show "(\<Sum>k<numeral n. f k * x^k) = f 0 + (\<Sum>k<pred_numeral n. f (k+1) * x^k) * x" 
eberlm@62049
    53
    by (simp add: numeral_eq_Suc)
eberlm@62049
    54
qed simp_all
eberlm@62049
    55
eberlm@62049
    56
eberlm@62049
    57
subsection \<open>The Harmonic numbers\<close>
eberlm@62049
    58
eberlm@62049
    59
definition harm :: "nat \<Rightarrow> 'a :: real_normed_field" where
eberlm@62049
    60
  "harm n = (\<Sum>k=1..n. inverse (of_nat k))"
eberlm@62049
    61
eberlm@62049
    62
lemma harm_altdef: "harm n = (\<Sum>k<n. inverse (of_nat (Suc k)))"
eberlm@62049
    63
  unfolding harm_def by (induction n) simp_all
eberlm@62049
    64
eberlm@62049
    65
lemma harm_Suc: "harm (Suc n) = harm n + inverse (of_nat (Suc n))"
eberlm@62049
    66
  by (simp add: harm_def)
eberlm@62049
    67
eberlm@62049
    68
lemma harm_nonneg: "harm n \<ge> (0 :: 'a :: {real_normed_field,linordered_field})"
eberlm@62049
    69
  unfolding harm_def by (intro setsum_nonneg) simp_all
eberlm@62049
    70
eberlm@62049
    71
lemma of_real_harm: "of_real (harm n) = harm n"
eberlm@62049
    72
  unfolding harm_def by simp
eberlm@62049
    73
  
eberlm@62049
    74
lemma norm_harm: "norm (harm n) = harm n"
eberlm@62049
    75
  by (subst of_real_harm [symmetric]) (simp add: harm_nonneg)
eberlm@62049
    76
eberlm@62049
    77
lemma harm_expand: 
eberlm@62049
    78
  "harm (Suc 0) = 1"
eberlm@62049
    79
  "harm (numeral n) = harm (pred_numeral n) + inverse (numeral n)"
eberlm@62049
    80
proof -
eberlm@62049
    81
  have "numeral n = Suc (pred_numeral n)" by simp
eberlm@62049
    82
  also have "harm \<dots> = harm (pred_numeral n) + inverse (numeral n)"
eberlm@62049
    83
    by (subst harm_Suc, subst numeral_eq_Suc[symmetric]) simp
eberlm@62049
    84
  finally show "harm (numeral n) = harm (pred_numeral n) + inverse (numeral n)" .
eberlm@62049
    85
qed (simp add: harm_def)
eberlm@62049
    86
eberlm@62049
    87
lemma not_convergent_harm: "\<not>convergent (harm :: nat \<Rightarrow> 'a :: real_normed_field)"
eberlm@62049
    88
proof -
eberlm@62049
    89
  have "convergent (\<lambda>n. norm (harm n :: 'a)) \<longleftrightarrow>
eberlm@62049
    90
            convergent (harm :: nat \<Rightarrow> real)" by (simp add: norm_harm)
eberlm@62049
    91
  also have "\<dots> \<longleftrightarrow> convergent (\<lambda>n. \<Sum>k=Suc 0..Suc n. inverse (of_nat k) :: real)"
eberlm@62049
    92
    unfolding harm_def[abs_def] by (subst convergent_Suc_iff) simp_all
eberlm@62049
    93
  also have "... \<longleftrightarrow> convergent (\<lambda>n. \<Sum>k\<le>n. inverse (of_nat (Suc k)) :: real)"
eberlm@62049
    94
    by (subst setsum_shift_bounds_cl_Suc_ivl) (simp add: atLeast0AtMost)
eberlm@62049
    95
  also have "... \<longleftrightarrow> summable (\<lambda>n. inverse (of_nat n) :: real)"
eberlm@62049
    96
    by (subst summable_Suc_iff [symmetric]) (simp add: summable_iff_convergent')
eberlm@62049
    97
  also have "\<not>..." by (rule not_summable_harmonic)
eberlm@62049
    98
  finally show ?thesis by (blast dest: convergent_norm)
eberlm@62049
    99
qed
eberlm@62049
   100
eberlm@62049
   101
eberlm@62049
   102
subsection \<open>The Euler–Mascheroni constant\<close>
eberlm@62049
   103
eberlm@62049
   104
text \<open>
eberlm@62049
   105
  The limit of the difference between the partial harmonic sum and the natural logarithm
eberlm@62049
   106
  (approximately 0.577216). This value occurs e.g. in the definition of the Gamma function.
eberlm@62049
   107
 \<close>
eberlm@62049
   108
definition euler_mascheroni :: "'a :: real_normed_algebra_1" where
eberlm@62049
   109
  "euler_mascheroni = of_real (lim (\<lambda>n. harm n - ln (of_nat n)))"
eberlm@62049
   110
eberlm@62049
   111
lemma of_real_euler_mascheroni [simp]: "of_real euler_mascheroni = euler_mascheroni"
eberlm@62049
   112
  by (simp add: euler_mascheroni_def)
eberlm@62049
   113
eberlm@62049
   114
interpretation euler_mascheroni: antimono_fun_sum_integral_diff "\<lambda>x. inverse (x + 1)"
eberlm@62049
   115
  by unfold_locales (auto intro!: continuous_intros)
eberlm@62049
   116
eberlm@62049
   117
lemma euler_mascheroni_sum_integral_diff_series:
eberlm@62049
   118
  "euler_mascheroni.sum_integral_diff_series n = harm (Suc n) - ln (of_nat (Suc n))"
eberlm@62049
   119
proof -
eberlm@62049
   120
  have "harm (Suc n) = (\<Sum>k=0..n. inverse (of_nat k + 1) :: real)" unfolding harm_def
eberlm@62049
   121
    unfolding One_nat_def by (subst setsum_shift_bounds_cl_Suc_ivl) (simp add: add_ac)
eberlm@62049
   122
  moreover have "((\<lambda>x. inverse (x + 1) :: real) has_integral ln (of_nat n + 1) - ln (0 + 1))
eberlm@62049
   123
                   {0..of_nat n}"
eberlm@62049
   124
    by (intro fundamental_theorem_of_calculus)
eberlm@62049
   125
       (auto intro!: derivative_eq_intros simp: divide_inverse
eberlm@62049
   126
           has_field_derivative_iff_has_vector_derivative[symmetric])
eberlm@62049
   127
  hence "integral {0..of_nat n} (\<lambda>x. inverse (x + 1) :: real) = ln (of_nat (Suc n))"
eberlm@62049
   128
    by (auto dest!: integral_unique)
eberlm@62049
   129
  ultimately show ?thesis 
eberlm@62049
   130
    by (simp add: euler_mascheroni.sum_integral_diff_series_def atLeast0AtMost)
eberlm@62049
   131
qed
eberlm@62049
   132
eberlm@62049
   133
lemma euler_mascheroni_sequence_decreasing:
eberlm@62049
   134
  "m > 0 \<Longrightarrow> m \<le> n \<Longrightarrow> harm n - ln (of_nat n) \<le> harm m - ln (of_nat m :: real)"
eberlm@62049
   135
  by (cases m, simp, cases n, simp, hypsubst,
eberlm@62049
   136
      subst (1 2) euler_mascheroni_sum_integral_diff_series [symmetric],
eberlm@62049
   137
      rule euler_mascheroni.sum_integral_diff_series_antimono, simp)
eberlm@62049
   138
eberlm@62049
   139
lemma euler_mascheroni_sequence_nonneg:
eberlm@62049
   140
  "n > 0 \<Longrightarrow> harm n - ln (of_nat n) \<ge> (0::real)"
eberlm@62049
   141
  by (cases n, simp, hypsubst, subst euler_mascheroni_sum_integral_diff_series [symmetric],
eberlm@62049
   142
      rule euler_mascheroni.sum_integral_diff_series_nonneg)
eberlm@62049
   143
eberlm@62049
   144
lemma euler_mascheroni_convergent: "convergent (\<lambda>n. harm n - ln (of_nat n) :: real)"
eberlm@62049
   145
proof -
eberlm@62049
   146
  have A: "(\<lambda>n. harm (Suc n) - ln (of_nat (Suc n))) = 
eberlm@62049
   147
             euler_mascheroni.sum_integral_diff_series"
eberlm@62049
   148
    by (subst euler_mascheroni_sum_integral_diff_series [symmetric]) (rule refl)
eberlm@62049
   149
  have "convergent (\<lambda>n. harm (Suc n) - ln (of_nat (Suc n) :: real))"
eberlm@62049
   150
    by (subst A) (fact euler_mascheroni.sum_integral_diff_series_convergent)
eberlm@62049
   151
  thus ?thesis by (subst (asm) convergent_Suc_iff)
eberlm@62049
   152
qed
eberlm@62049
   153
eberlm@62049
   154
lemma euler_mascheroni_LIMSEQ: 
eberlm@62049
   155
  "(\<lambda>n. harm n - ln (of_nat n) :: real) \<longlonglongrightarrow> euler_mascheroni"
eberlm@62049
   156
  unfolding euler_mascheroni_def
eberlm@62049
   157
  by (simp add: convergent_LIMSEQ_iff [symmetric] euler_mascheroni_convergent)
eberlm@62049
   158
eberlm@62049
   159
lemma euler_mascheroni_LIMSEQ_of_real: 
eberlm@62049
   160
  "(\<lambda>n. of_real (harm n - ln (of_nat n))) \<longlonglongrightarrow> 
eberlm@62049
   161
      (euler_mascheroni :: 'a :: {real_normed_algebra_1, topological_space})"
eberlm@62049
   162
proof -
eberlm@62049
   163
  have "(\<lambda>n. of_real (harm n - ln (of_nat n))) \<longlonglongrightarrow> (of_real (euler_mascheroni) :: 'a)"
eberlm@62049
   164
    by (intro tendsto_of_real euler_mascheroni_LIMSEQ)
eberlm@62049
   165
  thus ?thesis by simp
eberlm@62049
   166
qed
eberlm@62049
   167
eberlm@62049
   168
lemma euler_mascheroni_sum:
eberlm@62049
   169
  "(\<lambda>n. inverse (of_nat (n+1)) + ln (of_nat (n+1)) - ln (of_nat (n+2)) :: real)
eberlm@62049
   170
       sums euler_mascheroni" 
eberlm@62049
   171
 using sums_add[OF telescope_sums[OF LIMSEQ_Suc[OF euler_mascheroni_LIMSEQ]]
eberlm@62049
   172
                   telescope_sums'[OF LIMSEQ_inverse_real_of_nat]]
eberlm@62049
   173
  by (simp_all add: harm_def algebra_simps)
eberlm@62049
   174
eberlm@62049
   175
lemma alternating_harmonic_series_sums: "(\<lambda>k. (-1)^k / real_of_nat (Suc k)) sums ln 2"
eberlm@62049
   176
proof -
eberlm@62049
   177
  let ?f = "\<lambda>n. harm n - ln (real_of_nat n)"
eberlm@62049
   178
  let ?g = "\<lambda>n. if even n then 0 else (2::real)"
eberlm@62049
   179
  let ?em = "\<lambda>n. harm n - ln (real_of_nat n)"
eberlm@62049
   180
  have "eventually (\<lambda>n. ?em (2*n) - ?em n + ln 2 = (\<Sum>k<2*n. (-1)^k / real_of_nat (Suc k))) at_top"
eberlm@62049
   181
    using eventually_gt_at_top[of "0::nat"]
eberlm@62049
   182
  proof eventually_elim
eberlm@62049
   183
    fix n :: nat assume n: "n > 0"
eberlm@62049
   184
    have "(\<Sum>k<2*n. (-1)^k / real_of_nat (Suc k)) =
eberlm@62049
   185
              (\<Sum>k<2*n. ((-1)^k + ?g k) / of_nat (Suc k)) - (\<Sum>k<2*n. ?g k / of_nat (Suc k))"
eberlm@62049
   186
      by (simp add: setsum.distrib algebra_simps divide_inverse)
eberlm@62049
   187
    also have "(\<Sum>k<2*n. ((-1)^k + ?g k) / real_of_nat (Suc k)) = harm (2*n)"
eberlm@62049
   188
      unfolding harm_altdef by (intro setsum.cong) (auto simp: field_simps)
eberlm@62049
   189
    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))"
eberlm@62049
   190
      by (intro setsum.mono_neutral_right) auto
eberlm@62049
   191
    also have "\<dots> = (\<Sum>k|k<2*n \<and> odd k. 2 / (real_of_nat (Suc k)))"
eberlm@62049
   192
      by (intro setsum.cong) auto
eberlm@62049
   193
    also have "(\<Sum>k|k<2*n \<and> odd k. 2 / (real_of_nat (Suc k))) = harm n" 
eberlm@62049
   194
      unfolding harm_altdef
eberlm@62049
   195
      by (intro setsum.reindex_cong[of "\<lambda>n. 2*n+1"]) (auto simp: inj_on_def field_simps elim!: oddE)
eberlm@62049
   196
    also have "harm (2*n) - harm n = ?em (2*n) - ?em n + ln 2" using n
eberlm@62049
   197
      by (simp_all add: algebra_simps ln_mult)
eberlm@62049
   198
    finally show "?em (2*n) - ?em n + ln 2 = (\<Sum>k<2*n. (-1)^k / real_of_nat (Suc k))" ..
eberlm@62049
   199
  qed
eberlm@62049
   200
  moreover have "(\<lambda>n. ?em (2*n) - ?em n + ln (2::real)) 
eberlm@62049
   201
                     \<longlonglongrightarrow> euler_mascheroni - euler_mascheroni + ln 2"
eberlm@62049
   202
    by (intro tendsto_intros euler_mascheroni_LIMSEQ filterlim_compose[OF euler_mascheroni_LIMSEQ]
eberlm@62049
   203
              filterlim_subseq) (auto simp: subseq_def)
eberlm@62049
   204
  hence "(\<lambda>n. ?em (2*n) - ?em n + ln (2::real)) \<longlonglongrightarrow> ln 2" by simp
eberlm@62049
   205
  ultimately have "(\<lambda>n. (\<Sum>k<2*n. (-1)^k / real_of_nat (Suc k))) \<longlonglongrightarrow> ln 2"
eberlm@62049
   206
    by (rule Lim_transform_eventually)
eberlm@62049
   207
  
eberlm@62049
   208
  moreover have "summable (\<lambda>k. (-1)^k * inverse (real_of_nat (Suc k)))"
eberlm@62049
   209
    using LIMSEQ_inverse_real_of_nat
eberlm@62049
   210
    by (intro summable_Leibniz(1) decseq_imp_monoseq decseq_SucI) simp_all
eberlm@62049
   211
  hence A: "(\<lambda>n. \<Sum>k<n. (-1)^k / real_of_nat (Suc k)) \<longlonglongrightarrow> (\<Sum>k. (-1)^k / real_of_nat (Suc k))"
eberlm@62049
   212
    by (simp add: summable_sums_iff divide_inverse sums_def)
eberlm@62049
   213
  from filterlim_compose[OF this filterlim_subseq[of "op * (2::nat)"]]
eberlm@62049
   214
    have "(\<lambda>n. \<Sum>k<2*n. (-1)^k / real_of_nat (Suc k)) \<longlonglongrightarrow> (\<Sum>k. (-1)^k / real_of_nat (Suc k))"
eberlm@62049
   215
    by (simp add: subseq_def)
eberlm@62049
   216
  ultimately have "(\<Sum>k. (- 1) ^ k / real_of_nat (Suc k)) = ln 2" by (intro LIMSEQ_unique)
eberlm@62049
   217
  with A show ?thesis by (simp add: sums_def)
eberlm@62049
   218
qed
eberlm@62049
   219
eberlm@62049
   220
lemma alternating_harmonic_series_sums': 
eberlm@62049
   221
  "(\<lambda>k. inverse (real_of_nat (2*k+1)) - inverse (real_of_nat (2*k+2))) sums ln 2"
eberlm@62049
   222
unfolding sums_def
eberlm@62049
   223
proof (rule Lim_transform_eventually)
eberlm@62049
   224
  show "(\<lambda>n. \<Sum>k<2*n. (-1)^k / (real_of_nat (Suc k))) \<longlonglongrightarrow> ln 2"
eberlm@62049
   225
    using alternating_harmonic_series_sums unfolding sums_def 
eberlm@62049
   226
    by (rule filterlim_compose) (rule mult_nat_left_at_top, simp)
eberlm@62049
   227
  show "eventually (\<lambda>n. (\<Sum>k<2*n. (-1)^k / (real_of_nat (Suc k))) =
eberlm@62049
   228
            (\<Sum>k<n. inverse (real_of_nat (2*k+1)) - inverse (real_of_nat (2*k+2)))) sequentially"
eberlm@62049
   229
  proof (intro always_eventually allI)
eberlm@62049
   230
    fix n :: nat
eberlm@62049
   231
    show "(\<Sum>k<2*n. (-1)^k / (real_of_nat (Suc k))) =
eberlm@62049
   232
              (\<Sum>k<n. inverse (real_of_nat (2*k+1)) - inverse (real_of_nat (2*k+2)))"
eberlm@62049
   233
      by (induction n) (simp_all add: inverse_eq_divide)
eberlm@62049
   234
  qed
eberlm@62049
   235
qed               
eberlm@62049
   236
eberlm@62049
   237
eberlm@62049
   238
subsection \<open>Approximation of the Euler--Mascheroni constant\<close>
eberlm@62049
   239
eberlm@62049
   240
(* FIXME: ugly *)
eberlm@62049
   241
(* TODO: Move ? *)
eberlm@62049
   242
lemma ln_inverse_approx_le:
eberlm@62049
   243
  assumes "(x::real) > 0" "a > 0"
eberlm@62049
   244
  shows   "ln (x + a) - ln x \<le> a * (inverse x + inverse (x + a))/2" (is "_ \<le> ?A")
eberlm@62049
   245
proof -
eberlm@62049
   246
  def f' \<equiv> "(inverse (x + a) - inverse x)/a"
eberlm@62049
   247
  have f'_nonpos: "f' \<le> 0" using assms by (simp add: f'_def divide_simps)
eberlm@62049
   248
  let ?f = "\<lambda>t. (t - x) * f' + inverse x"
eberlm@62049
   249
  let ?F = "\<lambda>t. (t - x)^2 * f' / 2 + t * inverse x"
eberlm@62049
   250
  have diff: "\<forall>t\<in>{x..x+a}. (?F has_vector_derivative ?f t) 
eberlm@62049
   251
                               (at t within {x..x+a})" using assms
eberlm@62049
   252
    by (auto intro!: derivative_eq_intros 
eberlm@62049
   253
             simp: has_field_derivative_iff_has_vector_derivative[symmetric])
eberlm@62049
   254
  from assms have "(?f has_integral (?F (x+a) - ?F x)) {x..x+a}"
eberlm@62049
   255
    by (intro fundamental_theorem_of_calculus[OF _ diff])
eberlm@62049
   256
       (auto simp: has_field_derivative_iff_has_vector_derivative[symmetric] field_simps
eberlm@62049
   257
             intro!: derivative_eq_intros)
eberlm@62049
   258
  also have "?F (x+a) - ?F x = (a*2 + f'*a\<^sup>2*x) / (2*x)" using assms by (simp add: field_simps)
eberlm@62049
   259
  also have "f'*a^2 = - (a^2) / (x*(x + a))" using assms 
eberlm@62049
   260
    by (simp add: divide_simps f'_def power2_eq_square)
eberlm@62049
   261
  also have "(a*2 + - a\<^sup>2/(x*(x+a))*x) / (2*x) = ?A" using assms
eberlm@62049
   262
    by (simp add: divide_simps power2_eq_square) (simp add: algebra_simps)
eberlm@62049
   263
  finally have int1: "((\<lambda>t. (t - x) * f' + inverse x) has_integral ?A) {x..x + a}" .
eberlm@62049
   264
eberlm@62049
   265
  from assms have int2: "(inverse has_integral (ln (x + a) - ln x)) {x..x+a}"
eberlm@62049
   266
    by (intro fundamental_theorem_of_calculus)
eberlm@62049
   267
       (auto simp: has_field_derivative_iff_has_vector_derivative[symmetric] divide_simps
eberlm@62049
   268
             intro!: derivative_eq_intros)
eberlm@62049
   269
  hence "ln (x + a) - ln x = integral {x..x+a} inverse" by (simp add: integral_unique)
eberlm@62049
   270
  also have ineq: "\<forall>xa\<in>{x..x + a}. inverse xa \<le> (xa - x) * f' + inverse x"
eberlm@62049
   271
  proof
eberlm@62049
   272
    fix t assume t': "t \<in> {x..x+a}"
eberlm@62049
   273
    with assms have t: "0 \<le> (t - x) / a" "(t - x) / a \<le> 1" by simp_all
eberlm@62049
   274
    have "inverse t = inverse ((1 - (t - x) / a) *\<^sub>R x + ((t - x) / a) *\<^sub>R (x + a))" (is "_ = ?A")
eberlm@62049
   275
      using assms t' by (simp add: field_simps)
eberlm@62049
   276
    also from assms have "convex_on {x..x+a} inverse" by (intro convex_on_inverse) auto
eberlm@62049
   277
    from convex_onD_Icc[OF this _ t] assms 
eberlm@62049
   278
      have "?A \<le> (1 - (t - x) / a) * inverse x + (t - x) / a * inverse (x + a)" by simp
eberlm@62049
   279
    also have "\<dots> = (t - x) * f' + inverse x" using assms
eberlm@62049
   280
      by (simp add: f'_def divide_simps) (simp add: f'_def field_simps)
eberlm@62049
   281
    finally show "inverse t \<le> (t - x) * f' + inverse x" .
eberlm@62049
   282
  qed
eberlm@62049
   283
  hence "integral {x..x+a} inverse \<le> integral {x..x+a} ?f" using f'_nonpos assms
eberlm@62049
   284
    by (intro integral_le has_integral_integrable[OF int1] has_integral_integrable[OF int2] ineq)
eberlm@62049
   285
  also have "\<dots> = ?A" using int1 by (rule integral_unique)
eberlm@62049
   286
  finally show ?thesis .
eberlm@62049
   287
qed
eberlm@62049
   288
eberlm@62049
   289
lemma ln_inverse_approx_ge:
eberlm@62049
   290
  assumes "(x::real) > 0" "x < y"
eberlm@62049
   291
  shows   "ln y - ln x \<ge> 2 * (y - x) / (x + y)" (is "_ \<ge> ?A")
eberlm@62049
   292
proof -
eberlm@62049
   293
  def m \<equiv> "(x+y)/2"
eberlm@62049
   294
  def f' \<equiv> "-inverse (m^2)"
eberlm@62049
   295
  from assms have m: "m > 0" by (simp add: m_def)
eberlm@62049
   296
  let ?F = "\<lambda>t. (t - m)^2 * f' / 2 + t / m"
eberlm@62049
   297
  from assms have "((\<lambda>t. (t - m) * f' + inverse m) has_integral (?F y - ?F x)) {x..y}"
eberlm@62049
   298
    by (intro fundamental_theorem_of_calculus)
eberlm@62049
   299
       (auto simp: has_field_derivative_iff_has_vector_derivative[symmetric] divide_simps
eberlm@62049
   300
             intro!: derivative_eq_intros)
eberlm@62049
   301
  also from m have "?F y - ?F x = ((y - m)^2 - (x - m)^2) * f' / 2 + (y - x) / m" 
eberlm@62049
   302
    by (simp add: field_simps)
eberlm@62049
   303
  also have "((y - m)^2 - (x - m)^2) = 0" by (simp add: m_def power2_eq_square field_simps)
eberlm@62049
   304
  also have "0 * f' / 2 + (y - x) / m = ?A" by (simp add: m_def)
eberlm@62049
   305
  finally have int1: "((\<lambda>t. (t - m) * f' + inverse m) has_integral ?A) {x..y}" .
eberlm@62049
   306
eberlm@62049
   307
  from assms have int2: "(inverse has_integral (ln y - ln x)) {x..y}"
eberlm@62049
   308
    by (intro fundamental_theorem_of_calculus)
eberlm@62049
   309
       (auto simp: has_field_derivative_iff_has_vector_derivative[symmetric] divide_simps
eberlm@62049
   310
             intro!: derivative_eq_intros)
eberlm@62049
   311
  hence "ln y - ln x = integral {x..y} inverse" by (simp add: integral_unique)
eberlm@62049
   312
  also have ineq: "\<forall>xa\<in>{x..y}. inverse xa \<ge> (xa - m) * f' + inverse m"
eberlm@62049
   313
  proof
eberlm@62049
   314
    fix t assume t: "t \<in> {x..y}"
eberlm@62049
   315
    from t assms have "inverse t - inverse m \<ge> f' * (t - m)"
eberlm@62049
   316
      by (intro convex_on_imp_above_tangent[of "{0<..}"] convex_on_inverse)
eberlm@62049
   317
         (auto simp: m_def interior_open f'_def power2_eq_square intro!: derivative_eq_intros)
eberlm@62049
   318
    thus "(t - m) * f' + inverse m \<le> inverse t" by (simp add: algebra_simps)
eberlm@62049
   319
  qed
eberlm@62049
   320
  hence "integral {x..y} inverse \<ge> integral {x..y} (\<lambda>t. (t - m) * f' + inverse m)"
eberlm@62049
   321
    using int1 int2 by (intro integral_le has_integral_integrable)
eberlm@62049
   322
  also have "integral {x..y} (\<lambda>t. (t - m) * f' + inverse m) = ?A"
eberlm@62049
   323
    using integral_unique[OF int1] by simp
eberlm@62049
   324
  finally show ?thesis .
eberlm@62049
   325
qed
eberlm@62049
   326
eberlm@62049
   327
eberlm@62049
   328
lemma euler_mascheroni_lower: 
eberlm@62049
   329
        "euler_mascheroni \<ge> harm (Suc n) - ln (real_of_nat (n + 2)) + 1/real_of_nat (2 * (n + 2))"
eberlm@62049
   330
  and euler_mascheroni_upper:
eberlm@62049
   331
        "euler_mascheroni \<le> harm (Suc n) - ln (real_of_nat (n + 2)) + 1/real_of_nat (2 * (n + 1))"
eberlm@62049
   332
proof -
eberlm@62049
   333
  def D \<equiv> "\<lambda>n. inverse (of_nat (n+1)) + ln (of_nat (n+1)) - ln (of_nat (n+2)) :: real"
eberlm@62049
   334
  let ?g = "\<lambda>n. ln (of_nat (n+2)) - ln (of_nat (n+1)) - inverse (of_nat (n+1)) :: real"
eberlm@62049
   335
  def inv \<equiv> "\<lambda>n. inverse (real_of_nat n)"
eberlm@62049
   336
  fix n :: nat
eberlm@62049
   337
  note summable = sums_summable[OF euler_mascheroni_sum, folded D_def]
eberlm@62049
   338
  have sums: "(\<lambda>k. (inv (Suc (k + (n+1))) - inv (Suc (Suc k + (n+1))))/2) sums ((inv (Suc (0 + (n+1))) - 0)/2)"
eberlm@62049
   339
    unfolding inv_def
eberlm@62049
   340
    by (intro sums_divide telescope_sums' LIMSEQ_ignore_initial_segment LIMSEQ_inverse_real_of_nat)
eberlm@62049
   341
  have sums': "(\<lambda>k. (inv (Suc (k + n)) - inv (Suc (Suc k + n)))/2) sums ((inv (Suc (0 + n)) - 0)/2)"
eberlm@62049
   342
    unfolding inv_def
eberlm@62049
   343
    by (intro sums_divide telescope_sums' LIMSEQ_ignore_initial_segment LIMSEQ_inverse_real_of_nat)
eberlm@62049
   344
  from euler_mascheroni_sum have "euler_mascheroni = (\<Sum>k. D k)"
eberlm@62049
   345
    by (simp add: sums_iff D_def)
eberlm@62049
   346
  also have "\<dots> = (\<Sum>k. D (k + Suc n)) + (\<Sum>k\<le>n. D k)"
eberlm@62049
   347
    by (subst suminf_split_initial_segment[OF summable, of "Suc n"], subst lessThan_Suc_atMost) simp
eberlm@62049
   348
  finally have sum: "(\<Sum>k\<le>n. D k) - euler_mascheroni = -(\<Sum>k. D (k + Suc n))" by simp
eberlm@62049
   349
eberlm@62049
   350
  note sum
eberlm@62049
   351
  also have "\<dots> \<le> -(\<Sum>k. (inv (k + Suc n + 1) - inv (k + Suc n + 2)) / 2)"
eberlm@62049
   352
  proof (intro le_imp_neg_le suminf_le allI summable_ignore_initial_segment[OF summable])
eberlm@62049
   353
    fix k' :: nat
eberlm@62049
   354
    def k \<equiv> "k' + Suc n"
eberlm@62049
   355
    hence k: "k > 0" by (simp add: k_def)
eberlm@62049
   356
    have "real_of_nat (k+1) > 0" by (simp add: k_def)
eberlm@62049
   357
    with ln_inverse_approx_le[OF this zero_less_one]
eberlm@62049
   358
      have "ln (of_nat k + 2) - ln (of_nat k + 1) \<le> (inv (k+1) + inv (k+2))/2"
eberlm@62049
   359
      by (simp add: inv_def add_ac)
eberlm@62049
   360
    hence "(inv (k+1) - inv (k+2))/2 \<le> inv (k+1) + ln (of_nat (k+1)) - ln (of_nat (k+2))"
eberlm@62049
   361
      by (simp add: field_simps)
eberlm@62049
   362
    also have "\<dots> = D k" unfolding D_def inv_def ..
eberlm@62049
   363
    finally show "D (k' + Suc n) \<ge> (inv (k' + Suc n + 1) - inv (k' + Suc n + 2)) / 2"
eberlm@62049
   364
      by (simp add: k_def)
eberlm@62049
   365
    from sums_summable[OF sums] 
eberlm@62049
   366
      show "summable (\<lambda>k. (inv (k + Suc n + 1) - inv (k + Suc n + 2))/2)" by simp
eberlm@62049
   367
  qed
eberlm@62049
   368
  also from sums have "\<dots> = -inv (n+2) / 2" by (simp add: sums_iff)
eberlm@62049
   369
  finally have "euler_mascheroni \<ge> (\<Sum>k\<le>n. D k) + 1 / (of_nat (2 * (n+2)))" 
eberlm@62049
   370
    by (simp add: inv_def field_simps of_nat_mult)
eberlm@62049
   371
  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)))"
eberlm@62049
   372
    unfolding harm_altdef D_def by (subst lessThan_Suc_atMost) (simp add:  setsum.distrib setsum_subtractf)
eberlm@62049
   373
  also have "(\<Sum>k\<le>n. ln (real_of_nat (Suc k+1)) - ln (of_nat (k+1))) = ln (of_nat (n+2))"
eberlm@62049
   374
    by (subst atLeast0AtMost [symmetric], subst setsum_Suc_diff) simp_all
eberlm@62049
   375
  finally show "euler_mascheroni \<ge> harm (Suc n) - ln (real_of_nat (n + 2)) + 1/real_of_nat (2 * (n + 2))"
eberlm@62049
   376
    by simp
eberlm@62049
   377
  
eberlm@62049
   378
  note sum
eberlm@62049
   379
  also have "-(\<Sum>k. D (k + Suc n)) \<ge> -(\<Sum>k. (inv (Suc (k + n)) - inv (Suc (Suc k + n)))/2)"
eberlm@62049
   380
  proof (intro le_imp_neg_le suminf_le allI summable_ignore_initial_segment[OF summable])
eberlm@62049
   381
    fix k' :: nat
eberlm@62049
   382
    def k \<equiv> "k' + Suc n"
eberlm@62049
   383
    hence k: "k > 0" by (simp add: k_def)
eberlm@62049
   384
    have "real_of_nat (k+1) > 0" by (simp add: k_def)
eberlm@62049
   385
    from ln_inverse_approx_ge[of "of_nat k + 1" "of_nat k + 2"]
eberlm@62049
   386
      have "2 / (2 * real_of_nat k + 3) \<le> ln (of_nat (k+2)) - ln (real_of_nat (k+1))"
eberlm@62049
   387
      by (simp add: add_ac)
eberlm@62049
   388
    hence "D k \<le> 1 / real_of_nat (k+1) - 2 / (2 * real_of_nat k + 3)" 
eberlm@62049
   389
      by (simp add: D_def inverse_eq_divide inv_def)
eberlm@62049
   390
    also have "\<dots> = inv ((k+1)*(2*k+3))" unfolding inv_def by (simp add: field_simps)
eberlm@62049
   391
    also have "\<dots> \<le> inv (2*k*(k+1))" unfolding inv_def using k
eberlm@62049
   392
      by (intro le_imp_inverse_le) 
eberlm@62049
   393
         (simp add: algebra_simps, simp del: of_nat_add)
eberlm@62049
   394
    also have "\<dots> = (inv k - inv (k+1))/2" unfolding inv_def using k
eberlm@62049
   395
      by (simp add: divide_simps del: of_nat_mult) (simp add: algebra_simps)
eberlm@62049
   396
    finally show "D k \<le> (inv (Suc (k' + n)) - inv (Suc (Suc k' + n)))/2" unfolding k_def by simp
eberlm@62049
   397
  next
eberlm@62049
   398
    from sums_summable[OF sums'] 
eberlm@62049
   399
      show "summable (\<lambda>k. (inv (Suc (k + n)) - inv (Suc (Suc k + n)))/2)" by simp
eberlm@62049
   400
  qed
eberlm@62049
   401
  also from sums' have "(\<Sum>k. (inv (Suc (k + n)) - inv (Suc (Suc k + n)))/2) = inv (n+1)/2"
eberlm@62049
   402
    by (simp add: sums_iff)
eberlm@62049
   403
  finally have "euler_mascheroni \<le> (\<Sum>k\<le>n. D k) + 1 / of_nat (2 * (n+1))" 
eberlm@62049
   404
    by (simp add: inv_def field_simps)
eberlm@62049
   405
  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)))"
eberlm@62049
   406
    unfolding harm_altdef D_def by (subst lessThan_Suc_atMost) (simp add:  setsum.distrib setsum_subtractf)
eberlm@62049
   407
  also have "(\<Sum>k\<le>n. ln (real_of_nat (Suc k+1)) - ln (of_nat (k+1))) = ln (of_nat (n+2))"
eberlm@62049
   408
    by (subst atLeast0AtMost [symmetric], subst setsum_Suc_diff) simp_all
eberlm@62049
   409
  finally show "euler_mascheroni \<le> harm (Suc n) - ln (real_of_nat (n + 2)) + 1/real_of_nat (2 * (n + 1))"
eberlm@62049
   410
    by simp
eberlm@62049
   411
qed
eberlm@62049
   412
eberlm@62049
   413
lemma euler_mascheroni_pos: "euler_mascheroni > (0::real)"
eberlm@62049
   414
  using euler_mascheroni_lower[of 0] ln_2_less_1 by (simp add: harm_def)
eberlm@62049
   415
eberlm@62049
   416
lemma ln_approx_aux:
eberlm@62049
   417
  fixes n :: nat and x :: real
eberlm@62049
   418
  defines "y \<equiv> (x-1)/(x+1)"
eberlm@62049
   419
  assumes x: "x > 0" "x \<noteq> 1"
eberlm@62049
   420
  shows "inverse (2*y^(2*n+1)) * (ln x - (\<Sum>k<n. 2*y^(2*k+1) / of_nat (2*k+1))) \<in> 
eberlm@62049
   421
            {0..(1 / (1 - y^2) / of_nat (2*n+1))}"
eberlm@62049
   422
proof -
eberlm@62049
   423
  from x have norm_y: "norm y < 1" unfolding y_def by simp
eberlm@62049
   424
  from power_strict_mono[OF this, of 2] have norm_y': "norm y^2 < 1" by simp
eberlm@62049
   425
eberlm@62049
   426
  let ?f = "\<lambda>k. 2 * y ^ (2*k+1) / of_nat (2*k+1)"
eberlm@62049
   427
  note sums = ln_series_quadratic[OF x(1)]
eberlm@62049
   428
  def c \<equiv> "inverse (2*y^(2*n+1))"
eberlm@62049
   429
  let ?d = "c * (ln x - (\<Sum>k<n. ?f k))"
eberlm@62049
   430
  have "\<forall>k. y\<^sup>2^k / of_nat (2*(k+n)+1) \<le> y\<^sup>2 ^ k / of_nat (2*n+1)"
eberlm@62049
   431
    by (intro allI divide_left_mono mult_right_mono mult_pos_pos zero_le_power[of "y^2"]) simp_all
eberlm@62049
   432
  moreover {
eberlm@62049
   433
    have "(\<lambda>k. ?f (k + n)) sums (ln x - (\<Sum>k<n. ?f k))"
eberlm@62049
   434
      using sums_split_initial_segment[OF sums] by (simp add: y_def)
eberlm@62049
   435
    hence "(\<lambda>k. c * ?f (k + n)) sums ?d" by (rule sums_mult)
eberlm@62049
   436
    also have "(\<lambda>k. c * (2*y^(2*(k+n)+1) / of_nat (2*(k+n)+1))) =
eberlm@62049
   437
                   (\<lambda>k. (c * (2*y^(2*n+1))) * ((y^2)^k / of_nat (2*(k+n)+1)))"
eberlm@62049
   438
      by (simp only: ring_distribs power_add power_mult) (simp add: mult_ac)
eberlm@62049
   439
    also from x have "c * (2*y^(2*n+1)) = 1" by (simp add: c_def y_def)
eberlm@62049
   440
    finally have "(\<lambda>k. (y^2)^k / of_nat (2*(k+n)+1)) sums ?d" by simp
eberlm@62049
   441
  } note sums' = this
eberlm@62049
   442
  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))"
eberlm@62049
   443
    by (intro sums_divide geometric_sums) (simp_all add: norm_power)
eberlm@62049
   444
  ultimately have "?d \<le> (1 / (1 - y^2) / of_nat (2*n+1))" by (rule sums_le)
eberlm@62049
   445
  moreover have "c * (ln x - (\<Sum>k<n. 2 * y ^ (2 * k + 1) / real_of_nat (2 * k + 1))) \<ge> 0"
eberlm@62049
   446
    by (intro sums_le[OF _ sums_zero sums']) simp_all
eberlm@62049
   447
  ultimately show ?thesis unfolding c_def by simp
eberlm@62049
   448
qed
eberlm@62049
   449
eberlm@62049
   450
lemma
eberlm@62049
   451
  fixes n :: nat and x :: real
eberlm@62049
   452
  defines "y \<equiv> (x-1)/(x+1)"
eberlm@62049
   453
  defines "approx \<equiv> (\<Sum>k<n. 2*y^(2*k+1) / of_nat (2*k+1))"
eberlm@62049
   454
  defines "d \<equiv> y^(2*n+1) / (1 - y^2) / of_nat (2*n+1)"
eberlm@62049
   455
  assumes x: "x > 1"
eberlm@62049
   456
  shows   ln_approx_bounds: "ln x \<in> {approx..approx + 2*d}"
eberlm@62049
   457
  and     ln_approx_abs:    "abs (ln x - (approx + d)) \<le> d"
eberlm@62049
   458
proof -
eberlm@62049
   459
  def c \<equiv> "2*y^(2*n+1)"
eberlm@62049
   460
  from x have c_pos: "c > 0" unfolding c_def y_def 
eberlm@62049
   461
    by (intro mult_pos_pos zero_less_power) simp_all
eberlm@62049
   462
  have A: "inverse c * (ln x - (\<Sum>k<n. 2*y^(2*k+1) / of_nat (2*k+1))) \<in>
eberlm@62049
   463
              {0.. (1 / (1 - y^2) / of_nat (2*n+1))}" using assms unfolding y_def c_def
eberlm@62049
   464
    by (intro ln_approx_aux) simp_all
eberlm@62049
   465
  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))"
eberlm@62049
   466
    by simp
eberlm@62049
   467
  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))" 
eberlm@62049
   468
    by (auto simp add: divide_simps)
eberlm@62049
   469
  with c_pos have "ln x \<le> c / (1 - y^2) / of_nat (2*n+1) + approx"
eberlm@62049
   470
    by (subst (asm) pos_divide_le_eq) (simp_all add: mult_ac approx_def)
eberlm@62049
   471
  moreover {
eberlm@62049
   472
    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))))"
eberlm@62049
   473
      by (intro mult_nonneg_nonneg[of c]) simp_all
eberlm@62049
   474
    also have "\<dots> = (c * inverse c) * (ln x - (\<Sum>k<n. 2*y^(2*k+1) / of_nat (2*k+1)))"  
eberlm@62049
   475
      by (simp add: mult_ac)
eberlm@62049
   476
    also from c_pos have "c * inverse c = 1" by simp
eberlm@62049
   477
    finally have "ln x \<ge> approx" by (simp add: approx_def)
eberlm@62049
   478
  }
eberlm@62049
   479
  ultimately show "ln x \<in> {approx..approx + 2*d}" by (simp add: c_def d_def)
eberlm@62049
   480
  thus "abs (ln x - (approx + d)) \<le> d" by auto
eberlm@62049
   481
qed
eberlm@62049
   482
eberlm@62049
   483
context
eberlm@62049
   484
begin
eberlm@62049
   485
eberlm@62049
   486
qualified lemma ln_approx_abs': 
eberlm@62049
   487
  assumes "x > (1::real)"
eberlm@62049
   488
  assumes "(x-1)/(x+1) = y"
eberlm@62049
   489
  assumes "y^2 = ysqr"
eberlm@62049
   490
  assumes "(\<Sum>k<n. inverse (of_nat (2*k+1)) * ysqr^k) = approx"
eberlm@62049
   491
  assumes "y*ysqr^n / (1 - ysqr) / of_nat (2*n+1) = d"
eberlm@62049
   492
  assumes "d \<le> e"
eberlm@62049
   493
  shows   "abs (ln x - (2*y*approx + d)) \<le> e"
eberlm@62049
   494
proof -
eberlm@62049
   495
  note ln_approx_abs[OF assms(1), of n]
eberlm@62049
   496
  also note assms(2)
eberlm@62049
   497
  also have "y^(2*n+1) = y*ysqr^n" by (simp add: assms(3)[symmetric] power_mult)
eberlm@62049
   498
  also note assms(3)
eberlm@62049
   499
  also note assms(5)
eberlm@62049
   500
  also note assms(5)
eberlm@62049
   501
  also note assms(6)
eberlm@62049
   502
  also have "(\<Sum>k<n. 2*y^(2*k+1) / real_of_nat (2 * k + 1)) = (2*y) * approx"
eberlm@62049
   503
    apply (subst assms(4)[symmetric], subst setsum_right_distrib)
eberlm@62049
   504
    apply (simp add: assms(3)[symmetric] power_mult)
eberlm@62049
   505
    apply (simp add: mult_ac divide_simps)?
eberlm@62049
   506
    done
eberlm@62049
   507
  finally show ?thesis .
eberlm@62049
   508
qed
eberlm@62049
   509
eberlm@62049
   510
lemma ln_2_approx: "\<bar>ln 2 - 0.69314718055\<bar> < inverse (2 ^ 36 :: real)" (is ?thesis1)
eberlm@62049
   511
  and ln_2_bounds: "ln (2::real) \<in> {0.693147180549..0.693147180561}" (is ?thesis2)
eberlm@62049
   512
proof -
eberlm@62049
   513
  def approx \<equiv> "0.69314718055 :: real" and approx' \<equiv> "4465284211343447 / 6442043387911560 :: real"
eberlm@62049
   514
  def d \<equiv> "inverse (195259926456::real)"
eberlm@62049
   515
  have "dist (ln 2) approx \<le> dist (ln 2) approx' + dist approx' approx" by (rule dist_triangle)
eberlm@62049
   516
  also have "\<bar>ln (2::real) - (2 * (1/3) * (651187280816108 / 626309773824735) +
eberlm@62049
   517
                 inverse 195259926456)\<bar> \<le> inverse 195259926456"
eberlm@62049
   518
  proof (rule ln_approx_abs'[where n = 10])
eberlm@62049
   519
    show "(1/3::real)^2 = 1/9" by (simp add: power2_eq_square)
eberlm@62049
   520
  qed (simp_all add: eval_nat_numeral)
eberlm@62049
   521
  hence A: "dist (ln 2) approx' \<le> d" by (simp add: dist_real_def approx'_def d_def)
eberlm@62049
   522
  hence "dist (ln 2) approx' + dist approx' approx \<le> \<dots> + dist approx' approx"
eberlm@62049
   523
    by (rule add_right_mono)
eberlm@62049
   524
  also have "\<dots> < inverse (2 ^ 36)" by (simp add: dist_real_def approx'_def approx_def d_def)
eberlm@62049
   525
  finally show ?thesis1 unfolding dist_real_def approx_def .
eberlm@62049
   526
  
eberlm@62049
   527
  from A have "ln 2 \<in> {approx' - d..approx' + d}" 
eberlm@62049
   528
    by (simp add: dist_real_def abs_real_def split: split_if_asm)
eberlm@62049
   529
  also have "\<dots> \<subseteq> {0.693147180549..0.693147180561}"
eberlm@62049
   530
    by (subst atLeastatMost_subset_iff, rule disjI2) (simp add: approx'_def d_def)
eberlm@62049
   531
  finally show ?thesis2 .
eberlm@62049
   532
qed
eberlm@62049
   533
eberlm@62049
   534
end
eberlm@62049
   535
eberlm@62049
   536
eberlm@62049
   537
lemma euler_mascheroni_bounds:
eberlm@62049
   538
  fixes n :: nat assumes "n \<ge> 1" defines "t \<equiv> harm n - ln (of_nat (Suc n)) :: real"
eberlm@62049
   539
  shows "euler_mascheroni \<in> {t + inverse (of_nat (2*(n+1)))..t + inverse (of_nat (2*n))}"
eberlm@62049
   540
  using assms euler_mascheroni_upper[of "n-1"] euler_mascheroni_lower[of "n-1"]
eberlm@62049
   541
  unfolding t_def by (cases n) (simp_all add: harm_Suc t_def inverse_eq_divide of_nat_mult)
eberlm@62049
   542
eberlm@62049
   543
lemma euler_mascheroni_bounds':
eberlm@62049
   544
  fixes n :: nat assumes "n \<ge> 1" "ln (real_of_nat (Suc n)) \<in> {l<..<u}"
eberlm@62049
   545
  shows "euler_mascheroni \<in> 
eberlm@62049
   546
           {harm n - u + inverse (of_nat (2*(n+1)))<..<harm n - l + inverse (of_nat (2*n))}"
eberlm@62049
   547
  using euler_mascheroni_bounds[OF assms(1)] assms(2) by auto
eberlm@62049
   548
eberlm@62049
   549
lemma euler_mascheroni_approx: 
eberlm@62049
   550
  defines "approx \<equiv> 0.577257 :: real" and "e \<equiv> 0.000063 :: real"
eberlm@62049
   551
  shows   "abs (euler_mascheroni - approx :: real) < e"
eberlm@62049
   552
  (is "abs (_ - ?approx) < ?e")
eberlm@62049
   553
proof -
eberlm@62049
   554
  def l \<equiv> "47388813395531028639296492901910937/82101866951584879688289000000000000 :: real"
eberlm@62049
   555
  def u \<equiv> "142196984054132045946501548559032969 / 246305600854754639064867000000000000 :: real"
eberlm@62049
   556
  have impI: "P \<longrightarrow> Q" if Q for P Q using that by blast
eberlm@62049
   557
  have hsum_63: "harm 63 = (310559566510213034489743057 / 65681493561267903750631200 ::real)"
eberlm@62049
   558
    by (simp add: harm_expand)
eberlm@62049
   559
  from harm_Suc[of 63] have hsum_64: "harm 64 = 
eberlm@62049
   560
          623171679694215690971693339 / (131362987122535807501262400::real)" 
eberlm@62049
   561
    by (subst (asm) hsum_63) simp
eberlm@62049
   562
  have "ln (64::real) = real (6::nat) * ln 2" by (subst ln_realpow[symmetric]) simp_all
eberlm@62049
   563
  hence "ln (real_of_nat (Suc 63)) \<in> {4.158883083293<..<4.158883083367}" using ln_2_bounds by simp
eberlm@62049
   564
  from euler_mascheroni_bounds'[OF _ this]
eberlm@62049
   565
    have "(euler_mascheroni :: real) \<in> {l<..<u}" 
eberlm@62049
   566
    by (simp add: hsum_63 del: greaterThanLessThan_iff) (simp only: l_def u_def)
eberlm@62049
   567
  also have "\<dots> \<subseteq> {approx - e<..<approx + e}"
eberlm@62049
   568
    by (subst greaterThanLessThan_subseteq_greaterThanLessThan, rule impI) 
eberlm@62049
   569
       (simp add: approx_def e_def u_def l_def)
eberlm@62049
   570
  finally show ?thesis by (simp add: abs_real_def)
eberlm@62049
   571
qed
eberlm@62049
   572
eberlm@62049
   573
end